Merge branch 'feature-runtime-selection-geometry' into 'develop'

Feature runtime selection geometry

See merge request Development/openfoam!411

applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C

@@ -704,7 +704,7 @@ int main(int argc, char *argv[])
                 regionName,
                 runTimeExtruded.constant(),
                 runTimeExtruded,
-                IOobject::NO_READ,
+                IOobject::READ_IF_PRESENT,  // Read fv* if present
                 IOobject::AUTO_WRITE,
                 false
             ),

applications/utilities/mesh/manipulation/checkMesh/checkMesh.C

@@ -154,6 +154,7 @@ int main(int argc, char *argv[])
         "cellShapes",
         "cellVolume",
         "cellVolumeRatio",
+        "cellAspectRatio",
         "minTetVolume",
         "minPyrVolume",
         "cellRegion",

applications/utilities/mesh/manipulation/checkMesh/writeFields.C

@@ -7,6 +7,7 @@
 #include "tetPointRef.H"
 #include "regionSplit.H"
 #include "wallDist.H"
+#include "cellAspectRatio.H"
 
 using namespace Foam;
 
@@ -318,6 +319,32 @@ void Foam::writeFields
         aspectRatio.write();
     }
 
+    if (selectedFields.found("cellAspectRatio"))
+    {
+        volScalarField aspectRatio
+        (
+            IOobject
+            (
+                "cellAspectRatio",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            mesh,
+            dimensionedScalar(dimless, Zero),
+            zeroGradientFvPatchScalarField::typeName
+        );
+
+        aspectRatio.ref().field() = cellAspectRatio(mesh);
+
+        aspectRatio.correctBoundaryConditions();
+        Info<< "    Writing approximate aspect ratio to "
+            << aspectRatio.name() << endl;
+        aspectRatio.write();
+    }
+
 
     // cell type
     // ~~~~~~~~~

applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C

@@ -680,7 +680,7 @@ autoPtr<mapPolyMesh> createRegionMesh
             regionName,
             mesh.time().timeName(),
             mesh.time(),
-            IOobject::NO_READ,
+            IOobject::READ_IF_PRESENT,  // read fv* if present
             IOobject::AUTO_WRITE
         ),
         mesh

applications/utilities/parallelProcessing/redistributePar/loadOrCreateMesh.C

@@ -36,6 +36,7 @@ License
 #include "pointSet.H"
 #include "faceSet.H"
 #include "cellSet.H"
+#include "basicFvGeometryScheme.H"
 
 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 
@@ -254,6 +255,19 @@ Foam::autoPtr<Foam::fvMesh> Foam::loadOrCreateMesh
     auto meshPtr = autoPtr<fvMesh>::New(io);
     fvMesh& mesh = *meshPtr;
 
+    // Make sure to use a non-parallel geometry calculation method
+    {
+        tmp<fvGeometryScheme> basicGeometry
+        (
+            fvGeometryScheme::New
+            (
+                mesh,
+                dictionary(),
+                basicFvGeometryScheme::typeName
+            )
+        );
+        mesh.geometry(basicGeometry);
+    }
 
 
     // Sync patches

applications/utilities/parallelProcessing/redistributePar/redistributePar.C

@@ -92,6 +92,7 @@ Usage
 #include "zeroGradientFvPatchFields.H"
 #include "topoSet.H"
 #include "regionProperties.H"
+#include "basicFvGeometryScheme.H"
 
 #include "parFvFieldReconstructor.H"
 #include "parLagrangianRedistributor.H"
@@ -152,6 +153,30 @@ scalar getMergeDistance
 }
 
 
+void setBasicGeometry(fvMesh& mesh)
+{
+    // Set the fvGeometryScheme to basic since it does not require
+    // any parallel communication to construct the geometry. During
+    // redistributePar one might temporarily end up with processors
+    // with zero procBoundaries. Normally procBoundaries trigger geometry
+    // calculation (e.g. send over cellCentres) so on the processors without
+    // procBoundaries this will not happen. The call to the geometry calculation
+    // is not synchronised and this might lead to a hang for geometry schemes
+    // that do require synchronisation
+
+    tmp<fvGeometryScheme> basicGeometry
+    (
+        fvGeometryScheme::New
+        (
+            mesh,
+            dictionary(),
+            basicFvGeometryScheme::typeName
+        )
+    );
+    mesh.geometry(basicGeometry);
+}
+
+
 void printMeshData(const polyMesh& mesh)
 {
     // Collect all data on master
@@ -2675,6 +2700,10 @@ int main(int argc, char *argv[])
                     );
                     fvMesh& mesh = meshPtr();
 
+                    // Use basic geometry calculation to avoid synchronisation
+                    // problems. See comment in routine
+                    setBasicGeometry(mesh);
+
                     // Global matching tolerance
                     const scalar tolDim = getMergeDistance
                     (
@@ -2736,6 +2765,8 @@ int main(int argc, char *argv[])
                 ),
                 true            // read on master only
             );
+            setBasicGeometry(baseMeshPtr());
+
 
             Info<< "Reading local, decomposed mesh" << endl;
             autoPtr<fvMesh> meshPtr = loadOrCreateMesh

src/OpenFOAM/db/IOobject/IOobject.C

@@ -438,6 +438,20 @@ Foam::IOobject::IOobject
 {}
 
 
+Foam::IOobject::IOobject
+(
+    const IOobject& io,
+    readOption ro,
+    writeOption wo
+)
+:
+    IOobject(io)
+{
+    rOpt_ = ro;
+    wOpt_ = wo;
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 const Foam::objectRegistry& Foam::IOobject::db() const

src/OpenFOAM/db/IOobject/IOobject.H

@@ -340,6 +340,14 @@ public:
             const word& name
         );
 
+        //- Copy construct, resetting io options
+        IOobject
+        (
+            const IOobject& io,
+            readOption,
+            writeOption
+        );
+
         //- Clone
         autoPtr<IOobject> clone() const
         {

src/OpenFOAM/include/createMesh.H

@@ -63,9 +63,11 @@ else
                 runTime.timeName(),
                 runTime,
                 Foam::IOobject::MUST_READ
-            )
+            ),
+            false
         )
     );
+    meshPtr().init(true);   // initialise all (lower levels and current)
 }
 
 Foam::fvMesh& mesh = meshPtr();

src/OpenFOAM/include/createNamedMesh.H

@@ -21,5 +21,7 @@ Foam::fvMesh mesh
         runTime.timeName(),
         runTime,
         Foam::IOobject::MUST_READ
-    )
+    ),
+    false
 );
+mesh.init(true);   // initialise all (lower levels and current)

src/OpenFOAM/matrices/solution/solution.C

@@ -155,6 +155,50 @@ Foam::solution::solution
 }
 
 
+Foam::solution::solution
+(
+    const objectRegistry& obr,
+    const fileName& dictName,
+    const dictionary& dict
+)
+:
+    IOdictionary
+    (
+        IOobject
+        (
+            dictName,
+            obr.time().system(),
+            obr,
+            (
+                obr.readOpt() == IOobject::MUST_READ
+             || obr.readOpt() == IOobject::READ_IF_PRESENT
+              ? IOobject::MUST_READ_IF_MODIFIED
+              : obr.readOpt()
+            ),
+            IOobject::NO_WRITE
+        ),
+        dict
+    ),
+    cache_(dictionary::null),
+    caching_(false),
+    fieldRelaxDict_(dictionary::null),
+    eqnRelaxDict_(dictionary::null),
+    fieldRelaxDefault_(0),
+    eqnRelaxDefault_(0),
+    solvers_(dictionary::null)
+{
+    if
+    (
+        readOpt() == IOobject::MUST_READ
+     || readOpt() == IOobject::MUST_READ_IF_MODIFIED
+     || (readOpt() == IOobject::READ_IF_PRESENT && headerOk())
+    )
+    {
+        read(solutionDict());
+    }
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::label Foam::solution::upgradeSolverDict

src/OpenFOAM/matrices/solution/solution.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -101,13 +102,22 @@ public:
 
     // Constructors
 
-        //- Construct for given objectRegistry and dictionary
+        //- Construct for given objectRegistry and dictionary name
         solution
         (
             const objectRegistry& obr,
             const fileName& dictName
         );
 
+        //- Construct for given objectRegistry, dictionary name and (optional)
+        //  content (gets used in case of NO_READ or dictionary cannot be read)
+        solution
+        (
+            const objectRegistry& obr,
+            const fileName& dictName,
+            const dictionary& dict
+        );
+
 
     // Member Functions
 

src/OpenFOAM/meshes/data/data.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -54,6 +55,24 @@ Foam::data::data(const objectRegistry& obr)
 }
 
 
+Foam::data::data(const objectRegistry& obr, const dictionary& dict)
+:
+    IOdictionary
+    (
+        IOobject
+        (
+            "data",
+            obr.time().system(),
+            obr,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        dict
+    ),
+    prevTimeIndex_(0)
+{}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 const Foam::dictionary& Foam::data::solverPerformanceDict() const

src/OpenFOAM/meshes/data/data.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -82,6 +83,9 @@ public:
         //- Construct for objectRegistry
         data(const objectRegistry& obr);
 
+        //- Construct for objectRegistry and initial contents
+        data(const objectRegistry& obr, const dictionary& dict);
+
 
     // Member Functions
 

src/OpenFOAM/meshes/polyMesh/polyMesh.C

@@ -70,20 +70,29 @@ void Foam::polyMesh::calcDirections() const
 
     forAll(boundaryMesh(), patchi)
     {
-        if (boundaryMesh()[patchi].size())
+        const polyPatch& pp = boundaryMesh()[patchi];
+        if (isA<emptyPolyPatch>(pp))
         {
-            if (isA<emptyPolyPatch>(boundaryMesh()[patchi]))
+            // Force calculation of geometric properties, independent of
+            // size. This avoids parallel synchronisation problems.
+            const vectorField::subField fa(pp.faceAreas());
+
+            if (pp.size())
             {
                 nEmptyPatches++;
-                emptyDirVec += sum(cmptMag(boundaryMesh()[patchi].faceAreas()));
+                emptyDirVec += sum(cmptMag(fa));
             }
-            else if (isA<wedgePolyPatch>(boundaryMesh()[patchi]))
-            {
-                const wedgePolyPatch& wpp = refCast<const wedgePolyPatch>
-                (
-                    boundaryMesh()[patchi]
-                );
+        }
+        else if (isA<wedgePolyPatch>(pp))
+        {
+            const wedgePolyPatch& wpp = refCast<const wedgePolyPatch>(pp);
 
+            // Force calculation of geometric properties, independent of
+            // size. This avoids parallel synchronisation problems.
+            (void)wpp.faceNormals();
+
+            if (pp.size())
+            {
                 nWedgePatches++;
                 wedgeDirVec += cmptMag(wpp.centreNormal());
             }
@@ -161,7 +170,7 @@ Foam::autoPtr<Foam::labelIOList> Foam::polyMesh::readTetBasePtIs() const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::polyMesh::polyMesh(const IOobject& io)
+Foam::polyMesh::polyMesh(const IOobject& io, const bool doInit)
 :
     objectRegistry(io),
     primitiveMesh(),
@@ -328,12 +337,6 @@ Foam::polyMesh::polyMesh(const IOobject& io)
         neighbour_.write();
     }
 
-    // Calculate topology for the patches (processor-processor comms etc.)
-    boundary_.updateMesh();
-
-    // Calculate the geometry for the patches (transformation tensors etc.)
-    boundary_.calcGeometry();
-
     // Warn if global empty mesh
     if (returnReduce(boundary_.empty(), orOp<bool>()))
     {
@@ -352,8 +355,30 @@ Foam::polyMesh::polyMesh(const IOobject& io)
         }
     }
 
+    if (doInit)
+    {
+        polyMesh::init(false);  // do not init lower levels
+    }
+}
+
+
+bool Foam::polyMesh::init(const bool doInit)
+{
+    if (doInit)
+    {
+        primitiveMesh::init(doInit);
+    }
+
+    // Calculate topology for the patches (processor-processor comms etc.)
+    boundary_.updateMesh();
+
+    // Calculate the geometry for the patches (transformation tensors etc.)
+    boundary_.calcGeometry();
+
     // Initialise demand-driven data
     calcDirections();
+
+    return false;
 }
 
 
@@ -377,7 +402,7 @@ Foam::polyMesh::polyMesh
             instance(),
             meshSubDir,
             *this,
-            io.readOpt(),
+            IOobject::NO_READ,  //io.readOpt(),
             io.writeOpt()
         ),
         std::move(points)
@@ -390,7 +415,7 @@ Foam::polyMesh::polyMesh
             instance(),
             meshSubDir,
             *this,
-            io.readOpt(),
+            IOobject::NO_READ,  //io.readOpt(),
             io.writeOpt()
         ),
         std::move(faces)
@@ -403,7 +428,7 @@ Foam::polyMesh::polyMesh
             instance(),
             meshSubDir,
             *this,
-            io.readOpt(),
+            IOobject::NO_READ,  //io.readOpt(),
             io.writeOpt()
         ),
         std::move(owner)
@@ -416,7 +441,7 @@ Foam::polyMesh::polyMesh
             instance(),
             meshSubDir,
             *this,
-            io.readOpt(),
+            IOobject::NO_READ,  //io.readOpt(),
             io.writeOpt()
         ),
         std::move(neighbour)
@@ -430,7 +455,7 @@ Foam::polyMesh::polyMesh
             instance(),
             meshSubDir,
             *this,
-            io.readOpt(),
+            IOobject::NO_READ,  //io.readOpt(),
             io.writeOpt()
         ),
         *this,
@@ -440,7 +465,7 @@ Foam::polyMesh::polyMesh
     comm_(UPstream::worldComm),
     geometricD_(Zero),
     solutionD_(Zero),
-    tetBasePtIsPtr_(readTetBasePtIs()),
+    tetBasePtIsPtr_(nullptr),
     cellTreePtr_(nullptr),
     pointZones_
     (
@@ -450,7 +475,7 @@ Foam::polyMesh::polyMesh
             instance(),
             meshSubDir,
             *this,
-            io.readOpt(),
+            IOobject::NO_READ,  //io.readOpt(),
             IOobject::NO_WRITE
         ),
         *this,
@@ -464,7 +489,7 @@ Foam::polyMesh::polyMesh
             instance(),
             meshSubDir,
             *this,
-            io.readOpt(),
+            IOobject::NO_READ,  //io.readOpt(),
             IOobject::NO_WRITE
         ),
         *this,
@@ -478,7 +503,7 @@ Foam::polyMesh::polyMesh
             instance(),
             meshSubDir,
             *this,
-            io.readOpt(),
+            IOobject::NO_READ,  //io.readOpt(),
             IOobject::NO_WRITE
         ),
         *this,
@@ -591,7 +616,7 @@ Foam::polyMesh::polyMesh
     comm_(UPstream::worldComm),
     geometricD_(Zero),
     solutionD_(Zero),
-    tetBasePtIsPtr_(readTetBasePtIs()),
+    tetBasePtIsPtr_(nullptr),
     cellTreePtr_(nullptr),
     pointZones_
     (

src/OpenFOAM/meshes/polyMesh/polyMesh.H

@@ -319,11 +319,11 @@ public:
     // Constructors
 
         //- Read construct from IOobject
-        explicit polyMesh(const IOobject& io);
+        polyMesh(const IOobject& io, const bool doInit = true);
 
         //- Construct from IOobject or as zero-sized mesh
         //  Boundary is added using addPatches() member function
-        polyMesh(const IOobject& io, const zero, bool syncPar=true);
+        polyMesh(const IOobject& io, const zero, const bool syncPar=true);
 
         //- Construct from IOobject or from components.
         //  Boundary is added using addPatches() member function
@@ -593,6 +593,9 @@ public:
                 const List<cellZone*>& cz
             );
 
+            //- Initialise all non-demand-driven data
+            virtual bool init(const bool doInit);
+
             //- Update the mesh based on the mesh files saved in
             //  time directories
             virtual readUpdateState readUpdate();

src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshCheck.C

@@ -753,13 +753,26 @@ bool Foam::polyMesh::checkMeshMotion
     vectorField fCtrs(nFaces());
     vectorField fAreas(nFaces());
 
-    makeFaceCentresAndAreas(newPoints, fCtrs, fAreas);
+    primitiveMeshTools::makeFaceCentresAndAreas
+    (
+        *this,
+        newPoints,
+        fCtrs,
+        fAreas
+    );
 
     // Check cell volumes and calculate new cell centres
     vectorField cellCtrs(nCells());
     scalarField cellVols(nCells());
 
-    makeCellCentresAndVols(fCtrs, fAreas, cellCtrs, cellVols);
+    primitiveMeshTools::makeCellCentresAndVols
+    (
+        *this,
+        fCtrs,
+        fAreas,
+        cellCtrs,
+        cellVols
+    );
 
     // Check cell volumes
     bool error = checkCellVolumes

src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C

@@ -281,6 +281,44 @@ void Foam::primitiveMesh::reset
 }
 
 
+void Foam::primitiveMesh::resetGeometry
+(
+    pointField&& faceCentres,
+    pointField&& faceAreas,
+    pointField&& cellCentres,
+    scalarField&& cellVolumes
+)
+{
+    if
+    (
+        faceCentres.size() != nFaces_
+     || faceAreas.size() != nFaces_
+     || cellCentres.size() != nCells_
+     || cellVolumes.size() != nCells_
+    )
+    {
+        FatalErrorInFunction
+            << "Wrong sizes of passed in face/cell data"
+            << abort(FatalError);
+    }
+
+    // Remove old geometry
+    clearGeom();
+
+    faceCentresPtr_ = new pointField(std::move(faceCentres));
+    faceAreasPtr_ = new pointField(std::move(faceAreas));
+    cellCentresPtr_ = new pointField(std::move(cellCentres));
+    cellVolumesPtr_ = new scalarField(std::move(cellVolumes));
+
+    if (debug)
+    {
+        Pout<< "primitiveMesh::resetGeometry : geometry reset for"
+            << " nFaces:" << faceCentresPtr_->size()
+            << " nCells:" << cellCentresPtr_->size() << endl;
+    }
+}
+
+
 Foam::tmp<Foam::scalarField> Foam::primitiveMesh::movePoints
 (
     const pointField& newPoints,
@@ -324,4 +362,26 @@ const Foam::cellShapeList& Foam::primitiveMesh::cellShapes() const
 }
 
 
+
+void Foam::primitiveMesh::updateGeom()
+{
+    if (!faceCentresPtr_)
+    {
+        calcFaceCentresAndAreas();
+    }
+    if (!faceAreasPtr_)
+    {
+        calcFaceCentresAndAreas();
+    }
+    if (!cellCentresPtr_)
+    {
+        calcCellCentresAndVols();
+    }
+    if (!cellVolumesPtr_)
+    {
+        calcCellCentresAndVols();
+    }
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H

@@ -255,22 +255,9 @@ protected:
 
             //- Calculate face centres and areas
             void calcFaceCentresAndAreas() const;
-            void makeFaceCentresAndAreas
-            (
-                const pointField& p,
-                vectorField& fCtrs,
-                vectorField& fAreas
-            ) const;
 
             //- Calculate cell centres and volumes
             void calcCellCentresAndVols() const;
-            void makeCellCentresAndVols
-            (
-                const vectorField& fCtrs,
-                const vectorField& fAreas,
-                vectorField& cellCtrs,
-                scalarField& cellVols
-            ) const;
 
             //- Calculate edge vectors
             void calcEdgeVectors() const;
@@ -476,6 +463,21 @@ public:
             cellList& cells
         );
 
+        //- Reset the local geometry
+        void resetGeometry
+        (
+            pointField&& faceCentres,
+            pointField&& faceAreas,
+            pointField&& cellCentres,
+            scalarField&& cellVolumes
+        );
+
+        //- Initialise all non-demand-driven data
+        virtual bool init(const bool doInit)
+        {
+            return false;
+        }
+
 
         // Access
 
@@ -891,6 +893,8 @@ public:
 
             const labelList& cellEdges(const label celli) const;
 
+            //- Update all geometric data
+            virtual void updateGeom();
 
             //- Clear geometry
             void clearGeom();

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCentresAndVols.C

@@ -30,7 +30,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "primitiveMesh.H"
-#include "PrecisionAdaptor.H"
+#include "primitiveMeshTools.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -61,7 +61,14 @@ void Foam::primitiveMesh::calcCellCentresAndVols() const
     scalarField& cellVols = *cellVolumesPtr_;
 
     // Make centres and volumes
-    makeCellCentresAndVols(faceCentres(), faceAreas(), cellCtrs, cellVols);
+    primitiveMeshTools::makeCellCentresAndVols
+    (
+        *this,
+        faceCentres(),
+        faceAreas(),
+        cellCtrs,
+        cellVols
+    );
 
     if (debug)
     {
@@ -72,111 +79,14 @@ void Foam::primitiveMesh::calcCellCentresAndVols() const
 }
 
 
