From 9ba219d0f3af3f897772afef92e3b7c53929d559 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 6 Jun 2013 06:04:00 +0100
Subject: [PATCH 1/2] ENH: MGridGen: enforce all coarse cells to be single
 connected region

---
 .../GAMGAgglomeration/GAMGAgglomeration.C     | 119 ++++++++++++++++++
 .../GAMGAgglomeration/GAMGAgglomeration.H     |  10 ++
 .../MGridGenGAMGAgglomerate.C                 |  22 +++-
 3 files changed, 150 insertions(+), 1 deletion(-)

diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C
index dd51eb645b..66a495f6b2 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C
@@ -535,4 +535,123 @@ const Foam::labelListListList& Foam::GAMGAgglomeration::boundaryFaceMap
 }
 
 
+bool Foam::GAMGAgglomeration::checkRestriction
+(
+    labelList& newRestrict,
+    label& nNewCoarse,
+    const lduAddressing& fineAddressing,
+    const labelUList& restrict,
+    const label nCoarse
+)
+{
+    if (fineAddressing.size() != restrict.size())
+    {
+        FatalErrorIn
+        (
+            "checkRestriction(..)"
+        )   << "nCells:" << fineAddressing.size()
+            << " agglom:" << restrict.size()
+            << abort(FatalError);
+    }
+
+    // Seed (master) for every region
+    labelList master(identity(fineAddressing.size()));
+
+    // Now loop and transport master through region
+    const labelUList& lower = fineAddressing.lowerAddr();
+    const labelUList& upper = fineAddressing.upperAddr();
+
+    while (true)
+    {
+        label nChanged = 0;
+
+        forAll(lower, faceI)
+        {
+            label own = lower[faceI];
+            label nei = upper[faceI];
+
+            if (restrict[own] == restrict[nei])
+            {
+                // coarse-mesh-internal face
+
+                if (master[own] < master[nei])
+                {
+                    master[nei] = master[own];
+                    nChanged++;
+                }
+                else if (master[own] > master[nei])
+                {
+                    master[own] = master[nei];
+                    nChanged++;
+                }
+            }
+        }
+
+        reduce(nChanged, sumOp<label>());
+
+        if (nChanged == 0)
+        {
+            break;
+        }
+    }
+
+
+    // Count number of regions/masters per coarse cell
+    labelListList coarseToMasters(nCoarse);
+    nNewCoarse = 0;
+    forAll(restrict, cellI)
+    {
+        labelList& masters = coarseToMasters[restrict[cellI]];
+
+        if (findIndex(masters, master[cellI]) == -1)
+        {
+            masters.append(master[cellI]);
+            nNewCoarse++;
+        }
+    }
+
+    if (nNewCoarse > nCoarse)
+    {
+        //WarningIn("GAMGAgglomeration::checkRestriction(..)")
+        //    << "Have " << nCoarse
+        //    << " agglomerated cells but " << nNewCoarse
+        //    << " disconnected regions" << endl;
+
+        // Keep coarseToMasters[0] the original coarse, allocate new ones
+        // for the others
+        labelListList coarseToNewCoarse(coarseToMasters.size());
+
+        nNewCoarse = nCoarse;
+
+        forAll(coarseToMasters, coarseI)
+        {
+            const labelList& masters = coarseToMasters[coarseI];
+
+            labelList& newCoarse = coarseToNewCoarse[coarseI];
+            newCoarse.setSize(masters.size());
+            newCoarse[0] = coarseI;
+            for (label i = 1; i < newCoarse.size(); i++)
+            {
+                newCoarse[i] = nNewCoarse++;
+            }
+        }
+
+        newRestrict.setSize(fineAddressing.size());
+        forAll(restrict, cellI)
+        {
+            label coarseI = restrict[cellI];
+
+            label index = findIndex(coarseToMasters[coarseI], master[cellI]);
+            newRestrict[cellI] = coarseToNewCoarse[coarseI][index];
+        }
+
+        return false;
+    }
+    else
+    {
+        return true;
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.H b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.H
index 4f6ad8af07..37cf4a9197 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.H
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.H
@@ -474,6 +474,16 @@ public:
             const labelListListList& boundaryFaceMap(const label fineLeveli)
             const;
 
+        //- Given restriction determines if coarse cells are connected.
+        //  Return ok is so, otherwise creates new restriction that is
+        static bool checkRestriction
+        (
+            labelList& newRestrict,
+            label& nNewCoarse,
+            const lduAddressing& fineAddressing,
+            const labelUList& restrict,
+            const label nCoarse
+        );
 };
 
 
diff --git a/src/fvAgglomerationMethods/MGridGenGamgAgglomeration/MGridGenGAMGAgglomerate.C b/src/fvAgglomerationMethods/MGridGenGamgAgglomeration/MGridGenGAMGAgglomerate.C
index 75b7e5084c..e3663315aa 100644
--- a/src/fvAgglomerationMethods/MGridGenGamgAgglomeration/MGridGenGAMGAgglomerate.C
+++ b/src/fvAgglomerationMethods/MGridGenGamgAgglomeration/MGridGenGAMGAgglomerate.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) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -156,6 +156,26 @@ Foam::tmp<Foam::labelField> Foam::MGridGenGAMGAgglomeration::agglomerate
         finalAgglom.begin()
     );
 
