From 01b785f23a6dceb90eff6f981a26e4cca1c6b8f9 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 21 Dec 2017 13:12:22 +0000
Subject: [PATCH] ENH: 2:1 constraint

---
 .../polyTopoChange/hexRef8/hexRef8.C          |  44 ++--
 .../polyTopoChange/hexRef8/hexRef8.H          |  22 +-
 .../meshRefinement/meshRefinement.H           |  10 +
 .../meshRefinement/meshRefinementRefine.C     | 102 +++++++-
 .../snappyHexMeshDriver/snappyRefineDriver.C  | 242 ++----------------
 5 files changed, 176 insertions(+), 244 deletions(-)

diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C
index 74dcafca1d..79c7e7446a 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C
@@ -1564,6 +1564,7 @@ void Foam::hexRef8::walkFaceFromMid
 Foam::label Foam::hexRef8::faceConsistentRefinement
 (
     const bool maxSet,
+    const labelUList& cellLevel,
     bitSet& refineCell
 ) const
 {
@@ -1573,10 +1574,10 @@ Foam::label Foam::hexRef8::faceConsistentRefinement
     for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
     {
         label own = mesh_.faceOwner()[facei];
-        label nei = mesh_.faceNeighbour()[facei];
+        label ownLevel = cellLevel[own] + refineCell.get(own);
 
-        label ownLevel = cellLevel_[own] + refineCell.get(own);
-        label neiLevel = cellLevel_[nei] + refineCell.get(nei);
+        label nei = mesh_.faceNeighbour()[facei];
+        label neiLevel = cellLevel[nei] + refineCell.get(nei);
 
         if (ownLevel > (neiLevel+1))
         {
@@ -1613,7 +1614,7 @@ Foam::label Foam::hexRef8::faceConsistentRefinement
     {
         label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
 
-        neiLevel[i] = cellLevel_[own] + refineCell.get(own);
+        neiLevel[i] = cellLevel[own] + refineCell.get(own);
     }
 
     // Swap to neighbour
@@ -1623,7 +1624,7 @@ Foam::label Foam::hexRef8::faceConsistentRefinement
     forAll(neiLevel, i)
     {
         label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
-        label ownLevel = cellLevel_[own] + refineCell.get(own);
+        label ownLevel = cellLevel[own] + refineCell.get(own);
 
         if (ownLevel > (neiLevel[i]+1))
         {
@@ -1650,6 +1651,7 @@ Foam::label Foam::hexRef8::faceConsistentRefinement
 // Debug: check if wanted refinement is compatible with 2:1
 void Foam::hexRef8::checkWantedRefinementLevels
 (
+    const labelUList& cellLevel,
     const labelList& cellsToRefine
 ) const
 {
@@ -1658,10 +1660,10 @@ void Foam::hexRef8::checkWantedRefinementLevels
     for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
     {
         label own = mesh_.faceOwner()[facei];
-        label nei = mesh_.faceNeighbour()[facei];
+        label ownLevel = cellLevel[own] + refineCell.get(own);
 
-        label ownLevel = cellLevel_[own] + refineCell.get(own);
-        label neiLevel = cellLevel_[nei] + refineCell.get(nei);
+        label nei = mesh_.faceNeighbour()[facei];
+        label neiLevel = cellLevel[nei] + refineCell.get(nei);
 
         if (mag(ownLevel-neiLevel) > 1)
         {
@@ -1669,11 +1671,11 @@ void Foam::hexRef8::checkWantedRefinementLevels
             dumpCell(nei);
             FatalErrorInFunction
                 << "cell:" << own
-                << " current level:" << cellLevel_[own]
+                << " current level:" << cellLevel[own]
                 << " level after refinement:" << ownLevel
                 << nl
                 << "neighbour cell:" << nei
-                << " current level:" << cellLevel_[nei]
+                << " current level:" << cellLevel[nei]
                 << " level after refinement:" << neiLevel
                 << nl
                 << "which does not satisfy 2:1 constraints anymore."
@@ -1689,7 +1691,7 @@ void Foam::hexRef8::checkWantedRefinementLevels
     {
         label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
 
-        neiLevel[i] = cellLevel_[own] + refineCell.get(own);
+        neiLevel[i] = cellLevel[own] + refineCell.get(own);
     }
 
     // Swap to neighbour
@@ -1701,7 +1703,7 @@ void Foam::hexRef8::checkWantedRefinementLevels
         label facei = i + mesh_.nInternalFaces();
 
         label own = mesh_.faceOwner()[facei];
-        label ownLevel = cellLevel_[own] + refineCell.get(own);
+        label ownLevel = cellLevel[own] + refineCell.get(own);
 
         if (mag(ownLevel - neiLevel[i]) > 1)
         {
@@ -1715,7 +1717,7 @@ void Foam::hexRef8::checkWantedRefinementLevels
                 << " on patch " << patchi << " "
                 << mesh_.boundaryMesh()[patchi].name()
                 << " owner cell " << own
-                << " current level:" << cellLevel_[own]
+                << " current level:" << cellLevel[own]
                 << " level after refinement:" << ownLevel
                 << nl
                 << " (coupled) neighbour cell will get refinement "
@@ -2251,6 +2253,7 @@ Foam::hexRef8::hexRef8
 
 Foam::labelList Foam::hexRef8::consistentRefinement
 (
+    const labelUList& cellLevel,
     const labelList& cellsToRefine,
     const bool maxSet
 ) const
@@ -2264,7 +2267,12 @@ Foam::labelList Foam::hexRef8::consistentRefinement
 
     while (true)
     {
-        label nChanged = faceConsistentRefinement(maxSet, refineCell);
+        label nChanged = faceConsistentRefinement
+        (
+            maxSet,
+            cellLevel,
+            refineCell
+        );
 
         reduce(nChanged, sumOp<label>());
 
@@ -2286,7 +2294,7 @@ Foam::labelList Foam::hexRef8::consistentRefinement
 
     if (debug)
     {
-        checkWantedRefinementLevels(newCellsToRefine);
+        checkWantedRefinementLevels(cellLevel, newCellsToRefine);
     }
 
     return newCellsToRefine;
@@ -3089,11 +3097,11 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2
             refineCell.set(celli);
         }
     }
-    faceConsistentRefinement(true, refineCell);
+    faceConsistentRefinement(true, cellLevel_, refineCell);
 
     while (true)
     {
-        label nChanged = faceConsistentRefinement(true, refineCell);
+        label nChanged = faceConsistentRefinement(true, cellLevel_, refineCell);
 
         reduce(nChanged, sumOp<label>());
 
@@ -3141,7 +3149,7 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2
 
         const bitSet savedRefineCell(refineCell);
 
-        label nChanged = faceConsistentRefinement(true, refineCell);
+        label nChanged = faceConsistentRefinement(true, cellLevel_, refineCell);
 
         {
             cellSet cellsOut2
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.H b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.H
index a2388f3b60..3f4292d527 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.H
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.H
@@ -300,11 +300,16 @@ class hexRef8
         label faceConsistentRefinement
         (
             const bool maxSet,
+            const labelUList& cellLevel,
             bitSet& refineCell
         ) const;
 
         //- Check wanted refinement for 2:1 consistency
-        void checkWantedRefinementLevels(const labelList&) const;
+        void checkWantedRefinementLevels
+        (
+            const labelUList& cellLevel,
+            const labelList&
+        ) const;
 
 
         // Cellshape recognition
@@ -420,10 +425,25 @@ public:
             //  removes cells to refine (maxSet = false)
             labelList consistentRefinement
             (
+                const labelUList& cellLevel,
                 const labelList& cellsToRefine,
                 const bool maxSet
             ) const;
 
+            //- Given valid mesh and current cell level and proposed
+            //  cells to refine calculate any clashes (due to 2:1) and return
+            //  ok list of cells to refine.
+            //  Either adds cells to refine to set (maxSet = true) or
+            //  removes cells to refine (maxSet = false)
+            labelList consistentRefinement
+            (
+                const labelList& cellsToRefine,
+                const bool maxSet
+            ) const
+            {
+                return consistentRefinement(cellLevel_, cellsToRefine, maxSet);
+            }
+
             //- Like consistentRefinement but slower:
             //
             //  - specify number of cells between consecutive refinement levels
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H
index 1d9e724825..2ae2f5634c 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H
@@ -1037,6 +1037,16 @@ public:
                 const scalar maxLoadUnbalance
             );
 
+
+            //- Calculate list of cells to directionally refine
+            labelList directionalRefineCandidates
+            (
+                const label maxGlobalCells,
+                const label maxLocalCells,
+                const labelList& currentLevel,
+                const direction dir
+            ) const;
+
             //- Directionally refine in direction cmpt
             autoPtr<mapPolyMesh> directionalRefine
             (
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C
index 53728b952a..9954987315 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C
@@ -2078,7 +2078,7 @@ Foam::labelList Foam::meshRefinement::refineCandidates
         label nAllowRefine = labelMax / Pstream::nProcs();
 
         // Marked for refinement (>= 0) or not (-1). Actual value is the
-        // index of the surface it intersects.
+        // index of the surface it intersects / shell it is inside.
         labelList refineCell(mesh_.nCells(), -1);
         label nRefine = 0;
 
@@ -2629,6 +2629,106 @@ Foam::meshRefinement::balanceAndRefine
 }
 
 
+Foam::labelList Foam::meshRefinement::directionalRefineCandidates
+(
+    const label maxGlobalCells,
+    const label maxLocalCells,
+    const labelList& currentLevel,
+    const direction dir
+) const
+{
+    const labelList& cellLevel = meshCutter_.cellLevel();
+    const pointField& cellCentres = mesh_.cellCentres();
+
+    label totNCells = mesh_.globalData().nTotalCells();
+
+    labelList cellsToRefine;
+
+    if (totNCells >= maxGlobalCells)
+    {
+        Info<< "No cells marked for refinement since reached limit "
+            << maxGlobalCells << '.' << endl;
+    }
+    else
+    {
+        // Disable refinement shortcut. nAllowRefine is per processor limit.
+        label nAllowRefine = labelMax / Pstream::nProcs();
+
+        // Marked for refinement (>= 0) or not (-1). Actual value is the
+        // index of the surface it intersects / shell it is inside
+        labelList refineCell(mesh_.nCells(), -1);
+        label nRefine = 0;
+
+        // Find cells inside the shells with directional levels
+        labelList insideShell;
+        shells_.findDirectionalLevel
+        (
+            cellCentres,
+            cellLevel,
+            currentLevel,  // current directional level
+            dir,
+            insideShell
+        );
+
+        // Mark for refinement
+        forAll(insideShell, celli)
+        {
+            if (insideShell[celli] >= 0)
+            {
+                bool reachedLimit = !markForRefine
+                (
+                    insideShell[celli],    // mark with any positive value
+                    nAllowRefine,
+                    refineCell[celli],
+                    nRefine
+                );
+
+                if (reachedLimit)
+                {
+                    if (debug)
+                    {
+                        Pout<< "Stopped refining cells"
+                            << " since reaching my cell limit of "
+                            << mesh_.nCells()+nAllowRefine << endl;
+                    }
+                    break;
+                }
+            }
+        }
+
+        // Limit refinement
+        // ~~~~~~~~~~~~~~~~
+
+        {
+            label nUnmarked = unmarkInternalRefinement(refineCell, nRefine);
+            if (nUnmarked > 0)
+            {
+                Info<< "Unmarked for refinement due to limit shells"
+                    << "                : " << nUnmarked << " cells."  << endl;
+            }
+        }
+
+
+
+        // Pack cells-to-refine
+        // ~~~~~~~~~~~~~~~~~~~~
+
+        cellsToRefine.setSize(nRefine);
+        nRefine = 0;
+
+        forAll(refineCell, cellI)
+        {
+            if (refineCell[cellI] != -1)
+            {
+                cellsToRefine[nRefine++] = cellI;
+            }
+        }
+    }
+
+    return cellsToRefine;
+}
+
+
 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::directionalRefine
 (
     const string& msg,
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C
index 68df871bd5..15fe6edf80 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C
@@ -1572,190 +1572,6 @@ Foam::label Foam::snappyRefineDriver::shellRefine
 
     return iter;
 }
-//XXXXXXXXXX
-Foam::List<Foam::direction> Foam::snappyRefineDriver::faceDirection() const
-{
-    const fvMesh& mesh = meshRefiner_.mesh();
-
-    List<direction> faceDir(mesh.nFaces());
-
-    vectorField nf(mesh.faceAreas());
-    nf /= mag(nf);
-
-    forAll(nf, facei)
-    {
-        const vector& n = nf[facei];
-        if (mag(n[0]) > 0.5)
-        {
-            faceDir[facei] = 0;
-        }
-        else if (mag(n[1]) > 0.5)
-        {
-            faceDir[facei] = 1;
-        }
-        else if (mag(n[2]) > 0.5)
-        {
-            faceDir[facei] = 2;
-        }
-    }
-
-    return faceDir;
-}
-Foam::label Foam::snappyRefineDriver::faceConsistentRefinement
-(
-    const PackedBoolList& isDirFace,
-    const List<labelVector>& dirCellLevel,
-    const direction dir,
-    const bool maxSet,
-    PackedBoolList& refineCell
-) const
-{
-    const fvMesh& mesh = meshRefiner_.mesh();
-
-    label nChanged = 0;
-
-    // Internal faces.
-    for (label facei = 0; facei < mesh.nInternalFaces(); facei++)
-    {
-        if (isDirFace[facei])
-        {
-            label own = mesh.faceOwner()[facei];
-            label ownLevel = dirCellLevel[own][dir] + refineCell.get(own);
-
-            label nei = mesh.faceNeighbour()[facei];
-            label neiLevel = dirCellLevel[nei][dir] + refineCell.get(nei);
-
-            if (ownLevel > (neiLevel+1))
-            {
-                if (maxSet)
-                {
-                    refineCell.set(nei);
-                }
-                else
-                {
-                    refineCell.unset(own);
-                }
-                nChanged++;
-            }
-            else if (neiLevel > (ownLevel+1))
-            {
-                if (maxSet)
-                {
-                    refineCell.set(own);
-                }
-                else
-                {
-                    refineCell.unset(nei);
-                }
-                nChanged++;
-            }
-        }
-    }
-
-
-    // Coupled faces. Swap owner level to get neighbouring cell level.
-    // (only boundary faces of neiLevel used)
-    labelList neiLevel(mesh.nFaces()-mesh.nInternalFaces());
-
-    forAll(neiLevel, i)
-    {
-        label own = mesh.faceOwner()[i+mesh.nInternalFaces()];
-
-        neiLevel[i] = dirCellLevel[own][dir] + refineCell.get(own);
-    }
-
-    // Swap to neighbour
-    syncTools::swapBoundaryFaceList(mesh, neiLevel);
-
-    // Now we have neighbour value see which cells need refinement
-    forAll(neiLevel, i)
-    {
-        label facei = i+mesh.nInternalFaces();
-        if (isDirFace[facei])
-        {
-            label own = mesh.faceOwner()[facei];
-            label ownLevel = dirCellLevel[own][dir] + refineCell.get(own);
-
-            if (ownLevel > (neiLevel[i]+1))
-            {
-                if (!maxSet)
-                {
-                    refineCell.unset(own);
-                    nChanged++;
-                }
-            }
-            else if (neiLevel[i] > (ownLevel+1))
-            {
-                if (maxSet)
-                {
-                    refineCell.set(own);
-                    nChanged++;
-                }
-            }
-        }
-    }
-
-    return nChanged;
-}
-Foam::labelList Foam::snappyRefineDriver::consistentDirectionalRefinement
-(
-    const direction dir,
-    const List<labelVector>& dirCellLevel,
-    const labelUList& candidateCells,
-    const bool maxSet
-) const
-{
-    const fvMesh& mesh = meshRefiner_.mesh();
-
-    const List<direction> faceDir(faceDirection());
-
-    PackedBoolList isDirFace(mesh.nFaces());
-    forAll(faceDir, facei)
-    {
-        if (faceDir[facei] == dir)
-        {
-            isDirFace[facei] = true;
-        }
-    }
-
-    // Loop, modifying cellsToRefine, until no more changes to due to 2:1
-    // conflicts.
-    // maxSet = false : unselect cells to refine
-    // maxSet = true  : select cells to refine
-
-    // Go to straight boolList.
-    PackedBoolList isRefineCell(mesh.nCells());
-    isRefineCell.set(candidateCells);
-
-    while (true)
-    {
-        label nChanged = faceConsistentRefinement
-        (
-            isDirFace,
-            dirCellLevel,
-            dir,
-            maxSet,
-            isRefineCell
-        );
-
-        reduce(nChanged, sumOp<label>());
-
-        if (debug)
-        {
-            Pout<< "snappyRefineDriver::consistentDirectionalRefinement"
-                << " : Changed " << nChanged
-                << " refinement levels due to 2:1 conflicts."
-                << endl;
-        }
-
-        if (nChanged == 0)
-        {
-            break;
-        }
-    }
-
-    return isRefineCell.used();
-}
 
 
 Foam::label Foam::snappyRefineDriver::directionalShellRefine
@@ -1814,53 +1630,31 @@ Foam::label Foam::snappyRefineDriver::directionalShellRefine
             // - original cellLevel (using mapping) mentioned in levelIncrement
             // - dirCellLevel not yet up to cellLevel+levelIncrement
 
-            labelList candidateCells;
+
+            // Extract component of directional level
+            labelList currentLevel(dirCellLevel.size());
+            forAll(dirCellLevel, celli)
             {
-                // Extract component of directional level
-                labelList currentLevel(dirCellLevel.size());
-                forAll(dirCellLevel, celli)
-                {
-                    currentLevel[celli] = dirCellLevel[celli][dir];
-                }
-
-                // Find cells inside the shells
-                labelList insideShell;
-                shells.findDirectionalLevel
-                (
-                    mesh.cellCentres(),
-                    cellLevel,
-                    currentLevel,  // directional level
-                    dir,
-                    insideShell
-                );
-
-                // Extract
-                label n = 0;
-                forAll(insideShell, celli)
-                {
-                    if (insideShell[celli] >= 0)
-                    {
-                        n++;
-                    }
-                }
-                candidateCells.setSize(n);
-                n = 0;
-                forAll(insideShell, celli)
-                {
-                    if (insideShell[celli] >= 0)
-                    {
-                        candidateCells[n++] = celli;
-                    }
-                }
+                currentLevel[celli] = dirCellLevel[celli][dir];
             }
 
+            labelList candidateCells
+            (
+                meshRefiner_.directionalRefineCandidates
+                (
+                    refineParams.maxGlobalCells(),
+                    refineParams.maxLocalCells(),
+                    currentLevel,
+                    dir
+                )
+            );
+
             // Extend to keep 2:1 ratio
             labelList cellsToRefine
             (
-                consistentDirectionalRefinement
+                meshRefiner_.meshCutter().consistentRefinement
                 (
-                    dir,
-                    dirCellLevel,
+                    currentLevel,
                     candidateCells,
                     true
                 )