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.

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

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 List<FixedList<label, 8>> blockCells = b.cells();
-            const word& zoneName = b.blockDef().zoneName();
+            const word& zoneName = b.zoneName();
 
             if (zoneName.size())
             {

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);
 

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();
 

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)

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);
 

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);
 

src/mesh/blockMesh/Make/files

@@ -1,14 +1,16 @@
-blockEdges/BSpline.C
+blockEdges/blockEdge/blockEdge.C
-blockEdges/CatmullRomSpline.C
+blockEdges/lineDivide/lineDivide.C
-blockEdges/polyLine.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
+blockFaces/blockFace/blockFace.C
-blockEdges/blockEdge.C
+blockFaces/projectFace/projectFace.C
-blockEdges/lineEdge.C
-blockEdges/polyLineEdge.C
-blockEdges/lineDivide.C
-blockEdges/BSplineEdge.C
-blockEdges/splineEdge.C
 
 gradingDescriptor/gradingDescriptor.C
 gradingDescriptor/gradingDescriptors.C

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

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),
+    blockDescriptor(vertices, edges, faces, is)
-    points_(0),
+{
-    cells_(0),
+    createPoints();
-    boundaryPatches_(0)
+    createBoundary();
-{}
+}
 
 
 Foam::block::block(const blockDescriptor& blockDesc)
 :
