diff --git a/src/mesh/blockMesh/Make/files b/src/mesh/blockMesh/Make/files
index 123b6510..37bbe7ba 100644
--- a/src/mesh/blockMesh/Make/files
+++ b/src/mesh/blockMesh/Make/files
@@ -24,5 +24,6 @@ blockMesh/blockMeshCreate.C
 blockMesh/blockMeshTopology.C
 blockMesh/blockMeshCheck.C
 blockMesh/blockMeshMerge.C
+blockMesh/blockMeshMergeFast.C
 
 LIB = $(FOAM_LIBBIN)/libblockMesh
diff --git a/src/mesh/blockMesh/block/block.C b/src/mesh/blockMesh/block/block.C
index 74b7eb38..3130bcef 100644
--- a/src/mesh/blockMesh/block/block.C
+++ b/src/mesh/blockMesh/block/block.C
@@ -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-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/mesh/blockMesh/block/block.H b/src/mesh/blockMesh/block/block.H
index 05f41b92..fba8517c 100644
--- a/src/mesh/blockMesh/block/block.H
+++ b/src/mesh/blockMesh/block/block.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -126,13 +126,10 @@ public:
         // Access
 
             //- Return the block definition
-            inline const blockDescriptor& blockDef() const
-            {
-                return *this;
-            }
+            inline const blockDescriptor& blockDef() const;
 
             //- Vertex label offset for a particular i,j,k position
-            label vtxLabel(label i, label j, label k) const;
+            inline label vtxLabel(label i, label j, label k) const;
 
             //- Return the points for filling the block
             const pointField& points() const;
@@ -161,6 +158,10 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#include "blockI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //
diff --git a/src/mesh/blockMesh/block/blockCreate.C b/src/mesh/blockMesh/block/blockCreate.C
index 038f278c..85c6c583 100644
--- a/src/mesh/blockMesh/block/blockCreate.C
+++ b/src/mesh/blockMesh/block/blockCreate.C
@@ -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-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -28,20 +28,9 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-Foam::label Foam::block::vtxLabel(label i, label j, label k) const
-{
-    return
-    (
-        i
-      + j * (meshDensity().x() + 1)
-      + k * (meshDensity().x() + 1) * (meshDensity().y() + 1)
-    );
-}
-
-
 void Foam::block::createPoints() const
 {
-    // set local variables for mesh specification
+    // Set local variables for mesh specification
     const label ni = meshDensity().x();
     const label nj = meshDensity().y();
     const label nk = meshDensity().z();
@@ -57,12 +46,12 @@ void Foam::block::createPoints() const
     const point& p011 = blockPoint(7);
 
 
-    // list of edge point and weighting factors
+    // List of edge point and weighting factors
     const List< List<point> >& p = blockEdgePoints();
     const scalarListList& w = blockEdgeWeights();
 
     //
-    // generate vertices
+    // Generate vertices
     //
     vertices_.clear();
     vertices_.setSize(nPoints());
@@ -75,7 +64,7 @@ void Foam::block::createPoints() const
             {
                 const label vertexNo = vtxLabel(i, j, k);
 
-                // points on edges
+                // Points on edges
                 vector edgex1 = p000 + (p100 - p000)*w[0][i];
                 vector edgex2 = p010 + (p110 - p010)*w[1][i];
                 vector edgex3 = p011 + (p111 - p011)*w[2][i];
@@ -92,7 +81,7 @@ void Foam::block::createPoints() const
                 vector edgez4 = p010 + (p011 - p010)*w[11][k];
 
 
-                // calculate the importance factors for all edges
+                // Calculate the importance factors for all edges
 
                 // x-direction
                 scalar impx1 =
@@ -190,7 +179,7 @@ void Foam::block::createPoints() const
                 impz4 /= magImpz;
 
 
-                // calculate the correction vectors
+                // Calculate the correction vectors
                 vector corx1 = impx1*(p[0][i] - edgex1);
                 vector corx2 = impx2*(p[1][i] - edgex2);
                 vector corx3 = impx3*(p[2][i] - edgex3);
@@ -207,7 +196,7 @@ void Foam::block::createPoints() const
                 vector corz4 = impz4*(p[11][k] - edgez4);
 
 
-                // multiply by the importance factor
+                // Multiply by the importance factor
 
                 // x-direction
                 edgex1 *= impx1;
@@ -228,7 +217,7 @@ void Foam::block::createPoints() const
                 edgez4 *= impz4;
 
 
-                // add the contributions
+                // Add the contributions
                 vertices_[vertexNo] =
                 (
                     edgex1 + edgex2 + edgex3 + edgex4
@@ -255,7 +244,7 @@ void Foam::block::createCells() const
     const label nk = meshDensity().z();
 
     //
-    // generate cells
+    // Generate cells
     //
     cells_.clear();
     cells_.setSize(nCells());
@@ -292,7 +281,7 @@ void Foam::block::createBoundary() const
     const label nk = meshDensity().z();
 
     //
-    // generate boundaries on each side of the hex
+    // Generate boundaries on each side of the hex
     //
     boundaryPatches_.clear();
     boundaryPatches_.setSize(6);
@@ -311,7 +300,7 @@ void Foam::block::createBoundary() const
         {
             boundaryPatches_[wallLabel][wallCellLabel].setSize(4);
 
-            // set the points
+            // Set the points
             boundaryPatches_[wallLabel][wallCellLabel][0] =
                 vtxLabel(0, j, k);
             boundaryPatches_[wallLabel][wallCellLabel][1] =
@@ -321,7 +310,7 @@ void Foam::block::createBoundary() const
             boundaryPatches_[wallLabel][wallCellLabel][3] =
                 vtxLabel(0, j + 1, k);
 
-            // update the counter
+            // Update the counter
             wallCellLabel++;
         }
     }
@@ -338,7 +327,7 @@ void Foam::block::createBoundary() const
         {
             boundaryPatches_[wallLabel][wallCellLabel].setSize(4);
 
-            // set the points
+            // Set the points
             boundaryPatches_[wallLabel][wallCellLabel][0] =
                 vtxLabel(ni, j, k);
             boundaryPatches_[wallLabel][wallCellLabel][1] =
@@ -348,7 +337,7 @@ void Foam::block::createBoundary() const
             boundaryPatches_[wallLabel][wallCellLabel][3] =
                 vtxLabel(ni, j, k+1);
 
-            // update the counter
+            // Update the counter
             wallCellLabel++;
         }
     }
@@ -366,7 +355,7 @@ void Foam::block::createBoundary() const
         {
             boundaryPatches_[wallLabel][wallCellLabel].setSize(4);
 
-            // set the points
+            // Set the points
             boundaryPatches_[wallLabel][wallCellLabel][0] =
                 vtxLabel(i, 0, k);
             boundaryPatches_[wallLabel][wallCellLabel][1] =
@@ -376,7 +365,7 @@ void Foam::block::createBoundary() const
             boundaryPatches_[wallLabel][wallCellLabel][3] =
                 vtxLabel(i, 0, k + 1);
 
-            // update the counter
+            // Update the counter
             wallCellLabel++;
         }
     }
