/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \    /   O peration     | Website:  https://openfoam.org
    \  /    A nd           | Version:  9
     \/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

// Reduced case in two ways:
// 1. meshDensity 0.5: half the cells in every direction, so a cell that was
//    1 by 1 is now 2 by 2.
// 2. The domain was cut, xOutlet from 50 to 24 and yMax from 18 to 12 rInner.
//    simpleGrading is the ratio between the last cell of a block and its first,
//    so a shorter block needs a smaller ratio to keep the same growth per cell.

convertToMeters        1.0;
meshDensity  0.5;

rInner       0.5;                   // Cylinder radius
rOuter       #calc "2.5*$rInner";
xOutlet      #calc "24*$rInner";    // reference: 50, i.e. 25D behind the cylinder
xInlet       #calc "-10*$rInner";
yMax         #calc "12*$rInner";    // reference: 18, i.e. 9D above the cylinder
rInnerSinPi4 #calc "$rInner*sin(degToRad(45.0))";
rOuterSinPi4 #calc "$rOuter*sin(degToRad(45.0))";

thetaJet     10.0;                  //degrees
xJetInner    #calc "$rInner*sin(degToRad(0.5*$thetaJet))";
yJetInner    #calc "$rInner*cos(degToRad(0.5*$thetaJet))";
xJetOuter    #calc "$rOuter*sin(degToRad(0.5*$thetaJet))";
yJetOuter    #calc "$rOuter*cos(degToRad(0.5*$thetaJet))";

zmin         -0.5;
zmax          0.5;

nRadial      #calc "round(23 * $meshDensity)";
nTheta1      #calc "round(7  * $meshDensity)";
nTheta2      #calc "round(11 * $meshDensity)";
nJet         #calc "round(6  * $meshDensity)";
nxoutlet     #calc "round(38 * $meshDensity)";  // reference: 52 over 25D
nxinlet      #calc "round(14 * $meshDensity)";
ny           #calc "round(19 * $meshDensity)";  // reference: 23 over 9D
nz           1;

vertices #codeStream
{
    codeInclude
    #{
        #include "pointField.H"
    #};

    code
    #{
        pointField points(22);
        points[0]  = point($rInner, 0, $zmin);
        points[1]  = point($rOuter, 0, $zmin);
        points[2]  = point($rInnerSinPi4, $rInnerSinPi4, $zmin);
        points[3]  = point($rOuterSinPi4, $rOuterSinPi4, $zmin);
        points[4]  = point($xJetInner, $yJetInner, $zmin);
        points[5]  = point($xJetOuter, $yJetOuter, $zmin);
        points[6]  = point(-$xJetInner, $yJetInner, $zmin);
        points[7]  = point(-$xJetOuter, $yJetOuter, $zmin);
        points[8]  = point(-$rInnerSinPi4, $rInnerSinPi4, $zmin);
        points[9]  = point(-$rOuterSinPi4, $rOuterSinPi4, $zmin);
        points[10] = point(-$rInner, 0, $zmin);
        points[11] = point(-$rOuter, 0, $zmin);
        points[12] = point($xOutlet, 0, $zmin);
        points[13] = point($xOutlet, $rOuterSinPi4, $zmin);
        points[14] = point($xOutlet, $yMax, $zmin);
        points[15] = point($rOuterSinPi4, $yMax, $zmin);
        points[16] = point($xJetOuter, $yMax, $zmin);
        points[17] = point(-$xJetOuter, $yMax, $zmin);
        points[18] = point(-$rOuterSinPi4, $yMax, $zmin);
        points[19] = point($xInlet, $yMax, $zmin);
        points[20] = point($xInlet, $rOuterSinPi4, $zmin);
        points[21] = point($xInlet, 0, $zmin);
     
        // Mirror -z points to +z side
        label sz = points.size();
        points.setSize(2*sz);
        for (label i = 0; i < sz; ++i)
        {
            const point& pt = points[i];
            points[i + sz] = point(pt.x(), pt.y(), $zmax);
        }

        os  << points;
    #};
};


