diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
index 22b905f42e..b9fedffe61 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2020 OpenCFD Ltd.
+    Copyright (C) 2015-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -1658,7 +1658,7 @@ int main(int argc, char *argv[])
     }
 
     // Mesh distribution engine (uses tolerance to reconstruct meshes)
-    fvMeshDistribute distributor(mesh, mergeDist);
+    fvMeshDistribute distributor(mesh);
 
 
 
diff --git a/applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C b/applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C
index e4f717f682..25ddd33ca1 100644
--- a/applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C
+++ b/applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -36,11 +36,25 @@ Description
     Writes point/face/cell procAddressing so afterwards reconstructPar can be
     used to reconstruct fields.
 
-    Note:
-    - uses geometric matching tolerance (set with -mergeTol (at your option)
+Usage
+    \b reconstructParMesh [OPTION]
 
-    If the parallel case does not have correct procBoundaries use the
-    -fullMatch option which will check all boundary faces (bit slower).
+    Options:
+      - \par -fullMatch
+        Does geometric matching on all boundary faces. Assumes no point
+        ordering
+
+      - \par -procMatch
+        Assumes processor patches already in face order but not point order.
+        This is the pre v2106 default behaviour but might be removed if the new
+        topological method works well
+
+      - \par -mergeTol \<tol\>
+        Specifies non-default merge tolerance (fraction of mesh bounding box)
+        for above options
+
+    The default is to assume all processor boundaries are correctly ordered
+    (both faces and points) in which case no merge tolerance is needed.
 
 \*---------------------------------------------------------------------------*/
 
@@ -57,6 +71,7 @@ Description
 #include "polyTopoChange.H"
 #include "extrapolatedCalculatedFvPatchFields.H"
 #include "topoSet.H"
+#include "fvMeshTools.H"
 
 using namespace Foam;
 
@@ -67,22 +82,6 @@ using namespace Foam;
 static const scalar defaultMergeTol = 1e-7;
 
 
-static void renumber
-(
-    const labelList& map,
-    labelList& elems
-)
-{
-    forAll(elems, i)
-    {
-        if (elems[i] >= 0)
-        {
-            elems[i] = map[elems[i]];
-        }
-    }
-}
-
-
 // Determine which faces are coupled. Uses geometric merge distance.
 // Looks either at all boundaryFaces (fullMatch) or only at the
 // procBoundaries for proci. Assumes that masterMesh contains already merged
@@ -366,12 +365,11 @@ boundBox procBounds
 }
 
 
-void writeCellDistance
+void writeDistribution
 (
     Time& runTime,
     const fvMesh& masterMesh,
     const labelListList& cellProcAddressing
-
 )
 {
     // Write the decomposition as labelList for use with 'manual'
@@ -446,6 +444,175 @@ void writeCellDistance
 }
 
 
+void writeMesh
+(
+    const fvMesh& mesh,
+    const labelListList& cellProcAddressing
+)
+{
+    Info<< "\nWriting merged mesh to "
+        << mesh.time().path()/mesh.time().timeName()
+        << nl << endl;
+
+    if (!mesh.write())
+    {
+        FatalErrorInFunction
+            << "Failed writing polyMesh."
+            << exit(FatalError);
+    }
+    topoSet::removeFiles(mesh);
+}
+
+
+void writeMaps
+(
+    const label masterInternalFaces,
+    const labelUList& masterOwner,
+    const polyMesh& procMesh,
+    const labelUList& cellProcAddressing,
+    const labelUList& faceProcAddressing,
+    const labelUList& pointProcAddressing,
+    const labelUList& boundaryProcAddressing
+)
+{
+    // From processor point to reconstructed mesh point
+
+    Info<< "Writing pointProcAddressing to "
+        << procMesh.time().caseName()
+          /procMesh.facesInstance()
+          /polyMesh::meshSubDir
+        << endl;
+
+    labelIOList
+    (
+        IOobject
+        (
+            "pointProcAddressing",
+            procMesh.facesInstance(),
+            polyMesh::meshSubDir,
+            procMesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false                       // Do not register
+        ),
+        pointProcAddressing
+    ).write();
+
+
+    // From processor face to reconstructed mesh face
+
+    Info<< "Writing faceProcAddressing to "
+        << procMesh.time().caseName()
+          /procMesh.facesInstance()
+          /polyMesh::meshSubDir
+        << endl;
+
+    labelIOList faceProcAddr
+    (
+        IOobject
+        (
+            "faceProcAddressing",
+            procMesh.facesInstance(),
+            polyMesh::meshSubDir,
+            procMesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false                       // Do not register
+        ),
+        faceProcAddressing
+    );
+
+    // Now add turning index to faceProcAddressing.
+    // See reconstructPar for meaning of turning index.
+    forAll(faceProcAddr, procFacei)
+    {
+        const label masterFacei = faceProcAddr[procFacei];
+
+        if
+        (
+           !procMesh.isInternalFace(procFacei)
+         && masterFacei < masterInternalFaces
+        )
+        {
+            // proc face is now external but used to be internal face.
+            // Check if we have owner or neighbour.
+
+            label procOwn = procMesh.faceOwner()[procFacei];
+            label masterOwn = masterOwner[masterFacei];
+
+            if (cellProcAddressing[procOwn] == masterOwn)
+            {
+                // No turning. Offset by 1.
+                faceProcAddr[procFacei]++;
+            }
+            else
+            {
+                // Turned face.
+                faceProcAddr[procFacei] = -1 - faceProcAddr[procFacei];
+            }
+        }
+        else
+        {
+            // No turning. Offset by 1.
+            faceProcAddr[procFacei]++;
+        }
+    }
+
+    faceProcAddr.write();
+
+
+    // From processor cell to reconstructed mesh cell
+
+    Info<< "Writing cellProcAddressing to "
+        << procMesh.time().caseName()
+          /procMesh.facesInstance()
+          /polyMesh::meshSubDir
+        << endl;
+
+    labelIOList
+    (
+        IOobject
+        (
+            "cellProcAddressing",
+            procMesh.facesInstance(),
+            polyMesh::meshSubDir,
+            procMesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false                       // Do not register
+        ),
+        cellProcAddressing
+    ).write();
+
+
+
+    // From processor patch to reconstructed mesh patch
+
+    Info<< "Writing boundaryProcAddressing to "
+        << procMesh.time().caseName()
+          /procMesh.facesInstance()
+          /polyMesh::meshSubDir
+        << endl;
+
+    labelIOList
+    (
+        IOobject
+        (
+            "boundaryProcAddressing",
+            procMesh.facesInstance(),
+            polyMesh::meshSubDir,
+            procMesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false                       // Do not register
+        ),
+        boundaryProcAddressing
+    ).write();
+
+    Info<< endl;
+}
+
+
 int main(int argc, char *argv[])
 {
     argList::addNote
@@ -470,6 +637,11 @@ int main(int argc, char *argv[])
         "Do (slower) geometric matching on all boundary faces"
     );
     argList::addBoolOption
+    (
+        "procMatch",
+        "Do matching on processor faces only"
+    );
+    argList::addBoolOption
     (
         "cellDist",
         "Write cell distribution as a labelList - for use with 'manual' "
@@ -509,40 +681,46 @@ int main(int argc, char *argv[])
         Info<< "Operating on region " << regionName << nl << endl;
     }
 
-    const scalar mergeTol =
-        args.getOrDefault<scalar>("mergeTol", defaultMergeTol);
-
-    scalar writeTol = Foam::pow(10.0, -scalar(IOstream::defaultPrecision()));
-
-    Info<< "Merge tolerance : " << mergeTol << nl
-        << "Write tolerance : " << writeTol << endl;
-
-    if (runTime.writeFormat() == IOstream::ASCII && mergeTol < writeTol)
-    {
-        FatalErrorInFunction
-            << "Your current settings specify ASCII writing with "
-            << IOstream::defaultPrecision() << " digits precision." << endl
-            << "Your merging tolerance (" << mergeTol << ") is finer than this."
-            << endl
-            << "Please change your writeFormat to binary"
-            << " or increase the writePrecision" << endl
-            << "or adjust the merge tolerance (-mergeTol)."
-            << exit(FatalError);
-    }
-
-
     const bool fullMatch = args.found("fullMatch");
+    const bool procMatch = args.found("procMatch");
     const bool writeCellDist = args.found("cellDist");
 
     if (fullMatch)
     {
         Info<< "Doing geometric matching on all boundary faces." << nl << endl;
     }
-    else
+    else if (procMatch)
     {
         Info<< "Doing geometric matching on correct procBoundaries only."
             << nl << "This assumes a correct decomposition." << endl;
     }
+    else
+    {
+        Info<< "Assuming correct, fully matched procBoundaries." << nl << endl;
+    }
+
+    scalar mergeTol = args.getOrDefault<scalar>("mergeTol", defaultMergeTol);
+    if (fullMatch || procMatch)
+    {
+        scalar writeTol =
+            Foam::pow(10.0, -scalar(IOstream::defaultPrecision()));
+
+        Info<< "Merge tolerance : " << mergeTol << nl
+            << "Write tolerance : " << writeTol << endl;
+
+        if (runTime.writeFormat() == IOstream::ASCII && mergeTol < writeTol)
+        {
+            FatalErrorInFunction
+                << "Your current settings specify ASCII writing with "
+                << IOstream::defaultPrecision() << " digits precision." << endl
+                << "Your merging tolerance (" << mergeTol << ")"
+                << " is finer than this." << endl
+                << "Please change your writeFormat to binary"
+                << " or increase the writePrecision" << endl
+                << "or adjust the merge tolerance (-mergeTol)."
+                << exit(FatalError);
+        }
+    }
 
     label nProcs = fileHandler().nProcs(args.path());
 
