From 009203188fabc28f3eda77dd1cb31ad589f95cfb Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Thu, 13 Oct 2016 15:05:24 +0100
Subject: [PATCH] blockMesh: New experimental support for projecting block face
 point to geometric surfaces

For example, to mesh a sphere with a single block the geometry is defined in the
blockMeshDict as a searchableSurface:

    geometry
    {
        sphere
        {
            type searchableSphere;
            centre (0 0 0);
            radius 1;
        }
    }

The vertices, block topology and curved edges are defined in the usual
way, for example

    v 0.5773502;
    mv -0.5773502;

    a 0.7071067;
    ma -0.7071067;

    vertices
    (
        ($mv $mv $mv)
        ( $v $mv $mv)
        ( $v  $v $mv)
        ($mv  $v $mv)
        ($mv $mv  $v)
        ( $v $mv  $v)
        ( $v  $v  $v)
        ($mv  $v  $v)
    );

    blocks
    (
        hex (0 1 2 3 4 5 6 7) (10 10 10) simpleGrading (1 1 1)
    );

    edges
    (
        arc 0 1 (0 $ma $ma)
        arc 2 3 (0 $a $ma)
        arc 6 7 (0 $a $a)
        arc 4 5 (0 $ma $a)

        arc 0 3 ($ma 0 $ma)
        arc 1 2 ($a 0 $ma)
        arc 5 6 ($a 0 $a)
        arc 4 7 ($ma 0 $a)

        arc 0 4 ($ma $ma 0)
        arc 1 5 ($a $ma 0)
        arc 2 6 ($a $a 0)
        arc 3 7 ($ma $a 0)
    );

which will produce a mesh in which the block edges conform to the sphere
but the faces of the block lie somewhere between the original cube and
the spherical surface which is a consequence of the edge-based
transfinite interpolation.

Now the projection of the block faces to the geometry specified above
can also be specified:

    faces
    (
        project (0 4 7 3) sphere
        project (2 6 5 1) sphere
        project (1 5 4 0) sphere
        project (3 7 6 2) sphere
        project (0 3 2 1) sphere
        project (4 5 6 7) sphere
    );

which produces a mesh that actually conforms to the sphere.

See OpenFOAM-dev/tutorials/mesh/blockMesh/sphere

This functionality is experimental and will undergo further development
and generalization in the future to support more complex surfaces,
feature edge specification and extraction etc.  Please get involved if
you would like to see blockMesh become a more flexible block-structured
mesher.

Henry G. Weller, CFD Direct.
---
 .../mesh/generation/blockMesh/Make/options    |   2 +
 .../mesh/generation/blockMesh/blockMesh.C     |   6 +-
 .../vtkPV3blockMesh/vtkPV3blockMesh.C         |   2 +-
 .../vtkPV3blockMesh/vtkPV3blockMeshConvert.C  |   4 +-
 .../vtkPVblockMesh/Make/options               |   2 +
 .../vtkPVblockMesh/vtkPVblockMesh.C           |   2 +-
 .../vtkPVblockMesh/vtkPVblockMeshConvert.C    |   6 +-
 src/mesh/blockMesh/Make/files                 |  22 +-
 src/mesh/blockMesh/Make/options               |   4 +-
 src/mesh/blockMesh/block/block.C              |  58 +-
 src/mesh/blockMesh/block/block.H              |  52 +-
 src/mesh/blockMesh/block/blockCreate.C        | 620 ++++++++++--------
 src/mesh/blockMesh/block/blockI.H             |  18 +-
 .../blockDescriptor/blockDescriptor.C         | 122 +++-
 .../blockDescriptor/blockDescriptor.H         |  88 ++-
 .../blockDescriptor/blockDescriptorEdges.C    |  48 +-
 .../blockDescriptor/blockDescriptorI.H        | 115 ++++
 .../blockDescriptor/blockDescriptorTest.C     |  24 +
 .../blockEdges/{ => BSplineEdge}/BSpline.C    |   1 -
 .../blockEdges/{ => BSplineEdge}/BSpline.H    |   0
 .../{ => BSplineEdge}/BSplineEdge.C           |   7 +-
 .../{ => BSplineEdge}/BSplineEdge.H           |   7 +-
 .../blockEdges/{ => arcEdge}/arcEdge.C        |   7 +-
 .../blockEdges/{ => arcEdge}/arcEdge.H        |  11 +-
 .../blockEdges/{ => blockEdge}/blockEdge.C    |  24 +-
 .../blockEdges/{ => blockEdge}/blockEdge.H    |  41 +-
 .../blockEdges/{ => blockEdge}/blockEdgeI.H   |   0
 .../{ => blockEdge}/blockEdgeList.H           |   0
 .../blockEdges/{ => lineDivide}/lineDivide.C  |   0
 .../blockEdges/{ => lineDivide}/lineDivide.H  |   0
 .../blockEdges/{ => lineEdge}/lineEdge.C      |   8 +-
 .../blockEdges/{ => lineEdge}/lineEdge.H      |  12 +-
 .../blockEdges/{ => polyLineEdge}/polyLine.C  |   0
 .../blockEdges/{ => polyLineEdge}/polyLine.H  |   0
 .../{ => polyLineEdge}/polyLineEdge.C         |   8 +-
 .../{ => polyLineEdge}/polyLineEdge.H         |   7 +-
 .../{ => splineEdge}/CatmullRomSpline.C       |   1 -
 .../{ => splineEdge}/CatmullRomSpline.H       |   0
 .../blockEdges/{ => splineEdge}/splineEdge.C  |   7 +-
 .../blockEdges/{ => splineEdge}/splineEdge.H  |   7 +-
 .../blockFaces/blockFace/blockFace.C          |  98 +++
 .../blockFaces/blockFace/blockFace.H          | 162 +++++
 .../blockFaces/blockFace/blockFaceI.H         |  46 ++
 .../blockFaces/blockFace/blockFaceList.H      |  55 ++
 .../blockFaces/projectFace/projectFace.C      | 100 +++
 .../blockFaces/projectFace/projectFace.H      | 111 ++++
 src/mesh/blockMesh/blockMesh/blockMesh.C      |  24 +-
 src/mesh/blockMesh/blockMesh/blockMesh.H      |  27 +-
 .../blockMesh/blockMesh/blockMeshCreate.C     |  92 +--
 src/mesh/blockMesh/blockMesh/blockMeshMerge.C |  46 +-
 .../blockMesh/blockMesh/blockMeshMergeFast.C  |  22 +-
 .../blockMesh/blockMesh/blockMeshTopology.C   |  32 +-
 .../searchableSurfaceCollection.H             |   8 +-
 .../searchableSurface/searchableSurfaces.H    |   3 +
 tutorials/mesh/blockMesh/sphere/Allrun        |   9 +
 .../blockMesh/sphere/system/blockMeshDict     |  97 +++
 .../mesh/blockMesh/sphere/system/controlDict  |  49 ++
 .../mesh/blockMesh/sphere/system/fvSchemes    |  37 ++
 .../mesh/blockMesh/sphere/system/fvSolution   |  18 +
 59 files changed, 1774 insertions(+), 605 deletions(-)
 create mode 100644 src/mesh/blockMesh/blockDescriptor/blockDescriptorTest.C
 rename src/mesh/blockMesh/blockEdges/{ => BSplineEdge}/BSpline.C (99%)
 rename src/mesh/blockMesh/blockEdges/{ => BSplineEdge}/BSpline.H (100%)
 rename src/mesh/blockMesh/blockEdges/{ => BSplineEdge}/BSplineEdge.C (95%)
 rename src/mesh/blockMesh/blockEdges/{ => BSplineEdge}/BSplineEdge.H (95%)
 rename src/mesh/blockMesh/blockEdges/{ => arcEdge}/arcEdge.C (97%)
 rename src/mesh/blockMesh/blockEdges/{ => arcEdge}/arcEdge.H (94%)
 rename src/mesh/blockMesh/blockEdges/{ => blockEdge}/blockEdge.C (90%)
 rename src/mesh/blockMesh/blockEdges/{ => blockEdge}/blockEdge.H (88%)
 rename src/mesh/blockMesh/blockEdges/{ => blockEdge}/blockEdgeI.H (100%)
 rename src/mesh/blockMesh/blockEdges/{ => blockEdge}/blockEdgeList.H (100%)
 rename src/mesh/blockMesh/blockEdges/{ => lineDivide}/lineDivide.C (100%)
 rename src/mesh/blockMesh/blockEdges/{ => lineDivide}/lineDivide.H (100%)
 rename src/mesh/blockMesh/blockEdges/{ => lineEdge}/lineEdge.C (95%)
 rename src/mesh/blockMesh/blockEdges/{ => lineEdge}/lineEdge.H (93%)
 rename src/mesh/blockMesh/blockEdges/{ => polyLineEdge}/polyLine.C (100%)
 rename src/mesh/blockMesh/blockEdges/{ => polyLineEdge}/polyLine.H (100%)
 rename src/mesh/blockMesh/blockEdges/{ => polyLineEdge}/polyLineEdge.C (95%)
 rename src/mesh/blockMesh/blockEdges/{ => polyLineEdge}/polyLineEdge.H (95%)
 rename src/mesh/blockMesh/blockEdges/{ => splineEdge}/CatmullRomSpline.C (99%)
 rename src/mesh/blockMesh/blockEdges/{ => splineEdge}/CatmullRomSpline.H (100%)
 rename src/mesh/blockMesh/blockEdges/{ => splineEdge}/splineEdge.C (95%)
 rename src/mesh/blockMesh/blockEdges/{ => splineEdge}/splineEdge.H (95%)
 create mode 100644 src/mesh/blockMesh/blockFaces/blockFace/blockFace.C
 create mode 100644 src/mesh/blockMesh/blockFaces/blockFace/blockFace.H
 create mode 100644 src/mesh/blockMesh/blockFaces/blockFace/blockFaceI.H
 create mode 100644 src/mesh/blockMesh/blockFaces/blockFace/blockFaceList.H
 create mode 100644 src/mesh/blockMesh/blockFaces/projectFace/projectFace.C
 create mode 100644 src/mesh/blockMesh/blockFaces/projectFace/projectFace.H
 create mode 100755 tutorials/mesh/blockMesh/sphere/Allrun
 create mode 100644 tutorials/mesh/blockMesh/sphere/system/blockMeshDict
 create mode 100644 tutorials/mesh/blockMesh/sphere/system/controlDict
 create mode 100644 tutorials/mesh/blockMesh/sphere/system/fvSchemes
 create mode 100644 tutorials/mesh/blockMesh/sphere/system/fvSolution

diff --git a/applications/utilities/mesh/generation/blockMesh/Make/options b/applications/utilities/mesh/generation/blockMesh/Make/options
index 2e07c477da..9d13a25de6 100644
--- a/applications/utilities/mesh/generation/blockMesh/Make/options
+++ b/applications/utilities/mesh/generation/blockMesh/Make/options
@@ -1,9 +1,11 @@
 EXE_INC = \
     -I$(LIB_SRC)/mesh/blockMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fileFormats/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude
 
 EXE_LIBS = \
     -lblockMesh \
     -lmeshTools \
+    -lfileFormats \
     -ldynamicMesh
diff --git a/applications/utilities/mesh/generation/blockMesh/blockMesh.C b/applications/utilities/mesh/generation/blockMesh/blockMesh.C
index 2be39caae9..a7860c8f53 100644
--- a/applications/utilities/mesh/generation/blockMesh/blockMesh.C
+++ b/applications/utilities/mesh/generation/blockMesh/blockMesh.C
@@ -176,8 +176,6 @@ int main(int argc, char *argv[])
     Info<< "Creating block mesh from\n    "
         << meshDictIO.objectPath() << endl;
 
-    blockMesh::verbose(true);
-
     IOdictionary meshDict(meshDictIO);
     blockMesh blocks(meshDict, regionName);
 
