Merge branch 'develop' of develop.openfoam.com:Development/OpenFOAM-plus into develop

applications/solvers/DNS/dnsFoam/Allwmake

@@ -1,12 +1,13 @@
 #!/bin/sh
 cd ${0%/*} || exit 1                         # Run from this directory
+. $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
 . $WM_PROJECT_DIR/wmake/scripts/have_fftw
 
 #------------------------------------------------------------------------------
 
 if have_fftw
 then
-    wmake
+    wmake $targetType
 else
     echo "==> skip dnsFoam solver (no FFTW)"
 fi

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict

@@ -306,6 +306,35 @@ castellatedMeshControls
         //    // - cellLevel will include any directional refinement
         //    //   (i.e. it will be the maximum of all three directions)
         //}
+
+        //wakeBox
+        //{
+        //    mode        inside;
+        //    // Dummy base level
+        //    levels      ((10000 0));
+        //
+        //    // Optional directional refinement (after all other refinement)
+        //    // Directional refinement
+        //    // for all cells according to 'mode' ('inside' or 'outside';
+        //    // 'distance' not supported) and within certain range. E.g.
+        //    //  - for all cells with level 2-5
+        //    //  - do one split in x direction
+        //    levelIncrement  (2 5 (1 0 0));
+        //
+        //    // Note
+        //    // - ignores 'levels' and gap* settings.
+        //    // - the cellLevel/pointLevels files are no longer consistent
+        //    //   with the mesh, the resulting mesh is no longer compatible
+        //    //   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/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C

@@ -3148,7 +3148,6 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2
         bitSet refineCell(mesh_.nCells(), newCellsToRefine);
 
         const bitSet savedRefineCell(refineCell);
-
         label nChanged = faceConsistentRefinement(true, cellLevel_, refineCell);
 
         {

src/functionObjects/field/fieldMinMax/fieldMinMax.C

@@ -123,12 +123,6 @@ Foam::functionObjects::fieldMinMax::fieldMinMax
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::functionObjects::fieldMinMax::~fieldMinMax()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 bool Foam::functionObjects::fieldMinMax::read(const dictionary& dict)

src/functionObjects/field/fieldMinMax/fieldMinMax.H

@@ -156,6 +156,14 @@ protected:
         //- No copy assignment
         void operator=(const fieldMinMax&) = delete;
 
+        //- Calculate the field min/max for a given field type
+        template<class Type>
+        void calcMinMaxFieldType
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& field,
+            const word& outputFieldName
+        );
+
         //- Calculate the field min/max
         template<class Type>
         void calcMinMaxFields
@@ -183,7 +191,7 @@ public:
 
 
     //- Destructor
-    virtual ~fieldMinMax();
+    virtual ~fieldMinMax() = default;
 
 
     // Member Functions

src/functionObjects/field/fieldMinMax/fieldMinMaxTemplates.C

@@ -111,139 +111,59 @@ void Foam::functionObjects::fieldMinMax::output
 
 
 template<class Type>
-void Foam::functionObjects::fieldMinMax::calcMinMaxFields
+void Foam::functionObjects::fieldMinMax::calcMinMaxFieldType
 (
-    const word& fieldName,
+    const GeometricField<Type, fvPatchField, volMesh>& field,
-    const modeType& mode
+    const word& outputFieldName
 )
-{
-    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
-
-    if (obr_.foundObject<fieldType>(fieldName))
 {
     const label proci = Pstream::myProcNo();
 
-        const fieldType& field = lookupObject<fieldType>(fieldName);
+    // Find min internal field value info
-
-        const volVectorField::Boundary& CfBoundary =
-            mesh_.C().boundaryField();
-
-        switch (mode)
-        {
-            case mdMag:
-            {
-                const volScalarField magField(mag(field));
-                const volScalarField::Boundary& magFieldBoundary =
-                    magField.boundaryField();
-
-                scalarList minVs(Pstream::nProcs());
-                labelList minCells(Pstream::nProcs());
-                List<vector> minCs(Pstream::nProcs());
-                label minProci = findMin(magField);
-                minVs[proci] = magField[minProci];
-                minCells[proci] = minProci;
-                minCs[proci] = mesh_.C()[minProci];
-
-                scalarList maxVs(Pstream::nProcs());
-                labelList maxCells(Pstream::nProcs());
-                List<vector> maxCs(Pstream::nProcs());
-                label maxProci = findMax(magField);
-                maxVs[proci] = magField[maxProci];
-                maxCells[proci] = maxProci;
-                maxCs[proci] = mesh_.C()[maxProci];
-
-                forAll(magFieldBoundary, patchi)
-                {
-                    const scalarField& mfp = magFieldBoundary[patchi];
-                    if (mfp.size())
-                    {
-                        const vectorField& Cfp = CfBoundary[patchi];
-
-                        const labelList& faceCells =
-                            magFieldBoundary[patchi].patch().faceCells();
-
-                        label minPi = findMin(mfp);
-                        if (mfp[minPi] < minVs[proci])
-                        {
-                            minVs[proci] = mfp[minPi];
-                            minCells[proci] = faceCells[minPi];
-                            minCs[proci] = Cfp[minPi];
-                        }
-
-                        label maxPi = findMax(mfp);
-                        if (mfp[maxPi] > maxVs[proci])
-                        {
-                            maxVs[proci] = mfp[maxPi];
-                            maxCells[proci] = faceCells[maxPi];
-                            maxCs[proci] = Cfp[maxPi];
-                        }
-                    }
-                }
-
-                Pstream::gatherList(minVs);
-                Pstream::scatterList(minVs);
-                Pstream::gatherList(minCells);
-                Pstream::scatterList(minCells);
-                Pstream::gatherList(minCs);
-                Pstream::scatterList(minCs);
-
-                Pstream::gatherList(maxVs);
-                Pstream::scatterList(maxVs);
-                Pstream::gatherList(maxCells);
-                Pstream::scatterList(maxCells);
-                Pstream::gatherList(maxCs);
-                Pstream::scatterList(maxCs);
-
-                label mini = findMin(minVs);
-                scalar minValue = minVs[mini];
-                const label minCell = minCells[mini];
-                const vector& minC = minCs[mini];
-
-                label maxi = findMax(maxVs);
-                scalar maxValue = maxVs[maxi];
-                const label maxCell = maxCells[maxi];
-                const vector& maxC = maxCs[maxi];
-
-                output
-                (
-                    fieldName,
-                    word("mag(" + fieldName + ")"),
-                    minCell,
-                    maxCell,
-                    minC,
-                    maxC,
-                    mini,
-                    maxi,
-                    minValue,
-                    maxValue
-                );
-                break;
-            }
-            case mdCmpt:
-            {
-                const typename fieldType::Boundary&
-                    fieldBoundary = field.boundaryField();
-
     List<Type> minVs(Pstream::nProcs());
     labelList minCells(Pstream::nProcs());
     List<vector> minCs(Pstream::nProcs());
+
     label minProci = findMin(field);
+    if (minProci != -1)
+    {
         minVs[proci] = field[minProci];
         minCells[proci] = minProci;
         minCs[proci] = mesh_.C()[minProci];
+    }
+    else
+    {
+        minVs[proci] = pTraits<Type>::max;
+        minCells[proci] = -1;
+        minCs[proci] = vector::max;
+    }
 
+    // Find max internal field value info
     List<Type> maxVs(Pstream::nProcs());
     labelList maxCells(Pstream::nProcs());
     List<vector> maxCs(Pstream::nProcs());
+
     label maxProci = findMax(field);
+    if (maxProci != -1)
+    {
         maxVs[proci] = field[maxProci];
         maxCells[proci] = maxProci;
         maxCs[proci] = mesh_.C()[maxProci];
+    }
+    else
+    {
+        maxVs[proci] = pTraits<Type>::min;
+        maxCells[proci] = -1;
+        maxCs[proci] = vector::max;
+    }
+
+    // Find min and max boundary field info
+    const auto& fieldBoundary = field.boundaryField();
+    const auto& CfBoundary = mesh_.C().boundaryField();
 
     forAll(fieldBoundary, patchi)
     {
         const Field<Type>& fp = fieldBoundary[patchi];
-
         if (fp.size())
         {
             const vectorField& Cfp = CfBoundary[patchi];
@@ -269,6 +189,7 @@ void Foam::functionObjects::fieldMinMax::calcMinMaxFields
         }
     }
 
+    // Collect info from all processors and output
     Pstream::gatherList(minVs);
     Pstream::scatterList(minVs);
     Pstream::gatherList(minCells);
@@ -284,19 +205,19 @@ void Foam::functionObjects::fieldMinMax::calcMinMaxFields
     Pstream::scatterList(maxCs);
 
     label mini = findMin(minVs);
-                Type minValue = minVs[mini];
+    const Type& minValue = minVs[mini];
     const label minCell = minCells[mini];
     const vector& minC = minCs[mini];
 
     label maxi = findMax(maxVs);
-                Type maxValue = maxVs[maxi];
+    const Type& maxValue = maxVs[maxi];
     const label maxCell = maxCells[maxi];
     const vector& maxC = maxCs[maxi];
 
     output
     (
-                    fieldName,
+        field.name(),
-                    fieldName,
+        outputFieldName,
         minCell,
         maxCell,
         minC,
@@ -306,6 +227,36 @@ void Foam::functionObjects::fieldMinMax::calcMinMaxFields
         minValue,
         maxValue
     );
+}
+
+
+template<class Type>
+void Foam::functionObjects::fieldMinMax::calcMinMaxFields
+(
+    const word& fieldName,
+    const modeType& mode
+)
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+
+    if (obr_.foundObject<fieldType>(fieldName))
+    {
+        const fieldType& field = lookupObject<fieldType>(fieldName);
+
+        switch (mode)
+        {
+            case mdMag:
+            {
+                calcMinMaxFieldType<scalar>
+                (
+                    mag(field),
+                    word("mag(" + fieldName + ")")
+                );
+                break;
+            }
+            case mdCmpt:
+            {
+                calcMinMaxFieldType(field, fieldName);
                 break;
             }
             default:

src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H

@@ -1053,7 +1053,6 @@ public:
                 const scalar maxLoadUnbalance
             );
 
-
             //- Calculate list of cells to directionally refine
             labelList directionalRefineCandidates
             (

src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C

@@ -595,6 +595,9 @@ Foam::shellSurfaces::shellSurfaces
     distances_.setSize(shellI);
     levels_.setSize(shellI);
     dirLevels_.setSize(shellI);
+    smoothDirection_.setSize(shellI);
+    nSmoothExpansion_.setSize(shellI);
+    nSmoothPosition_.setSize(shellI);
 
     extendedGapLevel_.setSize(shellI);
     extendedGapMode_.setSize(shellI);
@@ -671,6 +674,19 @@ Foam::shellSurfaces::shellSurfaces
                 }
             }
 
+            // Directional smoothing
+            // ~~~~~~~~~~~~~~~~~~~~~
+
+            nSmoothExpansion_[shellI] = 0;
+            nSmoothPosition_[shellI] = 0;
+            smoothDirection_[shellI] =
+                dict.lookupOrDefault("smoothDirection", vector::zero);
+
+            if (smoothDirection_[shellI] != vector::zero)
+            {
+                dict.lookup("nSmoothExpansion") >> nSmoothExpansion_[shellI];
+                dict.lookup("nSmoothPosition") >> nSmoothPosition_[shellI];
+            }
 
 
             // Gap specification
@@ -797,6 +813,24 @@ Foam::labelPairList Foam::shellSurfaces::directionalSelectLevel() const
 }
 
 
+const Foam::labelList& Foam::shellSurfaces::nSmoothExpansion() const
+{
+    return nSmoothExpansion_;
+}
+
+
+const Foam::vectorField& Foam::shellSurfaces::smoothDirection() const
+{
+    return smoothDirection_;
+}
+
+
+const Foam::labelList& Foam::shellSurfaces::nSmoothPosition() const
+{
+    return nSmoothPosition_;
+}
+
+
 void Foam::shellSurfaces::findHigherLevel
 (
     const pointField& pt,

src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.H

@@ -86,9 +86,21 @@ private:
         //- Per shell per distance the refinement level
         labelListList levels_;
 
+
+        // 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_;
+
+            //- Per shell the positional smoothing iterations
+            labelList nSmoothPosition_;
+
 
         // Gap level refinement
 
@@ -231,6 +243,15 @@ 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;
+
+            //- Per shell the positional smoothing iterations
+            const labelList& nSmoothPosition() 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/pointData/pointData.C

@@ -1,53 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     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 "pointData.H"
-
-// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
-
-Foam::Ostream& Foam::operator<<(Ostream& os, const pointData& wDist)
-{
-    if (os.format() == IOstream::ASCII)
-    {
-        return os
-            << static_cast<const pointEdgePoint&>(wDist)
-            << token::SPACE << wDist.s() << token::SPACE << wDist.v();
-    }
-    else
-    {
-        return os
-            << static_cast<const pointEdgePoint&>(wDist)
-            << wDist.s() << wDist.v();
-    }
-}
-
-
-Foam::Istream& Foam::operator>>(Istream& is, pointData& wDist)
-{
-    return is >> static_cast<pointEdgePoint&>(wDist) >> wDist.s_ >> wDist.v_;
-}
-
-
-// ************************************************************************* //

src/mesh/snappyHexMesh/snappyHexMeshDriver/pointData/pointData.H

@@ -1,193 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     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::pointData
-
-Description
-    Variant of pointEdgePoint with some transported additional data.
-    WIP - should be templated on data like wallDistData.
-    Passive vector v_ is not a coordinate (so no enterDomain/leaveDomain
-    transformation needed)
-
-SourceFiles
-    pointDataI.H
-    pointData.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef pointData_H
-#define pointData_H
-
-#include "pointEdgePoint.H"
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Forward declaration of friend functions and operators
-
-class pointData;
-
-Istream& operator>>(Istream&, pointData&);
-Ostream& operator<<(Ostream&, const pointData&);
-
-
-/*---------------------------------------------------------------------------*\
-                           Class pointData Declaration
-\*---------------------------------------------------------------------------*/
-
-class pointData
-:
-    public pointEdgePoint
-{
-    // Private data
-
-        //- Additional information.
-        scalar s_;
-
-        //- Additional information.
-        vector v_;
-
-public:
-
-    // Constructors
-
-        //- Construct null
-        inline pointData();
-
-        //- Construct from origin, distance
-        inline pointData
-        (
-            const point& origin,
-            const scalar distSqr,
-            const scalar s,
-            const vector& v
-        );
-
-        //- Construct as copy
-        inline pointData(const pointData&);
-
-
-    // Member Functions
-
-        // Access
-
-            inline scalar s() const;
-
-            inline const vector& v() const;
-
-
-        // Needed by meshWave
-
-            //- Apply rotation matrix to origin
-            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 pointData& 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 pointData& 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 pointData& 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 pointData& pointInfo,
-                const scalar tol,
-                TrackingData& td
-            );
-
-    // Member Operators
-
-        // Needed for List IO
-        inline bool operator==(const pointData&) const;
-        inline bool operator!=(const pointData&) const;
-
-
-    // IOstream Operators
-
-        friend Ostream& operator<<(Ostream&, const pointData&);
-        friend Istream& operator>>(Istream&, pointData&);
-};
-
-
-//- Data associated with pointData as contiguous as pointEdgePoint
-template<>
-inline bool contiguous<pointData>()
-{
-    return contiguous<pointEdgePoint>();
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "pointDataI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/mesh/snappyHexMesh/snappyHexMeshDriver/pointData/pointDataI.H

@@ -1,237 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     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  * * * * * * * * * * * * * * //
-
-// Null constructor
-inline Foam::pointData::pointData()
-:
-    pointEdgePoint(),
-    s_(GREAT),
-    v_(point::max)
-{}
-
-
-// Construct from origin, distance
-inline Foam::pointData::pointData
-(
-    const point& origin,
-    const scalar distSqr,
-    const scalar s,
-    const vector& v
-)
-:
-    pointEdgePoint(origin, distSqr),
-    s_(s),
-    v_(v)
-{}
-
-
-// Construct as copy
-inline Foam::pointData::pointData(const pointData& wpt)
-:
-    pointEdgePoint(wpt),
-    s_(wpt.s()),
-    v_(wpt.v())
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-inline Foam::scalar Foam::pointData::s() const
-{
-    return s_;
-}
-
-
-inline const Foam::vector& Foam::pointData::v() const
-{
-    return v_;
-}
-
-
-template<class TrackingData>
-inline void Foam::pointData::transform
-(
-    const tensor& rotTensor,
-    TrackingData& td
-)
-{
-    pointEdgePoint::transform(rotTensor, td);
-    v_ = Foam::transform(rotTensor, v_);
-}
-
-
-// Update this with information from connected edge
-template<class TrackingData>
-inline bool Foam::pointData::updatePoint
-(
-    const polyMesh& mesh,
-    const label pointi,
-    const label edgeI,
-    const pointData& edgeInfo,
-    const scalar tol,
-    TrackingData& td
-)
-{
-    if
-    (
-        pointEdgePoint::updatePoint
-        (
-            mesh,
-            pointi,
-            edgeI,
-            edgeInfo,
-            tol,
-            td
-        )
-    )
-    {
-        s_ = edgeInfo.s_;
-        v_ = edgeInfo.v_;
-        return true;
-    }
-    else
-    {
-        return false;
-    }
-}
-
-// Update this with new information on same point
-template<class TrackingData>
-inline bool Foam::pointData::updatePoint
-(
-    const polyMesh& mesh,
-    const label pointi,
-    const pointData& newPointInfo,
-    const scalar tol,
-    TrackingData& td
-)
-{
-    if
-    (
-        pointEdgePoint::updatePoint
-        (
-            mesh,
-            pointi,
-            newPointInfo,
-            tol,
-            td
-        )
-    )
-    {
-        s_ = newPointInfo.s_;
-        v_ = newPointInfo.v_;
-        return true;
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
-// Update this with new information on same point. No extra information.
-template<class TrackingData>
-inline bool Foam::pointData::updatePoint
-(
-    const pointData& newPointInfo,
-    const scalar tol,
-    TrackingData& td
-)
-{
-    if (pointEdgePoint::updatePoint(newPointInfo, tol, td))
-    {
-        s_ = newPointInfo.s_;
-        v_ = newPointInfo.v_;
-        return true;
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
-// Update this with information from connected point
-template<class TrackingData>
-inline bool Foam::pointData::updateEdge
-(
-    const polyMesh& mesh,
-    const label edgeI,
-    const label pointi,
-    const pointData& pointInfo,
-    const scalar tol,
-    TrackingData& td
-
-)
-{
-    if
-    (
-        pointEdgePoint::updateEdge
-        (
-            mesh,
-            edgeI,
-            pointi,
-            pointInfo,
-            tol,
-            td
-        )
-    )
-    {
-        s_ = pointInfo.s_;
-        v_ = pointInfo.v_;
-        return true;
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
-inline bool Foam::pointData::operator==(const Foam::pointData& rhs)
-const
-{
-    return
-        pointEdgePoint::operator==(rhs)
-     && (s() == rhs.s())
-     && (v() == rhs.v());
-}
-
-
-inline bool Foam::pointData::operator!=(const Foam::pointData& rhs)
-const
-{
-    return !(*this == rhs);
-}
-
-
-// ************************************************************************* //

src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C

@@ -41,6 +41,8 @@ License
 #include "labelVector.H"
 #include "profiling.H"
 #include "searchableSurfaces.H"
+#include "fvMeshSubset.H"
+#include "interpolationTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -1642,7 +1644,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
@@ -1800,6 +1801,517 @@ Foam::label Foam::snappyRefineDriver::directionalShellRefine
 }
 
 
+void Foam::snappyRefineDriver::mergeAndSmoothRatio
+(
+    const scalarList& allSeedPointDist,
+    const label nSmoothExpansion,
+    List<Tuple2<scalar, scalar>>&  keyAndValue
+)
+{
+    // Merge duplicate distance from coupled locations to get unique
+    // distances to operate on, do on master
+    SortableList<scalar> unmergedDist(allSeedPointDist);
+    DynamicList<scalar> mergedDist;
+
+    scalar prevDist = GREAT;
+    forAll(unmergedDist, i)
+    {
+        scalar curDist = unmergedDist[i];
+        scalar difference = mag(curDist - prevDist);
+        if (difference > meshRefiner_.mergeDistance())
+        //if (difference > 0.01)
+        {
+             mergedDist.append(curDist);
+             prevDist = curDist;
+        }
+    }
+
+    // Sort the unique distances
+    SortableList<scalar> sortedDist(mergedDist);
+    labelList indexSet = sortedDist.indices();
+
+    // Get updated position starting from original (undistorted) mesh
+    scalarList seedPointsNewLocation = sortedDist;
+
+    scalar initResidual = 0.0;
+    scalar prevIterResidual = GREAT;
+
+    for (label iter = 0; iter < nSmoothExpansion; iter++)
+    {
+
+        // Position based edge averaging algorithm operated on
+        // all seed plane locations in normalized form.
+        //
+        //   0   1   2   3   4   5   6  (edge numbers)
+        //  ---x---x---x---x---x---x---
+        //     0   1   2   3   4   5    (point numbers)
+        //
+        // Average of edge 1-3 in terms of position
+        //  = (point3 - point0)/3
+        // Keeping points 0-1 frozen, new position of point 2
+        //  = position2 + (average of edge 1-3 as above)
+        for(label i = 2; i<mergedDist.size()-1; i++)
+        {
+            scalar oldX00 = sortedDist[i-2];
+            scalar oldX1 = sortedDist[i+1];
+            scalar curX0 = seedPointsNewLocation[i-1];
+            seedPointsNewLocation[i] = curX0 + (oldX1 - oldX00)/3;
+        }
+
+        const scalarField residual(seedPointsNewLocation-sortedDist);
+        {
+            scalar res(sumMag(residual));
+
+            if (iter == 0)
+            {
+                initResidual = res;
+            }
+            res /= initResidual;
+
+            if (mag(prevIterResidual - res) < SMALL)
+            {
+                if (debug)
+                {
+                    Pout<< "Converged with iteration " << iter 
+                        << " initResidual: " << initResidual
+                        << " final residual : " << res << endl;
+                }
+                break;
+            }
+            else
+            {
+                prevIterResidual = res;
+            }
+        }
+
+        // Update the field for next iteration, avoid moving points
+        sortedDist = seedPointsNewLocation;
+
+    }
+
+    keyAndValue.setSize(mergedDist.size());
+
+    forAll(mergedDist, i)
+    {
+        keyAndValue[i].first() = mergedDist[i];
+        label index = indexSet[i];
+        keyAndValue[i].second() = seedPointsNewLocation[index];
+    }
+}
+
+
+Foam::label Foam::snappyRefineDriver::directionalSmooth
+(
+    const refinementParameters& refineParams
+)
+{
+    addProfiling(split, "snappyHexMesh::refine::smooth");
+    Info<< nl
+        << "Directional expansion ratio smoothing" << nl
+        << "-------------------------------------" << nl
+        << endl;
+
+    fvMesh& baseMesh = meshRefiner_.mesh();
+    const searchableSurfaces& geometry = meshRefiner_.surfaces().geometry();
+    const shellSurfaces& shells = meshRefiner_.shells();
+
+    label iter = 0;
+
+    forAll(shells.nSmoothExpansion(), shellI)
+    {
+        if 
+        (
+            shells.nSmoothExpansion()[shellI] > 0
+         || shells.nSmoothPosition()[shellI] > 0
+        )
+        {
+            label surfi = shells.shells()[shellI];
+            const vector& userDirection = shells.smoothDirection()[shellI];
+
+
+            // Extract inside points
+            labelList pointLabels;
+            {
+                // Get inside points
+                List<volumeType> volType;
+                geometry[surfi].getVolumeType(baseMesh.points(), volType);
+
+                label nInside = 0;
+                forAll(volType, pointi)
+                {
+                    if (volType[pointi] == volumeType::INSIDE)
+                    {
+                        nInside++;
+                    }
+                }
+                pointLabels.setSize(nInside);
+                nInside = 0;
+                forAll(volType, pointi)
+                {
+                    if (volType[pointi] == volumeType::INSIDE)
+                    {
+                        pointLabels[nInside++] = pointi;
+                    }
+                }
+
+                //bitSet isInsidePoint(baseMesh.nPoints());
+                //forAll(volType, pointi)
+                //{
+                //    if (volType[pointi] == volumeType::INSIDE)
+                //    {
+                //        isInsidePoint.set(pointi);
+                //    }
+                //}
+                //pointLabels = isInsidePoint.used();
+            }
+
+            // Mark all directed edges
+            bitSet isXEdge(baseMesh.edges().size());
+            {
+                const edgeList& edges = baseMesh.edges();
+                forAll(edges, edgei)
+                {
+                    const edge& e = edges[edgei];
+                    vector eVec(e.vec(baseMesh.points()));
+                    eVec /= mag(eVec);
+                    if (mag(eVec&userDirection) > 0.9)
+                    {
+                        isXEdge.set(edgei);
+                    }
+                }
+            }
+
+            // Get the extreme of smoothing region and 
+            // normalize all points within
+            const scalar totalLength = 
+                geometry[surfi].bounds().span()
+              & userDirection;
+            const scalar startPosition = 
+                geometry[surfi].bounds().min() 
+              & userDirection;
+
+            scalarField normalizedPosition(pointLabels.size(), 0);
+            forAll(pointLabels, i)
+            {
+                label pointi = pointLabels[i];
+                normalizedPosition[i] = 
+                  (
+                    ((baseMesh.points()[pointi]&userDirection) - startPosition)
+                  / totalLength
+                  );
+            }
+
+            // Sort the normalized position
+            labelList order;
+            sortedOrder(normalizedPosition, order);
+
+            DynamicList<scalar> seedPointDist;
+
+            // Select points from finest refinement (one point-per plane)
+            scalar prevDist = GREAT;
+            forAll(order, i)
+            {
+                label pointi = order[i];
+                scalar curDist = normalizedPosition[pointi];
+                if (mag(curDist - prevDist) > meshRefiner_.mergeDistance())
+                {
+                    seedPointDist.append(curDist);
+                    prevDist = curDist;
+                }
+            }
+
+            // Collect data from all processors
+            scalarList allSeedPointDist;
+            {
+                List<scalarList> gatheredDist(Pstream::nProcs());
+                gatheredDist[Pstream::myProcNo()] = seedPointDist;
+                Pstream::gatherList(gatheredDist);
+
+                // Combine processor lists into one big list.
+                allSeedPointDist =
+                    ListListOps::combine<scalarList>
+                    (
+                        gatheredDist, accessOp<scalarList>()
+                    );
+            }
+
+            // Pre-set the points not to smooth (after expansion)
+            bitSet isFrozenPoint(baseMesh.nPoints(), true);
+
+            {
+                scalar minSeed = min(allSeedPointDist);
+                Pstream::scatter(minSeed);
+                scalar maxSeed = max(allSeedPointDist);
+                Pstream::scatter(maxSeed);
+
+                forAll(normalizedPosition, posI)
+                {
+                    const scalar pos = normalizedPosition[posI];
+                    if
+                    (
+                        (mag(pos-minSeed) < meshRefiner_.mergeDistance())
+                     || (mag(pos-maxSeed) < meshRefiner_.mergeDistance())
+                    )
+                    {
+                        // Boundary point: freeze
+                        isFrozenPoint.set(pointLabels[posI]);
+                    }
+                    else
+                    {
+                        // Internal to moving region
+                        isFrozenPoint.unset(pointLabels[posI]);
+                    }
+                }
+            }
+
+            Info<< "Smoothing " << geometry[surfi].name() << ':' << nl
+                << "    Direction                   : " << userDirection << nl
+                << "    Number of points            : "
+                << returnReduce(pointLabels.size(), sumOp<label>())
+                << " (out of " << baseMesh.globalData().nTotalPoints()
+                << ")" << nl
+                << "    Smooth expansion iterations : "
+                << shells.nSmoothExpansion()[shellI] << nl
+                << "    Smooth position iterations  : "
+                << shells.nSmoothPosition()[shellI] << nl
+                << "    Number of planes            : "
+                << allSeedPointDist.size()
+                << endl;
+
+            // Make lookup from original normalized distance to new value
+            List<Tuple2<scalar, scalar>> keyAndValue(allSeedPointDist.size());
+
+            // Filter unique seed distances and iterate for user given steps
+            // or until convergence. Then get back map from old to new distance
+            if (Pstream::master())
+            {
+                mergeAndSmoothRatio
+                (
+                    allSeedPointDist,
+                    shells.nSmoothExpansion()[shellI],
+                    keyAndValue
+                );
+            }
+
+            Pstream::scatter(keyAndValue);
+
+            // Construct an iterpolation table for further queries
+            // - although normalized values are used for query,
+            //   it might flow out of bounds due to precision, hence clamped
+            const interpolationTable<scalar> table
+            (
+                keyAndValue,
+                bounds::repeatableBounding::CLAMP,
+                "undefined"
+            );
+
+            // Now move the points directly on the baseMesh.
+            pointField baseNewPoints(baseMesh.points());
+            forAll(pointLabels, i)
+            {
+                label pointi = pointLabels[i];
+                const point& curPoint = baseMesh.points()[pointi];
+                scalar curDist = normalizedPosition[i];
+                scalar newDist = table(curDist);
+                scalar newPosition = startPosition + newDist*totalLength;
+                baseNewPoints[pointi] +=
+                    userDirection * (newPosition - (curPoint &userDirection));
+            }
+
+            // Moving base mesh with expansion ratio smoothing
+            vectorField disp(baseNewPoints-baseMesh.points());
+            syncTools::syncPointList
+            (
+                baseMesh,
+                disp,
+                maxMagSqrEqOp<vector>(),
+                point::zero
+            );
+            baseMesh.movePoints(baseMesh.points()+disp);
+
+            if (debug&meshRefinement::MESH)
+            {
+                const_cast<Time&>(baseMesh.time())++;
+
+                Pout<< "Writing directional expansion ratio smoothed"
+                    << " mesh to time " << meshRefiner_.timeName() << endl;
+
+                meshRefiner_.write
+                (
+                    meshRefinement::debugType(debug),
+                    meshRefinement::writeType
+                    (
+                        meshRefinement::writeLevel()
+                      | meshRefinement::WRITEMESH
+                    ),
+                    baseMesh.time().path()/meshRefiner_.timeName()
+                );
+            }
+
+            // Now we have moved the points in user specified region. Smooth
+            // them with neighbour points to avoid skewed cells
+            // Instead of moving actual mesh, operate on copy
+            pointField baseMeshPoints(baseMesh.points());
+            scalar initResidual = 0.0;
+            scalar prevIterResidual = GREAT;
+            for (iter = 0; iter < shells.nSmoothPosition()[shellI]; iter++)
+            {
+                {
+                    const edgeList& edges = baseMesh.edges();
+                    const labelListList& pointEdges = baseMesh.pointEdges();
+
+                    pointField unsmoothedPoints(baseMeshPoints);
+
+                    scalarField sumOther(baseMesh.nPoints(), 0.0);
+                    labelList nSumOther(baseMesh.nPoints(), 0);
+                    labelList nSumXEdges(baseMesh.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]]++;
+                        if (isXEdge[edgei])
+                        {
+                            nSumXEdges[e[0]]++;
+                            nSumXEdges[e[1]]++;
+                        }
+                    }
+
+                    syncTools::syncPointList
+                    (
+                        baseMesh,
+                        nSumXEdges,
+                        plusEqOp<label>(),
+                        0
+                    );
+
+                    forAll(pointLabels, i)
+                    {
+                        label pointi = pointLabels[i];
+
+                        if (nSumXEdges[pointi] < 2)
+                        {
+                            // Hanging node. Remove the (single!) X edge so it
+                            // will follow points above or below instead
+                            const labelList& pEdges = pointEdges[pointi];
+                            forAll(pEdges, pE)
+                            {
+                                label edgei = pEdges[pE];
+                                if (isXEdge[edgei])
+                                {
+                                    const edge& e = edges[edgei];
+                                    label otherPt = e.otherVertex(pointi);
+                                    nSumOther[pointi]--;
+                                    sumOther[pointi] -=
+                                      (
+                                          unsmoothedPoints[otherPt]
+                                        & userDirection
+                                      );
+                                }
+                            }
+                        }
+                    }
+
+                    syncTools::syncPointList
+                    (
+                        baseMesh,
+                        sumOther,
+                        plusEqOp<scalar>(),
+                        0.0
+                    );
+                    syncTools::syncPointList
+                    (
+                        baseMesh,
+                        nSumOther,
+                        plusEqOp<label>(),
+                        0
+                    );
+
+                    forAll(pointLabels, i)
+                    {
+                        label pointi = pointLabels[i];
+
+                        if ((nSumOther[pointi] >= 2) && !isFrozenPoint[pointi])
+                        {
+                            scalar smoothPos =
+                                0.5
+                               *(
+                                    (unsmoothedPoints[pointi]&userDirection)
+                                   +sumOther[pointi]/nSumOther[pointi]
+                                );
+
+                            vector& v = baseNewPoints[pointi];
+                            v += (smoothPos-(v&userDirection))*userDirection;
+                        }
+                    }
+
+                    const vectorField residual(baseNewPoints - baseMeshPoints);
+                    {
+                        scalar res(gSum(mag(residual)));
+
+                        if (iter == 0)
+                        {
+                            initResidual = res;
+                        }
+                        res /= initResidual;
+
+                        if (mag(prevIterResidual - res) < SMALL)
+                        {
+                            Info<< "Converged smoothing in iteration " << iter
+                                << " initResidual: " << initResidual
+                                << " final residual : " << res << endl;
+                            break;
+                        }
+                        else
+                        {
+                            prevIterResidual = res;
+                        }
+                    }
+
+                    // Just copy new location instead of moving base mesh
+                    baseMeshPoints = baseNewPoints;
+                }
+            }
+
+            // Move mesh to new location
+            vectorField dispSmooth(baseMeshPoints-baseMesh.points());
+            syncTools::syncPointList
+            (
+                baseMesh,
+                dispSmooth,
+                maxMagSqrEqOp<vector>(),
+                point::zero
+            );
+            baseMesh.movePoints(baseMesh.points()+dispSmooth);
+
+            if (debug&meshRefinement::MESH)
+            {
+                const_cast<Time&>(baseMesh.time())++;
+
+                Pout<< "Writing positional smoothing iteration "
+                    << iter << " mesh to time " << meshRefiner_.timeName()
+                    << endl;
+                meshRefiner_.write
+                (
+                    meshRefinement::debugType(debug),
+                    meshRefinement::writeType
+                    (
+                        meshRefinement::writeLevel()
+                      | meshRefinement::WRITEMESH
+                    ),
+                    baseMesh.time().path()/meshRefiner_.timeName()
+                );
+            }
+        }
+    }
+    return iter;
+}
+
+
 void Foam::snappyRefineDriver::baffleAndSplitMesh
 (
     const refinementParameters& refineParams,
@@ -2307,6 +2819,16 @@ void Foam::snappyRefineDriver::doRefine
         100    // maxIter
     );
 
+    if 
+    (
+        max(meshRefiner_.shells().nSmoothExpansion()) > 0
+     || max(meshRefiner_.shells().nSmoothPosition()) > 0
+    )
+    {
+        directionalSmooth(refineParams);
+    }
+
+
     // Introduce baffles at surface intersections. Remove sections unreachable
     // from keepPoint.
     baffleAndSplitMesh

src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -36,6 +36,8 @@ SourceFiles
 
 #include "wordPairHashTable.H"
 #include "labelList.H"
+#include "scalarField.H"
+#include "Tuple2.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -49,6 +51,7 @@ class snapParameters;
 class meshRefinement;
 class decompositionMethod;
 class fvMeshDistribute;
+class fvMesh;
 
 /*---------------------------------------------------------------------------*\
                            Class snappyRefineDriver Declaration
@@ -142,7 +145,7 @@ class snappyRefineDriver
             const label maxIter
         );
 
-        // Directional refinement
+        // Directional refinement and smoothing
 
             //- Refine (directional) all cells inside/outside shell
             label directionalShellRefine
@@ -151,6 +154,18 @@ class snappyRefineDriver
                 const label maxIter
             );
 
+            //- Calculate local edge length from cell volumes
+            void mergeAndSmoothRatio
+            (
+                const scalarList& allSeedPointDist,
+                const label nSmoothExpansion,
+                List<Tuple2<scalar, scalar>>&  keyAndValue
+            );
+
+            //- Smooth the directional expansion ratio
+            label directionalSmooth(const refinementParameters& refineParams);
+
+
         //- Add baffles and remove unreachable cells
         void baffleAndSplitMesh
         (
@@ -229,7 +244,6 @@ public:
             const refinementParameters& refineParams,
             const HashTable<Pair<word>>& faceZoneToPatches
         );
-
 };
 
 

src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C

@@ -28,6 +28,7 @@ License
 #include "meshTools.H"
 #include "mapDistribute.H"
 #include "flipOp.H"
+#include "profiling.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -100,10 +101,9 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::projectPointsToSurface
     pointField& pts
 ) const
 {
-    if (debug)
+    addProfiling(ami, "AMIInterpolation::projectPointsToSurface");
-    {
+
-        Info<< "AMI: projecting points to surface" << endl;
+    DebugInfo<< "AMI: projecting points to surface" << endl;
-    }
 
     List<pointIndexHit> nearInfo;
 
@@ -120,8 +120,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::projectPointsToSurface
         }
         else
         {
-            pts[i] = pts[i];
+            // Point remains unchanged
-            nMiss++;
+            ++nMiss;
         }
     }
 
@@ -148,6 +148,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::normaliseWeights
     const scalar lowWeightTol
 )
 {
+    addProfiling(ami, "AMIInterpolation::normaliseWeights");
+
     // Normalise the weights
     wghtSum.setSize(wght.size(), 0.0);
     label nLowWeight = 0;
@@ -233,6 +235,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::agglomerate
     autoPtr<mapDistribute>& tgtMap
 )
 {
+    addProfiling(ami, "AMIInterpolation::agglomerate");
+
     label sourceCoarseSize =
     (
         sourceRestrictAddressing.size()
@@ -329,7 +333,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::agglomerate
                     {
                         oldToNew[coarseElem] = newi;
                         newSubMap[newi] = coarseElem;
-                        newi++;
+                        ++newi;
                     }
                 }
                 newSubMap.setSize(newi);
@@ -407,7 +411,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::agglomerate
                             oldToNew[coarseElem] = newi;
                             tgtCompactMap[fineElem] = compacti;
                             newConstructMap[newi] = compacti++;
-                            newi++;
+                            ++newi;
                         }
                         else
                         {
@@ -843,6 +847,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
     const TargetPatch& tgtPatch
 )
 {
+    addProfiling(ami, "AMIInterpolation::update");
+
     label srcTotalSize = returnReduce(srcPatch.size(), sumOp<label>());
     label tgtTotalSize = returnReduce(tgtPatch.size(), sumOp<label>());
 
@@ -1057,6 +1063,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::append
     const TargetPatch& tgtPatch
 )
 {
+    addProfiling(ami, "AMIInterpolation::append");
+
     // Create a new interpolation
     autoPtr<AMIInterpolation<SourcePatch, TargetPatch>> newPtr
     (
@@ -1283,6 +1291,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolateToTarget
     const UList<Type>& defaultValues
 ) const
 {
+    addProfiling(ami, "AMIInterpolation::interpolateToTarget");
+
     if (fld.size() != srcAddress_.size())
     {
         FatalErrorInFunction
@@ -1368,6 +1378,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolateToSource
     const UList<Type>& defaultValues
 ) const
 {
+    addProfiling(ami, "AMIInterpolation::interpolateToSource");
+
     if (fld.size() != tgtAddress_.size())
     {
         FatalErrorInFunction

src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/faceAreaWeightAMI/faceAreaWeightAMI.C

@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "faceAreaWeightAMI.H"
+#include "profiling.H"
 
 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
 
@@ -38,6 +39,8 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::calcAddressing
     label tgtFacei
 )
 {
+    addProfiling(ami, "faceAreaWeightAMI::calcAddressing");
+
     // construct weights and addressing
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -123,6 +126,8 @@ bool Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::processSourceFace
     List<DynamicList<scalar>>& tgtWght
 )
 {
+    addProfiling(ami, "faceAreaWeightAMI::processSourceFace");
+
     if (tgtStartFacei == -1)
     {
         return false;
@@ -197,6 +202,8 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::setNextFaces
     bool errorOnNotFound
 ) const
 {
+    addProfiling(ami, "faceAreaWeightAMI::setNextFaces");
+
     const labelList& srcNbrFaces = this->srcPatch_.faceFaces()[srcFacei];
 
     // initialise tgtFacei
@@ -303,6 +310,8 @@ Foam::scalar Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::interArea
     const label tgtFacei
 ) const
 {
+    addProfiling(ami, "faceAreaWeightAMI::interArea");
+
     scalar area = 0;
 
     const pointField& srcPoints = this->srcPatch_.points();
@@ -444,6 +453,8 @@ restartUncoveredSourceFace
     List<DynamicList<scalar>>& tgtWght
 )
 {
+    addProfiling(ami, "faceAreaWeightAMI::restartUncoveredSourceFace");
+
     // Collect all src faces with a low weight
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -581,6 +592,8 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::calculate
     label tgtFacei
 )
 {
+    addProfiling(ami, "faceAreaWeightAMI::calculate");
+
     bool ok =
         this->initialise
         (

src/meshTools/meshSearch/meshSearch.C

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright 2015 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright 2015-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -38,6 +38,85 @@ namespace Foam
     defineTypeNameAndDebug(meshSearch, 0);
 
     scalar meshSearch::tol_ = 1e-3;
+
+    // Intersection operation that checks previous successful hits so that they
+    // are not duplicated
+    class findUniqueIntersectOp
+    :
+        public treeDataFace::findIntersectOp
+    {
+    public:
+
+        const indexedOctree<treeDataFace>& tree_;
+
+        const List<pointIndexHit>& hits_;
+
+    public:
+
+        //- Construct from components
+        findUniqueIntersectOp
+        (
+            const indexedOctree<treeDataFace>& tree,
+            const List<pointIndexHit>& hits
+        )
+        :
+            treeDataFace::findIntersectOp(tree),
+            tree_(tree),
+            hits_(hits)
+        {}
+
+        //- Calculate intersection of triangle with ray. Sets result
+        //  accordingly
+        bool operator()
+        (
+            const label index,
+            const point& start,
+            const point& end,
+            point& intersectionPoint
+        ) const
+        {
+            const primitiveMesh& mesh = tree_.shapes().mesh();
+
+            // Check whether this hit has already happened. If the new face
+            // index is the same as an existing hit then return no new hit. If
+            // the new face shares a point with an existing hit face and the
+            // line passes through both faces in the same direction, then this
+            // is also assumed to be a duplicate hit.
+            const label newFacei = tree_.shapes().faceLabels()[index];
+            const face& newFace = mesh.faces()[newFacei];
+            const scalar newDot = mesh.faceAreas()[newFacei] & (end - start);
+            forAll(hits_, hiti)
+            {
+                const label oldFacei = hits_[hiti].index();
+                const face& oldFace = mesh.faces()[oldFacei];
+                const scalar oldDot =
+                    mesh.faceAreas()[oldFacei] & (end - start);
+
+                if
+                (
+                    hits_[hiti].index() == newFacei
+                 || (
+                        newDot*oldDot > 0
+                     && (labelHashSet(newFace) & labelHashSet(oldFace)).size()
+                    )
+                )
+                {
+                    return false;
+                }
+            }
+
+            const bool hit =
+                treeDataFace::findIntersectOp::operator()
+                (
+                    index,
+                    start,
+                    end,
+                    intersectionPoint
+                );
+
+            return hit;
+        }
+    };
 }
 
 
@@ -466,25 +545,6 @@ Foam::label Foam::meshSearch::findNearestBoundaryFaceWalk
 }
 
 
-Foam::vector Foam::meshSearch::offset
-(
-    const point& bPoint,
-    const label bFacei,
-    const vector& dir
-) const
-{
-    // Get the neighbouring cell
-    label ownerCelli = mesh_.faceOwner()[bFacei];
-
-    const point& c = mesh_.cellCentres()[ownerCelli];
-
-    // Typical dimension: distance from point on face to cell centre
-    scalar typDim = mag(c - bPoint);
-
-    return tol_*typDim*dir;
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::meshSearch::meshSearch
@@ -808,40 +868,20 @@ Foam::List<Foam::pointIndexHit> Foam::meshSearch::intersections
 ) const
 {
     DynamicList<pointIndexHit> hits;
+    findUniqueIntersectOp iop(boundaryTree(), hits);
 
-    const vector edgeVec = normalised(pEnd - pStart);
+    while (true)
-
-    point pt = pStart;
-
-    pointIndexHit bHit;
-    do
     {
-        bHit = intersection(pt, pEnd);
+        // Get the next hit, or quit
+        pointIndexHit curHit = boundaryTree().findLine(pStart, pEnd, iop);
+        if (!curHit.hit()) break;
 
-        if (bHit.hit())
+        // Change index into octreeData into face label
-        {
+        curHit.setIndex(boundaryTree().shapes().faceLabels()[curHit.index()]);
-            hits.append(bHit);
 
-            const vector& area = mesh_.faceAreas()[bHit.index()];
+        hits.append(curHit);
-
-            scalar typDim = Foam::sqrt(mag(area));
-
-            if ((mag(bHit.hitPoint() - pEnd)/typDim) < SMALL)
-            {
-                break;
     }
 
-            // Restart from hitPoint shifted a little bit in the direction
-            // of the destination
-
-            pt =
-                bHit.hitPoint()
-              + offset(bHit.hitPoint(), bHit.index(), edgeVec);
-        }
-
-    } while (bHit.hit());
-
-
     hits.shrink();
 
     return hits;

src/meshTools/meshSearch/meshSearch.H

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -135,15 +135,6 @@ class meshSearch
                 const label seedFacei
             ) const;
 
-            //- Calculate offset vector in direction dir with as length a
-            //  fraction of the cell size (of the cell straddling boundary face)
-            vector offset
-            (
-                const point& bPoint,
-                const label bFacei,
-                const vector& dir
-            ) const;
-
 
         //- No copy construct
         meshSearch(const meshSearch&) = delete;

src/overset/cellCellStencil/cellCellStencil/cellCellStencilObject.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2017-2018 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -55,12 +55,7 @@ typedef MeshObject
 
 class cellCellStencilObject
 :
-    public MeshObject
+    public Stencil,
-    <
-        fvMesh,
-        MoveableMeshObject,
-        cellCellStencilObject
-    >,
     public cellCellStencil
 {
     // Private data

src/overset/dynamicOversetFvMesh/dynamicOversetFvMesh.C

@@ -44,6 +44,15 @@ namespace Foam
 bool Foam::dynamicOversetFvMesh::updateAddressing() const
 {
     const cellCellStencilObject& overlap = Stencil::New(*this);
+    // Make sure that stencil is not still in the initial, non-uptodate
+    // state. This gets triggered by e.g. Poisson wall distance that
+    // triggers the stencil even though time has not been updated (and hence
+    // mesh.update() has not been called.
+    if (!stencilIsUpToDate_)
+    {
+        stencilIsUpToDate_ = true;
+        const_cast<cellCellStencilObject&>(overlap).update();
+    }
 
     // The (processor-local part of the) stencil determines the local
     // faces to add to the matrix. tbd: parallel
@@ -227,6 +236,7 @@ Foam::dynamicOversetFvMesh::dynamicOversetFvMesh(const IOobject& io)
 {
     // Load stencil (but do not update)
     (void)Stencil::New(*this, false);
+    stencilIsUpToDate_ = false;
 }
 
 
@@ -390,6 +400,59 @@ bool Foam::dynamicOversetFvMesh::writeObject
         volZoneID.correctBoundaryConditions();
         volZoneID.writeObject(fmt, ver, cmp, valid);
     }
+    if (debug)
+    {
+        const cellCellStencilObject& overlap = Stencil::New(*this);
+        const labelIOList& zoneID = overlap.zoneID();
+        const labelListList& cellStencil = overlap.cellStencil();
+
+        labelList donorZoneID(zoneID);
+        overlap.cellInterpolationMap().distribute(donorZoneID);
+
+        forAll(cellStencil, cellI)
+        {
+            const labelList& stencil = cellStencil[cellI];
+            if (stencil.size())
+            {
+                donorZoneID[cellI] = zoneID[stencil[0]];
+                for (label i = 1; i < stencil.size(); i++)
+                {
+                    if (zoneID[stencil[i]] != donorZoneID[cellI])
+                    {
+                        WarningInFunction << "Mixed donor meshes for cell "
+                            << cellI << " at " << C()[cellI]
+                            << " donors:" << UIndirectList<point>(C(), stencil)
+                            << endl;
+                        donorZoneID[cellI] = -2;
+                    }
+                }
+            }
+        }
+
+        volScalarField volDonorZoneID
+        (
+            IOobject
+            (
+                "donorZoneID",
+                this->time().timeName(),
+                *this,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            *this,
+            dimensionedScalar("minOne", dimless, scalar(-1)),
+            zeroGradientFvPatchScalarField::typeName
+        );
+        forAll(donorZoneID, celli)
+        {
+            volDonorZoneID[celli] = donorZoneID[celli];
+        }
+        //- Do not correctBoundaryConditions since re-interpolates!
+        //volDonorZoneID.correctBoundaryConditions();
+        volDonorZoneID.writeObject(fmt, ver, cmp, valid);
+    }
+
     return ok;
 }
 

src/overset/dynamicOversetFvMesh/dynamicOversetFvMesh.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015-2017 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2015-2018 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -65,6 +65,8 @@ protected:
         //  lduAddressing (true)
         mutable bool active_;
 
+        mutable bool stencilIsUpToDate_;
+
         //- Extended addressing (extended with local interpolation stencils)
         mutable autoPtr<fvMeshPrimitiveLduAddressing> lduPtr_;
 

tutorials/basic/overLaplacianDyMFoam/heatTransfer/0.orig/cellDisplacement

@@ -0,0 +1,59 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1806                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    location    "0.1";
+    object      cellDisplacement;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 0 0 0 0 0];
+
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    free
+    {
+        patchType       overset;
+        type            cellMotion;
+        value           uniform (0 0 0);
+    }
+
+    walls
+    {
+        type            cellMotion;
+        value           uniform (0 0 0);
+    }
+    hole
+    {
+        type            cellMotion;
+        value           uniform (0 0 0);
+    }
+    frontAndBack
+    {
+        type            empty;
+    }
+    left1
+    {
+        type            cellMotion;
+        value           uniform (0 0 0);
+    }
+    right1
+    {
+        type            cellMotion;
+        value           uniform (0 0 0);
+    }
+}
+
+
+// ************************************************************************* //

tutorials/incompressible/overPimpleDyMFoam/simpleRotor/Allrun.pre

@@ -10,9 +10,7 @@ runApplication topoSet
 runApplication subsetMesh box -patch hole -overwrite
 
 # Select cellSets
-rm log.topoSet
+runApplication -s box topoSet
-
-runApplication topoSet
 
 restore0Dir
 

tutorials/mesh/snappyHexMesh/directionalRefinement/Allclean

@@ -0,0 +1,9 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+cleanCase
+
+#------------------------------------------------------------------------------

tutorials/mesh/snappyHexMesh/directionalRefinement/Allrun

@@ -0,0 +1,15 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+runApplication blockMesh
+
+runApplication decomposePar
+
+runParallel snappyHexMesh -overwrite
+
+runApplication reconstructParMesh -constant
+
+#------------------------------------------------------------------------------

tutorials/mesh/snappyHexMesh/directionalRefinement/README.txt

@@ -0,0 +1,9 @@
+Testcase for directional and automatic gap refinement. The geometry
+is two nested boxes
+with a single small pipe between them.
+
+- 1 cell initial mesh
+- consistent normal orientation of surface so
+  specify a 'gapMode' to limit refinement only to
+  gaps on the 'outside' of the surface
+- directional refinement 

tutorials/mesh/snappyHexMesh/directionalRefinement/constant/triSurface/mech_test.obj

@@ -0,0 +1,124 @@
+# Wavefront OBJ file
+# Regions:
+#     0    patch0
+#
+# points    : 40
+# triangles : 76
+#
+v 80 20 -20
+v -80 20 -20
+v 100 0 0
+v -100 0 0
+v 80 20 -280
+v -80 20 -280
+v 80 280 -100
+v -80 280 -100
+v 100 0 -300
+v -100 0 -300
+v 100 300 -120
+v -100 300 -120
+v 80 280 -280
+v -80 280 -280
+v 100 300 -300
+v -100 300 -300
+v 0 445 0
+v 0 445 -20
+v 3.53538 446.464 0
+v -3.53538 446.464 0
+v -3.53538 446.464 -20
+v 3.53538 446.464 -20
+v 5 450 0
+v -5 450 0
+v -5 450 -20
+v 5 450 -20
+v -3.53538 453.536 0
+v 3.53538 453.536 0
+v -3.53538 453.536 -20
+v 3.53538 453.536 -20
+v 0 455 0
+v 0 455 -20
+v 80 580 -20
+v -80 580 -20
+v 80 580 -100
+v -80 580 -100
+v 100 600 0
+v -100 600 0
+v 100 600 -120
+v -100 600 -120
+g patch0
+f 39 37 3
+f 4 3 37
+f 11 39 3
+f 15 11 3
+f 9 15 3
+f 10 9 3
+f 10 3 4
+f 40 37 39
+f 19 4 37
+f 38 37 40
+f 27 37 38
+f 23 19 37
+f 28 23 37
+f 31 28 37
+f 27 31 37
+f 40 39 11
+f 16 11 15
+f 12 11 16
+f 40 11 12
+f 16 15 9
+f 16 9 10
+f 12 10 4
+f 38 12 4
+f 24 38 4
+f 20 24 4
+f 17 20 4
+f 19 17 4
+f 12 16 10
+f 38 40 12
+f 27 38 24
+f 25 24 20
+f 29 27 24
+f 29 24 25
+f 21 20 17
+f 21 25 20
+f 18 17 19
+f 21 17 18
+f 22 19 23
+f 18 19 22
+f 26 23 28
+f 22 23 26
+f 30 28 31
+f 26 28 30
+f 32 31 27
+f 30 31 32
+f 32 27 29
+f 7 5 1
+f 2 1 5
+f 33 7 1
+f 26 33 1
+f 34 1 2
+f 21 1 34
+f 22 26 1
+f 18 22 1
+f 21 18 1
+f 7 13 5
+f 6 5 13
+f 6 2 5
+f 14 13 7
+f 6 13 14
+f 33 35 7
+f 8 7 35
+f 14 7 8
+f 34 35 33
+f 36 35 34
+f 8 35 36
+f 30 34 33
+f 26 30 33
+f 36 34 2
+f 8 36 2
+f 14 8 2
+f 6 14 2
+f 29 25 34
+f 21 34 25
+f 32 29 34
+f 30 32 34

tutorials/mesh/snappyHexMesh/directionalRefinement/system/blockMeshDict

@@ -0,0 +1,102 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scale   1;
+
+vertices
+(
+    (-400 -100  -400)
+    ( 400 -100  -400)
+    ( 400  700  -400)
+    (-400  700  -400)
+    (-400 -100   400)
+    ( 400 -100   400)
+    ( 400  700   400)
+    (-400  700   400)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (1 1 1) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    maxY
+    {
+        type patch;
+        faces
+        (
+            (3 7 6 2)
+        );
+    }
+
+    minX
+    {
+        type patch;
+        faces
+        (
+            (0 4 7 3)
+        );
+    }
+
+    maxX
+    {
+        type patch;
+        faces
+        (
+            (2 6 5 1)
+        );
+    }
+
+    minY
+    {
+        type patch;
+        faces
+        (
+            (1 5 4 0)
+        );
+    }
+
+    ground
+    {
+        type patch;
+        faces
+        (
+            (0 3 2 1)
+        );
+    }
+
+    maxZ
+    {
+        type patch;
+        faces
+        (
+            (4 5 6 7)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/directionalRefinement/system/controlDict

@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus  s                               |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     snappyHexMesh;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         2000;
+
+deltaT          1;
+
+writeControl    timeStep;
+
+writeInterval   100;
+
+purgeWrite      0;
+
+writeFormat     binary;
+
+writePrecision  6;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable true;
+
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/directionalRefinement/system/decomposeParDict

@@ -0,0 +1,23 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 8;
+
+method          scotch;
+
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/directionalRefinement/system/fvSchemes

@@ -0,0 +1,57 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         steadyState;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+
+    div(phi,U)      bounded Gauss upwind;
+    div(phi,T)      bounded Gauss upwind;
+    div(phi,k)      bounded Gauss upwind;
+    div(phi,epsilon) bounded Gauss upwind;
+    div(phi,R)      bounded Gauss upwind;
+    div(R)          Gauss linear;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear limited corrected 0.333;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         limited corrected 0.333;
+}
+
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/directionalRefinement/system/fvSolution

@@ -0,0 +1,69 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p_rgh
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-08;
+        relTol          0.01;
+    }
+
+    "(U|T|k|epsilon)"
+    {
+        solver          PBiCGStab;
+        preconditioner  DILU;
+        tolerance       1e-07;
+        relTol          0.1;
+    }
+}
+
+SIMPLE
+{
+    nNonOrthogonalCorrectors 2;
+    pRefCell        0;
+    pRefValue       0;
+
+    residualControl
+    {
+        p_rgh           1e-2;
+        U               1e-4;
+        T               1e-3;
+
+        // possibly check turbulence fields
+        "(k|epsilon|omega)" 1e-3;
+    }
+}
+
+relaxationFactors
+{
+    fields
+    {
+        p_rgh           0.7;
+    }
+    equations
+    {
+        U               0.2;
+        T               0.5;
+        "(k|epsilon)"   0.7;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/directionalRefinement/system/meshQualityDict

@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      meshQualityDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Include defaults parameters from master dictionary
+#includeEtc "caseDicts/meshQualityDict"
+
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/directionalRefinement/system/snappyHexMeshDict

@@ -0,0 +1,326 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      snappyHexMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Which of the steps to run
+castellatedMesh true;
+snap            false;
+addLayers       false;
+
+
+// Geometry. Definition of all surfaces. All surfaces are of class
+// searchableSurface.
+// Surfaces are used
+// - to specify refinement for any mesh cell intersecting it
+// - to specify refinement for any mesh cell inside/outside/near
+// - to 'snap' the mesh boundary to the surface
+geometry
+{
+    mech_test.obj
+    {
+        type triSurfaceMesh;
+    }
+    all
+    {
+        type searchableBox;
+        min (-1000 -1000 -1000);
+        max (1000 1000 1000);
+    }
+
+    dirRefinement
+    {
+        type searchableBox;
+        min (-75 25 -275);
+        max ( 75 275 -25);
+    }
+};
+
+
+
+// Settings for the castellatedMesh generation.
+castellatedMeshControls
+{
+
+    // Refinement parameters
+    // ~~~~~~~~~~~~~~~~~~~~~
+
+    // If local number of cells is >= maxLocalCells on any processor
+    // switches from from refinement followed by balancing
+    // (current method) to (weighted) balancing before refinement.
+    maxLocalCells 100000;
+
+    // Overall cell limit (approximately). Refinement will stop immediately
+    // upon reaching this number so a refinement level might not complete.
+    // Note that this is the number of cells before removing the part which
+    // is not 'visible' from the keepPoint. The final number of cells might
+    // actually be a lot less.
+    maxGlobalCells 2000000;
+
+    // The surface refinement loop might spend lots of iterations refining just a
+    // few cells. This setting will cause refinement to stop if <= minimumRefine
+    // are selected for refinement. Note: it will at least do one iteration
+    // (unless the number of cells to refine is 0)
+    minRefinementCells 0;
+
+    // Number of buffer layers between different levels.
+    // 1 means normal 2:1 refinement restriction, larger means slower
+    // refinement.
+    nCellsBetweenLevels 1;
+
+
+
+    // Explicit feature edge refinement
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Specifies a level for any cell intersected by its edges.
+    // This is a featureEdgeMesh, read from constant/triSurface for now.
+    features
+    (
+    );
+
+
+
+    // Surface based refinement
+    // ~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Specifies two levels for every surface. The first is the minimum level,
+    // every cell intersecting a surface gets refined up to the minimum level.
+    // The second level is the maximum level. Cells that 'see' multiple
+    // intersections where the intersections make an
+    // angle > resolveFeatureAngle get refined up to the maximum level.
+
+    refinementSurfaces
+    {
+        mech_test.obj
+        {
+            // Surface-wise min and max refinement level
+            level (0 0);
+        }
+    }
+
+    // Resolve sharp angles
+    resolveFeatureAngle 60;
+
+
+    // Region-wise refinement
+    // ~~~~~~~~~~~~~~~~~~~~~~
+
+    // Specifies refinement level for cells in relation to a surface. One of
+    // three modes
+    // - distance. 'levels' specifies per distance to the surface the
+    //   wanted refinement level. The distances need to be specified in
+    //   descending order.
+    // - inside. 'levels' is only one entry and only the level is used. All
+    //   cells inside the surface get refined up to the level. The surface
+    //   needs to be closed for this to be possible.
+    // - outside. Same but cells outside.
+
+    refinementRegions
+    {
+        all
+        {
+            mode        inside;
+
+            // Dummy base level
+            levels      ((10000 0));
+
+            // Optional: automatic gap refinement
+            // If cells
+            // - have level 0..9
+            // - and are in a gap < 3 cell sizes across
+            // - with the gap on the inside ('inside'), outside ('outside')
+            //   or both ('mixed') of the surface
+            // refine them. Ignores 'levels' and 'mode' settings
+            gapLevel    (4 0 10);
+            gapMode     outside;
+        }
+
+
+        dirRefinement
+        {
+            mode        inside;
+
+            // Disable uniform refinement
+            levels      ((10000 0));
+
+            // Optional: directional refinement
+            // (after all other refinement). Directional refinement
+            // for all cells according to 'mode' ('inside' or 'outside';
+            // 'distance' not supported) and within certain range. E.g.
+            //  - for all cells with level 2-5
+            //  - do one split in x direction and three splits in z direction.
+            // Ignores 'levels' and gap* settings. The resulting mesh is
+            // no longer compatible with e.g. dynamic refinement/unrefinement.
+            levelIncrement  (2 5 (1 0 3));
+        }
+    }
+
+
+    // Mesh selection
+    // ~~~~~~~~~~~~~~
+
+    // After refinement patches get added for all refinementSurfaces and
+    // all cells intersecting the surfaces get put into these patches. The
+    // section reachable from the locationInMesh is kept.
+    // NOTE: This point should never be on a face, always inside a cell, even
+    // after refinement.
+    locationInMesh (-100 -5 -300);
+
+
+    // Whether any faceZones (as specified in the refinementSurfaces)
+    // are only on the boundary of corresponding cellZones or also allow
+    // free-standing zone faces. Not used if there are no faceZones.
+    allowFreeStandingZoneFaces false;
+}
+
+
+
+// Settings for the snapping.
+snapControls
+{
+    //- Number of patch smoothing iterations before finding correspondence
+    //  to surface
+    nSmoothPatch 3;
+
+    //- Relative distance for points to be attracted by surface feature point
+    //  or edge. True distance is this factor times local
+    //  maximum edge length.
+    tolerance 2.0;
+
+    //- Number of mesh displacement relaxation iterations.
+    nSolveIter 30;
+
+    //- Maximum number of snapping relaxation iterations. Should stop
+    //  before upon reaching a correct mesh.
+    nRelaxIter 5;
+
+
+    // Feature snapping
+
+        //- Number of feature edge snapping iterations.
+        //  Leave out altogether to disable.
+        nFeatureSnapIter 10;
+
+        //- Detect (geometric) features by sampling the surface (default=false)
+        implicitFeatureSnap true;
+
+        //- Use castellatedMeshControls::features (default = true)
+        explicitFeatureSnap false;
+}
+
+
+
+// Settings for the layer addition.
+addLayersControls
+{
+    // Are the thickness parameters below relative to the undistorted
+    // size of the refined cell outside layer (true) or absolute sizes (false).
+    relativeSizes true;
+
+    // Per final patch (so not geometry!) the layer information
+    layers
+    {
+    }
+
+    // Expansion factor for layer mesh
+    expansionRatio 1.0;
+
+    // Wanted thickness of final added cell layer. If multiple layers
+    // is the thickness of the layer furthest away from the wall.
+    // Relative to undistorted size of cell outside layer.
+    // is the thickness of the layer furthest away from the wall.
+    // See relativeSizes parameter.
+    finalLayerThickness 0.5;
+
+    // Minimum thickness of cell layer. If for any reason layer
+    // cannot be above minThickness do not add layer.
+    // Relative to undistorted size of cell outside layer.
+    // See relativeSizes parameter.
+    minThickness 0.25;
+
+    // If points get not extruded do nGrow layers of connected faces that are
+    // also not grown. This helps convergence of the layer addition process
+    // close to features.
+    // Note: changed(corrected) w.r.t 17x! (didn't do anything in 17x)
+    nGrow 0;
+
+
+    // Advanced settings
+
+    // When not to extrude surface. 0 is flat surface, 90 is when two faces
+    // are perpendicular
+    featureAngle 60;
+
+    // Maximum number of snapping relaxation iterations. Should stop
+    // before upon reaching a correct mesh.
+    nRelaxIter 5;
+
+    // Number of smoothing iterations of surface normals
+    nSmoothSurfaceNormals 1;
+
+    // Number of smoothing iterations of interior mesh movement direction
+    nSmoothNormals 3;
+
+    // Smooth layer thickness over surface patches
+    nSmoothThickness 10;
+
+    // Stop layer growth on highly warped cells
+    maxFaceThicknessRatio 0.5;
+
+    // Reduce layer growth where ratio thickness to medial
+    // distance is large
+    maxThicknessToMedialRatio 0.3;
+
+    // Angle used to pick up medial axis points
+    // Note: changed(corrected) w.r.t 16x! 90 degrees corresponds to 130 in 16x.
+    minMedianAxisAngle 90;
+
+    // Create buffer region for new layer terminations
+    nBufferCellsNoExtrude 0;
+
+
+    // Overall max number of layer addition iterations. The mesher will exit
+    // if it reaches this number of iterations; possibly with an illegal
+    // mesh.
+    nLayerIter 50;
+}
+
+
+
+// Generic mesh quality settings. At any undoable phase these determine
+// where to undo.
+meshQualityControls
+{
+    #include "meshQualityDict"
+
+    // Advanced
+
+    //- Number of error distribution iterations
+    nSmoothScale 4;
+    //- amount to scale back displacement at error points
+    errorReduction 0.75;
+}
+
+
+// Advanced
+
+// Merge tolerance. Is fraction of overall bounding box of initial mesh.
+// Note: the write tolerance needs to be higher than this.
+mergeTolerance 1e-6;
+
+//debugFlags (mesh);
+//writeFlags (scalarLevels);
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/iglooWithFridgesDirectionalRefinement/0/T

@@ -0,0 +1,43 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 265;
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    wall
+    {
+        type            fixedValue;
+        value           uniform 265;
+    }
+
+    twoFridgeFreezers_seal_0
+    {
+        type            fixedValue;
+        value           uniform 303;
+    }
+
+    twoFridgeFreezers_herring_1
+    {
+        $twoFridgeFreezers_seal_0;
+    }
+}
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/iglooWithFridgesDirectionalRefinement/0/U

@@ -0,0 +1,31 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    wall
+    {
+        type            noSlip;
+    }
+}
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/iglooWithFridgesDirectionalRefinement/0/alphat

@@ -0,0 +1,35 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alphat;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    wall
+    {
+        type            alphatJayatillekeWallFunction;
+        Prt             0.85;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/iglooWithFridgesDirectionalRefinement/0/epsilon

@@ -0,0 +1,34 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -3 0 0 0 0];
+
+internalField   uniform 0.01;
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    wall
+    {
+        type            epsilonWallFunction;
+        value           uniform 0.01;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/iglooWithFridgesDirectionalRefinement/0/k

@@ -0,0 +1,34 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0.1;
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    wall
+    {
+        type            kqRWallFunction;
+        value           uniform 0.1;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/iglooWithFridgesDirectionalRefinement/0/nut

@@ -0,0 +1,34 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      nut;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    wall
+    {
+        type            nutkWallFunction;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/iglooWithFridgesDirectionalRefinement/0/p

@@ -0,0 +1,32 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    wall
+    {
+        type            calculated;
+        value           $internalField;
+    }
+}
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/iglooWithFridgesDirectionalRefinement/0/p_rgh

@@ -0,0 +1,33 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    wall
+    {
+        type            fixedFluxPressure;
+        rho             rhok;
+        value           uniform 0;
+    }
+}
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/iglooWithFridgesDirectionalRefinement/Allrun

@@ -0,0 +1,16 @@
+#!/bin/sh
+cd ${0%/*} || exit 1                        # Run from this directory
+. $WM_PROJECT_DIR/bin/tools/RunFunctions    # Tutorial run functions
+
+runApplication blockMesh
+
+## Serial
+#runApplication snappyHexMesh -overwrite
+#runApplication $(getApplication)
+
+## Parallel
+runApplication decomposePar -fileHandler collated
+runParallel snappyHexMesh -overwrite -fileHandler collated
+runApplication reconstructParMesh -constant -fileHandler collated -mergeTol 1e-6
+
+#------------------------------------------------------------------------------

tutorials/mesh/snappyHexMesh/iglooWithFridgesDirectionalRefinement/constant/g

@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       uniformDimensionedVectorField;
+    location    "constant";
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           (0 0 -9.81);
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/iglooWithFridgesDirectionalRefinement/constant/transportProperties

@@ -0,0 +1,35 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      transportProperties;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+transportModel Newtonian;
+
+// Laminar viscosity
+nu              1e-05;
+
+// Thermal expansion coefficient
+beta            3e-03;
+
+// Reference temperature
+TRef            300;
+
+// Laminar Prandtl number
+Pr              0.7;
+
+// Turbulent Prandtl number
+Prt             0.85;
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/iglooWithFridgesDirectionalRefinement/constant/triSurface/fridgeA.eMesh

@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       featureEdgeMesh;
+    object      points;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Points
+(
+(1.99 1.99 0.01)     //0
+(3.01 1.99 0.01)     //1
+(3.01 3.01 0.01)     //2
+(1.99 3.01 0.01)     //3
+
+(1.99 1.99 2.01)     //4
+(3.01 1.99 2.01)     //5
+(3.01 3.01 2.01)     //6
+(1.99 3.01 2.01)     //7
+)
+
+// Edges
+(
+(0 1)
+(1 2)
+(2 3)
+(3 0)
+
+(4 5)
+(5 6)
+(6 7)
+(7 4)
+
+(0 4)
+(1 5)
+(2 6)
+(3 7)
+)
+
+
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/iglooWithFridgesDirectionalRefinement/constant/turbulenceProperties

@@ -0,0 +1,30 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType RAS;
+
+RAS
+{
+    RASModel        kEpsilon;
+
+    turbulence      on;
+
+    printCoeffs     on;
+}
+
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/iglooWithFridgesDirectionalRefinement/system/blockMeshDict

@@ -0,0 +1,102 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scale   1;
+
+vertices
+(
+    (-2.03 -2.0  0)
+    ( 8.03 -2.0  0)
+    ( 8.03  8.0  0)
+    (-2.03  8.0  0)
+    (-2.03 -2.0  5)
+    ( 8.03 -2.0  5)
+    ( 8.03  8.0  5)
+    (-2.03  8.0  5)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (20 20 10) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    maxY
+    {
+        type symmetryPlane;
+        faces
+        (
+            (3 7 6 2)
+        );
+    }
+
+    minX
+    {
+        type symmetryPlane;
+        faces
+        (
+            (0 4 7 3)
+        );
+    }
+
+    maxX
+    {
+        type symmetryPlane;
+        faces
+        (
+            (2 6 5 1)
+        );
+    }
+
+    minY
+    {
+        type symmetryPlane;
+        faces
+        (
+            (1 5 4 0)
+        );
+    }
+
+    ground
+    {
+        type        wall;
+        faces
+        (
+            (0 3 2 1)
+        );
+    }
+
+    maxZ
+    {
+        type symmetryPlane;
+        faces
+        (
+            (4 5 6 7)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/iglooWithFridgesDirectionalRefinement/system/controlDict

@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     buoyantBoussinesqSimpleFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         2000;
+
+deltaT          1;
+
+writeControl    timeStep;
+
+writeInterval   100;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable true;
+
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/iglooWithFridgesDirectionalRefinement/system/decomposeParDict

@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 2;
+
+method          scotch;
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/iglooWithFridgesDirectionalRefinement/system/fvSchemes

@@ -0,0 +1,58 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         steadyState;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+
+    div(phi,U)      bounded Gauss upwind;
+    div(phi,T)      bounded Gauss upwind;
+    div(phi,k)      bounded Gauss upwind;
+    div(phi,epsilon) bounded Gauss upwind;
+    div(phi,R)      bounded Gauss upwind;
+    div(R)          Gauss linear;
+    div((nuEff*dev2(T(grad(U))))) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear limited corrected 0.33;
+    laplacian(diffusivity,cellDisplacement) Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         limited corrected 0.33;
+}
+
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/iglooWithFridgesDirectionalRefinement/system/fvSolution

@@ -0,0 +1,75 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p_rgh
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-08;
+        relTol          0.01;
+    }
+
+    "(U|T|k|epsilon)"
+    {
+        solver          PBiCGStab;
+        preconditioner  DILU;
+        tolerance       1e-07;
+        relTol          0.1;
+    }
+
+    cellDisplacement
+    {
+        $p_rgh;
+        relTol          0;
+    }
+}
+
+SIMPLE
+{
+    nNonOrthogonalCorrectors 2;
+    pRefCell        0;
+    pRefValue       0;
+
+    residualControl
+    {
+        p_rgh           1e-2;
+        U               1e-4;
+        T               1e-3;
+
+        // possibly check turbulence fields
+        "(k|epsilon|omega)" 1e-3;
+    }
+}
+
+relaxationFactors
+{
+    fields
+    {
+        p_rgh           0.7;
+    }
+    equations
+    {
+        U               0.2;
+        T               0.5;
+        "(k|epsilon)"   0.7;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/iglooWithFridgesDirectionalRefinement/system/meshQualityDict

@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      meshQualityDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Include defaults parameters from master dictionary
+#includeEtc "caseDicts/meshQualityDict"
+
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/iglooWithFridgesDirectionalRefinement/system/snappyHexMeshDict

@@ -0,0 +1,420 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      snappyHexMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Which of the steps to run
+castellatedMesh true;
+snap            true;
+addLayers       true;
+
+
+// Geometry. Definition of all surfaces. All surfaces are of class
+// searchableSurface.
+// Surfaces are used
+// - to specify refinement for any mesh cell intersecting it
+// - to specify refinement for any mesh cell inside/outside/near
+// - to 'snap' the mesh boundary to the surface
+geometry
+{
+    igloo
+    {
+        type searchableSphere;
+        centre (3 3 0);
+        radius 4;
+    }
+
+    box1
+    {
+        type searchableBox;
+        min (0 0 0);
+        max (1 1 1);
+    }
+
+    twoFridgeFreezers
+    {
+        type searchableSurfaceCollection;
+
+        mergeSubRegions true;
+
+        seal
+        {
+            surface box1;
+            scale (1.0 1.0 2.1);
+            transform
+            {
+                coordinateSystem
+                {
+                    type    cartesian;
+                    origin  (2 2 0);
+                    coordinateRotation
+                    {
+                        type    axesRotation;
+                        e1      (1 0 0);
+                        e3      (0 0 1);
+                    }
+                }
+            }
+        }
+        herring
+        {
+            surface box1;
+            scale (1.0 1.0 2.1);
+            transform
+            {
+                coordinateSystem
+                {
+                    type    cartesian;
+                    origin  (3.5 3 0);
+                    coordinateRotation
+                    {
+                        type    axesRotation;
+                        e1      (1 0 0);
+                        e3      (0 0 1);
+                    }
+                }
+            }
+        }
+    }
+
+    // For directional refinement
+    dirRefineBox1
+    {
+        type searchableBox;
+        min (-10 -10 -10);
+        max ( 10  10  0.5);
+    }
+    dirRefineBox2
+    {
+        type searchableBox;
+        min (-10 -10 -10);
+        max ( 10  10  0.25);
+    }
+};
+
+
+
+// Settings for the castellatedMesh generation.
+castellatedMeshControls
+{
+
+    // Refinement parameters
+    // ~~~~~~~~~~~~~~~~~~~~~
+
+    // If local number of cells is >= maxLocalCells on any processor
+    // switches from from refinement followed by balancing
+    // (current method) to (weighted) balancing before refinement.
+    maxLocalCells 100000;
+
+    // Overall cell limit (approximately). Refinement will stop immediately
+    // upon reaching this number so a refinement level might not complete.
+    // Note that this is the number of cells before removing the part which
+    // is not 'visible' from the keepPoint. The final number of cells might
+    // actually be a lot less.
+    maxGlobalCells 2000000;
+
+    // The surface refinement loop might spend lots of iterations refining
+    // just a few cells. This setting will cause refinement to stop if <=
+    // minimumRefineare selected for refinement. Note: it will at least do one
+    // iteration (unless the number of cells to refine is 0)
+    minRefinementCells 100;
+
+    // Number of buffer layers between different levels.
+    // 1 means normal 2:1 refinement restriction, larger means slower
+    // refinement.
+    nCellsBetweenLevels 1;
+
+
+
+    // Explicit feature edge refinement
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Specifies a level for any cell intersected by its edges.
+    // This is a featureEdgeMesh, read from constant/triSurface for now.
+    features
+    (
+        {
+            file "fridgeA.eMesh";
+            levels ((0.01 3));
+        }
+    );
+
+
+
+    // Surface based refinement
+    // ~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Specifies two levels for every surface. The first is the minimum level,
+    // every cell intersecting a surface gets refined up to the minimum level.
+    // The second level is the maximum level. Cells that 'see' multiple
+    // intersections where the intersections make an
+    // angle > resolveFeatureAngle get refined up to the maximum level.
+
+    refinementSurfaces
+    {
+        twoFridgeFreezers
+        {
+            // Surface-wise min and max refinement level
+            level (2 2);
+
+            regions
+            {
+                // Region-wise override
+                "cook.*"
+                {
+                    level (3 3);
+                }
+            }
+
+            // Optional specification of patch type (default is wall). No
+            // constraint types (cyclic, symmetry) etc. are allowed.
+            patchInfo
+            {
+                type        wall;
+            }
+        }
+
+        "iglo.*"
+        {
+            // Surface-wise min and max refinement level
+            level (1 1);
+
+            // Optional specification of patch type (default is wall). No
+            // constraint types (cyclic, symmetry) etc. are allowed.
+            patchInfo
+            {
+                type        wall;
+            }
+        }
+    }
+
+    // Resolve sharp angles on fridges
+    resolveFeatureAngle 60;
+
+
+    // Region-wise refinement
+    // ~~~~~~~~~~~~~~~~~~~~~~
+
+    // Specifies refinement level for cells in relation to a surface. One of
+    // three modes
+    // - distance. 'levels' specifies per distance to the surface the
+    //   wanted refinement level. The distances need to be specified in
+    //   descending order.
+    // - inside. 'levels' is only one entry and only the level is used. All
+    //   cells inside the surface get refined up to the level. The surface
+    //   needs to be closed for this to be possible.
+    // - outside. Same but cells outside.
+
+    refinementRegions
+    {
+        dirRefineBox1
+        {
+            mode        inside;
+
+            // Disable uniform refinement
+            levels      ((10000 0));
+
+            // Optional: directional refinement
+            // (after all other refinement). Directional refinement
+            // for all cells according to 'mode' ('inside' or 'outside';
+            // 'distance' not supported) and within certain range. E.g.
+            //  - for all cells with level 0-1
+            //  - do one split in z direction. The resulting mesh is
+            // no longer compatible with e.g. dynamic refinement/unrefinement.
+            levelIncrement  (0 1 (0 0 1));
+        }
+        dirRefineBox2
+        {
+            mode        inside;
+
+            // Disable uniform refinement
+            levels      ((10000 0));
+
+            // Optional: directional refinement
+            // (after all other refinement). Directional refinement
+            // for all cells according to 'mode' ('inside' or 'outside';
+            // 'distance' not supported) and within certain range. E.g.
+            //  - for all cells with level 0-1
+            //  - do two splits in z direction. The resulting mesh is
+            // no longer compatible with e.g. dynamic refinement/unrefinement.
+            levelIncrement  (0 1 (0 0 2));
+        }
+    }
+
+
+    // Mesh selection
+    // ~~~~~~~~~~~~~~
+
+    // After refinement patches get added for all refinementSurfaces and
+    // all cells intersecting the surfaces get put into these patches. The
+    // section reachable from the locationInMesh is kept.
+    // NOTE: This point should never be on a face, always inside a cell, even
+    // after refinement.
+    locationInMesh (3 0.28 0.43);
+
+
+    // Whether any faceZones (as specified in the refinementSurfaces)
+    // are only on the boundary of corresponding cellZones or also allow
+    // free-standing zone faces. Not used if there are no faceZones.
+    allowFreeStandingZoneFaces true;
+}
+
+
+
+// Settings for the snapping.
+snapControls
+{
+    //- Number of patch smoothing iterations before finding correspondence
+    //  to surface
+    nSmoothPatch 3;
+
+    //- Relative distance for points to be attracted by surface feature point
+    //  or edge. True distance is this factor times local
+    //  maximum edge length.
+    tolerance 2.0;
+
+    //- Number of mesh displacement relaxation iterations.
+    nSolveIter 30;
+
+    //- Maximum number of snapping relaxation iterations. Should stop
+    //  before upon reaching a correct mesh.
+    nRelaxIter 5;
+
+
+    // Feature snapping
+
+        //- Number of feature edge snapping iterations.
+        //  Leave out altogether to disable.
+        nFeatureSnapIter 10;
+
+        //- Detect (geometric) features by sampling the surface (default=false)
+        implicitFeatureSnap true;
+
+        //- Use castellatedMeshControls::features (default = true)
+        explicitFeatureSnap false;
+}
+
+
+
+// Settings for the layer addition.
+addLayersControls
+{
+    // Are the thickness parameters below relative to the undistorted
+    // size of the refined cell outside layer (true) or absolute sizes (false).
+    relativeSizes true;
+
+    // Per final patch (so not geometry!) the layer information
+    layers
+    {
+        "two.*"
+        {
+            nSurfaceLayers 3;
+        }
+        ground
+        {
+            nSurfaceLayers 3;
+        }
+        igloo
+        {
+            nSurfaceLayers 1;
+        }
+    }
+
+    // Expansion factor for layer mesh
+    expansionRatio 1.0;
+
+    // Wanted thickness of final added cell layer. If multiple layers
+    // is the thickness of the layer furthest away from the wall.
+    // Relative to undistorted size of cell outside layer.
+    // See relativeSizes parameter.
+    finalLayerThickness 0.5;
+
+    // Minimum thickness of cell layer. If for any reason layer
+    // cannot be above minThickness do not add layer.
+    // Relative to undistorted size of cell outside layer.
+    // See relativeSizes parameter.
+    minThickness 0.05;
+
+    // If points get not extruded do nGrow layers of connected faces that are
+    // also not grown. This helps convergence of the layer addition process
+    // close to features.
+    // Note: changed(corrected) w.r.t 17x! (didn't do anything in 17x)
+    nGrow 0;
+
+
+    // Advanced settings
+
+    // When not to extrude surface. 0 is flat surface, 90 is when two faces
+    // are perpendicular
+    featureAngle 60;
+
+    // Maximum number of snapping relaxation iterations. Should stop
+    // before upon reaching a correct mesh.
+    nRelaxIter 5;
+
+    // Stop layer growth on highly warped cells
+    maxFaceThicknessRatio 0.5;
+
+    // Smooth layer thickness over surface patches
+    nSmoothThickness 10;
+
+
+
+    //- Use displacementMotionSolver to shrink mesh
+    meshShrinker displacementMotionSolver;
+
+    //- Use laplacian for shrinking
+    solver displacementLaplacian;
+
+    displacementLaplacianCoeffs
+    {
+        diffusivity quadratic inverseDistance ("two.*" igloo);
+    }
+
+
+    // Create buffer region for new layer terminations
+    nBufferCellsNoExtrude 0;
+
+    // Overall max number of layer addition iterations. The mesher will exit
+    // if it reaches this number of iterations; possibly with an illegal
+    // mesh.
+    nLayerIter 50;
+}
+
+
+
+// Generic mesh quality settings. At any undoable phase these determine
+// where to undo.
+meshQualityControls
+{
+    #include "meshQualityDict"
+
+    // Advanced
+
+    //- Number of error distribution iterations
+    nSmoothScale 4;
+    //- Amount to scale back displacement at error points
+    errorReduction 0.75;
+}
+
+
+// Advanced
+
+// Merge tolerance. Is fraction of overall bounding box of initial mesh.
+// Note: the write tolerance needs to be higher than this.
+mergeTolerance 1e-6;
+
+writeFlags (scalarLevels);
+
+// ************************************************************************* //