@@ -612,31 +790,33 @@ int main(int argc, char *argv[])
         }
 
 
-        // Read point on individual processors to determine merge tolerance
-        // (otherwise single cell domains might give problems)
-
-        const boundBox bb = procBounds(args, databases, regionDir);
-        const scalar mergeDist = mergeTol*bb.mag();
-
-        Info<< "Overall mesh bounding box  : " << bb << nl
-            << "Relative tolerance         : " << mergeTol << nl
-            << "Absolute matching distance : " << mergeDist << nl
-            << endl;
-
-
         // Addressing from processor to reconstructed case
         labelListList cellProcAddressing(nProcs);
         labelListList faceProcAddressing(nProcs);
         labelListList pointProcAddressing(nProcs);
         labelListList boundaryProcAddressing(nProcs);
 
+
         // Internal faces on the final reconstructed mesh
         label masterInternalFaces;
 
         // Owner addressing on the final reconstructed mesh
         labelList masterOwner;
 
+        if (procMatch)
         {
+            // Read point on individual processors to determine merge tolerance
+            // (otherwise single cell domains might give problems)
+
+            const boundBox bb = procBounds(args, databases, regionDir);
+            const scalar mergeDist = mergeTol*bb.mag();
+
+            Info<< "Overall mesh bounding box  : " << bb << nl
+                << "Relative tolerance         : " << mergeTol << nl
+                << "Absolute matching distance : " << mergeDist << nl
+                << endl;
+
+
             // Construct empty mesh.
             // fvMesh** masterMesh = new fvMesh*[nProcs];
             PtrList<fvMesh> masterMesh(nProcs);
@@ -818,29 +998,187 @@ int main(int argc, char *argv[])
             masterInternalFaces = masterMesh[0].nInternalFaces();
             masterOwner = masterMesh[0].faceOwner();
 
-
-            Info<< "\nWriting merged mesh to "
-                << runTime.path()/runTime.timeName()
-                << nl << endl;
-
-            if (!masterMesh[0].write())
-            {
-                FatalErrorInFunction
-                    << "Failed writing polyMesh."
-                    << exit(FatalError);
-            }
-            topoSet::removeFiles(masterMesh[0]);
-
+            // Write reconstructed mesh
+            writeMesh(masterMesh[0], cellProcAddressing);
             if (writeCellDist)
             {
-                writeCellDistance
-                (
-                    runTime,
-                    masterMesh[0],
-                    cellProcAddressing
-                );
+                writeDistribution(runTime, masterMesh[0], cellProcAddressing);
             }
         }
+        else
+        {
+            // Load all meshes
+            PtrList<fvMesh> fvMeshes(nProcs);
+            for (label proci=0; proci<nProcs; proci++)
+            {
+                fvMeshes.set
+                (
+                    proci,
+                    new fvMesh
+                    (
+                        IOobject
+                        (
+                            regionName,
+                            databases[proci].timeName(),
+                            databases[proci]
+                        )
+                    )
+                );
+            }
+
+            // Construct pointers to meshes
+            UPtrList<polyMesh> meshes(fvMeshes.size());
+            forAll(fvMeshes, proci)
+            {
+                meshes.set(proci, &fvMeshes[proci]);
+            }
+
+            // Collect statistics
+            label nCells = 0;
+            label nFaces = 0;
+            label nPoints = 0;
+            forAll(meshes, proci)
+            {
+                const polyMesh& mesh = meshes[proci];
+                nCells += mesh.nCells();
+                nFaces += mesh.nFaces();
+                nPoints += mesh.nPoints();
+            }
+
+            // Get pairs of patches to stitch. These pairs have to
+            // - have ordered, opposite faces (so one to one correspondence)
+            List<DynamicList<label>> localPatch;
+            List<DynamicList<label>> remoteProc;
+            List<DynamicList<label>> remotePatch;
+            const label nGlobalPatches = polyMeshAdder::procPatchPairs
+            (
+                meshes,
+                localPatch,
+                remoteProc,
+                remotePatch
+            );
+
+            // Collect matching boundary faces on patches-to-stitch
+            labelListList localBoundaryFace;
+            labelListList remoteFaceProc;
+            labelListList remoteBoundaryFace;
+            polyMeshAdder::patchFacePairs
+            (
+                meshes,
+                localPatch,
+                remoteProc,
+                remotePatch,
+                localBoundaryFace,
+                remoteFaceProc,
+                remoteBoundaryFace
+            );
+
+            // All matched faces assumed to have vertex0 matched
+            labelListList remoteFaceStart(meshes.size());
+            forAll(meshes, proci)
+            {
+                const labelList& procFaces = localBoundaryFace[proci];
+                remoteFaceStart[proci].setSize(procFaces.size(), 0);
+            }
+
+
+
+            labelListList patchMap(meshes.size());
+            labelListList pointZoneMap(meshes.size());
+            labelListList faceZoneMap(meshes.size());
+            labelListList cellZoneMap(meshes.size());
+            forAll(meshes, proci)
+            {
+                const polyMesh& mesh = meshes[proci];
+                patchMap[proci] = identity(mesh.boundaryMesh().size());
+
+                // Remove excess patches
+                patchMap[proci].setSize(nGlobalPatches);
+
+                pointZoneMap[proci] = identity(mesh.pointZones().size());
+                faceZoneMap[proci] = identity(mesh.faceZones().size());
+                cellZoneMap[proci] = identity(mesh.cellZones().size());
+            }
+
+
+            refPtr<fvMesh> masterMeshPtr;
+            {
+                // Do in-place addition on proc0.
+
+                const labelList oldFaceOwner(fvMeshes[0].faceOwner());
+
+                fvMeshAdder::add
+                (
+                    0,              // index of mesh to modify (== mesh_)
+                    fvMeshes,
+                    oldFaceOwner,
+
+                    // Coupling info
+                    localBoundaryFace,
+                    remoteFaceProc,
+                    remoteBoundaryFace,
+
+                    boundaryProcAddressing,
+                    cellProcAddressing,
+                    faceProcAddressing,
+                    pointProcAddressing
+                );
+
+                // Remove zero-faces processor patches
+                const polyBoundaryMesh& pbm = fvMeshes[0].boundaryMesh();
+                labelList oldToNew(pbm.size(), -1);
+                label newi = 0;
+                // Non processor patches first
+                forAll(pbm, patchi)
+                {
+                    const auto& pp = pbm[patchi];
+                    if (!isA<processorPolyPatch>(pp) || pp.size())
+                    {
+                        oldToNew[patchi] = newi++;
+                    }
+                }
+                const label nNonProcPatches = newi;
+
+                // Move all deletable patches to the end
+                forAll(oldToNew, patchi)
+                {
+                    if (oldToNew[patchi] == -1)
+                    {
+                        oldToNew[patchi] = newi++;
+                    }
+                }
+                fvMeshTools::reorderPatches
+                (
+                    fvMeshes[0],
+                    oldToNew,
+                    nNonProcPatches,
+                    false
+                );
+
+                masterMeshPtr = fvMeshes[0];
+            }
+
+
+            const fvMesh& masterMesh = masterMeshPtr();
+
+            // Number of internal faces on the final reconstructed mesh
+            masterInternalFaces = masterMesh.nInternalFaces();
+            // Owner addressing on the final reconstructed mesh
+            masterOwner = masterMesh.faceOwner();
+
+            // Write reconstructed mesh
+            const word oldCaseName = masterMesh.time().caseName();
+            const_cast<Time&>(masterMesh.time()).caseName() =
+                runTime.caseName();
+
+            writeMesh(masterMesh, cellProcAddressing);
+            if (writeCellDist)
+            {
+                writeDistribution(runTime, masterMesh, cellProcAddressing);
+            }
+            const_cast<Time&>(masterMesh.time()).caseName() = oldCaseName;
+        }
+
 
 
         // Write the addressing
@@ -863,143 +1201,16 @@ int main(int argc, char *argv[])
                 )
             );
 
