Point based consistency check

This is necessary when one would run a dynamic refinement simulation with more
than 2 refinement levels. Ordinary face checking can produce 8:1 point
inconsistency, even on hexahedral meshes. Therefore, if the refinemement level
is > 2, point based consistency is switched on.

src/dynamicMesh/dynamicFvMesh/dynamicRefinePolyFvMesh/dynamicRefinePolyFvMesh.C

@@ -611,15 +611,23 @@ labelList dynamicRefinePolyFvMesh::selectRefineCells
         }
     }
 
-    // Guarantee 2:1 refinement after refinement
-    labelList consistentSet
+    // Guarantee 2:1 refinement after refinement.
+
+    // Create a label list for consistent set
+    labelList consistentSet;
+
+    // Note: the return type of consistentRefinement is Xfer<labelList>, so
+    // we will transfer its contents into this list
+    labelList transferredConsistentSet
     (
         meshCutter_.consistentRefinement
         (
             candidates.shrink(),
-            true               // Add to set to guarantee 2:1
+            true,                 // Add to set to guarantee 2:1,
+            pointBasedRefinement_ // Whether to use point based refinement
         )
     );
+    consistentSet.transfer(transferredConsistentSet);
 
     Info<< "Selected " << returnReduce(consistentSet.size(), sumOp<label>())
         << " cells for refinement out of " << globalData().nTotalCells()
@@ -677,9 +685,11 @@ labelList dynamicRefinePolyFvMesh::selectUnrefinePoints
         meshCutter_.consistentUnrefinement
         (
             newSplitPoints,
-            false
+            false,                // Remove from the set
+            pointBasedRefinement_ // Whether to use point based unrefinement
         )
     );
+
     Info<< "Selected " << returnReduce(consistentSet.size(), sumOp<label>())
         << " split points out of a possible "
         << returnReduce(splitPoints.size(), sumOp<label>())
@@ -750,6 +760,7 @@ dynamicRefinePolyFvMesh::dynamicRefinePolyFvMesh(const IOobject& io)
         .lookupOrDefault<Switch>("singleMotionUpdate", true)
     ),
     curTimeIndex_(-1),
+    pointBasedRefinement_(false), // Set in the update() member function
     meshCutter_(*this),
     dumpLevel_(false),
     nRefinementIterations_(0)
@@ -848,6 +859,42 @@ bool dynamicRefinePolyFvMesh::update()
                 << exit(FatalError);
         }
 
+        // Determining preferred consistency refinement based on maxRefinement
+        if (maxRefinement > 2)
+        {
+            // Refinement level higher than 2 might produce 8:1 point
+            // inconsistency, make sure that the point based consistency check
+            // is used
+            pointBasedRefinement_ = true;
+        }
+
+        // Check whether the user insists on specifying refinement strategy
+        if (refineDict.found("pointBasedRefinement"))
+        {
+            const Switch userPointBasedRefinement =
+                refineDict.lookup("pointBasedRefinement");
+
+            // Check whether the user insisted on not using point based checking
+            // and in case we determined that it should be used, issue a
+            // warning and let the user have a try at it
+            if (!userPointBasedRefinement && pointBasedRefinement_)
+            {
+                WarningIn("dynamicRefinePolyFvMesh::update()")
+                    << "You are insisting on using face based consistency"
+                    << " check for dynamic refinement."
+                    << nl
+                    << "Since you are allowing more than two maximum"
+                    << " refinement levels, this might produce erroneous mesh"
+                    << " due to 8:1 point conflicts."
+                    << nl
+                    << "In order to supress this message and use point based"
+                    << " consistency checks, set pointBasedRefinement to on."
+                    << endl;
+            }
+
+            pointBasedRefinement_ = userPointBasedRefinement;
+        }
+
         const word fieldName(refineDict.lookup("field"));
 
         const volScalarField& vFld = lookupObject<volScalarField>(fieldName);
@@ -861,7 +908,7 @@ bool dynamicRefinePolyFvMesh::update()
         const label nBufferLayers =
             readLabel(refineDict.lookup("nBufferLayers"));
 