-    blockDescriptor(blockDesc),
+    blockDescriptor(blockDesc)
-    points_(0),
-    cells_(0),
-    boundaryPatches_(0)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::block::~block()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-const Foam::pointField& Foam::block::points() const
 {
-    if (points_.empty())
+    createPoints();
-    {
+    createBoundary();
-        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_;
 }
 
 

src/mesh/blockMesh/block/block.H

@@ -68,25 +68,19 @@ class block
     // Private data
 
         //- List of points
-        mutable pointField points_;
+        pointField points_;
-
-        //- List of cells
-        mutable labelListList cells_;
 
         //- 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;
+        void createPoints();
-
-        //- Creates cells for filling the block
-        void createCells() const;
 
         //- 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;
+            inline const FixedList<List<FixedList<label, 4>>, 6>&
-
+                boundaryPatches() const;
-
-        // Edit
-
-            //- Clear geometry (internal points, cells, boundaryPatches)
-            void clearGeom();
 
 
     // Ostream Operator

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++)
+    points_[pointLabel(0,  0,  0)] = p000;
-    {
+    points_[pointLabel(ni, 0,  0)] = p100;
-        for (label j = 0; j <= nj; j++)
+    points_[pointLabel(ni, nj, 0)] = p110;
-        {
+    points_[pointLabel(0,  nj, 0)] = p010;
-            for (label i = 0; i <= ni; i++)
+    points_[pointLabel(0,  0,  nk)] = p001;
-            {
+    points_[pointLabel(ni, 0,  nk)] = p101;
-                const label vertexNo = vtxLabel(i, j, k);
+    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 =
+                scalar wx1 = (1 - w0)*(1 - w4)*(1 - w8) + w0*(1 - w5)*(1 - w9);
-                (
+                scalar wx2 = (1 - w1)*w4*(1 - w11)      + w1*w5*(1 - w10);
-                    (1 - w[0][i])*(1 - w[4][j])*(1 - w[8][k])
+                scalar wx3 = (1 - w2)*w7*w11            + w2*w6*w10;
-                  + w[0][i]*(1 - w[5][j])*(1 - w[9][k])
+                scalar wx4 = (1 - w3)*(1 - w7)*w8       + w3*(1 - w6)*w9;
-                );
 
-                scalar impx2 =
+                const scalar sumWx = wx1 + wx2 + wx3 + wx4;
-                (
+                wx1 /= sumWx;
-                    (1 - w[1][i])*w[4][j]*(1 - w[11][k])
+                wx2 /= sumWx;
-                  + w[1][i]*w[5][j]*(1 - w[10][k])
+                wx3 /= sumWx;
-                );
+                wx4 /= sumWx;
-
-                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;
 
 
                 // y-direction
-                scalar impy1 =
+                scalar wy1 = (1 - w4)*(1 - w0)*(1 - w8) + w4*(1 - w1)*(1 - w11);
-                (
+                scalar wy2 = (1 - w5)*w0*(1 - w9)       + w5*w1*(1 - w10);
-                    (1 - w[4][j])*(1 - w[0][i])*(1 - w[8][k])
+                scalar wy3 = (1 - w6)*w3*w9             + w6*w2*w10;
-                  + w[4][j]*(1 - w[1][i])*(1 - w[11][k])
+                scalar wy4 = (1 - w7)*(1 - w3)*w8       + w7*(1 - w2)*w11;
-                );
 
-                scalar impy2 =
+                const scalar sumWy = wy1 + wy2 + wy3 + wy4;
-                (
+                wy1 /= sumWy;
-                    (1 - w[5][j])*w[0][i]*(1 - w[9][k])
+                wy2 /= sumWy;
-                  + w[5][j]*w[1][i]*(1 - w[10][k])
+                wy3 /= sumWy;
-                );
+                wy4 /= sumWy;
-
-                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;
 
 
                 // z-direction
-                scalar impz1 =
+                scalar wz1 = (1 - w8)*(1 - w0)*(1 - w4) + w8*(1 - w3)*(1 - w7);
-                (
+                scalar wz2 = (1 - w9)*w0*(1 - w5)       + w9*w3*(1 - w6);
-                    (1 - w[8][k])*(1 - w[0][i])*(1 - w[4][j])
+                scalar wz3 = (1 - w10)*w1*w5            + w10*w2*w6;
-                  + w[8][k]*(1 - w[3][i])*(1 - w[7][j])
+                scalar wz4 = (1 - w11)*(1 - w1)*w4      + w11*(1 - w2)*w7;
-                );
 
-                scalar impz2 =
+                const scalar sumWz = wz1 + wz2 + wz3 + wz4;
-                (
+                wz1 /= sumWz;
-                    (1 - w[9][k])*w[0][i]*(1 - w[5][j])
+                wz2 /= sumWz;
-                  + w[9][k]*w[3][i]*(1 - w[6][j])
+                wz3 /= sumWz;
-                );
+                wz4 /= sumWz;
-
-                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;
 
 
                 // Points on straight edges
-                const vector edgex1 = p000 + (p100 - p000)*w[0][i];
+                const vector edgex1 = p000 + (p100 - p000)*w0;
-                const vector edgex2 = p010 + (p110 - p010)*w[1][i];
+                const vector edgex2 = p010 + (p110 - p010)*w1;
-                const vector edgex3 = p011 + (p111 - p011)*w[2][i];
+                const vector edgex3 = p011 + (p111 - p011)*w2;
-                const vector edgex4 = p001 + (p101 - p001)*w[3][i];
+                const vector edgex4 = p001 + (p101 - p001)*w3;
 
-                const vector edgey1 = p000 + (p010 - p000)*w[4][j];
+                const vector edgey1 = p000 + (p010 - p000)*w4;
-                const vector edgey2 = p100 + (p110 - p100)*w[5][j];
+                const vector edgey2 = p100 + (p110 - p100)*w5;
-                const vector edgey3 = p101 + (p111 - p101)*w[6][j];
+                const vector edgey3 = p101 + (p111 - p101)*w6;
-                const vector edgey4 = p001 + (p011 - p001)*w[7][j];
+                const vector edgey4 = p001 + (p011 - p001)*w7;
 
-                const vector edgez1 = p000 + (p001 - p000)*w[8][k];
+                const vector edgez1 = p000 + (p001 - p000)*w8;
-                const vector edgez2 = p100 + (p101 - p100)*w[9][k];
+                const vector edgez2 = p100 + (p101 - p100)*w9;
-                const vector edgez3 = p110 + (p111 - p110)*w[10][k];
+                const vector edgez3 = p110 + (p111 - p110)*w10;
-                const vector edgez4 = p010 + (p011 - p010)*w[11][k];
+                const vector edgez4 = p010 + (p011 - p010)*w11;
 
                 // Add the contributions
-                points_[vertexNo] =
+                points_[vijk] =
                 (
-                    impx1*edgex1 + impx2*edgex2 + impx3*edgex3 + impx4*edgex4
+                    wx1*edgex1 + wx2*edgex2 + wx3*edgex3 + wx4*edgex4
-                  + impy1*edgey1 + impy2*edgey2 + impy3*edgey3 + impy4*edgey4
+                  + wy1*edgey1 + wy2*edgey2 + wy3*edgey3 + wy4*edgey4
-                  + impz1*edgez1 + impz2*edgez2 + impz3*edgez3 + impz4*edgez4
+                  + wz1*edgez1 + wz2*edgez2 + wz3*edgez3 + wz4*edgez4
-                )/3.0;
+                )/3;
 
 
-                // Calculate the correction vectors
+                // Apply curved-edge correction if block has curved edges
-                const vector corx1 = impx1*(p[0][i] - edgex1);
+                if (nCurvedEdges)
-                const vector corx2 = impx2*(p[1][i] - edgex2);
+                {
-                const vector corx3 = impx3*(p[2][i] - edgex3);
+                    // Calculate the correction vectors
-                const vector corx4 = impx4*(p[3][i] - edgex4);
+                    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 cory1 = wy1*(p[4][j] - edgey1);
-                const vector cory2 = impy2*(p[5][j] - edgey2);
+                    const vector cory2 = wy2*(p[5][j] - edgey2);
-                const vector cory3 = impy3*(p[6][j] - edgey3);
+                    const vector cory3 = wy3*(p[6][j] - edgey3);
-                const vector cory4 = impy4*(p[7][j] - edgey4);
+                    const vector cory4 = wy4*(p[7][j] - edgey4);
 
-                const vector corz1 = impz1*(p[8][k] - edgez1);
+                    const vector corz1 = wz1*(p[8][k] - edgez1);
-                const vector corz2 = impz2*(p[9][k] - edgez2);
+                    const vector corz2 = wz2*(p[9][k] - edgez2);
-                const vector corz3 = impz3*(p[10][k] - edgez3);
+                    const vector corz3 = wz3*(p[10][k] - edgez3);
-                const vector corz4 = impz4*(p[11][k] - edgez4);
+                    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
+                    (1 - w0)*(1 - w4)*(1 - w8)
-                  + cory1 + cory2 + cory3 + cory4
+                  + (1 - w1)*w4*(1 - w11)
-                  + corz1 + corz2 + corz3 + corz4
+                  + (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();
 
-    //
+    List<FixedList<label, 8>> cells(nCells());
-    // Generate cells
-    //
-    cells_.clear();
-    cells_.setSize(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);
+                celli++;
-                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++;
             }
         }
     }
+
+    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();
 
-    //
+    label patchi = 0;
-    // Generate boundaries on each side of the hex
+    label facei = 0;
-    //
-    boundaryPatches_.clear();
-    boundaryPatches_.setSize(6);
-
 
     // x-direction
 
-    label wallLabel = 0;
-    label wallCellLabel = 0;
-
     // x-min
-    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(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
+            facei++;
-            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++;
         }
     }
 
     // x-max
