ENH: add treatment for PDRblockMesh outer region expansions (#1906)

// Treatment of the outer region
   outer
   {
       type        sphere;
       onGround    true;
       expansion   relative;

       ratios      1.1;

       size        3;
       nCells      10;
   }

applications/utilities/mesh/generation/PDRblockMesh/PDRblockMesh.C

@@ -57,7 +57,6 @@ Usage
 #include "Time.H"
 #include "IOdictionary.H"
 #include "OSspecific.H"
-#include "OFstream.H"
 
 using namespace Foam;
 
@@ -97,6 +96,22 @@ int main(int argc, char *argv[])
     );
     argList::addOptionCompat("no-clean", {"noClean", -2006});
 
+    argList::addBoolOption
+    (
+        "no-outer",
+        "Create without any other region"
+    );
+    argList::addBoolOption
+    (
+        "print-dict",
+        "Print blockMeshDict equivalent and exit"
+    );
+    argList::addBoolOption
+    (
+        "write-dict",
+        "Write system/blockMeshDict.PDRblockMesh and exit"
+    );
+
     argList::addOption("dict", "file", "Alternative PDRblockMeshDict");
     argList::addOption
     (
@@ -111,6 +126,9 @@ int main(int argc, char *argv[])
     // Remove old files, unless disabled
     const bool removeOldFiles = !args.found("no-clean");
 
+    // Suppress creation of the outer region
+    const bool noOuterRegion = args.found("no-outer");
+
     const word regionName(polyMesh::defaultRegion);
     const word regionPath;
 
@@ -149,8 +167,38 @@ int main(int argc, char *argv[])
     Info<< "Creating PDRblockMesh from "
         << runTime.relativePath(dictIO.objectPath()) << endl;
 
+    // Always start from a PDRblock
     PDRblock blkMesh(meshDict, true);
 
+    if (args.found("print-dict"))
+    {
+        Info<< nl << "Equivalent blockMeshDict" << nl << nl;
+
+        blkMesh.blockMeshDict(Info, true);
+
+        Info<< "\nEnd\n" << endl;
+        return 0;
+    }
+
+    if (args.found("write-dict"))
+    {
+        // Generate system/blockMeshDict and exit
+        blkMesh.writeBlockMeshDict
+        (
+            IOobject
+            (
+                "blockMeshDict.PDRblockMesh",
+                runTime.system(),   // instance
+                runTime,            // registry
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false  // Do not register
+            )
+        );
+
+        Info<< "\nEnd\n" << endl;
+        return 0;
+    }
 
     // Instance for resulting mesh
     if (useTime)
@@ -170,14 +218,19 @@ int main(int argc, char *argv[])
 
 
     Info<< nl << "Creating polyMesh from PDRblockMesh" << endl;
+    if (noOuterRegion)
+    {
+        Info<< "Outer region disabled, using ijk generation" << nl;
+    }
 
     autoPtr<polyMesh> meshPtr =
-        blkMesh.mesh
     (
-            IOobject(regionName, meshInstance, runTime)
+        args.found("no-outer")
+      ? blkMesh.innerMesh(IOobject(regionName, meshInstance, runTime))
+      : blkMesh.mesh(IOobject(regionName, meshInstance, runTime))
     );
 
-    const polyMesh& mesh = *meshPtr;
+    polyMesh& mesh = *meshPtr;
 
     // Set the precision of the points data to 10
     IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));

applications/utilities/preProcessing/PDR/pdrFields/PDRmeshArrays.C

@@ -94,7 +94,7 @@ void Foam::PDRmeshArrays::classify
         << "  nPoints:" << pdrBlock.nPoints()
         << "  nCells:" << pdrBlock.nCells()
         << "  nFaces:" << pdrBlock.nFaces() << nl
-        << "  min-edge:" << pdrBlock.minEdgeLen() << nl;
+        << "  min-edge:" << pdrBlock.edgeLimits().min() << nl;
 
     Info<< "Classifying ijk indexing... " << nl;
 

etc/caseDicts/annotated/PDRblockMeshDict

@@ -0,0 +1,97 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2006                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      PDRblockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//- Point scaling (optional)
+scale   1;
+
+x
+{
+    points  ( -13.28 -0.10 6.0 19.19 );
+    nCells  ( 10 12 10 );
+    ratios  ( -5.16 1 5.16 );
+}
+
+y
+{
+    points  ( -12.98 0 5.50 18.48 );
+    nCells  ( 10 11 10 );
+    ratios  ( -5.16 1 5.16 );
+}
+
+z
+{
+    points  ( 0.00 4.80 17.26 );
+    nCells  ( 10 10 );
+    ratios  ( 1 5.16 );
+}
+
+
+defaultPatch
+{
+    name    walls;
+    type    wall;
+}
+
+
+// Faces: 0=x-min, 1=x-max, 2=y-min, 3=y-max, 4=z-min, 5=z-max
+
+boundary
+(
+    outer
+    {
+        type    patch;
+        faces   ( 0 1 2 3 5 );
+    }
+
+    mergingFaces
+    {
+        type    wall;
+        faces   ();
+    }
+
+    blockedFaces
+    {
+        type    wall;
+        faces   ();
+    }
+
+    ground
+    {
+        type    wall;
+        faces   ( 4 );
+    }
+);
+
+
+// Treatment of the outer region
+
+outer
+{
+    type        sphere;     // (none | extend | box | sphere)  [default: none]
+    onGround    true;       // Module on the ground?  [default: false]
+    expansion   relative;   // (uniform | ratio | relative)  [default: ratio]
+
+    ratios      1.1;
+
+    size        3;          // Overall outer/inner size
+    nCells      10;         // Number of cells for outer region
+
+    // size       (3 4);
+    // nCells     (10 12);
+}
+
+
+// ************************************************************************* //

src/mesh/blockMesh/Make/files

@@ -40,7 +40,9 @@ blockMesh/blockMeshMergeTopological.C
 blockMeshTools/blockMeshTools.C
 
 PDRblockMesh/PDRblock.C
+PDRblockMesh/PDRblockBlockMesh.C
 PDRblockMesh/PDRblockCreate.C
 PDRblockMesh/PDRblockLocation.C
+PDRblockMesh/PDRblockOuter.C
 
 LIB = $(FOAM_LIBBIN)/libblockMesh

src/mesh/blockMesh/PDRblockMesh/PDRblock.C

@@ -28,6 +28,7 @@ License
 #include "PDRblock.H"
 #include "ListOps.H"
 #include "emptyPolyPatch.H"
+#include "gradingDescriptors.H"
 
 // Output when verbosity is enabled
 #undef  Log
@@ -52,21 +53,22 @@ Foam::PDRblock::expansionNames_
 
 namespace Foam
 {
-    //- Calculate geometric expansion factor from expansion ratio
-    inline scalar calcGexp(const scalar expRatio, const label nDiv)
-    {
-        return nDiv > 1 ? pow(expRatio, 1.0/(nDiv - 1)) : 0.0;
-    }
 
-    //- Calculate geometric ratio from relative ratio
+//- Calculate geometric expansion factor from expansion ratio
-    inline scalar relativeToGeometricRatio
+inline scalar calcGexp(const scalar expRatio, const label nDiv)
-    (
+{
+    return nDiv > 1 ? pow(expRatio, 1.0/(nDiv - 1)) : 0.0;
+}
+
+//- Calculate geometric ratio from relative ratio
+inline scalar relativeToGeometricRatio
+(
     const scalar expRatio,
     const label nDiv
-    )
+)
-    {
+{
     return nDiv > 1 ? pow(expRatio, (nDiv - 1)) : 1.0;
-    }
+}
 
 } // End namespace Foam
 
