Merge branch 'develop' of develop.openfoam.com:Development/OpenFOAM-plus into develop

applications/solvers/combustion/XiFoam/XiEngineFoam/logSummary.H

@@ -1,15 +1,22 @@
-Info<< "Mean pressure:" << p.weightedAverage(mesh.V()).value() << endl;
+{
-Info<< "Mean temperature:" << thermo.T().weightedAverage(mesh.V()).value()
+    const scalar meanP = p.weightedAverage(mesh.V()).value();
-    << endl;
+    const scalar meanT = thermo.T().weightedAverage(mesh.V()).value();
-Info<< "Mean u':"
+    const scalar meanUp =
-    << (sqrt((2.0/3.0)*turbulence->k()))().weightedAverage(mesh.V()).value()
+        (sqrt((2.0/3.0)*turbulence->k()))().weightedAverage(mesh.V()).value();
-    << endl;
+    const scalar meanB = b.weightedAverage(mesh.V()).value();
 
-logSummaryFile()
+    Info<< "Mean pressure:" << meanP << nl
+        << "Mean temperature:" << meanT << nl
+        << "Mean u':" << meanUp << endl;
+
+    if (Pstream::master())
+    {
+        logSummaryFile()
             << runTime.theta() << tab
-    << p.weightedAverage(mesh.V()).value() << tab
+            << meanP << tab
-    << thermo.T().weightedAverage(mesh.V()).value() << tab
+            << meanT << tab
-    << (sqrt((2.0/3.0)*turbulence->k()))().weightedAverage(mesh.V()).value()
+            << meanUp << tab
-    << tab
+            << 1 - meanB
-    << 1 - b.weightedAverage(mesh.V()).value()
             << endl;
+    }
+}

etc/caseDicts/annotated/decomposeParDict

@@ -53,11 +53,10 @@ regions
 }
 
 
-// Coefficients for the decomposition method are either as a
+// Coefficients for the decomposition method are either as a general
-// general "coeffs" dictionary or method-specific "<method>Coeffs".
+// "coeffs" dictionary or method-specific "<method>Coeffs".
 // For multiLevel, using multiLevelCoeffs only.
 
