ENH: additional constructor and methods for PDRblock (#1216)

- construct from positions

- reset positions, minEdgeLen, find grid index

ENH: add accessor for obtaining the size of a single i-j-k dimension

- eg, obtain the 'i' dimension with any of these methods:

      ijkAddr.size(vector::X)
      ijkAddr.sizes().x()
      ijkAddr.sizes()[0]

src/OpenFOAM/meshes/ijkMesh/IjkField.H

@@ -104,6 +104,12 @@ public:
         //- Return i,j,k addressing sizes for modification
         inline labelVector& sizes();
 
+        //- The field size
+        using Field<Type>::size;
+
+        //- The addressing dimension in the given direction
+        inline const label& size(const vector::components cmpt) const;
+
 
     // Edit
 

src/OpenFOAM/meshes/ijkMesh/IjkFieldI.H

@@ -157,6 +157,16 @@ inline Foam::labelVector& Foam::IjkField<Type>::sizes()
 }
 
 
+template<class Type>
+inline const Foam::label& Foam::IjkField<Type>::size
+(
+    const vector::components cmpt
+) const
+{
+    return ijk_.size(cmpt);
+}
+
+
 template<class Type>
 inline void Foam::IjkField<Type>::clear()
 {

src/OpenFOAM/meshes/ijkMesh/ijkAddressing.H

@@ -36,6 +36,7 @@ SourceFiles
 #define ijkAddressing_H
 
 #include "labelVector.H"
+#include "vector.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -72,17 +73,20 @@ public:
 
     // Access
 
-        //- Return the i,j,k addressing sizes
+        //- Addressing is considered empty if any component is zero
+        inline bool empty() const;
+
+        //- The (i,j,k) addressing dimensions
         inline const labelVector& sizes() const;
 
-        //- Return the i,j,k addressing sizes for modification
+        //- Return the (i,j,k) dimensions for modification
         inline labelVector& sizes();
 
         //- Return the total i*j*k size
         inline label size() const;
 
-        //- Addressing is considered empty if any component is zero
+        //- The addressing dimension in the given direction
-        inline bool empty() const;
+        inline const label& size(const vector::components cmpt) const;
 
         //- Reset to (0,0,0) sizing
         inline void clear();
@@ -93,13 +97,13 @@ public:
         //- Change the sizing parameters
         inline void reset(const labelVector& newSizes);
 
-        //- Linear addressing index (offset) for an i,j,k position.
+        //- Linear addressing index (offset) for an (i,j,k) position.
         inline label index(const label i, const label j, const label k) const;
 
-        //- Linear addressing index (offset) for an i,j,k position.
+        //- 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.
+        //- The (i,j,k) indexing from linear addressing.
         inline labelVector index(const label idx) const;
 
 

src/OpenFOAM/meshes/ijkMesh/ijkAddressingI.H

@@ -58,6 +58,12 @@ inline Foam::ijkAddressing::ijkAddressing
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+inline bool Foam::ijkAddressing::empty() const
+{
+    return (!sizes_.x() || !sizes_.y() || !sizes_.z());
+}
+
+
 inline const Foam::labelVector& Foam::ijkAddressing::sizes() const
 {
     return sizes_;
@@ -77,9 +83,12 @@ inline Foam::label Foam::ijkAddressing::size() const
 }
 
 
-inline bool Foam::ijkAddressing::empty() const
+inline const Foam::label& Foam::ijkAddressing::size
+(
+    const vector::components cmpt
+) const
 {
-    return (!sizes_.x() || !sizes_.y() || !sizes_.z());
+    return sizes_[cmpt];
 }
 
 

src/OpenFOAM/meshes/ijkMesh/ijkMeshI.H

@@ -143,13 +143,13 @@ inline Foam::label Foam::ijkMesh::nBoundaryFaces
             return n.y()*n.z();
             break;
 
-        // Face 2, 3 == y-min, y-max
+        // Face 2,3 == y-min, y-max
         case 2:
         case 3:
             return n.z()*n.x();
             break;
 
-        // Face 4, 5 == z-min, z-max
+        // Face 4,5 == z-min, z-max
         case 4:
         case 5:
             return n.x()*n.y();

src/mesh/blockMesh/Make/files

@@ -41,5 +41,6 @@ blockMeshTools/blockMeshTools.C
 
 PDRblockMesh/PDRblock.C
 PDRblockMesh/PDRblockCreate.C
+PDRblockMesh/PDRblockLocation.C
 
 LIB = $(FOAM_LIBBIN)/libblockMesh

src/mesh/blockMesh/PDRblockMesh/PDRblock.C

@@ -39,8 +39,87 @@ namespace Foam
 }
 
 
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+bool Foam::PDRblock::checkMonotonic
+(
+    const direction cmpt,
+    const UList<scalar>& 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;
+}
+
+
 // * * * * * * * * * * * * * 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,
