From 24b2335f190bb94a687613ea2e8a21798f039223 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Mon, 31 Oct 2016 18:00:15 +0000
Subject: [PATCH] blockMesh: Added edge projection

New functionality contributed by Mattijs Janssens:
  - new edge projection: projectCurve for use with new geometry
    'searchableCurve'
  - new tutorial 'pipe'
  - naming of vertices and blocks (see pipe tutorial). Including back
    substitution for error messages.
---
 .../meshes/meshShapes/cellShape/cellShape.H   |  10 +-
 .../meshes/meshShapes/cellShape/cellShapeI.H  |  18 +
 src/mesh/blockMesh/Make/files                 |   9 +-
 src/mesh/blockMesh/Make/options               |   2 +
 .../blockDescriptor/blockDescriptor.C         |  99 +++-
 .../blockDescriptor/blockDescriptor.H         |  26 ++
 .../blockDescriptor/blockDescriptorI.H        |   6 +
 .../blockDescriptorTemplates.C                | 114 +++++
 .../blockEdges/BSplineEdge/BSplineEdge.C      |   4 +-
 .../blockEdges/BSplineEdge/BSplineEdge.H      |   2 +
 .../blockMesh/blockEdges/arcEdge/arcEdge.C    |   4 +-
 .../blockMesh/blockEdges/arcEdge/arcEdge.H    |   2 +
 .../blockEdges/blockEdge/blockEdge.C          |  30 +-
 .../blockEdges/blockEdge/blockEdge.H          |  26 +-
 .../blockMesh/blockEdges/lineEdge/lineEdge.C  |   4 +-
 .../blockMesh/blockEdges/lineEdge/lineEdge.H  |   2 +
 .../blockEdges/polyLineEdge/polyLineEdge.C    |   4 +-
 .../blockEdges/polyLineEdge/polyLineEdge.H    |   2 +
 .../projectCurveEdge/projectCurveEdge.C       | 255 +++++++++++
 .../projectCurveEdge/projectCurveEdge.H       | 128 ++++++
 .../blockEdges/projectEdge/projectEdge.C      | 131 ++++--
 .../blockEdges/projectEdge/projectEdge.H      |   2 +
 .../blockEdges/splineEdge/splineEdge.C        |   4 +-
 .../blockEdges/splineEdge/splineEdge.H        |   2 +
 .../blockFaces/blockFace/blockFace.C          |  50 +-
 .../blockFaces/blockFace/blockFace.H          |  26 +-
 .../blockFaces/projectFace/projectFace.C      | 278 ++++++++++-
 .../blockFaces/projectFace/projectFace.H      |  18 +
 src/mesh/blockMesh/blockMesh/blockMesh.C      |   2 +-
 src/mesh/blockMesh/blockMesh/blockMesh.H      |   2 +-
 src/mesh/blockMesh/blockMesh/blockMeshCheck.C |  26 +-
 .../blockMesh/blockMesh/blockMeshTopology.C   |  31 +-
 .../blockVertices/blockVertex/blockVertex.C   |   8 +-
 .../blockVertices/blockVertex/blockVertex.H   |  16 +-
 .../blockVertices/namedVertex/namedVertex.C   |  81 ++++
 .../blockVertices/namedVertex/namedVertex.H   | 104 +++++
 .../blockVertices/pointVertex/pointVertex.C   |   2 +
 .../blockVertices/pointVertex/pointVertex.H   |   2 +
 .../projectVertex/projectVertex.C             |   9 +-
 .../projectVertex/projectVertex.H             |   2 +
 src/mesh/blockMesh/{ => blocks}/block/block.C |  55 ++-
 src/mesh/blockMesh/{ => blocks}/block/block.H |  50 +-
 .../{ => blocks}/block/blockCreate.C          |   0
 .../blockMesh/{ => blocks}/block/blockI.H     |   0
 .../blockMesh/{ => blocks}/block/blockList.H  |   3 +-
 .../blockMesh/blocks/namedBlock/namedBlock.C  |  75 +++
 .../blockMesh/blocks/namedBlock/namedBlock.H  |  91 ++++
 .../searchableCurve/searchableCurve.C         | 432 ++++++++++++++++++
 .../searchableCurve/searchableCurve.H         | 245 ++++++++++
 tutorials/mesh/blockMesh/pipe/Allrun          |   9 +
 .../pipe/constant/triSurface/curve.obj        |   9 +
 .../pipe/constant/triSurface/curve2.vtk       |  13 +
 .../mesh/blockMesh/pipe/system/blockMeshDict  | 177 +++++++
 .../mesh/blockMesh/pipe/system/controlDict    |  58 +++
 .../mesh/blockMesh/pipe/system/fvSchemes      |  37 ++
 .../mesh/blockMesh/pipe/system/fvSolution     |  18 +
 56 files changed, 2701 insertions(+), 114 deletions(-)
 create mode 100644 src/mesh/blockMesh/blockDescriptor/blockDescriptorTemplates.C
 create mode 100644 src/mesh/blockMesh/blockEdges/projectCurveEdge/projectCurveEdge.C
 create mode 100644 src/mesh/blockMesh/blockEdges/projectCurveEdge/projectCurveEdge.H
 create mode 100644 src/mesh/blockMesh/blockVertices/namedVertex/namedVertex.C
 create mode 100644 src/mesh/blockMesh/blockVertices/namedVertex/namedVertex.H
 rename src/mesh/blockMesh/{ => blocks}/block/block.C (61%)
 rename src/mesh/blockMesh/{ => blocks}/block/block.H (78%)
 rename src/mesh/blockMesh/{ => blocks}/block/blockCreate.C (100%)
 rename src/mesh/blockMesh/{ => blocks}/block/blockI.H (100%)
 rename src/mesh/blockMesh/{ => blocks}/block/blockList.H (95%)
 create mode 100644 src/mesh/blockMesh/blocks/namedBlock/namedBlock.C
 create mode 100644 src/mesh/blockMesh/blocks/namedBlock/namedBlock.H
 create mode 100644 src/mesh/blockMesh/searchableCurve/searchableCurve.C
 create mode 100644 src/mesh/blockMesh/searchableCurve/searchableCurve.H
 create mode 100755 tutorials/mesh/blockMesh/pipe/Allrun
 create mode 100644 tutorials/mesh/blockMesh/pipe/constant/triSurface/curve.obj
 create mode 100644 tutorials/mesh/blockMesh/pipe/constant/triSurface/curve2.vtk
 create mode 100644 tutorials/mesh/blockMesh/pipe/system/blockMeshDict
 create mode 100644 tutorials/mesh/blockMesh/pipe/system/controlDict
 create mode 100644 tutorials/mesh/blockMesh/pipe/system/fvSchemes
 create mode 100644 tutorials/mesh/blockMesh/pipe/system/fvSolution

diff --git a/src/OpenFOAM/meshes/meshShapes/cellShape/cellShape.H b/src/OpenFOAM/meshes/meshShapes/cellShape/cellShape.H
index 371d95a063..653da9c46d 100644
--- a/src/OpenFOAM/meshes/meshShapes/cellShape/cellShape.H
+++ b/src/OpenFOAM/meshes/meshShapes/cellShape/cellShape.H
@@ -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
@@ -92,6 +92,14 @@ public:
             const bool doCollapse = false
         );
 
+        //- Construct from components
+        inline cellShape
+        (
+            const word& model,
+            const labelList&,
+            const bool doCollapse = false
+        );
+
         //- Construct from Istream
         inline cellShape(Istream& is);
 
diff --git a/src/OpenFOAM/meshes/meshShapes/cellShape/cellShapeI.H b/src/OpenFOAM/meshes/meshShapes/cellShape/cellShapeI.H
index 5120cd4dfa..42055e04c9 100644
--- a/src/OpenFOAM/meshes/meshShapes/cellShape/cellShapeI.H
+++ b/src/OpenFOAM/meshes/meshShapes/cellShape/cellShapeI.H
@@ -25,6 +25,7 @@ License
 
 #include "Istream.H"
 #include "cell.H"
+#include "cellModeller.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -51,6 +52,23 @@ inline Foam::cellShape::cellShape
 }
 
 
+inline Foam::cellShape::cellShape
+(
+    const word& model,
+    const labelList& l,
+    const bool doCollapse
+)
+:
+    labelList(l),
+    m(cellModeller::lookup(model))
+{
+    if (doCollapse)
+    {
+        collapse();
+    }
+}
+
+
 inline Foam::cellShape::cellShape(Istream& is)
 {
     is >> *this;
diff --git a/src/mesh/blockMesh/Make/files b/src/mesh/blockMesh/Make/files
index 4d1da72068..f552dec8f9 100644
--- a/src/mesh/blockMesh/Make/files
+++ b/src/mesh/blockMesh/Make/files
@@ -1,6 +1,7 @@
 blockVertices/blockVertex/blockVertex.C
 blockVertices/pointVertex/pointVertex.C
 blockVertices/projectVertex/projectVertex.C
+blockVertices/namedVertex/namedVertex.C
 
 blockEdges/blockEdge/blockEdge.C
 blockEdges/lineDivide/lineDivide.C
@@ -13,6 +14,7 @@ blockEdges/BSplineEdge/BSplineEdge.C
 blockEdges/splineEdge/CatmullRomSpline.C
 blockEdges/splineEdge/splineEdge.C
 blockEdges/projectEdge/projectEdge.C
+blockEdges/projectCurveEdge/projectCurveEdge.C
 
 blockFaces/blockFace/blockFace.C
 blockFaces/projectFace/projectFace.C
@@ -23,8 +25,9 @@ gradingDescriptor/gradingDescriptors.C
 blockDescriptor/blockDescriptor.C
 blockDescriptor/blockDescriptorEdges.C
 
-block/block.C
-block/blockCreate.C
+blocks/block/block.C
+blocks/block/blockCreate.C
+blocks/namedBlock/namedBlock.C
 
 blockMesh/blockMesh.C
 blockMesh/blockMeshCreate.C
@@ -33,4 +36,6 @@ blockMesh/blockMeshCheck.C
 blockMesh/blockMeshMerge.C
 blockMesh/blockMeshMergeFast.C
 
+searchableCurve/searchableCurve.C
+
 LIB = $(FOAM_LIBBIN)/libblockMesh
diff --git a/src/mesh/blockMesh/Make/options b/src/mesh/blockMesh/Make/options
index bf6ca71abf..5467dc3ce3 100644
--- a/src/mesh/blockMesh/Make/options
+++ b/src/mesh/blockMesh/Make/options
@@ -1,9 +1,11 @@
 EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/fileFormats/lnInclude \
+    -I$(LIB_SRC)/edgeMesh/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude
 
 LIB_LIBS = \
     -lmeshTools \
     -lfileFormats \
+    -ledgeMesh \
     -lsurfMesh
diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C
index fb9a413013..ddcb9466e0 100644
--- a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C
+++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C
@@ -128,6 +128,90 @@ void Foam::blockDescriptor::findCurvedFaces()
 }
 
 
