From 7bc0c8b3d8fdfd4ca6f5abb365fa835fc612301a Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Thu, 4 Jul 2019 11:27:14 +0100
Subject: [PATCH 01/11] Bugfix: addes support for dealing with zero-sized patch
 fields (immersed boundary)

---
 .../finiteVolume/fvFieldDecomposerDecomposeFields.C  | 12 ++++++++++--
 1 file changed, 10 insertions(+), 2 deletions(-)

diff --git a/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/fvFieldDecomposerDecomposeFields.C b/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/fvFieldDecomposerDecomposeFields.C
index 80d92197e..5611b07bd 100644
--- a/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/fvFieldDecomposerDecomposeFields.C
+++ b/src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/fvFieldDecomposerDecomposeFields.C
@@ -158,9 +158,17 @@ fvFieldDecomposer::decomposeField
 
         const label patchStart = field.mesh().boundaryMesh()[patchI].start();
 
-        forAll (p, i)
+        // Bugfix
+        // Special patches soch as overset and immersed boundary have
+        // their faces and values in a separate list; the list of
+        // values does not correspond to faces.  Skip such patches.
+        // HJ, 4/Jul/2019
+        if (p.size() == field.mesh().boundaryMesh()[patchI].size())
         {
-            allFaceField[patchStart + i] = p[i];
+            forAll (p, i)
+            {
+                allFaceField[patchStart + i] = p[i];
+            }
         }
     }
 

From d5030fc035ef6f9e97d3ec47cd9c8bcaf3f2f87d Mon Sep 17 00:00:00 2001
From: Vuko Vukcevic <v.vukcevic@wikki.co.uk>
Date: Tue, 28 May 2019 20:52:04 +0200
Subject: [PATCH 02/11] Removed temporary debug stuff

1. Writing down mesh during load balancing,
2. Left-over sleep statement,
3. Output formatting (some output wrapped around debug statements)
---
 .../dynamicPolyRefinementFvMesh.C                | 16 +++++++++-------
 .../topoChangerFvMeshLoadBalance.C               | 15 ++++++++++-----
 2 files changed, 19 insertions(+), 12 deletions(-)

diff --git a/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C b/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C
index 5a5eacf23..680e7186b 100644
--- a/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C
+++ b/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C
@@ -271,13 +271,13 @@ bool Foam::dynamicPolyRefinementFvMesh::update()
             (
                 refinementSelectionPtr_->refinementCellCandidates()()
             );
-            Pout<< "Selected " << refCandidates.size()
+            Info<< "Selected " << refCandidates.size()
                 << " refinement candidates."
                 << endl;
         }
         else
         {
-            Pout<< "Skipping refinement for this time-step..." << endl;
+            Info<< "Skipping refinement for this time-step..." << endl;
         }
 
         // Set cells to refine. Note: refinement needs to make sure that face
@@ -297,13 +297,13 @@ bool Foam::dynamicPolyRefinementFvMesh::update()
             (
                 refinementSelectionPtr_->unrefinementPointCandidates()()
             );
-            Pout<< "Selected " << unrefCandidates.size()
+            Info<< "Selected " << unrefCandidates.size()
                 << " unrefinement candidates."
                 << endl;
         }
         else
         {
-            Pout<< "Skipping unrefinement for this time-step..." << endl;
+            Info<< "Skipping unrefinement for this time-step..." << endl;
         }
 
         // Set split points to unrefine around.
@@ -344,13 +344,13 @@ bool Foam::dynamicPolyRefinementFvMesh::update()
         // some topo changes
         if (sizeCellMap)
         {
-            Pout<< "Successfully performed polyhedral refinement. "
+            Info<< "Successfully performed polyhedral refinement. "
                 << "Changed from " << nOldCells << " to " << sizeCellMap
                 << " cells." << endl;
         }
         else
         {
-            Pout<< "Refinement/unrefinement not performed in this time step "
+            Info<< "Refinement/unrefinement not performed in this time step "
                 << "since no cells were selected." << endl;
         }
 
@@ -362,7 +362,9 @@ 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;
 }
diff --git a/src/dynamicMesh/topoChangerFvMesh/topoChangerFvMesh/topoChangerFvMeshLoadBalance.C b/src/dynamicMesh/topoChangerFvMesh/topoChangerFvMesh/topoChangerFvMeshLoadBalance.C
index 63349bb8a..8285623ee 100644
--- a/src/dynamicMesh/topoChangerFvMesh/topoChangerFvMesh/topoChangerFvMeshLoadBalance.C
+++ b/src/dynamicMesh/topoChangerFvMesh/topoChangerFvMesh/topoChangerFvMeshLoadBalance.C
@@ -291,7 +291,7 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
                 true    // Create passive processor patches
             );
             fvMesh& procMesh = procMeshPtr();