@@ -286,8 +284,8 @@ int main(int argc, char *argv[])
         forAll(blocks, blockI)
         {
             const block& b = blocks[blockI];
-            const labelListList& blockCells = b.cells();
-            const word& zoneName = b.blockDef().zoneName();
+            const List<FixedList<label, 8>> blockCells = b.cells();
+            const word& zoneName = b.zoneName();
 
             if (zoneName.size())
             {
diff --git a/applications/utilities/postProcessing/graphics/PV3Readers/PV3blockMeshReader/vtkPV3blockMesh/vtkPV3blockMesh.C b/applications/utilities/postProcessing/graphics/PV3Readers/PV3blockMeshReader/vtkPV3blockMesh/vtkPV3blockMesh.C
index a3ccb75891..3aea81a887 100644
--- a/applications/utilities/postProcessing/graphics/PV3Readers/PV3blockMeshReader/vtkPV3blockMesh/vtkPV3blockMesh.C
+++ b/applications/utilities/postProcessing/graphics/PV3Readers/PV3blockMeshReader/vtkPV3blockMesh/vtkPV3blockMesh.C
@@ -75,7 +75,7 @@ void Foam::vtkPV3blockMesh::updateInfoBlocks
     const int nBlocks = blkMesh.size();
     for (int blockI = 0; blockI < nBlocks; ++blockI)
     {
-        const blockDescriptor& blockDef = blkMesh[blockI].blockDef();
+        const blockDescriptor& blockDef = blkMesh[blockI];
 
         word partName = Foam::name(blockI);
 
diff --git a/applications/utilities/postProcessing/graphics/PV3Readers/PV3blockMeshReader/vtkPV3blockMesh/vtkPV3blockMeshConvert.C b/applications/utilities/postProcessing/graphics/PV3Readers/PV3blockMeshReader/vtkPV3blockMesh/vtkPV3blockMeshConvert.C
index 876fc47a3e..8dce33a982 100644
--- a/applications/utilities/postProcessing/graphics/PV3Readers/PV3blockMeshReader/vtkPV3blockMesh/vtkPV3blockMeshConvert.C
+++ b/applications/utilities/postProcessing/graphics/PV3Readers/PV3blockMeshReader/vtkPV3blockMesh/vtkPV3blockMeshConvert.C
@@ -77,7 +77,7 @@ void Foam::vtkPV3blockMesh::convertMeshBlocks
             continue;
         }
 
-        const blockDescriptor& blockDef = blkMesh[blockI].blockDef();
+        const blockDescriptor& blockDef = blkMesh[blockI];
 
         vtkUnstructuredGrid* vtkmesh = vtkUnstructuredGrid::New();
 
@@ -168,7 +168,7 @@ void Foam::vtkPV3blockMesh::convertMeshEdges
         // search each block
         forAll(blkMesh, blockI)
         {
-            const blockDescriptor& blockDef = blkMesh[blockI].blockDef();
+            const blockDescriptor& blockDef = blkMesh[blockI];
 
             edgeList blkEdges = blockDef.blockShape().edges();
 
diff --git a/applications/utilities/postProcessing/graphics/PVReaders/PVblockMeshReader/vtkPVblockMesh/Make/options b/applications/utilities/postProcessing/graphics/PVReaders/PVblockMeshReader/vtkPVblockMesh/Make/options
index ccbea978f7..9dcba4b39b 100644
--- a/applications/utilities/postProcessing/graphics/PVReaders/PVblockMeshReader/vtkPVblockMesh/Make/options
+++ b/applications/utilities/postProcessing/graphics/PVReaders/PVblockMeshReader/vtkPVblockMesh/Make/options
@@ -1,5 +1,6 @@
 EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fileFormats/lnInclude \
     -I$(LIB_SRC)/mesh/blockMesh/lnInclude \
     -I$(ParaView_INCLUDE_DIR) \
     -I$(ParaView_INCLUDE_DIR)/vtkkwiml \
@@ -8,6 +9,7 @@ EXE_INC = \
 
 LIB_LIBS = \
     -lmeshTools \
+    -lfileFormats \
     -lblockMesh \
     -L$(FOAM_LIBBIN) -lvtkPVReaders \
     $(GLIBS)
diff --git a/applications/utilities/postProcessing/graphics/PVReaders/PVblockMeshReader/vtkPVblockMesh/vtkPVblockMesh.C b/applications/utilities/postProcessing/graphics/PVReaders/PVblockMeshReader/vtkPVblockMesh/vtkPVblockMesh.C
index 22cdf2db50..a9932ffcc9 100644
--- a/applications/utilities/postProcessing/graphics/PVReaders/PVblockMeshReader/vtkPVblockMesh/vtkPVblockMesh.C
+++ b/applications/utilities/postProcessing/graphics/PVReaders/PVblockMeshReader/vtkPVblockMesh/vtkPVblockMesh.C
@@ -75,7 +75,7 @@ void Foam::vtkPVblockMesh::updateInfoBlocks
     const int nBlocks = blkMesh.size();
     for (int blockI = 0; blockI < nBlocks; ++blockI)
     {
-        const blockDescriptor& blockDef = blkMesh[blockI].blockDef();
+        const blockDescriptor& blockDef = blkMesh[blockI];
 
         word partName = Foam::name(blockI);
 
diff --git a/applications/utilities/postProcessing/graphics/PVReaders/PVblockMeshReader/vtkPVblockMesh/vtkPVblockMeshConvert.C b/applications/utilities/postProcessing/graphics/PVReaders/PVblockMeshReader/vtkPVblockMesh/vtkPVblockMeshConvert.C
index a287be2dd9..9c7f8a06ba 100644
--- a/applications/utilities/postProcessing/graphics/PVReaders/PVblockMeshReader/vtkPVblockMesh/vtkPVblockMeshConvert.C
+++ b/applications/utilities/postProcessing/graphics/PVReaders/PVblockMeshReader/vtkPVblockMesh/vtkPVblockMeshConvert.C
@@ -77,7 +77,7 @@ void Foam::vtkPVblockMesh::convertMeshBlocks
             continue;
         }
 
-        const blockDescriptor& blockDef = blkMesh[blockI].blockDef();
+        const blockDescriptor& blockDef = blkMesh[blockI];
 
         vtkUnstructuredGrid* vtkmesh = vtkUnstructuredGrid::New();
 
@@ -168,12 +168,12 @@ void Foam::vtkPVblockMesh::convertMeshEdges
         // search each block
         forAll(blkMesh, blockI)
         {
-            const blockDescriptor& blockDef = blkMesh[blockI].blockDef();
+            const blockDescriptor& blockDef = blkMesh[blockI];
 
             edgeList blkEdges = blockDef.blockShape().edges();
 
             // List of edge point and weighting factors
-            List<point> edgesPoints[12];
+            pointField edgesPoints[12];
             scalarList edgesWeights[12];
             blockDef.edgesPointsWeights(edgesPoints, edgesWeights);
 
diff --git a/src/mesh/blockMesh/Make/files b/src/mesh/blockMesh/Make/files
index 27bed65597..98e1c9cf6e 100644
--- a/src/mesh/blockMesh/Make/files
+++ b/src/mesh/blockMesh/Make/files
@@ -1,14 +1,16 @@
-blockEdges/BSpline.C
-blockEdges/CatmullRomSpline.C
-blockEdges/polyLine.C
+blockEdges/blockEdge/blockEdge.C
+blockEdges/lineDivide/lineDivide.C
+blockEdges/lineEdge/lineEdge.C
+blockEdges/polyLineEdge/polyLine.C
+blockEdges/polyLineEdge/polyLineEdge.C
+blockEdges/arcEdge/arcEdge.C
+blockEdges/BSplineEdge/BSpline.C
+blockEdges/BSplineEdge/BSplineEdge.C
+blockEdges/splineEdge/CatmullRomSpline.C
+blockEdges/splineEdge/splineEdge.C
 
-blockEdges/arcEdge.C
-blockEdges/blockEdge.C
-blockEdges/lineEdge.C
-blockEdges/polyLineEdge.C
-blockEdges/lineDivide.C
-blockEdges/BSplineEdge.C
-blockEdges/splineEdge.C
+blockFaces/blockFace/blockFace.C
+blockFaces/projectFace/projectFace.C
 
 gradingDescriptor/gradingDescriptor.C
 gradingDescriptor/gradingDescriptors.C
diff --git a/src/mesh/blockMesh/Make/options b/src/mesh/blockMesh/Make/options
index 1618ab57ec..3be93bd21f 100644
--- a/src/mesh/blockMesh/Make/options
+++ b/src/mesh/blockMesh/Make/options
@@ -1,7 +1,7 @@
 EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/dynamicMesh/lnInclude
+    -I$(LIB_SRC)/fileFormats/lnInclude
 
 LIB_LIBS = \
     -lmeshTools \
-    -ldynamicMesh
+    -lfileFormats
diff --git a/src/mesh/blockMesh/block/block.C b/src/mesh/blockMesh/block/block.C
index 0b345511f6..1722b0095e 100644
--- a/src/mesh/blockMesh/block/block.C
+++ b/src/mesh/blockMesh/block/block.C
@@ -31,63 +31,23 @@ Foam::block::block
 (
     const pointField& vertices,
     const blockEdgeList& edges,
+    const blockFaceList& faces,
     Istream& is
 )
 :
-    blockDescriptor(vertices, edges, is),
-    points_(0),
-    cells_(0),
-    boundaryPatches_(0)
-{}
+    blockDescriptor(vertices, edges, faces, is)
+{
+    createPoints();
+    createBoundary();
+}
 
 
 Foam::block::block(const blockDescriptor& blockDesc)
 :
-    blockDescriptor(blockDesc),
-    points_(0),
-    cells_(0),
-    boundaryPatches_(0)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::block::~block()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-const Foam::pointField& Foam::block::points() const
+    blockDescriptor(blockDesc)
 {
-    if (points_.empty())
-    {
-        createPoints();
-    }
-
-    return points_;
-}
-
-
-const Foam::labelListList& Foam::block::cells() const
-{
-    if (cells_.empty())
-    {
-        createCells();
-    }
-
-    return cells_;
-}
-
-
-const Foam::labelListListList& Foam::block::boundaryPatches() const
-{
-    if (boundaryPatches_.empty())
-    {
-        createBoundary();
-    }
-
-    return boundaryPatches_;
+    createPoints();
+    createBoundary();
 }
 
 
diff --git a/src/mesh/blockMesh/block/block.H b/src/mesh/blockMesh/block/block.H
index 0f48b975e0..7b55c5f59f 100644
--- a/src/mesh/blockMesh/block/block.H
+++ b/src/mesh/blockMesh/block/block.H
@@ -68,25 +68,19 @@ class block
     // Private data
 
         //- List of points
-        mutable pointField points_;
-
-        //- List of cells
-        mutable labelListList cells_;
+        pointField points_;
 
         //- Boundary patches
-        mutable labelListListList boundaryPatches_;
+        FixedList<List<FixedList<label, 4>>, 6> boundaryPatches_;
 
 
     // Private Member Functions
 
         //- Creates vertices for cells filling the block
-        void createPoints() const;
-
-        //- Creates cells for filling the block
-        void createCells() const;
+        void createPoints();
 
         //- Creates boundary patch faces for the block
-        void createBoundary() const;
+        void createBoundary();
 
         //- Disallow default bitwise copy construct
         block(const block&);
@@ -94,6 +88,7 @@ class block
         //- Disallow default bitwise assignment
         void operator=(const block&);
 
+
 public:
 
     // Constructors
@@ -103,6 +98,7 @@ public:
         (
             const pointField& vertices,
             const blockEdgeList& edges,
+            const blockFaceList& faces,
             Istream& is
         );
 
@@ -122,50 +118,42 @@ public:
         {
             const pointField& points_;
             const blockEdgeList& edges_;
+            const blockFaceList& faces_;
 
         public:
 
-            iNew(const pointField& points, const blockEdgeList& edges)
+            iNew
+            (
+                const pointField& points,
+                const blockEdgeList& edges,
+                const blockFaceList& faces
+            )
             :
                 points_(points),
-                edges_(edges)
+                edges_(edges),
+                faces_(faces)
             {}
 
             autoPtr<block> operator()(Istream& is) const
             {
-                return autoPtr<block>(new block(points_, edges_, is));
+                return autoPtr<block>(new block(points_, edges_, faces_, is));
             }
         };
 
 
-    //- Destructor
-    ~block();
-
-
     // Member Functions
 
         // Access
 
-            //- Return the block definition
-            inline const blockDescriptor& blockDef() const;
-
-            //- Vertex label offset for a particular i,j,k position
-            inline label vtxLabel(label i, label j, label k) const;
-
             //- Return the points for filling the block
-            const pointField& points() const;
+            inline const pointField& points() const;
 
             //- Return the cells for filling the block
-            const labelListList& cells() const;
+            List<FixedList<label, 8>> cells() const;
 
             //- Return the boundary patch faces for the block
-            const labelListListList& boundaryPatches() const;
-
-
-        // Edit
-
-            //- Clear geometry (internal points, cells, boundaryPatches)
-            void clearGeom();
+            inline const FixedList<List<FixedList<label, 4>>, 6>&
+                boundaryPatches() const;
 
 
     // Ostream Operator
diff --git a/src/mesh/blockMesh/block/blockCreate.C b/src/mesh/blockMesh/block/blockCreate.C
index 0aac98f197..8f63d629fd 100644
--- a/src/mesh/blockMesh/block/blockCreate.C
+++ b/src/mesh/blockMesh/block/blockCreate.C
@@ -27,7 +27,22 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::block::createPoints() const
+#define w0 w[0][i]
+#define w1 w[1][i]
+#define w2 w[2][i]
+#define w3 w[3][i]
+
+#define w4 w[4][j]
+#define w5 w[5][j]
+#define w6 w[6][j]
+#define w7 w[7][j]
+
+#define w8 w[8][k]
+#define w9 w[9][k]
+#define w10 w[10][k]
+#define w11 w[11][k]
+
+void Foam::block::createPoints()
 {
     // Set local variables for mesh specification
     const label ni = density().x();
@@ -45,168 +60,285 @@ void Foam::block::createPoints() const
     const point& p011 = blockPoint(7);
 
     // List of edge point and weighting factors
-    List<point> p[12];
+    pointField p[12];
     scalarList w[12];
-    edgesPointsWeights(p, w);
+    label nCurvedEdges = edgesPointsWeights(p, w);
 
-    //
-    // Generate vertices
-    //
-    points_.clear();
     points_.setSize(nPoints());
 
-    for (label k = 0; k <= nk; k++)
-    {
-        for (label j = 0; j <= nj; j++)
-        {
-            for (label i = 0; i <= ni; i++)
-            {
-                const label vertexNo = vtxLabel(i, j, k);
+    points_[pointLabel(0,  0,  0)] = p000;
+    points_[pointLabel(ni, 0,  0)] = p100;
+    points_[pointLabel(ni, nj, 0)] = p110;
+    points_[pointLabel(0,  nj, 0)] = p010;
+    points_[pointLabel(0,  0,  nk)] = p001;
+    points_[pointLabel(ni, 0,  nk)] = p101;
+    points_[pointLabel(ni, nj, nk)] = p111;
+    points_[pointLabel(0,  nj, nk)] = p011;
 
-                // Calculate the importance factors for all edges
+    for (label k=0; k<=nk; k++)
+    {
+        for (label j=0; j<=nj; j++)
+        {
+            for (label i=0; i<=ni; i++)
+            {
+                // Skip block vertices
+                if (vertex(i, j, k)) continue;
+
+                const label vijk = pointLabel(i, j, k);
+
+                // Calculate the weighting factors for all edges
 
                 // x-direction
-                scalar impx1 =
-                (
-                    (1 - w[0][i])*(1 - w[4][j])*(1 - w[8][k])
-                  + w[0][i]*(1 - w[5][j])*(1 - w[9][k])
-                );
+                scalar wx1 = (1 - w0)*(1 - w4)*(1 - w8) + w0*(1 - w5)*(1 - w9);
+                scalar wx2 = (1 - w1)*w4*(1 - w11)      + w1*w5*(1 - w10);
+                scalar wx3 = (1 - w2)*w7*w11            + w2*w6*w10;
+                scalar wx4 = (1 - w3)*(1 - w7)*w8       + w3*(1 - w6)*w9;
 
-                scalar impx2 =
-                (
-                    (1 - w[1][i])*w[4][j]*(1 - w[11][k])
-                  + w[1][i]*w[5][j]*(1 - w[10][k])
-                );
-
-                scalar impx3 =
-                (
-                     (1 - w[2][i])*w[7][j]*w[11][k]
-                   + w[2][i]*w[6][j]*w[10][k]
-                );
-
-                scalar impx4 =
-                (
-                    (1 - w[3][i])*(1 - w[7][j])*w[8][k]
-                  + w[3][i]*(1 - w[6][j])*w[9][k]
-                );
-
-                const scalar magImpx = impx1 + impx2 + impx3 + impx4;
-                impx1 /= magImpx;
-                impx2 /= magImpx;
-                impx3 /= magImpx;
-                impx4 /= magImpx;
+                const scalar sumWx = wx1 + wx2 + wx3 + wx4;
+                wx1 /= sumWx;
+                wx2 /= sumWx;
+                wx3 /= sumWx;
+                wx4 /= sumWx;
 
 
                 // y-direction
-                scalar impy1 =
-                (
-                    (1 - w[4][j])*(1 - w[0][i])*(1 - w[8][k])
-                  + w[4][j]*(1 - w[1][i])*(1 - w[11][k])
-                );
+                scalar wy1 = (1 - w4)*(1 - w0)*(1 - w8) + w4*(1 - w1)*(1 - w11);
+                scalar wy2 = (1 - w5)*w0*(1 - w9)       + w5*w1*(1 - w10);
+                scalar wy3 = (1 - w6)*w3*w9             + w6*w2*w10;
+                scalar wy4 = (1 - w7)*(1 - w3)*w8       + w7*(1 - w2)*w11;
 
-                scalar impy2 =
-                (
-                    (1 - w[5][j])*w[0][i]*(1 - w[9][k])
-                  + w[5][j]*w[1][i]*(1 - w[10][k])
-                );
-
-                scalar impy3 =
-                (
-                    (1 - w[6][j])*w[3][i]*w[9][k]
-                  + w[6][j]*w[2][i]*w[10][k]
-                );
-
-                scalar impy4 =
-                (
-                    (1 - w[7][j])*(1 - w[3][i])*w[8][k]
-                  + w[7][j]*(1 - w[2][i])*w[11][k]
-                );
-
-                const scalar magImpy = impy1 + impy2 + impy3 + impy4;
-                impy1 /= magImpy;
-                impy2 /= magImpy;
-                impy3 /= magImpy;
-                impy4 /= magImpy;
+                const scalar sumWy = wy1 + wy2 + wy3 + wy4;
+                wy1 /= sumWy;
+                wy2 /= sumWy;
+                wy3 /= sumWy;
+                wy4 /= sumWy;
 
 
                 // z-direction
-                scalar impz1 =
-                (
-                    (1 - w[8][k])*(1 - w[0][i])*(1 - w[4][j])
-                  + w[8][k]*(1 - w[3][i])*(1 - w[7][j])
-                );
+                scalar wz1 = (1 - w8)*(1 - w0)*(1 - w4) + w8*(1 - w3)*(1 - w7);
+                scalar wz2 = (1 - w9)*w0*(1 - w5)       + w9*w3*(1 - w6);
+                scalar wz3 = (1 - w10)*w1*w5            + w10*w2*w6;
+                scalar wz4 = (1 - w11)*(1 - w1)*w4      + w11*(1 - w2)*w7;
 
-                scalar impz2 =
-                (
-                    (1 - w[9][k])*w[0][i]*(1 - w[5][j])
-                  + w[9][k]*w[3][i]*(1 - w[6][j])
-                );
-
-                scalar impz3 =
-                (
-                    (1 - w[10][k])*w[1][i]*w[5][j]
-                  + w[10][k]*w[2][i]*w[6][j]
-                );
-
-                scalar impz4 =
-                (
-                    (1 - w[11][k])*(1 - w[1][i])*w[4][j]
-                  + w[11][k]*(1 - w[2][i])*w[7][j]
-                );
-
-                const scalar magImpz = impz1 + impz2 + impz3 + impz4;
-                impz1 /= magImpz;
-                impz2 /= magImpz;
-                impz3 /= magImpz;
-                impz4 /= magImpz;
+                const scalar sumWz = wz1 + wz2 + wz3 + wz4;
+                wz1 /= sumWz;
+                wz2 /= sumWz;
+                wz3 /= sumWz;
+                wz4 /= sumWz;
 
 
                 // Points on straight edges
-                const vector edgex1 = p000 + (p100 - p000)*w[0][i];
-                const vector edgex2 = p010 + (p110 - p010)*w[1][i];
-                const vector edgex3 = p011 + (p111 - p011)*w[2][i];
-                const vector edgex4 = p001 + (p101 - p001)*w[3][i];
+                const vector edgex1 = p000 + (p100 - p000)*w0;
+                const vector edgex2 = p010 + (p110 - p010)*w1;
+                const vector edgex3 = p011 + (p111 - p011)*w2;
+                const vector edgex4 = p001 + (p101 - p001)*w3;
 
-                const vector edgey1 = p000 + (p010 - p000)*w[4][j];
-                const vector edgey2 = p100 + (p110 - p100)*w[5][j];
-                const vector edgey3 = p101 + (p111 - p101)*w[6][j];
-                const vector edgey4 = p001 + (p011 - p001)*w[7][j];
+                const vector edgey1 = p000 + (p010 - p000)*w4;
+                const vector edgey2 = p100 + (p110 - p100)*w5;
+                const vector edgey3 = p101 + (p111 - p101)*w6;
+                const vector edgey4 = p001 + (p011 - p001)*w7;
 
-                const vector edgez1 = p000 + (p001 - p000)*w[8][k];
-                const vector edgez2 = p100 + (p101 - p100)*w[9][k];
-                const vector edgez3 = p110 + (p111 - p110)*w[10][k];
-                const vector edgez4 = p010 + (p011 - p010)*w[11][k];
+                const vector edgez1 = p000 + (p001 - p000)*w8;
+                const vector edgez2 = p100 + (p101 - p100)*w9;
+                const vector edgez3 = p110 + (p111 - p110)*w10;
+                const vector edgez4 = p010 + (p011 - p010)*w11;
 
                 // Add the contributions
-                points_[vertexNo] =
+                points_[vijk] =
                 (
-                    impx1*edgex1 + impx2*edgex2 + impx3*edgex3 + impx4*edgex4
-                  + impy1*edgey1 + impy2*edgey2 + impy3*edgey3 + impy4*edgey4
-                  + impz1*edgez1 + impz2*edgez2 + impz3*edgez3 + impz4*edgez4
-                )/3.0;
+                    wx1*edgex1 + wx2*edgex2 + wx3*edgex3 + wx4*edgex4
+                  + wy1*edgey1 + wy2*edgey2 + wy3*edgey3 + wy4*edgey4
+                  + wz1*edgez1 + wz2*edgez2 + wz3*edgez3 + wz4*edgez4
+                )/3;
 
 
-                // Calculate the correction vectors
-                const vector corx1 = impx1*(p[0][i] - edgex1);
-                const vector corx2 = impx2*(p[1][i] - edgex2);
-                const vector corx3 = impx3*(p[2][i] - edgex3);
-                const vector corx4 = impx4*(p[3][i] - edgex4);
+                // Apply curved-edge correction if block has curved edges
+                if (nCurvedEdges)
+                {
+                    // Calculate the correction vectors
+                    const vector corx1 = wx1*(p[0][i] - edgex1);
+                    const vector corx2 = wx2*(p[1][i] - edgex2);
+                    const vector corx3 = wx3*(p[2][i] - edgex3);
+                    const vector corx4 = wx4*(p[3][i] - edgex4);
 
-                const vector cory1 = impy1*(p[4][j] - edgey1);
-                const vector cory2 = impy2*(p[5][j] - edgey2);
-                const vector cory3 = impy3*(p[6][j] - edgey3);
-                const vector cory4 = impy4*(p[7][j] - edgey4);
+                    const vector cory1 = wy1*(p[4][j] - edgey1);
+                    const vector cory2 = wy2*(p[5][j] - edgey2);
+                    const vector cory3 = wy3*(p[6][j] - edgey3);
+                    const vector cory4 = wy4*(p[7][j] - edgey4);
 
-                const vector corz1 = impz1*(p[8][k] - edgez1);
-                const vector corz2 = impz2*(p[9][k] - edgez2);
-                const vector corz3 = impz3*(p[10][k] - edgez3);
-                const vector corz4 = impz4*(p[11][k] - edgez4);
+                    const vector corz1 = wz1*(p[8][k] - edgez1);
+                    const vector corz2 = wz2*(p[9][k] - edgez2);
+                    const vector corz3 = wz3*(p[10][k] - edgez3);
+                    const vector corz4 = wz4*(p[11][k] - edgez4);
 
-                points_[vertexNo] +=
+                    points_[vijk] +=
+                    (
+                        corx1 + corx2 + corx3 + corx4
+                      + cory1 + cory2 + cory3 + cory4
+                      + corz1 + corz2 + corz3 + corz4
+                    );
+                }
+            }
+        }
+    }
+
+    if (!nCurvedFaces()) return;
+
+    // Apply curvature correction to face points
+    FixedList<pointField, 6> facePoints(this->facePoints(points_));
+    correctFacePoints(facePoints);
+
+    // Loop over points and apply the correction from the from the i-faces
+    for (label ii=0; ii<=ni; ii++)
+    {
+        // Process the points on the curved faces last
+        label i = (ii + 1)%(ni + 1);
+
+        for (label j=0; j<=nj; j++)
+        {
+            for (label k=0; k<=nk; k++)
+            {
+                // Skip non-curved faces and edges
+                if (flatFaceOrEdge(i, j, k)) continue;
+
+                const label vijk = pointLabel(i, j, k);
+
+                scalar wf0 =
                 (
-                    corx1 + corx2 + corx3 + corx4
-                  + cory1 + cory2 + cory3 + cory4
-                  + corz1 + corz2 + corz3 + corz4
+                    (1 - w0)*(1 - w4)*(1 - w8)
+                  + (1 - w1)*w4*(1 - w11)
+                  + (1 - w2)*w7*w11
+                  + (1 - w3)*(1 - w7)*w8
+                );
+
+                scalar wf1 =
+                (
+                    w0*(1 - w5)*(1 - w9)
+                  + w1*w5*(1 - w10)
+                  + w2*w5*w10
+                  + w3*(1 - w6)*w9
+                );
+
+                const scalar sumWf = wf0 + wf1;
+                wf0 /= sumWf;
+                wf1 /= sumWf;
+
+                points_[vijk] +=
+                (
+                    wf0
+                   *(
+                       facePoints[0][facePointLabel(0, j, k)]
+                     - points_[pointLabel(0, j, k)]
+                    )
+                  + wf1
+                   *(
+                       facePoints[1][facePointLabel(1, j, k)]
+                     - points_[pointLabel(ni, j, k)]
+                    )
+                );
+            }
+        }
+    }
+
+    // Loop over points and apply the correction from the from the j-faces
+    for (label jj=0; jj<=nj; jj++)
+    {
+        // Process the points on the curved faces last
+        label j = (jj + 1)%(nj + 1);
+
+        for (label i=0; i<=ni; i++)
+        {
+            for (label k=0; k<=nk; k++)
+            {
+                // Skip flat faces and edges
+                if (flatFaceOrEdge(i, j, k)) continue;
+
+                const label vijk = pointLabel(i, j, k);
+
+                scalar wf2 =
+                (
+                    (1 - w4)*(1 - w1)*(1 - w8)
+                  + (1 - w5)*w0*(1 - w9)
+                  + (1 - w6)*w3*w9
+                  + (1 - w7)*(1 - w3)*w8
+                );
+
+                scalar wf3 =
+                (
+                    w4*(1 - w1)*(1 - w11)
+                  + w5*w1*(1 - w10)
+                  + w6*w2*w10
+                  + w7*(1 - w2)*w11
+                );
+
+                const scalar sumWf = wf2 + wf3;
+                wf2 /= sumWf;
+                wf3 /= sumWf;
+
+                points_[vijk] +=
+                (
+                    wf2
+                   *(
+                       facePoints[2][facePointLabel(2, i, k)]
+                     - points_[pointLabel(i, 0, k)]
+                    )
+                  + wf3
+                   *(
+                       facePoints[3][facePointLabel(3, i, k)]
+                     - points_[pointLabel(i, nj, k)]
+                    )
+                );
+            }
+        }
+    }
+
+    // Loop over points and apply the correction from the from the k-faces
+    for (label kk=0; kk<=nk; kk++)
+    {
+        // Process the points on the curved faces last
+        label k = (kk + 1)%(nk + 1);
+
+        for (label i=0; i<=ni; i++)
+        {
+            for (label j=0; j<=nj; j++)
+            {
+                // Skip flat faces and edges
+                if (flatFaceOrEdge(i, j, k)) continue;
+
+                const label vijk = pointLabel(i, j, k);
+
+                scalar wf4 =
+                (
+                    (1 - w8)*(1 - w0)*(1 - w4)
+                  + (1 - w9)*w0*(1 - w5)
+                  + (1 - w10)*w1*w5
+                  + (1 - w11)*(1 - w1)*w4
+                );
+
+                scalar wf5 =
+                (
+                    w8*(1 - w3)*(1 - w7)
+                  + w9*w3*(1 - w6)
+                  + w10*w2*w6
+                  + w11*(1 - w2)*w7
+                );
+
+                const scalar sumWf = wf4 + wf5;
+                wf4 /= sumWf;
+                wf5 /= sumWf;
+
+                points_[vijk] +=
+                (
+                    wf4
+                   *(
+                       facePoints[4][facePointLabel(4, i, j)]
+                     - points_[pointLabel(i, j, 0)]
+                    )
+                  + wf5
+                   *(
+                       facePoints[5][facePointLabel(5, i, j)]
+                     - points_[pointLabel(i, j, nk)]
+                    )
                 );
             }
         }
@@ -214,230 +346,164 @@ void Foam::block::createPoints() const
 }
 
 
-void Foam::block::createCells() const
+Foam::List<Foam::FixedList<Foam::label, 8>> Foam::block::cells() const
 {
     const label ni = density().x();
     const label nj = density().y();
     const label nk = density().z();
 
-    //
-    // Generate cells
-    //
-    cells_.clear();
-    cells_.setSize(nCells());
+    List<FixedList<label, 8>> cells(nCells());
 
-    label cellNo = 0;
+    label celli = 0;
 
-    for (label k = 0; k < nk; k++)
+    for (label k=0; k<nk; k++)
     {
-        for (label j = 0; j < nj; j++)
+        for (label j=0; j<nj; j++)
         {
-            for (label i = 0; i < ni; i++)
+            for (label i=0; i<ni; i++)
             {
-                cells_[cellNo].setSize(8);
+                cells[celli][0] = pointLabel(i,   j,   k);
+                cells[celli][1] = pointLabel(i+1, j,   k);
+                cells[celli][2] = pointLabel(i+1, j+1, k);
+                cells[celli][3] = pointLabel(i,   j+1, k);
+                cells[celli][4] = pointLabel(i,   j,   k+1);
+                cells[celli][5] = pointLabel(i+1, j,   k+1);
+                cells[celli][6] = pointLabel(i+1, j+1, k+1);
+                cells[celli][7] = pointLabel(i,   j+1, k+1);
 
-                cells_[cellNo][0] =  vtxLabel(i, j, k);
-                cells_[cellNo][1] =  vtxLabel(i+1, j, k);
-                cells_[cellNo][2] =  vtxLabel(i+1, j+1, k);
-                cells_[cellNo][3] =  vtxLabel(i, j+1, k);
-                cells_[cellNo][4] =  vtxLabel(i, j, k+1);
-                cells_[cellNo][5] =  vtxLabel(i+1, j, k+1);
-                cells_[cellNo][6] =  vtxLabel(i+1, j+1, k+1);
-                cells_[cellNo][7] =  vtxLabel(i, j+1, k+1);
-                cellNo++;
+                celli++;
             }
         }
     }
+
+    return cells;
 }
 
 
-void Foam::block::createBoundary() const
+void Foam::block::createBoundary()
 {
     const label ni = density().x();
     const label nj = density().y();
     const label nk = density().z();
 
-    //
-    // Generate boundaries on each side of the hex
-    //
-    boundaryPatches_.clear();
-    boundaryPatches_.setSize(6);
-
+    label patchi = 0;
+    label facei = 0;
 
     // x-direction
 
-    label wallLabel = 0;
-    label wallCellLabel = 0;
-
     // x-min
-    boundaryPatches_[wallLabel].setSize(nj*nk);
-    for (label k = 0; k < nk; k++)
+    boundaryPatches_[patchi].setSize(nj*nk);
+    for (label k=0; k<nk; k++)
     {
-        for (label j = 0; j < nj; j++)
+        for (label j=0; j<nj; j++)
         {
-            boundaryPatches_[wallLabel][wallCellLabel].setSize(4);
+            boundaryPatches_[patchi][facei][0] = pointLabel(0, j,   k);
+            boundaryPatches_[patchi][facei][1] = pointLabel(0, j,   k+1);
+            boundaryPatches_[patchi][facei][2] = pointLabel(0, j+1, k+1);
+            boundaryPatches_[patchi][facei][3] = pointLabel(0, j+1, k);
 
-            // Set the points
-            boundaryPatches_[wallLabel][wallCellLabel][0] =
-                vtxLabel(0, j, k);
-            boundaryPatches_[wallLabel][wallCellLabel][1] =
-                vtxLabel(0, j, k + 1);
-            boundaryPatches_[wallLabel][wallCellLabel][2] =
-                vtxLabel(0, j + 1, k + 1);
-            boundaryPatches_[wallLabel][wallCellLabel][3] =
-                vtxLabel(0, j + 1, k);
-
-            // Update the counter
-            wallCellLabel++;
+            facei++;
         }
     }
 
     // x-max
-    wallLabel++;
-    wallCellLabel = 0;
+    patchi++;
+    facei = 0;
 
-    boundaryPatches_[wallLabel].setSize(nj*nk);
+    boundaryPatches_[patchi].setSize(nj*nk);
 
-    for (label k = 0; k < nk; k++)
+    for (label k=0; k<nk; k++)
     {
-        for (label j = 0; j < nj; j++)
+        for (label j=0; j<nj; j++)
         {
-            boundaryPatches_[wallLabel][wallCellLabel].setSize(4);
+            boundaryPatches_[patchi][facei][0] = pointLabel(ni, j,   k);
+            boundaryPatches_[patchi][facei][1] = pointLabel(ni, j+1, k);
+            boundaryPatches_[patchi][facei][2] = pointLabel(ni, j+1, k+1);
+            boundaryPatches_[patchi][facei][3] = pointLabel(ni, j,   k+1);
 
-            // Set the points
-            boundaryPatches_[wallLabel][wallCellLabel][0] =
-                vtxLabel(ni, j, k);
-            boundaryPatches_[wallLabel][wallCellLabel][1] =
-                vtxLabel(ni, j+1, k);
-            boundaryPatches_[wallLabel][wallCellLabel][2] =
-                vtxLabel(ni, j+1, k+1);
-            boundaryPatches_[wallLabel][wallCellLabel][3] =
-                vtxLabel(ni, j, k+1);
-
-            // Update the counter
-            wallCellLabel++;
+            facei++;
         }
     }
 
     // y-direction
 
     // y-min
-    wallLabel++;
-    wallCellLabel = 0;
+    patchi++;
+    facei = 0;
 
-    boundaryPatches_[wallLabel].setSize(ni*nk);
-    for (label i = 0; i < ni; i++)
+    boundaryPatches_[patchi].setSize(ni*nk);
+    for (label i=0; i<ni; i++)
     {
-        for (label k = 0; k < nk; k++)
+        for (label k=0; k<nk; k++)
         {
-            boundaryPatches_[wallLabel][wallCellLabel].setSize(4);
+            boundaryPatches_[patchi][facei][0] = pointLabel(i,   0, k);
+            boundaryPatches_[patchi][facei][1] = pointLabel(i+1, 0, k);
+            boundaryPatches_[patchi][facei][2] = pointLabel(i+1, 0, k+1);
+            boundaryPatches_[patchi][facei][3] = pointLabel(i,   0, k+1);
 
-            // Set the points
-            boundaryPatches_[wallLabel][wallCellLabel][0] =
-                vtxLabel(i, 0, k);
-            boundaryPatches_[wallLabel][wallCellLabel][1] =
-                vtxLabel(i + 1, 0, k);
-            boundaryPatches_[wallLabel][wallCellLabel][2] =
-                vtxLabel(i + 1, 0, k + 1);
-            boundaryPatches_[wallLabel][wallCellLabel][3] =
-                vtxLabel(i, 0, k + 1);
-
-            // Update the counter
-            wallCellLabel++;
+            facei++;
         }
     }
 
     // y-max
-    wallLabel++;
-    wallCellLabel = 0;
+    patchi++;
+    facei = 0;
 
-    boundaryPatches_[wallLabel].setSize(ni*nk);
+    boundaryPatches_[patchi].setSize(ni*nk);
 
-    for (label i = 0; i < ni; i++)
+    for (label i=0; i<ni; i++)
     {
-        for (label k = 0; k < nk; k++)
+        for (label k=0; k<nk; k++)
         {
-            boundaryPatches_[wallLabel][wallCellLabel].setSize(4);
+            boundaryPatches_[patchi][facei][0] = pointLabel(i,   nj, k);
+            boundaryPatches_[patchi][facei][1] = pointLabel(i,   nj, k+1);
+            boundaryPatches_[patchi][facei][2] = pointLabel(i+1, nj, k+1);
+            boundaryPatches_[patchi][facei][3] = pointLabel(i+1, nj, k);
 
-            // Set the points
-            boundaryPatches_[wallLabel][wallCellLabel][0] =
-                vtxLabel(i, nj, k);
-            boundaryPatches_[wallLabel][wallCellLabel][1] =
-                vtxLabel(i, nj, k + 1);
-            boundaryPatches_[wallLabel][wallCellLabel][2] =
-                vtxLabel(i + 1, nj, k + 1);
-            boundaryPatches_[wallLabel][wallCellLabel][3] =
-                vtxLabel(i + 1, nj, k);
-
-            // Update the counter
-            wallCellLabel++;
+            facei++;
         }
     }
 
     // z-direction
 
     // z-min
-    wallLabel++;
-    wallCellLabel = 0;
+    patchi++;
+    facei = 0;
 
-    boundaryPatches_[wallLabel].setSize(ni*nj);
+    boundaryPatches_[patchi].setSize(ni*nj);
 
-    for (label i = 0; i < ni; i++)
+    for (label i=0; i<ni; i++)
     {
-        for (label j = 0; j < nj; j++)
+        for (label j=0; j<nj; j++)
         {
-            boundaryPatches_[wallLabel][wallCellLabel].setSize(4);
+            boundaryPatches_[patchi][facei][0] = pointLabel(i,   j,   0);
+            boundaryPatches_[patchi][facei][1] = pointLabel(i,   j+1, 0);
+            boundaryPatches_[patchi][facei][2] = pointLabel(i+1, j+1, 0);
+            boundaryPatches_[patchi][facei][3] = pointLabel(i+1, j,   0);
 
-            // Set the points
-            boundaryPatches_[wallLabel][wallCellLabel][0] =
-                vtxLabel(i, j, 0);
-            boundaryPatches_[wallLabel][wallCellLabel][1] =
-                vtxLabel(i, j + 1, 0);
-            boundaryPatches_[wallLabel][wallCellLabel][2] =
-                vtxLabel(i + 1, j + 1, 0);
-            boundaryPatches_[wallLabel][wallCellLabel][3] =
-                vtxLabel(i + 1, j, 0);
-
-            // Update the counter
-            wallCellLabel++;
+            facei++;
         }
     }
 
     // z-max
-    wallLabel++;
-    wallCellLabel = 0;
+    patchi++;
+    facei = 0;
 
-    boundaryPatches_[wallLabel].setSize(ni*nj);
+    boundaryPatches_[patchi].setSize(ni*nj);
 
-    for (label i = 0; i < ni; i++)
+    for (label i=0; i<ni; i++)
     {
-        for (label j = 0; j < nj; j++)
+        for (label j=0; j<nj; j++)
         {
-            boundaryPatches_[wallLabel][wallCellLabel].setSize(4);
+            boundaryPatches_[patchi][facei][0] = pointLabel(i,   j,   nk);
+            boundaryPatches_[patchi][facei][1] = pointLabel(i+1, j,   nk);
+            boundaryPatches_[patchi][facei][2] = pointLabel(i+1, j+1, nk);
+            boundaryPatches_[patchi][facei][3] = pointLabel(i,   j+1, nk);
 
-            // Set the points
-            boundaryPatches_[wallLabel][wallCellLabel][0] =
-                vtxLabel(i, j, nk);
-            boundaryPatches_[wallLabel][wallCellLabel][1] =
-                vtxLabel(i + 1, j, nk);
-            boundaryPatches_[wallLabel][wallCellLabel][2] =
-                vtxLabel(i + 1, j + 1, nk);
-            boundaryPatches_[wallLabel][wallCellLabel][3] =
-                vtxLabel(i, j + 1, nk);
-
-            // Update the counter
-            wallCellLabel++;
+            facei++;
         }
     }
 }
 
 
-void Foam::block::clearGeom()
-{
-    points_.clear();
-    cells_.clear();
-    boundaryPatches_.clear();
-}
-
-
 // ************************************************************************* //
diff --git a/src/mesh/blockMesh/block/blockI.H b/src/mesh/blockMesh/block/blockI.H
index 10ac991eb8..27711dc587 100644
--- a/src/mesh/blockMesh/block/blockI.H
+++ b/src/mesh/blockMesh/block/blockI.H
@@ -23,24 +23,18 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline Foam::label Foam::block::vtxLabel(label i, label j, label k) const
+inline const Foam::pointField& Foam::block::points() const
 {
-    return
-    (
-        i
-      + j * (density().x() + 1)
-      + k * (density().x() + 1) * (density().y() + 1)
-    );
+    return points_;
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-inline const Foam::blockDescriptor& Foam::block::blockDef() const
+inline const Foam::FixedList<Foam::List<Foam::FixedList<Foam::label, 4>>, 6>&
+Foam::block::boundaryPatches() const
 {
-    return *this;
+    return boundaryPatches_;
 }
 
 
diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C
index 2d6b308ecd..7ccf24d422 100644
--- a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C
+++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C
@@ -102,6 +102,32 @@ void Foam::blockDescriptor::check(const Istream& is)
 }
 
 
+void Foam::blockDescriptor::findCurvedFaces()
+{
+    const faceList blockFaces(blockShape().faces());
+
+    forAll(blockFaces, blockFacei)
+    {
+        forAll(faces_, facei)
+        {
+            if
+            (
+                face::sameVertices
+                (
+                    faces_[facei].vertices(),
+                    blockFaces[blockFacei]
+                )
+            )
+            {
+                curvedFaces_[blockFacei] = facei;
+                nCurvedFaces_++;
+                break;
+            }
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::blockDescriptor::blockDescriptor
@@ -109,6 +135,7 @@ Foam::blockDescriptor::blockDescriptor
     const cellShape& bshape,
     const pointField& vertices,
     const blockEdgeList& edges,
+    const blockFaceList& faces,
     const Vector<label>& density,
     const UList<gradingDescriptors>& expand,
     const word& zoneName
@@ -116,10 +143,13 @@ Foam::blockDescriptor::blockDescriptor
 :
     vertices_(vertices),
     edges_(edges),
+    faces_(faces),
     blockShape_(bshape),
     density_(density),
     expand_(expand),
-    zoneName_(zoneName)
+    zoneName_(zoneName),
+    curvedFaces_(-1),
+    nCurvedFaces_(0)
 {
     if (expand_.size() != 12)
     {
@@ -127,6 +157,8 @@ Foam::blockDescriptor::blockDescriptor
             << "Unknown definition of expansion ratios"
             << exit(FatalError);
     }
+
+    findCurvedFaces();
 }
 
 
@@ -134,19 +166,19 @@ Foam::blockDescriptor::blockDescriptor
 (
     const pointField& vertices,
     const blockEdgeList& edges,
+    const blockFaceList& faces,
     Istream& is
 )
 :
     vertices_(vertices),
     edges_(edges),
+    faces_(faces),
     blockShape_(is),
     density_(),
-    expand_
-    (
-        12,
-        gradingDescriptors()
-    ),
-    zoneName_()
+    expand_(12, gradingDescriptors()),
+    zoneName_(),
+    curvedFaces_(-1),
+    nCurvedFaces_(0)
 {
     // Examine next token
     token t(is);
@@ -231,13 +263,83 @@ Foam::blockDescriptor::blockDescriptor
     }
 
     check(is);
+
+    findCurvedFaces();
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::blockDescriptor::~blockDescriptor()
-{}
+Foam::FixedList<Foam::pointField, 6>
+Foam::blockDescriptor::facePoints(const pointField& points) const
+{
+    // Caches points for curvature correction
+    FixedList<pointField, 6> facePoints;
+
+    // Set local variables for mesh specification
+    const label ni = density_.x();
+    const label nj = density_.y();
+    const label nk = density_.z();
+
+    facePoints[0].setSize((nj + 1)*(nk + 1));
+    facePoints[1].setSize((nj + 1)*(nk + 1));
+
+    for (label j=0; j<=nj; j++)
+    {
+        for (label k=0; k<=nk; k++)
+        {
+            facePoints[0][facePointLabel(0, j, k)] =
+                points[pointLabel(0, j, k)];
+            facePoints[1][facePointLabel(1, j, k)] =
+                points[pointLabel(ni, j, k)];
+        }
+    }
+
+    facePoints[2].setSize((ni + 1)*(nk + 1));
+    facePoints[3].setSize((ni + 1)*(nk + 1));
+
+    for (label i=0; i<=ni; i++)
+    {
+        for (label k=0; k<=nk; k++)
+        {
+            facePoints[2][facePointLabel(2, i, k)] =
+                points[pointLabel(i, 0, k)];
+            facePoints[3][facePointLabel(3, i, k)] =
+                points[pointLabel(i, nj, k)];
+        }
+    }
+
+    facePoints[4].setSize((ni + 1)*(nj + 1));
+    facePoints[5].setSize((ni + 1)*(nj + 1));
+
+    for (label i=0; i<=ni; i++)
+    {
+        for (label j=0; j<=nj; j++)
+        {
+            facePoints[4][facePointLabel(4, i, j)] =
+                points[pointLabel(i, j, 0)];
+            facePoints[5][facePointLabel(5, i, j)] =
+                points[pointLabel(i, j, nk)];
+        }
+    }
+
+    return facePoints;
+}
+
+
+void Foam::blockDescriptor::correctFacePoints
+(
+    FixedList<pointField, 6>& facePoints
+) const
+{
+    forAll(curvedFaces_, blockFacei)
+    {
+        if (curvedFaces_[blockFacei] != -1)
+        {
+            faces_[curvedFaces_[blockFacei]].project(facePoints[blockFacei]);
+        }
+    }
+}
 
 
 // * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H
index 0a95cacd07..fd7ac5d1eb 100644
--- a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H
+++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H
@@ -1,4 +1,4 @@
-/*---------------------------------------------------------------------------*\
+/*---------------------------------------------------------------------------* \
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
@@ -31,14 +31,14 @@ Description
     For a given block, the correspondence between the ordering of vertex labels
     and face labels is shown below.  For vertex numbering in the sequence 0 to 7
     (block, centre): faces 0 (f0) and 1 are left and right, respectively; faces
-    2 and 3 are bottom and top; and faces 4 and 5 are front the back:
+    2 and 3 are front the back; and faces 4 and 5 are bottom and top:
     \verbatim
              4 ---- 5
-        f3   |\     |\   f5
+        f5   |\     |\   f3
          |   | 7 ---- 6    \
          |   0 |--- 1 |     \
-         |    \|     \|      f4
-        f2     3 ---- 2
+         |    \|     \|      f2
+        f4     3 ---- 2
 
              f0 ----- f1
      \endverbatim
@@ -56,6 +56,7 @@ SourceFiles
 #include "pointField.H"
 #include "scalarList.H"
 #include "blockEdgeList.H"
+#include "blockFaceList.H"
 #include "gradingDescriptors.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -84,6 +85,9 @@ class blockDescriptor
         //- Reference to a list of block edges
         const blockEdgeList& edges_;
 
+        //- Reference to the list of curved faces
+        const blockFaceList& faces_;
+
         //- Block shape
         cellShape blockShape_;
 
@@ -96,16 +100,23 @@ class blockDescriptor
         //- Name of the zone (empty word if none)
         word zoneName_;
 
+        //- Curved-face labels for each block-face (-1 for flat faces)
+        FixedList<label, 6> curvedFaces_;
+
+        //- Number of curved faces in this block
+        label nCurvedFaces_;
+
 
     // Private Member Functions
 
         //- Check block has outward-pointing faces
         void check(const Istream& is);
 
-        //- Calculate the points and weights for the specified edge
-        void edgePointsWeights
+        //- Calculate the points and weights for the specified edge.
+        //  Return the number of curved edges
+        label edgePointsWeights
         (
-            List<point> (&edgePoints)[12],
+            pointField (&edgePoints)[12],
             scalarList (&edgeWeights)[12],
             const label edgei,
             const label start,
@@ -113,6 +124,8 @@ class blockDescriptor
             const label dim
         ) const;
 
+        void findCurvedFaces();
+
 
     // Private Member Functions
 
@@ -130,6 +143,7 @@ public:
             const cellShape&,
             const pointField& vertices,
             const blockEdgeList&,
+            const blockFaceList&,
             const Vector<label>& density,
             const UList<gradingDescriptors>& expand,
             const word& zoneName = ""
@@ -140,17 +154,15 @@ public:
         (
             const pointField& vertices,
             const blockEdgeList&,
+            const blockFaceList&,
             Istream&
         );
 
 
-    //- Destructor
-    ~blockDescriptor();
-
-
     // Member Functions
 
-        // Access
+        //- Return reference to the list of curved faces
+        inline const blockFaceList& faces() const;
 
         //- Return the block shape
         inline const cellShape& blockShape() const;
@@ -167,16 +179,60 @@ public:
         //- Return the number of cells
         inline label nCells() const;
 
+        //- Curved-face labels for each block-face (-1 for flat faces)
+        inline const FixedList<label, 6>& curvedFaces() const;
+
+        //- Number of curved faces in this block
+        inline label nCurvedFaces() const;
+
         //- Return block point for local label i
         inline const point& blockPoint(const label i) const;
 
-        //- Calculate the points and weights for all edges
-        void edgesPointsWeights
+        //- Vertex label offset for a particular i,j,k position
+        inline label pointLabel
         (
-            List<point> (&edgePoints)[12],
+            const label i,
+            const label j,
+            const label k
+        ) const;
+
+        //- Face vertex label offset for a particular i,j,k position
+        inline label facePointLabel
+        (
+            const label facei,
+            const label i,
+            const label j
+        ) const;
+
+        //- Return true if point i,j,k addresses a block vertex
+        inline bool vertex(const label i, const label j, const label k) const;
+
+        //- Return true if point i,j,k addresses a block edge
+        inline bool edge(const label i, const label j, const label k) const;
+
+        //- Calculate the points and weights for all edges.
+        //  Return the number of curved edges
+        label edgesPointsWeights
+        (
+            pointField (&edgePoints)[12],
             scalarList (&edgeWeights)[12]
         ) const;
 
+        //- Return true if point i,j,k addresses a block flat face or edge
+        inline bool flatFaceOrEdge
+        (
+            const label i,
+            const label j,
+            const label k
+        ) const;
+
+        //- Return the list of face-points for all of the faces of the block
+        FixedList<pointField, 6> facePoints(const pointField& points) const;
+
+        //- Correct the location of the given face-points
+        //  to lie on the faces of the block
+        void correctFacePoints(FixedList<pointField, 6>&) const;
+
 
     // IOstream Operators
 
diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptorEdges.C b/src/mesh/blockMesh/blockDescriptor/blockDescriptorEdges.C
index 55895c34a8..77ca51abdd 100644
--- a/src/mesh/blockMesh/blockDescriptor/blockDescriptorEdges.C
+++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptorEdges.C
@@ -29,9 +29,9 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::blockDescriptor::edgePointsWeights
+Foam::label Foam::blockDescriptor::edgePointsWeights
 (
-    List<point> (&edgePoints)[12],
+    pointField (&edgePoints)[12],
     scalarList (&edgeWeights)[12],
     const label edgei,
     const label start,
@@ -79,16 +79,16 @@ void Foam::blockDescriptor::edgePointsWeights
                 edgePoints[edgei].setSize(p.size());
                 edgeWeights[edgei].setSize(d.size());
 
-                label pMax = p.size() - 1;
-                forAll(p, pI)
+                label pn = p.size() - 1;
+                forAll(p, pi)
                 {
-                    edgePoints[edgei][pI]  = p[pMax - pI];
-                    edgeWeights[edgei][pI] = 1.0 - d[pMax - pI];
+                    edgePoints[edgei][pi]  = p[pn - pi];
+                    edgeWeights[edgei][pi] = 1 - d[pn - pi];
                 }
             }
 
             // Found curved-edge: done
-            return;
+            return 1;
         }
     }
 
@@ -103,39 +103,43 @@ void Foam::blockDescriptor::edgePointsWeights
 
     edgePoints[edgei]  = divEdge.points();
     edgeWeights[edgei] = divEdge.lambdaDivisions();
+
+    return 0;
 }
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::blockDescriptor::edgesPointsWeights
+Foam::label Foam::blockDescriptor::edgesPointsWeights
 (
-    List<point> (&edgePoints)[12],
+    pointField (&edgePoints)[12],
     scalarList (&edgeWeights)[12]
 ) const
 {
-    // These edges correspond to the "hex" cellModel
+    label nCurvedEdges = 0;
 
     // X-direction
     const label ni = density_.x();
-    edgePointsWeights(edgePoints, edgeWeights, 0,  0, 1, ni);
-    edgePointsWeights(edgePoints, edgeWeights, 1,  3, 2, ni);
-    edgePointsWeights(edgePoints, edgeWeights, 2,  7, 6, ni);
-    edgePointsWeights(edgePoints, edgeWeights, 3,  4, 5, ni);
+    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 0,  0, 1, ni);
+    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 1,  3, 2, ni);
+    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 2,  7, 6, ni);
+    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 3,  4, 5, ni);
 
     // Y-direction
     const label nj = density_.y();
-    edgePointsWeights(edgePoints, edgeWeights, 4,  0, 3, nj);
-    edgePointsWeights(edgePoints, edgeWeights, 5,  1, 2, nj);
-    edgePointsWeights(edgePoints, edgeWeights, 6,  5, 6, nj);
-    edgePointsWeights(edgePoints, edgeWeights, 7,  4, 7, nj);
+    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 4,  0, 3, nj);
+    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 5,  1, 2, nj);
+    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 6,  5, 6, nj);
+    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 7,  4, 7, nj);
 
     // Z-direction
     const label nk = density_.z();
-    edgePointsWeights(edgePoints, edgeWeights, 8,  0, 4, nk);
-    edgePointsWeights(edgePoints, edgeWeights, 9,  1, 5, nk);
-    edgePointsWeights(edgePoints, edgeWeights, 10, 2, 6, nk);
-    edgePointsWeights(edgePoints, edgeWeights, 11, 3, 7, nk);
+    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 8,  0, 4, nk);
+    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 9,  1, 5, nk);
+    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 10, 2, 6, nk);
+    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 11, 3, 7, nk);
+
+    return nCurvedEdges;
 }
 
 
diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptorI.H b/src/mesh/blockMesh/blockDescriptor/blockDescriptorI.H
index 37439b1a0b..823860a820 100644
--- a/src/mesh/blockMesh/blockDescriptor/blockDescriptorI.H
+++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptorI.H
@@ -25,6 +25,12 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+inline const Foam::blockFaceList& Foam::blockDescriptor::faces() const
+{
+    return faces_;
+}
+
+
 inline const Foam::cellShape& Foam::blockDescriptor::blockShape() const
 {
     return blockShape_;
@@ -65,10 +71,119 @@ inline Foam::label Foam::blockDescriptor::nCells() const
 }
 
 
+inline const Foam::FixedList<Foam::label, 6>&
+Foam::blockDescriptor::curvedFaces() const
+{
+    return curvedFaces_;
+}
+
+
+inline Foam::label Foam::blockDescriptor::nCurvedFaces() const
+{
+    return nCurvedFaces_;
+}
+
+
 inline const Foam::point& Foam::blockDescriptor::blockPoint(const label i) const
 {
     return vertices_[blockShape_[i]];
 }
 
 
+inline Foam::label Foam::blockDescriptor::pointLabel
+(
+    const label i,
+    const label j,
+    const label k
+) const
+{
+    return
+    (
+        i
+      + j*(density_.x() + 1)
+      + k*(density_.x() + 1)*(density_.y() + 1)
+    );
+}
+
+
+inline Foam::label Foam::blockDescriptor::facePointLabel
+(
+    const label facei,
+    const label i,
+    const label j
+) const
+{
+    if (facei == 0 || facei == 1)
+    {
+        return
+        (
+            i
+          + j*(density_.y() + 1)
+        );
+    }
+    else if (facei == 2 || facei == 3)
+    {
+        return
+        (
+            i
+          + j*(density_.x() + 1)
+        );
+    }
+    else
+    {
+        return
+        (
+            i
+          + j*(density_.x() + 1)
+        );
+    }
+}
+
+
+inline bool Foam::blockDescriptor::vertex
+(
+    const label i, const label j, const label k
+) const
+{
+    bool iEnd = (i == 0 || i == density_.x());
+    bool jEnd = (j == 0 || j == density_.y());
+    bool kEnd = (k == 0 || k == density_.z());
+
+    return (iEnd && jEnd && kEnd);
+}
+
+
+inline bool Foam::blockDescriptor::edge
+(
+    const label i, const label j, const label k
+) const
+{
+    bool iEnd = (i == 0 || i == density_.x());
+    bool jEnd = (j == 0 || j == density_.y());
+    bool kEnd = (k == 0 || k == density_.z());
+
+    return (iEnd && jEnd) || (iEnd && kEnd) || (jEnd && kEnd);
+}
+
+
+inline bool Foam::blockDescriptor::flatFaceOrEdge
+(
+    const label i, const label j, const label k
+) const
+{
+    if (i == 0 && curvedFaces_[0] == -1) return true;
+    if (i == density_.x() && curvedFaces_[1] == -1) return true;
+    if (j == 0 && curvedFaces_[2] == -1) return true;
+    if (j == density_.y() && curvedFaces_[3] == -1) return true;
+    if (k == 0 && curvedFaces_[4] == -1) return true;
+    if (k == density_.z() && curvedFaces_[5] == -1) return true;
+
+    bool iEnd = (i == 0 || i == density_.x());
+    bool jEnd = (j == 0 || j == density_.y());
+    bool kEnd = (k == 0 || k == density_.z());
+
+    return (iEnd && jEnd) || (iEnd && kEnd) || (jEnd && kEnd);
+}
+
+
 // ************************************************************************* //
diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptorTest.C b/src/mesh/blockMesh/blockDescriptor/blockDescriptorTest.C
new file mode 100644
index 0000000000..477e75d54a
--- /dev/null
+++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptorTest.C
@@ -0,0 +1,24 @@
+namespace Foam
+{
+    inline point faceCorr(const point& p)
+    {
+        return p/mag(p);
+    }
+}
+
+
+Foam::label Foam::blockDescriptor::correctFacePoints
+(
+    FixedList<List<point>, 6>& facePoints
+) const
+{
+    forAll(facePoints, facei)
+    {
+        forAll(facePoints[facei], pointi)
+        {
+            facePoints[facei][pointi] = faceCorr(facePoints[facei][pointi]);
+        }
+    }
+
+    return 6;
+}
diff --git a/src/mesh/blockMesh/blockEdges/BSpline.C b/src/mesh/blockMesh/blockEdges/BSplineEdge/BSpline.C
similarity index 99%
rename from src/mesh/blockMesh/blockEdges/BSpline.C
rename to src/mesh/blockMesh/blockEdges/BSplineEdge/BSpline.C
index 357b346b1c..c9e68e3c8c 100644
--- a/src/mesh/blockMesh/blockEdges/BSpline.C
+++ b/src/mesh/blockMesh/blockEdges/BSplineEdge/BSpline.C
@@ -23,7 +23,6 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "error.H"
 #include "BSpline.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
diff --git a/src/mesh/blockMesh/blockEdges/BSpline.H b/src/mesh/blockMesh/blockEdges/BSplineEdge/BSpline.H
similarity index 100%
rename from src/mesh/blockMesh/blockEdges/BSpline.H
rename to src/mesh/blockMesh/blockEdges/BSplineEdge/BSpline.H
diff --git a/src/mesh/blockMesh/blockEdges/BSplineEdge.C b/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.C
similarity index 95%
rename from src/mesh/blockMesh/blockEdges/BSplineEdge.C
rename to src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.C
index 8a89be7d51..206d5db957 100644
--- a/src/mesh/blockMesh/blockEdges/BSplineEdge.C
+++ b/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.C
@@ -57,7 +57,12 @@ Foam::BSplineEdge::BSplineEdge
 {}
 
 
-Foam::BSplineEdge::BSplineEdge(const pointField& points, Istream& is)
+Foam::BSplineEdge::BSplineEdge
+(
+    const searchableSurfaces& geometry,
+    const pointField& points,
+    Istream& is
+)
 :
     blockEdge(points, is),
     BSpline(appendEndPoints(points, start_, end_, pointField(is)))
diff --git a/src/mesh/blockMesh/blockEdges/BSplineEdge.H b/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.H
similarity index 95%
rename from src/mesh/blockMesh/blockEdges/BSplineEdge.H
rename to src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.H
index b860d9bd6d..8d11e15446 100644
--- a/src/mesh/blockMesh/blockEdges/BSplineEdge.H
+++ b/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.H
@@ -79,7 +79,12 @@ public:
         );
 
         //- Construct from Istream, setting pointsList
-        BSplineEdge(const pointField&, Istream&);
+        BSplineEdge
+        (
+            const searchableSurfaces& geometry,
+            const pointField&,
+            Istream&
+        );
 
 
     //- Destructor
diff --git a/src/mesh/blockMesh/blockEdges/arcEdge.C b/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.C
similarity index 97%
rename from src/mesh/blockMesh/blockEdges/arcEdge.C
rename to src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.C
index 5b8c327ee6..55597f2dea 100644
--- a/src/mesh/blockMesh/blockEdges/arcEdge.C
+++ b/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.C
@@ -118,7 +118,12 @@ Foam::arcEdge::arcEdge
 {}
 
 
-Foam::arcEdge::arcEdge(const pointField& points, Istream& is)
+Foam::arcEdge::arcEdge
+(
+    const searchableSurfaces& geometry,
+    const pointField& points,
+    Istream& is
+)
 :
     blockEdge(points, is),
     p1_(points_[start_]),
diff --git a/src/mesh/blockMesh/blockEdges/arcEdge.H b/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.H
similarity index 94%
rename from src/mesh/blockMesh/blockEdges/arcEdge.H
rename to src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.H
index ba89dbd5c7..c41ae669b8 100644
--- a/src/mesh/blockMesh/blockEdges/arcEdge.H
+++ b/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.H
@@ -69,6 +69,7 @@ class arcEdge
         //- Disallow default bitwise assignment
         void operator=(const arcEdge&);
 
+
 public:
 
     //- Runtime type information
@@ -86,11 +87,17 @@ public:
         );
 
         //- Construct from Istream setting pointsList
-        arcEdge(const pointField& points, Istream&);
+        arcEdge
+        (
+            const searchableSurfaces& geometry,
+            const pointField& points,
+            Istream&
+        );
 
 
     //- Destructor
-    virtual ~arcEdge(){}
+    virtual ~arcEdge()
+    {}
 
 
     // Member Functions
diff --git a/src/mesh/blockMesh/blockEdges/blockEdge.C b/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.C
similarity index 90%
rename from src/mesh/blockMesh/blockEdges/blockEdge.C
rename to src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.C
index d9f2f2059a..a2cc741822 100644
--- a/src/mesh/blockMesh/blockEdges/blockEdge.C
+++ b/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.C
@@ -23,7 +23,6 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "error.H"
 #include "blockEdge.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -50,7 +49,11 @@ Foam::blockEdge::blockEdge
 {}
 
 
-Foam::blockEdge::blockEdge(const pointField& points, Istream& is)
+Foam::blockEdge::blockEdge
+(
+    const pointField& points,
+    Istream& is
+)
 :
     points_(points),
     start_(readLabel(is)),
@@ -58,14 +61,6 @@ Foam::blockEdge::blockEdge(const pointField& points, Istream& is)
 {}
 
 
-Foam::blockEdge::blockEdge(const blockEdge& c)
-:
-    points_(c.points_),
-    start_(c.start_),
-    end_(c.end_)
-{}
-
-
 Foam::autoPtr<Foam::blockEdge> Foam::blockEdge::clone() const
 {
     NotImplemented;
@@ -75,6 +70,7 @@ Foam::autoPtr<Foam::blockEdge> Foam::blockEdge::clone() const
 
 Foam::autoPtr<Foam::blockEdge> Foam::blockEdge::New
 (
+    const searchableSurfaces& geometry,
     const pointField& points,
     Istream& is
 )
@@ -99,7 +95,7 @@ Foam::autoPtr<Foam::blockEdge> Foam::blockEdge::New
             << abort(FatalError);
     }
 
-    return autoPtr<blockEdge>(cstrIter()(points, is));
+    return autoPtr<blockEdge>(cstrIter()(geometry, points, is));
 }
 
 
@@ -131,12 +127,6 @@ Foam::pointField Foam::blockEdge::appendEndPoints
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
-void Foam::blockEdge::operator=(const blockEdge&)
-{
-    NotImplemented;
-}
-
-
 Foam::Ostream& Foam::operator<<(Ostream& os, const blockEdge& p)
 {
     os << p.start_ << tab << p.end_ << endl;
diff --git a/src/mesh/blockMesh/blockEdges/blockEdge.H b/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.H
similarity index 88%
rename from src/mesh/blockMesh/blockEdges/blockEdge.H
rename to src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.H
index 254eb09d83..24a64509ab 100644
--- a/src/mesh/blockMesh/blockEdges/blockEdge.H
+++ b/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.H
@@ -36,11 +36,7 @@ SourceFiles
 #ifndef blockEdges_H
 #define blockEdges_H
 
-#include "edge.H"
-#include "pointField.H"
-#include "typeInfo.H"
-#include "HashTable.H"
-#include "autoPtr.H"
+#include "searchableSurfaces.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -65,9 +61,11 @@ protected:
     // Protected data
 
         const pointField& points_;
+
         const label start_;
         const label end_;
 
+
     // Protected Member Functions
 
         //- Return a complete point field by appending the start/end points
@@ -94,10 +92,11 @@ public:
             blockEdge,
             Istream,
             (
+                const searchableSurfaces& geometry,
                 const pointField& points,
                 Istream& is
             ),
-            (points, is)
+            (geometry, points, is)
         );
 
 
@@ -112,39 +111,48 @@ public:
         );
 
         //- Construct from Istream setting pointsList
-        blockEdge(const pointField&, Istream&);
-
-        //- Copy construct
-        blockEdge(const blockEdge&);
+        blockEdge
+        (
+            const pointField&,
+            Istream&
+        );
 
         //- Clone function
         virtual autoPtr<blockEdge> clone() const;
 
         //- New function which constructs and returns pointer to a blockEdge
-        static autoPtr<blockEdge> New(const pointField&, Istream&);
+        static autoPtr<blockEdge> New
+        (
+            const searchableSurfaces& geometry,
+            const pointField&,
+            Istream&
+        );
 
         //- Class used for the read-construction of
         //  PtrLists of blockEdge
         class iNew
         {
+            const searchableSurfaces& geometry_;
             const pointField& points_;
 
         public:
 
-            iNew(const pointField& points)
+            iNew(const searchableSurfaces& geometry, const pointField& points)
             :
+                geometry_(geometry),
                 points_(points)
             {}
 
             autoPtr<blockEdge> operator()(Istream& is) const
             {
-                return blockEdge::New(points_, is);
+                return blockEdge::New(geometry_, points_, is);
             }
         };
 
 
     //- Destructor
-    virtual ~blockEdge(){}
+    virtual ~blockEdge()
+    {}
 
 
     // Member Functions
@@ -184,11 +192,6 @@ public:
         virtual scalar length() const = 0;
 
 
-    // Member operators
-
-        void operator=(const blockEdge&);
-
-
     // Ostream operator
 
         friend Ostream& operator<<(Ostream&, const blockEdge&);
diff --git a/src/mesh/blockMesh/blockEdges/blockEdgeI.H b/src/mesh/blockMesh/blockEdges/blockEdge/blockEdgeI.H
similarity index 100%
rename from src/mesh/blockMesh/blockEdges/blockEdgeI.H
rename to src/mesh/blockMesh/blockEdges/blockEdge/blockEdgeI.H
diff --git a/src/mesh/blockMesh/blockEdges/blockEdgeList.H b/src/mesh/blockMesh/blockEdges/blockEdge/blockEdgeList.H
similarity index 100%
rename from src/mesh/blockMesh/blockEdges/blockEdgeList.H
rename to src/mesh/blockMesh/blockEdges/blockEdge/blockEdgeList.H
diff --git a/src/mesh/blockMesh/blockEdges/lineDivide.C b/src/mesh/blockMesh/blockEdges/lineDivide/lineDivide.C
similarity index 100%
rename from src/mesh/blockMesh/blockEdges/lineDivide.C
rename to src/mesh/blockMesh/blockEdges/lineDivide/lineDivide.C
diff --git a/src/mesh/blockMesh/blockEdges/lineDivide.H b/src/mesh/blockMesh/blockEdges/lineDivide/lineDivide.H
similarity index 100%
rename from src/mesh/blockMesh/blockEdges/lineDivide.H
rename to src/mesh/blockMesh/blockEdges/lineDivide/lineDivide.H
diff --git a/src/mesh/blockMesh/blockEdges/lineEdge.C b/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.C
similarity index 95%
rename from src/mesh/blockMesh/blockEdges/lineEdge.C
rename to src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.C
index c43580d99d..5ceaafc5c5 100644
--- a/src/mesh/blockMesh/blockEdges/lineEdge.C
+++ b/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.C
@@ -23,7 +23,6 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "error.H"
 #include "lineEdge.H"
 #include "addToRunTimeSelectionTable.H"
 
@@ -49,7 +48,12 @@ Foam::lineEdge::lineEdge
 {}
 
 
-Foam::lineEdge::lineEdge(const pointField& points, Istream& is)
+Foam::lineEdge::lineEdge
+(
+    const searchableSurfaces& geometry,
+    const pointField& points,
+    Istream& is
+)
 :
     blockEdge(points, is)
 {}
diff --git a/src/mesh/blockMesh/blockEdges/lineEdge.H b/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.H
similarity index 93%
rename from src/mesh/blockMesh/blockEdges/lineEdge.H
rename to src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.H
index 90b0504feb..7eb5de7ec6 100644
--- a/src/mesh/blockMesh/blockEdges/lineEdge.H
+++ b/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.H
@@ -57,17 +57,23 @@ public:
     //- Runtime type information
     TypeName("line");
 
+
     // Constructors
 
         //- Construct from components
         lineEdge(const pointField&, const label start, const label end);
 
         //- Construct from Istream with a pointField
-        lineEdge(const pointField&, Istream&);
+        lineEdge
+        (
+            const searchableSurfaces& geometry,
+            const pointField&,
+            Istream&
+        );
 
 
-        //- Destructor
-        virtual ~lineEdge();
+    //- Destructor
+    virtual ~lineEdge();
 
 
     // Member Functions
diff --git a/src/mesh/blockMesh/blockEdges/polyLine.C b/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLine.C
similarity index 100%
rename from src/mesh/blockMesh/blockEdges/polyLine.C
rename to src/mesh/blockMesh/blockEdges/polyLineEdge/polyLine.C
diff --git a/src/mesh/blockMesh/blockEdges/polyLine.H b/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLine.H
similarity index 100%
rename from src/mesh/blockMesh/blockEdges/polyLine.H
rename to src/mesh/blockMesh/blockEdges/polyLineEdge/polyLine.H
diff --git a/src/mesh/blockMesh/blockEdges/polyLineEdge.C b/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.C
similarity index 95%
rename from src/mesh/blockMesh/blockEdges/polyLineEdge.C
rename to src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.C
index f8e7169f88..c79097ce20 100644
--- a/src/mesh/blockMesh/blockEdges/polyLineEdge.C
+++ b/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.C
@@ -23,7 +23,6 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "error.H"
 #include "polyLineEdge.H"
 #include "addToRunTimeSelectionTable.H"
 
@@ -51,7 +50,12 @@ Foam::polyLineEdge::polyLineEdge
 {}
 
 
-Foam::polyLineEdge::polyLineEdge(const pointField& ps, Istream& is)
+Foam::polyLineEdge::polyLineEdge
+(
+    const searchableSurfaces& geometry,
+    const pointField& ps,
+    Istream& is
+)
 :
     blockEdge(ps, is),
     polyLine(appendEndPoints(ps, start_, end_, pointField(is)))
diff --git a/src/mesh/blockMesh/blockEdges/polyLineEdge.H b/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.H
similarity index 95%
rename from src/mesh/blockMesh/blockEdges/polyLineEdge.H
rename to src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.H
index 94481d6d93..7458b393e8 100644
--- a/src/mesh/blockMesh/blockEdges/polyLineEdge.H
+++ b/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.H
@@ -79,7 +79,12 @@ public:
         );
 
         //- Construct from Istream
-        polyLineEdge(const pointField&, Istream&);
+        polyLineEdge
+        (
+            const searchableSurfaces& geometry,
+            const pointField&,
+            Istream&
+        );
 
 
         //- Destructor
diff --git a/src/mesh/blockMesh/blockEdges/CatmullRomSpline.C b/src/mesh/blockMesh/blockEdges/splineEdge/CatmullRomSpline.C
similarity index 99%
rename from src/mesh/blockMesh/blockEdges/CatmullRomSpline.C
rename to src/mesh/blockMesh/blockEdges/splineEdge/CatmullRomSpline.C
index 6c847abac8..7b1595d7b7 100644
--- a/src/mesh/blockMesh/blockEdges/CatmullRomSpline.C
+++ b/src/mesh/blockMesh/blockEdges/splineEdge/CatmullRomSpline.C
@@ -23,7 +23,6 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "error.H"
 #include "CatmullRomSpline.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
diff --git a/src/mesh/blockMesh/blockEdges/CatmullRomSpline.H b/src/mesh/blockMesh/blockEdges/splineEdge/CatmullRomSpline.H
similarity index 100%
rename from src/mesh/blockMesh/blockEdges/CatmullRomSpline.H
rename to src/mesh/blockMesh/blockEdges/splineEdge/CatmullRomSpline.H
diff --git a/src/mesh/blockMesh/blockEdges/splineEdge.C b/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.C
similarity index 95%
rename from src/mesh/blockMesh/blockEdges/splineEdge.C
rename to src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.C
index 7c9fa8142b..511220bffc 100644
--- a/src/mesh/blockMesh/blockEdges/splineEdge.C
+++ b/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.C
@@ -57,7 +57,12 @@ Foam::splineEdge::splineEdge
 {}
 
 
-Foam::splineEdge::splineEdge(const pointField& points, Istream& is)
+Foam::splineEdge::splineEdge
+(
+    const searchableSurfaces& geometry,
+    const pointField& points,
+    Istream& is
+)
 :
     blockEdge(points, is),
     CatmullRomSpline(appendEndPoints(points, start_, end_, pointField(is)))
diff --git a/src/mesh/blockMesh/blockEdges/splineEdge.H b/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.H
similarity index 95%
rename from src/mesh/blockMesh/blockEdges/splineEdge.H
rename to src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.H
index 4580d16241..cb9afb9032 100644
--- a/src/mesh/blockMesh/blockEdges/splineEdge.H
+++ b/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.H
@@ -79,7 +79,12 @@ public:
         );
 
         //- Construct from Istream, setting pointsList
-        splineEdge(const pointField&, Istream&);
+        splineEdge
+        (
+            const searchableSurfaces& geometry,
+            const pointField&,
+            Istream&
+        );
 
 
     //- Destructor
diff --git a/src/mesh/blockMesh/blockFaces/blockFace/blockFace.C b/src/mesh/blockMesh/blockFaces/blockFace/blockFace.C
new file mode 100644
index 0000000000..388c9fcc5d
--- /dev/null
+++ b/src/mesh/blockMesh/blockFaces/blockFace/blockFace.C
@@ -0,0 +1,98 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "blockFace.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(blockFace, 0);
+    defineRunTimeSelectionTable(blockFace, Istream);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::blockFace::blockFace(const face& vertices)
+:
+    vertices_(vertices)
+{}
+
+
+Foam::blockFace::blockFace(Istream& is)
+:
+    vertices_(is)
+{}
+
+
+Foam::autoPtr<Foam::blockFace> Foam::blockFace::clone() const
+{
+    NotImplemented;
+    return autoPtr<blockFace>(nullptr);
+}
+
+
+Foam::autoPtr<Foam::blockFace> Foam::blockFace::New
+(
+    const searchableSurfaces& geometry,
+    Istream& is
+)
+{
+    if (debug)
+    {
+        InfoInFunction << "Constructing blockFace" << endl;
+    }
+
+    const word faceType(is);
+
+    IstreamConstructorTable::iterator cstrIter =
+        IstreamConstructorTablePtr_->find(faceType);
+
+    if (cstrIter == IstreamConstructorTablePtr_->end())
+    {
+        FatalErrorInFunction
+            << "Unknown blockFace type "
+            << faceType << nl << nl
+            << "Valid blockFace types are" << endl
+            << IstreamConstructorTablePtr_->sortedToc()
+            << abort(FatalError);
+    }
+
+    return autoPtr<blockFace>(cstrIter()(geometry, is));
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const blockFace& p)
+{
+    os << p.vertices_ << endl;
+
+    return os;
+}
+
+
+// ************************************************************************* //
diff --git a/src/mesh/blockMesh/blockFaces/blockFace/blockFace.H b/src/mesh/blockMesh/blockFaces/blockFace/blockFace.H
new file mode 100644
index 0000000000..bbee0cb818
--- /dev/null
+++ b/src/mesh/blockMesh/blockFaces/blockFace/blockFace.H
@@ -0,0 +1,162 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::blockFace
+
+Description
+    Define a curved face.
+
+SourceFiles
+    blockFace.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef blockFaces_H
+#define blockFaces_H
+
+#include "searchableSurfaces.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of friend functions and operators
+
+class blockFace;
+
+Ostream& operator<<(Ostream&, const blockFace&);
+
+
+/*---------------------------------------------------------------------------*\
+                         Class blockFace Declaration
+\*---------------------------------------------------------------------------*/
+
+class blockFace
+{
+protected:
+
+    // Protected data
+
+        //- Block face vertices
+        const face vertices_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("blockFace");
+
+
+    // Declare run-time constructor selection tables
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            blockFace,
+            Istream,
+            (
+                const searchableSurfaces& geometry,
+                Istream& is
+            ),
+            (geometry, is)
+        );
+
+
+    // Constructors
+
+        //- Construct from face vertices
+        blockFace(const face& vertices);
+
+        //- Construct from Istream
+        blockFace(Istream&);
+
+        //- Clone function
+        virtual autoPtr<blockFace> clone() const;
+
+        //- New function which constructs and returns pointer to a blockFace
+        static autoPtr<blockFace> New
+        (
+            const searchableSurfaces& geometry,
+            Istream&
+        );
+
+        //- Class used for the read-construction of
+        //  PtrLists of blockFace
+        class iNew
+        {
+            const searchableSurfaces& geometry_;
+
+        public:
+
+            iNew(const searchableSurfaces& geometry)
+            :
+                geometry_(geometry)
+            {}
+
+            autoPtr<blockFace> operator()(Istream& is) const
+            {
+                return blockFace::New(geometry_, is);
+            }
+        };
+
+
+    //- Destructor
+    virtual ~blockFace()
+    {}
+
+
+    // Member Functions
+
+        //- Return block face vertices
+        inline const face& vertices() const;
+
+        //- Compare with given blockFace
+        inline bool compare(const blockFace&) const;
+
+        //- Compare with the given block and block face
+        inline bool compare(const face& vertices) const;
+
+        virtual void project(pointField& points) const = 0;
+
+
+    // Ostream operator
+
+        friend Ostream& operator<<(Ostream&, const blockFace&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "blockFaceI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/mesh/blockMesh/blockFaces/blockFace/blockFaceI.H b/src/mesh/blockMesh/blockFaces/blockFace/blockFaceI.H
new file mode 100644
index 0000000000..f96aeafb2a
--- /dev/null
+++ b/src/mesh/blockMesh/blockFaces/blockFace/blockFaceI.H
@@ -0,0 +1,46 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const Foam::face& Foam::blockFace::vertices() const
+{
+    return vertices_;
+}
+
+
+inline bool Foam::blockFace::compare(const blockFace& bf) const
+{
+    return compare(bf.vertices());
+}
+
+
+inline bool Foam::blockFace::compare(const face& vertices) const
+{
+    return face::sameVertices(vertices_, vertices);
+}
+
+
+// ************************************************************************* //
diff --git a/src/mesh/blockMesh/blockFaces/blockFace/blockFaceList.H b/src/mesh/blockMesh/blockFaces/blockFace/blockFaceList.H
new file mode 100644
index 0000000000..e9fb46909e
--- /dev/null
+++ b/src/mesh/blockMesh/blockFaces/blockFace/blockFaceList.H
@@ -0,0 +1,55 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Typedef
+    Foam::blockFaceList
+
+Description
+    A PtrList of blockFaces
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef blockFaceList_H
+#define blockFaceList_H
+
+#include "blockFace.H"
+#include "PtrList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+typedef PtrList<blockFace> blockFaceList;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/mesh/blockMesh/blockFaces/projectFace/projectFace.C b/src/mesh/blockMesh/blockFaces/projectFace/projectFace.C
new file mode 100644
index 0000000000..4876230c07
--- /dev/null
+++ b/src/mesh/blockMesh/blockFaces/projectFace/projectFace.C
@@ -0,0 +1,100 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "projectFace.H"
+#include "unitConversion.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(projectFace, 0);
+    addToRunTimeSelectionTable(blockFace, projectFace, Istream);
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+const Foam::searchableSurface& Foam::projectFace::lookupSurface
+(
+    const searchableSurfaces& geometry,
+    Istream& is
+) const
+{
+    word name(is);
+
+    forAll(geometry, i)
+    {
+        if (geometry[i].name() == name)
+        {
+            return geometry[i];
+        }
+    }
+
+    FatalIOErrorInFunction(is)
+        << "Cannot find surface " << name << " in geometry"
+        << exit(FatalIOError);
+
+    return geometry[0];
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::projectFace::projectFace
+(
+    const searchableSurfaces& geometry,
+    Istream& is
+)
+:
+    blockFace(is),
+    surface_(lookupSurface(geometry, is))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::projectFace::project(pointField& points) const
+{
+    List<pointIndexHit> hits;
+    scalarField nearestDistSqr
+    (
+        points.size(),
+        magSqr(points[0] - points[points.size()-1])
+    );
+    surface_.findNearest(points, nearestDistSqr, hits);
+
+    forAll(hits, i)
+    {
+        if (hits[i].hit())
+        {
+            points[i] = hits[i].hitPoint();
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/mesh/blockMesh/blockFaces/projectFace/projectFace.H b/src/mesh/blockMesh/blockFaces/projectFace/projectFace.H
new file mode 100644
index 0000000000..787f5fe152
--- /dev/null
+++ b/src/mesh/blockMesh/blockFaces/projectFace/projectFace.H
@@ -0,0 +1,111 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::projectFace
+
+Description
+    Projects the given set of face points onto the selected surface of the
+    geometry provided as a searchableSurfaces object.
+
+SourceFiles
+    projectFace.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef projectFace_H
+#define projectFace_H
+
+#include "blockFace.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class projectFace Declaration
+\*---------------------------------------------------------------------------*/
+
+class projectFace
+:
+    public blockFace
+{
+    // Private data
+
+        //- The surface onto which the points are projected
+        const searchableSurface& surface_;
+
+
+    // Private Member Functions
+
+        const searchableSurface& lookupSurface
+        (
+            const searchableSurfaces& geometry,
+            Istream& is
+        ) const;
+
+        //- Disallow default bitwise copy construct
+        projectFace(const projectFace&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const projectFace&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("project");
+
+
+    // Constructors
+
+        //- Construct from Istream setting pointsList
+        projectFace
+        (
+            const searchableSurfaces& geometry,
+            Istream&
+        );
+
+
+    //- Destructor
+    virtual ~projectFace()
+    {}
+
+
+    // Member Functions
+
+        //- Project the given points onto the surface
+        virtual void project(pointField& points) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/mesh/blockMesh/blockMesh/blockMesh.C b/src/mesh/blockMesh/blockMesh/blockMesh.C
index 90f75c8014..458b7ee2b6 100644
--- a/src/mesh/blockMesh/blockMesh/blockMesh.C
+++ b/src/mesh/blockMesh/blockMesh/blockMesh.C
@@ -24,12 +24,11 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "blockMesh.H"
+#include "Time.H"
 #include "Switch.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-bool Foam::blockMesh::blockMesh::verboseOutput(false);
-
 namespace Foam
 {
     defineDebugSwitch(blockMesh, 0);
@@ -40,6 +39,23 @@ namespace Foam
 
 Foam::blockMesh::blockMesh(const IOdictionary& dict, const word& regionName)
 :
+    verboseOutput(dict.lookupOrDefault<Switch>("verbose", true)),
+    geometry_
+    (
+        IOobject
+        (
+            "geometry",                 // dummy name
+            dict.time().constant(),     // instance
+            "geometry",                 // local
+            dict.time(),                // registry
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        ),
+        dict.found("geometry")
+      ? dict.subDict("geometry")
+      : dictionary(),
+        true
+    ),
     scaleFactor_(1.0),
     vertices_(dict.lookup("vertices")),
     topologyPtr_(createTopology(dict, regionName))
@@ -170,9 +186,9 @@ Foam::label Foam::blockMesh::numZonedBlocks() const
 {
     label num = 0;
 
-    forAll(*this, blockI)
+    forAll(*this, blocki)
     {
-        if (operator[](blockI).zoneName().size())
+        if (operator[](blocki).zoneName().size())
         {
             num++;
         }
diff --git a/src/mesh/blockMesh/blockMesh/blockMesh.H b/src/mesh/blockMesh/blockMesh/blockMesh.H
index 6416d2d58f..2b07036d3b 100644
--- a/src/mesh/blockMesh/blockMesh/blockMesh.H
+++ b/src/mesh/blockMesh/blockMesh/blockMesh.H
@@ -43,9 +43,11 @@ SourceFiles
 #define blockMesh_H
 
 #include "blockList.H"
+#include "searchableSurfaces.H"
 #include "polyMesh.H"
 #include "IOdictionary.H"
 #include "blockEdgeList.H"
+#include "blockFaceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -63,7 +65,10 @@ class blockMesh
     // Private data
 
         //- Switch for verbose output
-        static bool verboseOutput;
+        bool verboseOutput;
+
+        //- Optional searchable geometry to project face-points to
+        searchableSurfaces geometry_;
 
         //- The scaling factor to convert to metres
         scalar scaleFactor_;
@@ -74,6 +79,9 @@ class blockMesh
         //- The list of curved edges
         blockEdgeList edges_;
 
+        //- The list of curved faces
+        blockFaceList faces_;
+
         //- The blocks themselves (the topology) as a polyMesh
         polyMesh* topologyPtr_;
 
@@ -166,6 +174,12 @@ public:
 
         // Access
 
+            //- Optional searchable geometry to project face-points to
+            const searchableSurfaces& geometry() const
+            {
+                return geometry_;
+            }
+
             //- Reference to point field defining the blockMesh
             //  these points have not been scaled by scaleFactor
             const pointField& vertices() const;
@@ -179,6 +193,12 @@ public:
                 return edges_;
             }
 
+            //- Return the curved faces
+            const blockFaceList& faces() const
+            {
+                return faces_;
+            }
+
             //- The scaling factor used to convert to metres
             scalar scaleFactor() const;
 
@@ -204,11 +224,8 @@ public:
 
         // Edit
 
-            //- Clear geometry (internal points, cells, boundaryPatches)
-            void clearGeom();
-
             //- Enable/disable verbose information about the progress
-            static void verbose(const bool on=true);
+            void verbose(const bool on=true);
 
 
         // Write
diff --git a/src/mesh/blockMesh/blockMesh/blockMeshCreate.C b/src/mesh/blockMesh/blockMesh/blockMeshCreate.C
index d16f4d3674..4ace695ad1 100644
--- a/src/mesh/blockMesh/blockMesh/blockMeshCreate.C
+++ b/src/mesh/blockMesh/blockMesh/blockMeshCreate.C
@@ -37,41 +37,37 @@ void Foam::blockMesh::createPoints() const
         Info<< "Creating points with scale " << scaleFactor_ << endl;
     }
 
-    //
-    // generate points
-    //
-    points_.clear();
     points_.setSize(nPoints_);
 
-    forAll(blocks, blockI)
+    forAll(blocks, blocki)
     {
-        const pointField& blockPoints = blocks[blockI].points();
+        const pointField& blockPoints = blocks[blocki].points();
 
         if (verboseOutput)
         {
-            const Vector<label>& density = blocks[blockI].density();
+            const Vector<label>& density = blocks[blocki].density();
 
-            label v0 = blocks[blockI].vtxLabel(0, 0, 0);
-            label vi1 = blocks[blockI].vtxLabel(1, 0, 0);
-            scalar diStart = mag(blockPoints[vi1]-blockPoints[v0]);
+            label v0 = blocks[blocki].pointLabel(0, 0, 0);
+            label vi1 = blocks[blocki].pointLabel(1, 0, 0);
+            scalar diStart = mag(blockPoints[vi1] - blockPoints[v0]);
 
-            label vinM1 = blocks[blockI].vtxLabel(density.x()-1, 0, 0);
-            label vin = blocks[blockI].vtxLabel(density.x(), 0, 0);
-            scalar diFinal = mag(blockPoints[vin]-blockPoints[vinM1]);
+            label vinM1 = blocks[blocki].pointLabel(density.x()-1, 0, 0);
+            label vin = blocks[blocki].pointLabel(density.x(), 0, 0);
+            scalar diFinal = mag(blockPoints[vin] - blockPoints[vinM1]);
 
-            label vj1 = blocks[blockI].vtxLabel(0, 1, 0);
-            scalar djStart = mag(blockPoints[vj1]-blockPoints[v0]);
-            label vjnM1 = blocks[blockI].vtxLabel(0, density.y()-1, 0);
-            label vjn = blocks[blockI].vtxLabel(0, density.y(), 0);
-            scalar djFinal = mag(blockPoints[vjn]-blockPoints[vjnM1]);
+            label vj1 = blocks[blocki].pointLabel(0, 1, 0);
+            scalar djStart = mag(blockPoints[vj1] - blockPoints[v0]);
+            label vjnM1 = blocks[blocki].pointLabel(0, density.y()-1, 0);
+            label vjn = blocks[blocki].pointLabel(0, density.y(), 0);
+            scalar djFinal = mag(blockPoints[vjn] - blockPoints[vjnM1]);
 
-            label vk1 = blocks[blockI].vtxLabel(0, 0, 1);
-            scalar dkStart = mag(blockPoints[vk1]-blockPoints[v0]);
-            label vknM1 = blocks[blockI].vtxLabel(0, 0, density.z()-1);
-            label vkn = blocks[blockI].vtxLabel(0, 0, density.z());
-            scalar dkFinal = mag(blockPoints[vkn]-blockPoints[vknM1]);
+            label vk1 = blocks[blocki].pointLabel(0, 0, 1);
+            scalar dkStart = mag(blockPoints[vk1] - blockPoints[v0]);
+            label vknM1 = blocks[blocki].pointLabel(0, 0, density.z()-1);
+            label vkn = blocks[blocki].pointLabel(0, 0, density.z());
+            scalar dkFinal = mag(blockPoints[vkn] - blockPoints[vknM1]);
 
-            Info<< "    Block " << blockI << " cell size :" << nl
+            Info<< "    Block " << blocki << " cell size :" << nl
                 << "        i : " << scaleFactor_*diStart << " .. "
                 << scaleFactor_*diFinal << nl
                 << "        j : " << scaleFactor_*djStart << " .. "
@@ -87,7 +83,7 @@ void Foam::blockMesh::createPoints() const
             [
                 mergeList_
                 [
-                    blockOffsets_[blockI] + blockPointi
+                    blockOffsets_[blocki] + blockPointi
                 ]
             ] = scaleFactor_ * blockPoints[blockPointi];
         }
@@ -105,17 +101,13 @@ void Foam::blockMesh::createCells() const
         Info<< "Creating cells" << endl;
     }
 
-    //
-    // generate cells
-    //
-    cells_.clear();
     cells_.setSize(nCells_);
 
     label cellLabel = 0;
 
-    forAll(blocks, blockI)
+    forAll(blocks, blocki)
     {
-        const labelListList& blockCells = blocks[blockI].cells();
+        const List<FixedList<label, 8>> blockCells(blocks[blocki].cells());
 
         forAll(blockCells, blockCelli)
         {
@@ -127,7 +119,7 @@ void Foam::blockMesh::createCells() const
                     mergeList_
                     [
                         blockCells[blockCelli][cellPointi]
-                      + blockOffsets_[blockI]
+                      + blockOffsets_[blocki]
                     ];
             }
 
@@ -153,9 +145,9 @@ Foam::faceList Foam::blockMesh::createPatchFaces
 
     forAll(patchTopologyFaces, patchTopologyFaceLabel)
     {
-        const label blockI = blockLabels[patchTopologyFaceLabel];
+        const label blocki = blockLabels[patchTopologyFaceLabel];
 
-        faceList blockFaces = blocks[blockI].blockShape().faces();
+        faceList blockFaces = blocks[blocki].blockShape().faces();
 
         forAll(blockFaces, blockFaceLabel)
         {
@@ -166,7 +158,7 @@ Foam::faceList Foam::blockMesh::createPatchFaces
             )
             {
                 nFaces +=
-                    blocks[blockI].boundaryPatches()[blockFaceLabel].size();
+                    blocks[blocki].boundaryPatches()[blockFaceLabel].size();
             }
         }
     }
@@ -178,20 +170,20 @@ Foam::faceList Foam::blockMesh::createPatchFaces
 
     forAll(patchTopologyFaces, patchTopologyFaceLabel)
     {
-        const label blockI = blockLabels[patchTopologyFaceLabel];
+        const label blocki = blockLabels[patchTopologyFaceLabel];
 
-        faceList blockFaces = blocks[blockI].blockShape().faces();
+        faceList blockFaces = blocks[blocki].blockShape().faces();
 
         forAll(blockFaces, blockFaceLabel)
         {
             if
             (
                 blockFaces[blockFaceLabel]
-                == patchTopologyFaces[patchTopologyFaceLabel]
+             == patchTopologyFaces[patchTopologyFaceLabel]
             )
             {
-                const labelListList& blockPatchFaces =
-                    blocks[blockI].boundaryPatches()[blockFaceLabel];
+                const List<FixedList<label, 4>>& blockPatchFaces =
+                    blocks[blocki].boundaryPatches()[blockFaceLabel];
 
                 forAll(blockPatchFaces, blockFaceLabel)
                 {
@@ -202,7 +194,7 @@ Foam::faceList Foam::blockMesh::createPatchFaces
                         mergeList_
                         [
                             blockPatchFaces[blockFaceLabel][0]
-                          + blockOffsets_[blockI]
+                          + blockOffsets_[blocki]
                         ];
 
                     label nUnique = 1;
@@ -218,7 +210,7 @@ Foam::faceList Foam::blockMesh::createPatchFaces
                             mergeList_
                             [
                                 blockPatchFaces[blockFaceLabel][facePointLabel]
-                              + blockOffsets_[blockI]
+                              + blockOffsets_[blocki]
                             ];
 
                         if (quadFace[nUnique] != quadFace[nUnique-1])
@@ -264,29 +256,13 @@ void Foam::blockMesh::createPatches() const
         Info<< "Creating patches" << endl;
     }
 
-    //
-    // generate points
-    //
-
-    patches_.clear();
     patches_.setSize(topoPatches.size());
 
     forAll(topoPatches, patchi)
     {
         patches_[patchi] = createPatchFaces(topoPatches[patchi]);
     }
-
 }
 
 
-void Foam::blockMesh::clearGeom()
-{
-    blockList& blocks = *this;
-
-    forAll(blocks, blockI)
-    {
-        blocks[blockI].clearGeom();
-    }
-}
-
 // ************************************************************************* //
diff --git a/src/mesh/blockMesh/blockMesh/blockMeshMerge.C b/src/mesh/blockMesh/blockMesh/blockMeshMerge.C
index edac5b9e49..79a42caec2 100644
--- a/src/mesh/blockMesh/blockMesh/blockMeshMerge.C
+++ b/src/mesh/blockMesh/blockMesh/blockMeshMerge.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -41,12 +41,12 @@ void Foam::blockMesh::calcMergeInfo()
     nPoints_ = 0;
     nCells_  = 0;
 
-    forAll(blocks, blockI)
+    forAll(blocks, blocki)
     {
-        blockOffsets_[blockI] = nPoints_;
+        blockOffsets_[blocki] = nPoints_;
 
-        nPoints_ += blocks[blockI].nPoints();
-        nCells_  += blocks[blockI].nCells();
+        nPoints_ += blocks[blocki].nPoints();
+        nCells_  += blocks[blocki].nCells();
     }
 
 
@@ -100,7 +100,7 @@ void Foam::blockMesh::calcMergeInfo()
                 << exit(FatalError);
         }
 
-        const labelListList& blockPfaceFaces =
+        const List<FixedList<label, 4>>& blockPfaceFaces =
             blocks[blockPlabel].boundaryPatches()[blockPfaceLabel];
 
         labelListList& curPairs = glueMergePairs[blockFaceLabel];
@@ -122,7 +122,7 @@ void Foam::blockMesh::calcMergeInfo()
 
         forAll(blockPfaceFaces, blockPfaceFaceLabel)
         {
-            const labelList& blockPfaceFacePoints
+            const FixedList<label, 4>& blockPfaceFacePoints
                 = blockPfaceFaces[blockPfaceFaceLabel];
 
             forAll(blockPfaceFacePoints, blockPfaceFacePointLabel)
@@ -209,7 +209,7 @@ void Foam::blockMesh::calcMergeInfo()
                 << exit(FatalError);
         }
 
-        const labelListList& blockNfaceFaces =
+        const List<FixedList<label, 4>>& blockNfaceFaces =
             blocks[blockNlabel].boundaryPatches()[blockNfaceLabel];
 
         if (blockPfaceFaces.size() != blockNfaceFaces.size())
@@ -225,7 +225,7 @@ void Foam::blockMesh::calcMergeInfo()
         // master block over all point of all faces of slave block
         forAll(blockPfaceFaces, blockPfaceFaceLabel)
         {
-            const labelList& blockPfaceFacePoints
+            const FixedList<label, 4>& blockPfaceFacePoints
                 = blockPfaceFaces[blockPfaceFaceLabel];
 
             labelList& cp = curPairs[blockPfaceFaceLabel];
@@ -236,7 +236,7 @@ void Foam::blockMesh::calcMergeInfo()
             {
                 forAll(blockNfaceFaces, blockNfaceFaceLabel)
                 {
-                    const labelList& blockNfaceFacePoints
+                    const FixedList<label, 4>& blockNfaceFacePoints
                         = blockNfaceFaces[blockNfaceFaceLabel];
 
                     forAll(blockNfaceFacePoints, blockNfaceFacePointLabel)
@@ -299,7 +299,7 @@ void Foam::blockMesh::calcMergeInfo()
     }
 
 
-    const faceList::subList blockInternalFaces
+    const faceList::subList blockinternalFaces
     (
         blockFaces,
         topology().nInternalFaces()
@@ -313,7 +313,7 @@ void Foam::blockMesh::calcMergeInfo()
         changedPointMerge = false;
         nPasses++;
 
-        forAll(blockInternalFaces, blockFaceLabel)
+        forAll(blockinternalFaces, blockFaceLabel)
         {
             label blockPlabel = faceOwnerBlocks[blockFaceLabel];
             label blockNlabel = faceNeighbourBlocks[blockFaceLabel];
@@ -334,7 +334,7 @@ void Foam::blockMesh::calcMergeInfo()
                 if
                 (
                     blockFaces[blockPfaces[blockPfaceLabel]]
-                 == blockInternalFaces[blockFaceLabel]
+                 == blockinternalFaces[blockFaceLabel]
                 )
                 {
                     break;
@@ -353,7 +353,7 @@ void Foam::blockMesh::calcMergeInfo()
                 if
                 (
                     blockFaces[blockNfaces[blockNfaceLabel]]
-                 == blockInternalFaces[blockFaceLabel]
+                 == blockinternalFaces[blockFaceLabel]
                 )
                 {
                     break;
@@ -361,12 +361,12 @@ void Foam::blockMesh::calcMergeInfo()
             }
 
 
-            const labelListList& blockPfaceFaces =
+            const List<FixedList<label, 4>>& blockPfaceFaces =
                 blocks[blockPlabel].boundaryPatches()[blockPfaceLabel];
 
             forAll(blockPfaceFaces, blockPfaceFaceLabel)
             {
-                const labelList& blockPfaceFacePoints
+                const FixedList<label, 4>& blockPfaceFacePoints
                     = blockPfaceFaces[blockPfaceFaceLabel];
 
                 const labelList& cp = curPairs[blockPfaceFaceLabel];
@@ -419,7 +419,7 @@ void Foam::blockMesh::calcMergeInfo()
         Info<< endl;
     }
 
-    forAll(blockInternalFaces, blockFaceLabel)
+    forAll(blockinternalFaces, blockFaceLabel)
     {
         label blockPlabel = faceOwnerBlocks[blockFaceLabel];
         label blockNlabel = faceNeighbourBlocks[blockFaceLabel];
@@ -442,7 +442,7 @@ void Foam::blockMesh::calcMergeInfo()
             if
             (
                 blockFaces[blockPfaces[blockPfaceLabel]]
-             == blockInternalFaces[blockFaceLabel]
+             == blockinternalFaces[blockFaceLabel]
             )
             {
                 foundFace = true;
@@ -469,7 +469,7 @@ void Foam::blockMesh::calcMergeInfo()
             if
             (
                 blockFaces[blockNfaces[blockNfaceLabel]]
-             == blockInternalFaces[blockFaceLabel]
+             == blockinternalFaces[blockFaceLabel]
             )
             {
                 foundFace = true;
@@ -484,15 +484,15 @@ void Foam::blockMesh::calcMergeInfo()
                 << exit(FatalError);
         }
 
-        const labelListList& blockPfaceFaces =
+        const List<FixedList<label, 4>>& blockPfaceFaces =
             blocks[blockPlabel].boundaryPatches()[blockPfaceLabel];
 
-        const labelListList& blockNfaceFaces =
+        const List<FixedList<label, 4>>& blockNfaceFaces =
             blocks[blockNlabel].boundaryPatches()[blockNfaceLabel];
 
         forAll(blockPfaceFaces, blockPfaceFaceLabel)
         {
-            const labelList& blockPfaceFacePoints
+            const FixedList<label, 4>& blockPfaceFacePoints
                 = blockPfaceFaces[blockPfaceFaceLabel];
 
             forAll(blockPfaceFacePoints, blockPfaceFacePointLabel)
@@ -518,7 +518,7 @@ void Foam::blockMesh::calcMergeInfo()
 
         forAll(blockNfaceFaces, blockNfaceFaceLabel)
         {
-            const labelList& blockNfaceFacePoints
+            const FixedList<label, 4>& blockNfaceFacePoints
                 = blockNfaceFaces[blockNfaceFaceLabel];
 
             forAll(blockNfaceFacePoints, blockNfaceFacePointLabel)
diff --git a/src/mesh/blockMesh/blockMesh/blockMeshMergeFast.C b/src/mesh/blockMesh/blockMesh/blockMeshMergeFast.C
index 2ad5892f97..7ab1b1d167 100644
--- a/src/mesh/blockMesh/blockMesh/blockMeshMergeFast.C
+++ b/src/mesh/blockMesh/blockMesh/blockMeshMergeFast.C
@@ -209,17 +209,17 @@ inline label facePoint
     switch (facei)
     {
         case 0:
-            return block.vtxLabel(0, i, j);
+            return block.pointLabel(0, i, j);
         case 1:
-            return block.vtxLabel(block.density().x(), i, j);
+            return block.pointLabel(block.density().x(), i, j);
         case 2:
-            return block.vtxLabel(i, 0, j);
+            return block.pointLabel(i, 0, j);
         case 3:
-            return block.vtxLabel(i, block.density().y(), j);
+            return block.pointLabel(i, block.density().y(), j);
         case 4:
-            return block.vtxLabel(i, j, 0);
+            return block.pointLabel(i, j, 0);
         case 5:
-            return block.vtxLabel(i, j, block.density().z());
+            return block.pointLabel(i, j, block.density().z());
         default:
             return -1;
     }
@@ -234,7 +234,7 @@ inline label facePointN
     const label k
 )
 {
-    return block.vtxLabel
+    return block.pointLabel
     (
         unsignIndex(i, block.density().x()),
         unsignIndex(j, block.density().y()),
@@ -305,12 +305,12 @@ void Foam::blockMesh::calcMergeInfoFast()
     nPoints_ = 0;
     nCells_  = 0;
 
-    forAll(blocks, blockI)
+    forAll(blocks, blocki)
     {
-        blockOffsets_[blockI] = nPoints_;
+        blockOffsets_[blocki] = nPoints_;
 
-        nPoints_ += blocks[blockI].nPoints();
-        nCells_  += blocks[blockI].nCells();
+        nPoints_ += blocks[blocki].nPoints();
+        nCells_  += blocks[blocki].nCells();
     }
 
     if (verboseOutput)
diff --git a/src/mesh/blockMesh/blockMesh/blockMeshTopology.C b/src/mesh/blockMesh/blockMesh/blockMeshTopology.C
index 1044803811..aaa76ac22d 100644
--- a/src/mesh/blockMesh/blockMesh/blockMeshTopology.C
+++ b/src/mesh/blockMesh/blockMesh/blockMeshTopology.C
@@ -307,9 +307,9 @@ void Foam::blockMesh::createCellShapes
     const blockMesh& blocks = *this;
 
     tmpBlockCells.setSize(blocks.size());
-    forAll(blocks, blockI)
+    forAll(blocks, blocki)
     {
-        tmpBlockCells[blockI] = blocks[blockI].blockShape();
+        tmpBlockCells[blocki] = blocks[blocki].blockShape();
     }
 }
 
@@ -354,7 +354,7 @@ Foam::polyMesh* Foam::blockMesh::createTopology
         blockEdgeList edges
         (
             meshDescription.lookup("edges"),
-            blockEdge::iNew(vertices_)
+            blockEdge::iNew(geometry_, vertices_)
         );
 
         edges_.transfer(edges);
@@ -365,6 +365,28 @@ Foam::polyMesh* Foam::blockMesh::createTopology
     }
 
 
+    // Read the block faces
+    if (meshDescription.found("faces"))
+    {
+        if (verboseOutput)
+        {
+            Info<< "Creating block faces" << endl;
+        }
+
+        blockFaceList faces
+        (
+            meshDescription.lookup("faces"),
+            blockFace::iNew(geometry_)
+        );
+
+        faces_.transfer(faces);
+    }
+    else if (verboseOutput)
+    {
+        Info<< "No non-planar block faces defined" << endl;
+    }
+
+
     // Create the blocks
     if (verboseOutput)
     {
@@ -374,12 +396,14 @@ Foam::polyMesh* Foam::blockMesh::createTopology
         blockList blocks
         (
             meshDescription.lookup("blocks"),
-            block::iNew(vertices_, edges_)
+            block::iNew(vertices_, edges_, faces_)
         );
+
         transfer(blocks);
     }
 
 
+
     polyMesh* blockMeshPtr = nullptr;
 
     // Create the patches
diff --git a/src/meshTools/searchableSurface/searchableSurfaceCollection.H b/src/meshTools/searchableSurface/searchableSurfaceCollection.H
index 66a9e81228..6026ebc87c 100644
--- a/src/meshTools/searchableSurface/searchableSurfaceCollection.H
+++ b/src/meshTools/searchableSurface/searchableSurfaceCollection.H
@@ -47,8 +47,6 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declaration of classes
-
 /*---------------------------------------------------------------------------*\
                  Class searchableSurfaceCollection Declaration
 \*---------------------------------------------------------------------------*/
@@ -57,8 +55,6 @@ class searchableSurfaceCollection
 :
     public searchableSurface
 {
-private:
-
     // Private Member Data
 
         // Per instance data
@@ -82,6 +78,7 @@ private:
 
         //- Region names
         mutable wordList regions_;
+
         //- From individual regions to collection regions
         mutable labelList regionOffset_;
 
@@ -133,6 +130,7 @@ public:
             const dictionary& dict
         );
 
+
     //- Destructor
     virtual ~searchableSurfaceCollection();
 
@@ -163,7 +161,6 @@ public:
             return transform_;
         }
 
-
         virtual const wordList& regions() const;
 
         //- Whether supports volume type below
@@ -281,7 +278,6 @@ public:
                 NotImplemented;
                 return false;
             }
-
 };
 
 
diff --git a/src/meshTools/searchableSurface/searchableSurfaces.H b/src/meshTools/searchableSurface/searchableSurfaces.H
index 9c85e3f042..b367a2a0d2 100644
--- a/src/meshTools/searchableSurface/searchableSurfaces.H
+++ b/src/meshTools/searchableSurface/searchableSurfaces.H
@@ -92,6 +92,7 @@ public:
 
     ClassName("searchableSurfaces");
 
+
     // Constructors
 
         //- Construct with length specified. Fill later.
@@ -119,6 +120,7 @@ public:
         {
             return names_;
         }
+
         wordList& names()
         {
             return names_;
@@ -128,6 +130,7 @@ public:
         {
             return regionNames_;
         }
+
         List<wordList>& regionNames()
         {
             return regionNames_;
diff --git a/tutorials/mesh/blockMesh/sphere/Allrun b/tutorials/mesh/blockMesh/sphere/Allrun
new file mode 100755
index 0000000000..3c4ef3e7ff
--- /dev/null
+++ b/tutorials/mesh/blockMesh/sphere/Allrun
@@ -0,0 +1,9 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+runApplication blockMesh
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/mesh/blockMesh/sphere/system/blockMeshDict b/tutorials/mesh/blockMesh/sphere/system/blockMeshDict
new file mode 100644
index 0000000000..56a52367f9
--- /dev/null
+++ b/tutorials/mesh/blockMesh/sphere/system/blockMeshDict
@@ -0,0 +1,97 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+geometry
+{
+    sphere
+    {
+        type searchableSphere;
+        centre (0 0 0);
+        radius 1;
+    }
+}
+
+scale 1;
+
+v 0.5773502;
+mv -0.5773502;
+
+a 0.7071067;
+ma -0.7071067;
+
+vertices
+(
+    ($mv $mv $mv)
+    ( $v $mv $mv)
+    ( $v  $v $mv)
+    ($mv  $v $mv)
+    ($mv $mv  $v)
+    ( $v $mv  $v)
+    ( $v  $v  $v)
+    ($mv  $v  $v)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (10 10 10) simpleGrading (1 1 1)
+);
+
+edges
+(
+    arc 0 1 (0 $ma $ma)
+    arc 2 3 (0 $a $ma)
+    arc 6 7 (0 $a $a)
+    arc 4 5 (0 $ma $a)
+
+    arc 0 3 ($ma 0 $ma)
+    arc 1 2 ($a 0 $ma)
+    arc 5 6 ($a 0 $a)
+    arc 4 7 ($ma 0 $a)
+
+    arc 0 4 ($ma $ma 0)
+    arc 1 5 ($a $ma 0)
+    arc 2 6 ($a $a 0)
+    arc 3 7 ($ma $a 0)
+);
+
+faces
+(
+    project (0 4 7 3) sphere
+    project (2 6 5 1) sphere
+    project (1 5 4 0) sphere
+    project (3 7 6 2) sphere
+    project (0 3 2 1) sphere
+    project (4 5 6 7) sphere
+);
+
+boundary
+(
+    walls
+    {
+        type wall;
+        faces
+        (
+            (0 4 7 3)
+            (2 6 5 1)
+            (1 5 4 0)
+            (3 7 6 2)
+            (0 3 2 1)
+            (4 5 6 7)
+        );
+    }
+);
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/blockMesh/sphere/system/controlDict b/tutorials/mesh/blockMesh/sphere/system/controlDict
new file mode 100644
index 0000000000..6b63048145
--- /dev/null
+++ b/tutorials/mesh/blockMesh/sphere/system/controlDict
@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     blockMesh;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         0;
+
+deltaT          0;
+
+writeControl    timeStep;
+
+writeInterval   1;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable true;
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/blockMesh/sphere/system/fvSchemes b/tutorials/mesh/blockMesh/sphere/system/fvSchemes
new file mode 100644
index 0000000000..606b28ce04
--- /dev/null
+++ b/tutorials/mesh/blockMesh/sphere/system/fvSchemes
@@ -0,0 +1,37 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{}
+
+gradSchemes
+{}
+
+divSchemes
+{}
+
+laplacianSchemes
+{}
+
+interpolationSchemes
+{}
+
+snGradSchemes
+{}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/blockMesh/sphere/system/fvSolution b/tutorials/mesh/blockMesh/sphere/system/fvSolution
new file mode 100644
index 0000000000..6db4e1a6ba
--- /dev/null
+++ b/tutorials/mesh/blockMesh/sphere/system/fvSolution
@@ -0,0 +1,18 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// ************************************************************************* //