+void Foam::blockDescriptor::read
+(
+    Istream& is,
+    label& val,
+    const dictionary& dict
+)
+{
+    token t(is);
+    if (t.isLabel())
+    {
+        val = t.labelToken();
+    }
+    else if (t.isWord())
+    {
+        const word& varName = t.wordToken();
+        const entry* ePtr = dict.lookupScopedEntryPtr
+        (
+            varName,
+            true,
+            true
+        );
+        if (ePtr)
+        {
+            // Read as label
+            val = Foam::readLabel(ePtr->stream());
+        }
+        else
+        {
+            FatalIOErrorInFunction(is)
+                << "Undefined variable "
+                << varName << ". Valid variables are " << dict
+                << exit(FatalIOError);
+        }
+    }
+    else
+    {
+        FatalIOErrorInFunction(is)
+            << "Illegal token " << t.info()
+            << " when trying to read label"
+            << exit(FatalIOError);
+    }
+
+    is.fatalCheck
+    (
+        "operator>>(Istream&, List<T>&) : reading entry"
+    );
+}
+
+
+Foam::label Foam::blockDescriptor::read
+(
+    Istream& is,
+    const dictionary& dict
+)
+{
+    label val;
+    read(is, val, dict);
+    return val;
+}
+
+
+void Foam::blockDescriptor::write
+(
+    Ostream& os,
+    const label val,
+    const dictionary& dict
+)
+{
+    forAllConstIter(dictionary, dict, iter)
+    {
+        if (iter().isStream())
+        {
+            label keyVal(Foam::readLabel(iter().stream()));
+            if (keyVal == val)
+            {
+                os << iter().keyword();
+                return;
+            }
+        }
+    }
+    os << val;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::blockDescriptor::blockDescriptor
@@ -164,6 +248,8 @@ Foam::blockDescriptor::blockDescriptor
 
 Foam::blockDescriptor::blockDescriptor
 (
+    const dictionary& dict,
+    const label index,
     const pointField& vertices,
     const blockEdgeList& edges,
     const blockFaceList& faces,
@@ -173,13 +259,24 @@ Foam::blockDescriptor::blockDescriptor
     vertices_(vertices),
     edges_(edges),
     faces_(faces),
-    blockShape_(is),
     density_(),
     expand_(12, gradingDescriptors()),
     zoneName_(),
     curvedFaces_(-1),
     nCurvedFaces_(0)
 {
+    // Read cell model and list of vertices (potentially with variables)
+    word model(is);
+    blockShape_ = cellShape
+    (
+        model,
+        read<label>
+        (
+            is,
+            dict.subOrEmptyDict("namedVertices")
+        )
+    );
+
     // Examine next token
     token t(is);
 
diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H
index b0f7b1412c..45c5e75561 100644
--- a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H
+++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H
@@ -152,6 +152,8 @@ public:
         //- Construct from Istream
         blockDescriptor
         (
+            const dictionary& dict,
+            const label index,
             const pointField& vertices,
             const blockEdgeList&,
             const blockFaceList&,
@@ -161,6 +163,9 @@ public:
 
     // Member Functions
 
+        //- Reference to point field defining the block mesh
+        inline const pointField& vertices() const;
+
         //- Return reference to the list of curved faces
         inline const blockFaceList& faces() const;
 
@@ -233,6 +238,23 @@ public:
         //  to lie on the faces of the block
         void correctFacePoints(FixedList<pointField, 6>&) const;
 
+        //- In-place read with dictionary lookup
+        static void read(Istream&, label&, const dictionary&);
+
+        //- In-place read with dictionary lookup
+        template<class T>
+        static void read(Istream&, List<T>&, const dictionary&);
+
+        //- Return-read with dictionary lookup
+        static label read(Istream&, const dictionary&);
+
+        //- Return-read with dictionary lookup
+        template<class T>
+        static List<T> read(Istream& is, const dictionary&);
+
+        //- Write with dictionary lookup
+        static void write(Ostream&, const label, const dictionary&);
+
 
     // IOstream Operators
 
@@ -248,6 +270,10 @@ public:
 
 #include "blockDescriptorI.H"
 
+#ifdef NoRepository
+    #include "blockDescriptorTemplates.C"
+#endif
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif
diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptorI.H b/src/mesh/blockMesh/blockDescriptor/blockDescriptorI.H
index 823860a820..13ec39c2e6 100644
--- a/src/mesh/blockMesh/blockDescriptor/blockDescriptorI.H
+++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptorI.H
@@ -25,6 +25,12 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+inline const Foam::pointField& Foam::blockDescriptor::vertices() const
+{
+    return vertices_;
+}
+
+
 inline const Foam::blockFaceList& Foam::blockDescriptor::faces() const
 {
     return faces_;
diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptorTemplates.C b/src/mesh/blockMesh/blockDescriptor/blockDescriptorTemplates.C
new file mode 100644
index 0000000000..aacd20d8bd
--- /dev/null
+++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptorTemplates.C
@@ -0,0 +1,114 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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  * * * * * * * * * * * * * //
+
+template<class T>
+void Foam::blockDescriptor::read
+(
+    Istream& is,
+    List<T>& L,
+    const dictionary& dict
+)
+{
+    token firstToken(is);
+
+    if (firstToken.isLabel())
+    {
+        label s = firstToken.labelToken();
+
+        // Set list length to that read
+        L.setSize(s);
+
+        // Read list contents depending on data format
+
+        // Read beginning of contents
+        char delimiter = is.readBeginList("List");
+
+        if (s)
+        {
+            if (delimiter == token::BEGIN_LIST)
+            {
+                for (label i=0; i<s; i++)
+                {
+                    read(is, L[i], dict);
+                }
+            }
+        }
+
+        // Read end of contents
+        is.readEndList("List");
+    }
+    else if (firstToken.isPunctuation())
+    {
+        if (firstToken.pToken() != token::BEGIN_LIST)
+        {
+            FatalIOErrorInFunction(is)
+                << "incorrect first token, expected '(', found "
+                << firstToken.info()
+                << exit(FatalIOError);
+        }
+
+        SLList<T> sll;
+
+        while (true)
+        {
+            token t(is);
+            if (t.isPunctuation() && t.pToken() == token::END_LIST)
+            {
+                break;
+            }
+            is.putBack(t);
+            T elem;
+            read(is, elem, dict);
+            sll.append(elem);
+        }
+
+        // Convert the singly-linked list to this list
+        L = sll;
+    }
+    else
+    {
+        FatalIOErrorInFunction(is)
+            << "incorrect first token, expected <int> or '(', found "
+            << firstToken.info()
+            << exit(FatalIOError);
+    }
+}
+
+
+template<class T>
+Foam::List<T> Foam::blockDescriptor::read
+(
+    Istream& is,
+    const dictionary& dict
+)
+{
+    List<T> L;
+    read(is, L, dict);
+    return L;
+}
+
+// ************************************************************************* //
diff --git a/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.C b/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.C
index 4e16c18264..dc62dbc3e4 100644
--- a/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.C
+++ b/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.C
@@ -62,12 +62,14 @@ Foam::blockEdges::BSplineEdge::BSplineEdge
 
 Foam::blockEdges::BSplineEdge::BSplineEdge
 (
+    const dictionary& dict,
+    const label index,
     const searchableSurfaces& geometry,
     const pointField& points,
     Istream& is
 )
 :
-    blockEdge(points, is),
+    blockEdge(dict, index, points, is),
     BSpline(appendEndPoints(points, start_, end_, pointField(is)))
 {
     token t(is);
diff --git a/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.H b/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.H
index 81c06d63a7..72cbeeab7a 100644
--- a/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.H
+++ b/src/mesh/blockMesh/blockEdges/BSplineEdge/BSplineEdge.H
@@ -83,6 +83,8 @@ public:
         //- Construct from Istream, setting pointsList
         BSplineEdge
         (
+            const dictionary& dict,
+            const label index,
             const searchableSurfaces& geometry,
             const pointField&,
             Istream&
diff --git a/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.C b/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.C
index 5d012e0531..793cc27b2f 100644
--- a/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.C
+++ b/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.C
@@ -123,12 +123,14 @@ Foam::blockEdges::arcEdge::arcEdge
 
 Foam::blockEdges::arcEdge::arcEdge
 (
+    const dictionary& dict,
+    const label index,
     const searchableSurfaces& geometry,
     const pointField& points,
     Istream& is
 )
 :
-    blockEdge(points, is),
+    blockEdge(dict, index, points, is),
     p1_(points_[start_]),
     p2_(is),
     p3_(points_[end_]),
diff --git a/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.H b/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.H
index b7a71ff6d1..9d41e8cceb 100644
--- a/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.H
+++ b/src/mesh/blockMesh/blockEdges/arcEdge/arcEdge.H
@@ -91,6 +91,8 @@ public:
         //- Construct from Istream setting pointsList
         arcEdge
         (
+            const dictionary& dict,
+            const label index,
             const searchableSurfaces& geometry,
             const pointField& points,
             Istream&
diff --git a/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.C b/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.C
index 4ed2647a74..0ef5fd0408 100644
--- a/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.C
+++ b/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.C
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "blockEdge.H"
+#include "blockDescriptor.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -51,13 +52,15 @@ Foam::blockEdge::blockEdge
 
 Foam::blockEdge::blockEdge
 (
+    const dictionary& dict,
+    const label index,
     const pointField& points,
     Istream& is
 )
 :
     points_(points),
-    start_(readLabel(is)),
-    end_(readLabel(is))
+    start_(blockDescriptor::read(is, dict.subOrEmptyDict("namedVertices"))),
+    end_(blockDescriptor::read(is, dict.subOrEmptyDict("namedVertices")))
 {}
 
 
@@ -70,6 +73,8 @@ Foam::autoPtr<Foam::blockEdge> Foam::blockEdge::clone() const
 
 Foam::autoPtr<Foam::blockEdge> Foam::blockEdge::New
 (
+    const dictionary& dict,
+    const label index,
     const searchableSurfaces& geometry,
     const pointField& points,
     Istream& is
@@ -95,7 +100,7 @@ Foam::autoPtr<Foam::blockEdge> Foam::blockEdge::New
             << abort(FatalError);
     }
 
-    return autoPtr<blockEdge>(cstrIter()(geometry, points, is));
+    return autoPtr<blockEdge>(cstrIter()(dict, index, geometry, points, is));
 }
 
 
@@ -139,6 +144,25 @@ Foam::blockEdge::position(const scalarList& lambdas) const
 }
 
 
+void Foam::blockEdge::write(Ostream& os, const dictionary& d) const
+{
+    const dictionary* varDictPtr = d.subDictPtr("namedVertices");
+    if (varDictPtr)
+    {
+        const dictionary& varDict = *varDictPtr;
+
+        blockDescriptor::write(os, start_, varDict);
+        os << tab;
+        blockDescriptor::write(os, end_, varDict);
+        os << endl;
+    }
+    else
+    {
+        os << start_ << tab << end_ << endl;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 Foam::Ostream& Foam::operator<<(Ostream& os, const blockEdge& p)
diff --git a/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.H b/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.H
index e6e3196ea9..bde77adbc2 100644
--- a/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.H
+++ b/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.H
@@ -92,11 +92,13 @@ public:
             blockEdge,
             Istream,
             (
+                const dictionary& dict,
+                const label index,
                 const searchableSurfaces& geometry,
                 const pointField& points,
                 Istream& is
             ),
-            (geometry, points, is)
+            (dict, index, geometry, points, is)
         );
 
 
@@ -113,6 +115,8 @@ public:
         //- Construct from Istream setting pointsList
         blockEdge
         (
+            const dictionary& dict,
+            const label index,
             const pointField&,
             Istream&
         );
@@ -123,6 +127,8 @@ public:
         //- New function which constructs and returns pointer to a blockEdge
         static autoPtr<blockEdge> New
         (
+            const dictionary& dict,
+            const label index,
             const searchableSurfaces& geometry,
             const pointField&,
             Istream&
@@ -132,20 +138,29 @@ public:
         //  PtrLists of blockEdge
         class iNew
         {
+            const dictionary& dict_;
             const searchableSurfaces& geometry_;
             const pointField& points_;
+            mutable label index_;
 
         public:
 
-            iNew(const searchableSurfaces& geometry, const pointField& points)
+            iNew
+            (
+                const dictionary& dict,
+                const searchableSurfaces& geometry,
+                const pointField& points
+            )
             :
+                dict_(dict),
                 geometry_(geometry),
-                points_(points)
+                points_(points),
+                index_(0)
             {}
 
             autoPtr<blockEdge> operator()(Istream& is) const
             {
-                return blockEdge::New(geometry_, points_, is);
+                return blockEdge::New(dict_, index_++, geometry_, points_, is);
             }
         };
 
@@ -195,6 +210,9 @@ public:
         //- Return the length of the curve
         virtual scalar length() const = 0;
 
+        //- Write edge with variable backsubstitution
+        void write(Ostream&, const dictionary&) const;
+
 
     // Ostream operator
 
diff --git a/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.C b/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.C
index 45b4434850..9206151f56 100644
--- a/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.C
+++ b/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.C
@@ -53,12 +53,14 @@ Foam::blockEdges::lineEdge::lineEdge
 
 Foam::blockEdges::lineEdge::lineEdge
 (
+    const dictionary& dict,
+    const label index,
     const searchableSurfaces& geometry,
     const pointField& points,
     Istream& is
 )
 :
-    blockEdge(points, is)
+    blockEdge(dict, index, points, is)
 {}
 
 
diff --git a/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.H b/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.H
index e70ffd10bf..c5266b28c3 100644
--- a/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.H
+++ b/src/mesh/blockMesh/blockEdges/lineEdge/lineEdge.H
@@ -68,6 +68,8 @@ public:
         //- Construct from Istream with a pointField
         lineEdge
         (
+            const dictionary& dict,
+            const label index,
             const searchableSurfaces& geometry,
             const pointField&,
             Istream&
diff --git a/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.C b/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.C
index 1396636a66..63e8fb4d00 100644
--- a/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.C
+++ b/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.C
@@ -55,12 +55,14 @@ Foam::blockEdges::polyLineEdge::polyLineEdge
 
 Foam::blockEdges::polyLineEdge::polyLineEdge
 (
+    const dictionary& dict,
+    const label index,
     const searchableSurfaces& geometry,
     const pointField& ps,
     Istream& is
 )
 :
-    blockEdge(ps, is),
+    blockEdge(dict, index, ps, is),
     polyLine(appendEndPoints(ps, start_, end_, pointField(is)))
 {}
 
diff --git a/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.H b/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.H
index adcd7fd1b4..62e42ca31c 100644
--- a/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.H
+++ b/src/mesh/blockMesh/blockEdges/polyLineEdge/polyLineEdge.H
@@ -83,6 +83,8 @@ public:
         //- Construct from Istream
         polyLineEdge
         (
+            const dictionary& dict,
+            const label index,
             const searchableSurfaces& geometry,
             const pointField&,
             Istream&
diff --git a/src/mesh/blockMesh/blockEdges/projectCurveEdge/projectCurveEdge.C b/src/mesh/blockMesh/blockEdges/projectCurveEdge/projectCurveEdge.C
new file mode 100644
index 0000000000..0e0166fd44
--- /dev/null
+++ b/src/mesh/blockMesh/blockEdges/projectCurveEdge/projectCurveEdge.C
@@ -0,0 +1,255 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "searchableSurfacesQueries.H"
+#include "projectCurveEdge.H"
+#include "unitConversion.H"
+#include "addToRunTimeSelectionTable.H"
+#include "pointConstraint.H"
+#include "OBJstream.H"
+#include "linearInterpolationWeights.H"
+#include "searchableCurve.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(projectCurveEdge, 0);
+    addToRunTimeSelectionTable(blockEdge, projectCurveEdge, Istream);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::projectCurveEdge::projectCurveEdge
+(
+    const dictionary& dict,
+    const label index,
+    const searchableSurfaces& geometry,
+    const pointField& points,
+    Istream& is
+)
+:
+    blockEdge(dict, index, points, is),
+    geometry_(geometry)
+{
+    wordList names(is);
+    surfaces_.setSize(names.size());
+    forAll(names, i)
+    {
+        surfaces_[i] = geometry_.findSurfaceID(names[i]);
+
+        if (surfaces_[i] == -1)
+        {
+            FatalIOErrorInFunction(is)
+                << "Cannot find surface " << names[i] << " in geometry"
+                << exit(FatalIOError);
+        }
+
+        if (isA<searchableCurve>(geometry_[surfaces_[i]]))
+        {
+            Info<< type() << " : Using curved surface "
+                << geometry_[surfaces_[i]].name()
+                << " to predict starting points." << endl;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::pointField>
+Foam::projectCurveEdge::position(const scalarList& lambdas) const
+{
+    // For debugging to tag the output
+    static label eIter = 0;
+
+    autoPtr<OBJstream> debugStr;
+    if (debug)
+    {
+        debugStr.reset
+        (
+            new OBJstream("projectCurveEdge_" + Foam::name(eIter++) + ".obj")
+        );
+        Info<< "Writing lines from straight-line start points"
+            << " to projected points to " << debugStr().name() << endl;
+    }
+
+
+    tmp<pointField> tpoints(new pointField(lambdas.size()));
+    pointField& points = tpoints.ref();
+
+    const point& startPt = points_[start_];
+    const point& endPt = points_[end_];
+    const vector d = endPt-startPt;
+
+    // Initial guess
+    forAll(lambdas, i)
+    {
+        points[i] = startPt+lambdas[i]*d;
+    }
+
+    // Use special interpolation to keep initial guess on same position on
+    // surface
+    forAll(surfaces_, i)
+    {
+        if (isA<searchableCurve>(geometry_[surfaces_[i]]))
+        {
+            const searchableCurve& s =
+                refCast<const searchableCurve>(geometry_[surfaces_[i]]);
+            List<pointIndexHit> nearInfo;
+            s.findNearest
+            (
+                points[0],
+                points.last(),
+                scalarField(lambdas),
+                scalarField(points.size(), magSqr(d)),
+                nearInfo
+            );
+            forAll(nearInfo, i)
+            {
+                if (nearInfo[i].hit())
+                {
+                    points[i] = nearInfo[i].hitPoint();
+                }
+            }
+
+            break;
+        }
+    }
+
+
+
+    // Upper limit for number of iterations
+    const label maxIter = 10;
+    // Residual tolerance
+    const scalar relTol = 0.1;
+    const scalar absTol = 1e-4;
+
+    scalar initialResidual = 0.0;
+
+    for (label iter = 0; iter < maxIter; iter++)
+    {
+        // Do projection
+        {
+            List<pointConstraint> constraints(lambdas.size());
+            pointField start(points);
+            searchableSurfacesQueries::findNearest
+            (
+                geometry_,
+                surfaces_,
+                start,
+                scalarField(start.size(), magSqr(d)),
+                points,
+                constraints
+            );
+
+            // Reset start and end point
+            if (lambdas[0] < SMALL)
+            {
+                points[0] = startPt;
+            }
+            if (lambdas.last() > 1.0-SMALL)
+            {
+                points.last() = endPt;
+            }
+
+            if (debugStr.valid())
+            {
+                forAll(points, i)
+                {
+                    debugStr().write(linePointRef(start[i], points[i]));
+                }
+            }
+        }
+
+        // Calculate lambdas (normalised coordinate along edge)
+        scalarField projLambdas(points.size());
+        {
+            projLambdas[0] = 0.0;
+            for (label i = 1; i < points.size(); i++)
+            {
+                projLambdas[i] = projLambdas[i-1] + mag(points[i]-points[i-1]);
+            }
+            projLambdas /= projLambdas.last();
+        }
+        linearInterpolationWeights interpolator(projLambdas);
+
+        // Compare actual distances and move points (along straight line;
+        // not along surface)
+        vectorField residual(points.size(), vector::zero);
+        labelList indices;
+        scalarField weights;
+        for (label i = 1; i < points.size() - 1; i++)
+        {
+            interpolator.valueWeights(lambdas[i], indices, weights);
+
+            point predicted = vector::zero;
+            forAll(indices, indexi)
+            {
+                predicted += weights[indexi]*points[indices[indexi]];
+            }
+            residual[i] = predicted-points[i];
+        }
+
+        scalar scalarResidual = sum(mag(residual));
+
+        if (debug)
+        {
+            Pout<< "Iter:" << iter << " initialResidual:" << initialResidual
+                << " residual:" << scalarResidual << endl;
+        }
+
+        if (scalarResidual < absTol*0.5*lambdas.size())
+        {
+            break;
+        }
+        else if (iter == 0)
+        {
+            initialResidual = scalarResidual;
+        }
+        else if (scalarResidual/initialResidual < relTol)
+        {
+            break;
+        }
+
+
+        if (debugStr.valid())
+        {
+            forAll(points, i)
+            {
+                const linePointRef ln(points[i], points[i]+residual[i]);
+                debugStr().write(ln);
+            }
+        }
+
+        points += residual;
+    }
+
+    return tpoints;
+}
+
+
+// ************************************************************************* //
diff --git a/src/mesh/blockMesh/blockEdges/projectCurveEdge/projectCurveEdge.H b/src/mesh/blockMesh/blockEdges/projectCurveEdge/projectCurveEdge.H
new file mode 100644
index 0000000000..27d7f7ffca
--- /dev/null
+++ b/src/mesh/blockMesh/blockEdges/projectCurveEdge/projectCurveEdge.H
@@ -0,0 +1,128 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::projectCurveEdge
+
+Description
+    Defines the edge from the projection onto a surface (single surface)
+    or intersection of two surfaces.
+
+SourceFiles
+    projectCurveEdge.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef projectCurveEdge_H
+#define projectCurveEdge_H
+
+#include "blockEdge.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class pointConstraint;
+
+/*---------------------------------------------------------------------------*\
+                      Class projectCurveEdge Declaration
+\*---------------------------------------------------------------------------*/
+
+class projectCurveEdge
+:
+    public blockEdge
+{
+    // Private data
+
+        const searchableSurfaces& geometry_;
+
+        //- The indices of surfaces onto which the points are projected
+        labelList surfaces_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        projectCurveEdge(const projectCurveEdge&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const projectCurveEdge&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("projectCurve");
+
+
+    // Constructors
+
+        //- Construct from Istream setting pointsList
+        projectCurveEdge
+        (
+            const dictionary& dict,
+            const label index,
+            const searchableSurfaces& geometry,
+            const pointField& points,
+            Istream&
+        );
+
+
+    //- Destructor
+    virtual ~projectCurveEdge()
+    {}
+
+
+    // Member Functions
+
+        //- Return the point positions corresponding to the curve parameters
+        //  0 <= lambda <= 1
+        virtual point position(const scalar) const
+        {
+            NotImplemented;
+            return point::max;
+        }
+
+        //- Return the point positions corresponding to the curve parameters
+        //  0 <= lambda <= 1
+        virtual tmp<pointField> position(const scalarList&) const;
+
+        //- Return the length of the curve
+        virtual scalar length() const
+        {
+            NotImplemented;
+            return 1;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.C b/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.C
index 9f10f34191..eefdae7af8 100644
--- a/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.C
+++ b/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.C
@@ -29,6 +29,7 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "pointConstraint.H"
 #include "OBJstream.H"
+#include "linearInterpolationWeights.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -79,12 +80,14 @@ void Foam::projectEdge::findNearest
 
 Foam::projectEdge::projectEdge
 (
+    const dictionary& dict,
+    const label index,
     const searchableSurfaces& geometry,
     const pointField& points,
     Istream& is
 )
 :
-    blockEdge(points, is),
+    blockEdge(dict, index, points, is),
     geometry_(geometry)
 {
     wordList names(is);
@@ -125,7 +128,20 @@ Foam::point Foam::projectEdge::position(const scalar lambda) const
 Foam::tmp<Foam::pointField>
 Foam::projectEdge::position(const scalarList& lambdas) const
 {
-    static label iter = 0;
+    // For debugging to tag the output
+    static label eIter = 0;
+
+    autoPtr<OBJstream> debugStr;
+    if (debug)
+    {
+        debugStr.reset
+        (
+            new OBJstream("projectEdge_" + Foam::name(eIter++) + ".obj")
+        );
+        Info<< "Writing lines from straight-line start points"
+            << " to projected points to " << debugStr().name() << endl;
+    }
+
 
     tmp<pointField> tpoints(new pointField(lambdas.size()));
     pointField& points = tpoints.ref();
@@ -141,17 +157,26 @@ Foam::projectEdge::position(const scalarList& lambdas) const
     }
 
 
-    for (label i = 0; i < 3; i++)
+    // Upper limit for number of iterations
+    const label maxIter = 10;
+    // Residual tolerance
+    const scalar relTol = 0.1;
+    const scalar absTol = 1e-4;
+
+    scalar initialResidual = 0.0;
+
+    for (label iter = 0; iter < maxIter; iter++)
     {
         // Do projection
         {
             List<pointConstraint> constraints(lambdas.size());
+            pointField start(points);
             searchableSurfacesQueries::findNearest
             (
                 geometry_,
                 surfaces_,
-                pointField(points),
-                scalarField(points.size(), magSqr(d)),
+                start,
+                scalarField(start.size(), magSqr(d)),
                 points,
                 constraints
             );
@@ -165,63 +190,77 @@ Foam::projectEdge::position(const scalarList& lambdas) const
             {
                 points.last() = endPt;
             }
-        }
 
-        // Calculate distances
-        scalarField nearLength(points.size());
-        {
-            nearLength[0] = 0.0;
-            for(label i = 1; i < points.size(); i++)
+            if (debugStr.valid())
             {
-                nearLength[i] = nearLength[i-1] + mag(points[i]-points[i-1]);
+                forAll(points, i)
+                {
+                    debugStr().write(linePointRef(start[i], points[i]));
+                }
             }
         }
 
+        // Calculate lambdas (normalised coordinate along edge)
+        scalarField projLambdas(points.size());
+        {
+            projLambdas[0] = 0.0;
+            for (label i = 1; i < points.size(); i++)
+            {
+                projLambdas[i] = projLambdas[i-1] + mag(points[i]-points[i-1]);
+            }
+            projLambdas /= projLambdas.last();
+        }
+        linearInterpolationWeights interpolator(projLambdas);
+
         // Compare actual distances and move points (along straight line;
         // not along surface)
-        for(label i = 1; i < points.size() - 1; i++)
+        vectorField residual(points.size(), vector::zero);
+        labelList indices;
+        scalarField weights;
+        for (label i = 1; i < points.size() - 1; i++)
         {
-            scalar nearDelta = mag(points[i]-points[i-1])/nearLength.last();
-            scalar wantedDelta = lambdas[i]-lambdas[i-1];
+            interpolator.valueWeights(lambdas[i], indices, weights);
 
-            vector v(points[i]-points[i-1]);
-            points[i] = points[i-1]+wantedDelta/nearDelta*v;
-        }
-    }
-
-
-    if (debug)
-    {
-        OBJstream str("projectEdge_" + Foam::name(iter++) + ".obj");
-        Info<< "Writing lines from straight-line start points"
-            << " to projected points to " << str.name() << endl;
-
-        pointField startPts(lambdas.size());
-        forAll(lambdas, i)
-        {
-            startPts[i] = startPt+lambdas[i]*d;
+            point predicted = vector::zero;
+            forAll(indices, indexi)
+            {
+                predicted += weights[indexi]*points[indices[indexi]];
+            }
+            residual[i] = predicted-points[i];
         }
 
-        pointField nearPts(lambdas.size());
-        List<pointConstraint> nearConstraints(lambdas.size());
+        scalar scalarResidual = sum(mag(residual));
+
+        if (debug)
         {
-            const scalar distSqr = magSqr(d);
-            searchableSurfacesQueries::findNearest
-            (
-                geometry_,
-                surfaces_,
-                startPts,
-                scalarField(startPts.size(), distSqr),
-                nearPts,
-                nearConstraints
-            );
+            Pout<< "Iter:" << iter << " initialResidual:" << initialResidual
+                << " residual:" << scalarResidual << endl;
         }
 
-        forAll(startPts, i)
+        if (scalarResidual < absTol*0.5*lambdas.size())
         {
-            str.write(linePointRef(startPts[i], nearPts[i]));
-            str.write(linePointRef(nearPts[i], points[i]));
+            break;
         }
+        else if (iter == 0)
+        {
+            initialResidual = scalarResidual;
+        }
+        else if (scalarResidual/initialResidual < relTol)
+        {
+            break;
+        }
+
+
+        if (debugStr.valid())
+        {
+            forAll(points, i)
+            {
+                const linePointRef ln(points[i], points[i]+residual[i]);
+                debugStr().write(ln);
+            }
+        }
+
+        points += residual;
     }
 
     return tpoints;
diff --git a/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.H b/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.H
index cf8baebe75..320f658eb5 100644
--- a/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.H
+++ b/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.H
@@ -84,6 +84,8 @@ public:
         //- Construct from Istream setting pointsList
         projectEdge
         (
+            const dictionary& dict,
+            const label index,
             const searchableSurfaces& geometry,
             const pointField& points,
             Istream&
diff --git a/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.C b/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.C
index 5a6be74f4e..0e70e255c8 100644
--- a/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.C
+++ b/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.C
@@ -62,12 +62,14 @@ Foam::blockEdges::splineEdge::splineEdge
 
 Foam::blockEdges::splineEdge::splineEdge
 (
+    const dictionary& dict,
+    const label index,
     const searchableSurfaces& geometry,
     const pointField& points,
     Istream& is
 )
 :
-    blockEdge(points, is),
+    blockEdge(dict, index, points, is),
     CatmullRomSpline(appendEndPoints(points, start_, end_, pointField(is)))
 {
     token t(is);
diff --git a/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.H b/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.H
index 4080b606c8..e1457a2d87 100644
--- a/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.H
+++ b/src/mesh/blockMesh/blockEdges/splineEdge/splineEdge.H
@@ -83,6 +83,8 @@ public:
         //- Construct from Istream, setting pointsList
         splineEdge
         (
+            const dictionary& dict,
+            const label index,
             const searchableSurfaces& geometry,
             const pointField&,
             Istream&
diff --git a/src/mesh/blockMesh/blockFaces/blockFace/blockFace.C b/src/mesh/blockMesh/blockFaces/blockFace/blockFace.C
index 388c9fcc5d..4ee92176c7 100644
--- a/src/mesh/blockMesh/blockFaces/blockFace/blockFace.C
+++ b/src/mesh/blockMesh/blockFaces/blockFace/blockFace.C
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "blockFace.H"
+#include "blockDescriptor.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -42,9 +43,21 @@ Foam::blockFace::blockFace(const face& vertices)
 {}
 
 
-Foam::blockFace::blockFace(Istream& is)
+Foam::blockFace::blockFace
+(
+    const dictionary& dict,
+    const label index,
+    Istream& is
+)
 :
-    vertices_(is)
+    vertices_
+    (
+        blockDescriptor::read<label>
+        (
+            is,
+            dict.subOrEmptyDict("namedVertices")
+        )
+    )
 {}
 
 
@@ -57,6 +70,8 @@ Foam::autoPtr<Foam::blockFace> Foam::blockFace::clone() const
 
 Foam::autoPtr<Foam::blockFace> Foam::blockFace::New
 (
+    const dictionary& dict,
+    const label index,
     const searchableSurfaces& geometry,
     Istream& is
 )
@@ -81,7 +96,36 @@ Foam::autoPtr<Foam::blockFace> Foam::blockFace::New
             << abort(FatalError);
     }
 
-    return autoPtr<blockFace>(cstrIter()(geometry, is));
+    return autoPtr<blockFace>(cstrIter()(dict, index, geometry, is));
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::blockFace::write(Ostream& os, const dictionary& d) const
+{
+    const dictionary* varDictPtr = d.subDictPtr("namedVertices");
+    if (varDictPtr)
+    {
+        const dictionary& varDict = *varDictPtr;
+
+        // Write size and start delimiter
+        os << vertices_.size() << token::BEGIN_LIST;
+
+        // Write contents
+        forAll(vertices_, i)
+        {
+            if (i > 0) os << token::SPACE;
+            blockDescriptor::write(os, vertices_[i], varDict);
+        }
+
+        // Write end delimiter
+        os << token::END_LIST;
+    }
+    else
+    {
+        os << vertices_ << endl;
+    }
 }
 
 
diff --git a/src/mesh/blockMesh/blockFaces/blockFace/blockFace.H b/src/mesh/blockMesh/blockFaces/blockFace/blockFace.H
index fee6771d79..f890aa3ae7 100644
--- a/src/mesh/blockMesh/blockFaces/blockFace/blockFace.H
+++ b/src/mesh/blockMesh/blockFaces/blockFace/blockFace.H
@@ -77,10 +77,12 @@ public:
             blockFace,
             Istream,
             (
+                const dictionary& dict,
+                const label index,
                 const searchableSurfaces& geometry,
                 Istream& is
             ),
-            (geometry, is)
+            (dict, index, geometry, is)
         );
 
 
@@ -90,7 +92,12 @@ public:
         blockFace(const face& vertices);
 
         //- Construct from Istream
-        blockFace(Istream&);
+        blockFace
+        (
+            const dictionary& dict,
+            const label index,
+            Istream&
+        );
 
         //- Clone function
         virtual autoPtr<blockFace> clone() const;
@@ -98,6 +105,8 @@ public:
         //- New function which constructs and returns pointer to a blockFace
         static autoPtr<blockFace> New
         (
+            const dictionary& dict,
+            const label index,
             const searchableSurfaces& geometry,
             Istream&
         );
@@ -106,18 +115,22 @@ public:
         //  PtrLists of blockFace
         class iNew
         {
+            const dictionary& dict_;
             const searchableSurfaces& geometry_;
+            mutable label index_;
 
         public:
 
-            iNew(const searchableSurfaces& geometry)
+            iNew(const dictionary& dict, const searchableSurfaces& geometry)
             :
-                geometry_(geometry)
+                dict_(dict),
+                geometry_(geometry),
+                index_(0)
             {}
 
             autoPtr<blockFace> operator()(Istream& is) const
             {
-                return blockFace::New(geometry_, is);
+                return blockFace::New(dict_, index_++, geometry_, is);
             }
         };
 
@@ -145,6 +158,9 @@ public:
             pointField& points
         ) const = 0;
 
+        //- Write face with variable backsubstitution
+        void write(Ostream&, const dictionary&) const;
+
 
     // Ostream operator
 
diff --git a/src/mesh/blockMesh/blockFaces/projectFace/projectFace.C b/src/mesh/blockMesh/blockFaces/projectFace/projectFace.C
index 2cf9c97427..13ca420a1c 100644
--- a/src/mesh/blockMesh/blockFaces/projectFace/projectFace.C
+++ b/src/mesh/blockMesh/blockFaces/projectFace/projectFace.C
@@ -27,6 +27,8 @@ License
 #include "unitConversion.H"
 #include "addToRunTimeSelectionTable.H"
 #include "blockDescriptor.H"
+#include "OBJstream.H"
+#include "linearInterpolationWeights.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -66,15 +68,74 @@ const Foam::searchableSurface& Foam::blockFaces::projectFace::lookupSurface
 }
 
 
+Foam::label Foam::blockFaces::projectFace::index
+(
+    const labelPair& n,
+    const labelPair& coord
+) const
+{
+    return coord.first()+coord.second()*n.first();
+}
+
+
+void Foam::blockFaces::projectFace::calcLambdas
+(
+    const labelPair& n,
+    const pointField& points,
+    scalarField& lambdaI,
+    scalarField& lambdaJ
+) const
+{
+    lambdaI.setSize(points.size());
+    lambdaI = 0.0;
+    lambdaJ.setSize(points.size());
+    lambdaJ = 0.0;
+
+    for (label i = 1; i < n.first(); i++)
+    {
+        for (label j = 1; j < n.second(); j++)
+        {
+            label ij = index(n, labelPair(i, j));
+            label iMin1j = index(n, labelPair(i-1, j));
+            lambdaI[ij] = lambdaI[iMin1j] + mag(points[ij]-points[iMin1j]);
+
+            label ijMin1 = index(n, labelPair(i, j-1));
+            lambdaJ[ij] = lambdaJ[ijMin1] + mag(points[ij]-points[ijMin1]);
+        }
+    }
+
+    for (label i = 1; i < n.first(); i++)
+    {
+        label ijLast = index(n, labelPair(i, n.second()-1));
+        for (label j = 1; j < n.second(); j++)
+        {
+            label ij = index(n, labelPair(i, j));
+            lambdaJ[ij] /= lambdaJ[ijLast];
+        }
+    }
+    for (label j = 1; j < n.second(); j++)
+    {
+        label iLastj = index(n, labelPair(n.first()-1, j));
+        for (label i = 1; i < n.first(); i++)
+        {
+            label ij = index(n, labelPair(i, j));
+            lambdaI[ij] /= lambdaI[iLastj];
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::blockFaces::projectFace::projectFace
 (
+    const dictionary& dict,
+    const label index,
     const searchableSurfaces& geometry,
     Istream& is
 )
 :
-    blockFace(is),
+    blockFace(dict, index, is),
     surface_(lookupSurface(geometry, is))
 {}
 
@@ -88,20 +149,215 @@ void Foam::blockFaces::projectFace::project
     pointField& points
 ) const
 {
-    List<pointIndexHit> hits;
-    scalarField nearestDistSqr
-    (
-        points.size(),
-        magSqr(points[0] - points[points.size()-1])
-    );
-    surface_.findNearest(points, nearestDistSqr, hits);
+    // For debugging to tag the output
+    static label fIter = 0;
 
-    forAll(hits, i)
+    autoPtr<OBJstream> debugStr;
+    if (debug)
     {
-        if (hits[i].hit())
+        debugStr.reset
+        (
+            new OBJstream("projectFace_" + Foam::name(fIter++) + ".obj")
+        );
+        Info<< "Face:" << blockFacei << " on block:" << desc.blockShape()
+            << " with verts:" << desc.vertices()
+            << " writing lines from start points"
+            << " to projected points to " << debugStr().name() << endl;
+    }
+
+
+    // Find out range of vertices in face
+    labelPair n(-1, -1);
+    switch (blockFacei)
+    {
+        case 0:
+        case 1:
         {
-            points[i] = hits[i].hitPoint();
+            n.first() = desc.density()[1]+1;
+            n.second() = desc.density()[2]+1;
         }
+        break;
+
+        case 2:
+        case 3:
+        {
+            n.first() = desc.density()[0]+1;
+            n.second() = desc.density()[2]+1;
+        }
+        break;
+
+        case 4:
+        case 5:
+        {
+            n.first() = desc.density()[0]+1;
+            n.second() = desc.density()[1]+1;
+        }
+        break;
+    }
+
+
+    // Calculate initial normalised edge lengths (= u,v coordinates)
+    scalarField lambdaI(points.size(), 0.0);
+    scalarField lambdaJ(points.size(), 0.0);
+    calcLambdas(n, points, lambdaI, lambdaJ);
+
+
+    // Upper limit for number of iterations
+    const label maxIter = 10;
+    // Residual tolerance
+    const scalar relTol = 0.1;
+
+    scalar initialResidual = 0.0;
+    scalar iResidual = 0.0;
+    scalar jResidual = 0.0;
+
+    for (label iter = 0; iter < maxIter;  iter++)
+    {
+        // Do projection
+        {
+            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())
+                {
+                    const point& hitPt = hits[i].hitPoint();
+                    if (debugStr.valid())
+                    {
+                        debugStr().write(linePointRef(points[i], hitPt));
+                    }
+                    points[i] = hitPt;
+                }
+            }
+        }
+
+        if (debug)
+        {
+            Pout<< "Iter:" << iter << " initialResidual:" << initialResidual
+                << " iResidual+jResidual:" << iResidual+jResidual << endl;
+        }
+
+
+        if (iter > 0 && (iResidual+jResidual)/initialResidual < relTol)
+        {
+            break;
+        }
+
+
+        // Predict along i
+        vectorField residual(points.size(), vector::zero);
+
+        // Work arrays for interpolation
+        labelList indices;
+        scalarField weights;
+        for (label j = 1; j < n.second()-1; j++)
+        {
+            // Calculate actual lamdba along constant j line
+            scalarField projLambdas(n.first());
+            projLambdas[0] = 0.0;
+            for (label i = 1; i < n.first(); i++)
+            {
+                label ij = index(n, labelPair(i, j));
+                label iMin1j = index(n, labelPair(i-1, j));
+                projLambdas[i] =
+                    projLambdas[i-1]
+                   +mag(points[ij]-points[iMin1j]);
+            }
+            projLambdas /= projLambdas.last();
+
+            linearInterpolationWeights interpolator(projLambdas);
+
+            for (label i = 1; i < n.first()-1; i++)
+            {
+                label ij = index(n, labelPair(i, j));
+
+                interpolator.valueWeights(lambdaI[ij], indices, weights);
+
+                point predicted = vector::zero;
+                forAll(indices, indexi)
+                {
+                    label ptIndex = index(n, labelPair(indices[indexi], j));
+                    predicted += weights[indexi]*points[ptIndex];
+                }
+                residual[ij] = predicted-points[ij];
+            }
+        }
+
+        if (debugStr.valid())
+        {
+            forAll(points, i)
+            {
+                const linePointRef ln(points[i], points[i]+residual[i]);
+                debugStr().write(ln);
+            }
+        }
+
+        iResidual = sum(mag(residual));
+
+        // Update points before doing j. Note: is this needed? Complicates
+        // residual checking.
+        points += residual;
+
+
+        // Predict along j
+        residual = vector::zero;
+        for (label i = 1; i < n.first()-1; i++)
+        {
+            // Calculate actual lamdba along constant i line
+            scalarField projLambdas(n.second());
+            projLambdas[0] = 0.0;
+            for (label j = 1; j < n.second(); j++)
+            {
+                label ij = index(n, labelPair(i, j));
+                label ijMin1 = index(n, labelPair(i, j-1));
+                projLambdas[j] =
+                    projLambdas[j-1]
+                   +mag(points[ij]-points[ijMin1]);
+            }
+
+            projLambdas /= projLambdas.last();
+
+            linearInterpolationWeights interpolator(projLambdas);
+
+            for (label j = 1; j < n.second()-1; j++)
+            {
+                label ij = index(n, labelPair(i, j));
+
+                interpolator.valueWeights(lambdaJ[ij], indices, weights);
+
+                point predicted = vector::zero;
+                forAll(indices, indexi)
+                {
+                    label ptIndex = index(n, labelPair(i, indices[indexi]));
+                    predicted += weights[indexi]*points[ptIndex];
+                }
+                residual[ij] = predicted-points[ij];
+            }
+        }
+
+        if (debugStr.valid())
+        {
+            forAll(points, i)
+            {
+                const linePointRef ln(points[i], points[i]+residual[i]);
+                debugStr().write(ln);
+            }
+        }
+
+        jResidual = sum(mag(residual));
+
+        if (iter == 0)
+        {
+            initialResidual = iResidual + jResidual;
+        }
+
+        points += residual;
     }
 }
 
diff --git a/src/mesh/blockMesh/blockFaces/projectFace/projectFace.H b/src/mesh/blockMesh/blockFaces/projectFace/projectFace.H
index c1b56bbcdb..635318b05c 100644
--- a/src/mesh/blockMesh/blockFaces/projectFace/projectFace.H
+++ b/src/mesh/blockMesh/blockFaces/projectFace/projectFace.H
@@ -67,6 +67,22 @@ class projectFace
             Istream& is
         ) const;
 
+        //- Convert i,j to single index
+        label index
+        (
+            const labelPair& n,
+            const labelPair& coord
+        ) const;
+
+        //- Calulate lambdas (but unnormalised)
+        void calcLambdas
+        (
+            const labelPair& n,
+            const pointField& points,
+            scalarField& lambdaI,
+            scalarField& lambdaJ
+        ) const;
+
         //- Disallow default bitwise copy construct
         projectFace(const projectFace&);
 
@@ -85,6 +101,8 @@ public:
         //- Construct from Istream setting pointsList
         projectFace
         (
+            const dictionary& dict,
+            const label index,
             const searchableSurfaces& geometry,
             Istream&
         );
diff --git a/src/mesh/blockMesh/blockMesh/blockMesh.C b/src/mesh/blockMesh/blockMesh/blockMesh.C
index 366288e5ae..4a08adbfbe 100644
--- a/src/mesh/blockMesh/blockMesh/blockMesh.C
+++ b/src/mesh/blockMesh/blockMesh/blockMesh.C
@@ -59,7 +59,7 @@ Foam::blockMesh::blockMesh(const IOdictionary& dict, const word& regionName)
     blockVertices_
     (
         dict.lookup("vertices"),
-        blockVertex::iNew(geometry_)
+        blockVertex::iNew(dict, geometry_)
     ),
     vertices_(Foam::vertices(blockVertices_)),
     topologyPtr_(createTopology(dict, regionName))
diff --git a/src/mesh/blockMesh/blockMesh/blockMesh.H b/src/mesh/blockMesh/blockMesh/blockMesh.H
index fcd3698c34..093972d12d 100644
--- a/src/mesh/blockMesh/blockMesh/blockMesh.H
+++ b/src/mesh/blockMesh/blockMesh/blockMesh.H
@@ -139,7 +139,7 @@ class blockMesh
 
         polyMesh* createTopology(const IOdictionary&, const word& regionName);
 
-        void check(const polyMesh&) const;
+        void check(const polyMesh&, const dictionary&) const;
 
         //- Determine the merge info and the final number of cells/points
         void calcMergeInfo();
diff --git a/src/mesh/blockMesh/blockMesh/blockMeshCheck.C b/src/mesh/blockMesh/blockMesh/blockMeshCheck.C
index 10bd9e3116..fc0c08b862 100644
--- a/src/mesh/blockMesh/blockMesh/blockMeshCheck.C
+++ b/src/mesh/blockMesh/blockMesh/blockMeshCheck.C
@@ -27,7 +27,7 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::blockMesh::check(const polyMesh& bm) const
+void Foam::blockMesh::check(const polyMesh& bm, const dictionary& dict) const
 {
     Info<< nl << "Check topology" << endl;
 
@@ -40,8 +40,9 @@ void Foam::blockMesh::check(const polyMesh& bm) const
         {
             if (edges_[cei].compare(edges_[cej]) != 0)
             {
-                Info<< "    Curved edge " << edges_[cej]
-                    << "    is a duplicate of curved edge " << edges_[cei]
+                Info<< "    Curved edge ";
+                edges_[cej].write(Info, dict);
+                Info<< "    is a duplicate of curved edge " << edges_[cei]
                     << endl;
                 ok = false;
                 break;
@@ -74,8 +75,9 @@ void Foam::blockMesh::check(const polyMesh& bm) const
 
         if (!found)
         {
-            Info<< "    Curved edge " << edges_[cei]
-                << "    does not correspond to a block edge."
+            Info<< "    Curved edge ";
+            edges_[cei].write(Info, dict);
+            Info<< "    does not correspond to a block edge."
                 << endl;
             ok = false;
         }
@@ -90,9 +92,11 @@ void Foam::blockMesh::check(const polyMesh& bm) const
         {
             if (faces_[cfi].compare(faces_[cfj]) != 0)
             {
-                Info<< "    Curved face " << faces_[cfj]
-                    << "    is a duplicate of curved face " << faces_[cfi]
-                    << endl;
+                Info<< "    Curved face ";
+                faces_[cfj].write(Info, dict);
+                Info<< "    is a duplicate of curved face ";
+                faces_[cfi].write(Info, dict);
+                Info<< endl;
                 ok = false;
                 break;
             }
@@ -112,9 +116,9 @@ void Foam::blockMesh::check(const polyMesh& bm) const
 
         if (!found)
         {
-            Info<< "    Curved face " << faces_[cfi]
-                << "    does not correspond to a block face."
-                << endl;
+            Info<< "    Curved face ";
+            faces_[cfi].write(Info, dict);
+            Info<< "    does not correspond to a block face." << endl;
             ok = false;
         }
     }
diff --git a/src/mesh/blockMesh/blockMesh/blockMeshTopology.C b/src/mesh/blockMesh/blockMesh/blockMeshTopology.C
index ab99844205..a4c9acead3 100644
--- a/src/mesh/blockMesh/blockMesh/blockMeshTopology.C
+++ b/src/mesh/blockMesh/blockMesh/blockMeshTopology.C
@@ -110,6 +110,11 @@ void Foam::blockMesh::readPatches
     wordList& nbrPatchNames
 )
 {
+    // Collect all variables
+    dictionary varDict(meshDescription.subOrEmptyDict("namedVertices"));
+    varDict.merge(meshDescription.subOrEmptyDict("namedBlocks"));
+
+
     ITstream& patchStream(meshDescription.lookup("patches"));
 
     // Read number of patches in mesh
@@ -160,7 +165,11 @@ void Foam::blockMesh::readPatches
             >> patchNames[nPatches];
 
         // Read patch faces
-        patchStream >> tmpBlocksPatches[nPatches];
+        tmpBlocksPatches[nPatches] = blockDescriptor::read<face>
+        (
+            patchStream,
+            varDict
+        );
 
 
         // Check for multiple patches
@@ -258,6 +267,11 @@ void Foam::blockMesh::readBoundary
     PtrList<dictionary>& patchDicts
 )
 {
+    // Collect all variables
+    dictionary varDict(meshDescription.subOrEmptyDict("namedVertices"));
+    varDict.merge(meshDescription.subOrEmptyDict("namedBlocks"));
+
+
     // Read like boundary file
     const PtrList<entry> patchesInfo
     (
@@ -286,7 +300,12 @@ void Foam::blockMesh::readBoundary
         patchDicts.set(patchi, new dictionary(patchInfo.dict()));
 
         // Read block faces
-        patchDicts[patchi].lookup("faces") >> tmpBlocksPatches[patchi];
+        tmpBlocksPatches[patchi] = blockDescriptor::read<face>
+        (
+            patchDicts[patchi].lookup("faces"),
+            varDict
+        );
+
 
         checkPatchLabels
         (
@@ -353,7 +372,7 @@ Foam::polyMesh* Foam::blockMesh::createTopology
         blockEdgeList edges
         (
             meshDescription.lookup("edges"),
-            blockEdge::iNew(geometry_, vertices_)
+            blockEdge::iNew(meshDescription, geometry_, vertices_)
         );
 
         edges_.transfer(edges);
@@ -375,7 +394,7 @@ Foam::polyMesh* Foam::blockMesh::createTopology
         blockFaceList faces
         (
             meshDescription.lookup("faces"),
-            blockFace::iNew(geometry_)
+            blockFace::iNew(meshDescription, geometry_)
         );
 
         faces_.transfer(faces);
@@ -395,7 +414,7 @@ Foam::polyMesh* Foam::blockMesh::createTopology
         blockList blocks
         (
             meshDescription.lookup("blocks"),
-            block::iNew(vertices_, edges_, faces_)
+            block::iNew(meshDescription, vertices_, edges_, faces_)
         );
 
         transfer(blocks);
@@ -545,7 +564,7 @@ Foam::polyMesh* Foam::blockMesh::createTopology
         );
     }
 
-    check(*blockMeshPtr);
+    check(*blockMeshPtr, meshDescription);
 
     return blockMeshPtr;
 }
diff --git a/src/mesh/blockMesh/blockVertices/blockVertex/blockVertex.C b/src/mesh/blockMesh/blockVertices/blockVertex/blockVertex.C
index fe8ab41dc5..6bdf9c2ade 100644
--- a/src/mesh/blockMesh/blockVertices/blockVertex/blockVertex.C
+++ b/src/mesh/blockMesh/blockVertices/blockVertex/blockVertex.C
@@ -50,6 +50,8 @@ Foam::autoPtr<Foam::blockVertex> Foam::blockVertex::clone() const
 
 Foam::autoPtr<Foam::blockVertex> Foam::blockVertex::New
 (
+    const dictionary& dict,
+    const label index,
     const searchableSurfaces& geometry,
     Istream& is
 )
@@ -61,14 +63,14 @@ Foam::autoPtr<Foam::blockVertex> Foam::blockVertex::New
 
     token firstToken(is);
 
-    if (firstToken.pToken() == token::BEGIN_LIST)
+    if (firstToken.isPunctuation() && firstToken.pToken() == token::BEGIN_LIST)
     {
         // Putback the opening bracket
         is.putBack(firstToken);
 
         return autoPtr<blockVertex>
         (
-            new blockVertices::pointVertex(geometry, is)
+            new blockVertices::pointVertex(dict, index, geometry, is)
         );
     }
     else if (firstToken.isWord())
@@ -88,7 +90,7 @@ Foam::autoPtr<Foam::blockVertex> Foam::blockVertex::New
                 << abort(FatalError);
         }
 
-        return autoPtr<blockVertex>(cstrIter()(geometry, is));
+        return autoPtr<blockVertex>(cstrIter()(dict, index, geometry, is));
     }
     else
     {
diff --git a/src/mesh/blockMesh/blockVertices/blockVertex/blockVertex.H b/src/mesh/blockMesh/blockVertices/blockVertex/blockVertex.H
index 7d5bb9fbd6..5c6c8006a1 100644
--- a/src/mesh/blockMesh/blockVertices/blockVertex/blockVertex.H
+++ b/src/mesh/blockMesh/blockVertices/blockVertex/blockVertex.H
@@ -63,10 +63,12 @@ public:
             blockVertex,
             Istream,
             (
+                const dictionary& dict,
+                const label index,
                 const searchableSurfaces& geometry,
                 Istream& is
             ),
-            (geometry, is)
+            (dict, index, geometry, is)
         );
 
 
@@ -81,6 +83,8 @@ public:
         //- New function which constructs and returns pointer to a blockVertex
         static autoPtr<blockVertex> New
         (
+            const dictionary& dict,
+            const label index,
             const searchableSurfaces& geometry,
             Istream&
         );
@@ -89,18 +93,22 @@ public:
         //  PtrLists of blockVertex
         class iNew
         {
+            const dictionary& dict_;
             const searchableSurfaces& geometry_;
+            mutable label index_;
 
         public:
 
-            iNew(const searchableSurfaces& geometry)
+            iNew(const dictionary& dict, const searchableSurfaces& geometry)
             :
-                geometry_(geometry)
+                dict_(dict),
+                geometry_(geometry),
+                index_(0)
             {}
 
             autoPtr<blockVertex> operator()(Istream& is) const
             {
-                return blockVertex::New(geometry_, is);
+                return blockVertex::New(dict_, index_++, geometry_, is);
             }
         };
 
diff --git a/src/mesh/blockMesh/blockVertices/namedVertex/namedVertex.C b/src/mesh/blockMesh/blockVertices/namedVertex/namedVertex.C
new file mode 100644
index 0000000000..b4392fac26
--- /dev/null
+++ b/src/mesh/blockMesh/blockVertices/namedVertex/namedVertex.C
@@ -0,0 +1,81 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "namedVertex.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace blockVertices
+{
+    defineTypeNameAndDebug(namedVertex, 0);
+    addToRunTimeSelectionTable(blockVertex, namedVertex, Istream);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::blockVertices::namedVertex::namedVertex
+(
+    const dictionary& dict,
+    const label index,
+    const searchableSurfaces& geometry,
+    Istream& is
+)
+:
+    name_(is),
+    vertexPtr_(blockVertex::New(dict, index, geometry, is))
+{
+    //Info<< "Vertex " << name_ << " at " <<  vertexPtr_().operator point()
+    //    << " has index " << index << endl;
+
+    dictionary& d = const_cast<dictionary&>(dict);
+
+    const dictionary* varDictPtr = d.subDictPtr("namedVertices");
+    if (varDictPtr)
+    {
+        const_cast<dictionary&>(*varDictPtr).add(name_, index);
+    }
+    else
+    {
+        dictionary varDict;
+        varDict.add(name_, index);
+        d.add("namedVertices", varDict);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::blockVertices::namedVertex::operator point() const
+{
+    return vertexPtr_().operator point();
+}
+
+
+// ************************************************************************* //
diff --git a/src/mesh/blockMesh/blockVertices/namedVertex/namedVertex.H b/src/mesh/blockMesh/blockVertices/namedVertex/namedVertex.H
new file mode 100644
index 0000000000..117107ad5f
--- /dev/null
+++ b/src/mesh/blockMesh/blockVertices/namedVertex/namedVertex.H
@@ -0,0 +1,104 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::blockVertices::namedVertex
+
+Description
+    Gives name to a vertex.
+
+SourceFiles
+    namedVertex.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef blockVertices_namedVertex_H
+#define blockVertices_namedVertex_H
+
+#include "blockVertex.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace blockVertices
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class namedVertex Declaration
+\*---------------------------------------------------------------------------*/
+
+class namedVertex
+:
+    public blockVertex
+{
+protected:
+
+    // Protected member data
+
+        //- The dictionary variable name for the vertex number
+        word name_;
+
+        //- The vertex location
+        autoPtr<blockVertex> vertexPtr_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("named");
+
+
+    // Constructors
+
+        //- Construct from Istream setting pointsList
+        namedVertex
+        (
+            const dictionary&,
+            const label index,
+            const searchableSurfaces& geometry,
+            Istream&
+        );
+
+
+    //- Destructor
+    virtual ~namedVertex()
+    {}
+
+
+    // Member Functions
+
+        virtual operator point() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace blockVertices
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/mesh/blockMesh/blockVertices/pointVertex/pointVertex.C b/src/mesh/blockMesh/blockVertices/pointVertex/pointVertex.C
index 4b3120d98b..4af95df163 100644
--- a/src/mesh/blockMesh/blockVertices/pointVertex/pointVertex.C
+++ b/src/mesh/blockMesh/blockVertices/pointVertex/pointVertex.C
@@ -42,6 +42,8 @@ namespace blockVertices
 
 Foam::blockVertices::pointVertex::pointVertex
 (
+    const dictionary&,
+    const label index,
     const searchableSurfaces& geometry,
     Istream& is
 )
diff --git a/src/mesh/blockMesh/blockVertices/pointVertex/pointVertex.H b/src/mesh/blockMesh/blockVertices/pointVertex/pointVertex.H
index 1179c66bd5..929f92a25c 100644
--- a/src/mesh/blockMesh/blockVertices/pointVertex/pointVertex.H
+++ b/src/mesh/blockMesh/blockVertices/pointVertex/pointVertex.H
@@ -70,6 +70,8 @@ public:
         //- Construct from Istream setting pointsList
         pointVertex
         (
+            const dictionary&,
+            const label index,
             const searchableSurfaces& geometry,
             Istream&
         );
diff --git a/src/mesh/blockMesh/blockVertices/projectVertex/projectVertex.C b/src/mesh/blockMesh/blockVertices/projectVertex/projectVertex.C
index 04cff7d250..f5d23fc50c 100644
--- a/src/mesh/blockMesh/blockVertices/projectVertex/projectVertex.C
+++ b/src/mesh/blockMesh/blockVertices/projectVertex/projectVertex.C
@@ -45,11 +45,13 @@ namespace blockVertices
 
 Foam::blockVertices::projectVertex::projectVertex
 (
+    const dictionary& dict,
+    const label index,
     const searchableSurfaces& geometry,
     Istream& is
 )
 :
-    pointVertex(geometry, is),
+    pointVertex(dict, index, geometry, is),
     geometry_(geometry)
 {
     wordList names(is);
@@ -79,8 +81,11 @@ Foam::blockVertices::projectVertex::operator point() const
 
 
     // Note: how far do we need to search? Probably not further than
-    //       span of surfaces themselves.
+    //       span of surfaces themselves. Make sure to limit in case
+    //       of e.g. searchablePlane which has infinite bb.
     boundBox bb(searchableSurfacesQueries::bounds(geometry_, surfaces_));
+    bb.min() = max(bb.min(), point(-GREAT, -GREAT, -GREAT));
+    bb.max() = min(bb.max(), point(GREAT, GREAT, GREAT));
 
     searchableSurfacesQueries::findNearest
     (
diff --git a/src/mesh/blockMesh/blockVertices/projectVertex/projectVertex.H b/src/mesh/blockMesh/blockVertices/projectVertex/projectVertex.H
index 929b6a68eb..d47dce6e8c 100644
--- a/src/mesh/blockMesh/blockVertices/projectVertex/projectVertex.H
+++ b/src/mesh/blockMesh/blockVertices/projectVertex/projectVertex.H
@@ -81,6 +81,8 @@ public:
         //- Construct from Istream setting pointsList
         projectVertex
         (
+            const dictionary&,
+            const label index,
             const searchableSurfaces& geometry,
             Istream&
         );
diff --git a/src/mesh/blockMesh/block/block.C b/src/mesh/blockMesh/blocks/block/block.C
similarity index 61%
rename from src/mesh/blockMesh/block/block.C
rename to src/mesh/blockMesh/blocks/block/block.C
index 1722b0095e..5d75aa5574 100644
--- a/src/mesh/blockMesh/block/block.C
+++ b/src/mesh/blockMesh/blocks/block/block.C
@@ -25,17 +25,27 @@ License
 
 #include "block.H"
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(block, 0);
+    defineRunTimeSelectionTable(block, Istream);
+}
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::block::block
 (
+    const dictionary& dict,
+    const label index,
     const pointField& vertices,
     const blockEdgeList& edges,
     const blockFaceList& faces,
     Istream& is
 )
 :
-    blockDescriptor(vertices, edges, faces, is)
+    blockDescriptor(dict, index, vertices, edges, faces, is)
 {
     createPoints();
     createBoundary();
@@ -51,6 +61,49 @@ Foam::block::block(const blockDescriptor& blockDesc)
 }
 
 
+Foam::autoPtr<Foam::block> Foam::block::New
+(
+    const dictionary& dict,
+    const label index,
+    const pointField& points,
+    const blockEdgeList& edges,
+    const blockFaceList& faces,
+    Istream& is
+)
+{
+    if (debug)
+    {
+        InfoInFunction << "Constructing block" << endl;
+    }
+
+    const word blockOrCellShapeType(is);
+
+    IstreamConstructorTable::iterator cstrIter =
+        IstreamConstructorTablePtr_->find(blockOrCellShapeType);
+
+    if (cstrIter == IstreamConstructorTablePtr_->end())
+    {
+        is.putBack(blockOrCellShapeType);
+        return autoPtr<block>(new block(dict, index, points, edges, faces, is));
+    }
+    else
+    {
+        return autoPtr<block>
+        (
+            cstrIter()
+            (
+                dict,
+                index,
+                points,
+                edges,
+                faces,
+                is
+            )
+        );
+    }
+}
+
+
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 
 Foam::Ostream& Foam::operator<<(Ostream& os, const block& b)
diff --git a/src/mesh/blockMesh/block/block.H b/src/mesh/blockMesh/blocks/block/block.H
similarity index 78%
rename from src/mesh/blockMesh/block/block.H
rename to src/mesh/blockMesh/blocks/block/block.H
index 7b55c5f59f..fb3012016d 100644
--- a/src/mesh/blockMesh/block/block.H
+++ b/src/mesh/blockMesh/blocks/block/block.H
@@ -91,11 +91,36 @@ class block
 
 public:
 
+    //- Runtime type information
+    TypeName("block");
+
+
+    // Declare run-time constructor selection tables
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            block,
+            Istream,
+            (
+                const dictionary& dict,
+                const label index,
+                const pointField& vertices,
+                const blockEdgeList& edges,
+                const blockFaceList& faces,
+                Istream& is
+            ),
+            (dict, index, vertices, edges, faces, is)
+        );
+
+
     // Constructors
 
         //- Construct from components with Istream
         block
         (
+            const dictionary& dict,
+            const label index,
             const pointField& vertices,
             const blockEdgeList& edges,
             const blockFaceList& faces,
@@ -112,35 +137,56 @@ public:
             return autoPtr<block>(nullptr);
         }
 
+        //- New function which constructs and returns pointer to a block
+        static autoPtr<block> New
+        (
+            const dictionary& dict,
+            const label index,
+            const pointField& points,
+            const blockEdgeList& edges,
+            const blockFaceList& faces,
+            Istream&
+        );
+
         //- Class used for the read-construction of
         //  PtrLists of blocks
         class iNew
         {
+            const dictionary& dict_;
             const pointField& points_;
             const blockEdgeList& edges_;
             const blockFaceList& faces_;
+            mutable label index_;
 
         public:
 
             iNew
             (
+                const dictionary& dict,
                 const pointField& points,
                 const blockEdgeList& edges,
                 const blockFaceList& faces
             )
             :
+                dict_(dict),
                 points_(points),
                 edges_(edges),
-                faces_(faces)
+                faces_(faces),
+                index_(0)
             {}
 
             autoPtr<block> operator()(Istream& is) const
             {
-                return autoPtr<block>(new block(points_, edges_, faces_, is));
+                return block::New(dict_, index_++, points_, edges_, faces_, is);
             }
         };
 
 
+    //- Destructor
+    virtual ~block()
+    {}
+
+
     // Member Functions
 
         // Access
diff --git a/src/mesh/blockMesh/block/blockCreate.C b/src/mesh/blockMesh/blocks/block/blockCreate.C
similarity index 100%
rename from src/mesh/blockMesh/block/blockCreate.C
rename to src/mesh/blockMesh/blocks/block/blockCreate.C
diff --git a/src/mesh/blockMesh/block/blockI.H b/src/mesh/blockMesh/blocks/block/blockI.H
similarity index 100%
rename from src/mesh/blockMesh/block/blockI.H
rename to src/mesh/blockMesh/blocks/block/blockI.H
diff --git a/src/mesh/blockMesh/block/blockList.H b/src/mesh/blockMesh/blocks/block/blockList.H
similarity index 95%
rename from src/mesh/blockMesh/block/blockList.H
rename to src/mesh/blockMesh/blocks/block/blockList.H
index bfca6c9c60..80b8a328af 100644
--- a/src/mesh/blockMesh/block/blockList.H
+++ b/src/mesh/blockMesh/blocks/block/blockList.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -53,4 +53,3 @@ typedef PtrList<block> blockList;
 #endif
 
 // ************************************************************************* //
-
diff --git a/src/mesh/blockMesh/blocks/namedBlock/namedBlock.C b/src/mesh/blockMesh/blocks/namedBlock/namedBlock.C
new file mode 100644
index 0000000000..9408e4fd3e
--- /dev/null
+++ b/src/mesh/blockMesh/blocks/namedBlock/namedBlock.C
@@ -0,0 +1,75 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "namedBlock.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace blocks
+{
+    defineTypeNameAndDebug(namedBlock, 0);
+    addToRunTimeSelectionTable(block, namedBlock, Istream);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::blocks::namedBlock::namedBlock
+(
+    const dictionary& dict,
+    const label index,
+    const pointField& vertices,
+    const blockEdgeList& edges,
+    const blockFaceList& faces,
+    Istream& is
+)
+:
+    word(is),
+    block(dict, index, vertices, edges, faces, is)
+{
+    //Info<< "Block " << static_cast<word&>(*this)
+    //    << " has index " << index << endl;
+
+    dictionary& d = const_cast<dictionary&>(dict);
+
+    const dictionary* varDictPtr = d.subDictPtr("namedBlocks");
+    if (varDictPtr)
+    {
+        const_cast<dictionary&>(*varDictPtr).add(*this, index);
+    }
+    else
+    {
+        dictionary varDict;
+        varDict.add(*this, index);
+        d.add("namedBlocks", varDict);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/mesh/blockMesh/blocks/namedBlock/namedBlock.H b/src/mesh/blockMesh/blocks/namedBlock/namedBlock.H
new file mode 100644
index 0000000000..98968f60b0
--- /dev/null
+++ b/src/mesh/blockMesh/blocks/namedBlock/namedBlock.H
@@ -0,0 +1,91 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::blocks::namedBlock
+
+Description
+    Gives name to a block.
+
+SourceFiles
+    namedBlock.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef blocks_namedBlock_H
+#define blocks_namedBlock_H
+
+#include "block.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace blocks
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class namedBlock Declaration
+\*---------------------------------------------------------------------------*/
+
+class namedBlock
+:
+    public word,
+    public block
+{
+public:
+
+    //- Runtime type information
+    TypeName("named");
+
+
+    // Constructors
+
+        //- Construct from Istream setting pointsList
+        namedBlock
+        (
+            const dictionary& dict,
+            const label index,
+            const pointField& vertices,
+            const blockEdgeList& edges,
+            const blockFaceList& faces,
+            Istream& is
+        );
+
+
+    //- Destructor
+    virtual ~namedBlock()
+    {}
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace blocks
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/mesh/blockMesh/searchableCurve/searchableCurve.C b/src/mesh/blockMesh/searchableCurve/searchableCurve.C
new file mode 100644
index 0000000000..b8cc234d0e
--- /dev/null
+++ b/src/mesh/blockMesh/searchableCurve/searchableCurve.C
@@ -0,0 +1,432 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "searchableCurve.H"
+#include "addToRunTimeSelectionTable.H"
+#include "Time.H"
+#include "edgeMesh.H"
+#include "indexedOctree.H"
+#include "treeDataEdge.H"
+#include "linearInterpolationWeights.H"
+#include "quaternion.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(searchableCurve, 0);
+    addToRunTimeSelectionTable(searchableSurface, searchableCurve, dict);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::searchableCurve::searchableCurve
+(
+    const IOobject& io,
+    const dictionary& dict
+)
+:
+    searchableSurface(io),
+    eMeshPtr_
+    (
+        edgeMesh::New
+        (
+            IOobject
+            (
+                dict.lookup("file"),                // name
+                io.time().constant(),               // instance
+                "triSurface",                       // local
+                io.time(),                          // registry
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE,
+                false
+            ).objectPath()
+        )
+    ),
+    radius_(readScalar(dict.lookup("radius")))
+{
+    const edgeMesh& eMesh = eMeshPtr_();
+
+    const pointField& points = eMesh.points();
+    const edgeList& edges = eMesh.edges();
+    bounds() = boundBox(points, false);
+
+    vector halfSpan(0.5*bounds().span());
+    point ctr(bounds().midpoint());
+
+    bounds().min() = ctr - mag(halfSpan)*vector(1, 1, 1);
+    bounds().max() = ctr + mag(halfSpan)*vector(1, 1, 1);
+
+    // Calculate bb of all points
+    treeBoundBox bb(bounds());
+
+    // Slightly extended bb. Slightly off-centred just so on symmetric
+    // geometry there are less face/edge aligned items.
+    bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+    bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+
+    edgeTree_.reset
+    (
+        new indexedOctree<treeDataEdge>
+        (
+            treeDataEdge
+            (
+                false,                  // do not cache bb
+                edges,
+                points,
+                identity(edges.size())
+            ),
+            bb,     // overall search domain
+            8,      // maxLevel
+            10,     // leafsize
+            3.0     // duplicity
+        )
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::searchableCurve::~searchableCurve()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::wordList& Foam::searchableCurve::regions() const
+{
+    if (regions_.empty())
+    {
+        regions_.setSize(1);
+        regions_[0] = "region0";
+    }
+    return regions_;
+}
+
+
+Foam::label Foam::searchableCurve::size() const
+{
+    return eMeshPtr_().points().size();
+}
+
+
+Foam::tmp<Foam::pointField> Foam::searchableCurve::coordinates() const
+{
+    return eMeshPtr_().points();
+}
+
+
+void Foam::searchableCurve::boundingSpheres
+(
+    pointField& centres,
+    scalarField& radiusSqr
+) const
+{
+    centres = eMeshPtr_().points();
+    radiusSqr.setSize(centres.size());
+    radiusSqr = Foam::sqr(radius_);
+    // Add a bit to make sure all points are tested inside
+    radiusSqr += Foam::sqr(SMALL);
+}
+
+
+void Foam::searchableCurve::findNearest
+(
+    const pointField& samples,
+    const scalarField& nearestDistSqr,
+    List<pointIndexHit>& info
+) const
+{
+    const indexedOctree<treeDataEdge>& tree = edgeTree_();
+
+    info.setSize(samples.size());
+
+    forAll(samples, i)
+    {
+        info[i] = tree.findNearest(samples[i], nearestDistSqr[i]);
+
+        if (info[i].hit())
+        {
+            vector d(samples[i]-info[i].hitPoint());
+            info[i].setPoint(info[i].hitPoint() + d/mag(d)*radius_);
+        }
+    }
+}
+
+
+void Foam::searchableCurve::findNearest
+(
+    const point& start,
+    const point& end,
+    const scalarField& rawLambdas,
+    const scalarField& nearestDistSqr,
+    List<pointIndexHit>& info
+) const
+{
+    const edgeMesh& mesh = eMeshPtr_();
+    const indexedOctree<treeDataEdge>& tree = edgeTree_();
+    const edgeList& edges = mesh.edges();
+    const pointField& points = mesh.points();
+    const labelListList& pointEdges = mesh.pointEdges();
+
+    const scalar maxDistSqr(Foam::magSqr(bounds().span()));
+
+    // Normalise lambdas
+    const scalarField lambdas
+    (
+        (rawLambdas-rawLambdas[0])
+       /(rawLambdas.last()-rawLambdas[0])
+    );
+
+
+    // Nearest point on curve and local axis direction
+    pointField curvePoints(lambdas.size());
+    vectorField axialVecs(lambdas.size());
+
+    const pointIndexHit startInfo = tree.findNearest(start, maxDistSqr);
+    curvePoints[0] = startInfo.hitPoint();
+    axialVecs[0] = edges[startInfo.index()].vec(points);
+
+    const pointIndexHit endInfo = tree.findNearest(end, maxDistSqr);
+    curvePoints.last() = endInfo.hitPoint();
+    axialVecs.last() = edges[endInfo.index()].vec(points);
+
+
+
+    scalarField curveLambdas(points.size(), -1.0);
+
+    {
+        scalar endDistance = -1.0;
+
+        // Determine edge lengths walking from start to end.
+
+        const point& start = curvePoints[0];
+        const point& end = curvePoints.last();
+
+        label edgei = startInfo.index();
+        const edge& startE = edges[edgei];
+
+        label pointi = startE[0];
+        if ((startE.vec(points)&(end-start)) > 0)
+        {
+            pointi = startE[1];
+        }
+
+        curveLambdas[pointi] = mag(points[pointi]-curvePoints[0]);
+        label otherPointi = startE.otherVertex(pointi);
+        curveLambdas[otherPointi] = -mag(points[otherPointi]-curvePoints[0]);
+
+        //Pout<< "for point:" << points[pointi] << " have distance "
+        //    << curveLambdas[pointi] << endl;
+
+
+        while (true)
+        {
+            const labelList& pEdges = pointEdges[pointi];
+            if (pEdges.size() == 1)
+            {
+                break;
+            }
+            else if (pEdges.size() != 2)
+            {
+                FatalErrorInFunction << "Curve " << name()
+                    << " is not a single path. This is not supported"
+                    << exit(FatalError);
+                break;
+            }
+
+            label oldEdgei = edgei;
+            if (pEdges[0] == oldEdgei)
+            {
+                edgei = pEdges[1];
+            }
+            else
+            {
+                edgei = pEdges[0];
+            }
+
+            if (edgei == endInfo.index())
+            {
+                endDistance = curveLambdas[pointi] + mag(end-points[pointi]);
+
+                //Pout<< "Found end edge:" << edges[edgei].centre(points)
+                //    << " endPt:" << end
+                //    << " point before:" << points[pointi]
+                //    << " accumulated length:" << endDistance << endl;
+            }
+
+
+            label oldPointi = pointi;
+            pointi = edges[edgei].otherVertex(oldPointi);
+
+            if (curveLambdas[pointi] >= 0)
+            {
+                break;
+            }
+
+            curveLambdas[pointi] =
+                curveLambdas[oldPointi] + edges[edgei].mag(points);
+        }
+
+        // Normalise curveLambdas
+        forAll(curveLambdas, i)
+        {
+            if (curveLambdas[i] >= 0)
+            {
+                curveLambdas[i] /= endDistance;
+            }
+        }
+    }
+
+
+
+    // Interpolation engine
+    linearInterpolationWeights interpolator(curveLambdas);
+
+    // Find wanted location along curve
+    labelList indices;
+    scalarField weights;
+    for (label i = 1; i < curvePoints.size()-1; i++)
+    {
+        interpolator.valueWeights(lambdas[i], indices, weights);
+
+        if (indices.size() == 1)
+        {
+            // On outside of curve. Choose one of the connected edges.
+            label pointi = indices[0];
+            const point& p0 = points[pointi];
+            label edge0 = pointEdges[pointi][0];
+            const edge& e0 = edges[edge0];
+            axialVecs[i] = e0.vec(points);
+            curvePoints[i] = weights[0]*p0;
+        }
+        else if (indices.size() == 2)
+        {
+            const point& p0 = points[indices[0]];
+            const point& p1 = points[indices[1]];
+            axialVecs[i] = p1-p0;
+            curvePoints[i] = weights[0]*p0+weights[1]*p1;
+        }
+    }
+    axialVecs /= mag(axialVecs);
+
+
+    // Now we have the lambdas, curvePoints and axialVecs.
+
+
+
+    info.setSize(lambdas.size());
+    info = pointIndexHit();
+
+    // Given the current lambdas interpolate radial direction inbetween
+    // endpoints (all projected onto the starting coordinate system)
+    quaternion qStart;
+    vector radialStart;
+    {
+        radialStart = start-curvePoints[0];
+        radialStart -= (radialStart&axialVecs[0])*axialVecs[0];
+        radialStart /= mag(radialStart);
+        qStart = quaternion(radialStart, 0.0);
+
+        info[0] = pointIndexHit(true, start, 0);
+    }
+
+    quaternion qProjectedEnd;
+    {
+        vector radialEnd(end-curvePoints.last());
+        radialEnd -= (radialEnd&axialVecs.last())*axialVecs.last();
+        radialEnd /= mag(radialEnd);
+
+        vector projectedEnd = radialEnd;
+        projectedEnd -= (projectedEnd&axialVecs[0])*axialVecs[0];
+        projectedEnd /= mag(projectedEnd);
+        qProjectedEnd = quaternion(projectedEnd, 0.0);
+
+        info.last() = pointIndexHit(true, end, 0);
+    }
+
+    for (label i = 1; i < lambdas.size()-1; i++)
+    {
+        quaternion q(slerp(qStart, qProjectedEnd, lambdas[i]));
+        vector radialDir(q.transform(radialStart));
+
+        radialDir -= (radialDir&axialVecs[i])*axialVecs.last();
+        radialDir /= mag(radialDir);
+
+        info[i] = pointIndexHit(true, curvePoints[i]+radius_*radialDir, 0);
+    }
+}
+
+
+void Foam::searchableCurve::getRegion
+(
+    const List<pointIndexHit>& info,
+    labelList& region
+) const
+{
+    region.setSize(info.size());
+    region = 0;
+}
+
+
+void Foam::searchableCurve::getNormal
+(
+    const List<pointIndexHit>& info,
+    vectorField& normal
+) const
+{
+    const edgeMesh& mesh = eMeshPtr_();
+    const indexedOctree<treeDataEdge>& tree = edgeTree_();
+    const edgeList& edges = mesh.edges();
+    const pointField& points = mesh.points();
+
+    normal.setSize(info.size());
+    normal = Zero;
+
+    forAll(info, i)
+    {
+        if (info[i].hit())
+        {
+            // Find nearest on curve
+            pointIndexHit curvePt = tree.findNearest
+            (
+                info[i].hitPoint(),
+                Foam::magSqr(bounds().span())
+            );
+
+            normal[i] = info[i].hitPoint()-curvePt.hitPoint();
+
+            // Subtract axial direction
+            vector axialVec = edges[curvePt.index()].vec(points);
+            axialVec /= mag(axialVec);
+            normal[i] -= (normal[i]&axialVec)*axialVec;
+            normal[i] /= mag(normal[i]);
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/mesh/blockMesh/searchableCurve/searchableCurve.H b/src/mesh/blockMesh/searchableCurve/searchableCurve.H
new file mode 100644
index 0000000000..3c5f122cbf
--- /dev/null
+++ b/src/mesh/blockMesh/searchableCurve/searchableCurve.H
@@ -0,0 +1,245 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::searchableCurve
+
+Description
+    Searching on edgemesh with constant radius
+
+SourceFiles
+    searchableCurve.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef searchableCurve_H
+#define searchableCurve_H
+
+#include "treeBoundBox.H"
+#include "searchableSurface.H"
+//#include "edgeMesh.H"
+//#include "indexedOctree.H"
+//#include "treeDataEdge.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class edgeMesh;
+class treeDataEdge;
+template <class Type> class indexedOctree;
+
+/*---------------------------------------------------------------------------*\
+                           Class searchableCurve Declaration
+\*---------------------------------------------------------------------------*/
+
+class searchableCurve
+:
+    public searchableSurface
+{
+private:
+
+    // Private Member Data
+
+        //- Feature
+        autoPtr<edgeMesh> eMeshPtr_;
+
+        //- Search structure
+        autoPtr<indexedOctree<treeDataEdge>> edgeTree_;
+
+        //- Radius
+        const scalar radius_;
+
+        //- Names of regions
+        mutable wordList regions_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        searchableCurve(const searchableCurve&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const searchableCurve&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("searchableCurve");
+
+
+    // Constructors
+
+        //- Construct from dictionary (used by searchableSurface)
+        searchableCurve
+        (
+            const IOobject& io,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~searchableCurve();
+
+
+    // Member Functions
+
+        virtual const wordList& regions() const;
+
+        //- Whether supports volume type below
+        virtual bool hasVolumeType() const
+        {
+            return true;
+        }
+
+        //- Range of local indices that can be returned.
+        virtual label size() const;
+
+        //- Get representative set of element coordinates
+        //  Usually the element centres (should be of length size()).
+        virtual tmp<pointField> coordinates() const;
+
+        //- Get bounding spheres (centre and radius squared), one per element.
+        //  Any point on element is guaranteed to be inside.
+        virtual void boundingSpheres
+        (
+            pointField& centres,
+            scalarField& radiusSqr
+        ) const;
+
+        //- Get the points that define the surface.
+        virtual tmp<pointField> points() const
+        {
+            return coordinates();
+        }
+
+        //- Does any part of the surface overlap the supplied bound box?
+        virtual bool overlaps(const boundBox& bb) const
+        {
+            NotImplemented;
+            return false;
+        }
+
+
+        // Multiple point queries.
+
+            virtual void findNearest
+            (
+                const pointField& sample,
+                const scalarField& nearestDistSqr,
+                List<pointIndexHit>&
+            ) const;
+
+            //- Unique to parametric geometry: given points find
+            //  an interpolated (along the curve) point on the surface.
+            //  The lambdas[0] is equivalent for start, lambdas.last()
+            //  is equivalent for end.
+            virtual void findNearest
+            (
+                const point& start,
+                const point& end,
+                const scalarField& lambdas,
+                const scalarField& nearestDistSqr,
+                List<pointIndexHit>&
+            ) const;
+
+            virtual void findLine
+            (
+                const pointField& start,
+                const pointField& end,
+                List<pointIndexHit>&
+            ) const
+            {
+                NotImplemented;
+            }
+
+            virtual void findLineAny
+            (
+                const pointField& start,
+                const pointField& end,
+                List<pointIndexHit>&
+            ) const
+            {
+                NotImplemented;
+            }
+
+            //- Get all intersections in order from start to end.
+            virtual void findLineAll
+            (
+                const pointField& start,
+                const pointField& end,
+                List<List<pointIndexHit>>&
+            ) const
+            {
+                NotImplemented;
+            }
+
+            //- From a set of points and indices get the region
+            virtual void getRegion
+            (
+                const List<pointIndexHit>&,
+                labelList& region
+            ) const;
+
+            //- From a set of points and indices get the normal
+            virtual void getNormal
+            (
+                const List<pointIndexHit>&,
+                vectorField& normal
+            ) const;
+
+            //- Determine type (inside/outside/mixed) for point. unknown if
+            //  cannot be determined (e.g. non-manifold surface)
+            virtual void getVolumeType
+            (
+                const pointField&,
+                List<volumeType>&
+            ) const
+            {
+                NotImplemented;
+            }
+
+
+        // regIOobject implementation
+
+            bool writeData(Ostream&) const
+            {
+                NotImplemented;
+                return false;
+            }
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/blockMesh/pipe/Allrun b/tutorials/mesh/blockMesh/pipe/Allrun
new file mode 100755
index 0000000000..3c4ef3e7ff
--- /dev/null
+++ b/tutorials/mesh/blockMesh/pipe/Allrun
@@ -0,0 +1,9 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+runApplication blockMesh
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/mesh/blockMesh/pipe/constant/triSurface/curve.obj b/tutorials/mesh/blockMesh/pipe/constant/triSurface/curve.obj
new file mode 100644
index 0000000000..9d5e40f6d7
--- /dev/null
+++ b/tutorials/mesh/blockMesh/pipe/constant/triSurface/curve.obj
@@ -0,0 +1,9 @@
+v -5   0.4  0
+v -4   0.4  0
+v -3   0.4  0.7
+v -2.5 0.4  2
+v -2   0.4  2
+v -1   0.4  0.7
+v -0.6 0.4  0
+v 10   0.4  0
+l 1 2 3 4 5 6 7 8
diff --git a/tutorials/mesh/blockMesh/pipe/constant/triSurface/curve2.vtk b/tutorials/mesh/blockMesh/pipe/constant/triSurface/curve2.vtk
new file mode 100644
index 0000000000..02142a3426
--- /dev/null
+++ b/tutorials/mesh/blockMesh/pipe/constant/triSurface/curve2.vtk
@@ -0,0 +1,13 @@
+# vtk DataFile Version 4.0
+vtk output
+ASCII
+DATASET POLYDATA
+POINTS 17 double
+-5.0437 0.4 -0.0159845 -4.54848 0.4 0.00946291 -4.30524 0.4 0.0219617
+-4.11475 0.4 0.0317502 -3.65463 0.4 0.0664504 -3.42509 0.4 0.242198
+-3.26981 0.4 0.570689 -3.04354 0.4 0.986036 -2.80622 0.4 1.28924
+-2.45212 0.4 1.43367 -2.10187 0.4 1.42911 -1.8115 0.4 1.2018
+-1.52708 0.4 0.866397 -1.30229 0.4 0.49514 -1.04633 0.4 0.189424
+-0.5819 0.4 -5.75752e-05 4 0.4 0
+LINES 1 18
+17 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
diff --git a/tutorials/mesh/blockMesh/pipe/system/blockMeshDict b/tutorials/mesh/blockMesh/pipe/system/blockMeshDict
new file mode 100644
index 0000000000..39f45b9230
--- /dev/null
+++ b/tutorials/mesh/blockMesh/pipe/system/blockMeshDict
@@ -0,0 +1,177 @@
+/*--------------------------------*- 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
+{
+    cylinder
+    {
+        type searchableCylinder;
+        point1 (0 -4 0);
+        point2 (0 4 0);
+        radius 0.7;
+    }
+    cylinder3
+    {
+        type searchableCylinder;
+        point1 (-10 0.4 0);
+        point2 (10 0.4 0);
+        radius 0.5;
+    }
+    cylinder2
+    {
+        type    searchableCurve;
+        file    "curve2.vtk";
+        radius  0.5;
+    }
+    inletPlane
+    {
+        type    searchablePlate;
+        origin  (-4 -50 -50);
+        span    (0 100 100);
+    }
+}
+
+vertices
+(
+    // Vertical cylinder
+    named v0 project (-1 -0.1 -1) (cylinder cylinder2)
+    named v1 project ( 1 -0.1 -1) (cylinder)
+    named v2 project ( 1  0.9 -1) (cylinder)
+    named v3 project (-1  0.9 -1) (cylinder cylinder2)
+    named v4 project (-1 -0.1  1) (cylinder cylinder2)
+    named v5 project ( 1 -0.1  1) (cylinder)
+    named v6 project ( 1  0.9  1) (cylinder)
+    named v7 project (-1  0.9  1) (cylinder cylinder2)
+
+
+    // Horizontal cylinder
+    named v8 project (-4   0 -0.5) (cylinder2 inletPlane)
+    named v9 project (-4   1 -0.5) (cylinder2 inletPlane)
+    named v10 project (-4  0  0.5) (cylinder2 inletPlane)
+    named v11 project (-4  1  0.5) (cylinder2 inletPlane)
+
+
+    // On top of vertical cylinder
+    named v12 project (-1 2 -1) (cylinder)
+    named v13 project ( 1 2 -1) (cylinder)
+    named v14 project ( 1 2  1) (cylinder)
+    named v15 project (-1 2  1) (cylinder)
+
+    // Below vertical cylinder
+    named v16 project (-1 -1 -1) (cylinder)
+    named v17 project ( 1 -1 -1) (cylinder)
+    named v18 project ( 1 -1  1) (cylinder)
+    named v19 project (-1 -1  1) (cylinder)
+);
+
+blocks
+(
+    hex (v0 v1 v2 v3 v4 v5 v6 v7) (8 8 8) simpleGrading (1 1 1)
+    named sideBlock hex (v0 v3 v9 v8 v4 v7 v11 v10) (8 20 8)
+        simpleGrading (1 1 1)
+
+    hex ( v7 v6 v2 v3 v15 v14 v13 v12) (8 8 8) simpleGrading (1 1 1)
+    hex (v16 v19 v18 v17 v0 v4 v5 v1) (8 8 8) simpleGrading (1 1 1)
+);
+
+edges
+(
+    project v0 v1 (cylinder)
+    project v1 v2 (cylinder)
+    project v2 v3 (cylinder)
+
+    project v1 v5 (cylinder)
+    project v2 v6 (cylinder)
+
+    project v4 v5 (cylinder)
+    project v5 v6 (cylinder)
+    project v6 v7 (cylinder)
+
+    // Common face
+    project v3 v0 (cylinder cylinder2)
+    project v3 v7 (cylinder cylinder2)
+    project v7 v4 (cylinder cylinder2)
+    project v0 v4 (cylinder cylinder2)
+
+    // Inlet
+    project v8 v10 (cylinder2 inletPlane)
+    project v10 v11 (cylinder2 inletPlane)
+    project v11 v9 (cylinder2 inletPlane)
+    project v9 v8 (cylinder2 inletPlane)
+
+    // Sides of horizontal cylinder. Use projectCurve to do interpolation
+    // for radial direction to keep points along edges at constant radial
+    // direction.
+    projectCurve v8  v0 (cylinder2)
+    projectCurve v9  v3 (cylinder2)
+    projectCurve v11 v7 (cylinder2)
+    projectCurve v10 v4 (cylinder2)
+
+
+
+    // Top cylinder
+    project v12 v15 (cylinder)
+    project v15 v14 (cylinder)
+    project v14 v13 (cylinder)
+    project v13 v12 (cylinder)
+
+    // Bottom cylinder
+    project v16 v17 (cylinder)
+    project v17 v18 (cylinder)
+    project v18 v19 (cylinder)
+    project v19 v16 (cylinder)
+);
+
+faces
+(
+    // Common face
+    project (v0 v4 v7 v3) cylinder
+
+    project (v8 v0 v4 v10) cylinder2
+    project (v10 v4 v7 v11) cylinder2
+    project (v11 v7 v3 v9) cylinder2
+    project (v8 v9 v3 v0) cylinder2
+);
+
+
+defaultPatch
+{
+    name walls;
+    type wall;
+}
+
+boundary
+(
+    side
+    {
+        type    patch;
+        faces   ((sideBlock 3));   //((v8 v10 v11 v9));
+    }
+
+    inlet
+    {
+        type    patch;
+        faces   ((v17 v18 v19 v16));
+    }
+
+    outlet
+    {
+        type    patch;
+        faces   ((v12 v15 v14 v13));
+    }
+);
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/blockMesh/pipe/system/controlDict b/tutorials/mesh/blockMesh/pipe/system/controlDict
new file mode 100644
index 0000000000..728c8b0b1c
--- /dev/null
+++ b/tutorials/mesh/blockMesh/pipe/system/controlDict
@@ -0,0 +1,58 @@
+/*--------------------------------*- 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;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+libs            ("libblockMesh.so");
+
+DebugSwitches
+{
+//    project 1;
+//    searchableCurve 1;
+//    projectCurve 1;
+}
+
+application     blockMesh;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         0;
+
+deltaT          0;
+
+writeControl    timeStep;
+
+writeInterval   1;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable true;
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/blockMesh/pipe/system/fvSchemes b/tutorials/mesh/blockMesh/pipe/system/fvSchemes
new file mode 100644
index 0000000000..606b28ce04
--- /dev/null
+++ b/tutorials/mesh/blockMesh/pipe/system/fvSchemes
@@ -0,0 +1,37 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{}
+
+gradSchemes
+{}
+
+divSchemes
+{}
+
+laplacianSchemes
+{}
+
+interpolationSchemes
+{}
+
+snGradSchemes
+{}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/blockMesh/pipe/system/fvSolution b/tutorials/mesh/blockMesh/pipe/system/fvSolution
new file mode 100644
index 0000000000..6db4e1a6ba
--- /dev/null
+++ b/tutorials/mesh/blockMesh/pipe/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;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// ************************************************************************* //