-void Foam::primitiveMesh::makeCellCentresAndVols
-(
-    const vectorField& fCtrs,
-    const vectorField& fAreas,
-    vectorField& cellCtrs_s,
-    scalarField& cellVols_s
-) const
-{
-    typedef Vector<solveScalar> solveVector;
-
-    PrecisionAdaptor<solveVector, vector> tcellCtrs(cellCtrs_s, false);
-    Field<solveVector>& cellCtrs = tcellCtrs.ref();
-    PrecisionAdaptor<solveScalar, scalar> tcellVols(cellVols_s, false);
-    Field<solveScalar>& cellVols = tcellVols.ref();
-
-    cellCtrs = Zero;
-    cellVols = 0.0;
-
-    const labelList& own = faceOwner();
-    const labelList& nei = faceNeighbour();
-
-    // first estimate the approximate cell centre as the average of
-    // face centres
-
-    Field<solveVector> cEst(nCells(), Zero);
-    labelField nCellFaces(nCells(), Zero);
-
-    forAll(own, facei)
-    {
-        cEst[own[facei]] += solveVector(fCtrs[facei]);
-        ++nCellFaces[own[facei]];
-    }
-
-    forAll(nei, facei)
-    {
-        cEst[nei[facei]] += solveVector(fCtrs[facei]);
-        ++nCellFaces[nei[facei]];
-    }
-
-    forAll(cEst, celli)
-    {
-        cEst[celli] /= nCellFaces[celli];
-    }
-
-    forAll(own, facei)
-    {
-        const solveVector fc(fCtrs[facei]);
-        const solveVector fA(fAreas[facei]);
-
-        // Calculate 3*face-pyramid volume
-        solveScalar pyr3Vol =
-            fA & (fc - cEst[own[facei]]);
-
-        // Calculate face-pyramid centre
-        solveVector pc = (3.0/4.0)*fc + (1.0/4.0)*cEst[own[facei]];
-
-        // Accumulate volume-weighted face-pyramid centre
-        cellCtrs[own[facei]] += pyr3Vol*pc;
-
-        // Accumulate face-pyramid volume
-        cellVols[own[facei]] += pyr3Vol;
-    }
-
-    forAll(nei, facei)
-    {
-        const solveVector fc(fCtrs[facei]);
-        const solveVector fA(fAreas[facei]);
-
-        // Calculate 3*face-pyramid volume
-        solveScalar pyr3Vol =
-            fA & (cEst[nei[facei]] - fc);
-
-        // Calculate face-pyramid centre
-        solveVector pc = (3.0/4.0)*fc + (1.0/4.0)*cEst[nei[facei]];
-
-        // Accumulate volume-weighted face-pyramid centre
-        cellCtrs[nei[facei]] += pyr3Vol*pc;
-
-        // Accumulate face-pyramid volume
-        cellVols[nei[facei]] += pyr3Vol;
-    }
-
-    forAll(cellCtrs, celli)
-    {
-        if (mag(cellVols[celli]) > VSMALL)
-        {
-            cellCtrs[celli] /= cellVols[celli];
-        }
-        else
-        {
-            cellCtrs[celli] = cEst[celli];
-        }
-    }
-
-    cellVols *= (1.0/3.0);
-}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 const Foam::vectorField& Foam::primitiveMesh::cellCentres() const
 {
     if (!cellCentresPtr_)
     {
-        calcCellCentresAndVols();
+        //calcCellCentresAndVols();
+        const_cast<primitiveMesh&>(*this).updateGeom();
     }
 
     return *cellCentresPtr_;
@@ -187,7 +97,8 @@ const Foam::scalarField& Foam::primitiveMesh::cellVolumes() const
 {
     if (!cellVolumesPtr_)
     {
-        calcCellCentresAndVols();
+        //calcCellCentresAndVols();
+        const_cast<primitiveMesh&>(*this).updateGeom();
     }
 
     return *cellVolumesPtr_;

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C

@@ -27,11 +27,182 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "primitiveMeshTools.H"
+#include "primitiveMesh.H"
 #include "syncTools.H"
 #include "pyramidPointFaceRef.H"
+#include "PrecisionAdaptor.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+void Foam::primitiveMeshTools::makeFaceCentresAndAreas
+(
+    const primitiveMesh& mesh,
+    const pointField& p,
+    vectorField& fCtrs,
+    vectorField& fAreas
+)
+{
+    const faceList& fs = mesh.faces();
+
+    forAll(fs, facei)
+    {
+        const labelList& f = fs[facei];
+        const label nPoints = f.size();
+
+        // If the face is a triangle, do a direct calculation for efficiency
+        // and to avoid round-off error-related problems
+        if (nPoints == 3)
+        {
+            fCtrs[facei] = (1.0/3.0)*(p[f[0]] + p[f[1]] + p[f[2]]);
+            fAreas[facei] = 0.5*((p[f[1]] - p[f[0]])^(p[f[2]] - p[f[0]]));
+        }
+        else
+        {
+            typedef Vector<solveScalar> solveVector;
+
+            solveVector sumN = Zero;
+            solveScalar sumA = 0.0;
+            solveVector sumAc = Zero;
+
+            solveVector fCentre = p[f[0]];
+            for (label pi = 1; pi < nPoints; pi++)
+            {
+                fCentre += solveVector(p[f[pi]]);
+            }
+
+            fCentre /= nPoints;
+
+            for (label pi = 0; pi < nPoints; pi++)
+            {
+                const label nextPi(pi == nPoints-1 ? 0 : pi+1);
+                const solveVector nextPoint(p[f[nextPi]]);
+                const solveVector thisPoint(p[f[pi]]);
+
+                solveVector c = thisPoint + nextPoint + fCentre;
+                solveVector n = (nextPoint - thisPoint)^(fCentre - thisPoint);
+                solveScalar a = mag(n);
+                sumN += n;
+                sumA += a;
+                sumAc += a*c;
+            }
+
+            // This is to deal with zero-area faces. Mark very small faces
+            // to be detected in e.g., processorPolyPatch.
+            if (sumA < ROOTVSMALL)
+            {
+                fCtrs[facei] = fCentre;
+                fAreas[facei] = Zero;
+            }
+            else
+            {
+                fCtrs[facei] = (1.0/3.0)*sumAc/sumA;
+                fAreas[facei] = 0.5*sumN;
+            }
+        }
+    }
+}
+
+
+void Foam::primitiveMeshTools::makeCellCentresAndVols
+(
+    const primitiveMesh& mesh,
+    const vectorField& fCtrs,
+    const vectorField& fAreas,
+    vectorField& cellCtrs_s,
+    scalarField& cellVols_s
+)
+{
+    typedef Vector<solveScalar> solveVector;
+
+    PrecisionAdaptor<solveVector, vector> tcellCtrs(cellCtrs_s);
+    Field<solveVector>& cellCtrs = tcellCtrs.ref();
+    PrecisionAdaptor<solveScalar, scalar> tcellVols(cellVols_s);
+    Field<solveScalar>& cellVols = tcellVols.ref();
+
+    // Clear the fields for accumulation
+    cellCtrs = Zero;
+    cellVols = 0.0;
+
+    const labelList& own = mesh.faceOwner();
+    const labelList& nei = mesh.faceNeighbour();
+
+    // first estimate the approximate cell centre as the average of
+    // face centres
+
+    Field<solveVector> cEst(mesh.nCells(), Zero);
+    labelField nCellFaces(mesh.nCells(), Zero);
+
+    forAll(own, facei)
+    {
+        cEst[own[facei]] += solveVector(fCtrs[facei]);
+        ++nCellFaces[own[facei]];
+    }
+
+    forAll(nei, facei)
+    {
+        cEst[nei[facei]] += solveVector(fCtrs[facei]);
+        ++nCellFaces[nei[facei]];
+    }
+
+    forAll(cEst, celli)
+    {
+        cEst[celli] /= nCellFaces[celli];
+    }
+
+    forAll(own, facei)
+    {
+        const solveVector fc(fCtrs[facei]);
+        const solveVector fA(fAreas[facei]);
+
+        // Calculate 3*face-pyramid volume
+        solveScalar pyr3Vol =
+            fA & (fc - cEst[own[facei]]);
+
+        // Calculate face-pyramid centre
+        solveVector pc = (3.0/4.0)*fc + (1.0/4.0)*cEst[own[facei]];
+
+        // Accumulate volume-weighted face-pyramid centre
+        cellCtrs[own[facei]] += pyr3Vol*pc;
+
+        // Accumulate face-pyramid volume
+        cellVols[own[facei]] += pyr3Vol;
+    }
+
+    forAll(nei, facei)
+    {
+        const solveVector fc(fCtrs[facei]);
+        const solveVector fA(fAreas[facei]);
+
+        // Calculate 3*face-pyramid volume
+        solveScalar pyr3Vol =
+            fA & (cEst[nei[facei]] - fc);
+
+        // Calculate face-pyramid centre
+        solveVector pc = (3.0/4.0)*fc + (1.0/4.0)*cEst[nei[facei]];
+
+        // Accumulate volume-weighted face-pyramid centre
+        cellCtrs[nei[facei]] += pyr3Vol*pc;
+
+        // Accumulate face-pyramid volume
+        cellVols[nei[facei]] += pyr3Vol;
+    }
+
+    forAll(cellCtrs, celli)
+    {
+        if (mag(cellVols[celli]) > VSMALL)
+        {
+            cellCtrs[celli] /= cellVols[celli];
+        }
+        else
+        {
+            cellCtrs[celli] = cEst[celli];
+        }
+    }
+
+    cellVols *= (1.0/3.0);
+}
+
+
 Foam::scalar Foam::primitiveMeshTools::faceSkewness
 (
     const primitiveMesh& mesh,

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.H

@@ -38,22 +38,45 @@ SourceFiles
 #ifndef primitiveMeshTools_H
 #define primitiveMeshTools_H
 
-#include "primitiveMesh.H"
 #include "bitSet.H"
+#include "primitiveFields.H"
+#include "pointField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
+// Forward declarations
+class primitiveMesh;
+
 /*---------------------------------------------------------------------------*\
-                        Namespace primitiveMeshTools Declaration
+                   Namespace primitiveMeshTools Declaration
 \*---------------------------------------------------------------------------*/
 
 class primitiveMeshTools
 {
 public:
 
+    //- Calculate face centres and areas
+    static void makeFaceCentresAndAreas
+    (
+        const primitiveMesh& mesh,
+        const pointField& p,
+        vectorField& fCtrs,
+        vectorField& fAreas
+    );
+
+    //- Calculate cell centres and volumes from face properties
+    static void makeCellCentresAndVols
+    (
+        const primitiveMesh& mesh,
+        const vectorField& fCtrs,
+        const vectorField& fAreas,
+        vectorField& cellCtrs,
+        scalarField& cellVols
+    );
+
     //- Generate non-orthogonality field (internal faces only)
     static tmp<scalarField> faceOrthogonality
     (

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C

@@ -33,6 +33,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "primitiveMesh.H"
+#include "primitiveMeshTools.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -60,7 +61,7 @@ void Foam::primitiveMesh::calcFaceCentresAndAreas() const
     faceAreasPtr_ = new vectorField(nFaces());
     vectorField& fAreas = *faceAreasPtr_;
 
-    makeFaceCentresAndAreas(points(), fCtrs, fAreas);
+    primitiveMeshTools::makeFaceCentresAndAreas(*this, points(), fCtrs, fAreas);
 
     if (debug)
     {
@@ -71,81 +72,14 @@ void Foam::primitiveMesh::calcFaceCentresAndAreas() const
 }
 
 
-void Foam::primitiveMesh::makeFaceCentresAndAreas
-(
-    const pointField& p,
-    vectorField& fCtrs,
-    vectorField& fAreas
-) const
-{
-    const faceList& fs = faces();
-
-    forAll(fs, facei)
-    {
-        const labelList& f = fs[facei];
-        const label nPoints = f.size();
-
-        // If the face is a triangle, do a direct calculation for efficiency
-        // and to avoid round-off error-related problems
-        if (nPoints == 3)
-        {
-            fCtrs[facei] = (1.0/3.0)*(p[f[0]] + p[f[1]] + p[f[2]]);
-            fAreas[facei] = 0.5*((p[f[1]] - p[f[0]])^(p[f[2]] - p[f[0]]));
-        }
-        else
-        {
-            typedef Vector<solveScalar> solveVector;
-
-            solveVector sumN = Zero;
-            solveScalar sumA = 0.0;
-            solveVector sumAc = Zero;
-
-            solveVector fCentre = p[f[0]];
-            for (label pi = 1; pi < nPoints; pi++)
-            {
-                fCentre += solveVector(p[f[pi]]);
-            }
-
-            fCentre /= nPoints;
-
-            for (label pi = 0; pi < nPoints; pi++)
-            {
-                const label nextPi(pi == nPoints-1 ? 0 : pi+1);
-                const solveVector nextPoint(p[f[nextPi]]);
-                const solveVector thisPoint(p[f[pi]]);
-
-                solveVector c = thisPoint + nextPoint + fCentre;
-                solveVector n = (nextPoint - thisPoint)^(fCentre - thisPoint);
-                solveScalar a = mag(n);
-                sumN += n;
-                sumA += a;
-                sumAc += a*c;
-            }
-
-            // This is to deal with zero-area faces. Mark very small faces
-            // to be detected in e.g., processorPolyPatch.
-            if (sumA < ROOTVSMALL)
-            {
-                fCtrs[facei] = fCentre;
-                fAreas[facei] = Zero;
-            }
-            else
-            {
-                fCtrs[facei] = (1.0/3.0)*sumAc/sumA;
-                fAreas[facei] = 0.5*sumN;
-            }
-        }
-    }
-}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 const Foam::vectorField& Foam::primitiveMesh::faceCentres() const
 {
     if (!faceCentresPtr_)
     {
-        calcFaceCentresAndAreas();
+        //calcFaceCentresAndAreas();
+        const_cast<primitiveMesh&>(*this).updateGeom();
     }
 
     return *faceCentresPtr_;
@@ -156,7 +90,8 @@ const Foam::vectorField& Foam::primitiveMesh::faceAreas() const
 {
     if (!faceAreasPtr_)
     {
-        calcFaceCentresAndAreas();
+        //calcFaceCentresAndAreas();
+        const_cast<primitiveMesh&>(*this).updateGeom();
     }
 
     return *faceAreasPtr_;

src/dynamicFvMesh/dynamicFvMesh/dynamicFvMesh.H

@@ -104,7 +104,7 @@ public:
     // Constructors
 
         //- Construct from an IOobject
-        explicit dynamicFvMesh(const IOobject& io);
+        explicit dynamicFvMesh(const IOobject& io); //, const bool doInit=true);
 
         //- Construct from components without boundary.
         //  Boundary is added using addFvPatches() member function
@@ -139,7 +139,11 @@ public:
         //- Select, construct and return the dynamicFvMesh
         //  If the constant/dynamicMeshDict does not exist
         //  a staticFvMesh is returned
-        static autoPtr<dynamicFvMesh> New(const IOobject& io);
+        static autoPtr<dynamicFvMesh> New
+        (
+            const IOobject& io  //,
+            //const bool doInit=true
+        );
 
 
         //- Select, construct and return the dynamicFvMesh
@@ -148,7 +152,8 @@ public:
         static autoPtr<dynamicFvMesh> New
         (
             const argList& args,
-            const Time& runTime
+            const Time& runTime //,
+            //const bool doInit=true
         );
 
 
@@ -158,6 +163,9 @@ public:
 
     // Member Functions
 
+        //- Initialise all non-demand-driven data
+        //virtual bool init(const bool doInit);
+
         //- Is mesh dynamic
         virtual bool dynamic() const
         {

src/dynamicFvMesh/include/createDynamicFvMesh.H

@@ -1,6 +1,9 @@
 Info<< "Create mesh for time = "
     << runTime.timeName() << nl << endl;
 
-autoPtr<dynamicFvMesh> meshPtr(dynamicFvMesh::New(args, runTime));
+autoPtr<dynamicFvMesh> meshPtr(dynamicFvMesh::New(args, runTime));//, false));
 
 dynamicFvMesh& mesh = meshPtr();
+
+// initialise all (lower levels and current)
+mesh.init(true);

src/dynamicMesh/fvMeshSubset/fvMeshSubset.C

@@ -268,7 +268,7 @@ void Foam::fvMeshSubset::removeCellsImpl
             baseMesh().name(),
             baseMesh().time().timeName(),
             baseMesh().time(),
-            IOobject::NO_READ,
+            IOobject::READ_IF_PRESENT,  // read fv* if present
             IOobject::NO_WRITE
         ),
         baseMesh(),
@@ -980,9 +980,10 @@ void Foam::fvMeshSubset::setCellSubset
             baseMesh().name(),
             baseMesh().time().timeName(),
             baseMesh().time(),
-            IOobject::NO_READ,
+            IOobject::NO_READ,      // do not read any dictionaries
             IOobject::NO_WRITE
         ),
+        baseMesh(),                 // get dictionaries from base mesh
         std::move(newPoints),
         std::move(newFaces),
         std::move(newCells),

src/dynamicMesh/polyMeshFilter/polyMeshFilter.C

@@ -79,7 +79,7 @@ Foam::autoPtr<Foam::fvMesh> Foam::polyMeshFilter::copyMesh(const fvMesh& mesh)
             mesh.name(),
             mesh.polyMesh::instance(),
             mesh.time(),
-            IOobject::NO_READ,
+            IOobject::READ_IF_PRESENT,  // read fv* if present
             IOobject::NO_WRITE,
             false
         ),

src/finiteVolume/Make/files

@@ -1,6 +1,17 @@
 fvMesh/fvMeshGeometry.C
 fvMesh/fvMesh.C
 
+fvGeometryScheme = fvMesh/fvGeometryScheme
+$(fvGeometryScheme)/fvGeometryScheme/fvGeometryScheme.C
+$(fvGeometryScheme)/basic/basicFvGeometryScheme.C
+$(fvGeometryScheme)/highAspectRatio/highAspectRatioFvGeometryScheme.C
+$(fvGeometryScheme)/highAspectRatio/cellAspectRatio.C
+$(fvGeometryScheme)/averageNeighbour/averageNeighbourFvGeometryScheme.C
+$(fvGeometryScheme)/stabilised/stabilisedFvGeometryScheme.C
+
+surfaceInterpolation = interpolation/surfaceInterpolation
+$(surfaceInterpolation)/surfaceInterpolation/surfaceInterpolation.C
+
 fvMesh/singleCellFvMesh/singleCellFvMesh.C
 
 fvMesh/simplifiedFvMesh/simplifiedFvMesh/simplifiedFvMesh.C
@@ -53,7 +64,6 @@ fvMeshMapper = fvMesh/fvMeshMapper
 $(fvMeshMapper)/fvPatchMapper.C
 $(fvMeshMapper)/fvSurfaceMapper.C
 
-
 extendedStencil = fvMesh/extendedStencil
 
 cellToCell = $(extendedStencil)/cellToCell
@@ -316,8 +326,6 @@ volPointInterpolation = interpolation/volPointInterpolation
 $(volPointInterpolation)/volPointInterpolation.C
 $(volPointInterpolation)/pointConstraints.C
 
-surfaceInterpolation = interpolation/surfaceInterpolation
-$(surfaceInterpolation)/surfaceInterpolation/surfaceInterpolation.C
 $(surfaceInterpolation)/surfaceInterpolationScheme/surfaceInterpolationSchemes.C
 
 $(surfaceInterpolation)/blendedSchemeBase/blendedSchemeBaseName.C

src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -327,6 +327,140 @@ Foam::fvSchemes::fvSchemes(const objectRegistry& obr)
 }
 
 
+Foam::fvSchemes::fvSchemes(const objectRegistry& obr, const dictionary& dict)
+:
+    IOdictionary
+    (
+        IOobject
+        (
+            "fvSchemes",
+            obr.time().system(),
+            obr,
+            (
+                obr.readOpt() == IOobject::MUST_READ
+             || obr.readOpt() == IOobject::READ_IF_PRESENT
+              ? IOobject::MUST_READ_IF_MODIFIED
+              : obr.readOpt()
+            ),
+            IOobject::NO_WRITE
+        ),
+        dict
+    ),
+    ddtSchemes_
+    (
+        ITstream
+        (
+            objectPath() + ".ddtSchemes",
+            tokenList()
+        )()
+    ),
+    defaultDdtScheme_
+    (
+        ddtSchemes_.name() + ".default",
+        tokenList()
+    ),
+    d2dt2Schemes_
+    (
+        ITstream
+        (
+            objectPath() + ".d2dt2Schemes",
+            tokenList()
+        )()
+    ),
+    defaultD2dt2Scheme_
+    (
+        d2dt2Schemes_.name() + ".default",
+        tokenList()
+    ),
+    interpolationSchemes_
+    (
+        ITstream
+        (
+            objectPath() + ".interpolationSchemes",
+            tokenList()
+        )()
+    ),
+    defaultInterpolationScheme_
+    (
+        interpolationSchemes_.name() + ".default",
+        tokenList()
+    ),
+    divSchemes_
+    (
+        ITstream
+        (
+            objectPath() + ".divSchemes",
+            tokenList()
+        )()
+    ),
+    defaultDivScheme_
+    (
+        divSchemes_.name() + ".default",
+        tokenList()
+    ),
+    gradSchemes_
+    (
+        ITstream
+        (
+            objectPath() + ".gradSchemes",
+            tokenList()
+        )()
+    ),
+    defaultGradScheme_
+    (
+        gradSchemes_.name() + ".default",
+        tokenList()
+    ),
+    snGradSchemes_
+    (
+        ITstream
+        (
+            objectPath() + ".snGradSchemes",
+            tokenList()
+        )()
+    ),
+    defaultSnGradScheme_
+    (
+        snGradSchemes_.name() + ".default",
+        tokenList()
+    ),
+    laplacianSchemes_
+    (
+        ITstream
+        (
+            objectPath() + ".laplacianSchemes",
+            tokenList()
+        )()
+    ),
+    defaultLaplacianScheme_
+    (
+        laplacianSchemes_.name() + ".default",
+        tokenList()
+    ),
+    fluxRequired_
+    (
+        ITstream
+        (
+            objectPath() + ".fluxRequired",
+            tokenList()
+        )()
+    ),
+    defaultFluxRequired_(false),
+    steady_(false)
+{
+    if
+    (
+        readOpt() == IOobject::MUST_READ
+     || readOpt() == IOobject::MUST_READ_IF_MODIFIED
+     || (readOpt() == IOobject::READ_IF_PRESENT && headerOk())
+    )
+    {
+        read(schemesDict());
+    }
+}
+
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 bool Foam::fvSchemes::read()

src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -111,6 +112,11 @@ public:
         //- Construct for objectRegistry
         fvSchemes(const objectRegistry& obr);
 
+        //- Construct from objectRegistry and supplied (optional) content
+        //  (gets used in case of NO_READ or fvSchemes dictionary cannot be
+        //   read)
+        fvSchemes(const objectRegistry& obr, const dictionary& dict);
+
 
     // Member Functions
 

src/finiteVolume/finiteVolume/fvSolution/fvSolution.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -69,6 +70,12 @@ public:
         :
             solution(obr, "fvSolution")
         {}
+
+        //- Construct for objectRegistry and optional contents
+        fvSolution(const objectRegistry& obr, const dictionary& dict)
+        :
+            solution(obr, "fvSolution", dict)
+        {}
 };
 
 

src/finiteVolume/fvMesh/fvGeometryScheme/averageNeighbour/averageNeighbourFvGeometryScheme.C

@@ -0,0 +1,880 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "averageNeighbourFvGeometryScheme.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvMesh.H"
+#include "cellAspectRatio.H"
+#include "syncTools.H"
+#include "polyMeshTools.H"
+#include "unitConversion.H"
+#include "OBJstream.H"
+#include "surfaceWriter.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(averageNeighbourFvGeometryScheme, 0);
+    addToRunTimeSelectionTable
+    (
+        fvGeometryScheme,
+        averageNeighbourFvGeometryScheme,
+        dict
+    );
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::label Foam::averageNeighbourFvGeometryScheme::clipFaceTet
+(
+    const scalar minRatio,
+    const vectorField& faceCentres,
+    const vectorField& faceNormals,
+
+    vectorField& faceCorrection
+) const
+{
+    // Clip correction vector if any triangle becomes too small. Return number
+    // of correction vectors clipped
+
+    typedef Vector<solveScalar> solveVector;
+
+    const pointField& p = mesh_.points();
+
+    label nClipped = 0;
+    for (label facei = 0; facei < mesh_.nFaces(); facei++)
+    {
+        const vector& fcCorr = faceCorrection[facei];
+        if (fcCorr != vector::zero)
+        {
+            const vector& fn = faceNormals[facei];
+            const point& fc = faceCentres[facei];
+            const face& f = mesh_.faces()[facei];
+
+            forAll(f, fp)
+            {
+                const solveVector thisPt(p[f[fp]]);
+                const solveVector nextPt(p[f.fcValue(fp)]);
+                const solveVector d(nextPt-thisPt);
+
+                // Calculate triangle area with correction
+                const solveVector nCorr(d^(fc+fcCorr - thisPt));
+
+                if ((nCorr & fn) < 0)
+                {
+                    // Triangle points wrong way
+                    faceCorrection[facei] = vector::zero;
+                    nClipped++;
+                    break;
+                }
+                else
+                {
+                    // Calculate triangle area without correction
+                    const solveVector n(d^(fc - thisPt));
+                    if ((n & fn) < 0)
+                    {
+                        // Original triangle points the wrong way, new one is ok
+                    }
+                    else
+                    {
+                        // Both point correctly. Make sure triangle doesn't get
+                        // too small
+                        if (mag(nCorr) < minRatio*mag(n))
+                        {
+                            faceCorrection[facei] = vector::zero;
+                            nClipped++;
+                            break;
+                        }
+                    }
+                }
+            }
+        }
+    }
+    return returnReduce(nClipped, sumOp<label>());
+}
+
+
+void Foam::averageNeighbourFvGeometryScheme::makePyrHeights
+(
+    const pointField& cellCentres,
+    const vectorField& faceCentres,
+    const vectorField& faceNormals,
+
+    scalarField& ownHeight,
+    scalarField& neiHeight
+) const
+{
+    ownHeight.setSize(mesh_.nFaces());
+    neiHeight.setSize(mesh_.nInternalFaces());
+
+    typedef Vector<solveScalar> solveVector;
+
+    const labelList& own = mesh_.faceOwner();
+    const labelList& nei = mesh_.faceNeighbour();
+
+    for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
+    {
+        const solveVector n = faceNormals[facei];
+        const solveVector fc = faceCentres[facei];
+        ownHeight[facei] = ((fc-cellCentres[own[facei]])&n);
+        neiHeight[facei] = ((cellCentres[nei[facei]]-fc)&n);
+    }
+
+    for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); facei++)
+    {
+        const solveVector n = faceNormals[facei];
+        const solveVector fc = faceCentres[facei];
+        ownHeight[facei] = ((fc-cellCentres[own[facei]])&n);
+    }
+}
+
+
+Foam::label Foam::averageNeighbourFvGeometryScheme::clipPyramids
+(
+    const pointField& cellCentres,
+    const vectorField& faceCentres,
+    const vectorField& faceNormals,
+
+    const scalarField& minOwnHeight,
+    const scalarField& minNeiHeight,
+
+    vectorField& correction
+) const
+{
+    // Clip correction vector if any pyramid becomes too small. Return number of
+    // cells clipped
+
+    typedef Vector<solveScalar> solveVector;
+
+    const labelList& own = mesh_.faceOwner();
+    const labelList& nei = mesh_.faceNeighbour();
+
+    label nClipped = 0;
+    for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
+    {
+        const vector& n = faceNormals[facei];
+        const point& fc = faceCentres[facei];
+
+        const label ownCelli = own[facei];
+        if (correction[ownCelli] != vector::zero)
+        {
+            const solveVector ownCc(cellCentres[ownCelli]+correction[ownCelli]);
+            const scalar ownHeight = ((fc-ownCc)&n);
+            if (ownHeight < minOwnHeight[facei])
+            {
+                //Pout<< "    internalface:" << fc
+                //    << " own:" << ownCc
+                //    << " pyrHeight:" << ownHeight
+                //    << " minHeight:" << minOwnHeight[facei]
+                //    << endl;
+                correction[ownCelli] = vector::zero;
+                nClipped++;
+            }
+        }
+
+        const label neiCelli = nei[facei];
+        if (correction[neiCelli] != vector::zero)
+        {
+            const solveVector neiCc(cellCentres[neiCelli]+correction[neiCelli]);
+            const scalar neiHeight = ((neiCc-fc)&n);
+            if (neiHeight < minNeiHeight[facei])
+            {
+                //Pout<< "    internalface:" << fc
+                //    << " nei:" << neiCc
+                //    << " pyrHeight:" << neiHeight
+                //    << " minHeight:" << minNeiHeight[facei]
+                //    << endl;
+                correction[neiCelli] = vector::zero;
+                nClipped++;
+            }
+        }
+    }
+
+    for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); facei++)
+    {
+        const vector& n = faceNormals[facei];
+        const point& fc = faceCentres[facei];
+
+        const label ownCelli = own[facei];
+        if (correction[ownCelli] != vector::zero)
+        {
+            const solveVector ownCc(cellCentres[ownCelli]+correction[ownCelli]);
+            const scalar ownHeight = ((fc-ownCc)&n);
+            if (ownHeight < minOwnHeight[facei])
+            {
+                //Pout<< "    boundaryface:" << fc
+                //    << " own:" << ownCc
+                //    << " pyrHeight:" << ownHeight
+                //    << " minHeight:" << minOwnHeight[facei]
+                //    << endl;
+                correction[ownCelli] = vector::zero;
+                nClipped++;
+            }
+        }
+    }
+    return returnReduce(nClipped, sumOp<label>());
+}
+
+
+Foam::tmp<Foam::pointField>
+Foam::averageNeighbourFvGeometryScheme::averageNeighbourCentres
+(
+    const pointField& cellCentres,
+    const vectorField& faceNormals,
+    const scalarField& faceWeights
+) const
+{
+    typedef Vector<solveScalar> solveVector;
+
+    const labelList& own = mesh_.faceOwner();
+    const labelList& nei = mesh_.faceNeighbour();
+
+
+    tmp<pointField> tcc(new pointField(mesh_.nCells(), Zero));
+    pointField& cc = tcc.ref();
+
+    Field<solveScalar> cellWeights(mesh_.nCells(), Zero);
+
+    // Internal faces
+    for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
+    {
+        const vector& n = faceNormals[facei];
+        const point& ownCc = cellCentres[own[facei]];
+        const point& neiCc = cellCentres[nei[facei]];
+
+        solveVector d(neiCc-ownCc);
+
+        // 1. Normalise contribution. This increases actual non-ortho
+        // since it does not 'see' the tangential offset of neighbours
+        //neiCc = ownCc + (d&n)*n;
+
+        // 2. Remove normal contribution, i.e. get tangential vector
+        //    (= non-ortho correction vector?)
+        d -= (d&n)*n;
+
+        // Apply half to both sides (as a correction)
+        // Note: should this be linear weights instead of 0.5?
+        const scalar w = 0.5*faceWeights[facei];
+        cc[own[facei]] += w*d;
+        cellWeights[own[facei]] += w;
+
+        cc[nei[facei]] -= w*d;
+        cellWeights[nei[facei]] += w;
+    }
+
+
+    // Boundary faces. Bypass stored cell centres
+    pointField neiCellCentres;
+    syncTools::swapBoundaryCellPositions(mesh_, cellCentres, neiCellCentres);
+
+    const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
+    forAll(pbm, patchi)
+    {
+        const polyPatch& pp = pbm[patchi];
+        if (pp.coupled())
+        {
+            const labelUList& fc = pp.faceCells();
+
+            forAll(fc, i)
+            {
+                const label meshFacei = pp.start()+i;
+                const label bFacei = meshFacei-mesh_.nInternalFaces();
+
+                const vector& n = faceNormals[meshFacei];
+
+                const point& ownCc = cellCentres[fc[i]];
+                const point& neiCc = neiCellCentres[bFacei];
+
+                solveVector d(neiCc-ownCc);
+
+                // 1. Normalise contribution. This increases actual non-ortho
+                // since it does not 'see' the tangential offset of neighbours
+                //neiCc = ownCc + (d&n)*n;
+
+                // 2. Remove normal contribution, i.e. get tangential vector
+                //    (= non-ortho correction vector?)
+                d -= (d&n)*n;
+
+                // Apply half to both sides (as a correction)
+                const scalar w = 0.5*faceWeights[meshFacei];
+                cc[fc[i]] += w*d;
+                cellWeights[fc[i]] += w;
+            }
+        }
+    }
+
+    // Now cc is still the correction vector. Add to cell original centres.
+    forAll(cc, celli)
+    {
+        if (cellWeights[celli] > VSMALL)
+        {
+            cc[celli] = cellCentres[celli] + cc[celli]/cellWeights[celli];
+        }
+        else
+        {
+            cc[celli] = cellCentres[celli];
+        }
+    }
+
+    return tcc;
+}
+
+
+Foam::tmp<Foam::pointField>
+Foam::averageNeighbourFvGeometryScheme::averageCentres
+(
+//    const scalar ratio,             // Amount of change in face-triangles area
+    const pointField& cellCentres,
+    const pointField& faceCentres,
+    const vectorField& faceNormals
+) const
+{
+    typedef Vector<solveScalar> solveVector;
+
+    const labelList& own = mesh_.faceOwner();
+    const labelList& nei = mesh_.faceNeighbour();
+
+
+    tmp<pointField> tnewFc(new pointField(faceCentres));
+    pointField& newFc = tnewFc.ref();
+
+    // Internal faces
+    for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
+    {
+        const vector& n = faceNormals[facei];
+        const point& oldFc = faceCentres[facei];
+
+        const solveVector ownCc(cellCentres[own[facei]]);
+        const solveVector neiCc(cellCentres[nei[facei]]);
+
+        solveVector deltaCc(neiCc-ownCc);
+        solveVector deltaFc(oldFc-ownCc);
+
+        //solveVector d(neiCc-ownCc);
+        //// 1. Normalise contribution. This increases actual non-ortho
+        //// since it does not 'see' the tangential offset of neighbours
+        ////neiCc = ownCc + s*n;
+        //
+        //// 2. Remove normal contribution, i.e. get tangential vector
+        ////    (= non-ortho correction vector?)
+        //d -= s*n;
+        //newFc[facei] = faceCentres[facei]+d;
+
+        // Get linear weight (normal distance to face)
+        const solveScalar f = (deltaFc&n)/(deltaCc&n);
+        const solveVector avgCc((1.0-f)*ownCc + f*neiCc);
+
+        solveVector d(avgCc-oldFc);
+        // Remove normal contribution, i.e. get tangential vector
+        //    (= non-ortho correction vector?)
+        d -= (d&n)*n;
+
+//        // Clip to limit change in 
+//        d *= ratio;
+
+
+        newFc[facei] = oldFc + d;
+    }
+
+
+    // Boundary faces. Bypass stored cell centres
+    pointField neiCellCentres;
+    syncTools::swapBoundaryCellPositions(mesh_, cellCentres, neiCellCentres);
+
+    const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
+    forAll(pbm, patchi)
+    {
+        const polyPatch& pp = pbm[patchi];
+        const labelUList& fc = pp.faceCells();
+
+        if (pp.coupled())
+        {
+            forAll(fc, i)
+            {
+                // Same as internal faces
+                const label facei = pp.start()+i;
+                const label bFacei = facei-mesh_.nInternalFaces();
+
+                const vector& n = faceNormals[facei];
+                const point& oldFc = faceCentres[facei];
+
+                const solveVector ownCc(cellCentres[fc[i]]);
+                const solveVector neiCc(neiCellCentres[bFacei]);
+
+                solveVector deltaCc(neiCc-ownCc);
+                solveVector deltaFc(oldFc-ownCc);
+
+                // Get linear weight (normal distance to face)
+                const solveScalar f = (deltaFc&n)/(deltaCc&n);
+                const solveVector avgCc((1.0-f)*ownCc + f*neiCc);
+
+                solveVector d(avgCc-oldFc);
+                // Remove normal contribution, i.e. get tangential vector
+                //    (= non-ortho correction vector?)
+                d -= (d&n)*n;
+
+                newFc[facei] = oldFc + d;
+            }
+        }
+        else
+        {
+            // Zero-grad?
+            forAll(fc, i)
+            {
+                const label facei = pp.start()+i;
+
+                const vector& n = faceNormals[facei];
+                const point& oldFc = faceCentres[facei];
+                const solveVector ownCc(cellCentres[fc[i]]);
+
+                solveVector d(ownCc-oldFc);
+                // Remove normal contribution, i.e. get tangential vector
+                //    (= non-ortho correction vector?)
+                d -= (d&n)*n;
+
+                newFc[facei] = oldFc+d;
+            }
+        }
+    }
+
+    return tnewFc;
+}
+
+
+void Foam::averageNeighbourFvGeometryScheme::makeNonOrthoWeights
+(
+    const pointField& cellCentres,
+    const vectorField& faceNormals,
+
+    scalarField& cosAngles,
+    scalarField& faceWeights
+) const
+{
+    cosAngles =
+        max
+        (
+            scalar(0),
+            min
+            (
+                scalar(1),
+                polyMeshTools::faceOrthogonality
+                (
+                    mesh_,
+                    faceNormals,
+                    cellCentres
+                )
+            )
+        );
+
+
+    // Make weight: 0 for ortho faces, 1 for 90degrees non-ortho
+    //const scalarField faceWeights(scalar(1)-cosAngles);
+    faceWeights.setSize(cosAngles.size());
+    {
+        const scalar minCos = Foam::cos(degToRad(80));
+        const scalar maxCos = Foam::cos(degToRad(10));
+
+        forAll(cosAngles, facei)
+        {
+            const scalar cosAngle = cosAngles[facei];
+            if (cosAngle < minCos)
+            {
+                faceWeights[facei] = 1.0;
+            }
+            else if (cosAngle > maxCos)
+            {
+                faceWeights[facei] = 0.0;
+            }
+            else
+            {
+                faceWeights[facei] =
+                    1.0-(cosAngle-minCos)/(maxCos-minCos);
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::averageNeighbourFvGeometryScheme::averageNeighbourFvGeometryScheme
+(
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    highAspectRatioFvGeometryScheme(mesh, dict),
+    nIters_
+    (
+        dict.getCheckOrDefault<label>
+        (
+            "nIters",
+            1,
+            [&](const label& nIters)
+            {
+                return nIters >= 0;
+            }
+        )
+    ),
+    relax_
+    (
+        dict.getCheck<scalar>
+        (
+            "relax",
+            [&](const scalar& relax)
+            {
+                return relax > 0 && relax <= 1;
+            }
+        )
+    ),
+    minRatio_
+    (
+        dict.getCheckOrDefault<scalar>
+        (
+            "minRatio",
+            0.5,
+            [&](const scalar& minRatio)
+            {
+                return minRatio >= 0 && minRatio <= 1;
+            }
+        )
+    )
+{
+    if (debug)
+    {
+        Pout<< "averageNeighbourFvGeometryScheme :"
+            << " nIters:" << nIters_
+            << " relax:" << relax_
+            << " minRatio:" << minRatio_ << endl;
+    }
+
+    // Force local calculation
+    movePoints();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::averageNeighbourFvGeometryScheme::movePoints()
+{
+    if (debug)
+    {
+        Pout<< "averageNeighbourFvGeometryScheme::movePoints() : "
+            << "recalculating primitiveMesh centres" << endl;
+    }
+
+    //if
+    //(
+    //   !mesh_.hasCellCentres()
+    //&& !mesh_.hasFaceCentres()
+    //&& !mesh_.hasCellVolumes()
+    //&& !mesh_.hasFaceAreas()
+    //)
+    {
+        highAspectRatioFvGeometryScheme::movePoints();
+
+        // Note: at this point the highAspectRatioFvGeometryScheme constructor
+        //       will have already reset the primitive geometry!
+
+        vectorField faceAreas(mesh_.faceAreas());
+        const scalarField magFaceAreas(mag(faceAreas));
+        const vectorField faceNormals(faceAreas/magFaceAreas);
+
+
+        // Calculate aspectratio weights
+        // - 0 if aratio < minAspect_
+        // - 1 if aratio >= maxAspect_
+        scalarField cellWeight, faceWeight;
+        calcAspectRatioWeights(cellWeight, faceWeight);
+
+        // Relaxation
+        cellWeight *= relax_;
+        //faceWeight *= relax_;
+
+        // Calculate current pyramid heights
+        scalarField minOwnHeight;
+        scalarField minNeiHeight;
+        makePyrHeights
+        (
+            mesh_.cellCentres(),
+            mesh_.faceCentres(),
+            faceNormals,
+
+            minOwnHeight,
+            minNeiHeight
+        );
+
+        // How much is the cell centre to vary inside the cell.
+        minOwnHeight *= minRatio_;
+        minNeiHeight *= minRatio_;
+
+
+
+        autoPtr<OBJstream> osPtr;
+        autoPtr<surfaceWriter> writerPtr;
+        if (debug)
+        {
+            osPtr.set
+            (
+                new OBJstream
+                (
+                    mesh_.time().timePath()
+                  / "cellCentre_correction.obj"
+                )
+            );
+            Pout<< "averageNeighbourFvGeometryScheme::movePoints() :"
+                << " writing cell centre path to " << osPtr().name() << endl;
+
+
+            // Write current non-ortho
+            fileName outputDir =
+            (
+                mesh_.time().globalPath()
+              / functionObject::outputPrefix
+              / mesh_.pointsInstance()
+            );
+            outputDir.clean();
+            writerPtr = surfaceWriter::New
+            (
+                "ensight" //"vtk"
+                // options
+            );
+
+            // Use outputDir/TIME/surface-name
+            writerPtr->useTimeDir() = true;
+
+            writerPtr->beginTime(mesh_.time());
+
+            writerPtr->open
+            (
+                mesh_.points(),
+                mesh_.faces(),
+                (outputDir / "cosAngle"),
+                true  // merge parallel bits
+            );
+
+            writerPtr->endTime();
+        }
+
+
+        // Current cellCentres. These get adjusted to lower the
+        // non-orthogonality
+        pointField cellCentres(mesh_.cellCentres());
+
+        // Modify cell centres to be more in-line with the face normals
+        for (label iter = 0; iter < nIters_; iter++)
+        {
+            // Get neighbour average (weighted by face area). This gives
+            // preference to the dominant faces. However if the non-ortho
+            // is not caused by the dominant faces this moves to the wrong
+            // direction.
+            //tmp<pointField> tcc
+            //(
+            //    averageNeighbourCentres
+            //    (
+            //        cellCentres,
+            //        faceNormals,
+            //        magFaceAreas
+            //    )
+            //);
+
+            // Get neighbour average weighted by non-ortho. Question: how to
+            // weight boundary faces?
+            tmp<pointField> tcc;
+            {
+                scalarField cosAngles;
+                scalarField faceWeights;
+                makeNonOrthoWeights
+                (
+                    cellCentres,
+                    faceNormals,
+
+                    cosAngles,
+                    faceWeights
+                );
+
+                if (writerPtr.valid())
+                {
+                    writerPtr->beginTime(instant(scalar(iter)));
+                    writerPtr->write("cosAngles", cosAngles);
+                    writerPtr->endTime();
+                }
+
+                if (debug)
+                {
+                    forAll(cosAngles, facei)
+                    {
+                        if (cosAngles[facei] < Foam::cos(degToRad(85.0)))
+                        {
+                            Pout<< "    face:" << facei
+                                << " at:" << mesh_.faceCentres()[facei]
+                                << " cosAngle:" << cosAngles[facei]
+                                << " faceWeight:" << faceWeights[facei]
+                                << endl;
+                        }
+                    }
+                }
+
+                tcc = averageNeighbourCentres
+                (
+                    cellCentres,
+                    faceNormals,
+                    faceWeights
+                );
+            }
+
+
+            // Calculate correction for cell centres. Leave low-aspect
+            // ratio cells unaffected (i.e. correction = 0)
+            vectorField correction(cellWeight*(tcc-cellCentres));
+
+            // Clip correction vector if pyramid becomes too small
+            const label nClipped = clipPyramids
+            (
+                cellCentres,
+                mesh_.faceCentres(),
+                faceNormals,
+
+                minOwnHeight, // minimum owner pyramid height. Usually fraction
+                minNeiHeight, // of starting mesh
+
+                correction
+            );
+
+            cellCentres += correction;
+
+            if (debug)
+            {
+                forAll(cellCentres, celli)
+                {
+                    const point& cc = cellCentres[celli];
+                    osPtr().write(linePointRef(cc-correction[celli], cc));
+                }
+
+                const scalarField magCorrection(mag(correction));
+                const scalarField faceOrthogonality
+                (
+                    min
+                    (
+                        scalar(1),
+                        polyMeshTools::faceOrthogonality
+                        (
+                            mesh_,
+                            faceAreas,
+                            cellCentres
+                        )
+                    )
+                );
+                const scalarField nonOrthoAngle
+                (
+                    radToDeg
+                    (
+                        Foam::acos(faceOrthogonality)
+                    )
+                );
+                Pout<< "    iter:" << iter
+                    << " nClipped:" << nClipped
+                    << " average displacement:" << gAverage(magCorrection)
+                    << " non-ortho angle : average:" << gAverage(nonOrthoAngle)
+                    << " max:" << gMax(nonOrthoAngle) << endl;
+            }
+        }
+
+        tmp<pointField> tfc
+        (
+            averageCentres
+            (
+                cellCentres,
+                mesh_.faceCentres(),
+                faceNormals
+            )
+        );
+        vectorField faceCorrection(faceWeight*(tfc-mesh_.faceCentres()));
+        // Clip correction vector to not have min triangle shrink
+        // by more than minRatio
+        clipFaceTet
+        (
+            minRatio_,
+            mesh_.faceCentres(),
+            faceNormals,
+            faceCorrection
+        );
+        vectorField faceCentres(mesh_.faceCentres() + faceCorrection);
+
+        if (debug)
+        {
+            Pout<< "averageNeighbourFvGeometryScheme::movePoints() :"
+                << " averageNeighbour weight"
+                << " max:" << gMax(cellWeight) << " min:" << gMin(cellWeight)
+                << " average:" << gAverage(cellWeight) << endl;
+
+            // Dump lines from old to new location
+            const fileName tp(mesh_.time().timePath());
+            mkDir(tp);
+            OBJstream str(tp/"averageNeighbourCellCentres.obj");
+            Pout<< "Writing lines from old to new cell centre to " << str.name()
+                << endl;
+            forAll(mesh_.cellCentres(), celli)
+            {
+                const point& oldCc = mesh_.cellCentres()[celli];
+                const point& newCc = cellCentres[celli];
+                str.write(linePointRef(oldCc, newCc));
+            }
+        }
+        if (debug)
+        {
+            // Dump lines from old to new location
+            const fileName tp(mesh_.time().timePath());
+            OBJstream str(tp/"averageFaceCentres.obj");
+            Pout<< "Writing lines from old to new face centre to " << str.name()
+                << endl;
+            forAll(mesh_.faceCentres(), facei)
+            {
+                const point& oldFc = mesh_.faceCentres()[facei];
+                const point& newFc = faceCentres[facei];
+                str.write(linePointRef(oldFc, newFc));
+            }
+        }
+
+        scalarField cellVolumes(mesh_.cellVolumes());
+
+        // Store on primitiveMesh
+        //const_cast<fvMesh&>(mesh_).clearGeom();
+        const_cast<fvMesh&>(mesh_).primitiveMesh::resetGeometry
+        (
+            std::move(faceCentres),
+            std::move(faceAreas),
+            std::move(cellCentres),
+            std::move(cellVolumes)
+        );
+    }
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/fvGeometryScheme/averageNeighbour/averageNeighbourFvGeometryScheme.H

@@ -0,0 +1,178 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::averageNeighbourFvGeometryScheme
+
+Description
+    Geometry calculation scheme to minimise non-orthogonality/
+
+SourceFiles
+    averageNeighbourFvGeometryScheme.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef averageNeighbourFvGeometryScheme_H
+#define averageNeighbourFvGeometryScheme_H
+
+#include "highAspectRatioFvGeometryScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class primitiveMesh;
+
+/*---------------------------------------------------------------------------*\
+               Class averageNeighbourFvGeometryScheme Declaration
+\*---------------------------------------------------------------------------*/
+
+class averageNeighbourFvGeometryScheme
+:
+    public highAspectRatioFvGeometryScheme
+{
+private:
+
+        //- No copy construct
+        averageNeighbourFvGeometryScheme
+        (
+            const averageNeighbourFvGeometryScheme&
+        ) = delete;
+
+        //- No copy assignment
+        void operator=(const averageNeighbourFvGeometryScheme&) = delete;
+
+
+protected:
+
+        //- Number of averaging iterations
+        const label nIters_;
+
+        //- Blending between old-iteration cell centres and current average
+        const scalar relax_;
+
+        //- Clipping for pyramid heights - allowable shrinkage as fraction
+        //  of original
+        const scalar minRatio_;
+
+
+        //- Clip face-centre correction vector if new triangle area
+        //- would get below min. Return number of clipped faces.
+        label clipFaceTet
+        (
+            const scalar minRatio,
+            const vectorField& faceCentres,
+            const vectorField& faceNormals,
+            vectorField& faceCorrection
+        ) const;
+
+        //- Calculate pyramid heights
+        void makePyrHeights
+        (
+            const pointField& cellCentres,
+            const vectorField& faceCentres,
+            const vectorField& faceNormals,
+
+            scalarField& ownHeight,
+            scalarField& neiHeight
+        ) const;
+
+        //- Clip correction vector if new pyramid height would get below min.
+        //- Return number of clipped cells.
+        label clipPyramids
+        (
+            const pointField& cellCentres,
+            const vectorField& faceCentres,
+            const vectorField& faceNormals,
+
+            const scalarField& minOwnHeight,
+            const scalarField& minNeiHeight,
+
+            vectorField& correction
+        ) const;
+
+        //- Average neighbouring cell centres to minimise non-ortho
+        tmp<pointField> averageNeighbourCentres
+        (
+            const pointField& cellCentres,
+            const vectorField& faceNormals,
+            const scalarField& faceWeights
+        ) const;
+
+        tmp<pointField> averageCentres
+        (
+            const pointField& cellCentres,
+            const pointField& faceCentres,
+            const vectorField& faceNormals
+        ) const;
+
+        //- Make weights based on non-orthogonality
+        void makeNonOrthoWeights
+        (
+            const pointField& cellCentres,
+            const vectorField& faceNormals,
+
+            scalarField& cosAngles,
+            scalarField& faceWeights
+        ) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("averageNeighbour");
+
+
+    // Constructors
+
+        //- Construct from mesh
+        averageNeighbourFvGeometryScheme
+        (
+            const fvMesh& mesh,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~averageNeighbourFvGeometryScheme() = default;
+
+
+    // Member Functions
+
+        //- Do what is necessary if the mesh has moved
+        virtual void movePoints();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/fvGeometryScheme/basic/basicFvGeometryScheme.C

@@ -0,0 +1,388 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "basicFvGeometryScheme.H"
+#include "addToRunTimeSelectionTable.H"
+#include "surfaceFields.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(basicFvGeometryScheme, 0);
+    addToRunTimeSelectionTable(fvGeometryScheme, basicFvGeometryScheme, dict);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::basicFvGeometryScheme::basicFvGeometryScheme
+(
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    fvGeometryScheme(mesh, dict)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::basicFvGeometryScheme::movePoints()
+{
+    if (debug)
+    {
+        Pout<< "basicFvGeometryScheme::movePoints() : "
+            << "recalculating primitiveMesh centres" << endl;
+    }
+    // Use lower level to calculate the geometry
+    const_cast<fvMesh&>(mesh_).primitiveMesh::updateGeom();
+}
+
+
+Foam::tmp<Foam::surfaceScalarField> Foam::basicFvGeometryScheme::weights() const
+{
+    if (debug)
+    {
+        Pout<< "basicFvGeometryScheme::weights() : "
+            << "Constructing weighting factors for face interpolation"
+            << endl;
+    }
+
+    tmp<surfaceScalarField> tweights
+    (
+        new surfaceScalarField
+        (
+            IOobject
+            (
+                "weights",
+                mesh_.pointsInstance(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false // Do not register
+            ),
+            mesh_,
+            dimless
+        )
+    );
+    surfaceScalarField& weights = tweights.ref();
+    weights.setOriented();
+
+    // Set local references to mesh data
+    // Note that we should not use fvMesh sliced fields at this point yet
+    // since this causes a loop when generating weighting factors in
+    // coupledFvPatchField evaluation phase
+    const labelUList& owner = mesh_.owner();
+    const labelUList& neighbour = mesh_.neighbour();
+
+    const vectorField& Cf = mesh_.faceCentres();
+    const vectorField& C = mesh_.cellCentres();
+    const vectorField& Sf = mesh_.faceAreas();
+
+    // ... and reference to the internal field of the weighting factors
+    scalarField& w = weights.primitiveFieldRef();
+
+    forAll(owner, facei)
+    {
+        // Note: mag in the dot-product.
+        // For all valid meshes, the non-orthogonality will be less than
+        // 90 deg and the dot-product will be positive.  For invalid
+        // meshes (d & s <= 0), this will stabilise the calculation
+        // but the result will be poor.
+        scalar SfdOwn = mag(Sf[facei] & (Cf[facei] - C[owner[facei]]));
+        scalar SfdNei = mag(Sf[facei] & (C[neighbour[facei]] - Cf[facei]));
+        w[facei] = SfdNei/(SfdOwn + SfdNei);
+    }
+
+    surfaceScalarField::Boundary& wBf = weights.boundaryFieldRef();
+
+    forAll(mesh_.boundary(), patchi)
+    {
+        mesh_.boundary()[patchi].makeWeights(wBf[patchi]);
+    }
+
+    if (debug)
+    {
+        Pout<< "basicFvGeometryScheme::weights : "
+            << "Finished constructing weighting factors for face interpolation"
+            << endl;
+    }
+    return tweights;
+}
+
+
+Foam::tmp<Foam::surfaceScalarField>
+Foam::basicFvGeometryScheme::deltaCoeffs() const
+{
+    if (debug)
+    {
+        Pout<< "basicFvGeometryScheme::deltaCoeffs() : "
+            << "Constructing differencing factors array for face gradient"
+            << endl;
+    }
+
+    // Force the construction of the weighting factors
+    // needed to make sure deltaCoeffs are calculated for parallel runs.
+    (void)mesh_.weights();
+
+    tmp<surfaceScalarField> tdeltaCoeffs
+    (
+        new surfaceScalarField
+        (
+            IOobject
+            (
+                "deltaCoeffs",
+                mesh_.pointsInstance(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false // Do not register
+            ),
+            mesh_,
+            dimless/dimLength
+        )
+    );
+    surfaceScalarField& deltaCoeffs = tdeltaCoeffs.ref();
+    deltaCoeffs.setOriented();
+
+
+    // Set local references to mesh data
+    const volVectorField& C = mesh_.C();
+    const labelUList& owner = mesh_.owner();
+    const labelUList& neighbour = mesh_.neighbour();
+
+    forAll(owner, facei)
+    {
+        deltaCoeffs[facei] = 1.0/mag(C[neighbour[facei]] - C[owner[facei]]);
+    }
+
+    surfaceScalarField::Boundary& deltaCoeffsBf =
+        deltaCoeffs.boundaryFieldRef();
+
+    forAll(deltaCoeffsBf, patchi)
+    {
+        const fvPatch& p = mesh_.boundary()[patchi];
+        deltaCoeffsBf[patchi] = 1.0/mag(p.delta());
+
+        // Optionally correct
+        p.makeDeltaCoeffs(deltaCoeffsBf[patchi]);
+    }
+
+    return tdeltaCoeffs;
+}
+
+
+Foam::tmp<Foam::surfaceScalarField>
+Foam::basicFvGeometryScheme::nonOrthDeltaCoeffs() const
+{
+    if (debug)
+    {
+        Pout<< "basicFvGeometryScheme::nonOrthDeltaCoeffs() : "
+            << "Constructing differencing factors array for face gradient"
+            << endl;
+    }
+
+    // Force the construction of the weighting factors
+    // needed to make sure deltaCoeffs are calculated for parallel runs.
+    weights();
+
+    tmp<surfaceScalarField> tnonOrthDeltaCoeffs
+    (
+        new surfaceScalarField
+        (
+            IOobject
+            (
+                "nonOrthDeltaCoeffs",
+                mesh_.pointsInstance(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false // Do not register
+            ),
+            mesh_,
+            dimless/dimLength
+        )
+    );
+    surfaceScalarField& nonOrthDeltaCoeffs = tnonOrthDeltaCoeffs.ref();
+    nonOrthDeltaCoeffs.setOriented();
+
+
+    // Set local references to mesh data
+    const volVectorField& C = mesh_.C();
+    const labelUList& owner = mesh_.owner();
+    const labelUList& neighbour = mesh_.neighbour();
+    const surfaceVectorField& Sf = mesh_.Sf();
+    const surfaceScalarField& magSf = mesh_.magSf();
+
+    forAll(owner, facei)
+    {
+        vector delta = C[neighbour[facei]] - C[owner[facei]];
+        vector unitArea = Sf[facei]/magSf[facei];
+
+        // Standard cell-centre distance form
+        //NonOrthDeltaCoeffs[facei] = (unitArea & delta)/magSqr(delta);
+
+        // Slightly under-relaxed form
+        //NonOrthDeltaCoeffs[facei] = 1.0/mag(delta);
+
+        // More under-relaxed form
+        //NonOrthDeltaCoeffs[facei] = 1.0/(mag(unitArea & delta) + VSMALL);
+
+        // Stabilised form for bad meshes
+        nonOrthDeltaCoeffs[facei] = 1.0/max(unitArea & delta, 0.05*mag(delta));
+    }
+
+    surfaceScalarField::Boundary& nonOrthDeltaCoeffsBf =
+        nonOrthDeltaCoeffs.boundaryFieldRef();
+
+    forAll(nonOrthDeltaCoeffsBf, patchi)
+    {
+        fvsPatchScalarField& patchDeltaCoeffs = nonOrthDeltaCoeffsBf[patchi];
+
+        const fvPatch& p = patchDeltaCoeffs.patch();
+
+        const vectorField patchDeltas(mesh_.boundary()[patchi].delta());
+
+        forAll(p, patchFacei)
+        {
+            vector unitArea =
+                Sf.boundaryField()[patchi][patchFacei]
+               /magSf.boundaryField()[patchi][patchFacei];
+
+            const vector& delta = patchDeltas[patchFacei];
+
+            patchDeltaCoeffs[patchFacei] =
+                1.0/max(unitArea & delta, 0.05*mag(delta));
+        }
+
+        // Optionally correct
+        p.makeNonOrthoDeltaCoeffs(patchDeltaCoeffs);
+    }
+    return tnonOrthDeltaCoeffs;
+}
+
+
+Foam::tmp<Foam::surfaceVectorField>
+Foam::basicFvGeometryScheme::nonOrthCorrectionVectors() const
+{
+    if (debug)
+    {
+        Pout<< "surfaceInterpolation::makeNonOrthCorrectionVectors() : "
+            << "Constructing non-orthogonal correction vectors"
+            << endl;
+    }
+
+    tmp<surfaceVectorField> tnonOrthCorrectionVectors
+    (
+        new surfaceVectorField
+        (
+            IOobject
+            (
+                "nonOrthCorrectionVectors",
+                mesh_.pointsInstance(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false // Do not register
+            ),
+            mesh_,
+            dimless
+        )
+    );
+    surfaceVectorField& corrVecs = tnonOrthCorrectionVectors.ref();
+    corrVecs.setOriented();
+
+    // Set local references to mesh data
+    const volVectorField& C = mesh_.C();
+    const labelUList& owner = mesh_.owner();
+    const labelUList& neighbour = mesh_.neighbour();
+    const surfaceVectorField& Sf = mesh_.Sf();
+    const surfaceScalarField& magSf = mesh_.magSf();
+    tmp<surfaceScalarField> tNonOrthDeltaCoeffs(nonOrthDeltaCoeffs());
+    const surfaceScalarField& NonOrthDeltaCoeffs = tNonOrthDeltaCoeffs();
+
+    forAll(owner, facei)
+    {
+        vector unitArea = Sf[facei]/magSf[facei];
+        vector delta = C[neighbour[facei]] - C[owner[facei]];
+
+        corrVecs[facei] = unitArea - delta*NonOrthDeltaCoeffs[facei];
+    }
+
+    // Boundary correction vectors set to zero for boundary patches
+    // and calculated consistently with internal corrections for
+    // coupled patches
+
+    surfaceVectorField::Boundary& corrVecsBf = corrVecs.boundaryFieldRef();
+
+    forAll(corrVecsBf, patchi)
+    {
+        fvsPatchVectorField& patchCorrVecs = corrVecsBf[patchi];
+
+        const fvPatch& p = patchCorrVecs.patch();
+
+        if (!patchCorrVecs.coupled())
+        {
+            patchCorrVecs = Zero;
+        }
+        else
+        {
+            const fvsPatchScalarField& patchNonOrthDeltaCoeffs =
+                NonOrthDeltaCoeffs.boundaryField()[patchi];
+
+            const vectorField patchDeltas(mesh_.boundary()[patchi].delta());
+
+            forAll(p, patchFacei)
+            {
+                vector unitArea =
+                    Sf.boundaryField()[patchi][patchFacei]
+                   /magSf.boundaryField()[patchi][patchFacei];
+
+                const vector& delta = patchDeltas[patchFacei];
+
+                patchCorrVecs[patchFacei] =
+                    unitArea - delta*patchNonOrthDeltaCoeffs[patchFacei];
+            }
+        }
+
+        // Optionally correct
+        p.makeNonOrthoCorrVectors(patchCorrVecs);
+    }
+
+    if (debug)
+    {
+        Pout<< "surfaceInterpolation::makeNonOrthCorrectionVectors() : "
+            << "Finished constructing non-orthogonal correction vectors"
+            << endl;
+    }
+    return tnonOrthCorrectionVectors;
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/fvGeometryScheme/basic/basicFvGeometryScheme.H

@@ -0,0 +1,107 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::basicFvGeometryScheme
+
+Description
+    Default geometry calculation scheme. Slight stabilisation for bad meshes.
+
+SourceFiles
+    basicFvGeometryScheme.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef basicFvGeometryScheme_H
+#define basicFvGeometryScheme_H
+
+#include "fvGeometryScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class basicFvGeometryScheme Declaration
+\*---------------------------------------------------------------------------*/
+
+class basicFvGeometryScheme
+:
+    public fvGeometryScheme
+{
+    // Private Member Functions
+
+        //- No copy construct
+        basicFvGeometryScheme(const basicFvGeometryScheme&) = delete;
+
+        //- No copy assignment
+        void operator=(const basicFvGeometryScheme&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("basic");
+
+
+    // Constructors
+
+        //- Construct from mesh
+        basicFvGeometryScheme(const fvMesh& mesh, const dictionary& dict);
+
+
+    //- Destructor
+    virtual ~basicFvGeometryScheme() = default;
+
+
+    // Member Functions
+
+        //- Do what is necessary if the mesh has moved
+        virtual void movePoints();
+
+        //- Return linear difference weighting factors
+        virtual tmp<surfaceScalarField> weights() const;
+
+        //- Return cell-centre difference coefficients
+        virtual tmp<surfaceScalarField> deltaCoeffs() const;
+
+        //- Return non-orthogonal cell-centre difference coefficients
+        virtual tmp<surfaceScalarField> nonOrthDeltaCoeffs() const;
+
+        //- Return non-orthogonality correction vectors
+        virtual tmp<surfaceVectorField> nonOrthCorrectionVectors() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/fvGeometryScheme/fvGeometryScheme/fvGeometryScheme.C

@@ -0,0 +1,80 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "fvGeometryScheme.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(fvGeometryScheme, 0);
+
+    defineRunTimeSelectionTable(fvGeometryScheme, dict);
+}
+
+
+// * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::fvGeometryScheme>
+Foam::fvGeometryScheme::New
+(
+    const fvMesh& mesh,
+    const dictionary& dict,
+    const word& defaultScheme
+)
+{
+    const entry* ePtr = dict.findEntry("method");
+    const word schemeName
+    (
+        ePtr
+      ? word(ePtr->stream())
+      : dict.getOrDefault<word>("type", defaultScheme)
+    );
+
+    if (debug)
+    {
+        InfoInFunction << "Geometry scheme = " << schemeName << endl;
+    }
+
+    auto cstrIter = dictConstructorTablePtr_->cfind(schemeName);
+
+    if (!cstrIter.found())
+    {
+        FatalIOErrorInLookup
+        (
+            dict,
+            "fvGeometryScheme",
+            schemeName,
+            *dictConstructorTablePtr_
+        ) << exit(FatalIOError);
+    }
+
+    return cstrIter()(mesh, dict);
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/fvGeometryScheme/fvGeometryScheme/fvGeometryScheme.H

@@ -0,0 +1,164 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::fvGeometryScheme
+
+Description
+    Abstract base class for geometry calculation schemes.
+
+SourceFiles
+    fvGeometryScheme.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fvGeometryScheme_H
+#define fvGeometryScheme_H
+
+#include "tmp.H"
+#include "surfaceFieldsFwd.H"
+#include "typeInfo.H"
+#include "runTimeSelectionTables.H"
+#include "pointField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class fvMesh;
+
+/*---------------------------------------------------------------------------*\
+                      Class fvGeometryScheme Declaration
+\*---------------------------------------------------------------------------*/
+
+class fvGeometryScheme
+:
+    public refCount
+{
+    // Private Member Functions
+
+        //- No copy construct
+        fvGeometryScheme(const fvGeometryScheme&) = delete;
+
+        //- No copy assignment
+        void operator=(const fvGeometryScheme&) = delete;
+
+
+protected:
+
+        //- Hold reference to mesh
+        const fvMesh& mesh_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("fvGeometryScheme");
+
+
+    // Declare run-time constructor selection tables
+
+        declareRunTimeSelectionTable
+        (
+            tmp,
+            fvGeometryScheme,
+            dict,
+            (
+                const fvMesh& mesh,
+                const dictionary& dict
+            ),
+            (mesh, dict)
+        );
+
+
+    // Constructors
+
+        //- Construct from mesh
+        fvGeometryScheme(const fvMesh& mesh, const dictionary& dict)
+        :
+            mesh_(mesh)
+        {}
+
+
+    // Selectors
+
+        //- Return new tmp interpolation scheme
+        static tmp<fvGeometryScheme> New
+        (
+            const fvMesh& mesh,
+            const dictionary& dict,
+            const word& defaultScheme
+        );
+
+
+    //- Destructor
+    virtual ~fvGeometryScheme() = default;
+
+
+    // Member Functions
+
+        //- Return mesh reference
+        const fvMesh& mesh() const
+        {
+            return mesh_;
+        }
+
+        //- Update basic geometric properties from provided points
+        virtual void movePoints()
+        {}
+
+        //- Return linear difference weighting factors
+        virtual tmp<surfaceScalarField> weights() const = 0;
+
+        //- Return cell-centre difference coefficients
+        virtual tmp<surfaceScalarField> deltaCoeffs() const = 0;
+
+        //- Return non-orthogonal cell-centre difference coefficients
+        virtual tmp<surfaceScalarField> nonOrthDeltaCoeffs() const = 0;
+
+        //- Return non-orthogonality correction vectors
+        virtual tmp<surfaceVectorField> nonOrthCorrectionVectors() const = 0;
+
+        ////- Selector for wall distance method. WIP. Ideally return wall
+        ////  distance or meshObject?
+        //virtual autoPtr<patchDistMethod> newPatchDistMethod
+        //(
+        //    const dictionary& dict,
+        //    const labelHashSet& patchIDs,
+        //    const word& defaultPatchDistMethod
+        //) const = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/fvGeometryScheme/highAspectRatio/cellAspectRatio.C

@@ -0,0 +1,123 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "cellAspectRatio.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(cellAspectRatio, 0);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::cellAspectRatio::cellAspectRatio(const polyMesh& mesh)
+:
+    MeshObject<polyMesh, Foam::MoveableMeshObject, cellAspectRatio>(mesh)
+{
+    calcAspectRatio();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::cellAspectRatio::~cellAspectRatio()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::cellAspectRatio::calcAspectRatio()
+{
+    if (debug)
+    {
+        InfoInFunction << "Calculating cell aspect ratio" << endl;
+    }
+
+    const polyMesh& mesh = mesh_;
+    const pointField& cellCentres = mesh.cellCentres();
+    const scalarField& cellVolumes = mesh.cellVolumes();
+    const vectorField& faceAreas = mesh.faceAreas();
+    const vectorField& faceCentres = mesh.faceCentres();
+    const cellList& cells = mesh.cells();
+    //const faceList& faces = mesh.faces();
+    //const pointField& points = mesh.points();
+
+    scalarField& aRatio = *this;
+    aRatio.setSize(mesh.nCells());
+
+    forAll(cells, celli)
+    {
+        const point& cc = cellCentres[celli];
+        const cell& cFaces = cells[celli];
+
+        scalar sumA = Zero;
+        scalar maxMag = Zero;
+
+        for (const label facei : cFaces)
+        {
+            const vector& n = faceAreas[facei];
+
+            sumA += mag(n);
+
+            //// Max distance from point to cell centre
+            //const face& f = faces[facei];
+            //for (const label pointi : f)
+            //{
+            //    const point& pt = points[pointi];
+            //    const vector d(pt-cc);
+            //    maxMag = max(maxMag, magSqr(d));
+            //}
+
+            // Max distance from face centre to cell centre
+            const point& fc = faceCentres[facei];
+            maxMag = max(maxMag, magSqr(fc-cc));
+        }
+        sumA /= cFaces.size();
+
+        // Local length scale
+        const scalar length = cellVolumes[celli]/sumA;
+
+        // Max edge length
+        maxMag = Foam::sqrt(maxMag);
+
+        //aRatio[celli] = Foam::sqrt(4.0/3.0)*maxMag/length;
+        aRatio[celli] = 2.0*maxMag/length;
+    }
+
+    if (debug)
+    {
+        InfoInFunction << "Calculated cell aspect ratio min:" << gMin(aRatio)
+            << " max:" << gMax(aRatio) << " average:" << gAverage(aRatio)
+            << endl;
+    }
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/fvGeometryScheme/highAspectRatio/cellAspectRatio.H

@@ -0,0 +1,97 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::cellAspectRatio
+
+Description
+    (Rough approximation of) cell aspect ratio
+
+SourceFiles
+    cellAspectRatio.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef cellAspectRatio_H
+#define cellAspectRatio_H
+
+#include "MeshObject.H"
+#include "polyMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class cellAspectRatio Declaration
+\*---------------------------------------------------------------------------*/
+
+class cellAspectRatio
+:
+    public MeshObject<polyMesh, MoveableMeshObject, cellAspectRatio>,
+    public scalarField
+{
+    // Private Member Functions
+
+        //- Construct aspect ratio
+        void calcAspectRatio();
+
+
+public:
+
+    // Declare name of the class and its debug switch
+    TypeName("cellAspectRatio");
+
+
+    // Constructors
+
+        //- Construct given an polyMesh
+        explicit cellAspectRatio(const polyMesh&);
+
+
+    //- Destructor
+    virtual ~cellAspectRatio();
+
+
+    // Member functions
+
+        //- Ignore mesh motion for now
+        virtual bool movePoints()
+        {
+            return false;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/fvGeometryScheme/highAspectRatio/highAspectRatioFvGeometryScheme.C

@@ -0,0 +1,466 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "highAspectRatioFvGeometryScheme.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvMesh.H"
+#include "syncTools.H"
+#include "cellAspectRatio.H"
+#include "emptyPolyPatch.H"
+#include "wedgePolyPatch.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(highAspectRatioFvGeometryScheme, 0);
+    addToRunTimeSelectionTable
+    (
+        fvGeometryScheme,
+        highAspectRatioFvGeometryScheme,
+        dict
+    );
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+//void Foam::highAspectRatioFvGeometryScheme::cellClosedness
+//(
+//    const vectorField& areas,
+//    const scalarField& vols,
+//    const tensorField& cellCoords,
+//
+//    scalarField& aratio
+//) const
+//{
+//    // From primitiveMeshTools::cellClosedness:
+//    //  calculate aspect ratio in given direction
+//    const labelList& own = mesh_.faceOwner();
+//    const labelList& nei = mesh_.faceNeighbour();
+//
+//    // Loop through cell faces and sum up the face area vectors for each cell.
+//    // This should be zero in all vector components
+//
+//    vectorField sumMagClosed(mesh_.nCells(), Zero);
+//
+//    forAll(own, facei)
+//    {
+//        // Add to owner
+//        vector& v = sumMagClosed[own[facei]];
+//        v += cmptMag(cellCoords[own[facei]] & areas[facei]);
+//    }
+//
+//    forAll(nei, facei)
+//    {
+//        // Subtract from neighbour
+//        vector& v = sumMagClosed[nei[facei]];
+//        v += cmptMag(cellCoords[nei[facei]] & areas[facei]);
+//    }
+//
+//    // Check the sums
+//    aratio.setSize(mesh_.nCells());
+//
+//    forAll(sumMagClosed, celli)
+//    {
+//        // Calculate the aspect ration as the maximum of Cartesian component
+//        // aspect ratio to the total area hydraulic area aspect ratio
+//        scalar minCmpt = VGREAT;
+//        scalar maxCmpt = -VGREAT;
+//        for (direction dir = 0; dir < vector::nComponents; dir++)
+//        {
+//            minCmpt = min(minCmpt, sumMagClosed[celli][dir]);
+//            maxCmpt = max(maxCmpt, sumMagClosed[celli][dir]);
+//        }
+//
+//        scalar aspectRatio = maxCmpt/(minCmpt + ROOTVSMALL);
+//        const scalar v = max(ROOTVSMALL, vols[celli]);
+//
+//        aspectRatio = max
+//        (
+//            aspectRatio,
+//            1.0/6.0*cmptSum(sumMagClosed[celli])/Foam::pow(v, 2.0/3.0)
+//        );
+//
+//        aratio[celli] = aspectRatio;
+//    }
+//}
+//
+//
+//void Foam::highAspectRatioFvGeometryScheme::cellDirections
+//(
+//    tensorField& T,
+//    vectorField& lambda
+//) const
+//{
+//    // Calculate principal directions in increasing order
+//
+//    T.setSize(mesh_.nCells());
+//    lambda.setSize(mesh_.nCells());
+//
+//    forAll(T, celli)
+//    {
+//        tensor J = Zero;
+//        {
+//            const List<tetIndices> cellTets
+//            (
+//                polyMeshTetDecomposition::cellTetIndices
+//                (
+//                    mesh_,
+//                    celli
+//                )
+//            );
+//            triFaceList faces(cellTets.size());
+//            forAll(cellTets, cTI)
+//            {
+//                faces[cTI] = cellTets[cTI].faceTriIs(mesh_);
+//            }
+//
+//            scalar m = 0.0;
+//            vector cM = Zero;
+//            J = Zero;
+//            momentOfInertia::massPropertiesShell
+//            (
+//                mesh_.points(),
+//                faces,
+//                1.0,
+//                m,
+//                cM,
+//                J
+//            );
+//        }
+//
+//        lambda[celli] = Foam::eigenValues(J);
+//        T[celli] = Foam::eigenVectors(J, lambda[celli]);
+//    }
+//}
+
+void Foam::highAspectRatioFvGeometryScheme::calcAspectRatioWeights
+(
+    scalarField& cellWeight,
+    scalarField& faceWeight
+) const
+{
+    //scalarField aratio;
+    //{
+    //    tensorField principalDirections;
+    //    vectorField lambdas;
+    //    cellDirections(principalDirections, lambdas);
+    //
+    //    cellClosedness
+    //    (
+    //        mesh_.faceAreas(),
+    //        mesh_.cellVolumes(),
+    //        principalDirections,
+    //        aratio
+    //    );
+    //}
+    const cellAspectRatio aratio(mesh_);
+
+    // Weighting for correction
+    // - 0 if aratio < minAspect_
+    // - 1 if aratio >= maxAspect_
+
+    scalar delta(maxAspect_-minAspect_);
+    if (delta < ROOTVSMALL)
+    {
+        delta = SMALL;
+    }
+
+    cellWeight =
+    max
+    (
+        min
+        (
+            (aratio-minAspect_)/delta,
+            1.0
+        ),
+        0.0
+    );
+
+    faceWeight.setSize(mesh_.nFaces());
+
+    for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
+    {
+        const label own = mesh_.faceOwner()[facei];
+        const label nei = mesh_.faceNeighbour()[facei];
+        faceWeight[facei] = max(cellWeight[own], cellWeight[nei]);
+    }
+    scalarField nbrCellWeight;
+    syncTools::swapBoundaryCellList
+    (
+        mesh_,
+        cellWeight,
+        nbrCellWeight
+    );
+    for
+    (
+        label facei = mesh_.nInternalFaces();
+        facei < mesh_.nFaces();
+        facei++
+    )
+    {
+        const label own = mesh_.faceOwner()[facei];
+        const label bFacei = facei-mesh_.nInternalFaces();
+        faceWeight[facei] = max(cellWeight[own], nbrCellWeight[bFacei]);
+    }
+}
+
+
+void Foam::highAspectRatioFvGeometryScheme::makeAverageCentres
+(
+    const polyMesh& mesh,
+    const pointField& p,
+    const pointField& faceAreas,
+    const scalarField& magFaceAreas,
+    pointField& faceCentres,
+    pointField& cellCentres
+)
+{
+    if (debug)
+    {
+        Pout<< "highAspectRatioFvGeometryScheme::makeAverageCentres() : "
+            << "calculating weighted average face/cell centre" << endl;
+    }
+
+    typedef Vector<solveScalar> solveVector;
+
+    const faceList& fs = mesh.faces();
+
+    // Start off from primitiveMesh faceCentres (preserved for triangles)
+    faceCentres.setSize(mesh.nFaces());
+
+    forAll(fs, facei)
+    {
+        const labelList& f = fs[facei];
+        const label nPoints = f.size();
+
+        if (nPoints == 3)
+        {
+            faceCentres[facei] = (1.0/3.0)*(p[f[0]] + p[f[1]] + p[f[2]]);
+        }
+        else
+        {
+            solveScalar sumA = 0.0;
+            solveVector sumAc = Zero;
+
+            for (label pi = 0; pi < nPoints; pi++)
+            {
+                const label nextPi(pi == nPoints-1 ? 0 : pi+1);
+                const solveVector nextPoint(p[f[nextPi]]);
+                const solveVector thisPoint(p[f[pi]]);
+
+                const solveVector eMid = 0.5*(thisPoint+nextPoint);
+                const solveScalar a = mag(nextPoint-thisPoint);
+
+                sumAc += a*eMid;
+                sumA += a;
+            }
+            // This is to deal with zero-area faces. Mark very small faces
+            // to be detected in e.g. processorPolyPatch.
+            if (sumA >= ROOTVSMALL)
+            {
+                faceCentres[facei] = sumAc/sumA;
+            }
+            else
+            {
+                // Unweighted average of points
+                sumAc = Zero;
+                for (label pi = 0; pi < nPoints; pi++)
+                {
+                    sumAc += p[f[pi]];
+                }
+                faceCentres[facei] = sumAc/nPoints;
+            }
+        }
+    }
+
+
+    cellCentres.setSize(mesh.nCells());
+    cellCentres = Zero;
+    {
+        const labelList& own = mesh.faceOwner();
+        const labelList& nei = mesh.faceNeighbour();
+
+        Field<solveScalar> cellWeights(mesh.nCells(), Zero);
+        for (label facei = 0; facei < mesh.nInternalFaces(); facei++)
+        {
+            const solveScalar magfA(magFaceAreas[facei]);
+            const solveVector weightedFc(magfA*faceCentres[facei]);
+
+            // Accumulate area-weighted face-centre
+            cellCentres[own[facei]] += weightedFc;
+            cellCentres[nei[facei]] += weightedFc;
+
+            // Accumulate weights
+            cellWeights[own[facei]] += magfA;
+            cellWeights[nei[facei]] += magfA;
+        }
+
+        const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+        for (const polyPatch& pp : pbm)
+        {
+            if (!isA<emptyPolyPatch>(pp) && !isA<wedgePolyPatch>(pp))
+            {
+                for
+                (
+                    label facei = pp.start();
+                    facei < pp.start()+pp.size();
+                    facei++
+                )
+                {
+                    const solveScalar magfA(magFaceAreas[facei]);
+                    const solveVector weightedFc(magfA*faceCentres[facei]);
+
+                    cellCentres[own[facei]] += weightedFc;
+                    cellWeights[own[facei]] += magfA;
+                }
+            }
+        }
+
+        forAll(cellCentres, celli)
+        {
+            if (mag(cellWeights[celli]) > VSMALL)
+            {
+                cellCentres[celli] /= cellWeights[celli];
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::highAspectRatioFvGeometryScheme::highAspectRatioFvGeometryScheme
+(
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    basicFvGeometryScheme(mesh, dict),
+    minAspect_(dict.get<scalar>("minAspect")),
+    maxAspect_(dict.get<scalar>("maxAspect"))
+{
+    if (maxAspect_ < minAspect_)
+    {
+        FatalIOErrorInFunction(dict)
+            << "minAspect " << minAspect_
+            << " has to be less than maxAspect " << maxAspect_
+            << exit(FatalIOError);
+    }
+    if (minAspect_ < 0 || maxAspect_ < 0)
+    {
+        FatalIOErrorInFunction(dict)
+            << "Illegal aspect ratio : minAspect:" << minAspect_
+            << " maxAspect:" << maxAspect_
+            << exit(FatalIOError);
+    }
+
+    // Force local calculation
+    movePoints();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::highAspectRatioFvGeometryScheme::movePoints()
+{
+    if (debug)
+    {
+        Pout<< "highAspectRatioFvGeometryScheme::movePoints() : "
+            << "recalculating primitiveMesh centres" << endl;
+    }
+
+    if
+    (
+       !mesh_.hasCellCentres()
+    && !mesh_.hasFaceCentres()
+    && !mesh_.hasCellVolumes()
+    && !mesh_.hasFaceAreas()
+    )
+    {
+        // Use lower level to calculate the geometry
+        const_cast<fvMesh&>(mesh_).primitiveMesh::updateGeom();
+
+        pointField avgFaceCentres;
+        pointField avgCellCentres;
+        makeAverageCentres
+        (
+            mesh_,
+            mesh_.points(),
+            mesh_.faceAreas(),
+            mag(mesh_.faceAreas()),
+            avgFaceCentres,
+            avgCellCentres
+        );
+
+
+        // Calculate aspectratio weights
+        // - 0 if aratio < minAspect_
+        // - 1 if aratio >= maxAspect_
+        scalarField cellWeight, faceWeight;
+        calcAspectRatioWeights(cellWeight, faceWeight);
+
+
+        // Weight with average ones
+        vectorField faceCentres
+        (
+            (1.0-faceWeight)*mesh_.faceCentres()
+          + faceWeight*avgFaceCentres
+        );
+        vectorField cellCentres
+        (
+            (1.0-cellWeight)*mesh_.cellCentres()
+          + cellWeight*avgCellCentres
+        );
+
+
+        if (debug)
+        {
+            Pout<< "highAspectRatioFvGeometryScheme::movePoints() :"
+                << " highAspectRatio weight"
+                << " max:" << gMax(cellWeight) << " min:" << gMin(cellWeight)
+                << " average:" << gAverage(cellWeight) << endl;
+        }
+
+        vectorField faceAreas(mesh_.faceAreas());
+        scalarField cellVolumes(mesh_.cellVolumes());
+
+        // Store on primitiveMesh
+        //const_cast<fvMesh&>(mesh_).clearGeom();
+        const_cast<fvMesh&>(mesh_).primitiveMesh::resetGeometry
+        (
+            std::move(faceCentres),
+            std::move(faceAreas),
+            std::move(cellCentres),
+            std::move(cellVolumes)
+        );
+    }
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/fvGeometryScheme/highAspectRatio/highAspectRatioFvGeometryScheme.H

@@ -0,0 +1,137 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::highAspectRatioFvGeometryScheme
+
+Description
+    Geometry calculation scheme with automatic stabilisation for high-aspect
+    ratio cells.
+
+SourceFiles
+    highAspectRatioFvGeometryScheme.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef highAspectRatioFvGeometryScheme_H
+#define highAspectRatioFvGeometryScheme_H
+
+#include "basicFvGeometryScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class primitiveMesh;
+class polyMesh;
+
+/*---------------------------------------------------------------------------*\
+               Class highAspectRatioFvGeometryScheme Declaration
+\*---------------------------------------------------------------------------*/
+
+class highAspectRatioFvGeometryScheme
+:
+    public basicFvGeometryScheme
+{
+protected:
+
+        const scalar minAspect_;
+
+        const scalar maxAspect_;
+
+
+    // Protected Member Functions
+
+        //- Calculate cell and face weight. Is 0 for cell < minAspect, 1 for
+        //  cell > maxAspect
+        void calcAspectRatioWeights
+        (
+            scalarField& cellWeight,
+            scalarField& faceWeight
+        ) const;
+
+        //- Helper : calculate (weighted) average face and cell centres
+        static void makeAverageCentres
+        (
+            const polyMesh& mesh,
+            const pointField& points,
+            const pointField& faceAreas,
+            const scalarField& magFaceAreas,
+            pointField& faceCentres,
+            pointField& cellCentres
+        );
+
+
+private:
+
+    // Private Member Functions
+
+        //- No copy construct
+        highAspectRatioFvGeometryScheme
+        (
+            const highAspectRatioFvGeometryScheme&
+        ) = delete;
+
+        //- No copy assignment
+        void operator=(const highAspectRatioFvGeometryScheme&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("highAspectRatio");
+
+
+    // Constructors
+
+        //- Construct from mesh
+        highAspectRatioFvGeometryScheme
+        (
+            const fvMesh& mesh,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~highAspectRatioFvGeometryScheme() = default;
+
+
+    // Member Functions
+
+        //- Do what is necessary if the mesh has moved
+        virtual void movePoints();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/fvGeometryScheme/stabilised/stabilisedFvGeometryScheme.C

@@ -0,0 +1,210 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "stabilisedFvGeometryScheme.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvMesh.H"
+#include "PrecisionAdaptor.H"
+#include "primitiveMeshTools.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(stabilisedFvGeometryScheme, 0);
+    addToRunTimeSelectionTable
+    (
+        fvGeometryScheme,
+        stabilisedFvGeometryScheme,
+        dict
+    );
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::stabilisedFvGeometryScheme::makeFaceCentresAndAreas
+(
+    const polyMesh& mesh,
+    const pointField& p,
+    vectorField& fCtrs,
+    vectorField& fAreas
+)
+{
+    const faceList& fs = mesh.faces();
+
+    forAll(fs, facei)
+    {
+        const labelList& f = fs[facei];
+        label nPoints = f.size();
+
+        // If the face is a triangle, do a direct calculation for efficiency
+        // and to avoid round-off error-related problems
+        if (nPoints == 3)
+        {
+            fCtrs[facei] = (1.0/3.0)*(p[f[0]] + p[f[1]] + p[f[2]]);
+            fAreas[facei] = 0.5*((p[f[1]] - p[f[0]])^(p[f[2]] - p[f[0]]));
+        }
+
+        // For more complex faces, decompose into triangles
+        else
+        {
+            typedef Vector<solveScalar> solveVector;
+
+            // Compute an estimate of the centre as the average of the points
+            solveVector fCentre = p[f[0]];
+            for (label pi = 1; pi < nPoints; pi++)
+            {
+                fCentre += solveVector(p[f[pi]]);
+            }
+            fCentre /= nPoints;
+
+            // Compute the face area normal and unit normal by summing up the
+            // normals of the triangles formed by connecting each edge to the
+            // point average.
+            solveVector sumA = Zero;
+            for (label pi = 0; pi < nPoints; pi++)
+            {
+                const label nextPi(pi == nPoints-1 ? 0 : pi+1);
+                const solveVector nextPoint(p[f[nextPi]]);
+                const solveVector thisPoint(p[f[pi]]);
+
+                const solveVector a =
+                    (nextPoint - thisPoint)^(fCentre - thisPoint);
+
+                sumA += a;
+            }
+            const solveVector sumAHat = normalised(sumA);
+
+            // Compute the area-weighted sum of the triangle centres. Note use
+            // the triangle area projected in the direction of the face normal
+            // as the weight, *not* the triangle area magnitude. Only the
+            // former makes the calculation independent of the initial estimate.
+            solveScalar sumAn = 0.0;
+            solveVector sumAnc = Zero;
+            for (label pi = 0; pi < nPoints; pi++)
+            {
+                const label nextPi(pi == nPoints-1 ? 0 : pi+1);
+                const solveVector nextPoint(p[f[nextPi]]);
+                const solveVector thisPoint(p[f[pi]]);
+
+                const solveVector c = thisPoint + nextPoint + fCentre;
+                const solveVector a =
+                    (nextPoint - thisPoint)^(fCentre - thisPoint);
+
+                const scalar an = a & sumAHat;
+
+                sumAn += an;
+                sumAnc += an*c;
+            }
+
+            // Complete calculating centres and areas. If the face is too small
+            // for the sums to be reliably divided then just set the centre to
+            // the initial estimate.
+            if (sumAn > ROOTVSMALL)
+            {
+                fCtrs[facei] = (1.0/3.0)*sumAnc/sumAn;
+            }
+            else
+            {
+                fCtrs[facei] = fCentre;
+            }
+            fAreas[facei] = 0.5*sumA;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::stabilisedFvGeometryScheme::stabilisedFvGeometryScheme
+(
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    basicFvGeometryScheme(mesh, dict)
+{
+    // Force local calculation
+    movePoints();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::stabilisedFvGeometryScheme::movePoints()
+{
+    if (debug)
+    {
+        Pout<< "stabilisedFvGeometryScheme::movePoints() : "
+            << "recalculating primitiveMesh centres" << endl;
+    }
+
+    if
+    (
+       !mesh_.hasCellCentres()
+    && !mesh_.hasFaceCentres()
+    && !mesh_.hasCellVolumes()
+    && !mesh_.hasFaceAreas()
+    )
+    {
+        vectorField faceCentres(mesh_.nFaces());
+        vectorField faceAreas(mesh_.nFaces());
+
+        makeFaceCentresAndAreas
+        (
+            mesh_,
+            mesh_.points(),
+            faceCentres,
+            faceAreas
+        );
+
+        vectorField cellCentres(mesh_.nCells());
+        scalarField cellVolumes(mesh_.nCells());
+
+        primitiveMeshTools::makeCellCentresAndVols
+        (
+            mesh_,
+            faceCentres,
+            faceAreas,
+            cellCentres,
+            cellVolumes
+        );
+
+        const_cast<fvMesh&>(mesh_).primitiveMesh::resetGeometry
+        (
+            std::move(faceCentres),
+            std::move(faceAreas),
+            std::move(cellCentres),
+            std::move(cellVolumes)
+        );
+    }
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/fvGeometryScheme/stabilised/stabilisedFvGeometryScheme.H

@@ -0,0 +1,126 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::stabilisedFvGeometryScheme
+
+Description
+    Geometry calculation scheme that implements face geometry calculation
+    using normal-component-of-area weighted triangle contributions.
+
+    This implements the Foundation 'Corrected face-centre calculations'
+    as a separate geometry scheme. Only implements the primitiveMesh parts,
+    not the individual face calculation.
+
+SourceFiles
+    stabilisedFvGeometryScheme.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef stabilisedFvGeometryScheme_H
+#define stabilisedFvGeometryScheme_H
+
+#include "basicFvGeometryScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class primitiveMesh;
+class polyMesh;
+
+/*---------------------------------------------------------------------------*\
+               Class stabilisedFvGeometryScheme Declaration
+\*---------------------------------------------------------------------------*/
+
+class stabilisedFvGeometryScheme
+:
+    public basicFvGeometryScheme
+{
+protected:
+
+    // Protected Member Functions
+
+        //- Calculate face area and centre weighted using pyramid height
+        static void makeFaceCentresAndAreas
+        (
+            const polyMesh& mesh,
+            const pointField& p,
+            vectorField& fCtrs,
+            vectorField& fAreas
+        );
+
+
+private:
+
+    // Private Member Functions
+
+        //- No copy construct
+        stabilisedFvGeometryScheme
+        (
+            const stabilisedFvGeometryScheme&
+        ) = delete;
+
+        //- No copy assignment
+        void operator=(const stabilisedFvGeometryScheme&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("stabilised");
+
+
+    // Constructors
+
+        //- Construct from mesh
+        stabilisedFvGeometryScheme
+        (
+            const fvMesh& mesh,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~stabilisedFvGeometryScheme() = default;
+
+
+    // Member Functions
+
+        //- Do what is necessary if the mesh has moved
+        virtual void movePoints();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/fvMesh.C

@@ -240,11 +240,11 @@ void Foam::fvMesh::clearOut()
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::fvMesh::fvMesh(const IOobject& io)
+Foam::fvMesh::fvMesh(const IOobject& io, const bool doInit)
 :
-    polyMesh(io),
-    surfaceInterpolation(*this),
+    polyMesh(io, doInit),
     fvSchemes(static_cast<const objectRegistry&>(*this)),
+    surfaceInterpolation(*this),
     fvSolution(static_cast<const objectRegistry&>(*this)),
     data(static_cast<const objectRegistry&>(*this)),
     boundary_(*this, boundaryMesh()),
@@ -261,6 +261,25 @@ Foam::fvMesh::fvMesh(const IOobject& io)
 {
     DebugInFunction << "Constructing fvMesh from IOobject" << endl;
 
+    if (doInit)
+    {
+        fvMesh::init(false);    // do not initialise lower levels
+    }
+}
+
+
+bool Foam::fvMesh::init(const bool doInit)
+{
+    if (doInit)
+    {
+        // Construct basic geometry calculation engine. Note: do before
+        // doing anything with primitiveMesh::cellCentres etc.
+        (void)geometry();
+
+        // Intialise my data
+        polyMesh::init(doInit);
+    }
+
     // Check the existence of the cell volumes and read if present
     // and set the storage of V00
     if (fileHandler().isFile(time().timePath()/dbDir()/"V0"))
@@ -322,6 +341,7 @@ Foam::fvMesh::fvMesh(const IOobject& io)
 
         moving(true);
     }
+    return false;
 }
 
 
@@ -344,11 +364,11 @@ Foam::fvMesh::fvMesh
         std::move(allNeighbour),
         syncPar
     ),
-    surfaceInterpolation(*this),
     fvSchemes(static_cast<const objectRegistry&>(*this)),
+    surfaceInterpolation(*this),
     fvSolution(static_cast<const objectRegistry&>(*this)),
     data(static_cast<const objectRegistry&>(*this)),
-    boundary_(*this, boundaryMesh()),
+    boundary_(*this),
     lduPtr_(nullptr),
     curTimeIndex_(time().timeIndex()),
     VPtr_(nullptr),
@@ -381,8 +401,8 @@ Foam::fvMesh::fvMesh
         std::move(cells),
         syncPar
     ),
-    surfaceInterpolation(*this),
     fvSchemes(static_cast<const objectRegistry&>(*this)),
+    surfaceInterpolation(*this),
     fvSolution(static_cast<const objectRegistry&>(*this)),
     data(static_cast<const objectRegistry&>(*this)),
     boundary_(*this),
@@ -407,6 +427,108 @@ Foam::fvMesh::fvMesh(const IOobject& io, const zero, const bool syncPar)
 {}
 
 
+Foam::fvMesh::fvMesh
+(
+    const IOobject& io,
+    const fvMesh& baseMesh,
+    pointField&& points,
+    faceList&& faces,
+    labelList&& allOwner,
+    labelList&& allNeighbour,
+    const bool syncPar
+)
+:
+    polyMesh
+    (
+        io,
+        std::move(points),
+        std::move(faces),
+        std::move(allOwner),
+        std::move(allNeighbour),
+        syncPar
+    ),
+    fvSchemes
+    (
+        static_cast<const objectRegistry&>(*this),
+        static_cast<const fvSchemes&>(baseMesh)
+    ),
+    surfaceInterpolation(*this),
+    fvSolution
+    (
+        static_cast<const objectRegistry&>(*this),
+        static_cast<const fvSolution&>(baseMesh)
+    ),
+    data
+    (
+        static_cast<const objectRegistry&>(*this),
+        static_cast<const data&>(baseMesh)
+    ),
+    boundary_(*this),
+    lduPtr_(nullptr),
+    curTimeIndex_(time().timeIndex()),
+    VPtr_(nullptr),
+    V0Ptr_(nullptr),
+    V00Ptr_(nullptr),
+    SfPtr_(nullptr),
+    magSfPtr_(nullptr),
+    CPtr_(nullptr),
+    CfPtr_(nullptr),
+    phiPtr_(nullptr)
+{
+    DebugInFunction << "Constructing fvMesh as copy and primitives" << endl;
+}
+
+
+Foam::fvMesh::fvMesh
+(
+    const IOobject& io,
+    const fvMesh& baseMesh,
+    pointField&& points,
+    faceList&& faces,
+    cellList&& cells,
+    const bool syncPar
+)
+:
+    polyMesh
+    (
+        io,
+        std::move(points),
+        std::move(faces),
+        std::move(cells),
+        syncPar
+    ),
+    fvSchemes
+    (
+        static_cast<const objectRegistry&>(*this),
+        static_cast<const fvSchemes&>(baseMesh)
+    ),
+    surfaceInterpolation(*this),
+    fvSolution
+    (
+        static_cast<const objectRegistry&>(*this),
+        static_cast<const fvSolution&>(baseMesh)
+    ),
+    data
+    (
+        static_cast<const objectRegistry&>(*this),
+        static_cast<const data&>(baseMesh)
+    ),
+    boundary_(*this),
+    lduPtr_(nullptr),
+    curTimeIndex_(time().timeIndex()),
+    VPtr_(nullptr),
+    V0Ptr_(nullptr),
+    V00Ptr_(nullptr),
+    SfPtr_(nullptr),
+    magSfPtr_(nullptr),
+    CPtr_(nullptr),
+    CfPtr_(nullptr),
+    phiPtr_(nullptr)
+{
+    DebugInFunction << "Constructing fvMesh as copy and primitives" << endl;
+}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::fvMesh::~fvMesh()
@@ -804,6 +926,13 @@ Foam::tmp<Foam::scalarField> Foam::fvMesh::movePoints(const pointField& p)
 }
 
 
+void Foam::fvMesh::updateGeom()
+{
+    // Let surfaceInterpolation handle geometry calculation
+    surfaceInterpolation::updateGeom();
+}
+
+
 void Foam::fvMesh::updateMesh(const mapPolyMesh& mpm)
 {
     DebugInFunction << endl;

src/finiteVolume/fvMesh/fvMesh.H

@@ -77,7 +77,6 @@ class volMesh;
 template<class Type>
 class fvMatrix;
 
-
 /*---------------------------------------------------------------------------*\
                            Class fvMesh Declaration
 \*---------------------------------------------------------------------------*/
@@ -86,8 +85,8 @@ class fvMesh
 :
     public polyMesh,
     public lduMesh,
-    public surfaceInterpolation,
     public fvSchemes,
+    public surfaceInterpolation,    // needs input from fvSchemes
     public fvSolution,
     public data
 {
@@ -98,7 +97,6 @@ protected:
         //- Boundary mesh
         fvBoundaryMesh boundary_;
 
-
     // Demand-driven data
 
         mutable fvMeshLduAddressing* lduPtr_;
@@ -186,7 +184,7 @@ public:
     // Constructors
 
         //- Construct from IOobject
-        explicit fvMesh(const IOobject& io);
+        explicit fvMesh(const IOobject& io, const bool doInit=true);
 
         //- Construct from IOobject or as zero-sized mesh
         //  Boundary is added using addFvPatches() member function
@@ -215,6 +213,32 @@ public:
             const bool syncPar = true
         );
 
+        //- Construct as copy (for dictionaries) and components without
+        //  boundary. Boundary is added using addFvPatches() member function
+        fvMesh
+        (
+            const IOobject& io,
+            const fvMesh& baseMesh,
+            pointField&& points,
+            faceList&& faces,
+            labelList&& allOwner,
+            labelList&& allNeighbour,
+            const bool syncPar = true
+        );
+
+        //- Construct as copy (for dictionaries) without boundary from cells
+        //  rather than owner/neighbour. Boundary is added using addFvPatches()
+        //  member function
+        fvMesh
+        (
+            const IOobject& io,
+            const fvMesh& baseMesh,
+            pointField&& points,
+            faceList&& faces,
+            cellList&& cells,
+            const bool syncPar = true
+        );
+
 
     //- Destructor
     virtual ~fvMesh();
@@ -224,6 +248,9 @@ public:
 
         // Helpers
 
+            //- Initialise all non-demand-driven data
+            virtual bool init(const bool doInit);
+
             //- Add boundary patches. Constructor helper
             void addFvPatches
             (
@@ -438,6 +465,10 @@ public:
             //- Move points, returns volumes swept by faces in motion
             virtual tmp<scalarField> movePoints(const pointField&);
 
+            //- Update all geometric data. This gets redirected up from
+            //  primitiveMesh level
+            virtual void updateGeom();
+
             //- Map all fields in time using given map.
             virtual void mapFields(const mapPolyMesh& mpm);
 

src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.H

@@ -74,7 +74,16 @@ class fvPatch
         const fvBoundaryMesh& boundaryMesh_;
 
 
-protected:
+    // Private Member Functions
+
+        //- No copy construct
+        fvPatch(const fvPatch&) = delete;
+
+        //- No copy assignment
+        void operator=(const fvPatch&) = delete;
+
+
+public:
 
     // Protected Member Functions
 
@@ -96,11 +105,6 @@ protected:
         //- Correct patches after moving points
         virtual void movePoints();
 
-        //- No copy construct
-        fvPatch(const fvPatch&) = delete;
-
-        //- No copy assignment
-        void operator=(const fvPatch&) = delete;
 
 public:
 

src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2017 OpenCFD Ltd.
+    Copyright (C) 2017-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -29,11 +29,12 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
+#include "surfaceInterpolation.H"
 #include "fvMesh.H"
 #include "volFields.H"
 #include "surfaceFields.H"
-#include "demandDrivenData.H"
 #include "coupledFvPatch.H"
+#include "basicFvGeometryScheme.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -47,10 +48,10 @@ namespace Foam
 
 void Foam::surfaceInterpolation::clearOut()
 {
-    deleteDemandDrivenData(weights_);
-    deleteDemandDrivenData(deltaCoeffs_);
-    deleteDemandDrivenData(nonOrthDeltaCoeffs_);
-    deleteDemandDrivenData(nonOrthCorrectionVectors_);
+    weights_.clear();
+    deltaCoeffs_.clear();
+    nonOrthDeltaCoeffs_.clear();
+    nonOrthCorrectionVectors_.clear();
 }
 
 
@@ -76,358 +77,107 @@ Foam::surfaceInterpolation::~surfaceInterpolation()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-const Foam::surfaceScalarField& Foam::surfaceInterpolation::weights() const
+const Foam::fvGeometryScheme& Foam::surfaceInterpolation::geometry() const
 {
-    if (!weights_)
+    if (!geometryPtr_.valid())
     {
-        makeWeights();
+        geometryPtr_ = fvGeometryScheme::New
+        (
+            mesh_,
+            mesh_.schemesDict().subOrEmptyDict("geometry"),
+            basicFvGeometryScheme::typeName
+        );
     }
 
-    return (*weights_);
+    return geometryPtr_();
+}
+
+
+void Foam::surfaceInterpolation::geometry(tmp<fvGeometryScheme>& schemePtr)
+{
+    geometryPtr_ = schemePtr;
+}
+
+
+const Foam::surfaceScalarField& Foam::surfaceInterpolation::weights() const
+{
+    if (!weights_.valid())
+    {
+        weights_.set(geometry().weights().ptr());
+    }
+
+    return weights_();
 }
 
 
 const Foam::surfaceScalarField& Foam::surfaceInterpolation::deltaCoeffs() const
 {
-    if (!deltaCoeffs_)
+    if (!deltaCoeffs_.valid())
     {
-        makeDeltaCoeffs();
+        deltaCoeffs_.set(geometry().deltaCoeffs().ptr());
     }
 
-    return (*deltaCoeffs_);
+    return deltaCoeffs_();
 }
 
 
 const Foam::surfaceScalarField&
 Foam::surfaceInterpolation::nonOrthDeltaCoeffs() const
 {
-    if (!nonOrthDeltaCoeffs_)
+    if (!nonOrthDeltaCoeffs_.valid())
     {
-        makeNonOrthDeltaCoeffs();
+        nonOrthDeltaCoeffs_.set(geometry().nonOrthDeltaCoeffs().ptr());
     }
 
-    return (*nonOrthDeltaCoeffs_);
+    return nonOrthDeltaCoeffs_();
 }
 
 
 const Foam::surfaceVectorField&
 Foam::surfaceInterpolation::nonOrthCorrectionVectors() const
 {
-    if (!nonOrthCorrectionVectors_)
+    if (!nonOrthCorrectionVectors_.valid())
     {
-        makeNonOrthCorrectionVectors();
+        nonOrthCorrectionVectors_.set
+        (
+            geometry().nonOrthCorrectionVectors().ptr()
+        );
     }
 
-    return (*nonOrthCorrectionVectors_);
+    return nonOrthCorrectionVectors_();
 }
 
 
 bool Foam::surfaceInterpolation::movePoints()
 {
-    deleteDemandDrivenData(weights_);
-    deleteDemandDrivenData(deltaCoeffs_);
-    deleteDemandDrivenData(nonOrthDeltaCoeffs_);
-    deleteDemandDrivenData(nonOrthCorrectionVectors_);
+    if (debug)
+    {
+        Pout<< "surfaceInterpolation::movePoints() : "
+            << "Updating geometric properties using the fvGeometryScheme"
+            << endl;
+    }
+
+    // Do any primitive geometry calculation
+    const_cast<fvGeometryScheme&>(geometry()).movePoints();
+
+    weights_.clear();
+    deltaCoeffs_.clear();
+    nonOrthDeltaCoeffs_.clear();
+    nonOrthCorrectionVectors_.clear();
 
     return true;
 }
 
 
-void Foam::surfaceInterpolation::makeWeights() const
+void Foam::surfaceInterpolation::updateGeom()
 {
     if (debug)
     {
-        Pout<< "surfaceInterpolation::makeWeights() : "
-            << "Constructing weighting factors for face interpolation"
-            << endl;
+        Pout<< "surfaceInterpolation::updateGeom() : "
+            << "Updating geometric properties" << endl;
     }
 
-    weights_ = new surfaceScalarField
-    (
-        IOobject
-        (
-            "weights",
-            mesh_.pointsInstance(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false // Do not register
-        ),
-        mesh_,
-        dimless
-    );
-    surfaceScalarField& weights = *weights_;
-    weights.setOriented();
-
-    // Set local references to mesh data
-    // Note that we should not use fvMesh sliced fields at this point yet
-    // since this causes a loop when generating weighting factors in
-    // coupledFvPatchField evaluation phase
-    const labelUList& owner = mesh_.owner();
-    const labelUList& neighbour = mesh_.neighbour();
-
-    const vectorField& Cf = mesh_.faceCentres();
-    const vectorField& C = mesh_.cellCentres();
-    const vectorField& Sf = mesh_.faceAreas();
-
-    // ... and reference to the internal field of the weighting factors
-    scalarField& w = weights.primitiveFieldRef();
-    forAll(owner, facei)
-    {
-        // Note: mag in the dot-product.
-        // For all valid meshes, the non-orthogonality will be less than
-        // 90 deg and the dot-product will be positive.  For invalid
-        // meshes (d & s <= 0), this will stabilise the calculation
-        // but the result will be poor.
-        scalar SfdOwn = mag(Sf[facei] & (Cf[facei] - C[owner[facei]]));
-        scalar SfdNei = mag(Sf[facei] & (C[neighbour[facei]] - Cf[facei]));
-        w[facei] = SfdNei/(SfdOwn + SfdNei);
-    }
-
-    surfaceScalarField::Boundary& wBf = weights.boundaryFieldRef();
-
-    forAll(mesh_.boundary(), patchi)
-    {
-        mesh_.boundary()[patchi].makeWeights(wBf[patchi]);
-    }
-
-    if (debug)
-    {
-        Pout<< "surfaceInterpolation::makeWeights() : "
-            << "Finished constructing weighting factors for face interpolation"
-            << endl;
-    }
-}
-
-
-void Foam::surfaceInterpolation::makeDeltaCoeffs() const
-{
-    if (debug)
-    {
-        Pout<< "surfaceInterpolation::makeDeltaCoeffs() : "
-            << "Constructing differencing factors array for face gradient"
-            << endl;
-    }
-
-    // Force the construction of the weighting factors
-    // needed to make sure deltaCoeffs are calculated for parallel runs.
-    weights();
-
-    deltaCoeffs_ = new surfaceScalarField
-    (
-        IOobject
-        (
-            "deltaCoeffs",
-            mesh_.pointsInstance(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false // Do not register
-        ),
-        mesh_,
-        dimless/dimLength
-    );
-    surfaceScalarField& deltaCoeffs = *deltaCoeffs_;
-    deltaCoeffs.setOriented();
-
-
-    // Set local references to mesh data
-    const volVectorField& C = mesh_.C();
-    const labelUList& owner = mesh_.owner();
-    const labelUList& neighbour = mesh_.neighbour();
-
-    forAll(owner, facei)
-    {
-        deltaCoeffs[facei] = 1.0/mag(C[neighbour[facei]] - C[owner[facei]]);
-    }
-
-    surfaceScalarField::Boundary& deltaCoeffsBf =
-        deltaCoeffs.boundaryFieldRef();
-
-    forAll(deltaCoeffsBf, patchi)
-    {
-        const fvPatch& p = mesh_.boundary()[patchi];
-        deltaCoeffsBf[patchi] = 1.0/mag(p.delta());
-
-        // Optionally correct
-        p.makeDeltaCoeffs(deltaCoeffsBf[patchi]);
-    }
-}
-
-
-void Foam::surfaceInterpolation::makeNonOrthDeltaCoeffs() const
-{
-    if (debug)
-    {
-        Pout<< "surfaceInterpolation::makeNonOrthDeltaCoeffs() : "
-            << "Constructing differencing factors array for face gradient"
-            << endl;
-    }
-
-    // Force the construction of the weighting factors
-    // needed to make sure deltaCoeffs are calculated for parallel runs.
-    weights();
-
-    nonOrthDeltaCoeffs_ = new surfaceScalarField
-    (
-        IOobject
-        (
-            "nonOrthDeltaCoeffs",
-            mesh_.pointsInstance(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false // Do not register
-        ),
-        mesh_,
-        dimless/dimLength
-    );
-    surfaceScalarField& nonOrthDeltaCoeffs = *nonOrthDeltaCoeffs_;
-    nonOrthDeltaCoeffs.setOriented();
-
-
-    // Set local references to mesh data
-    const volVectorField& C = mesh_.C();
-    const labelUList& owner = mesh_.owner();
-    const labelUList& neighbour = mesh_.neighbour();
-    const surfaceVectorField& Sf = mesh_.Sf();
-    const surfaceScalarField& magSf = mesh_.magSf();
-
-    forAll(owner, facei)
-    {
-        vector delta = C[neighbour[facei]] - C[owner[facei]];
-        vector unitArea = Sf[facei]/magSf[facei];
-
-        // Standard cell-centre distance form
-        //NonOrthDeltaCoeffs[facei] = (unitArea & delta)/magSqr(delta);
-
-        // Slightly under-relaxed form
-        //NonOrthDeltaCoeffs[facei] = 1.0/mag(delta);
-
-        // More under-relaxed form
-        //NonOrthDeltaCoeffs[facei] = 1.0/(mag(unitArea & delta) + VSMALL);
-
-        // Stabilised form for bad meshes
-        nonOrthDeltaCoeffs[facei] = 1.0/max(unitArea & delta, 0.05*mag(delta));
-    }
-
-    surfaceScalarField::Boundary& nonOrthDeltaCoeffsBf =
-        nonOrthDeltaCoeffs.boundaryFieldRef();
-
-    forAll(nonOrthDeltaCoeffsBf, patchi)
-    {
-        fvsPatchScalarField& patchDeltaCoeffs = nonOrthDeltaCoeffsBf[patchi];
-
-        const fvPatch& p = patchDeltaCoeffs.patch();
-
-        const vectorField patchDeltas(mesh_.boundary()[patchi].delta());
-
-        forAll(p, patchFacei)
-        {
-            vector unitArea =
-                Sf.boundaryField()[patchi][patchFacei]
-               /magSf.boundaryField()[patchi][patchFacei];
-
-            const vector& delta = patchDeltas[patchFacei];
-
-            patchDeltaCoeffs[patchFacei] =
-                1.0/max(unitArea & delta, 0.05*mag(delta));
-        }
-
-        // Optionally correct
-        p.makeNonOrthoDeltaCoeffs(patchDeltaCoeffs);
-    }
-}
-
-
-void Foam::surfaceInterpolation::makeNonOrthCorrectionVectors() const
-{
-    if (debug)
-    {
-        Pout<< "surfaceInterpolation::makeNonOrthCorrectionVectors() : "
-            << "Constructing non-orthogonal correction vectors"
-            << endl;
-    }
-
-    nonOrthCorrectionVectors_ = new surfaceVectorField
-    (
-        IOobject
-        (
-            "nonOrthCorrectionVectors",
-            mesh_.pointsInstance(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false // Do not register
-        ),
-        mesh_,
-        dimless
-    );
-    surfaceVectorField& corrVecs = *nonOrthCorrectionVectors_;
-    corrVecs.setOriented();
-
-    // Set local references to mesh data
-    const volVectorField& C = mesh_.C();
-    const labelUList& owner = mesh_.owner();
-    const labelUList& neighbour = mesh_.neighbour();
-    const surfaceVectorField& Sf = mesh_.Sf();
-    const surfaceScalarField& magSf = mesh_.magSf();
-    const surfaceScalarField& NonOrthDeltaCoeffs = nonOrthDeltaCoeffs();
-
-    forAll(owner, facei)
-    {
-        vector unitArea = Sf[facei]/magSf[facei];
-        vector delta = C[neighbour[facei]] - C[owner[facei]];
-
-        corrVecs[facei] = unitArea - delta*NonOrthDeltaCoeffs[facei];
-    }
-
-    // Boundary correction vectors set to zero for boundary patches
-    // and calculated consistently with internal corrections for
-    // coupled patches
-
-    surfaceVectorField::Boundary& corrVecsBf = corrVecs.boundaryFieldRef();
-
-    forAll(corrVecsBf, patchi)
-    {
-        fvsPatchVectorField& patchCorrVecs = corrVecsBf[patchi];
-
-        const fvPatch& p = patchCorrVecs.patch();
-
-        if (!patchCorrVecs.coupled())
-        {
-            patchCorrVecs = Zero;
-        }
-        else
-        {
-            const fvsPatchScalarField& patchNonOrthDeltaCoeffs =
-                NonOrthDeltaCoeffs.boundaryField()[patchi];
-
-            const vectorField patchDeltas(mesh_.boundary()[patchi].delta());
-
-            forAll(p, patchFacei)
-            {
-                vector unitArea =
-                    Sf.boundaryField()[patchi][patchFacei]
-                   /magSf.boundaryField()[patchi][patchFacei];
-
-                const vector& delta = patchDeltas[patchFacei];
-
-                patchCorrVecs[patchFacei] =
-                    unitArea - delta*patchNonOrthDeltaCoeffs[patchFacei];
-            }
-        }
-
-        // Optionally correct
-        p.makeNonOrthoCorrVectors(patchCorrVecs);
-    }
-
-    if (debug)
-    {
-        Pout<< "surfaceInterpolation::makeNonOrthCorrectionVectors() : "
-            << "Finished constructing non-orthogonal correction vectors"
-            << endl;
-    }
+    const_cast<fvGeometryScheme&>(geometry()).movePoints();
 }
 
 

src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -48,8 +49,11 @@ SourceFiles
 namespace Foam
 {
 
+class fvMesh;
+class fvGeometryScheme;
+
 /*---------------------------------------------------------------------------*\
-                     Class surfaceInterpolation Declaration
+                    Class surfaceInterpolation Declaration
 \*---------------------------------------------------------------------------*/
 
 class surfaceInterpolation
@@ -59,34 +63,23 @@ class surfaceInterpolation
         // Reference to fvMesh
         const fvMesh& mesh_;
 
+
         // Demand-driven data
 
+            //- Geometry calculation
+            mutable tmp<fvGeometryScheme> geometryPtr_;
+
             //- Linear difference weighting factors
-            mutable surfaceScalarField* weights_;
+            mutable autoPtr<surfaceScalarField> weights_;
 
             //- Cell-centre difference coefficients
-            mutable surfaceScalarField* deltaCoeffs_;
+            mutable autoPtr<surfaceScalarField> deltaCoeffs_;
 
             //- Non-orthogonal cell-centre difference coefficients
-            mutable surfaceScalarField* nonOrthDeltaCoeffs_;
+            mutable autoPtr<surfaceScalarField> nonOrthDeltaCoeffs_;
 
             //- Non-orthogonality correction vectors
-            mutable surfaceVectorField* nonOrthCorrectionVectors_;
-
-
-    // Private Member Functions
-
-        //- Construct central-differencing weighting factors
-        void makeWeights() const;
-
-        //- Construct face-gradient difference factors
-        void makeDeltaCoeffs() const;
-
-        //- Construct face-gradient difference factors
-        void makeNonOrthDeltaCoeffs() const;
-
-        //- Construct non-orthogonality correction vectors
-        void makeNonOrthCorrectionVectors() const;
+            mutable autoPtr<surfaceVectorField> nonOrthCorrectionVectors_;
 
 
 protected:
@@ -112,26 +105,35 @@ public:
 
 
     //- Destructor
-    ~surfaceInterpolation();
+    virtual ~surfaceInterpolation();
 
 
     // Member functions
 
+        //- Return reference to geometry calculation scheme
+        virtual const fvGeometryScheme& geometry() const;
+
+        //- Set geometry calculation scheme
+        void geometry(tmp<fvGeometryScheme>&);
+
         //- Return reference to linear difference weighting factors
-        const surfaceScalarField& weights() const;
+        virtual const surfaceScalarField& weights() const;
 
         //- Return reference to cell-centre difference coefficients
-        const surfaceScalarField& deltaCoeffs() const;
+        virtual const surfaceScalarField& deltaCoeffs() const;
 
         //- Return reference to non-orthogonal cell-centre difference
         //  coefficients
-        const surfaceScalarField& nonOrthDeltaCoeffs() const;
+        virtual const surfaceScalarField& nonOrthDeltaCoeffs() const;
 
         //- Return reference to non-orthogonality correction vectors
-        const surfaceVectorField& nonOrthCorrectionVectors() const;
+        virtual const surfaceVectorField& nonOrthCorrectionVectors() const;
 
         //- Do what is necessary if the mesh has moved
-        bool movePoints();
+        virtual bool movePoints();
+
+        //- Update all geometric data
+        virtual void updateGeom();
 };
 
 

src/mesh/snappyHexMesh/snappyHexMeshDriver/layerParameters/layerParameters.C

@@ -35,6 +35,20 @@ License
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
+const Foam::Enum
+<
+    Foam::layerParameters::thicknessModelType
+>
+Foam::layerParameters::thicknessModelTypeNames_
+({
+    { thicknessModelType::FIRST_AND_TOTAL, "firstAndOverall" },
+    { thicknessModelType::FIRST_AND_EXPANSION, "firstAndExpansion" },
+    { thicknessModelType::FINAL_AND_TOTAL, "finalAndOverall" },
+    { thicknessModelType::FINAL_AND_EXPANSION, "finalAndExpansion" },
+    { thicknessModelType::TOTAL_AND_EXPANSION, "overallAndExpansion" },
+    { thicknessModelType::FIRST_AND_RELATIVE_FINAL, "firstAndRelativeFinal" },
+});
+
 const Foam::scalar Foam::layerParameters::defaultConcaveAngle = 90;
 
 
@@ -44,7 +58,7 @@ Foam::scalar Foam::layerParameters::layerExpansionRatio
 (
     const label n,
     const scalar totalOverFirst
-) const
+)
 {
     if (n <= 1)
     {
@@ -66,11 +80,11 @@ Foam::scalar Foam::layerParameters::layerExpansionRatio
     if (totalOverFirst < n)
     {
         minR = 0.0;
-        maxR = pow(totalOverFirst/n, 1/(n-1));
+        maxR = pow(totalOverFirst/n, scalar(1)/(n-1));
     }
     else
     {
-        minR = pow(totalOverFirst/n, 1/(n-1));
+        minR = pow(totalOverFirst/n, scalar(1)/(n-1));
         maxR = totalOverFirst/(n - 1);
     }
 
@@ -94,6 +108,253 @@ Foam::scalar Foam::layerParameters::layerExpansionRatio
 }
 
 
+void Foam::layerParameters::readLayerParameters
+(
+    const bool verbose,
+    const dictionary& dict,
+    const thicknessModelType& spec,
+    scalar& firstLayerThickness,
+    scalar& finalLayerThickness,
+    scalar& thickness,
+    scalar& expansionRatio
+)
+{
+    // Now we have determined the layer-specification read the actual fields
+    switch (spec)
+    {
+        case FIRST_AND_TOTAL:
+            if (verbose)
+            {
+                Info<< "Layer specification as" << nl
+                    << "- first layer thickness ('firstLayerThickness')" << nl
+                    << "- overall thickness ('thickness')" << endl;
+            }
+            firstLayerThickness = readScalar
+            (
+                dict.lookup("firstLayerThickness")
+            );
+            thickness = readScalar(dict.lookup("thickness"));
+        break;
+
+        case FIRST_AND_EXPANSION:
+            if (verbose)
+            {
+                Info<< "Layer specification as" << nl
+                    << "- first layer thickness ('firstLayerThickness')" << nl
+                    << "- expansion ratio ('expansionRatio')" << endl;
+            }
+            firstLayerThickness = readScalar
+            (
+                dict.lookup("firstLayerThickness")
+            );
+            expansionRatio = readScalar(dict.lookup("expansionRatio"));
+        break;
+
+        case FINAL_AND_TOTAL:
+            if (verbose)
+            {
+                Info<< "Layer specification as" << nl
+                    << "- final layer thickness ('finalLayerThickness')" << nl
+                    << "- overall thickness ('thickness')" << endl;
+            }
+            finalLayerThickness = readScalar
+            (
+                dict.lookup("finalLayerThickness")
+            );
+            thickness = readScalar(dict.lookup("thickness"));
+        break;
+
+        case FINAL_AND_EXPANSION:
+            if (verbose)
+            {
+                Info<< "Layer specification as" << nl
+                    << "- final layer thickness ('finalLayerThickness')" << nl
+                    << "- expansion ratio ('expansionRatio')" << endl;
+            }
+            finalLayerThickness = readScalar
+            (
+                dict.lookup("finalLayerThickness")
+            );
+            expansionRatio = readScalar(dict.lookup("expansionRatio"));
+        break;
+
+        case TOTAL_AND_EXPANSION:
+            if (verbose)
+            {
+                Info<< "Layer specification as" << nl
+                    << "- overall thickness ('thickness')" << nl
+                    << "- expansion ratio ('expansionRatio')" << endl;
+            }
+            thickness = readScalar(dict.lookup("thickness"));
+            expansionRatio = readScalar(dict.lookup("expansionRatio"));
+        break;
+
+        case FIRST_AND_RELATIVE_FINAL:
+            if (verbose)
+            {
+                Info<< "Layer specification as" << nl
+                    << "- absolute first layer thickness"
+                    << " ('firstLayerThickness')"
+                    << nl
+                    << "- and final layer thickness"
+                    << " ('finalLayerThickness')" << nl
+                    << endl;
+            }
+            firstLayerThickness = readScalar
+            (
+                dict.lookup("firstLayerThickness")
+            );
+            finalLayerThickness = readScalar
+            (
+                dict.lookup("finalLayerThickness")
+            );
+        break;
+
+        default:
+            FatalIOErrorIn
+            (
+                "layerParameters::layerParameters(..)",
+                dict
+            )   << "problem." << exit(FatalIOError);
+        break;
+    }
+}
+
+
+void Foam::layerParameters::calculateLayerParameters
+(
+    const thicknessModelType& spec,
+    const label nLayers,
+    scalar& firstThickness,
+    scalar& finalThickness,
+    scalar& thickness,
+    scalar& expansionRatio
+)
+{
+    // Calculate the non-read parameters
+    switch (spec)
+    {
+        case FIRST_AND_TOTAL:
+            expansionRatio = layerExpansionRatio
+            (
+                spec,
+                nLayers,
+                firstThickness,
+                VGREAT,
+                thickness,              //totalThickness
+                VGREAT                  //expansionRatio
+            );
+            finalThickness =
+                thickness
+               *finalLayerThicknessRatio(nLayers, expansionRatio);
+
+        break;
+
+        case FIRST_AND_EXPANSION:
+            thickness = layerThickness
+            (
+                spec,
+                nLayers,
+                firstThickness,         //firstThickness
+                VGREAT,                 //finalThickness
+                VGREAT,                 //totalThickness
+                expansionRatio          //expansionRatio
+            );
+            finalThickness =
+                thickness
+               *finalLayerThicknessRatio(nLayers, expansionRatio);
+
+        break;
+
+        case FINAL_AND_TOTAL:
+            firstThickness = firstLayerThickness
+            (
+                spec,
+                nLayers,
+                VGREAT,                 //firstThickness
+                VGREAT,                 //finalThickness
+                thickness,              //totalThickness
+                VGREAT                  //expansionRatio
+            );
+            expansionRatio = layerExpansionRatio
+            (
+                spec,
+                nLayers,
+                VGREAT,                 //firstThickness
+                finalThickness,         //finalThickness
+                thickness,              //totalThickness
+                VGREAT                  //expansionRatio
+            );
+
+        break;
+
+        case FINAL_AND_EXPANSION:
+            firstThickness = firstLayerThickness
+            (
+                spec,
+                nLayers,
+                VGREAT,                 //firstThickness
+                finalThickness,         //finalThickness
+                VGREAT,                 //thickness
+                expansionRatio          //expansionRatio
+            );
+            thickness = layerThickness
+            (
+                spec,
+                nLayers,
+                VGREAT,                 //firstThickness
+                finalThickness,         //finalThickness
+                VGREAT,                 //totalThickness
+                expansionRatio          //expansionRatio
+            );
+        break;
+
+        case TOTAL_AND_EXPANSION:
+            firstThickness = firstLayerThickness
+            (
+                spec,
+                nLayers,
+                VGREAT,                 //firstThickness
+                finalThickness,         //finalThickness
+                VGREAT,                 //thickness
+                expansionRatio          //expansionRatio
+            );
+            finalThickness =
+                thickness
+               *finalLayerThicknessRatio(nLayers, expansionRatio);
+
+        break;
+
+        case FIRST_AND_RELATIVE_FINAL:
+            thickness = layerThickness
+            (
+                spec,
+                nLayers,
+                firstThickness,         //firstThickness
+                finalThickness,         //finalThickness
+                VGREAT,                 //totalThickness
+                VGREAT                  //expansionRatio
+            );
+            expansionRatio = layerExpansionRatio
+            (
+                spec,
+                nLayers,
+                firstThickness,         //firstThickness
+                finalThickness,         //finalThickness
+                VGREAT,                 //totalThickness
+                VGREAT                  //expansionRatio
+            );
+
+        break;
+
+        default:
+            FatalErrorInFunction << "Illegal thicknessModel " << spec
+                << exit(FatalError);
+        break;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::layerParameters::layerParameters
@@ -105,10 +366,13 @@ Foam::layerParameters::layerParameters
 :
     dict_(dict),
     dryRun_(dryRun),
-    relativeSizes_(meshRefinement::get<bool>(dict, "relativeSizes", dryRun)),
-    additionalReporting_(dict.getOrDefault("additionalReporting", false)),
-    layerSpec_(ILLEGAL),
     numLayers_(boundaryMesh.size(), -1),
+    relativeSizes_
+    (
+        boundaryMesh.size(),
+        meshRefinement::get<bool>(dict, "relativeSizes", dryRun)
+    ),
+    layerModels_(boundaryMesh.size(), FIRST_AND_TOTAL),
     firstLayerThickness_(boundaryMesh.size(), -123),
     finalLayerThickness_(boundaryMesh.size(), -123),
     thickness_(boundaryMesh.size(), -123),
@@ -142,6 +406,7 @@ Foam::layerParameters::layerParameters
     ),
     nLayerIter_(meshRefinement::get<label>(dict, "nLayerIter", dryRun)),
     nRelaxedIter_(labelMax),
+    additionalReporting_(dict.getOrDefault("additionalReporting", false)),
     meshShrinker_
     (
         dict.getOrDefault
@@ -153,100 +418,112 @@ Foam::layerParameters::layerParameters
 {
     // Detect layer specification mode
 
-    label nSpec = 0;
+    word spec;
+    if (dict.readIfPresent("thicknessModel", spec))
+    {
+        layerModels_ = thicknessModelTypeNames_[spec];
+    }
+    else
+    {
+        // Count number of specifications
+        label nSpec = 0;
 
-    bool haveFirst = dict.found("firstLayerThickness");
-    if (haveFirst)
-    {
-        firstLayerThickness_ = scalarField
-        (
-            boundaryMesh.size(),
-            dict.get<scalar>("firstLayerThickness")
-        );
-        nSpec++;
-    }
-    bool haveFinal = dict.found("finalLayerThickness");
-    if (haveFinal)
-    {
-        finalLayerThickness_ = scalarField
-        (
-            boundaryMesh.size(),
-            dict.get<scalar>("finalLayerThickness")
-        );
-        nSpec++;
-    }
-    bool haveTotal = dict.found("thickness");
-    if (haveTotal)
-    {
-        thickness_ = scalarField
-        (
-            boundaryMesh.size(),
-            dict.get<scalar>("thickness")
-        );
-        nSpec++;
-    }
-    bool haveExp = dict.found("expansionRatio");
-    if (haveExp)
-    {
-        expansionRatio_ = scalarField
-        (
-            boundaryMesh.size(),
-            dict.get<scalar>("expansionRatio")
-        );
-        nSpec++;
+        bool haveFirst = dict.found("firstLayerThickness");
+        if (haveFirst)
+        {
+            nSpec++;
+        }
+        bool haveFinal = dict.found("finalLayerThickness");
+        if (haveFinal)
+        {
+            nSpec++;
+        }
+        bool haveTotal = dict.found("thickness");
+        if (haveTotal)
+        {
+            nSpec++;
+        }
+        bool haveExp = dict.found("expansionRatio");
+        if (haveExp)
+        {
+            nSpec++;
+        }
+
+        if (nSpec == 2 && haveFirst && haveTotal)
+        {
+            layerModels_ = FIRST_AND_TOTAL;
+            //Info<< "Layer thickness specified as first layer"
+            //    << " and overall thickness." << endl;
+        }
+        else if (nSpec == 2 && haveFirst && haveExp)
+        {
+            layerModels_ = FIRST_AND_EXPANSION;
+            //Info<< "Layer thickness specified as first layer"
+            //    << " and expansion ratio." << endl;
+        }
+        else if (nSpec == 2 && haveFinal && haveTotal)
+        {
+            layerModels_ = FINAL_AND_TOTAL;
+            //Info<< "Layer thickness specified as final layer"
+            //    << " and overall thickness." << endl;
+        }
+        else if (nSpec == 2 && haveFinal && haveExp)
+        {
+            layerModels_ = FINAL_AND_EXPANSION;
+            //Info<< "Layer thickness specified as final layer"
+            //    << " and expansion ratio." << endl;
+        }
+        else if (nSpec == 2 && haveTotal && haveExp)
+        {
+            layerModels_ = TOTAL_AND_EXPANSION;
+            //Info<< "Layer thickness specified as overall thickness"
+            //    << " and expansion ratio." << endl;
+        }
+        else if (nSpec == 2 && haveFirst && haveFinal)
+        {
+            layerModels_ = FIRST_AND_RELATIVE_FINAL;
+            //Info<< "Layer thickness specified as absolute first and"
+            //    << " relative final layer ratio." << endl;
+        }
+        else
+        {
+            FatalIOErrorInFunction(dict)
+                << "Over- or underspecified layer thickness."
+                << " Please specify" << nl
+                << "    first layer thickness ('firstLayerThickness')"
+                << " and overall thickness ('thickness') or" << nl
+                << "    first layer thickness ('firstLayerThickness')"
+                << " and expansion ratio ('expansionRatio') or" << nl
+                << "    final layer thickness ('finalLayerThickness')"
+                << " and expansion ratio ('expansionRatio') or" << nl
+                << "    final layer thickness ('finalLayerThickness')"
+                << " and overall thickness ('thickness') or" << nl
+                << "    overall thickness ('thickness')"
+                << " and expansion ratio ('expansionRatio'"
+                << exit(FatalIOError);
+        }
     }
 
 
-    if (haveFirst && haveTotal)
-    {
-        layerSpec_ = FIRST_AND_TOTAL;
-        Info<< "Layer thickness specified as first layer and overall thickness."
-            << endl;
-    }
-    else if (haveFirst && haveExp)
-    {
-        layerSpec_ = FIRST_AND_EXPANSION;
-        Info<< "Layer thickness specified as first layer and expansion ratio."
-            << endl;
-    }
-    else if (haveFinal && haveTotal)
-    {
-        layerSpec_ = FINAL_AND_TOTAL;
-        Info<< "Layer thickness specified as final layer and overall thickness."
-            << endl;
-    }
-    else if (haveFinal && haveExp)
-    {
-        layerSpec_ = FINAL_AND_EXPANSION;
-        Info<< "Layer thickness specified as final layer and expansion ratio."
-            << endl;
-    }
-    else if (haveTotal && haveExp)
-    {
-        layerSpec_ = TOTAL_AND_EXPANSION;
-        Info<< "Layer thickness specified as overall thickness"
-            << " and expansion ratio." << endl;
-    }
-
-
-    if (layerSpec_ == ILLEGAL || nSpec != 2)
-    {
-        FatalIOErrorInFunction(dict)
-            << "Over- or underspecified layer thickness."
-            << " Please specify" << nl
-            << "    first layer thickness ('firstLayerThickness')"
-            << " and overall thickness ('thickness') or" << nl
-            << "    first layer thickness ('firstLayerThickness')"
-            << " and expansion ratio ('expansionRatio') or" << nl
-            << "    final layer thickness ('finalLayerThickness')"
-            << " and expansion ratio ('expansionRatio') or" << nl
-            << "    final layer thickness ('finalLayerThickness')"
-            << " and overall thickness ('thickness') or" << nl
-            << "    overall thickness ('thickness')"
-            << " and expansion ratio ('expansionRatio'"
-            << exit(FatalIOError);
-    }
-
+    // Now we have determined the layer-specification read the actual fields
+    scalar firstThickness;
+    scalar finalThickness;
+    scalar thickness;
+    scalar expansionRatio;
+    readLayerParameters
+    (
+        true,               // verbose
+        dict,
+        layerModels_[0],    // spec
+        firstThickness,
+        finalThickness,
+        thickness,
+        expansionRatio
+    );
+    firstLayerThickness_ = firstThickness;
+    finalLayerThickness_ = finalThickness;
+    thickness_ = thickness;
+    expansionRatio_ = expansionRatio;
 
     dict.readIfPresent("nRelaxedIter", nRelaxedIter_);
 
@@ -293,88 +570,147 @@ Foam::layerParameters::layerParameters
                     numLayers_[patchi] =
                         layerDict.get<label>("nSurfaceLayers");
 
-                    switch (layerSpec_)
+                    word spec;
+                    if (layerDict.readIfPresent("thicknessModel", spec))
                     {
-                        case FIRST_AND_TOTAL:
-                            layerDict.readIfPresent
-                            (
-                                "firstLayerThickness",
-                                firstLayerThickness_[patchi]
-                            );
-                            layerDict.readIfPresent
-                            (
-                                "thickness",
-                                thickness_[patchi]
-                            );
-                        break;
+                        // If the thickness model is explicitly specified
+                        // we want full specification of all parameters
+                        layerModels_[patchi] = thicknessModelTypeNames_[spec];
+                        readLayerParameters
+                        (
+                            false,                  // verbose
+                            layerDict,
+                            layerModels_[patchi],   // spec
+                            firstLayerThickness_[patchi],
+                            finalLayerThickness_[patchi],
+                            thickness_[patchi],
+                            expansionRatio_[patchi]
+                        );
+                        minThickness_[patchi] = readScalar
+                        (
+                            layerDict.lookup("minThickness")
+                        );
+                    }
+                    else
+                    {
+                        // Optional override of thickness parameters
+                        switch (layerModels_[patchi])
+                        {
+                            case FIRST_AND_TOTAL:
+                                layerDict.readIfPresent
+                                (
+                                    "firstLayerThickness",
+                                    firstLayerThickness_[patchi]
+                                );
+                                layerDict.readIfPresent
+                                (
+                                    "thickness",
+                                    thickness_[patchi]
+                                );
+                            break;
 
-                        case FIRST_AND_EXPANSION:
-                            layerDict.readIfPresent
-                            (
-                                "firstLayerThickness",
-                                firstLayerThickness_[patchi]
-                            );
-                            layerDict.readIfPresent
-                            (
-                                "expansionRatio",
-                                expansionRatio_[patchi]
-                            );
-                        break;
+                            case FIRST_AND_EXPANSION:
+                                layerDict.readIfPresent
+                                (
+                                    "firstLayerThickness",
+                                 firstLayerThickness_[patchi]
+                                );
+                                layerDict.readIfPresent
+                                (
+                                    "expansionRatio",
+                                    expansionRatio_[patchi]
+                                );
+                            break;
 
-                        case FINAL_AND_TOTAL:
-                            layerDict.readIfPresent
-                            (
-                                "finalLayerThickness",
-                                finalLayerThickness_[patchi]
-                            );
-                            layerDict.readIfPresent
-                            (
-                                "thickness",
-                                thickness_[patchi]
-                            );
-                        break;
+                            case FINAL_AND_TOTAL:
+                                layerDict.readIfPresent
+                                (
+                                    "finalLayerThickness",
+                                    finalLayerThickness_[patchi]
+                                );
+                                layerDict.readIfPresent
+                                (
+                                    "thickness",
+                                    thickness_[patchi]
+                                );
+                            break;
 
-                        case FINAL_AND_EXPANSION:
-                            layerDict.readIfPresent
-                            (
-                                "finalLayerThickness",
-                                finalLayerThickness_[patchi]
-                            );
-                            layerDict.readIfPresent
-                            (
-                                "expansionRatio",
-                                expansionRatio_[patchi]
-                            );
-                        break;
+                            case FINAL_AND_EXPANSION:
+                                layerDict.readIfPresent
+                                (
+                                    "finalLayerThickness",
+                                    finalLayerThickness_[patchi]
+                                );
+                                layerDict.readIfPresent
+                                (
+                                    "expansionRatio",
+                                    expansionRatio_[patchi]
+                                );
+                            break;
 
-                        case TOTAL_AND_EXPANSION:
-                            layerDict.readIfPresent
-                            (
-                                "thickness",
-                                thickness_[patchi]
-                            );
-                            layerDict.readIfPresent
-                            (
-                                "expansionRatio",
-                                expansionRatio_[patchi]
-                            );
-                        break;
+                            case TOTAL_AND_EXPANSION:
+                                layerDict.readIfPresent
+                                (
+                                    "thickness",
+                                    thickness_[patchi]
+                                );
+                                layerDict.readIfPresent
+                                (
+                                    "expansionRatio",
+                                    expansionRatio_[patchi]
+                                );
+                            break;
 
-                        default:
-                            FatalIOErrorInFunction(dict)
-                                << "problem." << exit(FatalIOError);
-                        break;
+                            case FIRST_AND_RELATIVE_FINAL:
+                                layerDict.readIfPresent
+                                (
+                                    "firstLayerThickness",
+                                    firstLayerThickness_[patchi]
+                                );
+                                layerDict.readIfPresent
+                                (
+                                    "finalLayerThickness",
+                                    finalLayerThickness_[patchi]
+                                );
+                            break;
+
+                            default:
+                                FatalIOErrorInFunction(dict)
+                                    << "problem." << exit(FatalIOError);
+                            break;
+                        }
+
+                        layerDict.readIfPresent
+                        (
+                            "minThickness",
+                            minThickness_[patchi]
+                        );
                     }
 
                     layerDict.readIfPresent
                     (
-                        "minThickness",
-                        minThickness_[patchi]
+                        "relativeSizes",
+                        relativeSizes_[patchi]
                     );
                 }
             }
         }
     }
+
+
+    forAll(numLayers_, patchi)
+    {
+        // Calculate the remaining parameters
+        calculateLayerParameters
+        (
+            layerModels_[patchi],
+            numLayers_[patchi],
+            firstLayerThickness_[patchi],
+            finalLayerThickness_[patchi],
+            thickness_[patchi],
+            expansionRatio_[patchi]
+        );
+    }
 }
 
 
@@ -382,14 +718,15 @@ Foam::layerParameters::layerParameters
 
 Foam::scalar Foam::layerParameters::layerThickness
 (
+    const thicknessModelType layerSpec,
     const label nLayers,
-    const scalar firstLayerThickess,
-    const scalar finalLayerThickess,
+    const scalar firstLayerThickness,
+    const scalar finalLayerThickness,
     const scalar totalThickness,
     const scalar expansionRatio
-) const
+)
 {
-    switch (layerSpec_)
+    switch (layerSpec)
     {
         case FIRST_AND_TOTAL:
         case FINAL_AND_TOTAL:
@@ -403,11 +740,11 @@ Foam::scalar Foam::layerParameters::layerThickness
         {
             if (mag(expansionRatio-1) < SMALL)
             {
-                return firstLayerThickess * nLayers;
+                return firstLayerThickness * nLayers;
             }
             else
             {
-                return firstLayerThickess
+                return firstLayerThickness
                    *(1.0 - pow(expansionRatio, nLayers))
                    /(1.0 - expansionRatio);
             }
@@ -418,22 +755,54 @@ Foam::scalar Foam::layerParameters::layerThickness
         {
             if (mag(expansionRatio-1) < SMALL)
             {
-                return finalLayerThickess * nLayers;
+                return finalLayerThickness * nLayers;
             }
             else
             {
                 scalar invExpansion = 1.0 / expansionRatio;
-                return finalLayerThickess
+                return finalLayerThickness
                    *(1.0 - pow(invExpansion, nLayers))
                    /(1.0 - invExpansion);
             }
         }
         break;
 
+        case FIRST_AND_RELATIVE_FINAL:
+        {
+            if (mag(expansionRatio-1) < SMALL)
+            {
+                return firstLayerThickness * nLayers;
+            }
+            else
+            {
+                scalar ratio = layerExpansionRatio
+                (
+                    layerSpec,
+                    nLayers,
+                    firstLayerThickness,
+                    finalLayerThickness,
+                    totalThickness,
+                    expansionRatio
+                );
+
+                if (mag(ratio-1) < SMALL)
+                {
+                    return firstLayerThickness * nLayers;
+                }
+                else
+                {
+                    return firstLayerThickness *
+                        (1.0 - pow(ratio, nLayers))
+                      / (1.0 - ratio);
+                }
+            }
+        }
+        break;
+
         default:
         {
             FatalErrorInFunction
-                << exit(FatalError);
+                << layerSpec << exit(FatalError);
             return -VGREAT;
         }
     }
@@ -442,14 +811,15 @@ Foam::scalar Foam::layerParameters::layerThickness
 
 Foam::scalar Foam::layerParameters::layerExpansionRatio
 (
+    const thicknessModelType layerSpec,
     const label nLayers,
-    const scalar firstLayerThickess,
-    const scalar finalLayerThickess,
+    const scalar firstLayerThickness,
+    const scalar finalLayerThickness,
     const scalar totalThickness,
     const scalar expansionRatio
-) const
+)
 {
-    switch (layerSpec_)
+    switch (layerSpec)
     {
         case FIRST_AND_EXPANSION:
         case FINAL_AND_EXPANSION:
@@ -461,23 +831,58 @@ Foam::scalar Foam::layerParameters::layerExpansionRatio
 
         case FIRST_AND_TOTAL:
         {
-            return layerExpansionRatio
-            (
-                nLayers,
-                totalThickness/firstLayerThickess
-            );
+            if (firstLayerThickness < SMALL)
+            {
+                // Do what?
+                return 1;
+            }
+            else
+            {
+                return layerExpansionRatio
+                (
+                    nLayers,
+                    totalThickness/firstLayerThickness
+                );
+            }
         }
         break;
 
         case FINAL_AND_TOTAL:
         {
-            return
-                1.0
-               /layerExpansionRatio
+            if (finalLayerThickness < SMALL)
+            {
+                // Do what?
+                return 1;
+            }
+            else
+            {
+                return
+                    1.0
+                  / layerExpansionRatio
+                    (
+                        nLayers,
+                        totalThickness/finalLayerThickness
+                    );
+            }
+        }
+        break;
+
+        case FIRST_AND_RELATIVE_FINAL:
+        {
+            if (firstLayerThickness < SMALL || nLayers <= 1)
+            {
+                return 1.0;
+            }
+            else
+            {
+                // Note: at this point the finalLayerThickness is already
+                // absolute
+                return pow
                 (
-                    nLayers,
-                    totalThickness/finalLayerThickess
+                    finalLayerThickness/firstLayerThickness,
+                    1.0/(nLayers-1)
                 );
+            }
         }
         break;
 
@@ -493,24 +898,34 @@ Foam::scalar Foam::layerParameters::layerExpansionRatio
 
 Foam::scalar Foam::layerParameters::firstLayerThickness
 (
+    const thicknessModelType layerSpec,
     const label nLayers,
-    const scalar firstLayerThickess,
-    const scalar finalLayerThickess,
+    const scalar firstLayerThickness,
+    const scalar finalLayerThickness,
     const scalar totalThickness,
     const scalar expansionRatio
-) const
+)
 {
-    switch (layerSpec_)
+    switch (layerSpec)
     {
         case FIRST_AND_EXPANSION:
         case FIRST_AND_TOTAL:
+        case FIRST_AND_RELATIVE_FINAL:
         {
-            return firstLayerThickess;
+            return firstLayerThickness;
         }
 
         case FINAL_AND_EXPANSION:
         {
-            return finalLayerThickess*pow(1.0/expansionRatio, nLayers-1);
+            if (expansionRatio < SMALL)
+            {
+                // Do what?
+                return 0.0;
+            }
+            else
+            {
+                return finalLayerThickness*pow(1.0/expansionRatio, nLayers-1);
+            }
         }
         break;
 
@@ -518,13 +933,14 @@ Foam::scalar Foam::layerParameters::firstLayerThickness
         {
             scalar r = layerExpansionRatio
             (
+                layerSpec,
                 nLayers,
-                firstLayerThickess,
-                finalLayerThickess,
+                firstLayerThickness,
+                finalLayerThickness,
                 totalThickness,
                 expansionRatio
             );
-            return finalLayerThickess/pow(r, nLayers-1);
+            return finalLayerThickness/pow(r, nLayers-1);
         }
         break;
 
@@ -554,7 +970,7 @@ Foam::scalar Foam::layerParameters::finalLayerThicknessRatio
 (
     const label nLayers,
     const scalar expansionRatio
-) const
+)
 {
     if (nLayers > 0)
     {

src/mesh/snappyHexMesh/snappyHexMeshDriver/layerParameters/layerParameters.H

@@ -41,6 +41,7 @@ SourceFiles
 #include "scalarField.H"
 #include "labelList.H"
 #include "Switch.H"
+#include "Enum.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -67,14 +68,14 @@ public:
         //  - final and total thickness specified
         //  - final and expansion ratio specified
         //  - total thickness and expansion ratio specified
-        enum layerSpecification
+        enum thicknessModelType
         {
-            ILLEGAL,
             FIRST_AND_TOTAL,
             FIRST_AND_EXPANSION,
             FINAL_AND_TOTAL,
             FINAL_AND_EXPANSION,
-            TOTAL_AND_EXPANSION
+            TOTAL_AND_EXPANSION,
+            FIRST_AND_RELATIVE_FINAL
         };
 
 
@@ -85,6 +86,9 @@ private:
         //- Default angle for faces to be concave
         static const scalar defaultConcaveAngle;
 
+        //- thicknessModelType names
+        static const Enum<thicknessModelType> thicknessModelTypeNames_;
+
 
     // Private Data
 
@@ -93,28 +97,25 @@ private:
         //- In dry-run mode?
         const bool dryRun_;
 
-        //- Are sizes relative to local cell size
-        bool relativeSizes_;
-
-        //- Any additional reporting
-        bool additionalReporting_;
-
 
         // Per patch (not region!) information
 
-            //- How thickness is specified.
-            layerSpecification layerSpec_;
-
             //- How many layers to add.
             labelList numLayers_;
 
-                scalarField firstLayerThickness_;
+            //- Are sizes relative to local cell size
+            boolList relativeSizes_;
 
-                scalarField finalLayerThickness_;
+            //- How thickness is specified.
+            List<thicknessModelType> layerModels_;
 
-                scalarField thickness_;
+            scalarField firstLayerThickness_;
 
-                scalarField expansionRatio_;
+            scalarField finalLayerThickness_;
+
+            scalarField thickness_;
+
+            scalarField expansionRatio_;
 
             //- Minimum total thickness
             scalarField minThickness_;
@@ -136,6 +137,9 @@ private:
 
         label nRelaxedIter_;
 
+        //- Any additional reporting
+        bool additionalReporting_;
+
         word meshShrinker_;
 
 
@@ -143,11 +147,33 @@ private:
 
         //- Calculate expansion ratio from overall size v.s. thickness of
         //  first layer.
-        scalar layerExpansionRatio
+        static scalar layerExpansionRatio
         (
             const label n,
             const scalar totalOverFirst
-        ) const;
+        );
+
+        //- Read parameters according to thickness model
+        static void readLayerParameters
+        (
+            const bool verbose,
+            const dictionary&,
+            const thicknessModelType& spec,
+            scalar& firstLayerThickness,
+            scalar& finalLayerThickness,
+            scalar& thickness,
+            scalar& expansionRatio
+        );
+
+        void calculateLayerParameters
+        (
+            const thicknessModelType& spec,
+            const label nLayers,
+            scalar& firstLayerThickness,
+            scalar& finalLayerThickness,
+            scalar& thickness,
+            scalar& expansionRatio
+        );
 
         //- No copy construct
         layerParameters(const layerParameters&) = delete;
@@ -191,15 +217,15 @@ public:
 
             //- Are size parameters relative to inner cell size or
             //  absolute distances.
-            bool relativeSizes() const
+            const boolList& relativeSizes() const
             {
                 return relativeSizes_;
             }
 
-            //- Any additional reporting requested?
-            bool additionalReporting() const
+            //- Specification of layer thickness
+            const List<thicknessModelType>& layerModels() const
             {
-                return additionalReporting_;
+                return layerModels_;
             }
 
             // Expansion factor for layer mesh
@@ -296,6 +322,12 @@ public:
                 return nBufferCellsNoExtrude_;
             }
 
+            //- Any additional reporting requested?
+            bool additionalReporting() const
+            {
+                return additionalReporting_;
+            }
+
             //- Type of mesh shrinker
             const word& meshShrinker() const
             {
@@ -306,45 +338,48 @@ public:
         // Helper
 
             //- Determine overall thickness. Uses two of the four parameters
-            //  according to the layerSpecification
-            scalar layerThickness
+            //  according to the thicknessModel
+            static scalar layerThickness
             (
+                const thicknessModelType,
                 const label nLayers,
-                const scalar firstLayerThickess,
-                const scalar finalLayerThickess,
+                const scalar firstLayerThickness,
+                const scalar finalLayerThickness,
                 const scalar totalThickness,
                 const scalar expansionRatio
-            ) const;
+            );
 
             //- Determine expansion ratio. Uses two of the four parameters
-            //  according to the layerSpecification
-            scalar layerExpansionRatio
+            //  according to the thicknessModel
+            static scalar layerExpansionRatio
             (
+                const thicknessModelType,
                 const label nLayers,
-                const scalar firstLayerThickess,
-                const scalar finalLayerThickess,
+                const scalar firstLayerThickness,
+                const scalar finalLayerThickness,
                 const scalar totalThickness,
                 const scalar expansionRatio
-            ) const;
+            );
 
             //- Determine first layer (near-wall) thickness. Uses two of the
-            //  four parameters according to the layerSpecification
-            scalar firstLayerThickness
+            //  four parameters according to the thicknessModel
+            static scalar firstLayerThickness
             (
+                const thicknessModelType,
                 const label nLayers,
-                const scalar firstLayerThickess,
-                const scalar finalLayerThickess,
+                const scalar firstLayerThickness,
+                const scalar finalLayerThickness,
                 const scalar totalThickness,
                 const scalar expansionRatio
-            ) const;
+            );
 
             //- Determine ratio of final layer thickness to
             //  overall layer thickness
-            scalar finalLayerThicknessRatio
+            static scalar finalLayerThicknessRatio
             (
                 const label nLayers,
                 const scalar expansionRatio
-            ) const;
+            );
 };
 
 

src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C

@@ -820,6 +820,7 @@ void Foam::snappyLayerDriver::handleWarpedFaces
 (
     const indirectPrimitivePatch& pp,
     const scalar faceRatio,
+    const boolList& relativeSizes,
     const scalar edge0Len,
     const labelList& cellLevel,
     pointField& patchDisp,
@@ -828,6 +829,7 @@ void Foam::snappyLayerDriver::handleWarpedFaces
 ) const
 {
     const fvMesh& mesh = meshRefiner_.mesh();
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
     Info<< nl << "Handling cells with warped patch faces ..." << nl;
 
@@ -838,16 +840,19 @@ void Foam::snappyLayerDriver::handleWarpedFaces
     forAll(pp, i)
     {
         const face& f = pp[i];
+        label faceI = pp.addressing()[i];
+        label patchI = patches.patchID()[faceI-mesh.nInternalFaces()];
 
-        if (f.size() > 3)
+        // It is hard to calculate some length scale if not in relative
+        // mode so disable this check.
+
+        if (relativeSizes[patchI] && f.size() > 3)
         {
-            label facei = pp.addressing()[i];
-
-            label ownLevel = cellLevel[mesh.faceOwner()[facei]];
+            label ownLevel = cellLevel[mesh.faceOwner()[faceI]];
             scalar edgeLen = edge0Len/(1<<ownLevel);
 
             // Normal distance to face centre plane
-            const point& fc = mesh.faceCentres()[facei];
+            const point& fc = mesh.faceCentres()[faceI];
             const vector& fn = pp.faceNormals()[i];
 
             scalarField vProj(f.size());
@@ -1412,10 +1417,14 @@ void Foam::snappyLayerDriver::calculateLayerThickness
     minThickness.setSize(pp.nPoints());
     minThickness = GREAT;
 
-    forAll(patchIDs, i)
-    {
-        label patchi = patchIDs[i];
+    thickness.setSize(pp.nPoints());
+    thickness = GREAT;
 
+    expansionRatio.setSize(pp.nPoints());
+    expansionRatio = GREAT;
+
+    for (const label patchi : patchIDs)
+    {
         const labelList& meshPoints = patches[patchi].meshPoints();
 
         forAll(meshPoints, patchPointi)
@@ -1495,32 +1504,11 @@ void Foam::snappyLayerDriver::calculateLayerThickness
     // Now the thicknesses are set according to the minimum of connected
     // patches.
 
+    // Determine per point the max cell level of connected cells
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    // Rework relative thickness into absolute
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // by multiplying with the internal cell size.
-
-    if (layerParams.relativeSizes())
+    labelList maxPointLevel(pp.nPoints(), labelMin);
     {
-        if
-        (
-            min(layerParams.minThickness()) < 0
-         || max(layerParams.minThickness()) > 2
-        )
-        {
-            FatalErrorInFunction
-                << "Thickness should be factor of local undistorted cell size."
-                << " Valid values are [0..2]." << nl
-                << " minThickness:" << layerParams.minThickness()
-                << exit(FatalError);
-        }
-
-
-        // Determine per point the max cell level of connected cells
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-        labelList maxPointLevel(pp.nPoints(), labelMin);
-
         forAll(pp, i)
         {
             label ownLevel = cellLevel[mesh.faceOwner()[pp.addressing()[i]]];
@@ -1541,44 +1529,149 @@ void Foam::snappyLayerDriver::calculateLayerThickness
             maxEqOp<label>(),
             labelMin            // null value
         );
+    }
 
 
-        forAll(maxPointLevel, pointi)
+    // Rework relative thickness into absolute
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // by multiplying with the internal cell size.
+    // Note that we cannot loop over the patches since then points on
+    // multiple patches would get multiplied with edgeLen twice ..
+    {
+        // Multiplication factor for relative sizes
+        scalarField edgeLen(pp.nPoints(), GREAT);
+
+        labelList spec(pp.nPoints(), layerParameters::FIRST_AND_TOTAL);
+
+        bitSet isRelativePoint(mesh.nPoints());
+
+        for (const label patchi : patchIDs)
         {
-            // Find undistorted edge size for this level.
-            scalar edgeLen = edge0Len/(1<<maxPointLevel[pointi]);
-            firstLayerThickness[pointi] *= edgeLen;
-            finalLayerThickness[pointi] *= edgeLen;
-            totalThickness[pointi] *= edgeLen;
-            minThickness[pointi] *= edgeLen;
+            const labelList& meshPoints = patches[patchi].meshPoints();
+            const layerParameters::thicknessModelType patchSpec =
+                layerParams.layerModels()[patchi];
+            const bool relSize = layerParams.relativeSizes()[patchi];
+
+            for (const label meshPointi : meshPoints)
+            {
+                const label ppPointi = pp.meshPointMap()[meshPointi];
+
+                // Note: who wins if different specs?
+
+                // Calculate undistorted edge size for this level.
+                edgeLen[ppPointi] = min
+                (
+                    edgeLen[ppPointi],
+                    edge0Len/(1<<maxPointLevel[ppPointi])
+                );
+                spec[ppPointi] = max(spec[ppPointi], patchSpec);
+                isRelativePoint[meshPointi] =
+                    isRelativePoint[meshPointi]
+                 || relSize;
+            }
+        }
+
+        syncTools::syncPointList
+        (
+            mesh,
+            pp.meshPoints(),
+            edgeLen,
+            minEqOp<scalar>(),
+            GREAT                               // null value
+        );
+        syncTools::syncPointList
+        (
+            mesh,
+            pp.meshPoints(),
+            spec,
+            maxEqOp<label>(),
+            label(layerParameters::FIRST_AND_TOTAL) // null value
+        );
+        syncTools::syncPointList
+        (
+            mesh,
+            isRelativePoint,
+            orEqOp<unsigned int>(),
+            0
+        );
+
+
+
+
+        forAll(pp.meshPoints(), pointi)
+        {
+            const label meshPointi = pp.meshPoints()[pointi];
+            const layerParameters::thicknessModelType pointSpec =
+                static_cast<layerParameters::thicknessModelType>(spec[pointi]);
+
+            if (pointSpec == layerParameters::FIRST_AND_RELATIVE_FINAL)
+            {
+                // This overrules the relative sizes flag for
+                // first (always absolute) and final (always relative)
+                finalLayerThickness[pointi] *= edgeLen[pointi];
+                if (isRelativePoint[meshPointi])
+                {
+                    totalThickness[pointi] *= edgeLen[pointi];
+                    minThickness[pointi] *= edgeLen[pointi];
+                }
+            }
+            else if (isRelativePoint[meshPointi])
+            {
+                firstLayerThickness[pointi] *= edgeLen[pointi];
+                finalLayerThickness[pointi] *= edgeLen[pointi];
+                totalThickness[pointi] *= edgeLen[pointi];
+                minThickness[pointi] *= edgeLen[pointi];
+            }
+
+            thickness[pointi] = min
+            (
+                thickness[pointi],
+                layerParameters::layerThickness
+                (
+                    pointSpec,
+                    patchNLayers[pointi],
+                    firstLayerThickness[pointi],
+                    finalLayerThickness[pointi],
+                    totalThickness[pointi],
+                    expRatio[pointi]
+                )
+            );
+            expansionRatio[pointi] = min
+            (
+                expansionRatio[pointi],
+                layerParameters::layerExpansionRatio
+                (
+                    pointSpec,
+                    patchNLayers[pointi],
+                    firstLayerThickness[pointi],
+                    finalLayerThickness[pointi],
+                    totalThickness[pointi],
+                    expRatio[pointi]
+                )
+            );
         }
     }
 
+    // Synchronise the determined thicknes. Note that this should not be
+    // necessary since the inputs to the calls to layerThickness,
+    // layerExpansionRatio above are already parallel consistent
 
-
-    // Rework thickness parameters into overall thickness
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    forAll(firstLayerThickness, pointi)
-    {
-        thickness[pointi] = layerParams.layerThickness
-        (
-            patchNLayers[pointi],
-            firstLayerThickness[pointi],
-            finalLayerThickness[pointi],
-            totalThickness[pointi],
-            expRatio[pointi]
-        );
-
-        expansionRatio[pointi] = layerParams.layerExpansionRatio
-        (
-            patchNLayers[pointi],
-            firstLayerThickness[pointi],
-            finalLayerThickness[pointi],
-            totalThickness[pointi],
-            expRatio[pointi]
-        );
-    }
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        thickness,
+        minEqOp<scalar>(),
+        GREAT               // null value
+    );
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        expansionRatio,
+        minEqOp<scalar>(),
+        GREAT               // null value
+    );
 
     //Info<< "calculateLayerThickness : min:" << gMin(thickness)
     //    << " max:" << gMax(thickness) << endl;
@@ -1614,6 +1707,8 @@ void Foam::snappyLayerDriver::calculateLayerThickness
             label patchi = patchIDs[i];
 
             const labelList& meshPoints = patches[patchi].meshPoints();
+            const layerParameters::thicknessModelType spec =
+                layerParams.layerModels()[patchi];
 
             scalar sumThickness = 0;
             scalar sumNearWallThickness = 0;
@@ -1629,6 +1724,7 @@ void Foam::snappyLayerDriver::calculateLayerThickness
                     sumThickness += thickness[ppPointi];
                     sumNearWallThickness += layerParams.firstLayerThickness
                     (
+                        spec,
                         patchNLayers[ppPointi],
                         firstLayerThickness[ppPointi],
                         finalLayerThickness[ppPointi],
@@ -1806,6 +1902,8 @@ void Foam::snappyLayerDriver::getPatchDisplacement
     const indirectPrimitivePatch& pp,
     const scalarField& thickness,
     const scalarField& minThickness,
+    const scalarField& expansionRatio,
+
     pointField& patchDisp,
     labelList& patchNLayers,
     List<extrudeMode>& extrudeStatus
@@ -1836,6 +1934,98 @@ void Foam::snappyLayerDriver::getPatchDisplacement
     label nExtrudeRemove = 0;
 
 
+////XXXXXXXX
+//    {
+//        OBJstream twoStr
+//        (
+//            mesh.time().path()
+//          / "twoFacePoints_"
+//          + meshRefiner_.timeName()
+//          + ".obj"
+//        );
+//        OBJstream multiStr
+//        (
+//            mesh.time().path()
+//          / "multiFacePoints_"
+//          + meshRefiner_.timeName()
+//          + ".obj"
+//        );
+//        Pout<< "Writing points inbetween two faces on same cell to "
+//            << twoStr.name() << endl;
+//        Pout<< "Writing points inbetween three or more faces on same cell to "
+//            << multiStr.name() << endl;
+//        // Check whether inbetween patch faces on same cell
+//        Map<labelList> cellToFaces;
+//        forAll(pointNormals, patchPointi)
+//        {
+//            const labelList& pFaces = pointFaces[patchPointi];
+//
+//            cellToFaces.clear();
+//            forAll(pFaces, pFacei)
+//            {
+//                const label patchFacei = pFaces[pFacei];
+//                const label meshFacei = pp.addressing()[patchFacei];
+//                const label celli = mesh.faceOwner()[meshFacei];
+//                Map<labelList>::iterator faceFnd = cellToFaces.find(celli);
+//                if (faceFnd.found())
+//                {
+//                    labelList& faces = faceFnd();
+//                    if (!faces.found(patchFacei))
+//                    {
+//                        faces.append(patchFacei);
+//                    }
+//                }
+//                else
+//                {
+//                    cellToFaces.insert(celli, labelList(1, patchFacei));
+//                }
+//            }
+//
+//            forAllConstIter(Map<labelList>, cellToFaces, iter)
+//            {
+//                if (iter().size() == 2)
+//                {
+//                    twoStr.write(pp.localPoints()[patchPointi]);
+//                }
+//                else if (iter().size() > 2)
+//                {
+//                    multiStr.write(pp.localPoints()[patchPointi]);
+//
+//                    const scalar ratio =
+//                        layerParameters::finalLayerThicknessRatio
+//                        (
+//                            patchNLayers[patchPointi],
+//                            expansionRatio[patchPointi]
+//                        );
+//                    // Get thickness of cell next to bulk
+//                    const vector finalDisp
+//                    (
+//                        ratio*patchDisp[patchPointi]
+//                    );
+//
+//                    //Pout<< "** point:" << pp.localPoints()[patchPointi]
+//                    //    << " on cell:" << iter.key()
+//                    //    << " faces:" << iter()
+//                    //    << " displacement was:" << patchDisp[patchPointi]
+//                    //    << " ratio:" << ratio
+//                    //    << " finalDispl:" << finalDisp;
+//
+//                    // Half this thickness
+//                    patchDisp[patchPointi] -= 0.8*finalDisp;
+//
+//                    //Pout<< " new displacement:"
+//                    //    << patchDisp[patchPointi] << endl;
+//                }
+//            }
+//        }
+//
+//        Pout<< "Written " << multiStr.nVertices()
+//            << " points inbetween three or more faces on same cell to "
+//            << multiStr.name() << endl;
+//    }
+////XXXXXXXX
+
+
     // Check if no extrude possible.
     forAll(pointNormals, patchPointi)
     {
@@ -3724,7 +3914,7 @@ void Foam::snappyLayerDriver::addLayers
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
         // It is hard to calculate some length scale if not in relative
         // mode so disable this check.
-        if (layerParams.relativeSizes())
+        if (!layerParams.relativeSizes().found(false))
         {
             // Undistorted edge length
             const scalar edge0Len =
@@ -3735,6 +3925,7 @@ void Foam::snappyLayerDriver::addLayers
             (
                 *pp,
                 layerParams.maxFaceThicknessRatio(),
+                layerParams.relativeSizes(),
                 edge0Len,
                 cellLevel,
 
@@ -3948,6 +4139,8 @@ void Foam::snappyLayerDriver::addLayers
                 *pp,
                 thickness,
                 minThickness,
+                expansionRatio,
+
                 patchDisp,
                 patchNLayers,
                 extrudeStatus
@@ -4077,7 +4270,7 @@ void Foam::snappyLayerDriver::addLayers
 
             forAll(nPatchPointLayers, i)
             {
-                scalar ratio = layerParams.finalLayerThicknessRatio
+                scalar ratio = layerParameters::finalLayerThicknessRatio
                 (
                     nPatchPointLayers[i],
                     expansionRatio[i]
@@ -4133,7 +4326,7 @@ void Foam::snappyLayerDriver::addLayers
                     mesh.name(),
                     static_cast<polyMesh&>(mesh).instance(),
                     mesh.time(),  // register with runTime
-                    IOobject::NO_READ,
+                    IOobject::READ_IF_PRESENT,  // read fv* if present
                     static_cast<polyMesh&>(mesh).writeOpt()
                 ),              // io params from original mesh but new name
                 mesh,           // original mesh

src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.H

@@ -235,6 +235,7 @@ private:
                 (
                     const indirectPrimitivePatch& pp,
                     const scalar faceRatio,
+                    const boolList& relativeSizes,
                     const scalar edge0Len,
                     const labelList& cellLevel,
                     pointField& patchDisp,
@@ -328,6 +329,8 @@ private:
                     const indirectPrimitivePatch& pp,
                     const scalarField& thickness,
                     const scalarField& minThickness,
+                    const scalarField& expansionRatio,
+
                     pointField& patchDisp,
                     labelList& patchNLayers,
                     List<extrudeMode>& extrudeStatus

src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C

@@ -1476,6 +1476,261 @@ Foam::label Foam::snappyRefineDriver::refinementInterfaceRefine
     return iter;
 }
 
+bool Foam::snappyRefineDriver::usesHigherLevel
+(
+    const labelUList& boundaryPointLevel,
+    const labelUList& f,
+    const label cLevel
+) const
+{
+    for (const label pointi : f)
+    {
+        if (boundaryPointLevel[pointi] > cLevel)
+        {
+            return true;
+        }
+    }
+    return false;
+}
+
+
+Foam::label Foam::snappyRefineDriver::boundaryRefinementInterfaceRefine
+(
+    const refinementParameters& refineParams,
+    const label maxIter
+)
+{
+    if (refineParams.minRefineCells() == -1)
+    {
+        // Special setting to be able to restart shm on meshes with inconsistent
+        // cellLevel/pointLevel
+        return 0;
+    }
+
+    if (dryRun_)
+    {
+        return 0;
+    }
+
+    addProfiling(interface, "snappyHexMesh::refine::transition");
+    const fvMesh& mesh = meshRefiner_.mesh();
+
+    label iter = 0;
+
+    if (refineParams.interfaceRefine())
+    {
+        for (;iter < maxIter; iter++)
+        {
+            Info<< nl
+                << "Boundary refinement iteration " << iter << nl
+                << "-------------------------------" << nl
+                << endl;
+
+            const labelList& surfaceIndex = meshRefiner_.surfaceIndex();
+            const hexRef8& cutter = meshRefiner_.meshCutter();
+            const labelList& cellLevel = cutter.cellLevel();
+            const faceList& faces = mesh.faces();
+            const cellList& cells = mesh.cells();
+
+
+            // Determine cells to refine
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+            // Point/face on boundary
+            bitSet isBoundaryPoint(mesh.nPoints());
+            bitSet isBoundaryFace(mesh.nFaces());
+            {
+                forAll(surfaceIndex, facei)
+                {
+                    if (surfaceIndex[facei] != -1)
+                    {
+                        isBoundaryFace.set(facei);
+                        isBoundaryPoint.set(faces[facei]);
+                    }
+                }
+                const labelList meshPatchIDs(meshRefiner_.meshedPatches());
+                for (const label patchi : meshPatchIDs)
+                {
+                    const polyPatch& pp = mesh.boundaryMesh()[patchi];
+                    forAll(pp, i)
+                    {
+                        isBoundaryFace.set(pp.start()+i);
+                        isBoundaryPoint.set(pp[i]);
+                    }
+                }
+
+                syncTools::syncPointList
+                (
+                    mesh,
+                    isBoundaryPoint,
+                    orEqOp<unsigned int>(),
+                    0
+                );
+            }
+
+            // Mark max boundary face level onto boundary points. All points
+            // not on a boundary face stay 0.
+            labelList boundaryPointLevel(mesh.nPoints(), 0);
+            {
+                forAll(cells, celli)
+                {
+                    const cell& cFaces = cells[celli];
+                    const label cLevel = cellLevel[celli];
+
+                    for (const label facei : cFaces)
+                    {
+                        if (isBoundaryFace(facei))
+                        {
+                            const face& f = faces[facei];
+                            for (const label pointi : f)
+                            {
+                                boundaryPointLevel[pointi] =
+                                max
+                                (
+                                    boundaryPointLevel[pointi],
+                                    cLevel
+                                );
+                            }
+                        }
+                    }
+                }
+
+                syncTools::syncPointList
+                (
+                    mesh,
+                    boundaryPointLevel,
+                    maxEqOp<label>(),
+                    0
+                );
+            }
+
+
+            // Detect cells with a point but not face on the boundary
+            labelList candidateCells;
+            {
+                const cellList& cells = mesh.cells();
+
+                cellSet candidateCellSet
+                (
+                    mesh,
+                    "candidateCells",
+                    cells.size()/100
+                );
+
+                forAll(cells, celli)
+                {
+                    const cell& cFaces = cells[celli];
+                    const label cLevel = cellLevel[celli];
+
+                    bool isBoundaryCell = false;
+                    for (const label facei : cFaces)
+                    {
+                        if (isBoundaryFace(facei))
+                        {
+                            isBoundaryCell = true;
+                            break;
+                        }
+                    }
+
+                    if (!isBoundaryCell)
+                    {
+                        for (const label facei : cFaces)
+                        {
+                            const face& f = mesh.faces()[facei];
+                            if (usesHigherLevel(boundaryPointLevel, f, cLevel))
+                            {
+                                candidateCellSet.insert(celli);
+                                break;
+                            }
+                        }
+                    }
+                }
+                if (debug&meshRefinement::MESH)
+                {
+                    Pout<< "Dumping " << candidateCellSet.size()
+                        << " cells to cellSet candidateCellSet." << endl;
+                    candidateCellSet.instance() = meshRefiner_.timeName();
+                    candidateCellSet.write();
+                }
+                candidateCells = candidateCellSet.toc();
+            }
+
+            labelList cellsToRefine
+            (
+                meshRefiner_.meshCutter().consistentRefinement
+                (
+                    candidateCells,
+                    true
+                )
+            );
+            Info<< "Determined cells to refine in = "
+                << mesh.time().cpuTimeIncrement() << " s" << endl;
+
+
+            label nCellsToRefine = cellsToRefine.size();
+            reduce(nCellsToRefine, sumOp<label>());
+
+            Info<< "Selected for refinement : " << nCellsToRefine
+                << " cells (out of " << mesh.globalData().nTotalCells()
+                << ')' << endl;
+
+            // Stop when no cells to refine. After a few iterations check if too
+            // few cells
+            if
+            (
+                nCellsToRefine == 0
+                // || (
+                //        iter >= 1
+                //     && nCellsToRefine <= refineParams.minRefineCells()
+                //    )
+            )
+            {
+                Info<< "Stopping refining since too few cells selected."
+                    << nl << endl;
+                break;
+            }
+
+
+            if (debug)
+            {
+                const_cast<Time&>(mesh.time())++;
+            }
+
+
+            if
+            (
+                returnReduce
+                (
+                    (mesh.nCells() >= refineParams.maxLocalCells()),
+                    orOp<bool>()
+                )
+            )
+            {
+                meshRefiner_.balanceAndRefine
+                (
+                    "boundary cell refinement iteration " + name(iter),
+                    decomposer_,
+                    distributor_,
+                    cellsToRefine,
+                    refineParams.maxLoadUnbalance()
+                );
+            }
+            else
+            {
+                meshRefiner_.refineAndBalance
+                (
+                    "boundary cell refinement iteration " + name(iter),
+                    decomposer_,
+                    distributor_,
+                    cellsToRefine,
+                    refineParams.maxLoadUnbalance()
+                );
+            }
+        }
+    }
+    return iter;
+}
+
 
 void Foam::snappyRefineDriver::removeInsideCells
 (
@@ -3195,6 +3450,18 @@ void Foam::snappyRefineDriver::doRefine
         motionDict
     );
 
+
+    //- Commented out for now since causes zoning errors (sigsegv) on
+    //  case with faceZones. TBD.
+    //// Refine any cells are point/edge connected to the boundary and have a
+    //// lower refinement level than the neighbouring cells
+    //boundaryRefinementInterfaceRefine
+    //(
+    //    refineParams,
+    //    10      // maxIter
+    //);
+
+
     // Mesh is at its finest. Do optional zoning (cellZones and faceZones)
     wordPairHashTable zonesToFaceZone;
     zonify(refineParams, zonesToFaceZone);

src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.H

@@ -165,6 +165,21 @@ class snappyRefineDriver
             const label maxIter
         );
 
+        //- Helper: see if any element in f has higher level than cLevel
+        bool usesHigherLevel
+        (
+            const labelUList& boundaryPointLevel,
+            const labelUList& f,
+            const label cLevel
+        ) const;
+
+        //- Refine cells with a point/edge but not face on the boundary
+        label boundaryRefinementInterfaceRefine
+        (
+            const refinementParameters& refineParams,
+            const label maxIter
+        );
+
         //- Remove all cells within intersected region
         void removeInsideCells
         (

src/meshTools/polyTopoChange/polyTopoChange.H

@@ -617,7 +617,8 @@ public:
                 const bool orderPoints = false
             );
 
-            //- Create new mesh with old mesh patches
+            //- Create new mesh with old mesh patches. Additional dictionaries
+            //  (fv* etc) read according to IO flags
             template<class Type>
             autoPtr<mapPolyMesh> makeMesh
             (

src/meshTools/polyTopoChange/polyTopoChangeTemplates.C

@@ -182,14 +182,14 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::makeMesh
     // Create the mesh
     // ~~~~~~~~~~~~~~~
 
-    IOobject noReadIO(io);
-    noReadIO.readOpt() = IOobject::NO_READ;
-    noReadIO.writeOpt() = IOobject::AUTO_WRITE;
+    //IOobject noReadIO(io);
+    //noReadIO.readOpt() = IOobject::NO_READ;
+    //noReadIO.writeOpt() = IOobject::AUTO_WRITE;
     newMeshPtr.reset
     (
         new Type
         (
-            noReadIO,
+            io,         //noReadIO
             std::move(newPoints),
             std::move(faces_),
             std::move(faceOwner_),

tutorials/mesh/snappyHexMesh/airfoilWithLayers/Allclean

@@ -0,0 +1,9 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/CleanFunctions      # Tutorial clean functions
+#------------------------------------------------------------------------------
+
+cleanCase0
+rm -rf constant/extendedFeatureEdgeMesh
+
+#------------------------------------------------------------------------------

tutorials/mesh/snappyHexMesh/airfoilWithLayers/Allrun

@@ -0,0 +1,22 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+runApplication blockMesh
+runApplication surfaceFeatureExtract
+
+mkdir -p 0
+
+# Run with basic
+foamDictionary -entry geometry.type -set basic system/fvSchemes
+runApplication -s basic snappyHexMesh
+runApplication -s basic checkMesh -writeAllFields
+foamListTimes -rm
+
+# Run with highAspectRatio
+foamDictionary -entry geometry.type -set highAspectRatio system/fvSchemes
+runApplication -s highAspectRatio snappyHexMesh
+runApplication -s highAspectRatio checkMesh -writeAllFields
+
+#------------------------------------------------------------------------------

tutorials/mesh/snappyHexMesh/airfoilWithLayers/constant/transportProperties

@@ -0,0 +1,32 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2006                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Air temperature at standard atmospheric pressure = 540 R = 300 K
+// See README.md
+// c = 347.336 [m/s]
+// Ma = UInf/c = 0.15
+// UInf = 52.1004 [m/s]
+// ReChord = 6e6 (per chord)
+// Chord = 1 [m]
+// ReChord = UInf*Chord/nuFluid
+// nuFluid = UInf*Chord/ReChord = 8.6834e-06 [m2/s]
+
+transportModel  Newtonian;
+
+nu              8.6834e-06;
+
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/airfoilWithLayers/constant/triSurface/internalFace.stl

@@ -0,0 +1,72 @@
+solid PART_6
+  facet normal -0.000000000e+00 -1.767905329e-06 -1.000000000e+00
+    outer loop
+      vertex   2.000000000e+00 -1.000000000e+00 5.352202629e-08
+      vertex   1.000000000e+00 -1.000000000e+00 5.352202629e-08
+      vertex   1.000000125e+00 0.000000000e+00 -1.714383303e-06
+    endloop
+  endfacet
+  facet normal -2.117582368e-22 -1.767905329e-06 -1.000000000e+00
+    outer loop
+      vertex   1.000000125e+00 0.000000000e+00 -1.714383303e-06
+      vertex   2.000000125e+00 0.000000000e+00 -1.714383303e-06
+      vertex   2.000000000e+00 -1.000000000e+00 5.352202629e-08
+    endloop
+  endfacet
+  facet normal -0.000000000e+00 -1.767905329e-06 -1.000000000e+00
+    outer loop
+      vertex   3.000000000e+00 -1.000000000e+00 5.352202629e-08
+      vertex   2.000000000e+00 -1.000000000e+00 5.352202629e-08
+      vertex   2.000000125e+00 0.000000000e+00 -1.714383303e-06
+    endloop
+  endfacet
+  facet normal 2.117582368e-22 -1.767905329e-06 -1.000000000e+00
+    outer loop
+      vertex   2.000000125e+00 0.000000000e+00 -1.714383303e-06
+      vertex   3.000000125e+00 0.000000000e+00 -1.714383303e-06
+      vertex   3.000000000e+00 -1.000000000e+00 5.352202629e-08
+    endloop
+  endfacet
+  facet normal -0.000000000e+00 -1.767905329e-06 -1.000000000e+00
+    outer loop
+      vertex   4.000000000e+00 -1.000000000e+00 5.352202629e-08
+      vertex   3.000000000e+00 -1.000000000e+00 5.352202629e-08
+      vertex   3.000000125e+00 0.000000000e+00 -1.714383303e-06
+    endloop
+  endfacet
+  facet normal 0.000000000e+00 -1.767905329e-06 -1.000000000e+00
+    outer loop
+      vertex   3.000000125e+00 0.000000000e+00 -1.714383303e-06
+      vertex   4.000000125e+00 0.000000000e+00 -1.714383303e-06
+      vertex   4.000000000e+00 -1.000000000e+00 5.352202629e-08
+    endloop
+  endfacet
+  facet normal 2.117582368e-22 -1.767905329e-06 -1.000000000e+00
+    outer loop
+      vertex   5.000000125e+00 0.000000000e+00 -1.714383303e-06
+      vertex   6.000000125e+00 0.000000000e+00 -1.714383303e-06
+      vertex   6.000000000e+00 -1.000000000e+00 5.352202629e-08
+    endloop
+  endfacet
+  facet normal 0.000000000e+00 -1.767905329e-06 -1.000000000e+00
+    outer loop
+      vertex   5.000000125e+00 0.000000000e+00 -1.714383303e-06
+      vertex   6.000000000e+00 -1.000000000e+00 5.352202629e-08
+      vertex   5.000000000e+00 -1.000000000e+00 5.352202629e-08
+    endloop
+  endfacet
+  facet normal -2.117582368e-22 -1.767905329e-06 -1.000000000e+00
+    outer loop
+      vertex   4.000000125e+00 0.000000000e+00 -1.714383303e-06
+      vertex   5.000000125e+00 0.000000000e+00 -1.714383303e-06
+      vertex   5.000000000e+00 -1.000000000e+00 5.352202629e-08
+    endloop
+  endfacet
+  facet normal 0.000000000e+00 -1.767905329e-06 -1.000000000e+00
+    outer loop
+      vertex   4.000000000e+00 -1.000000000e+00 5.352202629e-08
+      vertex   4.000000125e+00 0.000000000e+00 -1.714383303e-06
+      vertex   5.000000000e+00 -1.000000000e+00 5.352202629e-08
+    endloop
+  endfacet
+endsolid

tutorials/mesh/snappyHexMesh/airfoilWithLayers/system/blockMeshDict

@@ -0,0 +1,56 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 1;
+
+vertices
+(
+    (-1 -1 -1)
+    ( 2 -1 -1)
+    ( 2  0 -1)
+    (-1  0 -1)
+    (-1 -1  1)
+    ( 2 -1  1)
+    ( 2  0  1)
+    (-1  0  1)
+);
+
+blocks
+(
+//    hex (0 1 2 3 4 5 6 7) (15 10 10) simpleGrading (1 1 1)
+    hex (0 1 2 3 4 5 6 7) (15 5 10) simpleGrading (1 1 1)
+);
+
+edges
+();
+boundary
+(
+    airFlow
+    {
+        type patch;
+        faces
+        (
+            (3 7 6 2)
+            (1 5 4 0)  //back
+            (2 6 5 1)  //outlet
+            (0 4 7 3)  //inlet
+            (0 3 2 1)  //lowerWall
+            (4 5 6 7)  //upperWall
+        );
+    }
+
+);

tutorials/mesh/snappyHexMesh/airfoilWithLayers/system/controlDict

@@ -0,0 +1,55 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2006                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//DebugSwitches
+//{
+//    fvGeometryScheme    1;
+//    highAspectRatio     1;
+//    basic               1;
+//}
+
+application     simpleFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         15000;
+
+deltaT          1;
+
+writeControl    timeStep;
+
+writeInterval   5000;
+
+purgeWrite      2;
+
+writeFormat     binary;
+
+writePrecision  15;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   8;
+
+runTimeModifiable false;
+
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/airfoilWithLayers/system/fvSchemes

@@ -0,0 +1,73 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2006                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+FoamFile
+{
+    version         2;
+    format          ascii;
+    class           dictionary;
+    object          fvSchemes;
+}
+
+ddtSchemes
+{
+    default         steadyState;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+    grad(p)         Gauss linear;
+    unlimitedGrad(U) Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+    div(phi,U)      bounded Gauss linearUpwindV unlimitedGrad(U);
+    turbulence      bounded Gauss limitedLinear 1;
+    div(phi,k)      bounded Gauss limitedLinear 1;
+    div(phi,omega)  bounded Gauss limitedLinear 1;
+    div(phi,nuTilda) bounded Gauss limitedLinear 1;
+    div(phi,epsilon) bounded Gauss limitedLinear 1;
+    div(phi,phit)   bounded Gauss limitedLinear 1;
+    div(phi,f)      bounded Gauss limitedLinear 1;
+    div(phi,gammaInt) bounded Gauss linearUpwind grad;
+    div(phi,ReThetat) bounded Gauss linearUpwind grad;
+    div((nuEff*dev2(T(grad(U))))) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear limited corrected 0.33;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         limited corrected 0.33;
+}
+
+wallDist
+{
+    method          meshWave;
+}
+
+geometry
+{
+    type            highAspectRatio;
+    minAspect       10;
+    maxAspect       100;
+}
+
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/airfoilWithLayers/system/fvSolution

@@ -0,0 +1,71 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2006                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p
+    {
+        solver          GAMG;
+        smoother        GaussSeidel;
+        tolerance       1e-06;
+        relTol          0.1;
+    }
+
+    "(U|k|omega|nuTilda|gammaInt|ReThetat)"
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-08;
+        relTol          0.1;
+        maxIter         50;
+    }
+
+    "(epsilon|phit)"
+    {
+         solver          PBiCGStab;
+         preconditioner  DILU;
+         tolerance       1e-8;
+         relTol          0;
+    }
+
+    f
+    {
+        solver          PBiCGStab;
+        preconditioner  DIC;
+        tolerance       1e-8;
+        relTol          0;
+    }
+}
+
+SIMPLE
+{
+    nNonOrthogonalCorrectors 0;
+    consistent          true;
+}
+
+relaxationFactors
+{
+    equations
+    {
+        ".*"            0.9;
+        "(gammaInt|ReThetat|k|nuTilda)"    0.8;
+        "(phit|f)"      0.7;
+        epsilon         0.5;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/airfoilWithLayers/system/meshQualityDict

@@ -0,0 +1,78 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      meshQualityDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//- Maximum non-orthogonality allowed. Set to 180 to disable.
+maxNonOrtho 65;
+
+//- Max skewness allowed. Set to <0 to disable.
+maxBoundarySkewness 20;
+maxInternalSkewness 4;
+
+//- Max concaveness allowed. Is angle (in degrees) below which concavity
+//  is allowed. 0 is straight face, <0 would be convex face.
+//  Set to 180 to disable.
+maxConcave 80;
+
+//- Minimum pyramid volume. Is absolute volume of cell pyramid.
+//  Set to a sensible fraction of the smallest cell volume expected.
+//  Set to very negative number (e.g. -1E30) to disable.
+minVol 1e-13;
+
+//- Minimum quality of the tet formed by the face-centre
+//  and variable base point minimum decomposition triangles and
+//  the cell centre. Set to very negative number (e.g. -1E30) to
+//  disable.
+//     <0 = inside out tet,
+//      0 = flat tet
+//      1 = regular tet
+minTetQuality 1e-15;
+
+//- Minimum face area. Set to <0 to disable.
+minArea -1;
+
+//- Minimum face twist. Set to <-1 to disable. dot product of face normal
+//  (itself the average of the triangle normals)
+//  and face centre triangles normal
+minTwist 0.02;
+
+//- Minimum normalised cell determinant. This is the determinant of all
+//  the areas of internal faces. It is a measure of how much of the
+//  outside area of the cell is to other cells. The idea is that if all
+//  outside faces of the cell are 'floating' (zeroGradient) the
+//  'fixedness' of the cell is determined by the area of the internal faces.
+//  1 = hex, <= 0 = folded or flattened illegal cell
+minDeterminant 0.001;
+
+//- Relative position of face in relation to cell centres (0.5 for orthogonal
+//  mesh) (0 -> 0.5)
+minFaceWeight 0.05;
+
+//- Volume ratio of neighbouring cells (0 -> 1)
+minVolRatio 0.01;
+
+//- Per triangle normal compared to that of preceding triangle. Must be >0
+//  for Fluent compatibility
+minTriangleTwist -1;
+
+
+//- If >0 : preserve cells with all points on the surface if the
+//  resulting volume after snapping (by approximation) is larger than
+//  minVolCollapseRatio times old volume (i.e. not collapsed to flat cell).
+//  If <0 : delete always.
+//minVolCollapseRatio 0.1;
+
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/airfoilWithLayers/system/snappyHexMeshDict

@@ -0,0 +1,753 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  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;
+
+
+// Optional: single region surfaces get patch names according to
+//           surface only. Multi-region surfaces get patch name
+//           surface "_ "region. Default is true
+// singleRegionName false;
+
+// Optional: avoid patch-face merging. Allows mesh to be used for
+//           refinement/unrefinement
+// mergePatchFaces    false; // default true
+
+// Optional: preserve all generated patches. Default is to remove
+//          zero-sized patches.
+// keepPatches 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
+{
+    refinement1
+    {
+        type searchableBox;
+        min (-1 -2 -1);
+        max ( 1  2  1);
+    }
+/*
+    // Shell for directional refinement
+    wakeBox
+    {
+        type searchableBox;
+        min (1.5 1 -0.5);
+        max (3.5 2 0.5);
+    }
+*/
+    internalFace.stl { name internalFace; type triSurfaceMesh;}
+
+    aerofoil.stl
+    {
+        name aerofoil;
+        type triSurfaceMesh;
+
+        //tolerance   1E-5;   // optional:non-default tolerance on intersections
+        //maxTreeDepth 10;    // optional:depth of octree. Decrease only in case
+                              // of memory limitations.
+
+        // Per region the patchname. If not provided will be <surface>_<region>.
+        // Note: this name cannot be used to identity this region in any
+        //       other part of this dictionary; it is only a name
+        //       for the combination of surface+region (which is only used
+        //       when creating patches)
+//        regions
+//        {
+//            secondSolid
+//            {
+//                name mySecondPatch;
+//            }
+//        }
+    }
+
+/*
+    sphere2
+    {
+        type searchableSphere;
+        centre  (1.5 1.5 1.5);
+        radius  1.03;
+    }
+*/
+};
+
+
+// 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 1000000;
+
+    // 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 20000000;
+
+    // The surface refinement loop might spend lots of iterations refining just
+    // a few cells. This setting will cause refinement to stop if
+    // <= minRefinementCells cells are selected for refinement. Note: it will
+    // at least do one iteration unless
+    //  a: the number of cells to refine is 0
+    //  b: minRefinementCells = -1. This is a special value indicating
+    //     no refinement.
+    minRefinementCells 0;
+
+    // Allow a certain level of imbalance during refining
+    // (since balancing is quite expensive)
+    // Expressed as fraction of perfect balance (= overall number of cells /
+    // nProcs). 0=balance always.
+    maxLoadUnbalance 0.10;
+
+    // Number of buffer layers between different levels.
+    // 1 means normal 2:1 refinement restriction, larger means slower
+    // refinement.
+    nCellsBetweenLevels 4;
+
+
+    // Explicit feature edge refinement
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Specifies a level for any cell intersected by explicitly provided
+    // edges.
+    // This is a featureEdgeMesh, read from constant/triSurface for now.
+    // Specify 'levels' in the same way as the 'distance' mode in the
+    // refinementRegions (see below). The old specification
+    //      level   2;
+    // is equivalent to
+    //      levels  ((0 2));
+
+    features
+    (
+        {
+            file "aerofoil.eMesh";
+            level 0;        // Have level 2 for all cells intersected
+        //                      // by feature.
+        //    levels ((0.1 2)); // Have level 2 for all cells within
+        //                      // 0.1 of feature.
+        }
+    );
+
+
+    // 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
+    {
+        internalFace { level (0 0); faceZone internalFace; }//faceType baffle;}
+        aerofoil
+        {
+            // Surface-wise min and max refinement level
+            level (2 2);
+
+
+
+
+            //- Optional increment (on top of max level) in small gaps
+            //gapLevelIncrement 2;
+
+            //- Optional angle to detect small-large cell situation
+            //  perpendicular to the surface. Is the angle of face w.r.t.
+            //  the local surface normal. Use on flat(ish) surfaces only.
+            //  Otherwise leave out or set to negative number.
+            //perpendicularAngle 10;
+
+            //- Optional faceZone and (for closed surface) cellZone with
+            //  how to select the cells that are in the cellZone
+            //  (inside / outside / specified insidePoint)
+            //  The orientation of the faceZone is
+            //  - if on cellZone(s) : point out of (minimum) cellZone
+            //  - if freestanding   : oriented according to surface
+
+            //faceZone sphere;
+            //cellZone sphere;
+            //cellZoneInside inside;    // outside/insidePoint
+            //insidePoint    (1 1 1);   // if (cellZoneInside == insidePoint)
+
+            //- Optional specification of what to do with faceZone faces:
+            //      internal : keep them as internal faces (default)
+            //      baffle   : create baffles from them. This gives more
+            //                 freedom in mesh motion
+            //      boundary : create free-standing boundary faces (baffles
+            //                 but without the shared points)
+            //faceType baffle;
+        }
+    }
+
+    // Feature angle:
+    // - used if min and max refinement level of a surface differ
+    // - used if feature snapping (see snapControls below) is used
+    resolveFeatureAngle 30;
+
+    //- Optional increment (on top of max level) in small gaps
+    //gapLevelIncrement 2;
+
+
+    // Planar angle:
+    // - used to determine if surface normals
+    //   are roughly the same or opposite. Used
+    //      - in proximity refinement
+    //      - to decide when to merge free-standing baffles
+    //        (if e.g. running in surfaceSimplify mode set this to 180 to
+    //         merge all baffles)
+    //      - in snapping to avoid snapping to nearest on 'wrong' side
+    //        of thin gap
+    //
+    // If not specified same as resolveFeatureAngle
+    planarAngle 30;
+
+
+    // 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
+    //   increasing 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
+    {
+//        refinement1
+//        {
+//            mode inside;
+//            levels ((1.0 1));
+//        }
+/*
+        //sphere.stl
+        //{
+        //    mode inside;
+        //    levels ((1.0 4));
+        //    // Optional override of uniform refinement level such
+        //    // that in small gaps we're getting more cells.
+        //    // The specification is
+        //    //  - numGapCells : minimum number of cells in the gap
+        //    //                  (usually >3; lower than this might not
+        //    //                   resolve correctly)
+        //    //  - minLevel    : min refinement level at which to kick in
+        //    //  - maxLevel    : upper refinement level (to avoid refinement
+        //    //                  continuing on a single extraneous feature)
+        //    // All three settings can be overridden on a surface by
+        //    // surface basis in the refinementSurfaces section.
+        //    gapLevel (<numGapCells> <minLevel> <maxlevel>);
+        //    // Optional: when doing the gapLevel refinement directly remove
+        //    //           based on orientation w.r.t. gap. This limits the
+        //    //           amount of cells before doing the 'locationInMesh'
+        //    //           cell selection. Default is 'mixed' i.e. keep cells
+        //    //           whilst doing the gap-level refinement.
+        //    //gapMode inside;  // inside/outside/mixed
+        //}
+
+        //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)
+        //}
+
+        //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
+        //    smoothDirection     (1 0 0);
+        //    // Smoothing of expansion ratio
+        //    nSmoothExpansion    100;
+        //    // Smoothing of positions
+        //    nSmoothPosition     100;
+        //}
+*/
+    }
+
+
+
+    // Optionally limit refinement in geometric region. This limits all
+    // refinement (from features, refinementSurfaces, refinementRegions)
+    // in a given geometric region. The syntax is exactly the same as for the
+    // refinementRegions; the cell level now specifies the upper limit
+    // for any cell. (a special setting is cell level -1 which will remove
+    // any cells inside the region). Note that it does not override the
+    // refinement constraints given by the nCellsBetweenLevels setting.
+    limitRegions
+    {
+    }
+
+
+    // 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 location(s)InMesh is kept.
+    // NOTE: This point should never be on a face, always inside a cell, even
+    // after refinement.
+    //
+    // There are two different ways of specifying the regions to keep:
+    // 1. a single locationInMesh. This is the unzoned part of the mesh.
+    //    All the 'zoned' surfaces are marked as such
+    //    in the refinementSurfaces with the faceZone and cellZone keywords.
+    //    It is illegal to have the locationInMesh inside a surface for which
+    //    a cellZone is specified.
+    //
+    // or
+    //
+    // 2. multiple locationsInMesh, with per location the name of the cellZone.
+    //    This uses walking to determine zones and automatically creates
+    //    faceZones on the outside of cellZones. The special name 'none' is
+    //    used to indicate the unzoned/background part of the mesh.
+
+
+        // Ad 1. Specify a single location and how to treat faces inbetween
+        //       cellZones
+        locationInMesh (-0.457 -0.5 0.43);
+
+        // Whether any faceZones (as specified in the refinementSurfaces)
+        // are only on the boundary of corresponding cellZones.
+        // Not used if there are no faceZones. The behaviour has changed
+        // with respect to previous versions:
+        //  true  : all intersections with surface are put in faceZone
+        //          (same behaviour as before)
+        //  false : depending on the type of surface intersected:
+        //      - if intersecting surface has faceZone only (so no cellZone)
+        //        leave in faceZone (so behave as if set to true) (= changed
+        //        behaviour)
+        //      - if intersecting surface has faceZone and cellZone
+        //        remove if inbetween same cellZone or if on boundary
+        //        (same behaviour as before)
+        allowFreeStandingZoneFaces true;
+
+
+
+        // 2. Specify multiple locations with optional cellZones for the
+        //    regions (use cellZone "none" to specify the unzoned cells)
+        //    FaceZones are automatically constructed from the
+        //    names of the cellZones: <cellZoneA> _to_ <cellZoneB>
+        //    where the cellZoneA is the lowest numbered cellZone (non-cellZone
+        //    cells are in a special region called "none" which is always
+        //    last).
+
+
+
+
+    // Optional locations that should not be reachable from
+    // location(s)InMesh
+//    locationsOutsideMesh ((100 100 100));
+
+    // Optional: do not remove cells likely to give snapping problems
+    // handleSnapProblems false;
+
+    // Optional: switch off topological test for cells to-be-squashed
+    //           and use geometric test instead
+    //useTopologicalSnapDetection false;
+
+    // Optional: do not refine surface cells with opposite faces of
+    //           differing refinement levels
+    //interfaceRefine false;
+
+    // Optional: use an erosion instead of region assignment to allocate
+    //           left-over cells to the background region (i.e. make cellZones
+    //           consistent with the intersections of the surface).
+    //           Erosion is specified as a number of erosion iterations.
+    //           Erosion has less chance of bleeding and changing the zone
+    //           for a complete region.
+    //nCellZoneErodeIter 2;
+}
+
+// Settings for the snapping.
+snapControls
+{
+    // Number of patch smoothing iterations before finding correspondence
+    // to surface
+    nSmoothPatch 3;
+
+    // Optional: number of smoothing iterations for internal points on
+    // refinement interfaces. This will reduce non-orthogonality on
+    // refinement interfaces.
+    //nSmoothInternal $nSmoothPatch;
+
+    // Maximum relative distance for points to be attracted by surface.
+    // True distance is this factor times local maximum edge length.
+    tolerance 1.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;
+
+    // (wip) disable snapping to opposite near surfaces (revert to 22x
+    //  behaviour)
+    // detectNearSurfacesSnap false;
+
+
+    // Feature snapping
+
+        // Number of feature edge snapping iterations.
+        // Leave out altogether to disable.
+        nFeatureSnapIter 10;
+
+        // Detect (geometric only) features by sampling the surface
+        // (default=false).
+        implicitFeatureSnap false;
+
+        // Use castellatedMeshControls::features (default = true)
+        explicitFeatureSnap true;
+
+        // Detect features between multiple surfaces
+        // (only for explicitFeatureSnap, default = false)
+        multiRegionFeatureSnap false;
+
+
+        //- When to run face splitting (never at first iteration, always
+        //  at last iteration). Is interval. Default -1 (disabled)
+        //nFaceSplitInterval 5;
+
+
+        // (wip) Optional for explicit feature snapping:
+        //- Detect baffle edges. Default is true.
+        //detectBaffles false;
+        //- On any faces where points are on multiple regions (see
+        //  multiRegionFeatureSnap) have the other points follow these points
+        //  instead of having their own independent movement, i.e. have snapping
+        //  to multi-region edges/points take priority. This might aid snapping
+        //  to sharp edges that are also region edges. The default is false.
+        //releasePoints true;
+        //- Walk along feature edges, adding missing ones. Default is true.
+        //stringFeatures false;
+        //- If diagonal attraction also attract other face points. Default is
+        //  false
+        //avoidDiagonal true;
+        //- When splitting what concave faces to leave intact. Default is 45
+        //  degrees.
+        //concaveAngle 30;
+        //- When splitting the minimum area ratio of faces. If face split
+        //  causes ratio of area less than this do not split. Default is 0.3
+        //minAreaRatio 0.3;
+        //- Attract points only to the surface they originate from. Default
+        //  false. This can improve snapping of intersecting surfaces.
+          strictRegionSnap true;
+}
+
+// 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;
+
+    // Layer thickness specification. This can be specified in one of following
+    // ways:
+    // - expansionRatio and finalLayerThickness (cell nearest internal mesh)
+    // - expansionRatio and firstLayerThickness (cell on surface)
+    // - overall thickness and firstLayerThickness
+    // - overall thickness and finalLayerThickness
+    // - overall thickness and expansionRatio
+    //
+    // Note: the mode thus selected is global, i.e. one cannot override the
+    //       mode on a per-patch basis (only the values can be overridden)
+
+        // Expansion factor for layer mesh
+        expansionRatio 1.5;
+
+        // Wanted thickness of the layer furthest away from the wall.
+        // If relativeSizes this is relative to undistorted size of cell
+        // outside layer.
+        finalLayerThickness 0.3;
+
+        // Wanted thickness of the layer next to the wall.
+        // If relativeSizes this is relative to undistorted size of cell
+        // outside layer.
+        //firstLayerThickness 0.3;
+
+        // Wanted overall thickness of layers.
+        // If relativeSizes this is relative to undistorted size of cell
+        // outside layer.
+        //thickness 0.5
+
+
+    // Minimum overall thickness of total layers. If for any reason layer
+    // cannot be above minThickness do not add layer.
+    // If relativeSizes this is relative to undistorted size of cell
+    // outside layer..
+    minThickness 0.1;
+
+
+    // Per final patch or faceZone (so not geometry!) the layer information
+    // Note: This behaviour changed after 21x. Any non-mentioned patches
+    //       now slide unless:
+    //          - nSurfaceLayers is explicitly mentioned to be 0.
+    //          - angle to nearest surface < slipFeatureAngle (see below)
+    layers
+    {
+        "internalFace.*" {nSurfaceLayers 20; }
+        aerofoil
+        {
+            nSurfaceLayers 20;
+
+        }
+    }
+
+    // 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 1.7.x! (didn't do anything in 1.7.x)
+    nGrow -1;
+
+    // Advanced settings
+
+
+    // Static analysis of starting mesh
+
+        // When not to extrude surface. 0 is flat surface, 90 is when two faces
+        // are perpendicular. Note: was not working correctly < 1806
+        featureAngle 180;
+
+        // When to merge patch faces. Default is featureAngle. Useful when
+        // featureAngle is large.
+        //mergePatchFacesAngle 45;
+
+        // Stop layer growth on highly warped cells
+        maxFaceThicknessRatio 1000;//0.5;
+
+
+    // Patch displacement
+
+        // Number of smoothing iterations of surface normals
+        nSmoothSurfaceNormals 1;
+
+        // Smooth layer thickness over surface patches
+        nSmoothThickness 10;
+
+
+
+    // Choice of mesh shrinking algorithm
+
+        // Optional mesh shrinking algorithm (default is displacementMedialAxis)
+        // The displacementMotionSolver is a wrapper around the displacement
+        // motion solvers. It needs specification of the solver to use and
+        // its control dictionary.
+        //meshShrinker displacementMotionSolver;
+        //solver displacementLaplacian;
+        //displacementLaplacianCoeffs
+        //{
+        //    diffusivity quadratic inverseDistance
+        //    (
+        //        sphere.stl_firstSolid
+        //        maxY
+        //    );
+        //}
+        // Note that e.g. displacementLaplacian needs entries in
+        // fvSchemes, fvSolution. Also specify a minIter > 1 when solving
+        // cellDisplacement since otherwise solution might not be sufficiently
+        // accurate on points.
+
+
+    // Medial axis analysis (for use with default displacementMedialAxis)
+
+        // Angle used to pick up medial axis points
+        // Note: changed(corrected) w.r.t 1.7.x! 90 degrees corresponds to 130
+        // in 1.7.x.
+        minMedialAxisAngle 90;
+
+        // Reduce layer growth where ratio thickness to medial
+        // distance is large
+        maxThicknessToMedialRatio 0.3;
+
+        // Number of smoothing iterations of interior mesh movement direction
+        nSmoothNormals 3;
+
+        // Optional: limit the number of steps walking away from the surface.
+        // Default is unlimited.
+        //nMedialAxisIter 10;
+
+        // Optional: smooth displacement after medial axis determination.
+        // default is 0.
+        //nSmoothDisplacement 90;
+
+        // (wip)Optional: do not extrude any point where
+        //   (false) : all surrounding faces are not fully extruded
+        //   (true)  : all surrounding points are not extruded
+        // Default is false.
+        //detectExtrusionIsland true;
+
+        // Optional: do not extrude around sharp edges if both faces are not
+        // fully extruded i.e. if one of the faces on either side would
+        // become a wedge.
+        // Default is 0.5*featureAngle. Set to -180 always attempt extrusion
+        //layerTerminationAngle 25;
+
+        // Optional: disable shrinking of edges that have one (or two) points
+        // on an extruded patch.
+        // Default is false to enable single/two cell thick channels to still
+        // have layers. In <=1806 this was true by default. On larger gaps it
+        // should have no effect.
+        //disableWallEdges true;
+
+        // Optional: at non-patched sides allow mesh to slip if extrusion
+        // direction makes angle larger than slipFeatureAngle. Default is
+        // 0.5*featureAngle.
+        slipFeatureAngle 10;
+
+        // Maximum number of snapping relaxation iterations. Should stop
+        // before upon reaching a correct mesh.
+        nRelaxIter 5;
+
+
+    // Mesh shrinking
+
+        // Create buffer region for new layer terminations, i.e. gradually
+        // step down number of layers. Set to <0 to terminate layer in one go.
+        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;
+
+        // Max number of iterations after which relaxed meshQuality controls
+        // get used. Up to nRelaxedIter it uses the settings in
+        // meshQualityControls,
+        // after nRelaxedIter it uses the values in
+        // meshQualityControls::relaxed.
+        nRelaxedIter 0;
+
+        // Additional reporting: if there are just a few faces where there
+        // are mesh errors (after adding the layers) print their face centres.
+        // This helps in tracking down problematic mesh areas.
+        //additionalReporting true;
+}
+
+// Generic mesh quality settings. At any undoable phase these determine
+// where to undo.
+meshQualityControls
+{
+    // Specify mesh quality constraints in separate dictionary so can
+    // be reused (e.g. checkMesh -meshQuality)
+    #include "meshQualityDict"
+
+    minDeterminant 1e-8;
+
+    // Optional : some meshing phases allow usage of relaxed rules.
+    // See e.g. addLayersControls::nRelaxedIter.
+    relaxed
+    {
+        // Maximum non-orthogonality allowed. Set to 180 to disable.
+        maxNonOrtho 75;
+        minTetQuality -1e30;
+        minTwist    -1;
+    }
+
+
+    // Advanced
+
+        // Number of error distribution iterations
+        nSmoothScale 4;
+        // amount to scale back displacement at error points
+        errorReduction 0.75;
+}
+
+// Advanced
+
+//// Debug flags
+//debugFlags
+//(
+//    mesh            // write intermediate meshes
+//    intersections   // write current mesh intersections as .obj files
+//    featureSeeds    // write information about explicit feature edge
+//                    // refinement
+//    attraction      // write attraction as .obj files
+//    layerInfo       // write information about layers
+//);
+//
+//// Write flags
+//writeFlags
+//(
+//    scalarLevels    // write volScalarField with cellLevel for postprocessing
+//    layerSets       // write cellSets, faceSets of faces in layer
+//    layerFields     // write volScalarField for layer coverage
+//);
+
+
+//// Format for writing lines. E.g. leak path. Default is vtk format.
+//setFormat ensight;
+
+// Merge tolerance. Is fraction of overall bounding box of initial mesh.
+// Note: the write tolerance needs to be higher than this.
+mergeTolerance 1e-6;
+
+// ************************************************************************* //

tutorials/mesh/snappyHexMesh/airfoilWithLayers/system/surfaceFeatureExtractDict

@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v1912                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      surfaceFeatureExtractDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+aerofoil.stl
+{
+    // Extract raw features (none | extractFromFile | extractFromSurface)
+    extractionMethod    extractFromSurface;
+
+    // Mark edges whose adjacent surface normals are at an angle less
+    // than includedAngle as features
+    // - 0  : selects no edges
+    // - 180: selects all edges
+    includedAngle       120;
+
+curvature true;
+    // Do not mark region edges
+    geometricTestOnly   yes;
+
+    // Generate additional intersection features (none | self | region)
+    intersectionMethod  none;
+
+    // Tolerance for surface intersections
+    // tolerance           1e-3;
+
+// Output options:
+
+    // Write features to obj format for postprocessing
+    writeObj            yes;
+}
+
+
+
+
+// ************************************************************************* //