Merge branch 'nextRelease' of ssh://git.code.sf.net/p/foam-extend/foam-extend-4.0 into nextRelease

applications/solvers/immersedBoundary/pimpleDyMIbFoam/pEqn.H

@@ -35,7 +35,8 @@
         }
     }
 
-    // Explicitly relax pressure for momentum corrector except for last corrector
+    // Explicitly relax pressure for momentum corrector
+    // except for last corrector
     if (!pimple.finalIter())
     {
         p.relax();
@@ -54,5 +55,8 @@
             )
         )
       + (1 - UUrf)*U.prevIter();
+
+    // U.correctBoundaryConditions() must be called with relative fluxes
+    // because of inletOutlet BC's
     U.correctBoundaryConditions();
 }

applications/utilities/mesh/manipulation/setSet/setSet.C

@@ -46,7 +46,7 @@ Description
 #include <stdio.h>
 
 
-#if READLINE != 0
+#if HAS_READLINE
 #include <readline/readline.h>
 #include <readline/history.h>
 #endif
@@ -56,7 +56,7 @@ using namespace Foam;
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 
-#if READLINE != 0
+#if HAS_READLINE
 static const char* historyFile = ".setSet";
 #endif
 
@@ -785,7 +785,7 @@ int main(int argc, char *argv[])
 
         fileStreamPtr = new std::ifstream(batchFile.c_str());
     }
