openfoam/src/mesh/blockMesh/PDRblockMesh/PDRblockBlockMesh.C

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2020-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "PDRblock.H"
#include "ListOps.H"
#include "gradingDescriptors.H"
#include "objectRegistry.H"
#include "Time.H"
#include "IOdictionary.H"
#include "Fstream.H"
#include "OTstream.H"
#include "edgeHashes.H"

#include "cellModel.H"
#include "blockMesh.H"
#include "polyMesh.H"
#include "searchableSphere.H"

// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //

namespace Foam
{

// Name for the projection geometry
static const word projKeyword("project");

// Name for the projection geometry
static const word projGeomName("sphere");


//- Calculate geometric ratio from relative ratio
inline scalar relativeToGeometricRatio
(
    const scalar expRatio,
    const label nDiv
)
{
    return nDiv > 1 ? pow(expRatio, (nDiv - 1)) : 1.0;
}


// Output space-separated flat list. No size prefix.
template<class T>
static Ostream& outputFlatList(Ostream& os, const UList<T>& list)
{
    os  << token::BEGIN_LIST;
    label i = 0;
    for (const label val : list)
    {
        if (i++) os << token::SPACE;
        os  << val;
    }
    os  << token::END_LIST;

    return os;
}


// Begin indent list
static inline Ostream& begIndentList(Ostream& os)
{
    os  << indent << incrIndent << token::BEGIN_LIST << nl;
    return os;
}

// End indent list
static inline Ostream& endIndentList(Ostream& os)
{
    os  << decrIndent << indent << token::END_LIST;
    return os;
}


// Output list contents (newline separated) indented.
template<class T>
static Ostream& outputIndent(Ostream& os, const UList<T>& list)
{
    for (const T& val : list)
    {
        os  << indent << val << nl;
    }
    return os;
}


//
static Ostream& serializeHex
(
    Ostream& os,
    const labelUList& hexVerts,
    const labelVector& hexCount,
    const Vector<gradingDescriptors> hexGrade,
    const word& zoneName = word::null
)
{
    os  << indent << cellModel::modelNames[cellModel::HEX] << token::SPACE;
    outputFlatList(os, hexVerts);

    if (!zoneName.empty())
    {
        os  << token::SPACE << zoneName;
    }

    os  << token::SPACE << hexCount << nl
        << indent << word("edgeGrading") << nl;

    begIndentList(os);

    // Grading (x/y/z)
    for (const gradingDescriptors& gds : hexGrade)
    {
        begIndentList(os);
        outputIndent(os, gds);
        endIndentList(os) << nl;
    }

    endIndentList(os) << nl;
    return os;
}


// Generate list with entries:
//
//    project (x y z) (geometry)
//
static Ostream& serializeProjectPoints
(
    Ostream& os,
    const UList<point>& list
)
{
    for (const point& p : list)
    {
        os  << indent << projKeyword << token::SPACE
            << p
            << token::SPACE
            << token::BEGIN_LIST << projGeomName << token::END_LIST << nl;
    }

    return os;
}


// Generate entry:
//
//    project (beg end) (geometry)
//
static Ostream& serializeProjectEdge
(
    Ostream& os,
    const edge& e
)
{
    os  << indent << projKeyword << token::SPACE;

    if (e.sorted())
    {
        os  << e.first() << token::SPACE << e.second();
    }
    else
    {
        os  << e.second() << token::SPACE << e.first();
    }

    os  << token::SPACE
        << token::BEGIN_LIST << projGeomName << token::END_LIST << nl;

    return os;
}


// Generate entry:
//
//    (0 1 2 ..)
//
static Ostream& serializeFace
(
    Ostream& os,
    const face& list
)
{
    os  << indent;
    outputFlatList(os, list);
    os  << nl;
    return os;
}


// Generate entry:
//
//    project (0 1 2 ..) geometry
//
static Ostream& serializeProjectFace
(
    Ostream& os,
    const face& list
)
{
    os  << indent << projKeyword << token::SPACE;
    outputFlatList(os, list);
    os  << token::SPACE << projGeomName << nl;

    return os;
}

} // namespace Foam


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