@@ -393,7 +382,7 @@ void Foam::block::createBoundary() const
         {
             boundaryPatches_[wallLabel][wallCellLabel].setSize(4);
 
-            // set the points
+            // Set the points
             boundaryPatches_[wallLabel][wallCellLabel][0] =
                 vtxLabel(i, nj, k);
             boundaryPatches_[wallLabel][wallCellLabel][1] =
@@ -403,7 +392,7 @@ void Foam::block::createBoundary() const
             boundaryPatches_[wallLabel][wallCellLabel][3] =
                 vtxLabel(i + 1, nj, k);
 
-            // update the counter
+            // Update the counter
             wallCellLabel++;
         }
     }
@@ -422,7 +411,7 @@ void Foam::block::createBoundary() const
         {
             boundaryPatches_[wallLabel][wallCellLabel].setSize(4);
 
-            // set the points
+            // Set the points
             boundaryPatches_[wallLabel][wallCellLabel][0] =
                 vtxLabel(i, j, 0);
             boundaryPatches_[wallLabel][wallCellLabel][1] =
@@ -432,7 +421,7 @@ void Foam::block::createBoundary() const
             boundaryPatches_[wallLabel][wallCellLabel][3] =
                 vtxLabel(i + 1, j, 0);
 
-            // update the counter
+            // Update the counter
             wallCellLabel++;
         }
     }
@@ -449,7 +438,7 @@ void Foam::block::createBoundary() const
         {
             boundaryPatches_[wallLabel][wallCellLabel].setSize(4);
 
-            // set the points
+            // Set the points
             boundaryPatches_[wallLabel][wallCellLabel][0] =
                 vtxLabel(i, j, nk);
             boundaryPatches_[wallLabel][wallCellLabel][1] =
@@ -459,7 +448,7 @@ void Foam::block::createBoundary() const
             boundaryPatches_[wallLabel][wallCellLabel][3] =
                 vtxLabel(i, j + 1, nk);
 
-            // update the counter
+            // Update the counter
             wallCellLabel++;
         }
     }
diff --git a/src/mesh/blockMesh/block/blockI.H b/src/mesh/blockMesh/block/blockI.H
new file mode 100644
index 00000000..2d9136b7
--- /dev/null
+++ b/src/mesh/blockMesh/block/blockI.H
@@ -0,0 +1,47 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2015 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+inline Foam::label Foam::block::vtxLabel(label i, label j, label k) const
+{
+    return
+    (
+        i
+      + j * (meshDensity().x() + 1)
+      + k * (meshDensity().x() + 1) * (meshDensity().y() + 1)
+    );
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const Foam::blockDescriptor& Foam::block::blockDef() const
+{
+    return *this;
+}
+
+
+// ************************************************************************* //
diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C
index bafbc065..fad60ea2 100644
--- a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C
+++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C
@@ -59,7 +59,7 @@ Foam::blockDescriptor::blockDescriptor
             << exit(FatalError);
     }
 
-    // create a list of edges
+    // Create a list of edges
     makeBlockEdges();
 }
 