-
 multiLevelCoeffs
 {
     // multiLevel decomposition methods to apply in turn.
@@ -107,6 +106,11 @@ multiLevelCoeffs
 
 // Other example coefficients
 
+coeffs
+{
+    n       (2 1 1);
+}
+
 simpleCoeffs
 {
     n       (2 1 1);

src/OpenFOAM/global/fileOperations/fileOperation/fileOperation.C

@@ -104,51 +104,51 @@ Foam::instantList Foam::fileOperation::sortTimes
 )
 {
     // Initialise instant list
-    instantList Times(dirEntries.size() + 1);
+    instantList times(dirEntries.size() + 1);
     label nTimes = 0;
 
     // Check for "constant"
     bool haveConstant = false;
-    forAll(dirEntries, i)
+    for (const fileName& dirName : dirEntries)
     {
-        if (dirEntries[i] == constantName)
+        if (dirName == constantName)
         {
-            Times[nTimes].value() = 0;
-            Times[nTimes].name() = dirEntries[i];
-            nTimes++;
             haveConstant = true;
+            times[nTimes].value() = 0;
+            times[nTimes].name() = constantName;
+            ++nTimes;
             break;
         }
     }
 
     // Read and parse all the entries in the directory
-    forAll(dirEntries, i)
+    for (const fileName& dirName : dirEntries)
     {
         scalar timeValue;
-        if (readScalar(dirEntries[i], timeValue))
+        if (readScalar(dirName, timeValue))
         {
-            Times[nTimes].value() = timeValue;
+            times[nTimes].value() = timeValue;
-            Times[nTimes].name() = dirEntries[i];
+            times[nTimes].name() = dirName;
-            nTimes++;
+            ++nTimes;
         }
     }
 
     // Reset the length of the times list
-    Times.setSize(nTimes);
+    times.setSize(nTimes);
 
     if (haveConstant)
     {
         if (nTimes > 2)
         {
-            std::sort(&Times[1], Times.end(), instant::less());
+            std::sort(&times[1], times.end(), instant::less());
         }
     }
     else if (nTimes > 1)
     {
-        std::sort(&Times[0], Times.end(), instant::less());
+        std::sort(&times[0], times.end(), instant::less());
     }
 
-    return Times;
+    return times;
 }
 
 
@@ -161,15 +161,15 @@ void Foam::fileOperation::mergeTimes
 {
     if (extraTimes.size())
     {
-        bool haveConstant =
+        const bool haveConstant =
         (
-            times.size() > 0
+            times.size()
          && times[0].name() == constantName
         );
 
-        bool haveExtraConstant =
+        const bool haveExtraConstant =
         (
-            extraTimes.size() > 0
+            extraTimes.size()
          && extraTimes[0].name() == constantName
         );
 
@@ -228,9 +228,7 @@ void Foam::fileOperation::mergeTimes
 
 bool Foam::fileOperation::isFileOrDir(const bool isFile, const fileName& f)
 {
-    return
+    return (isFile ? Foam::isFile(f) : Foam::isDir(f));
-        (isFile && Foam::isFile(f))
-     || (!isFile && Foam::isDir(f));
 }
 
 
@@ -431,10 +429,8 @@ Foam::autoPtr<Foam::fileOperation> Foam::fileOperation::New
     bool verbose
 )
 {
-    if (debug)
+    DebugInFunction
-    {
+        << "Constructing fileHandler" << endl;
-        InfoInFunction << "Constructing fileHandler" << endl;
-    }
 
     auto cstrIter = wordConstructorTablePtr_->cfind(handlerType);
 
@@ -570,14 +566,12 @@ Foam::fileName Foam::fileOperation::filePath(const fileName& fName) const
         }
         return fName;
     }
-    else
+
-    {
     if (debug)
     {
         Pout<< "fileOperation::filePath : Not found" << endl;
     }
     return fileName::null;
-    }
 }
 
 
@@ -755,13 +749,10 @@ Foam::IOobject Foam::fileOperation::findInstance
 
     if (exists(io))
     {
-        if (debug)
+        DebugInFunction
-        {
-            InfoInFunction
             << "Found exact match for \"" << io.name()
             << "\" in " << io.instance()/io.local()
             << endl;
-        }
 
         return io;
     }
@@ -770,9 +761,9 @@ Foam::IOobject Foam::fileOperation::findInstance
     // closest to and lower than current time
 
     instantList ts = time.times();
-    label instanceI;
+    label instanceI = ts.size()-1;
 
-    for (instanceI = ts.size()-1; instanceI >= 0; --instanceI)
+    for (; instanceI >= 0; --instanceI)
     {
         if (ts[instanceI].value() <= startValue)
         {
@@ -780,7 +771,7 @@ Foam::IOobject Foam::fileOperation::findInstance
         }
     }
 