-        // Cells marked for refinement or otherwise protected from unrefinement.
+        // Cells marked for refinement
         PackedBoolList refineCell(nCells());
 
         if (globalData().nTotalCells() < maxCells)

src/dynamicMesh/dynamicFvMesh/dynamicRefinePolyFvMesh/dynamicRefinePolyFvMesh.H

@@ -74,6 +74,11 @@ class dynamicRefinePolyFvMesh
             //- Helper varaible: current time index
             label curTimeIndex_;
 
+            //- Switch whether to use point based refinement (ensures 4:1 point
+            //  consistency, along with the 2:1 face consistency
+            //  consistency). Set based on max refinement level and user input
+            Switch pointBasedRefinement_;
+
 
 protected:
 

src/dynamicMesh/dynamicMesh/directTopoChange/directTopoChange/directActions/polyRef.C

@@ -1584,6 +1584,341 @@ Foam::label Foam::polyRef::faceConsistentRefinement
 }
 
 
+// Updates refineCell (cells marked for refinement) such that across all points
+// there will be 4:1 consistency after refinement.
+Foam::label Foam::polyRef::pointConsistentRefinement
+(
+    PackedList<1>& refineCell
+) const
+{
+    // Count number of changed cells
+    label nChanged = 0;
+
+    // Collect all points from cells to refine. Assume that 10% of mesh points
+    // are going to be affected to prevent excessive resizing.
+    labelHashSet pointsToConsider(mesh_.nPoints()/10);
+
+    // Get cell points
+    const labelListList& cellPoints = mesh_.cellPoints();
+
+    // Collect points
+    forAll (cellPoints, cellI)
+    {
+        // Get current points
+        const labelList& curPoints = cellPoints[cellI];
+
+        forAll (curPoints, pointI)
+        {
+            pointsToConsider.insert(curPoints[pointI]);
+        }
+    }
+
+    // Maximum cell refinement level for each point
+    labelList maxRefLevel(mesh_.nPoints(), 0);
+
+    // Get point cells
+    const labelListList& pointCells = mesh_.pointCells();
+
+    // Loop through all points and collect maximum point level for each point
+    forAllConstIter (labelHashSet, pointsToConsider, iter)
+    {
+        // Get point index
+        const label pointI = iter.key();
+
+        // Get the cells for this point
+        const labelList& curCells = pointCells[pointI];
+
+        // Find maximum refinement level for this points
+        forAll (curCells, cellI)
+        {
+            const label curCellI = curCells[cellI];
+            const label curLevel =
+                cellLevel_[curCellI] + refineCell.get(curCellI);
+
+            if (curLevel > maxRefLevel[pointI])
+            {
+                maxRefLevel[pointI] = curLevel;
+            }
+        }
+    }
+
+    // Sync maximum refinement level across coupled boundaries
+    syncTools::syncPointList
+    (
+        mesh_,
+        maxRefLevel,
+        maxEqOp<label>(),
+        0,   // Null value
+        true // Apply separation for parallel cyclics
+    );
+
+    // Now that the levels are synced, go through considered points and add
+    // cells to refine
+    forAllConstIter (labelHashSet, pointsToConsider, iter)
+    {
+        // Get point index
+        const label pointI = iter.key();
+
+        // Get the cells for this point
+        const labelList& curCells = pointCells[iter.key()];
+
+        // Loop through these point cells and set cells for refinement which
+        // would end up having refinement level smaller than maximum level - 1
+        forAll (curCells, cellI)
+        {
+            const label curCellI = curCells[cellI];
+            const label willBeRefined = refineCell.get(curCellI);
+            const label curLevel = cellLevel_[curCellI] + willBeRefined;
+
+            if (curLevel < maxRefLevel[pointI] - 1)
+            {
+                if (willBeRefined == 0)
+                {
+                    // Set the cell for refinement and increment the counter
+                    refineCell.set(curCellI, 1);
+                    ++nChanged;
+                }
+                else
+                {
+                    FatalErrorIn
+                    (
+                        "label polyRef::pointConsistentRefinement"
+                        "(PackedList<1>& refineCells) const"
+                    )   << "Cell is marked for refinement, but the 4:1 point"
+                        << " consistency cannot be ensured." << nl
+                        << "Something went wrong before this step."
+                        << endl;
+                }
+            }
+        }
+    }
+
+    return nChanged;
+}
+
+
+// Updates unrefineCell (cells marked for unrefinement) so across all faces
+// there will be 2:1 consistency after unrefinement.
+Foam::label Foam::polyRef::faceConsistentUnrefinement
+(
+    const bool maxSet,
+    PackedList<1>& unrefineCell
+) const
+{
+    label nChanged = 0;
+
+    // Internal faces.
+    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+    {
+        label own = mesh_.faceOwner()[faceI];
+        label ownLevel = cellLevel_[own] - unrefineCell.get(own);
+
+        label nei = mesh_.faceNeighbour()[faceI];
+        label neiLevel = cellLevel_[nei] - unrefineCell.get(nei);
+
+        if (ownLevel < (neiLevel - 1))
+        {
+            if (maxSet)
+            {
+                unrefineCell.set(nei, 1);
+            }
+            else
+            {
+                if (unrefineCell.get(own) == 0)
+                {
+                    FatalErrorIn("label polyRef::faceConsistentUnrefinement(..)")
+                        << "Unrefinement problem" << abort(FatalError);
+                }
+
+                unrefineCell.set(own, 0);
+            }
+
+            ++nChanged;
+        }
+        else if (neiLevel < (ownLevel - 1))
+        {
+            if (maxSet)
+            {
+                unrefineCell.set(own, 1);
+            }
+            else
+            {
+                if (unrefineCell.get(nei) == 0)
+                {
+                    FatalErrorIn("label polyRef::faceConsistentUnrefinement(..)")
+                        << "Unrefinement problem" << abort(FatalError);
+                }
+                unrefineCell.set(nei, 0);
+            }
+            ++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] = cellLevel_[own] - unrefineCell.get(own);
+    }
+
+    // Swap to neighbour
+    syncTools::swapBoundaryFaceList(mesh_, neiLevel, false);
+
+    // Now we have neighbour value see which cells need refinement
+    forAll(neiLevel, i)
+    {
+        label own = mesh_.faceOwner()[i + mesh_.nInternalFaces()];
+        label ownLevel = cellLevel_[own] - unrefineCell.get(own);
+
+        if (ownLevel < (neiLevel[i] - 1))
+        {
+            if (!maxSet)
+            {
+                if (unrefineCell.get(own) == 0)
+                {
+                    FatalErrorIn("label polyRef::faceConsistentUnrefinement(..)")
+                        << "Unrefinement problem" << abort(FatalError);
+                }
+
+                unrefineCell.set(own, 0);
+                ++nChanged;
+            }
+        }
+        else if (neiLevel[i] < (ownLevel - 1))
+        {
+            if (maxSet)
+            {
+                if (unrefineCell.get(own) == 1)
+                {
+                    FatalErrorIn("label polyRef::faceConsistentUnrefinement(..)")
+                        << "Unrefinement problem" << abort(FatalError);
+                }
+
+                unrefineCell.set(own, 1);
+                ++nChanged;
+            }
+        }
+    }
+
+    return nChanged;
+}
+
+
+// Updates unrefineCell (cells marked for unrefinement) such that across all
+// points there will be 4:1 consistency after unrefinement.
+Foam::label Foam::polyRef::pointConsistentUnrefinement
+(
+    const PackedList<1>& unrefinePoint,
+    PackedList<1>& unrefineCell
+) const
+{
+    // Count number of changed cells
+    label nChanged = 0;
+
+    // Get a dynamicList for all unrefine point candidates. Assume that 10% of
+    // mesh points are going to be affected to prevent excessive resizing
+    dynamicLabelList pointsToConsider(mesh_.nPoints()/10);
+
+    forAll (unrefinePoint, pointI)
+    {
+        if (unrefinePoint.get(pointI) == 1)
+        {
+            pointsToConsider.append(pointI);
+        }
+    }
+
+    // Minimum cell refinement level for each point
+    labelList minRefLevel(mesh_.nPoints(), 0);
+
+    // Get point cells
+    const labelListList& pointCells = mesh_.pointCells();
+
+    // Loop through all points and collect minimum point level for each point
+    forAll (pointsToConsider, i)
+    {
+        // Get point index
+        const label pointI = pointsToConsider[i];
+
+        // Get cells for this point
+        const labelList& curCells = pointCells[pointI];
+
+        // Find minimum refinement level for this points
+        forAll (curCells, cellI)
+        {
+            const label curCellI = curCells[cellI];
+            const label curLevel =
+                cellLevel_[curCellI] - unrefineCell.get(curCellI);
+
+            if (curLevel < minRefLevel[pointI])
+            {
+                minRefLevel[pointI] = curLevel;
+            }
+        }
+    }
+
+    // Sync minimum refinement level across coupled boundaries
+    syncTools::syncPointList
+    (
+        mesh_,
+        minRefLevel,
+        minEqOp<label>(),
+        0,   // Null value
+        true // Apply separation for parallel cyclics
+    );
+
+    // Now that the levels are synced, go through considered points and add
+    // cells to unrefine
+    forAll (pointsToConsider, i)
+    {
+        // Get point index
+        const label pointI = pointsToConsider[i];
+
+        // Get the cells for this point
+        const labelList& curCells = pointCells[pointI];
+
+        // Loop through these point cells and set cells for unrefinement which
+        // would end up having refinement level greater than level + 1
+        forAll (curCells, cellI)
+        {
+            const label curCellI = curCells[cellI];
+            const label willBeUnrefined = unrefineCell.get(curCellI);
+            const label curLevel = cellLevel_[curCellI] - willBeUnrefined;
+
+            if (curLevel > minRefLevel[pointI] + 1)
+            {
+                if (willBeUnrefined == 0)
+                {
+                    // Set the cell for unrefinement and increment the counter
+                    unrefineCell.set(curCellI, 1);
+                    ++nChanged;
+                }
+                else
+                {
+                    FatalErrorIn
+                    (
+                        "label polyRef::pointConsistentUnrefinement"
+                        "\n("
+                        "\n    const PackedList<1>& unrefinePoints,"
+                        "\n    PackedList<1>& unrefineCells"
+                        "\n) const"
+                    )   << "Cell is marked for unrefinement, but the 4:1 point"
+                        << " consistency cannot be ensured." << nl
+                        << "Something went wrong before this step."
+                        << endl;
+                }
+            }
+        }
+    }
+
+    return nChanged;
+}
+
+
 // Debug: check if wanted refinement is compatible with 2:1
 void Foam::polyRef::checkWantedRefinementLevels
 (
@@ -1981,14 +2316,15 @@ Foam::labelList Foam::polyRef::cellPoints(const label cellI) const
 }
 
 
