From 920ff7bcc466a5c15357098b5a7f437864463d6c Mon Sep 17 00:00:00 2001
From: Vuko Vukcevic <vuko.vukcevic@fsb.hr>
Date: Wed, 23 Jan 2019 13:03:44 +0100
Subject: [PATCH] Manually merged Henrik's DLB and some other stuff

---
 ...htRcThermalDiffusivityFvPatchScalarField.C |   2 +-
 .../finiteVolume/domainDecomposition.C        |  20 ++
 .../finiteVolume/domainDecomposition.H        |   7 +
 .../fvFieldReconstructorReconstructFields.C   |   8 +-
 .../finiteVolume/processorMeshesRebuild.C     |  28 +-
 .../processorMeshesReconstructor.C            |   2 +
 .../processorMeshesReconstructor.H            |   9 +
 .../finiteVolume/sharedPoints.C               | 171 +++++++++-
 .../finiteVolume/sharedPoints.H               |  15 +-
 .../dynamicPolyRefinementFvMesh.C             |  14 +-
 .../topoChangerFvMeshLoadBalance.C            |  58 +++-
 .../topoChangerFvMeshLoadBalanceTemplates.C   |  24 +-
 .../compressible/compressibleCourantNo.H      |  21 +-
 .../cfdTools/incompressible/CourantNo.H       |   2 +-
 .../fixedGradient/fixedGradientFvPatchField.C |   8 +-
 .../basic/mixed/mixedFvPatchField.C           |   8 +-
 .../inletOutlet/inletOutletFvPatchField.C     |   2 +-
 .../finiteVolume/fvSchemes/fvSchemes.C        |  29 +-
 src/finiteVolume/fvMesh/fvMesh.C              |   3 +
 .../fvMesh/wallDist/nearWallDist.C            |   2 +-
 src/foam/Make/files                           |   1 +
 .../GeometricField/GeometricField.C           |  10 +
 .../GeometricField/GeometricField.H           |   4 +
 src/foam/fields/cloud/cloud.C                 |   7 +
 src/foam/fields/cloud/cloud.H                 |   4 +-
 .../pointBoundaryMesh/pointBoundaryMesh.C     |  32 +-
 .../pointBoundaryMesh/pointBoundaryMesh.H     |   6 +-
 src/foam/meshes/pointMesh/pointMesh.C         |   5 +-
 .../globalMeshData/globalProcFaceIndex.C      |   6 +-
 .../globalMeshData/globalProcFaceIndex.H      |   2 +-
 .../globalMeshData/globalProcPointIndex.C     | 319 ++++++++++++++++++
 .../globalMeshData/globalProcPointIndex.H     | 135 ++++++++
 ...undaryDynamicRefineSolidBodyMotionFvMesh.C |   2 +-
 ...aryEpsilonWallFunctionFvPatchScalarField.C |   2 +-
 ...ersedBoundaryKqRWallFunctionFvPatchField.C |   2 +-
 ...ndaryOmegaWallFunctionFvPatchScalarField.C |   2 +-
 .../CloudDistributeTemplate.H                 |   2 +-
 .../basic/CloudTemplate/CloudTemplate.C       |  34 ++
 .../basic/CloudTemplate/CloudTemplate.H       |   3 +
 .../InjectionModel/InjectionModel.C           |   2 +-
 .../oversetMesh/donorAcceptor/donorAcceptor.H |   3 +-
 .../adaptiveOverlapFringe.C                   | 125 ++++---
 .../adaptiveOverlapFringe.H                   |  10 +-
 .../donorSuitability/donorSuitability.H       |  15 +
 .../overlapFringe/overlapFringe.C             |  63 ++--
 .../overlapFringe/overlapFringe.H             |   2 +-
 .../oversetMesh/oversetMeshAddressing.C       |  21 +-
 .../oversetMesh/oversetRegion/oversetRegion.C |  92 ++---
 src/sampling/sampledSet/curve/curveSet.C      |   4 +-
 src/sampling/sampledSet/face/faceOnlySet.C    |  16 +-
 src/sampling/sampledSet/face/faceOnlySet.H    |   7 +
 .../sampledSet/sampledSet/sampledSet.C        |  24 +-
 .../sampledSet/sampledSet/sampledSet.H        |   8 +-
 src/sampling/sampledSet/uniform/uniformSet.C  |   4 +-
 .../basicChemistryModel/basicChemistryModel.C |  17 +-
 .../basicChemistryModel/basicChemistryModel.H |   6 +-
 .../incompressible/RAS/RASModel/RASModel.H    |   2 +-
 .../turbulenceModel/turbulenceModel.H         |   2 +-
 58 files changed, 1225 insertions(+), 209 deletions(-)
 create mode 100644 src/foam/meshes/polyMesh/globalMeshData/globalProcPointIndex.C
 create mode 100644 src/foam/meshes/polyMesh/globalMeshData/globalProcPointIndex.H

diff --git a/src/conjugateHeatTransfer/fvPatchFields/chtRcThermalDiffusivity/chtRcThermalDiffusivityFvPatchScalarField.C b/src/conjugateHeatTransfer/fvPatchFields/chtRcThermalDiffusivity/chtRcThermalDiffusivityFvPatchScalarField.C
index 54497a652..9f15e1af4 100644
--- a/src/conjugateHeatTransfer/fvPatchFields/chtRcThermalDiffusivity/chtRcThermalDiffusivityFvPatchScalarField.C
+++ b/src/conjugateHeatTransfer/fvPatchFields/chtRcThermalDiffusivity/chtRcThermalDiffusivityFvPatchScalarField.C
@@ -261,7 +261,7 @@ Foam::chtRcThermalDiffusivityFvPatchScalarField::calcThermalDiffusivity
 
     const scalarField kOwn = fOwn/(1.0 - p.weights())/mld.magDelta(p.index());
     const scalarField kNei = fNei/p.weights()/mld.magDelta(p.index());
- 
+
     tmp<scalarField> kTmp(new scalarField(p.size()));
     scalarField& k = kTmp();
 
diff --git a/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/domainDecomposition.C b/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/domainDecomposition.C
index 3d9069908..372ec808d 100644
--- a/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/domainDecomposition.C
+++ b/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/domainDecomposition.C
@@ -34,6 +34,7 @@ License
 #include "OSspecific.H"
 #include "Map.H"
 #include "DynamicList.H"
+#include "globalMeshData.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -53,6 +54,7 @@ Foam::domainDecomposition::domainDecomposition
     nProcs_(readInt(decompositionDict_.lookup("numberOfSubdomains"))),
     distributed_(false),
     gfIndex_(mesh_),
+    gpIndex_(mesh_),
     cellToProc_(mesh_.nCells()),
     patchNbrCellToProc_(mesh_.boundaryMesh().size()),
     procPointAddressing_(nProcs_),
@@ -503,6 +505,24 @@ Foam::autoPtr<Foam::fvMesh> Foam::domainDecomposition::processorMesh
 }
 
 
+Foam::labelList Foam::domainDecomposition::globalPointIndex
+(
+    const label procI
+) const
+{
+    const labelList& gppi = gpIndex_.globalLabel();
+    const labelList& ppAddr = procPointAddressing_[procI];
+    labelList globalPointIndex(ppAddr.size());
+
+    forAll (globalPointIndex, gpI)
+    {
+        globalPointIndex[gpI] = gppi[ppAddr[gpI]];
+    }
+
+    return globalPointIndex;
+}
+
+
 bool Foam::domainDecomposition::writeDecomposition()
 {
     Info<< "\nConstructing processor meshes" << endl;
diff --git a/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/domainDecomposition.H b/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/domainDecomposition.H
index b4f02aa41..fc32d24e4 100644
--- a/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/domainDecomposition.H
+++ b/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/domainDecomposition.H
@@ -43,6 +43,7 @@ SourceFiles
 #include "PtrList.H"
 #include "point.H"
 #include "globalProcFaceIndex.H"
+#include "globalProcPointIndex.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -72,6 +73,9 @@ class domainDecomposition
         //- Global face index
         globalProcFaceIndex gfIndex_;
 
+        //- Global face index
+        globalProcPointIndex gpIndex_;
+
         //- Processor label for each cell
         labelList cellToProc_;
 
@@ -205,6 +209,9 @@ public:
                 const bool createPassiveProcPatches = false
             ) const;
 
