diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
index 7f74c3761d..9349455847 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -405,7 +405,7 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
                 icellWeights
             );
 
-            fvMeshDistribute distributor(mesh_, mergeDist_);
+            fvMeshDistribute distributor(mesh_);
 
             autoPtr<mapDistributePolyMesh> mapDist = distributor.distribute
             (
@@ -813,7 +813,6 @@ Foam::backgroundMeshDecomposition::backgroundMeshDecomposition
         )
     ),
     decomposerPtr_(decompositionMethod::New(decomposeDict_)),
-    mergeDist_(1e-6*mesh_.bounds().mag()),
     spanScale_(coeffsDict.lookup<scalar>("spanScale")),
     minCellSizeLimit_
     (
@@ -998,7 +997,7 @@ Foam::backgroundMeshDecomposition::distribute
 
     Info<< "    Redistributing background mesh cells" << endl;
 
-    fvMeshDistribute distributor(mesh_, mergeDist_);
+    fvMeshDistribute distributor(mesh_);
 
     autoPtr<mapDistributePolyMesh> mapDist = distributor.distribute(newDecomp);
 
diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
index b7264a2f64..ae7edf54ce 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -1247,7 +1247,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 0432e048e6..080717687a 100644
--- a/applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C
+++ b/applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,17 +25,11 @@ Application
     reconstructParMesh
 
 Description
-    Reconstructs a mesh using geometric information only.
+    Reconstructs a mesh.
 
     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)
-
-    If the parallel case does not have correct procBoundaries use the
-    -fullMatch option which will check all boundary faces (bit slower).
-
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
@@ -56,320 +50,111 @@ 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-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
-// all the processors < proci.
 autoPtr<faceCoupleInfo> determineCoupledFaces
 (
-    const bool fullMatch,
     const label masterMeshProcStart,
     const label masterMeshProcEnd,
     const polyMesh& masterMesh,
     const label meshToAddProcStart,
     const label meshToAddProcEnd,
-    const polyMesh& meshToAdd,
-    const scalar mergeDist
+    const polyMesh& meshToAdd
 )
 {
-    if (fullMatch || masterMesh.nCells() == 0)
-    {
-        return autoPtr<faceCoupleInfo>
-        (
-            new faceCoupleInfo
-            (
-                masterMesh,
-                meshToAdd,
-                mergeDist,      // Absolute merging distance
-                true            // Matching faces identical
-            )
-        );
-    }
-    else
-    {
-        // Pick up all patches on masterMesh ending in "toDDD" where DDD is
-        // the processor number proci.
+    const polyBoundaryMesh& masterPatches = masterMesh.boundaryMesh();
+    const polyBoundaryMesh& addPatches = meshToAdd.boundaryMesh();
 
-        const polyBoundaryMesh& masterPatches = masterMesh.boundaryMesh();
-
-
-        DynamicList<label> masterFaces
-        (
-            masterMesh.nFaces()
-          - masterMesh.nInternalFaces()
-        );
-
-
-        forAll(masterPatches, patchi)
-        {
-            const polyPatch& pp = masterPatches[patchi];
-
-            if (isA<processorPolyPatch>(pp))
-            {
-                for
-                (
-                    label proci=meshToAddProcStart;
-                    proci<meshToAddProcEnd;
-                    proci++
-                )
-                {
-                    const string toProcString("to" + name(proci));
-                    if (
-                        pp.name().rfind(toProcString)
-                     == (pp.name().size()-toProcString.size())
-                    )
-                    {
-                        label meshFacei = pp.start();
-                        forAll(pp, i)
-                        {
-                            masterFaces.append(meshFacei++);
-                        }
-                        break;
-                    }
-                }
-
-            }
-        }
-        masterFaces.shrink();
-
-
-        // Pick up all patches on meshToAdd ending in "procBoundaryDDDtoYYY"
-        // where DDD is the processor number proci and YYY is < proci.
-
-        const polyBoundaryMesh& addPatches = meshToAdd.boundaryMesh();
-
-        DynamicList<label> addFaces
-        (
-            meshToAdd.nFaces()
-          - meshToAdd.nInternalFaces()
-        );
-
-        forAll(addPatches, patchi)
-        {
-            const polyPatch& pp = addPatches[patchi];
-
-            if (isA<processorPolyPatch>(pp))
-            {
-                bool isConnected = false;
-
-                for
-                (
-                    label mergedProci=masterMeshProcStart;
-                    !isConnected && (mergedProci < masterMeshProcEnd);
-                    mergedProci++
-                )
-                {
-                    for
-                    (
-                        label proci = meshToAddProcStart;
-                        proci < meshToAddProcEnd;
-                        proci++
-                    )
-                    {
-                        const word fromProcString
-                        (
-                            processorPolyPatch::newName(proci, mergedProci)
-                        );
-
-                        if (pp.name() == fromProcString)
-                        {
-                            isConnected = true;
-                            break;
-                        }
-                    }
-                }
-
-                if (isConnected)
-                {
-                    label meshFacei = pp.start();
-                    forAll(pp, i)
-                    {
-                        addFaces.append(meshFacei++);
-                    }
-                }
-            }
-        }
-        addFaces.shrink();
-
-        return autoPtr<faceCoupleInfo>
-        (
-            new faceCoupleInfo
-            (
-                masterMesh,
-                masterFaces,
-                meshToAdd,
-                addFaces,
-                mergeDist,      // Absolute merging distance
-                true,           // Matching faces identical?
-                false,          // If perfect match are faces already ordered
-                                // (e.g. processor patches)
-                false           // are faces each on separate patch?
-            )
-        );
-    }
-}
-
-
-autoPtr<mapPolyMesh> mergeSharedPoints
-(
-    const scalar mergeDist,
-    polyMesh& mesh,
-    labelListList& pointProcAddressing
-)
-{
-    // Find out which sets of points get merged and create a map from
-    // mesh point to unique point.
-    Map<label> pointToMaster
+    DynamicList<label> masterFaces
     (
-        fvMeshAdder::findSharedPoints
-        (
-            mesh,
-            mergeDist
-        )
+        masterMesh.nFaces() - masterMesh.nInternalFaces()
+    );
+    DynamicList<label> addFaces
+    (
+        meshToAdd.nFaces() - meshToAdd.nInternalFaces()
     );
 
-    Info<< "mergeSharedPoints : detected " << pointToMaster.size()
-        << " points that are to be merged." << endl;
-
-    if (returnReduce(pointToMaster.size(), sumOp<label>()) == 0)
+    for
+    (
+        label masterProci = masterMeshProcStart;
+        masterProci < masterMeshProcEnd;
+        masterProci++
+    )
     {
-        return autoPtr<mapPolyMesh>(nullptr);
-    }
-
-    polyTopoChange meshMod(mesh);
-
-    fvMeshAdder::mergePoints(mesh, pointToMaster, meshMod);
-
-    // Change the mesh (no inflation). Note: parallel comms allowed.
-    autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false, true);
-
-    // Update fields. No inflation, parallel sync.
-    mesh.updateMesh(map);
-
-    // pointProcAddressing give indices into the master mesh so adapt them
-    // for changed point numbering.
-
-    // Adapt constructMaps for merged points.
-    forAll(pointProcAddressing, proci)
-    {
-        labelList& constructMap = pointProcAddressing[proci];
-
-        forAll(constructMap, i)
+        for
+        (
+            label addProci = meshToAddProcStart;
+            addProci < meshToAddProcEnd;
+            addProci++
+        )
         {
-            label oldPointi = constructMap[i];
+            const word masterToAddName
+            (
+                "procBoundary" + name(masterProci) + "to" + name(addProci)
+            );
+            const word addToMasterName
+            (
+                "procBoundary" + name(addProci) + "to" + name(masterProci)
+            );
 
-            // New label of point after changeMesh.
-            label newPointi = map().reversePointMap()[oldPointi];
+            const label masterToAddID =
+                masterPatches.findPatchID(masterToAddName);
+            const label addToMasterID =
+                addPatches.findPatchID(addToMasterName);
 
-            if (newPointi < -1)
+            if (masterToAddID != -1 && addToMasterID != -1)
             {
-                constructMap[i] = -newPointi-2;
+                const polyPatch& masterPp = masterPatches[masterToAddID];
+
+                forAll(masterPp, i)
+                {
+                    masterFaces.append(masterPp.start() + i);
+                }
+
+                const polyPatch& addPp = addPatches[addToMasterID];
+
+                forAll(addPp, i)
+                {
+                    addFaces.append(addPp.start() + i);
+                }
             }
-            else if (newPointi >= 0)
-            {
-                constructMap[i] = newPointi;
-            }
-            else
+
+            if ((masterToAddID != -1) != (addToMasterID != -1))
             {
+                const label foundProci =
+                    masterToAddID != -1 ? masterProci : addProci;
+                const word& foundName =
+                    masterToAddID != -1 ? masterToAddName : addToMasterName;
+
+                const label missingProci =
+                    masterToAddID != -1 ? addProci : masterProci;
+                const word& missingName =
+                    masterToAddID != -1 ? addToMasterName : masterToAddName;
+
                 FatalErrorInFunction
-                    << "Problem. oldPointi:" << oldPointi
-                    << " newPointi:" << newPointi << abort(FatalError);
+                    << "Patch " << foundName << " found on processor "
+                    << foundProci << " but corresponding patch "
+                    << missingName << " missing on processor "
+                    << missingProci << exit(FatalError);
             }
         }
     }
 
-    return map;
+    masterFaces.shrink();
+    addFaces.shrink();
+
+    return autoPtr<faceCoupleInfo>
+    (
+        new faceCoupleInfo
+        (
+            masterMesh,
+            masterFaces,
+            meshToAdd,
+            addFaces
+        )
+    );
 }
 
 
