diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
index c62a340a02..508bdcb2d8 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
@@ -1342,6 +1342,17 @@ int main(int argc, char *argv[])
     const snapParameters snapParams(snapDict, dryRun);
 
 
+    Info<< "Setting refinement level of surface to be consistent"
+        << " with curvature." << endl;
+    surfaces.setCurvatureMinLevelFields
+    (
+        refineParams.curvature(),
+        meshRefiner.meshCutter().level0EdgeLength()
+    );
+    Info<< "Checked curvature refinement in = "
+        << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
+
+
 
     // Add all the cellZones and faceZones
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
diff --git a/etc/caseDicts/annotated/snappyHexMeshDict b/etc/caseDicts/annotated/snappyHexMeshDict
index e4d319f8f0..b66214a74e 100644
--- a/etc/caseDicts/annotated/snappyHexMeshDict
+++ b/etc/caseDicts/annotated/snappyHexMeshDict
@@ -230,6 +230,21 @@ castellatedMeshControls
             //  outside locations. Default is only after all refinement has
             //  been done.
             //leakLevel 10;
+
+            // Additional refinement for regions of high curvature. Expressed
+            // (bit similar to gapLevel) as:
+            // - number of cells per radius of curvature. (usually 1 is
+            //   good enough)
+            // - starting cell level? Not used at the moment.
+            // - maximum cell level. This can be smaller or larger than the
+            //   max 'surface' level
+            // - minumum curvature radius to ignore (expressed as a cell level).
+            //   This can be used to avoid detecting small sharp surface
+            //   features. Set to -1 to ignore.
+            //   (sometimes you want more refinement than sharp features since
+            //    these can be done with feature edge snapping (so can leave
+            //    level (0 0))
+            //curvatureLevel (10 0 10 -1);
         }
     }
 
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H
index fb0d7305d0..2e985707b7 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H
@@ -471,6 +471,19 @@ private:
                 label& nRefine
             ) const;
 
+            //- Refine cells containing triangles with high refinement level
+            //  (currently due to high curvature or being inside shell with
+            //   high level)
+            label markSurfaceFieldRefinement
+            (
+                const label nAllowRefine,
+                const labelList& neiLevel,
+                const pointField& neiCc,
+
+                labelList& refineCell,
+                label& nRefine
+            ) const;
+
             //- Helper: count number of normals1 that are in normals2
             label countMatches
             (
@@ -479,6 +492,18 @@ private:
                 const scalar tol = 1e-6
             ) const;
 
+            //- Detect if two intersection points are high curvature (w.r.t.
+            //  lengthscale
+            bool highCurvature
+            (
+                const scalar minCosAngle,
+                const scalar lengthScale,
+                const point& p0,
+                const vector& n0,
+                const point& p1,
+                const vector& n1
+            ) const;
+
             //- Mark cells for surface curvature based refinement. Marks if
             //  local curvature > curvature or if on different regions
             //  (markDifferingRegions)
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementGapRefine.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementGapRefine.C
index 01af062032..a4e26a976b 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementGapRefine.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementGapRefine.C
@@ -1833,4 +1833,120 @@ Foam::label Foam::meshRefinement::markSmallFeatureRefinement
 }
 
 
+////XXXXXXXX
+Foam::label Foam::meshRefinement::markSurfaceFieldRefinement
+(
+    const label nAllowRefine,
+    const labelList& neiLevel,
+    const pointField& neiCc,
+
+    labelList& refineCell,
+    label& nRefine
+) const
+{
+    const labelList& cellLevel = meshCutter_.cellLevel();
+    const labelList& surfaceIndices = surfaces_.surfaces();
+
+    label oldNRefine = nRefine;
+
+    //- Force calculation of tetBasePt
+    (void)mesh_.tetBasePtIs();
+    (void)mesh_.cellTree();
+    const indexedOctree<treeDataCell>& tree = mesh_.cellTree();
+
+
+    forAll(surfaceIndices, surfI)
+    {
+        label geomI = surfaceIndices[surfI];
+        const searchableSurface& geom = surfaces_.geometry()[geomI];
+
+        // Get the element index in a roundabout way. Problem is e.g.
+        // distributed surface where local indices differ from global
+        // ones (needed for getRegion call)
+
+        pointField ctrs;
+        labelList region;
+        labelList minLevelField;
+        {
+            // Representative local coordinates and bounding sphere
+            scalarField radiusSqr;
+            geom.boundingSpheres(ctrs, radiusSqr);
+
+            List<pointIndexHit> info;
+            geom.findNearest(ctrs, radiusSqr, info);
+
+            forAll(info, i)
+            {
+                if (!info[i].hit())
+                {
+                    FatalErrorInFunction
+                        << "fc:" << ctrs[i]
+                        << " radius:" << radiusSqr[i]
+                        << exit(FatalError);
+                }
+            }
+
+            geom.getRegion(info, region);
+            geom.getField(info, minLevelField);
+        }
+
+        if (minLevelField.size() != geom.size())
+        {
+            Pout<< "** no minLevelField" << endl;
+            continue;
+        }
+
+
+        label nOldRefine = 0;
+
+        forAll(ctrs, i)
+        {
+            label cellI = -1;
+            if (tree.nodes().size() && tree.bb().contains(ctrs[i]))
+            {
+                cellI = tree.findInside(ctrs[i]);
+            }
+
+            if
+            (
+                cellI != -1
+             && refineCell[cellI] == -1
+             && minLevelField[i] > cellLevel[cellI]
+            )
+            {
+                if
+                (
+                   !markForRefine
+                    (
+                        surfI,
+                        nAllowRefine,
+                        refineCell[cellI],
+                        nRefine
+                    )
+                )
+                {
+                    break;
+                }
+            }
+        }
+
+        Info<< "For surface " << geom.name() << " found "
+            << returnReduce(nRefine-nOldRefine, sumOp<label>())
+            << " cells containing cached refinement field" << endl;
+
+        if
+        (
+            returnReduce(nRefine, sumOp<label>())
+          > returnReduce(nAllowRefine, sumOp<label>())
+        )
+        {
+            Info<< "Reached refinement limit." << endl;
+        }
+    }
+
+    return returnReduce(nRefine-oldNRefine, sumOp<label>());
+}
+////XXXXXXXXX
+
+
 // ************************************************************************* //
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C
index 2344db2dd9..e1873dca63 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C
@@ -1109,6 +1109,71 @@ Foam::label Foam::meshRefinement::countMatches
 }
 
 
