PDRobstacleLegacyRead.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) 2016 Shell Research Ltd.
    Copyright (C) 2019-2020 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;

#undef USE_ZERO_INSTANCE_GROUPS
// #define USE_ZERO_INSTANCE_GROUPS


// Counting
Foam::labelPair Foam::PDRlegacy::readObstacleFiles
(
    const fileName& obsFileDir,
    const wordList& obsFileNames,
    Map<obstacleGrouping>& groups
)
{
    // Default group with single instance and position (0,0,0)
    groups(0).clear();
    groups(0).append(point::zero);

    string buffer;

    if (!obsFileNames.empty())
    {
        Info<< "Counting groups in obstacle files" << nl;
    }
    for (const word& inputFile : obsFileNames)
    {
        Info<< "    file: " << inputFile << nl;

        fileName path = (obsFileDir / inputFile);

        IFstream is(path);

        while (is.good())
        {
            // Process each line of obstacle files
            is.getLine(buffer);

            const auto firstCh = buffer.find_first_not_of(" \t\n\v\f\r");

            if
            (
                firstCh == std::string::npos
             || buffer[firstCh] == '#'
            )
            {
                // Empty line or comment line
                continue;
            }

            int typeId;
            double x, y, z; // Double (not scalar) to match sscanf spec

            if
            (
                sscanf(buffer.c_str(), "%d %lf %lf %lf", &typeId, &x, &y, &z)<4
             || typeId == 0
             || typeId == PDRobstacle::MESH_PLANE
            )
            {
                continue;
            }

            x *= pars.scale;
            y *= pars.scale;
            z *= pars.scale;

            const label groupId = typeId / 100;
            typeId %= 100;

            if (typeId == PDRobstacle::OLD_INLET)
            {
                Info<< "Ignored old-inlet type" << nl;
                continue;
            }
            else if (typeId == PDRobstacle::GRATING && pars.ignoreGratings)
            {
                Info<< "Ignored grating" << nl;
                continue;
            }

            if (typeId == 0)
            {
                // Defining a group location
                groups(groupId).append(x, y, z);
            }
            else if (PDRobstacle::isCylinder(typeId))
            {
                // Increment cylinder count for the group
                groups(groupId).addCylinder();
            }
            else
            {
                // Increment obstacle count for the group
                groups(groupId).addObstacle();
            }
        }
    }


    label nTotalObs = 0;
    label nTotalCyl = 0;

    label nMissedObs = 0;
    label nMissedCyl = 0;

    forAllConstIters(groups, iter)
    {
        const auto& group = iter.val();

        nTotalObs += group.nTotalObstacle();
        nTotalCyl += group.nTotalCylinder();

        if (group.empty())
        {
            nMissedObs += group.nObstacle();
            nMissedCyl += group.nCylinder();
        }
    }

    for (const label groupId : groups.sortedToc())
    {
        const auto& group = groups[groupId];

        if (groupId)
        {
            if (group.size())
            {
                Info<< "Found " << group.size()
                    << " instances of group " << groupId << " ("
                    << group.nObstacle() << " obstacles "
                    << group.nCylinder() << " cylinders)"
                    << nl;
            }
        }
        else
        {
            // The group 0 is for ungrouped obstacles
            Info<< "Found "
                << group.nObstacle() << " obstacles "
                << group.nCylinder() << " cylinders not in groups" << nl;
        }
    }

    Info<< "Number of obstacles: "
        << (nTotalObs + nTotalCyl) << " ("
        << nTotalCyl << " cylinders)" << nl;

    if (nMissedObs + nMissedCyl)
    {
        #ifdef USE_ZERO_INSTANCE_GROUPS

        nTotalObs += nMissedObs;
        nTotalCyl += nMissedCyl;
        Info<< "Adding " << (nMissedObs + nMissedCyl)
            << " obstacles in groups without instances to default group" << nl;

        #else

        Warning
            << nl << "Found " << (nMissedObs + nMissedCyl)
            << " obstacles in groups without instances" << nl << nl;

        if (pars.debugLevel > 1)
        {
            for (const label groupId : groups.sortedToc())
            {
                const auto& group = groups[groupId];

                if
                (
                    groupId && group.empty()
                 && (group.nObstacle() || group.nCylinder())
                )
                {
                    Info<< "    Group " << groupId << " ("
                        << group.nObstacle() << " obstacles "
                        << group.nCylinder() << " cylinders)"
                        << nl;
                }
            }
        }
        #endif
    }

    return labelPair(nTotalObs, nTotalCyl);
}


