diff --git a/src/ODE/sixDOF/sixDOFODE/constraints/rotationalConstraints/rotationalConstraint/rotationalConstraint.H b/src/ODE/sixDOF/sixDOFODE/constraints/rotationalConstraints/rotationalConstraint/rotationalConstraint.H
index 4bbf61634..72081fbd4 100644
--- a/src/ODE/sixDOF/sixDOFODE/constraints/rotationalConstraints/rotationalConstraint/rotationalConstraint.H
+++ b/src/ODE/sixDOF/sixDOFODE/constraints/rotationalConstraints/rotationalConstraint/rotationalConstraint.H
@@ -72,8 +72,15 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class sixDOFODE;
 
+// Forward declaration of friend functions and operators
+class rotationalConstraint;
+
+Ostream& operator<<(Ostream&, const rotationalConstraint&);
+
+
 /*---------------------------------------------------------------------------*\
                    Class rotationalConstraint Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/ODE/sixDOF/sixDOFODE/constraints/translationalConstraints/translationalConstraint/translationalConstraint.H b/src/ODE/sixDOF/sixDOFODE/constraints/translationalConstraints/translationalConstraint/translationalConstraint.H
index 08104e909..f72760082 100644
--- a/src/ODE/sixDOF/sixDOFODE/constraints/translationalConstraints/translationalConstraint/translationalConstraint.H
+++ b/src/ODE/sixDOF/sixDOFODE/constraints/translationalConstraints/translationalConstraint/translationalConstraint.H
@@ -72,8 +72,15 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class sixDOFODE;
 
+// Forward declaration of friend functions and operators
+class translationalConstraint;
+
+Ostream& operator<<(Ostream&, const translationalConstraint&);
+
+
 /*---------------------------------------------------------------------------*\
                    Class translationalConstraint Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/ODE/sixDOF/sixDOFODE/restraints/combinedRestraints/combinedRestraint/combinedRestraint.H b/src/ODE/sixDOF/sixDOFODE/restraints/combinedRestraints/combinedRestraint/combinedRestraint.H
index 15519b7a3..45fce5fd8 100644
--- a/src/ODE/sixDOF/sixDOFODE/restraints/combinedRestraints/combinedRestraint/combinedRestraint.H
+++ b/src/ODE/sixDOF/sixDOFODE/restraints/combinedRestraints/combinedRestraint/combinedRestraint.H
@@ -54,8 +54,15 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class sixDOFODE;
 
+// Forward declaration of friend functions and operators
+class combinedRestraint;
+
+Ostream& operator<<(Ostream&, const combinedRestraint&);
+
+
 /*---------------------------------------------------------------------------*\
                    Class combinedRestraint Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/ODE/sixDOF/sixDOFODE/restraints/rotationalRestraints/rotationalRestraint/rotationalRestraint.H b/src/ODE/sixDOF/sixDOFODE/restraints/rotationalRestraints/rotationalRestraint/rotationalRestraint.H
index 6aed49eef..0ab11cd15 100644
--- a/src/ODE/sixDOF/sixDOFODE/restraints/rotationalRestraints/rotationalRestraint/rotationalRestraint.H
+++ b/src/ODE/sixDOF/sixDOFODE/restraints/rotationalRestraints/rotationalRestraint/rotationalRestraint.H
@@ -54,8 +54,15 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class sixDOFODE;
 
+// Forward declaration of friend functions and operators
+class rotationalRestraint;
+
+Ostream& operator<<(Ostream&, const rotationalRestraint&);
+
+
 /*---------------------------------------------------------------------------*\
                    Class rotationalRestraint Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/ODE/sixDOF/sixDOFODE/restraints/translationalRestraints/translationalRestraint/translationalRestraint.H b/src/ODE/sixDOF/sixDOFODE/restraints/translationalRestraints/translationalRestraint/translationalRestraint.H
index bfb1f434f..32ae439b4 100644
--- a/src/ODE/sixDOF/sixDOFODE/restraints/translationalRestraints/translationalRestraint/translationalRestraint.H
+++ b/src/ODE/sixDOF/sixDOFODE/restraints/translationalRestraints/translationalRestraint/translationalRestraint.H
@@ -54,8 +54,15 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class sixDOFODE;
 
+// Forward declaration of friend functions and operators
+class translationalRestraint;
+
+Ostream& operator<<(Ostream&, const translationalRestraint&);
+
+
 /*---------------------------------------------------------------------------*\
                    Class translationalRestraint Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H
index e84157b92..a945c6c93 100644
--- a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H
+++ b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H
@@ -52,6 +52,12 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+class sixDOFqODE;
+
+Ostream& operator<<(Ostream&, const sixDOFqODE&);
+
+
 /*---------------------------------------------------------------------------*\
                          Class sixDOFqODE Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/ODE/translationODE/translationODE.H b/src/ODE/translationODE/translationODE.H
index d7d32cefc..18838a1e9 100644
--- a/src/ODE/translationODE/translationODE.H
+++ b/src/ODE/translationODE/translationODE.H
@@ -51,6 +51,12 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+class translationODE;
+
+Ostream& operator<<(Ostream&, const translationODE&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class translationODE Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/OSspecific/POSIX/memInfo/memInfo.H b/src/OSspecific/POSIX/memInfo/memInfo.H
index cf660d34f..4f7344491 100644
--- a/src/OSspecific/POSIX/memInfo/memInfo.H
+++ b/src/OSspecific/POSIX/memInfo/memInfo.H
@@ -47,6 +47,17 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
+class Istream;
+class Ostream;
+
+// Forward declaration of friend functions and operators
+class memInfo;
+
+Istream& operator>>(Istream&, memInfo&);
+Ostream& operator<<(Ostream&, const memInfo&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class memInfo Declaration
 \*---------------------------------------------------------------------------*/
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/conversion/ensight/part/ensightPart.H b/src/conversion/ensight/part/ensightPart.H
index eb60906bd..f2c609f07 100644
--- a/src/conversion/ensight/part/ensightPart.H
+++ b/src/conversion/ensight/part/ensightPart.H
@@ -51,6 +51,13 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+class ensightPart;
+
+Ostream& operator<<(Ostream&, const ensightPart&);
+ensightGeoFile& operator<<(ensightGeoFile&, const ensightPart&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class ensightPart Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/conversion/ensight/part/ensightParts.H b/src/conversion/ensight/part/ensightParts.H
index bbd4bde03..7a4b0c2e6 100644
--- a/src/conversion/ensight/part/ensightParts.H
+++ b/src/conversion/ensight/part/ensightParts.H
@@ -46,6 +46,12 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+class ensightParts;
+
+ensightGeoFile& operator<<(ensightGeoFile&, const ensightParts&);
+
+
 /*---------------------------------------------------------------------------*\
                        Class ensightParts Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/decomposeMesh.C b/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/decomposeMesh.C
index 3d0fab6bf..45c459f87 100644
--- a/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/decomposeMesh.C
+++ b/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/decomposeMesh.C
@@ -38,9 +38,83 @@ Description
 #include "primitiveMesh.H"
 #include "processorPolyPatch.H"
 #include "cyclicPolyPatch.H"
+#include "SortableList.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+bool Foam::domainDecomposition::addInterProcessorBoundaryData
+(
+    const label& procI,
+    const label& procJ,
+    const label faceI,
+
+    List<SLList<label> >& interProcBoundaries,
+    List<SLList<SLList< label> > >& interProcBFaces
+) const
+{
+    // Algorithm: Syncronously loop through the two lists containing all
+    // neighbouring processor for a given processor and check whether this
+    // processor is already in the list of neighbours. If the processor is
+    // already there, simply append the face. If it's not there, append the
+    // processor and the first face.
+
+    // Get iterators to beginning of the two lists
+    SLList<label>::iterator curInterProcBdrsIter =
+        interProcBoundaries[procI].begin();
+
+    SLList<SLList<label> >::iterator curInterProcBFacesIter =
+        interProcBFaces[procI].begin();
+
+    // Helper flag to distinguish when a particular neighbour processor has been
+    // encountered
+    bool interProcBouFound = false;
+
+    // WARNING: Synchronous SLList iterators: Assuming interProcBoundaries and
+    // interProcBFaces are always in sync as they should be
+    for
+    (
+        ;
+
+        curInterProcBdrsIter != interProcBoundaries[procI].end()
+     && curInterProcBFacesIter != interProcBFaces[procI].end();
+
+      ++curInterProcBdrsIter,
+      ++curInterProcBFacesIter
+    )
+    {
+        if (curInterProcBdrsIter() == procJ)
+        {
+            // Inter-processor boundary exists. Mark that we found this
+            // neighbour processor
+            interProcBouFound = true;
+
+            // Append the face into the list
+            curInterProcBFacesIter().append(faceI);
+
+            // Break out of the loop
+            break;
+        }
+    }
+
+    if (!interProcBouFound)
+    {
+        // Inter-processor boundary does not exist and needs to be created
+
+        // Append procJ to procI list
+        interProcBoundaries[procI].append(procJ);
+
+        // Append face to procI. Note: SLList construct taking a single
+        // label creates a list with only that element
+        interProcBFaces[procI].append(SLList<label>(faceI));
+    }
+
+    // Return whether the inter processor boundary has been found in the list
+    return interProcBouFound;
+}
+
+
+// * * * * * * * * * * * * * Public Member Functions * * * * * * * * * * * * //
+
 void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
 {
     // Decide which cell goes to which processor
@@ -48,8 +122,11 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
 
     // Distribute the cells according to the given processor label
 
-    // calculate the addressing information for the original mesh
-    Info<< "\nCalculating original mesh data" << endl;
+    // Calculate the addressing information for the original mesh
+    if (debug)
+    {
+        Pout<< "\nCalculating original mesh data" << endl;
+    }
 
     // Set references to the original mesh
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
@@ -62,7 +139,10 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
     // loop through the list of processor labels for the cell and add the
     // cell shape to the list of cells for the appropriate processor
 
-    Info<< "\nDistributing cells to processors" << endl;
+    if (debug)
+    {
+        Pout<< "\nDistributing cells to processors" << endl;
+    }
 
     // Memory management
     {
@@ -90,7 +170,10 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
         }
     }
 
-    Info << "\nDistributing faces to processors" << endl;
+    if (debug)
+    {
+        Pout << "\nDistributing faces to processors" << endl;
+    }
 
     // Loop through internal faces and decide which processor they belong to
     // First visit all internal faces. If cells at both sides belong to the
@@ -124,8 +207,6 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
         // Face labels belonging to each inter-processor boundary
         List<SLList<SLList<label> > > interProcBFaces(nProcs_);
 
-        List<SLList<label> > procPatchIndex(nProcs_);
-
         // Rewrite:
         // Handling of coupled patches is changed.  HJ, 11/Apr/2018
 
@@ -153,10 +234,10 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
         // When running the decomposition in parallel, it is assumed that
         // the result will be used in dynamic load balancing
         // For load balancing, the processor patches need to match in order for
-        // the patch-to-patch matching to work propely on mesh reconstruction
+        // the patch-to-patch matching to work properly on mesh reconstruction
         // Therefore, all processor patches need to be split into the matching
-        // and non-matching part be examining the cellToProc_ data on
-        // the neighbour side.  This has been prepared in patchNbrCellToProc_
+        // and non-matching part by examining the cellToProc_ data on
+        // the neighbour side. This has been prepared in patchNbrCellToProc_
         // HJ, 11/Apr/2018
 
         // The correct order of dumping the faces is to dump the
@@ -178,7 +259,7 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
                 const processorPolyPatch& procPatch =
                     refCast<const processorPolyPatch>(patches[patchI]);
 
-                // DO ONLY SLAVE SIDE
+                // Note: DO ONLY SLAVE SIDE
                 if (!procPatch.master())
                 {
                     // Get patch start
@@ -200,191 +281,67 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
 
                         // Check change in processor type across the processor
                         // boundary
-                        // If ownerProc and neighbourProc are the same, the
-                        // processor face will be merged, meaning that it
-                        // remains in the (old) processor patch
-                        // If ownerProc and neighbourProc are different,
+
+                        // If procI and procJ are the same, the processor face
+                        // will be merged, meaning that it remains in the (old)
+                        // processor patch. If procI and procJ are different,
                         // this will be a new processor boundary created from
-                        //  the existing processor face and added afterwards
+                        // the existing processor face and added afterwards
                         if (ownerProc != neighbourProc)
                         {
-                            // Search algorithm repeated in processor patches.
-                            // HJ, 11/Apr/2018
-                            SLList<label>::iterator curInterProcBdrsOwnIter =
-                                interProcBoundaries[ownerProc].begin();
-
-                            SLList<SLList<label> >::iterator
-                                curInterProcBFacesOwnIter =
-                                interProcBFaces[ownerProc].begin();
-
-                            bool interProcBouFound = false;
-
-                            // WARNING: Synchronous SLList iterators
-
-                            for
+                            // Insert inter-processor data for ownerProc
+                            addInterProcessorBoundaryData
                             (
-                                ;
-                                curInterProcBdrsOwnIter
-                             != interProcBoundaries[ownerProc].end()
-                             && curInterProcBFacesOwnIter
-                             != interProcBFaces[ownerProc].end();
-                                ++curInterProcBdrsOwnIter,
-                                ++curInterProcBFacesOwnIter
-                            )
-                            {
-                                if (curInterProcBdrsOwnIter() == neighbourProc)
-                                {
-                                    // Inter - processor boundary exists.
-                                    // Add the face
-                                    interProcBouFound = true;
+                                ownerProc,     // Processor to append to
+                                neighbourProc, // Processor to append
+                                patchStart + patchFaceI, // Face index to append
 
-                                    // Add to owner only!
-                                    curInterProcBFacesOwnIter().append
-                                    (
-                                        patchStart + patchFaceI
-                                    );
-                                }
-
-                                if (interProcBouFound) break;
-                            }
-
-                            if (!interProcBouFound)
-                            {
-                                // inter - processor boundaries do not exist and
-                                // need to be created
-
-                                // set the new addressing information
-
-                                // Add to owner only!
-                                interProcBoundaries[ownerProc].append
-                                (
-                                    neighbourProc
-                                );
-                                interProcBFaces[ownerProc].append
-                                (
-                                    SLList<label>(patchStart + patchFaceI)
-                                );
-                            }
+                                interProcBoundaries,
+                                interProcBFaces
+                            );
                         }
-                        // Note: cannot insert regular faces here, because
-                        // they are out of sequence.
-                        // HJ, 24/Apr/2018
-                    }
-                }
-            }
-        }
+
+                        // Note: cannot insert regular faces here, because they
+                        // are out of sequence.  HJ, 24/Apr/2018
+
+                    } // End for all patch faces
+                } // End if this is slave processor
+            } // End if this is a processor patch
+        } // End for all patches
 
         // Internal mesh faces
         forAll (neighbour, faceI)
         {
-            if (cellToProc_[owner[faceI]] != cellToProc_[neighbour[faceI]])
+            const label& ownerProc = cellToProc_[owner[faceI]];
+            const label& neighbourProc = cellToProc_[neighbour[faceI]];
+
+            // Check whether we'll end up on the new processor boundary (if
+            // ownerProc and neighbourProc are different)
+            if (ownerProc != neighbourProc)
             {
-                // Inter - processor patch face found. Go through the list of
-                // inside boundaries for the owner processor and try to find
-                // this inter-processor patch.
-
-                const label ownerProc = cellToProc_[owner[faceI]];
-                const label neighbourProc = cellToProc_[neighbour[faceI]];
-
-                // Search algorithm repeated in processor patches.  Reconsider
-                // HJ, 11/Apr/2018
-                SLList<label>::iterator curInterProcBdrsOwnIter =
-                    interProcBoundaries[ownerProc].begin();
-
-                SLList<SLList<label> >::iterator curInterProcBFacesOwnIter =
-                    interProcBFaces[ownerProc].begin();
-
-                bool interProcBouFound = false;
-
-                // WARNING: Synchronous SLList iterators
-
-                for
+                // Insert inter-processor data for ownerProc
+                addInterProcessorBoundaryData
                 (
-                    ;
-                    curInterProcBdrsOwnIter
-                 != interProcBoundaries[ownerProc].end()
-                 && curInterProcBFacesOwnIter
-                 != interProcBFaces[ownerProc].end();
-                    ++curInterProcBdrsOwnIter, ++curInterProcBFacesOwnIter
-                )
-                {
-                    if (curInterProcBdrsOwnIter() == neighbourProc)
-                    {
-                        // Inter - processor boundary exists. Add the face
-                        interProcBouFound = true;
+                    ownerProc,     // Processor to append to
+                    neighbourProc, // Processor to append
+                    faceI,         // Face index to append
 
-                        curInterProcBFacesOwnIter().append(faceI);
+                    interProcBoundaries,
+                    interProcBFaces
+                );
 
-                        SLList<label>::iterator curInterProcBdrsNeiIter =
-                            interProcBoundaries[neighbourProc].begin();
+                // Insert inter-processor data for neighbourProc
+                // Note: ownerProc and neighbourProc swapped compared to the
+                // call above
+                addInterProcessorBoundaryData
+                (
+                    neighbourProc, // Processor to append to
+                    ownerProc,     // Processor to append
+                    faceI,         // Face index to append
 
-                        SLList<SLList<label> >::iterator
-                            curInterProcBFacesNeiIter =
-                            interProcBFaces[neighbourProc].begin();
-
-                        bool neighbourFound = false;
-
-                        // WARNING: Synchronous SLList iterators
-
-                        for
-                        (
-                            ;
-                            curInterProcBdrsNeiIter !=
-                            interProcBoundaries[neighbourProc].end()
-                         && curInterProcBFacesNeiIter !=
-                            interProcBFaces[neighbourProc].end();
-                            ++curInterProcBdrsNeiIter,
-                            ++curInterProcBFacesNeiIter
-                        )
-                        {
-                            if (curInterProcBdrsNeiIter() == ownerProc)
-                            {
-                                // boundary found. Add the face
-                                neighbourFound = true;
-
-                                curInterProcBFacesNeiIter().append(faceI);
-                            }
-
-                            if (neighbourFound) break;
-                        }
-
-                        if (!neighbourFound)
-                        {
-                            // Owner found (slave proc boundary), but neighbour
-                            // not found: add it
-                            interProcBoundaries[neighbourProc].append
-                            (
-                                ownerProc
-                            );
-
-                            interProcBFaces[neighbourProc].append
-                            (
-                                SLList<label>(faceI)
-                            );
-                        }
-                    }
-
-                    if (interProcBouFound) break;
-                }
-
-                if (!interProcBouFound)
-                {
-                    // inter - processor boundaries do not exist and need to
-                    // be created
-
-                    // set the new addressing information
-
-                    // Add owner side
-                    interProcBoundaries[ownerProc].append(neighbourProc);
-                    interProcBFaces[ownerProc].append(SLList<label>(faceI));
-
-                    // Add neighbour side
-                    interProcBoundaries[neighbourProc].append(ownerProc);
-                    interProcBFaces[neighbourProc].append
-                    (
-                        SLList<label>(faceI)
-                    );
-                }
+                    interProcBoundaries,
+                    interProcBFaces
+                );
             }
         }
 
@@ -402,7 +359,7 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
                 const processorPolyPatch& procPatch =
                     refCast<const processorPolyPatch>(patches[patchI]);
 
-                // DO ONLY MASTER SIDE
+                // Note: DO ONLY MASTER SIDE
                 if (procPatch.master())
                 {
                     // Get patch start
@@ -424,79 +381,31 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
 
                         // Check change in processor type across the processor
                         // boundary
+
                         // If ownerProc and neighbourProc are the same, the
                         // processor face will be merged, meaning that it
                         // remains in the (old) processor patch
                         // If ownerProc and neighbourProc are different,
                         // this will be a new processor boundary created from
-                        //  the existing processor face and added afterwards
+                        // the existing processor face and added afterwards
                         if (ownerProc != neighbourProc)
                         {
-                            // Search algorithm repeated in processor patches.
-                            // HJ, 11/Apr/2018
-                            SLList<label>::iterator curInterProcBdrsOwnIter =
-                                interProcBoundaries[ownerProc].begin();
-
-                            SLList<SLList<label> >::iterator
-                                curInterProcBFacesOwnIter =
-                                interProcBFaces[ownerProc].begin();
-
-                            bool interProcBouFound = false;
-
-                            // WARNING: Synchronous SLList iterators
-
-                            for
+                            // Insert inter-processor data for ownerProc
+                            addInterProcessorBoundaryData
                             (
-                                ;
-                                curInterProcBdrsOwnIter
-                             != interProcBoundaries[ownerProc].end()
-                             && curInterProcBFacesOwnIter
-                             != interProcBFaces[ownerProc].end();
-                                ++curInterProcBdrsOwnIter,
-                                ++curInterProcBFacesOwnIter
-                            )
-                            {
-                                if (curInterProcBdrsOwnIter() == neighbourProc)
-                                {
-                                    // Inter - processor boundary exists.
-                                    // Add the face
-                                    interProcBouFound = true;
+                                ownerProc,     // Processor to append to
+                                neighbourProc, // Processor to append
+                                patchStart + patchFaceI, // Face index to append
 
-                                    // Add to owner only!
-                                    curInterProcBFacesOwnIter().append
-                                    (
-                                        patchStart + patchFaceI
-                                    );
-                                }
-
-                                if (interProcBouFound) break;
-                            }
-
-                            if (!interProcBouFound)
-                            {
-                                // inter - processor boundaries do not exist and
-                                // need to be created
-
-                                // set the new addressing information
-
-                                // Add to owner only!
-                                interProcBoundaries[ownerProc].append
-                                (
-                                    neighbourProc
-                                );
-                                interProcBFaces[ownerProc].append
-                                (
-                                    SLList<label>(patchStart + patchFaceI)
-                                );
-                            }
+                                interProcBoundaries,
+                                interProcBFaces
+                            );
                         }
-                        // Note: cannot insert regular faces here, because
-                        // they are out of sequence.
-                        // HJ, 24/Apr/2018
-                    }
-                }
-            }
-        }
+
+                    } // End for all patch faces
+                } // End if this is slave processor
+            } // End if this is a processor patch
+        } // End for all patches
 
         // Loop through patches. For cyclic boundaries detect inter-processor
         // faces; for all other, add faces to the face list and remember start
@@ -513,8 +422,8 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
 
             const label patchStart = patches[patchI].start();
 
-            // Do normal patches.  Note that processor patches
-            // have already been partially done and need special treatment
+            // Do normal patches. Note: processor patches have already been
+            // partially done and need special treatment
             if
             (
                 isA<processorPolyPatch>(patches[patchI])
@@ -522,31 +431,170 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
             )
             {
                 // Only collect faces where the owner and neighbour processor
-                // index are the same.
-                // If owner and neighbour processor index are different,
-                // the face was already collected into a separate patch
-                // HJ, 23/Apr/2018
+                // index are the same. If owner and neighbour processor index
+                // are different, the face was already collected into a separate
+                // patch. HJ, 23/Apr/2018.
 
+                // Cast to processor patch
+                const processorPolyPatch& procPatch =
+                    refCast<const processorPolyPatch>(patches[patchI]);
+
+                // Is this slave processor?
+                const bool isSlave = !procPatch.master();
+
+                // Get face cells on this side
                 const unallocLabelList& fc = patches[patchI].faceCells();
 
+                // Get face cells on the other side (communicated during
+                // distributeCells() call)
+                const unallocLabelList& nfc = patchNbrFaceCells_[patchI];
+
                 // Get neighbour cellToProc addressing across the interface
                 const labelList& curNbrPtc = patchNbrCellToProc_[patchI];
 
+                // Note: in order to avoid upper triangular ordering errors for
+                // neighbours after the reconstruction, we need to make sure
+                // that the master processor inserts possibly multiple faces for
+                // a single owner cell in the correct order. For each owner
+                // cell with multiple faces, we need to collect neighbour cell
+                // indices on the other side and possibly swap the order of
+                // adding faces. The same "swapping" of insertion order needs to
+                // happend on the slave side, but now we sort on the remote
+                // (master) data and not on the local (slave) data as we do for
+                // master processor. VV, 16/Feb/2019.
+
+                // A map containing a list of patch faces and neighbour (owner)
+                // cells for master (slave) processor, for each owner cell (for
+                // master, neighbour cell for slave) (key). In the pair, first
+                // entry is the list of patch faces and the second one is the
+                // list of neigbhour (for master proc or owner for slave proc)
+                // cells. We'll sort the list of neighbouring cells and use the
+                // indices to add the patch faces in the correct order later
+                // on. VV, 16/Feb/2019.
+                Map<Pair<dynamicLabelList> > faceCellToNbrFaceCells
+                (
+                    // Reasonable size estimate (minimum = 10)
+                    max(10, fc.size()/10)
+                );
+
                 forAll (fc, patchFaceI)
                 {
+                    // Get cell on this side: make a copy for swapping if this
+                    // is slave processor later on
+                    label ownCellI = fc[patchFaceI];
+
                     // Local owner proc is looked up using faceCells
-                    const label ownerProc = cellToProc_[fc[patchFaceI]];
+                    const label ownerProc = cellToProc_[ownCellI];
 
                     // Neighbour proc is looked up directly
                     const label neighbourProc = curNbrPtc[patchFaceI];
 
                     // If the owner and neighbour processor index is the same,
-                    // the face remains in the processor patch
+                    // the face remains in the processor patch.
                     // In load balancing, it will be re-merged on reconstruction
-                    // HJ, 23/Apr/2018
+                    // HJ, 23/Apr/2018.
                     if (ownerProc == neighbourProc)
                     {
-                        // Add the face
+                        // Get the face cell on the other side: make a copy for
+                        // swapping if this is slave processor
+                        label neiCellI = nfc[patchFaceI];
+
+                        if (isSlave)
+                        {
+                            // This is slave, need to swap owner and neighbour
+                            const label tmpCellI = ownCellI;
+                            ownCellI = neiCellI;
+                            neiCellI = tmpCellI;
+                        }
+
+                        // Add into the map
+                        if (faceCellToNbrFaceCells.found(ownCellI))
+                        {
+                            // This owner cell has already been visited
+                            Pair<dynamicLabelList>& tll =
+                                faceCellToNbrFaceCells[ownCellI];
+
+                            // Insert patch face index into the first list
+                            tll.first().append(patchFaceI);
+
+                            // Insert neighbour into the second list
+                            tll.second().append(neiCellI);
+                        }
+                        else
+                        {
+                            // This owner cell has not been visited
+
+                            // Create dynamicList containing the patch face
+                            dynamicLabelList pfl(4);
+                            pfl.append(patchFaceI);
+
+                            // Create dynamicList containing the neighbour
+                            dynamicLabelList nl(4);
+                            nl.append(neiCellI);
+
+                            // Insert the pair of lists into the map
+                            faceCellToNbrFaceCells.insert
+                            (
+                                ownCellI,
+                                Pair<dynamicLabelList>(pfl, nl)
+                            );
+                        }
+                    } // End if owner and neighbour going to some processor
+                } // End for all patch faces
+
+                // Get sorted table of contents from the map to make sure that
+                // we insert faces in the correct order (on both sides, since
+                // for the slave processor, the map actually contains list of
+                // patch faces and owner cells (on my, slave side) for all
+                // neighbour cells (on the other, master side)
+                const labelList sortedOwn = faceCellToNbrFaceCells.sortedToc();
+
+                // Loop through ordered owner cells with faces on processor
+                // patch
+                forAll (sortedOwn, i)
+                {
+                    // Get current owner cell
+                    const label& ownCellI = sortedOwn[i];
+
+                    // The cell has not been handled yet and removed from
+                    // the map. Get the pair of lists
+                    Pair<dynamicLabelList>& patchFacesAndNeighbours =
+                        faceCellToNbrFaceCells[ownCellI];
+
+                    // Get the list of patch faces and list of neighbours
+                    dynamicLabelList& pfl = patchFacesAndNeighbours.first();
+                    dynamicLabelList& nl = patchFacesAndNeighbours.second();
+
+                    // Create sorted neighbour list.
+                    // Notes:
+                    // 1. Transfer the contents of the dynamic list,
+                    // 2. Sorted on construction
+                    SortableList<label> sortedNbrs(nl.xfer());
+
+                    // Get sorted indices for correct patch face insertion
+                    // order
+                    const labelList& sortedIndices = sortedNbrs.indices();
+
+                    // Loop through sorted indices
+                    forAll (sortedIndices, j)
+                    {
+                        // Get the sorted index into the patch face
+                        const label& sI = sortedIndices[j];
+
+                        // Get patch face index
+                        const label& patchFaceI = pfl[sI];
+
+                        // Get owner processor. Note: by definition/filtering
+                        // above, owner and neighbour proc are the same. In
+                        // order to avoid using cellToProc list for slave
+                        // processor (where ownCellI is actually found on the
+                        // other side), use curNbrPtc with patchFace indes
+                        const label& ownerProc = curNbrPtc[patchFaceI];
+
+                        // Add the patch face into the list in the correct
+                        // order for owner processor. Note: map contains
+                        // only faces where owner proc = neighbour proc, so
+                        // no need to double check
                         procFaceList[ownerProc].append(patchStart + patchFaceI);
 
                         // Increment the number of faces for this patch
@@ -556,8 +604,7 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
             }
             else if (isA<cyclicPolyPatch>(patches[patchI]))
             {
-                // Cyclic patch special treatment
-
+                // Cyclic patch requires special treatment
                 const polyPatch& cPatch = patches[patchI];
 
                 const label cycOffset = cPatch.size()/2;
@@ -578,11 +625,14 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
 
                 forAll (firstFaceCells, patchFaceI)
                 {
-                    if
-                    (
-                        cellToProc_[firstFaceCells[patchFaceI]]
-                     != cellToProc_[secondFaceCells[patchFaceI]]
-                    )
+                    // Get owner and neighbour processor labels
+                    const label& ownerProc =
+                        cellToProc_[firstFaceCells[patchFaceI]];
+                    const label& neighbourProc =
+                        cellToProc_[secondFaceCells[patchFaceI]];
+
+                    // Check whether ownerProc and neighbourProc are different
+                    if (ownerProc != neighbourProc)
                     {
                         // This face becomes an inter-processor boundary face
                         // inter - processor patch face found. Go through
@@ -590,139 +640,52 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
                         // processor and try to find this inter-processor
                         // patch.
 
-                        cyclicParallel_ = true;
+                        // Insert inter-processor data for ownerProc and return
+                        // whether the neighbour was already present in the list
+                        const bool ownerInterProcFound =
+                            addInterProcessorBoundaryData
+                            (
+                                ownerProc,     // Processor to append to
+                                neighbourProc, // Processor to append
+                                patchStart + patchFaceI, // Face index to append
 
-                        label ownerProc =
-                            cellToProc_[firstFaceCells[patchFaceI]];
+                                interProcBoundaries,
+                                interProcBFaces
+                            );
 
-                        label neighbourProc =
-                            cellToProc_[secondFaceCells[patchFaceI]];
+                        // Insert inter-processor data for neighbourProc and
+                        // return whether the neighbour was already present in
+                        // the list
+                        const bool neighbourInterProcFound =
+                            addInterProcessorBoundaryData
+                            (
+                                neighbourProc, // Processor to append to
+                                ownerProc,     // Processor to append
+                                patchStart + patchFaceI, // Face index to append
 
-                        SLList<label>::iterator curInterProcBdrsOwnIter =
-                            interProcBoundaries[ownerProc].begin();
+                                interProcBoundaries,
+                                interProcBFaces
+                            );
 
-                        SLList<SLList<label> >::iterator
-                            curInterProcBFacesOwnIter =
-                            interProcBFaces[ownerProc].begin();
-
-                        bool interProcBouFound = false;
-
-                        // WARNING: Synchronous SLList iterators
-
-                        for
-                        (
-                            ;
-                            curInterProcBdrsOwnIter !=
-                            interProcBoundaries[ownerProc].end()
-                         && curInterProcBFacesOwnIter !=
-                            interProcBFaces[ownerProc].end();
-                            ++curInterProcBdrsOwnIter,
-                            ++curInterProcBFacesOwnIter
-                        )
+                        if (ownerInterProcFound && !neighbourInterProcFound)
                         {
-                            if (curInterProcBdrsOwnIter() == neighbourProc)
-                            {
-                                // the inter - processor boundary exists.
-                                // Add the face
-                                interProcBouFound = true;
-
-                                curInterProcBFacesOwnIter().append
-                                    (patchStart + patchFaceI);
-
-                                SLList<label>::iterator curInterProcBdrsNeiIter
-                                   = interProcBoundaries[neighbourProc].begin();
-
-                                SLList<SLList<label> >::iterator
-                                    curInterProcBFacesNeiIter =
-                                    interProcBFaces[neighbourProc].begin();
-
-                                bool neighbourFound = false;
-
-                                // WARNING: Synchronous SLList iterators
-
-                                for
-                                (
-                                    ;
-                                    curInterProcBdrsNeiIter
-                                   != interProcBoundaries[neighbourProc].end()
-                                 && curInterProcBFacesNeiIter
-                                   != interProcBFaces[neighbourProc].end();
-                                    ++curInterProcBdrsNeiIter,
-                                    ++curInterProcBFacesNeiIter
-                                )
-                                {
-                                    if (curInterProcBdrsNeiIter() == ownerProc)
-                                    {
-                                        // boundary found. Add the face
-                                        neighbourFound = true;
-
-                                        curInterProcBFacesNeiIter()
-                                            .append
-                                            (
-                                                patchStart
-                                              + cycOffset
-                                              + patchFaceI
-                                            );
-                                    }
-
-                                    if (neighbourFound) break;
-                                }
-
-                                if (interProcBouFound && !neighbourFound)
-                                {
-                                    FatalErrorIn
-                                    (
-                                        "domainDecomposition::decomposeMesh()"
-                                    )   << "Inconsistency in inter-processor "
-                                        << "boundary lists for processors "
-                                        << ownerProc << " and "
-                                        << neighbourProc
-                                        << " in cyclic boundary matching"
-                                        << abort(FatalError);
-                                }
-                            }
-
-                            if (interProcBouFound) break;
-                        }
-
-                        if (!interProcBouFound)
-                        {
-                            // inter - processor boundaries do not exist
-                            // and need to be created
-
-                            // set the new addressing information
-
-                            // owner
-                            interProcBoundaries[ownerProc]
-                                .append(neighbourProc);
-                            interProcBFaces[ownerProc]
-                                .append(SLList<label>(patchStart + patchFaceI));
-
-                            // neighbour
-                            interProcBoundaries[neighbourProc]
-                                .append(ownerProc);
-                            interProcBFaces[neighbourProc]
-                                .append
-                                (
-                                    SLList<label>
-                                    (
-                                        patchStart
-                                      + cycOffset
-                                      + patchFaceI
-                                    )
-                                );
+                            FatalErrorIn
+                            (
+                                "domainDecomposition::decomposeMesh()"
+                            )   << "Inconsistency in inter-processor "
+                                << "boundary lists for processors "
+                                << ownerProc << " and "
+                                << neighbourProc
+                                << " in cyclic boundary matching"
+                                << abort(FatalError);
                         }
                     }
                     else
                     {
-                        // This cyclic face remains on the processor
-                        label ownerProc =
-                            cellToProc_[firstFaceCells[patchFaceI]];
-
-                        // add the face
+                        // Add the face
                         procFaceList[ownerProc].append(patchStart + patchFaceI);
 
-                        // increment the number of faces for this patch
+                        // Increment the number of faces for this patch
                         procPatchSize_[ownerProc][patchI]++;
 
                         // Note: I cannot add the other side of the cyclic
@@ -744,7 +707,7 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
                     )
                     {
                         // This cyclic face remains on the processor
-                        label ownerProc =
+                        const label ownerProc =
                             cellToProc_[firstFaceCells[patchFaceI]];
 
                         // Add the second face
@@ -766,10 +729,10 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
                 {
                     const label curProc = cellToProc_[fc[patchFaceI]];
 
-                    // add the face
+                    // Add the face
                     procFaceList[curProc].append(patchStart + patchFaceI);
 
-                    // increment the number of faces for this patch
+                    // Increment the number of faces for this patch
                     procPatchSize_[curProc][patchI]++;
                 }
             }
@@ -781,7 +744,11 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
 
         if (decompositionDict_.found("globalFaceZones"))
         {
-            wordList fzNames(decompositionDict_.lookup("globalFaceZones"));
+            // Get faze zone names and face zones
+            const wordList fzNames
+            (
+                decompositionDict_.lookup("globalFaceZones")
+            );
 
             const faceZoneMesh& fz = mesh_.faceZones();
 
@@ -797,8 +764,11 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
                         << exit(FatalError);
                 }
 
-                Info<< "Preserving global face zone " << fzNames[nameI]
-                    << endl;
+                if (debug)
+                {
+                    Pout<< "Preserving global face zone " << fzNames[nameI]
+                        << endl;
+                }
 
                 const faceZone& curFz = fz[zoneID];
 
@@ -807,14 +777,13 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
                 // otherwise, add the face to all processor face lists
                 forAll (curFz, faceI)
                 {
-                    const label curFaceID = curFz[faceI];
+                    const label& curFaceID = curFz[faceI];
 
                     if (curFaceID < owner.size())
                     {
                         // This is an active mesh face, and it already belongs
                         // to one CPU.  Find out which and add it to the others
-
-                        const label curProc = cellToProc_[owner[curFaceID]];
+                        const label& curProc = cellToProc_[owner[curFaceID]];
 
                         forAll (procZoneFaceList, procI)
                         {
@@ -836,10 +805,8 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
             }
         }
 
-
         // Convert linked lists into normal lists
         // Add inter-processor boundaries and remember start indices
-
         forAll (procFaceList, procI)
         {
             // Get internal and regular boundary processor faces
@@ -863,9 +830,9 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
             for
             (
                 SLList<SLList<label> >::const_iterator curInterProcBFacesIter =
-                    interProcBFaces[procI].begin();
-                curInterProcBFacesIter != interProcBFaces[procI].end();
-                ++curInterProcBFacesIter
+                    interProcBFaces[procI].cbegin();
+                curInterProcBFacesIter != interProcBFaces[procI].cend();
+              ++curInterProcBFacesIter
             )
             {
                 nFacesOnProcessor += curInterProcBFacesIter().size();
@@ -890,9 +857,9 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
             for
             (
                 SLList<label>::const_iterator curProcFacesIter =
-                    curProcFaces.begin();
-                curProcFacesIter != curProcFaces.end();
-                ++curProcFacesIter
+                    curProcFaces.cbegin();
+                curProcFacesIter != curProcFaces.cend();
+              ++curProcFacesIter
             )
             {
                 curProcFaceAddressing[nFaces] = curProcFacesIter() + 1;
@@ -928,9 +895,12 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
             for
             (
                 ;
+
                 curInterProcBdrsIter != interProcBoundaries[procI].end()
              && curInterProcBFacesIter != interProcBFaces[procI].end();
-                ++curInterProcBdrsIter, ++curInterProcBFacesIter
+
+              ++curInterProcBdrsIter,
+              ++curInterProcBFacesIter
             )
             {
                 curProcNeighbourProcessors[nProcPatches] =
@@ -951,7 +921,7 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
                     SLList<label>::iterator curFacesIter =
                         curInterProcBFacesIter().begin();
                     curFacesIter != curInterProcBFacesIter().end();
-                    ++curFacesIter
+                  ++curFacesIter
                 )
                 {
                     // Add the face
@@ -968,7 +938,7 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
                         curProcFaceAddressing[nFaces] = -(curFacesIter() + 1);
                     }
 
-                    // increment the size
+                    // Increment the size
                     curSize++;
 
                     nFaces++;
@@ -986,9 +956,9 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
             for
             (
                 SLList<label>::const_iterator curProcZoneFacesIter =
-                    curProcZoneFaces.begin();
-                curProcZoneFacesIter != curProcZoneFaces.end();
-                ++curProcZoneFacesIter
+                    curProcZoneFaces.cbegin();
+                curProcZoneFacesIter != curProcZoneFaces.cend();
+              ++curProcZoneFacesIter
             )
             {
                 curProcFaceAddressing[nFaces] = curProcZoneFacesIter() + 1;
@@ -997,7 +967,11 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
         } // End for all processors
     } // End of memory management
 
-    Info << "\nCalculating processor boundary addressing" << endl;
+    if (debug)
+    {
+        Pout << "\nCalculating processor boundary addressing" << endl;
+    }
+
     // For every patch of processor boundary, find the index of the original
     // patch. Mis-alignment is caused by the fact that patches with zero size
     // are omitted. For processor patches, set index to -1.
@@ -1066,7 +1040,11 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
         curBoundaryAddressing.setSize(nPatches);
     }
 
-    Info << "\nDistributing points to processors" << endl;
+    if (debug)
+    {
+        Pout << "\nDistributing points to processors" << endl;
+    }
+
     // For every processor, loop through the list of faces for the processor.
     // For every face, loop through the list of points and mark the point as
     // used for the processor. Collect the list of used points for the
diff --git a/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/distributeCells.C b/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/distributeCells.C
index edbb3309b..78297eb6b 100644
--- a/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/distributeCells.C
+++ b/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/distributeCells.C
@@ -184,6 +184,7 @@ void Foam::domainDecomposition::distributeCells()
     {
         const fvBoundaryMesh& patches = mesh_.boundary();
 
+        // Send cellToProc data
         forAll (patches, patchI)
         {
             if (isA<processorFvPatch>(patches[patchI]))
@@ -212,7 +213,43 @@ void Foam::domainDecomposition::distributeCells()
                     cpPatch.internalFieldTransfer
                     (
                         Pstream::blocking,
-                        cellToProc_
+                        cellToProc_ // Dummy argument
+                    );
+            }
+        }
+
+        // Send face cells for correct ordering of faces in parallel load
+        // balancing when certain processor faces become internal faces on a
+        // given processor
+        forAll (patches, patchI)
+        {
+            if (isA<processorFvPatch>(patches[patchI]))
+            // if (patches[patchI].coupled())
+            {
+                const lduInterface& cpPatch =
+                    refCast<const lduInterface>(patches[patchI]);
+
+                cpPatch.initTransfer
+                (
+                    Pstream::blocking,
+                    patches[patchI].faceCells()
+                );
+            }
+        }
+
+        forAll (patches, patchI)
+        {
+            if (isA<processorFvPatch>(patches[patchI]))
+            // if (patches[patchI].coupled())
+            {
+                const lduInterface& cpPatch =
+                    refCast<const lduInterface>(patches[patchI]);
+
+                patchNbrFaceCells_[patchI] =
+                    cpPatch.transfer
+                    (
+                        Pstream::blocking,
+                        patches[patchI].faceCells() // Dummy argument
                     );
             }
         }
diff --git a/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/domainDecomposition.C b/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/domainDecomposition.C
index 3d9069908..cdb9ea84a 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,8 +54,10 @@ Foam::domainDecomposition::domainDecomposition
     nProcs_(readInt(decompositionDict_.lookup("numberOfSubdomains"))),
     distributed_(false),
     gfIndex_(mesh_),
+    gpIndex_(mesh_),
     cellToProc_(mesh_.nCells()),
     patchNbrCellToProc_(mesh_.boundaryMesh().size()),
+    patchNbrFaceCells_(mesh_.boundaryMesh().size()),
     procPointAddressing_(nProcs_),
     procFaceAddressing_(nProcs_),
     nInternalProcFaces_(nProcs_),
@@ -122,8 +125,8 @@ Foam::autoPtr<Foam::fvMesh> Foam::domainDecomposition::processorMesh
     forAll (curFaceLabels, faceI)
     {
         // Mark the original face as used
-        // Remember to decrement the index by one (turning index)
-        // HJ, 5/Dec/2001
+        // Remember to decrement the index by one (turning index) (see
+        // procFaceAddressing_ data member comment). HJ, 5/Dec/2001
         label curF = mag(curFaceLabels[faceI]) - 1;
 
         faceLookup[curF] = faceI;
@@ -174,8 +177,8 @@ Foam::autoPtr<Foam::fvMesh> Foam::domainDecomposition::processorMesh
     // Note: loop over owner, not all faces: sizes are different
     forAll (procOwner, faceI)
     {
-        // Remember to decrement the index by one (turning index)
-        // HJ, 28/Mar/2009
+        // Remember to decrement the index by one (turning index) (see
+        // procFaceAddressing_ data member comment). HJ, 5/Dec/2001
         label curF = mag(curFaceLabels[faceI]) - 1;
 
         if (curFaceLabels[faceI] >= 0)
@@ -193,8 +196,8 @@ Foam::autoPtr<Foam::fvMesh> Foam::domainDecomposition::processorMesh
     // Note: loop over neighbour, not all faces: sizes are different
     forAll (procNeighbour, faceI)
     {
-        // Remember to decrement the index by one (turning index)
-        // HJ, 28/Mar/2009
+        // Remember to decrement the index by one (turning index) (see
+        // procFaceAddressing_ data member comment). HJ, 5/Dec/2001
         label curF = mag(curFaceLabels[faceI]) - 1;
 
         if (curFaceLabels[faceI] >= 0)
@@ -285,6 +288,8 @@ Foam::autoPtr<Foam::fvMesh> Foam::domainDecomposition::processorMesh
 
             forAll (patchGlobalIndex, fI)
             {
+                // Remember to decrement the index by one (turning index) (see
+                // procFaceAddressing_ data member comment). HJ, 5/Dec/2001
                 patchGlobalIndex[fI] =
                     globalIndex[mag(curFaceLabels[curPatchStart + fI]) - 1];
             }
@@ -503,6 +508,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..7fa550ac9 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_;
 
@@ -79,6 +83,11 @@ class domainDecomposition
         //  Data is used when running decomposition in parallel
         labelListList patchNbrCellToProc_;
 
+        //- Face cells for neighbour cell for each processor boundary
+        //  Data is used for correct upper triangular (neighbour list) ordering
+        //  when running decomposition in parallel
+        labelListList patchNbrFaceCells_;
+
         //- Labels of points for each processor
         labelListList procPointAddressing_;
 
@@ -87,9 +96,9 @@ class domainDecomposition
         //  Only the processor boundary faces are affected: if the sign of the
         //  index is negative, the processor face is the reverse of the
         //  original face. In order to do this properly, all face
-        //   indices will be incremented by 1 and the decremented as
-        //  necessary t avoid the problem of face number zero having no
-        //  sign.  HJ, 5/Dec/2001
+        //  indices will be incremented by 1 and then decremented as
+        //  necessary to avoid the problem of face number zero having no
+        //  sign. HJ, 5/Dec/2001
         labelListList procFaceAddressing_;
 
         //- Number of internal faces for each processor
@@ -133,12 +142,28 @@ class domainDecomposition
         //- Create cell-to processor list using domain decomposition tools
         void distributeCells();
 
+        //- Helper function for decomposeMesh. Given two lists of inter processor
+        //  boundaries and inter processor boundary faces, add neighbouring procJ
+        //  to list of the procI. Also adds faceI into the procI list and return
+        //  whether the procJ is already present in the list of visited
+        //  processors (needed for error handling on cyclic patches)
+        bool addInterProcessorBoundaryData
+        (
+            const label& procI,
+            const label& procJ,
+            const label faceI,
+
+            List<SLList<label> >& interProcBoundaries,
+            List<SLList<SLList<label> > >& interProcBFaces
+        ) const;
+
 
 public:
 
     // Declare name of the class and its debug switch
     ClassName("domainDecomposition");
 
+
     // Constructors
 
         //- Construct from mesh and dictionary
@@ -205,6 +230,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 e099977f8..df1afae56 100644
--- a/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesRebuild.C
+++ b/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesRebuild.C
@@ -437,15 +437,13 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
         }
     }
 
-    // Prepare patch reconstruction
-    HashTable<label, word> patchNameLookup
-    (
-        meshes_[firstValidMesh()].boundaryMesh().size()
-    );
-    DynamicList<word> patchTypeLookup
-    (
-        meshes_[firstValidMesh()].boundaryMesh().size()
-    );
+    // Get first boundary mesh size
+    const label firstBndryMeshSize =
+        meshes_[firstValidMesh()].boundaryMesh().size();
+
+    // Prepare patch reconstruction. Note: default key type is word in HashTable
+    HashTable<label> patchNameLookup(firstBndryMeshSize);
+    DynamicList<word> patchTypeLookup(firstBndryMeshSize);
 
     label nReconPatches = 0;
 
@@ -459,35 +457,37 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
             {
                 if (!isA<processorPolyPatch>(procPatches[patchI]))
                 {
-                    const word& patchName = procPatches[patchI].name();
+                    // Regular patch, get the patch and the name
+                    const polyPatch& patch = procPatches[patchI];
+                    const word& patchName = patch.name();
 
-                    // Regular patch.  Try to find it in a list
-                    if (!patchNameLookup.found(patchName))
+                    // Insert the pair into the hash table
+                    if (patchNameLookup.insert(patchName, nReconPatches))
                     {
-                        // Patch not found.  Add it
-                        patchNameLookup.insert(patchName, nReconPatches);
-                        patchTypeLookup.append(procPatches[patchI].type());
-                        nReconPatches++;
+                        // This is the first time we're inserting this patch,
+                        // add its type into the list and increment the counter
+                        patchTypeLookup.append(patch.type());
+                        ++nReconPatches;
                     }
+                    // else already present in the table
                 }
             }
         }
     }
 
-    // Fill in patch names and types
+    // Fill in patch names
     wordList reconPatchNames(patchNameLookup.size());
-    wordList reconPatchTypes(patchNameLookup.size());
 
-    wordList patchNameToc = patchNameLookup.toc();
-
-    forAll (patchNameToc, pnI)
+    // Loop through all the patch names and collect names and types
+    forAllConstIter(HashTable<label>, patchNameLookup, iter)
     {
-        const label pnIndex = patchNameLookup.find(patchNameToc[pnI])();
-
-        reconPatchNames[pnIndex] = patchNameToc[pnI];
-        reconPatchTypes[pnIndex] = patchTypeLookup[pnIndex];
+        reconPatchNames[iter()] = iter.key();
     }
 
+    // Transfer the contents of the patchTypeLookup dynamic list into the
+    // ordinary list. Note: the ordering remains the same
+    const wordList reconPatchTypes(patchTypeLookup.xfer());
+
     // Prepare point, face and patch reconstruction
     label nReconPoints = 0;
     label nReconFaces = 0;
@@ -498,30 +498,31 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
     {
         if (meshes_.set(procI))
         {
-            // Count total number of points and faces
-            nReconPoints += meshes_[procI].nPoints();
-            nReconFaces += meshes_[procI].allFaces().size();
-            nReconCells += meshes_[procI].nCells();
+            // Get current mesh
+            const fvMesh& procMesh = meshes_[procI];
 
-            const polyBoundaryMesh& procPatches = meshes_[procI].boundaryMesh();
+            // Count total number of points and faces
+            nReconPoints += procMesh.nPoints();
+            nReconFaces += procMesh.allFaces().size();
+            nReconCells += procMesh.nCells();
+
+            const polyBoundaryMesh& procPatches = procMesh.boundaryMesh();
 
             forAll (procPatches, patchI)
             {
                 if (!isA<processorPolyPatch>(procPatches[patchI]))
                 {
+                    // Get current patch
+                    const polyPatch& patch = procPatches[patchI];
+
                     // Find processor patch index in reconstructed boundary
-                    const label pnIndex = patchNameLookup.find
-                    (
-                        procPatches[patchI].name()
-                    )();
+                    const label pnIndex = patchNameLookup.find(patch.name())();
 
                     // Check patch name and type
                     if
                     (
-                        procPatches[patchI].name()
-                     != reconPatchNames[pnIndex]
-                     || procPatches[patchI].type()
-                     != reconPatchTypes[pnIndex]
+                        patch.name() != reconPatchNames[pnIndex]
+                     || patch.type() != reconPatchTypes[pnIndex]
                     )
                     {
                         FatalErrorIn
@@ -530,13 +531,12 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
                             "reconstructMesh(const Time& db)"
                         )   << "Patch name and type does not match "
                             << "across processors for patch "
-                            << procPatches[patchI].name() << " type: "
-                            << procPatches[patchI].type()
+                            << patch.name() << " type: " << patch.type()
                             << abort(FatalError);
                     }
 
                     // Record number of faces in patch
-                    reconPatchSizes[pnIndex] += procPatches[patchI].size();
+                    reconPatchSizes[pnIndex] += patch.size();
                 }
             }
         }
@@ -553,10 +553,11 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
 
     forAll (reconPatchFaces, patchI)
     {
+        // Set size of reconstructed patch faces
         reconPatchFaces[patchI].setSize(reconPatchSizes[patchI]);
 
-        reconPatchOwner[patchI].setSize(reconPatchSizes[patchI]);
-        reconPatchOwner[patchI] = -1;
+        // Set size of reconstructed patch face owners with default value of -1
+        reconPatchOwner[patchI].setSize(reconPatchSizes[patchI], -1);
     }
 
     // Size the mapping arrays
@@ -652,10 +653,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);
@@ -784,7 +805,11 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
     // Dump first valid mesh without checking
     {
         const label fvmId = firstValidMesh();
-        Pout<< "Dump mesh " << fvmId << endl;
+
+        if (debug)
+        {
+            Pout<< "Dump mesh " << fvmId << endl;
+        }
 
         cellOffset[fvmId] = 0;
 
@@ -874,14 +899,14 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
             }
         }
 
-        // Get sorting order.  Note make a copy of indices because
+        // Get sorting order. Note: make a copy of indices because
         // sortable list will be deleted
-        labelList procVisitOrder =
+        const labelList procVisitOrder =
             SortableList<label>(nrbProcIndex).indices();
 
         forAll (procVisitOrder, pvoI)
         {
-            const label patchI = procVisitOrder[pvoI];
+            const label& patchI = procVisitOrder[pvoI];
 
             if (isA<processorPolyPatch>(procPatches[patchI]))
             {
@@ -944,15 +969,14 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
 
                 forAll (indices, i)
                 {
-                    // Location in reconstructed patch where the face
-                    // is inserted
-                    const label insertSlot = indices[i];
+                    // Location in reconstructed patch where the face is
+                    // inserted
+                    const label& insertSlot = indices[i];
 
                     // Calculate face index depending on the ordering
                     const label faceI = patchStart + i;
 
                     face newFace = curFaces[faceI];
-
                     inplaceRenumber(ppAddr, newFace);
 
                     // Insert into correct slot
@@ -1015,7 +1039,9 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
             cpAddr[cellI] = cellI; // + cellOffset[firstValidMesh()];
         }
 
-        // Set cell offset for the next mesh
+        // Set cell offset for the next mesh. Note: if the next mesh is not
+        // valid, the cell offset is propagated to the next one in the processor
+        // loop below.
         if (cellOffset.size() > firstValidMesh() + 1)
         {
             cellOffset[firstValidMesh() + 1] =
@@ -1025,14 +1051,25 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
 
 
     // Dump all other meshes, merging the processor boundaries
-
     for (label procI = firstValidMesh() + 1; procI < meshes_.size(); procI++)
     {
-        if (meshes_.set(procI))
+        if (!meshes_.set(procI))
         {
-            Pout<< "Dump mesh " << procI
-                << " cell offset: " << cellOffset[procI]
-                << endl;
+            // Next mesh is not valid: simply propagate cell offset
+            if (cellOffset.size() > procI + 1)
+            {
+                cellOffset[procI + 1] = cellOffset[procI];
+            }
+        }
+        else
+        {
+            // Valid mesh, combine it
+            if (debug)
+            {
+                Pout<< "Dump mesh " << procI
+                    << " cell offset: " << cellOffset[procI]
+                    << endl;
+            }
 
             const polyMesh& curMesh = meshes_[procI];
             const polyBoundaryMesh& procPatches = curMesh.boundaryMesh();
@@ -1170,7 +1207,7 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
                         }
                     } // End of "is neighbour"
                 } // End of "is processor"
-            }
+            } // End for all patches
 
             // Dump unmarked points into the global point list
             label nMergedPoints = 0;
@@ -1274,7 +1311,7 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
 
             // Get sorting order.  Note make a copy of indices because
             // sortable list will be deleted
-            labelList procVisitOrder =
+            const labelList procVisitOrder =
                 SortableList<label>(nrbProcIndex).indices();
 
             forAll (procVisitOrder, pvoI)
@@ -1347,13 +1384,14 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
                             faceI++
                         )
                         {
-                            label faceInPatch = faceI - patchStart;
+                            const label faceInPatch = faceI - patchStart;
 
                             // Calculate master index
-                            label masterIndex = masterProcPatch.start()
-                                + faceInPatch;
+                            const label masterIndex =
+                                masterProcPatch.start() + faceInPatch;
 
-                            label masterFp = masterFaceAddr[masterIndex] - 1;
+                            const label masterFp =
+                                masterFaceAddr[masterIndex] - 1;
 
                             // Record face-cells for the neighbour
                             fpAddr[faceI] = -masterFaceAddr[masterIndex];
@@ -1470,14 +1508,6 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
                 cellOffset[procI + 1] = cellOffset[procI] + meshes_[procI].nCells();
             }
         }
-        else
-        {
-            // No valid mesh.  Propagate cell offset
-            if (cellOffset.size() > procI + 1)
-            {
-                cellOffset[procI + 1] = cellOffset[procI];
-            }
-        }
     }
 
     // Resize the lists
@@ -1519,12 +1549,24 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
 
     // Mesh assembly completed
 
-    Info<< "Global mesh size (final): " << nl
-        << "    nPoints = " << reconPoints.size() << nl
-        << "    nFaces = " << reconFaces.size() << nl
-        << "    nCells = " << nReconCells << nl
-        << "    nPatches = " << reconPatchSizes.size() << nl
-        << "    nPatchFaces = " << reconPatchSizes << endl;
+    if (!Pstream::parRun())
+    {
+        Info<< "Global mesh size (final): " << nl
+            << "    nPoints = " << reconPoints.size() << nl
+            << "    nFaces = " << reconFaces.size() << nl
+            << "    nCells = " << nReconCells << nl
+            << "    nPatches = " << reconPatchSizes.size() << nl
+            << "    nPatchFaces = " << reconPatchSizes << endl;
+    }
+    else if (debug)
+    {
+        Pout<< "Local mesh size (final): " << nl
+            << "    nPoints = " << reconPoints.size() << nl
+            << "    nFaces = " << reconFaces.size() << nl
+            << "    nCells = " << nReconCells << nl
+            << "    nPatches = " << reconPatchSizes.size() << nl
+            << "    nPatchFaces = " << reconPatchSizes << endl;
+    }
 
     // Renumber the face-processor addressing list for all pieces
     // now that the number of internal faces is known
diff --git a/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesReconstructor.C b/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesReconstructor.C
index 18ace98bb..289794ab5 100644
--- a/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesReconstructor.C
+++ b/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesReconstructor.C
@@ -25,6 +25,11 @@ License
 
 #include "processorMeshesReconstructor.H"
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(Foam::processorMeshesReconstructor, 0);
+
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void Foam::processorMeshesReconstructor::readMeshes(PtrList<Time>& databases)
@@ -71,6 +76,7 @@ Foam::processorMeshesReconstructor::processorMeshesReconstructor
 :
     meshName_(meshName),
     meshes_(),
+    globalPointIndex_(),
     pointProcAddressing_(),
     faceProcAddressing_(),
     cellProcAddressing_(),
@@ -86,6 +92,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..bab934041 100644
--- a/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesReconstructor.H
+++ b/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesReconstructor.H
@@ -29,8 +29,6 @@ Description
     matching processor boundaries and build a single combined mesh by
     matching processor patches to each other.
 
-    In order to
-
 Author
     Hrvoje Jasak, Wikki Ltd.  All rights reserved.
 
@@ -72,6 +70,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_;
 
@@ -121,6 +122,10 @@ class processorMeshesReconstructor
 
 public:
 
+    //- Declare name of the class and its debug switch
+    ClassName("processorMeshesReconstructor");
+
+
     // Constructors
 
         //- Construct given name.  Set meshes later
@@ -161,6 +166,12 @@ public:
             return meshes_;
         }
 
+        //- Return to global point index
+        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/dynamicMesh/directTopoChange/boundaryMesh/boundaryPatch.H b/src/dynamicMesh/dynamicMesh/directTopoChange/boundaryMesh/boundaryPatch.H
index 34e1106f2..3a31d0377 100644
--- a/src/dynamicMesh/dynamicMesh/directTopoChange/boundaryMesh/boundaryPatch.H
+++ b/src/dynamicMesh/dynamicMesh/directTopoChange/boundaryMesh/boundaryPatch.H
@@ -44,6 +44,12 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+class boundaryPatch;
+
+Ostream& operator<<(Ostream&, const boundaryPatch&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class boundaryPatch Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/dynamicMesh/dynamicMesh/directTopoChange/meshCut/directions/directionInfo/directionInfo.H b/src/dynamicMesh/dynamicMesh/directTopoChange/meshCut/directions/directionInfo/directionInfo.H
index 36a79cc36..1dd0df734 100644
--- a/src/dynamicMesh/dynamicMesh/directTopoChange/meshCut/directions/directionInfo/directionInfo.H
+++ b/src/dynamicMesh/dynamicMesh/directTopoChange/meshCut/directions/directionInfo/directionInfo.H
@@ -62,6 +62,8 @@ SourceFiles
 
 namespace Foam
 {
+
+// Forward declaration of classes
 class polyPatch;
 class polyMesh;
 class primitiveMesh;
@@ -69,6 +71,13 @@ class edge;
 class face;
 class polyMesh;
 
+// Forward declaration of friend functions and operators
+class directionInfo;
+
+Ostream& operator<<(Ostream&, const directionInfo&);
+Istream& operator>>(Istream&, directionInfo&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class directionInfo Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/dynamicMesh/dynamicMesh/directTopoChange/meshCut/refineCell/refineCell.H b/src/dynamicMesh/dynamicMesh/directTopoChange/meshCut/refineCell/refineCell.H
index ccfa90941..efc2574ab 100644
--- a/src/dynamicMesh/dynamicMesh/directTopoChange/meshCut/refineCell/refineCell.H
+++ b/src/dynamicMesh/dynamicMesh/directTopoChange/meshCut/refineCell/refineCell.H
@@ -43,9 +43,16 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class Istream;
 class Ostream;
 
+// Forward declaration of friend functions and operators
+class refineCell;
+
+Ostream& operator<<(Ostream&, const refineCell&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class refineCell Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/dynamicMesh/dynamicMesh/directTopoChange/meshCut/wallLayerCells/wallNormalInfo/wallNormalInfo.H b/src/dynamicMesh/dynamicMesh/directTopoChange/meshCut/wallLayerCells/wallNormalInfo/wallNormalInfo.H
index 72a9acd28..c55737a52 100644
--- a/src/dynamicMesh/dynamicMesh/directTopoChange/meshCut/wallLayerCells/wallNormalInfo/wallNormalInfo.H
+++ b/src/dynamicMesh/dynamicMesh/directTopoChange/meshCut/wallLayerCells/wallNormalInfo/wallNormalInfo.H
@@ -51,6 +51,13 @@ namespace Foam
 class polyPatch;
 class polyMesh;
 
+// Forward declaration of friend functions and operators
+class wallNormalInfo;
+
+Ostream& operator<<(Ostream&, const wallNormalInfo&);
+Istream& operator>>(Istream&, wallNormalInfo&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class wallNormalInfo Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyMeshModifier/polyMeshModifier.H b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyMeshModifier/polyMeshModifier.H
index 244eb570c..480f0dcad 100644
--- a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyMeshModifier/polyMeshModifier.H
+++ b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyMeshModifier/polyMeshModifier.H
@@ -55,6 +55,12 @@ class polyTopoChanger;
 class polyTopoChange;
 class mapPolyMesh;
 
+// Forward declaration of friend functions and operators
+class polyMeshModifier;
+
+Ostream& operator<<(Ostream&, const polyMeshModifier&);
+
+
 /*---------------------------------------------------------------------------*\
                        Class polyMeshModifier Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/polyhedralRefinement/polyhedralRefinement.C b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/polyhedralRefinement/polyhedralRefinement.C
index d3036301c..0d5f2db57 100644
--- a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/polyhedralRefinement/polyhedralRefinement.C
+++ b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/polyhedralRefinement/polyhedralRefinement.C
@@ -31,6 +31,7 @@ Author
 #include "faceSet.H"
 #include "pointSet.H"
 #include "addToRunTimeSelectionTable.H"
+#include "meshTools.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -2042,7 +2043,16 @@ void Foam::polyhedralRefinement::setCellsToRefine
     // layers
     for (label i = 0; i < nRefinementBufferLayers_; ++i)
     {
-        extendMarkedCellsAcrossFaces(refineCell);
+        meshTools::extendMarkedCellsAcrossFaces(mesh_, refineCell);
+    }
+
+    // Remove all cells that would exceed the maximum refinement level
+    forAll (refineCell, cellI)
+    {
+        if (refineCell[cellI] && (cellLevel_[cellI] + 1 > maxRefinementLevel_))
+        {
+            refineCell[cellI] = false;
+        }
     }
 
     // Remove all cells that would exceed the maximum refinement level
@@ -2086,6 +2096,12 @@ void Foam::polyhedralRefinement::setCellsToRefine
         ++nIters;
         nTotalAddCells += nAddCells;
 
+	if (debug)
+        {
+            Info<< "Added " << nAddCells << " cells in iteration "
+                << nIters << nl;
+        }
+
     } while (nAddCells > 0);
 
     Info<< "Added " << nTotalAddCells // nTotalAddCells already reduced
@@ -2205,7 +2221,7 @@ void Foam::polyhedralRefinement::setSplitPointsToUnrefine
 
     // Note: if there is no dynamic load balancing, points at the boundary
     // cannot be split points by definition. However, in dynamic load balancing
-    // runs, it is possible that a split point end on processor boundary, in
+    // runs, it is possible that a split point ends on processor boundary, in
     // which case we will simply avoid (actually delay) unrefining until this
     // becomes internal point again. VV, 4/Jul/2018.
     const label nInternalFaces = mesh_.nInternalFaces();
@@ -2213,7 +2229,7 @@ void Foam::polyhedralRefinement::setSplitPointsToUnrefine
 
     for (label faceI = nInternalFaces; faceI < nFaces; ++faceI)
     {
-        // Get the face and make sure that the points are unarked
+        // Get the face and make sure that the points are unmarked
         const face& f = meshFaces[faceI];
 
         forAll (f, fpI)
@@ -2260,7 +2276,7 @@ void Foam::polyhedralRefinement::setSplitPointsToUnrefine
     // unrefinement buffer layers
     for (label i = 0; i < nUnrefinementBufferLayers_; ++i)
     {
-        extendMarkedCellsAcrossPoints(protectedCell);
+        meshTools::extendMarkedCellsAcrossPoints(mesh_, protectedCell);
     }
 
     // Loop through all cells and if the cell should be protected, protect all
@@ -2368,6 +2384,12 @@ void Foam::polyhedralRefinement::setSplitPointsToUnrefine
         ++nIters;
         nTotalRemCells += nRemCells;
 
+        if (debug)
+        {
+            Info<< "Removed " << nRemCells << " cells in iteration "
+                << nIters << nl;
+        }
+
     } while (nRemCells > 0);
 
     Info<< "Removed " << nTotalRemCells // nTotalRemCells already reduced
diff --git a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/prismatic2DRefinement/prismatic2DRefinement.C b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/prismatic2DRefinement/prismatic2DRefinement.C
index 203263431..7fa88a836 100644
--- a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/prismatic2DRefinement/prismatic2DRefinement.C
+++ b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/prismatic2DRefinement/prismatic2DRefinement.C
@@ -33,6 +33,7 @@ Author
 #include "emptyPolyPatch.H"
 #include "wedgePolyPatch.H"
 #include "addToRunTimeSelectionTable.H"
+#include "meshTools.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -2601,7 +2602,16 @@ void Foam::prismatic2DRefinement::setCellsToRefine
     // layers
     for (label i = 0; i < nRefinementBufferLayers_; ++i)
     {
-        extendMarkedCellsAcrossFaces(refineCell);
+        meshTools::extendMarkedCellsAcrossFaces(mesh_, refineCell);
+    }
+
+    // Remove all cells that exceed the maximum refinement level
+    forAll (refineCell, cellI)
+    {
+        if (refineCell[cellI] && (cellLevel_[cellI] + 1 > maxRefinementLevel_))
+        {
+            refineCell[cellI] = false;
+        }
     }
 
     // Remove all cells that exceed the maximum refinement level
@@ -2695,7 +2705,7 @@ void Foam::prismatic2DRefinement::setSplitPointsToUnrefine
     // 1. Has pointLevel_ > 0 (obviously),
     // 2. A point that has the same pointLevel_ as ALL of the points of its
     //    edges. In other words, for each point, we will look through all the
-    //    edges of the point. For each edges, we will visit both points and
+    //    edges of the point. For each edge, we will visit both points and
     //    check point levels. All point levels must be the same for this point
     //    candidate to be a split point. This is quite useful since there is no
     //    need to store the refinement history
@@ -2779,11 +2789,14 @@ void Foam::prismatic2DRefinement::setSplitPointsToUnrefine
 
         if (isA<processorPolyPatch>(patch))
         {
+            // Get patch start
+            const label startIndex = patch.start();
+
             // Loop through all the faces
             forAll (patch, i)
             {
                 // Get global face index and face
-                const label faceI = patch.size() + i;
+                const label faceI = startIndex + i;
                 const face& f = meshFaces[faceI];
 
                 // Make sure that we don't split around point at all points of
@@ -2835,7 +2848,7 @@ void Foam::prismatic2DRefinement::setSplitPointsToUnrefine
     // unrefinement buffer layers
     for (label i = 0; i < nUnrefinementBufferLayers_; ++i)
     {
-        extendMarkedCellsAcrossPoints(protectedCell);
+        meshTools::extendMarkedCellsAcrossPoints(mesh_, protectedCell);
     }
 
     // Loop through all cells and if the cell should be protected, protect all
diff --git a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/refinement/refinement.C b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/refinement/refinement.C
index 9e156aa87..eaf3e075e 100644
--- a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/refinement/refinement.C
+++ b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/refinement/refinement.C
@@ -66,130 +66,6 @@ void Foam::refinement::setInstance(const fileName& inst) const
 }
 
 
-void Foam::refinement::extendMarkedCellsAcrossFaces
-(
-    boolList& markedCell
-) const
-{
-    // Mark all faces for all marked cells
-    const label nFaces = mesh_.nFaces();
-    boolList markedFace(nFaces, false);
-
-    // Get mesh cells
-    const cellList& meshCells = mesh_.cells();
-
-    // Loop through all cells
-    forAll (markedCell, cellI)
-    {
-        if (markedCell[cellI])
-        {
-            // This cell is marked, get its faces
-            const cell& cFaces = meshCells[cellI];
-
-            forAll (cFaces, i)
-            {
-                markedFace[cFaces[i]] = true;
-            }
-        }
-    }
-
-    // Snyc the face list across processor boundaries
-    syncTools::syncFaceList(mesh_, markedFace, orEqOp<bool>(), false);
-
-    // Get necessary mesh data
-    const label nInternalFaces = mesh_.nInternalFaces();
-    const labelList& owner = mesh_.faceOwner();
-    const labelList& neighbour = mesh_.faceNeighbour();
-
-    // Internal faces
-    for (label faceI = 0; faceI < nInternalFaces; ++faceI)
-    {
-        if (markedFace[faceI])
-        {
-            // Face is marked, mark both owner and neighbour if the maximum
-            // refinement level is not exceeded
-            const label& own = owner[faceI];
-            const label& nei = neighbour[faceI];
-
-            // Mark owner and neighbour cells
-            markedCell[own] = true;
-            markedCell[nei] = true;
-        }
-    }
-
-    // Boundary faces
-    for (label faceI = nInternalFaces; faceI < nFaces; ++faceI)
-    {
-        if (markedFace[faceI])
-        {
-            // Face is marked, mark owner if the maximum refinement level is not
-            // exceeded
-            const label& own = owner[faceI];
-
-            // Mark owner
-            markedCell[own] = true;
-        }
-    }
-}
-
-
-void Foam::refinement::extendMarkedCellsAcrossPoints
-(
-    boolList& markedCell
-) const
-{
-    // Mark all points for all marked cells
-    const label nPoints = mesh_.nPoints();
-    boolList markedPoint(nPoints, false);
-
-    // Get cell points
-    const labelListList& meshCellPoints = mesh_.cellPoints();
-
-    // Loop through all cells
-    forAll (markedCell, cellI)
-    {
-        if (markedCell[cellI])
-        {
-            // This cell is marked, get its points
-            const labelList& cPoints = meshCellPoints[cellI];
-
-            forAll (cPoints, i)
-            {
-                markedPoint[cPoints[i]] = true;
-            }
-        }
-    }
-
-    // Snyc point list across processor boundaries
-    syncTools::syncPointList
-    (
-        mesh_,
-        markedPoint,
-        orEqOp<bool>(),
-        true, // Default value
-        true  // Apply separation for parallel cyclics
-    );
-
-    // Get point cells
-    const labelListList& meshPointCells = mesh_.pointCells();
-
-    // Loop through all points
-    forAll (markedPoint, pointI)
-    {
-        if (markedPoint[pointI])
-        {
-            // This point is marked, mark all of its cells
-            const labelList& pCells = meshPointCells[pointI];
-
-            forAll (pCells, i)
-            {
-                markedCell[pCells[i]] = true;
-            }
-        }
-    }
-}
-
-
 Foam::label Foam::refinement::addFace
 (
     polyTopoChange& ref,
@@ -797,7 +673,14 @@ Foam::label Foam::refinement::faceConsistentRefinement
         // Note: we are using more stringent 1:1 consistency across coupled
         // boundaries in order to simplify handling of edge based consistency
         // checks for parallel runs
-        if (neiLevel[i] > curOwnLevel)
+        // Bugfix related to PLB: Check whether owner is already marked for
+        // refinement. Will allow 2:1 consistency across certain processor faces
+        // where we have a new processor boundary. VV, 23/Jan/2019.
+        if
+        (
+            (neiLevel[i] > curOwnLevel)
+         && !cellsToRefine[own]
+        )
         {
             // Neighbour level is higher than owner level, owner must be
             // marked for refinement
@@ -949,9 +832,11 @@ Foam::label Foam::refinement::faceConsistentUnrefinement
                     << "Try increasing nUnrefinementBufferLayers. "
                     << abort(FatalError);
             }
-
-            cellsToUnrefine[own] = false;
-            ++nRemCells;
+            else
+            {
+                cellsToUnrefine[own] = false;
+                ++nRemCells;
+            }
         }
         else if (neiLevel < (ownLevel - 1))
         {
@@ -977,9 +862,11 @@ Foam::label Foam::refinement::faceConsistentUnrefinement
                     << "Try increasing nUnrefinementBufferLayers. "
                     << abort(FatalError);
             }
-
-            cellsToUnrefine[nei] = false;
-            ++nRemCells;
+            else
+            {
+                cellsToUnrefine[nei] = false;
+                ++nRemCells;
+            }
         }
     }
 
@@ -1011,7 +898,7 @@ Foam::label Foam::refinement::faceConsistentUnrefinement
 
         // Note: we are using more stringent 1:1 consistency across coupled
         // boundaries in order to simplify handling of edge based consistency
-        // checkes for parallel runs
+        // checks for parallel runs
         if (curOwnLevel < neiLevel[i])
         {
             // Owner level is smaller than neighbour level, we must not
@@ -1046,7 +933,11 @@ Foam::label Foam::refinement::faceConsistentUnrefinement
                         << "This is probably because the refinement and "
                         << "unrefinement regions are very close." << nl
                         << "Try increasing nUnrefinementBufferLayers. "
-                        << abort(FatalError);
+                        << nl
+                        << "Another possibility is that you are running "
+                        << "with Dynamic Load Balancing, in which case "
+                        << "this should be fine."
+                        << endl;
                 }
             }
             else
@@ -1123,26 +1014,31 @@ Foam::label Foam::refinement::edgeConsistentUnrefinement
                     // unrefinement
                     if (!cellsToUnrefine[cellI])
                     {
-                        FatalErrorIn
-                        (
-                            "label refinement::"
-                            "edgeConsistentUnrefinement"
-                            "(boolList& cellsToUnrefine)"
-                        )   << "Cell not marked for unrefinement, indicating a"
-                            << " previous unnoticed problem with unrefinement."
-                            << nl
-                            << "cellI: " << cellI << ", cellJ: " << cellJ
-                            << nl
-                            << "Level of cellI: " << cellILevel
-                            << ", level of cellJ: " << cellJLevel << nl
-                            << "This is probably because the refinement and "
-                            << "unrefinement regions are very close." << nl
-                            << "Try increasing nUnrefinementBufferLayers. "
-                            << abort(FatalError);
+                        if (debug)
+                        {
+                            WarningIn
+                            (
+                                "label refinement::"
+                                "edgeConsistentUnrefinement"
+                                "(boolList& cellsToUnrefine)"
+                            )   << "Cell not marked for unrefinement, indicating a"
+                                << " previous unnoticed problem with unrefinement."
+                                << nl
+                                << "cellI: " << cellI << ", cellJ: " << cellJ
+                                << nl
+                                << "Level of cellI: " << cellILevel
+                                << ", level of cellJ: " << cellJLevel << nl
+                                << "This is probably because the refinement and "
+                                << "unrefinement regions are very close." << nl
+                                << "Try increasing nUnrefinementBufferLayers. "
+                                << endl;
+                        }
+                    }
+                    else
+                    {
+                        cellsToUnrefine[cellI] = false;
+                        ++nRemCells;
                     }
-
-                    cellsToUnrefine[cellI] = false;
-                    ++nRemCells;
                 }
                 else if (cellJLevel < cellILevel - 1)
                 {
@@ -1153,26 +1049,31 @@ Foam::label Foam::refinement::edgeConsistentUnrefinement
                     // unrefinement
                     if (!cellsToUnrefine[cellJ])
                     {
-                        FatalErrorIn
-                        (
-                            "label refinement::"
-                            "edgeConsistentUnrefinement"
-                            "(boolList& cellsToUnrefine)"
-                        )   << "Cell not marked for unrefinement, indicating a"
-                            << " previous unnoticed problem with unrefinement."
-                            << nl
-                            << "cellI: " << cellI << ", cellJ: " << cellJ
-                            << nl
-                            << "Level of cellI: " << cellILevel
-                            << ", level of cellJ: " << cellJLevel << nl
-                            << "This is probably because the refinement and "
-                            << "unrefinement regions are very close." << nl
-                            << "Try increasing nUnrefinementBufferLayers. "
-                            << abort(FatalError);
+                        if (debug)
+                        {
+                            WarningIn
+                            (
+                                "label refinement::"
+                                "edgeConsistentUnrefinement"
+                                "(boolList& cellsToUnrefine)"
+                            )   << "Cell not marked for unrefinement, indicating a"
+                                << " previous unnoticed problem with unrefinement."
+                                << nl
+                                << "cellI: " << cellI << ", cellJ: " << cellJ
+                                << nl
+                                << "Level of cellI: " << cellILevel
+                                << ", level of cellJ: " << cellJLevel << nl
+                                << "This is probably because the refinement and "
+                                << "unrefinement regions are very close." << nl
+                                << "Try increasing nUnrefinementBufferLayers. "
+                                << endl;
+                        }
+                    }
+                    else
+                    {
+                        cellsToUnrefine[cellJ] = false;
+                        ++nRemCells;
                     }
-
-                    cellsToUnrefine[cellJ] = false;
-                    ++nRemCells;
                 }
             }
         }
@@ -1442,7 +1343,7 @@ void Foam::refinement::setRefinement(polyTopoChange& ref) const
         mesh_,
         pointLevel_,
         maxEqOp<label>(),
-        label(0), // Null value
+        0,   // Null value
         true // Apply separation for parallel cyclics
     );
 
@@ -1599,10 +1500,21 @@ void Foam::refinement::updateMesh(const mapPolyMesh& map)
         }
     }
 
-    // Note: new point level is going to be synced at processor boundaries just
-    // before the next step in setRefinement. Need to investigate why the sync
-    // is not done properly if I put it here. Something is not updated yet.
-    // VV, 31/Jan/2018.
+    // Sync the new point level. Note: here, we assume that the call to
+    // updateMesh happened after polyBoundaryMesh::updateMesh where the
+    // processor data is fully rebuilt. If this is not the case, the point
+    // level will remain unsynced and will cause all kinds of trouble that
+    // are extremely difficult to spot. See the change in
+    // polyTopoChanger::changeMesh order of calling polyMesh::updateMesh and
+    // polyTopoChanger::update. VV, 19/Feb/2019.
+    syncTools::syncPointList
+    (
+        mesh_,
+        newPointLevel,
+        maxEqOp<label>(),
+        0,   // Null value
+        true // Apply separation for parallel cyclics
+    );
 
     // Transfer the new point level into the data member
     pointLevel_.transfer(newPointLevel);
diff --git a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/refinement/refinement.H b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/refinement/refinement.H
index 74b79cf92..f361ca707 100644
--- a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/refinement/refinement.H
+++ b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/refinement/refinement.H
@@ -159,12 +159,6 @@ protected:
             //- Set file instance for cellLevel_ and pointLevel_
             void setInstance(const fileName& inst) const;
 
-            //- Extend marked cells across faces
-            void extendMarkedCellsAcrossFaces(boolList& markedCell) const;
-
-            //- Extend marked cells across points
-            void extendMarkedCellsAcrossPoints(boolList& markedCell) const;
-
 
         // Local topology modification functions (operate on cells/faces)
 
@@ -189,8 +183,8 @@ protected:
                 const label nei
             ) const;
 
-            //- Modifies existing faceI for either new owner/neighbour or new face
-            //  points. Reverses if nessecary
+            //- Modifies existing faceI for either new owner/neighbour or new
+            //  face points. Reverses if nessecary
             void modifyFace
             (
                 polyTopoChange& ref,
diff --git a/src/dynamicMesh/dynamicMesh/polyTopoChange/polyTopoChanger/polyTopoChanger.H b/src/dynamicMesh/dynamicMesh/polyTopoChange/polyTopoChanger/polyTopoChanger.H
index 7926a4343..ba55603bd 100644
--- a/src/dynamicMesh/dynamicMesh/polyTopoChange/polyTopoChanger/polyTopoChanger.H
+++ b/src/dynamicMesh/dynamicMesh/polyTopoChange/polyTopoChanger/polyTopoChanger.H
@@ -44,12 +44,17 @@ SourceFiles
 namespace Foam
 {
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
+// Forward declaration of classes
 class polyMesh;
 class mapPolyMesh;
 class polyBoundaryMesh;
 
+// Forward declaration of friend functions and operators
+class polyTopoChanger;
+
+Ostream& operator<<(Ostream&, const polyTopoChanger&);
+
+
 /*---------------------------------------------------------------------------*\
                       Class polyTopoChanger Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/dynamicMesh/dynamicMesh/polyTopoChange/polyTopoChanger/polyTopoChangerChangeMesh.C b/src/dynamicMesh/dynamicMesh/polyTopoChange/polyTopoChanger/polyTopoChangerChangeMesh.C
index 325845f55..48573620e 100644
--- a/src/dynamicMesh/dynamicMesh/polyTopoChange/polyTopoChanger/polyTopoChangerChangeMesh.C
+++ b/src/dynamicMesh/dynamicMesh/polyTopoChange/polyTopoChanger/polyTopoChangerChangeMesh.C
@@ -2490,10 +2490,17 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChanger::changeMesh()
             topoChangeRequest()()
         );
 
-        update(topoChangeMap());
-
+        // Mesh data needs to be updated before the update of polyMeshModifiers
+        // because polyMeshModifiers might need all the new polyMesh data (see
+        // below for further comments). VV, 19/Feb/2019
         mesh_.updateMesh(topoChangeMap());
 
+        // Bugfix: call to polyTopoChanger::update must happen after
+        // polyMesh::updateMesh where all the relevant mesh bits for parallel
+        // comms are updated. First noticed when the syncying of pointLevel in
+        // refinement::updateMesh was not syncying properly. VV, 19/Feb/2019
+        update(topoChangeMap());
+
         // Increment the morph index
         morphIndex_++;
 
diff --git a/src/dynamicMesh/dynamicMesh/refinementData/refinementData.H b/src/dynamicMesh/dynamicMesh/refinementData/refinementData.H
index 8f25a31af..ff34e916b 100644
--- a/src/dynamicMesh/dynamicMesh/refinementData/refinementData.H
+++ b/src/dynamicMesh/dynamicMesh/refinementData/refinementData.H
@@ -45,9 +45,17 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class polyPatch;
 class polyMesh;
 
+// Forward declaration of friend functions and operators
+class refinementData;
+
+Ostream& operator<<(Ostream&, const refinementData&);
+Istream& operator>>(Istream&, refinementData&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class refinementData Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/dynamicMesh/dynamicMesh/refinementData/refinementDistanceData.H b/src/dynamicMesh/dynamicMesh/refinementData/refinementDistanceData.H
index 1145fa55a..da2167396 100644
--- a/src/dynamicMesh/dynamicMesh/refinementData/refinementDistanceData.H
+++ b/src/dynamicMesh/dynamicMesh/refinementData/refinementDistanceData.H
@@ -45,9 +45,17 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class polyPatch;
 class polyMesh;
 
+// Forward declaration of friend functions and operators
+class refinementDistanceData;
+
+Ostream& operator<<(Ostream&, const refinementDistanceData&);
+Istream& operator>>(Istream&, refinementDistanceData&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class refinementDistanceData Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/dynamicMesh/dynamicMesh/refinementData/refinementHistory.H b/src/dynamicMesh/dynamicMesh/refinementData/refinementHistory.H
index aaf5a6d45..a1a25bb73 100644
--- a/src/dynamicMesh/dynamicMesh/refinementData/refinementHistory.H
+++ b/src/dynamicMesh/dynamicMesh/refinementData/refinementHistory.H
@@ -87,6 +87,14 @@ namespace Foam
 // Forward declaration of classes
 class mapPolyMesh;
 class mapDistributePolyMesh;
+class polyMesh;
+
+// Forward declaration of friend functions and operators
+class refinementHistory;
+
+Ostream& operator<<(Ostream&, const refinementHistory&);
+Istream& operator>>(Istream&, refinementHistory&);
+
 
 /*---------------------------------------------------------------------------*\
                            Class refinementHistory Declaration
@@ -360,6 +368,10 @@ public:
 };
 
 
+Istream& operator>>(Istream&, refinementHistory::splitCell8&);
+Ostream& operator<<(Ostream&, const refinementHistory::splitCell8&);
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C b/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C
index cec584173..5a5eacf23 100644
--- a/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C
+++ b/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C
@@ -271,15 +271,11 @@ bool Foam::dynamicPolyRefinementFvMesh::update()
             (
                 refinementSelectionPtr_->refinementCellCandidates()()
             );
-
-            if (debug)
-            {
-                Pout<< "Selected " << refCandidates.size()
-                    << " refinement candidates."
-                    << endl;
-            }
+            Pout<< "Selected " << refCandidates.size()
+                << " refinement candidates."
+                << endl;
         }
-        else if (debug)
+        else
         {
             Pout<< "Skipping refinement for this time-step..." << endl;
         }
@@ -301,15 +297,11 @@ bool Foam::dynamicPolyRefinementFvMesh::update()
             (
                 refinementSelectionPtr_->unrefinementPointCandidates()()
             );
-
-            if (debug)
-            {
-                Pout<< "Selected " << unrefCandidates.size()
-                    << " unrefinement candidates."
-                    << endl;
-            }
+            Pout<< "Selected " << unrefCandidates.size()
+                << " unrefinement candidates."
+                << endl;
         }
-        else if (debug)
+        else
         {
             Pout<< "Skipping unrefinement for this time-step..." << endl;
         }
@@ -352,13 +344,13 @@ bool Foam::dynamicPolyRefinementFvMesh::update()
         // some topo changes
         if (sizeCellMap)
         {
-            Info<< "Successfully performed polyhedral refinement. "
+            Pout<< "Successfully performed polyhedral refinement. "
                 << "Changed from " << nOldCells << " to " << sizeCellMap
                 << " cells." << endl;
         }
         else
         {
-            Info<< "Refinement/unrefinement not performed in this time step "
+            Pout<< "Refinement/unrefinement not performed in this time step "
                 << "since no cells were selected." << endl;
         }
 
@@ -370,9 +362,7 @@ bool Foam::dynamicPolyRefinementFvMesh::update()
         // per time step
         curTimeIndex_ = time().timeIndex();
     }
-
-    Info<< "No refinement/unrefinement" << endl;
-
+    Pout<< "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 4be8fed12..1d0ebb281 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"
 
 #include "IFstream.H"
 #include "OFstream.H"
@@ -107,7 +108,11 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
     // Now that each processor has filled in its own part, combine the data
     Pstream::gatherList(migratedCells);
     Pstream::scatterList(migratedCells);
-    Info<< "Migrated cells per processor: " << migratedCells << endl;
+
+    if (debug)
+    {
+        Info<< "Migrated cells per processor: " << migratedCells << endl;
+    }
 
     // Reading through second index now tells how many cells will arrive
     // from which processor
@@ -141,10 +146,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();
@@ -245,6 +254,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
@@ -256,8 +288,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();
@@ -289,6 +322,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);
@@ -345,6 +379,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
@@ -413,7 +450,9 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
                     );
                 }
 
-                //HJ Insert clouds missing.  HJ, 12/Oct/2018
+                // HR, 18.11.2018. Distribution of clouds is trivial and is
+                // treated in Cloud<ParticleType>::split, which is called in
+                // the constructor of CloudDistribute.
             }
         }
     }
@@ -450,9 +489,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
                         ),
@@ -461,6 +500,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
@@ -524,34 +567,40 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
         }
     }
 
-    forAll (procMeshes, procI)
+    if (debug)
     {
-        if (procMeshes.set(procI))
+        forAll (procMeshes, procI)
         {
-            Pout<< "procMesh " << procI
-                << " points " << procMeshes[procI].nPoints()
-                << " faces: " << procMeshes[procI].nFaces()
-                << " internal: " << procMeshes[procI].nInternalFaces()
-                << " cells: " << procMeshes[procI].nCells()
-                << " patches: " << procMeshes[procI].boundary().size()
-                << endl;
+            if (procMeshes.set(procI))
+            {
+                Pout<< "procMesh " << procI
+                    << " points " << procMeshes[procI].nPoints()
+                    << " faces: " << procMeshes[procI].nFaces()
+                    << " internal: " << procMeshes[procI].nInternalFaces()
+                    << " cells: " << procMeshes[procI].nCells()
+                    << " patches: " << procMeshes[procI].boundary().size()
+                    << endl;
+            }
         }
     }
 
     // Create the reconstructed mesh
     autoPtr<fvMesh> reconstructedMeshPtr =
-        meshRecon.reconstructMesh(time());
+        meshRecon.reconstructMesh(dummyTime);
     fvMesh& reconMesh = reconstructedMeshPtr();
 
-    Pout<< "Reconstructed mesh stats: "
-        << " nCells: " << reconMesh.nCells()
-        << " nFaces: " << reconMesh.nFaces()
-        << " nIntFaces: " << reconMesh.nInternalFaces()
-        << " polyPatches: "
-        << reconMesh.boundaryMesh().size()
-        << " patches: "
-        << reconMesh.boundary().size()
-        << endl;
+    if (debug)
+    {
+        Pout<< "Reconstructed mesh stats: "
+            << " nCells: " << reconMesh.nCells()
+            << " nFaces: " << reconMesh.nFaces()
+            << " nIntFaces: " << reconMesh.nInternalFaces()
+            << " polyPatches: "
+            << reconMesh.boundaryMesh().size()
+            << " patches: "
+            << reconMesh.boundary().size()
+            << endl;
+    }
 
     // Apply changes to the local mesh:
     // - refactor the boundary to match new patches.  Note: processor
@@ -920,8 +969,14 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
         );
     }
 
-    // Debug: remove?  HJ, 22/Oct/2018
-    // checkMesh(true);
+    // HR 13.12.18: Update the mesh objects
+    meshObjectBase::allUpdateTopology<polyMesh>(*this, meshMap);
+
+    if (debug)
+    {
+        Info<< "Checking reconstructed mesh after load balancing..." << endl;
+        checkMesh(true);
+    }
 
     return true;
 }
diff --git a/src/dynamicMesh/topoChangerFvMesh/topoChangerFvMesh/topoChangerFvMeshLoadBalanceTemplates.C b/src/dynamicMesh/topoChangerFvMesh/topoChangerFvMesh/topoChangerFvMeshLoadBalanceTemplates.C
index 406b2eff3..238d81f97 100644
--- a/src/dynamicMesh/topoChangerFvMesh/topoChangerFvMesh/topoChangerFvMeshLoadBalanceTemplates.C
+++ b/src/dynamicMesh/topoChangerFvMesh/topoChangerFvMesh/topoChangerFvMeshLoadBalanceTemplates.C
@@ -89,7 +89,6 @@ void Foam::topoChangerFvMesh::insertFields
         iter
     )
     {
-
         localFields[fI].set
         (
             Pstream::myProcNo(),
@@ -126,7 +125,7 @@ void Foam::topoChangerFvMesh::receiveFields
     receivedFields.setSize(nScalarFields);
 
     fromProc.readBegin("topoChangerFvMeshReceiveFields");
-    
+
     forAll (receivedFields, fI)
     {
         word fieldName(fromProc);
@@ -181,7 +180,10 @@ void Foam::topoChangerFvMesh::rebuildFields
 
         const PtrList<GeoField>& partFields = receivedFields[fieldI];
 
-        Pout<< "Rebuilding field " << masterField.name() << endl;
+        if (debug)
+        {
+            Pout<< "Rebuilding field " << masterField.name() << endl;
+        }
 
         // Check name match.  Note: there may be holes
         word partName;
@@ -263,15 +265,19 @@ 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()
-                << " new size = "
-                << GeoField::GeoMeshType::size(masterField.mesh())
-                << endl;
+            if (debug)
+            {
+                Pout<< "Resizing internal field: old size = "
+                    << masterField.size()
+                    << " new size = "
+                    << GeoField::GeoMeshType::size(masterField.mesh())
+                    << endl;
+            }
+
             masterField.setSize
             (
                 GeoField::GeoMeshType::size(masterField.mesh())
@@ -280,13 +286,16 @@ void Foam::topoChangerFvMesh::rebuildFields
 
         // Resize boundary (number of patches)
         typename GeoField::GeometricBoundaryField& patchFields =
-                masterField.boundaryField();
+            masterField.boundaryFieldNoStoreOldTimes();
 
         if (patchFields.size() != masterField.mesh().boundary().size())
         {
-            Pout<< "Resizing boundary field: "
+            if (debug)
+            {
+                Pout<< "Resizing boundary field: "
                     << masterField.mesh().boundary().size()
                     << endl;
+            }
 
             patchFields.setSize(masterField.mesh().boundary().size());
         }
@@ -297,9 +306,12 @@ void Foam::topoChangerFvMesh::rebuildFields
             if (meshMap.resetPatchFlag()[patchI])
             {
                 // Create a new constrained patch field
-                Pout<< "Inserting constrained patch field for patch "
-                    << masterField.mesh().boundary()[patchI].name()
-                    << endl;
+                if (debug)
+                {
+                    Pout<< "Inserting constrained patch field for patch "
+                        << masterField.mesh().boundary()[patchI].name()
+                        << endl;
+                }
 
                 patchFields.set
                 (
@@ -326,12 +338,15 @@ void Foam::topoChangerFvMesh::rebuildFields
             )
             {
                 // Resize patch field
-                Pout<< "Resizing patch field for patch "
-                    << masterField.mesh().boundary()[patchI].name()
-                    << " old size: " << patchFields[patchI].size()
-                    << " new size: "
-                    << masterField.mesh().boundary()[patchI].size()
-                    << endl;
+                if (debug)
+                {
+                    Pout<< "Resizing patch field for patch "
+                        << masterField.mesh().boundary()[patchI].name()
+                        << " old size: " << patchFields[patchI].size()
+                        << " new size: "
+                        << masterField.mesh().boundary()[patchI].size()
+                        << endl;
+                }
 
                 // Reset patch field size
                 patchFields[patchI].autoMap
@@ -350,7 +365,27 @@ void Foam::topoChangerFvMesh::rebuildFields
 
         // Increment field counter
         fieldI++;
-        Pout<< "... done" << endl;
+
+        if (debug)
+        {
+            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/edgeMesh/edgeMesh.H b/src/edgeMesh/edgeMesh.H
index f3560ab41..eabb780ae 100644
--- a/src/edgeMesh/edgeMesh.H
+++ b/src/edgeMesh/edgeMesh.H
@@ -45,6 +45,13 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+class edgeMesh;
+
+Ostream& operator<<(Ostream&, const edgeMesh&);
+Istream& operator>>(Istream&, edgeMesh&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class edgeMesh Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/equationReader/equationReader/equationReader.H b/src/equationReader/equationReader/equationReader.H
index a91822dd4..926ca99e9 100644
--- a/src/equationReader/equationReader/equationReader.H
+++ b/src/equationReader/equationReader/equationReader.H
@@ -107,16 +107,18 @@ Author
 namespace Foam
 {
 
-// Forward declaration of friend functions and operators
-
+// Forward declaration of classes
 class equation;
 class equationOperation;
 
+// Forward declaration of friend functions and operators
 // *** Located in equationReaderIO.C ***
 class equationReader;
+
 Istream& operator>>(Istream&, equationReader&);
 Ostream& operator<<(Ostream&, const equationReader&);
 
+
 /*---------------------------------------------------------------------------*\
                         Class equationReader Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/finiteArea/faMesh/faPatches/faPatch/faPatch.H b/src/finiteArea/faMesh/faPatches/faPatch/faPatch.H
index 20b19e59a..8afaa60b9 100644
--- a/src/finiteArea/faMesh/faPatches/faPatch/faPatch.H
+++ b/src/finiteArea/faMesh/faPatches/faPatch/faPatch.H
@@ -50,7 +50,12 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
 class faBoundaryMesh;
+class faPatch;
+
+Ostream& operator<<(Ostream&, const faPatch&);
+
 
 /*---------------------------------------------------------------------------*\
                            Class faPatch Declaration
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/general/findRefCell/findRefCell.C b/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C
index eb879b516..dfa1a002e 100644
--- a/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C
+++ b/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C
@@ -127,7 +127,7 @@ Foam::scalar Foam::getRefCellValue
 )
 {
     scalar refCellValue = (refCelli >= 0 ? field[refCelli] : 0.0);
-    return returnReduce<label>(refCellValue, sumOp<scalar>());
+    return returnReduce<scalar>(refCellValue, sumOp<scalar>());
 }
 
 
diff --git a/src/finiteVolume/cfdTools/general/porousMedia/porousZone.H b/src/finiteVolume/cfdTools/general/porousMedia/porousZone.H
index 36e184cc4..a19520601 100644
--- a/src/finiteVolume/cfdTools/general/porousMedia/porousZone.H
+++ b/src/finiteVolume/cfdTools/general/porousMedia/porousZone.H
@@ -80,8 +80,15 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class fvMesh;
 
+// Forward declaration of friend functions and operators
+class porousZone;
+
+Ostream& operator<<(Ostream&, const porousZone&);
+
+
 /*---------------------------------------------------------------------------*\
                         Class porousZone Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C
index 2c7ddbbf5..57dc90f7e 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C
+++ b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C
@@ -662,8 +662,7 @@ void Foam::solutionControl::calcSteadyConsistentFlux
                 "void solutionControl::calcSteadyConsistentFlux"
                 "\n("
                 "\n    surfaceScalarField& phi,"
-                "\n    const volVectorField& U,"
-                "\n    const volScalarField& rAU"
+                "\n    const volVectorField& U"
                 "\n)"
             )   << "Either aCoeff or faceU is allocated for field " << UName
                 << " while the other is not."
@@ -684,8 +683,7 @@ void Foam::solutionControl::calcSteadyConsistentFlux
                 "void solutionControl::calcSteadyConsistentFlux"
                 "\n("
                 "\n    surfaceScalarField& phi,"
-                "\n    const volVectorField& U,"
-                "\n    const volScalarField& rAU"
+                "\n    const volVectorField& U"
                 "\n)"
             )   << "Index is set, but the aCoeff field is not allocated for "
                 << UName << "." << nl
@@ -701,8 +699,7 @@ void Foam::solutionControl::calcSteadyConsistentFlux
                 "void solutionControl::calcSteadyConsistentFlux"
                 "\n("
                 "\n    surfaceScalarField& phi,"
-                "\n    const volVectorField& U,"
-                "\n    const volScalarField& rAU"
+                "\n    const volVectorField& U"
                 "\n)"
             )   << "Index is set, but the faceU field is allocated for "
                 << UName << "." << nl
@@ -794,8 +791,7 @@ void Foam::solutionControl::calcSteadyConsistentFlux
                 "void solutionControl::calcSteadyConsistentFlux"
                 "\n("
                 "\n    surfaceScalarField& phi,"
-                "\n    const volVectorField& U,"
-                "\n    const volScalarField& rAU"
+                "\n    const volVectorField& U"
                 "\n)"
             )   << "Either aCoeff or faceU is allocated for field " << UName
                 << " while the other is not."
@@ -816,8 +812,7 @@ void Foam::solutionControl::calcSteadyConsistentFlux
                 "void solutionControl::calcSteadyConsistentFlux"
                 "\n("
                 "\n    surfaceScalarField& phi,"
-                "\n    const volVectorField& U,"
-                "\n    const volScalarField& rAU"
+                "\n    const volVectorField& U"
                 "\n)"
             )   << "Index is set, but the aCoeff field is not allocated for "
                 << UName << "." << nl
@@ -833,8 +828,7 @@ void Foam::solutionControl::calcSteadyConsistentFlux
                 "void solutionControl::calcSteadyConsistentFlux"
                 "\n("
                 "\n    surfaceScalarField& phi,"
-                "\n    const volVectorField& U,"
-                "\n    const volScalarField& rAU"
+                "\n    const volVectorField& U"
                 "\n)"
             )   << "Index is set, but the faceU field is allocated for "
                 << UName << "." << nl
@@ -928,11 +922,11 @@ void Foam::solutionControl::calcSteadyMRFConsistentFlux
         {
             FatalErrorIn
             (
-                "void solutionControl::calcSteadyConsistentFlux"
+                "void solutionControl::calcSteadyMRFConsistentFlux"
                 "\n("
                 "\n    surfaceScalarField& phi,"
                 "\n    const volVectorField& U,"
-                "\n    const volScalarField& rAU"
+                "\n    const MRFZones& mrfZones"
                 "\n)"
             )   << "Either aCoeff or faceU is allocated for field " << UName
                 << " while the other is not."
@@ -954,7 +948,7 @@ void Foam::solutionControl::calcSteadyMRFConsistentFlux
                 "\n("
                 "\n    surfaceScalarField& phi,"
                 "\n    const volVectorField& U,"
-                "\n    const volScalarField& rAU"
+                "\n    const MRFZones& mrfZones"
                 "\n)"
             )   << "Index is set, but the aCoeff field is not allocated for "
                 << UName << "." << nl
@@ -971,7 +965,7 @@ void Foam::solutionControl::calcSteadyMRFConsistentFlux
                 "\n("
                 "\n    surfaceScalarField& phi,"
                 "\n    const volVectorField& U,"
-                "\n    const volScalarField& rAU"
+                "\n    const MRFZones& mrfZones"
                 "\n)"
             )   << "Index is set, but the faceU field is allocated for "
                 << UName << "." << nl
@@ -1044,10 +1038,10 @@ void Foam::solutionControl::reconstructTransientVelocity
             "void solutionControl::reconstructTransientVelocity"
             "\n("
             "\n    volVectorField& U,"
+            "\n    const surfaceScalarField& phi,"
             "\n    const volVectorField& ddtUEqn,"
             "\n    const volScalarField& rAU,"
             "\n    const volScalarField& p"
-            "\n    const surfaceScalarField& phi"
             "\n) const"
         )   << "faceU not calculated for field " << UName
             << ". Make sure you have called"
@@ -1110,7 +1104,8 @@ const Foam::surfaceScalarField& Foam::solutionControl::aCoeff
     {
         FatalErrorIn
         (
-            "const surfaceScalarField& solutionControl::aCoeff() const"
+            "const surfaceScalarField& solutionControl::aCoeff"
+            "(const word& UName) const"
         )   << "aCoeff not calculated for field " << UName
             << ". Make sure you have called"
             << " calcTransientConsistentFlux(...) or "
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.H b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.H
index ca9dcb720..0fe9b8af7 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.H
+++ b/src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.H
@@ -104,6 +104,20 @@ protected:
             label corrNonOrtho_;
 
 
+        // Fields necessary for time and under-relaxation consistency
+
+            //- A coeff (A^~) arising from consistency formulation. Note: we can
+            //  have multiple pressure/velocity systems, hence the PtrList
+            mutable PtrList<surfaceScalarField> aCoeffPtrs_;
+
+            //- Face velocity needed for consistent formulation. Note: we can
+            //  have multiple pressure/velocity systems, hence the PtrList
+            mutable PtrList<surfaceVectorField> faceUPtrs_;
+
+            //- Hash Table containing indices of PtrLists for given names
+            mutable HashTable<label> indices_;
+
+
     // Protected Member Functions
 
         //- Read controls from fvSolution dictionary
@@ -194,22 +208,6 @@ protected:
 
 private:
 
-    // Private data
-
-        // Fields necessary for time and under-relaxation consistency
-
-            //- A coeff (A^~) arising from consistency formulation. Note: we can
-            //  have multiple pressure/velocity systems, hence the PtrList
-            mutable PtrList<surfaceScalarField> aCoeffPtrs_;
-
-            //- Face velocity needed for consistent formulation. Note: we can
-            //  have multiple pressure/velocity systems, hence the PtrList
-            mutable PtrList<surfaceVectorField> faceUPtrs_;
-
-            //- Hash Table containing indices of PtrLists for given names
-            mutable HashTable<label> indices_;
-
-
     // Private member functions
 
         //- Disallow default bitwise copy construct
diff --git a/src/finiteVolume/cfdTools/incompressible/CourantNo.H b/src/finiteVolume/cfdTools/incompressible/CourantNo.H
index b1e05b72a..397351773 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
@@ -61,7 +61,11 @@ else
 
     const scalar deltaX = Foam::cbrt(mesh.V()[0]);
 
-    velMag = mag(U[0]);
+    // recover velocity field in a more general way
+    const volVectorField& URef
+	    = mesh.db().lookupObject<const volVectorField>("U");
+
+    velMag = mag(URef[0]);
 
     CoNum = velMag*deltaT/deltaX;
 
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 c112a8102..5967aad1a 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 80420e0bb..c6c143d97 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/algorithms/octree/octree/octreeDataPoint.H b/src/foam/algorithms/octree/octree/octreeDataPoint.H
index d6abf0209..bdf816acf 100644
--- a/src/foam/algorithms/octree/octree/octreeDataPoint.H
+++ b/src/foam/algorithms/octree/octree/octreeDataPoint.H
@@ -51,8 +51,15 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
 template<class Type> class octree;
 
+// Forward declaration of friend functions and operators
+class octreeDataPoint;
+
+Ostream& operator<<(Ostream&, const octreeDataPoint&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class octreeDataPoint Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/foam/algorithms/octree/octree/treeBoundBox.H b/src/foam/algorithms/octree/octree/treeBoundBox.H
index 84794793d..2ddc1a2e0 100644
--- a/src/foam/algorithms/octree/octree/treeBoundBox.H
+++ b/src/foam/algorithms/octree/octree/treeBoundBox.H
@@ -66,8 +66,19 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class Random;
 
+// Forward declaration of friend functions and operators
+class treeBoundBox;
+
+bool operator==(const treeBoundBox&, const treeBoundBox&);
+bool operator!=(const treeBoundBox&, const treeBoundBox&);
+
+Istream& operator>>(Istream& is, treeBoundBox&);
+Ostream& operator<<(Ostream& os, const treeBoundBox&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class treeBoundBox Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/foam/algorithms/octree/octree/treeLeaf.H b/src/foam/algorithms/octree/octree/treeLeaf.H
index 4eb474d31..0eb2626b9 100644
--- a/src/foam/algorithms/octree/octree/treeLeaf.H
+++ b/src/foam/algorithms/octree/octree/treeLeaf.H
@@ -50,13 +50,14 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class treeBoundBox;
 class Ostream;
 
 template<class Type> class octree;
-template<class Type> class treeLeaf;
 
 // Forward declaration of friend functions and operators
+template<class Type> class treeLeaf;
 
 template<class Type> Istream& operator>>(Istream&, treeLeaf<Type>&);
 template<class Type> Ostream& operator<<(Ostream&, const treeLeaf<Type>&);
diff --git a/src/foam/algorithms/octree/octree/treeNode.H b/src/foam/algorithms/octree/octree/treeNode.H
index 73338497b..17c6cde18 100644
--- a/src/foam/algorithms/octree/octree/treeNode.H
+++ b/src/foam/algorithms/octree/octree/treeNode.H
@@ -58,12 +58,12 @@ namespace Foam
 {
 
 // class intersection;
-
 template<class Type> class octree;
 template<class Type> class treeLeaf;
-template<class Type> class treeNode;
+
 
 // Forward declaration of friend functions and operators
+template<class Type> class treeNode;
 
 template<class Type> Istream& operator>>(Istream&, treeNode<Type>&);
 template<class Type> Ostream& operator<<(Ostream&, const treeNode<Type>&);
diff --git a/src/foam/containers/HashTables/HashPtrTable/HashPtrTable.H b/src/foam/containers/HashTables/HashPtrTable/HashPtrTable.H
index 14c5bae88..ecfaa3edd 100644
--- a/src/foam/containers/HashTables/HashPtrTable/HashPtrTable.H
+++ b/src/foam/containers/HashTables/HashPtrTable/HashPtrTable.H
@@ -43,6 +43,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class Istream;
 class Ostream;
 
diff --git a/src/foam/containers/HashTables/HashTable/HashTable.H b/src/foam/containers/HashTables/HashTable/HashTable.H
index 43967edde..6f864e8a3 100644
--- a/src/foam/containers/HashTables/HashTable/HashTable.H
+++ b/src/foam/containers/HashTables/HashTable/HashTable.H
@@ -55,10 +55,11 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declaration of friend functions and operators
-
+// Forward declaration of classes
 template<class T> class List;
 template<class T> class UList;
+
+// Forward declaration of friend functions and operators
 template<class T, class Key, class Hash> class HashTable;
 template<class T, class Key, class Hash> class HashPtrTable;
 
diff --git a/src/foam/containers/HashTables/StaticHashTable/StaticHashTable.H b/src/foam/containers/HashTables/StaticHashTable/StaticHashTable.H
index 1f1f166f6..ebb4f1c3f 100644
--- a/src/foam/containers/HashTables/StaticHashTable/StaticHashTable.H
+++ b/src/foam/containers/HashTables/StaticHashTable/StaticHashTable.H
@@ -53,9 +53,10 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declaration of friend functions and operators
-
+// Forward declaration of classes
 template<class T> class List;
+
+// Forward declaration of friend functions and operators
 template<class T, class Key, class Hash> class StaticHashTable;
 
 template<class T, class Key, class Hash> Istream& operator>>
diff --git a/src/foam/containers/LinkedLists/accessTypes/ILList/ILList.H b/src/foam/containers/LinkedLists/accessTypes/ILList/ILList.H
index 1cdb76dbc..f17cf99d5 100644
--- a/src/foam/containers/LinkedLists/accessTypes/ILList/ILList.H
+++ b/src/foam/containers/LinkedLists/accessTypes/ILList/ILList.H
@@ -43,6 +43,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class Istream;
 class Ostream;
 
diff --git a/src/foam/containers/LinkedLists/accessTypes/LList/LList.H b/src/foam/containers/LinkedLists/accessTypes/LList/LList.H
index 785a025c3..a8ae51865 100644
--- a/src/foam/containers/LinkedLists/accessTypes/LList/LList.H
+++ b/src/foam/containers/LinkedLists/accessTypes/LList/LList.H
@@ -44,6 +44,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class Istream;
 class Ostream;
 
diff --git a/src/foam/containers/LinkedLists/accessTypes/UILList/UILList.H b/src/foam/containers/LinkedLists/accessTypes/UILList/UILList.H
index 1beef3930..7fc266980 100644
--- a/src/foam/containers/LinkedLists/accessTypes/UILList/UILList.H
+++ b/src/foam/containers/LinkedLists/accessTypes/UILList/UILList.H
@@ -44,6 +44,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class Ostream;
 
 // Forward declaration of friend functions and operators
diff --git a/src/foam/containers/Lists/FixedList/FixedList.H b/src/foam/containers/Lists/FixedList/FixedList.H
index 062e906a3..edcc4a6fd 100644
--- a/src/foam/containers/Lists/FixedList/FixedList.H
+++ b/src/foam/containers/Lists/FixedList/FixedList.H
@@ -49,8 +49,12 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declaration of friend functions and operators
+// Forward declaration of classes
+template<class T> class UList;
+template<class T> class SLList;
 
+
+// Forward declaration of friend functions and operators
 template<class T, unsigned Size> class FixedList;
 
 template<class T, unsigned Size>
@@ -59,9 +63,6 @@ Istream& operator>>(Istream&, FixedList<T, Size>&);
 template<class T, unsigned Size>
 Ostream& operator<<(Ostream&, const FixedList<T, Size>&);
 
-template<class T> class UList;
-template<class T> class SLList;
-
 
 /*---------------------------------------------------------------------------*\
                            Class FixedList Declaration
diff --git a/src/foam/coordinateSystems/coordinateSystem.H b/src/foam/coordinateSystems/coordinateSystem.H
index a6f603f21..ea6c459be 100644
--- a/src/foam/coordinateSystems/coordinateSystem.H
+++ b/src/foam/coordinateSystems/coordinateSystem.H
@@ -143,9 +143,16 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declaration of friend functions and operators
+// Forward declaration of classes
 class boundBox;
 
+// Forward declaration of friend functions and operators
+class coordinateSystem;
+
+bool operator!=(const coordinateSystem&, const coordinateSystem&);
+Ostream& operator<<(Ostream&, const coordinateSystem&);
+
+
 /*---------------------------------------------------------------------------*\
                      Class coordinateSystem Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/foam/db/IOobjects/CompactIOList/CompactIOList.H b/src/foam/db/IOobjects/CompactIOList/CompactIOList.H
index 078cf0a24..7077f8c49 100644
--- a/src/foam/db/IOobjects/CompactIOList/CompactIOList.H
+++ b/src/foam/db/IOobjects/CompactIOList/CompactIOList.H
@@ -47,6 +47,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class Istream;
 class Ostream;
 
@@ -64,6 +65,7 @@ template<class T, class BaseType> Ostream& operator<<
     const CompactIOList<T, BaseType>&
 );
 
+
 /*---------------------------------------------------------------------------*\
                            Class CompactIOList Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/foam/db/dictionary/dictionary.H b/src/foam/db/dictionary/dictionary.H
index 2a0b01da7..984f5fec0 100644
--- a/src/foam/db/dictionary/dictionary.H
+++ b/src/foam/db/dictionary/dictionary.H
@@ -65,14 +65,17 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declaration of friend functions and operators
+// Forward declaration of classes
 class regExp;
-class dictionary;
 class SHA1Digest;
 
+// Forward declaration of friend functions and operators
+class dictionary;
+
 Istream& operator>>(Istream&, dictionary&);
 Ostream& operator<<(Ostream&, const dictionary&);
 
+
 /*---------------------------------------------------------------------------*\
                         Class dictionaryName Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/foam/db/dictionary/dictionaryEntry/dictionaryEntry.H b/src/foam/db/dictionary/dictionaryEntry/dictionaryEntry.H
index 4a3fad529..bd0e60df9 100644
--- a/src/foam/db/dictionary/dictionaryEntry/dictionaryEntry.H
+++ b/src/foam/db/dictionary/dictionaryEntry/dictionaryEntry.H
@@ -51,6 +51,12 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+class dictionaryEntry;
+
+Ostream& operator<<(Ostream&, const dictionaryEntry&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class dictionaryEntry Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/foam/dimensionedTypes/dimensionedType/dimensionedType.H b/src/foam/dimensionedTypes/dimensionedType/dimensionedType.H
index 50421ff6f..d4d728981 100644
--- a/src/foam/dimensionedTypes/dimensionedType/dimensionedType.H
+++ b/src/foam/dimensionedTypes/dimensionedType/dimensionedType.H
@@ -45,8 +45,10 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declaration of friend functions and operators
+// Forward declaration of classes
+class dictionary;
 
+// Forward declaration of friend functions and operators
 template<class Type> class dimensioned;
 
 template<class Type>
@@ -55,7 +57,6 @@ Istream& operator>>(Istream&, dimensioned<Type>&);
 template<class Type>
 Ostream& operator<<(Ostream&, const dimensioned<Type>&);
 
-class dictionary;
 
 /*---------------------------------------------------------------------------*\
                            Class dimensioned Declaration
diff --git a/src/foam/fields/CoeffField/diagTensorCoeffField.H b/src/foam/fields/CoeffField/diagTensorCoeffField.H
index d58d8713b..8e1189ce4 100644
--- a/src/foam/fields/CoeffField/diagTensorCoeffField.H
+++ b/src/foam/fields/CoeffField/diagTensorCoeffField.H
@@ -47,6 +47,16 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+
+Ostream& operator<<(Ostream&, const CoeffField<diagTensor>&);
+Ostream& operator<<(Ostream&, const tmp<CoeffField<diagTensor> >&);
+
+
+/*---------------------------------------------------------------------------*\
+              Class symmTensor4thOrderCoeffField Declaration
+\*---------------------------------------------------------------------------*/
+
 template<>
 class CoeffField<diagTensor>
 :
diff --git a/src/foam/fields/CoeffField/scalarCoeffField.H b/src/foam/fields/CoeffField/scalarCoeffField.H
index 060b43994..13357ec73 100644
--- a/src/foam/fields/CoeffField/scalarCoeffField.H
+++ b/src/foam/fields/CoeffField/scalarCoeffField.H
@@ -45,6 +45,16 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+
+Ostream& operator<<(Ostream&, const CoeffField<scalar>&);
+Ostream& operator<<(Ostream&, const tmp<CoeffField<scalar> >&);
+
+
+/*---------------------------------------------------------------------------*\
+                     Class scalarCoeffField Declaration
+\*---------------------------------------------------------------------------*/
+
 template<>
 class CoeffField<scalar>
 :
diff --git a/src/foam/fields/CoeffField/sphericalTensorCoeffField.H b/src/foam/fields/CoeffField/sphericalTensorCoeffField.H
index 773fa89f0..fec8a0f73 100644
--- a/src/foam/fields/CoeffField/sphericalTensorCoeffField.H
+++ b/src/foam/fields/CoeffField/sphericalTensorCoeffField.H
@@ -47,6 +47,16 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+
+Ostream& operator<<(Ostream&, const CoeffField<sphericalTensor>&);
+Ostream& operator<<(Ostream&, const tmp<CoeffField<sphericalTensor> >&);
+
+
+/*---------------------------------------------------------------------------*\
+               Class sphericalTensorCoeffField Declaration
+\*---------------------------------------------------------------------------*/
+
 template<>
 class CoeffField<sphericalTensor>
 :
diff --git a/src/foam/fields/CoeffField/symmTensor4thOrderCoeffField.H b/src/foam/fields/CoeffField/symmTensor4thOrderCoeffField.H
index 58ead4310..e7303216b 100644
--- a/src/foam/fields/CoeffField/symmTensor4thOrderCoeffField.H
+++ b/src/foam/fields/CoeffField/symmTensor4thOrderCoeffField.H
@@ -47,6 +47,16 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+
+Ostream& operator<<(Ostream&, const CoeffField<symmTensor4thOrder>&);
+Ostream& operator<<(Ostream&, const tmp<CoeffField<symmTensor4thOrder> >&);
+
+
+/*---------------------------------------------------------------------------*\
+              Class symmTensor4thOrderCoeffField Declaration
+\*---------------------------------------------------------------------------*/
+
 template<>
 class CoeffField<symmTensor4thOrder>
 :
diff --git a/src/foam/fields/CoeffField/symmTensorCoeffField.H b/src/foam/fields/CoeffField/symmTensorCoeffField.H
index 79f97250c..4a424e246 100644
--- a/src/foam/fields/CoeffField/symmTensorCoeffField.H
+++ b/src/foam/fields/CoeffField/symmTensorCoeffField.H
@@ -47,6 +47,16 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+
+Ostream& operator<<(Ostream&, const CoeffField<symmTensor>&);
+Ostream& operator<<(Ostream&, const tmp<CoeffField<symmTensor> >&);
+
+
+/*---------------------------------------------------------------------------*\
+                    Class symmTensorCoeffField Declaration
+\*---------------------------------------------------------------------------*/
+
 template<>
 class CoeffField<symmTensor>
 :
diff --git a/src/foam/fields/CoeffField/tensorCoeffField.H b/src/foam/fields/CoeffField/tensorCoeffField.H
index 6b837b826..8fc9a7c84 100644
--- a/src/foam/fields/CoeffField/tensorCoeffField.H
+++ b/src/foam/fields/CoeffField/tensorCoeffField.H
@@ -47,6 +47,16 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+
+Ostream& operator<<(Ostream&, const CoeffField<tensor>&);
+Ostream& operator<<(Ostream&, const tmp<CoeffField<tensor> >&);
+
+
+/*---------------------------------------------------------------------------*\
+                       Class tensorCoeffField Declaration
+\*---------------------------------------------------------------------------*/
+
 template<>
 class CoeffField<tensor>
 :
diff --git a/src/foam/fields/GeometricFields/GeometricField/GeometricField.C b/src/foam/fields/GeometricFields/GeometricField/GeometricField.C
index e61a7883f..70bd583c7 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..4ecb6528a 100644
--- a/src/foam/fields/GeometricFields/GeometricField/GeometricField.H
+++ b/src/foam/fields/GeometricFields/GeometricField/GeometricField.H
@@ -51,6 +51,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class dictionary;
 
 // Forward declaration of friend functions and operators
@@ -447,6 +448,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 03fee6e59..7e9842ae0 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/matrices/LUscalarMatrix/procLduInterface.H b/src/foam/matrices/LUscalarMatrix/procLduInterface.H
index 4edec793e..d5470eae7 100644
--- a/src/foam/matrices/LUscalarMatrix/procLduInterface.H
+++ b/src/foam/matrices/LUscalarMatrix/procLduInterface.H
@@ -43,8 +43,15 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class lduInterfaceField;
 
+// Forward declaration of friend functions and operators
+class procLduInterface;
+
+Ostream& operator<<(Ostream&, const procLduInterface&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class procLduInterface Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/foam/matrices/LUscalarMatrix/procLduMatrix.H b/src/foam/matrices/LUscalarMatrix/procLduMatrix.H
index 999e1262f..b39bc6a9b 100644
--- a/src/foam/matrices/LUscalarMatrix/procLduMatrix.H
+++ b/src/foam/matrices/LUscalarMatrix/procLduMatrix.H
@@ -45,6 +45,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class procLduInterface;
 class lduMatrix;
 
diff --git a/src/foam/matrices/crMatrix/crAddressing.H b/src/foam/matrices/crMatrix/crAddressing.H
index 1682a3521..605a22f5d 100644
--- a/src/foam/matrices/crMatrix/crAddressing.H
+++ b/src/foam/matrices/crMatrix/crAddressing.H
@@ -46,6 +46,12 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+class crAddressing;
+
+Ostream& operator<<(Ostream&, const crAddressing&);
+
+
 /*---------------------------------------------------------------------------*\
                         Class crAddressing Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/foam/matrices/crMatrix/crMatrix.H b/src/foam/matrices/crMatrix/crMatrix.H
index a377fc997..5a5ba50a3 100644
--- a/src/foam/matrices/crMatrix/crMatrix.H
+++ b/src/foam/matrices/crMatrix/crMatrix.H
@@ -46,6 +46,12 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+class crMatrix;
+
+Ostream& operator<<(Ostream&, const crMatrix&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class crMatrix Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/foam/matrices/lduMatrix/lduAddressing/lduInterfaces/lduInterface/lduInterface.C b/src/foam/matrices/lduMatrix/lduAddressing/lduInterfaces/lduInterface/lduInterface.C
index ac309dbba..11e870e9a 100644
--- a/src/foam/matrices/lduMatrix/lduAddressing/lduInterfaces/lduInterface/lduInterface.C
+++ b/src/foam/matrices/lduMatrix/lduAddressing/lduInterfaces/lduInterface/lduInterface.C
@@ -59,7 +59,7 @@ Foam::autoPtr<Foam::crMatrix> Foam::lduInterface::prolongationTransfer
     );
 
     // Dummy return to make the compiler happy
-    return autoPtr<crMatrix>(NULL);
+    return autoPtr<crMatrix>(nullptr);
 }
 
 
diff --git a/src/foam/meshes/boundBox/boundBox.H b/src/foam/meshes/boundBox/boundBox.H
index c23256ede..8d2742abb 100644
--- a/src/foam/meshes/boundBox/boundBox.H
+++ b/src/foam/meshes/boundBox/boundBox.H
@@ -40,11 +40,11 @@ namespace Foam
 {
 
 // Forward declaration of friend functions and operators
-
 class boundBox;
 template<class T> class tmp;
 
-Ostream& operator<<(Ostream& os, const boundBox& b);
+Istream& operator>>(Istream& is, boundBox&);
+Ostream& operator<<(Ostream& os, const boundBox&);
 
 
 /*---------------------------------------------------------------------------*\
diff --git a/src/foam/meshes/meshShapes/face/face.H b/src/foam/meshes/meshShapes/face/face.H
index 854eb08ce..b38a6dd6a 100644
--- a/src/foam/meshes/meshShapes/face/face.H
+++ b/src/foam/meshes/meshShapes/face/face.H
@@ -57,18 +57,20 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declaration of friend functions and operators
-
-class face;
+// Forward declaration of classes
 class triFace;
 
 template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
 class DynamicList;
 
+// Forward declaration of friend functions and operators
+class face;
+
 inline bool operator==(const face& a, const face& b);
 inline bool operator!=(const face& a, const face& b);
 inline Istream& operator>>(Istream&, face&);
 
+
 /*---------------------------------------------------------------------------*\
                            Class face Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/foam/meshes/meshShapes/triFace/triFace.H b/src/foam/meshes/meshShapes/triFace/triFace.H
index b90c948d9..f9b73755f 100644
--- a/src/foam/meshes/meshShapes/triFace/triFace.H
+++ b/src/foam/meshes/meshShapes/triFace/triFace.H
@@ -53,11 +53,12 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class face;
 
 // Forward declaration of friend functions and operators
-
 class triFace;
+
 inline bool operator==(const triFace&, const triFace&);
 inline bool operator!=(const triFace&, const triFace&);
 
diff --git a/src/foam/meshes/meshTools/meshTools.C b/src/foam/meshes/meshTools/meshTools.C
index 5745a09c2..b99542297 100644
--- a/src/foam/meshes/meshTools/meshTools.C
+++ b/src/foam/meshes/meshTools/meshTools.C
@@ -27,6 +27,7 @@ License
 #include "polyMesh.H"
 #include "hexMatcher.H"
 #include "faceZone.H"
+#include "syncTools.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
@@ -927,4 +928,128 @@ void Foam::meshTools::setFaceInfo
 }
 
 
+void Foam::meshTools::extendMarkedCellsAcrossFaces
+(
+    const polyMesh& mesh,
+    boolList& markedCell
+)
+{
+    // Mark all faces for all marked cells
+    const label nFaces = mesh.nFaces();
+    boolList markedFace(nFaces, false);
+
+    // Get mesh cells
+    const cellList& meshCells = mesh.cells();
+
+    // Loop through all cells
+    forAll (markedCell, cellI)
+    {
+        if (markedCell[cellI])
+        {
+            // This cell is marked, get its faces
+            const cell& cFaces = meshCells[cellI];
+
+            forAll (cFaces, i)
+            {
+                markedFace[cFaces[i]] = true;
+            }
+        }
+    }
+
+    // Snyc the face list across processor boundaries
+    syncTools::syncFaceList(mesh, markedFace, orEqOp<bool>(), false);
+
+    // Get necessary mesh data
+    const label nInternalFaces = mesh.nInternalFaces();
+    const labelList& owner = mesh.faceOwner();
+    const labelList& neighbour = mesh.faceNeighbour();
+
+    // Internal faces
+    for (label faceI = 0; faceI < nInternalFaces; ++faceI)
+    {
+        if (markedFace[faceI])
+        {
+            // Face is marked, mark both owner and neighbour
+            const label& own = owner[faceI];
+            const label& nei = neighbour[faceI];
+
+            // Mark owner and neighbour cells
+            markedCell[own] = true;
+            markedCell[nei] = true;
+        }
+    }
+
+    // Boundary faces
+    for (label faceI = nInternalFaces; faceI < nFaces; ++faceI)
+    {
+        if (markedFace[faceI])
+        {
+            // Face is marked, mark owner
+            const label& own = owner[faceI];
+
+            // Mark owner
+            markedCell[own] = true;
+        }
+    }
+}
+
+
+void Foam::meshTools::extendMarkedCellsAcrossPoints
+(
+    const polyMesh& mesh,
+    boolList& markedCell
+)
+{
+    // Mark all points for all marked cells
+    const label nPoints = mesh.nPoints();
+    boolList markedPoint(nPoints, false);
+
+    // Get cell points
+    const labelListList& meshCellPoints = mesh.cellPoints();
+
+    // Loop through all cells
+    forAll (markedCell, cellI)
+    {
+        if (markedCell[cellI])
+        {
+            // This cell is marked, get its points
+            const labelList& cPoints = meshCellPoints[cellI];
+
+            forAll (cPoints, i)
+            {
+                markedPoint[cPoints[i]] = true;
+            }
+        }
+    }
+
+    // Snyc point list across processor boundaries
+    syncTools::syncPointList
+    (
+        mesh,
+        markedPoint,
+        orEqOp<bool>(),
+        true, // Default value
+        true  // Apply separation for parallel cyclics
+    );
+
+    // Get point cells
+    const labelListList& meshPointCells = mesh.pointCells();
+
+    // Loop through all points
+    forAll (markedPoint, pointI)
+    {
+        if (markedPoint[pointI])
+        {
+            // This point is marked, mark all of its cells
+            const labelList& pCells = meshPointCells[pointI];
+
+            forAll (pCells, i)
+            {
+                markedCell[pCells[i]] = true;
+            }
+        }
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/src/foam/meshes/meshTools/meshTools.H b/src/foam/meshes/meshTools/meshTools.H
index f977f8fca..10ae04a43 100644
--- a/src/foam/meshes/meshTools/meshTools.H
+++ b/src/foam/meshes/meshTools/meshTools.H
@@ -318,6 +318,25 @@ namespace meshTools
             label& zoneFlip
         );
 
+
+    // Mark-up of mesh bits. Relocated from refinement polyMeshModifier
+
+        //- Extend marked cells across faces given a bool list of already marked
+        //  cells
+        void extendMarkedCellsAcrossFaces
+        (
+            const polyMesh& mesh,
+            boolList& markedCell
+        );
+
+        //- Extend marked cells across points given a bool list of already
+        //  marked cells
+        void extendMarkedCellsAcrossPoints
+        (
+            const polyMesh& mesh,
+            boolList& markedCell
+        );
+
 } // End namespace meshTools
 
 
diff --git a/src/foam/meshes/patchIdentifier/patchIdentifier.H b/src/foam/meshes/patchIdentifier/patchIdentifier.H
index 27d878a69..674a75e9e 100644
--- a/src/foam/meshes/patchIdentifier/patchIdentifier.H
+++ b/src/foam/meshes/patchIdentifier/patchIdentifier.H
@@ -43,11 +43,12 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class dictionary;
 
 // Forward declaration of friend functions and operators
-
 class patchIdentifier;
+
 Ostream& operator<<(Ostream&, const patchIdentifier&);
 
 
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..0b221932f 100644
--- a/src/foam/meshes/polyMesh/globalMeshData/globalProcFaceIndex.C
+++ b/src/foam/meshes/polyMesh/globalMeshData/globalProcFaceIndex.C
@@ -35,7 +35,7 @@ void Foam::globalProcFaceIndex::calcFaceIndex()
 
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
-    // Assing unique face label to all master processor faces
+    // Assign unique face label to all master processor faces
 
     // Count faces and processor faces per processor
     forAll (patches, patchI)
@@ -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..6655e3b53 100644
--- a/src/foam/meshes/polyMesh/globalMeshData/globalProcFaceIndex.H
+++ b/src/foam/meshes/polyMesh/globalMeshData/globalProcFaceIndex.H
@@ -33,7 +33,7 @@ Description
 
     Face offsets counts a number of unique faces on each processor,
     excluding slave processor patch faces, which are given master face index.
-    If needed, global face indes from all faces can be derived from this data
+    If needed, global face index from all faces can be derived from this data
 
     Currently, faces are ordered with internal faces first, followed by patch
     faces in patch order, excluding slave processor patches.
@@ -120,9 +120,9 @@ 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
+        const labelList& globalLabel() const
         {
             return globalLabel_;
         }
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/foam/meshes/polyMesh/zones/ZoneID/ZoneID.H b/src/foam/meshes/polyMesh/zones/ZoneID/ZoneID.H
index 6b0147582..af3e74bac 100644
--- a/src/foam/meshes/polyMesh/zones/ZoneID/ZoneID.H
+++ b/src/foam/meshes/polyMesh/zones/ZoneID/ZoneID.H
@@ -44,6 +44,7 @@ Description
 namespace Foam
 {
 
+// Forward declaration of classes
 template<class ZoneType, class MeshType> class ZoneMesh;
 
 // Forward declaration of friend functions and operators
diff --git a/src/foam/meshes/polyMesh/zones/faceZone/faceZone.H b/src/foam/meshes/polyMesh/zones/faceZone/faceZone.H
index 199e5c497..822a1e9d8 100644
--- a/src/foam/meshes/polyMesh/zones/faceZone/faceZone.H
+++ b/src/foam/meshes/polyMesh/zones/faceZone/faceZone.H
@@ -53,11 +53,12 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class mapPolyMesh;
 
 // Forward declaration of friend functions and operators
-
 class faceZone;
+
 Ostream& operator<<(Ostream&, const faceZone&);
 
 
diff --git a/src/foam/meshes/primitiveShapes/tetrahedron/tetrahedron.H b/src/foam/meshes/primitiveShapes/tetrahedron/tetrahedron.H
index 7213381b3..9232f818a 100644
--- a/src/foam/meshes/primitiveShapes/tetrahedron/tetrahedron.H
+++ b/src/foam/meshes/primitiveShapes/tetrahedron/tetrahedron.H
@@ -48,6 +48,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class Istream;
 class Ostream;
 
diff --git a/src/foam/meshes/primitiveShapes/triangle/triangle.H b/src/foam/meshes/primitiveShapes/triangle/triangle.H
index 0b3e3ff67..947e72e7a 100644
--- a/src/foam/meshes/primitiveShapes/triangle/triangle.H
+++ b/src/foam/meshes/primitiveShapes/triangle/triangle.H
@@ -45,6 +45,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class Istream;
 class Ostream;
 
diff --git a/src/foam/primitives/BlockCoeff/scalarBlockCoeff.H b/src/foam/primitives/BlockCoeff/scalarBlockCoeff.H
index cefd63035..d7385e44a 100644
--- a/src/foam/primitives/BlockCoeff/scalarBlockCoeff.H
+++ b/src/foam/primitives/BlockCoeff/scalarBlockCoeff.H
@@ -42,6 +42,11 @@ Author
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+
+Ostream& operator<<(Ostream&, const BlockCoeff<scalar>&);
+
+
 /*---------------------------------------------------------------------------*\
                           Class BlockCoeff Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/foam/primitives/Lists/hashedWordList.H b/src/foam/primitives/Lists/hashedWordList.H
index d06a05d9a..3a1e38c6e 100644
--- a/src/foam/primitives/Lists/hashedWordList.H
+++ b/src/foam/primitives/Lists/hashedWordList.H
@@ -44,9 +44,9 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
 class hashedWordList;
 
-// Forward declaration of friend functions and operators
 Istream& operator>>(Istream&, hashedWordList&);
 Ostream& operator<<(Ostream&, const hashedWordList&);
 
diff --git a/src/foam/primitives/bools/Switch/Switch.H b/src/foam/primitives/bools/Switch/Switch.H
index a568f7a81..6ee713940 100644
--- a/src/foam/primitives/bools/Switch/Switch.H
+++ b/src/foam/primitives/bools/Switch/Switch.H
@@ -45,11 +45,12 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declaration of friend functions and operators
-
-class Switch;
+// Forward declaration of classes
 class dictionary;
 
+// Forward declaration of friend functions and operators
+class Switch;
+
 Istream& operator>>(Istream&, Switch&);
 Ostream& operator<<(Ostream&, const Switch&);
 
diff --git a/src/foam/primitives/strings/fileName/fileName.H b/src/foam/primitives/strings/fileName/fileName.H
index 1644f8734..a75347991 100644
--- a/src/foam/primitives/strings/fileName/fileName.H
+++ b/src/foam/primitives/strings/fileName/fileName.H
@@ -53,6 +53,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 template<class T> class List;
 typedef List<word> wordList;
 
diff --git a/src/foam/primitives/strings/keyType/keyType.H b/src/foam/primitives/strings/keyType/keyType.H
index 3e7defe4e..32a3b52ab 100644
--- a/src/foam/primitives/strings/keyType/keyType.H
+++ b/src/foam/primitives/strings/keyType/keyType.H
@@ -49,6 +49,12 @@ namespace Foam
 class Istream;
 class Ostream;
 
+// Forward declaration of friend functions and operators
+class keyType;
+
+Istream& operator>>(Istream&, keyType&);
+Ostream& operator<<(Ostream&, const keyType&);
+
 
 /*---------------------------------------------------------------------------*\
                           Class keyType Declaration
diff --git a/src/foam/primitives/strings/wordRe/wordRe.H b/src/foam/primitives/strings/wordRe/wordRe.H
index bba4afc26..7774d631c 100644
--- a/src/foam/primitives/strings/wordRe/wordRe.H
+++ b/src/foam/primitives/strings/wordRe/wordRe.H
@@ -59,11 +59,13 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declaration of friend functions and operators
-class wordRe;
+// Forward declaration of classes
 class Istream;
 class Ostream;
 
+// Forward declaration of friend functions and operators
+class wordRe;
+
 Istream& operator>>(Istream&, wordRe&);
 Ostream& operator<<(Ostream&, const wordRe&);
 
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 5dd5add9a..29eb8f5a3 100644
--- a/src/immersedBoundary/immersedBoundaryTurbulence/wallFunctions/immersedBoundaryKqRWallFunction/immersedBoundaryKqRWallFunctionFvPatchField.C
+++ b/src/immersedBoundary/immersedBoundaryTurbulence/wallFunctions/immersedBoundaryKqRWallFunction/immersedBoundaryKqRWallFunctionFvPatchField.C
@@ -153,7 +153,7 @@ immersedBoundaryKqRWallFunctionFvPatchField
     (
         ptf.ibPatch(),
         ptf.setDeadValue(),
-        ptf.deadValue()        
+        ptf.deadValue()
     )
 {
     this->setPatchType(ptf);
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/basic/Particle/Particle.H b/src/lagrangian/basic/Particle/Particle.H
index fa94334b1..0c7e0fef3 100644
--- a/src/lagrangian/basic/Particle/Particle.H
+++ b/src/lagrangian/basic/Particle/Particle.H
@@ -44,6 +44,7 @@ Description
 namespace Foam
 {
 
+// Forward declaration of classes
 template<class Particle>
 class Cloud;
 
diff --git a/src/lagrangian/dieselSpray/parcel/parcel.H b/src/lagrangian/dieselSpray/parcel/parcel.H
index 24f035c07..65def467c 100644
--- a/src/lagrangian/dieselSpray/parcel/parcel.H
+++ b/src/lagrangian/dieselSpray/parcel/parcel.H
@@ -40,8 +40,15 @@ Description
 namespace Foam
 {
 
+// Forward declaration of classes
 class spray;
 
+// Forward declaration of friend functions and operators
+class parcel;
+
+Ostream& operator<<(Ostream&, const parcel&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class parcel Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/lagrangian/intermediate/phaseProperties/phaseProperties/phaseProperties.H b/src/lagrangian/intermediate/phaseProperties/phaseProperties/phaseProperties.H
index db83c51a1..b0e0c6a3b 100644
--- a/src/lagrangian/intermediate/phaseProperties/phaseProperties/phaseProperties.H
+++ b/src/lagrangian/intermediate/phaseProperties/phaseProperties/phaseProperties.H
@@ -46,6 +46,13 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+class phaseProperties;
+
+Istream& operator>>(Istream&, phaseProperties&);
+Ostream& operator<<(Ostream&, const phaseProperties&);
+
+
 /*---------------------------------------------------------------------------*\
                       Class phaseProperties Declaration
 \*---------------------------------------------------------------------------*/
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/lagrangian/molecularDynamics/molecule/interactionLists/referralLists/receivingReferralList.H b/src/lagrangian/molecularDynamics/molecule/interactionLists/referralLists/receivingReferralList.H
index b077d8dd9..8f5a360ac 100644
--- a/src/lagrangian/molecularDynamics/molecule/interactionLists/referralLists/receivingReferralList.H
+++ b/src/lagrangian/molecularDynamics/molecule/interactionLists/referralLists/receivingReferralList.H
@@ -43,6 +43,13 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+class receivingReferralList;
+
+Istream& operator>>(Istream&, receivingReferralList&);
+Ostream& operator<<(Ostream&, const receivingReferralList&);
+
+
 /*---------------------------------------------------------------------------*\
                     Class receivingReferralList Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/lagrangian/molecularDynamics/molecule/interactionLists/referralLists/sendingReferralList.H b/src/lagrangian/molecularDynamics/molecule/interactionLists/referralLists/sendingReferralList.H
index d7315eb8b..9de42a022 100644
--- a/src/lagrangian/molecularDynamics/molecule/interactionLists/referralLists/sendingReferralList.H
+++ b/src/lagrangian/molecularDynamics/molecule/interactionLists/referralLists/sendingReferralList.H
@@ -43,6 +43,13 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+class sendingReferralList;
+
+Istream& operator>>(Istream&, sendingReferralList&);
+Ostream& operator<<(Ostream&, const sendingReferralList&);
+
+
 /*---------------------------------------------------------------------------*\
                       Class sendingReferralList Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/lagrangian/molecularDynamics/molecule/interactionLists/referredMolecule/referredMolecule.H b/src/lagrangian/molecularDynamics/molecule/interactionLists/referredMolecule/referredMolecule.H
index a3ee23462..ba5d39af8 100644
--- a/src/lagrangian/molecularDynamics/molecule/interactionLists/referredMolecule/referredMolecule.H
+++ b/src/lagrangian/molecularDynamics/molecule/interactionLists/referredMolecule/referredMolecule.H
@@ -43,6 +43,16 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+class referredMolecule;
+
+inline bool operator==(const referredMolecule& a, const referredMolecule& b);
+inline bool operator!=(const referredMolecule& a, const referredMolecule& b);
+
+Istream& operator>>(Istream&, referredMolecule&);
+Ostream& operator<<(Ostream&, const referredMolecule&);
+
+
 /*---------------------------------------------------------------------------*\
                       Class referredMolecule Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
index f68ec313f..9fb06a631 100644
--- a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
+++ b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.H
@@ -47,8 +47,12 @@ SourceFiles
 namespace Foam
 {
 
-// Class forward declarations
+// Forward declaration of friend functions and operators
 class moleculeCloud;
+class molecule;
+
+Ostream& operator<<(Ostream&, const molecule&);
+
 
 /*---------------------------------------------------------------------------*\
                            Class molecule Declaration
diff --git a/src/lagrangian/molecularDynamics/molecule/reducedUnits/reducedUnits.H b/src/lagrangian/molecularDynamics/molecule/reducedUnits/reducedUnits.H
index 1e4f36f0e..4e5ec6a45 100644
--- a/src/lagrangian/molecularDynamics/molecule/reducedUnits/reducedUnits.H
+++ b/src/lagrangian/molecularDynamics/molecule/reducedUnits/reducedUnits.H
@@ -44,6 +44,12 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+class reducedUnits;
+
+Ostream& operator<<(Ostream&, const reducedUnits&);
+
+
 /*---------------------------------------------------------------------------*\
                          Class reducedUnits Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/lagrangian/solidParticle/solidParticle.H b/src/lagrangian/solidParticle/solidParticle.H
index 0331cf3f1..945c67d70 100644
--- a/src/lagrangian/solidParticle/solidParticle.H
+++ b/src/lagrangian/solidParticle/solidParticle.H
@@ -49,7 +49,12 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
 class solidParticleCloud;
+class solidParticle;
+
+Ostream& operator<<(Ostream&, const solidParticle&);
+
 
 /*---------------------------------------------------------------------------*\
                            Class solidParticle Declaration
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointData.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointData.H
index 54283f2f4..14d54abb3 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointData.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointData.H
@@ -49,10 +49,17 @@ SourceFiles
 namespace Foam
 {
 
-// Class forward declarations
+// Forward declaration of classes
 class polyPatch;
 class polyMesh;
 
+// Forward declaration of friend functions and operators
+class pointData;
+
+Ostream& operator<<(Ostream&, const pointData&);
+Istream& operator>>(Istream&, pointData&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class pointData Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/mesh/autoMesh/autoHexMesh/trackedParticle/ExactParticle.H b/src/mesh/autoMesh/autoHexMesh/trackedParticle/ExactParticle.H
index 1bb359827..c415b1e3c 100644
--- a/src/mesh/autoMesh/autoHexMesh/trackedParticle/ExactParticle.H
+++ b/src/mesh/autoMesh/autoHexMesh/trackedParticle/ExactParticle.H
@@ -40,6 +40,7 @@ Description
 namespace Foam
 {
 
+// Forward declaration of classes
 template<class ExactParticle>
 class Cloud;
 
diff --git a/src/mesh/autoMesh/autoHexMesh/trackedParticle/trackedParticle.H b/src/mesh/autoMesh/autoHexMesh/trackedParticle/trackedParticle.H
index 67402983e..13a1fd4cc 100644
--- a/src/mesh/autoMesh/autoHexMesh/trackedParticle/trackedParticle.H
+++ b/src/mesh/autoMesh/autoHexMesh/trackedParticle/trackedParticle.H
@@ -45,8 +45,15 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class trackedParticleCloud;
 
+// Forward declaration of friend functions and operators
+class trackedParticle;
+
+Ostream& operator<<(Ostream&, const trackedParticle&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class trackedParticle Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/mesh/blockMesh/block/block.H b/src/mesh/blockMesh/block/block.H
index 1a9564b99..4ca927882 100644
--- a/src/mesh/blockMesh/block/block.H
+++ b/src/mesh/blockMesh/block/block.H
@@ -50,13 +50,16 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of class
 class Istream;
 class Ostream;
 
 // Forward declaration of friend functions and operators
 class block;
+
 Ostream& operator<<(Ostream&, const block&);
 
+
 /*---------------------------------------------------------------------------*\
                            Class block Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H
index 1fe113120..2535d8266 100644
--- a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H
+++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H
@@ -47,14 +47,17 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class Istream;
 class Ostream;
 
 // Forward declaration of friend functions and operators
 class blockMesh;
 class blockDescriptor;
+
 Ostream& operator<<(Ostream&, const blockDescriptor&);
 
+
 /*---------------------------------------------------------------------------*\
                        Class blockDescriptor Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/mesh/blockMesh/curvedEdges/curvedEdge.H b/src/mesh/blockMesh/curvedEdges/curvedEdge.H
index 4f3294d52..cca97ef81 100644
--- a/src/mesh/blockMesh/curvedEdges/curvedEdge.H
+++ b/src/mesh/blockMesh/curvedEdges/curvedEdge.H
@@ -47,6 +47,12 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+class curvedEdge;
+
+Ostream& operator<<(Ostream&, const curvedEdge&);
+
+
 /*---------------------------------------------------------------------------*\
                          Class curvedEdge Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/mesh/cfMesh/utilities/containers/DynList/DynListI.H b/src/mesh/cfMesh/utilities/containers/DynList/DynListI.H
index 3692a12cd..a3728f895 100644
--- a/src/mesh/cfMesh/utilities/containers/DynList/DynListI.H
+++ b/src/mesh/cfMesh/utilities/containers/DynList/DynListI.H
@@ -150,6 +150,25 @@ inline Foam::DynList<T, staticSize>::DynList(const label s)
     # endif
 }
 
+
+#if WM_LABEL_SIZE == 64
+template<class T, Foam::label staticSize>
+inline Foam::DynList<T, staticSize>::DynList(const int s)
+:
+    dataPtr_(nullptr),
+    nAllocated_(0),
+    staticData_(),
+    nextFree_(0)
+{
+    setSize(s);
+
+    # ifdef DEBUG
+    checkAllocation();
+    # endif
+}
+#endif
+
+
 template<class T, Foam::label staticSize>
 inline Foam::DynList<T, staticSize>::DynList(const label s, const T& val)
 :
diff --git a/src/mesh/cfMesh/utilities/decomposeCells/decomposeCellsPyramids.C b/src/mesh/cfMesh/utilities/decomposeCells/decomposeCellsPyramids.C
index eff6ef4bf..16daf21ac 100644
--- a/src/mesh/cfMesh/utilities/decomposeCells/decomposeCellsPyramids.C
+++ b/src/mesh/cfMesh/utilities/decomposeCells/decomposeCellsPyramids.C
@@ -58,8 +58,10 @@ void decomposeCells::findAddressingForCell
     const faceListPMG& faces = mesh_.faces();
     forAll(faceEdges, feI)
     {
-        faceEdges[feI].setSize(faces[c[feI]].size());
-        faceEdges[feI] = -1;
+        DynList<label, 8>& fEdges = faceEdges[feI];
+
+        fEdges.setSize(faces[c[feI]].size());
+        fEdges = label(-1);
     }
 
     forAll(c, fI)
diff --git a/src/mesh/cfMesh/utilities/helperFunctions/helperFunctionsGeometryQueriesI.H b/src/mesh/cfMesh/utilities/helperFunctions/helperFunctionsGeometryQueriesI.H
index fb7a7eda7..170876362 100644
--- a/src/mesh/cfMesh/utilities/helperFunctions/helperFunctionsGeometryQueriesI.H
+++ b/src/mesh/cfMesh/utilities/helperFunctions/helperFunctionsGeometryQueriesI.H
@@ -1221,7 +1221,7 @@ inline bool doTrianglesOverlap
         x /= (mag(x) + VSMALL);
         vector y = vec ^ x;
 
-        DynList<point2D, 6> poly2D(3);
+	DynList<point2D, 6> poly2D(label(3));
         poly2D[0] = point2D((tri0.a() - origin) & x, (tri0.a() - origin) & y);
         poly2D[1] = point2D((tri0.b() - origin) & x, (tri0.b() - origin) & y);
         poly2D[2] = point2D((tri0.c() - origin) & x, (tri0.c() - origin) & y);
diff --git a/src/mesh/cfMesh/utilities/meshes/partTriMesh/partTriMeshParallelAddressing.C b/src/mesh/cfMesh/utilities/meshes/partTriMesh/partTriMeshParallelAddressing.C
index 14075472c..b4199e149 100644
--- a/src/mesh/cfMesh/utilities/meshes/partTriMesh/partTriMeshParallelAddressing.C
+++ b/src/mesh/cfMesh/utilities/meshes/partTriMesh/partTriMeshParallelAddressing.C
@@ -354,7 +354,7 @@ void partTriMesh::createBufferLayers()
         {
             const parTriFace& tri = receivedTrias[i];
 
-            DynList<label, 3> triPointLabels(3);
+	    DynList<label, 3> triPointLabels(label(3));
             for(label j=0;j<3;++j)
             {
                 const label gpI = tri.globalLabelOfPoint(j);
diff --git a/src/mesh/cfMesh/utilities/octrees/meshOctree/meshOctree.C b/src/mesh/cfMesh/utilities/octrees/meshOctree/meshOctree.C
index e34146e20..498f98d75 100644
--- a/src/mesh/cfMesh/utilities/octrees/meshOctree/meshOctree.C
+++ b/src/mesh/cfMesh/utilities/octrees/meshOctree/meshOctree.C
@@ -136,7 +136,7 @@ void meshOctree::setOctantVectorsAndPositions()
     //- set vrtLeavesPos_
     for(label vrtI=0;vrtI<8;++vrtI)
     {
-        FixedList<label, 3> vc(0);
+	FixedList<label, 3> vc(label(0));
 
         if( vrtI & 1 )
             vc[0] += 1;
diff --git a/src/mesh/cfMesh/utilities/octrees/meshOctree/meshOctreeAutomaticRefinement/meshOctreeAutomaticRefinement.C b/src/mesh/cfMesh/utilities/octrees/meshOctree/meshOctreeAutomaticRefinement/meshOctreeAutomaticRefinement.C
index 28153aefa..c9ee5a5f5 100644
--- a/src/mesh/cfMesh/utilities/octrees/meshOctree/meshOctreeAutomaticRefinement/meshOctreeAutomaticRefinement.C
+++ b/src/mesh/cfMesh/utilities/octrees/meshOctree/meshOctreeAutomaticRefinement/meshOctreeAutomaticRefinement.C
@@ -138,7 +138,7 @@ void meshOctreeAutomaticRefinement::setMaxRefLevel()
         {
             finished = false;
 
-            const scalar lSize = size / pow(2, label(maxRefLevel_));
+	    const scalar lSize = size / pow(2.0, label(maxRefLevel_));
 
             if( lSize < cs )
             {
diff --git a/src/mesh/cfMesh/utilities/octrees/meshOctree/meshOctreeCreator/meshOctreeCreatorCreateOctreeBoxes.C b/src/mesh/cfMesh/utilities/octrees/meshOctree/meshOctreeCreator/meshOctreeCreatorCreateOctreeBoxes.C
index 45af48b60..257120620 100644
--- a/src/mesh/cfMesh/utilities/octrees/meshOctree/meshOctreeCreator/meshOctreeCreatorCreateOctreeBoxes.C
+++ b/src/mesh/cfMesh/utilities/octrees/meshOctree/meshOctreeCreator/meshOctreeCreatorCreateOctreeBoxes.C
@@ -87,12 +87,12 @@ void meshOctreeCreator::setRootCubeSizeAndRefParameters()
     {
         finished = false;
 
-        const scalar lSize = size / Foam::pow(2, label(globalRefLevel_));
+	const scalar lSize = size / pow(label(2), label(globalRefLevel_));
 
         if( lSize < (maxSize * (1.0-SMALL)) )
         {
             const scalar bbSize =
-                0.5 * maxSize * Foam::pow(2, label(globalRefLevel_));
+		0.5 * maxSize * pow(label(2), label(globalRefLevel_));
             rootBox.max() = c + point(bbSize, bbSize, bbSize);
             rootBox.min() = c - point(bbSize, bbSize, bbSize);
             finished = true;
@@ -136,7 +136,7 @@ void meshOctreeCreator::setRootCubeSizeAndRefParameters()
         {
             finished = false;
 
-            const scalar lSize = maxSize / Foam::pow(2, addLevel);
+	    const scalar lSize = maxSize / Foam::pow(label(2), addLevel);
 
             if( lSize <= cs )
             {
@@ -233,7 +233,7 @@ void meshOctreeCreator::setRootCubeSizeAndRefParameters()
             {
                 finished = false;
 
-                const scalar lSize = maxSize / Foam::pow(2, addLevel);
+		const scalar lSize = maxSize / Foam::pow(label(2), addLevel);
 
                 if( lSize <= cs )
                 {
@@ -311,7 +311,7 @@ void meshOctreeCreator::setRootCubeSizeAndRefParameters()
             {
                 finished = false;
 
-                const scalar lSize = maxSize / Foam::pow(2, addLevel);
+		const scalar lSize = maxSize / Foam::pow(label(2), addLevel);
 
                 if( lSize <= cs )
                 {
@@ -383,7 +383,7 @@ void meshOctreeCreator::setRootCubeSizeAndRefParameters()
                     {
                         finished = false;
 
-                        const scalar lSize = maxSize / Foam::pow(2, nLevel);
+			const scalar lSize = maxSize / Foam::pow(label(2), nLevel);
 
                         if( lSize <= cs )
                         {
diff --git a/src/mesh/cfMesh/utilities/smoothers/geometry/meshOptimizer/tetMeshOptimisation/advancedSmoothers/meshUntangler/meshUntanglerCutRegionFaces.C b/src/mesh/cfMesh/utilities/smoothers/geometry/meshOptimizer/tetMeshOptimisation/advancedSmoothers/meshUntangler/meshUntanglerCutRegionFaces.C
index 83964742d..823fec41f 100644
--- a/src/mesh/cfMesh/utilities/smoothers/geometry/meshOptimizer/tetMeshOptimisation/advancedSmoothers/meshUntangler/meshUntanglerCutRegionFaces.C
+++ b/src/mesh/cfMesh/utilities/smoothers/geometry/meshOptimizer/tetMeshOptimisation/advancedSmoothers/meshUntangler/meshUntanglerCutRegionFaces.C
@@ -75,7 +75,7 @@ void meshUntangler::cutRegion::findNewFaces()
                 << " consisting of edges " << f << endl;
             # endif
 
-            pointUsage = 0;
+	    pointUsage = label(0);
 
             DynList<label, 8> newFace;
 
@@ -147,7 +147,7 @@ void meshUntangler::cutRegion::findNewFaces()
         //- find edges which form the faceInPlane
         DynList<label, 128> edgeUsage;
         edgeUsage.setSize(cEdges.size());
-        edgeUsage = 0;
+	edgeUsage = label(0);
         forAll(cFaces, fI)
         {
             const DynList<label, 8>& f = cFaces[fI];
diff --git a/src/mesh/cfMesh/utilities/surfaceTools/edgeExtraction/edgeExtractor/edgeExtractor.C b/src/mesh/cfMesh/utilities/surfaceTools/edgeExtraction/edgeExtractor/edgeExtractor.C
index a97e3f6c4..957eb556c 100644
--- a/src/mesh/cfMesh/utilities/surfaceTools/edgeExtraction/edgeExtractor/edgeExtractor.C
+++ b/src/mesh/cfMesh/utilities/surfaceTools/edgeExtraction/edgeExtractor/edgeExtractor.C
@@ -29,8 +29,6 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "error.H"
-#include "objectRegistry.H"
-#include "foamTime.H"
 #include "polyMeshGenModifier.H"
 #include "edgeExtractor.H"
 #include "meshSurfaceEngine.H"
@@ -1529,7 +1527,7 @@ bool edgeExtractor::checkConcaveEdgeCells()
 
                         DynList<label, 2> nFacesInPatch;
                         nFacesInPatch.setSize(2);
-                        nFacesInPatch = 0;
+			nFacesInPatch = label(0);
 
                         DynList<bool, 2> hasPatchPoints;
                         hasPatchPoints.setSize(2);
diff --git a/src/meshTools/PointEdgeWave/pointEdgePoint.H b/src/meshTools/PointEdgeWave/pointEdgePoint.H
index ca9d8fb05..e5959474d 100644
--- a/src/meshTools/PointEdgeWave/pointEdgePoint.H
+++ b/src/meshTools/PointEdgeWave/pointEdgePoint.H
@@ -54,6 +54,13 @@ namespace Foam
 class polyPatch;
 class polyMesh;
 
+// Forward declaration of friend functions and operators
+class pointEdgePoint;
+
+Ostream& operator<<(Ostream&, const pointEdgePoint&);
+Istream& operator>>(Istream&, pointEdgePoint&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class pointEdgePoint Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/meshTools/cellClassification/cellInfo.H b/src/meshTools/cellClassification/cellInfo.H
index 67d0f56fc..1d7ac6b1f 100644
--- a/src/meshTools/cellClassification/cellInfo.H
+++ b/src/meshTools/cellClassification/cellInfo.H
@@ -45,9 +45,18 @@ SourceFiles
 
 namespace Foam
 {
+
+// Forward declaration of classes
 class polyPatch;
 class polyMesh;
 
+// Forward declaration of friend functions and operators
+class cellInfo;
+
+Ostream& operator<<(Ostream&, const cellInfo&);
+Istream& operator>>(Istream&, cellInfo&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class cellInfo Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/meshTools/cellDist/wallPoint/wallPoint.H b/src/meshTools/cellDist/wallPoint/wallPoint.H
index dcdabd890..a86df4809 100644
--- a/src/meshTools/cellDist/wallPoint/wallPoint.H
+++ b/src/meshTools/cellDist/wallPoint/wallPoint.H
@@ -50,9 +50,10 @@ namespace Foam
 // Forward declaration of classes
 class polyPatch;
 class polyMesh;
-class wallPoint;
 
 // Forward declaration of friend functions and operators
+class wallPoint;
+
 Ostream& operator<<(Ostream&, const wallPoint&);
 Istream& operator>>(Istream&, wallPoint&);
 
diff --git a/src/meshTools/cellDist/wallPoint/wallPointData.H b/src/meshTools/cellDist/wallPoint/wallPointData.H
index f29286570..991c9cf7b 100644
--- a/src/meshTools/cellDist/wallPoint/wallPointData.H
+++ b/src/meshTools/cellDist/wallPoint/wallPointData.H
@@ -46,6 +46,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 template<class Type> class wallPointData;
 
 // Forward declaration of friend functions and operators
diff --git a/src/meshTools/triSurface/surfaceLocation/surfaceLocation.H b/src/meshTools/triSurface/surfaceLocation/surfaceLocation.H
index f3995c90d..b40e4795e 100644
--- a/src/meshTools/triSurface/surfaceLocation/surfaceLocation.H
+++ b/src/meshTools/triSurface/surfaceLocation/surfaceLocation.H
@@ -55,6 +55,14 @@ namespace Foam
 // Forward declaration of classes
 class triSurface;
 
+// Forward declaration of friend functions and operators
+class surfaceLocation;
+
+Istream& operator>>(Istream& is, surfaceLocation& sl);
+Ostream& operator<<(Ostream& os, const surfaceLocation& sl);
+Ostream& operator<<(Ostream&, const InfoProxy<surfaceLocation>&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class surfaceLocation Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/meshTools/triSurface/triSurface/interfaces/STL/STLtriangle.H b/src/meshTools/triSurface/triSurface/interfaces/STL/STLtriangle.H
index 09a341369..002c227ae 100644
--- a/src/meshTools/triSurface/triSurface/interfaces/STL/STLtriangle.H
+++ b/src/meshTools/triSurface/triSurface/interfaces/STL/STLtriangle.H
@@ -42,6 +42,12 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+class STLtriangle;
+
+inline Ostream& operator<<(Ostream&, const STLtriangle&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class STLtriangle Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/meshTools/triSurface/triSurface/surfacePatch/surfacePatch.H b/src/meshTools/triSurface/triSurface/surfacePatch/surfacePatch.H
index 48799b05a..53eaac03c 100644
--- a/src/meshTools/triSurface/triSurface/surfacePatch/surfacePatch.H
+++ b/src/meshTools/triSurface/triSurface/surfacePatch/surfacePatch.H
@@ -45,6 +45,12 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+class surfacePatch;
+
+Ostream& operator<<(Ostream&, const surfacePatch&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class surfacePatch Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/meshTools/triSurface/triSurface/surfacePatch/surfacePatchIOList.H b/src/meshTools/triSurface/triSurface/surfacePatch/surfacePatchIOList.H
index 6b153d21d..f31cc93ce 100644
--- a/src/meshTools/triSurface/triSurface/surfacePatch/surfacePatchIOList.H
+++ b/src/meshTools/triSurface/triSurface/surfacePatch/surfacePatchIOList.H
@@ -46,6 +46,10 @@ namespace Foam
 {
 
 // Forward declaration of classes
+// Forward declaration of friend functions and operators
+class surfacePatchIOList;
+
+Ostream& operator<<(Ostream&, const surfacePatchIOList&);
 
 /*---------------------------------------------------------------------------*\
                            Class surfacePatchIOList Declaration
diff --git a/src/meshTools/triSurface/triSurface/triSurface.H b/src/meshTools/triSurface/triSurface/triSurface.H
index 7d3a9cf73..eac31a00c 100644
--- a/src/meshTools/triSurface/triSurface/triSurface.H
+++ b/src/meshTools/triSurface/triSurface/triSurface.H
@@ -48,9 +48,17 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class Time;
 class IFstream;
 
+// Forward declaration of friend functions and operators
+
+class triSurface;
+
+Ostream& operator<<(Ostream&, const triSurface&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class triSurface Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/overset/oversetMesh/oversetFringe/overlapFringe/adaptiveOverlapFringe/adaptiveOverlapFringe.C b/src/overset/oversetMesh/oversetFringe/overlapFringe/adaptiveOverlapFringe/adaptiveOverlapFringe.C
index 2c679f70e..98282dee7 100644
--- a/src/overset/oversetMesh/oversetFringe/overlapFringe/adaptiveOverlapFringe/adaptiveOverlapFringe.C
+++ b/src/overset/oversetMesh/oversetFringe/overlapFringe/adaptiveOverlapFringe/adaptiveOverlapFringe.C
@@ -95,8 +95,8 @@ void Foam::adaptiveOverlapFringe::suitabilityFractionSlope
     forAllConstIter(FIFOStack<iterationData>, iterHist, it)
     {
         n += (it().iteration() - iterMean)*
-            (it().suitability() - suitabilityMean);
-        d += sqr(it().iteration() - iterMean);
+             (it().suitability() - suitabilityMean);
+        d += Foam::sqr((it().iteration() - iterMean));
     }
 
     reduce(n, sumOp<scalar>());
@@ -389,9 +389,9 @@ Foam::adaptiveOverlapFringe::adaptiveOverlapFringe
 )
 :
     oversetFringe(mesh, region, dict),
-    fringeHolesPtr_(nullptr),
-    acceptorsPtr_(nullptr),
-    finalDonorAcceptorsPtr_(nullptr),
+    fringeHolesPtr_(NULL),
+    acceptorsPtr_(NULL),
+    finalDonorAcceptorsPtr_(NULL),
 
     holesZoneName_(dict.lookupOrDefault<word>("holes", word())),
     initPatchNames_
diff --git a/src/sampling/sampledSurface/sampledSurface/sampledSurface.H b/src/sampling/sampledSurface/sampledSurface/sampledSurface.H
index bafbee892..66c0d934a 100644
--- a/src/sampling/sampledSurface/sampledSurface/sampledSurface.H
+++ b/src/sampling/sampledSurface/sampledSurface/sampledSurface.H
@@ -68,6 +68,12 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+class sampledSurface;
+
+Ostream& operator<<(Ostream&, const sampledSurface&);
+
+
 /*---------------------------------------------------------------------------*\
                       Class sampledSurface Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/surfMesh/surfZone/surfZone/surfZoneIOList.H b/src/surfMesh/surfZone/surfZone/surfZoneIOList.H
index e66008ca5..2de1e0571 100644
--- a/src/surfMesh/surfZone/surfZone/surfZoneIOList.H
+++ b/src/surfMesh/surfZone/surfZone/surfZoneIOList.H
@@ -45,7 +45,11 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declaration of classes
+// Forward declaration of friend functions and operators
+class surfZoneIOList;
+
+Ostream& operator<<(Ostream&, const surfZoneIOList&);
+
 
 /*---------------------------------------------------------------------------*\
                       Class surfZoneIOList Declaration
diff --git a/src/surfMesh/surfZone/surfZoneIdentifier/surfZoneIdentifier.H b/src/surfMesh/surfZone/surfZoneIdentifier/surfZoneIdentifier.H
index 0ee086ece..19df11aec 100644
--- a/src/surfMesh/surfZone/surfZoneIdentifier/surfZoneIdentifier.H
+++ b/src/surfMesh/surfZone/surfZoneIdentifier/surfZoneIdentifier.H
@@ -47,13 +47,15 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class dictionary;
 
 // Forward declaration of friend functions and operators
-
 class surfZoneIdentifier;
+
 Ostream& operator<<(Ostream&, const surfZoneIdentifier&);
 
+
 /*---------------------------------------------------------------------------*\
                     Class surfZoneIdentifier Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/surfMesh/surfaceFormats/stl/STLtriangle.H b/src/surfMesh/surfaceFormats/stl/STLtriangle.H
index bf4f8f9b1..aef7de706 100644
--- a/src/surfMesh/surfaceFormats/stl/STLtriangle.H
+++ b/src/surfMesh/surfaceFormats/stl/STLtriangle.H
@@ -44,6 +44,12 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+class STLtriangle;
+
+Ostream& operator<<(Ostream&, const STLtriangle&);
+
+
 /*---------------------------------------------------------------------------*\
                          Class STLtriangle Declaration
 \*---------------------------------------------------------------------------*/
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/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.H b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.H
index 526cfb036..1b36140c5 100755
--- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.H
+++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.H
@@ -55,6 +55,12 @@ Germany
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+class aungierRedlichKwong;
+
+Ostream& operator<<(Ostream&, const aungierRedlichKwong&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class aungierRedlichKwong Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinson.H b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinson.H
index fe43dbb2d..711c223b2 100755
--- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinson.H
+++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinson.H
@@ -57,6 +57,12 @@ Germany
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+class pengRobinson;
+
+Ostream& operator<<(Ostream&, const pengRobinson&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class Peng Robinson Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwong.H b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwong.H
index 826e70161..b5ab2ccbc 100755
--- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwong.H
+++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwong.H
@@ -50,6 +50,12 @@ Germany
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+class redlichKwong;
+
+Ostream& operator<<(Ostream&, const redlichKwong&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class redlichKwong Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.H b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.H
index 81d487902..c19a471a8 100755
--- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.H
+++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.H
@@ -55,6 +55,12 @@ Germany
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+class soaveRedlichKwong;
+
+Ostream& operator<<(Ostream&, const soaveRedlichKwong&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class Soave Readlich Kwong Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/thermophysicalModels/specie/specie/specie.H b/src/thermophysicalModels/specie/specie/specie.H
index b070825b8..0e9b8bfe0 100644
--- a/src/thermophysicalModels/specie/specie/specie.H
+++ b/src/thermophysicalModels/specie/specie/specie.H
@@ -46,9 +46,16 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class Istream;
 class Ostream;
 
+// Forward declaration of friend functions and operators
+class specie;
+
+Ostream& operator<<(Ostream&, const specie&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class specie Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/thermophysicalModels/specie/speciesTable/speciesTable.H b/src/thermophysicalModels/specie/speciesTable/speciesTable.H
index f3a7d0e4d..24d71ea6a 100644
--- a/src/thermophysicalModels/specie/speciesTable/speciesTable.H
+++ b/src/thermophysicalModels/specie/speciesTable/speciesTable.H
@@ -44,6 +44,12 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+class speciesTable;
+
+Istream& operator>>(Istream&, speciesTable&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class speciesTable Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/src/turbulenceModels/incompressible/LES/derivedFvPatchFields/decayingTurbulence/decayingVorton.H b/src/turbulenceModels/incompressible/LES/derivedFvPatchFields/decayingTurbulence/decayingVorton.H
index f9f0325e7..3f90ec738 100644
--- a/src/turbulenceModels/incompressible/LES/derivedFvPatchFields/decayingTurbulence/decayingVorton.H
+++ b/src/turbulenceModels/incompressible/LES/derivedFvPatchFields/decayingTurbulence/decayingVorton.H
@@ -49,6 +49,16 @@ Author
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+class decayingVorton;
+
+bool operator==(const decayingVorton&, const decayingVorton&);
+bool operator!=(const decayingVorton&, const decayingVorton&);
+
+Ostream& operator<<(Ostream& s, const decayingVorton& vt);
+Istream& operator>>(Istream& s, decayingVorton& vt);
+
+
 /*---------------------------------------------------------------------------*\
                         Class decayingVorton Declaration
 \*---------------------------------------------------------------------------*/
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;