+            //- Create and return global point index
+            labelList globalPointIndex(const label procI) const;
+
             //- Return processor point addressing
             const labelListList& procPointAddressing() const
             {
diff --git a/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/fvFieldReconstructorReconstructFields.C b/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/fvFieldReconstructorReconstructFields.C
index 47567fa5a..dcf8dc227 100644
--- a/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/fvFieldReconstructorReconstructFields.C
+++ b/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/fvFieldReconstructorReconstructFields.C
@@ -44,10 +44,10 @@ void Foam::fvFieldReconstructor::reconstructField
 ) const
 {
     // Get references to internal and boundary field
-    Field<Type>& iField = reconField.internalField();
+    Field<Type>& iField = reconField;
 
     typename GeometricField<Type, fvPatchField, volMesh>::
-        GeometricBoundaryField& bouField = reconField.boundaryField();
+        GeometricBoundaryField& bouField = reconField.boundaryFieldNoStoreOldTimes();
 
     forAll (procFields, procI)
     {
@@ -143,10 +143,10 @@ void Foam::fvFieldReconstructor::reconstructField
 ) const
 {
     // Create the internalField
-    Field<Type>& iField = reconField.internalField();
+    Field<Type>& iField = reconField;
 
     typename GeometricField<Type, fvsPatchField, surfaceMesh>::
-        GeometricBoundaryField& bouField = reconField.boundaryField();
+        GeometricBoundaryField& bouField = reconField.boundaryFieldNoStoreOldTimes();
 
     forAll (procFields, procI)
     {
diff --git a/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesRebuild.C b/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesRebuild.C
index fbe520e7b..9d9655059 100644
--- a/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesRebuild.C
+++ b/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesRebuild.C
@@ -652,10 +652,30 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
 
     reconPatchSizes = 0;
 
-    // Prepare handling for globally shared points.  This is equivalent
-    // to parallel processor points, but working on a PtrList of meshes
-    // on the same processor
-    sharedPoints sharedData(meshes_);
+    // HR 21.12.18 : A bit of a hack for the moment in order to support the new
+    // method (using global point addressing) and old method (syncing across
+    // processor BCs) of constructing the shared points. The old method is still
+    // needed since global point addressing does not exist when the sharedPoint
+    // object is constructed from reconstructParMesh.
+    //
+    // TODO: Unify the methods by constructing global point addressing from the
+    // mesh pieces when. Or even better, calculate procPointAddressing directly
+    // which will simplify the mesh merging immensely.
+    autoPtr<sharedPoints> sharedDataPtr;
+    if(globalPointIndex_.size())
+    {
+        // Determine globally shared points using global point
+        // adressing
+        sharedDataPtr.set(new sharedPoints(meshes_, globalPointIndex_));
+    }
+    else
+    {
+        // Prepare handling for globally shared points.  This is equivalent
+        // to parallel processor points, but working on a PtrList of meshes
+        // on the same processor
+        sharedDataPtr.set(new sharedPoints(meshes_));
+    }
+    sharedPoints& sharedData = sharedDataPtr();
 
     // Record global point index for shared points
     labelList globalPointMapping(sharedData.nGlobalPoints(), -1);
diff --git a/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesReconstructor.C b/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesReconstructor.C
index 18ace98bb..55a754bca 100644
--- a/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesReconstructor.C
+++ b/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesReconstructor.C
@@ -71,6 +71,7 @@ Foam::processorMeshesReconstructor::processorMeshesReconstructor
 :
     meshName_(meshName),
     meshes_(),
+    globalPointIndex_(),
     pointProcAddressing_(),
     faceProcAddressing_(),
     cellProcAddressing_(),
@@ -86,6 +87,7 @@ Foam::processorMeshesReconstructor::processorMeshesReconstructor
 :
     meshName_(meshName),
     meshes_(databases.size()),
+    globalPointIndex_(),
     pointProcAddressing_(),
     faceProcAddressing_(),
     cellProcAddressing_(),
diff --git a/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesReconstructor.H b/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesReconstructor.H
index 17e39cacf..a8e4af0c0 100644
--- a/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesReconstructor.H
+++ b/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesReconstructor.H
@@ -72,6 +72,9 @@ class processorMeshesReconstructor
         //- List of processor meshes
         PtrList<fvMesh> meshes_;
 
+        //- global point index per sub-mesh
+        labelListList globalPointIndex_;
+
         //- List of processor point addressing lists
         PtrList<labelIOList> pointProcAddressing_;
 
@@ -161,6 +164,12 @@ public:
             return meshes_;
         }
 
+        //- Return access to meshes
+        labelListList& globalPointIndex()
+        {
+            return globalPointIndex_;
+        }
+
         //- Return point-processor addressing arrays
         const PtrList<labelIOList>& pointProcAddressing() const;
 
diff --git a/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/sharedPoints.C b/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/sharedPoints.C
index bc8a97d40..dc8b1d220 100644
--- a/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/sharedPoints.C
+++ b/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/sharedPoints.C
@@ -451,9 +451,160 @@ void Foam::sharedPoints::calcSharedPoints()
 }
 
 
+void Foam::sharedPoints::calcSharedPoints
+(
+    const labelListList& globalPointIndex
+)
+{
+    typedef HashTable<Pair<label>, label, Hash<label> > markTable;
+    markTable marker;
+
+    List<labelHashSet> procSharedPoints(meshes_.size());
+
+    // Mark up points for the first time
+    forAll (meshes_, meshI)
+    {
+        if (meshes_.set(meshI))
+        {
+            const labelList& curGlobalPointIndex = globalPointIndex[meshI];
+
+            labelHashSet& sharedPoints = procSharedPoints[meshI];
+
+            forAll (curGlobalPointIndex, pI)
+            {
+                const label gpi = curGlobalPointIndex[pI];
+
+                markTable::iterator iter = marker.find(gpi);
+                if (iter == Map<label>::end())
+                {
+                    // Record meshI and pI
+                    marker.insert(gpi, Pair<label>(-meshI - 1, pI));
+                }
+                else
+                {
+                    if (iter().first() < 0)
+                    {
+                        if(iter().first() != -meshI - 1)
+                        {
+                            // Shared point detected. Add for both meshes
+
+                            // Add for first mesh
+                            const label firstMesh = -iter().first() - 1;
+                            const label firstPointI = iter().second();
+                            procSharedPoints[firstMesh].insert(firstPointI);
+
+                            // Add for current mesh
+                            sharedPoints.insert(pI);
+
+                            // Count shared points and update bookkeeping
+                            iter().first() = nGlobalPoints_;
+                            nGlobalPoints_++;
+
+                        }
+                    }
+                    else
+                    {
+                        // Existing shared point. Add for current mesh
+                        sharedPoints.insert(pI);
+                    }
+                }
+            }
+        }
+    }
+
+    forAll (meshes_, meshI)
+    {
+        if (meshes_.set(meshI))
+        {
+            const labelList& curGlobalPointIndex = globalPointIndex[meshI];
+
+            labelList& curSharedPoints = sharedPointLabels_[meshI];
+            curSharedPoints = procSharedPoints[meshI].toc();
+
+            labelList& curSharedAddr = sharedPointAddr_[meshI];
+            curSharedAddr.setSize(curSharedPoints.size());
+
+            forAll(curSharedPoints, i)
+            {
+                const label gpi = curSharedPoints[i];
+                curSharedAddr[i] = marker.find(curGlobalPointIndex[gpi])().first();
+            }
+        }
+    }
+
+    // debug
+    {
+        pointField gp(nGlobalPoints_);
+        boolList gpSet(nGlobalPoints_, false);
+
+        forAll (meshes_, meshI)
+        {
+            if (meshes_.set(meshI))
+            {
+                // Get points
+                const pointField& P = meshes_[meshI].points();
+
+                // Get list of local point labels that are globally shared
+                const labelList& curShared = sharedPointLabels_[meshI];
+
+                // Get index in global point list
+                const labelList& curAddr = sharedPointAddr_[meshI];
+
+                // Loop through all local points
+                forAll (curShared, i)
+                {
+                    if (!gpSet[curAddr[i]])
+                    {
+                        // Point not set: set it
+                        gp[curAddr[i]] = P[curShared[i]];
+                        gpSet[curAddr[i]] = true;
+                    }
+                    else
+                    {
+                        // Point already set: check location
+                        if (mag(gp[curAddr[i]] - P[curShared[i]]) > SMALL)
+                        {
+                            Info<< "MERGE MISMATCH: mesh" << meshI
+                                << " point: " << curShared[i]
+                                << " dist: " << gp[curAddr[i]] << " "
+                                << P[curShared[i]]
+                                << endl;
+                        }
+                    }
+                }
+            }
+        }
+
+    /*
+            // Grab marked points
+            OFstream ppp("points.vtk");
+            ppp << "# vtk DataFile Version 2.0" << nl
+                << "points.vtk" << nl
+                << "ASCII" << nl
+                << "DATASET POLYDATA" << nl
+                << "POINTS " << nGlobalPoints_ << " float" << nl;
+
+            Pout << "Global marking " << nGlobalPoints_ << endl;
+            forAll (gp, i)
+            {
+                ppp << float(gp[i].x()) << ' '
+                    << float(gp[i].y()) << ' '
+                    << float(gp[i].z())
+                    << nl;
+
+                Pout << gp[i] << endl;
+            }
+            */
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::sharedPoints::sharedPoints(const PtrList<fvMesh>& meshes)
+Foam::sharedPoints::sharedPoints
+(
+    const PtrList<fvMesh>& meshes
+)
 :
     meshes_(meshes),
     sharedPointAddr_(meshes_.size()),
@@ -463,11 +614,19 @@ Foam::sharedPoints::sharedPoints(const PtrList<fvMesh>& meshes)
     calcSharedPoints();
 }
 
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-// Foam::sharedPoints::~sharedPoints()
-// {}
+Foam::sharedPoints::sharedPoints
+(
+    const PtrList<fvMesh>& meshes,
+    const labelListList& globalPointIndex
+)
+:
+    meshes_(meshes),
+    sharedPointAddr_(meshes_.size()),
+    sharedPointLabels_(meshes_.size()),
+    nGlobalPoints_(0)
+{
+    calcSharedPoints(globalPointIndex);
+}
 
 
 // ************************************************************************* //
diff --git a/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/sharedPoints.H b/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/sharedPoints.H
index 555931114..5b31e9376 100644
--- a/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/sharedPoints.H
+++ b/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/sharedPoints.H
@@ -117,13 +117,26 @@ class sharedPoints
         //- Calculate shared points
         void calcSharedPoints();
 
+        //- Calculate shared points
+        void calcSharedPoints(const labelListList& globalPointIndex);
+
 
 public:
 
     // Constructors
 
         //- Construct from the list of meshes
-        sharedPoints(const PtrList<fvMesh>& meshes);
+        sharedPoints
+        (
+            const PtrList<fvMesh>& meshes
+        );
+
+        //- Construct from the list of meshes and global point addressing
+        sharedPoints
+        (
+            const PtrList<fvMesh>& meshes,
+            const labelListList& globalPointIndex
+        );
 
 
     //- Destructor
diff --git a/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C b/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C
index 5ae37fd8e..680e7186b 100644
--- a/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C
+++ b/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C
@@ -271,13 +271,13 @@ bool Foam::dynamicPolyRefinementFvMesh::update()
             (
                 refinementSelectionPtr_->refinementCellCandidates()()
             );
-            Pout<< "Selected " << refCandidates.size()
+            Info<< "Selected " << refCandidates.size()
                 << " refinement candidates."
                 << endl;
         }
         else
         {
-            Pout<< "Skipping refinement for this time-step..." << endl;
+            Info<< "Skipping refinement for this time-step..." << endl;
         }
 
         // Set cells to refine. Note: refinement needs to make sure that face
@@ -297,13 +297,13 @@ bool Foam::dynamicPolyRefinementFvMesh::update()
             (
                 refinementSelectionPtr_->unrefinementPointCandidates()()
             );
-            Pout<< "Selected " << unrefCandidates.size()
+            Info<< "Selected " << unrefCandidates.size()
                 << " unrefinement candidates."
                 << endl;
         }
         else
         {
-            Pout<< "Skipping unrefinement for this time-step..." << endl;
+            Info<< "Skipping unrefinement for this time-step..." << endl;
         }
 
         // Set split points to unrefine around.
@@ -344,13 +344,13 @@ bool Foam::dynamicPolyRefinementFvMesh::update()
         // some topo changes
         if (sizeCellMap)
         {
-            Pout<< "Successfully performed polyhedral refinement. "
+            Info<< "Successfully performed polyhedral refinement. "
                 << "Changed from " << nOldCells << " to " << sizeCellMap
                 << " cells." << endl;
         }
         else
         {
-            Pout<< "Refinement/unrefinement not performed in this time step "
+            Info<< "Refinement/unrefinement not performed in this time step "
                 << "since no cells were selected." << endl;
         }
 
@@ -363,7 +363,7 @@ bool Foam::dynamicPolyRefinementFvMesh::update()
         curTimeIndex_ = time().timeIndex();
     }
 
-    Pout<< "No refinement/unrefinement" << endl;
+    Info<< "No refinement/unrefinement" << endl;
 
     // No refinement/unrefinement at this time step. Return false
     return false;
diff --git a/src/dynamicMesh/topoChangerFvMesh/topoChangerFvMesh/topoChangerFvMeshLoadBalance.C b/src/dynamicMesh/topoChangerFvMesh/topoChangerFvMesh/topoChangerFvMeshLoadBalance.C
index 1cde2cc69..163ee34b6 100644
--- a/src/dynamicMesh/topoChangerFvMesh/topoChangerFvMesh/topoChangerFvMeshLoadBalance.C
+++ b/src/dynamicMesh/topoChangerFvMesh/topoChangerFvMesh/topoChangerFvMeshLoadBalance.C
@@ -33,6 +33,7 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "cloud.H"
 #include "cloudDistribute.H"
+#include "meshObjectBase.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
@@ -138,10 +139,14 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
     // Prepare receiving side
 
     // Create the reconstructor
-    processorMeshesReconstructor meshRecon("reconstructed");
+    //
+    // HR 21.12.18 : Use empty domainname to avoid auto-created of
+    // fvSchemes/fvSolution
+    processorMeshesReconstructor meshRecon("");
 
     PtrList<fvMesh>& procMeshes = meshRecon.meshes();
     procMeshes.setSize(meshDecomp.nProcs());
+    meshRecon.globalPointIndex().setSize(meshDecomp.nProcs());
 
     // Collect local fields for decomposition
     clearOut();
@@ -242,6 +247,29 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
         }
     }
 
+    // HR 14.12.18: Create dummy database pointing into the non-parallel case
+    // directory and disable the runTimeModifiable switch.  dummyTime is used
+    // for decomposed, received and reconstructed fvMeshes (ie. before its data
+    // is moved into *this).
+    //
+    // The pointing into the non-parallel case directory is somewhat a hack to
+    // prevent auto-creation of fvSchemes and fvSolution (they exist in the root
+    // case). This is potentially dangerous if we write anything, but fvMeshes
+    // derived from this database are NO_WRITE so we should be fine. Making if
+    // non-runTimeModifiable prevents registration of fvSolution with the
+    // fileMonitor (addWatch). Not all processors will necessarily receive a
+    // mesh and watches will then cause dead-locks!
+    Time dummyTime
+    (
+        time().rootPath(),
+        time().globalCaseName(),
+        "system",
+        "constant",
+        false
+    );
+
+    const_cast<Switch&>(dummyTime.runTimeModifiable()) = false;
+
     for (label procI = 0; procI < meshDecomp.nProcs(); procI++)
     {
         // Check if there is a mesh to send
@@ -253,8 +281,9 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
             autoPtr<fvMesh> procMeshPtr = meshDecomp.processorMesh
             (
                 procI,
-                time(),
-                "processorPart" + Foam::name(procI),
+                dummyTime,
+                "",     // HR 21.12.18 : Use empty domainname to avoid
+                        // auto-created offvSchemes/fvSolution
                 true    // Create passive processor patches
             );
             fvMesh& procMesh = procMeshPtr();
@@ -286,6 +315,7 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
 
                 // Send the mesh and fields to target processor
                 toProc << procMesh << nl;
+                toProc << meshDecomp.globalPointIndex(procI) << nl;
 
                 // Send fields
                 sendFields(volScalarFields, fieldDecomposer, toProc);
@@ -342,6 +372,9 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
                     procMeshPtr
                 );
 
