ENH: snappyHexMesh: directional smoothing support

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict

@@ -299,6 +299,12 @@ castellatedMeshControls
         //    //   with e.g. dynamic refinement/unrefinement.
         //    // - cellLevel will include any directional refinement
         //    //   (i.e. it will be the maximum of all three directions)
+        //
+        //    // Optional directional expansion-ratio smoothing (after all
+        //    // refinement). This will try to smooth out edge/cell size jumps
+        //    // Specify smoothing direction and number of iterations
+        //    nSmoothExpansion    100;
+        //    smoothDirection     (1 0 0);
         //}
     }
 

src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C

@@ -585,6 +585,8 @@ Foam::shellSurfaces::shellSurfaces
     distances_.setSize(shellI);
     levels_.setSize(shellI);
     dirLevels_.setSize(shellI);
+    smoothDirection_.setSize(shellI);
+    nSmoothExpansion_.setSize(shellI);
 
     extendedGapLevel_.setSize(shellI);
     extendedGapMode_.setSize(shellI);
@@ -661,6 +663,13 @@ Foam::shellSurfaces::shellSurfaces
                 }
             }
 
+            // Directional smoothing
+            // ~~~~~~~~~~~~~~~~~~~~~
+
+            smoothDirection_[shellI] =
+                dict.lookupOrDefault("smoothDirection", vector::zero);
+            nSmoothExpansion_[shellI] =
+                dict.lookupOrDefault("nSmoothExpansion", 0);
 
 
             // Gap specification
@@ -803,6 +812,18 @@ Foam::labelPairList Foam::shellSurfaces::directionalSelectLevel() const
 }
 
 
+const Foam::labelList& Foam::shellSurfaces::nSmoothExpansion() const
+{
+    return nSmoothExpansion_;
+}
+
+
+const Foam::vectorField& Foam::shellSurfaces::smoothDirection() const
+{
+    return smoothDirection_;
+}
+
+
 void Foam::shellSurfaces::findHigherLevel
 (
     const pointField& pt,

src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.H

@@ -86,8 +86,17 @@ private:
         //- Per shell per distance the refinement level
         labelListList levels_;
 
-        //- Per shell any additional directional refinement
+
-        List<Tuple2<labelPair,labelVector>> dirLevels_;
+        // Directional
+
+            //- Per shell any additional directional refinement
+            List<Tuple2<labelPair,labelVector>> dirLevels_;
+
+            //- Per shell the smoothing direction
+            vectorField smoothDirection_;
+
+            //- Per shell the directional smoothing iterations
+            labelList nSmoothExpansion_;
 
 
         // Gap level refinement
@@ -231,6 +240,12 @@ public:
                 const direction dir,
                 labelList& shell
             ) const;
+
+            //- Per shell the smoothing direction
+            const vectorField& smoothDirection() const;
+
+            //- Per shell the directional smoothing iterations
+            const labelList& nSmoothExpansion() const;
 };
 
 

src/mesh/snappyHexMesh/snappyHexMeshDriver/PointIntegrateData/PointIntegrateData.H

