ENH: snappyHexMesh: layer coverage volfields

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

@@ -1253,6 +1253,9 @@ int main(int argc, char *argv[])
     // Snap parameters
     const snapParameters snapParams(snapDict);
 
+    // Layer addition parameters
+    const layerParameters layerParams(layerDict, mesh.boundaryMesh());
+
 
     if (wantRefine)
     {
@@ -1346,9 +1349,6 @@ int main(int argc, char *argv[])
             globalToSlavePatch
         );
 
-        // Layer addition parameters
-        layerParameters layerParams(layerDict, mesh.boundaryMesh());
-
         // Use the maxLocalCells from the refinement parameters
         bool preBalance = returnReduce
         (

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict

@@ -304,11 +304,13 @@ addLayersControls
     // size of the refined cell outside layer (true) or absolute sizes (false).
     relativeSizes true;
 
-    // Layer thickness specification. This can be specified in one of four ways
+    // Layer thickness specification. This can be specified in one of following
+    // ways:
     // - expansionRatio and finalLayerThickness (cell nearest internal mesh)
     // - expansionRatio and firstLayerThickness (cell on surface)
     // - overall thickness and firstLayerThickness
     // - overall thickness and finalLayerThickness
+    // - overall thickness and expansionRatio
 
         // Expansion factor for layer mesh
         expansionRatio 1.0;

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C

@@ -52,6 +52,7 @@ Description
 #include "fixedValuePointPatchFields.H"
 #include "calculatedPointPatchFields.H"
 #include "cyclicSlipPointPatchFields.H"
+#include "fixedValueFvPatchFields.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -110,6 +111,31 @@ void Foam::autoLayerDriver::dumpDisplacement
 }
 
 
+Foam::tmp<Foam::scalarField> Foam::autoLayerDriver::avgPointData
+(
+    const indirectPrimitivePatch& pp,
+    const scalarField& pointFld
+)
+{
+    tmp<scalarField> tfaceFld(new scalarField(pp.size(), 0.0));
+    scalarField& faceFld = tfaceFld();
+
+    forAll(pp.localFaces(), faceI)
+    {
+        const face& f = pp.localFaces()[faceI];
+        if (f.size())
+        {
+            forAll(f, fp)
+            {
+                faceFld[faceI] += pointFld[f[fp]];
+            }
+            faceFld[faceI] /= f.size();
+        }
+    }
+    return tfaceFld;
+}
+
+
 // Check that primitivePatch is not multiply connected. Collect non-manifold
 // points in pointSet.
 void Foam::autoLayerDriver::checkManifold
@@ -2391,19 +2417,21 @@ Foam::label Foam::autoLayerDriver::countExtrusion
 }
 
 
-// Collect layer faces and layer cells into bools for ease of handling
+// Collect layer faces and layer cells into mesh fields for ease of handling
 void Foam::autoLayerDriver::getLayerCellsFaces
 (
     const polyMesh& mesh,
     const addPatchCellLayer& addLayer,
-    boolList& flaggedCells,
+    const scalarField& oldRealThickness,
-    boolList& flaggedFaces
+
+    labelList& cellNLayers,
+    scalarField& faceRealThickness
 )
 {
-    flaggedCells.setSize(mesh.nCells());
+    cellNLayers.setSize(mesh.nCells());
-    flaggedCells = false;
+    cellNLayers = 0;
-    flaggedFaces.setSize(mesh.nFaces());
+    faceRealThickness.setSize(mesh.nFaces());
-    flaggedFaces = false;
+    faceRealThickness = 0;
 
     // Mark all faces in the layer
     const labelListList& layerFaces = addLayer.layerFaces();
@@ -2415,27 +2443,195 @@ void Foam::autoLayerDriver::getLayerCellsFaces
     {
         const labelList& added = addedCells[oldPatchFaceI];
 
-        forAll(added, i)
+        const labelList& layer = layerFaces[oldPatchFaceI];
+
+        if (layer.size())
         {
-            flaggedCells[added[i]] = true;
+            forAll(added, i)
+            {
+                cellNLayers[added[i]] = layer.size()-1;
+            }
         }
     }
 
     forAll(layerFaces, oldPatchFaceI)
     {
         const labelList& layer = layerFaces[oldPatchFaceI];
+        const scalar realThickness = oldRealThickness[oldPatchFaceI];
 
         if (layer.size())
         {
-            for (label i = 1; i < layer.size()-1; i++)
+            // Layer contains both original boundary face and new boundary
+            // face so is nLayers+1
+            forAll(layer, i)
             {
-                flaggedFaces[layer[i]] = true;
+                faceRealThickness[layer[i]] = realThickness;
             }
         }
     }
 }
 
 