-
-            // From processor point to reconstructed mesh point
-
-            Info<< "Writing pointProcAddressing to "
-                << databases[proci].caseName()
-                  /procMesh.facesInstance()
-                  /polyMesh::meshSubDir
-                << endl;
-
-            labelIOList
+            writeMaps
             (
-                IOobject
-                (
-                    "pointProcAddressing",
-                    procMesh.facesInstance(),
-                    polyMesh::meshSubDir,
-                    procMesh,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE,
-                    false                       // Do not register
-                ),
-                pointProcAddressing[proci]
-            ).write();
-
-
-            // From processor face to reconstructed mesh face
-
-            Info<< "Writing faceProcAddressing to "
-                << databases[proci].caseName()
-                  /procMesh.facesInstance()
-                  /polyMesh::meshSubDir
-                << endl;
-
-            labelIOList faceProcAddr
-            (
-                IOobject
-                (
-                    "faceProcAddressing",
-                    procMesh.facesInstance(),
-                    polyMesh::meshSubDir,
-                    procMesh,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE,
-                    false                       // Do not register
-                ),
-                faceProcAddressing[proci]
-            );
-
-            // Now add turning index to faceProcAddressing.
-            // See reconstructPar for meaning of turning index.
-            forAll(faceProcAddr, procFacei)
-            {
-                const label masterFacei = faceProcAddr[procFacei];
-
-                if
-                (
-                   !procMesh.isInternalFace(procFacei)
-                 && masterFacei < masterInternalFaces
-                )
-                {
-                    // proc face is now external but used to be internal face.
-                    // Check if we have owner or neighbour.
-
-                    label procOwn = procMesh.faceOwner()[procFacei];
-                    label masterOwn = masterOwner[masterFacei];
-
-                    if (cellProcAddressing[proci][procOwn] == masterOwn)
-                    {
-                        // No turning. Offset by 1.
-                        faceProcAddr[procFacei]++;
-                    }
-                    else
-                    {
-                        // Turned face.
-                        faceProcAddr[procFacei] =
-                            -1 - faceProcAddr[procFacei];
-                    }
-                }
-                else
-                {
-                    // No turning. Offset by 1.
-                    faceProcAddr[procFacei]++;
-                }
-            }
-
-            faceProcAddr.write();
-
-
-            // From processor cell to reconstructed mesh cell
-
-            Info<< "Writing cellProcAddressing to "
-                << databases[proci].caseName()
-                  /procMesh.facesInstance()
-                  /polyMesh::meshSubDir
-                << endl;
-
-            labelIOList
-            (
-                IOobject
-                (
-                    "cellProcAddressing",
-                    procMesh.facesInstance(),
-                    polyMesh::meshSubDir,
-                    procMesh,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE,
-                    false                       // Do not register
-                ),
-                cellProcAddressing[proci]
-            ).write();
-
-
-
-            // From processor patch to reconstructed mesh patch
-
-            Info<< "Writing boundaryProcAddressing to "
-                << databases[proci].caseName()
-                  /procMesh.facesInstance()
-                  /polyMesh::meshSubDir
-                << endl;
-
-            labelIOList
-            (
-                IOobject
-                (
-                    "boundaryProcAddressing",
-                    procMesh.facesInstance(),
-                    polyMesh::meshSubDir,
-                    procMesh,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE,
-                    false                       // Do not register
-                ),
+                masterInternalFaces,
+                masterOwner,
+                procMesh,
+                cellProcAddressing[proci],
+                faceProcAddressing[proci],
+                pointProcAddressing[proci],
                 boundaryProcAddressing[proci]
-            ).write();
-
-            Info<< endl;
+            );
         }
     }
 