-Foam::labelList Foam::polyRef::consistentRefinement
+Foam::Xfer<Foam::labelList> Foam::polyRef::consistentRefinement
 (
     const labelList& cellsToRefine,
-    const bool maxSet
+    const bool maxSet,
+    const bool pointBasedRefinement
 ) const
 {
-    // Loop, modifying cellsToRefine, until no more changes to due to 2:1
-    // conflicts.
+    // Loop, modifying cellsToRefine, until no more changes to due to 2:1 face
+    // conflicts and optionally 4:1 point conflicts.
     // maxSet = false : unselect cells to refine
     // maxSet = true  : select cells to refine
 
@@ -2001,7 +2337,16 @@ Foam::labelList Foam::polyRef::consistentRefinement
 
     while (true)
     {
-        label nChanged = faceConsistentRefinement(maxSet, refineCell);
+        label nChanged = 0;
+
+        if (pointBasedRefinement)
+        {
+            // Check for 4:1 point based consistent unrefinement
+            nChanged += pointConsistentRefinement(refineCell);
+        }
+
+        // Check for 2:1 face based consistent unrefinement
+        nChanged += faceConsistentRefinement(maxSet, refineCell);
 
         reduce(nChanged, sumOp<label>());
 
@@ -2046,7 +2391,7 @@ Foam::labelList Foam::polyRef::consistentRefinement
         checkWantedRefinementLevels(newCellsToRefine);
     }
 
-    return newCellsToRefine;
+    return xferMove<labelList>(newCellsToRefine);
 }
 
 
