decomposePar: Speed-up decomposing with constraints

Patch contributed by Mattijs Janssens
Resolves bug-report http://bugs.openfoam.org/view.php?id=2159

src/meshTools/algorithms/MeshWave/FaceCellWave.C

@@ -49,10 +49,11 @@ int Foam::FaceCellWave<Type, TrackingData>::dummyTrackData_ = 12345;
 
 namespace Foam
 {
-    //- Combine operator for AMIInterpolation
     template<class Type, class TrackingData>
     class combine
     {
+        //- Combine operator for AMIInterpolation
+
         FaceCellWave<Type, TrackingData>& solver_;
 
         const cyclicAMIPolyPatch& patch_;
@@ -103,14 +104,8 @@ namespace Foam
 }
 
 
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-// Update info for celli, at position pt, with information from
-// neighbouring face/cell.
-// Updates:
-//      - changedCell_, changedCells_, nChangedCells_,
-//      - statistics: nEvals_, nUnvisitedCells_
 template<class Type, class TrackingData>
 bool Foam::FaceCellWave<Type, TrackingData>::updateCell
 (
@@ -121,6 +116,12 @@ bool Foam::FaceCellWave<Type, TrackingData>::updateCell
     Type& cellInfo
 )
 {
+    // Update info for celli, at position pt, with information from
+    // neighbouring face/cell.
+    // Updates:
+    //      - changedCell_, changedCells_
+    //      - statistics: nEvals_, nUnvisitedCells_
+
     nEvals_++;
 
     bool wasValid = cellInfo.valid(td_);
@@ -141,7 +142,7 @@ bool Foam::FaceCellWave<Type, TrackingData>::updateCell
         if (!changedCell_[celli])
         {
             changedCell_[celli] = true;
-            changedCells_[nChangedCells_++] = celli;
+            changedCells_.append(celli);
         }
     }
 
@@ -154,11 +155,6 @@ bool Foam::FaceCellWave<Type, TrackingData>::updateCell
 }
 
 
-// Update info for facei, at position pt, with information from
-// neighbouring face/cell.
-// Updates:
-//      - changedFace_, changedFaces_, nChangedFaces_,
-//      - statistics: nEvals_, nUnvisitedFaces_
 template<class Type, class TrackingData>
 bool Foam::FaceCellWave<Type, TrackingData>::updateFace
 (
@@ -169,6 +165,12 @@ bool Foam::FaceCellWave<Type, TrackingData>::updateFace
     Type& faceInfo
 )
 {
+    // Update info for facei, at position pt, with information from
+    // neighbouring face/cell.
+    // Updates:
+    //      - changedFace_, changedFaces_,
+    //      - statistics: nEvals_, nUnvisitedFaces_
+
     nEvals_++;
 
     bool wasValid = faceInfo.valid(td_);
@@ -189,7 +191,7 @@ bool Foam::FaceCellWave<Type, TrackingData>::updateFace
         if (!changedFace_[facei])
         {
             changedFace_[facei] = true;
-            changedFaces_[nChangedFaces_++] = facei;
+            changedFaces_.append(facei);
         }
     }
 
@@ -202,11 +204,6 @@ bool Foam::FaceCellWave<Type, TrackingData>::updateFace
 }
 
 
-// Update info for facei, at position pt, with information from
-// same face.
-// Updates:
-//      - changedFace_, changedFaces_, nChangedFaces_,
-//      - statistics: nEvals_, nUnvisitedFaces_
 template<class Type, class TrackingData>
 bool Foam::FaceCellWave<Type, TrackingData>::updateFace
 (
@@ -216,6 +213,12 @@ bool Foam::FaceCellWave<Type, TrackingData>::updateFace
     Type& faceInfo
 )
 {
+    // Update info for facei, at position pt, with information from
+    // same face.
+    // Updates:
+    //      - changedFace_, changedFaces_,
+    //      - statistics: nEvals_, nUnvisitedFaces_
+
     nEvals_++;
 
     bool wasValid = faceInfo.valid(td_);
@@ -235,7 +238,7 @@ bool Foam::FaceCellWave<Type, TrackingData>::updateFace
         if (!changedFace_[facei])
         {
             changedFace_[facei] = true;
-            changedFaces_[nChangedFaces_++] = facei;
+            changedFaces_.append(facei);
         }
     }
 
@@ -248,13 +251,14 @@ bool Foam::FaceCellWave<Type, TrackingData>::updateFace
 }
 
 
-// For debugging: check status on both sides of cyclic
 template<class Type, class TrackingData>
 void Foam::FaceCellWave<Type, TrackingData>::checkCyclic
 (
     const polyPatch& patch
 ) const
 {
+    // For debugging: check status on both sides of cyclic
+
     const cyclicPolyPatch& nbrPatch =
         refCast<const cyclicPolyPatch>(patch).neighbPatch();
 
@@ -293,7 +297,6 @@ void Foam::FaceCellWave<Type, TrackingData>::checkCyclic
 }
 
 
-// Check if has cyclic patches
 template<class Type, class TrackingData>
 template<class PatchType>
 bool Foam::FaceCellWave<Type, TrackingData>::hasPatch() const
@@ -309,7 +312,6 @@ bool Foam::FaceCellWave<Type, TrackingData>::hasPatch() const
 }
 
 
-// Copy face information into member data
 template<class Type, class TrackingData>
 void Foam::FaceCellWave<Type, TrackingData>::setFaceInfo
 (
@@ -335,12 +337,11 @@ void Foam::FaceCellWave<Type, TrackingData>::setFaceInfo
         // Mark facei as changed, both on list and on face itself.
 
         changedFace_[facei] = true;
-        changedFaces_[nChangedFaces_++] = facei;
+        changedFaces_.append(facei);
     }
 }
 
 