@@ -0,0 +1,258 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::PointIntegrateData
+
+Description
+    Integrate along selected edges using PointEdgeWave.
+
+SourceFiles
+    PointIntegrateDataI.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PointIntegrateData_H
+#define PointIntegrateData_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class Istream;
+class Ostream;
+template<class DataType>
+class PointIntegrateData;
+
+// Forward declaration of friend functions and operators
+template<class DataType>
+Ostream& operator<<(Ostream&, const PointIntegrateData<DataType>&);
+template<class DataType>
+Istream& operator>>(Istream&, PointIntegrateData<DataType>&);
+
+/*---------------------------------------------------------------------------*\
+                     Class PointIntegrateData declaration
+\*---------------------------------------------------------------------------*/
+
+template<class DataType>
+class PointIntegrateData
+{
+private:
+
+    // Private data
+
+        //- Valid flag
+        bool valid_;
+
+        //- Integrated data
+        DataType data_;
+
+
+public:
+
+    //- Class used to pass extra data
+    class trackingData
+    {
+    public:
+
+        UList<DataType>& edgeData_;
+
+        trackingData(UList<DataType>& edgeData)
+        :
+            edgeData_(edgeData)
+        {}
+    };
+
+
+
+    // Constructors
+
+        //- Construct null
+        inline PointIntegrateData();
+
+        //- Construct from data
+        inline PointIntegrateData(const DataType& data);
+
+
+    // Member Functions
+
+        // Access
+
+            //- Const access the data
+            inline const DataType& data() const
+            {
+                return data_;
+            };
+
+
+        // Needed by PointEdgeWave
+
+            //- Check whether original (invalid) value.
+            template<class TrackingData>
+            inline bool valid(TrackingData& td) const;
+
+            //- Check for identical geometrical data. Used for cyclics checking.
+            template<class TrackingData>
+            inline bool sameGeometry
+            (
+                const PointIntegrateData<DataType>&,
+                const scalar tol,
+                TrackingData& td
+            ) const;
+
+            //- Convert origin to relative vector to leaving point
+            //  (= point coordinate)
+            template<class TrackingData>
+            inline void leaveDomain
+            (
+                const polyPatch& patch,
+                const label patchPointi,
+                const point& pos,
+                TrackingData& td
+            );
+
+            //- Convert relative origin to absolute by adding entering point
+            template<class TrackingData>
+            inline void enterDomain
+            (
+                const polyPatch& patch,
+                const label patchPointi,
+                const point& pos,
+                TrackingData& td
+            );
+
+            //- Apply rotation matrix to the data
+            template<class TrackingData>
+            inline void transform
+            (
+                const tensor& rotTensor,
+                TrackingData& td
+            );
+
+            //- Influence of edge on point
+            template<class TrackingData>
+            inline bool updatePoint
+            (
+                const polyMesh& mesh,
+                const label pointI,
+                const label edgeI,
+                const PointIntegrateData<DataType>& edgeInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Influence of different value on same point.
+            //  Merge new and old info.
+            template<class TrackingData>
+            inline bool updatePoint
+            (
+                const polyMesh& mesh,
+                const label pointI,
+                const PointIntegrateData<DataType>& newPointInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Influence of different value on same point.
+            //  No information about current position whatsoever.
+            template<class TrackingData>
+            inline bool updatePoint
+            (
+                const PointIntegrateData<DataType>& newPointInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Influence of point on edge.
+            template<class TrackingData>
+            inline bool updateEdge
+            (
+                const polyMesh& mesh,
+                const label edgeI,
+                const label pointI,
+                const PointIntegrateData<DataType>& pointInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Same (like operator==)
+            template<class TrackingData>
+            inline bool equal
+            (
+                const PointIntegrateData<DataType>&,
+                TrackingData& td
+            ) const;
+
+
+    // Member Operators
+
+        // Needed for List IO
+        inline bool operator==(const PointIntegrateData<DataType>&) const;
+        inline bool operator!=(const PointIntegrateData<DataType>&) const;
+
+
+    // IOstream Operators
+
+        friend Ostream& operator<< <DataType>
+        (
+            Ostream&,
+            const PointIntegrateData<DataType>&
+        );
+        friend Istream& operator>> <DataType>
+        (
+            Istream&,
+            PointIntegrateData<DataType>&
+        );
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//- Data associated with PointIntegrateData types is contiguous
+
+template<>
+inline bool contiguous<PointIntegrateData<scalar>>()
+{
+    return true;
+}
+
+template<>
+inline bool contiguous<PointIntegrateData<vector>>()
+{
+    return true;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "PointIntegrateDataI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/mesh/snappyHexMesh/snappyHexMeshDriver/PointIntegrateData/PointIntegrateDataI.H

@@ -0,0 +1,304 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "polyMesh.H"
+#include "transform.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class DataType>
+inline Foam::PointIntegrateData<DataType>::PointIntegrateData()
+:
+    valid_(false)
+{}
+
+
+template<class DataType>
+inline Foam::PointIntegrateData<DataType>::PointIntegrateData
+(
+    const DataType& data
+)
+:
+    valid_(true),
+    data_(data)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class DataType>
+template<class TrackingData>
+inline bool Foam::PointIntegrateData<DataType>::valid(TrackingData& td) const
+{
+    return valid_;
+}
+
+
+template<class DataType>
+template<class TrackingData>
+inline bool Foam::PointIntegrateData<DataType>::sameGeometry
+(
+    const PointIntegrateData<DataType>&,
+    const scalar tol,
+    TrackingData& td
+) const
+{
+    return true;
+}
+
+
+template<class DataType>
+template<class TrackingData>
+inline void Foam::PointIntegrateData<DataType>::leaveDomain
+(
+    const polyPatch& patch,
+    const label patchPointi,
+    const point& pos,
+    TrackingData& td
+)
+{}
+
+
+template<class DataType>
+template<class TrackingData>
+inline void Foam::PointIntegrateData<DataType>::enterDomain
+(
+    const polyPatch& patch,
+    const label patchPointi,
+    const point& pos,
+    TrackingData& td
+)
+{}
+
+
+template<class DataType>
+template<class TrackingData>
+inline void Foam::PointIntegrateData<DataType>::transform
+(
+    const tensor& rotTensor,
+    TrackingData& td
+)
+{
+    this->data_ = Foam::transform(rotTensor, this->data_);
+}
+
+
+template<class DataType>
+template<class TrackingData>
+inline bool Foam::PointIntegrateData<DataType>::updatePoint
+(
+    const polyMesh& mesh,
+    const label pointI,
+    const label edgeI,
+    const PointIntegrateData<DataType>& edgeInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    // Update point from an edge
+    if (!valid_)
+    {
+        if (!edgeInfo.valid_)
+        {
+            FatalErrorInFunction<< "edgeInfo:" << edgeInfo << exit(FatalError);
+        }
+        this->operator=(edgeInfo);
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template<class DataType>
+template<class TrackingData>
+inline bool Foam::PointIntegrateData<DataType>::updatePoint
+(
+    const polyMesh& mesh,
+    const label pointI,
+    const PointIntegrateData<DataType>& newPointInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    // Update point from coupled point
+    if (!valid_)
+    {
+        if (!newPointInfo.valid_)
+        {
+            FatalErrorInFunction<< "newPointInfo:" << newPointInfo
+                << exit(FatalError);
+        }
+        this->operator=(newPointInfo);
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template<class DataType>
+template<class TrackingData>
+inline bool Foam::PointIntegrateData<DataType>::updatePoint
+(
+    const PointIntegrateData<DataType>& newPointInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    if (!valid_)
+    {
+        if (!newPointInfo.valid_)
+        {
+            FatalErrorInFunction<< "newPointInfo:" << newPointInfo
+                << exit(FatalError);
+        }
+        this->operator=(newPointInfo);
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template<class DataType>
+template<class TrackingData>
+inline bool Foam::PointIntegrateData<DataType>::updateEdge
+(
+    const polyMesh& mesh,
+    const label edgeI,
+    const label pointI,
+    const PointIntegrateData<DataType>& pointInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    if (!valid_)
+    {
+        if (!pointInfo.valid(td))
+        {
+            FatalErrorInFunction<< "problem: invalid point:"
+                << mesh.points()[pointI]
+                << " data:" << pointInfo
+                << exit(FatalError);
+        }
+        this->data_ = pointInfo.data_ + td.edgeData_[edgeI];
+        this->valid_ = true;
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+//- Same (like operator==)
+template<class DataType>
+template<class TrackingData>
+inline bool Foam::PointIntegrateData<DataType>::equal
+(
+    const PointIntegrateData<DataType>& pi,
+    TrackingData& td
+) const
+{
+    if (!valid_)
+    {
+        return false;
+    }
+    else if (!pi.valid_)
+    {
+        FatalErrorInFunction << "pi:" << pi
+            << exit(FatalError);
+        return false;
+    }
+    else
+    {
+        return this->data_ == pi.data_;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+template<class DataType>
+inline bool Foam::PointIntegrateData<DataType>::operator==
+(
+    const Foam::PointIntegrateData<DataType>& rhs
+)
+const
+{
+    return this->data_ == rhs.data_;
+}
+
+
+template<class DataType>
+inline bool Foam::PointIntegrateData<DataType>::operator!=
+(
+    const Foam::PointIntegrateData<DataType>& rhs
+)
+const
+{
+    return !(*this == rhs);
+}
+
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+template<class DataType>
+inline Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const PointIntegrateData<DataType>& pd
+)
+{
+    if (os.format() == IOstream::ASCII)
+    {
+        return os << pd.valid_ << token::SPACE << pd.data();
+    }
+    else
+    {
+        return os << pd.valid_ << pd.data();
+    }
+}
+
+
+template<class DataType>
+inline Foam::Istream& Foam::operator>>
+(
+    Istream& is,
+    PointIntegrateData<DataType>& pd
+)
+{
+    return is >> pd.valid_ >> pd.data_;
+}
+
+
+// ************************************************************************* //

src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C

@@ -668,7 +668,7 @@ void Foam::snappyLayerDriver::handleFeatureAngle
 (
     const indirectPrimitivePatch& pp,
     const labelList& meshEdges,
-    const scalar minCos,
+    const scalar minAngle,
     pointField& patchDisp,
     labelList& patchNLayers,
     List<extrudeMode>& extrudeStatus
@@ -676,7 +676,10 @@ void Foam::snappyLayerDriver::handleFeatureAngle
 {
     const fvMesh& mesh = meshRefiner_.mesh();
 
-    Info<< nl << "Handling feature edges ..." << endl;
+    const scalar minCos = Foam::cos(degToRad(minAngle));
+
+    Info<< nl << "Handling feature edges (angle < " << minAngle
+        << ") ..." << endl;
 
     if (minCos < 1-SMALL)
     {
@@ -3647,7 +3650,7 @@ void Foam::snappyLayerDriver::addLayers
         (
             pp,
             meshEdges,
-            degToRad(layerParams.featureAngle()),
+            layerParams.featureAngle(),
 
             patchDisp,
             patchNLayers,

src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.H

@@ -218,7 +218,7 @@ private:
                 (
                     const indirectPrimitivePatch& pp,
                     const labelList& meshEdges,
-                    const scalar minCos,
+                    const scalar minAngle,
                     pointField& patchDisp,
                     labelList& patchNLayers,
                     List<extrudeMode>& extrudeStatus

src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015-2017 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -41,6 +41,9 @@ License
 #include "labelVector.H"
 #include "profiling.H"
 #include "searchableSurfaces.H"
+#include "PointIntegrateData.H"
+#include "PointEdgeWave.H"
+#include "OBJstream.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -1627,7 +1630,6 @@ Foam::label Foam::snappyRefineDriver::directionalShellRefine
         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
@@ -1786,6 +1788,576 @@ Foam::label Foam::snappyRefineDriver::directionalShellRefine
 }
 
 
+void Foam::snappyRefineDriver::getVirtualEdgeLength
+(
+    const vector& userDirection,
+    const labelList& pointLabels,
+    const PackedBoolList& isXFace,
+    scalarField& maxUserSize
+) const
+{
+    const fvMesh& mesh = meshRefiner_.mesh();
+    const labelList& own = mesh.faceOwner();
+    const labelList& nei = mesh.faceNeighbour();
+    const vectorField& fA = mesh.faceAreas();
+    const pointField& faceCentres = mesh.faceCentres();
+    const pointField& cellCentres = mesh.cellCentres();
+
+    scalarField sumWeights(mesh.nCells(), 0.0);
+    scalarField sumArea(mesh.nCells(), 0.0);
+
+    forAll(own, facei)
+    {
+        const scalar fc = (faceCentres[facei]&userDirection);
+        const scalar cc = (cellCentres[own[facei]]&userDirection);
+        // Calculate 'pyramid volume'
+        scalar pyr3Vol = (fA[facei]&userDirection)*(fc-cc);
+        sumWeights[own[facei]] += pyr3Vol;
+        sumArea[own[facei]] += mag(fA[facei]&userDirection);
+    }
+
+    forAll(nei, facei)
+    {
+        const scalar fc = (faceCentres[facei]&userDirection);
+        const scalar cc = (cellCentres[nei[facei]]&userDirection);
+        // Calculate 'pyramid volume'
+        scalar pyr3Vol = (fA[facei]&userDirection)*(cc-fc);
+        sumWeights[nei[facei]] += pyr3Vol;
+        sumArea[nei[facei]] += mag(fA[facei]&userDirection);
+    }
+
+    const scalarField cellSize(2*sumWeights/sumArea);
+
+    maxUserSize.setSize(mesh.nPoints());
+    maxUserSize = -GREAT;
+
+    forAll(pointLabels, i)
+    {
+        label pointi = pointLabels[i];
+        const labelList& pFaces = mesh.pointFaces()[pointi];
+        scalar& maxSz = maxUserSize[pointi];
+        forAll(pFaces, pFacei)
+        {
+            label facei = pFaces[pFacei];
+            if (isXFace[facei])
+            {
+                const point& fc = faceCentres[facei];
+                if (((cellCentres[own[facei]]-fc)&userDirection) > 0.0)
+                {
+                    maxSz = max(maxSz, cellSize[own[facei]]);
+                }
+                if (mesh.isInternalFace(facei))
+                {
+                    if (((cellCentres[nei[facei]]-fc)&userDirection) > 0.0)
+                    {
+                        maxSz = max(maxSz, cellSize[nei[facei]]);
+                    }
+                }
+            }
+        }
+    }
+    syncTools::syncPointList
+    (
+        mesh,
+        maxUserSize,
+        maxEqOp<scalar>(),
+        -GREAT
+    );
+}
+
+
+Foam::label Foam::snappyRefineDriver::directionalSmooth
+(
+    const refinementParameters& refineParams
+)
+{
+    addProfiling(split, "snappyHexMesh::refine::smooth");
+    Info<< nl
+        << "Directional expansion ratio smoothing" << nl
+        << "-------------------------------------" << nl
+        << endl;
+
+    fvMesh& mesh = meshRefiner_.mesh();
+    const pointField& points = mesh.points();
+    const searchableSurfaces& geometry = meshRefiner_.surfaces().geometry();
+    const shellSurfaces& shells = meshRefiner_.shells();
+
+    label iter = 0;
+
+    forAll(shells.nSmoothExpansion(), shellI)
+    {
+        if (shells.nSmoothExpansion()[shellI] > 0)
+        {
+            label surfi = shells.shells()[shellI];
+            const vector& userDirection = shells.smoothDirection()[shellI];
+
+
+            // Extract inside points
+            labelList pointLabels;
+            PackedBoolList isInsidePoint(mesh.nPoints());
+            {
+                // Get inside points
+                List<volumeType> volType;
+                geometry[surfi].getVolumeType(points, volType);
+
+                forAll(volType, pointi)
+                {
+                    if (volType[pointi] == volumeType::INSIDE)
+                    {
+                        isInsidePoint.set(pointi);
+                    }
+                }
+                pointLabels = isInsidePoint.used();
+            }
+
+
+            const edgeList& edges = mesh.edges();
+            const labelListList& pointEdges = mesh.pointEdges();
+
+            // Mark all directed edges
+            PackedBoolList isXEdge(edges.size());
+            {
+                forAll(edges, edgei)
+                {
+                    const edge& e = edges[edgei];
+                    vector eVec(e.vec(points));
+                    eVec /= mag(eVec);
+                    if (mag(eVec&userDirection) > 0.9)
+                    {
+                        isXEdge.set(edgei);
+                    }
+                }
+            }
+
+
+            // Mark all directed faces (with all points inside the region)
+            PackedBoolList isXFace(mesh.nFaces());
+            {
+                const vectorField& faceAreas = mesh.faceAreas();
+                const faceList& faces = mesh.faces();
+
+                forAll(faces, facei)
+                {
+                    const face& f = faces[facei];
+
+                    bool allInRegion = true;
+                    for (label pointi : f)
+                    {
+                        if (!isInsidePoint[pointi])
+                        {
+                            allInRegion = false;
+                            break;
+                        }
+                    }
+
+                    if (allInRegion)
+                    {
+                        vector n(faceAreas[facei]);
+                        n /= mag(n);
+                        if (mag(n&userDirection) > 0.9)
+                        {
+                            isXFace[facei] = true;
+                        }
+                    }
+                }
+            }
+
+
+            // Get left and right local edge
+            labelList leftEdge(points.size(), -1);
+            labelList rightEdge(points.size(), -1);
+
+            forAll(pointLabels, i)
+            {
+                label pointi = pointLabels[i];
+                const labelList& pEdges = pointEdges[pointi];
+                forAll(pEdges, pEdgei)
+                {
+                    label edgei = pEdges[pEdgei];
+                    const edge& e = edges[edgei];
+                    label otherPointi = e.otherVertex(pointi);
+                    vector eVec(points[otherPointi]-points[pointi]);
+                    eVec /= mag(eVec);
+                    scalar cosAngle = (eVec&userDirection);
+
+                    if (cosAngle < -0.9)
+                    {
+                        leftEdge[pointi] = edgei;
+                    }
+                    else if (cosAngle > 0.9)
+                    {
+                        rightEdge[pointi] = edgei;
+                    }
+                }
+            }
+
+
+            pointField newPoints(points);
+
+            Info<< "Smoothing " << geometry[surfi].name() << ':' << nl
+                << "    Direction                   : " << userDirection << nl
+                << "    Number of points            : "
+                << returnReduce(pointLabels.size(), sumOp<label>())
+                << " (out of " << mesh.globalData().nTotalPoints()
+                << ")" << nl
+                << "    Number of directed edges    : "
+                << returnReduce(isXEdge.count(), sumOp<label>())
+                << " (out of " << returnReduce(edges.size(), sumOp<label>())
+                << ")" << nl
+                << "    Number of directed faces    : "
+                << returnReduce(isXFace.count(), sumOp<label>())
+                << " (out of " << mesh.globalData().nTotalFaces()
+                << ")" << nl
+                << "    Iterations                  : "
+                << shells.nSmoothExpansion()[shellI] << nl
+                << endl;
+
+            for (iter = 0; iter < shells.nSmoothExpansion()[shellI]; iter++)
+            {
+                if (debug)
+                {
+                    const_cast<Time&>(mesh.time())++;
+                }
+
+                // Determine size of connected edge
+                scalarField leftEdgeMag(points.size(), -GREAT);
+                scalarField rightEdgeMag(points.size(), -GREAT);
+                {
+                    forAll(edges, edgei)
+                    {
+                        if (isXEdge[edgei])
+                        {
+                            const edge& e = edges[edgei];
+                            forAll(e, fp)
+                            {
+                                label pointi = e[fp];
+                                if (leftEdge[pointi] == edgei)
+                                {
+                                    leftEdgeMag[pointi] = e.mag(points);
+                                }
+                                else if (rightEdge[pointi] == edgei)
+                                {
+                                    rightEdgeMag[pointi] = e.mag(points);
+                                }
+                            }
+                        }
+                    }
+                    syncTools::syncPointList
+                    (
+                        mesh,
+                        leftEdgeMag,
+                        maxEqOp<scalar>(),
+                        -GREAT
+                    );
+                    syncTools::syncPointList
+                    (
+                        mesh,
+                        rightEdgeMag,
+                        maxEqOp<scalar>(),
+                        -GREAT
+                    );
+                }
+
+
+                // Get local characteristic size (in provided direction) on all
+                // points
+                scalarField maxUserSize;
+                getVirtualEdgeLength
+                (
+                    userDirection,
+                    pointLabels,
+                    isXFace,
+                    maxUserSize
+                );
+
+
+                // Determine average edge
+                scalarField avgEdgeMag(edges.size(), 0.0);
+                forAll(edges, edgei)
+                {
+                    if (isXEdge[edgei])
+                    {
+                        const edge& e = edges[edgei];
+
+                        label nSum = 0;
+
+                        // Add both contributions from e[0]. One will be the
+                        // edge itself.
+                        if (leftEdgeMag[e[0]] > 0.0)
+                        {
+                            avgEdgeMag[edgei] += leftEdgeMag[e[0]];
+                            nSum++;
+                        }
+                        if (rightEdgeMag[e[0]] > 0.0)
+                        {
+                            avgEdgeMag[edgei] += rightEdgeMag[e[0]];
+                            nSum++;
+                        }
+                        else if (maxUserSize[e[0]] > 0.0)
+                        {
+                            // Dangling point. Characteristic length from cell
+                            // size
+                            avgEdgeMag[edgei] += maxUserSize[e[0]];
+                            nSum++;
+                        }
+
+                        // Make sure we don't add the edge itself from e[1]
+                        if (leftEdge[e[1]] == edgei)
+                        {
+                            // leftEdge is edgei. Use right contribtion
+                            if (rightEdgeMag[e[1]] > 0.0)
+                            {
+                                avgEdgeMag[edgei] += rightEdgeMag[e[1]];
+                                nSum++;
+                            }
+                            else if (maxUserSize[e[1]] > 0.0)
+                            {
+                                avgEdgeMag[edgei] += maxUserSize[e[1]];
+                                nSum++;
+                            }
+                        }
+                        else if (leftEdgeMag[e[1]] > 0.0)
+                        {
+                            avgEdgeMag[edgei] += leftEdgeMag[e[1]];
+                            nSum++;
+                        }
+
+                        if (nSum > 0)
+                        {
+                            avgEdgeMag[edgei] /= nSum;
+                        }
+                    }
+                }
+
+
+                // Edge integration
+                // ~~~~~~~~~~~~~~~~
+
+                // Using edge mag
+                List<PointIntegrateData<scalar>> pointData(mesh.nPoints());
+                List<PointIntegrateData<scalar>> edgeData(mesh.nEdges());
+
+                // Using average edge mag
+                List<PointIntegrateData<scalar>> avgPointData(mesh.nPoints());
+                List<PointIntegrateData<scalar>> avgEdgeData(mesh.nEdges());
+
+                {
+                    DynamicList<label> seedPoints;
+                    DynamicList<PointIntegrateData<scalar>> seedData;
+                    DynamicList<PointIntegrateData<scalar>> avgSeedData;
+                    {
+                        forAll(pointLabels, i)
+                        {
+                            label pointi = pointLabels[i];
+                            label leftEdgei = leftEdge[pointi];
+                            if (leftEdgeMag[pointi] <= 0.0)
+                            {
+                                // Boundary to the left.
+                                seedPoints.append(pointi);
+                                seedData.append(0.0);
+                                avgSeedData.append(0.0);
+                            }
+                            else if (leftEdgei != -1)
+                            {
+                                const edge& leftE = edges[leftEdgei];
+                                label otherPointi = leftE.otherVertex(pointi);
+                                if (!isInsidePoint[otherPointi])
+                                {
+                                    seedPoints.append(pointi);
+                                    seedData.append(leftEdgeMag[pointi]);
+                                    avgSeedData.append(avgEdgeMag[leftEdgei]);
+                                }
+                            }
+                        }
+                    }
+
+                    // Prevent walking on all non-isXEdge edges by setting valid
+                    // flag.
+                    const PointIntegrateData<scalar> dummyData(-123);
+                    forAll(edges, edgei)
+                    {
+                        const edge& e = edges[edgei];
+                        if (!isXEdge[edgei])
+                        {
+                            edgeData[edgei] = dummyData;
+                            avgEdgeData[edgei] = dummyData;
+                        }
+                        else if (!isInsidePoint[e[0]] || !isInsidePoint[e[1]])
+                        {
+                            edgeData[edgei] = dummyData;
+                            avgEdgeData[edgei] = dummyData;
+                        }
+                    }
+                    forAll(isInsidePoint, pointi)
+                    {
+                        if (!isInsidePoint[pointi])
+                        {
+                            pointData[pointi] = dummyData;
+                            avgPointData[pointi] = dummyData;
+                        }
+                    }
+
+                    // Integrate edgeMag
+                    {
+                        scalarField edgeMag(edges.size());
+                        forAll(edges, edgei)
+                        {
+                            edgeMag[edgei] = edges[edgei].mag(points);
+                        }
+                        PointIntegrateData<scalar>::trackingData td(edgeMag);
+
+                        // Do all calculations
+                        PointEdgeWave
+                        <
+                            PointIntegrateData<scalar>,
+                            PointIntegrateData<scalar>::trackingData
+                        > calc
+                        (
+                            mesh,
+                            seedPoints,
+                            seedData,
+                            pointData,
+                            edgeData,
+                            returnReduce(edges.size(), sumOp<label>()),
+                            td
+                        );
+                    }
+                    // Integrate avgEdgeMag
+                    {
+                        PointIntegrateData<scalar>::trackingData td(avgEdgeMag);
+
+                        // Do all calculations
+                        PointEdgeWave
+                        <
+                            PointIntegrateData<scalar>,
+                            PointIntegrateData<scalar>::trackingData
+                        > calc
+                        (
+                            mesh,
+                            seedPoints,
+                            avgSeedData,
+                            avgPointData,
+                            avgEdgeData,
+                            returnReduce(edges.size(), sumOp<label>()),
+                            td
+                        );
+                    }
+                }
+
+
+                // Update position
+                forAll(pointLabels, zonePointi)
+                {
+                    label meshPointi = pointLabels[zonePointi];
+                    if (pointData[meshPointi].data() >= 0.0)
+                    {
+                        newPoints[meshPointi] +=
+                            userDirection
+                          * (
+                                avgPointData[meshPointi].data()
+                              - pointData[meshPointi].data()
+                            );
+                    }
+                }
+
+
+                // Smooth new position with average of non-moving points
+                {
+                    pointField unsmoothedPoints(newPoints);
+                    scalarField sumOther(mesh.nPoints(), 0.0);
+                    labelList nSumOther(mesh.nPoints(), 0);
+                    forAll(edges, edgei)
+                    {
+                        const edge& e = edges[edgei];
+                        sumOther[e[0]] +=
+                            (unsmoothedPoints[e[1]]&userDirection);
+                        nSumOther[e[0]]++;
+                        sumOther[e[1]] +=
+                            (unsmoothedPoints[e[0]]&userDirection);
+                        nSumOther[e[1]]++;
+                    }
+                    syncTools::syncPointList
+                    (
+                        mesh,
+                        sumOther,
+                        plusEqOp<scalar>(),
+                        0.0
+                    );
+                    syncTools::syncPointList
+                    (
+                        mesh,
+                        nSumOther,
+                        plusEqOp<label>(),
+                        0
+                    );
+                    for (auto pointi : pointLabels)
+                    {
+                        label leftEdgei = leftEdge[pointi];
+                        bool isSeedPoint = false;
+                        if (leftEdgeMag[pointi] <= 0.0)
+                        {
+                            isSeedPoint = true;
+                        }
+                        else if (leftEdgei != -1)
+                        {
+                            const edge& e = edges[leftEdgei];
+                            label otherPointi = e.otherVertex(pointi);
+                            if (!isInsidePoint[otherPointi])
+                            {
+                                isSeedPoint = true;
+                            }
+                        }
+
+                        if (!isSeedPoint)
+                        {
+                            scalar smoothPos =
+                                0.5
+                               *(
+                                    (unsmoothedPoints[pointi]&userDirection)
+                                   +sumOther[pointi]/nSumOther[pointi]
+                                );
+
+                            vector& v = newPoints[pointi];
+                            v += (smoothPos-(v&userDirection))*userDirection;
+                        }
+                    }
+                }
+
+
+                // Move mesh to new location
+                vectorField disp(newPoints-mesh.points());
+                syncTools::syncPointList
+                (
+                    mesh,
+                    disp,
+                    maxMagSqrEqOp<vector>(),
+                    point::zero
+                );
+                mesh.movePoints(mesh.points()+disp);
+
+                if (debug&meshRefinement::MESH)
+                {
+                    Pout<< "Writing directional smooting iteration "
+                        << iter << " mesh to time " << meshRefiner_.timeName()
+                        << endl;
+                    meshRefiner_.write
+                    (
+                        meshRefinement::debugType(debug),
+                        meshRefinement::writeType
+                        (
+                            meshRefinement::writeLevel()
+                          | meshRefinement::WRITEMESH
+                        ),
+                        mesh.time().path()/meshRefiner_.timeName()
+                    );
+                }
+            }
+        }
+    }
+    return iter;
+}
+
+
 void Foam::snappyRefineDriver::baffleAndSplitMesh
 (
     const refinementParameters& refineParams,
@@ -2285,6 +2857,12 @@ void Foam::snappyRefineDriver::doRefine
         100    // maxIter
     );
 
+    if (max(meshRefiner_.shells().nSmoothExpansion()) > 0)
+    {
+        directionalSmooth(refineParams);
+    }
+
+
     // Introduce baffles at surface intersections. Remove sections unreachable
     // from keepPoint.
     baffleAndSplitMesh

src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.H

@@ -38,6 +38,8 @@ SourceFiles
 #include "labelList.H"
 #include "PackedBoolList.H"
 #include "labelVector.H"
+#include "vector.H"
+#include "scalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -144,7 +146,7 @@ class snappyRefineDriver
             const label maxIter
         );
 
-        // Directional refinement
+        // Directional refinement and smoothing
 
             //- Refine (directional) all cells inside/outside shell
             label directionalShellRefine
@@ -153,6 +155,21 @@ class snappyRefineDriver
                 const label maxIter
             );
 
+            //- Calculate local edge length from cell volumes
+            void getVirtualEdgeLength
+            (
+                const vector& userDirection,
+                const labelList& pointLabels,
+                const PackedBoolList& isXFace,
+                scalarField& maxUserSize
+            ) const;
+
+            //- Smooth the directional expansion ratio
+            label directionalSmooth
+            (
+                const refinementParameters& refineParams
+            );
+
         //- Add baffles and remove unreachable cells
         void baffleAndSplitMesh
         (