diff --git a/applications/utilities/parallelProcessing/redistributePar/redistributePar.C b/applications/utilities/parallelProcessing/redistributePar/redistributePar.C
index b2d6547ded..e1d2d5d26d 100644
--- a/applications/utilities/parallelProcessing/redistributePar/redistributePar.C
+++ b/applications/utilities/parallelProcessing/redistributePar/redistributePar.C
@@ -107,52 +107,6 @@ using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// Tolerance (as fraction of the bounding box). Needs to be fairly lax since
-// usually meshes get written with limited precision (6 digits)
-static const scalar defaultMergeTol = 1e-6;
-
-
-// Get merging distance when matching face centres
-scalar getMergeDistance
-(
-    const argList& args,
-    const Time& runTime,
-    const boundBox& bb
-)
-{
-    const scalar mergeTol =
-        args.getOrDefault<scalar>("mergeTol", defaultMergeTol);
-
-    const scalar writeTol =
-        Foam::pow(scalar(10), -scalar(IOstream::defaultPrecision()));
-
-    Info<< "Merge tolerance : " << mergeTol << nl
-        << "Write tolerance : " << writeTol << endl;
-
-    if (runTime.writeFormat() == IOstream::ASCII && mergeTol < writeTol)
-    {
-        FatalErrorInFunction
-            << "Your current settings specify ASCII writing with "
-            << IOstream::defaultPrecision() << " digits precision." << nl
-            << "Your merging tolerance (" << mergeTol << ") is finer than this."
-            << nl
-            << "Please change your writeFormat to binary"
-            << " or increase the writePrecision" << endl
-            << "or adjust the merge tolerance (-mergeTol)."
-            << exit(FatalError);
-    }
-
-    const scalar mergeDist = mergeTol * bb.mag();
-
-    Info<< "Overall meshes bounding box : " << bb << nl
-        << "Relative tolerance          : " << mergeTol << nl
-        << "Absolute matching distance  : " << mergeDist << nl
-        << endl;
-
-    return mergeDist;
-}
-
-
 void setBasicGeometry(fvMesh& mesh)
 {
     // Set the fvGeometryScheme to basic since it does not require
@@ -860,7 +814,6 @@ void correctCoupledBoundaryConditions(fvMesh& mesh)
 autoPtr<mapDistributePolyMesh> redistributeAndWrite
 (
     const Time& baseRunTime,
-    const scalar tolDim,
     const boolList& haveMesh,
     const fileName& meshSubDir,
     const bool doReadFields,
@@ -1150,7 +1103,7 @@ autoPtr<mapDistributePolyMesh> redistributeAndWrite
 
 
     // Mesh distribution engine
-    fvMeshDistribute distributor(mesh, tolDim);
+    fvMeshDistribute distributor(mesh);
 
     // Do all the distribution of mesh and fields
     autoPtr<mapDistributePolyMesh> rawMap = distributor.distribute(decomp);
@@ -2317,12 +2270,6 @@ int main(int argc, char *argv[])
         "Test without writing the decomposition. "
         "Changes -cellDist to only write volScalarField."
     );
-    argList::addOption
-    (
-        "mergeTol",
-        "scalar",
-        "The merge distance relative to the bounding box size (default 1e-6)"
-    );
     argList::addBoolOption
     (
         "cellDist",
@@ -2707,14 +2654,6 @@ int main(int argc, char *argv[])
                     // problems. See comment in routine
                     setBasicGeometry(mesh);
 
-                    // Global matching tolerance
-                    const scalar tolDim = getMergeDistance
-                    (
-                        args,
-                        runTime,
-                        mesh.bounds()
-                    );
-
 
                     // Determine decomposition
                     // ~~~~~~~~~~~~~~~~~~~~~~~
@@ -2728,7 +2667,6 @@ int main(int argc, char *argv[])
                     redistributeAndWrite
                     (
                         baseRunTime,
-                        tolDim,
                         haveMesh,
                         meshSubDir,
                         false,      // do not read fields
@@ -3060,14 +2998,6 @@ int main(int argc, char *argv[])
             //    << " nPatches:" << mesh.boundaryMesh().size() << endl;
 
 
-            // Global matching tolerance
-            const scalar tolDim = getMergeDistance
-            (
-                args,
-                runTime,
-                mesh.bounds()
-            );
-
             // Determine decomposition
             // ~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -3139,7 +3069,6 @@ int main(int argc, char *argv[])
             autoPtr<mapDistributePolyMesh> distMap = redistributeAndWrite
             (
                 baseRunTime,
-                tolDim,
                 haveMesh,
                 meshSubDir,
                 true,           // read fields
diff --git a/src/dynamicMesh/fvMeshAdder/fvMeshAdder.C b/src/dynamicMesh/fvMeshAdder/fvMeshAdder.C
index db8400b0e4..2b0925c30d 100644
--- a/src/dynamicMesh/fvMeshAdder/fvMeshAdder.C
+++ b/src/dynamicMesh/fvMeshAdder/fvMeshAdder.C
@@ -317,9 +317,15 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshAdder::add
         if (meshes.set(meshi))
         {
             constructPatchMap[meshi] = patchMap[meshi];
+            constructPatchMap[meshi].setSize
+            (
+                meshes[meshi].boundaryMesh().size(),
+                -1
+            );
         }
     }
 
+
     // Map all fields
     fvMeshAdder::MapVolFields<scalar>
     (
diff --git a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C
index d72d4cac10..23c5799651 100644
--- a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C
+++ b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C
@@ -704,8 +704,15 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::mergeSharedPoints
         }
     }
 
+    if (debug)
+    {
+        Pout<< "mergeSharedPoints : found " << nShared
+            << " points on processor boundaries" << nl << endl;
+    }
+
     Map<label> globalMasterToLocalMaster(2*nShared);
     Map<label> pointToMaster(2*nShared);
+    label nMatch = 0;
 
     forAll(pointToGlobalMaster, pointi)
     {
@@ -716,7 +723,9 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::mergeSharedPoints
 
             if (iter.found())
             {
+                // Matched to existing master
                 pointToMaster.insert(pointi, *iter);
+                nMatch++;
             }
             else
             {
@@ -727,7 +736,16 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::mergeSharedPoints
         }
     }
 
-    if (returnReduce(pointToMaster.size(), sumOp<label>()) == 0)
+    reduce(nMatch, sumOp<label>());
+
+    if (debug)
+    {
+        Pout<< "mergeSharedPoints : found "
+            << nMatch << " mergeable points" << nl << endl;
+    }
+
+
+    if (nMatch == 0)
     {
         return nullptr;
     }
@@ -1120,6 +1138,141 @@ void Foam::fvMeshDistribute::findCouples
 }
 
 
+void Foam::fvMeshDistribute::findCouples
+(
+    const UPtrList<polyMesh>& meshes,
+    const PtrList<labelList>& domainSourceFaces,
+    const PtrList<labelList>& domainSourceProcs,
+    const PtrList<labelList>& domainSourcePatchs,
+
+    labelListList& localBoundaryFace,
+    labelListList& remoteFaceProc,
+    labelListList& remoteBoundaryFace
+)
+{
+    // Collect all matching processor face pairs. These are all the
+    // faces for which we have both sides
+
+    // Pass 0: count number of inter-processor faces
+    labelList nProcFaces(meshes.size(), 0);
+    forAll(meshes, meshi)
+    {
+        if (meshes.set(meshi))
+        {
+            const labelList& domainProc = domainSourceProcs[meshi];
+            const labelList& domainPatch = domainSourcePatchs[meshi];
+
+            forAll(domainProc, bFacei)
+            {
+                if (domainProc[bFacei] != -1 && domainPatch[bFacei] == -1)
+                {
+                    nProcFaces[meshi]++;
+                }
+            }
+        }
+    }
+
+    if (debug)
+    {
+        Pout<< "fvMeshDistribute::findCouples : nProcFaces:"
+            << flatOutput(nProcFaces) << endl;
+    }
+
+
+    // Size
+    List<DynamicList<label>> dynLocalFace(Pstream::nProcs());
+    List<DynamicList<label>> dynRemoteProc(Pstream::nProcs());
+    List<DynamicList<label>> dynRemoteFace(Pstream::nProcs());
+
+    forAll(meshes, meshi)
+    {
+        if (meshes.set(meshi))
+        {
+            dynLocalFace[meshi].setCapacity(nProcFaces[meshi]);
+            dynRemoteProc[meshi].setCapacity(nProcFaces[meshi]);
+            dynRemoteFace[meshi].setCapacity(nProcFaces[meshi]);
+        }
+    }
+
+
+    // Insert all into big map. Find matches
+    HashTable<labelPair, labelPair, labelPair::Hash<>> map(2*sum(nProcFaces));
+
+    nProcFaces = 0;
+
+    forAll(meshes, meshi)
+    {
+        if (meshes.set(meshi))
+        {
+            const labelList& domainFace = domainSourceFaces[meshi];
+            const labelList& domainProc = domainSourceProcs[meshi];
+            const labelList& domainPatch = domainSourcePatchs[meshi];
+
+            forAll(domainProc, bFacei)
+            {
+                if (domainProc[bFacei] != -1 && domainPatch[bFacei] == -1)
+                {
+                    const labelPair key
+                    (
+                        domainProc[bFacei],
+                        domainFace[bFacei]
+                    );
+                    auto fnd = map.find(key);
+
+                    if (!fnd.found())
+                    {
+                        // Insert
+                        map.emplace(key, meshi, bFacei);
+                    }
+                    else
+                    {
+                        // Second occurence. Found match.
+                        const label matchProci = fnd().first();
+                        const label matchFacei = fnd().second();
+
+                        dynLocalFace[meshi].append(bFacei);
+                        dynRemoteProc[meshi].append(matchProci);
+                        dynRemoteFace[meshi].append(matchFacei);
+                        nProcFaces[meshi]++;
+
+                        dynLocalFace[matchProci].append(matchFacei);
+                        dynRemoteProc[matchProci].append(meshi);
+                        dynRemoteFace[matchProci].append(bFacei);
+                        nProcFaces[matchProci]++;
+                    }
+                }
+            }
+        }
+    }
+
+    if (debug)
+    {
+        Pout<< "fvMeshDistribute::findCouples : stored procFaces:"
+            << map.size() << endl;
+    }
+
+    localBoundaryFace.setSize(Pstream::nProcs());
+    remoteFaceProc.setSize(Pstream::nProcs());
+    remoteBoundaryFace.setSize(Pstream::nProcs());
+    forAll(meshes, meshi)
+    {
+        if (meshes.set(meshi))
+        {
+            localBoundaryFace[meshi] = std::move(dynLocalFace[meshi]);
+            remoteFaceProc[meshi] = std::move(dynRemoteProc[meshi]);
+            remoteBoundaryFace[meshi] = std::move(dynRemoteFace[meshi]);
+        }
+    }
+
+
+    if (debug)
+    {
+        Pout<< "fvMeshDistribute::findCouples : found matches:"
+            << flatOutput(nProcFaces) << endl;
+    }
+}
+
+
 // Map data on boundary faces to new mesh (resulting from adding two meshes)
 Foam::labelList Foam::fvMeshDistribute::mapBoundaryData
 (
@@ -1676,10 +1829,10 @@ Foam::autoPtr<Foam::fvMesh> Foam::fvMeshDistribute::receiveMesh
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::fvMeshDistribute::fvMeshDistribute(fvMesh& mesh, const scalar mergeTol)
+Foam::fvMeshDistribute::fvMeshDistribute(fvMesh& mesh)//, const scalar mergeTol)
 :
-    mesh_(mesh),
-    mergeTol_(mergeTol)
+    mesh_(mesh)
+    //mergeTol_(mergeTol)
 {}
 
 
@@ -2165,7 +2318,21 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
 
 
     // Start sending&receiving from buffers
-    pBufs.finishedSends();
+    {
+        if (debug)
+        {
+            Pout<< "Starting sending" << endl;
+        }
+
+        labelList recvSizes;
+        pBufs.finishedSends(recvSizes);
+
+        if (debug)
+        {
+            Pout<< "Finished sending and receiving : " << flatOutput(recvSizes)
+                << endl;
+        }
+    }
 
 
     // Subset the part that stays
@@ -2208,14 +2375,7 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
         subPointMap[Pstream::myProcNo()] = subMap().pointMap();
         subPatchMap[Pstream::myProcNo()] = identity(patches.size());
 
-        // Initialize all addressing into current mesh
-        constructCellMap[Pstream::myProcNo()] = identity(mesh_.nCells());
-        constructFaceMap[Pstream::myProcNo()] = identity(mesh_.nFaces(), 1);
-        constructPointMap[Pstream::myProcNo()] = identity(mesh_.nPoints());
-        constructPatchMap[Pstream::myProcNo()] = identity(patches.size());
-
-        // Subset the mesh data: neighbourCell/neighbourProc
-        // fields
+        // Subset the mesh data: neighbourCell/neighbourProc fields
         labelList domainSourceFace;
         labelList domainSourceProc;
         labelList domainSourcePatch;
@@ -2281,6 +2441,41 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
     // Disable parallel. Original state already known.
     UPstream::parRun(false);
 
+    PtrList<labelList> domainSourceFaces(Pstream::nProcs());
+    PtrList<labelList> domainSourceProcs(Pstream::nProcs());
+    PtrList<labelList> domainSourcePatchs(Pstream::nProcs());
+    PtrList<labelList> domainSourceNewNbrProcs(Pstream::nProcs());
+    PtrList<labelList> domainSourcePointMasters(Pstream::nProcs());
+
+    PtrList<fvMesh> domainMeshPtrs(Pstream::nProcs());
+
+    PtrList<PtrList<volScalarField>> vsfs(Pstream::nProcs());
+    PtrList<PtrList<volVectorField>> vvfs(Pstream::nProcs());
+    PtrList<PtrList<volSphericalTensorField>> vsptfs(Pstream::nProcs());
+    PtrList<PtrList<volSymmTensorField>> vsytfs(Pstream::nProcs());
+    PtrList<PtrList<volTensorField>> vtfs(Pstream::nProcs());
+
+    PtrList<PtrList<surfaceScalarField>> ssfs(Pstream::nProcs());
+    PtrList<PtrList<surfaceVectorField>> svfs(Pstream::nProcs());
+    PtrList<PtrList<surfaceSphericalTensorField>> ssptfs
+    (
+        Pstream::nProcs()
+    );
+    PtrList<PtrList<surfaceSymmTensorField>> ssytfs(Pstream::nProcs());
+    PtrList<PtrList<surfaceTensorField>> stfs(Pstream::nProcs());
+
+    PtrList<PtrList<volScalarField::Internal>> dsfs(Pstream::nProcs());
+    PtrList<PtrList<volVectorField::Internal>> dvfs(Pstream::nProcs());
+    PtrList<PtrList<volSphericalTensorField::Internal>> dstfs
+    (
+        Pstream::nProcs()
+    );
+    PtrList<PtrList<volSymmTensorField::Internal>> dsytfs
+    (
+        Pstream::nProcs()
+    );
+    PtrList<PtrList<volTensorField::Internal>> dtfs(Pstream::nProcs());
+
     forAll(nRevcCells, sendProc)
     {
         // Did processor sendProc send anything to me?
@@ -2301,36 +2496,26 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
 
 
             // Receive from sendProc
-            labelList domainSourceFace;
-            labelList domainSourceProc;
-            labelList domainSourcePatch;
-            labelList domainSourceNewNbrProc;
-            labelList domainSourcePointMaster;
+            domainSourceFaces.set(sendProc, new labelList(0));
+            labelList& domainSourceFace = domainSourceFaces[sendProc];
 
-            autoPtr<fvMesh> domainMeshPtr;
+            domainSourceProcs.set(sendProc, new labelList(0));
+            labelList& domainSourceProc = domainSourceProcs[sendProc];
 
-            PtrList<volScalarField> vsf;
-            PtrList<volVectorField> vvf;
-            PtrList<volSphericalTensorField> vsptf;
-            PtrList<volSymmTensorField> vsytf;
-            PtrList<volTensorField> vtf;
+            domainSourcePatchs.set(sendProc, new labelList(0));
+            labelList& domainSourcePatch = domainSourcePatchs[sendProc];
 
-            PtrList<surfaceScalarField> ssf;
-            PtrList<surfaceVectorField> svf;
-            PtrList<surfaceSphericalTensorField> ssptf;
-            PtrList<surfaceSymmTensorField> ssytf;
-            PtrList<surfaceTensorField> stf;
-
-            PtrList<volScalarField::Internal> dsf;
-            PtrList<volVectorField::Internal> dvf;
-            PtrList<volSphericalTensorField::Internal> dstf;
-            PtrList<volSymmTensorField::Internal> dsytf;
-            PtrList<volTensorField::Internal> dtf;
+            domainSourceNewNbrProcs.set(sendProc, new labelList(0));
+            labelList& domainSourceNewNbrProc =
+                domainSourceNewNbrProcs[sendProc];
 
+            domainSourcePointMasters.set(sendProc, new labelList(0));
+            labelList& domainSourcePointMaster =
+                domainSourcePointMasters[sendProc];
 
             // Opposite of sendMesh
             {
-                domainMeshPtr = receiveMesh
+                autoPtr<fvMesh> domainMeshPtr = receiveMesh
                 (
                     sendProc,
                     pointZoneNames,
@@ -2345,7 +2530,8 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
                     domainSourcePointMaster,
                     str
                 );
-                fvMesh& domainMesh = domainMeshPtr();
+                domainMeshPtrs.set(sendProc, domainMeshPtr.ptr());
+                fvMesh& domainMesh = domainMeshPtrs[sendProc];
                 // Force construction of various on mesh.
                 //(void)domainMesh.globalData();
 
@@ -2355,347 +2541,362 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
                 dictionary fieldDicts(str);
 
                 // Vol fields
+                vsfs.set(sendProc, new PtrList<volScalarField>(0));
                 receiveFields<volScalarField>
                 (
                     sendProc,
                     allFieldNames,
                     domainMesh,
-                    vsf,
+                    vsfs[sendProc],
                     fieldDicts
                 );
+                vvfs.set(sendProc, new PtrList<volVectorField>(0));
                 receiveFields<volVectorField>
                 (
                     sendProc,
                     allFieldNames,
                     domainMesh,
-                    vvf,
+                    vvfs[sendProc],
                     fieldDicts
                 );
+                vsptfs.set
+                (
+                    sendProc,
+                    new PtrList<volSphericalTensorField>(0)
+                );
                 receiveFields<volSphericalTensorField>
                 (
                     sendProc,
                     allFieldNames,
                     domainMesh,
-                    vsptf,
+                    vsptfs[sendProc],
                     fieldDicts
                 );
+                vsytfs.set(sendProc, new PtrList<volSymmTensorField>(0));
                 receiveFields<volSymmTensorField>
                 (
                     sendProc,
                     allFieldNames,
                     domainMesh,
-                    vsytf,
+                    vsytfs[sendProc],
                     fieldDicts
                 );
+                vtfs.set(sendProc, new PtrList<volTensorField>(0));
                 receiveFields<volTensorField>
                 (
                     sendProc,
                     allFieldNames,
                     domainMesh,
-                    vtf,
+                    vtfs[sendProc],
                     fieldDicts
                 );
 
                 // Surface fields
+                ssfs.set(sendProc, new PtrList<surfaceScalarField>(0));
                 receiveFields<surfaceScalarField>
                 (
                     sendProc,
                     allFieldNames,
                     domainMesh,
-                    ssf,
+                    ssfs[sendProc],
                     fieldDicts
                 );
+                svfs.set(sendProc, new PtrList<surfaceVectorField>(0));
                 receiveFields<surfaceVectorField>
                 (
                     sendProc,
                     allFieldNames,
                     domainMesh,
-                    svf,
+                    svfs[sendProc],
                     fieldDicts
                 );
+                ssptfs.set
+                (
+                    sendProc,
+                    new PtrList<surfaceSphericalTensorField>(0)
+                );
                 receiveFields<surfaceSphericalTensorField>
                 (
                     sendProc,
                     allFieldNames,
                     domainMesh,
-                    ssptf,
+                    ssptfs[sendProc],
                     fieldDicts
                 );
+                ssytfs.set(sendProc, new PtrList<surfaceSymmTensorField>(0));
                 receiveFields<surfaceSymmTensorField>
                 (
                     sendProc,
                     allFieldNames,
                     domainMesh,
-                    ssytf,
+                    ssytfs[sendProc],
                     fieldDicts
                 );
+                stfs.set(sendProc, new PtrList<surfaceTensorField>(0));
                 receiveFields<surfaceTensorField>
                 (
                     sendProc,
                     allFieldNames,
                     domainMesh,
-                    stf,
+                    stfs[sendProc],
                     fieldDicts
                 );
 
                 // Dimensioned fields
+                dsfs.set
+                (
+                    sendProc,
+                    new PtrList<volScalarField::Internal>(0)
+                );
                 receiveFields<volScalarField::Internal>
                 (
                     sendProc,
                     allFieldNames,
                     domainMesh,
-                    dsf,
+                    dsfs[sendProc],
                     fieldDicts
                 );
+                dvfs.set
+                (
+                    sendProc,
+                    new PtrList<volVectorField::Internal>(0)
+                );
                 receiveFields<volVectorField::Internal>
                 (
                     sendProc,
                     allFieldNames,
                     domainMesh,
-                    dvf,
+                    dvfs[sendProc],
                     fieldDicts
                 );
+                dstfs.set
+                (
+                    sendProc,
+                    new PtrList<volSphericalTensorField::Internal>(0)
+                );
                 receiveFields<volSphericalTensorField::Internal>
                 (
                     sendProc,
                     allFieldNames,
                     domainMesh,
-                    dstf,
+                    dstfs[sendProc],
                     fieldDicts
                 );
+                dsytfs.set
+                (
+                    sendProc,
+                    new PtrList<volSymmTensorField::Internal>(0)
+                );
                 receiveFields<volSymmTensorField::Internal>
                 (
                     sendProc,
                     allFieldNames,
                     domainMesh,
-                    dsytf,
+                    dsytfs[sendProc],
                     fieldDicts
                 );
+                dtfs.set
+                (
+                    sendProc,
+                    new PtrList<volTensorField::Internal>(0)
+                );
                 receiveFields<volTensorField::Internal>
                 (
                     sendProc,
                     allFieldNames,
                     domainMesh,
-                    dtf,
+                    dtfs[sendProc],
                     fieldDicts
                 );
             }
-            const fvMesh& domainMesh = domainMeshPtr();
+        }
+    }
+
+    // Clear out storage
+    pBufs.clear();
 
 
-            constructCellMap[sendProc] = identity(domainMesh.nCells());
-            constructFaceMap[sendProc] = identity(domainMesh.nFaces(), 1);
-            constructPointMap[sendProc] = identity(domainMesh.nPoints());
-            constructPatchMap[sendProc] =
-                identity(domainMesh.boundaryMesh().size());
+    // Set up pointers to meshes so we can include our mesh_
+    UPtrList<polyMesh> meshes(domainMeshPtrs.size());
+    UPtrList<fvMesh> fvMeshes(domainMeshPtrs.size());
+    forAll(domainMeshPtrs, proci)
+    {
+        if (domainMeshPtrs.set(proci))
+        {
+            meshes.set(proci, &domainMeshPtrs[proci]);
+            fvMeshes.set(proci, &domainMeshPtrs[proci]);
+        }
+    }
+
+    // 'Receive' from myself
+    {
+        meshes.set(Pstream::myProcNo(), &mesh_);
+        fvMeshes.set(Pstream::myProcNo(), &mesh_);
+
+        //domainSourceFaces.set(Pstream::myProcNo(), std::move(sourceFace));
+        domainSourceFaces.set(Pstream::myProcNo(), new labelList(0));
+        domainSourceFaces[Pstream::myProcNo()] = sourceFace;
+
+        domainSourceProcs.set(Pstream::myProcNo(), new labelList(0));
+        //std::move(sourceProc));
+        domainSourceProcs[Pstream::myProcNo()] = sourceProc;
+
+        domainSourcePatchs.set(Pstream::myProcNo(), new labelList(0));
+        //, std::move(sourcePatch));
+        domainSourcePatchs[Pstream::myProcNo()] = sourcePatch;
+
+        domainSourceNewNbrProcs.set(Pstream::myProcNo(), new labelList(0));
+        domainSourceNewNbrProcs[Pstream::myProcNo()] = sourceNewNbrProc;
+
+        domainSourcePointMasters.set(Pstream::myProcNo(), new labelList(0));
+        domainSourcePointMasters[Pstream::myProcNo()] = sourcePointMaster;
+    }
 
 
-            // Print a bit.
-            if (debug)
+    // Find matching faces that need to be stitched
+    labelListList localBoundaryFace(Pstream::nProcs());
+    labelListList remoteFaceProc(Pstream::nProcs());
+    labelListList remoteBoundaryFace(Pstream::nProcs());
+    findCouples
+    (
+        meshes,
+        domainSourceFaces,
+        domainSourceProcs,
+        domainSourcePatchs,
+
+        localBoundaryFace,
+        remoteFaceProc,
+        remoteBoundaryFace
+    );
+
+
+    const label nOldInternalFaces = mesh_.nInternalFaces();
+    const labelList oldFaceOwner(mesh_.faceOwner());
+
+    fvMeshAdder::add
+    (
+        Pstream::myProcNo(),    // index of mesh to modify (== mesh_)
+        fvMeshes,
+        oldFaceOwner,
+
+        // Coupling info
+        localBoundaryFace,
+        remoteFaceProc,
+        remoteBoundaryFace,
+
+        constructPatchMap,
+        constructCellMap,
+        constructFaceMap,
+        constructPointMap
+    );
+
+    if (debug)
+    {
+        Pout<< nl << "ADDED REMOTE MESHES:" << endl;
+        printMeshInfo(mesh_);
+        printFieldInfo<volScalarField>(mesh_);
+        printFieldInfo<volVectorField>(mesh_);
+        printFieldInfo<volSphericalTensorField>(mesh_);
+        printFieldInfo<volSymmTensorField>(mesh_);
+        printFieldInfo<volTensorField>(mesh_);
+        printFieldInfo<surfaceScalarField>(mesh_);
+        printFieldInfo<surfaceVectorField>(mesh_);
+        printFieldInfo<surfaceSphericalTensorField>(mesh_);
+        printFieldInfo<surfaceSymmTensorField>(mesh_);
+        printFieldInfo<surfaceTensorField>(mesh_);
+        Pout<< nl << endl;
+    }
+
+    {
+        //- Combine sourceProc, sourcePatch, sourceFace
+        sourceProc.setSize(mesh_.nBoundaryFaces());
+        sourceProc = -1;
+        sourcePatch.setSize(mesh_.nBoundaryFaces());
+        sourcePatch = -1;
+        sourceFace.setSize(mesh_.nBoundaryFaces());
+        sourceFace = -1;
+        sourceNewNbrProc.setSize(mesh_.nBoundaryFaces());
+        sourceNewNbrProc = -1;
+        sourcePointMaster.setSize(mesh_.nPoints());
+        sourcePointMaster = -1;
+
+        if (mesh_.nPoints() > 0)
+        {
+            forAll(meshes, meshi)
             {
-                Pout<< nl << "RECEIVED MESH FROM:" << sendProc << endl;
-                printMeshInfo(domainMesh);
-                printFieldInfo<volScalarField>(domainMesh);
-                printFieldInfo<volVectorField>(domainMesh);
-                printFieldInfo<volSphericalTensorField>(domainMesh);
-                printFieldInfo<volSymmTensorField>(domainMesh);
-                printFieldInfo<volTensorField>(domainMesh);
-                printFieldInfo<surfaceScalarField>(domainMesh);
-                printFieldInfo<surfaceVectorField>(domainMesh);
-                printFieldInfo<surfaceSphericalTensorField>(domainMesh);
-                printFieldInfo<surfaceSymmTensorField>(domainMesh);
-                printFieldInfo<surfaceTensorField>(domainMesh);
-            }
-
-
-            // Now this mesh we received (from sendProc) needs to be merged
-            // with the current mesh. On the current mesh we have for all
-            // boundaryfaces the original face and processor. See if we can
-            // match these up to the received domainSourceFace and
-            // domainSourceProc.
-            labelList masterCoupledFaces;
-            labelList slaveCoupledFaces;
-            findCouples
-            (
-                mesh_,
-
-                sourceFace,
-                sourceProc,
-                sourcePatch,
-
-                sendProc,
-                domainMesh,
-                domainSourceFace,
-                domainSourceProc,
-                domainSourcePatch,
-
-                masterCoupledFaces,
-                slaveCoupledFaces
-            );
-
-            // Generate additional coupling info (points, edges) from
-            // faces-that-match
-            faceCoupleInfo couples
-            (
-                mesh_,
-                masterCoupledFaces,
-                domainMesh,
-                slaveCoupledFaces,
-                mergeTol_,              // merge tolerance
-                true,                   // faces align
-                true,                   // couples are ordered already
-                false
-            );
-
-
-            // Add domainMesh to mesh
-            // ~~~~~~~~~~~~~~~~~~~~~~
-
-            autoPtr<mapAddedPolyMesh> map = fvMeshAdder::add
-            (
-                mesh_,
-                domainMesh,
-                couples,
-                false,          // no parallel comms
-                true            // fake complete mapping
-            );
-
-            // Update mesh data: sourceFace,sourceProc for added
-            // mesh.
-
-            sourceFace = mapBoundaryData
-            (
-                mesh_,
-                map(),
-                sourceFace,
-                domainMesh.nInternalFaces(),
-                domainSourceFace
-            );
-            sourceProc = mapBoundaryData
-            (
-                mesh_,
-                map(),
-                sourceProc,
-                domainMesh.nInternalFaces(),
-                domainSourceProc
-            );
-            sourcePatch = mapBoundaryData
-            (
-                mesh_,
-                map(),
-                sourcePatch,
-                domainMesh.nInternalFaces(),
-                domainSourcePatch
-            );
-            sourceNewNbrProc = mapBoundaryData
-            (
-                mesh_,
-                map(),
-                sourceNewNbrProc,
-                domainMesh.nInternalFaces(),
-                domainSourceNewNbrProc
-            );
-            // Update pointMaster data
-            sourcePointMaster = mapPointData
-            (
-                mesh_,
-                map(),
-                sourcePointMaster,
-                domainSourcePointMaster
-            );
-
-
-            // Update all addressing so xxProcAddressing points to correct
-            // item in masterMesh.
-            const labelList& oldCellMap = map().oldCellMap();
-            const labelList& oldFaceMap = map().oldFaceMap();
-            const labelList& oldPointMap = map().oldPointMap();
-            const labelList& oldPatchMap = map().oldPatchMap();
-
-            //Note: old mesh faces never flipped!
-            forAll(constructPatchMap, proci)
-            {
-                if (proci != sendProc && constructPatchMap[proci].size())
+                if (domainSourceFaces.set(meshi))
                 {
-                    // Processor already in mesh (either myProcNo or received)
-                    inplaceRenumber(oldCellMap, constructCellMap[proci]);
-                    inplaceRenumberWithFlip
+                    const label nIntFaces =
                     (
-                        oldFaceMap,
-                        false,
-                        true,
-                        constructFaceMap[proci]
+                        meshi == Pstream::myProcNo()
+                      ? nOldInternalFaces
+                      : meshes[meshi].nInternalFaces()
                     );
-                    inplaceRenumber(oldPointMap, constructPointMap[proci]);
-                    inplaceRenumber(oldPatchMap, constructPatchMap[proci]);
-                }
-            }
+                    const labelList& faceOwner
+                    (
+                        meshi == Pstream::myProcNo()
+                      ? oldFaceOwner
+                      : meshes[meshi].faceOwner()
+                    );
+
+                    labelList& faceMap = constructFaceMap[meshi];
+                    const labelList& cellMap = constructCellMap[meshi];
+
+                    const labelList& domainSourceFace =
+                        domainSourceFaces[meshi];
+                    const labelList& domainSourceProc =
+                        domainSourceProcs[meshi];
+                    const labelList& domainSourcePatch =
+                        domainSourcePatchs[meshi];
+                    const labelList& domainSourceNewNbr =
+                        domainSourceNewNbrProcs[meshi];
+                    UIndirectList<label>
+                    (
+                        sourcePointMaster,
+                        constructPointMap[meshi]
+                    ) = domainSourcePointMasters[meshi];
 
 
-            labelHashSet flippedAddedFaces;
-            {
-                // Find out if any faces of domain mesh were flipped (boundary
-                // faces becoming internal)
-                const label nBnd = domainMesh.nBoundaryFaces();
-                flippedAddedFaces.resize(nBnd/4);
-
-                for
-                (
-                    label domainFaceI = domainMesh.nInternalFaces();
-                    domainFaceI < domainMesh.nFaces();
-                    domainFaceI++
-                )
-                {
-                    label newFaceI = map().addedFaceMap()[domainFaceI];
-                    label newCellI = mesh_.faceOwner()[newFaceI];
-
-                    label domainCellI = domainMesh.faceOwner()[domainFaceI];
-
-                    if (newCellI != map().addedCellMap()[domainCellI])
+                    forAll(domainSourceFace, bFacei)
                     {
-                        flippedAddedFaces.insert(domainFaceI);
+                        const label oldFacei = bFacei+nIntFaces;
+                        const label allFacei = faceMap[oldFacei];
+                        const label allbFacei = allFacei-mesh_.nInternalFaces();
+
+                        if (allbFacei >= 0)
+                        {
+                            sourceProc[allbFacei] = domainSourceProc[bFacei];
+                            sourcePatch[allbFacei] = domainSourcePatch[bFacei];
+                            sourceFace[allbFacei] = domainSourceFace[bFacei];
+                            sourceNewNbrProc[allbFacei] =
+                                domainSourceNewNbr[bFacei];
+                        }
+                    }
+
+
+                    // Add flip to constructFaceMap
+                    forAll(faceMap, oldFacei)
+                    {
+                        const label allFacei = faceMap[oldFacei];
+                        const label allOwn = mesh_.faceOwner()[allFacei];
+
+                        if (cellMap[faceOwner[oldFacei]] == allOwn)
+                        {
+                            // Master face
+                            faceMap[oldFacei] += 1;
+                        }
+                        else
+                        {
+                            // Slave face. Flip.
+                            faceMap[oldFacei] = -faceMap[oldFacei] - 1;
+                        }
                     }
                 }
             }
-
-
-            // Added processor
-            inplaceRenumber(map().addedCellMap(), constructCellMap[sendProc]);
-            // Add flip
-            for (const label domainFaceI : flippedAddedFaces)
-            {
-                label& val = constructFaceMap[sendProc][domainFaceI];
-                val = -val;
-            }
-            inplaceRenumberWithFlip
-            (
-                map().addedFaceMap(),
-                false,
-                true,           // constructFaceMap has flip sign
-                constructFaceMap[sendProc]
-            );
-            inplaceRenumber(map().addedPointMap(), constructPointMap[sendProc]);
-            inplaceRenumber(map().addedPatchMap(), constructPatchMap[sendProc]);
-
-            if (debug)
-            {
-                Pout<< nl << "MERGED MESH FROM:" << sendProc << endl;
-                printMeshInfo(mesh_);
-                printFieldInfo<volScalarField>(mesh_);
-                printFieldInfo<volVectorField>(mesh_);
-                printFieldInfo<volSphericalTensorField>(mesh_);
-                printFieldInfo<volSymmTensorField>(mesh_);
-                printFieldInfo<volTensorField>(mesh_);
-                printFieldInfo<surfaceScalarField>(mesh_);
-                printFieldInfo<surfaceVectorField>(mesh_);
-                printFieldInfo<surfaceSphericalTensorField>(mesh_);
-                printFieldInfo<surfaceSymmTensorField>(mesh_);
-                printFieldInfo<surfaceTensorField>(mesh_);
-                Pout<< nl << endl;
-            }
         }
     }
 
+
     UPstream::parRun(oldParRun);  // Restore parallel state
 
 
diff --git a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.H b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.H
index edff18f8c5..3ea3774b27 100644
--- a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.H
+++ b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2018 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -42,9 +42,6 @@ Description
       expects procPatches on all)
     - initial mesh has to have procPatches last and all normal patches common
       to all processors and in the same order. This is checked.