-// Merge face information into member data
 template<class Type, class TrackingData>
 void Foam::FaceCellWave<Type, TrackingData>::mergeFaceInfo
 (
@@ -350,6 +351,8 @@ void Foam::FaceCellWave<Type, TrackingData>::mergeFaceInfo
     const List<Type>& changedFacesInfo
 )
 {
+    // Merge face information into member data
+
     for (label changedFacei = 0; changedFacei < nFaces; changedFacei++)
     {
         const Type& neighbourWallInfo = changedFacesInfo[changedFacei];
@@ -373,9 +376,6 @@ void Foam::FaceCellWave<Type, TrackingData>::mergeFaceInfo
 }
 
 
-// Construct compact patchFace change arrays for a (slice of a) single patch.
-// changedPatchFaces in local patch numbering.
-// Return length of arrays.
 template<class Type, class TrackingData>
 Foam::label Foam::FaceCellWave<Type, TrackingData>::getChangedPatchFaces
 (
@@ -386,6 +386,9 @@ Foam::label Foam::FaceCellWave<Type, TrackingData>::getChangedPatchFaces
     List<Type>& changedPatchFacesInfo
 ) const
 {
+    // Construct compact patchFace change arrays for a (slice of a) single
+    // patch. changedPatchFaces in local patch numbering.
+    // Return length of arrays.
     label nChangedPatchFaces = 0;
 
     for (label i = 0; i < nFaces; i++)
@@ -405,7 +408,6 @@ Foam::label Foam::FaceCellWave<Type, TrackingData>::getChangedPatchFaces
 }
 
 
-// Handle leaving domain. Implementation referred to Type
 template<class Type, class TrackingData>
 void Foam::FaceCellWave<Type, TrackingData>::leaveDomain
 (
@@ -415,6 +417,8 @@ void Foam::FaceCellWave<Type, TrackingData>::leaveDomain
     List<Type>& faceInfo
 ) const
 {
+    // Handle leaving domain. Implementation referred to Type
+
     const vectorField& fc = mesh_.faceCentres();
 
     for (label i = 0; i < nFaces; i++)
@@ -427,7 +431,6 @@ void Foam::FaceCellWave<Type, TrackingData>::leaveDomain
 }
 
 
-// Handle entering domain. Implementation referred to Type
 template<class Type, class TrackingData>
 void Foam::FaceCellWave<Type, TrackingData>::enterDomain
 (
@@ -437,6 +440,8 @@ void Foam::FaceCellWave<Type, TrackingData>::enterDomain
     List<Type>& faceInfo
 ) const
 {
+    // Handle entering domain. Implementation referred to Type
+
     const vectorField& fc = mesh_.faceCentres();
 
     for (label i = 0; i < nFaces; i++)
@@ -449,7 +454,6 @@ void Foam::FaceCellWave<Type, TrackingData>::enterDomain
 }
 
 
-// Transform. Implementation referred to Type
 template<class Type, class TrackingData>
 void Foam::FaceCellWave<Type, TrackingData>::transform
 (
@@ -458,6 +462,8 @@ void Foam::FaceCellWave<Type, TrackingData>::transform
     List<Type>& faceInfo
 )
 {
+    // Transform. Implementation referred to Type
+
     if (rotTensor.size() == 1)
     {
         const tensor& T = rotTensor[0];
@@ -477,7 +483,6 @@ void Foam::FaceCellWave<Type, TrackingData>::transform
 }
 
 
-// Offset mesh face. Used for transferring from one cyclic half to the other.
 template<class Type, class TrackingData>
 void Foam::FaceCellWave<Type, TrackingData>::offset
 (
@@ -487,6 +492,9 @@ void Foam::FaceCellWave<Type, TrackingData>::offset
     labelList& faces
 )
 {
+    // Offset mesh face. Used for transferring from one cyclic half to the
+    // other.
+
     for (label facei = 0; facei < nFaces; facei++)
     {
         faces[facei] += cycOffset;
@@ -494,10 +502,11 @@ void Foam::FaceCellWave<Type, TrackingData>::offset
 }
 
 
-// Tranfer all the information to/from neighbouring processors
 template<class Type, class TrackingData>
 void Foam::FaceCellWave<Type, TrackingData>::handleProcPatches()
 {
+    // Tranfer all the information to/from neighbouring processors
+
     const globalMeshData& pData = mesh_.globalData();
 
     // Which patches are processor patches
@@ -613,10 +622,11 @@ void Foam::FaceCellWave<Type, TrackingData>::handleProcPatches()
 }
 
 
-// Transfer information across cyclic halves.
 template<class Type, class TrackingData>
 void Foam::FaceCellWave<Type, TrackingData>::handleCyclicPatches()
 {
+    // Transfer information across cyclic halves.
+
     forAll(mesh_.boundaryMesh(), patchi)
     {
         const polyPatch& patch = mesh_.boundaryMesh()[patchi];
@@ -698,10 +708,11 @@ void Foam::FaceCellWave<Type, TrackingData>::handleCyclicPatches()
 }
 
 
-// Transfer information across cyclic halves.
 template<class Type, class TrackingData>
 void Foam::FaceCellWave<Type, TrackingData>::handleAMICyclicPatches()
 {
+    // Transfer information across cyclicAMI halves.
+
     forAll(mesh_.boundaryMesh(), patchi)
     {
         const polyPatch& patch = mesh_.boundaryMesh()[patchi];
@@ -802,9 +813,81 @@ void Foam::FaceCellWave<Type, TrackingData>::handleAMICyclicPatches()
 }
 
 
+template<class Type, class TrackingData>
+void Foam::FaceCellWave<Type, TrackingData>::handleExplicitConnections()
+{
+    // Collect changed information
+
+    DynamicList<label> f0Baffle(explicitConnections_.size());
+    DynamicList<Type> f0Info(explicitConnections_.size());
+
+    DynamicList<label> f1Baffle(explicitConnections_.size());
+    DynamicList<Type> f1Info(explicitConnections_.size());
+
+    forAll(explicitConnections_, connI)
+    {
+        const labelPair& baffle = explicitConnections_[connI];
+
+        label f0 = baffle[0];
+        if (changedFace_[f0])
+        {
+            f0Baffle.append(connI);
+            f0Info.append(allFaceInfo_[f0]);
+        }
+
+        label f1 = baffle[1];
+        if (changedFace_[f1])
+        {
+            f1Baffle.append(connI);
+            f1Info.append(allFaceInfo_[f1]);
+        }
+    }
+
+
+    // Update other side with changed information
+
+    forAll(f1Info, index)
+    {
+        const labelPair& baffle = explicitConnections_[f1Baffle[index]];
+
+        label f0 = baffle[0];
+        Type& currentWallInfo = allFaceInfo_[f0];
+
+        if (!currentWallInfo.equal(f1Info[index], td_))
+        {
+            updateFace
+            (
+                f0,
+                f1Info[index],
+                propagationTol_,
+                currentWallInfo
+            );
+        }
+    }
+
+    forAll(f0Info, index)
+    {
+        const labelPair& baffle = explicitConnections_[f0Baffle[index]];
+
+        label f1 = baffle[1];
+        Type& currentWallInfo = allFaceInfo_[f1];
+
+        if (!currentWallInfo.equal(f0Info[index], td_))
+        {
+            updateFace
+            (
+                f1,
+                f0Info[index],
+                propagationTol_,
+                currentWallInfo
+            );
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Set up only. Use setFaceInfo and iterate() to do actual calculation.
 template<class Type, class TrackingData>
 Foam::FaceCellWave<Type, TrackingData>::FaceCellWave
 (
@@ -815,15 +898,14 @@ Foam::FaceCellWave<Type, TrackingData>::FaceCellWave
 )
 :
     mesh_(mesh),
+    explicitConnections_(0),
     allFaceInfo_(allFaceInfo),
     allCellInfo_(allCellInfo),
     td_(td),
     changedFace_(mesh_.nFaces(), false),
     changedFaces_(mesh_.nFaces()),
-    nChangedFaces_(0),
     changedCell_(mesh_.nCells(), false),
     changedCells_(mesh_.nCells()),
-    nChangedCells_(0),
     hasCyclicPatches_(hasPatch<cyclicPolyPatch>()),
     hasCyclicAMIPatches_
     (
@@ -851,8 +933,6 @@ Foam::FaceCellWave<Type, TrackingData>::FaceCellWave
 }
 
 
-// Iterate, propagating changedFacesInfo across mesh, until no change (or
-// maxIter reached). Initial cell values specified.
 template<class Type, class TrackingData>
 Foam::FaceCellWave<Type, TrackingData>::FaceCellWave
 (
@@ -866,15 +946,14 @@ Foam::FaceCellWave<Type, TrackingData>::FaceCellWave
 )
 :
     mesh_(mesh),
+    explicitConnections_(0),
     allFaceInfo_(allFaceInfo),
     allCellInfo_(allCellInfo),
     td_(td),
     changedFace_(mesh_.nFaces(), false),
     changedFaces_(mesh_.nFaces()),
-    nChangedFaces_(0),
     changedCell_(mesh_.nCells(), false),
     changedCells_(mesh_.nCells()),
-    nChangedCells_(0),
     hasCyclicPatches_(hasPatch<cyclicPolyPatch>()),
     hasCyclicAMIPatches_
     (
@@ -911,8 +990,75 @@ Foam::FaceCellWave<Type, TrackingData>::FaceCellWave
         FatalErrorInFunction
             << "Maximum number of iterations reached. Increase maxIter." << endl
             << "    maxIter:" << maxIter << endl
-            << "    nChangedCells:" << nChangedCells_ << endl
+            << "    nChangedCells:" << changedCells_.size() << endl
-            << "    nChangedFaces:" << nChangedFaces_ << endl
+            << "    nChangedFaces:" << changedFaces_.size() << endl
+            << exit(FatalError);
+    }
+}
+
+
+template<class Type, class TrackingData>
+Foam::FaceCellWave<Type, TrackingData>::FaceCellWave
+(
+    const polyMesh& mesh,
+    const labelPairList& explicitConnections,
+    const bool handleCyclicAMI,
+    const labelList& changedFaces,
+    const List<Type>& changedFacesInfo,
+    UList<Type>& allFaceInfo,
+    UList<Type>& allCellInfo,
+    const label maxIter,
+    TrackingData& td
+)
+:
+    mesh_(mesh),
+    explicitConnections_(explicitConnections),
+    allFaceInfo_(allFaceInfo),
+    allCellInfo_(allCellInfo),
+    td_(td),
+    changedFace_(mesh_.nFaces(), false),
+    changedFaces_(mesh_.nFaces()),
+    changedCell_(mesh_.nCells(), false),
+    changedCells_(mesh_.nCells()),
+    hasCyclicPatches_(hasPatch<cyclicPolyPatch>()),
+    hasCyclicAMIPatches_
+    (
+        handleCyclicAMI
+     && returnReduce(hasPatch<cyclicAMIPolyPatch>(), orOp<bool>())
+    ),
+    nEvals_(0),
+    nUnvisitedCells_(mesh_.nCells()),
+    nUnvisitedFaces_(mesh_.nFaces())
+{
+    if
+    (
+        allFaceInfo.size() != mesh_.nFaces()
+     || allCellInfo.size() != mesh_.nCells()
+    )
+    {
+        FatalErrorInFunction
+            << "face and cell storage not the size of mesh faces, cells:"
+            << endl
+            << "    allFaceInfo   :" << allFaceInfo.size() << endl
+            << "    mesh_.nFaces():" << mesh_.nFaces() << endl
+            << "    allCellInfo   :" << allCellInfo.size() << endl
+            << "    mesh_.nCells():" << mesh_.nCells()
+            << exit(FatalError);
+    }
+
+    // Copy initial changed faces data
+    setFaceInfo(changedFaces, changedFacesInfo);
+
+    // Iterate until nothing changes
+    label iter = iterate(maxIter);
+
+    if ((maxIter > 0) && (iter >= maxIter))
+    {
+        FatalErrorInFunction
+            << "Maximum number of iterations reached. Increase maxIter." << endl
+            << "    maxIter:" << maxIter << endl
+            << "    nChangedCells:" << changedCells_.size() << endl
+            << "    nChangedFaces:" << changedFaces_.size() << endl
             << exit(FatalError);
     }
 }
@@ -920,7 +1066,6 @@ Foam::FaceCellWave<Type, TrackingData>::FaceCellWave
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-
 template<class Type, class TrackingData>
 Foam::label Foam::FaceCellWave<Type, TrackingData>::getUnsetCells() const
 {
@@ -935,21 +1080,16 @@ Foam::label Foam::FaceCellWave<Type, TrackingData>::getUnsetFaces() const
 }
 
 
-
-// Propagate cell to face
 template<class Type, class TrackingData>
 Foam::label Foam::FaceCellWave<Type, TrackingData>::faceToCell()
 {
+    // Propagate face to cell
+
     const labelList& owner = mesh_.faceOwner();
     const labelList& neighbour = mesh_.faceNeighbour();
     label nInternalFaces = mesh_.nInternalFaces();
 
-    for
+    forAll(changedFaces_, changedFacei)
-    (
-        label changedFacei = 0;
-        changedFacei < nChangedFaces_;
-        changedFacei++
-    )
     {
         label facei = changedFaces_[changedFacei];
         if (!changedFace_[facei])
@@ -1005,15 +1145,15 @@ Foam::label Foam::FaceCellWave<Type, TrackingData>::faceToCell()
     }
 
     // Handled all changed faces by now
-    nChangedFaces_ = 0;
+    changedFaces_.clear();
 
     if (debug & 2)
     {
-        Pout<< " Changed cells            : " << nChangedCells_ << endl;
+        Pout<< " Changed cells            : " << changedCells_.size() << endl;
     }
 
-    // Sum nChangedCells over all procs
+    // Sum changedCells over all procs
-    label totNChanged = nChangedCells_;
+    label totNChanged = changedCells_.size();
 
     reduce(totNChanged, sumOp<label>());
 
@@ -1021,18 +1161,14 @@ Foam::label Foam::FaceCellWave<Type, TrackingData>::faceToCell()
 }
 
 
-// Propagate cell to face
 template<class Type, class TrackingData>
 Foam::label Foam::FaceCellWave<Type, TrackingData>::cellToFace()
 {
+    // Propagate cell to face
+
     const cellList& cells = mesh_.cells();
 
-    for
+    forAll(changedCells_, changedCelli)
-    (
-        label changedCelli = 0;
-        changedCelli < nChangedCells_;
-        changedCelli++
-    )
     {
         label celli = changedCells_[changedCelli];
         if (!changedCell_[celli])
@@ -1070,7 +1206,11 @@ Foam::label Foam::FaceCellWave<Type, TrackingData>::cellToFace()
     }
 
     // Handled all changed cells by now
-    nChangedCells_ = 0;
+    changedCells_.clear();
+
+
+    // Transfer across any explicitly provided internal connections
+    handleExplicitConnections();
 
     if (hasCyclicPatches_)
     {
@@ -1091,11 +1231,11 @@ Foam::label Foam::FaceCellWave<Type, TrackingData>::cellToFace()
 
     if (debug & 2)
     {
-        Pout<< " Changed faces            : " << nChangedFaces_ << endl;
+        Pout<< " Changed faces            : " << changedFaces_.size() << endl;
     }
 
     // Sum nChangedFaces over all procs
-    label totNChanged = nChangedFaces_;
+    label totNChanged = changedFaces_.size();
 
     reduce(totNChanged, sumOp<label>());
 

src/meshTools/algorithms/MeshWave/FaceCellWave.H

@@ -47,9 +47,10 @@ SourceFiles
 #ifndef FaceCellWave_H
 #define FaceCellWave_H
 
-#include "boolList.H"
+#include "PackedBoolList.H"
-#include "labelList.H"
+#include "DynamicList.H"
 #include "primitiveFieldsFwd.H"
+#include "labelPair.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -76,11 +77,27 @@ class FaceCellWave
 :
     public FaceCellWaveName
 {
+private:
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        FaceCellWave(const FaceCellWave&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const FaceCellWave&);
+
+
+protected:
+
     // Private data
 
         //- Reference to mesh
         const polyMesh& mesh_;
 
+        //- Optional boundary faces that information should travel through
+        const labelPairList explicitConnections_;
+
         //- Information for all faces
         UList<Type>& allFaceInfo_;
 
@@ -91,18 +108,17 @@ class FaceCellWave
         TrackingData& td_;
 
         //- Has face changed
-        boolList changedFace_;
+        PackedBoolList changedFace_;
 
         //- List of changed faces
-        labelList changedFaces_;
+        DynamicList<label> changedFaces_;
 
-        //- Number of changed faces
+        //- Has cell changed
-        label nChangedFaces_;
+        PackedBoolList changedCell_;
 
         // Cells that have changed
-        boolList changedCell_;
+        DynamicList<label> changedCells_;
-        labelList changedCells_;
+
-        label nChangedCells_;
 
         //- Contains cyclics
         const bool hasCyclicPatches_;
@@ -118,15 +134,6 @@ class FaceCellWave
         label nUnvisitedFaces_;
 
 
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        FaceCellWave(const FaceCellWave&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const FaceCellWave&);
-
-
         //- Updates cellInfo with information from neighbour. Updates all
         //  statistics.
         bool updateCell
@@ -232,6 +239,10 @@ class FaceCellWave
             //- Merge data from across AMI cyclics
             void handleAMICyclicPatches();
 
+            //- Merge data across explicitly provided local connections (usually
+            //  baffles)
+            void handleExplicitConnections();
+
 
       // Private static data
 
@@ -286,6 +297,28 @@ public:
             TrackingData& td = dummyTrackData_
         );
 
+        //- Construct from mesh and explicitly connected boundary faces
+        //  and list of changed faces with the Type
+        //  for these faces. Iterates until nothing changes or maxIter reached.
+        //  (maxIter can be 0)
+        FaceCellWave
+        (
+            const polyMesh&,
+            const labelPairList& explicitConnections,
+            const bool handleCyclicAMI,
+            const labelList& initialChangedFaces,
+            const List<Type>& changedFacesInfo,
+            UList<Type>& allFaceInfo,
+            UList<Type>& allCellInfo,
+            const label maxIter,
+            TrackingData& td = dummyTrackData_
+        );
+
+
+    //- Destructor
+    virtual ~FaceCellWave()
+    {}
+
 
     // Member Functions
 
@@ -337,17 +370,16 @@ public:
 
             //- Propagate from face to cell. Returns total number of cells
             //  (over all processors) changed.
-            label faceToCell();
+           virtual label faceToCell();
 
             //- Propagate from cell to face. Returns total number of faces
             //  (over all processors) changed. (Faces on processorpatches are
             //  counted double)
-            label cellToFace();
+            virtual label cellToFace();
 
             //- Iterate until no changes or maxIter reached.  Returns actual
             //  number of iterations.
-            label iterate(const label maxIter);
+            virtual label iterate(const label maxIter);
-
 };
 
 

src/meshTools/regionSplit/regionSplit.C

@@ -24,10 +24,12 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "regionSplit.H"
+#include "FaceCellWave.H"
 #include "cyclicPolyPatch.H"
 #include "processorPolyPatch.H"
 #include "globalIndex.H"
 #include "syncTools.H"
+#include "minData.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -40,335 +42,93 @@ defineTypeNameAndDebug(regionSplit, 0);
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-// Handle (non-processor) coupled faces.
+void Foam::regionSplit::calcNonCompactRegionSplit
-void Foam::regionSplit::transferCoupledFaceRegion
-(
-    const label facei,
-    const label otherFacei,
-
-    labelList& faceRegion,
-    DynamicList<label>& newChangedFaces
-) const
-{
-    if (faceRegion[facei] >= 0)
-    {
-        if (faceRegion[otherFacei] == -1)
-        {
-            faceRegion[otherFacei] = faceRegion[facei];
-            newChangedFaces.append(otherFacei);
-        }
-        else if (faceRegion[otherFacei] == -2)
-        {
-            // otherFacei blocked but facei is not. Is illegal for coupled
-            // faces, not for explicit connections.
-        }
-        else if (faceRegion[otherFacei] != faceRegion[facei])
-        {
-            FatalErrorInFunction
-                  << "Problem : coupled face " << facei
-                  << " on patch " << mesh().boundaryMesh().whichPatch(facei)
-                  << " has region " << faceRegion[facei]
-                  << " but coupled face " << otherFacei
-                  << " has region " << faceRegion[otherFacei]
-                  << endl
-                  << "Is your blocked faces specification"
-                  << " synchronized across coupled boundaries?"
-                  << abort(FatalError);
-        }
-    }
-    else if (faceRegion[facei] == -1)
-    {
-        if (faceRegion[otherFacei] >= 0)
-        {
-            faceRegion[facei] = faceRegion[otherFacei];
-            newChangedFaces.append(facei);
-        }
-        else if (faceRegion[otherFacei] == -2)
-        {
-            // otherFacei blocked but facei is not. Is illegal for coupled
-            // faces, not for explicit connections.
-        }
-    }
-}
-
-
-void Foam::regionSplit::fillSeedMask
-(
-    const List<labelPair>& explicitConnections,
-    labelList& cellRegion,
-    labelList& faceRegion,
-    const label seedCellID,
-    const label markValue
-) const
-{
-    // Do seed cell
-    cellRegion[seedCellID] = markValue;
-
-
-    // Collect faces on seed cell
-    const cell& cFaces = mesh().cells()[seedCellID];
-
-    label nFaces = 0;
-
-    labelList changedFaces(cFaces.size());
-
-    forAll(cFaces, i)
-    {
-        label facei = cFaces[i];
-
-        if (faceRegion[facei] == -1)
-        {
-            faceRegion[facei] = markValue;
-            changedFaces[nFaces++] = facei;
-        }
-    }
-    changedFaces.setSize(nFaces);
-
-
-    // Loop over changed faces. FaceCellWave in small.
-
-    while (changedFaces.size())
-    {
-        //if (debug)
-        //{
-        //    Pout<< "regionSplit::fillSeedMask : changedFaces:"
-        //        << changedFaces.size() << endl;
-        //}
-
-        DynamicList<label> changedCells(changedFaces.size());
-
-        forAll(changedFaces, i)
-        {
-            label facei = changedFaces[i];
-
-            label own = mesh().faceOwner()[facei];
-
-            if (cellRegion[own] == -1)
-            {
-                cellRegion[own] = markValue;
-                changedCells.append(own);
-            }
-
-            if (mesh().isInternalFace(facei))
-            {
-                label nei = mesh().faceNeighbour()[facei];
-
-                if (cellRegion[nei] == -1)
-                {
-                    cellRegion[nei] = markValue;
-                    changedCells.append(nei);
-                }
-            }
-        }
-
-
-        //if (debug)
-        //{
-        //    Pout<< "regionSplit::fillSeedMask : changedCells:"
-        //        << changedCells.size() << endl;
-        //}
-
-        // Loop over changedCells and collect faces
-        DynamicList<label> newChangedFaces(changedCells.size());
-
-        forAll(changedCells, i)
-        {
-            label celli = changedCells[i];
-
-            const cell& cFaces = mesh().cells()[celli];
-
-            forAll(cFaces, cFacei)
-            {
-                label facei = cFaces[cFacei];
-
-                if (faceRegion[facei] == -1)
-                {
-                    faceRegion[facei] = markValue;
-                    newChangedFaces.append(facei);
-                }
-            }
-        }
-
-
-        //if (debug)
-        //{
-        //    Pout<< "regionSplit::fillSeedMask : changedFaces before sync:"
-        //        << changedFaces.size() << endl;
-        //}
-
-
-        // Check for changes to any locally coupled face.
-        // Global connections are done later.
-
-        const polyBoundaryMesh& patches = mesh().boundaryMesh();
-
-        forAll(patches, patchi)
-        {
-            const polyPatch& pp = patches[patchi];
-
-            if
-            (
-                isA<cyclicPolyPatch>(pp)
-             && refCast<const cyclicPolyPatch>(pp).owner()
-            )
-            {
-                // Transfer from neighbourPatch to here or vice versa.
-
-                const cyclicPolyPatch& cycPatch =
-                    refCast<const cyclicPolyPatch>(pp);
-
-                label facei = cycPatch.start();
-
-                forAll(cycPatch, i)
-                {
-                    label otherFacei = cycPatch.transformGlobalFace(facei);
-
-                    transferCoupledFaceRegion
-                    (
-                        facei,
-                        otherFacei,
-                        faceRegion,
-                        newChangedFaces
-                    );
-
-                    facei++;
-                }
-            }
-        }
-        forAll(explicitConnections, i)
-        {
-            transferCoupledFaceRegion
-            (
-                explicitConnections[i][0],
-                explicitConnections[i][1],
-                faceRegion,
-                newChangedFaces
-            );
-        }
-
-        //if (debug)
-        //{
-        //    Pout<< "regionSplit::fillSeedMask : changedFaces after sync:"
-        //        << newChangedFaces.size() << endl;
-        //}
-
-        changedFaces.transfer(newChangedFaces);
-    }
-}
-
-
-Foam::label Foam::regionSplit::calcLocalRegionSplit
 (
+    const globalIndex& globalFaces,
     const boolList& blockedFace,
     const List<labelPair>& explicitConnections,
 
     labelList& cellRegion
 ) const
 {
-    if (debug)
+    // Field on cells and faces.
-    {
+    List<minData> cellData(mesh().nCells());
-        if (blockedFace.size())
+    List<minData> faceData(mesh().nFaces());
-        {
-            // Check that blockedFace is synced.
-            boolList syncBlockedFace(blockedFace);
-            syncTools::swapFaceList(mesh(), syncBlockedFace);
 
-            forAll(syncBlockedFace, facei)
+    // Take over blockedFaces by seeding a negative number
+    // (so is always less than the decomposition)
+    label nUnblocked = 0;
+    forAll(faceData, facei)
     {
-                if (syncBlockedFace[facei] != blockedFace[facei])
+        if (blockedFace.size() && blockedFace[facei])
         {
-                    FatalErrorInFunction
+            faceData[facei] = minData(-2);
-                        << "Face " << facei << " not synchronised. My value:"
-                        << blockedFace[facei] << "  coupled value:"
-                        << syncBlockedFace[facei]
-                        << abort(FatalError);
-                }
         }
+        else
+        {
+            nUnblocked++;
         }
     }
 
-    // Region per face.
+    // Seed unblocked faces
-    // -1 unassigned
+    labelList seedFaces(nUnblocked);
-    // -2 blocked
+    List<minData> seedData(nUnblocked);
-    labelList faceRegion(mesh().nFaces(), -1);
+    nUnblocked = 0;
 
-    if (blockedFace.size())
+
+    forAll(faceData, facei)
     {
-        forAll(blockedFace, facei)
+        if (blockedFace.empty() || !blockedFace[facei])
         {
-            if (blockedFace[facei])
+            seedFaces[nUnblocked] = facei;
-            {
+            // Seed face with globally unique number
-                faceRegion[facei] = -2;
+            seedData[nUnblocked] = minData(globalFaces.toGlobal(facei));
-            }
+            nUnblocked++;
-        }
-    }
-
-
-    // Assign local regions
-    // ~~~~~~~~~~~~~~~~~~~~
-
-    // Start with region 0
-    label nLocalRegions = 0;
-
-    label unsetCelli = 0;
-
-    do
-    {
-        // Find first unset cell
-
-        for (; unsetCelli < mesh().nCells(); unsetCelli++)
-        {
-            if (cellRegion[unsetCelli] == -1)
-            {
-                break;
         }
     }
 
-        if (unsetCelli >= mesh().nCells())
-        {
-            break;
-        }
 
-        fillSeedMask
+    // Propagate information inwards
+    FaceCellWave<minData> deltaCalc
     (
+        mesh(),
         explicitConnections,
-            cellRegion,
+        false,  // disable walking through cyclicAMI for backwards compatibility
-            faceRegion,
+        seedFaces,
-            unsetCelli,
+        seedData,
-            nLocalRegions
+        faceData,
+        cellData,
+        mesh().globalData().nTotalCells()+1
     );
 
-        // Current unsetCell has now been handled. Go to next region.
-        nLocalRegions++;
-        unsetCelli++;
-    }
-    while (true);
 
-
+    // And extract
-    if (debug)
+    cellRegion.setSize(mesh().nCells());
-    {
     forAll(cellRegion, celli)
     {
-            if (cellRegion[celli] < 0)
+        if (cellData[celli].valid(deltaCalc.data()))
         {
-                FatalErrorInFunction
+            cellRegion[celli] = cellData[celli].data();
-                    << "cell:" << celli << " region:" << cellRegion[celli]
-                    << abort(FatalError);
-            }
         }
+        else
+        {
+            // Unvisited cell -> only possible if surrounded by blocked faces.
+            // If so make up region from any of the faces
+            const cell& cFaces = mesh().cells()[celli];
+            label facei = cFaces[0];
 
-        forAll(faceRegion, facei)
+            if (blockedFace.size() && !blockedFace[facei])
             {
-            if (faceRegion[facei] == -1)
+                FatalErrorIn("regionSplit::calcNonCompactRegionSplit(..)")
-            {
+                    << "Problem: unblocked face " << facei
-                FatalErrorInFunction
+                    << " at " << mesh().faceCentres()[facei]
-                    << "face:" << facei << " region:" << faceRegion[facei]
+                    << " on unassigned cell " << celli
-                    << abort(FatalError);
+                    << mesh().cellCentres()[celli]
+                    << exit(FatalError);
+            }
+            cellRegion[celli] = globalFaces.toGlobal(facei);
         }
     }
-    }
-
-    return nLocalRegions;
 }
 
 
@@ -383,176 +143,142 @@ Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
 {
     // See header in regionSplit.H
 
-    // 1. Do local analysis
+
-    label nLocalRegions = calcLocalRegionSplit
+    if (!doGlobalRegions)
+    {
+        // Block all parallel faces to avoid comms across
+        boolList coupledOrBlockedFace(blockedFace);
+        const polyBoundaryMesh& pbm = mesh().boundaryMesh();
+
+        if (coupledOrBlockedFace.size())
+        {
+            forAll(pbm, patchi)
+            {
+                const polyPatch& pp = pbm[patchi];
+                if (isA<processorPolyPatch>(pp))
+                {
+                    label facei = pp.start();
+                    forAll(pp, i)
+                    {
+                        coupledOrBlockedFace[facei++] = true;
+                    }
+                }
+            }
+        }
+
+        // Create dummy (local only) globalIndex
+        labelList offsets(Pstream::nProcs()+1, 0);
+        for (label i = Pstream::myProcNo()+1; i < offsets.size(); i++)
+        {
+            offsets[i] = mesh().nFaces();
+        }
+        const globalIndex globalRegions(offsets.xfer());
+
+        // Minimise regions across connected cells
+        // Note: still uses global decisions so all processors are running
+        //       in lock-step, i.e. slowest determines overall time.
+        //       To avoid this we could switch off Pstream::parRun.
+        calcNonCompactRegionSplit
         (
+            globalRegions,
+            coupledOrBlockedFace,
+            explicitConnections,
+            cellRegion
+        );
+
+        // Compact
+        Map<label> globalToCompact(mesh().nCells()/8);
+        forAll(cellRegion, celli)
+        {
+            label region = cellRegion[celli];
+
+            label globalRegion;
+
+            Map<label>::const_iterator fnd = globalToCompact.find(region);
+            if (fnd == globalToCompact.end())
+            {
+                globalRegion = globalRegions.toGlobal(globalToCompact.size());
+                globalToCompact.insert(region, globalRegion);
+            }
+            else
+            {
+                globalRegion = fnd();
+            }
+            cellRegion[celli] = globalRegion;
+        }
+
+
+        // Return globalIndex with size = localSize and all regions local
+        labelList compactOffsets(Pstream::nProcs()+1, 0);
+        for (label i = Pstream::myProcNo()+1; i < compactOffsets.size(); i++)
+        {
+            compactOffsets[i] = globalToCompact.size();
+        }
+
+        return autoPtr<globalIndex>(new globalIndex(compactOffsets.xfer()));
+    }
+
+
+
+    // Initial global region numbers
+    const globalIndex globalRegions(mesh().nFaces());
+
+    // Minimise regions across connected cells (including parallel)
+    calcNonCompactRegionSplit
+    (
+        globalRegions,
         blockedFace,
         explicitConnections,
         cellRegion
     );
 
-    if (!doGlobalRegions)
+
+    // Now our cellRegion will have
+    // - non-local regions (i.e. originating from other processors)
+    // - non-compact locally originating regions
+    // so we'll need to compact
+
+    // 4a: count per originating processor the number of regions
+    labelList nOriginating(Pstream::nProcs(), 0);
     {
-        return autoPtr<globalIndex>(new globalIndex(nLocalRegions));
+        labelHashSet haveRegion(mesh().nCells()/8);
-    }
 
-
-    // 2. Assign global regions
-    // ~~~~~~~~~~~~~~~~~~~~~~~~
-    // Offset local regions to create unique global regions.
-
-    globalIndex globalRegions(nLocalRegions);
-
-
-    // Convert regions to global ones
         forAll(cellRegion, celli)
         {
-        cellRegion[celli] = globalRegions.toGlobal(cellRegion[celli]);
+            label region = cellRegion[celli];
-    }
 
-
+            // Count originating processor. Use isLocal as efficiency since
-    // 3. Merge global regions
+            // most cells are locally originating.
-    // ~~~~~~~~~~~~~~~~~~~~~~~
+            if (globalRegions.isLocal(region))
-    // Regions across non-blocked proc patches get merged.
-    // This will set merged global regions to be the min of both.
-    // (this will create gaps in the global region list so they will get
-    // merged later on)
-
-    while (true)
             {
-        if (debug)
+                if (haveRegion.insert(region))
                 {
-            Pout<< nl << "-- Starting Iteration --" << endl;
+                    nOriginating[Pstream::myProcNo()]++;
-        }
-
-
-        const polyBoundaryMesh& patches = mesh().boundaryMesh();
-
-        labelList nbrRegion(mesh().nFaces()-mesh().nInternalFaces(), -1);
-        forAll(patches, patchi)
-        {
-            const polyPatch& pp = patches[patchi];
-
-            if (pp.coupled())
-            {
-                const labelList& patchCells = pp.faceCells();
-                SubList<label> patchNbrRegion
-                (
-                    nbrRegion,
-                    pp.size(),
-                    pp.start()-mesh().nInternalFaces()
-                );
-
-                forAll(patchCells, i)
-                {
-                    label facei = pp.start()+i;
-                    if (!blockedFace.size() || !blockedFace[facei])
-                    {
-                        patchNbrRegion[i] = cellRegion[patchCells[i]];
                 }
             }
-            }
+            else
-        }
-        syncTools::swapBoundaryFaceList(mesh(), nbrRegion);
-
-        Map<label> globalToMerged(mesh().nFaces()-mesh().nInternalFaces());
-
-        forAll(patches, patchi)
             {
-            const polyPatch& pp = patches[patchi];
+                label proci = globalRegions.whichProcID(region);
-
+                if (haveRegion.insert(region))
-            if (pp.coupled())
                 {
-                const labelList& patchCells = pp.faceCells();
+                    nOriginating[proci]++;
-                SubList<label> patchNbrRegion
-                (
-                    nbrRegion,
-                    pp.size(),
-                    pp.start()-mesh().nInternalFaces()
-                );
-
-                forAll(patchCells, i)
-                {
-                    label facei = pp.start()+i;
-
-                    if (!blockedFace.size() || !blockedFace[facei])
-                    {
-                        if (patchNbrRegion[i] < cellRegion[patchCells[i]])
-                        {
-                            //Pout<< "on patch:" << pp.name()
-                            //    << " cell:" << patchCells[i]
-                            //    << " at:"
-                            // << mesh().cellCentres()[patchCells[i]]
-                            //    << " was:" << cellRegion[patchCells[i]]
-                            //    << " nbr:" << patchNbrRegion[i]
-                            //    << endl;
-
-                            globalToMerged.insert
-                            (
-                                cellRegion[patchCells[i]],
-                                patchNbrRegion[i]
-                            );
                 }
             }
         }
     }
-        }
-
-
-        label nMerged = returnReduce(globalToMerged.size(), sumOp<label>());
-
-        if (debug)
-        {
-            Pout<< "nMerged:" << nMerged << endl;
-        }
-
-        if (nMerged == 0)
-        {
-            break;
-        }
-
-        // Renumber the regions according to the globalToMerged
-        forAll(cellRegion, celli)
-        {
-            label regionI = cellRegion[celli];
-            Map<label>::const_iterator iter = globalToMerged.find(regionI);
-            if (iter != globalToMerged.end())
-            {
-                 cellRegion[celli] = iter();
-            }
-        }
-    }
-
-
-    // Now our cellRegion will have non-local elements in it. So compact
-    // it.
-
-    // 4a: count. Use a labelHashSet to count regions only once.
-    label nCompact = 0;
-    {
-        labelHashSet localRegion(mesh().nFaces()-mesh().nInternalFaces());
-        forAll(cellRegion, celli)
-        {
-            if
-            (
-                globalRegions.isLocal(cellRegion[celli])
-             && localRegion.insert(cellRegion[celli])
-            )
-            {
-                nCompact++;
-            }
-        }
-    }
 
     if (debug)
     {
-        Pout<< "Compacted from " << nLocalRegions
+        Pout<< "Counted " << nOriginating[Pstream::myProcNo()]
-            << " down to " << nCompact << " local regions." << endl;
+            << " local regions." << endl;
     }
 
 
     // Global numbering for compacted local regions
-    autoPtr<globalIndex> globalCompactPtr(new globalIndex(nCompact));
+    autoPtr<globalIndex> globalCompactPtr
+    (
+        new globalIndex(nOriginating[Pstream::myProcNo()])
+    );
     const globalIndex& globalCompact = globalCompactPtr();
 
 
@@ -563,12 +289,15 @@ Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
     // labelList.
 
     // Local compaction map
-    Map<label> globalToCompact(2*nCompact);
+    Map<label> globalToCompact(2*nOriginating[Pstream::myProcNo()]);
     // Remote regions we want the compact number for
     List<labelHashSet> nonLocal(Pstream::nProcs());
     forAll(nonLocal, proci)
     {
-        nonLocal[proci].resize((nLocalRegions-nCompact)/Pstream::nProcs());
+        if (proci != Pstream::myProcNo())
+        {
+            nonLocal[proci].resize(2*nOriginating[proci]);
+        }
     }
 
     forAll(cellRegion, celli)
@@ -576,15 +305,12 @@ Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
         label region = cellRegion[celli];
         if (globalRegions.isLocal(region))
         {
-            Map<label>::const_iterator iter = globalToCompact.find(region);
+            // Insert new compact region (if not yet present)
-            if (iter == globalToCompact.end())
+            globalToCompact.insert
-            {
-                label compactRegion = globalCompact.toGlobal
             (
-                    globalToCompact.size()
+                region,
+                globalCompact.toGlobal(globalToCompact.size())
             );
-                globalToCompact.insert(region, compactRegion);
-            }
         }
         else
         {
@@ -598,22 +324,17 @@ Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
     // Convert the nonLocal (labelHashSets) to labelLists.
 
     labelListList sendNonLocal(Pstream::nProcs());
-    labelList nNonLocal(Pstream::nProcs(), 0);
     forAll(sendNonLocal, proci)
     {
-        sendNonLocal[proci].setSize(nonLocal[proci].size());
+        sendNonLocal[proci] = nonLocal[proci].toc();
-        forAllConstIter(labelHashSet, nonLocal[proci], iter)
-        {
-            sendNonLocal[proci][nNonLocal[proci]++] = iter.key();
-        }
     }
 
     if (debug)
     {
-        forAll(nNonLocal, proci)
+        forAll(sendNonLocal, proci)
         {
             Pout<< "    from processor " << proci
-                << " want " << nNonLocal[proci]
+                << " want " << sendNonLocal[proci].size()
                 << " region numbers."
                 << endl;
         }

src/meshTools/regionSplit/regionSplit.H

@@ -89,6 +89,9 @@ Description
 
     Can optionally keep all regions local to the processor.
 
+    Note: does not walk across cyclicAMI/cyclicACMI - since not 'coupled()'
+          at the patch level.
+
 
 SourceFiles
     regionSplit.C
@@ -126,30 +129,10 @@ class regionSplit
 
     // Private Member Functions
 
-        //- Transfer faceRegion data from one face to the other (or vice versa)
+        //- Calculate region split in non-compact (global) numbering.
-        void transferCoupledFaceRegion
+        void calcNonCompactRegionSplit
-        (
-            const label facei,
-            const label otherFacei,
-
-            labelList& faceRegion,
-            DynamicList<label>& newChangedFaces
-        ) const;
-
-        //- Given a seed cell label, fill cellRegion/faceRegion with markValue
-        //  for contiguous region around it
-        void fillSeedMask
-        (
-            const List<labelPair>& explicitConnections,
-            labelList& cellRegion,
-            labelList& faceRegion,
-            const label seedCellID,
-            const label markValue
-        ) const;
-
-        //- Calculate local region split. Return number of regions.
-        label calcLocalRegionSplit
         (
+            const globalIndex& globalFaces,
             const boolList& blockedFace,
             const List<labelPair>& explicitConnections,
             labelList& cellRegion