-    // continue searching from here
+    // Continue searching from here
     for (; instanceI >= 0; --instanceI)
     {
         // Shortcut: if actual directory is the timeName we've already tested it
@@ -796,13 +787,10 @@ Foam::IOobject Foam::fileOperation::findInstance
         io.instance() = ts[instanceI].name();
         if (exists(io))
         {
-            if (debug)
+            DebugInFunction
-            {
-                InfoInFunction
                 << "Found exact match for \"" << io.name()
                 << "\" in " << io.instance()/io.local()
                 << endl;
-            }
 
             return io;
         }
@@ -810,11 +798,8 @@ Foam::IOobject Foam::fileOperation::findInstance
         // Check if hit minimum instance
         if (ts[instanceI].name() == stopInstance)
         {
-            if (debug)
+            DebugInFunction
-            {
-                InfoInFunction
                 << "Hit stopInstance " << stopInstance << endl;
-            }
 
             if
             (
@@ -845,27 +830,30 @@ Foam::IOobject Foam::fileOperation::findInstance
         }
     }
 
-    // times() usually already includes the constant() so would have been
+    // Times usually already includes 'constant' so would have been checked
-    // checked above. Re-test if
+    // above.
-    // - times() is empty. Sometimes this can happen (e.g. decomposePar with
+    // However, re-test under these conditions:
-    //   collated)
+    // - Times is empty.
-    // - times()[0] is not constant
+    //   Sometimes this can happen (eg, decomposePar with collated)
-    if (!ts.size() || ts[0].name() != time.constant())
+    // - Times[0] is not constant
-    {
+    // - The startValue is negative (eg, kivaTest).
-        // Note. This needs to be a hard-coded constant, rather than the
+    //   This plays havoc with the reverse search, causing it to miss 'constant'
-        // constant function of the time, because the latter points to
-        // the case constant directory in parallel cases
 
+    if
+    (
+        ts.empty()
+     || ts.first().name() != time.constant()
+     || startValue < 0
+    )
+    {
         io.instance() = time.constant();
         if (exists(io))
         {
-            if (debug)
+            DebugInFunction
-            {
-                InfoInFunction
                 << "Found constant match for \"" << io.name()
                 << "\" in " << io.instance()/io.local()
                 << endl;
-            }
+
             return io;
         }
     }

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/OpenFOAM/primitives/functions/Function1/Table/TableBase.H

@@ -62,7 +62,7 @@ class TableBase
 {
 protected:
 
-    // Protected data
+    // Protected Data
 
         //- Table name
         const word name_;
@@ -93,9 +93,6 @@ protected:
         //- Return (demand driven) interpolator
         const interpolationWeights& interpolator() const;
 
-
-private:
-
         //- No copy assignment
         void operator=(const TableBase<Type>&) = delete;
 

src/OpenFOAM/primitives/functions/Function1/TableFile/TableFile.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016-2017 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2016-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2017 OpenFOAM Foundation
@@ -37,12 +37,12 @@ Foam::Function1Types::TableFile<Type>::TableFile
 )
 :
     TableBase<Type>(entryName, dict),
-    fName_("none")
+    fName_()
 {
     dict.readEntry("file", fName_);
 
     fileName expandedFile(fName_);
-    //IFstream is(expandedFile.expand());
+
     autoPtr<ISstream> isPtr(fileHandler().NewIFstream(expandedFile.expand()));
     ISstream& is = isPtr();
 
@@ -66,13 +66,6 @@ Foam::Function1Types::TableFile<Type>::TableFile(const TableFile<Type>& tbl)
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class Type>
-Foam::Function1Types::TableFile<Type>::~TableFile()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>

src/OpenFOAM/primitives/functions/Function1/TableFile/TableFile.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2017 OpenFOAM Foundation
@@ -27,7 +27,7 @@ Class
     Foam::Function1Types::TableFile
 
 Description
-    Templated table container function where data is read from file.
+    Templated table container function where data are read from file.
 
     Usage:
     \verbatim
@@ -40,8 +40,8 @@ Description
         }
     \endverbatim
 
-    Data is stored as a list of Tuple2's. First column is always stored as
+    Data are stored as a list of Tuple2's. First column is always stored as
-    scalar entries.  Data is read in the form, e.g. for an entry \<entryName\>
+    scalar entries.  Data are read in the form, e.g. for an entry \<entryName\>
     that is (scalar, vector):
     \verbatim
         (
@@ -77,7 +77,7 @@ class TableFile
 :
     public TableBase<Type>
 {
-    // Private data
+    // Private Data
 
         //- File name for csv table (optional)
         fileName fName_;
@@ -97,7 +97,7 @@ public:
 
     // Constructors
 
-        //- Construct from entry name and Istream
+        //- Construct from entry name and "file" found in dictionary
         TableFile(const word& entryName, const dictionary& dict);
 
         //- Copy constructor
@@ -105,7 +105,7 @@ public:
 
 
     //- Destructor
-    virtual ~TableFile();
+    virtual ~TableFile() = default;
 
 
     // I/O

src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H

@@ -87,7 +87,7 @@ Usage
             componentColumns 1(1);
             separator       ",";
             mergeSeparators no;
-            file            "<constant>/pressureVsU";
+            file            "<constant>/UvsPressure";
         }
         value           uniform 0;
     }

src/finiteVolume/fields/fvPatchFields/derived/swirlFanVelocity/swirlFanVelocityFvPatchField.C

@@ -107,7 +107,6 @@ void Foam::swirlFanVelocityFvPatchField::calcFanJump()
 }
 
 
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::swirlFanVelocityFvPatchField::swirlFanVelocityFvPatchField
@@ -156,7 +155,6 @@ Foam::swirlFanVelocityFvPatchField::swirlFanVelocityFvPatchField
         this->cyclicPatch().owner()
       ? Function1<scalar>::New("rpm", dict)
       : nullptr
-
     ),
     rEff_(dict.lookupOrDefault<scalar>("rEff", 0)),
     fanEff_(dict.lookupOrDefault<scalar>("fanEff", 1)),