-    - faces are matched topologically but points on the faces are not. So
-      expect problems -on separated patches (cyclics?) -on zero sized processor
-      edges.
 
 SourceFiles
     fvMeshDistribute.C
@@ -78,10 +75,6 @@ class fvMeshDistribute
         //- Underlying fvMesh
         fvMesh& mesh_;
 
-        //- Absolute merging tolerance (constructing meshes gets done using
-        //  geometric matching)
-        const scalar mergeTol_;
-
 
     // Private Member Functions
 
@@ -230,6 +223,20 @@ class fvMeshDistribute
                 labelList& slaveCoupledFaces
             );
 
+            //- Find cells on mesh whose faceID/procID match the neighbour
+            //  cell/proc of domainMesh. Store the matching face pairs
+            static void findCouples
+            (
+                const UPtrList<polyMesh>& meshes,
+                const PtrList<labelList>& domainSourceFaces,
+                const PtrList<labelList>& domainSourceProcs,
+                const PtrList<labelList>& domainSourcePatchs,
+
+                labelListList& localBoundaryFace,
+                labelListList& remoteFaceProc,
+                labelListList& remoteBoundaryFace
+            );
+
             //- Map data on boundary faces to new mesh (resulting from adding
             //  two meshes)
             static labelList mapBoundaryData