-            procMesh.write();
+
             // Create a field decomposer
             fvFieldDecomposer fieldDecomposer
             (
@@ -304,7 +304,11 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
 
             if (procI != Pstream::myProcNo())
             {
-                Pout<< "Send mesh and fields to processor " << procI << endl;
+                if (debug)
+                {
+                    Pout<< "Send mesh and fields to processor " << procI
+                        << endl;
+                }
 
                 OPstream toProc
                 (
@@ -449,8 +453,6 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
         }
     }
 
-    sleep(2);
-    
     // Collect pieces of mesh and fields from other processors
     for (label procI = 0; procI < meshDecomp.nProcs(); procI++)
     {
@@ -459,7 +461,10 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
             // Check if there is a mesh to send
             if (migratedCells[procI][Pstream::myProcNo()] > 0)
             {
-                Pout<< "Receive mesh and fields from " << procI << endl;
+                if (debug)
+                {
+                    Pout<< "Receive mesh and fields from " << procI << endl;
+                }
 
                 // Note: communication can be optimised.  HJ, 27/Feb/2018
                 IPstream fromProc

From 44cef61dbe0f31256e2f33c6c01d7b691a3e88ab Mon Sep 17 00:00:00 2001
From: Vuko Vukcevic <v.vukcevic@wikki.co.uk>
Date: Thu, 30 May 2019 15:46:02 +0200
Subject: [PATCH 03/11] Clean-up after merge (duplicate stuff)

---
 .../polyhedralRefinement/polyhedralRefinement.C    |  9 ---------
 .../prismatic2DRefinement/prismatic2DRefinement.C  |  9 ---------
 .../refinement/refinement/refinement.C             | 14 --------------
 3 files changed, 32 deletions(-)

diff --git a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/polyhedralRefinement/polyhedralRefinement.C b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/polyhedralRefinement/polyhedralRefinement.C
index 0d5f2db57..64862038f 100644
--- a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/polyhedralRefinement/polyhedralRefinement.C
+++ b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/polyhedralRefinement/polyhedralRefinement.C
@@ -2055,15 +2055,6 @@ void Foam::polyhedralRefinement::setCellsToRefine
         }
     }
 
-    // Remove all cells that would exceed the maximum refinement level
-    forAll (refineCell, cellI)
-    {
-        if (refineCell[cellI] && (cellLevel_[cellI] + 1 > maxRefinementLevel_))
-        {
-            refineCell[cellI] = false;
-        }
-    }
-
     // Make sure that the refinement is face consistent (2:1 consistency) and
     // point consistent (4:1 consistency) if necessary
 
diff --git a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/prismatic2DRefinement/prismatic2DRefinement.C b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/prismatic2DRefinement/prismatic2DRefinement.C
index 7fa88a836..e6e3677d5 100644
--- a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/prismatic2DRefinement/prismatic2DRefinement.C
+++ b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/prismatic2DRefinement/prismatic2DRefinement.C
@@ -2614,15 +2614,6 @@ void Foam::prismatic2DRefinement::setCellsToRefine
         }
     }
 
-    // Remove all cells that exceed the maximum refinement level
-    forAll (refineCell, cellI)
-    {
-        if (refineCell[cellI] && (cellLevel_[cellI] + 1 > maxRefinementLevel_))
-        {
-            refineCell[cellI] = false;
-        }
-    }
-
     // Make sure that the refinement is face consistent (2:1 consistency) and
     // point consistent (4:1 consistency) if necessary
 