+bool Foam::autoLayerDriver::writeLayerData
+(
+    const fvMesh& mesh,
+    const labelList& patchIDs,
+    const labelList& cellNLayers,
+    const scalarField& faceWantedThickness,
+    const scalarField& faceRealThickness
+) const
+{
+    bool allOk = true;
+
+    {
+        label nAdded = 0;
+        forAll(cellNLayers, cellI)
+        {
+            if (cellNLayers[cellI] > 0)
+            {
+                nAdded++;
+            }
+        }
+        cellSet addedCellSet(mesh, "addedCells", nAdded);
+        forAll(cellNLayers, cellI)
+        {
+            if (cellNLayers[cellI] > 0)
+            {
+                addedCellSet.insert(cellI);
+            }
+        }
+        addedCellSet.instance() = meshRefiner_.timeName();
+        Info<< "Writing "
+            << returnReduce(addedCellSet.size(), sumOp<label>())
+            << " added cells to cellSet "
+            << addedCellSet.name() << endl;
+        bool ok = addedCellSet.write();
+        allOk = allOk & ok;
+    }
+    {
+        label nAdded = 0;
+        forAll(faceRealThickness, faceI)
+        {
+            if (faceRealThickness[faceI] > 0)
+            {
+                nAdded++;
+            }
+        }
+
+        faceSet layerFacesSet(mesh, "layerFaces", nAdded);
+        forAll(faceRealThickness, faceI)
+        {
+            if (faceRealThickness[faceI] > 0)
+            {
+                layerFacesSet.insert(faceI);
+            }
+        }
+        layerFacesSet.instance() = meshRefiner_.timeName();
+        Info<< "Writing "
+            << returnReduce(layerFacesSet.size(), sumOp<label>())
+            << " faces inside added layer to faceSet "
+            << layerFacesSet.name() << endl;
+        bool ok = layerFacesSet.write();
+        allOk = allOk & ok;
+    }
+    {
+        volScalarField fld
+        (
+            IOobject
+            (
+                "nSurfaceLayers",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            mesh,
+            dimensionedScalar("zero", dimless, 0),
+            fixedValueFvPatchScalarField::typeName
+        );
+        const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+        forAll(patchIDs, i)
+        {
+            label patchI = patchIDs[i];
+            const polyPatch& pp = pbm[patchI];
+            const labelList& faceCells = pp.faceCells();
+            scalarField pfld(faceCells.size());
+            forAll(faceCells, i)
+            {
+                pfld[i] = cellNLayers[faceCells[i]];
+            }
+            fld.boundaryField()[patchI] == pfld;
+        }
+        Info<< "Writing volScalarField " << fld.name()
+            << " with actual number of layers" << endl;
+        bool ok = fld.write();
+        allOk = allOk & ok;
+    }
+    {
+        volScalarField fld
+        (
+            IOobject
+            (
+                "wantedThickness",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            mesh,
+            dimensionedScalar("zero", dimless, 0),
+            fixedValueFvPatchScalarField::typeName
+        );
+        const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+        forAll(patchIDs, i)
+        {
+            label patchI = patchIDs[i];
+            fld.boundaryField()[patchI] == pbm[patchI].patchSlice
+            (
+                faceWantedThickness
+            );
+        }
+        Info<< "Writing volScalarField " << fld.name()
+            << " with wanted thickness" << endl;
+        bool ok = fld.write();
+        allOk = allOk & ok;
+    }
+    {
+        volScalarField fld
+        (
+            IOobject
+            (
+                "thickness",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            mesh,
+            dimensionedScalar("zero", dimless, 0),
+            fixedValueFvPatchScalarField::typeName
+        );
+        const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+        forAll(patchIDs, i)
+        {
+            label patchI = patchIDs[i];
+            fld.boundaryField()[patchI] == pbm[patchI].patchSlice
+            (
+                faceRealThickness
+            );
+        }
+        Info<< "Writing volScalarField " << fld.name()
+            << " with layer thickness" << endl;
+        bool ok = fld.write();
+        allOk = allOk & ok;
+    }
+    return allOk;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::autoLayerDriver::autoLayerDriver
@@ -2618,7 +2814,7 @@ void Foam::autoLayerDriver::addLayers
     labelList patchNLayers(pp().nPoints(), 0);
 
     // Ideal number of cells added
-    label nIdealAddedCells = 0;
+    label nIdealTotAddedCells = 0;
 
     // Whether to add edge for all pp.localPoints.
     List<extrudeMode> extrudeStatus(pp().nPoints(), EXTRUDE);
@@ -2637,7 +2833,7 @@ void Foam::autoLayerDriver::addLayers
             patchDisp,
             patchNLayers,
             extrudeStatus,
-            nIdealAddedCells
+            nIdealTotAddedCells
         );
 
         // Precalculate mesh edge labels for patch edges
@@ -2835,11 +3031,20 @@ void Foam::autoLayerDriver::addLayers
     // Saved old points
     pointField oldPoints(mesh.points());
 
-    // Last set of topology changes. (changing mesh clears out polyTopoChange)
+    // Current set of topology changes. (changing mesh clears out
+    // polyTopoChange)
     polyTopoChange savedMeshMod(mesh.boundaryMesh().size());
+    // Per cell 0 or number of layers in the cell column it is part of
+    labelList cellNLayers;
+    // Per face actual overall layer thickness
+    scalarField faceRealThickness;
+    // Per face wanted overall layer thickness
+    scalarField faceWantedThickness(mesh.nFaces(), 0.0);
+    {
+        UIndirectList<scalar>(faceWantedThickness, pp().addressing()) =
+            avgPointData(pp, thickness);
+    }
 
-    boolList flaggedCells;
-    boolList flaggedFaces;
 
     for (label iteration = 0; iteration < layerParams.nLayerIter(); iteration++)
     {
@@ -3088,35 +3293,52 @@ void Foam::autoLayerDriver::addLayers
         (
             newMesh,
             addLayer,
-            flaggedCells,
+            avgPointData(pp, mag(patchDisp))(), // current thickness
-            flaggedFaces
+
+            cellNLayers,
+            faceRealThickness
         );
 
 
+        // Count number of added cells
+        label nAddedCells = 0;
+        forAll(cellNLayers, cellI)
+        {
+            if (cellNLayers[cellI] > 0)
+            {
+                nAddedCells++;
+            }
+        }
+
+
         if (debug&meshRefinement::MESH)
         {
             Info<< "Writing layer mesh to time " << meshRefiner_.timeName()
                 << endl;
             newMesh.write();
-            cellSet addedCellSet
+
-            (
+            cellSet addedCellSet(newMesh," addedCells", nAddedCells);
-                newMesh,
+            forAll(cellNLayers, cellI)
-                "addedCells",
+            {
-                findIndices(flaggedCells, true)
+                if (cellNLayers[cellI] > 0)
-            );
+                {
+                    addedCellSet.insert(cellI);
+                }
+            }
             addedCellSet.instance() = meshRefiner_.timeName();
             Info<< "Writing "
                 << returnReduce(addedCellSet.size(), sumOp<label>())
-                << " added cells to cellSet "
+                << " added cells to cellSet " << addedCellSet.name() << endl;
-                << addedCellSet.name() << endl;
             addedCellSet.write();
 
-            faceSet layerFacesSet
+            faceSet layerFacesSet(newMesh, "layerFaces", newMesh.nFaces()/100);
-            (
+            forAll(faceRealThickness, faceI)
-                newMesh,
+            {
-                "layerFaces",
+                if (faceRealThickness[faceI] > 0)
-                findIndices(flaggedCells, true)
+                {
-            );
+                    layerFacesSet.insert(faceI);
+                }
+            }
             layerFacesSet.instance() = meshRefiner_.timeName();
             Info<< "Writing "
                 << returnReduce(layerFacesSet.size(), sumOp<label>())
@@ -3140,26 +3362,17 @@ void Foam::autoLayerDriver::addLayers
             extrudeStatus
         );
 
-        label nExtruded = countExtrusion(pp, extrudeStatus);
+        label nTotExtruded = countExtrusion(pp, extrudeStatus);
         label nTotFaces = returnReduce(pp().size(), sumOp<label>());
-        label nAddedCells = 0;
+        label nTotAddedCells = returnReduce(nAddedCells, sumOp<label>());
-        {
+
-            forAll(flaggedCells, cellI)
+        Info<< "Extruding " << nTotExtruded
-            {
-                if (flaggedCells[cellI])
-                {
-                    nAddedCells++;
-                }
-            }
-            reduce(nAddedCells, sumOp<label>());
-        }
-        Info<< "Extruding " << nExtruded
             << " out of " << nTotFaces
-            << " faces (" << 100.0*nExtruded/nTotFaces << "%)."
+            << " faces (" << 100.0*nTotExtruded/nTotFaces << "%)."
             << " Removed extrusion at " << nTotChanged << " faces."
             << endl
-            << "Added " << nAddedCells << " out of " << nIdealAddedCells
+            << "Added " << nTotAddedCells << " out of " << nIdealTotAddedCells
-            << " cells (" << 100.0*nAddedCells/nIdealAddedCells << "%)."
+            << " cells (" << 100.0*nTotAddedCells/nIdealTotAddedCells << "%)."
             << endl;
 
         if (nTotChanged == 0)
@@ -3217,6 +3430,8 @@ void Foam::autoLayerDriver::addLayers
 
     meshRefiner_.updateMesh(map, labelList(0));
 
+    // Update numbering of faceWantedThickness
+    meshRefinement::updateList(map().faceMap(), 0.0, faceWantedThickness);
 
     // Update numbering on baffles
     forAll(baffles, i)
@@ -3237,8 +3452,9 @@ void Foam::autoLayerDriver::addLayers
 
         autoPtr<mapPolyMesh> map = meshRefiner_.mergeBaffles(baffles);
 
-        inplaceReorder(map().reverseCellMap(), flaggedCells);
+        inplaceReorder(map().reverseCellMap(), cellNLayers);
-        inplaceReorder(map().reverseFaceMap(), flaggedFaces);
+        inplaceReorder(map().reverseFaceMap(), faceWantedThickness);
+        inplaceReorder(map().reverseFaceMap(), faceRealThickness);
 
         Info<< "Converted baffles in = "
             << meshRefiner_.mesh().time().cpuTimeIncrement()
@@ -3272,29 +3488,23 @@ void Foam::autoLayerDriver::addLayers
         );
 
         // Re-distribute flag of layer faces and cells
-        map().distributeCellData(flaggedCells);
+        map().distributeCellData(cellNLayers);
-        map().distributeFaceData(flaggedFaces);
+        map().distributeFaceData(faceWantedThickness);
+        map().distributeFaceData(faceRealThickness);
     }
 
 
-    // Write mesh
+    // Write mesh data
-    // ~~~~~~~~~~
+    // ~~~~~~~~~~~~~~~
 
-    cellSet addedCellSet(mesh, "addedCells", findIndices(flaggedCells, true));
+    writeLayerData
-    addedCellSet.instance() = meshRefiner_.timeName();
+    (
-    Info<< "Writing "
+        mesh,
-        << returnReduce(addedCellSet.size(), sumOp<label>())
+        patchIDs,
-        << " added cells to cellSet "
+        cellNLayers,
-        << addedCellSet.name() << endl;
+        faceWantedThickness,
-    addedCellSet.write();
+        faceRealThickness
-
+    );
-    faceSet layerFacesSet(mesh, "layerFaces", findIndices(flaggedFaces, true));
-    layerFacesSet.instance() = meshRefiner_.timeName();
-    Info<< "Writing "
-        << returnReduce(layerFacesSet.size(), sumOp<label>())
-        << " faces inside added layer to faceSet "
-        << layerFacesSet.name() << endl;
-    layerFacesSet.write();
 }
 
 

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H

@@ -120,6 +120,13 @@ class autoLayerDriver
                 const List<extrudeMode>&
             );
 