+                meshRecon.globalPointIndex()[procI] =
+                    meshDecomp.globalPointIndex(procI);
+
                 // Set local fields
                 // Note: first index is field index and second index is procI
                 insertFields
@@ -411,6 +444,10 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
                 }
 
                 //HJ Insert clouds missing.  HJ, 12/Oct/2018
+                //
+                // HR 18.11.18 - Not missing. Step is trivial and is treated in
+                // Cloud<ParticleType>::split which is called in the constructor
+                // of CloudDistibute
             }
         }
     }
@@ -445,9 +482,9 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
                     (
                         IOobject
                         (
-                            "processorPart" + Foam::name(procI),
-                            time().timeName(),
-                            time(),
+                            "", // "processorPart" + Foam::name(procI),
+                            dummyTime.timeName(),
+                            dummyTime,
                             IOobject::NO_READ,
                             IOobject::NO_WRITE
                         ),
@@ -456,6 +493,10 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
                     )
                 );
 
+                // Receive the global points addr
+                labelList gppi(fromProc);
+                meshRecon.globalPointIndex()[procI] = gppi;
+
                 // Receive fields
                 // Note: first index is field index and second index is procI
                 receiveFields
@@ -535,7 +576,7 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
 
     // Create the reconstructed mesh
     autoPtr<fvMesh> reconstructedMeshPtr =
-        meshRecon.reconstructMesh(time());
+        meshRecon.reconstructMesh(dummyTime);
     fvMesh& reconMesh = reconstructedMeshPtr();
 
     Pout<< "Reconstructed mesh stats: "
@@ -915,6 +956,9 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
         );
     }
 
+    // HR 13.12.18: Update the mesh objects
+    meshObjectBase::allUpdateTopology<polyMesh>(*this, meshMap);
+
     // Debug: remove?  HJ, 22/Oct/2018
     // checkMesh(true);
 
diff --git a/src/dynamicMesh/topoChangerFvMesh/topoChangerFvMesh/topoChangerFvMeshLoadBalanceTemplates.C b/src/dynamicMesh/topoChangerFvMesh/topoChangerFvMesh/topoChangerFvMeshLoadBalanceTemplates.C
index 406b2eff3..ee2db7c26 100644
--- a/src/dynamicMesh/topoChangerFvMesh/topoChangerFvMesh/topoChangerFvMeshLoadBalanceTemplates.C
+++ b/src/dynamicMesh/topoChangerFvMesh/topoChangerFvMesh/topoChangerFvMeshLoadBalanceTemplates.C
@@ -126,7 +126,7 @@ void Foam::topoChangerFvMesh::receiveFields
     receivedFields.setSize(nScalarFields);
 
     fromProc.readBegin("topoChangerFvMeshReceiveFields");
-    
+
     forAll (receivedFields, fI)
     {
         word fieldName(fromProc);
@@ -263,12 +263,12 @@ void Foam::topoChangerFvMesh::rebuildFields
         // Resize internal field
         if
         (
-            masterField.internalField().size()
+            masterField.size()
          != GeoField::GeoMeshType::size(masterField.mesh())
         )
         {
             Pout<< "Resizing internal field: old size = "
-                << masterField.internalField().size()
+                << masterField.size()
                 << " new size = "
                 << GeoField::GeoMeshType::size(masterField.mesh())
                 << endl;
@@ -280,7 +280,7 @@ void Foam::topoChangerFvMesh::rebuildFields
 
         // Resize boundary (number of patches)
         typename GeoField::GeometricBoundaryField& patchFields =
-                masterField.boundaryField();
+            masterField.boundaryFieldNoStoreOldTimes();
 
         if (patchFields.size() != masterField.mesh().boundary().size())
         {
@@ -352,6 +352,22 @@ void Foam::topoChangerFvMesh::rebuildFields
         fieldI++;
         Pout<< "... done" << endl;
     }
+
+    // HR 14.12.18: We create new processor boundary faces from internal
+    // faces. The values on these faces could be initialised by interpolation.
+    // Instead we choose to fix the values by evaluating the boundaries.
+    // I tried to execute evaluateCoupled() at the end of
+    // fvFieldReconstructor::reconstructField, but this fails in a strange way.
+    forAllConstIter
+    (
+        typename HashTable<const GeoField*>,
+        geoFields,
+        iter
+    )
+    {
+        GeoField& masterField = const_cast<GeoField&>(*iter());
+        masterField.boundaryField().evaluateCoupled();
+    }
 }
 
 
diff --git a/src/finiteVolume/cfdTools/compressible/compressibleCourantNo.H b/src/finiteVolume/cfdTools/compressible/compressibleCourantNo.H
index 6c0ef2b54..3c972a10e 100644
--- a/src/finiteVolume/cfdTools/compressible/compressibleCourantNo.H
+++ b/src/finiteVolume/cfdTools/compressible/compressibleCourantNo.H
@@ -33,7 +33,12 @@ scalar CoNum = 0.0;
 scalar meanCoNum = 0.0;
 scalar velMag = 0.0;
 
-if (mesh.nInternalFaces())
+// HR 26.06.18: A parallel run has at least two cells and therefore at least
+// one internal face in the global mesh. It may be a processor boundary, but
+// this is captured by max(mag(phi)).
+// Old formulation hangs on parallel cases where one partition is degenerated
+// to a single cell.
+if (mesh.nInternalFaces() || Pstream::parRun())
 {
     surfaceScalarField phiOverRho = mag(phi)/fvc::interpolate(rho);
 
@@ -47,6 +52,20 @@ if (mesh.nInternalFaces())
 
     velMag = max(phiOverRho/mesh.magSf()).value();
 }
+else
+{
+    // Single cell mesh: Co is still defined; use cell formulation
+
+    const scalar deltaT = runTime.deltaT().value();
+
+    const scalar deltaX = Foam::cbrt(mesh.V()[0]);
+
+    velMag = mag(U[0]);
+
+    CoNum = velMag*deltaT/deltaX;
+
+    meanCoNum = CoNum;
+}
 
 Info<< "Courant Number mean: " << meanCoNum
     << " max: " << CoNum
diff --git a/src/finiteVolume/cfdTools/incompressible/CourantNo.H b/src/finiteVolume/cfdTools/incompressible/CourantNo.H
index b1e05b72a..d94892406 100644
--- a/src/finiteVolume/cfdTools/incompressible/CourantNo.H
+++ b/src/finiteVolume/cfdTools/incompressible/CourantNo.H
@@ -33,7 +33,7 @@ scalar CoNum = 0.0;
 scalar meanCoNum = 0.0;
 scalar velMag = 0.0;
 
-// HR 26.06.28: A parallel run has at least two cells and therefore at least
+// HR 26.06.18: A parallel run has at least two cells and therefore at least
 // one internal face in the global mesh. It may be a processor boundary, but
 // this is captured by max(mag(phi)).
 // Old formulation hangs on parallel cases where one partition is degenerated
diff --git a/src/finiteVolume/fields/fvPatchFields/basic/fixedGradient/fixedGradientFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/basic/fixedGradient/fixedGradientFvPatchField.C
index 429548c2d..fa1a73466 100644
--- a/src/finiteVolume/fields/fvPatchFields/basic/fixedGradient/fixedGradientFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/basic/fixedGradient/fixedGradientFvPatchField.C
@@ -67,11 +67,11 @@ fixedGradientFvPatchField<Type>::fixedGradientFvPatchField
     const dictionary& dict
 )
 :
-    fvPatchField<Type>(p, iF, dict),
+    // HR 15.12.18: Must not call evaluate during construction. Read value
+    // instead. This is needed for PLB.
+    fvPatchField<Type>(p, iF, dict, true),
     gradient_("gradient", dict, p.size())
-{
-    evaluate();
-}
+{}
 
 
 template<class Type>
diff --git a/src/finiteVolume/fields/fvPatchFields/basic/mixed/mixedFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/basic/mixed/mixedFvPatchField.C
index 72abd4f07..e542db57b 100644
--- a/src/finiteVolume/fields/fvPatchFields/basic/mixed/mixedFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/basic/mixed/mixedFvPatchField.C
@@ -54,13 +54,13 @@ mixedFvPatchField<Type>::mixedFvPatchField
     const dictionary& dict
 )
 :
-    fvPatchField<Type>(p, iF, dict),
+    // HR 15.12.18: Must not call evaluate during construction. Read value
+    // instead. This is needed for PLB.
+    fvPatchField<Type>(p, iF, dict, true),
     refValue_("refValue", dict, p.size()),
     refGrad_("refGradient", dict, p.size()),
     valueFraction_("valueFraction", dict, p.size())
-{
-    evaluate();
-}
+{}
 
 
 template<class Type>
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/inletOutlet/inletOutletFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/inletOutlet/inletOutletFvPatchField.C
index 9f0d8bbbb..147d141d1 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/inletOutlet/inletOutletFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/inletOutlet/inletOutletFvPatchField.C
@@ -62,7 +62,7 @@ inletOutletFvPatchField<Type>::inletOutletFvPatchField
 {
     // Read patch type
     this->readPatchType(dict);
-    
+
     this->refValue() = Field<Type>("inletValue", dict, p.size());
 
     if (dict.found("value"))
diff --git a/src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.C b/src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.C
index c9050cd4d..b0533cd44 100644
--- a/src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.C
+++ b/src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.C
@@ -70,7 +70,12 @@ Foam::fvSchemes::fvSchemes(const objectRegistry& obr)
             "fvSchemes",
             obr.time().system(),
             obr,
-            IOobject::READ_IF_PRESENT_IF_MODIFIED,  // Allow default dictionary creation
+             (
+                obr.readOpt() == IOobject::MUST_READ
+             || obr.readOpt() == IOobject::READ_IF_PRESENT
+              ? IOobject::MUST_READ_IF_MODIFIED
+              : obr.readOpt()
+            ),
             IOobject::NO_WRITE
         )
     ),
@@ -168,21 +173,17 @@ Foam::fvSchemes::fvSchemes(const objectRegistry& obr)
     fluxRequired_(),  // Do not read flux required option
     defaultFluxRequired_(false)
 {
-    if (!headerOk())
+    // HR 21.12.18 : Writing a default fvSchemes is not useful in PLB. Please
+    // specify MUST_READ on obr if you need this.
+    if
+    (
+        readOpt() == IOobject::MUST_READ
+     || readOpt() == IOobject::MUST_READ_IF_MODIFIED
+     || (readOpt() == IOobject::READ_IF_PRESENT && headerOk())
+    )
     {
-        if (debug)
-        {
-            InfoIn
-            (
-                "fvSchemes::fvSchemes(const objectRegistry& obr)"
-            )   << "fvSchemes dictionary not found.  Creating default."
-                << endl;
-        }
-
-        regIOobject::write();
+        read();
     }
-
-    read();
 }
 
 
diff --git a/src/finiteVolume/fvMesh/fvMesh.C b/src/finiteVolume/fvMesh/fvMesh.C
index 26a7a8e74..d4865544a 100644
--- a/src/finiteVolume/fvMesh/fvMesh.C
+++ b/src/finiteVolume/fvMesh/fvMesh.C
@@ -424,6 +424,9 @@ void Foam::fvMesh::resetFvPrimitives
     // Reset fvPatches  HJ, 16/Apr/2018
     boundary_.resetFvPatches(resetFvPatchFlag);
 
+    // HR 14.12.18: Indicate that the mesh is changing to e.g. nearWallDist
+    polyMesh::changing(true);
+
     // Clear all mesh data
     clearOut();
 }
diff --git a/src/finiteVolume/fvMesh/wallDist/nearWallDist.C b/src/finiteVolume/fvMesh/wallDist/nearWallDist.C
index a3de7479c..76f99052b 100644
--- a/src/finiteVolume/fvMesh/wallDist/nearWallDist.C
+++ b/src/finiteVolume/fvMesh/wallDist/nearWallDist.C
@@ -84,7 +84,7 @@ void Foam::nearWallDist::correct()
         {
             setSize(mesh_.boundary().size());
         }
