From a8f61f210f984c4f1e591f9d1eb975d704658757 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Sun, 12 Jul 2015 20:32:25 +0100
Subject: [PATCH] blockMesh: added experimental fast-merge algorithm

The standard merge-algorithm is N^2 over the face-points and uses a
geometric proximity test for the merge.  These are both choices for
implementation simplicity and are rather inefficient for large meshes.
I have now implemented an experimental linear topological merge
algorithm which is VERY fast and effective for meshes of any size.
Currently it will merge internal faces on meshes of arbitrary complexity
but does not yet handle edge or face collapse needed for wedges and
other degenerate blocks.

The new fast-merge algorithm may be selected using the optional
"fastMerge" entry:

fastMerge yes;

and if not present the standard N^2 algorithm will be used.

Henry G. Weller
CFD Direct
---
 src/mesh/blockMesh/Make/files                 |   1 +
 src/mesh/blockMesh/block/block.C              |   2 +-
 src/mesh/blockMesh/block/block.H              |  13 +-
 src/mesh/blockMesh/block/blockCreate.C        |  13 +-
 src/mesh/blockMesh/block/blockI.H             |  47 ++
 .../blockDescriptor/blockDescriptor.H         |   6 +
 .../blockDescriptor/blockDescriptorEdges.C    |   2 +-
 src/mesh/blockMesh/blockMesh/blockMesh.C      |  19 +-
 src/mesh/blockMesh/blockMesh/blockMesh.H      |  22 +-
 src/mesh/blockMesh/blockMesh/blockMeshMerge.C |   2 +-
 .../blockMesh/blockMesh/blockMeshMergeFast.C  | 590 ++++++++++++++++++
 11 files changed, 688 insertions(+), 29 deletions(-)
 create mode 100644 src/mesh/blockMesh/block/blockI.H
 create mode 100644 src/mesh/blockMesh/blockMesh/blockMeshMergeFast.C

diff --git a/src/mesh/blockMesh/Make/files b/src/mesh/blockMesh/Make/files
index 123b6510b..37bbe7ba1 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 74b7eb38b..3130bcef9 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 05f41b92a..fba8517c0 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 038f278ce..c5d1c5356 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 000000000..2d9136b72
--- /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 323e61160..ac5964ecc 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 e35e9ee76..dbfa505a7 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 ff7b1f902..a7d46b040 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 2ab190179..b4cad2a87 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 9f8e87469..f8d05d0f2 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 000000000..d6d2a2d5f
--- /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;
+}
+
+
+// ************************************************************************* //