@@ -2056,7 +2401,7 @@ Foam::labelList Foam::polyRef::consistentRefinement
 // - satisfies maxPointDiff (e.g. 4:1) across selected point connected
 //   cells. This is used to ensure that e.g. cells on the surface are not
 //   point connected to cells which are 8 times smaller.
-Foam::labelList Foam::polyRef::consistentSlowRefinement
+Foam::Xfer<Foam::labelList> Foam::polyRef::consistentSlowRefinement
 (
     const label maxFaceDiff,
     const labelList& cellsToRefine,
@@ -2542,11 +2887,11 @@ Foam::labelList Foam::polyRef::consistentSlowRefinement
             << " cells to refine." << endl;
     }
 
-    return newCellsToRefine;
+    return xferMove<labelList>(newCellsToRefine);
 }
 
 
-Foam::labelList Foam::polyRef::consistentSlowRefinement2
+Foam::Xfer<Foam::labelList> Foam::polyRef::consistentSlowRefinement2
 (
     const label maxFaceDiff,
     const labelList& cellsToRefine,
@@ -2974,7 +3319,7 @@ Foam::labelList Foam::polyRef::consistentSlowRefinement2
         }
     }
 
-    return newCellsToRefine;
+    return xferMove<labelList>(newCellsToRefine);
 }
 
 