+    {
+        label nNewCoarseCells = 0;
+        labelList newRestrictAddr;
+        bool ok = checkRestriction
+        (
+            newRestrictAddr,
+            nNewCoarseCells
+,
+            fineAddressing,
+            finalAgglom,
+            nCoarseCells
+        );
+
+        if (!ok)
+        {
+            nCoarseCells = nNewCoarseCells;
+            finalAgglom.transfer(newRestrictAddr);
+        }
+    }
+
     return tmp<labelField>(new labelField(finalAgglom));
 }
 

From 9426345a093eea9bd390e6c16af69870518c2c9b Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 6 Jun 2013 08:59:23 +0100
Subject: [PATCH 2/2] ENH: Test-GAMGAgglomeration: moved test into
 GAMGAgglomeration

---
 .../Test-GAMGAgglomeration.C                  | 130 +++---------------
 1 file changed, 18 insertions(+), 112 deletions(-)

diff --git a/applications/test/GAMGAgglomeration/Test-GAMGAgglomeration.C b/applications/test/GAMGAgglomeration/Test-GAMGAgglomeration.C
index a11b78772e..2dd7905294 100644
--- a/applications/test/GAMGAgglomeration/Test-GAMGAgglomeration.C
+++ b/applications/test/GAMGAgglomeration/Test-GAMGAgglomeration.C
@@ -35,116 +35,6 @@ Description
 #include "meshTools.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-void checkConnectedAgglomeration
-(
-    const lduMesh& mesh,
-    const labelUList& restrict,
-    const label nCoarse
-)
-{
-    if (mesh.lduAddr().size() != restrict.size())
-    {
-        FatalErrorIn
-        (
-            "checkConnectedAgglomeration(const lduMesh&, const labelList&)"
-        )   << "nCells:" << mesh.lduAddr().size()
-            << " agglom:" << restrict.size()
-            << abort(FatalError);
-    }
-
-    // Seed (master) for every region
-    labelList regionToMaster(nCoarse, -1);
-    labelList master(mesh.lduAddr().size(), -1);
-    forAll(restrict, cellI)
-    {
-        label region = restrict[cellI];
-        if (regionToMaster[region] == -1)
-        {
-            // Set cell to be master for region
-            //Pout<< "For region " << region
-            //    << " allocating local master " << cellI
-            //    << endl;
-            regionToMaster[region] = cellI;
-            master[cellI] = cellI;
-        }
-    }
-
-    // Now loop and transport master through region
-    const labelUList& lower = mesh.lduAddr().lowerAddr();
-    const labelUList& upper = mesh.lduAddr().upperAddr();
-
-    while (true)
-    {
-        label nChanged = 0;
-
-        forAll(lower, faceI)
-        {
-            label own = lower[faceI];
-            label nei = upper[faceI];
-
-            if (restrict[own] == restrict[nei])
-            {
-                // Region-internal face
-
-                if (master[own] != -1)
-                {
-                    if (master[nei] == -1)
-                    {
-                        master[nei] = master[own];
-                        nChanged++;
-                    }
-                    else if (master[nei] != master[own])
-                    {
-                        FatalErrorIn("checkConnectedAgglomeration(..)")
-                            << "problem" << abort(FatalError);
-                    }
-                }
-                else if (master[nei] != -1)
-                {
-                    master[own] = master[nei];
-                    nChanged++;
-                }
-            }
-        }
-
-        reduce(nChanged, sumOp<label>());
-
-        if (nChanged == 0)
-        {
-            break;
-        }
-    }
-
-    // Check that master is set for all cells
-    boolList singleRegion(nCoarse, true);
-    label nSet = nCoarse;
-    forAll(master, cellI)
-    {
-        if (master[cellI] == -1)
-        {
-            label region = restrict[cellI];
-            if (singleRegion[region] == true)
-            {
-                singleRegion[region] = false;
-                nSet--;
-            }
-        }
-    }
-
-    label totalNCoarse = returnReduce(nCoarse, sumOp<label>());
-    label totalNVisited = returnReduce(nSet, sumOp<label>());
-
-    if (totalNVisited < totalNCoarse)
-    {
-        WarningIn("checkConnectedAgglomeration(..)")
-            << "out of " << totalNCoarse
-            << " agglomerated cells have " << totalNCoarse-totalNVisited
-            << " cells that are not a single connected region" << endl;
-    }
-}
-
-
 // Main program:
 
 int main(int argc, char *argv[])
@@ -227,12 +117,28 @@ int main(int argc, char *argv[])
             << "    agglomerated size : "
             << returnReduce(coarseSize, sumOp<label>()) << endl;
 
-        checkConnectedAgglomeration
+        labelList newAddr;
+        label newCoarseSize = 0;
+        bool ok = GAMGAgglomeration::checkRestriction
         (
-            agglom.meshLevel(level),
+            newAddr,
+            newCoarseSize,
+
+            agglom.meshLevel(level).lduAddr(),
             addr,
             coarseSize
         );
+        if (!ok)
+        {
+            WarningIn(args.executable())
+                << "At level " << level
+                << " there are " << coarseSize
+                << " agglomerated cells but " << newCoarseSize
+                << " disconnected regions" << endl
+                << "    This means that some agglomerations (coarse cells)"
+                << "    consist of multiple disconnected regions."
+                << endl;
+        }
 
 
         forAll(addr, fineI)