Foam::Ostream& Foam::PDRblock::blockMeshDict
(
    Ostream& os,
    const bool withHeader
) const
{
    if (withHeader)
    {
        // Use dummy time for fake objectRegistry
        autoPtr<Time> dummyTimePtr(Time::New());

        IOdictionary iodict
        (
            IOobject
            (
                "blockMeshDict",
                dummyTimePtr->system(),
                *dummyTimePtr,
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false  // no register
            )
        );

        iodict.writeHeader(os);
    }

    const cellModel& hex = cellModel::ref(cellModel::HEX);

    // The mesh topology will normally be an O-grid with a central (inner)
    // block and 6 outer blocks.
    // The inner block is described by (0 1 2 3 4 5 6 7).
    // The additional points for the outer region: (8 9 19 11 12 13 14 15)

    // The outer blocks will be addressed according to their
    // placement w.r.t. the inner block, and defined such that they retain
    // the same global orientation as the inner block.
    // For example, face 0 of all blocks will be on the logical x-min side
    // of the mesh.

    // List of hex vertices, ordered with outer blocks first to allow
    // direct addressing by their logical position.
    const FixedList<labelList, 8>
        hexVerts
        ({
            {8, 0,  3, 11, 12,  4,  7, 15},     // x-min block (face 0)
            {1, 9, 10,  2,  5, 13, 14,  6},     // x-max block (face 1)
            {8, 9,  1,  0, 12, 13,  5,  4},     // y-min block (face 2)
            {3, 2, 10, 11,  7,  6, 14, 15},     // y-max block (face 3)
            {8, 9, 10, 11,  0,  1,  2,  3},     // z-min block (face 4)
            {4, 5,  6,  7, 12, 13, 14, 15},     // z-max block (face 5)
            {0, 1,  2,  3,  4,  5,  6,  7},     // Inner box description
            {8, 9, 10, 11, 12, 13, 14, 15},     // Outer box description
        });

    // The face or the logical block index (for the O-grid)
    enum faceIndex
    {
        X_Min = 0,
        X_Max = 1,
        Y_Min = 2,
        Y_Max = 3,
        Z_Min = 4,
        Z_Max = 5,
        Inner_Block = 6,    // Inner block description
        Outer_Block = 7,    // Outer bounding box description
    };

    // Lists of block/face for outside and ground faces
    DynamicList<labelPair> outerFaces(8);
    DynamicList<labelPair> groundFaces(8);

    // We handle a few fixed topology configurations

    enum blockTopologyType
    {
        INNER_ONLY = 0,   // No outer region
        EXTENDED   = 1,   // Outer created by extending inner region
        CLIP_BOTTOM = 5,  // Outer O-grid on 5 sides
        FULL_OUTER = 6    // Outer O-grid on 6 sides
    };



    // Expansion ratios need conversion from relative to geometric
    const bool useRelToGeom =
        (expansionType::EXPAND_RATIO == outer_.expandType_);


    // Physical dimensions

    Vector<gridControl> ctrl(control_);
    boundBox innerCorners(bounds(ctrl.x(), ctrl.y(), ctrl.z()));

    boundBox outerCorners;


    point radialCentre(innerCorners.centre());
    vector radialSizes(0.5*innerCorners.span());

    blockTopologyType outerTopology = INNER_ONLY;


    if (outer_.active())
    {
        outerTopology = FULL_OUTER;

        // Convert from relative size
        radialSizes.x() *= outer_.relSize_.x();
        radialSizes.y() *= outer_.relSize_.y();
        radialSizes.z() *= min(outer_.relSize_.x(), outer_.relSize_.y());

        if (outer_.onGround())
        {
            outerTopology = CLIP_BOTTOM;
            radialCentre.z() = innerCorners.min().z();
            radialSizes.z() *= 2;
        }

        // Corners for a box
        outerCorners.min() = radialCentre - radialSizes;
        outerCorners.max() = radialCentre + radialSizes;

        if (outer_.onGround())
        {
            outerCorners.min().z() = innerCorners.min().z();
        }


        if (outer_.isSphere())
        {
            // For spheroid projection, don't trust that blockMesh does it
            // properly. Give some reasonable estimates of the corners.

            // Use dummy Time for objectRegistry
            autoPtr<Time> dummyTimePtr(Time::New());

            const IOobject io
            (
                "sphere",
                *dummyTimePtr,
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false  // do not register
            );

            searchableSphere sphere(io, radialCentre, radialSizes);

            pointField queries(2);
            queries[0] = outerCorners.min();
            queries[1] = outerCorners.max();

            List<pointIndexHit> hits;
            sphere.findNearest
            (
                queries,
                scalarField(2, GREAT),
                hits
            );

            outerCorners.min() = hits[0].hitPoint();
            outerCorners.max() = hits[1].hitPoint();
        }
        else if (outerControl::OUTER_EXTEND == outer_.type_)
        {
            outerTopology = EXTENDED;

            // Extend the inner block
            label outerCount;
            scalar expRatio;

            outerCount = outer_.nCells_.x();
            expRatio = outer_.expansion_.x();
            if (useRelToGeom)
            {
                expRatio = relativeToGeometricRatio(expRatio, outerCount);
            }

            ctrl.x().prepend(outerCorners.min().x(), outerCount, -expRatio);
            ctrl.x().append(outerCorners.max().x(), outerCount, expRatio);


            outerCount = outer_.nCells_.y();
            expRatio = outer_.expansion_.y();
            if (useRelToGeom)
            {
                expRatio = relativeToGeometricRatio(expRatio, outerCount);
            }

            ctrl.y().prepend(outerCorners.min().y(), outerCount, -expRatio);
            ctrl.y().append(outerCorners.max().y(), outerCount, expRatio);

            outerCount = max(outer_.nCells_.x(), outer_.nCells_.y());
            expRatio = min(outer_.expansion_.x(), outer_.expansion_.y());
            if (useRelToGeom)
            {
                expRatio = relativeToGeometricRatio(expRatio, outerCount);
            }

            if (!outer_.onGround())
            {
                ctrl.z().prepend(outerCorners.min().z(), outerCount, -expRatio);
            }
            ctrl.z().append(outerCorners.max().z(), outerCount, expRatio);

            // Update corners
            innerCorners = bounds(ctrl.x(), ctrl.y(), ctrl.z());
            outerCorners = innerCorners;
        }
    }


    const Vector<gradingDescriptors> innerGrading(grading(ctrl));
    const labelVector innerCount(sizes(ctrl));

    labelVector hexCount;
    Vector<gradingDescriptors> hexGrade;


    const label radialCount = outer_.nCells_.x();
    scalar expRatio = outer_.expansion_.x();

    if (useRelToGeom)
    {
        expRatio = relativeToGeometricRatio(expRatio, radialCount);
    }

    const gradingDescriptors radialInward
    (
        gradingDescriptor{-expRatio}
    );

    const gradingDescriptors radialOutward
    (
        gradingDescriptor{expRatio}
    );


    if (EXTENDED == outerTopology)
    {
        // The inner block is extended to become the outer faces
        outerFaces.append
        ({
            labelPair(Inner_Block, X_Min),
            labelPair(Inner_Block, X_Max),
            labelPair(Inner_Block, Y_Min),
            labelPair(Inner_Block, Y_Max),
            labelPair(Inner_Block, Z_Max)
        });

        // The ground faces vs outside faces
        if (outer_.onGround())
        {
            groundFaces.append
            (
                labelPair(Inner_Block, Z_Min)
            );
        }
        else
        {
            outerFaces.append
            (
                labelPair(Inner_Block, Z_Min)
            );
        }
    }
    else if (CLIP_BOTTOM == outerTopology || FULL_OUTER == outerTopology)
    {
        // The outside faces
        outerFaces.append
        ({
            labelPair(X_Min, X_Min),
            labelPair(X_Max, X_Max),
            labelPair(Y_Min, Y_Min),
            labelPair(Y_Max, Y_Max),
            labelPair(Z_Max, Z_Max)
        });

        // The ground faces
        if (CLIP_BOTTOM == outerTopology)
        {
            groundFaces.append
            ({
                labelPair(X_Min, Z_Min),
                labelPair(X_Max, Z_Min),
                labelPair(Y_Min, Z_Min),
                labelPair(Y_Max, Z_Min),
                // Note: {Z_Min, Z_Min} will not exist
                labelPair(Inner_Block, Z_Min)
            });
        }
        else
        {
            outerFaces.append
            (
                labelPair(Z_Min, Z_Min)
            );
        }
    }


    if (outer_.isSphere())
    {
        os.beginBlock("geometry");
        {
            os.beginBlock(projGeomName);
            {
                os.writeEntry("type", "sphere");
                os.writeEntry("origin", radialCentre);
                os.writeEntry("radius", radialSizes);
            }
            os.endBlock();
        }
        os.endBlock();
    }

    // vertices
    {
        os  << nl << word("vertices") << nl;
        begIndentList(os);

        pointField corners(innerCorners.points());

        // inner
        outputIndent(os, corners);

        // outer
        if (CLIP_BOTTOM == outerTopology || FULL_OUTER == outerTopology)
        {
            corners = outerCorners.points();

            if (outer_.isSphere())
            {
                serializeProjectPoints(os, corners);
            }
            else
            {
                outputIndent(os, corners);
            }
        }

        endIndentList(os) << token::END_STATEMENT << nl;
    }


    // blocks
    {
        word innerZoneName = "inner";
        if (INNER_ONLY == outerTopology || EXTENDED == outerTopology)
        {
            innerZoneName.clear();
        }

        os  << nl << word("blocks") << nl;
        begIndentList(os);

        // Inner block
        hexCount = innerCount;

        serializeHex
        (
            os,
            hexVerts[Inner_Block],
            hexCount,
            innerGrading,
            innerZoneName
        );

        // outer
        if (CLIP_BOTTOM == outerTopology || FULL_OUTER == outerTopology)
        {
            // Radial direction = X
            hexCount = innerCount;
            hexGrade = innerGrading;

            hexCount.x() = radialCount;

            // Face 0: x-min
            {
                hexGrade.x() = radialInward;
                serializeHex(os, hexVerts[X_Min], hexCount, hexGrade);
            }
            // Face 1: x-max
            {
                hexGrade.x() = radialOutward;
                serializeHex(os, hexVerts[X_Max], hexCount, hexGrade);
            }


            // Radial direction = Y
            hexCount = innerCount;
            hexGrade = innerGrading;

            hexCount.y() = radialCount;

            // Face 2: y-min
            {
                hexGrade.y() = radialInward;
                serializeHex(os, hexVerts[Y_Min], hexCount, hexGrade);
            }
            // Face 3: y-max
            {
                hexGrade.y() = radialOutward;
                serializeHex(os, hexVerts[Y_Max], hexCount, hexGrade);
            }


            // Radial direction = Z
            hexCount = innerCount;
            hexGrade = innerGrading;

            hexCount.z() = radialCount;

            // Face 4: z-min
            if (!outer_.onGround())
            {
                hexGrade.z() = radialInward;
                serializeHex(os, hexVerts[Z_Min], hexCount, hexGrade);
            }
            // Face 5: z-max
            {
                hexGrade.z() = radialOutward;
                serializeHex(os, hexVerts[Z_Max], hexCount, hexGrade);
            }
        }

        endIndentList(os) << token::END_STATEMENT << nl;
    }


    // edges
    {
        os  << nl << word("edges") << nl;
        begIndentList(os);

        if (outer_.isSphere() && outerFaces.size())
        {
            // Edges for the outer face of the block
            edgeHashSet projEdges(32);

            for (const labelPair& pr : outerFaces)
            {
                projEdges.insert
                (
                    hex.face(pr.second(), hexVerts[pr.first()]).edges()
                );
            }

            for (const edge& e : projEdges.sortedToc())
            {
                serializeProjectEdge(os, e);
            }
        }

        endIndentList(os) << token::END_STATEMENT << nl;
    }


    // faces
    {
        os  << nl << word("faces") << nl;
        begIndentList(os);

        if (outer_.isSphere() && outerFaces.size())
        {
            for (const labelPair& pr : outerFaces)
            {
                serializeProjectFace
                (
                    os,
                    hex.face(pr.second(), hexVerts[pr.first()])
                );
            }
        }

        endIndentList(os) << token::END_STATEMENT << nl;
    }


    // boundary
    {
        os  << nl << word("boundary") << nl;
        begIndentList(os);

        // outer
        {
            os.beginBlock("outer");
            os.writeEntry("type", word("patch"));

            os  << indent << word("faces") << nl;
            begIndentList(os);

            for (const labelPair& pr : outerFaces)
            {
                serializeFace
                (
                    os,
                    hex.face(pr.second(), hexVerts[pr.first()])
                );
            }

            endIndentList(os) << token::END_STATEMENT << nl;

            os.endBlock();
        }

        if (outer_.onGround())
        {
            os.beginBlock("ground");
            os.writeEntry("type", word("wall"));

            os  << indent << word("faces") << nl;
            begIndentList(os);

            for (const labelPair& pr : groundFaces)
            {
                serializeFace
                (
                    os,
                    hex.face(pr.second(), hexVerts[pr.first()])
                );
            }

            endIndentList(os) << token::END_STATEMENT << nl;
            os.endBlock();
        }

        endIndentList(os) << token::END_STATEMENT << nl;
    }


    if (withHeader)
    {
        IOobject::writeEndDivider(os);
    }

    return os;
}


