DLB and AMR fixes and additional developments. Vuko Vukcevic

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

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

@@ -1728,8 +1728,337 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
     // due to the presence of old/new patches
     globalMesh.addFvPatches(reconPatches, false);
 
-    // TODO: 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];
+
+            // PART 1: Cell zones. Scope for clarity and safety
+            {
+                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.
+                        // 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 - 1;
+
+                            // 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
+                            ++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;
 }
 

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
 

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
 

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

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
+    // Check whether we read polyMeshModifiers from
-    // HJ, 16/Oct/2018
+    // constant/polyMesh/meshModifiers file in the base class
-    if (topoChanger_.empty())
+    if (!topoChanger_.empty())
     {
-        // Only one topo changer engine
+        // Already initialized, warn the user that we'll neglect it
-        topoChanger_.setSize(1);
+        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
+       // Clear the list
-        const label nGeometricDirs = this->nGeometricD();
+       topoChanger_.clear();
-
-        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();
     }
 
+    // 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_);
 }
@@ -271,13 +287,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 +313,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 +360,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 +378,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;
 }

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
@@ -839,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
     (
@@ -852,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
     (

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
+                << "bad parallel decomposition." << nl << nl
-                << "Please check your geometry and empty patches"
+                << "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++)
         {

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 mesh data
-    static word meshSubDir;
+
+        //- 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

src/foam/meshes/primitiveMesh/PrimitivePatchTemplate/PrimitivePatchMeshData.C

@@ -73,53 +73,27 @@ calcMeshData() const
     // number of faces in the patch
     Map<label> markedPoints(4*this->size());
 
-
+    // Note: foam-extend was relying on strong ordering of meshPoints on
-    // Important:
+    // processor patches which causes point ordering errors when doing Dynamic
-    // ~~~~~~~~~~
+    // Load Balancing. Because of this ordering error on two sides of processor
-    // In <= 1.5 the meshPoints would be in increasing order but this gives
+    // boundaries, point fields ended up having wrong values near processor
-    // problems in processor point synchronisation where we have to find out
+    // boundaries. Revert to unsorted version. VV, 11/July/2019.
-    // how the opposite side would have allocated points.
+    dynamicLabelList meshPoints(2*this->size());
-
-    // 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
     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
+    // Transfer to straight list (reuses storage)
-    // the ones marked by the usage loop above
+    meshPointsPtr_ = new labelList(meshPoints, true);
-    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);
-        }
-    }
 
     // Create local faces. Note that we start off from copy of original face
     // list (even though vertices are overwritten below). This is done so