-#if READLINE != 0
+#if HAS_READLINE
     else if (!read_history(historyFile))
     {
         Info<< "Successfully read history from " << historyFile << endl;
@@ -829,7 +829,7 @@ int main(int argc, char *argv[])
         }
         else
         {
-#           if READLINE != 0
+#           if HAS_READLINE
             {
                 char* linePtr = readline("readline>");
 

etc/cshrc

@@ -172,7 +172,7 @@ if ( ! $?WM_NCOMPPROCS ) setenv WM_NCOMPPROCS `$WM_PROJECT_DIR/bin/foamGetSystem
 # Run options (floating-point signal handling and memory initialisation)
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 setenv FOAM_SIGFPE
-# setenv FOAM_SETNAN
+setenv FOAM_SETNAN
 
 
 # Detect system type and set environment variables

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
 \*---------------------------------------------------------------------------*/

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
 \*---------------------------------------------------------------------------*/

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
 \*---------------------------------------------------------------------------*/

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
 \*---------------------------------------------------------------------------*/

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
 \*---------------------------------------------------------------------------*/

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
 \*---------------------------------------------------------------------------*/

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
 \*---------------------------------------------------------------------------*/

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
 \*---------------------------------------------------------------------------*/

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();
 

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
 \*---------------------------------------------------------------------------*/

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
 \*---------------------------------------------------------------------------*/

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
                     );
             }
         }

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)
+        // Remember to decrement the index by one (turning index) (see
-        // HJ, 5/Dec/2001
+        // 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)
+        // Remember to decrement the index by one (turning index) (see
-        // HJ, 28/Mar/2009
+        // 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)
+        // Remember to decrement the index by one (turning index) (see
-        // HJ, 28/Mar/2009
+        // 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;

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
+        //  indices will be incremented by 1 and then decremented as
-        //  necessary t avoid the problem of face number zero having no
+        //  necessary to avoid the problem of face number zero having no
-        //  sign.  HJ, 5/Dec/2001
+        //  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
             {

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

src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesRebuild.C

@@ -437,15 +437,13 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
         }
     }
 
-    // Prepare patch reconstruction
+    // Get first boundary mesh size
-    HashTable<label, word> patchNameLookup
+    const label firstBndryMeshSize =
-    (
+        meshes_[firstValidMesh()].boundaryMesh().size();
-        meshes_[firstValidMesh()].boundaryMesh().size()
+
-    );
+    // Prepare patch reconstruction. Note: default key type is word in HashTable
-    DynamicList<word> patchTypeLookup
+    HashTable<label> patchNameLookup(firstBndryMeshSize);
-    (
+    DynamicList<word> patchTypeLookup(firstBndryMeshSize);
-        meshes_[firstValidMesh()].boundaryMesh().size()
-    );
 
     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
+                    // Insert the pair into the hash table
-                    if (!patchNameLookup.found(patchName))
+                    if (patchNameLookup.insert(patchName, nReconPatches))
                     {
-                        // Patch not found.  Add it
+                        // This is the first time we're inserting this patch,
-                        patchNameLookup.insert(patchName, nReconPatches);
+                        // add its type into the list and increment the counter
-                        patchTypeLookup.append(procPatches[patchI].type());
+                        patchTypeLookup.append(patch.type());
-                        nReconPatches++;
+                        ++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();
+    // Loop through all the patch names and collect names and types
-
+    forAllConstIter(HashTable<label>, patchNameLookup, iter)
-    forAll (patchNameToc, pnI)
     {
-        const label pnIndex = patchNameLookup.find(patchNameToc[pnI])();
+        reconPatchNames[iter()] = iter.key();
-
-        reconPatchNames[pnIndex] = patchNameToc[pnI];
-        reconPatchTypes[pnIndex] = patchTypeLookup[pnIndex];
     }
 
+    // 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
+            // Get current mesh
-            nReconPoints += meshes_[procI].nPoints();
+            const fvMesh& procMesh = meshes_[procI];
-            nReconFaces += meshes_[procI].allFaces().size();
-            nReconCells += meshes_[procI].nCells();
 
-            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
+                    const label pnIndex = patchNameLookup.find(patch.name())();
-                    (
-                        procPatches[patchI].name()
-                    )();
 
                     // Check patch name and type
                     if
                     (
-                        procPatches[patchI].name()
+                        patch.name() != reconPatchNames[pnIndex]
-                     != reconPatchNames[pnIndex]
+                     || patch.type() != reconPatchTypes[pnIndex]
-                     || procPatches[patchI].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: "
+                            << patch.name() << " type: " << patch.type()
-                            << procPatches[patchI].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]);
+        // Set size of reconstructed patch face owners with default value of -1
-        reconPatchOwner[patchI] = -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
+    // HR 21.12.18 : A bit of a hack for the moment in order to support the new
-    // to parallel processor points, but working on a PtrList of meshes
+    // method (using global point addressing) and old method (syncing across
-    // on the same processor
+    // processor BCs) of constructing the shared points. The old method is still
-    sharedPoints sharedData(meshes_);
+    // 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
+                    // Location in reconstructed patch where the face is
-                    // is inserted
+                    // inserted
-                    const label insertSlot = indices[i];
+                    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
+            // Next mesh is not valid: simply propagate cell offset
-                << " cell offset: " << cellOffset[procI]
+            if (cellOffset.size() > procI + 1)
-                << endl;
+            {
+                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()
+                            const label masterIndex =
-                                + faceInPatch;
+                                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
+    if (!Pstream::parRun())
-        << "    nPoints = " << reconPoints.size() << nl
+    {
-        << "    nFaces = " << reconFaces.size() << nl
+        Info<< "Global mesh size (final): " << nl
-        << "    nCells = " << nReconCells << nl
+            << "    nPoints = " << reconPoints.size() << nl
-        << "    nPatches = " << reconPatchSizes.size() << nl
+            << "    nFaces = " << reconFaces.size() << nl
-        << "    nPatchFaces = " << reconPatchSizes << endl;
+            << "    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

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_(),

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;
 

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();
 }
 
-
+Foam::sharedPoints::sharedPoints
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+(
-
+    const PtrList<fvMesh>& meshes,
-// Foam::sharedPoints::~sharedPoints()
+    const labelListList& globalPointIndex
-// {}
+)
+:
+    meshes_(meshes),
+    sharedPointAddr_(meshes_.size()),
+    sharedPointLabels_(meshes_.size()),
+    nGlobalPoints_(0)
+{
+    calcSharedPoints(globalPointIndex);
+}
 
 
 // ************************************************************************* //

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

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
 \*---------------------------------------------------------------------------*/

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
 \*---------------------------------------------------------------------------*/

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
 \*---------------------------------------------------------------------------*/

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
 \*---------------------------------------------------------------------------*/

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
 \*---------------------------------------------------------------------------*/

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

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

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);
             }
-
+            else
-            cellsToUnrefine[own] = false;
+            {
-            ++nRemCells;
+                cellsToUnrefine[own] = false;
+                ++nRemCells;
+            }
         }
         else if (neiLevel < (ownLevel - 1))
         {
@@ -977,9 +862,11 @@ Foam::label Foam::refinement::faceConsistentUnrefinement
                     << "Try increasing nUnrefinementBufferLayers. "
                     << abort(FatalError);
             }
-
+            else
-            cellsToUnrefine[nei] = false;
+            {
-            ++nRemCells;
+                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
+                        if (debug)
-                        (
+                        {
-                            "label refinement::"
+                            WarningIn
-                            "edgeConsistentUnrefinement"
+                            (
-                            "(boolList& cellsToUnrefine)"
+                                "label refinement::"
-                        )   << "Cell not marked for unrefinement, indicating a"
+                                "edgeConsistentUnrefinement"
-                            << " previous unnoticed problem with unrefinement."
+                                "(boolList& cellsToUnrefine)"
-                            << nl
+                            )   << "Cell not marked for unrefinement, indicating a"
-                            << "cellI: " << cellI << ", cellJ: " << cellJ
+                                << " previous unnoticed problem with unrefinement."
-                            << nl
+                                << nl
-                            << "Level of cellI: " << cellILevel
+                                << "cellI: " << cellI << ", cellJ: " << cellJ
-                            << ", level of cellJ: " << cellJLevel << nl
+                                << nl
-                            << "This is probably because the refinement and "
+                                << "Level of cellI: " << cellILevel
-                            << "unrefinement regions are very close." << nl
+                                << ", level of cellJ: " << cellJLevel << nl
-                            << "Try increasing nUnrefinementBufferLayers. "
+                                << "This is probably because the refinement and "
-                            << abort(FatalError);
+                                << "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
+                        if (debug)
-                        (
+                        {
-                            "label refinement::"
+                            WarningIn
-                            "edgeConsistentUnrefinement"
+                            (
-                            "(boolList& cellsToUnrefine)"
+                                "label refinement::"
-                        )   << "Cell not marked for unrefinement, indicating a"
+                                "edgeConsistentUnrefinement"
-                            << " previous unnoticed problem with unrefinement."
+                                "(boolList& cellsToUnrefine)"
-                            << nl
+                            )   << "Cell not marked for unrefinement, indicating a"
-                            << "cellI: " << cellI << ", cellJ: " << cellJ
+                                << " previous unnoticed problem with unrefinement."
-                            << nl
+                                << nl
-                            << "Level of cellI: " << cellILevel
+                                << "cellI: " << cellI << ", cellJ: " << cellJ
-                            << ", level of cellJ: " << cellJLevel << nl
+                                << nl
-                            << "This is probably because the refinement and "
+                                << "Level of cellI: " << cellILevel
-                            << "unrefinement regions are very close." << nl
+                                << ", level of cellJ: " << cellJLevel << nl
-                            << "Try increasing nUnrefinementBufferLayers. "
+                                << "This is probably because the refinement and "
-                            << abort(FatalError);
+                                << "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
+    // Sync the new point level. Note: here, we assume that the call to
-    // before the next step in setRefinement. Need to investigate why the sync
+    // updateMesh happened after polyBoundaryMesh::updateMesh where the
-    // is not done properly if I put it here. Something is not updated yet.
+    // processor data is fully rebuilt. If this is not the case, the point
-    // VV, 31/Jan/2018.
+    // 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);

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
+            //- Modifies existing faceI for either new owner/neighbour or new
-            //  points. Reverses if nessecary
+            //  face points. Reverses if nessecary
             void modifyFace
             (
                 polyTopoChange& ref,

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
 \*---------------------------------------------------------------------------*/

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_++;
 

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
 \*---------------------------------------------------------------------------*/

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
 \*---------------------------------------------------------------------------*/

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

src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C

@@ -271,15 +271,11 @@ bool Foam::dynamicPolyRefinementFvMesh::update()
             (
                 refinementSelectionPtr_->refinementCellCandidates()()
             );
-
+            Pout<< "Selected " << refCandidates.size()
-            if (debug)
+                << " 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()()
             );
-
+            Pout<< "Selected " << unrefCandidates.size()
-            if (debug)
+                << " 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();
     }
-
+    Pout<< "No refinement/unrefinement" << endl;
-    Info<< "No refinement/unrefinement" << endl;
-
     // No refinement/unrefinement at this time step. Return false
     return false;
 }

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(),
+                dummyTime,
-                "processorPart" + Foam::name(procI),
+                "",     // 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),
+                            "", // "processorPart" + Foam::name(procI),
-                            time().timeName(),
+                            dummyTime.timeName(),
-                            time(),
+                            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
+            if (procMeshes.set(procI))
-                << " points " << procMeshes[procI].nPoints()
+            {
-                << " faces: " << procMeshes[procI].nFaces()
+                Pout<< "procMesh " << procI
-                << " internal: " << procMeshes[procI].nInternalFaces()
+                    << " points " << procMeshes[procI].nPoints()
-                << " cells: " << procMeshes[procI].nCells()
+                    << " faces: " << procMeshes[procI].nFaces()
-                << " patches: " << procMeshes[procI].boundary().size()
+                    << " internal: " << procMeshes[procI].nInternalFaces()
-                << endl;
+                    << " 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: "
+    if (debug)
-        << " nCells: " << reconMesh.nCells()
+    {
-        << " nFaces: " << reconMesh.nFaces()
+        Pout<< "Reconstructed mesh stats: "
-        << " nIntFaces: " << reconMesh.nInternalFaces()
+            << " nCells: " << reconMesh.nCells()
-        << " polyPatches: "
+            << " nFaces: " << reconMesh.nFaces()
-        << reconMesh.boundaryMesh().size()
+            << " nIntFaces: " << reconMesh.nInternalFaces()
-        << " patches: "
+            << " polyPatches: "
-        << reconMesh.boundary().size()
+            << reconMesh.boundaryMesh().size()
-        << endl;
+            << " 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
+    // HR 13.12.18: Update the mesh objects
-    // checkMesh(true);
+    meshObjectBase::allUpdateTopology<polyMesh>(*this, meshMap);
+
+    if (debug)
+    {
+        Info<< "Checking reconstructed mesh after load balancing..." << endl;
+        checkMesh(true);
+    }
 
     return true;
 }

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 = "
+            if (debug)
-                << masterField.internalField().size()
+            {
-                << " new size = "
+                Pout<< "Resizing internal field: old size = "
-                << GeoField::GeoMeshType::size(masterField.mesh())
+                    << masterField.size()
-                << endl;
+                    << " 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 "
+                if (debug)
-                    << masterField.mesh().boundary()[patchI].name()
+                {
-                    << endl;
+                    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 "
+                if (debug)
-                    << masterField.mesh().boundary()[patchI].name()
+                {
-                    << " old size: " << patchFields[patchI].size()
+                    Pout<< "Resizing patch field for patch "
-                    << " new size: "
+                        << masterField.mesh().boundary()[patchI].name()
-                    << masterField.mesh().boundary()[patchI].size()
+                        << " old size: " << patchFields[patchI].size()
-                    << endl;
+                        << " 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();
     }
 }
 

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
 \*---------------------------------------------------------------------------*/

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
 \*---------------------------------------------------------------------------*/

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

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

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>());
 }
 
 

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
 \*---------------------------------------------------------------------------*/

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 volVectorField& U"
-                "\n    const volScalarField& rAU"
                 "\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 volVectorField& U"
-                "\n    const volScalarField& rAU"
                 "\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 volVectorField& U"
-                "\n    const volScalarField& rAU"
                 "\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 volVectorField& U"
-                "\n    const volScalarField& rAU"
                 "\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 volVectorField& U"
-                "\n    const volScalarField& rAU"
                 "\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 volVectorField& U"
-                "\n    const volScalarField& rAU"
                 "\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 "

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

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;
 

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>

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>

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"))

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)
+        read();
-        {
-            InfoIn
-            (
-                "fvSchemes::fvSchemes(const objectRegistry& obr)"
-            )   << "fvSchemes dictionary not found.  Creating default."
-                << endl;
-        }
-
-        regIOobject::write();
     }
-
-    read();
 }
 
 

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();
 }

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

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
 

src/foam/algorithms/octree/octree/octreeDataPoint.H

@@ -51,8 +51,14 @@ 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
 \*---------------------------------------------------------------------------*/

src/foam/algorithms/octree/octree/treeBoundBox.H

@@ -66,8 +66,18 @@ 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
 \*---------------------------------------------------------------------------*/

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>&);

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>&);

src/foam/containers/HashTables/HashPtrTable/HashPtrTable.H

@@ -43,6 +43,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class Istream;
 class Ostream;
 

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;
 

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>>

src/foam/containers/LinkedLists/accessTypes/ILList/ILList.H

@@ -43,6 +43,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class Istream;
 class Ostream;
 

src/foam/containers/LinkedLists/accessTypes/LList/LList.H

@@ -44,6 +44,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class Istream;
 class Ostream;
 

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

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

src/foam/coordinateSystems/coordinateSystem.H

@@ -143,11 +143,17 @@ 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
+                       Class coordinateSystem Declaration
 \*---------------------------------------------------------------------------*/
 
 class coordinateSystem

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
 \*---------------------------------------------------------------------------*/

src/foam/db/IOstreams/Pstreams/Pstream.C

@@ -264,9 +264,9 @@ void Foam::Pstream::allocatePstreamCommunicator
     if (index == PstreamGlobals::MPIGroups_.size())
     {
         // Extend storage with dummy values
-        MPI_Group newGroup;
+        MPI_Group newGroup = MPI_GROUP_NULL;
         PstreamGlobals::MPIGroups_.append(newGroup);
-        MPI_Comm newComm;
+        MPI_Comm newComm = MPI_COMM_NULL;
         PstreamGlobals::MPICommunicators_.append(newComm);
     }
     else if (index > PstreamGlobals::MPIGroups_.size())

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
 \*---------------------------------------------------------------------------*/

src/foam/db/dictionary/dictionaryEntry/dictionaryEntry.H

@@ -51,6 +51,11 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+class dictionaryEntry;
+Ostream& operator<<(Ostream&, const dictionaryEntry&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class dictionaryEntry Declaration
 \*---------------------------------------------------------------------------*/

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

src/foam/fields/CoeffField/diagTensorCoeffField.H

@@ -47,6 +47,15 @@ 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>
 :

src/foam/fields/CoeffField/scalarCoeffField.H

@@ -45,6 +45,15 @@ 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>
 :

src/foam/fields/CoeffField/sphericalTensorCoeffField.H

@@ -47,6 +47,15 @@ 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>
 :

src/foam/fields/CoeffField/symmTensor4thOrderCoeffField.H

@@ -47,6 +47,15 @@ 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>
 :

src/foam/fields/CoeffField/symmTensorCoeffField.H

@@ -47,6 +47,15 @@ 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>
 :

src/foam/fields/CoeffField/tensorCoeffField.H

@@ -47,6 +47,15 @@ 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>
 :

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

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;
 

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()

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

src/foam/interpolations/interpolationTable/interpolationTable.C

@@ -273,6 +273,7 @@ Type Foam::interpolationTable<Type>::rateOfChange(const scalar value) const
                     << "    Zero rate of change."
                     << endl;
                 // fall-through to 'CLAMP'
+                [[fallthrough]];
             }
             case interpolationTable::CLAMP:
             {
@@ -306,6 +307,7 @@ Type Foam::interpolationTable<Type>::rateOfChange(const scalar value) const
                     << "    Zero rate of change."
                     << endl;
                 // fall-through to 'CLAMP'
+                [[fallthrough]];
             }
             case interpolationTable::CLAMP:
             {
@@ -415,6 +417,7 @@ Foam::interpolationTable<Type>::operator[](const label i) const
                     << "    Continuing with the first entry"
                     << endl;
                 // fall-through to 'CLAMP'
+                [[fallthrough]];
             }
             case interpolationTable::CLAMP:
             {
@@ -449,6 +452,7 @@ Foam::interpolationTable<Type>::operator[](const label i) const
                     << "    Continuing with the last entry"
                     << endl;
                 // fall-through to 'CLAMP'
+                [[fallthrough]];
             }
             case interpolationTable::CLAMP:
             {
@@ -502,6 +506,7 @@ Type Foam::interpolationTable<Type>::operator()(const scalar value) const
                     << "    Continuing with the first entry"
                     << endl;
                 // fall-through to 'CLAMP'
+                [[fallthrough]];
             }
             case interpolationTable::CLAMP:
             {
@@ -535,6 +540,7 @@ Type Foam::interpolationTable<Type>::operator()(const scalar value) const
                     << "    Continuing with the last entry"
                     << endl;
                 // fall-through to 'CLAMP'
+                [[fallthrough]];
             }
             case interpolationTable::CLAMP:
             {

src/foam/matrices/LUscalarMatrix/procLduInterface.H

@@ -43,8 +43,14 @@ 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
 \*---------------------------------------------------------------------------*/

src/foam/matrices/LUscalarMatrix/procLduMatrix.H

@@ -45,6 +45,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class procLduInterface;
 class lduMatrix;
 

src/foam/matrices/crMatrix/crAddressing.H

@@ -46,6 +46,11 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+class crAddressing;
+Ostream& operator<<(Ostream&, const crAddressing&);
+
+
 /*---------------------------------------------------------------------------*\
                         Class crAddressing Declaration
 \*---------------------------------------------------------------------------*/

src/foam/matrices/crMatrix/crMatrix.H

@@ -46,6 +46,11 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+class crMatrix;
+Ostream& operator<<(Ostream&, const crMatrix&);
+
+
 /*---------------------------------------------------------------------------*\
                            Class crMatrix Declaration
 \*---------------------------------------------------------------------------*/

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);
 }
 
 

src/foam/meshes/boundBox/boundBox.H

@@ -39,12 +39,13 @@ Description
 namespace Foam
 {
 
+template<class T> class tmp;
+
 // Forward declaration of friend functions and operators
 
 class boundBox;
-template<class T> class tmp;
+Istream& operator>>(Istream& is, boundBox&);
-
+Ostream& operator<<(Ostream& os, const boundBox&);
-Ostream& operator<<(Ostream& os, const boundBox& b);
 
 
 /*---------------------------------------------------------------------------*\

src/foam/meshes/meshShapes/face/face.H

@@ -57,18 +57,20 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declaration of friend functions and operators
+// Forward declaration of classes
-
-class face;
 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
 \*---------------------------------------------------------------------------*/

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&);
 

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

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
 
 

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&);
 
 

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();

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
+        );
 };
 
 

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);

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()
                 );