Foam::dictionary Foam::PDRblock::blockMeshDict() const
{
    OTstream os;
    blockMeshDict(os);

    ITstream is;
    is.transfer(os.tokens());

    return dictionary(is);
}


void Foam::PDRblock::writeBlockMeshDict(const IOobject& io) const
{
    IOdictionary iodict
    (
        IOobject
        (
            io.name(),
            io.db().time().system(),
            io.local(),
            io.db(),
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            false  // no register
        )
    );

    OFstream os(iodict.objectPath());

    Info<< nl
        << "Generate blockMeshDict: "
        << iodict.db().time().relativePath(os.name()) << endl;

    // Set precision for points to 10
    os.precision(max(10u, IOstream::defaultPrecision()));

    iodict.writeHeader(os);

    // Just like writeData, but without copying beforehand
    this->blockMeshDict(os);

    IOobject::writeEndDivider(os);
}


Foam::autoPtr<Foam::blockMesh>
Foam::PDRblock::createBlockMesh(const IOobject& io) const
{
    IOdictionary iodict
    (
        IOobject
        (
            "blockMeshDict.PDRblockMesh",
            io.db().time().system(),
            io.local(),
            io.db(),
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            false  // no register
        ),
        blockMeshDict()
    );

    return autoPtr<blockMesh>::New(iodict);
}