diff --git a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/refinement/refinement.C b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/refinement/refinement.C
index 63e1dd1fe..43f709ef5 100644
--- a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/refinement/refinement.C
+++ b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/refinement/refinement.C
@@ -1333,20 +1333,6 @@ bool Foam::refinement::changeTopology() const
 
 void Foam::refinement::setRefinement(polyTopoChange& ref) const
 {
-    // Make sure that the point levels are updated across coupled patches before
-    // setting refinement and unrefinement. Note: not sure why the sync is not
-    // performed correctly if I do it in updateMesh. This is a temporary
-    // solution, need to investigate in detail, but I assume something is not
-    // updated yet in that case. VV, 31/Jan/2018.
-    syncTools::syncPointList
-    (
-        mesh_,
-        pointLevel_,
-        maxEqOp<label>(),
-        label(0),   // Null value
-        true        // Apply separation for parallel cyclics
-    );
-
     // Set refinement and unrefinement
     this->setRefinementInstruction(ref);
     this->setUnrefinementInstruction(ref);

From d5948f3fb76d699b4fe91676a58489fa1b275fed Mon Sep 17 00:00:00 2001
From: Vuko Vukcevic <v.vukcevic@wikki.co.uk>
Date: Fri, 31 May 2019 10:39:45 +0200
Subject: [PATCH 04/11] Start/restart handling in dynamicPolyRefinementFvMesh

If the constant/polyMesh/meshModifiers file is present, we simply neglect it and
use the controls from dynamicMeshDict to avoid confusion.
---
 .../dynamicPolyRefinementFvMesh.C             | 174 ++++++++++--------
 1 file changed, 95 insertions(+), 79 deletions(-)

diff --git a/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C b/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C
index 680e7186b..4013e0681 100644
--- a/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C
+++ b/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C
@@ -117,89 +117,105 @@ Foam::dynamicPolyRefinementFvMesh::dynamicPolyRefinementFvMesh
     // protectedInitialRefinement)
     refinementSelectionPtr_()
 {
-    // Only create topo modifiers if they haven't been read in
-    // HJ, 16/Oct/2018
-    if (topoChanger_.empty())
+    // Check whether we read polyMeshModifiers from
+    // constant/polyMesh/meshModifiers file in the base class
+    if (!topoChanger_.empty())
     {
-        // Only one topo changer engine
-        topoChanger_.setSize(1);
+        // Already initialized, warn the user that we'll neglect it
+        WarningIn
+        (
+            "dynamicPolyRefinementFvMesh::dynamicPolyRefinementFvMesh"
+            "\n("
+            "\n    const IOobject& io"
+            "\n    const word subDictName"
+            "\n)"
+        )  << "Using controls from constant/dynamicMeshDict instead of"
+           << " constant/polyMesh/meshModifiers."
+           << nl
+           << "To supress this warning, delete meshModifiers file."
+           << endl;
 
-        // Get number of valid geometric dimensions
-        const label nGeometricDirs = this->nGeometricD();
-
-        switch(nGeometricDirs)
-        {
-            case 3:
-                // Add the polyhedralRefinement engine for
-                // 3D isotropic refinement
-                Info<< "3D case detected. "
-                    << "Adding polyhedralRefinement topology modifier" << endl;
-                topoChanger_.set
-                (
-                    0,
-                    new polyhedralRefinement
-                    (
-                        "polyhedralRefinement",
-                        refinementDict_,
-                        0,
-                        topoChanger_
-                    )
-                );
-                break;
-
-            case 2:
-                // Add the prismatic2DRefinement engine for
-                // 2D isotropic refinement
-                Info<< "2D case detected. "
-                    << "Adding prismatic2DRefinement topology modifier" << endl;
-                topoChanger_.set
-                (
-                    0,
-                    new prismatic2DRefinement
-                    (
-                        "prismatic2DRefinement",
-                        refinementDict_,
-                        0,
-                        topoChanger_
-                    )
-                );
-                break;
-
-            case 1:
-                FatalErrorIn
-                (
-                    "dynamicPolyRefinementFvMesh::dynamicPolyRefinementFvMesh"
-                    "\n("
-                    "\n    const IOobject& io,"
-                    "\n    const word subDictName"
-                    "\n)"
-                )   << "1D case detected. No valid refinement strategy is"
-                    <<  " available for 1D cases."
-                    << abort(FatalError);
-                break;
-
-            default:
-                FatalErrorIn
-                (
-                    "dynamicPolyRefinementFvMesh::dynamicPolyRefinementFvMesh"
-                    "\n("
-                    "\n    const IOobject& io,"
-                    "\n    const word subDictName"
-                    "\n)"
-                )   << "Invalid number of geometric meshes detected: "
-                    << nGeometricDirs
-                    << nl << "It appears that this mesh is neither 1D, 2D or 3D."
-                    << nl << " the mesh."
-                    << abort(FatalError);
-
-        }
-
-        // Write mesh and modifiers
-        topoChanger_.writeOpt() = IOobject::AUTO_WRITE;
-        topoChanger_.write();
-        write();
+       // Clear the list
+       topoChanger_.clear();
     }
 
+    // Only one topo changer engine
+    topoChanger_.setSize(1);
+
+    // Get number of valid geometric dimensions
+    const label nGeometricDirs = this->nGeometricD();
+
+    switch(nGeometricDirs)
+    {
+        case 3:
+            // Add the polyhedralRefinement engine for
+            // 3D isotropic refinement
+            Info<< "3D case detected. "
+                << "Adding polyhedralRefinement topology modifier" << endl;
+            topoChanger_.set
+            (
+                0,
+                new polyhedralRefinement
+                (
+                    "polyhedralRefinement",
+                    refinementDict_,
+                    0,
+                    topoChanger_
+                )
+            );
+            break;
+
+        case 2:
+            // Add the prismatic2DRefinement engine for
+            // 2D isotropic refinement
+            Info<< "2D case detected. "
+                << "Adding prismatic2DRefinement topology modifier" << endl;
+            topoChanger_.set
+            (
+                0,
+                new prismatic2DRefinement
+                (
+                    "prismatic2DRefinement",
+                    refinementDict_,
+                    0,
+                    topoChanger_
+                )
+            );
+            break;
+
+        case 1:
+            FatalErrorIn
+            (
+                "dynamicPolyRefinementFvMesh::dynamicPolyRefinementFvMesh"
+                "\n("
+                "\n    const IOobject& io,"
+                "\n    const word subDictName"
+                "\n)"
+            )   << "1D case detected. No valid refinement strategy is"
+                <<  " available for 1D cases."
+                << abort(FatalError);
+            break;
+
+        default:
+            FatalErrorIn
+            (
+                "dynamicPolyRefinementFvMesh::dynamicPolyRefinementFvMesh"
+                "\n("
+                "\n    const IOobject& io,"
+                "\n    const word subDictName"
+                "\n)"
+            )   << "Invalid number of geometric meshes detected: "
+                << nGeometricDirs
+                << nl << "It appears that this mesh is neither 1D, 2D or 3D."
+                << abort(FatalError);
+
+    }
+
+    // Write mesh and modifiers
+    topoChanger_.writeOpt() = IOobject::AUTO_WRITE;
+    topoChanger_.write();
+    write();
+
     // Initialize refinement selection algorithm after modifiers
     refinementSelectionPtr_ = refinementSelection::New(*this, refinementDict_);
 }

From fd18eaffa13c2f8d83ec7aebd296eedff10a731d Mon Sep 17 00:00:00 2001
From: Vuko Vukcevic <v.vukcevic@wikki.co.uk>
Date: Mon, 3 Jun 2019 09:36:41 +0200
Subject: [PATCH 05/11] Temporarily issuing an error if we have zones while
 doing DLB

---
 .../finiteVolume/processorMeshesRebuild.C     | 29 ++++++++++++++++++-
 1 file changed, 28 insertions(+), 1 deletion(-)

diff --git a/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesRebuild.C b/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesRebuild.C
index df1afae56..d1373341f 100644
--- a/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesRebuild.C
+++ b/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesRebuild.C
@@ -1728,7 +1728,34 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
     // due to the presence of old/new patches
     globalMesh.addFvPatches(reconPatches, false);
 
-    // TODO: point, face and cell zones
+    // TODO: point, face and cell zone
+    // Temporarily issue an error if we have point, face and cell zones
+    forAll (meshes_, procI)
+    {
+        if (meshes_.set(procI))
+        {
+            const polyMesh& curMesh = meshes_[procI];
+
+            if
+            (
+                !curMesh.cellZones().empty()
+             || !curMesh.faceZones().empty()
+             || !curMesh.pointZones().empty()
+            )
+            {
+                FatalErrorIn
+                (
+                    "autoPtr<fvMesh> processorMeshesReconstructor::"
+                    "reconstructMesh(const Time& db)"
+                )   << "Reconstructing cellZones, faceZones and pointZones is"
+                    << " currently not supported."
+                    << nl
+                    << "In order to get past this error, you can delete all the"
+                    << " zones in your mesh."
+                    << abort(FatalError);
+            }
+        }
+    }
 
     return globalMeshPtr;
 }

From 13f349b7691bea536392f615c3c89f6f589db96c Mon Sep 17 00:00:00 2001
From: Vuko Vukcevic <v.vukcevic@wikki.co.uk>
Date: Thu, 6 Jun 2019 13:13:33 +0200
Subject: [PATCH 06/11] Added empty and wedge direction tolerances in polyMesh