@@ -129,32 +131,20 @@ void Foam::PDRblock::adjustSizes()
         grid_.z().clear();
 
         bounds_ = boundBox::invertedBox;
-        minEdgeLen_ = Zero;
+        edgeLimits_.min() = 0;
+        edgeLimits_.max() = 0;
         return;
     }
 
     // Adjust boundBox
-    bounds_.min().x() = grid_.x().first();
+    bounds_ = bounds(grid_.x(), grid_.y(), grid_.z());
-    bounds_.min().y() = grid_.y().first();
-    bounds_.min().z() = grid_.z().first();
 
-    bounds_.max().x() = grid_.x().last();
+    // Min/max edge lengths
-    bounds_.max().y() = grid_.y().last();
+    edgeLimits_.clear();
-    bounds_.max().z() = grid_.z().last();
 
-    // Min edge length
+    edgeLimits_.add(grid_.x().edgeLimits());
-    minEdgeLen_ = GREAT;
+    edgeLimits_.add(grid_.y().edgeLimits());
-
+    edgeLimits_.add(grid_.z().edgeLimits());
-    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);
-        }
-    }
 }
 
 
@@ -166,14 +156,18 @@ void Foam::PDRblock::readGridControl
     expansionType expandType
 )
 {
-    //- The begin/end nodes for each segment
+    gridControl& ctrl = control_[cmpt];
-    scalarList knots;
+
+    // Begin/end nodes for each segment
+    scalarList& knots = static_cast<scalarList&>(ctrl);
 
     // The number of division per segment
-    labelList divisions;
+    labelList& divisions = ctrl.divisions_;
 
     // The expansion ratio per segment
-    scalarList expansion;
+    scalarList& expansion = ctrl.expansion_;
+
+    expansion.clear();  // expansion is optional
 
     Log << "Reading grid control for "
         << vector::componentNames[cmpt] << " direction" << nl;
@@ -514,6 +508,12 @@ void Foam::PDRblock::readBoundary(const dictionary& dict)
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+Foam::PDRblock::PDRblock()
+:
+    PDRblock(dictionary::null, false)
+{}
+
+
 Foam::PDRblock::PDRblock
 (
     const UList<scalar>& xgrid,
@@ -521,7 +521,7 @@ Foam::PDRblock::PDRblock
     const UList<scalar>& zgrid
 )
 :
-    PDRblock()
+    PDRblock(dictionary::null, false)
 {
     // Default boundaries with patchi == shapeFacei
     patches_.resize(6);
@@ -541,13 +541,21 @@ Foam::PDRblock::PDRblock
 }
 
 
-Foam::PDRblock::PDRblock(const dictionary& dict, bool verbose)
+Foam::PDRblock::PDRblock(const dictionary& dict, bool verboseOutput)
 :
-    PDRblock()
+    ijkMesh(),
+    meshDict_(dict),
+    grid_(),
+    outer_(),
+    bounds_(),
+    patches_(),
+    edgeLimits_(0,0),
+    verbose_(verboseOutput)
 {
-    verbose_ = verbose;
+    if (&dict != &dictionary::null)
-
+    {
         read(dict);
+    }
 }
 
 
@@ -575,6 +583,16 @@ bool Foam::PDRblock::read(const dictionary& dict)
 
     readBoundary(dict);
 
+    // Outer treatment: (none | extend | box | sphere)
+    outer_.clear();
+
+    const dictionary* outerDictPtr = dict.findDict("outer");
+    if (outerDictPtr)
+    {
+        outer_.read(*outerDictPtr);
+    }
+    outer_.report(Info);
+
     return true;
 }
 
