Merge branch 'master' of ssh://dm/home/dm4/OpenFOAM/OpenFOAM-dev

applications/solvers/lagrangian/sprayFoam/Make/options

@@ -1,4 +1,5 @@
 EXE_INC = \
+    -I$(FOAM_SOLVERS)/lagrangian/reactingParcelFoam \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I${LIB_SRC}/meshTools/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
@@ -20,8 +21,7 @@ EXE_INC = \
     -I$(LIB_SRC)/ODE/lnInclude \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
     -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
-    -I$(LIB_SRC)/combustionModels/lnInclude \
+    -I$(LIB_SRC)/combustionModels/lnInclude
-    -I$(FOAM_SOLVERS)/lagrangian/reactingParcelFoam
 
 
 EXE_LIBS = \

src/meshTools/regionSplit/regionSplit.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -128,7 +128,7 @@ void Foam::regionSplit::fillSeedMask
     changedFaces.setSize(nFaces);
 
 
-    // Loop over changed faces. MeshWave in small.
+    // Loop over changed faces. FaceCellWave in small.
 
     while (changedFaces.size())
     {
@@ -260,7 +260,7 @@ void Foam::regionSplit::fillSeedMask
 }
 
 
-Foam::label Foam::regionSplit::calcRegionSplit
+Foam::label Foam::regionSplit::calcLocalRegionSplit
 (
     const boolList& blockedFace,
     const List<labelPair>& explicitConnections,
@@ -282,7 +282,7 @@ Foam::label Foam::regionSplit::calcRegionSplit
                 {
                     FatalErrorIn
                     (
-                        "regionSplit::calcRegionSplit(..)"
+                        "regionSplit::calcLocalRegionSplit(..)"
                     )   << "Face " << faceI << " not synchronised. My value:"
                         << blockedFace[faceI] << "  coupled value:"
                         << syncBlockedFace[faceI]
@@ -313,7 +313,7 @@ Foam::label Foam::regionSplit::calcRegionSplit
     // ~~~~~~~~~~~~~~~~~~~~
 
     // Start with region 0
-    label nRegions = 0;
+    label nLocalRegions = 0;
 
     label unsetCellI = 0;
 
@@ -340,11 +340,11 @@ Foam::label Foam::regionSplit::calcRegionSplit
             cellRegion,
             faceRegion,
             unsetCellI,
-            nRegions
+            nLocalRegions
         );
 
         // Current unsetCell has now been handled. Go to next region.
-        nRegions++;
+        nLocalRegions++;
         unsetCellI++;
     }
     while (true);
@@ -356,7 +356,7 @@ Foam::label Foam::regionSplit::calcRegionSplit
         {
             if (cellRegion[cellI] < 0)
             {
-                FatalErrorIn("regionSplit::calcRegionSplit(..)")
+                FatalErrorIn("regionSplit::calcLocalRegionSplit(..)")
                     << "cell:" << cellI << " region:" << cellRegion[cellI]
                     << abort(FatalError);
             }
@@ -366,34 +366,62 @@ Foam::label Foam::regionSplit::calcRegionSplit
         {
             if (faceRegion[faceI] == -1)
             {
-                FatalErrorIn("regionSplit::calcRegionSplit(..)")
+                FatalErrorIn("regionSplit::calcLocalRegionSplit(..)")
                     << "face:" << faceI << " region:" << faceRegion[faceI]
                     << abort(FatalError);
             }
         }
     }
 