@@ -99,7 +99,7 @@ Foam::blockDescriptor::blockDescriptor
 
     if (t.isPunctuation())
     {
-        // new-style: read a list of 3 values
+        // New-style: read a list of 3 values
         if (t.pToken() == token::BEGIN_LIST)
         {
             is >> meshDensity_;
@@ -118,7 +118,7 @@ Foam::blockDescriptor::blockDescriptor
     }
     else
     {
-        // old-style: read three labels
+        // Old-style: read three labels
         is  >> meshDensity_.x()
             >> meshDensity_.y()
             >> meshDensity_.z();
@@ -134,7 +134,7 @@ Foam::blockDescriptor::blockDescriptor
 
     if (expRatios.size() == 1)
     {
-        // identical in x/y/z-directions
+        // Identical in x/y/z-directions
         expand_ = expRatios[0];
     }
     else if (expRatios.size() == 3)
@@ -171,7 +171,7 @@ Foam::blockDescriptor::blockDescriptor
             << exit(FatalError);
     }
 
-    // create a list of edges
+    // Create a list of edges
     makeBlockEdges();
 }
 
@@ -279,7 +279,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const blockDescriptor& bd)
 
     const List<gradingDescriptors>& expand = bd.expand_;
 
-    // can we use a compact notation?
+    // Can we use a compact notation?
     if
     (
         // x-direction
diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H
index 323e6116..ac5964ec 100644
--- a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H
+++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H
@@ -142,6 +142,12 @@ public:
 
         // Access
 
+        //- Return the number of cells in the i,j,k directions
+        const Vector<label>& density() const
+        {
+            return meshDensity_;
+        }
+
         //- Reference to point field defining the block mesh
         const pointField& blockPointField() const;
 
diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptorEdges.C b/src/mesh/blockMesh/blockDescriptor/blockDescriptorEdges.C
index e35e9ee7..dbfa505a 100644
--- a/src/mesh/blockMesh/blockDescriptor/blockDescriptorEdges.C
+++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptorEdges.C
@@ -68,7 +68,7 @@ void Foam::blockDescriptor::setEdge
     // Set reference to the list of labels defining the block
     const labelList& blockLabels = blockShape_;
 
-    // Set reference to global list of points
+    // Get list of points for this block
     const pointField blockPoints = blockShape_.points(blockPointField_);
 
     // Set the edge points/weights
diff --git a/src/mesh/blockMesh/blockMesh/blockMesh.C b/src/mesh/blockMesh/blockMesh/blockMesh.C
index ff7b1f90..a7d46b04 100644
--- a/src/mesh/blockMesh/blockMesh/blockMesh.C
+++ b/src/mesh/blockMesh/blockMesh/blockMesh.C
@@ -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-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,11 +24,17 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "blockMesh.H"
+#include "Switch.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 bool Foam::blockMesh::blockMesh::verboseOutput(false);
 
+namespace Foam
+{
+    defineDebugSwitch(blockMesh, 0);
+}
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -38,7 +44,16 @@ Foam::blockMesh::blockMesh(const IOdictionary& dict, const word& regionName)
     scaleFactor_(1.0),
     topologyPtr_(createTopology(dict, regionName))
 {
-    calcMergeInfo();
+    Switch fastMerge(dict.lookupOrDefault<Switch>("fastMerge", false));
+
+    if (fastMerge)
+    {
+         calcMergeInfoFast();
+    }
+    else
+    {
+        calcMergeInfo();
+    }
 }
 
 
diff --git a/src/mesh/blockMesh/blockMesh/blockMesh.H b/src/mesh/blockMesh/blockMesh/blockMesh.H
index e3a6619c..b4cad2a8 100644
--- a/src/mesh/blockMesh/blockMesh/blockMesh.H
+++ b/src/mesh/blockMesh/blockMesh/blockMesh.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -61,6 +61,8 @@ class blockMesh
     public blockList
 {
     // Private data
+
+        //- Switch for verbose output
         static bool verboseOutput;
 
         //- Point field defining the block mesh (corners)
@@ -134,18 +136,26 @@ class blockMesh
         //- Determine the merge info and the final number of cells/points
         void calcMergeInfo();
 
+        //- Determine the merge info and the final number of cells/points
+        void calcMergeInfoFast();
+
         faceList createPatchFaces(const polyPatch& patchTopologyFaces) const;
 
         void createPoints() const;
         void createCells() const;
         void createPatches() const;
 
-        //- as copy (not implemented)
+        //- As copy (not implemented)
         blockMesh(const blockMesh&);
 
 
 public:
 
+    // Static data members
+
+        ClassName("blockMesh");
+
+
     // Constructors
 
         //- Construct from IOdictionary
@@ -164,8 +174,10 @@ public:
             //  these points have not been scaled by scaleFactor
             const pointField& blockPointField() const;
 
+            //- Return the blockMesh topology as a polyMesh
             const polyMesh& topology() const;
 
+            //- Return the curved edges
             const curvedEdgeList& edges() const
             {
                 return edges_;
@@ -178,20 +190,18 @@ public:
             //  these points have been scaled by scaleFactor
             const pointField& points() const;
 
+            //- Return cell shapes list
             const cellShapeList& cells() const;
 
+            //- Return the patch face lists
             const faceListList& patches() const;
 
-
             //- Get patch information from the topology mesh
             PtrList<dictionary> patchDicts() const;
 
+            //- Return patch names
             wordList patchNames() const;
 
-//            wordList patchTypes() const;
-//
-//            wordList patchPhysicalTypes() const;
-
             //- Number of blocks with specified zones
             label numZonedBlocks() const;
 
@@ -204,6 +214,7 @@ public:
             //- Enable/disable verbose information about the progress
             static void verbose(const bool on=true);
 
+
         // Write
 
             //- Writes edges of blockMesh in OBJ format.
@@ -220,4 +231,3 @@ public:
 #endif
 
 // ************************************************************************* //
-
diff --git a/src/mesh/blockMesh/blockMesh/blockMeshMerge.C b/src/mesh/blockMesh/blockMesh/blockMeshMerge.C
index c5263c53..f8d05d0f 100644
--- a/src/mesh/blockMesh/blockMesh/blockMeshMerge.C
+++ b/src/mesh/blockMesh/blockMesh/blockMeshMerge.C
@@ -86,11 +86,7 @@ void Foam::blockMesh::calcMergeInfo()
             blockPfaceLabel++
         )
         {
-            if
-            (
-                blockFaces[blockPfaces[blockPfaceLabel]]
-             == blockFaces[blockFaceLabel]
-            )
+            if (blockPfaces[blockPfaceLabel] == blockFaceLabel)
             {
                 foundFace = true;
                 break;
@@ -345,7 +341,6 @@ void Foam::blockMesh::calcMergeInfo()
                 }
             }
 
-// FIXME? - there seems to be some logic missing here
 
             label blockNfaceLabel;
             for
@@ -365,8 +360,6 @@ void Foam::blockMesh::calcMergeInfo()
                 }
             }
 