+//XXXXXX
+//bool Foam::meshRefinement::highCurvature
+//(
+//    const scalar minCosAngle,
+//    const point& p0,
+//    const vector& n0,
+//    const point& p1,
+//    const vector& n1
+//) const
+//{
+//    return ((n0&n1) < minCosAngle);
+//}
+bool Foam::meshRefinement::highCurvature
+(
+    const scalar minCosAngle,
+    const scalar lengthScale,
+    const point& p0,
+    const vector& n0,
+    const point& p1,
+    const vector& n1
+) const
+{
+    const scalar cosAngle = (n0&n1);
+
+    if (cosAngle < minCosAngle)
+    {
+        // Sharp feature, independent of intersection points
+        return true;
+    }
+    else if (cosAngle > 1-1e-6)
+    {
+        // Co-planar
+        return false;
+    }
+    else
+    {
+        // Calculate radius of curvature
+
+        vector axis(n0 ^ n1);
+        const plane pl0(p0, (n0 ^ axis));
+        const scalar r1 = pl0.normalIntersect(p1, n1);
+
+        //const point radiusPoint(p1+r1*n1);
+        //DebugVar(radiusPoint);
+        const plane pl1(p1, (n1 ^ axis));
+        const scalar r0 = pl1.normalIntersect(p0, n0);
+
+        //const point radiusPoint(p0+r0*n0);
+        //DebugVar(radiusPoint);
+        //- Take average radius. Bit ad hoc
+        const scalar radius = 0.5*(mag(r1)+mag(r0));
+
+        if (radius < lengthScale)
+        {
+            return true;
+        }
+        else
+        {
+            return false;
+        }
+    }
+}
+//XXXXX
+
+
 // Mark cells for surface curvature based refinement.
 Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
 (
@@ -1187,9 +1252,12 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
     // Test for all intersections (with surfaces of higher max level than
     // minLevel) and cache per cell the interesting inter
     labelListList cellSurfLevels(mesh_.nCells());
+    List<pointList> cellSurfLocations(mesh_.nCells());
     List<vectorList> cellSurfNormals(mesh_.nCells());
 
     {
+        // Per segment the intersection point of the surfaces
+        List<pointList> surfaceLocation;
         // Per segment the normals of the surfaces
         List<vectorList> surfaceNormal;
         // Per segment the list of levels of the surfaces
@@ -1205,6 +1273,7 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
             labelList(surfaces_.maxLevel().size(), 0),  // min level
             surfaces_.maxLevel(),
 
+            surfaceLocation,
             surfaceNormal,
             surfaceLevel
         );
@@ -1216,12 +1285,14 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
         labelList visitOrder;
         forAll(surfaceNormal, pointI)
         {
+            pointList& pLocations = surfaceLocation[pointI];
             vectorList& pNormals = surfaceNormal[pointI];
             labelList& pLevel = surfaceLevel[pointI];
 
-            sortedOrder(pNormals, visitOrder, normalLess(pNormals));
+            sortedOrder(pNormals, visitOrder, normalLess(pLocations));
 
-            pNormals = List<point>(pNormals, visitOrder);
+            pLocations = List<point>(pLocations, visitOrder);
+            pNormals = List<vector>(pNormals, visitOrder);
             pLevel = labelUIndList(pLevel, visitOrder);
         }
 
@@ -1236,6 +1307,7 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
             label faceI = testFaces[i];
             label own = mesh_.faceOwner()[faceI];
 
+            const pointList& fPoints = surfaceLocation[i];
             const vectorList& fNormals = surfaceNormal[i];
             const labelList& fLevels = surfaceLevel[i];
 
@@ -1244,6 +1316,7 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
                 if (fLevels[hitI] > cellLevel[own])
                 {
                     cellSurfLevels[own].append(fLevels[hitI]);
+                    cellSurfLocations[own].append(fPoints[hitI]);
                     cellSurfNormals[own].append(fNormals[hitI]);
                 }
 
@@ -1253,6 +1326,7 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
                     if (fLevels[hitI] > cellLevel[nei])
                     {
                         cellSurfLevels[nei].append(fLevels[hitI]);
+                        cellSurfLocations[nei].append(fPoints[hitI]);
                         cellSurfNormals[nei].append(fNormals[hitI]);
                     }
                 }