+    return nLocalRegions;
+}
 
 
-    // Assign global regions
+Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
-    // ~~~~~~~~~~~~~~~~~~~~~
+(
+    const boolList& blockedFace,
+    const List<labelPair>& explicitConnections,
+
+    labelList& cellRegion
+) const
+{
+    // See header in regionSplit.H
+
+    // 1. Do local analysis
+    label nLocalRegions = calcLocalRegionSplit
+    (
+        blockedFace,
+        explicitConnections,
+        cellRegion
+    );
+
+    if (!Pstream::parRun())
+    {
+        return autoPtr<globalIndex>(new globalIndex(nLocalRegions));
+    }
+
+
+    // 2. Assign global regions
+    // ~~~~~~~~~~~~~~~~~~~~~~~~
     // Offset local regions to create unique global regions.
 
-    globalIndex globalRegions(nRegions);
+    globalIndex globalRegions(nLocalRegions);
 
 
-    // Merge global regions
+    // Convert regions to global ones
-    // ~~~~~~~~~~~~~~~~~~~~
+    forAll(cellRegion, cellI)
+    {
+        cellRegion[cellI] = globalRegions.toGlobal(cellRegion[cellI]);
+    }
+
+
+    // 3. Merge global regions
+    // ~~~~~~~~~~~~~~~~~~~~~~~
     // 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)
 
-    // Map from global to merged global
-    labelList mergedGlobal(identity(globalRegions.size()));
-
-
-    // See if any regions get merged. Only nessecary for parallel
     while (Pstream::parRun())
     {
         if (debug)
@@ -401,123 +429,81 @@ Foam::label Foam::regionSplit::calcRegionSplit
             Pout<< nl << "-- Starting Iteration --" << endl;
         }
 
+
         const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
-        // Send global regions across (or -2 if blocked face)
+        labelList nbrRegion(mesh_.nFaces()-mesh_.nInternalFaces(), -1);
         forAll(patches, patchI)
         {
             const polyPatch& pp = patches[patchI];
 
-            if (isA<processorPolyPatch>(pp))
+            if (pp.coupled())
             {
-                labelList myGlobalRegions(pp.size());
+                const labelList& patchCells = pp.faceCells();
-
+                SubList<label> patchNbrRegion
-                label faceI = pp.start();
-
-                forAll(pp, i)
-                {
-                    if (faceRegion[faceI] < 0)
-                    {
-                        myGlobalRegions[i] = faceRegion[faceI];
-                    }
-                    else
-                    {
-                        myGlobalRegions[i] = mergedGlobal
-                        [globalRegions.toGlobal(faceRegion[faceI])];
-                    }
-
-                    faceI++;
-                }
-
-                OPstream toProcNbr
                 (
-                    Pstream::blocking,
+                    nbrRegion,
-                    refCast<const processorPolyPatch>(pp).neighbProcNo()
+                    pp.size(),
+                    pp.start()-mesh_.nInternalFaces()
                 );
 
-                toProcNbr << myGlobalRegions;
+                forAll(patchCells, i)
+                {
+                    label faceI = pp.start()+i;
+                    if (!blockedFace.size() || !blockedFace[faceI])
+                    {
+                        patchNbrRegion[i] = cellRegion[patchCells[i]];
                     }
                 }
+            }
+        }
+        syncTools::swapBoundaryFaceList(mesh_, nbrRegion);
 
-
+        Map<label> globalToMerged(mesh_.nFaces()-mesh_.nInternalFaces());
-        // Receive global regions
-
-        label nMerged = 0;
 
         forAll(patches, patchI)
         {
             const polyPatch& pp = patches[patchI];
 
-            if (isA<processorPolyPatch>(pp))
+            if (pp.coupled())
             {
-                const processorPolyPatch& procPp =
+                const labelList& patchCells = pp.faceCells();
-                    refCast<const processorPolyPatch>(pp);
+                SubList<label> patchNbrRegion
-
-                IPstream fromProcNbr(Pstream::blocking, procPp.neighbProcNo());
-
-                labelList nbrRegions(fromProcNbr);
-
-
-                // Compare with my regions to see which get merged.
-
-                label faceI = pp.start();
-
-                forAll(pp, i)
-                {
-                    if
                 (
-                        faceRegion[faceI] < 0
+                    nbrRegion,
-                     || nbrRegions[i] < 0
+                    pp.size(),
-                    )
+                    pp.start()-mesh_.nInternalFaces()
+                );
+
+                forAll(patchCells, i)
                 {
-                        if (faceRegion[faceI] != nbrRegions[i])
+                    label faceI = pp.start()+i;
+
+                    if (!blockedFace.size() || !blockedFace[faceI])
                     {
-                            FatalErrorIn("regionSplit::calcRegionSplit(..)")
+                        if (patchNbrRegion[i] < cellRegion[patchCells[i]])
-                                << "On patch:" << pp.name()
-                                << " face:" << faceI
-                                << " my local region:" << faceRegion[faceI]
-                                << " neighbouring region:"
-                                << nbrRegions[i] << nl
-                                << "Maybe your blockedFaces are not"
-                                << " synchronized across coupled faces?"
-                                << abort(FatalError);
-                        }
-                    }
-                    else
                         {
-                        label uncompactGlobal =
+                            //Pout<< "on patch:" << pp.name()
-                            globalRegions.toGlobal(faceRegion[faceI]);
+                            //    << " cell:" << patchCells[i]
+                            //    << " at:"
+                            // << mesh_.cellCentres()[patchCells[i]]
+                            //    << " was:" << cellRegion[patchCells[i]]
+                            //    << " nbr:" << patchNbrRegion[i]
+                            //    << endl;
 
-                        label myGlobal = mergedGlobal[uncompactGlobal];
+                            globalToMerged.insert
-
+                            (
-                        if (myGlobal != nbrRegions[i])
+                                cellRegion[patchCells[i]],
-                        {
+                                patchNbrRegion[i]
-                            label minRegion = min(myGlobal, nbrRegions[i]);
+                            );
-
+                        }
-                            if (debug)
+                    }
-                            {
-                                Pout<< "Merging region " << myGlobal
-                                    << " (on proc " << Pstream::myProcNo()
-                                    << ") and region " << nbrRegions[i]
-                                    << " (on proc " << procPp.neighbProcNo()
-                                    << ") into region " << minRegion << endl;
-                            }
-
-                            mergedGlobal[uncompactGlobal] = minRegion;
-                            mergedGlobal[myGlobal] = minRegion;
-                            mergedGlobal[nbrRegions[i]] = minRegion;
-
-                            nMerged++;
-                        }
-                    }
-
-                    faceI++;
                 }
             }
         }
 
 