-// FIXME? - there seems to be some logic missing here
-
 
             const labelListList& blockPfaceFaces =
                 blocks[blockPlabel].boundaryPatches()[blockPfaceLabel];
@@ -416,7 +409,7 @@ void Foam::blockMesh::calcMergeInfo()
         {
             FatalErrorIn("blockMesh::calcMergeInfo()")
                 << "Point merging failed after max number of passes."
-                << abort(FatalError);
+                << exit(FatalError);
         }
     }
     while (changedPointMerge);
@@ -551,7 +544,7 @@ void Foam::blockMesh::calcMergeInfo()
     }
 
 
-    // sort merge list to return new point label (in new shorter list)
+    // Sort merge list to return new point label (in new shorter list)
     // given old point label
     label newPointLabel = 0;
 
@@ -560,7 +553,8 @@ void Foam::blockMesh::calcMergeInfo()
         if (mergeList_[pointLabel] > pointLabel)
         {
             FatalErrorIn("blockMesh::calcMergeInfo()")
-                << "ouch" << exit(FatalError);
+                << "Merge list contains point index out of range"
+                << exit(FatalError);
         }
 
         if
@@ -579,7 +573,6 @@ void Foam::blockMesh::calcMergeInfo()
     }
 
     nPoints_ = newPointLabel;
-
 }
 
 