-    wallLabel++;
+    patchi++;
-    wallCellLabel = 0;
+    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
+            facei++;
-            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++;
         }
     }
 
     // y-direction
 
     // y-min
-    wallLabel++;
+    patchi++;
-    wallCellLabel = 0;
+    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,   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
+            facei++;
-            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++;
         }
     }
 
     // y-max
-    wallLabel++;
+    patchi++;
-    wallCellLabel = 0;
+    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
+            facei++;
-            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++;
         }
     }
 
     // z-direction
 
     // z-min
-    wallLabel++;
+    patchi++;
-    wallCellLabel = 0;
+    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
+            facei++;
-            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++;
         }
     }
 
     // z-max
-    wallLabel++;
+    patchi++;
-    wallCellLabel = 0;
+    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
+            facei++;
-            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++;
         }
     }
 }
 
 
-void Foam::block::clearGeom()
-{
-    points_.clear();
-    cells_.clear();
-    boundaryPatches_.clear();
-}
-
-
 // ************************************************************************* //

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
+    return points_;
-    (
-        i
-      + j * (density().x() + 1)
-      + k * (density().x() + 1) * (density().y() + 1)
-    );
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+inline const Foam::FixedList<Foam::List<Foam::FixedList<Foam::label, 4>>, 6>&
-
+Foam::block::boundaryPatches() const
-inline const Foam::blockDescriptor& Foam::block::blockDef() const
 {
-    return *this;
+    return boundaryPatches_;
 }
 
 

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_
+    expand_(12, gradingDescriptors()),
-    (
+    zoneName_(),
-        12,
+    curvedFaces_(-1),
-        gradingDescriptors()
+    nCurvedFaces_(0)
-    ),
-    zoneName_()
 {
     // 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 * * * * * * * * * * * * * * //

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
-        f2     3 ---- 2
+        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
+        //- Calculate the points and weights for the specified edge.
-        void edgePointsWeights
+        //  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
+        //- Vertex label offset for a particular i,j,k position
-        void edgesPointsWeights
+        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
 

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;
+                label pn = p.size() - 1;
-                forAll(p, pI)
+                forAll(p, pi)
                 {
-                    edgePoints[edgei][pI]  = p[pMax - pI];
+                    edgePoints[edgei][pi]  = p[pn - pi];
-                    edgeWeights[edgei][pI] = 1.0 - d[pMax - 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);
+    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 0,  0, 1, ni);
-    edgePointsWeights(edgePoints, edgeWeights, 1,  3, 2, ni);
+    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 1,  3, 2, ni);
-    edgePointsWeights(edgePoints, edgeWeights, 2,  7, 6, ni);
+    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 2,  7, 6, ni);
-    edgePointsWeights(edgePoints, edgeWeights, 3,  4, 5, ni);
+    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 3,  4, 5, ni);
 
     // Y-direction
     const label nj = density_.y();