+            //- Average point wise data to face wise
+            static tmp<scalarField> avgPointData
+            (
+                const indirectPrimitivePatch&,
+                const scalarField& pointFld
+            );
+
             //- Check that primitivePatch is not multiply connected.
             //  Collect non-manifold points in pointSet.
             static void checkManifold
@@ -367,10 +374,23 @@ class autoLayerDriver
                 (
                     const polyMesh&,
                     const addPatchCellLayer&,
-                    boolList&,
+                    const scalarField& oldRealThickness,
-                    boolList&
+
+                    labelList& cellStatus,
+                    scalarField& faceRealThickness
                 );
 
+                //- Write cellSet,faceSet for layers
+                bool writeLayerData
+                (
+                    const fvMesh& mesh,
+                    const labelList& patchIDs,
+                    const labelList& cellNLayers,
+                    const scalarField& faceWantedThickness,
+                    const scalarField& faceRealThickness
+                ) const;
+
+
             // Mesh shrinking (to create space for layers)
 
                 //- Average field (over all subset of mesh points) by

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C

@@ -235,6 +235,12 @@ Foam::layerParameters::layerParameters
         Info<< "Layer thickness specified as final layer and expansion ratio."
             << endl;
     }
+    else if (haveTotal && haveExp)
+    {
+        layerSpec_ = TOTAL_AND_EXPANSION;
+        Info<< "Layer thickness specified as overall thickness"
+            << " and expansion ratio." << endl;
+    }
 
 
     if (layerSpec_ == ILLEGAL || nSpec != 2)