-        reduce(nMerged, sumOp<label>());
+        label nMerged = returnReduce(globalToMerged.size(), sumOp<label>());
 
         if (debug)
         {
@@ -529,51 +515,176 @@ Foam::label Foam::regionSplit::calcRegionSplit
             break;
         }
 
-        // Merge the compacted regions.
+        // Renumber the regions according to the globalToMerged
-        Pstream::listCombineGather(mergedGlobal, minEqOp<label>());
+        forAll(cellRegion, cellI)
-        Pstream::listCombineScatter(mergedGlobal);
+        {
+            label regionI = cellRegion[cellI];
+            Map<label>::const_iterator iter = globalToMerged.find(regionI);
+            if (iter != globalToMerged.end())
+            {
+                 cellRegion[cellI] = iter();
+            }
+        }
     }
 
 
-    // Compact global regions
+    // Now our cellRegion will have non-local elements in it. So compact
-    // ~~~~~~~~~~~~~~~~~~~~~~
+    // it.
 
-    // All procs will have the same global mergedGlobal region.
+    // 4a: count. Use a labelHashSet to count regions only once.
-    // There might be gaps in it however so compact.
+    label nCompact = 0;
-
-    labelList mergedToCompacted(globalRegions.size(), -1);
-
-    label compactI = 0;
-
-    forAll(mergedGlobal, i)
     {
-        label merged = mergedGlobal[i];
+        labelHashSet localRegion(mesh_.nFaces()-mesh_.nInternalFaces());
-
+        forAll(cellRegion, cellI)
-        if (mergedToCompacted[merged] == -1)
         {
-            mergedToCompacted[merged] = compactI++;
+            if
+            (
+                globalRegions.isLocal(cellRegion[cellI])
+             && localRegion.insert(cellRegion[cellI])
+            )
+            {
+                nCompact++;
+            }
         }
     }
 
     if (debug)
     {
-        Pout<< "Compacted down to " << compactI << " regions." << endl;
+        Pout<< "Compacted from " << nLocalRegions
+            << " down to " << nCompact << " local regions." << endl;
+    }
+
+
+    // Global numbering for compacted local regions
+    autoPtr<globalIndex> globalCompactPtr(new globalIndex(nCompact));
+    const globalIndex& globalCompact = globalCompactPtr();
+
+
+    // 4b: renumber
+    // Renumber into compact indices. Note that since we've already made
+    // all regions global we now need a Map to store the compacting information
+    // instead of a labelList - otherwise we could have used a straight
+    // labelList.
+
+    // Local compaction map
+    Map<label> globalToCompact(2*nCompact);
+    // Remote regions we want the compact number for
+    List<labelHashSet> nonLocal(Pstream::nProcs());
+    forAll(nonLocal, procI)
+    {
+        nonLocal[procI].resize((nLocalRegions-nCompact)/Pstream::nProcs());
     }
 
