PDRobstacleIO.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2019-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 "PDRsetFields.H"
#include "PDRobstacle.H"
#include "volumeType.H"

using namespace Foam;
using namespace Foam::constant;

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

bool Foam::PDRobstacle::read(Istream& is)
{
    this->clear();

    const word obsType(is);
    const dictionary dict(is);

    auto* mfuncPtr = readdictionaryMemberFunctionTable(obsType);

    if (!mfuncPtr)
    {
        FatalIOErrorInLookup
        (
            is,
            "obstacle",
            obsType,
            *readdictionaryMemberFunctionTablePtr_
        ) << exit(FatalIOError);
    }

    mfuncPtr(*this, dict);

    return true;
}


Foam::scalar Foam::PDRobstacle::readFiles
(
    const fileName& obsFileDir,
    const wordList& obsFileNames,
    const boundBox& meshBb,

    DynamicList<PDRobstacle>& blocks,
    DynamicList<PDRobstacle>& cylinders
)
{
    blocks.clear();
    cylinders.clear();

    scalar totVolume = 0;
    label nOutside = 0;
    label nProtruding = 0;

    scalar shift = pars.obs_expand;

    if (!obsFileNames.empty())
    {
        Info<< "Reading obstacle files" << nl;
    }

    label maxGroup = -1;

    for (const word& inputFile : obsFileNames)
    {
        Info<< "    file: " << inputFile << nl;

        fileName path = (obsFileDir / inputFile);

        IFstream is(path);
        dictionary inputDict(is);

        const scalar scaleFactor = inputDict.getOrDefault<scalar>("scale", 0);

        const label verbose = inputDict.getOrDefault<label>("verbose", 0);

        for (const entry& dEntry : inputDict)
        {
            if (!dEntry.isDict())
            {
                // ignore non-dictionary entry
                continue;
            }

            const dictionary& dict = dEntry.dict();

            if (!dict.getOrDefault("enabled", true))
            {
                continue;
            }

            label obsGroupId = 0;
            if (dict.readIfPresent("groupId", obsGroupId))
            {
                maxGroup = max(maxGroup, obsGroupId);
            }
            else
            {
                obsGroupId = ++maxGroup;
            }


            pointField pts;
            dict.readIfPresent("locations", pts);
            if (pts.empty())
            {
                pts.resize(1, Zero);
            }

            List<PDRobstacle> obsInput;
            dict.readEntry("obstacles", obsInput);

            label nCyl = 0; // The number of cylinders vs blocks

            for (PDRobstacle& obs : obsInput)
            {
                obs.groupId = obsGroupId;
                obs.scale(scaleFactor);

                if (obs.isCylinder())
                {
                    ++nCyl;
                }
            }

            const label nBlock = (obsInput.size() - nCyl);

            blocks.reserve(blocks.size() + nBlock*pts.size());
            cylinders.reserve(cylinders.size() + nCyl*pts.size());

            if (verbose)
            {
                Info<< "Read " << obsInput.size() << " obstacles ("
                    << nCyl << " cylinders) with "
                    << pts.size() << " locations" << nl;

                if (verbose > 1)
                {
                    Info<< "locations " << pts << nl
                        << "obstacles " << obsInput << nl;
                }
            }

            for (const PDRobstacle& scanObs : obsInput)
            {
                // Reject anything below minimum width
                if (scanObs.tooSmall(pars.min_width))
                {
                    continue;
                }

                for (const point& origin : pts)
                {
                    // A different (very small) shift for each obstacle
                    // so that faces cannot be coincident

                    shift += floatSMALL;
                    const scalar shift2 = shift * 2.0;


                switch (scanObs.typeId)
                {
                    case PDRobstacle::CYLINDER:
                    {
                        // Make a copy
                        PDRobstacle obs(scanObs);

                        // Offset for the group position
                        obs.pt += origin;

                        // Shift the end outwards so, if exactly on
                        // cell boundary, now overlap cell.
                        // So included in Aw.
                        obs.pt -= point::uniform(shift);
                        obs.len() += shift2;
                        obs.dia() -= floatSMALL;


                        // Trim against the mesh bounds.
                        // Ignore if it doesn't overlap, or bounds are invalid
                        const volumeType vt = obs.trim(meshBb);

                        switch (vt)
                        {
                            case volumeType::OUTSIDE:
                                ++nOutside;
                                continue; // Can ignore the rest
                                break;

                            case volumeType::MIXED:
                                ++nProtruding;
                                break;

                            default:
                                break;
                        }

                        // Later for position sorting
                        switch (obs.orient)
                        {
                            case vector::X:
                                obs.sortBias = obs.len();
                                break;
                            case vector::Y:
                                obs.sortBias = 0.5*obs.dia();
                                break;
                            case vector::Z:
                                obs.sortBias = 0.5*obs.dia();
                                break;
                        }

                        totVolume += obs.volume();
                        cylinders.append(obs);

                        break;
                    }

                    case PDRobstacle::DIAG_BEAM:
                    {
                        // Make a copy
                        PDRobstacle obs(scanObs);

                        // Offset for the group position
                        obs.pt += origin;

                        // Shift the end outwards so, if exactly on
                        // cell boundary, now overlap cell.
                        // So included in Aw.
                        obs.pt -= point::uniform(shift);
                        obs.len() += shift2;
                        obs.wa += shift2;
                        obs.wb += shift2;

                        totVolume += obs.volume();
                        cylinders.append(obs);

                        break;
                    }

                    case PDRobstacle::CUBOID_1:
                    case PDRobstacle::LOUVRE_BLOWOFF:
                    case PDRobstacle::CUBOID:
                    case PDRobstacle::WALL_BEAM:
                    case PDRobstacle::GRATING:
                    case PDRobstacle::RECT_PATCH:
                    {
                        // Make a copy
                        PDRobstacle obs(scanObs);

                        // Offset for the position of the group
                        obs.pt += origin;

                        if (obs.typeId == PDRobstacle::GRATING)
                        {
                            if (obs.slat_width <= 0)
                            {
                                obs.slat_width = pars.def_grating_slat_w;
                            }
                        }

                        // Shift the end outwards so, if exactly on
                        // cell boundary, now overlap cell.
                        // So included in Aw.
                        obs.pt -= point::uniform(shift);
                        obs.span += point::uniform(shift2);


                        // Trim against the mesh bounds.
                        // Ignore if it doesn't overlap, or bounds are invalid
                        const volumeType vt = obs.trim(meshBb);

                        switch (vt)
                        {
                            case volumeType::OUTSIDE:
                                ++nOutside;
                                continue; // Can ignore the rest
                                break;

                            case volumeType::MIXED:
                                ++nProtruding;
                                break;

                            default:
                                break;
                        }

                        totVolume += obs.volume();

                        blocks.append(obs);

                        break;
                    }
                }
            }
            }

            // Info<< "Cylinders: " << cylinders << nl;
        }

        if (nOutside || nProtruding)
        {
            Info<< "Warning: " << nOutside << " obstacles outside domain, "
                << nProtruding << " obstacles partly outside domain" << nl;
        }
    }

    // #ifdef FULLDEBUG
    // Info<< blocks << nl << cylinders << nl;
    // #endif


    Info<< "Number of obstacles: "
        << (blocks.size() + cylinders.size()) << " ("
        << cylinders.size() << " cylinders)" << nl;

    return totVolume;
}


// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //

Foam::Istream& Foam::operator>>(Istream& is, PDRobstacle& obs)
{
    obs.read(is);

    return is;
}


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