diff --git a/applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRarraysCalc.C b/applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRarraysCalc.C
index 3b4c27ccd4b1c022b7f1553f803c6cbbf163cd81..ada22eb6a0530025e8c6f3d75e1194b8370f414e 100644
--- a/applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRarraysCalc.C
+++ b/applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRarraysCalc.C
@@ -52,9 +52,25 @@ License
#endif
#include <cassert>
+// * * * * * * * * * * * * * * * 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
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
using namespace Foam;
-HashTable<string> fieldNotes
+static Foam::HashTable<Foam::string> fieldNotes
({
{ "Lobs", "" },
{ "Aw", "surface area per unit volume" },
@@ -1255,7 +1271,7 @@ void write_scalarField
{
const labelVector& cellIdx = meshIndexing.cellIndex[celli];
- if (cmptMin(cellIdx) < 0)
+ if (!isGoodIndex(cellIdx))
{
os << deflt << nl;
continue;
@@ -1529,7 +1545,7 @@ void write_symmTensorField
{
const labelVector& cellIdx = meshIndexing.cellIndex[celli];
- if (cmptMin(cellIdx) < 0)
+ if (!isGoodIndex(cellIdx))
{
os << deflt << nl;
continue;
@@ -1591,7 +1607,7 @@ void write_symmTensorFieldV
{
const labelVector& cellIdx = meshIndexing.cellIndex[celli];
- if (cmptMin(cellIdx) < 0)
+ if (!isGoodIndex(cellIdx))
{
os << deflt << nl;
continue;
@@ -1657,7 +1673,7 @@ void write_blocked_face_list
// The related i-j-k face index for the mesh face
const labelVector& faceIdx = meshIndexing.faceIndex[facei];
- if (cmptMin(faceIdx) < 0)
+ if (!isGoodIndex(faceIdx))
{
continue;
}
@@ -1923,7 +1939,7 @@ void write_blockedCellsSet
{
const labelVector& cellIdx = meshIndexing.cellIndex[celli];
- if (cmptMin(cellIdx) < 0)
+ if (!isGoodIndex(cellIdx))
{
continue;
}
diff --git a/applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRmeshArrays.C b/applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRmeshArrays.C
index 8e6a6c1fd95496c471d8ec8c4f2087f443e1e9b1..230638132291cbafc15377df04ccaf8f880b5f78 100644
--- a/applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRmeshArrays.C
+++ b/applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRmeshArrays.C
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2016 Shell Research Ltd.
- Copyright (C) 2019 OpenCFD Ltd.
+ Copyright (C) 2019-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -57,6 +57,20 @@ License
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
@@ -77,16 +91,73 @@ void Foam::PDRmeshArrays::classify
<< " nFaces:" << mesh.nFaces() << nl;
Info<< "PDRblock" << nl
- << " minEdgeLen:" << pdrBlock.minEdgeLen() << nl;
+ << " nPoints:" << pdrBlock.nPoints()
+ << " nCells:" << pdrBlock.nCells()
+ << " nFaces:" << pdrBlock.nFaces() << nl
+ << " min-edge:" << pdrBlock.minEdgeLen() << 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;
+ }
+ }
- // Bin points into i-j-k locations
+ 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)
+ for (label pointi = 0; pointi < mesh.nPoints(); ++pointi)
{
- const point& pt = mesh.points()[pointi];
- pointIndex[pointi] = pdrBlock.gridIndex(pt, gridPointRelTol);
+ const point& p = mesh.points()[pointi];
+ pointIndex[pointi] = pdrBlock.gridIndex(p, gridPointRelTol);
}
// Min x,y,z index
@@ -105,6 +176,26 @@ void Foam::PDRmeshArrays::classify
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])
@@ -135,6 +226,8 @@ void Foam::PDRmeshArrays::classify
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())
@@ -145,8 +238,8 @@ void Foam::PDRmeshArrays::classify
else if (limits.min() + 1 != limits.max())
{
FatalErrorInFunction
- << "Face " << facei
- << " not in " << vector::componentNames[cmpt] << "-plane" << nl
+ << "Face " << facei << " not in "
+ << vector::componentNames[cmpt] << "-plane" << nl
<< exit(FatalError);
}
}
@@ -172,78 +265,9 @@ void Foam::PDRmeshArrays::classify
break;
}
- faceIndex[facei] =
- labelVector
- (
- faceLimits.x().min(),
- faceLimits.y().min(),
- faceLimits.z().min()
- );
- }
-
-
- // Bin cells into i-j-k locations
- cellIndex = std::move(pointIndex);
- cellIndex = labelVector::uniform(maxPointId);
- cellIndex.resize(mesh.nCells(), labelVector::uniform(maxPointId));
-
- // Option 1: use PDRblock.findCell() method
- if (true)
- {
- const pointField& cc = mesh.cellCentres();
-
- for (label celli=0; celli < mesh.nCells(); ++celli)
- {
- cellIndex[celli] = pdrBlock.findCell(cc[celli]);
- }
- }
-
- // Option 2: walk cell faces and use faceIndex information
- if (false)
- {
- for (label celli=0; celli < mesh.nCells(); ++celli)
- {
- labelVector& cellIdx = cellIndex[celli];
-
- for (const label facei : mesh.cells()[celli])
- {
- cellIdx.x() = min(cellIdx.x(), faceIndex[facei].x());
- cellIdx.y() = min(cellIdx.y(), faceIndex[facei].y());
- cellIdx.z() = min(cellIdx.z(), faceIndex[facei].z());
- }
-
- if (cmptMin(cellIdx) < 0)
- {
- cellIdx = labelVector(-1,-1,-1);
- }
- }
- }
-
-
- // Verify that all i-j-k cells were 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 (cmptMin(cellIdx) >= 0)
- {
- cellFound(cellIdx) = true;
- }
- }
-
- label firstMissing = cellFound.find(false);
-
- if (firstMissing >= 0)
- {
- FatalErrorInFunction
- << "No cell found for " << pdrBlock.index(firstMissing)
- << " indexing"
- << exit(FatalError);
- }
+ faceIdx.x() = faceLimits.x().min();
+ faceIdx.y() = faceLimits.y().min();
+ faceIdx.z() = faceLimits.z().min();
}
}
diff --git a/applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRmeshArrays.H b/applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRmeshArrays.H
index 26142e54fb85b7784363d469bd0049f2937b6d64..6c2e329721478806319369f863df895604f8dc07 100644
--- a/applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRmeshArrays.H
+++ b/applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRmeshArrays.H
@@ -88,7 +88,7 @@ public:
// Constructors
- //- Construct null
+ //- Default construct
PDRmeshArrays() = default;
@@ -116,7 +116,6 @@ public:
//- Read OpenFOAM mesh and determine i-j-k indices for faces/cells
void read(const Time& runTime, const PDRblock& pdrBlock);
-
};
diff --git a/src/mesh/blockMesh/PDRblockMesh/PDRblockCreate.C b/src/mesh/blockMesh/PDRblockMesh/PDRblockCreate.C
index e12cff26c386f2645ff66381226e722c6a3ef791..61f07e7582d8ee56d3a8455e65090b66efd6def7 100644
--- a/src/mesh/blockMesh/PDRblockMesh/PDRblockCreate.C
+++ b/src/mesh/blockMesh/PDRblockMesh/PDRblockCreate.C
@@ -328,14 +328,19 @@ Foam::autoPtr<Foam::polyMesh> Foam::PDRblock::mesh(const IOobject& io) const
}
}
+
+ IOobject iomesh(io);
+ iomesh.writeOpt() = IOobject::AUTO_WRITE;
+
auto meshPtr = autoPtr<polyMesh>::New
(
- io,
+ iomesh,
std::move(pts),
std::move(faces),
std::move(own),
std::move(nei)
);
+ polyMesh& pmesh = *meshPtr;
PtrList<polyPatch> patches(patches_.size());
@@ -355,7 +360,7 @@ Foam::autoPtr<Foam::polyMesh> Foam::PDRblock::mesh(const IOobject& io) const
bentry.size_,
startFace,
patchi, // index
- meshPtr->boundaryMesh()
+ pmesh.boundaryMesh()
)
);
@@ -365,7 +370,7 @@ Foam::autoPtr<Foam::polyMesh> Foam::PDRblock::mesh(const IOobject& io) const
++patchi;
}
- meshPtr->addPatches(patches);
+ pmesh.addPatches(patches);
return meshPtr;
}