Foam::autoPtr<Foam::polyMesh>
Foam::PDRblock::meshBlockMesh(const IOobject& io) const
{
    const bool oldVerbose = blockMesh::verboseOutput;
    blockMesh::verboseOutput = false;

    autoPtr<polyMesh> meshPtr(createBlockMesh(io)->mesh(io));

    blockMesh::verboseOutput = oldVerbose;

    // This is a bit ugly.
    // For extend, we still wish to have an 'inner' cellZone,
    // but we meshed the entirety.

    if
    (
        outerControl::OUTER_EXTEND == outer_.type_
     && meshPtr->cellZones().empty()
    )
    {
        const boundBox innerBox
        (
            bounds(control_.x(), control_.y(), control_.z())
        );

        const label nZoneCellsMax =
        (
            control_.x().nCells()
          * control_.y().nCells()
          * control_.z().nCells()
        );


        polyMesh& pmesh = *meshPtr;

        List<cellZone*> cz(1);
        cz[0] = new cellZone
        (
            "inner",
            labelList(nZoneCellsMax),
            0, // zonei
            pmesh.cellZones()
        );

        cellZone& innerZone = *(cz[0]);

        const vectorField& cc = pmesh.cellCentres();

        label nZoneCells = 0;

        for
        (
            label celli = 0;
            celli < cc.size() && nZoneCells < nZoneCellsMax;
            ++celli
        )
        {
            if (innerBox.contains(cc[celli]))
            {
                innerZone[nZoneCells] = celli;
                ++nZoneCells;
            }
        }

        innerZone.resize(nZoneCells);

        pmesh.pointZones().clear();
        pmesh.faceZones().clear();
        pmesh.cellZones().clear();
        pmesh.addZones(List<pointZone*>(), List<faceZone*>(), cz);
    }

    return meshPtr;
}


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