@@ -89,22 +168,7 @@ void Foam::PDRblock::readGridControl
     }
 
     // Do points increase monotonically?
-    {
+    checkMonotonic(cmpt, knots);
-        const scalar& minVal = knots.first();
-
-        for (label pointi = 1; pointi < knots.size(); ++pointi)
-        {
-            if (knots[pointi] <= minVal)
-            {
-                FatalIOErrorInFunction(dict)
-                    << "Points are not monotonically increasing:"
-                    << " in " << vector::componentNames[cmpt]
-                    << " direction" << nl
-                    << exit(FatalIOError);
-            }
-        }
-    }
-
 
     if (verbose_)
     {
@@ -394,6 +458,33 @@ void Foam::PDRblock::readBoundary(const dictionary& dict)
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+Foam::PDRblock::PDRblock
+(
+    const UList<scalar>& xgrid,
+    const UList<scalar>& ygrid,
+    const UList<scalar>& 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()
@@ -411,50 +502,124 @@ bool Foam::PDRblock::read(const dictionary& dict)
     // Mark no scaling with invalid value
     const scalar scaleFactor = dict.lookupOrDefault<scalar>("scale", -1);
 
-    // Grid controls
     readGridControl(0, dict.subDict("x"), scaleFactor);
     readGridControl(1, dict.subDict("y"), scaleFactor);
     readGridControl(2, dict.subDict("z"), scaleFactor);
 
-    // Adjust i-j-k addressing
+    adjustSizes();
-    sizes().x() = grid_.x().nCells();
-    sizes().y() = grid_.y().nCells();
-    sizes().z() = grid_.z().nCells();
 
-    // 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();
-
-
-    // Boundaries
     readBoundary(dict);
 
     return true;
 }
 
 
+void Foam::PDRblock::reset
+(
+    const UList<scalar>& xgrid,
+    const UList<scalar>& ygrid,
+    const UList<scalar>& zgrid
+)
+{
+    static_cast<scalarList&>(grid_.x()) = xgrid;
+    static_cast<scalarList&>(grid_.y()) = ygrid;
+    static_cast<scalarList&>(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
 {
-    // Out-of-bounds is handled explicitly, for efficiency and consistency,
+    labelVector pos;
-    // 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 (findCell(pt, pos))
-    if (bounds_.contains(pt))
     {
-        for (direction cmpt=0; cmpt < labelVector::nComponents; ++cmpt)
+        return pos;
-        {
-            // Binary search
-            pos[cmpt] = findLower(grid_[cmpt], pt[cmpt]);
-        }
     }
 
-    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);
 }
 
 

src/mesh/blockMesh/PDRblockMesh/PDRblock.H

@@ -44,10 +44,10 @@ Description
 
     Grid coordinate controls
     \table
-        Property| Description                          | Required | Default
+        Property| Description                               | Required | Default
-        points  | Locations defining the mesh segment  | yes |
+        points  | Locations defining the mesh segment       | yes |
-        nCells  | Divisions per mesh segment           | yes |
+        nCells  | Divisions per mesh segment                | yes |
-        ratios  | Expansion ratio (end/start) per segment | no | uniform
+        ratios  | Expansion ratio (end/start) per segment   | no  | uniform
     \endtable
 
     The expansion ratios are defined as per blockMesh and represent the ratio
@@ -98,14 +98,20 @@ public:
         :
             public scalarList
         {
+            //- The locations are valid if they contain 2 or more points
+            inline bool valid() const;
+
             //- The number of cells in this direction.
             inline label nCells() const;
 
             //- The number of points in this direction.
             inline label nPoints() const;
 
-            //- The number of divisions (nCells) in this direction.
+            //- True if the location is within the range
-            inline label size() const;
+            inline bool contains(const scalar p) const;
+
+            //- Mid-point location, zero for an empty list.
+            inline scalar centre() const;
 
             //- Check that element index is within valid range.
             inline void checkIndex(const label i) const;
@@ -114,7 +120,22 @@ public:
             inline scalar width(const label i) const;
 
             //- Cell centre at element position.
+            //  Treats -1 and nCells positions like a halo cell.
             inline scalar C(const label i) const;
+
+            //- Find the cell index enclosing this location
+            //  \return -1 for out-of-bounds
+            label findCell(const scalar p) const;
+
+            //- Find the grid index, within the given tolerance
+            //  Return -1 for out-of-bounds and -2 for a point that is
+            //  within bounds, but not aligned with a grid point.
+            label findIndex(const scalar p, const scalar tol) const;
+
+            //- If out of range, return the respective min/max limits,
+            //- otherwise return the value itself.
+            //  If the range is invalid, always return the value.
+            inline const scalar& clip(const scalar& val) const;
         };
 
 
@@ -153,9 +174,18 @@ private:
         //- The boundary patch information
         PtrList<boundaryEntry> patches_;
 
+        //- The min edge length
+        scalar minEdgeLen_;
+
 
     // Private Member Functions
 
+        //- Check that points increase monotonically
+        static bool checkMonotonic(const direction cmpt, const UList<scalar>& pts);
+
+        //- Adjust sizing for updated grid points
+        void adjustSizes();
+
         //- Read and define grid points in given direction
         void readGridControl
         (
@@ -167,7 +197,6 @@ private:
         //- Read "boundary" information
         void readBoundary(const dictionary& dict);
 
-
         //- Populate point field for the block
         void createPoints(pointField& pts) const;
 
@@ -181,7 +210,6 @@ private:
             labelList::iterator& neiIter
         ) const;
 
-
         //- Add boundary faces for the shape face to lists
         //  Lists must be properly sized!
         //  \return the number of faces added
@@ -192,6 +220,19 @@ private:
             labelList::iterator& ownIter
         ) const;
 
+        //- Obtain i,j,k index for cell enclosing this location
+        //  \return false for out-of-bounds
+        bool findCell(const point& pt, labelVector& pos) const;
+
+        //- Obtain i,j,k grid index for point location
+        //  \return false for out-of-bounds and off-grid
+        bool gridIndex
+        (
+            const point& pt,
+            labelVector& pos,
+            const scalar tol
+        ) const;
+
 
 public:
 
@@ -200,6 +241,14 @@ public:
         //- Construct zero-size
         inline PDRblock();
 
+        //- Construct from components
+        PDRblock
+        (
+            const UList<scalar>& xgrid,
+            const UList<scalar>& ygrid,
+            const UList<scalar>& zgrid
+        );
+
         //- Construct from dictionary
         explicit PDRblock(const dictionary& dict, bool verbose=false);
 
@@ -209,6 +258,14 @@ public:
         //- Read dictionary
         bool read(const dictionary& dict);
 
+        //- Reset grid locations and mesh i-j-k sizing
+        void reset
+        (
+            const UList<scalar>& xgrid,
+            const UList<scalar>& ygrid,
+            const UList<scalar>& zgrid
+        );
+
 
     // Access
 
@@ -221,6 +278,9 @@ public:
         //- The mesh bounding box
         inline const boundBox& bounds() const;
 
+        //- The min edge length
+        inline const scalar& minEdgeLen() const;
+
         //- Cell size in x-direction at i position.
         inline scalar dx(const label i) const;
 
@@ -272,6 +332,12 @@ public:
         //  The value (-1,-1,-1) is returned for out-of-bounds (not found).
         labelVector findCell(const point& pt) const;
 
+        //- Obtain i,j,k grid index for point location within specified
+        //  relative tolerance of the minEdgeLen.
+        //  The value (-1,-1,-1) is returned for out-of-bounds (not found).
+        //  and off-grid
+        labelVector gridIndex(const point& pt, const scalar relTol=0.01) const;
+
 
     // Mesh generation
 

src/mesh/blockMesh/PDRblockMesh/PDRblockI.H

@@ -31,12 +31,19 @@ inline Foam::PDRblock::PDRblock()
     verbose_(false),
     grid_(),
     bounds_(),
-    patches_()
+    patches_(),
+    minEdgeLen_(Zero)
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+inline bool Foam::PDRblock::location::valid() const
+{
+    return (scalarList::size() > 1);
+}
+
+
 inline Foam::label Foam::PDRblock::location::nCells() const
 {
     return (scalarList::size()-1);
@@ -49,19 +56,25 @@ inline Foam::label Foam::PDRblock::location::nPoints() const
 }
 
 
-inline Foam::label Foam::PDRblock::location::size() const
+inline bool Foam::PDRblock::location::contains(const scalar p) const
 {
-    return (scalarList::size()-1);
+    return (scalarList::size() > 1 && first() <= p && p <= last());
+}
+
+
+inline Foam::scalar Foam::PDRblock::location::centre() const
+{
+    return scalarList::empty() ? 0 : (0.5*first() + 0.5*last());
 }
 
 
 inline void Foam::PDRblock::location::checkIndex(const label i) const
 {
-    if (i < 0 || i >= size())
+    if (i < 0 || i >= nCells())
     {
         FatalErrorInFunction
             << "The index " << i
-            << " is out of range [0," << size() << ']' << nl
+            << " is out of range [0," << nCells() << ']' << nl
             << abort(FatalError);
     }
 }
@@ -79,11 +92,45 @@ inline Foam::scalar Foam::PDRblock::location::width(const label i) const
 
 inline Foam::scalar Foam::PDRblock::location::C(const label i) const
 {
+    if (i == -1)
+    {
+        #ifdef FULLDEBUG
+        checkIndex(0);
+        #endif
+
+        // "Halo" centre [-1] == x0 - 1/2 (x1 - x0)
+        return first() - 0.5*(operator[](1) - first());
+    }
+    else if (i > 1 && i == scalarList::size()-1)
+    {
+        // "Halo" centre [nCells] == xN + 1/2 (xN - xN1)
+        return last() + 0.5*(last() - operator[](scalarList::size()-2));
+    }
+
     #ifdef FULLDEBUG
     checkIndex(i);
     #endif
 
-    return (operator[](i+1) + operator[](i));
+    return 0.5*(operator[](i+1) + operator[](i));
+}
+
+
+inline const Foam::scalar&
+Foam::PDRblock::location::clip(const scalar& val) const
+{
+    if (scalarList::size())
+    {
+        if (val < first())
+        {
+            return first();
+        }
+        else if (last() < val)
+        {
+            return last();
+        }
+    }
+
+    return val; // Pass-through
 }
 
 
@@ -94,6 +141,12 @@ Foam::PDRblock::grid() const
 }
 
 
