BUG: PDRsetFields fails with non-orthogonal outer region (fixes #1907)

- first sort the cells into their ijk bins, and restrict testing for
  face orientation to those faces with an owner or neighbour that has
  an ijk bin.

ENH: ensure polyMesh from PDRblockMesh is marked as AUTO_WRITE

- the particular polyMesh constructor inherits the writeOpt, which
  makes is dependent on the caller and thus somewhat fragile for the
  top level caller.

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

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 points into i-j-k locations
+    // 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& pt = mesh.points()[pointi];
+        const point& p = mesh.points()[pointi];
-        pointIndex[pointi] = pdrBlock.gridIndex(pt, gridPointRelTol);
+        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
+                    << "Face " << facei << " not in "
-                    << " not in " << vector::componentNames[cmpt] << "-plane" << nl
+                    << vector::componentNames[cmpt] << "-plane" << nl
                     << exit(FatalError);
             }
         }
@@ -172,78 +265,9 @@ void Foam::PDRmeshArrays::classify
                 break;
         }
 
-        faceIndex[facei] =
+        faceIdx.x() = faceLimits.x().min();
-            labelVector
+        faceIdx.y() = faceLimits.y().min();
-            (
+        faceIdx.z() = faceLimits.z().min();
-                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);
-        }
     }
 }
 

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);
-
 };
 
 

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