/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
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 .
\*---------------------------------------------------------------------------*/
#include "PDRblock.H"
#include "ListOps.H"
#include "emptyPolyPatch.H"
// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
namespace Foam
{
//- Calculate geometric expansion factor from expansion ratio
inline scalar calcGexp(const scalar expRatio, const label nDiv)
{
return nDiv > 1 ? pow(expRatio, 1.0/(nDiv - 1)) : 0.0;
}
}
// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
bool Foam::PDRblock::checkMonotonic
(
const direction cmpt,
const UList& pts
)
{
const label len = pts.size();
if (!len)
{
return false;
}
const scalar& minVal = pts[0];
for (label i=1; i < len; ++i)
{
if (pts[i] <= minVal)
{
FatalErrorInFunction
<< "Points in " << vector::componentNames[cmpt]
<< " direction do not increase monotonically" << nl
<< flatOutput(pts) << nl << nl
<< exit(FatalError);
}
}
return true;
}
const Foam::PDRblock& Foam::PDRblock::null()
{
return NullObjectRef();
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::PDRblock::adjustSizes()
{
// Adjust i-j-k addressing
sizes().x() = grid_.x().nCells();
sizes().y() = grid_.y().nCells();
sizes().z() = grid_.z().nCells();
if (sizes().x() <= 0 || sizes().y() <= 0 || sizes().z() <= 0)
{
// Sanity check. Silently disallow bad sizing
ijkMesh::clear();
grid_.x().clear();
grid_.y().clear();
grid_.z().clear();
bounds_ = boundBox::invertedBox;
minEdgeLen_ = Zero;
return;
}
// Adjust boundBox
bounds_.min().x() = grid_.x().first();
bounds_.min().y() = grid_.y().first();
bounds_.min().z() = grid_.z().first();
bounds_.max().x() = grid_.x().last();
bounds_.max().y() = grid_.y().last();
bounds_.max().z() = grid_.z().last();
// Min edge length
minEdgeLen_ = GREAT;
for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt)
{
const label nEdge = grid_[cmpt].nCells();
for (label edgei=0; edgei < nEdge; ++edgei)
{
const scalar len = grid_[cmpt].width(edgei);
minEdgeLen_ = min(minEdgeLen_, len);
}
}
}
void Foam::PDRblock::readGridControl
(
const direction cmpt,
const dictionary& dict,
const scalar scaleFactor
)
{
//- The begin/end nodes for each segment
scalarList knots;
// The number of division per segment
labelList divisions;
// The expansion ratio per segment
scalarList expansion;
if (verbose_)
{
Info<< "Reading grid control for "
<< vector::componentNames[cmpt] << " direction" << nl;
}
// Points
// ~~~~~~
dict.readEntry("points", knots);
const label nSegments = (knots.size()-1);
if (nSegments < 1)
{
FatalIOErrorInFunction(dict)
<< "Must define at least two control points:"
<< " in " << vector::componentNames[cmpt]
<< " direction" << nl
<< exit(FatalIOError);
}
// Apply scaling
if (scaleFactor > 0)
{
for (scalar& pt : knots)
{
pt *= scaleFactor;
}
}
// Do points increase monotonically?
checkMonotonic(cmpt, knots);
if (verbose_)
{
Info<< " points : " << flatOutput(knots) << nl;
}
// Divisions
// ~~~~~~~~~
dict.readEntry("nCells", divisions);
label nTotalDivs = 0;
for (const label ndiv : divisions)
{
nTotalDivs += ndiv;
if (ndiv < 1)
{
FatalIOErrorInFunction(dict)
<< "Negative or zero divisions:"
<< " in " << vector::componentNames[cmpt]
<< " direction" << nl
<< exit(FatalIOError);
}
}
if (divisions.size() != nSegments)
{
FatalIOErrorInFunction(dict)
<< "Mismatch in number of divisions and number of points:"
<< " in " << vector::componentNames[cmpt]
<< " direction" << nl
<< exit(FatalIOError);
}
else if (!nTotalDivs)
{
FatalIOErrorInFunction(dict)
<< "No grid divisions at all:"
<< " in " << vector::componentNames[cmpt]
<< " direction" << nl
<< exit(FatalIOError);
}
if (verbose_)
{
Info<< " nCells : " << flatOutput(divisions) << nl;
}
// Expansion ratios
// ~~~~~~~~~~~~~~~~
dict.readIfPresent("ratios", expansion);
// Correct expansion ratios - negative is the same as inverse.
for (scalar& expRatio : expansion)
{
if (expRatio < 0)
{
expRatio = 1.0/(-expRatio);
}
}
if (expansion.size() && expansion.size() != nSegments)
{
FatalIOErrorInFunction(dict)
<< "Mismatch in number of expansion ratios and divisions:"
<< " in " << vector::componentNames[cmpt]
<< " direction" << nl
<< exit(FatalIOError);
}
if (expansion.empty())
{
expansion.resize(nSegments, scalar(1));
if (verbose_)
{
Info<< "Warning: 'ratios' not specified, using constant spacing"
<< nl;
}
}
else
{
if (verbose_)
{
Info<< " ratios : " << flatOutput(expansion) << nl;
}
}
// Define the grid points
scalarList& gridPoint = grid_[cmpt];
gridPoint.resize(nTotalDivs+1);
label start = 0;
for (label segmenti=0; segmenti < nSegments; ++segmenti)
{
const label nDiv = divisions[segmenti];
const scalar expRatio = expansion[segmenti];
SubList 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(expRatio, 1.0))
{
for (label i=1; i < nDiv; ++i)
{
subPoint[i] = (subPoint[0] + (dist * i)/nDiv);
}
}
else
{
// Geometric expansion factor from the expansion ratio
const scalar expFact = calcGexp(expRatio, nDiv);
for (label i=1; i < nDiv; ++i)
{
subPoint[i] =
(
subPoint[0]
+ dist * (1.0 - pow(expFact, i))/(1.0 - pow(expFact, nDiv))
);
}
}
}
}
void Foam::PDRblock::readBoundary(const dictionary& dict)
{
if (verbose_)
{
Info<< "Reading boundary entries" << nl;
}
PtrList patchEntries;
if (dict.found("boundary"))
{
PtrList newEntries(dict.lookup("boundary"));
patchEntries.transfer(newEntries);
}
patches_.clear();
patches_.resize(patchEntries.size());
// Hex cell has 6 sides
const labelRange validFaceId(0, 6);
// Track which sides have been handled
labelHashSet handled;
forAll(patchEntries, patchi)
{
if (!patchEntries.set(patchi))
{
continue;
}
const entry& e = patchEntries[patchi];
if (!e.isDict())
{
FatalIOErrorInFunction(dict)
<< "Entry " << e
<< " in boundary section is not a valid dictionary."
<< exit(FatalIOError);
}
const dictionary& patchDict = e.dict();
const word& patchName = e.keyword();
const word patchType(patchDict.get("type"));
labelList faceLabels(patchDict.get("faces"));
labelHashSet patchFaces(faceLabels);
if (faceLabels.size() != patchFaces.size())
{
WarningInFunction
<< "Patch: " << patchName
<< ": Ignoring duplicate face ids" << nl;
}
label nErased = patchFaces.filterKeys(validFaceId);
if (nErased)
{
WarningInFunction
<< "Patch: " << patchName << ": Ignoring " << nErased
<< " faces with invalid identifiers" << nl;
}
nErased = patchFaces.erase(handled);
if (nErased)
{
WarningInFunction
<< "Patch: " << patchName << ": Ignoring " << nErased
<< " faces, which were already handled" << nl;
}
// Mark these faces as having been handled
handled += patchFaces;
// Save information for later access during mesh creation.
patches_.set(patchi, new boundaryEntry());
boundaryEntry& bentry = patches_[patchi];
bentry.name_ = patchName;
bentry.type_ = patchType;
bentry.size_ = 0;
bentry.faces_ = patchFaces.sortedToc();
// Count patch faces
for (const label shapeFacei : bentry.faces_)
{
bentry.size_ += nBoundaryFaces(shapeFacei);
}
}
// Deal with unhandled faces
labelHashSet missed(identity(6));
missed.erase(handled);
if (missed.size())
{
patches_.append(new boundaryEntry());
boundaryEntry& bentry = patches_.last();
bentry.name_ = "defaultFaces";
bentry.type_ = emptyPolyPatch::typeName;
bentry.size_ = 0;
bentry.faces_ = missed.sortedToc();
// Count patch faces
for (const label shapeFacei : bentry.faces_)
{
bentry.size_ += nBoundaryFaces(shapeFacei);
}
// Use name/type from "defaultPatch" entry if available
// - be generous with error handling
const dictionary* defaultEntry = dict.findDict("defaultPatch");
if (defaultEntry)
{
defaultEntry->readIfPresent("name", bentry.name_);
defaultEntry->readIfPresent("type", bentry.type_);
}
}
if (verbose_)
{
label patchi = 0;
for (const boundaryEntry& bentry : patches_)
{
Info<< " patch: " << patchi
<< " (size: " << bentry.size_
<< " type: " << bentry.type_
<< ") name: " << bentry.name_
<< " faces: " << flatOutput(bentry.faces_) << nl;
++patchi;
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::PDRblock::PDRblock
(
const UList& xgrid,
const UList& ygrid,
const UList& zgrid
)
:
PDRblock()
{
// Default boundaries with patchi == shapeFacei
patches_.resize(6);
for (label patchi=0; patchi < 6; ++patchi)
{
patches_.set(patchi, new boundaryEntry());
boundaryEntry& bentry = patches_[patchi];
bentry.name_ = "patch" + Foam::name(patchi);
bentry.type_ = "patch";
bentry.size_ = 0;
bentry.faces_ = labelList(one(), patchi);
}
reset(xgrid, ygrid, zgrid);
}
Foam::PDRblock::PDRblock(const dictionary& dict, bool verbose)
:
PDRblock()
{
verbose_ = verbose;
read(dict);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::PDRblock::read(const dictionary& dict)
{
// Mark no scaling with invalid value
const scalar scaleFactor = dict.lookupOrDefault("scale", -1);
readGridControl(0, dict.subDict("x"), scaleFactor);
readGridControl(1, dict.subDict("y"), scaleFactor);
readGridControl(2, dict.subDict("z"), scaleFactor);
adjustSizes();
readBoundary(dict);
return true;
}
void Foam::PDRblock::reset
(
const UList& xgrid,
const UList& ygrid,
const UList& zgrid
)
{
static_cast(grid_.x()) = xgrid;
static_cast(grid_.y()) = ygrid;
static_cast(grid_.z()) = zgrid;
#ifdef FULLDEBUG
for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt)
{
checkMonotonic(cmpt, grid_[cmpt]);
}
#endif
adjustSizes();
// Adjust boundaries
for (boundaryEntry& bentry : patches_)
{
bentry.size_ = 0;
// Count patch faces
for (const label shapeFacei : bentry.faces_)
{
bentry.size_ += nBoundaryFaces(shapeFacei);
}
}
}
bool Foam::PDRblock::findCell(const point& pt, labelVector& pos) const
{
// Out-of-bounds is handled explicitly, for efficiency and consistency,
// but principally to ensure that findLower() returns a valid
// result when the point is to the right of the bounds.
// Since findLower returns the lower index, it corresponds to the
// cell in which the point is found
if (!bounds_.contains(pt))
{
return false;
}
for (direction cmpt=0; cmpt < labelVector::nComponents; ++cmpt)
{
// Binary search
pos[cmpt] = findLower(grid_[cmpt], pt[cmpt]);
}
return true;
}
bool Foam::PDRblock::gridIndex
(
const point& pt,
labelVector& pos,
const scalar relTol
) const
{
const scalar tol = relTol * minEdgeLen_;
for (direction cmpt=0; cmpt < labelVector::nComponents; ++cmpt)
{
// Linear search
pos[cmpt] = grid_[cmpt].findIndex(pt[cmpt], tol);
if (pos[cmpt] < 0) return false;
}
return true;
}
Foam::labelVector Foam::PDRblock::findCell(const point& pt) const
{
labelVector pos;
if (findCell(pt, pos))
{
return pos;
}
return labelVector(-1,-1,-1);
}
Foam::labelVector Foam::PDRblock::gridIndex
(
const point& pt,
const scalar relTol
) const
{
labelVector pos;
if (gridIndex(pt, pos, relTol))
{
return pos;
}
return labelVector(-1,-1,-1);
}
// ************************************************************************* //