@@ -644,7 +662,7 @@ bool Foam::PDRblock::gridIndex
     const scalar relTol
 ) const
 {
-    const scalar tol = relTol * minEdgeLen_;
+    const scalar tol = relTol * edgeLimits_.min();
 
     for (direction cmpt=0; cmpt < labelVector::nComponents; ++cmpt)
     {
@@ -688,4 +706,33 @@ Foam::labelVector Foam::PDRblock::gridIndex
 }
 
 
+Foam::Vector<Foam::gradingDescriptors> Foam::PDRblock::grading() const
+{
+    return grading(control_);
+}
+
+
+Foam::gradingDescriptors Foam::PDRblock::grading(const direction cmpt) const
+{
+    switch (cmpt)
+    {
+        case vector::X :
+        case vector::Y :
+        case vector::Z :
+        {
+            return control_[cmpt].grading();
+            break;
+        }
+
+        default :
+            FatalErrorInFunction
+                << "Not gridControl for direction " << label(cmpt) << endl
+                << exit(FatalError);
+            break;
+    }
+
+    return gradingDescriptors();
+}
+
+
 // ************************************************************************* //

src/mesh/blockMesh/PDRblockMesh/PDRblock.H

@@ -73,6 +73,8 @@ SourceFiles
 #include "pointField.H"
 #include "faceList.H"
 #include "Enum.H"
+#include "vector2D.H"
+#include "labelVector2D.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -81,7 +83,9 @@ namespace Foam
 
 // Forward Declarations
 class IOobject;
+class blockMesh;
 class polyMesh;
+class gradingDescriptors;
 
 /*---------------------------------------------------------------------------*\
                            Class PDRblock Declaration
@@ -96,7 +100,7 @@ public:
     // Data Types
 
         //- The expansion type
-        enum expansionType
+        enum expansionType : uint8_t
         {
             EXPAND_UNIFORM,     //!< Uniform expansion (ie, no expansion)
             EXPAND_RATIO,       //!< End/start ratio
@@ -115,7 +119,7 @@ public:
         {
         public:
 
-            //- The locations are valid if they contain 2 or more points
+            //- The location list is valid if it contains 2 or more points
             inline bool valid() const;
 
             //- The number of cells in this direction.
@@ -136,6 +140,9 @@ public:
             //- Mid-point location, zero for an empty list.
             inline scalar centre() const;
 
+            //- The difference between min/max values, zero for an empty list.
+            inline scalar length() const;
+
             //- Check that element index is within valid range.
             inline void checkIndex(const label i) const;
 
@@ -146,6 +153,13 @@ public:
             //  Treats -1 and nCells positions like a halo cell.
             inline scalar C(const label i) const;
 
+            //- Return min/max edge lengths
+            scalarMinMax edgeLimits() const;
+
+            //- Return edge grading descriptors for the locations
+            //  \see Foam::gradingDescriptor
+            gradingDescriptors grading() const;
+
             //- Find the cell index enclosing this location
             //  \return -1 for out-of-bounds
             label findCell(const scalar p) const;
@@ -172,6 +186,33 @@ private:
 
     // Private Classes
 
+        //- The begin/end nodes for each segment,
+        //- with divisions and expansion for each segment
+        struct gridControl
+        :
+            public scalarList
+        {
+            //- The number of division per segment
+            labelList divisions_;
+
+            //- The expansion ratio per segment
+            scalarList expansion_;
+
+            //- Total number of cells in this direction
+            label nCells() const;
+
+            //- Return edge grading descriptors for the locations
+            //  \see Foam::gradingDescriptor
+            gradingDescriptors grading() const;
+
+            //- Add point/divisions/expand to end of list (push_back)
+            void append(const scalar p, label nDiv, scalar expRatio=1);
+
+            //- Add point/divisions/expand to front of list (push_front)
+            void prepend(const scalar p, label nDiv, scalar expRatio=1);
+        };
+
+
         //- Extracted patch settings
         struct boundaryEntry
         {
@@ -188,23 +229,100 @@ private:
             labelList faces_;
         };
 
+        //- The begin/end nodes for each segment,
+        //- with divisions and expansion for each segment
+        struct outerControl
+        {
+            //- The control type
+            enum controlType : uint8_t
+            {
+                OUTER_NONE = 0,     //!< No outer region
+                OUTER_EXTEND,       //!< Extend inner region (orthogonal)
+                OUTER_BOX,          //!< Cuboid
+                OUTER_SPHERE        //!< Spherical
+            };
+
+            //- Named enumerations for the control type
+            const static Enum<controlType> controlNames_;
+
+            //- The control type
+            controlType type_;
+
+            //- The expansion type
+            expansionType expandType_;
+
+            //- True if on the ground
+            bool onGround_;
+
+            //- Relative size(s) for the outer region
+            vector2D relSize_;
+
+            //- Number of cells in outer region
+            //  Generally only single component is used
+            labelVector2D nCells_;
+
+            //- Expansion ratio(s) for the outer region
+            vector2D expansion_;
+
+
+        // Constructors
+
+            //- Default construct. NONE
+            outerControl();
+
+
+        // Member Functions
+
+            //- Reset to default (NONE) values
+            void clear();
+
+            //- Is enabled (not NONE)
+            bool active() const;
+
+            //- Project on to sphere (is SPHERE)
+            bool isSphere() const;
+
+            //- Is the outer region on the ground?
+            bool onGround() const;
+
+            //- Define that the outer region is on the ground or not
+            //  \return the old value
+            bool onGround(const bool on);
+
+            //- Read content from dictionary
+            void read(const dictionary& dict);
+
+            //- Report information about outer region
+            void report(Ostream& os) const;
+        };
+
 
     // Private Data
 
-        //- Verbosity
+        //- Reference to mesh dictionary
-        bool verbose_;
+        const dictionary& meshDict_;
 
-        //- The grid points in all (i,j,k / x,y,z) directions.
+        //- The grid controls in (i,j,k / x,y,z) directions.
+        Vector<gridControl> control_;
+
+        //- The grid points in all (i,j,k / x,y,z) directions,
+        //- after applying the internal subdivisions.
         Vector<location> grid_;
 
+        //- Control for the outer-region (if any)
+        outerControl outer_;
+
         //- The mesh bounding box
         boundBox bounds_;
 
         //- The boundary patch information
         PtrList<boundaryEntry> patches_;
 
-        //- The min edge length
+        //- The min/max edge lengths
-        scalar minEdgeLen_;
+        scalarMinMax edgeLimits_;
+
+        //- Verbosity
+        bool verbose_;
 
 
     // Private Member Functions
@@ -267,6 +385,35 @@ private:
             const scalar tol
         ) const;
 
+        //- The bounding box of the grid points
+        static boundBox bounds
+        (
+            const scalarList& x, //!< X-points, monotonically increasing
+            const scalarList& y, //!< Y-points, monotonically increasing
+            const scalarList& z  //!< T-points, monotonically increasing
+        );
+
+        //- Equivalent edge grading descriptors in (x,y,z) directions.
+        static Vector<gradingDescriptors> grading
+        (
+            const Vector<gridControl>& ctrl
+        );
+
+        //- Mesh sizes based on the controls
+        static labelVector sizes
+        (
+            const Vector<gridControl>& ctrl
+        );
+
+
+    // Mesh Generation
+
+        //- Create a blockMesh
+        autoPtr<blockMesh> createBlockMesh(const IOobject& io) const;
+
+        //- Create polyMesh via blockMesh
+        autoPtr<polyMesh> meshBlockMesh(const IOobject& io) const;
+
 
 public:
 
@@ -278,8 +425,8 @@ public:
 
     // Constructors
 
-        //- Construct zero-size
+        //- Default construct, zero-size, inverted bounds etc
-        inline PDRblock();
+        PDRblock();
 
         //- Construct from components
         PDRblock
@@ -290,7 +437,7 @@ public:
         );
 
         //- Construct from dictionary
-        explicit PDRblock(const dictionary& dict, bool verbose=false);
+        explicit PDRblock(const dictionary& dict, bool verboseOutput=false);
 
 
     // Member Functions
@@ -312,14 +459,23 @@ public:
         //- The grid point locations in the i,j,k (x,y,z) directions.
         inline const Vector<location>& grid() const;
 
+        //- Equivalent edge grading descriptors in (x,y,z) directions.
+        Vector<gradingDescriptors> grading() const;
 
-    // Mesh information
+        //- Equivalent edge grading descriptors in specified (x,y,z) direction.
+        gradingDescriptors grading(const direction cmpt) const;
+
+
+    // Mesh Information
+
+        //- Mesh sizing as per ijkMesh
+        using ijkMesh::sizes;
 
         //- The mesh bounding box
         inline const boundBox& bounds() const;
 
-        //- The min edge length
+        //- The min/max edge length
-        inline const scalar& minEdgeLen() const;
+        inline const scalarMinMax& edgeLimits() const;
 
         //- Cell size in x-direction at i position.
         inline scalar dx(const label i) const;
@@ -379,17 +535,30 @@ public:
         labelVector findCell(const point& pt) const;
 
         //- Obtain i,j,k grid index for point location within specified
-        //  relative tolerance of the minEdgeLen.
+        //  relative tolerance of the min edge length
         //  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
+    // Mesh Generation
+
+        //- Output content for an equivalent blockMeshDict
+        //  Optionally generate header/footer content
+        Ostream& blockMeshDict(Ostream& os, const bool withHeader=false) const;
+
+        //- Content for an equivalent blockMeshDict
+        dictionary blockMeshDict() const;
+
+        //- Write an equivalent blockMeshDict
+        void writeBlockMeshDict(const IOobject& io) const;
 
         //- Create polyMesh for grid definition and patch information
         autoPtr<polyMesh> mesh(const IOobject& io) const;
 
+        //- Create polyMesh for inner-mesh only,
+        //- ignore any outer block definitions
+        autoPtr<polyMesh> innerMesh(const IOobject& io) const;
 };
 
 

src/mesh/blockMesh/PDRblockMesh/PDRblockBlockMesh.C

@@ -0,0 +1,920 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 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 "PDRblock.H"
+#include "ListOps.H"
+#include "cellModeller.H"
+#include "gradingDescriptors.H"
+#include "objectRegistry.H"
+#include "Time.H"
+#include "IOdictionary.H"
+#include "Fstream.H"
+#include "OTstream.H"
+#include "edgeHashes.H"
+
+#include "cellModel.H"
+#include "blockMesh.H"
+#include "polyMesh.H"
+#include "searchableSphere.H"
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Name for the projection geometry
+static const word projKeyword("project");
+
+// Name for the projection geometry
+static const word projGeomName("sphere");
+
+
+//- Calculate geometric ratio from relative ratio
+inline scalar relativeToGeometricRatio
+(
+    const scalar expRatio,
+    const label nDiv
+)
+{
+    return nDiv > 1 ? pow(expRatio, (nDiv - 1)) : 1.0;
+}
+
+
+// Output space-separated flat list. No size prefix.
+template<class T>
+static Ostream& outputFlatList(Ostream& os, const UList<T>& list)
+{
+    os  << token::BEGIN_LIST;
+    label i = 0;
+    for (const label val : list)
+    {
+        if (i++) os << token::SPACE;
+        os  << val;
+    }
+    os  << token::END_LIST;
+
+    return os;
+}
+
+
+// Begin indent list
+static inline Ostream& begIndentList(Ostream& os)
+{
+    os  << indent << incrIndent << token::BEGIN_LIST << nl;
+    return os;
+}
+
+// End indent list
+static inline Ostream& endIndentList(Ostream& os)
+{
+    os  << decrIndent << indent << token::END_LIST;
+    return os;
+}
+
+
+// Output list contents (newline separated) indented.
+template<class T>
+static Ostream& outputIndent(Ostream& os, const UList<T>& list)
+{
+    for (const T& val : list)
+    {
+        os  << indent << val << nl;
+    }
+    return os;
+}
+
+
+//
+static Ostream& serializeHex
+(
+    Ostream& os,
+    const labelUList& hexVerts,
+    const labelVector& hexCount,
+    const Vector<gradingDescriptors> hexGrade,
+    const word& zoneName = word::null
+)
+{
+    os  << indent << cellModel::modelNames[cellModel::HEX] << token::SPACE;
+    outputFlatList(os, hexVerts);
+
+    if (!zoneName.empty())
+    {
+        os  << token::SPACE << zoneName;
+    }
+
+    os  << token::SPACE << hexCount << nl
+        << indent << word("edgeGrading") << nl;
+
+    begIndentList(os);
+
+    // Grading (x/y/z)
+    for (const gradingDescriptors& gds : hexGrade)
+    {
+        begIndentList(os);
+        outputIndent(os, gds);
+        endIndentList(os) << nl;
+    }
+
+    endIndentList(os) << nl;
+    return os;
+}
+
+
+// Generate list with entries:
+//
+//    project (x y z) (geometry)
+//
+static Ostream& serializeProjectPoints
+(
+    Ostream& os,
+    const UList<point>& list
+)
+{
+    for (const point& p : list)
+    {
+        os  << indent << projKeyword << token::SPACE
+            << p
+            << token::SPACE
+            << token::BEGIN_LIST << projGeomName << token::END_LIST << nl;
+    }
+
+    return os;
+}
+
+
+// Generate entry:
+//
+//    project (beg end) (geometry)
+//
+static Ostream& serializeProjectEdge
+(
+    Ostream& os,
+    const edge& e
+)
+{
+    os  << indent << projKeyword << token::SPACE;
+
+    if (e.sorted())
+    {
+        os  << e.first() << token::SPACE << e.second();
+    }
+    else
+    {
+        os  << e.second() << token::SPACE << e.first();
+    }
+
+    os  << token::SPACE
+        << token::BEGIN_LIST << projGeomName << token::END_LIST << nl;
+
+    return os;
+}
+
+
+// Generate entry:
+//
+//    (0 1 2 ..)
+//
+static Ostream& serializeFace
+(
+    Ostream& os,
+    const face& list
+)
+{
+    os  << indent;
+    outputFlatList(os, list);
+    os  << nl;
+    return os;
+}
+
+
+// Generate entry:
+//
+//    project (0 1 2 ..) geometry
+//
+static Ostream& serializeProjectFace
+(
+    Ostream& os,
+    const face& list
+)
+{
+    os  << indent << projKeyword << token::SPACE;
+    outputFlatList(os, list);
+    os  << token::SPACE << projGeomName << nl;
+
+    return os;
+}
+
+} // namespace Foam
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::PDRblock::blockMeshDict
+(
+    Ostream& os,
+    const bool withHeader
+) const
+{
+    if (withHeader)
+    {
+        // Use dummy time for fake objectRegistry
+        autoPtr<Time> dummyTimePtr(Time::New());
+
+        IOdictionary iodict
+        (
+            IOobject
+            (
+                "blockMeshDict",
+                dummyTimePtr->system(),
+                *dummyTimePtr,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false  // no register
+            )
+        );
+
+        iodict.writeHeader(os);
+    }
+
+    const cellModel& hex = cellModel::ref(cellModel::HEX);
+
+    // The mesh topology will normally be an O-grid with a central (inner)
+    // block and 6 outer blocks.
+    // The inner block is described by (0 1 2 3 4 5 6 7).
+    // The additional points for the outer region: (8 9 19 11 12 13 14 15)
+
+    // The outer blocks will be addressed according to their
+    // placement w.r.t. the inner block, and defined such that they retain
+    // the same global orientation as the inner block.
+    // For example, face 0 of all blocks will be on the logical x-min side
+    // of the mesh.
+
+    // List of hex vertices, ordered with outer blocks first to allow
+    // direct addressing by their logical position.
+    const FixedList<labelList, 8>
+        hexVerts
+        ({
+            {8, 0,  3, 11, 12,  4,  7, 15},     // x-min block (face 0)
+            {1, 9, 10,  2,  5, 13, 14,  6},     // x-max block (face 1)
+            {8, 9,  1,  0, 12, 13,  5,  4},     // y-min block (face 2)
+            {3, 2, 10, 11,  7,  6, 14, 15},     // y-max block (face 3)
+            {8, 9, 10, 11,  0,  1,  2,  3},     // z-min block (face 4)
+            {4, 5,  6,  7, 12, 13, 14, 15},     // z-max block (face 5)
+            {0, 1,  2,  3,  4,  5,  6,  7},     // Inner box description
+            {8, 9, 10, 11, 12, 13, 14, 15},     // Outer box description
+        });
+
+    // The face or the logical block index (for the O-grid)
+    enum faceIndex
+    {
+        X_Min = 0,
+        X_Max = 1,
+        Y_Min = 2,
+        Y_Max = 3,
+        Z_Min = 4,
+        Z_Max = 5,
+        Inner_Block = 6,    // Inner block description
+        Outer_Block = 7,    // Outer bounding box description
+    };
+
+    // Lists of block/face for outside and ground faces
+    DynamicList<labelPair> outerFaces(8);
+    DynamicList<labelPair> groundFaces(8);
+
+    // We handle a few fixed topology configurations
+
+    enum blockTopologyType
+    {
+        INNER_ONLY = 0,   // No outer region
+        EXTENDED   = 1,   // Outer created by extending inner region
+        CLIP_BOTTOM = 5,  // Outer O-grid on 5 sides
+        FULL_OUTER = 6    // Outer O-grid on 6 sides
+    };
+
+
+
+    // Expansion ratios need conversion from relative to geometric
+    const bool useRelToGeom =
+        (expansionType::EXPAND_RATIO == outer_.expandType_);
+
+
+    // Physical dimensions
+
+    Vector<gridControl> ctrl(control_);
+    boundBox innerCorners(bounds(ctrl.x(), ctrl.y(), ctrl.z()));
+
+    boundBox outerCorners;
+
+
+    point radialCentre(innerCorners.centre());
+    vector radialSizes(0.5*innerCorners.span());
+
+    blockTopologyType outerTopology = INNER_ONLY;
+
+
+    if (outer_.active())
+    {
+        outerTopology = FULL_OUTER;
+
+        // Convert from relative size
+        radialSizes.x() *= outer_.relSize_.x();
+        radialSizes.y() *= outer_.relSize_.y();
+        radialSizes.z() *= min(outer_.relSize_.x(), outer_.relSize_.y());
+
+        if (outer_.onGround())
+        {
+            outerTopology = CLIP_BOTTOM;
+            radialCentre.z() = innerCorners.min().z();
+            radialSizes.z() *= 2;
+        }
+
+        // Corners for a box
+        outerCorners.min() = radialCentre - radialSizes;
+        outerCorners.max() = radialCentre + radialSizes;
+
+        if (outer_.onGround())
+        {
+            outerCorners.min().z() = innerCorners.min().z();
+        }
+
+
+        if (outer_.isSphere())
+        {
+            // For spheroid projection, don't trust that blockMesh does it
+            // properly. Give some reasonable estimates of the corners.
+
+            // Use dummy Time for objectRegistry
+            autoPtr<Time> dummyTimePtr(Time::New());
+
+            const IOobject io
+            (
+                "sphere",
+                *dummyTimePtr,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false  // do not register
+            );
+
+            searchableSphere sphere(io, radialCentre, radialSizes);
+
+            pointField queries(2);
+            queries[0] = outerCorners.min();
+            queries[1] = outerCorners.max();
+
+            List<pointIndexHit> hits;
+            sphere.findNearest
+            (
+                queries,
+                scalarField(2, GREAT),
+                hits
+            );
+
+            outerCorners.min() = hits[0].hitPoint();
+            outerCorners.max() = hits[1].hitPoint();
+        }
+        else if (outerControl::OUTER_EXTEND == outer_.type_)
+        {
+            outerTopology = EXTENDED;
+
+            // Extend the inner block
+            label outerCount;
+            scalar expRatio;
+
+            outerCount = outer_.nCells_.x();
+            expRatio = outer_.expansion_.x();
+            if (useRelToGeom)
+            {
+                expRatio = relativeToGeometricRatio(expRatio, outerCount);
+            }
+
+            ctrl.x().prepend(outerCorners.min().x(), outerCount, -expRatio);
+            ctrl.x().append(outerCorners.max().x(), outerCount, expRatio);
+
+
+            outerCount = outer_.nCells_.y();
+            expRatio = outer_.expansion_.y();
+            if (useRelToGeom)
+            {
+                expRatio = relativeToGeometricRatio(expRatio, outerCount);
+            }
+
+            ctrl.y().prepend(outerCorners.min().y(), outerCount, -expRatio);
+            ctrl.y().append(outerCorners.max().y(), outerCount, expRatio);
+
+            outerCount = max(outer_.nCells_.x(), outer_.nCells_.y());
+            expRatio = min(outer_.expansion_.x(), outer_.expansion_.y());
+            if (useRelToGeom)
+            {
+                expRatio = relativeToGeometricRatio(expRatio, outerCount);
+            }
+
+            if (!outer_.onGround())
+            {
+                ctrl.z().prepend(outerCorners.min().z(), outerCount, -expRatio);
+            }
+            ctrl.z().append(outerCorners.max().z(), outerCount, expRatio);
+
+            // Update corners
+            innerCorners = bounds(ctrl.x(), ctrl.y(), ctrl.z());
+            outerCorners = innerCorners;
+        }
+    }
+
+
+    const Vector<gradingDescriptors> innerGrading(grading(ctrl));
+    const labelVector innerCount(sizes(ctrl));
+
+    labelVector hexCount;
+    Vector<gradingDescriptors> hexGrade;
+
+
+    const label radialCount = outer_.nCells_.x();
+    scalar expRatio = outer_.expansion_.x();
+
+    if (useRelToGeom)
+    {
+        expRatio = relativeToGeometricRatio(expRatio, radialCount);
+    }
+
+    const gradingDescriptors radialInward
+    (
+        gradingDescriptor{-expRatio}
+    );
+
+    const gradingDescriptors radialOutward
+    (
+        gradingDescriptor{expRatio}
+    );
+
+
+    if (EXTENDED == outerTopology)
+    {
+        // The inner block is extended to become the outer faces
+        outerFaces.append
+        ({
+            labelPair(Inner_Block, X_Min),
+            labelPair(Inner_Block, X_Max),
+            labelPair(Inner_Block, Y_Min),
+            labelPair(Inner_Block, Y_Max),
+            labelPair(Inner_Block, Z_Max)
+        });
+
+        // The ground faces vs outside faces
+        if (outer_.onGround())
+        {
+            groundFaces.append
+            (
+                labelPair(Inner_Block, Z_Min)
+            );
+        }
+        else
+        {
+            outerFaces.append
+            (
+                labelPair(Inner_Block, Z_Min)
+            );
+        }
+    }
+    else if (CLIP_BOTTOM == outerTopology || FULL_OUTER == outerTopology)
+    {
+        // The outside faces
+        outerFaces.append
+        ({
+            labelPair(X_Min, X_Min),
+            labelPair(X_Max, X_Max),
+            labelPair(Y_Min, Y_Min),
+            labelPair(Y_Max, Y_Max),
+            labelPair(Z_Max, Z_Max)
+        });
+
+        // The ground faces
+        if (CLIP_BOTTOM == outerTopology)
+        {
+            groundFaces.append
+            ({
+                labelPair(X_Min, Z_Min),
+                labelPair(X_Max, Z_Min),
+                labelPair(Y_Min, Z_Min),
+                labelPair(Y_Max, Z_Min),
+                // Note: {Z_Min, Z_Min} will not exist
+                labelPair(Inner_Block, Z_Min)
+            });
+        }
+        else
+        {
+            outerFaces.append
+            (
+                labelPair(Z_Min, Z_Min)
+            );
+        }
+    }
+
+
+    if (outer_.isSphere())
+    {
+        os.beginBlock("geometry");
+        {
+            os.beginBlock(projGeomName);
+            {
+                os.writeEntry("type", "sphere");
+                os.writeEntry("origin", radialCentre);
+                os.writeEntry("radius", radialSizes);
+            }
+            os.endBlock();
+        }
+        os.endBlock();
+    }
+
+    // vertices
+    {
+        os  << nl << word("vertices") << nl;
+        begIndentList(os);
+
+        pointField corners(innerCorners.points());
+
+        // inner
+        outputIndent(os, corners);
+
+        // outer
+        if (CLIP_BOTTOM == outerTopology || FULL_OUTER == outerTopology)
+        {
+            corners = outerCorners.points();
+
+            if (outer_.isSphere())
+            {
+                serializeProjectPoints(os, corners);
+            }
+            else
+            {
+                outputIndent(os, corners);
+            }
+        }
+
+        endIndentList(os) << token::END_STATEMENT << nl;
+    }
+
+
+    // blocks
+    {
+        word innerZoneName = "inner";
+        if (INNER_ONLY == outerTopology || EXTENDED == outerTopology)
+        {
+            innerZoneName.clear();
+        }
+
+        os  << nl << word("blocks") << nl;
+        begIndentList(os);
+
+        // Inner block
+        hexCount = innerCount;
+
+        serializeHex
+        (
+            os,
+            hexVerts[Inner_Block],
+            hexCount,
+            innerGrading,
+            innerZoneName
+        );
+
+        // outer
+        if (CLIP_BOTTOM == outerTopology || FULL_OUTER == outerTopology)
+        {
+            // Radial direction = X
+            hexCount = innerCount;
+            hexGrade = innerGrading;
+
+            hexCount.x() = radialCount;
+
+            // Face 0: x-min
+            {
+                hexGrade.x() = radialInward;
+                serializeHex(os, hexVerts[X_Min], hexCount, hexGrade);
+            }
+            // Face 1: x-max
+            {
+                hexGrade.x() = radialOutward;
+                serializeHex(os, hexVerts[X_Max], hexCount, hexGrade);
+            }
+
+
+            // Radial direction = Y
+            hexCount = innerCount;
+            hexGrade = innerGrading;
+
+            hexCount.y() = radialCount;
+
+            // Face 2: y-min
+            {
+                hexGrade.y() = radialInward;
+                serializeHex(os, hexVerts[Y_Min], hexCount, hexGrade);
+            }
+            // Face 3: y-max
+            {
+                hexGrade.y() = radialOutward;
+                serializeHex(os, hexVerts[Y_Max], hexCount, hexGrade);
+            }
+
+
+            // Radial direction = Z
+            hexCount = innerCount;
+            hexGrade = innerGrading;
+
+            hexCount.z() = radialCount;
+
+            // Face 4: z-min
+            if (!outer_.onGround())
+            {
+                hexGrade.z() = radialInward;
+                serializeHex(os, hexVerts[Z_Min], hexCount, hexGrade);
+            }
+            // Face 5: z-max
+            {
+                hexGrade.z() = radialOutward;
+                serializeHex(os, hexVerts[Z_Max], hexCount, hexGrade);
+            }
+        }
+
+        endIndentList(os) << token::END_STATEMENT << nl;
+    }
+
+
+    // edges
+    {
+        os  << nl << word("edges") << nl;
+        begIndentList(os);
+
+        if (outer_.isSphere() && outerFaces.size())
+        {
+            // Edges for the outer face of the block
+            edgeHashSet projEdges(32);
+
+            for (const labelPair& pr : outerFaces)
+            {
+                projEdges.insert
+                (
+                    hex.face(pr.second(), hexVerts[pr.first()]).edges()
+                );
+            }
+
+            for (const edge& e : projEdges.sortedToc())
+            {
+                serializeProjectEdge(os, e);
+            }
+        }
+
+        endIndentList(os) << token::END_STATEMENT << nl;
+    }
+
+
+    // faces
+    {
+        os  << nl << word("faces") << nl;
+        begIndentList(os);
+
+        if (outer_.isSphere() && outerFaces.size())
+        {
+            for (const labelPair& pr : outerFaces)
+            {
+                serializeProjectFace
+                (
+                    os,
+                    hex.face(pr.second(), hexVerts[pr.first()])
+                );
+            }
+        }
+
+        endIndentList(os) << token::END_STATEMENT << nl;
+    }
+
+
+    // boundary
+    {
+        os  << nl << word("boundary") << nl;
+        begIndentList(os);
+
+        // outer
+        {
+            os.beginBlock("outer");
+            os.writeEntry("type", word("patch"));
+
+            os  << indent << word("faces") << nl;
+            begIndentList(os);
+
+            for (const labelPair& pr : outerFaces)
+            {
+                serializeFace
+                (
+                    os,
+                    hex.face(pr.second(), hexVerts[pr.first()])
+                );
+            }
+
+            endIndentList(os) << token::END_STATEMENT << nl;
+
+            os.endBlock();
+        }
+
+        if (outer_.onGround())
+        {
+            os.beginBlock("ground");
+            os.writeEntry("type", word("wall"));
+
+            os  << indent << word("faces") << nl;
+            begIndentList(os);
+
+            for (const labelPair& pr : groundFaces)
+            {
+                serializeFace
+                (
+                    os,
+                    hex.face(pr.second(), hexVerts[pr.first()])
+                );
+            }
+
+            endIndentList(os) << token::END_STATEMENT << nl;
+            os.endBlock();
+        }
+
+        endIndentList(os) << token::END_STATEMENT << nl;
+    }
+
+
+    if (withHeader)
+    {
+        IOobject::writeEndDivider(os);
+    }
+
+    return os;
+}
+
+
+Foam::dictionary Foam::PDRblock::blockMeshDict() const
+{
+    OTstream os;
+    blockMeshDict(os);
+
+    ITstream is("blockMeshDict", tokenList());
+    is.transfer(os.tokens());
+
+    return dictionary(is);
+}
+
+
+void Foam::PDRblock::writeBlockMeshDict(const IOobject& io) const
+{
+    IOdictionary iodict
+    (
+        IOobject
+        (
+            io.name(),
+            io.db().time().system(),
+            io.local(),
+            io.db(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false  // no register
+        )
+    );
+
+    OFstream os(iodict.objectPath());
+
+    Info<< nl
+        << "Generate blockMeshDict: "
+        << iodict.db().time().relativePath(os.name()) << endl;
+
+    // Set precision for points to 10
+    os.precision(max(10u, IOstream::defaultPrecision()));
+
+    iodict.writeHeader(os);
+
+    // Just like writeData, but without copying beforehand
+    this->blockMeshDict(os);
+
+    iodict.writeEndDivider(os);
+}
+
+
+Foam::autoPtr<Foam::blockMesh>
+Foam::PDRblock::createBlockMesh(const IOobject& io) const
+{
+    IOdictionary iodict
+    (
+        IOobject
+        (
+            "blockMeshDict.PDRblockMesh",
+            io.db().time().system(),
+            io.local(),
+            io.db(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false  // no register
+        ),
+        blockMeshDict()
+    );
+
+    return autoPtr<blockMesh>::New(iodict);
+}
+
+
+Foam::autoPtr<Foam::polyMesh>
+Foam::PDRblock::meshBlockMesh(const IOobject& io) const
+{
+    const bool oldVerbose = blockMesh::verboseOutput;
+    blockMesh::verboseOutput = false;
+
+    autoPtr<polyMesh> meshPtr(createBlockMesh(io)->mesh(io));
+
+    blockMesh::verboseOutput = oldVerbose;
+
+    // This is a bit ugly.
+    // For extend, we still wish to have an 'inner' cellZone,
+    // but we meshed the entirety.
+
+    if
+    (
+        outerControl::OUTER_EXTEND == outer_.type_
+     && meshPtr->cellZones().empty()
+    )
+    {
+        const boundBox innerBox
+        (
+            bounds(control_.x(), control_.y(), control_.z())
+        );
+
+        const label nZoneCellsMax =
+        (
+            control_.x().nCells()
+          * control_.y().nCells()
+          * control_.z().nCells()
+        );
+
+
+        polyMesh& pmesh = *meshPtr;
+
+        List<cellZone*> cz(1);
+        cz[0] = new cellZone
+        (
+            "inner",
+            labelList(nZoneCellsMax),
+            0, // zonei
+            pmesh.cellZones()
+        );
+
+        cellZone& innerZone = *(cz[0]);
+
+        const vectorField& cc = pmesh.cellCentres();
+
+        label nZoneCells = 0;
+
+        for
+        (
+            label celli = 0;
+            celli < cc.size() && nZoneCells < nZoneCellsMax;
+            ++celli
+        )
+        {
+            if (innerBox.contains(cc[celli]))
+            {
+                innerZone[nZoneCells] = celli;
+                ++nZoneCells;
+            }
+        }
+
+        innerZone.resize(nZoneCells);
+
+        pmesh.pointZones().clear();
+        pmesh.faceZones().clear();
+        pmesh.cellZones().clear();
+        pmesh.addZones(List<pointZone*>(), List<faceZone*>(), cz);
+    }
+
+    return meshPtr;
+}
+
+
+// ************************************************************************* //

src/mesh/blockMesh/PDRblockMesh/PDRblockCreate.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -300,7 +300,8 @@ Foam::label Foam::PDRblock::addBoundaryFaces
 }
 
 
-Foam::autoPtr<Foam::polyMesh> Foam::PDRblock::mesh(const IOobject& io) const
+Foam::autoPtr<Foam::polyMesh>
+Foam::PDRblock::innerMesh(const IOobject& io) const
 {
     pointField pts(nPoints());
 
@@ -376,4 +377,20 @@ Foam::autoPtr<Foam::polyMesh> Foam::PDRblock::mesh(const IOobject& io) const
 }
 
 
+Foam::autoPtr<Foam::polyMesh>
+Foam::PDRblock::mesh(const IOobject& io) const
+{
+    if (outer_.active())
+    {
+        Info<< "Outer region is active, using blockMesh generation" << nl;
+        return meshBlockMesh(io);
+    }
+    else
+    {
+        Info<< "Outer region is inactive, using ijk generation" << nl;
+        return innerMesh(io);
+    }
+}
+
+
 // ************************************************************************* //

src/mesh/blockMesh/PDRblockMesh/PDRblockI.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -25,20 +25,8 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-inline Foam::PDRblock::PDRblock()
-:
-    ijkMesh(),
-    verbose_(false),
-    grid_(),
-    bounds_(),
-    patches_(),
-    minEdgeLen_(Zero)
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+// Location
 
 inline bool Foam::PDRblock::location::valid() const
 {
@@ -82,6 +70,12 @@ inline Foam::scalar Foam::PDRblock::location::centre() const
 }
 
 
+inline Foam::scalar Foam::PDRblock::location::length() const
+{
+    return scalarList::empty() ? 0 : mag(last() - first());
+}
+
+
 inline void Foam::PDRblock::location::checkIndex(const label i) const
 {
     if (i < 0 || i >= nCells())
@@ -148,6 +142,8 @@ Foam::PDRblock::location::clip(const scalar& val) const
 }
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
 inline const Foam::Vector<Foam::PDRblock::location>&
 Foam::PDRblock::grid() const
 {
@@ -155,9 +151,9 @@ Foam::PDRblock::grid() const
 }
 
 
-inline const Foam::scalar& Foam::PDRblock::minEdgeLen() const
+inline const Foam::scalarMinMax& Foam::PDRblock::edgeLimits() const
 {
-    return minEdgeLen_;
+    return edgeLimits_;
 }
 
 

src/mesh/blockMesh/PDRblockMesh/PDRblockLocation.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -26,13 +26,211 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "PDRblock.H"
-#include "ListOps.H"
+#include "gradingDescriptors.H"
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Prepend a value by shifting contents
+template<class T>
+static void prependList(List<T>& list, const T& val)
+{
+    const label oldLen = list.size();
+    list.resize(oldLen + 1);
+
+    for (label i = oldLen; i > 0; --i)
+    {
+        list[i] = std::move(list[i-1]);
+    }
+
+    list[0] = val;
+}
+
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+Foam::boundBox Foam::PDRblock::bounds
+(
+    const scalarList& x,
+    const scalarList& y,
+    const scalarList& z
+)
+{
+    return boundBox
+    (
+        point(x.first(), y.first(), z.first()),
+        point(x.last(),  y.last(),  z.last())
+    );
+}
+
+
+Foam::Vector<Foam::gradingDescriptors>
+Foam::PDRblock::grading(const Vector<gridControl>& ctrl)
+{
+    return Vector<gradingDescriptors>
+    (
+        ctrl.x().grading(),
+        ctrl.y().grading(),
+        ctrl.z().grading()
+    );
+}
+
+Foam::labelVector
+Foam::PDRblock::sizes(const Vector<gridControl>& ctrl)
+{
+    return labelVector
+    (
+        ctrl.x().nCells(),
+        ctrl.y().nCells(),
+        ctrl.z().nCells()
+    );
+}
+
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+Foam::label Foam::PDRblock::gridControl::nCells() const
+{
+    label nTotal = 0;
+    for (const label nDiv : divisions_)
+    {
+        nTotal += nDiv;
+    }
+
+    return nTotal;
+}
+
+
+Foam::gradingDescriptors Foam::PDRblock::gridControl::grading() const
+{
+    // Begin/end nodes for each segment
+    const scalarList& knots = *this;
+
+    gradingDescriptors gds(divisions_.size());
+
+    forAll(gds, i)
+    {
+        //- Construct from components
+        gds[i] = gradingDescriptor
+        (
+            knots[i+1] - knots[i],  // blockFraction from delta
+            divisions_[i],          // nDivFraction  from nDivs
+            expansion_[i]
+        );
+    }
+
+    gds.normalise();
+
+    return gds;
+}
+
+
+void Foam::PDRblock::gridControl::append
+(
+    const scalar p,
+    const label nDiv,
+    scalar expRatio
+)
+{
+    // Begin/end nodes for each segment
+    scalarList& knots = *this;
+
+    // Is monotonic?
+    if (knots.size() && (p <= knots.last()))
+    {
+        WarningInFunction
+            << "Cannot append point " << p
+            << " which is <= last value " << knots.last() << endl;
+        return;
+    }
+
+    if (nDiv < 1)
+    {
+        WarningInFunction
+            << "Negative or zero divisions " << nDiv << endl;
+        return;
+    }
+
+    // Correct expansion ratios - negative is the same as inverse.
+    if (expRatio < 0)
+    {
+        expRatio = 1.0/(-expRatio);
+    }
+    else if (equal(expRatio, 0))
+    {
+        expRatio = 1;
+    }
+
+    // Now append (push_back)
+    knots.append(p);
+    divisions_.append(nDiv);
+    expansion_.append(expRatio);
+}
+
+
+void Foam::PDRblock::gridControl::prepend
+(
+    const scalar p,
+    const label nDiv,
+    scalar expRatio
+)
+{
+    // Begin/end nodes for each segment
+    scalarList& knots = static_cast<scalarList&>(*this);
+
+    // Is monotonic?
+    if (knots.size() && (p >= knots.first()))
+    {
+        WarningInFunction
+            << "Cannot prepend point " << p
+            << " which is >= first value " << knots.first() << endl;
+        return;
+    }
+
+    if (nDiv < 1)
+    {
+        WarningInFunction
+            << "Negative or zero divisions " << nDiv << endl;
+        return;
+    }
+
+    // Correct expansion ratios - negative is the same as inverse.
+    if (expRatio < 0)
+    {
+        expRatio = 1.0/(-expRatio);
+    }
+    else if (equal(expRatio, 0))
+    {
+        expRatio = 1;
+    }
+
+    // Now prepend (push_front)
+    prependList(knots, p);
+    prependList(divisions_, nDiv);
+    prependList(expansion_, expRatio);
+}
+
+
+Foam::scalarMinMax Foam::PDRblock::location::edgeLimits() const
+{
+    scalarMinMax limits;
+
+    for (label edgei = 0; edgei < this->nCells(); ++edgei)
+    {
+        limits.add(width(edgei));
+    }
+
+    return limits;
+}
+
+
 Foam::label Foam::PDRblock::location::findCell(const scalar p) const
 {
-    if (scalarList::empty())
+    if (scalarList::empty() || p < first() || p > last())
     {
         return -1;
     }

src/mesh/blockMesh/PDRblockMesh/PDRblockOuter.C

@@ -0,0 +1,245 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 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 "PDRblock.H"
+#include "dictionary.H"
+#include "Switch.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+const Foam::Enum
+<
+    Foam::PDRblock::outerControl::controlType
+>
+Foam::PDRblock::outerControl::controlNames_
+({
+    { controlType::OUTER_NONE,  "none" },
+    { controlType::OUTER_EXTEND, "extend" },
+    { controlType::OUTER_BOX, "box" },
+    { controlType::OUTER_SPHERE, "sphere" },
+});
+
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Get a single or a pair of values
+template<class T>
+static Vector2D<T> getLazyPair(const word& name, const dictionary& dict)
+{
+    if (token(dict.lookup(name)).isNumber())
+    {
+        return Vector2D<T>::uniform(dict.get<T>(name));
+    }
+
+    return dict.get<Vector2D<T>>(name);
+}
+
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::PDRblock::outerControl::outerControl()
+:
+    type_(controlType::OUTER_NONE),
+    expandType_(expansionType::EXPAND_RATIO),
+    onGround_(false),
+    relSize_(0,0),
+    nCells_(0,0),
+    expansion_(1,1)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::PDRblock::outerControl::clear()
+{
+    type_ = controlType::OUTER_NONE;
+    expandType_ = expansionType::EXPAND_RATIO;
+    onGround_ = false;
+    relSize_ = Zero;
+    nCells_ = Zero;
+    expansion_ = vector2D::uniform(1);
+}
+
+
+void Foam::PDRblock::outerControl::report(Ostream& os) const
+{
+    if (active())
+    {
+        os  << "Has outer region: " << controlNames_[type_] << nl
+            << " onGround : " << Switch::name(onGround_) << nl
+            << "    sizes : " << relSize_ << nl
+            << "   nCells : " << nCells_ << nl;
+    }
+    else
+    {
+        os  << "No outer region" << nl;
+    }
+}
+
+
+bool Foam::PDRblock::outerControl::active() const
+{
+    return (controlType::OUTER_NONE != type_);
+}
+
+
+bool Foam::PDRblock::outerControl::isSphere() const
+{
+    return (controlType::OUTER_SPHERE == type_);
+}
+
+
+bool Foam::PDRblock::outerControl::onGround() const
+{
+    return onGround_;
+}
+
+
+bool Foam::PDRblock::outerControl::onGround(const bool on)
+{
+    bool old(onGround_);
+    onGround_ = on;
+    return old;
+}
+
+
+void Foam::PDRblock::outerControl::read(const dictionary& dict)
+{
+    clear();
+
+    type_ = controlNames_.getOrDefault("type", dict, controlType::OUTER_NONE);
+    onGround_ = dict.getOrDefault("onGround", false);
+
+    if (controlType::OUTER_NONE == type_)
+    {
+        return;
+    }
+
+    // Everything else
+
+    nCells_ = getLazyPair<label>("nCells", dict);
+    relSize_ = getLazyPair<scalar>("size", dict);
+
+    expandType_ =
+        expansionNames_.getOrDefault
+        (
+            "expansion",
+            dict,
+            expansionType::EXPAND_RATIO
+        );
+
+
+    if (dict.found("ratios"))
+    {
+        expansion_ = getLazyPair<scalar>("ratios", dict);
+    }
+    else
+    {
+        if (expandType_ != expansionType::EXPAND_UNIFORM)
+        {
+            expandType_ = expansionType::EXPAND_UNIFORM;
+            // Info << "Warning: no 'ratios', use uniform spacing" << nl;
+        }
+    }
+
+    if (expandType_ == expansionType::EXPAND_UNIFORM)
+    {
+        expansion_ = vector2D::uniform(1);
+    }
+
+
+    // Errors
+    int die = 0;
+
+    if (nCells_.x() <= 1 || nCells_.y() <= 1)
+    {
+        if (!die++)
+        {
+            FatalIOErrorInFunction(dict);
+        }
+        FatalIOError
+            << "Too few outer cells: " << nCells_ << nl;
+    }
+
+    if (relSize_.x() <= 1 || relSize_.y() <= 1)
+    {
+        if (!die++)
+        {
+            FatalIOErrorInFunction(dict);
+        }
+        FatalIOError
+            << "Outer dimensions must be > 1. Had " << relSize_ << nl;
+    }
+
+    if (die)
+    {
+        FatalIOError << nl << exit(FatalIOError);
+    }
+
+
+    // Warnings
+
+    if
+    (
+        controlType::OUTER_BOX == type_
+     || controlType::OUTER_SPHERE == type_
+    )
+    {
+        if (relSize_.x() < 2 || relSize_.y() < 2)
+        {
+            WarningInFunction
+                << "Outer dimensions "
+                << relSize_ << " too small for "
+                << controlNames_[type_] << " - switching to "
+                << controlNames_[controlType::OUTER_EXTEND] << nl;
+
+            type_ = controlType::OUTER_EXTEND;
+        }
+    }
+
+    if (controlType::OUTER_SPHERE == type_)
+    {
+        if (relSize_.x() < 3 || relSize_.y() < 3)
+        {
+            WarningInFunction
+                << "Outer dimensions "
+                << relSize_ << " too small for "
+                << controlNames_[type_] << " - switching to "
+                << controlNames_[controlType::OUTER_BOX] << nl;
+
+            type_ = controlType::OUTER_EXTEND;
+        }
+    }
+}
+
+
+// ************************************************************************* //

tutorials/mesh/PDRblockMesh/box0/system/PDRblockMeshDict

@@ -75,4 +75,23 @@ boundary
     }
 );
 
+
+// Treatment of the outer region
+
+outer
+{
+    type        sphere;     // (none | extend | box | sphere)  [default: none]
+    onGround    true;       // Module on the ground?  [default: false]
+    expansion   relative;   // (uniform | ratio | relative)  [default: ratio]
+
+    ratios      1.1;
+
+    size        3;          // Overall outer/inner size
+    nCells      10;         // Number of cells for outer region
+
+    // size       (3 4);
+    // nCells     (10 12);
+}
+
+
 // ************************************************************************* //

tutorials/preProcessing/PDRsetFields/simplePipeCage/system/PDRblockMeshDict

@@ -74,4 +74,23 @@ boundary
     }
 );
 
+
+// Treatment of the outer region
+
+outer
+{
+    type        sphere;     // (none | extend | box | sphere)  [default: none]
+    onGround    true;       // Module on the ground?  [default: false]
+    expansion   relative;   // (uniform | ratio | relative)  [default: ratio]
+
+    ratios      1.1;
+
+    size        3;          // Overall outer/inner size
+    nCells      10;         // Number of cells for outer region
+
+    // size       (3 4);
+    // nCells     (10 12);
+}
+
+
 //***************************************************************************//