blocks
(
    hex ( 0 1 3 2 22 23 25 24 ) ($nRadial $nTheta1 $nz) simpleGrading (18 1 1)
    hex ( 2 3 5 4 24 25 27 26 ) ($nRadial $nTheta2 $nz) simpleGrading (18 0.25 1)
    hex ( 4 5 7 6 26 27 29 28 ) ($nRadial $nJet $nz) simpleGrading (18 1 1)
    hex ( 6 7 9 8 28 29 31 30 ) ($nRadial $nTheta2 $nz) simpleGrading (18 4 1)
    hex ( 8 9 11 10 30 31 33 32 ) ($nRadial $nTheta1 $nz) simpleGrading (18 1 1)
    hex ( 1 12 13 3 23 34 35 25 ) ($nxoutlet $nTheta1 $nz) simpleGrading (6 1 1)
    hex ( 3 13 14 15 25 35 36 37 ) ($nxoutlet $ny $nz) simpleGrading (6 5.5 1)
    hex ( 5 3 15 16 27 25 37 38 ) ($nTheta2 $ny $nz) simpleGrading (4 5.5 1)
    hex ( 7 5 16 17 29 27 38 39 ) ($nJet $ny $nz) simpleGrading (1 5.5 1)
    hex ( 9 7 17 18 31 29 39 40 ) ($nTheta2 $ny $nz) simpleGrading (0.25 5.5 1)
    hex ( 20 9 18 19 42 31 40 41 ) ($nxinlet $ny $nz) simpleGrading (0.2 5.5 1)
    hex ( 21 11 9 20 43 33 31 42 ) ($nxinlet $nTheta1 $nz) simpleGrading (0.2 1 1)

);

edges #codeStream
{
    codeInclude
    #{
        #include "unitConversion.H"
    #};

    code
    #{

        const scalar half = 0.5*$thetaJet;

        // vertex angles around the upper half, in degrees, in vertex order
        const scalar angle[6] = {0, 45, 90 - half, 90 + half, 135, 180};

        const scalar radius[2] = {$rInner, $rOuter};
        const scalar zPlane[2] = {$zmin, $zmax};

        const label nLower = 22;    // vertices before the mirrored copies

        os << "(" << nl;

        for (label ring = 0; ring < 2; ++ring)
        {
            for (label side = 0; side < 2; ++side)
            {
                const label first = ring + side*nLower;

                for (label i = 0; i < 5; ++i)
                {
                    const scalar mid = degToRad(0.5*(angle[i] + angle[i + 1]));

                    os  << "    arc " << first + 2*i
                        << " " << first + 2*(i + 1)
                        << " (" << radius[ring]*cos(mid)
                        << " " << radius[ring]*sin(mid)
                        << " " << zPlane[side] << ")" << nl;
                }
            }
        }

        os << ")" << nl;
    #};
};

boundary
(
    cylinder
    {
        type wall;
        faces
        (
            (0 2 22 24)
            (2 4 24 26)
            (6 8 28 30)
            (8 10 30 32)
        );
    }
    jet
    {
        type patch;
        faces
        (
            (4 6 26 28)
        );
    }    
    inlet
    {
        type patch;
        faces
        (
            (19 20 41 42)
            (20 21 42 43)
        );
    }
    outlet
    {
        type patch;
        faces
        (
            (12 13 34 35)
            (13 14 35 36)
        );
    }
    topAndBottom
    {
        type patch;
        faces
        (
            (14 15 36 37)
            (15 16 37 38)
            (16 17 38 39)
            (17 18 39 40)
            (18 19 40 41)
        );
    }
    symmetry
    {
        type symmetryPlane;
        faces
        (
            (0 1 22 23)
            (1 12 23 34)
            (10 11 32 33)
            (11 21 33 43)            
        );
    }    
    frontAndBack
    {
        type empty;
        faces
        (
            (0 1 3 2)
            (2 3 5 4)
            (4 5 7 6)
            (6 7 9 8)
            (8 9 11 10)
            (1 12 13 3)
            (3 13 14 15)
            (3 15 16 5)
            (5 16 17 7)
            (7 17 18 9)
            (9 18 19 20)
            (11 9 20 21)
            (22 23 25 24)
            (24 25 27 26)
            (26 27 29 28)
            (28 29 31 30)
            (30 31 33 32)
            (23 34 35 25)
            (25 35 36 37)
            (25 37 38 27)
            (27 38 39 29)
            (29 39 40 31)
            (31 40 41 42)
            (33 31 42 43)
        );
    }
);

// ************************************************************************* //