@@ -253,7 +259,9 @@ Foam::layerParameters::layerParameters
             << "    final layer thickness ('finalLayerThickness')"
             << " and expansion ratio ('expansionRatio') or" << nl
             << "    final layer thickness ('finalLayerThickness')"
-            << " and overall thickness ('thickness')"
+            << " and overall thickness ('thickness') or" << nl
+            << "    overall thickness ('thickness')"
+            << " and expansion ratio ('expansionRatio'"
             << exit(FatalIOError);
     }
 
@@ -354,6 +362,19 @@ Foam::layerParameters::layerParameters
                             );
                         break;
 
+                        case TOTAL_AND_EXPANSION:
+                            layerDict.readIfPresent
+                            (
+                                "thickness",
+                                thickness_[patchI]
+                            );
+                            layerDict.readIfPresent
+                            (
+                                "expansionRatio",
+                                expansionRatio_[patchI]
+                            );
+                        break;
+
                         default:
                             FatalIOErrorIn
                             (
@@ -390,6 +411,7 @@ Foam::scalar Foam::layerParameters::layerThickness
     {
         case FIRST_AND_TOTAL:
         case FINAL_AND_TOTAL:
+        case TOTAL_AND_EXPANSION:
         {
             return totalThickness;
         }
@@ -450,6 +472,7 @@ Foam::scalar Foam::layerParameters::layerExpansionRatio
     {
         case FIRST_AND_EXPANSION:
         case FINAL_AND_EXPANSION:
+        case TOTAL_AND_EXPANSION:
         {
             return expansionRatio;
         }
@@ -524,6 +547,18 @@ Foam::scalar Foam::layerParameters::firstLayerThickness
         }
         break;
 
+        case TOTAL_AND_EXPANSION:
+        {
+            scalar r = finalLayerThicknessRatio
+            (
+                nLayers,
+                expansionRatio
+            );
+            scalar finalThickness = r*totalThickness;
+            return finalThickness/pow(expansionRatio, nLayers-1);
+        }
+        break;
+
         default:
         {
             FatalErrorIn("layerParameters::layerThickness(..)")

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H

@@ -64,13 +64,15 @@ public:
         //  - first and expansion ratio specified
         //  - final and total thickness specified
         //  - final and expansion ratio specified
+        //  - total thickness and expansion ratio specified
         enum layerSpecification
         {
             ILLEGAL,
             FIRST_AND_TOTAL,
             FIRST_AND_EXPANSION,
             FINAL_AND_TOTAL,
-            FINAL_AND_EXPANSION
+            FINAL_AND_EXPANSION,
+            TOTAL_AND_EXPANSION
         };