-        
+
         // Update size of GeometricBoundaryField
         forAll (mesh_.boundary(), patchI)
         {
diff --git a/src/foam/Make/files b/src/foam/Make/files
index afee07ad2..c982c574d 100644
--- a/src/foam/Make/files
+++ b/src/foam/Make/files
@@ -447,6 +447,7 @@ $(globalMeshData)/globalMeshData.C
 $(globalMeshData)/globalPoints.C
 $(globalMeshData)/globalIndex.C
 $(globalMeshData)/globalProcFaceIndex.C
+$(globalMeshData)/globalProcPointIndex.C
 
 $(polyMesh)/syncTools/syncTools.C
 
diff --git a/src/foam/fields/GeometricFields/GeometricField/GeometricField.C b/src/foam/fields/GeometricFields/GeometricField/GeometricField.C
index 3ecdc54bb..6d80f868a 100644
--- a/src/foam/fields/GeometricFields/GeometricField/GeometricField.C
+++ b/src/foam/fields/GeometricFields/GeometricField/GeometricField.C
@@ -728,6 +728,16 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::boundaryField()
 }
 
 
+// Return reference to GeometricBoundaryField
+template<class Type, template<class> class PatchField, class GeoMesh>
+typename
+Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricBoundaryField&
+Foam::GeometricField<Type, PatchField, GeoMesh>::boundaryFieldNoStoreOldTimes()
+{
+    return boundaryField_;
+}
+
+
 // Store old-time field
 template<class Type, template<class> class PatchField, class GeoMesh>
 void Foam::GeometricField<Type, PatchField, GeoMesh>::storeOldTimes() const
diff --git a/src/foam/fields/GeometricFields/GeometricField/GeometricField.H b/src/foam/fields/GeometricFields/GeometricField/GeometricField.H
index 3f05f6419..e8028ae8c 100644
--- a/src/foam/fields/GeometricFields/GeometricField/GeometricField.H
+++ b/src/foam/fields/GeometricFields/GeometricField/GeometricField.H
@@ -447,6 +447,10 @@ public:
         //- Return reference to GeometricBoundaryField for const field
         inline const GeometricBoundaryField& boundaryField() const;
 
+        //- Return reference to GeometricBoundaryField without storing old
+        //- times
+        GeometricBoundaryField& boundaryFieldNoStoreOldTimes();
+
         //- Return the time index of the field
         inline label timeIndex() const;
 
diff --git a/src/foam/fields/cloud/cloud.C b/src/foam/fields/cloud/cloud.C
index 4c7681f41..0310f5442 100644
--- a/src/foam/fields/cloud/cloud.C
+++ b/src/foam/fields/cloud/cloud.C
@@ -64,6 +64,13 @@ Foam::autoPtr<Foam::cloudDistribute> Foam::cloud::cloudDist
 }
 
 
+Foam::labelList Foam::cloud::nParticlesPerCell() const
+{
+    NotImplemented;
+    return labelList(0);
+}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::cloud::~cloud()
diff --git a/src/foam/fields/cloud/cloud.H b/src/foam/fields/cloud/cloud.H
index 613bed946..01e8d7a08 100644
--- a/src/foam/fields/cloud/cloud.H
+++ b/src/foam/fields/cloud/cloud.H
@@ -101,7 +101,9 @@ public:
             //- Return size of the cloud
             virtual label size() const;
 
-    
+            //- Count and return number of particles per cell
+            virtual labelList nParticlesPerCell() const;
+
         // Edit
 
             //- Remap the cells of particles corresponding to the
diff --git a/src/foam/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.C b/src/foam/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.C
index 882a3cae5..50a994d7f 100644
--- a/src/foam/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.C
+++ b/src/foam/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.C
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "pointBoundaryMesh.H"
+#include "polyMesh.H"
 #include "polyBoundaryMesh.H"
 #include "facePointPatch.H"
 #include "globalPointPatch.H"
@@ -105,10 +106,39 @@ void Foam::pointBoundaryMesh::movePoints()
 }
 
 
