ENH: add findCell for PDRblockMesh

applications/test/IjkField/Test-IjkField.C

@@ -49,6 +49,36 @@ Ostream& print(const IjkField<T>& fld)
 
 int main(int argc, char *argv[])
 {
+    // Basic addressing checks
+    #if 0
+    {
+        ijkAddressing addr1(3, 4, 5);
+
+        Info<< "addressing: " << addr1.sizes() << nl;
+        Info<< "index of (2,2,2) " << addr1.index(2,2,2) << nl;
+
+        for (const label idx : labelRange(0, addr1.size()))
+        {
+            Info<< "index of " << idx << " => " << addr1.index(idx) << nl;
+        }
+
+        for (label k=0; k < addr1.sizes().z(); ++k)
+        {
+            for (label j=0; j < addr1.sizes().y(); ++j)
+            {
+                for (label i=0; i < addr1.sizes().x(); ++i)
+                {
+                    labelVector ijk(i,j,k);
+
+                    Info<< "index of " << addr1.index(ijk)
+                        << " <= " << ijk << nl;
+                }
+            }
+        }
+    }
+    #endif
+
+
     // Create with inconsistent sizes
     IjkField<label> field0({1, 2, 3}, identity(10));
 

applications/test/PDRblockMesh/Test-PDRblockMesh.C

@@ -157,6 +157,28 @@ int main(int argc, char *argv[])
             << "  length = " << mesh.width(mid) << nl;
     }
 
+    // Test findCell
+    {
+        Info<< nl << "findCell:" << nl;
+
+        for
+        (
+            const point& pt
+          : {
+                mesh.bounds().centre(),
+                mesh.bounds().min() - 0.1 * mesh.bounds().span(),
+                mesh.bounds().max() + 0.1 * mesh.bounds().span()
+            }
+        )
+        {
+            labelVector ijk = mesh.findCell(pt);
+
+            Info<< "    " << pt << " = " << ijk;
+            if (cmptMin(ijk) < 0) Info<< " [not found]";
+            Info<< nl;
+        }
+    }
+
     Info<< nl;
 
     // Fatal with FULLDEBUG

src/OpenFOAM/meshes/ijkMesh/ijkAddressing.H

@@ -50,7 +50,7 @@ class ijkAddressing
 {
     // Private Data
 
-        //- The number of cells in the i,j,k directions.
+        //- The number of items in the i,j,k directions.
         labelVector sizes_;
 
 
@@ -99,6 +99,9 @@ public:
         //- Linear addressing index (offset) for an i,j,k position.
         inline label index(const labelVector& ijk) const;
 
+        //- The i,j,k indexing from linear addressing.
+        inline labelVector index(const label idx) const;
+
 
     // Checks
 

src/OpenFOAM/meshes/ijkMesh/ijkAddressingI.H

@@ -141,6 +141,19 @@ inline Foam::label Foam::ijkAddressing::index(const labelVector& ijk) const
 }
 
 
+inline Foam::labelVector Foam::ijkAddressing::index(const label idx) const
+{
+    // No checkIndex
+
+    return labelVector
+    (
+        idx % (sizes_.x()),
+        idx / (sizes_.x()),
+        idx / (sizes_.x() * sizes_.y())
+    );
+}
+
+
 inline void Foam::ijkAddressing::checkIndex
 (
     const label i,

src/mesh/blockMesh/PDRblockMesh/PDRblock.C

@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "PDRblock.H"
+#include "ListOps.H"
 #include "emptyPolyPatch.H"
 
 // * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
@@ -437,4 +438,24 @@ bool Foam::PDRblock::read(const dictionary& dict)
 }
 
 
+Foam::labelVector Foam::PDRblock::findCell(const point& pt) 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.
+
+    labelVector pos(-1,-1,-1);
+    if (bounds_.contains(pt))
+    {
+        for (direction cmpt=0; cmpt < labelVector::nComponents; ++cmpt)
+        {
+            // Binary search
+            pos[cmpt] = findLower(grid_[cmpt], pt[cmpt]);
+        }
+    }
+
+    return pos;
+}
+
+
 // ************************************************************************* //

src/mesh/blockMesh/PDRblockMesh/PDRblock.H

@@ -266,6 +266,13 @@ public:
         inline scalar width(const labelVector& ijk) const;
 
 
+    // Searching
+
+        //- Return i,j,k index for cell enclosing this location
+        //  The value (-1,-1,-1) is returned for out-of-bounds (not found).
+        labelVector findCell(const point& pt) const;
+
+
     // Mesh generation
 
         //- Create polyMesh for grid definition and patch information