-    edgePointsWeights(edgePoints, edgeWeights, 4,  0, 3, nj);
+    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 4,  0, 3, nj);
-    edgePointsWeights(edgePoints, edgeWeights, 5,  1, 2, nj);
+    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 5,  1, 2, nj);
-    edgePointsWeights(edgePoints, edgeWeights, 6,  5, 6, nj);
+    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 6,  5, 6, nj);
-    edgePointsWeights(edgePoints, edgeWeights, 7,  4, 7, nj);
+    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 7,  4, 7, nj);
 
     // Z-direction
     const label nk = density_.z();
-    edgePointsWeights(edgePoints, edgeWeights, 8,  0, 4, nk);
+    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 8,  0, 4, nk);
-    edgePointsWeights(edgePoints, edgeWeights, 9,  1, 5, nk);
+    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 9,  1, 5, nk);
-    edgePointsWeights(edgePoints, edgeWeights, 10, 2, 6, nk);
+    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 10, 2, 6, nk);
-    edgePointsWeights(edgePoints, edgeWeights, 11, 3, 7, nk);
+    nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 11, 3, 7, nk);
+
+    return nCurvedEdges;
 }
 
 

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);
+}
+
+
 // ************************************************************************* //

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;
+}

src/mesh/blockMesh/blockEdges/BSplineEdge/BSpline.C

@@ -23,7 +23,6 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "error.H"
 #include "BSpline.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

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)))

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

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_]),

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

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;

src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.H

@@ -36,11 +36,7 @@ SourceFiles
 #ifndef blockEdges_H
 #define blockEdges_H
 
-#include "edge.H"
+#include "searchableSurfaces.H"
-#include "pointField.H"
-#include "typeInfo.H"
-#include "HashTable.H"
-#include "autoPtr.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&);
+        blockEdge
-
+        (
-        //- Copy construct
+            const pointField&,
-        blockEdge(const blockEdge&);
+            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&);

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)
 {}

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
+    //- Destructor
-        virtual ~lineEdge();
+    virtual ~lineEdge();
 
 
     // Member Functions

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)))

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

