779c3fe11e
- avoids naming ambiguity between querying sorted state vs returning a sorted list (for example) ENH: add 'good()' method to a few more classes
287 lines
8.6 KiB
C
287 lines
8.6 KiB
C
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / 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 "PDRmeshArrays.H"
|
|
#include "PDRblock.H"
|
|
#include "polyMesh.H"
|
|
#include "Time.H"
|
|
#include "IjkField.H"
|
|
|
|
// Notes
|
|
//
|
|
// Determines the face and cell numbers of all faces and cells in the
|
|
// central rectangular region where CAD_PDR operates. First,
|
|
// "points" is read and the coordinates (by which I mean here the
|
|
// indices in the x, y and z coordinate arrays) are determined. Then
|
|
// "faces" is read and for each the coordinates of the lower- x,y,z
|
|
// corner are determioned, also the orientation (X, Y or Z).
|
|
// (Orientation in the sense of e.g. + or -x is not noted.) The files
|
|
// "owner" and "neighbour" specify the six faces around each cell, so
|
|
// from these the coordinates of the cells are determined.
|
|
//
|
|
// Full checks are made that the mesh in the central region is consistent
|
|
// with CAD_PDR's mesh specified by the PDRmeshSpec file.
|
|
//
|
|
// Eventually, when writing out results, we shall work through the
|
|
// full list of cells, writing default values for any cells that are
|
|
// not in the central regtion.
|
|
|
|
|
|
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
|
|
|
Foam::scalar Foam::PDRmeshArrays::gridPointRelTol = 0.02;
|
|
|
|
|
|
// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
|
|
|
|
namespace
|
|
{
|
|
|
|
// A good ijk index has all components >= 0
|
|
static inline bool isGoodIndex(const Foam::labelVector& idx)
|
|
{
|
|
return (idx.x() >= 0 && idx.y() >= 0 && idx.z() >= 0);
|
|
}
|
|
|
|
} // End anonymous namespace
|
|
|
|
|
|
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
|
|
|
void Foam::PDRmeshArrays::classify
|
|
(
|
|
const polyMesh& mesh,
|
|
const PDRblock& pdrBlock
|
|
)
|
|
{
|
|
// Additional copy of i-j-k addressing
|
|
cellDims = pdrBlock.sizes();
|
|
faceDims = (cellDims + labelVector::one);
|
|
|
|
const label maxPointId = cmptMax(pdrBlock.sizes())+1;
|
|
|
|
Info<< "Mesh" << nl
|
|
<< " nPoints:" << mesh.nPoints()
|
|
<< " nCells:" << mesh.nCells()
|
|
<< " nFaces:" << mesh.nFaces() << nl;
|
|
|
|
Info<< "PDRblock" << nl
|
|
<< " nPoints:" << pdrBlock.nPoints()
|
|
<< " nCells:" << pdrBlock.nCells()
|
|
<< " nFaces:" << pdrBlock.nFaces() << nl
|
|
<< " min-edge:" << pdrBlock.edgeLimits().min() << nl;
|
|
|
|
Info<< "Classifying ijk indexing... " << nl;
|
|
|
|
|
|
// Bin cells into i-j-k locations with the PDRblock::findCell()
|
|
// method, which combines a bounding box rejection and binary
|
|
// search in the three directions.
|
|
|
|
cellIndex.resize(mesh.nCells());
|
|
{
|
|
const pointField& cc = mesh.cellCentres();
|
|
|
|
for (label celli = 0; celli < mesh.nCells(); ++celli)
|
|
{
|
|
cellIndex[celli] = pdrBlock.findCell(cc[celli]);
|
|
}
|
|
}
|
|
|
|
// Verify that all i-j-k cells were indeed found
|
|
{
|
|
// This could be more efficient - but we want to be picky
|
|
IjkField<bool> cellFound(pdrBlock.sizes(), false);
|
|
|
|
for (label celli=0; celli < cellIndex.size(); ++celli)
|
|
{
|
|
const labelVector& cellIdx = cellIndex[celli];
|
|
|
|
if (isGoodIndex(cellIdx))
|
|
{
|
|
cellFound(cellIdx) = true;
|
|
}
|
|
}
|
|
|
|
const label firstMiss = cellFound.find(false);
|
|
|
|
if (firstMiss >= 0)
|
|
{
|
|
label nMissing = 0;
|
|
for (label celli = firstMiss; celli < cellFound.size(); ++celli)
|
|
{
|
|
if (!cellFound[celli])
|
|
{
|
|
++nMissing;
|
|
}
|
|
}
|
|
|
|
FatalErrorInFunction
|
|
<< "No ijk location found for "
|
|
<< nMissing << " cells.\nFirst miss at: "
|
|
<< pdrBlock.index(firstMiss)
|
|
<< " indexing" << nl
|
|
<< exit(FatalError);
|
|
}
|
|
}
|
|
|
|
|
|
// Bin all mesh points into i-j-k locations
|
|
List<labelVector> pointIndex(mesh.nPoints());
|
|
|
|
for (label pointi = 0; pointi < mesh.nPoints(); ++pointi)
|
|
{
|
|
const point& p = mesh.points()[pointi];
|
|
pointIndex[pointi] = pdrBlock.gridIndex(p, gridPointRelTol);
|
|
}
|
|
|
|
// Min x,y,z index
|
|
const labelMinMax invertedLimits(maxPointId, -maxPointId);
|
|
Vector<labelMinMax> faceLimits;
|
|
|
|
const Vector<direction> faceBits
|
|
(
|
|
boundBox::XDIR,
|
|
boundBox::YDIR,
|
|
boundBox::ZDIR
|
|
);
|
|
|
|
faceIndex.resize(mesh.nFaces());
|
|
faceOrient.resize(mesh.nFaces());
|
|
|
|
for (label facei=0; facei < mesh.nFaces(); ++facei)
|
|
{
|
|
labelVector& faceIdx = faceIndex[facei];
|
|
|
|
// Check faces that are associated with i-j-k cells
|
|
const label own = mesh.faceOwner()[facei];
|
|
const label nei =
|
|
(
|
|
facei < mesh.nInternalFaces()
|
|
? mesh.faceNeighbour()[facei]
|
|
: own
|
|
);
|
|
|
|
if (!isGoodIndex(cellIndex[own]) && !isGoodIndex(cellIndex[nei]))
|
|
{
|
|
// Invalid
|
|
faceIdx.x() = faceIdx.y() = faceIdx.z() = -1;
|
|
faceOrient[facei] = vector::X;
|
|
continue;
|
|
}
|
|
|
|
|
|
faceLimits.x() = faceLimits.y() = faceLimits.z() = invertedLimits;
|
|
|
|
for (const label pointi : mesh.faces()[facei])
|
|
{
|
|
for (direction cmpt=0; cmpt < labelVector::nComponents; ++cmpt)
|
|
{
|
|
faceLimits[cmpt].add(pointIndex[pointi][cmpt]);
|
|
}
|
|
}
|
|
|
|
direction inPlane(0u);
|
|
|
|
for (direction cmpt=0; cmpt < labelVector::nComponents; ++cmpt)
|
|
{
|
|
const auto& limits = faceLimits[cmpt];
|
|
|
|
if (!limits.good())
|
|
{
|
|
// This should be impossible
|
|
FatalErrorInFunction
|
|
<< "Unexpected search failure for " << facei << " in "
|
|
<< vector::componentNames[cmpt] << "-direction" << nl
|
|
<< exit(FatalError);
|
|
}
|
|
|
|
if (limits.min() < 0)
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Face " << facei << " contains non-grid point in "
|
|
<< vector::componentNames[cmpt] << "-direction" << nl
|
|
<< mesh.faces()[facei] << ' '
|
|
<< mesh.faces()[facei].points(mesh.points())
|
|
<< exit(FatalError);
|
|
}
|
|
else if (limits.min() == limits.max())
|
|
{
|
|
// In plane
|
|
inPlane |= faceBits[cmpt];
|
|
}
|
|
else if (limits.min() + 1 != limits.max())
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Face " << facei << " not in "
|
|
<< vector::componentNames[cmpt] << "-plane" << nl
|
|
<< exit(FatalError);
|
|
}
|
|
}
|
|
|
|
switch (inPlane)
|
|
{
|
|
case boundBox::XDIR:
|
|
faceOrient[facei] = vector::X;
|
|
break;
|
|
|
|
case boundBox::YDIR:
|
|
faceOrient[facei] = vector::Y;
|
|
break;
|
|
|
|
case boundBox::ZDIR:
|
|
faceOrient[facei] = vector::Z;
|
|
break;
|
|
|
|
default:
|
|
FatalErrorInFunction
|
|
<< "Face " << facei << " not in an x/y/z plane?" << nl
|
|
<< exit(FatalError);
|
|
break;
|
|
}
|
|
|
|
faceIdx.x() = faceLimits.x().min();
|
|
faceIdx.y() = faceLimits.y().min();
|
|
faceIdx.z() = faceLimits.z().min();
|
|
}
|
|
}
|
|
|
|
|
|
void Foam::PDRmeshArrays::read
|
|
(
|
|
const Time& runTime,
|
|
const PDRblock& pdrBlock
|
|
)
|
|
{
|
|
#include "createPolyMesh.H"
|
|
classify(mesh, pdrBlock);
|
|
}
|
|
|
|
|
|
// ************************************************************************* //
|