@@ -1266,7 +1340,7 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
     if (debug)
     {
         label nSet = 0;
-        label nNormals = 9;
+        label nNormals = 0;
         forAll(cellSurfNormals, cellI)
         {
             const vectorList& normals = cellSurfNormals[cellI];
@@ -1296,21 +1370,38 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
     for
     (
         label cellI = 0;
-        !reachedLimit && cellI < cellSurfNormals.size();
+        !reachedLimit && cellI < cellSurfLocations.size();
         cellI++
     )
     {
+        const pointList& points = cellSurfLocations[cellI];
         const vectorList& normals = cellSurfNormals[cellI];
         const labelList& levels = cellSurfLevels[cellI];
 
+        // Current cell size
+        const scalar cellSize =
+            meshCutter_.level0EdgeLength()/pow(2.0, cellLevel[cellI]);
+
         // n^2 comparison of all normals in a cell
         for (label i = 0; !reachedLimit && i < normals.size(); i++)
         {
             for (label j = i+1; !reachedLimit && j < normals.size(); j++)
             {
-                if ((normals[i] & normals[j]) < curvature)
+                //if ((normals[i] & normals[j]) < curvature)
+                if
+                (
+                    highCurvature
+                    (
+                        curvature,
+                        cellSize,
+                        points[i],
+                        normals[i],
+                        points[j],
+                        normals[j]
+                    )
+                )
                 {
-                    label maxLevel = max(levels[i], levels[j]);
+                    const label maxLevel = max(levels[i], levels[j]);
 
                     if (cellLevel[cellI] < maxLevel)
                     {
@@ -1358,8 +1449,10 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
         label own = mesh_.faceOwner()[faceI];
         label nei = mesh_.faceNeighbour()[faceI];
 
+        const pointList& ownPoints = cellSurfLocations[own];
         const vectorList& ownNormals = cellSurfNormals[own];
         const labelList& ownLevels = cellSurfLevels[own];
+        const pointList& neiPoints = cellSurfLocations[nei];
         const vectorList& neiNormals = cellSurfNormals[nei];
         const labelList& neiLevels = cellSurfLevels[nei];
 
@@ -1379,13 +1472,30 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
 
         if (!ownIsSubset && !neiIsSubset)
         {
+            // Current average cell size. Is min? max? average?
+            const scalar cellSize =
+                meshCutter_.level0EdgeLength()
+              / pow(2.0, min(cellLevel[own], cellLevel[nei]));
+
             // n^2 comparison of between ownNormals and neiNormals
             for (label i = 0; !reachedLimit &&  i < ownNormals.size(); i++)
             {
                 for (label j = 0; !reachedLimit && j < neiNormals.size(); j++)
                 {
                     // Have valid data on both sides. Check curvature.
-                    if ((ownNormals[i] & neiNormals[j]) < curvature)
+                    //if ((ownNormals[i] & neiNormals[j]) < curvature)
+                    if
+                    (
+                        highCurvature
+                        (
+                            curvature,
+                            cellSize,
+                            ownPoints[i],
+                            ownNormals[i],
+                            neiPoints[j],
+                            neiNormals[j]
+                        )
+                    )
                     {
                         // See which side to refine.
                         if (cellLevel[own] < ownLevels[i])
@@ -1441,7 +1551,11 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
     }
 
 
-    // Send over surface normal to neighbour cell.
+    // Send over surface point/normal to neighbour cell.
+//    labelListList neiSurfaceLevel;
+//    syncTools::swapBoundaryCellList(mesh_, cellSurfLevels, neiSurfaceLevel);
+    List<vectorList> neiSurfacePoints;
+    syncTools::swapBoundaryCellList(mesh_, cellSurfLocations, neiSurfacePoints);
     List<vectorList> neiSurfaceNormals;
     syncTools::swapBoundaryCellList(mesh_, cellSurfNormals, neiSurfaceNormals);
 
@@ -1456,9 +1570,13 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
         label own = mesh_.faceOwner()[faceI];
         label bFaceI = faceI - mesh_.nInternalFaces();
 
+        const pointList& ownPoints = cellSurfLocations[own];
         const vectorList& ownNormals = cellSurfNormals[own];
         const labelList& ownLevels = cellSurfLevels[own];
+
+        const pointList& neiPoints = neiSurfacePoints[bFaceI];
         const vectorList& neiNormals = neiSurfaceNormals[bFaceI];
+        //const labelList& neiLevels = neiSurfaceLevel[bFacei];
 
         // Special case: owner normals set is a subset of the neighbour
         // normals. Do not do curvature refinement since other cell's normals
@@ -1475,13 +1593,30 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
 
         if (!ownIsSubset && !neiIsSubset)
         {
+            // Current average cell size. Is min? max? average?
+            const scalar cellSize =
+                meshCutter_.level0EdgeLength()
+              / pow(2.0, min(cellLevel[own], neiLevel[bFaceI]));
+
             // n^2 comparison of between ownNormals and neiNormals
             for (label i = 0; !reachedLimit &&  i < ownNormals.size(); i++)
             {
                 for (label j = 0; !reachedLimit && j < neiNormals.size(); j++)
                 {
                     // Have valid data on both sides. Check curvature.
-                    if ((ownNormals[i] & neiNormals[j]) < curvature)
+                    //if ((ownNormals[i] & neiNormals[j]) < curvature)
+                    if
+                    (
+                        highCurvature
+                        (
+                            curvature,
+                            cellSize,
+                            ownPoints[i],
+                            ownNormals[i],
+                            neiPoints[j],
+                            neiNormals[j]
+                        )
+                    )
                     {
                         if (cellLevel[own] < ownLevels[i])
                         {
@@ -2205,25 +2340,46 @@ Foam::labelList Foam::meshRefinement::refineCandidates
         // Refinement based on features smaller than cell size
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-        if
-        (
-            smallFeatureRefinement
-         && (planarCos >= -1 && planarCos <= 1)
-         && max(shells_.maxGapLevel()) > 0
-        )
+        if (smallFeatureRefinement)
         {
-            label nGap = markSmallFeatureRefinement
+            label nGap = 0;
+            if
             (
-                planarCos,
-                nAllowRefine,
-                neiLevel,
-                neiCc,
+                planarCos >= -1
+             && planarCos <= 1
+             && max(shells_.maxGapLevel()) > 0
+            )
+            {
+                nGap = markSmallFeatureRefinement
+                (
+                    planarCos,
+                    nAllowRefine,
+                    neiLevel,
+                    neiCc,
 
-                refineCell,
-                nRefine
-            );
+                    refineCell,
+                    nRefine
+                );
+            }
             Info<< "Marked for refinement due to close opposite surfaces       "
                 << ": " << nGap << " cells."  << endl;
+
+            label nCurv = 0;
+            if (max(surfaces_.maxCurvatureLevel()) > 0)
+            {
+                nCurv = markSurfaceFieldRefinement
+                (
+                    nAllowRefine,
+                    neiLevel,
+                    neiCc,
+
+                    refineCell,
+                    nRefine
+                );
+                Info<< "Marked for refinement"
+                    << " due to curvature                     "
+                    << ": " << nCurv << " cells."  << endl;
+            }
         }
 
 
diff --git a/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.C b/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.C
index c98a71adab..1e0de1362d 100644
--- a/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.C
+++ b/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2015-2020 OpenCFD Ltd.
+    Copyright (C) 2015-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -38,6 +38,8 @@ License
 // For dictionary::get wrapper
 #include "meshRefinement.H"
 
+#include "OBJstream.H"
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 Foam::labelList Foam::refinementSurfaces::findHigherLevel
@@ -194,15 +196,14 @@ Foam::refinementSurfaces::refinementSurfaces
     labelList globalMaxLevel(surfI, Zero);
     labelList globalLevelIncr(surfI, Zero);
 
-    FixedList<label, 3> nullGapLevel;
-    nullGapLevel[0] = 0;
-    nullGapLevel[1] = 0;
-    nullGapLevel[2] = 0;
-
+    const FixedList<label, 3> nullGapLevel({0, 0, 0});
     List<FixedList<label, 3>> globalGapLevel(surfI);
     List<volumeType> globalGapMode(surfI);
     boolList globalGapSelf(surfI);
 
+    const FixedList<label, 4> nullCurvLevel({0, 0, 0, -1});
+    List<FixedList<label, 4>> globalCurvLevel(surfI);
+
     scalarField globalAngle(surfI, -GREAT);
     PtrList<dictionary> globalPatchInfo(surfI);
 
@@ -216,6 +217,7 @@ Foam::refinementSurfaces::refinementSurfaces
     List<Map<FixedList<label, 3>>> regionGapLevel(surfI);
     List<Map<volumeType>> regionGapMode(surfI);
     List<Map<bool>> regionGapSelf(surfI);
+    List<Map<FixedList<label, 4>>> regionCurvLevel(surfI);
     List<Map<scalar>> regionAngle(surfI);
     List<Map<autoPtr<dictionary>>> regionPatchInfo(surfI);
     List<Map<label>> regionBlockLevel(surfI);
@@ -307,6 +309,19 @@ Foam::refinementSurfaces::refinementSurfaces
             globalGapSelf[surfI] =
                 dict.getOrDefault<bool>("gapSelf", true);
 
+            globalCurvLevel[surfI] = nullCurvLevel;
+            if (dict.readIfPresent("curvatureLevel", globalCurvLevel[surfI]))
+            {
+                if (globalCurvLevel[surfI][0] <= 0)
+                {
+                    FatalIOErrorInFunction(dict)
+                        << "Illegal curvatureLevel specification for surface "
+                        << names_[surfI]
+                        << " : curvatureLevel:" << globalCurvLevel[surfI]
+                        << exit(FatalIOError);
+                }
+            }
+
             const searchableSurface& surface = allGeometry_[surfaces_[surfI]];
 
             // Surface zones
@@ -435,6 +450,20 @@ Foam::refinementSurfaces::refinementSurfaces
                             )
                         );
 
+                        FixedList<label, 4> curvSpec(nullCurvLevel);
+                        if (regionDict.readIfPresent("curvatureLevel", curvSpec))
+                        {
+                            if (curvSpec[0] <= 0)
+                            {
+                                FatalIOErrorInFunction(dict)
+                                    << "Illegal curvatureLevel specification for surface "
+                                    << names_[surfI]
+                                    << " : curvatureLevel:" << curvSpec
+                                    << exit(FatalIOError);
+                            }
+                        }
+                        regionCurvLevel[surfI].insert(regionI, curvSpec);
+
                         if (regionDict.found("perpendicularAngle"))
                         {
                             regionAngle[surfI].insert
@@ -505,6 +534,8 @@ Foam::refinementSurfaces::refinementSurfaces
     extendedGapMode_ = volumeType::UNKNOWN;
     selfProximity_.setSize(nRegions);
     selfProximity_ = true;
+    extendedCurvatureLevel_.setSize(nRegions);
+    extendedCurvatureLevel_ = nullCurvLevel;
     perpendicularAngle_.setSize(nRegions);
     perpendicularAngle_ = -GREAT;
     patchInfo_.setSize(nRegions);
@@ -531,6 +562,8 @@ Foam::refinementSurfaces::refinementSurfaces
             extendedGapLevel_[globalRegionI] = globalGapLevel[surfI];
             extendedGapMode_[globalRegionI] = globalGapMode[surfI];
             selfProximity_[globalRegionI] = globalGapSelf[surfI];
+            extendedCurvatureLevel_[globalRegionI] =
+                globalCurvLevel[surfI];
             perpendicularAngle_[globalRegionI] = globalAngle[surfI];
             if (globalPatchInfo.set(surfI))
             {
@@ -560,6 +593,8 @@ Foam::refinementSurfaces::refinementSurfaces
                 regionGapMode[surfI][iter.key()];
             selfProximity_[globalRegionI] =
                 regionGapSelf[surfI][iter.key()];
+            extendedCurvatureLevel_[globalRegionI] =
+                regionCurvLevel[surfI][iter.key()];
         }
         forAllConstIters(regionAngle[surfI], iter)
         {
@@ -713,6 +748,26 @@ Foam::labelList Foam::refinementSurfaces::maxGapLevel() const
 }
 
 
+Foam::labelList Foam::refinementSurfaces::maxCurvatureLevel() const
+{
+    labelList surfaceMax(surfaces_.size(), Zero);
+
+    forAll(surfaces_, surfI)
+    {
+        const wordList& regionNames = allGeometry_[surfaces_[surfI]].regions();
+
+        forAll(regionNames, regionI)
+        {
+            label globalI = globalRegion(surfI, regionI);
+            const FixedList<label, 4>& gapInfo =
+                extendedCurvatureLevel_[globalI];
+            surfaceMax[surfI] = max(surfaceMax[surfI], gapInfo[2]);
+        }
+    }
+    return surfaceMax;
+}
+
+
 // Precalculate the refinement level for every element of the searchable
 // surface.
 void Foam::refinementSurfaces::setMinLevelFields(const shellSurfaces& shells)
@@ -763,13 +818,15 @@ void Foam::refinementSurfaces::setMinLevelFields(const shellSurfaces& shells)
 
             // In case of triangulated surfaces only cache value if triangle
             // centre and vertices are in same shell
-            if (isA<triSurface>(geom))
+            const auto* tsPtr = isA<const triSurface>(geom);
+
+            if (tsPtr)
             {
                 label nUncached = 0;
 
                 // Check if points differing from ctr level
 
-                const triSurface& ts = refCast<const triSurface>(geom);
+                const triSurface& ts = *tsPtr;
                 const pointField& points = ts.points();
 
                 // Determine minimum expected level to avoid having to
@@ -845,6 +902,312 @@ void Foam::refinementSurfaces::setMinLevelFields(const shellSurfaces& shells)
 }
 
 
+// Precalculate the refinement level for every element of the searchable
+// surface.
+void Foam::refinementSurfaces::setCurvatureMinLevelFields
+(
+    const scalar cosAngle,
+    const scalar level0EdgeLength
+)
+{
+    const labelList maxCurvLevel(maxCurvatureLevel());
+
+
+    forAll(surfaces_, surfI)
+    {
+        // Check if there is a specification of the extended curvature
+        if (maxCurvLevel[surfI] <= 0)
+        {
+            continue;
+        } 
+
+        const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
+
+        const auto* tsPtr = isA<const triSurface>(geom);
+
+        // Cache the refinement level (max of surface level and shell level)
+        // on a per-element basis. Only makes sense if there are lots of
+        // elements. Possibly should have 'enough' elements to have fine
+        // enough resolution but for now just make sure we don't catch e.g.
+        // searchableBox (size=6)
+        if (tsPtr)
+        {
+            // Representative local coordinates and bounding sphere
+            pointField ctrs;
+            scalarField radiusSqr;
+            geom.boundingSpheres(ctrs, radiusSqr);
+
+            labelList minLevelField(ctrs.size(), Zero);
+            //labelList surfMin(ctrs.size(), Zero);
+            labelList surfMax(ctrs.size(), Zero);
+            labelList nCurvCells(ctrs.size(), Zero);
+            labelList curvIgnore(ctrs.size(), -1);
+            {
+                // Get the element index in a roundabout way. Problem is e.g.
+                // distributed surface where local indices differ from global
+                // ones (needed for getRegion call)
+                List<pointIndexHit> info;
+                geom.findNearest(ctrs, radiusSqr, info);
+
+                // Get per element the region
+                labelList region;
+                geom.getRegion(info, region);
+
+                // See if a cached level field available (from e.g. shells)
+                labelList cachedField;
+                geom.getField(info, cachedField);
+
+                // From the region get the surface-wise refinement level
+                forAll(minLevelField, i)
+                {
+                    if (info[i].hit()) //Note: should not be necessary
+                    {
+                        const label globali = globalRegion(surfI, region[i]);
+                        curvIgnore[i] = extendedCurvatureLevel_[globali][3];
+                        nCurvCells[i] = extendedCurvatureLevel_[globali][0];
+                        //surfMin[i] = extendedCurvatureLevel_[globali][1];
+                        surfMax[i] = extendedCurvatureLevel_[globali][2];
+
+                        minLevelField[i] = minLevel(surfI, region[i]);
+                        if (cachedField.size() > i)
+                        {
+                            minLevelField[i] =
+                                max(minLevelField[i], cachedField[i]);
+                        }
+                    }
+                }
+            }
+
+            // Calculate per-triangle curvature. This is the max of the
+            // measured point-based curvature + some constraints.
+            scalarField cellCurv(ctrs.size(), Zero);
+            {
+                // Walk surface and detect sharp features. Returns maximum
+                // curvature (per surface point. Note: returns per point, not
+                // per localpoint)
+                const auto& ts = *tsPtr;
+                auto tcurv(triSurfaceTools::curvatures(ts));
+                auto& curv = tcurv.ref();
+
+                // Reset curvature on sharp edges (and neighbours since
+                // curvature uses extended stencil)
+                {
+                    const auto& edgeFaces = ts.edgeFaces();
+                    const auto& edges = ts.edges();
+                    const auto& points = ts.points();
+                    const auto& mp = ts.meshPoints();
+
+                    bitSet isOnSharpEdge(points.size());
+                    forAll(edgeFaces, edgei)
+                    {
+                        const auto& eFaces = edgeFaces[edgei];
+                        const edge meshE(mp, edges[edgei]);
+
+                        if (eFaces.size() == 2)
+                        {
+                            const auto& f0 = ts[eFaces[0]];
+                            const auto& f1 = ts[eFaces[1]];
+
+                            const vector n0 = f0.unitNormal(points);
+
+                            const int dir0 = f0.edgeDirection(meshE);
+                            const int dir1 = f1.edgeDirection(meshE);
+                            vector n1 = f1.unitNormal(points);
+                            if (dir0 == dir1)
+                            {
+                                // Flip since use edge in same direction
+                                // (should not be the case for 'proper'
+                                //  surfaces)
+                                n1 = -n1;
+                            }
+
+                            if ((n0&n1) < cosAngle)
+                            {
+                                isOnSharpEdge.set(meshE[0]);
+                                isOnSharpEdge.set(meshE[1]);
+                            }
+                        }
+                    }
+
+                    // Extend by one layer
+                    {
+                        bitSet oldOnSharpEdge(isOnSharpEdge);
+                        isOnSharpEdge = false;
+                        for (const auto& f : ts)
+                        {
+                            for (const label pointi : f)
+                            {
+                                if (oldOnSharpEdge[pointi])
+                                {
+                                    // Mark all points on triangle
+                                    isOnSharpEdge.set(f);
+                                    break;
+                                }
+                            }
+                        }
+                    }
+
+
+                    // Unmark curvature
+                    autoPtr<OBJstream> str;
+                    //if (debug)
+                    //{
+                    //    str.reset
+                    //    (
+                    //        new OBJstream
+                    //        (
+                    //            "sharpEdgePoints_"
+                    //          + geom.name()
+                    //          + ".obj"
+                    //        )
+                    //    );
+                    //    Info<< "Writing sharp edge points to "
+                    //        << str().name() << endl;
+                    //}
+
+                    for (const label pointi : isOnSharpEdge)
+                    {
+                        // Reset measured point-based curvature
+                        curv[pointi] = 0.0;
+                        if (str)
+                        {
+                            str().write(points[pointi]);
+                        }
+                    }
+                }
+
+                // Reset curvature on -almost- sharp edges.
+                // This resets the point-based curvature if the edge
+                // is considered to be a sharp edge based on its actual
+                // curvature. This is only used if the 'ignore' level is
+                // given.
+                {
+                    // Pass 1: accumulate constraints on the points - get
+                    //         the minimum of curvature constraints on the
+                    //         connected triangles. Looks at user-specified
+                    //         min curvature - does not use actual measured
+                    //         curvature
+                    scalarField pointMinCurv(ts.nPoints(), VGREAT);
+
+                    forAll(ts, i)
+                    {
+                        // Is ignore level given for surface
+                        const label level = curvIgnore[i];
+                        if (level >= 0)
+                        {
+                            // Convert to (inv) size
+                            const scalar length = level0EdgeLength/(2<<level);
+                            const scalar invLength = 1.0/length;
+                            for (const label pointi : ts[i])
+                            {
+                                if
+                                (
+                                    invLength < pointMinCurv[pointi]
+                                 && curv[pointi] > SMALL
+                                )
+                                {
+                                    //Pout<< "** at location:"
+                                    //    << ts.points()[pointi]
+                                    //    << " measured curv:" << curv[pointi]
+                                    //    << " radius:" << 1.0/curv[pointi]
+                                    //    << " ignore level:" << level
+                                    //    << " ignore radius:" << length
+                                    //    << " resetting minCurv to "
+                                    //    << invLength
+                                    //    << endl;
+                                }
+
+                                pointMinCurv[pointi] =
+                                    min(pointMinCurv[pointi], invLength);
+                            }
+                        }
+                    }
+
+                    // Clip curvature (should do nothing for most points unless
+                    // ignore-level is triggered)
+                    forAll(pointMinCurv, pointi)
+                    {
+                        if (pointMinCurv[pointi] < curv[pointi])
+                        {
+                            //Pout<< "** at location:" << ts.points()[pointi]
+                            //    << " measured curv:" << curv[pointi]
+                            //    << " radius:" << 1.0/curv[pointi]
+                            //    << " cellLimit:" << pointMinCurv[pointi]
+                            //    << endl;
+
+                            // Set up to ignore point
+                            //curv[pointi] = pointMinCurv[pointi];
+                            curv[pointi] = 0.0;
+                        }
+                    }
+                }
+
+
+                forAll(ts, i)
+                {
+                    const auto& f = ts[i];
+                    // Take max curvature (= min radius of curvature)
+                    cellCurv[i] = max(curv[f[0]], max(curv[f[1]], curv[f[2]]));
+                }
+            }
+
+
+            //if(debug)
+            //{
+            //    const scalar maxCurv = gMax(cellCurv);
+            //    if (maxCurv > SMALL)
+            //    {
+            //        const scalar r = scalar(1.0)/maxCurv;
+            //
+            //        Pout<< "For geometry " << geom.name()
+            //            << " have curvature max " << maxCurv
+            //            << " which equates to radius:" << r
+            //            << " which equates to refinement level "
+            //            << log2(level0EdgeLength/r)
+            //            << endl;
+            //    }
+            //}
+
+            forAll(cellCurv, i)
+            {
+                if (cellCurv[i] > SMALL && nCurvCells[i] > 0)
+                {
+                    //- ?If locally have a cached field override the
+                    //  surface-based level ignore any curvature?
+                    //if (minLevelField[i] > surfMin[i])
+                    //{
+                    //    // Ignore curvature
+                    //}
+                    //else
+                    //if (surfMin[i] == surfMax[i])
+                    //{
+                    //    // Ignore curvature. Bypass calculation below.
+                    //}
+                    //else
+                    {
+                        // Re-work the curvature into a radius and into a
+                        // number of cells
+                        const scalar r = scalar(1.0)/cellCurv[i];
+                        const scalar level =
+                            log2(nCurvCells[i]*level0EdgeLength/r);
+                        const label l = round(level);
+
+                        if (l > minLevelField[i] && l <= surfMax[i])
+                        {
+                            minLevelField[i] = l;
+                        }
+                    }
+                }
+            }
+
+
+            // Store minLevelField on surface
+            const_cast<searchableSurface&>(geom).setField(minLevelField);
+        }
+    }
+}
+
+
 // Find intersections of edge. Return -1 or first surface with higher minLevel
 // number.
 void Foam::refinementSurfaces::findHigherIntersection
diff --git a/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.H b/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.H
index 16d8582c09..2c4ff38380 100644
--- a/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.H
+++ b/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.H
@@ -110,6 +110,9 @@ class refinementSurfaces
         //  (in gap refinement)
         boolList selfProximity_;
 
+        //- From global region number to curvature specification
+        List<FixedList<label, 4>> extendedCurvatureLevel_;
+
         //- From global region number to perpendicular angle
         scalarField perpendicularAngle_;
 
@@ -263,6 +266,19 @@ public:
                 return selfProximity_;
             }
 
+            //- From global region number to specification of curvature
+            //  refinement: 4 labels specifying
+            //  - minimum wanted number of cells in the curvature radius
+            //  - ?minimum cell level when to start trying to detect gaps
+            //  - maximum cell level to refine to (so do not detect curvature if
+            //    celllevel >= maximum level)
+            //  - minimum radius to ignore (expressed as refinement level).
+            //    This can be used to ignore feature-edges. Set to -1 to ignore.
+            const List<FixedList<label, 4>>& extendedCurvatureLevel() const
+            {
+                return extendedCurvatureLevel_;
+            }
+
             //- From global region number to perpendicular angle
             const scalarField& perpendicularAngle() const
             {
@@ -307,12 +323,23 @@ public:
             //- Per surface the maximum extendedGapLevel over all its regions
             labelList maxGapLevel() const;
 
-            //- Calculate minLevelFields
+            //- Per surface the maximum curvatureLevel over all its regions
+            labelList maxCurvatureLevel() const;
+
+            //- Calculate minLevelFields according to both surface- and
+            //  shell-based levels
             void setMinLevelFields
             (
                 const shellSurfaces& shells
             );
 
+            //- Update minLevelFields according to (triSurface-only) curvature
+            void setCurvatureMinLevelFields
+            (
+                const scalar cosAngle,
+                const scalar level0EdgeLength
+            );
+
             ////- Helper: count number of triangles per region
             //static labelList countRegions(const triSurface&);
 
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C
index 5b336ed18f..7fc5cccc86 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C
@@ -230,10 +230,16 @@ Foam::label Foam::snappyRefineDriver::smallFeatureRefine
     label iter = 0;
 
     // See if any surface has an extendedGapLevel
-    labelList surfaceMaxLevel(meshRefiner_.surfaces().maxGapLevel());
-    labelList shellMaxLevel(meshRefiner_.shells().maxGapLevel());
+    const labelList surfaceMaxLevel(meshRefiner_.surfaces().maxGapLevel());
+    const labelList shellMaxLevel(meshRefiner_.shells().maxGapLevel());
+    const labelList curvMaxLevel(meshRefiner_.surfaces().maxCurvatureLevel());
 
-    if (max(surfaceMaxLevel) == 0 && max(shellMaxLevel) == 0)
+    if
+    (
+        max(surfaceMaxLevel) == 0
+     && max(shellMaxLevel) == 0
+     && max(curvMaxLevel) == 0
+    )
     {
         return iter;
     }
@@ -3382,6 +3388,7 @@ void Foam::snappyRefineDriver::doRefine
     (
         max(meshRefiner_.surfaces().maxGapLevel()) > 0
      || max(meshRefiner_.shells().maxGapLevel()) > 0
+     || max(meshRefiner_.surfaces().maxCurvatureLevel()) > 0
     )
     {
         // In case we use automatic gap level refinement do some pre-refinement
diff --git a/tutorials/mesh/snappyHexMesh/block_with_curvature/Allclean b/tutorials/mesh/snappyHexMesh/block_with_curvature/Allclean
new file mode 100755
index 0000000000..21a3ef143c
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/block_with_curvature/Allclean
@@ -0,0 +1,12 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/CleanFunctions      # Tutorial clean functions
+#------------------------------------------------------------------------------
+
+cleanCase0
+
+# Remove surface and features
+rm -rf constant/triSurface
+rm -rf constant/extendedFeatureEdgeMesh
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/mesh/snappyHexMesh/block_with_curvature/Allrun b/tutorials/mesh/snappyHexMesh/block_with_curvature/Allrun
new file mode 100755
index 0000000000..02e4d82318
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/block_with_curvature/Allrun
@@ -0,0 +1,18 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+mkdir -p constant/triSurface
+
+cp -f \
+    resources/geometry/curvature-box.stl.gz \
+    constant/triSurface
+
+runApplication surfaceFeatureExtract
+
+runApplication blockMesh
+
+runApplication snappyHexMesh -overwrite
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/mesh/snappyHexMesh/block_with_curvature/README.txt b/tutorials/mesh/snappyHexMesh/block_with_curvature/README.txt
new file mode 100644
index 0000000000..6a02065b59
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/block_with_curvature/README.txt
@@ -0,0 +1,6 @@
+testcase to demonstrate surface-curvature-based refinement.
+- starts from a single cell
+- has uniform surface refinement 0
+- but has curvature detection to refine up to 10
+- also uses limitRegions to not refine half the domain (in x direction)
+- note that there is still refinement bleeding due to the 2:1 limitation
diff --git a/tutorials/mesh/snappyHexMesh/block_with_curvature/resources/geometry/curvature-box.stl.gz b/tutorials/mesh/snappyHexMesh/block_with_curvature/resources/geometry/curvature-box.stl.gz
new file mode 100644
index 0000000000..925f4bdeaf
Binary files /dev/null and b/tutorials/mesh/snappyHexMesh/block_with_curvature/resources/geometry/curvature-box.stl.gz differ
diff --git a/tutorials/mesh/snappyHexMesh/block_with_curvature/system/blockMeshDict b/tutorials/mesh/snappyHexMesh/block_with_curvature/system/blockMeshDict
new file mode 100644
index 0000000000..388d4d4238
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/block_with_curvature/system/blockMeshDict
@@ -0,0 +1,58 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2112                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scale   1;
+
+vertices
+(
+    ( 0     -0.2  0)
+    ( 0.5   -0.2  0)
+    ( 0.5    0.3  0)
+    ( 0      0.3  0)
+    ( 0     -0.2  0.5)
+    ( 0.5   -0.2  0.5)
+    ( 0.5    0.3  0.5)
+    ( 0      0.3  0.5)
+);
+
+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) (1 1 1) simpleGrading (1 1 1)
+);
+
+edges
+();
+
+boundary
+(
+    allBoundary
+    {
+        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
+        );
+    }
+);
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/snappyHexMesh/block_with_curvature/system/controlDict b/tutorials/mesh/snappyHexMesh/block_with_curvature/system/controlDict
new file mode 100644
index 0000000000..48f8127e14
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/block_with_curvature/system/controlDict
@@ -0,0 +1,55 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2112                                 |
+|   \\  /    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;
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/snappyHexMesh/block_with_curvature/system/fvSchemes b/tutorials/mesh/snappyHexMesh/block_with_curvature/system/fvSchemes
new file mode 100644
index 0000000000..bd9bab5f29
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/block_with_curvature/system/fvSchemes
@@ -0,0 +1,53 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2112                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+}
+
+gradSchemes
+{
+}
+
+divSchemes
+{
+}
+
+laplacianSchemes
+{
+}
+
+interpolationSchemes
+{
+}
+
+snGradSchemes
+{
+}
+
+wallDist
+{
+}
+
+geometry
+{
+    type            highAspectRatio;
+    minAspect       10;
+    maxAspect       100;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/snappyHexMesh/block_with_curvature/system/fvSolution b/tutorials/mesh/snappyHexMesh/block_with_curvature/system/fvSolution
new file mode 100644
index 0000000000..3316c1cc74
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/block_with_curvature/system/fvSolution
@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2112                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/snappyHexMesh/block_with_curvature/system/meshQualityDict b/tutorials/mesh/snappyHexMesh/block_with_curvature/system/meshQualityDict
new file mode 100644
index 0000000000..e06a307879
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/block_with_curvature/system/meshQualityDict
@@ -0,0 +1,78 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2112                                 |
+|   \\  /    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;
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/snappyHexMesh/block_with_curvature/system/snappyHexMeshDict b/tutorials/mesh/snappyHexMesh/block_with_curvature/system/snappyHexMeshDict
new file mode 100644
index 0000000000..b1f8b166e0
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/block_with_curvature/system/snappyHexMeshDict
@@ -0,0 +1,587 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2112                                 |
+|   \\  /    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            false;
+addLayers       false;
+
+
+// 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
+{
+    box_limit
+    {
+        type box;
+        min  (-1000 -1000 -1000);
+        max  (0.125 1000 1000);
+    }
+    cylinder_all
+    {
+        //- Radius 0.1m
+        type triSurfaceMesh;
+        file curvature-box.stl;
+    }
+};
+
+
+// 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
+    (
+    );
+
+
+    // 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
+    {
+        cylinder_all
+        {
+            // Surface-wise min and max refinement level. max level is only used
+            // for sharp angles (see 'resolveFeatureAngle')
+            level           (0 0);
+
+            // Additional refinement for regions of high curvature. Expressed
+            // (bit similar to gapLevel) as:
+            // - number of cells per radius of curvature. (usually 1 is
+            //   good enough)
+            // - starting cell level? Not used at the moment.
+            // - maximum cell level. This can be smaller or larger than the
+            //   max 'surface' level
+            // - minumum curvature radius to ignore (expressed as a cell level).
+            //   This can be used to avoid detecting small sharp surface
+            //   features. Set to -1 to ignore.
+            //
+            // Sometimes you want more refinement than sharp features since
+            // these can be done with feature edge snapping (so can leave
+            // 'level (0 0)')
+            curvatureLevel (10 0 10 -1);
+
+            // To trigger small feature refinement
+            //gapLevel        (1 0 1);
+        }
+    }
+
+    // Feature angle:
+    // - used if min and max refinement level of a surface differ
+    // - used if feature snapping (see snapControls below) is used
+    resolveFeatureAngle 45;
+
+
+    // 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
+    {
+    }
+
+
+
+    // 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
+    {
+        box_limit
+        {
+            mode        inside;
+
+            // Don't refine at all inside 'box_limit'
+            levels      ((10000 0));
+        }
+    }
+
+
+    // 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.101 0.101 0.101);
+
+        // 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;
+}
+
+
+// Debug flags
+//debugFlags
+//(
+//    mesh            // write intermediate meshes
+//);
+
+
+//// 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;
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/snappyHexMesh/block_with_curvature/system/snappyHexMeshDict.orig b/tutorials/mesh/snappyHexMesh/block_with_curvature/system/snappyHexMeshDict.orig
new file mode 100644
index 0000000000..727fc6b77d
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/block_with_curvature/system/snappyHexMeshDict.orig
@@ -0,0 +1,622 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2112                                 |
+|   \\  /    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            false;
+addLayers       false;
+
+
+// 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
+{
+    all
+    {
+        type box;
+        min  (-1000 -1000 -1000);
+        max  (1000 1000 1000);
+    }
+    cylinder_all
+    {
+        //- Radius 0.1m
+        type triSurfaceMesh;
+        file curvature-box.stl;
+    }
+};
+
+
+// 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
+    (
+    );
+
+
+    // 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
+    {
+        cylinder_all
+        {
+            // Surface-wise min and max refinement level. max level is only used
+            // for sharp angles (see 'resolveFeatureAngle')
+            level           (0 0);
+
+            // Additional refinement for regions of high curvature. Expressed
+            // (bit similar to gapLevel) as:
+            // - number of cells per radius of curvature. (usually 1 is
+            //   good enough)
+            // - starting cell level? Not used at the moment.
+            // - maximum cell level. This can be smaller or larger than the
+            //   max 'surface' level
+            // - minumum curvature radius to ignore (expressed as a cell level).
+            //   This can be used to avoid detecting small sharp surface
+            //   features. Set to -1 to ignore.
+            //
+            // Sometimes you want more refinement than sharp features since
+            // these can be done with feature edge snapping (so can leave
+            // 'level (0 0)')
+            curvatureLevel (10 0 10 -1);
+
+            // To trigger small feature refinement
+            //gapLevel        (1 0 1);
+        }
+    }
+
+//- looking at triangle mesh itself? E.g. cables with 3 triangles do not need
+//refinement.
+//- what about 90 degree angles. Do NOT refine.
+//Refine always. With chamfer: currently not
+//  detected. So curvature refinement is just extension of resolveFeatureAngle
+//  mechanism.
+//- disconnected surfaces. Currently not needed. Can also only be on single patch.
+//- triangle clouds. Currently not needed.
+//- limited to volume as well (like gapLevel). Not currently. Happy with either
+//    - global (like resolveFeatureAngle)
+//    - per patch (like resolveFeatureAngle)
+
+    // Feature angle:
+    // - used if min and max refinement level of a surface differ
+    // - used if feature snapping (see snapControls below) is used
+    resolveFeatureAngle 45;
+
+
+    // 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
+    {
+        all
+        {
+            mode        inside;
+
+            // Dummy base level
+            levels      ((10000 0));
+
+            // If cells
+            // - have level 0..9
+            // - and are in a gap < 3 cell sizes across
+            // - with the gap on the inside ('inside'), outside ('outside')
+            //   or both ('mixed') of the surface
+            // refine them
+            //gapLevel    (4 0 10);
+            //gapMode     outside;
+        }
+    }
+
+
+
+    // 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.101 0.101 0.101);
+
+        // 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;
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/snappyHexMesh/block_with_curvature/system/surfaceFeatureExtractDict b/tutorials/mesh/snappyHexMesh/block_with_curvature/system/surfaceFeatureExtractDict
new file mode 100644
index 0000000000..2504325d35
--- /dev/null
+++ b/tutorials/mesh/snappyHexMesh/block_with_curvature/system/surfaceFeatureExtractDict
@@ -0,0 +1,60 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2112                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      surfaceFeatureExtractDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//- temp dictionary with common settings
+_surfaceExtract
+{
+    // 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;
+
+    // Output surface curvature
+    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;
+
+    // Write closeness/curvature/proximity fields as VTK for postprocessing
+    writeVTK            yes;
+}
+
+
+curvature-box.stl
+{
+    $_surfaceExtract;
+}
+
+
+//- Remove temp from dictionary so it doesn't think it is a surface
+#remove _surfaceExtract
+
+// ************************************************************************* //