src/mesh/blockMesh/blockEdges/splineEdge/CatmullRomSpline.C

@@ -23,7 +23,6 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "error.H"
 #include "CatmullRomSpline.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

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)))

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

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;
+}
+
+
+// ************************************************************************* //

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
+
+// ************************************************************************* //

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);
+}
+
+
+// ************************************************************************* //

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
+
+// ************************************************************************* //

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();
+        }
+    }
+}
+
+
+// ************************************************************************* //

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
+
+// ************************************************************************* //

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++;
         }

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

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 v0 = blocks[blocki].pointLabel(0, 0, 0);
-            label vi1 = blocks[blockI].vtxLabel(1, 0, 0);
+            label vi1 = blocks[blocki].pointLabel(1, 0, 0);
-            scalar diStart = mag(blockPoints[vi1]-blockPoints[v0]);
+            scalar diStart = mag(blockPoints[vi1] - blockPoints[v0]);
 
-            label vinM1 = blocks[blockI].vtxLabel(density.x()-1, 0, 0);
+            label vinM1 = blocks[blocki].pointLabel(density.x()-1, 0, 0);
-            label vin = blocks[blockI].vtxLabel(density.x(), 0, 0);
+            label vin = blocks[blocki].pointLabel(density.x(), 0, 0);
-            scalar diFinal = mag(blockPoints[vin]-blockPoints[vinM1]);
+            scalar diFinal = mag(blockPoints[vin] - blockPoints[vinM1]);
 
-            label vj1 = blocks[blockI].vtxLabel(0, 1, 0);
+            label vj1 = blocks[blocki].pointLabel(0, 1, 0);
-            scalar djStart = mag(blockPoints[vj1]-blockPoints[v0]);
+            scalar djStart = mag(blockPoints[vj1] - blockPoints[v0]);
-            label vjnM1 = blocks[blockI].vtxLabel(0, density.y()-1, 0);
+            label vjnM1 = blocks[blocki].pointLabel(0, density.y()-1, 0);
-            label vjn = blocks[blockI].vtxLabel(0, density.y(), 0);
+            label vjn = blocks[blocki].pointLabel(0, density.y(), 0);
-            scalar djFinal = mag(blockPoints[vjn]-blockPoints[vjnM1]);
+            scalar djFinal = mag(blockPoints[vjn] - blockPoints[vjnM1]);
 
-            label vk1 = blocks[blockI].vtxLabel(0, 0, 1);
+            label vk1 = blocks[blocki].pointLabel(0, 0, 1);
-            scalar dkStart = mag(blockPoints[vk1]-blockPoints[v0]);
+            scalar dkStart = mag(blockPoints[vk1] - blockPoints[v0]);
-            label vknM1 = blocks[blockI].vtxLabel(0, 0, density.z()-1);
+            label vknM1 = blocks[blocki].pointLabel(0, 0, density.z()-1);
-            label vkn = blocks[blockI].vtxLabel(0, 0, density.z());
+            label vkn = blocks[blocki].pointLabel(0, 0, density.z());
-            scalar dkFinal = mag(blockPoints[vkn]-blockPoints[vknM1]);
+            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 =
+                const List<FixedList<label, 4>>& blockPatchFaces =
-                    blocks[blockI].boundaryPatches()[blockFaceLabel];
+                    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();
-    }
-}
-
 // ************************************************************************* //

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();
+        nPoints_ += blocks[blocki].nPoints();
-        nCells_  += blocks[blockI].nCells();
+        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)

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();
+        nPoints_ += blocks[blocki].nPoints();
-        nCells_  += blocks[blockI].nCells();
+        nCells_  += blocks[blocki].nCells();
     }
 
     if (verboseOutput)

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

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;
             }
-
 };
 
 

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_;

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
+
+#------------------------------------------------------------------------------

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)
+        );
+    }
+);
+
+// ************************************************************************* //

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;
+
+
+// ************************************************************************* //

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
+{}
+
+
+// ************************************************************************* //

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;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// ************************************************************************* //