On same DLB+AMR cases, the round-off error caused the difference in global and
local empty direction to be slightly larger than SMALL.
---
 src/foam/meshes/polyMesh/polyMesh.C | 41 ++++++++++++++++++++++++-----
 src/foam/meshes/polyMesh/polyMesh.H | 20 +++++++++++---
 2 files changed, 51 insertions(+), 10 deletions(-)

diff --git a/src/foam/meshes/polyMesh/polyMesh.C b/src/foam/meshes/polyMesh/polyMesh.C
index 27eee495a..38a004391 100644
--- a/src/foam/meshes/polyMesh/polyMesh.C
+++ b/src/foam/meshes/polyMesh/polyMesh.C
@@ -43,6 +43,20 @@ defineTypeNameAndDebug(Foam::polyMesh, 0);
 Foam::word Foam::polyMesh::defaultRegion = "region0";
 Foam::word Foam::polyMesh::meshSubDir = "polyMesh";
 
+const Foam::debug::tolerancesSwitch
+Foam::polyMesh::emptyDirTol_
+(
+    "emptyDirectionTolerance",
+    1e-10
+);
+
+const Foam::debug::tolerancesSwitch
+Foam::polyMesh::wedgeDirTol_
+(
+    "wedgeDirectionTolerance",
+    1e-10
+);
+
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -93,20 +107,27 @@ void Foam::polyMesh::calcDirections() const
         globalEmptyDirVec /= mag(globalEmptyDirVec);
 
         // Check if all processors see the same 2-D from empty patches
-        if (mag(globalEmptyDirVec - emptyDirVec) > SMALL)
+        if (mag(globalEmptyDirVec - emptyDirVec) > emptyDirTol_())
         {
             FatalErrorIn
             (
                 "void polyMesh::calcDirections() const"
             )   << "Some processors detect different empty (2-D) "
                 << "directions.  Probably using empty patches on a "
-                << "bad parallel decomposition." << nl
-                << "Please check your geometry and empty patches"
+                << "bad parallel decomposition." << nl << nl
+                << "Global empty direction vector: " << globalEmptyDirVec << nl
+                << "Local empty direction vector: " << emptyDirVec << nl
+                << "If global and local empty direction vectors are different"
+                << " to a round-off error, try increasing the"
+                << " emptyDirectionTolerance in controlDict::Tolerances" << nl
+                << "The emptyDirectionTolerance for this run is: "
+                << emptyDirTol_() << nl << nl
+                << "Otherwise please check your geometry and empty patches."
                 << abort(FatalError);
         }
     }
 