@@ -4990,7 +5335,8 @@ Foam::labelList Foam::polyRef::getSplitPoints() const
 Foam::labelList Foam::polyRef::consistentUnrefinement
 (
     const labelList& pointsToUnrefine,
-    const bool maxSet
+    const bool maxSet,
+    const bool pointBasedUnrefinement
 ) const
 {
     if (debug)
@@ -5047,109 +5393,15 @@ Foam::labelList Foam::polyRef::consistentUnrefinement
 
         label nChanged = 0;
 
-
-        // Check 2:1 consistency taking refinement into account
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-        // Internal faces.
-        for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+        if (pointBasedUnrefinement)
         {
-            label own = mesh_.faceOwner()[faceI];
-            label ownLevel = cellLevel_[own] - unrefineCell.get(own);
-
-            label nei = mesh_.faceNeighbour()[faceI];
-            label neiLevel = cellLevel_[nei] - unrefineCell.get(nei);
-
-            if (ownLevel < (neiLevel-1))
-            {
-                // Since was 2:1 this can only occur if own is marked for
-                // unrefinement.
-
-                if (maxSet)
-                {
-                    unrefineCell.set(nei, 1);
-                }
-                else
-                {
-                    if (unrefineCell.get(own) == 0)
-                    {
-                        FatalErrorIn("polyRef::consistentUnrefinement(..)")
-                            << "problem" << abort(FatalError);
-                    }
-
-                    unrefineCell.set(own, 0);
-                }
-                nChanged++;
-            }
-            else if (neiLevel < (ownLevel-1))
-            {
-                if (maxSet)
-                {
-                    unrefineCell.set(own, 1);
-                }
-                else
-                {
-                    if (unrefineCell.get(nei) == 0)
-                    {
-                        FatalErrorIn("polyRef::consistentUnrefinement(..)")
-                            << "problem" << abort(FatalError);
-                    }
-
-                    unrefineCell.set(nei, 0);
-                }
-                nChanged++;
-            }
+            // Check for 4:1 point based consistent unrefinement
+            nChanged +=
+                pointConsistentUnrefinement(unrefinePoint, unrefineCell);
         }
 
-
-        // Coupled faces. Swap owner level to get neighbouring cell level.
-        labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
-
-        forAll(neiLevel, i)
-        {
-            label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
-
-            neiLevel[i] = cellLevel_[own] - unrefineCell.get(own);
-        }
-
-        // Swap to neighbour
-        syncTools::swapBoundaryFaceList(mesh_, neiLevel, false);
-
-        forAll(neiLevel, i)
-        {
-            label faceI = i+mesh_.nInternalFaces();
-            label own = mesh_.faceOwner()[faceI];
-            label ownLevel = cellLevel_[own] - unrefineCell.get(own);
-
-            if (ownLevel < (neiLevel[i]-1))
-            {
-                if (!maxSet)
-                {
-                    if (unrefineCell.get(own) == 0)
-                    {
-                        FatalErrorIn("polyRef::consistentUnrefinement(..)")
-                            << "problem" << abort(FatalError);
-                    }
-
-                    unrefineCell.set(own, 0);
-                    nChanged++;
-                }
-            }
-            else if (neiLevel[i] < (ownLevel-1))
-            {
-                if (maxSet)
-                {
-                    if (unrefineCell.get(own) == 1)
-                    {
-                        FatalErrorIn("polyRef::consistentUnrefinement(..)")
-                            << "problem" << abort(FatalError);
-                    }
-
-                    unrefineCell.set(own, 1);
-                    nChanged++;
-                }
-            }
-        }
+        // Check for 2:1 face based consistent unrefinement
+        nChanged += faceConsistentUnrefinement(maxSet, unrefineCell);
 
         reduce(nChanged, sumOp<label>());
 

src/dynamicMesh/dynamicMesh/directTopoChange/directTopoChange/directActions/polyRef.H

@@ -308,14 +308,37 @@ class polyRef
             dynamicLabelList& faceVerts
         ) const;
 