-    // Renumber cellRegion to be global regions
     forAll(cellRegion, cellI)
     {
         label region = cellRegion[cellI];
-
+        if (globalRegions.isLocal(region))
-        if (region >= 0)
         {
-            label merged = mergedGlobal[globalRegions.toGlobal(region)];
+            Map<label>::const_iterator iter = globalToCompact.find(region);
-
+            if (iter == globalToCompact.end())
-            cellRegion[cellI] = mergedToCompacted[merged];
+            {
+                label compactRegion = globalCompact.toGlobal
+                (
+                    globalToCompact.size()
+                );
+                globalToCompact.insert(region, compactRegion);
+            }
+        }
+        else
+        {
+            nonLocal[globalRegions.whichProcID(region)].insert(region);
         }
     }
 
-    return compactI;
+
+    // Now we have all the local regions compacted. Now we need to get the
+    // non-local ones from the processors to whom they are local.
+    // Convert the nonLocal (labelHashSets) to labelLists.
+
+    labelListList sendNonLocal(Pstream::nProcs());
+    labelList nNonLocal(Pstream::nProcs(), 0);
+    forAll(sendNonLocal, procI)
+    {
+        sendNonLocal[procI].setSize(nonLocal[procI].size());
+        forAllConstIter(labelHashSet, nonLocal[procI], iter)
+        {
+            sendNonLocal[procI][nNonLocal[procI]++] = iter.key();
+        }
+    }
+
+    if (debug)
+    {
+        forAll(nNonLocal, procI)
+        {
+            Pout<< "    from processor " << procI
+                << " want " << nNonLocal[procI]
+                << " region numbers."
+                << endl;
+        }
+        Pout<< endl;
+    }
+
+
+    // Get the wanted region labels into recvNonLocal
+    labelListList recvNonLocal;
+    labelListList sizes;
+    Pstream::exchange<labelList, label>
+    (
+        sendNonLocal,
+        recvNonLocal,
+        sizes
+    );
+
+    // Now we have the wanted compact region labels that procI wants in
+    // recvNonLocal[procI]. Construct corresponding list of compact
+    // region labels to send back.
+
+    labelListList sendWantedLocal(Pstream::nProcs());
+    forAll(recvNonLocal, procI)
+    {
+        const labelList& nonLocal = recvNonLocal[procI];
+        sendWantedLocal[procI].setSize(nonLocal.size());
+
+        forAll(nonLocal, i)
+        {
+            sendWantedLocal[procI][i] = globalToCompact[nonLocal[i]];
+        }
+    }
+
+
+    // Send back (into recvNonLocal)
+    recvNonLocal.clear();
+    sizes.clear();
+    Pstream::exchange<labelList, label>
+    (
+        sendWantedLocal,
+        recvNonLocal,
+        sizes
+    );
+    sendWantedLocal.clear();
+
+    // Now recvNonLocal contains for every element in setNonLocal the
+    // corresponding compact number. Insert these into the local compaction
+    // map.
+
+    forAll(recvNonLocal, procI)
+    {
+        const labelList& wantedRegions = sendNonLocal[procI];
+        const labelList& compactRegions = recvNonLocal[procI];
+
+        forAll(wantedRegions, i)
+        {
+            globalToCompact.insert(wantedRegions[i], compactRegions[i]);
+        }
+    }
+
+    // Finally renumber the regions
+    forAll(cellRegion, cellI)
+    {
+        cellRegion[cellI] = globalToCompact[cellRegion[cellI]];
+    }
+
+    return globalCompactPtr;
 }
 
 
@@ -582,9 +693,15 @@ Foam::label Foam::regionSplit::calcRegionSplit
 Foam::regionSplit::regionSplit(const polyMesh& mesh)
 :
     labelList(mesh.nCells(), -1),
-    mesh_(mesh),
+    mesh_(mesh)
-    nRegions_(calcRegionSplit(boolList(0, false), List<labelPair>(0), *this))
+{
-{}
+    globalNumberingPtr_ = calcRegionSplit
+    (
+        boolList(0, false), //blockedFaces
+        List<labelPair>(0), //explicitConnections,
+        *this
+    );
+}
 
 
 Foam::regionSplit::regionSplit
@@ -594,9 +711,15 @@ Foam::regionSplit::regionSplit
 )
 :
     labelList(mesh.nCells(), -1),
-    mesh_(mesh),
+    mesh_(mesh)
-    nRegions_(calcRegionSplit(blockedFace, List<labelPair>(0), *this))
+{
-{}
+    globalNumberingPtr_ = calcRegionSplit
+    (
+        blockedFace,        //blockedFaces
+        List<labelPair>(0), //explicitConnections,
+        *this
+    );
+}
 
 
 Foam::regionSplit::regionSplit
