ENH: snappyHexMesh: initial version of directional refinement

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C

@@ -1082,7 +1082,29 @@ int main(int argc, char *argv[])
 
     if (!limitDict.empty())
     {
-        Info<< "Read refinement shells in = "
+        Info<< "Read limit shells in = "
+            << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
+    }
+
+
+    // Optionally read directional refinement shells
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    const dictionary dirRefDict
+    (
+        refineDict.subOrEmptyDict("directionalRefinementRegions")
+    );
+
+    if (!dirRefDict.empty())
+    {
+        Info<< "Reading directional refinement shells." << endl;
+    }
+
+    shellSurfaces dirShells(allGeometry, dirRefDict);
+
+    if (!dirRefDict.empty())
+    {
+        Info<< "Read directional refinement shells in = "
             << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
     }
 
@@ -1119,6 +1141,7 @@ int main(int argc, char *argv[])
         surfaces,           // for surface intersection refinement
         features,           // for feature edges/point based refinement
         shells,             // for volume (inside/outside) refinement
+        dirShells,          // vol volume directional refinement
         limitShells         // limit of volume refinement
     );
     Info<< "Calculated surface intersections in = "

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict

@@ -49,6 +49,14 @@ geometry
         max (3.5 2 0.5);
     }
 
+    // Shell for directional refinement
+    refinementBox
+    {
+        type searchableBox;
+        min (1.5 1 -0.5);
+        max (3.5 2 0.5);
+    }
+
     sphere.stl
     {
         type triSurfaceMesh;
@@ -278,6 +286,21 @@ castellatedMeshControls
     }
 
 
+    directionalRefinementRegions
+    {
+        refinementBox       // Closed surface
+        {
+            mode inside;
+            levelIncrement  // Specification of additional refinement
+            (
+                (0 (1 0 0)) // For level 0 cells: add one level refinement in x
+                (1 (1 0 0)) // For level 1 cells: add one level refinement in x
+            );
+        }
+    }
+
+
+
     // Optionally limit refinement in geometric region. This limits all
     // refinement (from features, refinementSurfaces, refinementRegions)
     // in a given geometric region. The syntax is exactly the same as for the

src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C

@@ -1222,6 +1222,7 @@ Foam::meshRefinement::meshRefinement
     const refinementSurfaces& surfaces,
     const refinementFeatures& features,
     const shellSurfaces& shells,
+    const shellSurfaces& dirShells,
     const shellSurfaces& limitShells
 )
 :
@@ -1232,6 +1233,7 @@ Foam::meshRefinement::meshRefinement
     surfaces_(surfaces),
     features_(features),
     shells_(shells),
+    dirShells_(dirShells),
     limitShells_(limitShells),
     meshCutter_
     (

src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H

@@ -161,6 +161,9 @@ private:
         //- All shell-refinement interaction
         const shellSurfaces& shells_;
 
+        //- All directional shell-refinement interaction
+        const shellSurfaces& dirShells_;
+
         //- All limit-refinement interaction
         const shellSurfaces& limitShells_;
 
@@ -807,8 +810,9 @@ public:
             const bool overwrite,
             const refinementSurfaces&,
             const refinementFeatures&,
-            const shellSurfaces&,
+            const shellSurfaces&,   // omnidirectional refinement
-            const shellSurfaces&
+            const shellSurfaces&,   // directional refinement
+            const shellSurfaces&    // limit refinement
         );
 
 
@@ -861,6 +865,12 @@ public:
                 return shells_;
             }
 
+            //- Reference to directional shell-refinement shells
+            const shellSurfaces& dirShells() const
+            {
+                return dirShells_;
+            }
+
             //- Reference to meshcutting engine
             const hexRef8& meshCutter() const
             {
@@ -1037,6 +1047,14 @@ public:
                 const scalar maxLoadUnbalance
             );
 
+            //- Directionally refine in direction cmpt
+            autoPtr<mapPolyMesh> directionalRefine
+            (
+                const string& msg,
+                const direction cmpt,
+                const labelList& cellsToRefine
+            );
+
 
         // Baffle handling
 

src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C

@@ -41,6 +41,11 @@ License
 #include "cellSet.H"
 #include "treeDataCell.H"
 
+#include "cellCuts.H"
+#include "refineCell.H"
+#include "hexCellLooper.H"
+#include "meshCutter.H"
+
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
@@ -2624,4 +2629,63 @@ Foam::meshRefinement::balanceAndRefine
 }
 
 
+Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::directionalRefine
+(
+    const string& msg,
+    const direction cmpt,
+    const labelList& cellsToRefine
+)
+{
+    // Set splitting direction
+    vector refDir(Zero);
+    refDir[cmpt] = 1;
+    List<refineCell> refCells(cellsToRefine.size());
+    forAll(cellsToRefine, i)
+    {
+        refCells[i] = refineCell(cellsToRefine[i], refDir);
+    }
+
+    // How to walk circumference of cells
+    hexCellLooper cellWalker(mesh_);
+
+    // Analyse cuts
+    cellCuts cuts(mesh_, cellWalker, refCells);
+
+    // Cell cutter
+    Foam::meshCutter meshRefiner(mesh_);
+
+    polyTopoChange meshMod(mesh_);
+
+    // Insert mesh refinement into polyTopoChange.
+    meshRefiner.setRefinement(cuts, meshMod);
+
+    autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh_, false);
+
+    // Update fields
+    mesh_.updateMesh(morphMap);
+
+    // Move mesh (since morphing does not do this)
+    if (morphMap().hasMotionPoints())
+    {
+        mesh_.movePoints(morphMap().preMotionPoints());
+    }
+    else
+    {
+        // Delete mesh volumes.
+        mesh_.clearOut();
+    }
+
+    // Reset the instance for if in overwrite mode
+    mesh_.setInstance(timeName());
+
+    // Update stored refinement pattern
+    meshRefiner.updateMesh(morphMap);
+
+    // Update intersection info
+    updateMesh(morphMap, getChangedFaces(morphMap, cellsToRefine));
+
+    return morphMap;
+}
+
+
 // ************************************************************************* //

src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C

@@ -584,6 +584,7 @@ Foam::shellSurfaces::shellSurfaces
     modes_.setSize(shellI);
     distances_.setSize(shellI);
     levels_.setSize(shellI);
+    dirLevels_.setSize(shellI);
 
     extendedGapLevel_.setSize(shellI);
     extendedGapMode_.setSize(shellI);
@@ -615,6 +616,31 @@ Foam::shellSurfaces::shellSurfaces
             setAndCheckLevels(shellI, dict.lookup("levels"));
 
 
+            // Directional refinement
+            // ~~~~~~~~~~~~~~~~~~~~~~
+
+            if (dict.readIfPresent("levelIncrement", dirLevels_[shellI]))
+            {
+                if (modes_[shellI] == INSIDE)
+                {
+                    Info<< "Additional directional refinement level"
+                        << " for all cells inside " << geomName << endl;
+                }
+                else if (modes_[shellI] == OUTSIDE)
+                {
+                    Info<< "Additional directional refinement level"
+                        << " for all cells outside " << geomName << endl;
+                }
+                else
+                {
+                    FatalIOErrorInFunction(shellsDict)
+                        << "Unsupported mode "
+                        << refineModeNames_[modes_[shellI]]
+                        << exit(FatalIOError);
+                }
+            }
+
+
 
             // Gap specification
             // ~~~~~~~~~~~~~~~~~
@@ -824,4 +850,74 @@ void Foam::shellSurfaces::findLevel
 }
 
 
+void Foam::shellSurfaces::findDirectionalLevel
+(
+    const pointField& pt,
+    const labelList& ptLevel,
+    const labelList& dirLevel,  // directional level
+    const direction dir,
+    labelList& shell
+) const
+{
+    shell.setSize(pt.size());
+    shell = -1;
+
+    List<volumeType> volType;
+
+    // Current back to original
+    DynamicList<label> candidateMap(pt.size());
+
+    forAll(shells_, shelli)
+    {
+        if (modes_[shelli] == INSIDE || modes_[shelli] == OUTSIDE)
+        {
+            const LevelAndDirList& shellLevels = dirLevels_[shelli];
+
+            // Collect the cells that are of the right original level
+            candidateMap.clear();
+            forAll(ptLevel, celli)
+            {
+                label level = ptLevel[celli];
+                forAll(shellLevels, leveli)
+                {
+                    label selectLevel = shellLevels[leveli].first();
+                    label addLevel = shellLevels[leveli].second()[dir];
+
+                    if
+                    (
+                        level == selectLevel
+                     && dirLevel[celli] < level+addLevel
+                    )
+                    {
+                        candidateMap.append(celli);
+                        break;
+                    }
+                }
+            }
+
+            pointField candidatePt(pt, candidateMap);
+            allGeometry_[shells_[shelli]].getVolumeType(candidatePt, volType);
+
+            forAll(candidateMap, i)
+            {
+                if
+                (
+                    (
+                        modes_[shelli] == INSIDE
+                     && volType[i] == volumeType::INSIDE
+                    )
+                 || (
+                        modes_[shelli] == OUTSIDE
+                     && volType[i] == volumeType::OUTSIDE
+                    )
+                )
+                {
+                    shell[candidateMap[i]] = shelli;
+                }
+            }
+        }
+    }
+}
+
+
 // ************************************************************************* //

