PDRlegacyMeshSpec.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 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/>.
 
\*---------------------------------------------------------------------------*/
 
// Read spec that resemble this:
//
//   xmesh
//   (
//       ( -8.18  14  0.83 )
//       ( -0.69  11  1.00 )
//       (  0.69  14  1.20 )
//       (  8.18  0 )
//   )
//
//   ymesh
//   (
//       ...
//   )
 
#include "PDRlegacy.H"
 
// OpenFOAM includes
#include "error.H"
#include "IFstream.H"
#include "stringOps.H"
#include "OSspecific.H"
 
#define XMESH_TAG "xmesh"
#define YMESH_TAG "ymesh"
#define ZMESH_TAG "zmesh"
 
 
// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
 
namespace Foam
{
namespace PDRlegacy
{
namespace Detail
{
 
//- Simple structure for processing a mesh spec entry
//  Handles an entry with (float), (float int), or (float int float)
//
//  This is for handling a sub-spec that resembles this:
//
//  (
//      ( -8.18  14  0.83 )
//      ( -0.69  11  1.00 )
//      (  0.69  14  1.20 )
//      (  8.18  0 )
//  )
struct pdrMeshSpecLine
{
    scalar knot;
    label  ndiv;
    scalar factor;
 
    pdrMeshSpecLine() : knot(0), ndiv(0), factor(0) {}
 
    //- Cheap means to avoid uniform list output
    bool operator!=(const pdrMeshSpecLine&) const
    {
        return false;
    }
};
 
 
// Read mesh-spec entry
Istream& operator>>(Istream& is, pdrMeshSpecLine& spec)
{
    spec.knot = 0;
    spec.ndiv = 0;
    spec.factor = 0;
 
    is.readBegin("pdrMeshSpecLine");
 
    // Must have a point
    is >> spec.knot;
 
    token tok(is);
    if (tok.isLabel())
    {
        spec.ndiv = tok.labelToken();
 
        if (spec.ndiv)
        {
            is >> spec.factor;
        }
    }
    else
    {
        is.putBack(tok);
    }
 
    is.readEnd("pdrMeshSpecLine");
 
    is.check(FUNCTION_NAME);
    return is;
}
 
#ifdef FULLDEBUG
// Write mesh-spec entry
Ostream& operator<<(Ostream& os, const pdrMeshSpecLine& spec)
{
    os << token::BEGIN_LIST << spec.knot;
 
    if (spec.ndiv)
    {
        os  << token::SPACE << spec.ndiv
            << token::SPACE << spec.factor;
    }
 
    os << token::END_LIST;
 
    return os;
}
#endif
 
 
void read_spec(ISstream& is, const direction cmpt, List<scalar>& gridPoint)
{
    if (!gridPoint.empty())
    {
        FatalErrorInFunction
            << "Duplicate specification of "
            << vector::componentNames[cmpt]
            << " grid"
            << exit(FatalError);
    }
 
    List<pdrMeshSpecLine> specs(is);
 
    if (specs.size() < 2)
    {
        FatalErrorInFunction
            << "Grid specification for " << vector::componentNames[cmpt]
            << " is too small. Need at least two points!" << nl
            << exit(FatalError);
    }
 
    specs.last().ndiv = 0; // safety
 
 
    DynamicList<scalar> knots;
    DynamicList<label>  divisions;
    DynamicList<scalar> factors;
 
    for (const auto& spec : specs)
    {
        knots.append(spec.knot);
 
        if (spec.ndiv < 1)
        {
            break;
        }
        divisions.append(spec.ndiv);
        factors.append(spec.factor);
    }
 
    label nPoints = 1;
    for (const label nDiv : divisions)
    {
        nPoints += nDiv;
    }
 
    if (nPoints < 2)
    {
        FatalErrorInFunction
            << "No cells defined for direction "
            << vector::componentNames[cmpt] << nl
            << exit(FatalError);
    }
 
 
    // Define the grid points
    gridPoint.resize(nPoints);
 
    const label nSegments = divisions.size();
 
    label start = 0;
 
    for (label segmenti=0; segmenti < nSegments; ++segmenti)
    {
        const label nDiv = divisions[segmenti];
        const scalar factor = factors[segmenti];
 
        SubList<scalar> subPoint(gridPoint, nDiv+1, start);
        start += nDiv;
 
        subPoint[0] = knots[segmenti];
        subPoint[nDiv] = knots[segmenti+1];
 
        const scalar dist = (subPoint.last() - subPoint.first());
 
        if (equal(factor, scalar(1)))
        {
            for (label i=1; i < nDiv; ++i)
            {
                subPoint[i] = (subPoint[0] + (dist * i)/nDiv);
            }
        }
        else
        {
            scalar delta = dist * (1.0 - factor) / (1.0 - ::pow(factor, nDiv));
 
            scalar xyz = subPoint[0];
 
            for (label i=0; i < nDiv; ++i)
            {
                subPoint[i] = xyz;
                xyz += delta;
                delta *= factor;
            }
        }
    }
}
 
 
} // End namespace Detail
} // End namespace PDRlegacy
} // End namespace Foam
 
 
// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 
void Foam::PDRlegacy::print_info(const PDRblock& block)
{
    Info<< "PDRblock" << nl
        << "    nCells: " << block.sizes() << nl
        << "    Box: " << block.bounds() << nl
        << "x " << flatOutput(block.grid().x()) << nl
        << "y " << flatOutput(block.grid().y()) << nl
        << "z " << flatOutput(block.grid().z()) << nl
        << endl;
}
 
 
void Foam::PDRlegacy::read_mesh_spec(const fileName& casepath, PDRblock& block)
{
    Info<< "Reading pdrMeshSpec (legacy format)" << nl;
 
    bool processed = false;
 
    for (const fileName dirName : { "system", "constant/polyMesh" })
    {
        fileName path
        (
            casepath / dirName / "pdrMeshSpec"
        );
 
        if (Foam::isFile(path))
        {
            IFstream is(path);
 
            read_mesh_spec(is, block);
            processed = true;
            break;
        }
    }
 
    if (!processed)
    {
        FatalErrorInFunction
            << "Did not process pdrMeshSpec" << nl
            << exit(FatalError);
    }
}
 
 
void Foam::PDRlegacy::read_mesh_spec(ISstream& is, PDRblock& block)
{
    Vector<scalarList> grid;
 
    string line;
 
    while (is.good())
    {
        is.getLine(line);
        stringOps::inplaceTrim(line);
 
        if (line == XMESH_TAG)
        {
            Detail::read_spec(is, vector::X, grid.x());
        }
        else if (line == YMESH_TAG)
        {
            Detail::read_spec(is, vector::Y, grid.y());
        }
        else if (line == ZMESH_TAG)
        {
            Detail::read_spec(is, vector::Z, grid.z());
        }
    }
 
    for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt)
    {
        if (grid[cmpt].empty())
        {
            FatalErrorInFunction
                << "No specification for "
                << vector::componentNames[cmpt]
                << " grid" << nl
                << exit(FatalError);
        }
    }
 
    block.reset(grid.x(), grid.y(), grid.z());
 
    #ifdef FULLDEBUG
    print_info(block);
    #endif
}
 
 
// ************************************************************************* //