-boundBox procBounds
-(
-    const argList& args,
-    const PtrList<Time>& databases,
-    const word& regionDir
-)
-{
-    boundBox bb = boundBox::invertedBox;
-
-    forAll(databases, proci)
-    {
-        fileName pointsInstance
-        (
-            databases[proci].findInstance
-            (
-                regionDir/polyMesh::meshSubDir,
-                "points"
-            )
-        );
-
-        if (pointsInstance != databases[proci].timeName())
-        {
-            FatalErrorInFunction
-                << "Your time was specified as " << databases[proci].timeName()
-                << " but there is no polyMesh/points in that time." << endl
-                << "(there is a points file in " << pointsInstance
-                << ")" << endl
-                << "Please rerun with the correct time specified"
-                << " (through the -constant, -time or -latestTime "
-                << "(at your option)."
-                << endl << exit(FatalError);
-        }
-
-        Info<< "Reading points from "
-            << databases[proci].caseName()
-            << " for time = " << databases[proci].timeName()
-            << nl << endl;
-
-        pointIOField points
-        (
-            IOobject
-            (
-                "points",
-                databases[proci].findInstance
-                (
-                    regionDir/polyMesh::meshSubDir,
-                    "points"
-                ),
-                regionDir/polyMesh::meshSubDir,
-                databases[proci],
-                IOobject::MUST_READ,
-                IOobject::NO_WRITE,
-                false
-            )
-        );
-
-        boundBox domainBb(points, false);
-
-        bb.min() = min(bb.min(), domainBb.min());
-        bb.max() = max(bb.max(), domainBb.max());
-    }
-
-    return bb;
-}
-
-
-void writeCellDistance
+void writeCellDistribution
 (
     Time& runTime,
     const fvMesh& masterMesh,
@@ -451,27 +236,9 @@ void writeCellDistance
 
 int main(int argc, char *argv[])
 {
-    argList::addNote
-    (
-        "reconstruct a mesh using geometric information only"
-    );
-
-    // Enable -constant ... if someone really wants it
-    // Enable -withZero to prevent accidentally trashing the initial fields
+    argList::addNote("reconstruct a mesh");
     timeSelector::addOptions(true, true);
     argList::noParallel();
-    argList::addOption
-    (
-        "mergeTol",
-        "scalar",
-        "specify the merge distance relative to the bounding box size "
-        "(default 1e-7)"
-    );
-    argList::addBoolOption
-    (
-        "fullMatch",
-        "do (slower) geometric matching on all boundary faces"
-    );
     argList::addBoolOption
     (
         "cellDist",
@@ -484,22 +251,6 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
 
-    Info<< "This is an experimental tool which tries to merge"
-        << " individual processor" << nl
-        << "meshes back into one master mesh. Use it if the original"
-        << " master mesh has" << nl
-        << "been deleted or if the processor meshes have been modified"
-        << " (topology change)." << nl
-        << "This tool will write the resulting mesh to a new time step"
-        << " and construct" << nl
-        << "xxxxProcAddressing files in the processor meshes so"
-        << " reconstructPar can be" << nl
-        << "used to regenerate the fields on the master mesh." << nl
-        << nl
-        << "Not well tested & use at your own risk!" << nl
-        << endl;
-
-
     const wordList regionNames(selectRegionNames(args, runTime));
     if (regionNames.size() > 1)
     {
@@ -515,52 +266,11 @@ int main(int argc, char *argv[])
         Info<< "Operating on region " << regionNames[0] << nl << endl;
     }
 
-
-    scalar mergeTol = defaultMergeTol;
-    args.optionReadIfPresent("mergeTol", mergeTol);
-
-    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.optionFound("fullMatch");
-
-    if (fullMatch)
-    {
-        Info<< "Doing geometric matching on all boundary faces." << nl << endl;
-    }
-    else
-    {
-        Info<< "Doing geometric matching on correct procBoundaries only."
-            << nl << "This assumes a correct decomposition." << endl;
-    }
-
-    bool writeCellDist = args.optionFound("cellDist");
-
-
     label nProcs = fileHandler().nProcs(args.path());
-
     Info<< "Found " << nProcs << " processor directories" << nl << endl;
 
-
     // Read all time databases
     PtrList<Time> databases(nProcs);
-
     forAll(databases, proci)
     {
         Info<< "Reading database "
@@ -630,18 +340,6 @@ 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);
@@ -656,9 +354,9 @@ int main(int argc, char *argv[])
 
             {
                 // Construct empty mesh.
-                // fvMesh** masterMesh = new fvMesh*[nProcs];
                 PtrList<fvMesh> masterMesh(nProcs);
 
+                // Read all the meshes
                 for (label proci=0; proci<nProcs; proci++)
                 {
                     masterMesh.set
@@ -696,17 +394,15 @@ int main(int argc, char *argv[])
                     boundaryProcAddressing[proci] =
                         identity(meshToAdd.boundaryMesh().size());
 
-                    // Find geometrically shared points/faces.
+                    // Find shared points/faces
                     autoPtr<faceCoupleInfo> couples = determineCoupledFaces
                     (
-                        fullMatch,
                         proci,
                         proci,
                         masterMesh[proci],
                         proci,
                         proci,
-                        meshToAdd,
-                        mergeDist
+                        meshToAdd
                     );
 
                     // Add elements to mesh
@@ -718,15 +414,29 @@ int main(int argc, char *argv[])
                     );
 
                     // Added processor
-                    renumber(map().addedCellMap(), cellProcAddressing[proci]);
-                    renumber(map().addedFaceMap(), faceProcAddressing[proci]);
-                    renumber(map().addedPointMap(), pointProcAddressing[proci]);
-                    renumber
+                    inplaceRenumber
+                    (
+                        map().addedCellMap(),
+                        cellProcAddressing[proci]
+                    );
+                    inplaceRenumber
+                    (
+                        map().addedFaceMap(),
+                        faceProcAddressing[proci]
+                    );
+                    inplaceRenumber
+                    (
+                        map().addedPointMap(),
+                        pointProcAddressing[proci]
+                    );
+                    inplaceRenumber
                     (
                         map().addedPatchMap(),
                         boundaryProcAddressing[proci]
                     );
                 }
+
+                // Merge the meshes
                 for (label step=2; step<nProcs*2; step*=2)
                 {
                     for (label proci=0; proci<nProcs; proci+=step)
@@ -740,17 +450,15 @@ int main(int argc, char *argv[])
                         Info<< "Merging mesh " << proci << " with " << next
                             << endl;
 
-                        // Find geometrically shared points/faces.
+                        // Find shared points/faces
                         autoPtr<faceCoupleInfo> couples = determineCoupledFaces
                         (
-                            fullMatch,
                             proci,
                             next,
                             masterMesh[proci],
                             next,
                             proci+step,
-                            masterMesh[next],
-                            mergeDist
+                            masterMesh[next]
                         );
 
                         // Add elements to mesh
@@ -764,26 +472,22 @@ int main(int argc, char *argv[])
                         // Processors that were already in masterMesh
                         for (label mergedI=proci; mergedI<next; mergedI++)
                         {
-                            renumber
+                            inplaceRenumber
                             (
                                 map().oldCellMap(),
                                 cellProcAddressing[mergedI]
                             );
-
-                            renumber
+                            inplaceRenumber
                             (
                                 map().oldFaceMap(),
                                 faceProcAddressing[mergedI]
                             );
-
-                            renumber
+                            inplaceRenumber
                             (
                                 map().oldPointMap(),
                                 pointProcAddressing[mergedI]
                             );
-
-                            // Note: boundary is special since can contain -1.
-                            renumber
+                            inplaceRenumber
                             (
                                 map().oldPatchMap(),
                                 boundaryProcAddressing[mergedI]
@@ -798,25 +502,22 @@ int main(int argc, char *argv[])
                             addedI++
                         )
                         {
-                            renumber
+                            inplaceRenumber
                             (
                                 map().addedCellMap(),
                                 cellProcAddressing[addedI]
                             );
-
-                            renumber
+                            inplaceRenumber
                             (
                                 map().addedFaceMap(),
                                 faceProcAddressing[addedI]
                             );
-
-                            renumber
+                            inplaceRenumber
                             (
                                 map().addedPointMap(),
                                 pointProcAddressing[addedI]
                             );
-
-                            renumber
+                            inplaceRenumber
                             (
                                 map().addedPatchMap(),
                                 boundaryProcAddressing[addedI]
@@ -835,14 +536,6 @@ int main(int argc, char *argv[])
                         << nl << nl << endl;
                 }
 
-                // See if any points on the mastermesh have become connected
-                // because of connections through processor meshes.
-                mergeSharedPoints
-                (
-                    mergeDist,
-                    masterMesh[0],
-                    pointProcAddressing
-                );
 
                 // Save some properties on the reconstructed mesh
                 masterInternalFaces = masterMesh[0].nInternalFaces();
@@ -860,9 +553,9 @@ int main(int argc, char *argv[])
                         << exit(FatalError);
                 }
 
-                if (writeCellDist)
+                if (args.optionFound("cellDist"))
                 {
-                    writeCellDistance
+                    writeCellDistribution
                     (
                         runTime,
                         masterMesh[0],
diff --git a/applications/utilities/parallelProcessing/redistributePar/redistributePar.C b/applications/utilities/parallelProcessing/redistributePar/redistributePar.C
index 73b9f6a4d1..562da7029a 100644
--- a/applications/utilities/parallelProcessing/redistributePar/redistributePar.C
+++ b/applications/utilities/parallelProcessing/redistributePar/redistributePar.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -58,65 +58,6 @@ Description
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// 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
-)
-{
-    scalar mergeTol = defaultMergeTol;
-    args.optionReadIfPresent("mergeTol", mergeTol);
-
-    scalar writeTol =
-        Foam::pow(scalar(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);
-    }
-
-    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 printMeshData(Ostream& os, const polyMesh& mesh)
-//{
-//    os  << "Number of points:           " << mesh.points().size() << nl
-//        << "          faces:            " << mesh.faces().size() << nl
-//        << "          internal faces:   " << mesh.faceNeighbour().size() << nl
-//        << "          cells:            " << mesh.cells().size() << nl
-//        << "          boundary patches: " << mesh.boundaryMesh().size() << nl
-//        << "          point zones:      " << mesh.pointZones().size() << nl
-//        << "          face zones:       " << mesh.faceZones().size() << nl
-//        << "          cell zones:       " << mesh.cellZones().size() << nl;
-//}
-
-
 void printMeshData(const polyMesh& mesh)
 {
     // Collect all data on master
@@ -365,72 +306,11 @@ void readFields
 }
 
 
-// Debugging: compare two fields.
-void compareFields
-(
-    const scalar tolDim,
-    const volVectorField& a,
-    const volVectorField& b
-)
-{
-    forAll(a, celli)
-    {
-        if (mag(b[celli] - a[celli]) > tolDim)
-        {
-            FatalErrorInFunction
-                << "Did not map volVectorField correctly:" << nl
-                << "cell:" << celli
-                << " transfer b:" << b[celli]
-                << " real cc:" << a[celli]
-                << abort(FatalError);
-        }
-    }
-    forAll(a.boundaryField(), patchi)
-    {
-        // We have real mesh cellcentre and
-        // mapped original cell centre.
-
-        const fvPatchVectorField& aBoundary =
-            a.boundaryField()[patchi];
-
-        const fvPatchVectorField& bBoundary =
-            b.boundaryField()[patchi];
-
-        if (!bBoundary.coupled())
-        {
-            forAll(aBoundary, i)
-            {
-                if (mag(aBoundary[i] - bBoundary[i]) > tolDim)
-                {
-                    WarningInFunction
-                        << "Did not map volVectorField correctly:"
-                        << endl
-                        << "patch:" << patchi << " patchFace:" << i
-                        << " cc:" << endl
-                        << "    real    :" << aBoundary[i] << endl
-                        << "    mapped  :" << bBoundary[i] << endl
-                        << "This might be just a precision entry"
-                        << " on writing the mesh." << endl;
-                        //<< abort(FatalError);
-                }
-            }
-        }
-    }
-}
-
-
-
 int main(int argc, char *argv[])
 {
     #include "addRegionOption.H"
     #include "addOverwriteOption.H"
-    argList::addOption
-    (
-        "mergeTol",
-        "scalar",
-        "specify the merge distance relative to the bounding box size "
-        "(default 1e-6)"
-    );
+
     // Include explicit constant options, have zero from time range
     timeSelector::addOptions();
 
@@ -781,16 +661,8 @@ int main(int argc, char *argv[])
     // Used to test correctness of mapping
     // volVectorField mapCc("mapCc", 1*mesh.C());
 
-    // Global matching tolerance
-    const scalar tolDim = getMergeDistance
-    (
-        args,
-        runTime,
-        mesh.bounds()
-    );
-
     // Mesh distribution engine
-    fvMeshDistribute distributor(mesh, tolDim);
+    fvMeshDistribute distributor(mesh);
 
     // Pout<< "Wanted distribution:"
     //    << distributor.countCells(finalDecomp) << nl << endl;
@@ -819,9 +691,6 @@ int main(int argc, char *argv[])
     mesh.write();
 
 
-    // Debugging: test mapped cellcentre field.
-    // compareFields(tolDim, mesh.C(), mapCc);
-
     // Print nice message
     // ~~~~~~~~~~~~~~~~~~
 
diff --git a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C
index 7ae0d3e512..79691713e3 100644
--- a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C
+++ b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -656,10 +656,95 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::mergeSharedPoints
         return autoPtr<mapPolyMesh>(nullptr);
     }
 
-
+    // Create the mesh change engine to merge the points
     polyTopoChange meshMod(mesh_);
+    {
+        // Remove all non-master points.
+        forAll(mesh_.points(), pointi)
+        {
+            Map<label>::const_iterator iter = pointToMaster.find(pointi);
 
-    fvMeshAdder::mergePoints(mesh_, pointToMaster, meshMod);
+            if (iter != pointToMaster.end())
+            {
+                if (iter() != pointi)
+                {
+                    meshMod.removePoint(pointi, iter());
+                }
+            }
+        }
+
+        // Modify faces for points. Note: could use pointFaces here but want to
+        // avoid addressing calculation.
+        const faceList& faces = mesh_.faces();
+
+        forAll(faces, facei)
+        {
+            const face& f = faces[facei];
+
+            bool hasMerged = false;
+
+            forAll(f, fp)
+            {
+                label pointi = f[fp];
+
+                Map<label>::const_iterator iter = pointToMaster.find(pointi);
+
+                if (iter != pointToMaster.end())
+                {
+                    if (iter() != pointi)
+                    {
+                        hasMerged = true;
+                        break;
+                    }
+                }
+            }
+
+            if (hasMerged)
+            {
+                face newF(f);
+
+                forAll(f, fp)
+                {
+                    label pointi = f[fp];
+
+                    Map<label>::const_iterator iter =
+                        pointToMaster.find(pointi);
+
+                    if (iter != pointToMaster.end())
+                    {
+                        newF[fp] = iter();
+                    }
+                }
+
+                label patchID = mesh_.boundaryMesh().whichPatch(facei);
+                label nei = (patchID == -1 ? mesh_.faceNeighbour()[facei] : -1);
+                label zoneID = mesh_.faceZones().whichZone(facei);
+                bool zoneFlip = false;
+
+                if (zoneID >= 0)
+                {
+                    const faceZone& fZone = mesh_.faceZones()[zoneID];
+                    zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
+                }
+
+                meshMod.setAction
+                (
+                    polyModifyFace
+                    (
+                        newF,                       // modified face
+                        facei,                      // label of face
+                        mesh_.faceOwner()[facei],   // owner
+                        nei,                        // neighbour
+                        false,                      // face flip
+                        patchID,                    // patch for face
+                        false,                      // remove from zone
+                        zoneID,                     // zone for face
+                        zoneFlip                    // face flip in zone
+                    )
+                );
+            }
+        }
+    }
 
     // Change the mesh (no inflation). Note: parallel comms allowed.
     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
@@ -1624,10 +1709,9 @@ Foam::autoPtr<Foam::fvMesh> Foam::fvMeshDistribute::receiveMesh
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from components
-Foam::fvMeshDistribute::fvMeshDistribute(fvMesh& mesh, const scalar mergeTol)
+Foam::fvMeshDistribute::fvMeshDistribute(fvMesh& mesh)
 :
-    mesh_(mesh),
-    mergeTol_(mergeTol)
+    mesh_(mesh)
 {}
 
 
@@ -2676,11 +2760,7 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
                 mesh_,
                 masterCoupledFaces,
                 domainMesh,
-                slaveCoupledFaces,
-                mergeTol_,              // merge tolerance
-                true,                   // faces align
-                true,                   // couples are ordered already
-                false
+                slaveCoupledFaces
             );
 
 
diff --git a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.H b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.H
index a50cf3f4f1..ead216fa2a 100644
--- a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.H
+++ b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -75,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
 
@@ -344,8 +340,8 @@ public:
 
     // Constructors
 
-        //- Construct from mesh and absolute merge tolerance
-        fvMeshDistribute(fvMesh& mesh, const scalar mergeTol);
+        //- Construct from mesh
+        fvMeshDistribute(fvMesh& mesh);
 
         //- Disallow default bitwise copy construction
         fvMeshDistribute(const fvMeshDistribute&) = delete;
diff --git a/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C b/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C
index 19e12133cc..1deeb657a8 100644
--- a/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C
+++ b/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,7 +37,6 @@ License
 namespace Foam
 {
     defineTypeNameAndDebug(faceCoupleInfo, 0);
-    const scalar faceCoupleInfo::angleTol_ = 1e-3;
 }
 
 
@@ -47,44 +46,31 @@ void Foam::faceCoupleInfo::writeOBJ
 (
     const fileName& fName,
     const edgeList& edges,
-    const pointField& points,
-    const bool compact
+    const pointField& points
 )
 {
     OFstream str(fName);
 
     labelList pointMap(points.size(), -1);
 
-    if (compact)
+    label newPointi = 0;
+
+    forAll(edges, edgeI)
     {
-        label newPointi = 0;
+        const edge& e = edges[edgeI];
 
-        forAll(edges, edgeI)
+        forAll(e, eI)
         {
-            const edge& e = edges[edgeI];
+            label pointi = e[eI];
 
-            forAll(e, eI)
+            if (pointMap[pointi] == -1)
             {
-                label pointi = e[eI];
+                pointMap[pointi] = newPointi++;
 
-                if (pointMap[pointi] == -1)
-                {
-                    pointMap[pointi] = newPointi++;
-
-                    meshTools::writeOBJ(str, points[pointi]);
-                }
+                meshTools::writeOBJ(str, points[pointi]);
             }
         }
     }
-    else
-    {
-        forAll(points, pointi)
-        {
-            meshTools::writeOBJ(str, points[pointi]);
-        }
-
-        pointMap = identity(points.size());
-    }
 
     forAll(edges, edgeI)
     {
@@ -123,7 +109,6 @@ void Foam::faceCoupleInfo::writePointsFaces() const
 {
     const indirectPrimitivePatch& m = masterPatch();
     const indirectPrimitivePatch& s = slavePatch();
-    const primitiveFacePatch& c = cutFaces();
 
     // Patches
     {
@@ -136,83 +121,40 @@ void Foam::faceCoupleInfo::writePointsFaces() const
         Pout<< "Writing slavePatch to " << str.name() << endl;
         meshTools::writeOBJ(str, s.localFaces(), s.localPoints());
     }
-    {
-        OFstream str("cutFaces.obj");
-        Pout<< "Writing cutFaces to " << str.name() << endl;
-        meshTools::writeOBJ(str, c.localFaces(), c.localPoints());
-    }
 
     // Point connectivity
     {
-        Pout<< "Writing cutToMasterPoints to cutToMasterPoints.obj" << endl;
+        Pout<< "Writing masterToSlavePoints to masterToSlavePoints.obj" << endl;
+
+        const labelListList coupleToMasterPoints(this->coupleToMasterPoints());
+        const labelListList coupleToSlavePoints(this->coupleToSlavePoints());
+        pointField coupleMasterPoints(coupleToMasterPoints.size());
+        pointField coupleSlavePoints(coupleToSlavePoints.size());
+        forAll(coupleToMasterPoints, couplePointi)
+        {
+            const label masterPointi = coupleToMasterPoints[couplePointi][0];
+            const label slavePointi = coupleToSlavePoints[couplePointi][0];
+            coupleMasterPoints[couplePointi] = m.localPoints()[masterPointi];
+            coupleSlavePoints[couplePointi] = m.localPoints()[slavePointi];
+        }
 
         writeOBJ
         (
-            "cutToMasterPoints.obj",
-            m.localPoints(),
-            pointField(c.localPoints(), masterToCutPoints_));
-    }
-    {
-        Pout<< "Writing cutToSlavePoints to cutToSlavePoints.obj" << endl;
-
-        writeOBJ
-        (
-            "cutToSlavePoints.obj",
-            s.localPoints(),
-            pointField(c.localPoints(), slaveToCutPoints_)
+            "masterToSlavePoints.obj",
+            coupleMasterPoints,
+            coupleSlavePoints
         );
     }
 
     // Face connectivity
     {
-        Pout<< "Writing cutToMasterFaces to cutToMasterFaces.obj" << endl;
-
-        pointField equivMasterFaces(c.size());
-
-        forAll(cutToMasterFaces(), cutFacei)
-        {
-            label masterFacei = cutToMasterFaces()[cutFacei];
-
-            if (masterFacei != -1)
-            {
-                equivMasterFaces[cutFacei] = m[masterFacei].centre(m.points());
-            }
-            else
-            {
-                WarningInFunction
-                    << "No master face for cut face " << cutFacei
-                    << " at position " << c[cutFacei].centre(c.points())
-                    << endl;
-
-                equivMasterFaces[cutFacei] = Zero;
-            }
-        }
+        Pout<< "Writing masterToSlaveFaces to masterToSlaveFaces.obj" << endl;
 
         writeOBJ
         (
-            "cutToMasterFaces.obj",
-            calcFaceCentres<List>(c, cutPoints(), 0, c.size()),
-            equivMasterFaces
-        );
-    }
-
-    {
-        Pout<< "Writing cutToSlaveFaces to cutToSlaveFaces.obj" << endl;
-
-        pointField equivSlaveFaces(c.size());
-
-        forAll(cutToSlaveFaces(), cutFacei)
-        {
-            label slaveFacei = cutToSlaveFaces()[cutFacei];
-
-            equivSlaveFaces[cutFacei] = s[slaveFacei].centre(s.points());
-        }
-
-        writeOBJ
-        (
-            "cutToSlaveFaces.obj",
-            calcFaceCentres<List>(c, cutPoints(), 0, c.size()),
-            equivSlaveFaces
+            "masterToSlaveFaces.obj",
+            calcFaceCentres<IndirectList>(m, m.points(), 0, m.size()),
+            calcFaceCentres<IndirectList>(s, s.points(), 0, s.size())
         );
     }
 
@@ -220,1804 +162,14 @@ void Foam::faceCoupleInfo::writePointsFaces() const
 }
 
 
-void Foam::faceCoupleInfo::writeEdges
-(
-    const labelList& cutToMasterEdges,
-    const labelList& cutToSlaveEdges
-) const
-{
-    const indirectPrimitivePatch& m = masterPatch();
-    const indirectPrimitivePatch& s = slavePatch();
-    const primitiveFacePatch& c = cutFaces();
-
-    // Edge connectivity
-    {
-        OFstream str("cutToMasterEdges.obj");
-        Pout<< "Writing cutToMasterEdges to " << str.name() << endl;
-
-        label vertI = 0;
-
-        forAll(cutToMasterEdges, cutEdgeI)
-        {
-            if (cutToMasterEdges[cutEdgeI] != -1)
-            {
-                const edge& masterEdge =
-                    m.edges()[cutToMasterEdges[cutEdgeI]];
-                const edge& cutEdge = c.edges()[cutEdgeI];
-
-                meshTools::writeOBJ(str, m.localPoints()[masterEdge[0]]);
-                vertI++;
-                meshTools::writeOBJ(str, m.localPoints()[masterEdge[1]]);
-                vertI++;
-                meshTools::writeOBJ(str, c.localPoints()[cutEdge[0]]);
-                vertI++;
-                meshTools::writeOBJ(str, c.localPoints()[cutEdge[1]]);
-                vertI++;
-                str << "l " << vertI-3 << ' ' << vertI-2 << nl;
-                str << "l " << vertI-3 << ' ' << vertI-1 << nl;
-                str << "l " << vertI-3 << ' ' << vertI << nl;
-                str << "l " << vertI-2 << ' ' << vertI-1 << nl;
-                str << "l " << vertI-2 << ' ' << vertI << nl;
-                str << "l " << vertI-1 << ' ' << vertI << nl;
-            }
-        }
-    }
-    {
-        OFstream str("cutToSlaveEdges.obj");
-        Pout<< "Writing cutToSlaveEdges to " << str.name() << endl;
-
-        label vertI = 0;
-
-        labelList slaveToCut(invert(s.nEdges(), cutToSlaveEdges));
-
-        forAll(slaveToCut, edgeI)
-        {
-            if (slaveToCut[edgeI] != -1)
-            {
-                const edge& slaveEdge = s.edges()[edgeI];
-                const edge& cutEdge = c.edges()[slaveToCut[edgeI]];
-
-                meshTools::writeOBJ(str, s.localPoints()[slaveEdge[0]]);
-                vertI++;
-                meshTools::writeOBJ(str, s.localPoints()[slaveEdge[1]]);
-                vertI++;
-                meshTools::writeOBJ(str, c.localPoints()[cutEdge[0]]);
-                vertI++;
-                meshTools::writeOBJ(str, c.localPoints()[cutEdge[1]]);
-                vertI++;
-                str << "l " << vertI-3 << ' ' << vertI-2 << nl;
-                str << "l " << vertI-3 << ' ' << vertI-1 << nl;
-                str << "l " << vertI-3 << ' ' << vertI << nl;
-                str << "l " << vertI-2 << ' ' << vertI-1 << nl;
-                str << "l " << vertI-2 << ' ' << vertI << nl;
-                str << "l " << vertI-1 << ' ' << vertI << nl;
-            }
-        }
-    }
-
-    Pout<< endl;
-}
-
-
-Foam::labelList Foam::faceCoupleInfo::findMappedEdges
-(
-    const edgeList& edges,
-    const labelList& pointMap,
-    const indirectPrimitivePatch& patch
-)
-{
-    labelList toPatchEdges(edges.size());
-
-    forAll(toPatchEdges, edgeI)
-    {
-        const edge& e = edges[edgeI];
-
-        label v0 = pointMap[e[0]];
-        label v1 = pointMap[e[1]];
-
-        toPatchEdges[edgeI] =
-            meshTools::findEdge
-            (
-                patch.edges(),
-                patch.pointEdges()[v0],
-                v0,
-                v1
-            );
-    }
-    return toPatchEdges;
-}
-
-
-bool Foam::faceCoupleInfo::regionEdge
-(
-    const polyMesh& slaveMesh,
-    const label slaveEdgeI
-) const
-{
-    const labelList& eFaces = slavePatch().edgeFaces()[slaveEdgeI];
-
-    if (eFaces.size() == 1)
-    {
-        return true;
-    }
-    else
-    {
-        // Count how many different patches connected to this edge.
-
-        label patch0 = -1;
-
-        forAll(eFaces, i)
-        {
-            label facei = eFaces[i];
-
-            label meshFacei = slavePatch().addressing()[facei];
-
-            label patchi = slaveMesh.boundaryMesh().whichPatch(meshFacei);
-
-            if (patch0 == -1)
-            {
-                patch0 = patchi;
-            }
-            else if (patchi != patch0)
-            {
-                // Found two different patches connected to this edge.
-                return true;
-            }
-        }
-        return false;
-    }
-}
-
-
-Foam::label Foam::faceCoupleInfo::mostAlignedCutEdge
-(
-    const bool report,
-    const polyMesh& slaveMesh,
-    const bool patchDivision,
-    const labelList& cutToMasterEdges,
-    const labelList& cutToSlaveEdges,
-    const label pointi,
-    const label edgeStart,
-    const label edgeEnd
-) const
-{
-    // Find edge using pointi that is most aligned with vector between master
-    // points. Patchdivision tells us whether or not to use patch information to
-    // match edges.
-
-    const pointField& localPoints = cutFaces().localPoints();
-
-    const labelList& pEdges = cutFaces().pointEdges()[pointi];
-
-    if (report)
-    {
-        Pout<< "mostAlignedEdge : finding nearest edge among "
-            << UIndirectList<edge>(cutFaces().edges(), pEdges)()
-            << " connected to point " << pointi
-            << " coord:" << localPoints[pointi]
-            << " running between " << edgeStart << " coord:"
-            << localPoints[edgeStart]
-            << " and " << edgeEnd << " coord:"
-            << localPoints[edgeEnd]
-            << endl;
-    }
-
-    // Find the edge that gets us nearest end.
-
-    label maxEdgeI = -1;
-    scalar maxCos = -great;
-
-    forAll(pEdges, i)
-    {
-        label edgeI = pEdges[i];
-
-        if
-        (
-           !(
-                patchDivision
-             && cutToMasterEdges[edgeI] == -1
-            )
-         || (
-                patchDivision
-             && regionEdge(slaveMesh, cutToSlaveEdges[edgeI])
-            )
-        )
-        {
-            const edge& e = cutFaces().edges()[edgeI];
-
-            label otherPointi = e.otherVertex(pointi);
-
-            if (otherPointi == edgeEnd)
-            {
-                // Shortcut: found edge end point.
-                if (report)
-                {
-                    Pout<< "    mostAlignedEdge : found end point " << edgeEnd
-                        << endl;
-                }
-                return edgeI;
-            }
-
-            // Get angle between edge and edge to masterEnd
-
-            vector eVec(localPoints[otherPointi] - localPoints[pointi]);
-
-            scalar magEVec = mag(eVec);
-
-            if (magEVec < vSmall)
-            {
-                WarningInFunction
-                    << "Crossing zero sized edge " << edgeI
-                    << " coords:" << localPoints[otherPointi]
-                    << localPoints[pointi]
-                    << " when walking from " << localPoints[edgeStart]
-                    << " to " << localPoints[edgeEnd]
-                    << endl;
-                return edgeI;
-            }
-
-            eVec /= magEVec;
-
-            vector eToEndPoint(localPoints[edgeEnd] - localPoints[otherPointi]);
-            eToEndPoint /= mag(eToEndPoint);
-
-            scalar cosAngle = eVec & eToEndPoint;
-
-            if (report)
-            {
-                Pout<< "    edge:" << e << " points:" << localPoints[pointi]
-                    << localPoints[otherPointi]
-                    << "  vec:" << eVec
-                    << "  vecToEnd:" << eToEndPoint
-                    << " cosAngle:" << cosAngle
-                    << endl;
-            }
-
-            if (cosAngle > maxCos)
-            {
-                maxCos = cosAngle;
-                maxEdgeI = edgeI;
-            }
-        }
-    }
-
-    if (maxCos > 1 - angleTol_)
-    {
-        return maxEdgeI;
-    }
-    else
-    {
-        return -1;
-    }
-}
-
-
-void Foam::faceCoupleInfo::setCutEdgeToPoints(const labelList& cutToMasterEdges)
-{
-    labelListList masterToCutEdges
-    (
-        invertOneToMany
-        (
-            masterPatch().nEdges(),
-            cutToMasterEdges
-        )
-    );
-
-    const edgeList& cutEdges = cutFaces().edges();
-
-    // Size extra big so searching is faster
-    cutEdgeToPoints_.resize
-    (
-        masterPatch().nEdges()
-      + slavePatch().nEdges()
-      + cutEdges.size()
-    );
-
-    forAll(masterToCutEdges, masterEdgeI)
-    {
-        const edge& masterE = masterPatch().edges()[masterEdgeI];
-
-        // Pout<< "Master:" << masterPatch().localPoints()[masterE[0]] << ' '
-        //    << masterPatch().localPoints()[masterE[1]] << endl;
-
-        const labelList& stringedEdges = masterToCutEdges[masterEdgeI];
-
-        if (stringedEdges.empty())
-        {
-            FatalErrorInFunction
-                << "Did not match all of master edges to cutFace edges"
-                << nl
-                << "First unmatched edge:" << masterEdgeI << " endPoints:"
-                << masterPatch().localPoints()[masterE[0]]
-                << masterPatch().localPoints()[masterE[1]]
-                << endl
-                << "This usually means that the slave patch is not a"
-                << " subdivision of the master patch"
-                << abort(FatalError);
-        }
-        else if (stringedEdges.size() > 1)
-        {
-            // String up the edges between e[0] and e[1]. Store the points
-            // in between e[0] and e[1] (all in cutFaces() labels)
-
-            DynamicList<label> splitPoints(stringedEdges.size()-1);
-
-            // Unsplit edge endpoints
-            const edge unsplitEdge
-            (
-                masterToCutPoints_[masterE[0]],
-                masterToCutPoints_[masterE[1]]
-            );
-
-            label startVertI = unsplitEdge[0];
-            label startEdgeI = -1;
-
-            while (startVertI != unsplitEdge[1])
-            {
-                // Loop over all string of edges. Update
-                // - startVertI : previous vertex
-                // - startEdgeI : previous edge
-                // and insert any points into splitPoints
-
-                // For checking
-                label oldStart = startVertI;
-
-                forAll(stringedEdges, i)
-                {
-                    label edgeI = stringedEdges[i];
-
-                    if (edgeI != startEdgeI)
-                    {
-                        const edge& e = cutEdges[edgeI];
-
-                        // Pout<< "    cut:" << e << " at:"
-                        //    << cutFaces().localPoints()[e[0]]
-                        //    << ' ' << cutFaces().localPoints()[e[1]] << endl;
-
-                        if (e[0] == startVertI)
-                        {
-                            startEdgeI = edgeI;
-                            startVertI = e[1];
-                            if (e[1] != unsplitEdge[1])
-                            {
-                                splitPoints.append(e[1]);
-                            }
-                            break;
-                        }
-                        else if (e[1] == startVertI)
-                        {
-                            startEdgeI = edgeI;
-                            startVertI = e[0];
-                            if (e[0] != unsplitEdge[1])
-                            {
-                                splitPoints.append(e[0]);
-                            }
-                            break;
-                        }
-                    }
-                }
-
-                // Check
-                if (oldStart == startVertI)
-                {
-                    FatalErrorInFunction
-                        << " unsplitEdge:" << unsplitEdge
-                        << " does not correspond to split edges "
-                        << UIndirectList<edge>(cutEdges, stringedEdges)()
-                        << abort(FatalError);
-                }
-            }
-
-            // Pout<< "For master edge:"
-            //    << unsplitEdge
-            //    << " Found stringed points "
-            //    <<  UIndirectList<point>
-            //        (
-            //            cutFaces().localPoints(),
-            //            splitPoints.shrink()
-            //        )()
-            //    << endl;
-
-            cutEdgeToPoints_.insert(unsplitEdge, splitPoints.shrink());
-        }
-    }
-}
-
-
-Foam::label Foam::faceCoupleInfo::matchFaces
-(
-    const scalar absTol,
-    const pointField& points0,
-    const face& f0,
-    const pointField& points1,
-    const face& f1,
-    const bool sameOrientation
-)
-{
-    if (f0.size() != f1.size())
-    {
-        FatalErrorInFunction
-            << "Different sizes for supposedly matching faces." << nl
-            << "f0:" << f0 << " coords:" << UIndirectList<point>(points0, f0)()
-            << nl
-            << "f1:" << f1 << " coords:" << UIndirectList<point>(points1, f1)()
-            << abort(FatalError);
-    }
-
-    const scalar absTolSqr = sign(absTol)*sqr(absTol);
-
-
-    label matchFp = -1;
-
-    forAll(f0, startFp)
-    {
-        // See -if starting from startFp on f0- the two faces match.
-        bool fullMatch = true;
-
-        label fp0 = startFp;
-        label fp1 = 0;
-
-        forAll(f1, i)
-        {
-            scalar distSqr = Foam::magSqr(points0[f0[fp0]] - points1[f1[fp1]]);
-
-            if (distSqr > absTolSqr)
-            {
-                fullMatch = false;
-                break;
-            }
-
-            fp0 = f0.fcIndex(fp0);  // walk forward
-
-            if (sameOrientation)
-            {
-                fp1 = f1.fcIndex(fp1);
-            }
-            else
-            {
-                fp1 = f1.rcIndex(fp1);
-            }
-        }
-
-        if (fullMatch)
-        {
-            matchFp = startFp;
-            break;
-        }
-    }
-
-    if (matchFp == -1)
-    {
-        FatalErrorInFunction
-            << "No unique match between two faces" << nl
-            << "Face " << f0 << " coords "
-            << UIndirectList<point>(points0, f0)() << nl
-            << "Face " << f1 << " coords "
-            << UIndirectList<point>(points1, f1)()
-            << "when using tolerance " << absTol
-            << " and forwardMatching:" << sameOrientation
-            << abort(FatalError);
-    }
-
-    return matchFp;
-}
-
-
-bool Foam::faceCoupleInfo::matchPointsThroughFaces
-(
-    const scalar absTol,
-    const pointField& cutPoints,
-    const faceList& cutFaces,
-    const pointField& patchPoints,
-    const faceList& patchFaces,
-    const bool sameOrientation,
-
-    labelList& patchToCutPoints,    // patch to (uncompacted) cut points
-    labelList& cutToCompact,        // compaction list for cut points
-    labelList& compactToCut         // inverse ,,
-)
-{
-    // Find correspondence from patch points to cut points. This might detect
-    // shared points so the output is a patch-to-cut point list and a compaction
-    // list for the cut points (which will always be equal or more connected
-    // than the patch).  Returns true if there are any duplicates.
-
-    // From slave to cut point
-    patchToCutPoints.setSize(patchPoints.size());
-    patchToCutPoints = -1;
-
-    // Compaction list for cut points: either -1 or index into master which
-    // gives the point to compact to.
-    labelList cutPointRegion(cutPoints.size(), -1);
-    DynamicList<label> cutPointRegionMaster;
-
-    forAll(patchFaces, patchFacei)
-    {
-        const face& patchF = patchFaces[patchFacei];
-
-        // const face& cutF = cutFaces[patchToCutFaces[patchFacei]];
-        const face& cutF = cutFaces[patchFacei];
-
-        // Do geometric matching to get position of cutF[0] in patchF
-        label patchFp = matchFaces
-        (
-            absTol,
-            patchPoints,
-            patchF,
-            cutPoints,
-            cutF,
-            sameOrientation        // orientation
-        );
-
-        forAll(cutF, cutFp)
-        {
-            label cutPointi = cutF[cutFp];
-            label patchPointi = patchF[patchFp];
-
-            // const point& cutPt = cutPoints[cutPointi];
-            // const point& patchPt = patchPoints[patchPointi];
-            // if (mag(cutPt - patchPt) > small)
-            //{
-            //    FatalErrorInFunction
-            //    << "cutP:" << cutPt
-            //    << " patchP:" << patchPt
-            //    << abort(FatalError);
-            //}
-
-            if (patchToCutPoints[patchPointi] == -1)
-            {
-                patchToCutPoints[patchPointi] = cutPointi;
-            }
-            else if (patchToCutPoints[patchPointi] != cutPointi)
-            {
-                // Multiple cut points connecting to same patch.
-                // Check if already have region & region master for this set
-                label otherCutPointi = patchToCutPoints[patchPointi];
-
-                // Pout<< "PatchPoint:" << patchPt
-                //    << " matches to:" << cutPointi
-                //    << " coord:" << cutPoints[cutPointi]
-                //    << " and to:" << otherCutPointi
-                //    << " coord:" << cutPoints[otherCutPointi]
-                //    << endl;
-
-                if (cutPointRegion[otherCutPointi] != -1)
-                {
-                    // Have region for this set. Copy.
-                    label region = cutPointRegion[otherCutPointi];
-                    cutPointRegion[cutPointi] = region;
-
-                    // Update region master with min point label
-                    cutPointRegionMaster[region] = min
-                    (
-                        cutPointRegionMaster[region],
-                        cutPointi
-                    );
-                }
-                else
-                {
-                    // Create new region.
-                    label region = cutPointRegionMaster.size();
-                    cutPointRegionMaster.append
-                    (
-                        min(cutPointi, otherCutPointi)
-                    );
-                    cutPointRegion[cutPointi] = region;
-                    cutPointRegion[otherCutPointi] = region;
-                }
-            }
-
-            if (sameOrientation)
-            {
-                patchFp = patchF.fcIndex(patchFp);
-            }
-            else
-            {
-                patchFp = patchF.rcIndex(patchFp);
-            }
-        }
-    }
-
-    // Rework region&master into compaction array
-    compactToCut.setSize(cutPointRegion.size());
-    cutToCompact.setSize(cutPointRegion.size());
-    cutToCompact = -1;
-    label compactPointi = 0;
-
-    forAll(cutPointRegion, i)
-    {
-        if (cutPointRegion[i] == -1)
-        {
-            // Unduplicated point. Allocate new compacted point.
-            cutToCompact[i] = compactPointi;
-            compactToCut[compactPointi] = i;
-            compactPointi++;
-        }
-        else
-        {
-            // Duplicate point. Get master.
-
-            label masterPointi = cutPointRegionMaster[cutPointRegion[i]];
-
-            if (cutToCompact[masterPointi] == -1)
-            {
-                cutToCompact[masterPointi] = compactPointi;
-                compactToCut[compactPointi] = masterPointi;
-                compactPointi++;
-            }
-            cutToCompact[i] = cutToCompact[masterPointi];
-        }
-    }
-    compactToCut.setSize(compactPointi);
-
-    return compactToCut.size() != cutToCompact.size();
-}
-
-
-Foam::scalar Foam::faceCoupleInfo::maxDistance
-(
-    const face& cutF,
-    const pointField& cutPoints,
-    const face& masterF,
-    const pointField& masterPoints
-)
-{
-    scalar maxDist = -great;
-
-    forAll(cutF, fp)
-    {
-        const point& cutPt = cutPoints[cutF[fp]];
-
-        pointHit pHit = masterF.nearestPoint(cutPt, masterPoints);
-
-        maxDist = max(maxDist, pHit.distance());
-    }
-    return maxDist;
-}
-
-
-void Foam::faceCoupleInfo::findPerfectMatchingFaces
-(
-    const primitiveMesh& mesh0,
-    const primitiveMesh& mesh1,
-    const scalar absTol,
-
-    labelList& mesh0Faces,
-    labelList& mesh1Faces
-)
-{
-    // Face centres of external faces (without invoking
-    // mesh.faceCentres since mesh might have been clearedOut)
-
-    pointField fc0
-    (
-        calcFaceCentres<List>
-        (
-            mesh0.faces(),
-            mesh0.points(),
-            mesh0.nInternalFaces(),
-            mesh0.nFaces() - mesh0.nInternalFaces()
-        )
-    );
-
-    pointField fc1
-    (
-        calcFaceCentres<List>
-        (
-            mesh1.faces(),
-            mesh1.points(),
-            mesh1.nInternalFaces(),
-            mesh1.nFaces() - mesh1.nInternalFaces()
-        )
-    );
-
-
-    if (debug)
-    {
-        Pout<< "Face matching tolerance : " << absTol << endl;
-    }
-
-
-    // Match geometrically
-    labelList from1To0;
-    bool matchedAllFaces = matchPoints
-    (
-        fc1,
-        fc0,
-        scalarField(fc1.size(), absTol),
-        false,
-        from1To0
-    );
-
-    if (matchedAllFaces)
-    {
-        WarningInFunction
-            << "Matched ALL " << fc1.size()
-            << " boundary faces of mesh0 to boundary faces of mesh1." << endl
-            << "This is only valid if the mesh to add is fully"
-            << " enclosed by the mesh it is added to." << endl;
-    }
-
-
-    // Collect matches.
-    label nMatched = 0;
-
-    mesh0Faces.setSize(fc0.size());
-    mesh1Faces.setSize(fc1.size());
-
-    forAll(from1To0, i)
-    {
-        if (from1To0[i] != -1)
-        {
-            mesh1Faces[nMatched] = i + mesh1.nInternalFaces();
-            mesh0Faces[nMatched] = from1To0[i] + mesh0.nInternalFaces();
-
-            nMatched++;
-        }
-    }
-
-    mesh0Faces.setSize(nMatched);
-    mesh1Faces.setSize(nMatched);
-}
-
-
-void Foam::faceCoupleInfo::findSlavesCoveringMaster
-(
-    const primitiveMesh& mesh0,
-    const primitiveMesh& mesh1,
-    const scalar absTol,
-
-    labelList& mesh0Faces,
-    labelList& mesh1Faces
-)
-{
-    // Construct octree from all mesh0 boundary faces
-    labelList bndFaces(mesh0.nFaces()-mesh0.nInternalFaces());
-    forAll(bndFaces, i)
-    {
-        bndFaces[i] = mesh0.nInternalFaces() + i;
-    }
-
-    treeBoundBox overallBb(mesh0.points());
-
-    indexedOctree<treeDataFace> tree
-    (
-        treeDataFace    // all information needed to search faces
-        (
-            false,                      // do not cache bb
-            mesh0,
-            bndFaces                    // boundary faces only
-        ),
-        overallBb.extend(1e-4),         // overall search domain
-        8,                              // maxLevel
-        10,                             // leafsize
-        3.0                             // duplicity
-    );
-
-    if (debug)
-    {
-        Pout<< "findSlavesCoveringMaster :"
-            << " constructed octree for mesh0 boundary faces" << endl;
-    }
-
-    // Search nearest face for every mesh1 boundary face.
-
-    labelHashSet mesh0Set(mesh0.nFaces() - mesh0.nInternalFaces());
-    labelHashSet mesh1Set(mesh1.nFaces() - mesh1.nInternalFaces());
-
-    for
-    (
-        label mesh1Facei = mesh1.nInternalFaces();
-        mesh1Facei < mesh1.nFaces();
-        mesh1Facei++
-    )
-    {
-        const face& f1 = mesh1.faces()[mesh1Facei];
-
-        // Generate face centre (prevent cellCentres() reconstruction)
-        point fc(f1.centre(mesh1.points()));
-
-        pointIndexHit nearInfo = tree.findNearest(fc, sign(absTol)*sqr(absTol));
-
-        if (nearInfo.hit())
-        {
-            label mesh0Facei = tree.shapes().faceLabels()[nearInfo.index()];
-
-            // Check if points of f1 actually lie on top of mesh0 face
-            // This is the bit that might fail since if f0 is severely warped
-            // and f1's points are not present on f0 (since subdivision)
-            // then the distance of the points to f0 might be quite large
-            // and the test will fail. This all should in fact be some kind
-            // of walk from known corresponding points/edges/faces.
-            scalar dist =
-                maxDistance
-                (
-                    f1,
-                    mesh1.points(),
-                    mesh0.faces()[mesh0Facei],
-                    mesh0.points()
-                );
-
-            if (dist < absTol)
-            {
-                mesh0Set.insert(mesh0Facei);
-                mesh1Set.insert(mesh1Facei);
-            }
-        }
-    }
-
-    if (debug)
-    {
-        Pout<< "findSlavesCoveringMaster :"
-            << " matched " << mesh1Set.size() << " mesh1 faces to "
-            << mesh0Set.size() << " mesh0 faces" << endl;
-    }
-
-    mesh0Faces = mesh0Set.toc();
-    mesh1Faces = mesh1Set.toc();
-}
-
-
-Foam::label Foam::faceCoupleInfo::growCutFaces
-(
-    const labelList& cutToMasterEdges,
-    Map<labelList>& candidates
-)
-{
-    if (debug)
-    {
-        Pout<< "growCutFaces :"
-            << " growing cut faces to masterPatch" << endl;
-    }
-
-    label nTotChanged = 0;
-
-    while (true)
-    {
-        const labelListList& cutFaceEdges = cutFaces().faceEdges();
-        const labelListList& cutEdgeFaces = cutFaces().edgeFaces();
-
-        label nChanged = 0;
-
-        forAll(cutToMasterFaces_, cutFacei)
-        {
-            const label masterFacei = cutToMasterFaces_[cutFacei];
-
-            if (masterFacei != -1)
-            {
-                // We now have a cutFace for which we have already found a
-                // master face. Grow this masterface across any internal edge
-                // (internal: no corresponding master edge)
-
-                const labelList& fEdges = cutFaceEdges[cutFacei];
-
-                forAll(fEdges, i)
-                {
-                    const label cutEdgeI = fEdges[i];
-
-                    if (cutToMasterEdges[cutEdgeI] == -1)
-                    {
-                        // So this edge:
-                        // - internal to the cutPatch (since all region edges
-                        //   marked before)
-                        // - on cutFacei which corresponds to masterFace.
-                        // Mark all connected faces with this masterFace.
-
-                        const labelList& eFaces = cutEdgeFaces[cutEdgeI];
-
-                        forAll(eFaces, j)
-                        {
-                            const label facei = eFaces[j];
-
-                            if (cutToMasterFaces_[facei] == -1)
-                            {
-                                cutToMasterFaces_[facei] = masterFacei;
-                                candidates.erase(facei);
-                                nChanged++;
-                            }
-                            else if (cutToMasterFaces_[facei] != masterFacei)
-                            {
-                                const pointField& cutPoints =
-                                    cutFaces().points();
-                                const pointField& masterPoints =
-                                    masterPatch().points();
-
-                                const edge& e = cutFaces().edges()[cutEdgeI];
-
-                                label myMaster = cutToMasterFaces_[facei];
-                                const face& myF = masterPatch()[myMaster];
-
-                                const face& nbrF = masterPatch()[masterFacei];
-
-                                FatalErrorInFunction
-                                    << "Cut face "
-                                    << cutFaces()[facei].points(cutPoints)
-                                    << " has master "
-                                    << myMaster
-                                    << " but also connects to nbr face "
-                                    << cutFaces()[cutFacei].points(cutPoints)
-                                    << " with master " << masterFacei
-                                    << nl
-                                    << "myMasterFace:"
-                                    << myF.points(masterPoints)
-                                    << "  nbrMasterFace:"
-                                    << nbrF.points(masterPoints) << nl
-                                    << "Across cut edge " << cutEdgeI
-                                    << " coord:"
-                                    << cutFaces().localPoints()[e[0]]
-                                    << cutFaces().localPoints()[e[1]]
-                                    << abort(FatalError);
-                            }
-                        }
-                    }
-                }
-            }
-        }
-
-        if (debug)
-        {
-            Pout<< "growCutFaces : Grown an additional " << nChanged
-                << " cut to master face" << " correspondence" << endl;
-        }
-
-        nTotChanged += nChanged;
-
-        if (nChanged == 0)
-        {
-            break;
-        }
-    }
-
-    return nTotChanged;
-}
-
-
-void Foam::faceCoupleInfo::checkMatch(const labelList& cutToMasterEdges) const
-{
-    const pointField& cutLocalPoints = cutFaces().localPoints();
-
-    const pointField& masterLocalPoints = masterPatch().localPoints();
-    const faceList& masterLocalFaces = masterPatch().localFaces();
-
-    forAll(cutToMasterEdges, cutEdgeI)
-    {
-        const edge& e = cutFaces().edges()[cutEdgeI];
-
-        if (cutToMasterEdges[cutEdgeI] == -1)
-        {
-            // Internal edge. Check that master face is same on both sides.
-            const labelList& cutEFaces = cutFaces().edgeFaces()[cutEdgeI];
-
-            label masterFacei = -1;
-
-            forAll(cutEFaces, i)
-            {
-                label cutFacei = cutEFaces[i];
-
-                if (cutToMasterFaces_[cutFacei] != -1)
-                {
-                    if (masterFacei == -1)
-                    {
-                        masterFacei = cutToMasterFaces_[cutFacei];
-                    }
-                    else if (masterFacei != cutToMasterFaces_[cutFacei])
-                    {
-                        label myMaster = cutToMasterFaces_[cutFacei];
-                        const face& myF = masterLocalFaces[myMaster];
-
-                        const face& nbrF = masterLocalFaces[masterFacei];
-
-                        FatalErrorInFunction
-                            << "Internal CutEdge " << e
-                            << " coord:"
-                            << cutLocalPoints[e[0]]
-                            << cutLocalPoints[e[1]]
-                            << " connects to master " << myMaster
-                            << " and to master " << masterFacei << nl
-                            << "myMasterFace:"
-                            << myF.points(masterLocalPoints)
-                            << "  nbrMasterFace:"
-                            << nbrF.points(masterLocalPoints)
-                            << abort(FatalError);
-                    }
-                }
-            }
-        }
-    }
-}
-
-
-Foam::label Foam::faceCoupleInfo::matchEdgeFaces
-(
-    const labelList& cutToMasterEdges,
-    Map<labelList>& candidates
-)
-{
-    // Extends matching information by elimination across cutFaces using more
-    // than one region edge. Updates cutToMasterFaces_ and sets candidates which
-    // is for every cutface on a region edge the possible master faces.
-
-    // For every unassigned cutFacei the possible list of master faces.
-    candidates.clear();
-    candidates.resize(cutFaces().size());
-
-    label nChanged = 0;
-
-    forAll(cutToMasterEdges, cutEdgeI)
-    {
-        label masterEdgeI = cutToMasterEdges[cutEdgeI];
-
-        if (masterEdgeI != -1)
-        {
-            // So cutEdgeI is matched to masterEdgeI. For all cut faces using
-            // the cut edge store the master face as a candidate match.
-
-            const labelList& cutEFaces = cutFaces().edgeFaces()[cutEdgeI];
-            const labelList& masterEFaces =
-                masterPatch().edgeFaces()[masterEdgeI];
-
-            forAll(cutEFaces, i)
-            {
-                label cutFacei = cutEFaces[i];
-
-                if (cutToMasterFaces_[cutFacei] == -1)
-                {
-                    // So this face (cutFacei) is on an edge (cutEdgeI) that
-                    // is used by some master faces (masterEFaces). Check if
-                    // the cutFacei and some of these masterEFaces share more
-                    // than one edge (which uniquely defines face)
-
-                    // Combine master faces with current set of candidate
-                    // master faces.
-                    Map<labelList>::iterator fnd = candidates.find(cutFacei);
-
-                    if (fnd == candidates.end())
-                    {
-                        // No info yet for cutFacei. Add all master faces as
-                        // candidates
-                        candidates.insert(cutFacei, masterEFaces);
-                    }
-                    else
-                    {
-                        // From some other cutEdgeI there are already some
-                        // candidate master faces. Check the overlap with
-                        // the current set of master faces.
-                        const labelList& masterFaces = fnd();
-
-                        DynamicList<label> newCandidates(masterFaces.size());
-
-                        forAll(masterEFaces, j)
-                        {
-                            if (findIndex(masterFaces, masterEFaces[j]) != -1)
-                            {
-                                newCandidates.append(masterEFaces[j]);
-                            }
-                        }
-
-                        if (newCandidates.size() == 1)
-                        {
-                            // We found a perfect match. Delete entry from
-                            // candidates map.
-                            cutToMasterFaces_[cutFacei] = newCandidates[0];
-                            candidates.erase(cutFacei);
-                            nChanged++;
-                        }
-                        else
-                        {
-                            // Should not happen?
-                            // Pout<< "On edge:" << cutEdgeI
-                            //    << " have connected masterFaces:"
-                            //    << masterEFaces
-                            //    << " and from previous edge we have"
-                            //    << " connected masterFaces" << masterFaces
-                            //    << " . Overlap:" << newCandidates << endl;
-
-                            fnd() = newCandidates.shrink();
-                        }
-                    }
-                }
-
-            }
-        }
-    }
-
-    if (debug)
-    {
-        Pout<< "matchEdgeFaces : Found " << nChanged
-            << " faces where there was"
-            << " only one remaining choice for cut-master correspondence"
-            << endl;
-    }
-
-    return nChanged;
-}
-
-
-Foam::label Foam::faceCoupleInfo::geometricMatchEdgeFaces
-(
-    Map<labelList>& candidates
-)
-{
-    const pointField& cutPoints = cutFaces().points();
-
-    label nChanged = 0;
-
-    // Mark all master faces that have been matched so far.
-
-    labelListList masterToCutFaces
-    (
-        invertOneToMany
-        (
-            masterPatch().size(),
-            cutToMasterFaces_
-        )
-    );
-
-    forAllConstIter(Map<labelList>, candidates, iter)
-    {
-        label cutFacei = iter.key();
-
-        const face& cutF = cutFaces()[cutFacei];
-
-        if (cutToMasterFaces_[cutFacei] == -1)
-        {
-            const labelList& masterFaces = iter();
-
-            // Find the best matching master face.
-            scalar minDist = great;
-            label minMasterFacei = -1;
-
-            forAll(masterFaces, i)
-            {
-                label masterFacei = masterFaces[i];
-
-                if (masterToCutFaces[masterFaces[i]].empty())
-                {
-                    scalar dist = maxDistance
-                    (
-                        cutF,
-                        cutPoints,
-                        masterPatch()[masterFacei],
-                        masterPatch().points()
-                    );
-
-                    if (dist < minDist)
-                    {
-                        minDist = dist;
-                        minMasterFacei = masterFacei;
-                    }
-                }
-            }
-
-            if (minMasterFacei != -1)
-            {
-                cutToMasterFaces_[cutFacei] = minMasterFacei;
-                masterToCutFaces[minMasterFacei] = cutFacei;
-                nChanged++;
-            }
-        }
-    }
-
-    // (inefficiently) force candidates to be uptodate.
-    forAll(cutToMasterFaces_, cutFacei)
-    {
-        if (cutToMasterFaces_[cutFacei] != -1)
-        {
-            candidates.erase(cutFacei);
-        }
-    }
-
-
-    if (debug)
-    {
-        Pout<< "geometricMatchEdgeFaces : Found " << nChanged
-            << " faces where there was"
-            << " only one remaining choice for cut-master correspondence"
-            << endl;
-    }
-
-    return nChanged;
-}
-
-
-void Foam::faceCoupleInfo::perfectOrderedPointMatch()
-{
-    // Calculate the set of cut faces in between master and slave patch assuming
-    // perfect match (and optional face ordering on slave)
-
-    if (debug)
-    {
-        Pout<< "perfectOrderedPointMatch :"
-            << " Matching master and slave to cut."
-            << " Master and slave faces are identical and ordered" << nl;
-    }
-
-    // Faces correspond exactly
-    cutToMasterFaces_ = identity(masterPatch().size());
-    cutToSlaveFaces_ = identity(masterPatch().size());
-
-    // Cut patch is the same as the master patch
-    cutPoints_ = masterPatch().localPoints();
-    cutFacesPtr_.reset
-    (
-        new primitiveFacePatch
-        (
-            masterPatch().localFaces(),
-            cutPoints_
-        )
-    );
-
-    // Cut points are the same as the master points
-    masterToCutPoints_ = identity(cutPoints_.size());
-
-    // Slave points number around the faces in the reverse direction
-    slaveToCutPoints_.setSize(slavePatch().nPoints());
-    forAll(cutFaces(), cutFacei)
-    {
-        const label slaveFacei = cutToSlaveFaces_[cutFacei];
-
-        const face& cutF = cutFaces().localFaces()[cutFacei];
-        const face& slaveF = slavePatch().localFaces()[slaveFacei];
-
-        label slaveFp = 0;
-
-        forAll(cutF, cutFp)
-        {
-            const label cutPointi = cutF[cutFp];
-            const label slavePointi = slaveF[slaveFp];
-
-            slaveToCutPoints_[slavePointi] = cutPointi;
-
-            slaveFp = slaveF.rcIndex(slaveFp);
-        }
-    }
-}
-
-
-void Foam::faceCoupleInfo::perfectPointMatch(const scalar absTol)
-{
-    // Calculate the set of cut faces in between master and slave patch assuming
-    // perfect match (and optional face ordering on slave)
-
-    if (debug)
-    {
-        Pout<< "perfectPointMatch :"
-            << " Matching master and slave to cut."
-            << " Master and slave faces are identical and not ordered" << nl;
-    }
-
-    cutToMasterFaces_ = identity(masterPatch().size());
-    cutPoints_ = masterPatch().localPoints();
-    cutFacesPtr_.reset
-    (
-        new primitiveFacePatch
-        (
-            masterPatch().localFaces(),
-            cutPoints_
-        )
-    );
-    masterToCutPoints_ = identity(cutPoints_.size());
-
-    // Faces do not have to be ordered, but all have to match
-
-    bool matchedAllFaces = matchPoints
-    (
-        calcFaceCentres<List>
-        (
-            cutFaces(),
-            cutFaces().points(),
-            0,
-            cutFaces().size()
-        ),
-        calcFaceCentres<IndirectList>
-        (
-            slavePatch(),
-            slavePatch().points(),
-            0,
-            slavePatch().size()
-        ),
-        scalarField(slavePatch().size(), absTol),
-        false,
-        cutToSlaveFaces_
-    );
-
-    // If some of the face centres did not match, then try to match the
-    // point averages instead. There is no division by the face area in
-    // calculating the point average, so this is more stable when faces
-    // collapse onto a line or point.
-    if (!matchedAllFaces)
-    {
-        labelList cutToSlaveFacesTemp(cutToSlaveFaces_.size(), -1);
-
-        matchPoints
-        (
-            calcFacePointAverages<List>
-            (
-                cutFaces(),
-                cutFaces().points(),
-                0,
-                cutFaces().size()
-            ),
-            calcFacePointAverages<IndirectList>
-            (
-                slavePatch(),
-                slavePatch().points(),
-                0,
-                slavePatch().size()
-            ),
-            scalarField(slavePatch().size(), absTol),
-            true,
-            cutToSlaveFacesTemp
-        );
-
-        cutToSlaveFaces_ = max(cutToSlaveFaces_, cutToSlaveFacesTemp);
-
-        matchedAllFaces = min(cutToSlaveFaces_) != -1;
-    }
-
-    if (!matchedAllFaces)
-    {
-        FatalErrorInFunction
-            << "Did not match all of the master faces to the slave faces"
-            << endl
-            << "This usually means that the slave patch and master patch"
-            << " do not align to within " << absTol << " metre."
-            << abort(FatalError);
-    }
-
-    // Find correspondence from slave points to cut points. This might
-    // detect shared points so the output is a slave-to-cut point list
-    // and a compaction list.
-
-    labelList cutToCompact, compactToCut;
-    matchPointsThroughFaces
-    (
-        absTol,
-        cutFaces().localPoints(),
-        reorder(cutToSlaveFaces_, cutFaces().localFaces()),
-        slavePatch().localPoints(),
-        slavePatch().localFaces(),
-        false,                      // slave and cut have opposite
-                                    // orientation
-        slaveToCutPoints_,          // slave to (uncompacted) cut points
-        cutToCompact,               // compaction map: from cut to compacted
-        compactToCut                // compaction map: from compacted to cut
-    );
-
-
-    // Use compaction lists to renumber cutPoints.
-    cutPoints_ = UIndirectList<point>(cutPoints_, compactToCut)();
-    {
-        const faceList& cutLocalFaces = cutFaces().localFaces();
-
-        faceList compactFaces(cutLocalFaces.size());
-        forAll(cutLocalFaces, i)
-        {
-            compactFaces[i] = renumber(cutToCompact, cutLocalFaces[i]);
-        }
-        cutFacesPtr_.reset
-        (
-            new primitiveFacePatch
-            (
-                compactFaces,
-                cutPoints_
-            )
-        );
-    }
-    inplaceRenumber(cutToCompact, slaveToCutPoints_);
-    inplaceRenumber(cutToCompact, masterToCutPoints_);
-}
-
-
-void Foam::faceCoupleInfo::subDivisionMatch
-(
-    const polyMesh& slaveMesh,
-    const bool patchDivision,
-    const scalar absTol
-)
-{
-    if (debug)
-    {
-        Pout<< "subDivisionMatch :"
-            << " Matching master and slave to cut."
-            << " Slave can be subdivision of master but all master edges"
-            << " have to be covered by slave edges." << endl;
-    }
-
-    // Assume slave patch is subdivision of the master patch so cutFaces is a
-    // copy of the slave (but reversed since orientation has to be like master).
-    cutPoints_ = slavePatch().localPoints();
-    {
-        faceList cutFaces(slavePatch().size());
-
-        forAll(cutFaces, i)
-        {
-            cutFaces[i] = slavePatch().localFaces()[i].reverseFace();
-        }
-        cutFacesPtr_.reset(new primitiveFacePatch(cutFaces, cutPoints_));
-    }
-
-    // Cut is copy of slave so addressing to slave is trivial.
-    cutToSlaveFaces_ = identity(cutFaces().size());
-    slaveToCutPoints_ = identity(slavePatch().nPoints());
-
-    // Cut edges to slave patch
-    labelList cutToSlaveEdges
-    (
-        findMappedEdges
-        (
-            cutFaces().edges(),
-            slaveToCutPoints_,
-            slavePatch()
-        )
-    );
-
-
-    if (debug)
-    {
-        OFstream str("regionEdges.obj");
-
-        Pout<< "subDivisionMatch :"
-            << " addressing for slave patch fully done."
-            << " Dumping region edges to " << str.name() << endl;
-
-        label vertI = 0;
-
-        forAll(slavePatch().edges(), slaveEdgeI)
-        {
-            if (regionEdge(slaveMesh, slaveEdgeI))
-            {
-                const edge& e = slavePatch().edges()[slaveEdgeI];
-
-                meshTools::writeOBJ(str, slavePatch().localPoints()[e[0]]);
-                vertI++;
-                meshTools::writeOBJ(str, slavePatch().localPoints()[e[1]]);
-                vertI++;
-                str<< "l " << vertI-1 << ' ' << vertI << nl;
-            }
-        }
-    }
-
-
-    // Addressing from cut to master.
-
-    // Cut points to master patch
-    // - determine master points to cut points. Has to be full!
-    // - invert to get cut to master
-    if (debug)
-    {
-        Pout<< "subDivisionMatch :"
-            << " matching master points to cut points" << endl;
-    }
-
-    bool matchedAllPoints = matchPoints
-    (
-        masterPatch().localPoints(),
-        cutPoints_,
-        scalarField(masterPatch().nPoints(), absTol),
-        true,
-        masterToCutPoints_
-    );
-
-    if (!matchedAllPoints)
-    {
-        FatalErrorInFunction
-            << "Did not match all of the master points to the slave points"
-            << endl
-            << "This usually means that the slave patch is not a"
-            << " subdivision of the master patch"
-            << abort(FatalError);
-    }
-
-
-    // Do masterEdges to cutEdges :
-    // - mark all edges between two masterEdge endpoints. (geometric test since
-    //   this is the only distinction)
-    // - this gives the borders in between which all cutfaces come from
-    //   a single master face.
-    if (debug)
-    {
-        Pout<< "subDivisionMatch :"
-            << " matching cut edges to master edges" << endl;
-    }
-
-    const edgeList& masterEdges = masterPatch().edges();
-    const pointField& masterPoints = masterPatch().localPoints();
-
-    const edgeList& cutEdges = cutFaces().edges();
-
-    labelList cutToMasterEdges(cutFaces().nEdges(), -1);
-
-    forAll(masterEdges, masterEdgeI)
-    {
-        const edge& masterEdge = masterEdges[masterEdgeI];
-
-        label cutPoint0 = masterToCutPoints_[masterEdge[0]];
-        label cutPoint1 = masterToCutPoints_[masterEdge[1]];
-
-        // Find edges between cutPoint0 and cutPoint1.
-
-        label cutPointi = cutPoint0;
-        do
-        {
-            // Find edge (starting at pointi on cut), aligned with master
-            // edge.
-            label cutEdgeI =
-                mostAlignedCutEdge
-                (
-                    false,
-                    slaveMesh,
-                    patchDivision,
-                    cutToMasterEdges,
-                    cutToSlaveEdges,
-                    cutPointi,
-                    cutPoint0,
-                    cutPoint1
-                );
-
-            if (cutEdgeI == -1)
-            {
-                // Problem. Report while matching to produce nice error message
-                mostAlignedCutEdge
-                (
-                    true,
-                    slaveMesh,
-                    patchDivision,
-                    cutToMasterEdges,
-                    cutToSlaveEdges,
-                    cutPointi,
-                    cutPoint0,
-                    cutPoint1
-                );
-
-                Pout<< "Dumping unmatched pointEdges to errorEdges.obj"
-                    << endl;
-
-                writeOBJ
-                (
-                    "errorEdges.obj",
-                    edgeList
-                    (
-                        UIndirectList<edge>
-                        (
-                            cutFaces().edges(),
-                            cutFaces().pointEdges()[cutPointi]
-                        )
-                    ),
-                    cutFaces().localPoints(),
-                    false
-                );
-
-                FatalErrorInFunction
-                    << "Problem in finding cut edges corresponding to"
-                    << " master edge " << masterEdge
-                    << " points:" << masterPoints[masterEdge[0]]
-                    << ' ' << masterPoints[masterEdge[1]]
-                    << " corresponding cut edge: (" << cutPoint0
-                    << ") " << cutPoint1
-                    << abort(FatalError);
-            }
-
-            cutToMasterEdges[cutEdgeI] = masterEdgeI;
-
-            cutPointi = cutEdges[cutEdgeI].otherVertex(cutPointi);
-
-        } while (cutPointi != cutPoint1);
-    }
-
-    if (debug)
-    {
-        // Write all matched edges.
-        writeEdges(cutToMasterEdges, cutToSlaveEdges);
-    }
-
-    // Rework cutToMasterEdges into list of points in between two endpoints
-    // (cutEdgeToPoints_)
-    setCutEdgeToPoints(cutToMasterEdges);
-
-
-    // Now that we have marked the cut edges that lie on top of master edges
-    // we can find cut faces that have two (or more) of these edges.
-    // Note that there can be multiple faces having two or more matched edges
-    // since the cut faces can form a non-manifold surface(?)
-    // So the matching is done as an elimination process where for every
-    // cutFace on a matched edge we store the possible master faces and
-    // eliminate more and more until we only have one possible master face
-    // left.
-
-    if (debug)
-    {
-        Pout<< "subDivisionMatch :"
-            << " matching (topological) cut faces to masterPatch" << endl;
-    }
-
-    // For every unassigned cutFacei the possible list of master faces.
-    Map<labelList> candidates(cutFaces().size());
-
-    cutToMasterFaces_.setSize(cutFaces().size());
-    cutToMasterFaces_ = -1;
-
-    while (true)
-    {
-        // See if there are any edges left where there is only one remaining
-        // candidate.
-        label nChanged = matchEdgeFaces(cutToMasterEdges, candidates);
-
-        checkMatch(cutToMasterEdges);
-
-
-        // Now we should have found a face correspondence for every cutFace
-        // that uses two or more cutEdges. Floodfill from these across
-        // 'internal' cutedges (i.e. edges that do not have a master
-        // equivalent) to determine some more cutToMasterFaces
-        nChanged += growCutFaces(cutToMasterEdges, candidates);
-
-        checkMatch(cutToMasterEdges);
-
-        if (nChanged == 0)
-        {
-            break;
-        }
-    }
-
-    if (debug)
-    {
-        Pout<< "subDivisionMatch :"
-            << " matching (geometric) cut faces to masterPatch" << endl;
-    }
-
-    // Above should have matched all cases fully. Cannot prove this yet so
-    // do any remaining unmatched edges by a geometric test.
-    while (true)
-    {
-        // See if there are any edges left where there is only one remaining
-        // candidate.
-        label nChanged = geometricMatchEdgeFaces(candidates);
-
-        checkMatch(cutToMasterEdges);
-
-        nChanged += growCutFaces(cutToMasterEdges, candidates);
-
-        checkMatch(cutToMasterEdges);
-
-        if (nChanged == 0)
-        {
-            break;
-        }
-    }
-
-
-    // All cut faces matched?
-    forAll(cutToMasterFaces_, cutFacei)
-    {
-        if (cutToMasterFaces_[cutFacei] == -1)
-        {
-            const face& cutF = cutFaces()[cutFacei];
-
-            FatalErrorInFunction
-                << "Did not match all of cutFaces to a master face" << nl
-                << "First unmatched cut face:" << cutFacei << " with points:"
-                << UIndirectList<point>(cutFaces().points(), cutF)()
-                << nl
-                << "This usually means that the slave patch is not a"
-                << " subdivision of the master patch"
-                << abort(FatalError);
-        }
-    }
-
-    if (debug)
-    {
-        Pout<< "subDivisionMatch :"
-            << " finished matching master and slave to cut" << endl;
-    }
-
-    if (debug)
-    {
-        writePointsFaces();
-    }
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::faceCoupleInfo::faceCoupleInfo
-(
-    const polyMesh& masterMesh,
-    const polyMesh& slaveMesh,
-    const scalar absTol,
-    const bool perfectMatch
-)
-:
-    masterPatchPtr_(nullptr),
-    slavePatchPtr_(nullptr),
-    cutPoints_(0),
-    cutFacesPtr_(nullptr),
-    cutToMasterFaces_(0),
-    masterToCutPoints_(0),
-    cutToSlaveFaces_(0),
-    slaveToCutPoints_(0),
-    cutEdgeToPoints_(0)
-{
-    // Get faces on both meshes that are aligned.
-    // (not ordered i.e. masterToMesh[0] does
-    // not couple to slaveToMesh[0]. This ordering is done later on)
-    labelList masterToMesh;
-    labelList slaveToMesh;
-
-    if (perfectMatch)
-    {
-        // Identical faces so use tight face-centre comparison.
-        findPerfectMatchingFaces
-        (
-            masterMesh,
-            slaveMesh,
-            absTol,
-            masterToMesh,
-            slaveToMesh
-        );
-    }
-    else
-    {
-        // Slave subdivision of master so use 'nearest'. Bit dodgy, especially
-        // with using absTol (which is quite small)
-        findSlavesCoveringMaster
-        (
-            masterMesh,
-            slaveMesh,
-            absTol,
-            masterToMesh,
-            slaveToMesh
-        );
-    }
-
-    // Construct addressing engines for both sides
-    masterPatchPtr_.reset
-    (
-        new indirectPrimitivePatch
-        (
-            IndirectList<face>(masterMesh.faces(), masterToMesh),
-            masterMesh.points()
-        )
-    );
-
-    slavePatchPtr_.reset
-    (
-        new indirectPrimitivePatch
-        (
-            IndirectList<face>(slaveMesh.faces(), slaveToMesh),
-            slaveMesh.points()
-        )
-    );
-
-
-    if (perfectMatch)
-    {
-        // Faces are perfectly aligned but probably not ordered.
-        perfectPointMatch(absTol);
-    }
-    else
-    {
-        // Slave faces are subdivision of master face. Faces not ordered.
-        subDivisionMatch(slaveMesh, false, absTol);
-    }
-
-    if (debug)
-    {
-        writePointsFaces();
-    }
-}
-
-
 Foam::faceCoupleInfo::faceCoupleInfo
 (
     const polyMesh& masterMesh,
     const labelList& masterAddressing,
     const polyMesh& slaveMesh,
-    const labelList& slaveAddressing,
-    const scalar absTol,
-    const bool perfectMatch,
-    const bool orderedFaces,
-    const bool patchDivision
+    const labelList& slaveAddressing
 )
 :
     masterPatchPtr_
@@ -2036,21 +188,16 @@ Foam::faceCoupleInfo::faceCoupleInfo
             slaveMesh.points()
         )
     ),
-    cutPoints_(0),
-    cutFacesPtr_(nullptr),
-    cutToMasterFaces_(0),
-    masterToCutPoints_(0),
-    cutToSlaveFaces_(0),
-    slaveToCutPoints_(0),
-    cutEdgeToPoints_(0)
+    nCouplePoints_(0),
+    masterToCouplePoints_(),
+    slaveToCouplePoints_()
 {
-    if (perfectMatch && (masterAddressing.size() != slaveAddressing.size()))
+    if (masterAddressing.size() != slaveAddressing.size())
     {
         FatalErrorInFunction
-            << "Perfect match specified but number of master and slave faces"
-            << " differ." << endl
+            << "Number of master and slave faces differ." << endl
             << "master:" << masterAddressing.size()
-            << "  slave:" << slaveAddressing.size()
+            << " slave:" << slaveAddressing.size()
             << abort(FatalError);
     }
 
@@ -2065,6 +212,7 @@ Foam::faceCoupleInfo::faceCoupleInfo
             << "Faces to be coupled have to be boundary faces."
             << abort(FatalError);
     }
+
     if
     (
         slaveAddressing.size()
@@ -2077,19 +225,46 @@ Foam::faceCoupleInfo::faceCoupleInfo
             << abort(FatalError);
     }
 
+    // Initialise point addressing
+    nCouplePoints_ = 0;
+    masterToCouplePoints_ = labelList(masterPatch().nPoints(), -1);
+    slaveToCouplePoints_ = labelList(slavePatch().nPoints(), -1);
 
-    if (perfectMatch && orderedFaces)
+    // Slave points number around the master faces in the reverse direction and
+    // have the same first point
+    forAll(masterPatch(), coupleFacei)
     {
-        perfectOrderedPointMatch();
-    }
-    else if (perfectMatch)
-    {
-        perfectPointMatch(absTol);
-    }
-    else
-    {
-        // Slave faces are subdivision of master face. Faces not ordered.
-        subDivisionMatch(slaveMesh, patchDivision, absTol);
+        const face& masterF = masterPatch().localFaces()[coupleFacei];
+        const face& slaveF = slavePatch().localFaces()[coupleFacei];
+
+        label slaveFp = 0;
+
+        forAll(masterF, masterFp)
+        {
+            const label masterPointi = masterF[masterFp];
+            const label slavePointi = slaveF[slaveFp];
+
+            // If this point already has a coupled point index then use it,
+            // else create a new one
+            label couplePointi = -1;
+            if (masterToCouplePoints_[masterPointi] != -1)
+            {
+                couplePointi = masterToCouplePoints_[masterPointi];
+            }
+            else if (slaveToCouplePoints_[slavePointi] != -1)
+            {
+                couplePointi = slaveToCouplePoints_[slavePointi];
+            }
+            else
+            {
+                couplePointi = nCouplePoints_ ++;
+            }
+
+            masterToCouplePoints_[masterPointi] = couplePointi;
+            slaveToCouplePoints_[slavePointi] = couplePointi;
+
+            slaveFp = slaveF.rcIndex(slaveFp);
+        }
     }
 
     if (debug)
@@ -2105,53 +280,4 @@ Foam::faceCoupleInfo::~faceCoupleInfo()
 {}
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-Foam::labelList Foam::faceCoupleInfo::faceLabels(const polyPatch& pp)
-{
-    labelList faces(pp.size());
-
-    label facei = pp.start();
-
-    forAll(pp, i)
-    {
-        faces[i] = facei++;
-    }
-    return faces;
-}
-
-
-Foam::Map<Foam::label> Foam::faceCoupleInfo::makeMap(const labelList& lst)
-{
-    Map<label> map(lst.size());
-
-    forAll(lst, i)
-    {
-        if (lst[i] != -1)
-        {
-            map.insert(i, lst[i]);
-        }
-    }
-    return map;
-}
-
-
-Foam::Map<Foam::labelList> Foam::faceCoupleInfo::makeMap
-(
-    const labelListList& lst
-)
-{
-    Map<labelList> map(lst.size());
-
-    forAll(lst, i)
-    {
-        if (lst[i].size())
-        {
-            map.insert(i, lst[i]);
-        }
-    }
-    return map;
-}
-
-
 // ************************************************************************* //
diff --git a/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.H b/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.H
index 49c687643f..9ae9ae59c0 100644
--- a/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.H
+++ b/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,108 +26,12 @@ Class
 
 Description
     Container for information needed to couple to meshes. When constructed
-    from two meshes and a geometric tolerance finds the corresponding
-    boundary faces.
-
-    The information it keeps is the set of faces&points (cutFaces,
-    cutPoints) that should replace a set of faces on the master
-    (masterPatch) and a set of faces on the slave (slavePatch)
-
-
-    Uses same tolerance to match faces and points on matched faces since
-    they both originate from the same points and the tolerance usually
-    comes from writing these points with limited precision (6 by default)
-
-    -# Perfect match:
-       - one-to-one match for faces and points.
-       - the cut is always the 'most connected' of the master and slave so
-         multiple master or slave points might point to the same cut point.
-
-    \verbatim
-    e.g. master:
-
-        +--+
-        |  |
-        |  |
-        +--+
-            +--+
-            |  |
-            |  |
-            +--+
-    slave:
-            +--+
-            |  |
-            |  |
-            +--+
-        +--+
-        |  |
-        |  |
-        +--+
-    \endverbatim
-    adding both together creates a singly connected 2x2 cavity so suddenly
-    the duplicate master points and the duplicate slave points all become
-    a single cut point.
-
-
-    -# Subdivision match:
-       - Can be constructed from slave being subdivision of master with the
-         polyPatch constructor.
-       - Does not include above shared-point detection!
-
-    Notes on multiple slave faces per master:
-
-    As long as
-    - all master edges are present in slave
-    - slave can have extra edges/points/faces BUT all subfaces have to have
-      at least one point on a maste face.
-
-    \verbatim
-    So master:
-    +-------+
-    |       |
-    |       |
-    |       |
-    |       |
-    |       |
-    |       |
-    |       |
-    +-------+
-
-    slave:
-    +---+---+
-    |\  |  /|
-    | \ | / |
-    |  \|/  |
-    +---+---+
-    |  /|\  |
-    | / | \ |
-    |/  |  \|
-    +---+---+
-    is ok.
-    \endverbatim
-
-    For this kind of matching the order is:
-    - match cutpoint to masterpoint
-    - find those cutEdges that align with a master edge. This gives two sets
-      of cut edges: those that have a master equivalent ('border edges') and
-      those that don't ('internal edges'). The border edges now divide the
-      cutFaces into regions with the same masterFace correspondence.
-    - find cutFaces that are fully determined by the border edges they use.
-    - all cutFaces that are connected through an internal edge have the same
-      master face.
-
-
-    Note: matching refined faces onto master is a bit dodgy and will probably
-    only work for unwarped faces. Also it will fail if e.g. face is split
-    into 3x3 since then middle face has no point/edge in common with master.
-    (problem is in face matching (findSlavesCoveringMaster), probably
-     point/edge matching might just work)
-
+    from two meshes and a list of coupled faces returns the mapping between
+    points.
 
 SourceFiles
     faceCoupleInfo.C
 
-
 \*---------------------------------------------------------------------------*/
 
 #ifndef faceCoupleInfo_H
@@ -135,20 +39,13 @@ SourceFiles
 
 #include "pointField.H"
 #include "indirectPrimitivePatch.H"
-#include "primitiveFacePatch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-typedef HashTable<labelList, edge, Hash<edge>> edgeLookup;
-
-
 // Forward declaration of classes
-class face;
-class primitiveMesh;
-class polyPatch;
 class polyMesh;
 
 /*---------------------------------------------------------------------------*\
@@ -159,49 +56,20 @@ class faceCoupleInfo
 {
     // Private Data
 
-        //- Angle matching tolerance.
-        static const scalar angleTol_;
-
         //- Master patch
         autoPtr<indirectPrimitivePatch> masterPatchPtr_;
 
         //- Slave patch
         autoPtr<indirectPrimitivePatch> slavePatchPtr_;
 
+        //- The number of coupled points
+        label nCouplePoints_;
 
-        //- Description of cut.
-        //  - Cut is the matching area between the slave
-        //  and the master.
-        //  - cut is the finest of master and slave. It can never be
-        //  coarser than either one of them. (so face addressing we keep is
-        //  cut-to-master and cut-to-slave)
-        //  - multiple master or slave points can end up becoming one cut point
-        //  (so point addressing we keep is master-to-cut and slave-to-cut)
+        //- Master-to-couple point addressing
+        labelList masterToCouplePoints_;
 
-        //  Cut consists of faces and points (note: could be expressed as some
-        //  kind of PrimitivePatch which holds points instead of reference to
-        //  them)
-        //  Orientation of cutFaces should be same as masterFaces!
-        pointField cutPoints_;
-        autoPtr<primitiveFacePatch> cutFacesPtr_;
-
-        //- Additional point coupling information. Is between points on
-        //  boundary of both meshes.
-
-        // Addressing to/from cut
-
-            //- master
-            labelList cutToMasterFaces_;
-            labelList masterToCutPoints_;
-
-            //- slave
-            labelList cutToSlaveFaces_;
-            labelList slaveToCutPoints_;
-
-            //- For edges originating from splitting of edges:
-            //  given the two endpoints of the unsplit edge give the list
-            //  of in between vertices
-            edgeLookup cutEdgeToPoints_;
+        //- Slave-to-couple point addressing
+        labelList slaveToCouplePoints_;
 
 
     // Private Member Functions
@@ -233,8 +101,7 @@ class faceCoupleInfo
             (
                 const fileName& fName,
                 const edgeList& edges,
-                const pointField& points,
-                const bool compact = true
+                const pointField& points
             );
 
             //- Write edges
@@ -249,136 +116,6 @@ class faceCoupleInfo
             //  as .obj files.
             void writePointsFaces() const;
 
-            //- Write connections between corresponding edges as .obj files.
-            void writeEdges(const labelList&, const labelList&) const;
-
-
-        // Edge handling/matching
-
-            //- Find corresponding edges on patch when having only a map for
-            //  the points.
-            labelList findMappedEdges
-            (
-                const edgeList& edges,
-                const labelList& pointMap,
-                const indirectPrimitivePatch&
-            );
-
-            //- Check if edge on slavePatch corresponds to an edge between faces
-            //  in two different polyPatches on the mesh.
-            bool regionEdge(const polyMesh&, const label slaveEdgeI) const;
-
-            //- Finds edge connected to point most aligned with master edge.
-            label mostAlignedCutEdge
-            (
-                const bool report,
-                const polyMesh& slaveMesh,
-                const bool patchDivision,
-                const labelList& cutToMasterEdges,
-                const labelList& cutToSlaveEdges,
-                const label pointi,
-                const label edgeStart,
-                const label edgeEnd
-            ) const;
-
-            //- From (many-to-one) map of cut edges to master edges determine
-            //  points in between. I.e. just string up the edges. Stores this
-            //  all on cutEdgeToPoints_
-            void setCutEdgeToPoints(const labelList& cutToMasterEdges);
-
-        // Face matching
-
-            //- Matches two faces.
-            //  Determines rotation for f1 to match up with f0,
-            //  i.e. the index in f0 of the first point of f1.
-            static label matchFaces
-            (
-                const scalar absTol,
-                const pointField& points0,
-                const face& f0,
-                const pointField& points1,
-                const face& f1,
-                const bool sameOrientation
-            );
-
-            //- Matches points on patch to points on cut.
-            static bool matchPointsThroughFaces
-            (
-                const scalar absTol,
-                const pointField& cutPoints,
-                const faceList& cutFaces,
-                const pointField& patchPoints,
-                const faceList& patchFaces,
-                const bool sameOrientation,
-
-                labelList& patchToCutPoints,// patch to (uncompacted) cut points
-                labelList& cutToCompact,    // compaction list
-                labelList& compactToCut     // compaction list
-            );
-
-            //- Returns max distance to masterF of any point on cutF.
-            static scalar maxDistance
-            (
-                const face& cutF,
-                const pointField& cutPoints,
-                const face& masterF,
-                const pointField& masterPoints
-            );
-
-            //- Finds matching (boundary)face centres.
-            //  Since faces identical uses geometric match on face centres.
-            static void findPerfectMatchingFaces
-            (
-                const primitiveMesh& mesh0,
-                const primitiveMesh& mesh1,
-                const scalar absTol,
-
-                labelList& mesh0Faces,
-                labelList& mesh1Faces
-            );
-
-            //- Find matching (boundary)faces. Matching if slave is on top of
-            //  master face (slaves is subdivision of master)
-            static void findSlavesCoveringMaster
-            (
-                const primitiveMesh& mesh0,
-                const primitiveMesh& mesh1,
-                const scalar absTol,
-
-                labelList& mesh0Faces,
-                labelList& mesh1Faces
-            );
-
-            //- Grow cutToMasterFace across 'internal' edges.
-            label growCutFaces(const labelList&, Map<labelList>&);
-
-            void checkMatch(const labelList& cutToMasterEdges) const;
-
-            //- Gets a list of cutFaces (that use a master edge) and the
-            // candidate master faces.
-            // Checks among these master faces if there is only one remaining
-            // unmatched one.
-            label matchEdgeFaces(const labelList&, Map<labelList>& candidates);
-
-            //- Gets a list of cutFaces (that use a master edge) and the
-            //  candidate master faces.
-            //  Finds most aligned master face.
-            label geometricMatchEdgeFaces(Map<labelList>& candidates);
-
-        //- Find point and edge correspondence for perfect matching faces
-        void perfectOrderedPointMatch();
-
-        //- Find point and edge correspondence for perfect matching faces
-        void perfectPointMatch(const scalar absTol);
-
-        //- Find point and edge correspondence for slaves being subdivision of
-        //  master.
-        void subDivisionMatch
-        (
-            const polyMesh& slaveMesh,
-            const bool patchDivision,
-            const scalar absTol
-        );
 
 public:
 
@@ -388,43 +125,16 @@ public:
 
     // Constructors
 
-        //- Construct from two meshes and absolute tolerance.
-        //  Finds out matches geometrically. No checking for nonsense match.
-        //  Tolerance is absolute one so use with care.
-        //  perfectMatch : each point/edge/face has corresponding point on other
-        //                 side
-        //                 if this is false then assumes slave is subdivision.
-        //                 Matching then will work only for non-warped faces
-        //                 since does nearest-to-face comparison with absTol.
-        faceCoupleInfo
-        (
-            const polyMesh& mesh0,
-            const polyMesh& mesh1,
-            const scalar absTol,
-            const bool perfectMatch
-        );
-
-        //- Construct from meshes and subset of mesh faces
-        //  (i.e. indirectPrimitivePatch addressing)
-        //  All faces in patch are considered matched (but don't have to be
-        //  ordered)
-        //  perfectMatch : each point/edge/face has corresponding point on other
-        //                 side
-        //  orderedFaces : faces in patch are ordered (so masterAddressing[i]
-        //  matches slaveAddressing[i])
-        //  patchDivision: faces in slave mesh that originate from the
-        //  same master face have the same patch. Used by some triangulation
-        //  methods.
+        //- Construct from meshes and subset of mesh faces (i.e.,
+        //  indirectPrimitivePatch addressing). Faces should be supplied in
+        //  order and the points within them should be ordered (i.e., coupled
+        //  faces should have the same zero-point and opposite orientations).
         faceCoupleInfo
         (
             const polyMesh& masterMesh,
             const labelList& masterAddressing,
             const polyMesh& slaveMesh,
-            const labelList& slaveAddressing,
-            const scalar absTol,
-            const bool perfectMatch,
-            const bool orderedFaces,
-            const bool patchDivision
+            const labelList& slaveAddressing
         );
 
 
@@ -432,19 +142,8 @@ public:
     ~faceCoupleInfo();
 
 
-
     // Member Functions
 
-        //- Utility functions
-
-            //- Get patch face labels
-            static labelList faceLabels(const polyPatch&);
-
-            //- Create Map from List
-            static Map<label> makeMap(const labelList&);
-            static Map<labelList> makeMap(const labelListList&);
-
-
         // Access
 
             //- Addressing engine for coupled faces on mesh0
@@ -459,49 +158,40 @@ public:
                 return slavePatchPtr_();
             }
 
-            //- Addressing engine for combined set of faces.
-            const primitiveFacePatch& cutFaces() const
+            //- Access master-to-couple point addressing
+            const labelList& masterToCouplePoints() const
             {
-                return cutFacesPtr_();
+                return masterToCouplePoints_;
             }
 
-            //- Points for combined set of faces.
-            const pointField& cutPoints() const
+            //- Access slave-to-couple point addressing
+            const labelList& slaveToCouplePoints() const
             {
-                return cutPoints_;
+                return slaveToCouplePoints_;
             }
 
-
-        // Addressing from meshes to cut and vice versa.
-
-            //- Master face for every face on cut. Will always be at least
-            //  one but there might be multiple cut faces pointing to the same
-            //  master
-            const labelList& cutToMasterFaces() const
+            //- Generate couple-to-master point addressing
+            labelListList coupleToMasterPoints() const
             {
-                return cutToMasterFaces_;
-            }
-            const labelList& masterToCutPoints() const
-            {
-                return masterToCutPoints_;
+                return invertOneToMany(nCouplePoints_, masterToCouplePoints_);
             }
 
-            const labelList& cutToSlaveFaces() const
+            //- Generate couple-to-slave point addressing
+            labelListList coupleToSlavePoints() const
             {
-                return cutToSlaveFaces_;
-            }
-            const labelList& slaveToCutPoints() const
-            {
-                return slaveToCutPoints_;
+                return invertOneToMany(nCouplePoints_, slaveToCouplePoints_);
             }
 
-            //- From two cut points (original edge) to list of inserted
-            //  points
-            const edgeLookup& cutEdgeToPoints() const
+            //- Return a coupled face
+            const face coupleFace(const label coupleFacei) const
             {
-                return cutEdgeToPoints_;
+                return
+                    renumber
+                    (
+                        masterToCouplePoints(),
+                        masterPatch().localFaces()[coupleFacei]
+                    );
             }
-
  };
 
 
diff --git a/src/dynamicMesh/polyMeshAdder/polyMeshAdder.C b/src/dynamicMesh/polyMeshAdder/polyMeshAdder.C
index 48cfe65159..4c423f0ef9 100644
--- a/src/dynamicMesh/polyMeshAdder/polyMeshAdder.C
+++ b/src/dynamicMesh/polyMeshAdder/polyMeshAdder.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,7 +31,6 @@ License
 #include "SortableList.H"
 #include "Time.H"
 #include "globalMeshData.H"
-#include "mergePoints.H"
 #include "polyModifyFace.H"
 #include "polyRemovePoint.H"
 #include "polyTopoChange.H"
@@ -395,81 +394,6 @@ Foam::labelList Foam::polyMeshAdder::getFaceOrder
 }
 
 
-// Extends face f with split points. cutEdgeToPoints gives for every
-// edge the points introduced in between the endpoints.
-void Foam::polyMeshAdder::insertVertices
-(
-    const edgeLookup& cutEdgeToPoints,
-    const Map<label>& meshToMaster,
-    const labelList& masterToCutPoints,
-    const face& masterF,
-
-    DynamicList<label>& workFace,
-    face& allF
-)
-{
-    workFace.clear();
-
-    // Check any edge for being cut (check on the cut so takes account
-    // for any point merging on the cut)
-
-    forAll(masterF, fp)
-    {
-        label v0 = masterF[fp];
-        label v1 = masterF.nextLabel(fp);
-
-        // Copy existing face point
-        workFace.append(allF[fp]);
-
-        // See if any edge between v0,v1
-
-        Map<label>::const_iterator v0Fnd = meshToMaster.find(v0);
-        if (v0Fnd != meshToMaster.end())
-        {
-            Map<label>::const_iterator v1Fnd = meshToMaster.find(v1);
-            if (v1Fnd != meshToMaster.end())
-            {
-                // Get edge in cutPoint numbering
-                edge cutEdge
-                (
-                    masterToCutPoints[v0Fnd()],
-                    masterToCutPoints[v1Fnd()]
-                );
-
-                edgeLookup::const_iterator iter = cutEdgeToPoints.find(cutEdge);
-
-                if (iter != cutEdgeToPoints.end())
-                {
-                    const edge& e = iter.key();
-                    const labelList& addedPoints = iter();
-
-                    // cutPoints first in allPoints so no need for renumbering
-                    if (e[0] == cutEdge[0])
-                    {
-                        forAll(addedPoints, i)
-                        {
-                            workFace.append(addedPoints[i]);
-                        }
-                    }
-                    else
-                    {
-                        forAllReverse(addedPoints, i)
-                        {
-                            workFace.append(addedPoints[i]);
-                        }
-                    }
-                }
-            }
-        }
-    }
-
-    if (workFace.size() != allF.size())
-    {
-        allF.transfer(workFace);
-    }
-}
-
-
 // Adds primitives (cells, faces, points)
 // Cells:
 //  - all of mesh0
@@ -511,7 +435,6 @@ void Foam::polyMeshAdder::mergePrimitives
     const polyBoundaryMesh& patches0 = mesh0.boundaryMesh();
     const polyBoundaryMesh& patches1 = mesh1.boundaryMesh();
 
-    const primitiveFacePatch& cutFaces = coupleInfo.cutFaces();
     const indirectPrimitivePatch& masterPatch = coupleInfo.masterPatch();
     const indirectPrimitivePatch& slavePatch = coupleInfo.slavePatch();
 
@@ -528,53 +451,37 @@ void Foam::polyMeshAdder::mergePrimitives
     from1ToAllPoints.setSize(mesh1.nPoints());
     from1ToAllPoints = -1;
 
-    // Copy coupled points (on cut)
+    // Copy coupled points
     {
-        const pointField& cutPoints = coupleInfo.cutPoints();
-
-        // const labelListList& cutToMasterPoints =
-        //   coupleInfo.cutToMasterPoints();
-        labelListList cutToMasterPoints
+        const labelListList coupleToMasterPoints
         (
-            invertOneToMany
-            (
-                cutPoints.size(),
-                coupleInfo.masterToCutPoints()
-            )
+            coupleInfo.coupleToMasterPoints()
+        );
+        const labelListList coupleToSlavePoints
+        (
+            coupleInfo.coupleToSlavePoints()
         );
 
-        // const labelListList& cutToSlavePoints =
-        //    coupleInfo.cutToSlavePoints();
-        labelListList cutToSlavePoints
-        (
-            invertOneToMany
-            (
-                cutPoints.size(),
-                coupleInfo.slaveToCutPoints()
-            )
-        );
-
-        forAll(cutPoints, i)
+        forAll(coupleToMasterPoints, couplePointi)
         {
-            allPoints[allPointi] = cutPoints[i];
-
-            // Mark all master and slave points referring to this point.
-
-            const labelList& masterPoints = cutToMasterPoints[i];
+            const labelList& masterPoints = coupleToMasterPoints[couplePointi];
 
             forAll(masterPoints, j)
             {
                 label mesh0Pointi = masterPatch.meshPoints()[masterPoints[j]];
                 from0ToAllPoints[mesh0Pointi] = allPointi;
+                allPoints[allPointi] = mesh0.points()[mesh0Pointi];
             }
 
-            const labelList& slavePoints = cutToSlavePoints[i];
+            const labelList& slavePoints = coupleToSlavePoints[couplePointi];
 
             forAll(slavePoints, j)
             {
                 label mesh1Pointi = slavePatch.meshPoints()[slavePoints[j]];
                 from1ToAllPoints[mesh1Pointi] = allPointi;
+                allPoints[allPointi] = mesh1.points()[mesh1Pointi];
             }
+
             allPointi++;
         }
     }
@@ -635,14 +542,9 @@ void Foam::polyMeshAdder::mergePrimitives
 
     // Copy coupled faces. Every coupled face has an equivalent master and
     // slave. Also uncount as boundary faces all the newly coupled faces.
-    const labelList& cutToMasterFaces = coupleInfo.cutToMasterFaces();
-    const labelList& cutToSlaveFaces = coupleInfo.cutToSlaveFaces();
-
-    forAll(cutFaces, i)
+    forAll(masterPatch, coupleFacei)
     {
-        label masterFacei = cutToMasterFaces[i];
-
-        label mesh0Facei = masterPatch.addressing()[masterFacei];
+        const label mesh0Facei = masterPatch.addressing()[coupleFacei];
 
         if (from0ToAllFaces[mesh0Facei] == -1)
         {
@@ -654,9 +556,7 @@ void Foam::polyMeshAdder::mergePrimitives
             nFacesPerPatch[patch0]--;
         }
 
-        label slaveFacei = cutToSlaveFaces[i];
-
-        label mesh1Facei = slavePatch.addressing()[slaveFacei];
+        const label mesh1Facei = slavePatch.addressing()[coupleFacei];
 
         if (from1ToAllFaces[mesh1Facei] == -1)
         {
@@ -668,7 +568,7 @@ void Foam::polyMeshAdder::mergePrimitives
 
         // Copy cut face (since cutPoints are copied first no renumbering
         // necessary)
-        allFaces[allFacei] = cutFaces[i];
+        allFaces[allFacei] = coupleInfo.coupleFace(coupleFacei);
         allOwner[allFacei] = mesh0.faceOwner()[mesh0Facei];
         allNeighbour[allFacei] = mesh1.faceOwner()[mesh1Facei] + mesh0.nCells();
 
@@ -750,85 +650,6 @@ void Foam::polyMeshAdder::mergePrimitives
     allOwner.setSize(allFacei);
     allNeighbour.setSize(allFacei);
 
-
-    // So now we have all ok for one-to-one mapping.
-    // For split slace faces:
-    // - mesh consistent with slave side
-    // - mesh not consistent with owner side. It is not zipped up, the
-    //   original faces need edges split.
-
-    // Use brute force to prevent having to calculate addressing:
-    // - get map from master edge to split edges.
-    // - check all faces to find any edge that is split.
-    {
-        // From two cut-points to labels of cut-points in between.
-        // (in order: from e[0] to e[1]
-        const edgeLookup& cutEdgeToPoints = coupleInfo.cutEdgeToPoints();
-
-        // Get map of master face (in mesh labels) that are in cut. These faces
-        // do not need to be renumbered.
-        labelHashSet masterCutFaces(cutToMasterFaces.size());
-        forAll(cutToMasterFaces, i)
-        {
-            label meshFacei = masterPatch.addressing()[cutToMasterFaces[i]];
-
-            masterCutFaces.insert(meshFacei);
-        }
-
-        DynamicList<label> workFace(100);
-
-        forAll(from0ToAllFaces, face0)
-        {
-            if (!masterCutFaces.found(face0))
-            {
-                label allFacei = from0ToAllFaces[face0];
-
-                insertVertices
-                (
-                    cutEdgeToPoints,
-                    masterPatch.meshPointMap(),
-                    coupleInfo.masterToCutPoints(),
-                    mesh0.faces()[face0],
-
-                    workFace,
-                    allFaces[allFacei]
-                );
-            }
-        }
-
-        // Same for slave face
-
-        labelHashSet slaveCutFaces(cutToSlaveFaces.size());
-        forAll(cutToSlaveFaces, i)
-        {
-            label meshFacei = slavePatch.addressing()[cutToSlaveFaces[i]];
-
-            slaveCutFaces.insert(meshFacei);
-        }
-
-        forAll(from1ToAllFaces, face1)
-        {
-            if (!slaveCutFaces.found(face1))
-            {
-                label allFacei = from1ToAllFaces[face1];
-
-                insertVertices
-                (
-                    cutEdgeToPoints,
-                    slavePatch.meshPointMap(),
-                    coupleInfo.slaveToCutPoints(),
-                    mesh1.faces()[face1],
-
-                    workFace,
-                    allFaces[allFacei]
-                );
-            }
-        }
-    }
-
-    // Now we have a full facelist and owner/neighbour addressing.
-
-
     // Cells
     // ~~~~~
 
@@ -1603,7 +1424,7 @@ Foam::autoPtr<Foam::polyMesh> Foam::polyMeshAdder::add
             fromAllTo1Patches,
             mesh0.nInternalFaces()
           + mesh1.nInternalFaces()
-          + coupleInfo.cutFaces().size(),
+          + coupleInfo.masterPatch().size(),
             nFaces,
 
             from0ToAllPatches,
@@ -1998,310 +1819,4 @@ Foam::autoPtr<Foam::mapAddedPolyMesh> Foam::polyMeshAdder::add
 }
 
 
-Foam::Map<Foam::label> Foam::polyMeshAdder::findSharedPoints
-(
-    const polyMesh& mesh,
-    const scalar mergeDist
-)
-{
-    const labelList& sharedPointLabels = mesh.globalData().sharedPointLabels();
-    const labelList& sharedPointAddr = mesh.globalData().sharedPointAddr();
-
-    // Because of adding the missing pieces e.g. when redistributing a mesh
-    // it can be that there are multiple points on the same processor that
-    // refer to the same shared point.
-
-    // Invert point-to-shared addressing
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    Map<labelList> sharedToMesh(sharedPointLabels.size());
-
-    label nMultiple = 0;
-
-    forAll(sharedPointLabels, i)
-    {
-        label pointi = sharedPointLabels[i];
-
-        label sharedI = sharedPointAddr[i];
-
-        Map<labelList>::iterator iter = sharedToMesh.find(sharedI);
-
-        if (iter != sharedToMesh.end())
-        {
-            // sharedI already used by other point. Add this one.
-
-            nMultiple++;
-
-            labelList& connectedPointLabels = iter();
-
-            label sz = connectedPointLabels.size();
-
-            // Check just to make sure.
-            if (findIndex(connectedPointLabels, pointi) != -1)
-            {
-                FatalErrorInFunction
-                    << "Duplicate point in sharedPoint addressing." << endl
-                    << "When trying to add point " << pointi << " on shared "
-                    << sharedI  << " with connected points "
-                    << connectedPointLabels
-                    << abort(FatalError);
-            }
-
-            connectedPointLabels.setSize(sz+1);
-            connectedPointLabels[sz] = pointi;
-        }
-        else
-        {
-            sharedToMesh.insert(sharedI, labelList(1, pointi));
-        }
-    }
-
-
-    // Assign single master for every shared with multiple geometric points
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    Map<label> pointToMaster(nMultiple);
-
-    forAllConstIter(Map<labelList>, sharedToMesh, iter)
-    {
-        const labelList& connectedPointLabels = iter();
-
-        // Pout<< "For shared:" << iter.key()
-        //    << " found points:" << connectedPointLabels
-        //    << " at coords:"
-        //    <<  pointField(mesh.points(), connectedPointLabels) << endl;
-
-        if (connectedPointLabels.size() > 1)
-        {
-            const pointField connectedPoints
-            (
-                mesh.points(),
-                connectedPointLabels
-            );
-
-            labelList toMergedPoints;
-            label nUnique = Foam::mergePoints
-            (
-                connectedPoints,
-                mergeDist,
-                false,
-                toMergedPoints
-            );
-
-            if (nUnique < connectedPoints.size())
-            {
-                // Invert toMergedPoints
-                const labelListList mergeSets
-                (
-                    invertOneToMany
-                    (
-                        nUnique,
-                        toMergedPoints
-                    )
-                );
-
-                // Find master for valid merges
-                forAll(mergeSets, setI)
-                {
-                    const labelList& mergeSet = mergeSets[setI];
-
-                    if (mergeSet.size() > 1)
-                    {
-                        // Pick lowest numbered point
-                        label masterPointi = labelMax;
-
-                        forAll(mergeSet, i)
-                        {
-                            label pointi = connectedPointLabels[mergeSet[i]];
-
-                            masterPointi = min(masterPointi, pointi);
-                        }
-
-                        forAll(mergeSet, i)
-                        {
-                            label pointi = connectedPointLabels[mergeSet[i]];
-
-                            // Pout<< "Merging point " << pointi
-                            //    << " at " << mesh.points()[pointi]
-                            //    << " into master point "
-                            //    << masterPointi
-                            //    << " at " << mesh.points()[masterPointi]
-                            //    << endl;
-
-                            pointToMaster.insert(pointi, masterPointi);
-                        }
-                    }
-                }
-            }
-        }
-    }
-
-    //- Old: geometric merging. Causes problems for two close shared points.
-    // labelList sharedToMerged;
-    // label nUnique = Foam::mergePoints
-    //(
-    //    pointField
-    //    (
-    //        mesh.points(),
-    //        sharedPointLabels
-    //    ),
-    //    mergeDist,
-    //    false,
-    //    sharedToMerged
-    //);
-    //
-    //// Find out which sets of points get merged and create a map from
-    //// mesh point to unique point.
-    //
-    // Map<label> pointToMaster(10*sharedToMerged.size());
-    //
-    // if (nUnique < sharedPointLabels.size())
-    //{
-    //    labelListList mergeSets
-    //    (
-    //        invertOneToMany
-    //        (
-    //            sharedToMerged.size(),
-    //            sharedToMerged
-    //        )
-    //    );
-    //
-    //    label nMergeSets = 0;
-    //
-    //    forAll(mergeSets, setI)
-    //    {
-    //        const labelList& mergeSet = mergeSets[setI];
-    //
-    //        if (mergeSet.size() > 1)
-    //        {
-    //            // Take as master the shared point with the lowest mesh
-    //            // point label. (rather arbitrarily - could use max or
-    //            // any other one of the points)
-    //
-    //            nMergeSets++;
-    //
-    //            label masterI = labelMax;
-    //
-    //            forAll(mergeSet, i)
-    //            {
-    //                label sharedI = mergeSet[i];
-    //
-    //                masterI = min(masterI, sharedPointLabels[sharedI]);
-    //            }
-    //
-    //            forAll(mergeSet, i)
-    //            {
-    //                label sharedI = mergeSet[i];
-    //
-    //                pointToMaster.insert(sharedPointLabels[sharedI], masterI);
-    //            }
-    //        }
-    //    }
-    //
-    //    // if (debug)
-    //    //{
-    //    //    Pout<< "polyMeshAdder : merging:"
-    //    //        << pointToMaster.size() << " into " << nMergeSets
-    //    //        << " sets." << endl;
-    //    //}
-    //}
-
-    return pointToMaster;
-}
-
-
-void Foam::polyMeshAdder::mergePoints
-(
-    const polyMesh& mesh,
-    const Map<label>& pointToMaster,
-    polyTopoChange& meshMod
-)
-{
-    // Remove all non-master points.
-    forAll(mesh.points(), pointi)
-    {
-        Map<label>::const_iterator iter = pointToMaster.find(pointi);
-
-        if (iter != pointToMaster.end())
-        {
-            if (iter() != pointi)
-            {
-                meshMod.removePoint(pointi, iter());
-            }
-        }
-    }
-
-    // Modify faces for points. Note: could use pointFaces here but want to
-    // avoid addressing calculation.
-    const faceList& faces = mesh.faces();
-
-    forAll(faces, facei)
-    {
-        const face& f = faces[facei];
-
-        bool hasMerged = false;
-
-        forAll(f, fp)
-        {
-            label pointi = f[fp];
-
-            Map<label>::const_iterator iter = pointToMaster.find(pointi);
-
-            if (iter != pointToMaster.end())
-            {
-                if (iter() != pointi)
-                {
-                    hasMerged = true;
-                    break;
-                }
-            }
-        }
-
-        if (hasMerged)
-        {
-            face newF(f);
-
-            forAll(f, fp)
-            {
-                label pointi = f[fp];
-
-                Map<label>::const_iterator iter = pointToMaster.find(pointi);
-
-                if (iter != pointToMaster.end())
-                {
-                    newF[fp] = iter();
-                }
-            }
-
-            label patchID = mesh.boundaryMesh().whichPatch(facei);
-            label nei = (patchID == -1 ? mesh.faceNeighbour()[facei] : -1);
-            label zoneID = mesh.faceZones().whichZone(facei);
-            bool zoneFlip = false;
-
-            if (zoneID >= 0)
-            {
-                const faceZone& fZone = mesh.faceZones()[zoneID];
-                zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
-            }
-
-            meshMod.setAction
-            (
-                polyModifyFace
-                (
-                    newF,                       // modified face
-                    facei,                      // label of face
-                    mesh.faceOwner()[facei],    // owner
-                    nei,                        // neighbour
-                    false,                      // face flip
-                    patchID,                    // patch for face
-                    false,                      // remove from zone
-                    zoneID,                     // zone for face
-                    zoneFlip                    // face flip in zone
-                )
-            );
-        }
-    }
-}
-
-
 // ************************************************************************* //
diff --git a/src/dynamicMesh/polyMeshAdder/polyMeshAdder.H b/src/dynamicMesh/polyMeshAdder/polyMeshAdder.H
index aa5ec39c7c..932a70d687 100644
--- a/src/dynamicMesh/polyMeshAdder/polyMeshAdder.H
+++ b/src/dynamicMesh/polyMeshAdder/polyMeshAdder.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -59,7 +59,6 @@ class polyTopoChange;
 
 class polyMeshAdder
 {
-
 private:
 
     // Private Member Functions
@@ -89,6 +88,7 @@ private:
 
         //- Get starts of patches
         static labelList getPatchStarts(const polyBoundaryMesh&);
+
         //- Get sizes of patches
         static labelList getPatchSizes(const polyBoundaryMesh&);
 
@@ -117,19 +117,6 @@ private:
             const labelList& neighbour
         );
 
-        //- Extends face f with split points. cutEdgeToPoints gives for every
-        // edge the points introduced in between the endpoints.
-        static void insertVertices
-        (
-            const edgeLookup& cutEdgeToPoints,
-            const Map<label>& meshToMaster,
-            const labelList& masterToCutPoints,
-            const face& masterFace,
-
-            DynamicList<label>& workFace,
-            face& f
-        );
-
         //- Merges primitives of two meshes.
         static void mergePrimitives
         (
@@ -250,7 +237,6 @@ public:
 
     // Member Functions
 
-
         //- Add two polyMeshes. Returns new polyMesh and map construct.
         static autoPtr<polyMesh> add
         (
@@ -269,32 +255,6 @@ public:
             const faceCoupleInfo& coupleInfo,
             const bool validBoundary = true
         );
-
-        // Point merging
-
-            //- Find topologically and geometrically shared points.
-            //
-            //  - should only be called for parallel correct mesh
-            //    (since uses mesh.globalData)
-            //  - returns Map from point to master point (all in mesh point
-            //    labels) for any sets of points that need to be merged.
-            static Map<label> findSharedPoints
-            (
-                const polyMesh&,
-                const scalar mergeTol
-            );
-
-            //- Helper: Merge points.
-            //
-            //  - Gets map from point to destination point
-            //  - Removes all points that don't map to themselves
-            //  - Modifies all faces that use the points to be removed.
-            static void mergePoints
-            (
-                const polyMesh&,
-                const Map<label>& pointToMaster,
-                polyTopoChange& meshMod
-            );
 };