diff --git a/src/mesh/blockMesh/Make/files b/src/mesh/blockMesh/Make/files
index 123b6510bf..37bbe7ba17 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 74b7eb38b7..3130bcef9a 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 05f41b92ad..fba8517c03 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 038f278ce9..c5d1c53568 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,17 +28,6 @@ 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
diff --git a/src/mesh/blockMesh/block/blockI.H b/src/mesh/blockMesh/block/blockI.H
new file mode 100644
index 0000000000..2d9136b726
--- /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.H b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H
index 323e611608..ac5964ecc5 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 e35e9ee764..dbfa505a7e 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 ff7b1f902b..a7d46b0401 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 2ab190179f..b4cad2a87d 100644
--- a/src/mesh/blockMesh/blockMesh/blockMesh.H
+++ b/src/mesh/blockMesh/blockMesh/blockMesh.H
@@ -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,6 +136,9 @@ 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;
@@ -146,6 +151,11 @@ class 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 9f8e87469e..f8d05d0f2d 100644
--- a/src/mesh/blockMesh/blockMesh/blockMeshMerge.C
+++ b/src/mesh/blockMesh/blockMesh/blockMeshMerge.C
@@ -544,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;
 
diff --git a/src/mesh/blockMesh/blockMesh/blockMeshMergeFast.C b/src/mesh/blockMesh/blockMesh/blockMeshMergeFast.C
new file mode 100644
index 0000000000..d6d2a2d5f4
--- /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;
+}
+
+
+// ************************************************************************* //