@@ -248,7 +246,11 @@ void Foam::swirlFanVelocityFvPatchField::write(Ostream& os) const
         os.writeEntryIfDifferent<word>("p", "p", pName_);
         os.writeEntryIfDifferent<word>("rho", "rho", rhoName_);
         os.writeEntry("origin", origin_);
+
+        if (rpm_)
+        {
             rpm_->writeData(os);
+        }
 
         os.writeEntryIfDifferent<scalar>("rEff", 0.0, rEff_);
         os.writeEntryIfDifferent<bool>("useRealRadius", false, useRealRadius_);

src/finiteVolume/fields/fvPatchFields/derived/uniformJump/uniformJumpFvPatchField.C

@@ -131,6 +131,7 @@ template<class Type>
 void Foam::uniformJumpFvPatchField<Type>::write(Ostream& os) const
 {
     fixedJumpFvPatchField<Type>::write(os);
+
     if (this->cyclicPatch().owner())
     {
         jumpTable_->writeData(os);

src/finiteVolume/fields/fvPatchFields/derived/uniformJump/uniformJumpFvPatchField.H

@@ -94,7 +94,7 @@ protected:
 
     // Protected data
 
-        //- "jump" table
+        //- The "jump" table
         autoPtr<Function1<Type>> jumpTable_;
 
 

src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H

@@ -211,7 +211,7 @@ class surfaceFieldValue
 {
 public:
 
-    // Public data types
+    // Public Data Types
 
         //- Region type enumeration
         enum regionTypes

src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueTemplates.C

@@ -98,25 +98,8 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::getFieldValues
             if (sampledPtr_().interpolate())
             {
                 const interpolationCellPoint<Type> interp(fld);
-                tmp<Field<Type>> tintFld(sampledPtr_().interpolate(interp));
-                const Field<Type>& intFld = tintFld();
 
-                // Average
+                return sampledPtr_().sample(interp);
-                const faceList& faces = sampledPtr_().faces();
-                auto tavg = tmp<Field<Type>>::New(faces.size(), Zero);
-                auto& avg = tavg.ref();
-
-                forAll(faces, facei)
-                {
-                    const face& f = faces[facei];
-                    for (const label labi : f)
-                    {
-                        avg[facei] += intFld[labi];
-                    }
-                    avg[facei] /= f.size();
-                }
-
-                return tavg;
             }
             else
             {

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)
-        {
-            // 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

@@ -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;
+}
+
+
+// ************************************************************************* //

tutorials/combustion/XiEngineFoam/kivaTest/Allrun-parallel

@@ -0,0 +1,13 @@
+#!/bin/sh
+cd ${0%/*} || exit 1                        # Run from this directory
+. $WM_PROJECT_DIR/bin/tools/RunFunctions    # Tutorial run functions
+
+runApplication kivaToFoam -file otape17
+
+runApplication decomposePar
+
+runParallel $(getApplication)
+
+runApplication reconstructPar
+
+#------------------------------------------------------------------------------

tutorials/combustion/XiEngineFoam/kivaTest/system/decomposeParDict

@@ -0,0 +1,27 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1812                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains  4;
+
+method  hierarchical;
+
+coeffs
+{
+    n       (2 2 1);
+}
+
+
+// ************************************************************************* //