src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.H

@@ -39,12 +39,17 @@ SourceFiles
 #include "searchableSurface.H"
 #include "Enum.H"
 #include "Tuple2.H"
+#include "labelVector.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
+typedef Tuple2<label,labelVector> LevelAndDir;
+typedef List<LevelAndDir> LevelAndDirList;
+typedef List<LevelAndDirList> LevelAndDirListList;
+
 class searchableSurfaces;
 
 /*---------------------------------------------------------------------------*\
@@ -85,6 +90,9 @@ private:
         //- Per shell per distance the refinement level
         labelListList levels_;
 
+        //- Per shell, per refinement level additional directional refinement
+        LevelAndDirListList dirLevels_;
+
 
         // Gap level refinement
 
@@ -171,6 +179,13 @@ public:
                 return shells_;
             }
 
+            //- Raw access to directional refinement
+            //  Per shell a list of (level + additional level)
+            const LevelAndDirListList& dirLevels() const
+            {
+                return dirLevels_;
+            }
+
 
         // Query
 
@@ -215,6 +230,16 @@ public:
                 const labelList& ptLevel,
                 labelList& shell
             ) const;
+
+            //- Find any shell (or -1) with higher wanted directional level
+            void findDirectionalLevel
+            (
+                const pointField& pt,
+                const labelList& ptLevel,   // omnidirectional level
+                const labelList& dirLevel,  // directional level
+                const direction dir,
+                labelList& shell
+            ) const;
 };
 
 

src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C

@@ -40,6 +40,7 @@ License
 #include "IOmanip.H"
 #include "labelVector.H"
 #include "profiling.H"
+#include "searchableSurfaces.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -1571,6 +1572,400 @@ 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
+(
+    const refinementParameters& refineParams,
+    const label maxIter
+)
+{
+    addProfiling(shell, "snappyHexMesh::refine::directionalShell");
+    const fvMesh& mesh = meshRefiner_.mesh();
+    const shellSurfaces& shells = meshRefiner_.dirShells();
+
+    labelList& cellLevel =
+        const_cast<labelIOList&>(meshRefiner_.meshCutter().cellLevel());
+
+
+    // Determine the minimum and maximum cell levels that are candidates for
+    // directional refinement
+    label overallMinLevel = labelMax;
+    label overallMaxLevel = labelMin;
+    {
+        const LevelAndDirListList& dirLevels = shells.dirLevels();
+        forAll(dirLevels, shelli)
+        {
+            const LevelAndDirList& dirLevel = dirLevels[shelli];
+            forAll(dirLevel, i)
+            {
+                overallMinLevel = min(dirLevel[i].first(), overallMinLevel);
+                overallMaxLevel = max(dirLevel[i].first(), overallMaxLevel);
+            }
+        }
+    }
+
+    if (overallMinLevel > overallMaxLevel)
+    {
+        return 0;
+    }
+
+    // Maintain directional refinement levels
+    List<labelVector> dirCellLevel(cellLevel.size());
+    forAll(cellLevel, celli)
+    {
+        label l = cellLevel[celli];
+        dirCellLevel[celli] = labelVector(l, l, l);
+    }
+
+    label iter;
+    for (iter = 0; iter < maxIter; iter++)
+    {
+        Info<< nl
+            << "Directional shell refinement iteration " << iter << nl
+            << "--------------------------------------" << nl
+            << endl;
+
+        label nAllRefine = 0;
+
+        for (direction dir = 0; dir < vector::nComponents; dir++)
+        {
+            // Select the cells that need to be refined in certain direction:
+            // 
+            // - cell inside/outside shell
+            // - 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)
+                {
+                    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;
+                    }
+                }
+            }
+
+            // Extend to keep 2:1 ratio
+            labelList cellsToRefine
+            (
+                consistentDirectionalRefinement
+                (
+                    dir,
+                    dirCellLevel,
+                    candidateCells,
+                    true
+                )
+            );
+
+            Info<< "Determined cells to refine in = "
+                << mesh.time().cpuTimeIncrement() << " s" << endl;
+
+            label nCellsToRefine = cellsToRefine.size();
+            reduce(nCellsToRefine, sumOp<label>());
+
+            Info<< "Selected for direction " << vector::componentNames[dir]
+                << " refinement : " << nCellsToRefine
+                << " cells (out of " << mesh.globalData().nTotalCells()
+                << ')' << endl;
+
+            nAllRefine += nCellsToRefine;
+
+            // Stop when no cells to refine or have done minimum necessary
+            // iterations and not enough cells to refine.
+            if (nCellsToRefine == 0)
+            {
+                Info<< "Not refining direction " << dir
+                    << " since too few cells selected."
+                    << nl << endl;
+                break;
+            }
+
+
+            if (debug)
+            {
+                const_cast<Time&>(mesh.time())++;
+            }
+
+            PackedBoolList isRefineCell(mesh.nCells());
+            isRefineCell.set(cellsToRefine);
+
+            autoPtr<mapPolyMesh> map
+            (
+                meshRefiner_.directionalRefine
+                (
+                    "directional refinement iteration " + name(iter),
+                    dir,
+                    cellsToRefine
+                )
+            );
+
+            meshRefinement::updateList
+            (
+                map().cellMap(),
+                labelVector(0, 0, 0),
+                dirCellLevel
+            );
+
+            forAll(map().cellMap(), celli)
+            {
+                if (isRefineCell[map().cellMap()[celli]])
+                {
+                    dirCellLevel[celli][dir]++;
+                }
+            }
+        }
+
+
+        if (nAllRefine == 0)
+        {
+            Info<< "Stopping refining since no cells selected."
+                << nl << endl;
+            break;
+        }
+
+        meshRefiner_.printMeshInfo
+        (
+            debug,
+            "After directional refinement iteration " + name(iter)
+        );
+
+        if (debug&meshRefinement::MESH)
+        {
+            Pout<< "Writing directional refinement iteration " + name(iter)
+                << " mesh to time " << meshRefiner_.timeName() << endl;
+            meshRefiner_.write
+            (
+                meshRefinement::debugType(debug),
+                meshRefinement::writeType
+                (
+                    meshRefinement::writeLevel()
+                  | meshRefinement::WRITEMESH
+                ),
+                mesh.time().path()/meshRefiner_.timeName()
+            );
+        }
+    }
+
+    // Adjust cellLevel from dirLevel?
+    // Is cellLevel the max of dirLevel? Or the min?
+
+    return iter;
+}
 
 
 void Foam::snappyRefineDriver::baffleAndSplitMesh
@@ -2065,6 +2460,15 @@ void Foam::snappyRefineDriver::doRefine
         10      // maxIter
     );
 
+    // Directional refinement
+//XXXXX
+    directionalShellRefine
+    (
+        refineParams,
+        100    // maxIter
+    );
+//XXXXX
+
     // Introduce baffles at surface intersections. Remove sections unreachable
     // from keepPoint.
     baffleAndSplitMesh

src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.H

@@ -36,6 +36,8 @@ SourceFiles
 
 #include "wordPairHashTable.H"
 #include "labelList.H"
+#include "PackedBoolList.H"
+#include "labelVector.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -135,13 +137,42 @@ class snappyRefineDriver
             const label nBufferLayers
         );
 
-        //- Remove all cells inside/outside shell
+        //- Refine all cells inside/outside shell
         label shellRefine
         (
             const refinementParameters& refineParams,
             const label maxIter
         );
 
+        // Directional refinement
+
+            //- Calculate direction of face (using geometry)
+            List<direction> faceDirection() const;
+
+            //- Check if level across (directed) face obeys 2:1
+            label faceConsistentRefinement
+            (
+                const PackedBoolList& isDirFace,
+                const List<labelVector>& dirCellLevel,
+                const direction dir,
+                const bool maxSet,
+                PackedBoolList& refineCell
+            ) const;
+            //- Make sure 2:1 refinement in given direction
+            labelList consistentDirectionalRefinement
+            (
+                const direction dir,
+                const List<labelVector>& dirCellLevel,
+                const labelUList& candidateCells,
+                const bool maxSet
+            ) const;
+            //- Refine (directional) all cells inside/outside shell
+            label directionalShellRefine
+            (
+                const refinementParameters& refineParams,
+                const label maxIter
+            );
+
         //- Add baffles and remove unreachable cells
         void baffleAndSplitMesh
         (