-void Foam::pointBoundaryMesh::updateMesh()
+void Foam::pointBoundaryMesh::updateMesh
+(
+    const polyMesh& pMesh
+)
 {
+    const polyBoundaryMesh& pBoundary = pMesh.boundaryMesh();
     pointPatchList& patches = *this;
 
+    // 21.1.19 HR : Patches need to be recreated in PLB
+    patches.clear();
+    patches.setSize(pBoundary.size());
+
+    forAll(patches, patchI)
+    {
+        patches.set
+        (
+            patchI,
+            facePointPatch::New(pBoundary[patchI], *this).ptr()
+        );
+    }
+
+    // Add the globalPointPatch
+    if(pMesh.globalData().nGlobalPoints())
+    {
+        patches.setSize(pBoundary.size() + 1);
+
+        patches.set
+        (
+            patches.size() - 1,
+            new globalPointPatch(*this, patches.size() - 1)
+        );
+    }
+
     forAll(patches, patchi)
     {
         patches[patchi].initUpdateMesh();
diff --git a/src/foam/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.H b/src/foam/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.H
index 5e9fe7443..ed1199d9c 100644
--- a/src/foam/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.H
+++ b/src/foam/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.H
@@ -46,6 +46,7 @@ namespace Foam
 class pointMesh;
 class polyBoundaryMesh;
 class globalPointPatch;
+class polyMesh;
 
 /*---------------------------------------------------------------------------*\
                        Class pointBoundaryMesh Declaration
@@ -105,7 +106,10 @@ public:
         void movePoints();
 
         //- Correct polyBoundaryMesh after topology update
-        void updateMesh();
+        void updateMesh
+        (
+            const polyMesh& pMesh
+        );
 };
 
 
diff --git a/src/foam/meshes/pointMesh/pointMesh.C b/src/foam/meshes/pointMesh/pointMesh.C
index 2acfa7044..21367477a 100644
--- a/src/foam/meshes/pointMesh/pointMesh.C
+++ b/src/foam/meshes/pointMesh/pointMesh.C
@@ -122,7 +122,10 @@ bool Foam::pointMesh::updateMesh(const mapPolyMesh& mpm) const
 {
     // Casting const-ness to answer the interface of meshObject
     // HJ, 30/Aug/2010
-    const_cast<pointBoundaryMesh&>(boundary_).updateMesh();
+    const_cast<pointBoundaryMesh&>(boundary_).updateMesh
+    (
+        GeoMesh<polyMesh>::mesh_
+    );
 
     // Map all registered point fields
     mapFields(mpm);
diff --git a/src/foam/meshes/polyMesh/globalMeshData/globalProcFaceIndex.C b/src/foam/meshes/polyMesh/globalMeshData/globalProcFaceIndex.C
index 4f7b1cc79..149ff962a 100644
--- a/src/foam/meshes/polyMesh/globalMeshData/globalProcFaceIndex.C
+++ b/src/foam/meshes/polyMesh/globalMeshData/globalProcFaceIndex.C
@@ -118,7 +118,8 @@ void Foam::globalProcFaceIndex::calcFaceIndex()
 
                 OPstream toProc
                 (
-                    Pstream::nonBlocking,
+                    // HR 12.12.18: nonBlocking fails on PLB of Aachen bomb
+                    Pstream::blocking,
                     procPatch.neighbProcNo()
                 );
                 toProc<< curFaceLabels;
@@ -154,7 +155,8 @@ void Foam::globalProcFaceIndex::calcFaceIndex()
                 // Receive the data from master and insert into the list
                 IPstream fromProc
                 (
-                    Pstream::nonBlocking,
+                    // HR 12.12.18: nonBlocking fails on PLB of Aachen bomb
+                    Pstream::blocking,
                     procPatch.neighbProcNo()
                 );
 
diff --git a/src/foam/meshes/polyMesh/globalMeshData/globalProcFaceIndex.H b/src/foam/meshes/polyMesh/globalMeshData/globalProcFaceIndex.H
index f0703ae0c..2da075d12 100644
--- a/src/foam/meshes/polyMesh/globalMeshData/globalProcFaceIndex.H
+++ b/src/foam/meshes/polyMesh/globalMeshData/globalProcFaceIndex.H
@@ -120,7 +120,7 @@ public:
             return procFaceOffset_;
         }
 
-        //- Return lobal face label for all faces of current mesh
+        //- Return global face label for all faces of current mesh
         //  Sized to number of live faces in the mesh
         const labelList globalLabel() const
         {
diff --git a/src/foam/meshes/polyMesh/globalMeshData/globalProcPointIndex.C b/src/foam/meshes/polyMesh/globalMeshData/globalProcPointIndex.C
new file mode 100644
index 000000000..cc1953f0d
--- /dev/null
+++ b/src/foam/meshes/polyMesh/globalMeshData/globalProcPointIndex.C
@@ -0,0 +1,319 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.1
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "globalProcPointIndex.H"
+#include "processorPolyPatch.H"
+#include "globalMeshData.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::globalProcPointIndex::calcPointIndex()
+{
+    // Count and mark points in the following order
+    // 1. global points
+    // 2. processor patches
+    // 3. regular patches
+    // 4. anything else ("internal points")
+    //
+    // The marking is a follows
+    // -4 : Not visited yet
+    // -3 : Found in a regular patch
+    // -2 : Found in slave processor patch
+    // -1 : Found in a master processor patch
+    //  0 - nGlobalPoints-1 : Global point ID
+
+    // 1. Mark the global points
+    const labelList& spl = mesh_.globalData().sharedPointLabels();
+    const labelList& spa = mesh_.globalData().sharedPointAddr();
+
+    forAll(spl, spI)
+    {
+        globalLabel_[spl[spI]] = spa[spI];
+    }
+
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+    labelList procPointCounter(Pstream::nProcs(), 0);
+    labelList patchPointCounter(patches.size(), 0);
+    label& nProcPoints = procPointCounter[Pstream::myProcNo()];
+
+    // Master processor patches - first pass
+    forAll (patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+        const labelList& patchMeshPoints = pp.meshPoints();
+
+        if (isA<processorPolyPatch>(pp))
+        {
+            const processorPolyPatch& procPatch =
+                refCast<const processorPolyPatch>(pp);
+
+            label& nPatchPoints = patchPointCounter[patchI];
+
+            if (procPatch.master())
+            {
+                forAll (patchMeshPoints, mpI)
+                {
+                    label& pl = globalLabel_[patchMeshPoints[mpI]];
+
+                    if (pl < 0)
+                    {
+                        nProcPoints++;
+                        nPatchPoints++;
+                        pl = -1;
+                    }
+                }
+            }
+        }
+    }
+
+    // Slave processor patches - first pass
+    forAll (patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+        const labelList& patchMeshPoints = pp.meshPoints();
+
+        if (isA<processorPolyPatch>(pp))
+        {
+            const processorPolyPatch& procPatch =
+                refCast<const processorPolyPatch>(pp);
+
+            if (!procPatch.master())
+            {
+                forAll (patchMeshPoints, mpI)
+                {
+                    label& pl = globalLabel_[patchMeshPoints[mpI]];
+
+                    if (pl == -4)
+                    {
+                        pl = -2;
+                    }
+                }
+            }
+        }
+    }
+
+    // Regular patches - first pass
+    forAll (patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+        const labelList& patchMeshPoints = pp.meshPoints();
+
+        if (!isA<processorPolyPatch>(pp))
+        {
+            label& nPatchPoints = patchPointCounter[patchI];
+
+            forAll (patchMeshPoints, mpI)
+            {
+                label& pl = globalLabel_[patchMeshPoints[mpI]];
+
+                if (pl == -4)
+                {
+                    nProcPoints++;
+                    nPatchPoints++;
+                    pl = -3;
+                }
+            }
+        }
+    }
+
+    // "Internal" points - first pass
+    label nInternalPoints = 0;
+    forAll (globalLabel_, pI)
+    {
+        if (globalLabel_[pI] == -4)
+        {
+            nInternalPoints++;
+            nProcPoints++;
+        }
+    }
+
+    // Gather-Scatter counter data
+    Pstream::gatherList(procPointCounter);
+    Pstream::scatterList(procPointCounter);
+
+    // Convert counters to offsets
+    procPointOffset_[0] = mesh_.globalData().nGlobalPoints();
+    for (label procI = 1; procI < procPointOffset_.size(); procI++)
+    {
+        procPointOffset_[procI] =
+            procPointCounter[procI-1] + procPointOffset_[procI-1];
+    }
+
+    patchPointOffset_[0] =
+        procPointOffset_[Pstream::myProcNo()] + nInternalPoints;
+    for (label patchI = 1; patchI < patchPointOffset_.size(); patchI++)
+    {
+        patchPointOffset_[patchI] =
+            patchPointCounter[patchI-1] + patchPointOffset_[patchI-1];
+    }
+
+    const pointField& points = mesh_.points();
+
+    // Master processor patches - second pass
+    // Send patch offset to slave and assign global labels to points on
+    // master processor patches and regular patches
+    forAll (patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+        const labelList& patchMeshPoints = pp.meshPoints();
+        const label patchPO = patchPointOffset_[patchI];
+
+        if (isA<processorPolyPatch>(pp))
+        {
+            const processorPolyPatch& procPatch =
+                refCast<const processorPolyPatch>(pp);
+
+            if (procPatch.master())
+            {
+                // Master processor patch
+                // Send the offset to the slave and mark the points through a
+                // face loop to establish an order that can be untangled on
+                // slave side
+
+                OPstream toProc
+                (
+                    Pstream::nonBlocking,
+                    procPatch.neighbProcNo()
+                );
+                toProc<< patchPO;
+
+                pointField pointLocs(patchMeshPoints.size());
+                label glPointLabelsI = 0;
+                forAll (pp, faceI)
+                {
+                    const face& curFace = pp[faceI];
+
+                    forAll (curFace, fpI)
+                    {
+                        const label pointI = curFace[fpI];
+                        label& pl = globalLabel_[pointI];
+
+                        if (pl == -1)
+                        {
+                            pl = patchPO + glPointLabelsI;
+                            pointLocs[glPointLabelsI] = points[pointI];
+                            glPointLabelsI++;
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    // Slave processor and regular patches - second pass
+    // Receive data on slave and assign global labels to points on
+    // master processor patches and regular patches
+    forAll (patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (isA<processorPolyPatch>(pp))
+        {
+            const processorPolyPatch& procPatch =
+                refCast<const processorPolyPatch>(pp);
+
+            if (!procPatch.master())
+            {
+                // Slave processor patch - second pass
+                // Receive the offset from master and mark the points by taking
+                // into account that the points in the faces are in reverse
+                // order
+
+                IPstream fromProc
+                (
+                    Pstream::nonBlocking,
+                    procPatch.neighbProcNo()
+                );
+
+                label masterPatchPO(readLabel(fromProc));
+
+                label glPointLabelsI = 0;
+                forAll (pp, faceI)
+                {
+                    const face curFace = pp[faceI].reverseFace();
+
+                    forAll (curFace, fpI)
+                    {
+                        const label pointI = curFace[fpI];
+                        label& pl = globalLabel_[pointI];
+
+                        if (pl == -2)
+                        {
+                            pl = masterPatchPO + glPointLabelsI;
+                            glPointLabelsI++;
+                        }
+                    }
+                }
+            }
+        }
+        else
+        {
+            // Regular patch - second pass
+            const labelList& patchMeshPoints = pp.meshPoints();
+            const label patchPO = patchPointOffset_[patchI];
+            label glPointLabelsI = 0;
+            forAll (patchMeshPoints, mpI)
+            {
+                label& pl = globalLabel_[patchMeshPoints[mpI]];
+
+                if (pl == -3)
+                {
+                    pl = patchPO + glPointLabelsI;
+                    glPointLabelsI++;
+                }
+            }
+        }
+    }
+
+    // "Internal" points - second pass
+    nInternalPoints = 0;
+    forAll (globalLabel_, pI)
+    {
+        label& pl = globalLabel_[pI];
+
+        if (pl == -4)
+        {
+            pl = procPointOffset_[Pstream::myProcNo()] + nInternalPoints;
+            nInternalPoints++;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::globalProcPointIndex::globalProcPointIndex(const polyMesh& mesh)
+:
+    mesh_(mesh),
+    procPointOffset_(Pstream::nProcs()+1),
+    patchPointOffset_(mesh_.boundaryMesh().size()+1),
+    globalLabel_(mesh_.nPoints(), -4)
+{
+    calcPointIndex();
+}
+
+
+// ************************************************************************* //
diff --git a/src/foam/meshes/polyMesh/globalMeshData/globalProcPointIndex.H b/src/foam/meshes/polyMesh/globalMeshData/globalProcPointIndex.H
new file mode 100644
index 000000000..d72fa04c6
--- /dev/null
+++ b/src/foam/meshes/polyMesh/globalMeshData/globalProcPointIndex.H
@@ -0,0 +1,135 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.1
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::globalProcPointIndex
+
+Description
+    The class creates a unique global point index for each processor points
+    (pair) in the mesh. Master and slave processor points carry the same index.
+
+    Point offsets counts a number of unique points on each processor,
+    excluding slave processor patch points and global points.
+
+    Currently, faces are ordered with internal faces first, followed by patch
+    faces in patch order, excluding slave processor patches.
+    If needed, this can be changed to group processor faces with internal faces
+    to facilitate parallel I/O.  HJ, 4/May/2018
+
+Author
+    Henrik Rusche, Wikki GmbH
+
+SourceFiles
+    globalProcPointIndex.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef globalProcPointIndex_H
+#define globalProcPointIndex_H
+
+#include "polyMesh.H"
+#include "labelList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+
+/*---------------------------------------------------------------------------*\
+                      Class globalProcPointIndex Declaration
+\*---------------------------------------------------------------------------*/
+
+class globalProcPointIndex
+{
+    // Private data
+
+        //- Mesh reference
+        const polyMesh& mesh_;
+
+        //- Processor point index offset
+        labelList procPointOffset_;
+
+        //- point index offset for patches in local mesh
+        labelList patchPointOffset_;
+
+        //- Global poins label for all points of current mesh
+        //  Sized to number of live points in the mesh
+        labelList globalLabel_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        globalProcPointIndex(const globalProcPointIndex&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const globalProcPointIndex&);
+
+        //- Calculate point index
+        void calcPointIndex();
+
+
+public:
+
+    // Constructors
+
+        //- Construct from mesh
+        globalProcPointIndex(const polyMesh&);
+
+
+    //- Destructor - default
+
+
+    // Member Functions
+
+        //- Return point index offset per processor
+        inline const labelList& procPointOffset() const
+        {
+            return procPointOffset_;
+        }
+
+        inline const labelList& patchPointOffset() const
+        {
+            return patchPointOffset_;
+        }
+
+        //- Return global point label for all points of current mesh
+        //  Sized to number of live points in the mesh
+        const labelList globalLabel() const
+        {
+            return globalLabel_;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/immersedBoundary/immersedBoundaryDynamicMesh/immersedBoundaryDynamicRefineSolidBodyMotionFvMesh/immersedBoundaryDynamicRefineSolidBodyMotionFvMesh.C b/src/immersedBoundary/immersedBoundaryDynamicMesh/immersedBoundaryDynamicRefineSolidBodyMotionFvMesh/immersedBoundaryDynamicRefineSolidBodyMotionFvMesh.C
index 2b57b7cd4..bef23df19 100644
--- a/src/immersedBoundary/immersedBoundaryDynamicMesh/immersedBoundaryDynamicRefineSolidBodyMotionFvMesh/immersedBoundaryDynamicRefineSolidBodyMotionFvMesh.C
+++ b/src/immersedBoundary/immersedBoundaryDynamicMesh/immersedBoundaryDynamicRefineSolidBodyMotionFvMesh/immersedBoundaryDynamicRefineSolidBodyMotionFvMesh.C
@@ -127,7 +127,7 @@ bool immersedBoundaryDynamicRefineSolidBodyMotionFvMesh::update()
     {
         hasChanged = loadBalance(refinementDict());
     }
-    
+
     // Execute dummy mesh motion for the background mesh
     const pointField oldPoints = allPoints();
     fvMesh::movePoints(oldPoints);
diff --git a/src/immersedBoundary/immersedBoundaryTurbulence/wallFunctions/immersedBoundaryEpsilonWallFunction/immersedBoundaryEpsilonWallFunctionFvPatchScalarField.C b/src/immersedBoundary/immersedBoundaryTurbulence/wallFunctions/immersedBoundaryEpsilonWallFunction/immersedBoundaryEpsilonWallFunctionFvPatchScalarField.C
index 1d0d64fdf..164aa35f7 100644
--- a/src/immersedBoundary/immersedBoundaryTurbulence/wallFunctions/immersedBoundaryEpsilonWallFunction/immersedBoundaryEpsilonWallFunctionFvPatchScalarField.C
+++ b/src/immersedBoundary/immersedBoundaryTurbulence/wallFunctions/immersedBoundaryEpsilonWallFunction/immersedBoundaryEpsilonWallFunctionFvPatchScalarField.C
@@ -68,7 +68,7 @@ immersedBoundaryEpsilonWallFunctionFvPatchScalarField
         p,
         Switch(dict.lookup("setDeadValue")),
         readScalar(dict.lookup("deadValue"))
-    )    
+    )
 {}
 
 
diff --git a/src/immersedBoundary/immersedBoundaryTurbulence/wallFunctions/immersedBoundaryKqRWallFunction/immersedBoundaryKqRWallFunctionFvPatchField.C b/src/immersedBoundary/immersedBoundaryTurbulence/wallFunctions/immersedBoundaryKqRWallFunction/immersedBoundaryKqRWallFunctionFvPatchField.C
index 4fc53997a..b2d580158 100644
--- a/src/immersedBoundary/immersedBoundaryTurbulence/wallFunctions/immersedBoundaryKqRWallFunction/immersedBoundaryKqRWallFunctionFvPatchField.C
+++ b/src/immersedBoundary/immersedBoundaryTurbulence/wallFunctions/immersedBoundaryKqRWallFunction/immersedBoundaryKqRWallFunctionFvPatchField.C
@@ -150,7 +150,7 @@ immersedBoundaryKqRWallFunctionFvPatchField
     (
         ptf.ibPatch(),
         ptf.setDeadValue(),
-        ptf.deadValue()        
+        ptf.deadValue()
     )
 {
     this->checkType();
diff --git a/src/immersedBoundary/immersedBoundaryTurbulence/wallFunctions/immersedBoundaryOmegaWallFunctions/immersedBoundaryOmegaWallFunctionFvPatchScalarField.C b/src/immersedBoundary/immersedBoundaryTurbulence/wallFunctions/immersedBoundaryOmegaWallFunctions/immersedBoundaryOmegaWallFunctionFvPatchScalarField.C
index 2e8fc4c1e..eddb71c93 100644
--- a/src/immersedBoundary/immersedBoundaryTurbulence/wallFunctions/immersedBoundaryOmegaWallFunctions/immersedBoundaryOmegaWallFunctionFvPatchScalarField.C
+++ b/src/immersedBoundary/immersedBoundaryTurbulence/wallFunctions/immersedBoundaryOmegaWallFunctions/immersedBoundaryOmegaWallFunctionFvPatchScalarField.C
@@ -66,7 +66,7 @@ immersedBoundaryOmegaWallFunctionFvPatchScalarField
         p,
         Switch(dict.lookup("setDeadValue")),
         readScalar(dict.lookup("deadValue"))
-    )    
+    )
 {}
 
 
diff --git a/src/lagrangian/basic/CloudDistributeTemplate/CloudDistributeTemplate.H b/src/lagrangian/basic/CloudDistributeTemplate/CloudDistributeTemplate.H
index be366e58e..600bbb674 100644
--- a/src/lagrangian/basic/CloudDistributeTemplate/CloudDistributeTemplate.H
+++ b/src/lagrangian/basic/CloudDistributeTemplate/CloudDistributeTemplate.H
@@ -47,7 +47,7 @@ class Cloud;
 
 class polyMesh;
 class cloud;
-    
+
 /*---------------------------------------------------------------------------*\
                            Class Cloud Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/lagrangian/basic/CloudTemplate/CloudTemplate.C b/src/lagrangian/basic/CloudTemplate/CloudTemplate.C
index c8bf946d2..848d48158 100644
--- a/src/lagrangian/basic/CloudTemplate/CloudTemplate.C
+++ b/src/lagrangian/basic/CloudTemplate/CloudTemplate.C
@@ -506,6 +506,40 @@ void Foam::Cloud<ParticleType>::rebuild
 }
 
 
+template<class ParticleType>
+Foam::labelList Foam::Cloud<ParticleType>::nParticlesPerCell() const
+{
+    Pout << "In nParticlesPerCell" << endl;
+    labelList nppc (pMesh().nCells(), 0);
+
+    forAllConstIter(typename Cloud<ParticleType>, *this, pIter)
+    {
+        const ParticleType& p = pIter();
+        const label celli = p.cell();
+
+        // Check
+        if (celli < 0 || celli >= pMesh().nCells())
+        {
+            FatalErrorIn
+            (
+                "Foam::Cloud<ParticleType>::nParticlesPerCell()"
+            )
+                << "Illegal cell number " << celli
+                << " at position " << p.position() << nl
+                << "Cell number should be between 0 and "
+                << pMesh().nCells()-1 << nl
+                << "On this mesh the particle should be in cell "
+                << pMesh().findCell(p.position())
+                << exit(FatalError);
+        }
+
+        nppc[celli]++;
+    }
+
+    return nppc;
+}
+
+
 template<class ParticleType>
 void Foam::Cloud<ParticleType>::writePositions() const
 {
diff --git a/src/lagrangian/basic/CloudTemplate/CloudTemplate.H b/src/lagrangian/basic/CloudTemplate/CloudTemplate.H
index e5716ed85..ae5700216 100644
--- a/src/lagrangian/basic/CloudTemplate/CloudTemplate.H
+++ b/src/lagrangian/basic/CloudTemplate/CloudTemplate.H
@@ -290,6 +290,9 @@ public:
                 const PtrList<labelIOList>& faceProcAddressing
             );
 
+            //- Count and return number of particles per cell
+            virtual labelList nParticlesPerCell() const;
+
 
         // Read
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
index aaf119dc5..542c4a8eb 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
@@ -370,7 +370,7 @@ void Foam::InjectionModel<CloudType>::inject(TrackData& td)
     // constrainToMeshCentre is inside of an if-statement
     mesh.geometricD();
     mesh.bounds();
-    
+
     // Prepare for next time step
     label parcelsAdded = 0;
     scalar massAdded = 0.0;
diff --git a/src/overset/oversetMesh/donorAcceptor/donorAcceptor.H b/src/overset/oversetMesh/donorAcceptor/donorAcceptor.H
index 153b41335..e298792ff 100644
--- a/src/overset/oversetMesh/donorAcceptor/donorAcceptor.H
+++ b/src/overset/oversetMesh/donorAcceptor/donorAcceptor.H
@@ -134,7 +134,8 @@ public:
             donorProcNo_(-1),
             donorPoint_(vector::zero),
             extendedDonorCells_(),
-            extendedDonorPoints_()
+            extendedDonorPoints_(),
+            withinBB_(false)
         {}
 
         //- Construct from acceptor data
diff --git a/src/overset/oversetMesh/oversetFringe/overlapFringe/adaptiveOverlapFringe/adaptiveOverlapFringe.C b/src/overset/oversetMesh/oversetFringe/overlapFringe/adaptiveOverlapFringe/adaptiveOverlapFringe.C
index 8ae0fcc87..98282dee7 100644
--- a/src/overset/oversetMesh/oversetFringe/overlapFringe/adaptiveOverlapFringe/adaptiveOverlapFringe.C
+++ b/src/overset/oversetMesh/oversetFringe/overlapFringe/adaptiveOverlapFringe/adaptiveOverlapFringe.C
@@ -67,54 +67,44 @@ void Foam::adaptiveOverlapFringe::suitabilityFractionSlope
     scalar& beta
 ) const
 {
-    // Linear regression coefficients equations implemented down there:
-    // beta = (n*sum(x_i*y_i)-sum(x_i)*sum(y_i)) / (n*sum(x_i^2)-(sum(x_i))^2)
-    // alpha = sum(y_i)/n - beta* (sum(x_i)/n
+    // Linear regression coefficients y = alpha + beta*x
+    // beta = sum((x_i - x_mean)*(y_i - y_mean))/sum(x_i - x_mean)^2
+    // alpha = y_mean - beta*x_mean
+    // x -> iteration number
+    // y -> suitability
 
-    // Iteration ordinal numbers sum (iteration ordinal number represents x
-    // value in x-y coordinate system)
-    scalar iterSum = 0;
+    // Calculate mean values first
+    scalar iterMean = 0;
+    scalar suitabilityMean = 0;
 
-    // Sum of squared iteration ordinal numbers
-    scalar iterSquaredSum = 0;
-
-    // Average donor/acceptor suitability fraction values sum (average
-    // suitability represents y value in x-y coordinate system)
-    scalar suitabilitySum = 0;
-
-    // Sum of iteration ordinal number multiplied by donor/acceptor suitability
-    // fraction value from one iteration
-    scalar iterXSuitabilitySum = 0;
-
-    // Loop through all iteration data objects and calculate values needed for
-    // linear regression calculation
     forAllConstIter(FIFOStack<iterationData>, iterHist, it)
     {
-        iterSum += it().iteration();
-        iterSquaredSum += pow(it().iteration(), 2);
-        suitabilitySum += it().suitability();
-        iterXSuitabilitySum += it().iteration()*it().suitability();
+        iterMean += it().iteration();
+        suitabilityMean += it().suitability();
     }
 
-    // Number of iterations is equal to iterationDataHistory_ size
-    const label n = iterationDataHistory_.size();
+    reduce(iterMean, sumOp<scalar>());
+    iterMean /= iterHist.size();
+    reduce(suitabilityMean, sumOp<scalar>());
+    suitabilityMean /= iterHist.size();
 
-    // Calculate slope (i.e. if y = alpha + beta*x, slope is beta)
-    beta =
-        (n*iterXSuitabilitySum - iterSum*suitabilitySum)/
-        (n*iterSquaredSum - pow(iterSum, 2));
+    // Calculate denominator and numerator for beta
+    scalar n = 0;
+    scalar d = 0;
 
-    // Calculate mean iteration number - actually we don't need to calculate
-    // it beacuse we are only interested in trend, i.e we need only slope
-    // (beta). Maybe I will need it for diagrams in my master thesis, so
-    // calculate it for now.
-    const scalar meanIter = iterSum/n;
+    forAllConstIter(FIFOStack<iterationData>, iterHist, it)
+    {
+        n += (it().iteration() - iterMean)*
+             (it().suitability() - suitabilityMean);
+        d += Foam::sqr((it().iteration() - iterMean));
+    }
 
-    // Calculate mean suitability from all iterations
-    const scalar meanSuitability = suitabilitySum/n;
+    reduce(n, sumOp<scalar>());
+    reduce(d, sumOp<scalar>());
 
-    // Calculate intercept
-    alpha = meanSuitability - beta*meanIter;
+    // Calculate regression coefficients
+    beta = n/d;
+    alpha = suitabilityMean - beta*iterMean;
 }
 
 
@@ -179,20 +169,30 @@ void Foam::adaptiveOverlapFringe::calcAddressing() const
     // holes)
     boolList eligibleAcceptors(mesh.nCells(), true);
 
-    // Read user specified holes into allHoles list. Note: if the cell set
-    // is not found, the list will be empty
-    cellSet allHoles
-    (
-        mesh,
-        holesSetName_,
-        IOobject::READ_IF_PRESENT,
-        IOobject::NO_WRITE
-    );
+    // Read user specified holes into allHoles list. Note: use cellZone rather
+    // than cellSet to have correct behaviour on dynamic mesh simulations
+    // We will silently proceed if the zone is not found since this option is
+    // not mandatory but is useful in certain cases
+
+    // Get zone index
+    const label zoneID = mesh.cellZones().findZoneID(holesZoneName_);
+
+    // Create a hash set for allHoles
+    labelHashSet allHoles;
+
+    if (zoneID > -1)
+    {
+        // Get the zone for holes and append them to set
+        const labelList& specifiedHoles = mesh.cellZones()[zoneID];
+
+        allHoles.insert(specifiedHoles);
+    }
+    // else silently proceed without user-specified holes
 
     // Extend allHoles with cutHoles
     forAll (cutHoles, chI)
     {
-        // Note: cellSet is a hashSet so it automatically removes duplicates
+        // Note: duplicated are removed because we're using hash set
         allHoles.insert(cutHoles[chI]);
     }
 
@@ -393,7 +393,7 @@ Foam::adaptiveOverlapFringe::adaptiveOverlapFringe
     acceptorsPtr_(NULL),
     finalDonorAcceptorsPtr_(NULL),
 
-    holesSetName_(dict.lookupOrDefault<word>("holes", word())),
+    holesZoneName_(dict.lookupOrDefault<word>("holes", word())),
     initPatchNames_
     (
         dict.lookupOrDefault<wordList>("initPatchNames", wordList())
@@ -408,6 +408,14 @@ Foam::adaptiveOverlapFringe::adaptiveOverlapFringe
     (
         dict.lookupOrDefault<label>("specifiedIterationsNumber", 5)
     ),
+    minSuitabilityRate_
+    (
+        dict.lookupOrDefault<scalar>("minSuitabilityRate", 0.0)
+    ),
+    maxIter_
+    (
+        dict.lookupOrDefault<label>("maximumIterations", 100)
+    ),
     relativeCounter_(specifiedIterationsNumber_),
     additionalIterations_
     (
@@ -667,6 +675,7 @@ bool Foam::adaptiveOverlapFringe::updateIteration
         (additionalIterations_)
      && (fringeIter_ == relativeCounter_)
      && (specifiedIterationsNumber_ > 1)
+     && (fringeIter_ < maxIter_)
     )
     {
         // Slope (i.e. if y = alpha + beta*x, slope is beta)
@@ -695,7 +704,7 @@ bool Foam::adaptiveOverlapFringe::updateIteration
 
         // If the donor/acceptor suitability gradient is positive,
         // make 1 additional iteration
-        if (beta > 0)
+        if (beta > minSuitabilityRate_)
         {
             // Pop the bottom element of the stack
             iterationDataHistory_.pop();
@@ -710,7 +719,7 @@ bool Foam::adaptiveOverlapFringe::updateIteration
                     << " will be performed." << endl;
             }
         }
-     }
+    }
 
     // Go through unsuitable donor/acceptor pairs from previous iteration and
     // find a new batch of acceptors and holes for the next iteration if:
@@ -1076,16 +1085,26 @@ bool Foam::adaptiveOverlapFringe::updateIteration
         reduce(nAccToHoles, sumOp<label>());
         Info<< "Converted " << nAccToHoles << " acceptors to holes." << endl;
 
+        // Bugfix: Although we have found suitable overlap, we need to update
+        // acceptors as well because eligible donors for acceptors of other
+        // regions are calculated based on these acceptors (and holes)
+        labelList& acceptors = *acceptorsPtr_;
+        acceptors.setSize(finalDAPairs.size());
+        forAll (acceptors, aI)
+        {
+            acceptors[aI] = finalDAPairs[aI].acceptorCell();
+        }
+
         // Construct final donor/acceptor list
         finalDonorAcceptorsPtr_ = new donorAcceptorList
         (
-            finalDAPairs
+            finalDAPairs.xfer()
         );
 
         // Construct final fringe holes list
         fringeHolesPtr_ = new labelList
         (
-            fringeHoles
+            fringeHoles.xfer()
         );
 
         // Print information
diff --git a/src/overset/oversetMesh/oversetFringe/overlapFringe/adaptiveOverlapFringe/adaptiveOverlapFringe.H b/src/overset/oversetMesh/oversetFringe/overlapFringe/adaptiveOverlapFringe/adaptiveOverlapFringe.H
index f14b5830a..b33b7d548 100644
--- a/src/overset/oversetMesh/oversetFringe/overlapFringe/adaptiveOverlapFringe/adaptiveOverlapFringe.H
+++ b/src/overset/oversetMesh/oversetFringe/overlapFringe/adaptiveOverlapFringe/adaptiveOverlapFringe.H
@@ -105,7 +105,7 @@ class adaptiveOverlapFringe
             //- Name of the cell set defining initial holes (empty by default).
             //  Useful when the resolution of the background mesh is much
             //  coarser than the front mesh and no hole is found
-            const word holesSetName_;
+            const word holesZoneName_;
 
             //- Optional list of patches to start the iterative fringe assembly
             //  process (empty list by default). Useful when we actually have
@@ -125,6 +125,14 @@ class adaptiveOverlapFringe
             //  be made
             const label specifiedIterationsNumber_;
 
+            //- User specified suitability trend rate. If the trend is larger
+            //  than the specified trend rate, additional iterations are
+            //  performed
+            const scalar minSuitabilityRate_;
+
+            //- Maximum number of iterations
+            const label maxIter_;
+
             //- Relative iteration counter
             mutable label relativeCounter_;
 
diff --git a/src/overset/oversetMesh/oversetFringe/overlapFringe/donorSuitability/donorSuitability/donorSuitability.H b/src/overset/oversetMesh/oversetFringe/overlapFringe/donorSuitability/donorSuitability/donorSuitability.H
index c60210047..39a2b8cd8 100644
--- a/src/overset/oversetMesh/oversetFringe/overlapFringe/donorSuitability/donorSuitability/donorSuitability.H
+++ b/src/overset/oversetMesh/oversetFringe/overlapFringe/donorSuitability/donorSuitability/donorSuitability.H
@@ -224,6 +224,21 @@ public:
                         << " on processor: " << daPair.acceptorProcNo()
                         << " did not find donor candidate."
                         << nl
+                        << "Additional information: "
+                        << nl
+                        <<  (
+                                daPair.withinBB()
+                              ? "   Within bounding box"
+                              : "   Not within bounding box"
+                            )
+                        << nl << tab
+                        << "Donor index: " << daPair.donorCell()
+                        << nl << tab
+                        << "Donor processor index: " << daPair.donorProcNo()
+                        << nl << tab
+                        << "Number of extended donors: "
+                        << daPair.extendedDonorCells().size()
+                        << nl << nl
                         << "Please review your fringe assembly settings"
                         << " (or try using adaptiveOverlap fringe algorithm)."
                         << abort(FatalError);
diff --git a/src/overset/oversetMesh/oversetFringe/overlapFringe/overlapFringe/overlapFringe.C b/src/overset/oversetMesh/oversetFringe/overlapFringe/overlapFringe/overlapFringe.C
index f714b7d5f..f40f0a762 100644
--- a/src/overset/oversetMesh/oversetFringe/overlapFringe/overlapFringe/overlapFringe.C
+++ b/src/overset/oversetMesh/oversetFringe/overlapFringe/overlapFringe/overlapFringe.C
@@ -188,26 +188,37 @@ void Foam::overlapFringe::calcAddressing() const
     // holes)
     boolList eligibleAcceptors(mesh.nCells(), true);
 
-    // Read user specified holes into allHoles list. Note: if the cell set
-    // is not found, the list will be empty
-    labelList allHoles
-    (
-        cellSet
-        (
-            mesh,
-            holesSetName_,
-            IOobject::READ_IF_PRESENT,
-            IOobject::NO_WRITE
-        ).toc()
-    );
+    // Read user specified holes into allHoles list. Note: use cellZone rather
+    // than cellSet to have correct behaviour on dynamic mesh simulations
+    // We will silently proceed if the zone is not found since this option is
+    // not mandatory but is useful in certain cases
+
+    // Get zone index
+    const label zoneID = mesh.cellZones().findZoneID(holesZoneName_);
+
+    // Create a hash set for allHoles
+    labelHashSet allHoles;
+
+    if (zoneID > -1)
+    {
+        // Get the zone for holes and append them to set
+        const labelList& specifiedHoles = mesh.cellZones()[zoneID];
+
+        allHoles.insert(specifiedHoles);
+    }
+    // else silently proceed without user-specified holes
 
     // Extend allHoles with cutHoles
-    allHoles.append(cutHoles);
+    forAll (cutHoles, chI)
+    {
+        // Note: duplicated are removed because we're using hash set
+        allHoles.insert(cutHoles[chI]);
+    }
 
     // Mark all holes
-    forAll (allHoles, hI)
+    forAllConstIter (labelHashSet, allHoles, iter)
     {
-        const label& holeCellI = allHoles[hI];
+        const label& holeCellI = iter.key();
 
         // Mask eligible acceptors
         eligibleAcceptors[holeCellI] = false;
@@ -225,10 +236,10 @@ void Foam::overlapFringe::calcAddressing() const
     dynamicLabelList candidateAcceptors(mesh.nCells());
 
     // Loop through all holes and find acceptor candidates
-    forAll (allHoles, hI)
+    forAllConstIter (labelHashSet, allHoles, iter)
     {
         // Get neighbours of this hole cell
-        const labelList& hNbrs = cc[allHoles[hI]];
+        const labelList& hNbrs = cc[iter.key()];
 
         // Loop through neighbours of this hole cell
         forAll (hNbrs, nbrI)
@@ -370,7 +381,7 @@ void Foam::overlapFringe::calcAddressing() const
     // Transfer the acceptor list and allocate empty fringeHoles list, which
     // may be populated in updateIteration member function
     acceptorsPtr_ = new labelList(candidateAcceptors.xfer());
-    fringeHolesPtr_ = new labelList(cutHoles);
+    fringeHolesPtr_ = new labelList(allHoles.toc().xfer());
 }
 
 
@@ -398,7 +409,7 @@ Foam::overlapFringe::overlapFringe
     acceptorsPtr_(NULL),
     finalDonorAcceptorsPtr_(NULL),
 
-    holesSetName_(dict.lookupOrDefault<word>("holes", word())),
+    holesZoneName_(dict.lookupOrDefault<word>("holes", word())),
     initPatchNames_
     (
         dict.lookupOrDefault<wordList>("initPatchNames", wordList())
@@ -687,8 +698,15 @@ bool Foam::overlapFringe::updateIteration
             }
         }
 
-        // Tranfer back the allFringeHoles dynamic list into member data
-        fringeHolesPtr_->transfer(allFringeHoles);
+        // Bugfix: Although we have found suitable overlap, we need to update
+        // acceptors as well because eligible donors for acceptors of other
+        // regions are calculated based on these acceptors (and holes)
+        labelList& acceptors = *acceptorsPtr_;
+        acceptors.setSize(finalDAPairs.size());
+        forAll (acceptors, aI)
+        {
+            acceptors[aI] = finalDAPairs[aI].acceptorCell();
+        }
 
         // Transfer ownership of the final donor/acceptor list to the
         // finalDonorAcceptorsPtr_
@@ -697,6 +715,9 @@ bool Foam::overlapFringe::updateIteration
             finalDAPairs.xfer()
         );
 
+        // Tranfer back the allFringeHoles dynamic list into member data
+        fringeHolesPtr_->transfer(allFringeHoles);
+
         // At least 100*minGlobalFraction_ % of suitable donor/acceptor pairs
         // have been found.
         Info<< "Converted " << nAccToHoles << " acceptors to holes."
diff --git a/src/overset/oversetMesh/oversetFringe/overlapFringe/overlapFringe/overlapFringe.H b/src/overset/oversetMesh/oversetFringe/overlapFringe/overlapFringe/overlapFringe.H
index 40bb502f1..fe0b1ada3 100644
--- a/src/overset/oversetMesh/oversetFringe/overlapFringe/overlapFringe/overlapFringe.H
+++ b/src/overset/oversetMesh/oversetFringe/overlapFringe/overlapFringe/overlapFringe.H
@@ -96,7 +96,7 @@ class overlapFringe
             //- Name of the cell set defining initial holes (empty by default).
             //  Useful when the resolution of the background mesh is much
             //  coarser than the front mesh and no hole is found
-            const word holesSetName_;
+            const word holesZoneName_;
 
             //- Optional list of patches to start the iterative fringe assembly
             //  process (empty list by default). Useful when we actually have
diff --git a/src/overset/oversetMesh/oversetMesh/oversetMeshAddressing.C b/src/overset/oversetMesh/oversetMesh/oversetMeshAddressing.C
index bd922b8de..6e8cdfe80 100644
--- a/src/overset/oversetMesh/oversetMesh/oversetMeshAddressing.C
+++ b/src/overset/oversetMesh/oversetMesh/oversetMeshAddressing.C
@@ -55,18 +55,33 @@ void Foam::oversetMesh::calcCellClassification() const
         nHoleCells += regions_[regionI].holes().size();
     }
 
-    Pout<< "Number of acceptor cells: " << nAcceptorCells << endl;
     acceptorCellsPtr_ = new labelList(nAcceptorCells);
     labelList& acceptor = *acceptorCellsPtr_;
 
-    Pout<< "Number of donor cells: " << nDonorCells << endl;
     donorCellsPtr_ = new labelList(nDonorCells);
     labelList& donor = *donorCellsPtr_;
 
-    Pout<< "Number of hole cells: " << nHoleCells << endl;
     holeCellsPtr_ = new labelList(nHoleCells);
     labelList& hole = *holeCellsPtr_;
 
+    // Print out processor specific information in debug
+    if (oversetMesh::debug)
+    {
+        Pout<< "Number of acceptor cells: " << nAcceptorCells << endl;
+        Pout<< "Number of donor cells: " << nDonorCells << endl;
+        Pout<< "Number of hole cells: " << nHoleCells << endl;
+    }
+    // Else print global information
+    else
+    {
+        Info<< "Number of acceptor cells: "
+            << returnReduce(nAcceptorCells, sumOp<label>()) << endl;
+        Info<< "Number of donor cells: "
+            << returnReduce(nDonorCells, sumOp<label>()) << endl;
+        Info<< "Number of hole cells: "
+            << returnReduce(nHoleCells, sumOp<label>()) << endl;
+    }
+
     // Reset counters
     nAcceptorCells = 0;
     nDonorCells = 0;
diff --git a/src/overset/oversetMesh/oversetRegion/oversetRegion.C b/src/overset/oversetMesh/oversetRegion/oversetRegion.C
index f261ec605..674e0a9c9 100644
--- a/src/overset/oversetMesh/oversetRegion/oversetRegion.C
+++ b/src/overset/oversetMesh/oversetRegion/oversetRegion.C
@@ -190,18 +190,9 @@ void Foam::oversetRegion::calcDonorAcceptorCells() const
             // Update global flag
             foundGlobalOverlap &= regionFoundSuitableOverlap;
 
-            // If the overlap has not been found for this region, we need to
-            // reset:
-            //  - holeCells (depend on fringe holes)
-            //  - eligibleDonors (depend on fringe holes and acceptors),
-            //  - cellSearch (depends on eligible donors).
-            if (!regionFoundSuitableOverlap)
-            {
-                deleteDemandDrivenData(curRegion.cellSearchPtr_);
-            }
-
-            deleteDemandDrivenData(curRegion.holeCellsPtr_);
-            deleteDemandDrivenData(curRegion.eligibleDonorCellsPtr_);
+            // Deletion of demand driven data relocated to the end of
+            // updateDonorAcceptors() member function. Makes more sense to have
+            // it there. VV, 13/Jan/2019.
         }
     } while (!foundGlobalOverlap);
 
@@ -1317,6 +1308,20 @@ bool Foam::oversetRegion::updateDonorAcceptors() const
 
     // STAGE 11: Filter possibly multiple remote donors
 
+    // Sanity check before moving on. For each initial acceptor, we must
+    // receive at least 1 corresponding donor (even if it is not found,
+    // indicating an orphan cell) from different processors.
+    if (a.size() > completeDonorAcceptorList.size())
+    {
+        FatalErrorIn("void oversetRegion::updateDonorAcceptors() const")
+            << "Size of initial acceptor set: " << a.size()
+            << " is larger than the size of the distributed "
+            << " donor/acceptor list: " << completeDonorAcceptorList.size()
+            << nl
+            << "This should not have happened..."
+            << abort(FatalError);
+    }
+
     // Create a masking field indicating that a certain acceptor has been
     // visited
     boolList isVisited(a.size(), false);
@@ -1379,35 +1384,16 @@ bool Foam::oversetRegion::updateDonorAcceptors() const
                     curDA.donorPoint(),
                     curDA.withinBB()
                 );
+
+                // Bugfix: also need to reset extended donors since a better
+                // candidate has been found. VV, 1/Jan/2019
+                curDACombined.setExtendedDonors(curDA);
             }
         }
     }
 
-    // Update withinBB flag if the donor is within bounding box of acceptor
-    // (previously we checked whether the acceptor is within bounding box of
-    // donor)
-    forAll (combinedDonorAcceptorList, daI)
-    {
-        donorAcceptor& curDA = combinedDonorAcceptorList[daI];
-
-        // If the acceptor is not within bounding box of donor, set the flag
-        // other way around
-        if (!curDA.withinBB())
-        {
-            curDA.setWithinBB
-            (
-                mesh_.pointInCellBB
-                (
-                    curDA.donorPoint(),
-                    curDA.acceptorCell()
-                )
-            );
-        }
-    }
-
-    // Check whether all acceptors have been visited. Used for testing/debugging
-    // parallel comms
-    if (oversetMesh::debug)
+    // Check whether all acceptors have been visited. Useful check if in
+    // no-debug mode
     {
         bool allVisited = true;
 
@@ -1435,12 +1421,35 @@ bool Foam::oversetRegion::updateDonorAcceptors() const
                     << nl
                     << "Try switching off useLocalBoundingBoxes for all regions"
                     << nl
-                    << "(this optimisation is switched on by default)."
+                    << "(this optimisation is switched off by default)."
                     << abort(FatalError);
             }
         }
     }
 
+    // Update withinBB flag if the donor is within bounding box of acceptor
+    // (previously we checked whether the acceptor is within bounding box of
+    // donor)
+    forAll (combinedDonorAcceptorList, daI)
+    {
+        donorAcceptor& curDA = combinedDonorAcceptorList[daI];
+
+        // If the acceptor is not within bounding box of donor, set the flag
+        // other way around
+        if (!curDA.withinBB())
+        {
+            curDA.setWithinBB
+            (
+                mesh_.pointInCellBB
+                (
+                    curDA.donorPoint(),
+                    curDA.acceptorCell()
+                )
+            );
+        }
+    }
+
+
     // STAGE 12: Finish the iteration by updating the fringe, which will
     // actually hold final and some intermediate steps for donor/acceptor
     // assembly process
@@ -1450,6 +1459,13 @@ bool Foam::oversetRegion::updateDonorAcceptors() const
     bool suitableOverlapFound =
         fringePtr_->updateIteration(combinedDonorAcceptorList);
 
+    // Delete all necessary demand driven data for this region since we have
+    // just updated the iteration. Therefore, cellSearch, eligibleDonors and
+    // holes need to be updated
+    deleteDemandDrivenData(cellSearchPtr_);
+    deleteDemandDrivenData(eligibleDonorCellsPtr_);
+    deleteDemandDrivenData(holeCellsPtr_);
+
     return suitableOverlapFound;
 }
 
diff --git a/src/sampling/sampledSet/curve/curveSet.C b/src/sampling/sampledSet/curve/curveSet.C
index cfdbd5133..f8f9cd745 100644
--- a/src/sampling/sampledSet/curve/curveSet.C
+++ b/src/sampling/sampledSet/curve/curveSet.C
@@ -65,7 +65,7 @@ bool Foam::curveSet::trackToBoundary
         const vector offset =
             sampleCoords_[sampleI + 1] - sampleCoords_[sampleI];
 
-        const scalar smallDist = mag(tol*offset);
+        const scalar smallDist = mag(tol_()*offset);
 
         point oldPos = trackPt;
         label facei = -1;
@@ -174,7 +174,7 @@ void Foam::curveSet::calcSamples
         {
             const vector offset =
                 sampleCoords_[sampleI+1] - sampleCoords_[sampleI];
-            const scalar smallDist = mag(tol*offset);
+            const scalar smallDist = mag(tol_()*offset);
 
 
             // Get all boundary intersections
diff --git a/src/sampling/sampledSet/face/faceOnlySet.C b/src/sampling/sampledSet/face/faceOnlySet.C
index 3ddc58bf1..629871f2b 100644
--- a/src/sampling/sampledSet/face/faceOnlySet.C
+++ b/src/sampling/sampledSet/face/faceOnlySet.C
@@ -43,6 +43,16 @@ namespace Foam
 }
 
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+const Foam::debug::optimisationSwitch
+Foam::faceOnlySet::maxNSteps_
+(
+    "maximumFaceOnlySetStepNumber",
+    100
+);
+
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 
@@ -58,7 +68,7 @@ bool Foam::faceOnlySet::trackToBoundary
 {
     // distance vector between sampling points
     const vector offset = end_ - start_;
-    const vector smallVec = tol*offset;
+    const vector smallVec = tol_()*offset;
     const scalar smallDist = mag(smallVec);
 
     // Alias
@@ -115,7 +125,7 @@ void Foam::faceOnlySet::calcSamples
 
     const vector offset = (end_ - start_);
     const vector normOffset = offset/mag(offset);
-    const vector smallVec = tol*offset;
+    const vector smallVec = tol_()*offset;
     const scalar smallDist = mag(smallVec);
 
 
@@ -198,7 +208,7 @@ void Foam::faceOnlySet::calcSamples
     // index in bHits; current boundary intersection
     label bHitI = 1;
 
-    while(true)
+    while(segmentI <= maxNSteps_())
     {
         if (trackFaceI != -1)
         {
diff --git a/src/sampling/sampledSet/face/faceOnlySet.H b/src/sampling/sampledSet/face/faceOnlySet.H
index 712039b8d..0b7b7fd2f 100644
--- a/src/sampling/sampledSet/face/faceOnlySet.H
+++ b/src/sampling/sampledSet/face/faceOnlySet.H
@@ -63,6 +63,13 @@ class faceOnlySet
         point end_;
 
 
+    // Static data
+
+        //- Maximum number of steps for intersection search. Introduced to avoid
+        //  inifinite loops (IG 4/Dec/2018)
+        static const debug::optimisationSwitch maxNSteps_;
+
+
     // Private Member Functions
 
         //- Samples from startTrackPt/CellI. Updates particle/samplePt/sampleI
diff --git a/src/sampling/sampledSet/sampledSet/sampledSet.C b/src/sampling/sampledSet/sampledSet/sampledSet.C
index 12628065f..d6aa7d05f 100644
--- a/src/sampling/sampledSet/sampledSet/sampledSet.C
+++ b/src/sampling/sampledSet/sampledSet/sampledSet.C
@@ -33,13 +33,29 @@ License
 
 namespace Foam
 {
-    const scalar sampledSet::tol = 1e-6;
-
     defineTypeNameAndDebug(sampledSet, 0);
     defineRunTimeSelectionTable(sampledSet, word);
 }
 
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+const Foam::debug::tolerancesSwitch
+Foam::sampledSet::tol_
+(
+    "sampleSetTolerance",
+    1e-5
+);
+
+
+const Foam::debug::tolerancesSwitch
+Foam::sampledSet::pushTol_
+(
+    "sampleSetPushTolerance",
+    0.1
+);
+
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 Foam::label Foam::sampledSet::getBoundaryCell(const label faceI) const
@@ -189,7 +205,7 @@ Foam::point Foam::sampledSet::pushIn
     const point& cellCtr = mesh().cellCentres()[cellI];
 
     point newSample =
-        facePt + tol*(cellCtr - facePt);
+        facePt + pushTol_()*(cellCtr - facePt);
 
     if (!searchEngine().pointInCell(newSample, cellI))
     {
@@ -221,7 +237,7 @@ bool Foam::sampledSet::getTrackingPoint
     label& trackFaceI
 ) const
 {
-    const scalar smallDist = mag(tol*offset);
+    const scalar smallDist = mag(tol_()*offset);
 
     bool isGoodSample = false;
 
diff --git a/src/sampling/sampledSet/sampledSet/sampledSet.H b/src/sampling/sampledSet/sampledSet/sampledSet.H
index f7ba35aba..c81c279eb 100644
--- a/src/sampling/sampledSet/sampledSet/sampledSet.H
+++ b/src/sampling/sampledSet/sampledSet/sampledSet.H
@@ -48,6 +48,8 @@ SourceFiles
 #include "runTimeSelectionTables.H"
 #include "autoPtr.H"
 #include "coordSet.H"
+#include "debug.H"
+#include "tolerancesSwitch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -201,7 +203,11 @@ public:
 
         //- Tolerance when comparing points. Usually relative to difference
         //  between start_ and end_
-        const static scalar tol;
+        static const debug::tolerancesSwitch tol_;
+
+        //- Tolerance used to push points out of the face. It is relative to the
+        //  distance between the cell face centre and cell centre.
+        static const debug::tolerancesSwitch pushTol_;
 
 
     // Constructors
diff --git a/src/sampling/sampledSet/uniform/uniformSet.C b/src/sampling/sampledSet/uniform/uniformSet.C
index 0ce568d6b..7d4c0b0f6 100644
--- a/src/sampling/sampledSet/uniform/uniformSet.C
+++ b/src/sampling/sampledSet/uniform/uniformSet.C
@@ -96,7 +96,7 @@ bool Foam::uniformSet::trackToBoundary
 {
     // distance vector between sampling points
     const vector offset = (end_ - start_)/(nPoints_ - 1);
-    const vector smallVec = tol*offset;
+    const vector smallVec = tol_()*offset;
     const scalar smallDist = mag(smallVec);
 
     // Alias
@@ -233,7 +233,7 @@ void Foam::uniformSet::calcSamples
 
     const vector offset = (end_ - start_)/(nPoints_ - 1);
     const vector normOffset = offset/mag(offset);
-    const vector smallVec = tol*offset;
+    const vector smallVec = tol_()*offset;
     const scalar smallDist = mag(smallVec);
 
     // Get all boundary intersections
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.C
index bcc190429..bdfbad2b0 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.C
@@ -57,8 +57,21 @@ Foam::basicChemistryModel::basicChemistryModel
     chemistry_(lookup("chemistry")),
     deltaTChem_
     (
-        mesh.nCells(),
-        readScalar(lookup("initialChemicalTimeStep"))
+        IOobject
+        (
+            "deltaTChem",
+            mesh.time().constant(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar
+        (
+            "deltaTChem0",
+            dimTime,
+            readScalar(lookup("initialChemicalTimeStep"))
+        )
     )
 {}
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H
index 0f6e20bc9..67ee7fd9d 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H
@@ -38,7 +38,7 @@ SourceFiles
 
 #include "IOdictionary.H"
 #include "Switch.H"
-#include "scalarField.H"
+#include "volFields.H"
 #include "volFieldsFwd.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -77,7 +77,9 @@ protected:
         Switch chemistry_;
 
         //- Latest estimation of integration step
-        scalarField deltaTChem_;
+        // HR 14.12.18: Upgraded to volScalarField in order to enable automatic
+        // mapping/resizing after a topo change. Also needed for PLB.
+        volScalarField deltaTChem_;
 
 
     // Protected member functions
diff --git a/src/turbulenceModels/incompressible/RAS/RASModel/RASModel.H b/src/turbulenceModels/incompressible/RAS/RASModel/RASModel.H
index c1be1f43e..20d3b948c 100644
--- a/src/turbulenceModels/incompressible/RAS/RASModel/RASModel.H
+++ b/src/turbulenceModels/incompressible/RAS/RASModel/RASModel.H
@@ -298,7 +298,7 @@ public:
         ) const;
 
         //- Solve the turbulence equations and correct the turbulence viscosity
-        virtual void correct() = 0;
+        virtual void correct();
 
         //- Read RASProperties dictionary
         virtual bool read() = 0;
diff --git a/src/turbulenceModels/incompressible/turbulenceModel/turbulenceModel.H b/src/turbulenceModels/incompressible/turbulenceModel/turbulenceModel.H
index 8508129bc..1460becd3 100644
--- a/src/turbulenceModels/incompressible/turbulenceModel/turbulenceModel.H
+++ b/src/turbulenceModels/incompressible/turbulenceModel/turbulenceModel.H
@@ -207,7 +207,7 @@ public:
         virtual tmp<fvVectorMatrix> divDevReff() const = 0;
 
         //- Solve the turbulence equations and correct the turbulence viscosity
-        virtual void correct() = 0;
+        virtual void correct();
 
         //- Read turbulenceProperties dictionary
         virtual bool read() = 0;