@@ -607,12 +730,20 @@ Foam::regionSplit::regionSplit
 )
 :
     labelList(mesh.nCells(), -1),
-    mesh_(mesh),
+    mesh_(mesh)
-    nRegions_(calcRegionSplit(blockedFace, explicitConnections, *this))
+{
-{}
+    globalNumberingPtr_ = calcRegionSplit
+    (
+        blockedFace,            //blockedFaces
+        explicitConnections,    //explicitConnections,
+        *this
+    );
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 
+
+
 // ************************************************************************* //

src/meshTools/regionSplit/regionSplit.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,7 +26,66 @@ Class
 
 Description
     This class separates the mesh into distinct unconnected regions,
-    each of which is then given a label.
+    each of which is then given a label according to globalNumbering().
+
+
+    Say 6 cells, 3 processors, with single baffle on proc1.
+
+              baffle
+                |
+    +---+---+---+---+---+---+
+    |   |   |   |   |   |   |
+    +---+---+---+---+---+---+
+      proc0 | proc1 | proc2
+
+
+
+    1: determine local regions (uncoupled)
+
+    +---+---+---+---+---+---+
+    | 0 | 0 | 0 | 1 | 0 | 0 |
+    +---+---+---+---+---+---+
+      proc0 | proc1 | proc2
+
+
+
+    2: make global
+
+    +---+---+---+---+---+---+
+    | 0 | 0 | 1 | 2 | 3 | 3 |
+    +---+---+---+---+---+---+
+      proc0 | proc1 | proc2
+
+
+
+    3: merge connected across procs
+
+    +---+---+---+---+---+---+
+    | 0 | 0 | 0 | 2 | 2 | 2 |
+    +---+---+---+---+---+---+
+      proc0 | proc1 | proc2
+
+
+
+    4. determine locally owner regions. determine compact numbering for the
+    local regions and send these to all processors that need them:
+
+    proc0 uses regions:
+        - 0 which is local to it.
+    proc1 uses regions
+        - 0 which originates from proc0
+        - 2 which is local to it
+    proc2 uses regions
+        - 2 which originates from proc1
+
+    So proc1 needs to get the compact number for region 0 from proc0 and proc2
+    needs to get the compact number for region 2 from proc1:
+
+    +---+---+---+---+---+---+
+    | 0 | 0 | 0 | 1 | 1 | 1 |
+    +---+---+---+---+---+---+
+      proc0 | proc1 | proc2
+
 
 SourceFiles
     regionSplit.C
@@ -36,14 +95,17 @@ SourceFiles
 #ifndef regionSplit_H
 #define regionSplit_H
 
-#include "polyMesh.H"
+#include "globalIndex.H"
-#include "demandDrivenData.H"
+#include "labelPair.H"
+#include "boolList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
+class polyMesh;
+
 /*---------------------------------------------------------------------------*\
                            Class regionSplit Declaration
 \*---------------------------------------------------------------------------*/
@@ -57,8 +119,7 @@ class regionSplit
         //- Reference to mesh
         const polyMesh& mesh_;
 
-        //- Number of regions
+        mutable autoPtr<globalIndex> globalNumberingPtr_;
-        label nRegions_;
 
 
     // Private Member Functions
@@ -84,8 +145,16 @@ class regionSplit
             const label markValue
         ) const;
 
-        //- Calculate region split. Return number of regions.
+        //- Calculate local region split. Return number of regions.
-        label calcRegionSplit
+        label calcLocalRegionSplit
+        (
+            const boolList& blockedFace,
+            const List<labelPair>& explicitConnections,
+            labelList& cellRegion
+        ) const;
+
+        //- Calculate global region split. Return globalIndex.
+        autoPtr<globalIndex> calcRegionSplit
         (
             const boolList& blockedFace,
             const List<labelPair>& explicitConnections,
@@ -118,11 +187,18 @@ public:
 
     // Member Functions
 
-        //- Return number of regions
+        //- Return global region numbering
+        const globalIndex& globalNumbering() const
+        {
+            return globalNumberingPtr_();
+        }
+
+        //- Return total number of regions
         label nRegions() const
         {
-            return nRegions_;
+            return globalNumbering().size();
         }
+
 };
 
 

src/regionModels/regionModel/regionProperties/regionProperties.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -52,4 +52,5 @@ Foam::regionProperties::regionProperties(const Time& runTime)
 Foam::regionProperties::~regionProperties()
 {}
 
+
 // ************************************************************************* //