Foam::scalar Foam::PDRlegacy::readObstacleFiles
(
    const fileName& obsFileDir,
    const wordList& obsFileNames,
    const Map<obstacleGrouping>& groups,
    const boundBox& meshBb,

    DynamicList<PDRobstacle>& blocks,
    DynamicList<PDRobstacle>& cylinders
)
{
    // Catch programming errors
    if (!groups.found(0))
    {
        FatalErrorInFunction
            << "No default group 0 defined!" << nl
            << exit(FatalError);
    }

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

    scalar shift = pars.obs_expand;

    string buffer;

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

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

        fileName path = (obsFileDir / inputFile);

        IFstream is(path);

        label lineNo = 0;
        while (is.good())
        {
            // Process each line of obstacle files
            ++lineNo;
            is.getLine(buffer);

            const auto firstCh = buffer.find_first_not_of(" \t\n\v\f\r");

            if
            (
                firstCh == std::string::npos
             || buffer[firstCh] == '#'
            )
            {
                // Empty line or comment line
                continue;
            }

            // Quick reject

            int typeId;     // Int (not label) to match sscanf spec
            double x, y, z; // Double (not scalar) to match sscanf spec

            if
            (
                sscanf(buffer.c_str(), "%d %lf %lf %lf", &typeId, &x, &y, &z) < 4
             || typeId == 0
             || typeId == PDRobstacle::MESH_PLANE
            )
            {
                continue;
            }

            int groupId = typeId / 100;
            typeId %= 100;

            if
            (
                typeId == PDRobstacle::OLD_INLET
             || (typeId == PDRobstacle::GRATING && pars.ignoreGratings)
            )
            {
                // Silent - already warned during counting
                continue;
            }

            if (typeId == 0)
            {
                // Group location - not an obstacle
                continue;
            }

            if (!groups.found(groupId))
            {
                // Catch programming errors.
                // - group should be there after the previous read
                Warning
                    << "Encountered undefined group: " << groupId << nl;
                continue;
            }

            #ifdef USE_ZERO_INSTANCE_GROUPS
            const obstacleGrouping& group =
            (
                groups[groups[groupId].size() ? groupId : 0]
            );
            #else
            const obstacleGrouping& group = groups[groupId];
            #endif

            // Add the obstacle to the list with different position
            // offsets according to its group. Non-group obstacles
            // are treated as group 0, which has a single instance
            // with position (0,0,0) and are added only once.

            PDRobstacle scanObs;

            if
            (
                !scanObs.setFromLegacy
                (
                    (groupId * 100) + typeId,
                    buffer,
                    lineNo,
                    inputFile
                )
            )
            {
                continue;
            }

            scanObs.scale(pars.scale);

            // Ignore anything below minimum width
            if (scanObs.tooSmall(pars.min_width))
            {
                continue;
            }


            for (const point& origin : group)
            {
                // A different (very small) shift for each obstacle
                // so that faces cannot be coincident
                shift += floatSMALL;
                const scalar shift2 = shift * 2.0;

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

                        // Offset for the position of the group
                        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 position of the group
                        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:       // Cuboid "Type 1"
                    case PDRobstacle::LOUVRE_BLOWOFF: // Louvred wall or blow-off panel
                    case PDRobstacle::CUBOID:         // Cuboid
                    case PDRobstacle::WALL_BEAM:      // Beam against wall (treated here as normal cuboid)
                    case PDRobstacle::GRATING:        // Grating
                    case PDRobstacle::RECT_PATCH:     // Inlet, outlet or ather b.c. (rectangular)
                    {
                        // Make a copy
                        PDRobstacle obs(scanObs);

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

                        if (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;
                    }
                }
            }
        }

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


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

    return totVolume;
}


Foam::scalar Foam::PDRobstacle::legacyReadFiles
(
    const fileName& obsFileDir,
    const wordList& obsFileNames,
    const boundBox& meshBb,
    DynamicList<PDRobstacle>& blocks,
    DynamicList<PDRobstacle>& cylinders
)
{
    // Still just with legacy reading

    // Count the obstacles and get the group locations
    Map<PDRlegacy::obstacleGrouping> groups;

    const labelPair obsCounts =
        PDRlegacy::readObstacleFiles(obsFileDir, obsFileNames, groups);

    const label nObstacle = obsCounts.first();
    const label nCylinder = obsCounts.second();

    // Info<< "grouping: " << groups << endl;

    if (!nObstacle && !nCylinder)
    {
        FatalErrorInFunction
            << "No obstacles in domain" << nl
            << exit(FatalError);
    }

    blocks.clear();
    blocks.reserve(4 * max(4, nObstacle));

    cylinders.clear();
    cylinders.reserve(4 * max(4, nCylinder));

    return PDRlegacy::readObstacleFiles
    (
        obsFileDir, obsFileNames, groups,
        meshBb,
        blocks,
        cylinders
    );
}


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