-        //- Updates refineCell so consistent 2:1 refinement. Returns local
-        //  number of cells changed
+        //- Updates refineCell such that a face consistent 2:1 refinement is
+        //  obtained. Returns local number of cells changed
         label faceConsistentRefinement
         (
             const bool maxSet,
             PackedList<1>& refineCell
         ) const;
 
+        //- Updates refineCell such that a point consistent 4:1 refinement is
+        //  obtained. Returns local number of cells changed
+        label pointConsistentRefinement
+        (
+            PackedList<1>& refineCell
+        ) const;
+
+        //- Updates unrefineCell such that a face consistent 2:1 unrefinement is
+        //  obtained. Returns local number of cells changed
+        label faceConsistentUnrefinement
+        (
+            const bool maxSet,
+            PackedList<1>& unrefineCell
+        ) const;
+
+        //- Updates unrefineCell such that a point consistent 4:1 unrefinement
+        //  is obtained. Returns local number of cells changed
+        label pointConsistentUnrefinement
+        (
+            const PackedList<1>& unrefinePoints,
+            PackedList<1>& unrefineCell
+        ) const;
+
         //- Check wanted refinement for 2:1 consistency
         void checkWantedRefinementLevels(const labelList&) const;
 
@@ -398,10 +421,11 @@ public:
             //  ok list of cells to refine.
             //  Either adds cells to refine to set (maxSet = true) or
             //  removes cells to refine (maxSet = false)
-            labelList consistentRefinement
+            Xfer<labelList> consistentRefinement
             (
                 const labelList& cellsToRefine,
-                const bool maxSet
+                const bool maxSet,
+                const bool pointBasedRefinement
             ) const;
 
             //- Like consistentRefinement but slower:
@@ -417,7 +441,7 @@ public:
             //  facesToCheck  : additional faces where to implement the
             //                  maxFaceDiff thickness (usually only boundary
             //                  faces)
-            labelList consistentSlowRefinement
+            Xfer<labelList> consistentSlowRefinement
             (
                 const label maxFaceDiff,
                 const labelList& cellsToRefine,
@@ -429,7 +453,7 @@ public:
             //- Like consistentSlowRefinement but uses different meshWave
             //  (proper distance instead of toplogical count). No point checks
             //  yet
-            labelList consistentSlowRefinement2
+            Xfer<labelList> consistentSlowRefinement2
             (
                 const label maxFaceDiff,
                 const labelList& cellsToRefine,
@@ -523,7 +547,8 @@ public:
             labelList consistentUnrefinement
             (
                 const labelList& pointsToUnrefine,
-                const bool maxSet
+                const bool maxSet,
+                const bool pointBasedUnrefinement
             ) const;
 
             //- Remove some refinement. Needs to be supplied output of