+inline const Foam::scalar& Foam::PDRblock::minEdgeLen() const
+{
+    return minEdgeLen_;
+}
+
+
 inline const Foam::boundBox& Foam::PDRblock::bounds() const
 {
     return bounds_;

src/mesh/blockMesh/PDRblockMesh/PDRblockLocation.C

@@ -0,0 +1,108 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PDRblock.H"
+#include "ListOps.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::label Foam::PDRblock::location::findCell(const scalar p) const
+{
+    if (scalarList::empty())
+    {
+        return -1;
+    }
+    else if (equal(p, first()))
+    {
+        return 0;
+    }
+    else if (equal(p, last()))
+    {
+        return nCells()-1;
+    }
+    else if (p < first() || p > last())
+    {
+        // The point is out-of-bounds
+        return -1;
+    }
+
+    // Binary search, finds lower index and thus corresponds to the
+    // cell in which the point is found
+    return findLower(*this, p);
+}
+
+
+Foam::label Foam::PDRblock::location::findIndex
+(
+    const scalar p,
+    const scalar tol
+) const
+{
+    if (scalarList::empty())
+    {
+        return -1;
+    }
+    else if (Foam::mag(p - first()) <= tol)
+    {
+        return 0;
+    }
+    else if (Foam::mag(p - last()) <= tol)
+    {
+        return scalarList::size()-1;
+    }
+    else if (p < first() || p > last())
+    {
+        // The point is out-of-bounds
+        return -1;
+    }
+
+    // Linear search
+    label i = 0;
+    scalar delta = GREAT;
+
+    for (const scalar& val : *this)
+    {
+        const scalar diff = mag(p - val);
+
+        if (diff <= tol)
+        {
+            return i;
+        }
+        else if (delta < diff)
+        {
+            // Moving further away
+            break;
+        }
+
+        delta = diff;
+        ++i;
+    }
+
+    // This point is within bounds, but not near a grid-point
+    return -2;
+}
+
+
+// ************************************************************************* //