@@ -352,8 +359,8 @@ public:
 
     // Constructors
 
-        //- Construct from mesh and absolute merge tolerance
-        fvMeshDistribute(fvMesh& mesh, const scalar mergeTol);
+        //- Construct from mesh
+        fvMeshDistribute(fvMesh& mesh);
 
 
     // Member Functions
diff --git a/src/meshTools/polyTopoChange/polyTopoChange.C b/src/meshTools/polyTopoChange/polyTopoChange.C
index 9d7f32735a..675341571d 100644
--- a/src/meshTools/polyTopoChange/polyTopoChange.C
+++ b/src/meshTools/polyTopoChange/polyTopoChange.C
@@ -42,6 +42,7 @@ License
 #include "processorPolyPatch.H"
 #include "CompactListList.H"
 #include "ListOps.H"
+#include "mapPolyMesh.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -1298,36 +1299,43 @@ Foam::labelList Foam::polyTopoChange::selectFaces
 void Foam::polyTopoChange::calcPatchPointMap
 (
     const UList<Map<label>>& oldPatchMeshPointMaps,
+    const labelUList& patchMap,
     const polyBoundaryMesh& boundary,
     labelListList& patchPointMap
 ) const
 {
-    patchPointMap.setSize(boundary.size());
+    patchPointMap.setSize(patchMap.size());
 
-    forAll(boundary, patchi)
+    forAll(patchMap, patchi)
     {
-        const labelList& meshPoints = boundary[patchi].meshPoints();
+        const label oldPatchi = patchMap[patchi];
 
-        const Map<label>& oldMeshPointMap = oldPatchMeshPointMaps[patchi];
-
-        labelList& curPatchPointRnb = patchPointMap[patchi];
-
-        curPatchPointRnb.setSize(meshPoints.size());
-
-        forAll(meshPoints, i)
+        if (oldPatchi != -1)
         {
-            if (meshPoints[i] < pointMap_.size())
+            const labelList& meshPoints = boundary[patchi].meshPoints();
+
+            const Map<label>& oldMeshPointMap =
+                oldPatchMeshPointMaps[oldPatchi];
+
+            labelList& curPatchPointRnb = patchPointMap[patchi];
+
+            curPatchPointRnb.setSize(meshPoints.size());
+
+            forAll(meshPoints, i)
             {
-                // Check if old point was part of same patch
-                curPatchPointRnb[i] = oldMeshPointMap.lookup
-                (
-                    pointMap_[meshPoints[i]],
-                    -1
-                );
-            }
-            else
-            {
-                curPatchPointRnb[i] = -1;
+                if (meshPoints[i] < pointMap_.size())
+                {
+                    // Check if old point was part of same patch
+                    curPatchPointRnb[i] = oldMeshPointMap.lookup
+                    (
+                        pointMap_[meshPoints[i]],
+                        -1
+                    );
+                }
+                else
+                {
+                    curPatchPointRnb[i] = -1;
+                }
             }
         }
     }
@@ -1865,7 +1873,8 @@ void Foam::polyTopoChange::calcFaceZonePointMap
 void Foam::polyTopoChange::reorderCoupledFaces
 (
     const bool syncParallel,
-    const polyBoundaryMesh& boundary,
+    const polyBoundaryMesh& oldBoundary,
+    const labelUList& patchMap,     // new to old patches
     const labelUList& patchStarts,
     const labelUList& patchSizes,
     const pointField& points
@@ -1881,11 +1890,17 @@ void Foam::polyTopoChange::reorderCoupledFaces
     PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
 
     // Send ordering
-    forAll(boundary, patchi)
+    forAll(patchMap, patchi)
     {
-        if (syncParallel || !isA<processorPolyPatch>(boundary[patchi]))
+        const label oldPatchi = patchMap[patchi];
+
+        if
+        (
+            oldPatchi != -1
+         && (syncParallel || !isA<processorPolyPatch>(oldBoundary[oldPatchi]))
+        )
         {
-            boundary[patchi].initOrder
+            oldBoundary[oldPatchi].initOrder
             (
                 pBufs,
                 primitivePatch
@@ -1911,14 +1926,20 @@ void Foam::polyTopoChange::reorderCoupledFaces
 
     bool anyChanged = false;
 
-    forAll(boundary, patchi)
+    forAll(patchMap, patchi)
     {
-        if (syncParallel || !isA<processorPolyPatch>(boundary[patchi]))
+        const label oldPatchi = patchMap[patchi];
+
+        if
+        (
+            oldPatchi != -1
+         && (syncParallel || !isA<processorPolyPatch>(oldBoundary[oldPatchi]))
+        )
         {
             labelList patchFaceMap(patchSizes[patchi], -1);
             labelList patchFaceRotation(patchSizes[patchi], Zero);
 
-            const bool changed = boundary[patchi].order
+            const bool changed = oldBoundary[oldPatchi].order
             (
                 pBufs,
                 primitivePatch
@@ -1982,6 +2003,7 @@ void Foam::polyTopoChange::reorderCoupledFaces
 void Foam::polyTopoChange::compactAndReorder
 (
     const polyMesh& mesh,
+    const labelUList& patchMap,     // from new to old patch
     const bool syncParallel,
     const bool orderCells,
     const bool orderPoints,
@@ -2004,13 +2026,13 @@ void Foam::polyTopoChange::compactAndReorder
     List<Map<label>>& oldFaceZoneMeshPointMaps
 )
 {
-    if (mesh.boundaryMesh().size() != nPatches_)
+    if (patchMap.size() != nPatches_)
     {
         FatalErrorInFunction
             << "polyTopoChange was constructed with a mesh with "
             << nPatches_ << " patches." << endl
             << "The mesh now provided has a different number of patches "
-            << mesh.boundaryMesh().size()
+            << patchMap.size()
             << " which is illegal" << endl
             << abort(FatalError);
     }
@@ -2028,6 +2050,7 @@ void Foam::polyTopoChange::compactAndReorder
     (
         syncParallel,
         mesh.boundaryMesh(),
+        patchMap,
         patchStarts,
         patchSizes,
         newPoints
@@ -2245,7 +2268,7 @@ void Foam::polyTopoChange::addMesh
         points_.setCapacity(points_.size() + points.size());
         pointMap_.setCapacity(pointMap_.size() + points.size());
         reversePointMap_.setCapacity(reversePointMap_.size() + points.size());
-        pointZone_.resize(pointZone_.size() + points.size()/100);
+        pointZone_.resize(pointZone_.size() + points.size()/128);
 
         // Precalc offset zones
         labelList newZoneID(points.size(), -1);
@@ -2285,9 +2308,9 @@ void Foam::polyTopoChange::addMesh
 
         cellMap_.setCapacity(cellMap_.size() + nAllCells);
         reverseCellMap_.setCapacity(reverseCellMap_.size() + nAllCells);
-        cellFromPoint_.resize(cellFromPoint_.size() + nAllCells/100);
-        cellFromEdge_.resize(cellFromEdge_.size() + nAllCells/100);
-        cellFromFace_.resize(cellFromFace_.size() + nAllCells/100);
+        cellFromPoint_.resize(cellFromPoint_.size() + nAllCells/128);
+        cellFromEdge_.resize(cellFromEdge_.size() + nAllCells/128);
+        cellFromFace_.resize(cellFromFace_.size() + nAllCells/128);
         cellZone_.setCapacity(cellZone_.size() + nAllCells);
 
 
@@ -2343,10 +2366,10 @@ void Foam::polyTopoChange::addMesh
         faceNeighbour_.setCapacity(faceNeighbour_.size() + nAllFaces);
         faceMap_.setCapacity(faceMap_.size() + nAllFaces);
         reverseFaceMap_.setCapacity(reverseFaceMap_.size() + nAllFaces);
-        faceFromPoint_.resize(faceFromPoint_.size() + nAllFaces/100);
-        faceFromEdge_.resize(faceFromEdge_.size() + nAllFaces/100);
+        faceFromPoint_.resize(faceFromPoint_.size() + nAllFaces/128);
+        faceFromEdge_.resize(faceFromEdge_.size() + nAllFaces/128);
         flipFaceFlux_.setCapacity(faces_.size() + nAllFaces);
-        faceZone_.resize(faceZone_.size() + nAllFaces/100);
+        faceZone_.resize(faceZone_.size() + nAllFaces/128);
         faceZoneFlip_.setCapacity(faces_.size() + nAllFaces);
 
 
@@ -2434,7 +2457,7 @@ void Foam::polyTopoChange::setCapacity
     points_.setCapacity(nPoints);
     pointMap_.setCapacity(nPoints);
     reversePointMap_.setCapacity(nPoints);
-    pointZone_.resize(pointZone_.size() + nPoints/100);
+    pointZone_.resize(pointZone_.size() + nPoints/128);
 
     faces_.setCapacity(nFaces);
     region_.setCapacity(nFaces);
@@ -2442,17 +2465,17 @@ void Foam::polyTopoChange::setCapacity
     faceNeighbour_.setCapacity(nFaces);
     faceMap_.setCapacity(nFaces);
     reverseFaceMap_.setCapacity(nFaces);
-    faceFromPoint_.resize(faceFromPoint_.size() + nFaces/100);
-    faceFromEdge_.resize(faceFromEdge_.size() + nFaces/100);
+    faceFromPoint_.resize(faceFromPoint_.size() + nFaces/128);
+    faceFromEdge_.resize(faceFromEdge_.size() + nFaces/128);
     flipFaceFlux_.setCapacity(nFaces);
-    faceZone_.resize(faceZone_.size() + nFaces/100);
+    faceZone_.resize(faceZone_.size() + nFaces/128);
     faceZoneFlip_.setCapacity(nFaces);
 
     cellMap_.setCapacity(nCells);
     reverseCellMap_.setCapacity(nCells);
-    cellFromPoint_.resize(cellFromPoint_.size() + nCells/100);
-    cellFromEdge_.resize(cellFromEdge_.size() + nCells/100);
-    cellFromFace_.resize(cellFromFace_.size() + nCells/100);
+    cellFromPoint_.resize(cellFromPoint_.size() + nCells/128);
+    cellFromEdge_.resize(cellFromEdge_.size() + nCells/128);
+    cellFromFace_.resize(cellFromFace_.size() + nCells/128);
     cellZone_.setCapacity(nCells);
 }
 
@@ -2958,6 +2981,7 @@ void Foam::polyTopoChange::removeCell
 Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::changeMesh
 (
     polyMesh& mesh,
+    const labelUList& patchMap,
     const bool inflate,
     const bool syncParallel,
     const bool orderCells,
@@ -3003,6 +3027,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::changeMesh
     compactAndReorder
     (
         mesh,
+        patchMap,
         syncParallel,
         orderCells,
         orderPoints,
@@ -3159,6 +3184,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::changeMesh
     calcPatchPointMap
     (
         oldPatchMeshPointMaps,
+        patchMap,
         mesh.boundaryMesh(),
         patchPointMap
     );
@@ -3218,4 +3244,25 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::changeMesh
 }
 
 
+Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::changeMesh
+(
+    polyMesh& mesh,
+    const bool inflate,
+    const bool syncParallel,
+    const bool orderCells,
+    const bool orderPoints
+)
+{
+    return changeMesh
+    (
+        mesh,
+        identity(mesh.boundaryMesh().size()),
+        inflate,
+        syncParallel,
+        orderCells,
+        orderPoints
+    );
+}
+
+
 // ************************************************************************* //
diff --git a/src/meshTools/polyTopoChange/polyTopoChange.H b/src/meshTools/polyTopoChange/polyTopoChange.H
index 00fd3603b3..d7661e4c43 100644
--- a/src/meshTools/polyTopoChange/polyTopoChange.H
+++ b/src/meshTools/polyTopoChange/polyTopoChange.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -71,7 +71,6 @@ SourceFiles
 #include "pointField.H"
 #include "Map.H"
 #include "HashSet.H"
-#include "mapPolyMesh.H"
 #include "bitSet.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -91,6 +90,7 @@ class dictionary;
 class topoAction;
 class objectMap;
 class IOobject;
+class mapPolyMesh;
 template<class T, class Container> class CompactListList;
 
 /*---------------------------------------------------------------------------*\
@@ -349,6 +349,7 @@ class polyTopoChange
         void calcPatchPointMap
         (
             const UList<Map<label>>& oldPatchMeshPointMaps,
+            const labelUList& patchMap,
             const polyBoundaryMesh& boundary,
             labelListList& patchPointMap
         ) const;
@@ -394,6 +395,7 @@ class polyTopoChange
             (
                 const bool syncParallel,
                 const polyBoundaryMesh& boundary,
+                const labelUList& patchMap,     // new back to old patch labels
                 const labelUList& patchStarts,
                 const labelUList& patchSizes,
                 const pointField& points
@@ -402,6 +404,7 @@ class polyTopoChange
         void compactAndReorder
         (
             const polyMesh& mesh,
+            const labelUList& patchMap,      // from new to old patch
             const bool syncParallel,
             const bool orderCells,
             const bool orderPoints,
@@ -611,12 +614,37 @@ public:
             autoPtr<mapPolyMesh> changeMesh
             (
                 polyMesh& mesh,
+                const labelUList& patchMap,
                 const bool inflate,
                 const bool syncParallel = true,
                 const bool orderCells = false,
                 const bool orderPoints = false
             );
 
+            //- Without patch mapping
+            autoPtr<mapPolyMesh> changeMesh
+            (
+                polyMesh& mesh,
+                const bool inflate,
+                const bool syncParallel = true,
+                const bool orderCells = false,
+                const bool orderPoints = false
+            );
+
+            //- Create new mesh with old mesh patches. Additional dictionaries
+            //  (fv* etc) read according to IO flags
+            template<class Type>
+            autoPtr<mapPolyMesh> makeMesh
+            (
+                autoPtr<Type>& newMesh,
+                const IOobject& io,
+                const polyMesh& mesh,
+                const labelUList& patchMap,
+                const bool syncParallel = true,
+                const bool orderCells = false,
+                const bool orderPoints = false
+            );
+
             //- Create new mesh with old mesh patches. Additional dictionaries
             //  (fv* etc) read according to IO flags
             template<class Type>
diff --git a/src/meshTools/polyTopoChange/polyTopoChangeTemplates.C b/src/meshTools/polyTopoChange/polyTopoChangeTemplates.C
index 8f02476777..eea465b276 100644
--- a/src/meshTools/polyTopoChange/polyTopoChangeTemplates.C
+++ b/src/meshTools/polyTopoChange/polyTopoChangeTemplates.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2017-2020 OpenCFD Ltd.
+    Copyright (C) 2017-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -26,9 +26,10 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "polyTopoChange.H"
 #include "polyMesh.H"
 #include "HashOps.H"
+#include "emptyPolyPatch.H"
+#include "mapPolyMesh.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
@@ -105,6 +106,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::makeMesh
     autoPtr<Type>& newMeshPtr,
     const IOobject& io,
     const polyMesh& mesh,
+    const labelUList& patchMap,
     const bool syncParallel,
     const bool orderCells,
     const bool orderPoints
@@ -151,6 +153,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::makeMesh
     compactAndReorder
     (
         mesh,
+        patchMap,      // from new to old patch
         syncParallel,
         orderCells,
         orderPoints,
@@ -239,17 +242,35 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::makeMesh
     {
         const polyBoundaryMesh& oldPatches = mesh.boundaryMesh();
 
-        List<polyPatch*> newBoundary(oldPatches.size());
+        List<polyPatch*> newBoundary(patchMap.size());
 
-        forAll(oldPatches, patchi)
+        forAll(patchMap, patchi)
         {
-            newBoundary[patchi] = oldPatches[patchi].clone
-            (
-                newMesh.boundaryMesh(),
-                patchi,
-                patchSizes[patchi],
-                patchStarts[patchi]
-            ).ptr();
+            const label oldPatchi = patchMap[patchi];
+
+            if (oldPatchi != -1)
+            {
+                newBoundary[patchi] = oldPatches[oldPatchi].clone
+                (
+                    newMesh.boundaryMesh(),
+                    patchi,
+                    patchSizes[patchi],
+                    patchStarts[patchi]
+                ).ptr();
+            }
+            else
+            {
+                // Added patch
+                newBoundary[patchi] = new emptyPolyPatch
+                (
+                    "patch" + Foam::name(patchi),
+                    patchSizes[patchi],
+                    patchStarts[patchi],
+                    patchi,
+                    newMesh.boundaryMesh(),
+                    word::null
+                );
+            }
         }
         newMesh.addFvPatches(newBoundary);
     }
@@ -327,6 +348,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::makeMesh
     calcPatchPointMap
     (
         oldPatchMeshPointMaps,
+        patchMap,
         newMesh.boundaryMesh(),
         patchPointMap
     );
@@ -338,7 +360,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::makeMesh
     if (debug)
     {
         Pout<< "New mesh:" << nl;
-        writeMeshStats(mesh, Pout);
+        writeMeshStats(newMesh, Pout);
     }
 
     labelHashSet flipFaceFluxSet(HashSetOps::used(flipFaceFlux_));
@@ -390,4 +412,28 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::makeMesh
 }
 
 
+template<class Type>
+Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::makeMesh
+(
+    autoPtr<Type>& newMeshPtr,
+    const IOobject& io,
+    const polyMesh& mesh,
+    const bool syncParallel,
+    const bool orderCells,
+    const bool orderPoints
+)
+{
+    return makeMesh
+    (
+        newMeshPtr,
+        io,
+        mesh,
+        identity(mesh.boundaryMesh().size()),
+        syncParallel,
+        orderCells,
+        orderPoints
+    );
+}
+
+
 // ************************************************************************* //