-    if (mag(emptyDirVec) > SMALL)
+    if (mag(emptyDirVec) > emptyDirTol_())
     {
         for (direction cmpt = 0; cmpt < vector::nComponents; cmpt++)
         {
@@ -139,7 +160,7 @@ void Foam::polyMesh::calcDirections() const
         globalWedgeDirVec /= mag(globalWedgeDirVec);
 
         // Check if all processors see the same 2-D from wedge patches
-        if (mag(globalWedgeDirVec - wedgeDirVec) > SMALL)
+        if (mag(globalWedgeDirVec - wedgeDirVec) > wedgeDirTol_())
         {
             FatalErrorIn
             (
@@ -147,12 +168,20 @@ void Foam::polyMesh::calcDirections() const
             )   << "Some processors detect different wedge (2-D) "
                 << "directions.  Probably using wedge patches on a "
                 << "bad parallel decomposition." << nl
+                << "Global wedge direction vector: " << globalWedgeDirVec << nl
+                << "Local wedge direction vector: " << wedgeDirVec << nl
+                << "If global and local wedge direction vectors are different"
+                << " to a round-off error, try increasing the"
+                << " wedgeDirectionTolerance in controlDict::Tolerances" << nl
+                << "The wedgeDirectionTolerance for this run is: "
+                << wedgeDirTol_() << nl << nl
+                << "Otherwise please check your geometry and empty patches."
                 << "Please check your geometry and wedge patches"
                 << abort(FatalError);
         }
     }
 
-    if (mag(wedgeDirVec) > SMALL)
+    if (mag(wedgeDirVec) > wedgeDirTol_())
     {
         for (direction cmpt = 0; cmpt < vector::nComponents; cmpt++)
         {
diff --git a/src/foam/meshes/polyMesh/polyMesh.H b/src/foam/meshes/polyMesh/polyMesh.H
index f15acc99c..27a6e461b 100644
--- a/src/foam/meshes/polyMesh/polyMesh.H
+++ b/src/foam/meshes/polyMesh/polyMesh.H
@@ -235,11 +235,23 @@ public:
     //- Runtime type information
     TypeName("polyMesh");
 
-    //- Return the default region name
-    static word defaultRegion;
 
-    //- Return the mesh sub-directory name (usually "polyMesh")
-    static word meshSubDir;
+    //- Static mesh data
+
+        //- Return the default region name
+        static word defaultRegion;
+
+        //- Return the mesh sub-directory name (usually "polyMesh")
+        static word meshSubDir;
+
+
+    //- Static data to control empty and wedge directions
+
+        //- Empty direction tolerance
+        static const debug::tolerancesSwitch emptyDirTol_;
+
+        //- Wedge direction tolerance
+        static const debug::tolerancesSwitch wedgeDirTol_;
 
 
     // Constructors

From d9cfdea84f6c8c520030d098f31d69d833713c1e Mon Sep 17 00:00:00 2001
From: Vuko Vukcevic <v.vukcevic@wikki.co.uk>
Date: Fri, 7 Jun 2019 08:36:46 +0200
Subject: [PATCH 07/11] Reverted to inv in fvcReconstruct

---
 src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C b/src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C
index ed2e8fcab..c39989e68 100644
--- a/src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C
+++ b/src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C
@@ -90,10 +90,9 @@ reconstruct
     // Note: hinv inverse must be used to stabilise the inverse on bad meshes
     // but it gives strange failures.  Fixed hinv.  HJ, 22/Mar/2019
     // HJ, 19/Aug/2015
-    // Temporarily reverting hinv: Vuko Vukcevic, 6/Jun/2019
+    // Nope, still problematic, VV, 7/Jun/2019
     reconField.internalField() =
     (
-        // hinv
         inv
         (
             surfaceSum(sqr(mesh.Sf())/mesh.magSf())().internalField()

From 57d3af8960c4c82ebb1c3e9e5e982201a81bac6b Mon Sep 17 00:00:00 2001
From: Vuko Vukcevic <v.vukcevic@wikki.co.uk>
Date: Fri, 14 Jun 2019 14:36:09 +0200
Subject: [PATCH 08/11] Bugfix on creating PrimitivePatch::meshPointsPtr_

We must not enforce the order of points on processor patches when collecting
patch points from faces because this occasionally causes ordering errors on two
sides (e.g. in DLB).
---
 .../PrimitivePatchMeshData.C                  | 43 ++++---------------
 1 file changed, 9 insertions(+), 34 deletions(-)

diff --git a/src/foam/meshes/primitiveMesh/PrimitivePatchTemplate/PrimitivePatchMeshData.C b/src/foam/meshes/primitiveMesh/PrimitivePatchTemplate/PrimitivePatchMeshData.C
index ed7d4364e..b9ad4a1d5 100644
--- a/src/foam/meshes/primitiveMesh/PrimitivePatchTemplate/PrimitivePatchMeshData.C
+++ b/src/foam/meshes/primitiveMesh/PrimitivePatchTemplate/PrimitivePatchMeshData.C
@@ -73,53 +73,28 @@ calcMeshData() const
     // number of faces in the patch
     Map<label> markedPoints(4*this->size());
 
-
     // Important:
     // ~~~~~~~~~~
     // In <= 1.5 the meshPoints would be in increasing order but this gives
     // problems in processor point synchronisation where we have to find out
-    // how the opposite side would have allocated points.
-
-    // Note:
-    // ~~~~~
-    // This is all garbage.  All -ext versions will preserve strong ordering
-    // HJ, 17/Aug/2010
-
-    //- 1.5 code:
-    // If the point is used, set the mark to 1
+    // how the opposite side would have allocated points. We'll use unsorted
+    // version to avoid certain ordering problems
+    dynamicLabelList meshPoints(2*this->size());
     forAll(*this, facei)
     {
         const Face& curPoints = this->operator[](facei);
 
         forAll(curPoints, pointi)
         {
-            markedPoints.insert(curPoints[pointi], -1);
+            if (markedPoints.insert(curPoints[pointi], meshPoints.size()))
+            {
+                meshPoints.append(curPoints[pointi]);
+            }
         }
     }
 
-    // Create the storage and store the meshPoints.  Mesh points are
-    // the ones marked by the usage loop above
-    meshPointsPtr_ = new labelList(markedPoints.toc());
-    labelList& pointPatch = *meshPointsPtr_;
-
-    // Sort the list to preserve compatibility with the old ordering
-    sort(pointPatch);
-
-    // For every point in map give it its label in mesh points
-    forAll(pointPatch, pointi)
-    {
-        markedPoints.find(pointPatch[pointi])() = pointi;
-    }
-
-    forAll(*this, faceI)
-    {
-        const Face& curPoints = this->operator[](faceI);
-
-        forAll (curPoints, pointI)
-        {
-            markedPoints.insert(curPoints[pointI], -1);
-        }
-    }
+    // Transfer to straight list (reuses storage)
+    meshPointsPtr_ = new labelList(meshPoints, true);
 
     // Create local faces. Note that we start off from copy of original face
     // list (even though vertices are overwritten below). This is done so

From 49769f03979928522fcc5a5c9f0c55a2112128b8 Mon Sep 17 00:00:00 2001
From: Vuko Vukcevic <v.vukcevic@wikki.co.uk>
Date: Wed, 19 Jun 2019 16:32:10 +0200
Subject: [PATCH 09/11] Reconstruction of cell/face/point zones during DLB

Wrote a missing piece of code that maps cell/face/point zones from different
mesh bits to the new, reconstructed mesh when doing Dynamic Load Balancing.
---
 .../finiteVolume/processorMeshesRebuild.C     | 335 +++++++++++++++++-
 .../topoChangerFvMeshLoadBalance.C            |  52 +++
 2 files changed, 368 insertions(+), 19 deletions(-)

diff --git a/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesRebuild.C b/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesRebuild.C
index d1373341f..d5592f624 100644
--- a/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesRebuild.C
+++ b/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesRebuild.C
@@ -1728,35 +1728,332 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
     // due to the presence of old/new patches
     globalMesh.addFvPatches(reconPatches, false);
 
-    // TODO: point, face and cell zone
-    // Temporarily issue an error if we have point, face and cell zones
+    // Recombine cell, face and point zones.
+    // Note 1: all zones have to be present on same processors in the same
+    // order. This is the result of the decomposition. See
+    // domainDecomposition::processorMesh member function
+    // Note 2: the code below could be written in a generic way by using a
+    // template helper member function, but it's not straightforward since we
+    // don't have a list of ZoneMeshes for all processors
+
+    // Get index for the first valid mesh
+    const label fvmID = firstValidMesh();
+
+    // First pass: count maximum number of cells, faces and points in zones
+    labelList nCellsPerZone(meshes_[fvmID].cellZones().size(), 0);
+    labelList nFacesPerZone(meshes_[fvmID].faceZones().size(), 0);
+    labelList nPointsPerZone(meshes_[fvmID].pointZones().size(), 0);
+
     forAll (meshes_, procI)
     {
         if (meshes_.set(procI))
         {
+            // Grab the current mesh
             const polyMesh& curMesh = meshes_[procI];
 
-            if
-            (
-                !curMesh.cellZones().empty()
-             || !curMesh.faceZones().empty()
-             || !curMesh.pointZones().empty()
-            )
+            // PART 1: Cell zones. Scope for clarity and safety
             {
-                FatalErrorIn
-                (
-                    "autoPtr<fvMesh> processorMeshesReconstructor::"
-                    "reconstructMesh(const Time& db)"
-                )   << "Reconstructing cellZones, faceZones and pointZones is"
-                    << " currently not supported."
-                    << nl
-                    << "In order to get past this error, you can delete all the"
-                    << " zones in your mesh."
-                    << abort(FatalError);
+                const cellZoneMesh& cz = curMesh.cellZones();
+
+                forAll (cz, zoneI)
+                {
+                    // Count number of cells in the zone
+                    nCellsPerZone[zoneI] += cz[zoneI].size();
+                }
             }
-        }
+
+            // PART 2: Face zones. Scope for clarity and safety
+            {
+                const faceZoneMesh& fz = curMesh.faceZones();
+
+                forAll (fz, zoneI)
+                {
+                    // Count number of faces in the zone
+                    nFacesPerZone[zoneI] += fz[zoneI].size();
+                }
+            }
+
+            // PART 3: Point zones. Scope for clarity and safety
+            {
+                const pointZoneMesh& pz = curMesh.pointZones();
+
+                forAll (pz, zoneI)
+                {
+                    // Count number of points in the zone
+                    nPointsPerZone[zoneI] += pz[zoneI].size();
+                }
+            }
+        } // End if processor mesh set
+    } // End for all processor meshes
+
+    // Second pass: redistribute cells, faces and points in zones
+
+    // Create lists that contain labels of all cells/faces/points in a given
+    // zone coming from different processor meshes. Index is the zoneID, which
+    // is the same for all processor bits, while the other list collects all the
+    // cells/faces/points in a given zone.
+    labelListList reconCellZones(nCellsPerZone.size());
+    labelListList reconFaceZones(nFacesPerZone.size());
+    List<boolList> reconFaceZoneFlips(nFacesPerZone.size());
+    labelListList reconPointZones(nPointsPerZone.size());
+
+    // Size the lists appropriatelly for each zone
+    forAll (reconCellZones, zoneI)
+    {
+        reconCellZones[zoneI].setSize(nCellsPerZone[zoneI], -1);
+    }
+    forAll (reconFaceZones, zoneI)
+    {
+        reconFaceZones[zoneI].setSize(nFacesPerZone[zoneI], -1);
+        reconFaceZoneFlips[zoneI].setSize(nFacesPerZone[zoneI], false);
+    }
+    forAll (reconPointZones, zoneI)
+    {
+        reconPointZones[zoneI].setSize(nPointsPerZone[zoneI], -1);
     }
 
+    // Reset counting lists for indexing during list item assignement and for
+    // collecting the final number of faces/points in face/pointZones (since
+    // these can be actually fewer if e.g. a processor face becomes an internal
+    // face).
+    nCellsPerZone = 0;
+    nFacesPerZone = 0;
+    nPointsPerZone = 0;
+
+    // Loop through all the meshes and collect cells/faces/points in the zones
+    forAll (meshes_, procI)
+    {
+        if (meshes_.set(procI))
+        {
+            // Grab the current mesh
+            const polyMesh& curMesh = meshes_[procI];
+
+            // PART 1: Cell zones. Scope for clarity and safety
+            {
+                const cellZoneMesh& cz = curMesh.cellZones();
+
+                // Get old-to-new cell addressing for this mesh
+                const labelList& curCellProcAddr = cellProcAddressing_[procI];
+
+                forAll (cz, zoneI)
+                {
+                    // Get "new" zone cell index
+                    label& nCells = nCellsPerZone[zoneI];
+
+                    // Reference to the new recon zone
+                    labelList& zoneReconCells = reconCellZones[zoneI];
+
+                    // Get all the cells in this zone
+                    const labelList& zoneCells = cz[zoneI];
+
+                    // Loop through the cells
+                    forAll (zoneCells, i)
+                    {
+                        // Get cell index in the processor mesh
+                        const label& oldCellI = zoneCells[i];
+
+                        // Get cell index in the new mesh
+                        const label& newCellI = curCellProcAddr[oldCellI];
+
+                        // Redundancy check: check if the the newCellI is
+                        // -1. This should not happen because cells have perfect
+                        // 1-to-1 mapping
+                        if (newCellI != -1)
+                        {
+                            // Insert the cell in the new recon zone and
+                            // increment the counter
+                            zoneReconCells[nCells++] = newCellI;
+                        }
+                        else
+                        {
+                            FatalErrorIn
+                            (
+                                "autoPtr<fvMesh>"
+                                "\n processorMeshesReconstructor::"
+                                "reconstructMesh(const Time& db)"
+                            )   << "Found unmapped cell while reconstructing"
+                                << " cell zones."
+                                << nl
+                                << "Cell from processor: " << procI << nl
+                                << "Cell zone name: " << cz[zoneI].name() << nl
+                                << "Index in the cell zone: " << i << nl
+                                << "Old cell index: " << oldCellI << nl
+                                << "New cell index: " << newCellI
+                                << abort(FatalError);
+                        }
+
+                    } // End for all cells in the zone
+                } // End for all cell zones
+            } // End scope for cell zone handling
+
+            // PART 2: Face zones. Scope for clarity and safety
+            {
+                const faceZoneMesh& fz = curMesh.faceZones();
+
+                // Get old-to-new face addressing for this mesh
+                const labelList& curFaceProcAddr = faceProcAddressing_[procI];
+
+                forAll (fz, zoneI)
+                {
+                    // Get "new" zone face index
+                    label& nFaces = nFacesPerZone[zoneI];
+
+                    // Reference to the new recon zone and flips in the zone
+                    labelList& zoneReconFaces = reconFaceZones[zoneI];
+                    boolList& zoneReconFaceFlips = reconFaceZoneFlips[zoneI];
+
+                    // Get all the faces in this zone and also their flips
+                    const labelList& zoneFaces = fz[zoneI];
+                    const boolList& zoneFlipMap = fz[zoneI].flipMap();
+
+                    // Loop through the faces
+                    forAll (zoneFaces, i)
+                    {
+                        // Get the face index in the processor mesh
+                        const label& oldFaceI = zoneFaces[i];
+
+                        // Get the face index in the new, reconstructed mesh
+                        const label& newFaceI = curFaceProcAddr[oldFaceI];
+
+                        // Check if the face is mapped (if the index is >= 0)
+                        if (newFaceI > -1)
+                        {
+                            // This is a face that's been correctly
+                            // mapped, insert the face in the new zone
+                            zoneReconFaces[nFaces] = newFaceI;
+
+                            // Also store the flip map of the face. Note: I'm
+                            // pretty sure that we don't need to check whether
+                            // the flip map has been preserved because the
+                            // combined faces are inserted from master side
+                            // always.
+                            zoneReconFaceFlips[nFaces] = zoneFlipMap[i];
+
+                            // Increment the number of faces for this zone
+                            ++nFaces;
+                        }
+                    } // End for all faces in the zone
+                } // End for all face zones
+            } // End scope for face zone handling
+
+            // PART 3: Point zones. Scope for clarity and safety
+            {
+                const pointZoneMesh& fz = curMesh.pointZones();
+
+                // Get old-to-new point addressing for this mesh
+                const labelList& curPointProcAddr = pointProcAddressing_[procI];
+
+                forAll (fz, zoneI)
+                {
+                    // Get "new" zone point index
+                    label& nPoints = nPointsPerZone[zoneI];
+
+                    // Reference to the new recon zone
+                    labelList& zoneReconPoints = reconPointZones[zoneI];
+
+                    // Get all the points in this zone
+                    const labelList& zonePoints = fz[zoneI];
+
+                    // Loop through the points
+                    forAll (zonePoints, i)
+                    {
+                        // Get point index in the processor mesh
+                        const label& oldPointI = zonePoints[i];
+
+                        // Get point index in the new mesh
+                        const label& newPointI = curPointProcAddr[oldPointI];
+
+                        // Check if the point is mapped
+                        if (newPointI != -1)
+                        {
+                            // Insert the point in the new recon zone and
+                            // increment the counter
+                            zoneReconPoints[nPoints++] = newPointI;
+                        }
+
+                    } // End for all points in the zone
+                } // End for all point zones
+            } // End scope for point zone handling
+        } // End if the processor mesh is set
+    } // End for all processor meshes
+
+    // We need to resize the face and point zones to number of inserted
+    // faces/points because not all faces and points need to be
+    // inserted. There's nothing to do for cell zones because these are always
+    // mapped uniquely one-to-one
+    forAll (reconFaceZones, zoneI)
+    {
+        reconFaceZones[zoneI].setSize(nFacesPerZone[zoneI]);
+        reconFaceZoneFlips[zoneI].setSize(nFacesPerZone[zoneI]);
+    }
+    forAll (reconPointZones, zoneI)
+    {
+        reconPointZones[zoneI].setSize(nPointsPerZone[zoneI]);
+    }
+
+    // Now we have all the zones as ordinary lists without possible duplicate
+    // faces and points due to merging of processor boundaries. Create zone
+    // meshes
+
+    // PART 1: Cell zones
+    List<cellZone*> reconCz(reconCellZones.size());
+
+    // Loop through all the cell zones and create them
+    forAll (reconCz, zoneI)
+    {
+        // Notes:
+        // 1. Grab the name from the respective zone in the first valid mesh
+        // 2. Transfer the list of cell IDs, invalidating reconCellZones[zoneI]
+        reconCz[zoneI] = new cellZone
+        (
+            meshes_[fvmID].cellZones()[zoneI].name(),
+            reconCellZones[zoneI].xfer(),
+            zoneI,
+            globalMesh.cellZones()
+        );
+    }
+
+    // PART 2: Face zones
+    List<faceZone*> reconFz(reconFaceZones.size());
+
+    // Loop through all the face zones and create them
+    forAll (reconFz, zoneI)
+    {
+        // Notes:
+        // 1. Grab the name from the respective zone in the first valid mesh
+        // 2. Transfer the list of face IDs, invalidating reconFaceZones[zoneI]
+        reconFz[zoneI] = new faceZone
+        (
+            meshes_[fvmID].faceZones()[zoneI].name(),
+            reconFaceZones[zoneI].xfer(),
+            reconFaceZoneFlips[zoneI].xfer(),
+            zoneI,
+            globalMesh.faceZones()
+        );
+    }
+
+    // PART 3: Point zones
+    List<pointZone*> reconPz(reconPointZones.size());
+
+    // Loop through all the point zones and create them
+    forAll (reconPz, zoneI)
+    {
+        // Notes:
+        // 1. Grab the name from the respective zone in the first valid mesh
+        // 2. Transfer the list of point IDs, invalidating reconPointZones[zoneI]
+        reconPz[zoneI] = new pointZone
+        (
+            meshes_[fvmID].pointZones()[zoneI].name(),
+            reconPointZones[zoneI].xfer(),
+            zoneI,
+            globalMesh.pointZones()
+        );
+    }
+
+    // Add the zones into the mesh
+    globalMesh.addZones(reconPz, reconFz, reconCz);
+
+    // All done, return the global mesh pointer
     return globalMeshPtr;
 }
 
diff --git a/src/dynamicMesh/topoChangerFvMesh/topoChangerFvMesh/topoChangerFvMeshLoadBalance.C b/src/dynamicMesh/topoChangerFvMesh/topoChangerFvMesh/topoChangerFvMeshLoadBalance.C
index 8285623ee..ed95986d0 100644
--- a/src/dynamicMesh/topoChangerFvMesh/topoChangerFvMesh/topoChangerFvMeshLoadBalance.C
+++ b/src/dynamicMesh/topoChangerFvMesh/topoChangerFvMesh/topoChangerFvMeshLoadBalance.C
@@ -844,6 +844,9 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
         oldPatchNMeshPoints         // oldPatchNMeshPoints
     );
 
+    // Remove point, face and cell zones from the original mesh
+    removeZones();
+
     // Reset fvMesh and patches
     resetFvPrimitives
     (
@@ -857,6 +860,55 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
         true                       // Valid boundary
     );
 
+    // Get pointers to cell/face/point zones from reconstructed mesh and "link"
+    // them to the new mesh
+
+    // Cell zones
+    const cellZoneMesh& reconCellZones = reconMesh.cellZones();
+    List<cellZone*> czs(reconCellZones.size());
+    forAll (czs, zoneI)
+    {
+        czs[zoneI] = new cellZone
+        (
+            reconCellZones[zoneI].name(),
+            reconCellZones[zoneI],
+            zoneI,
+            this->cellZones()
+        );
+    }
+
+    // Face zones
+    const faceZoneMesh& reconFaceZones = reconMesh.faceZones();
+    List<faceZone*> fzs(reconFaceZones.size());
+    forAll (fzs, zoneI)
+    {
+        fzs[zoneI] = new faceZone
+        (
+            reconFaceZones[zoneI].name(),
+            reconFaceZones[zoneI],
+            reconFaceZones[zoneI].flipMap(),
+            zoneI,
+            this->faceZones()
+        );
+    }
+
+    // Point zones
+    const pointZoneMesh& reconPointZones = reconMesh.pointZones();
+    List<pointZone*> pzs(reconPointZones.size());
+    forAll (pzs, zoneI)
+    {
+        pzs[zoneI] = new pointZone
+        (
+            reconPointZones[zoneI].name(),
+            reconPointZones[zoneI],
+            zoneI,
+            this->pointZones()
+        );
+    }
+
+    // Add the zones to the mesh
+    addZones(pzs, fzs, czs);
+
     // Create field reconstructor
     fvFieldReconstructor fieldReconstructor
     (

From 15af49e74a3a0340835d1ad7451bca437c2bd766 Mon Sep 17 00:00:00 2001
From: Vuko Vukcevic <v.vukcevic@wikki.co.uk>
Date: Mon, 24 Jun 2019 09:51:40 +0200
Subject: [PATCH 10/11] Bugfix in faceZone reconstruction procedure

Forgot to take into account the "turning face" index and decrement by one
---
 .../finiteVolume/processorMeshesRebuild.C     | 21 ++++++++++++-------
 1 file changed, 13 insertions(+), 8 deletions(-)

diff --git a/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesRebuild.C b/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesRebuild.C
index d5592f624..19fd200fa 100644
--- a/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesRebuild.C
+++ b/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesRebuild.C
@@ -1915,18 +1915,23 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
                         // Get the face index in the new, reconstructed mesh
                         const label& newFaceI = curFaceProcAddr[oldFaceI];
 
-                        // Check if the face is mapped (if the index is >= 0)
-                        if (newFaceI > -1)
+                        // Check if the face is mapped.
+                        // Note:
+                        // 1. Need to decrement by 1 because of the face turning
+                        // 2. No need to handle negative new indices coming from
+                        //    slave processor because we'd end up with
+                        //    duplicate entries (two faces on two processors
+                        //    merged into a single one)
+                        if (newFaceI > 0)
                         {
                             // This is a face that's been correctly
                             // mapped, insert the face in the new zone
-                            zoneReconFaces[nFaces] = newFaceI;
+                            zoneReconFaces[nFaces] = newFaceI - 1;
 
-                            // Also store the flip map of the face. Note: I'm
-                            // pretty sure that we don't need to check whether
-                            // the flip map has been preserved because the
-                            // combined faces are inserted from master side
-                            // always.
+                            // Also store the flip map of the face. We don't
+                            // need to check whether the flip map has been
+                            // preserved because we only get the combined faces
+                            // that are inserted from master side.
                             zoneReconFaceFlips[nFaces] = zoneFlipMap[i];
 
                             // Increment the number of faces for this zone

From 907c80509df249d492303da488f4b3f4ea6c2439 Mon Sep 17 00:00:00 2001
From: Vuko Vukcevic <v.vukcevic@wikki.co.uk>
Date: Thu, 11 Jul 2019 12:44:40 +0200
Subject: [PATCH 11/11] Added comment in PrimitivePatchMeshData

Regarding usage of unsorted variant of meshPoints.
---
 .../PrimitivePatchTemplate/PrimitivePatchMeshData.C   | 11 +++++------
 1 file changed, 5 insertions(+), 6 deletions(-)

diff --git a/src/foam/meshes/primitiveMesh/PrimitivePatchTemplate/PrimitivePatchMeshData.C b/src/foam/meshes/primitiveMesh/PrimitivePatchTemplate/PrimitivePatchMeshData.C
index b9ad4a1d5..8993cf124 100644
--- a/src/foam/meshes/primitiveMesh/PrimitivePatchTemplate/PrimitivePatchMeshData.C
+++ b/src/foam/meshes/primitiveMesh/PrimitivePatchTemplate/PrimitivePatchMeshData.C
@@ -73,12 +73,11 @@ calcMeshData() const
     // number of faces in the patch
     Map<label> markedPoints(4*this->size());
 
-    // Important:
-    // ~~~~~~~~~~
-    // In <= 1.5 the meshPoints would be in increasing order but this gives
-    // problems in processor point synchronisation where we have to find out
-    // how the opposite side would have allocated points. We'll use unsorted
-    // version to avoid certain ordering problems
+    // Note: foam-extend was relying on strong ordering of meshPoints on
+    // processor patches which causes point ordering errors when doing Dynamic
+    // Load Balancing. Because of this ordering error on two sides of processor
+    // boundaries, point fields ended up having wrong values near processor
+    // boundaries. Revert to unsorted version. VV, 11/July/2019.
     dynamicLabelList meshPoints(2*this->size());
     forAll(*this, facei)
     {