diff --git a/src/mesh/blockMesh/blockMesh/blockMeshMergeFast.C b/src/mesh/blockMesh/blockMesh/blockMeshMergeFast.C
new file mode 100644
index 00000000..d6d2a2d5
--- /dev/null
+++ b/src/mesh/blockMesh/blockMesh/blockMeshMergeFast.C
@@ -0,0 +1,590 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 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 "blockMesh.H"
+
+// * * * * * * * * * * * * * * * Private Functions * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// faces
+// 6
+// (
+//     4(0 4 7 3)   // x-min
+//     4(1 2 6 5)   // x-max
+//     4(0 1 5 4)   // y-min
+//     4(3 7 6 2)   // y-max
+//     4(0 3 2 1)   // z-min
+//     4(4 5 6 7)   // z-max
+// );
+
+// Face-edge directions
+static const int faceEdgeDirs[6][4] =
+{
+    {2, 1, -2, -1},
+    {1, 2, -1, -2},
+    {1, 2, -1, -2},
+    {2, 1, -2, -1},
+    {2, 1, -2, -1},
+    {1, 2, -1, -2}
+};
+
+// The face-face-rotation direction correspondence map
+static Pair<int> faceFaceRotMap[6][6][4];
+
+// Generate the face-face-rotation direction correspondence map
+void genFaceFaceRotMap()
+{
+    for(int facePi=0; facePi<6; facePi++)
+    {
+        for(int faceNi=0; faceNi<6; faceNi++)
+        {
+            for(int rot=0; rot<4; rot++)
+            {
+                Pair<int>& map = faceFaceRotMap[facePi][faceNi][rot];
+
+                for(int Pp=0; Pp<2; Pp++)
+                {
+                    int Pdir = faceEdgeDirs[facePi][Pp];
+                    int Np = (3 - Pp + rot)%4;
+                    int Ndir = faceEdgeDirs[faceNi][Np];
+                    map[Pdir-1] = -Ndir;
+                }
+
+                // Handle sign change due to match-face transpose
+                if (mag(map[0]) == 2 && map[0]*map[1] < 0)
+                {
+                    map[0] = -map[0];
+                    map[1] = -map[1];
+                }
+            }
+        }
+    }
+}
+
+// Return the direction map for the merge-faces
+Pair<int> faceMap
+(
+    const label facePi,
+    const face& faceP,
+    const label faceNi,
+    const face& faceN
+)
+{
+    // Search for the point on faceN corresponding to the 0-point on faceP
+    for(int rot=0; rot<4; rot++)
+    {
+        if (faceN[rot] == faceP[0])
+        {
+            return faceFaceRotMap[facePi][faceNi][rot];
+        }
+    }
+
+    FatalErrorIn
+    (
+        "faceMap(const label facePi, const face& faceP, "
+        "const label faceNi, const face& faceN)"
+    )   << "Cannot find point correspondance for faces "
+        << faceP << " and " << faceN
+        << exit(FatalError);
+
+    return Pair<int>(0, 0);
+}
+
+// Set the block and face indices for all the merge faces
+void setBlockFaceCorrespondence
+(
+    const cellList& topoCells,
+    const faceList::subList& topoInternalFaces,
+    const labelList& topoFaceCell,
+    List<Pair<label> >& mergeBlock
+)
+{
+    forAll(topoInternalFaces, topoFacei)
+    {
+        label topoPi = topoFaceCell[topoFacei];
+        const labelList& topoPfaces = topoCells[topoPi];
+
+        bool foundFace = false;
+        label topoPfacei;
+        for
+        (
+            topoPfacei = 0;
+            topoPfacei < topoPfaces.size();
+            topoPfacei++
+        )
+        {
+            if (topoPfaces[topoPfacei] == topoFacei)
+            {
+                foundFace = true;
+                break;
+            }
+        }
+
+        if (!foundFace)
+        {
+            FatalErrorIn("setBlockFaceCorrespondence()")
+                << "Cannot find merge face for block " << topoPi
+                << exit(FatalError);
+        }
+
+        mergeBlock[topoFacei].first() = topoPi;
+        mergeBlock[topoFacei].second() = topoPfacei;
+    }
+}
+
+// Return the number of divisions in each direction for the face
+Pair<label> faceNij(const label facei, const block& block)
+{
+    Pair<label> fnij;
+
+    int i = facei/2;
+
+    if (i == 0)
+    {
+        fnij.first() = block.meshDensity().y() + 1;
+        fnij.second() = block.meshDensity().z() + 1;
+    }
+    else if (i == 1)
+    {
+        fnij.first() = block.meshDensity().x() + 1;
+        fnij.second() = block.meshDensity().z() + 1;
+    }
+    else if (i == 2)
+    {
+        fnij.first() = block.meshDensity().x() + 1;
+        fnij.second() = block.meshDensity().y() + 1;
+    }
+
+    return fnij;
+}
+
+// Sign the index corresponding to the map
+inline label signIndex(const int map, const label i)
+{
+    return map < 0 ? -i-1 : i;
+}
+
+// Reverse a signed index with the number of divisions
+inline label unsignIndex(const label i, const label ni)
+{
+    return i >= 0 ? i : ni + i + 1;
+}
+
+// Return the mapped index
+inline label mapij(const int map, const label i, const label j)
+{
+    return signIndex(map, mag(map) == 1 ? i : j);
+}
+
+// Return the face point index
+inline label facePoint
+(
+    const int facei,
+    const block& block,
+    const label i,
+    const label j
+)
+{
+    switch (facei)
+    {
+        case 0:
+            return block.vtxLabel(0, i, j);
+        case 1:
+            return block.vtxLabel(block.meshDensity().x(), i, j);
+        case 2:
+            return block.vtxLabel(i, 0, j);
+        case 3:
+            return block.vtxLabel(i, block.meshDensity().y(), j);
+        case 4:
+            return block.vtxLabel(i, j, 0);
+        case 5:
+            return block.vtxLabel(i, j, block.meshDensity().z());
+        default:
+            return -1;
+    }
+}
+
+// Return the neighbour face point from the signed indices
+inline label facePointN
+(
+    const block& block,
+    const label i,
+    const label j,
+    const label k
+)
+{
+    return block.vtxLabel
+    (
+        unsignIndex(i, block.meshDensity().x()),
+        unsignIndex(j, block.meshDensity().y()),
+        unsignIndex(k, block.meshDensity().z())
+    );
+}
+
+// Return the neighbour face point from the mapped indices
+inline label facePointN
+(
+    const int facei,
+    const block& block,
+    const label i,
+    const label j
+)
+{
+    switch (facei)
+    {
+        case 0:
+            return facePointN(block, 0, i, j);
+        case 1:
+            return facePointN(block, block.meshDensity().x(), i, j);
+        case 2:
+            return facePointN(block, i, 0, j);
+        case 3:
+            return facePointN(block, i, block.meshDensity().y(), j);
+        case 4:
+            return facePointN(block, i, j, 0);
+        case 5:
+            return facePointN(block, i, j, block.meshDensity().z());
+        default:
+            return -1;
+    }
+}
+
+// Return the neighbour face point using the map
+inline label facePointN
+(
+    const int facei,
+    const Pair<int>& fmap,
+    const block& block,
+    const label i,
+    const label j
+)
+{
+    return facePointN(facei, block, mapij(fmap[0], i, j), mapij(fmap[1], i, j));
+}
+
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::blockMesh::calcMergeInfoFast()
+{
+    // Generate the static face-face map
+    genFaceFaceRotMap();
+
+    const blockList& blocks = *this;
+
+    if (verboseOutput)
+    {
+        Info<< "Creating block offsets" << endl;
+    }
+
+    blockOffsets_.setSize(blocks.size());
+
+    nPoints_ = 0;
+    nCells_  = 0;
+
+    forAll(blocks, blockI)
+    {
+        blockOffsets_[blockI] = nPoints_;
+
+        nPoints_ += blocks[blockI].nPoints();
+        nCells_  += blocks[blockI].nCells();
+    }
+
+    if (verboseOutput)
+    {
+        Info<< "Creating merge list using the fast topological search"
+            << flush;
+    }
+
+    // Size merge list and initialize to -1
+    mergeList_.setSize(nPoints_, -1);
+
+    // Block mesh topology
+    const pointField& topoPoints = topology().points();
+    const cellList& topoCells = topology().cells();
+    const faceList& topoFaces = topology().faces();
+    const labelList& topoFaceOwn = topology().faceOwner();
+    const labelList& topoFaceNei = topology().faceNeighbour();
+
+    // Topological merging is only necessary for internal block faces
+    // Note edge and face collapse may apply to boundary faces
+    // but is not yet supported in the "fast" algorithm
+    const faceList::subList topoInternalFaces
+    (
+        topoFaces,
+        topology().nInternalFaces()
+    );
+
+    List<Pair<label> > mergeBlockP(topoInternalFaces.size());
+    setBlockFaceCorrespondence
+    (
+        topoCells,
+        topoInternalFaces,
+        topoFaceOwn,
+        mergeBlockP
+    );
+
+    List<Pair<label> > mergeBlockN(topoInternalFaces.size());
+    setBlockFaceCorrespondence
+    (
+        topoCells,
+        topoInternalFaces,
+        topoFaceNei,
+        mergeBlockN
+    );
+
+    if (debug)
+    {
+        Info<< endl;
+    }
+
+    forAll(topoInternalFaces, topoFacei)
+    {
+        if (debug)
+        {
+            Info<< "Processing face " << topoFacei << endl;
+        }
+
+        label blockPi = mergeBlockP[topoFacei].first();
+        label blockPfacei = mergeBlockP[topoFacei].second();
+
+        label blockNi = mergeBlockN[topoFacei].first();
+        label blockNfacei = mergeBlockN[topoFacei].second();
+
+        Pair<int> fmap
+        (
+            faceMap
+            (
+                blockPfacei,
+                blocks[blockPi].blockShape().faces()[blockPfacei],
+                blockNfacei,
+                blocks[blockNi].blockShape().faces()[blockNfacei]
+            )
+        );
+
+        if (debug)
+        {
+            Info<< "    Face map for faces "
+                << blocks[blockPi].blockShape().faces()[blockPfacei] << " "
+                << blocks[blockNi].blockShape().faces()[blockNfacei] << ": "
+                << fmap << endl;
+        }
+
+        const pointField& blockPpoints = blocks[blockPi].points();
+        const pointField& blockNpoints = blocks[blockNi].points();
+
+        Pair<label> Pnij(faceNij(blockPfacei, blocks[blockPi]));
+
+        // Check block subdivision correspondence
+        {
+            Pair<label> Nnij(faceNij(blockNfacei, blocks[blockNi]));
+            Pair<label> NPnij;
+            NPnij[0] = Nnij[mag(fmap[0]) - 1];
+            NPnij[1] = Nnij[mag(fmap[1]) - 1];
+
+            if (Pnij != NPnij)
+            {
+                FatalErrorIn("blockMesh::calcMergeInfoFast()")
+                    << "Sub-division mismatch between face "
+                    << blockPfacei << " of block " << blockPi << Pnij
+                    << " and face "
+                    << blockNfacei << " of block " << blockNi << Nnij
+                    << exit(FatalError);
+            }
+        }
+
+        // Calculate a suitable test distance from the bounding box of the face.
+        // Note this is used only as a sanity check and for diagnostics in
+        // case there is a grading inconsistency.
+        const boundBox bb(topoCells[blockPi].points(topoFaces, topoPoints));
+        const scalar testSqrDist = magSqr(1e-6*bb.span());
+
+        // Accumulate the maximum merge distance for diagnostics
+        scalar maxSqrDist = 0;
+
+        for (label j=0; j<Pnij.second(); j++)
+        {
+            for (label i=0; i<Pnij.first(); i++)
+            {
+                label blockPpointi =
+                    facePoint(blockPfacei, blocks[blockPi], i, j);
+
+                label blockNpointi =
+                    facePointN(blockNfacei, fmap, blocks[blockNi], i, j);
+
+                scalar sqrDist
+                (
+                    magSqr
+                    (
+                        blockPpoints[blockPpointi]
+                      - blockNpoints[blockNpointi]
+                    )
+                );
+
+                if (sqrDist > testSqrDist)
+                {
+                    FatalErrorIn("blockMesh::calcMergeInfoFast()")
+                        << "Point merge failure between face "
+                        << blockPfacei << " of block " << blockPi
+                        << " and face "
+                        << blockNfacei << " of block " << blockNi
+                        << endl
+                        << "    This may be due to inconsistent grading."
+                        << exit(FatalError);
+                }
+
+                maxSqrDist = max(maxSqrDist, sqrDist);
+
+                label Ppointi = blockPpointi + blockOffsets_[blockPi];
+                label Npointi = blockNpointi + blockOffsets_[blockNi];
+
+                label minPNi = min(Ppointi, Npointi);
+
+                if (mergeList_[Ppointi] != -1)
+                {
+                    minPNi = min(minPNi, mergeList_[Ppointi]);
+                }
+
+                if (mergeList_[Npointi] != -1)
+                {
+                    minPNi = min(minPNi, mergeList_[Npointi]);
+                }
+
+                mergeList_[Ppointi] = mergeList_[Npointi] = minPNi;
+            }
+        }
+
+        if (debug)
+        {
+            Info<< "    Max distance between merge points: "
+                << sqrt(maxSqrDist) << endl;
+        }
+    }
+
+
+    bool changedPointMerge = false;
+    label nPasses = 0;
+
+    do
+    {
+        changedPointMerge = false;
+        nPasses++;
+
+        forAll(topoInternalFaces, topoFacei)
+        {
+            label blockPi = mergeBlockP[topoFacei].first();
+            label blockPfacei = mergeBlockP[topoFacei].second();
+
+            label blockNi = mergeBlockN[topoFacei].first();
+            label blockNfacei = mergeBlockN[topoFacei].second();
+
+            Pair<int> fmap
+            (
+                faceMap
+                (
+                    blockPfacei,
+                    blocks[blockPi].blockShape().faces()[blockPfacei],
+                    blockNfacei,
+                    blocks[blockNi].blockShape().faces()[blockNfacei]
+                )
+            );
+
+            Pair<label> Pnij(faceNij(blockPfacei, blocks[blockPi]));
+
+            for (label j=0; j<Pnij.second(); j++)
+            {
+                for (label i=0; i<Pnij.first(); i++)
+                {
+                    label blockPpointi =
+                        facePoint(blockPfacei, blocks[blockPi], i, j);
+
+                    label blockNpointi =
+                        facePointN(blockNfacei, fmap, blocks[blockNi], i, j);
+
+                    label Ppointi =
+                        blockPpointi + blockOffsets_[blockPi];
+
+                    label Npointi =
+                        blockNpointi + blockOffsets_[blockNi];
+
+                    if (mergeList_[Ppointi] != mergeList_[Npointi])
+                    {
+                        changedPointMerge = true;
+
+                        mergeList_[Ppointi]
+                          = mergeList_[Npointi]
+                          = min(mergeList_[Ppointi], mergeList_[Npointi]);
+                    }
+                }
+            }
+        }
+
+        if (verboseOutput)
+        {
+            Info<< "." << flush;
+        }
+
+        if (nPasses > 100)
+        {
+            FatalErrorIn("blockMesh::calcMergeInfoFast()")
+                << "Point merging failed after 100 passes."
+                << exit(FatalError);
+        }
+
+    } while (changedPointMerge);
+
+
+    // Sort merge list and count number of unique points
+    label nUniqPoints = 0;
+
+    forAll(mergeList_, pointi)
+    {
+        if (mergeList_[pointi] > pointi)
+        {
+            FatalErrorIn("blockMesh::calcMergeInfoFast()")
+                << "Merge list contains point index out of range"
+                << exit(FatalError);
+        }
+
+        if (mergeList_[pointi] == -1 || mergeList_[pointi] == pointi)
+        {
+            mergeList_[pointi] = nUniqPoints++;
+        }
+        else
+        {
+            mergeList_[pointi] = mergeList_[mergeList_[pointi]];
+        }
+    }
+
+    // Correct number of points in mesh
+    nPoints_ = nUniqPoints;
+}
+
+
+// ************************************************************************* //
diff --git a/src/mesh/blockMesh/curvedEdges/curvedEdge.C b/src/mesh/blockMesh/curvedEdges/curvedEdge.C
index a764cc85..38a6637b 100644
--- a/src/mesh/blockMesh/curvedEdges/curvedEdge.C
+++ b/src/mesh/blockMesh/curvedEdges/curvedEdge.C
@@ -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-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -117,11 +117,11 @@ Foam::pointField Foam::curvedEdge::appendEndPoints
 {
     pointField allKnots(otherKnots.size() + 2);
 
-    // start/end knots
+    // Start/end knots
     allKnots[0] = points[start];
     allKnots[otherKnots.size() + 1] = points[end];
 
-    // intermediate knots
+    // Intermediate knots
     forAll(otherKnots, knotI)
     {
         allKnots[knotI+1] = otherKnots[knotI];
diff --git a/src/mesh/blockMesh/curvedEdges/lineDivide.C b/src/mesh/blockMesh/curvedEdges/lineDivide.C
index 7ac372c4..5169b9d2 100644
--- a/src/mesh/blockMesh/curvedEdges/lineDivide.C
+++ b/src/mesh/blockMesh/curvedEdges/lineDivide.C
@@ -59,17 +59,38 @@ Foam::lineDivide::lineDivide
     // Check that there are more divisions than sections
     if (nDiv >= gd.size())
     {
+        // Calculate distribution of divisions to be independent
+        // of the order of the sections
+        labelList secnDivs(gd.size());
+        label sumSecnDivs = 0;
+        label secnMaxDivs = 0;
+
+        forAll(gd, sectioni)
+        {
+            scalar nDivFrac = gd[sectioni].nDivFraction();
+            secnDivs[sectioni] = label(nDivFrac*nDiv + 0.5);
+            sumSecnDivs += secnDivs[sectioni];
+
+            // Find the section with the largest number of divisions
+            if (nDivFrac > gd[secnMaxDivs].nDivFraction())
+            {
+                secnMaxDivs = sectioni;
+            }
+        }
+
+        // Adjust the number of divisions on the section with the largest
+        // number so that the total is nDiv
+        if (sumSecnDivs != nDiv)
+        {
+            secnDivs[secnMaxDivs] += (nDiv - sumSecnDivs);
+        }
+
         forAll(gd, sectioni)
         {
             scalar blockFrac = gd[sectioni].blockFraction();
-            scalar nDivFrac = gd[sectioni].nDivFraction();
             scalar expRatio = gd[sectioni].expansionRatio();
 
-            label secnDiv = label(nDivFrac*nDiv + 0.5);
-            if (sectioni == gd.size() - 1)
-            {
-                secnDiv = nDiv - secnStart + 1;
-            }
+            label secnDiv = secnDivs[sectioni];
             label secnEnd = secnStart + secnDiv;
 
             // Calculate the spacing
@@ -109,7 +130,6 @@ Foam::lineDivide::lineDivide
         }
     }
 
-
     // Calculate the points
     for (label i = 0; i <= nDiv; i++)
     {
diff --git a/src/mesh/blockMesh/curvedEdges/polyLine.C b/src/mesh/blockMesh/curvedEdges/polyLine.C
index a3c05d34..1d0e39d4 100644
--- a/src/mesh/blockMesh/curvedEdges/polyLine.C
+++ b/src/mesh/blockMesh/curvedEdges/polyLine.C
@@ -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-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,7 +23,6 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "error.H"
 #include "polyLine.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
@@ -41,7 +40,7 @@ void Foam::polyLine::calcParam()
             param_[i] = param_[i-1] + mag(points_[i] - points_[i-1]);
         }
 
-        // normalize on the interval 0-1
+        // Normalize on the interval 0-1
         lineLength_ = param_.last();
         for (label i=1; i < param_.size() - 1; i++)
         {
@@ -85,7 +84,7 @@ Foam::label Foam::polyLine::nSegments() const
 
 Foam::label Foam::polyLine::localParameter(scalar& lambda) const
 {
-    // check endpoints
+    // Check endpoints
     if (lambda < SMALL)
     {
         lambda = 0;
@@ -97,22 +96,20 @@ Foam::label Foam::polyLine::localParameter(scalar& lambda) const
         return nSegments();
     }
 
-    // search table of cumulative distances to find which line-segment
-    // we are on. Check the upper bound.
+    // Search table of cumulative distances to find which line-segment
+    // we are on.
+    // Check the upper bound.
 
     label segmentI = 1;
     while (param_[segmentI] < lambda)
     {
         segmentI++;
     }
-    segmentI--;   // we want the corresponding lower bound
+    segmentI--;   // We want the corresponding lower bound
 
-    // the local parameter [0-1] on this line segment
+    // The local parameter [0-1] on this line segment
     lambda =
-    (
-        ( lambda - param_[segmentI] )
-      / ( param_[segmentI+1] - param_[segmentI] )
-    );
+        (lambda - param_[segmentI])/(param_[segmentI+1] - param_[segmentI]);
 
     return segmentI;
 }
@@ -120,7 +117,7 @@ Foam::label Foam::polyLine::localParameter(scalar& lambda) const
 
 Foam::point Foam::polyLine::position(const scalar mu) const
 {
-    // check endpoints
+    // Check end-points
     if (mu < SMALL)
     {
         return points_.first();
@@ -130,7 +127,6 @@ Foam::point Foam::polyLine::position(const scalar mu) const
         return points_.last();
     }
 
-
     scalar lambda = mu;
     label segment = localParameter(lambda);
     return position(segment, lambda);
@@ -143,7 +139,7 @@ Foam::point Foam::polyLine::position
     const scalar mu
 ) const
 {
-    // out-of-bounds
+    // Out-of-bounds
     if (segment < 0)
     {
         return points_.first();
@@ -156,7 +152,7 @@ Foam::point Foam::polyLine::position
     const point& p0 = points()[segment];
     const point& p1 = points()[segment+1];
 
-    // special cases - no calculation needed
+    // Special cases - no calculation needed
     if (mu <= 0.0)
     {
         return p0;
@@ -167,8 +163,8 @@ Foam::point Foam::polyLine::position
     }
     else
     {
-        // linear interpolation
-        return points_[segment] + mu * (p1 - p0);
+        // Linear interpolation
+        return points_[segment] + mu*(p1 - p0);
     }
 }
 
diff --git a/src/mesh/blockMesh/gradingDescriptor/gradingDescriptors.C b/src/mesh/blockMesh/gradingDescriptor/gradingDescriptors.C
index 5284deac..6e6418af 100644
--- a/src/mesh/blockMesh/gradingDescriptor/gradingDescriptors.C
+++ b/src/mesh/blockMesh/gradingDescriptor/gradingDescriptors.C
@@ -47,7 +47,7 @@ Foam::gradingDescriptors Foam::gradingDescriptors::inv() const
 
     forAll(ret, i)
     {
-        ret[i] = operator[](i).inv();
+        ret[i] = operator[](ret.size() - i - 1).inv();
     }
 
     return ret;