diff --git a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
index a413e8037e..81bf9400b3 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
@@ -561,6 +561,17 @@ Foam::labelList Foam::polyBoundaryMesh::patchSizes() const
 }
 
 
+Foam::List<Foam::labelRange> Foam::polyBoundaryMesh::patchRanges() const
+{
+    return
+        PtrListOps::get<labelRange>
+        (
+            *this,
+            [](const polyPatch& p) { return p.range(); }
+        );
+}
+
+
 Foam::label Foam::polyBoundaryMesh::start() const
 {
     return mesh_.nInternalFaces();
diff --git a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H
index 8f7e9c86b8..c3d0b0f4b5 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -174,6 +174,9 @@ public:
         //- Return a list of patch sizes
         labelList patchSizes() const;
 
+        //- Return a list of patch ranges
+        List<labelRange> patchRanges() const;
+
         //- The start label of the boundary faces in the polyMesh face list
         //  Same as mesh.nInternalFaces()
         label start() const;
diff --git a/src/dynamicMesh/fvMeshAdder/fvMeshAdder.C b/src/dynamicMesh/fvMeshAdder/fvMeshAdder.C
index cdb2f25b55..db8400b0e4 100644
--- a/src/dynamicMesh/fvMeshAdder/fvMeshAdder.C
+++ b/src/dynamicMesh/fvMeshAdder/fvMeshAdder.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,6 +29,7 @@ License
 #include "fvMesh.H"
 #include "fvMeshAdder.H"
 #include "faceCoupleInfo.H"
+#include "polyTopoChange.H"
 
 /* * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * * */
 
@@ -50,8 +52,8 @@ Foam::labelList Foam::fvMeshAdder::calcPatchMap
 {
     labelList newToOld(newPatch.size(), unmappedValue);
 
-    label newStart = newPatch.start();
-    label newSize = newPatch.size();
+    const label newStart = newPatch.start();
+    const label newSize = newPatch.size();
 
     for (label i = 0; i < oldSize; i++)
     {
@@ -129,4 +131,335 @@ Foam::autoPtr<Foam::mapAddedPolyMesh> Foam::fvMeshAdder::add
 }
 
 
+Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshAdder::add
+(
+    const label myProci,            // index of mesh to modify
+    UPtrList<fvMesh>& fvMeshes,
+    const labelList& oldFaceOwner,  // face owner for myProci mesh
+
+    // Coupling info
+    const labelListList& localBoundaryFace,
+    const labelListList& remoteFaceProc,
+    const labelListList& remoteBoundaryFace,
+
+    labelListList& constructPatchMap,
+    labelListList& constructCellMap,
+    labelListList& constructFaceMap,
+    labelListList& constructPointMap
+)
+{
+    // Do in-place addition. Modifies fvMeshes[myProci]
+
+    UPtrList<polyMesh> meshes(fvMeshes.size());
+    forAll(fvMeshes, proci)
+    {
+        if (fvMeshes.set(proci))
+        {
+            meshes.set(proci, &fvMeshes[proci]);
+        }
+    }
+
+
+    // All matched faces assumed to have vertex0 matched
+    labelListList remoteFaceStart(meshes.size());
+    forAll(localBoundaryFace, proci)
+    {
+        const labelList& procFaces = localBoundaryFace[proci];
+        remoteFaceStart[proci].setSize(procFaces.size(), 0);
+    }
+
+    // Assume all meshes have same global patches
+    labelListList patchMap(meshes.size());
+    labelListList pointZoneMap(meshes.size());
+    labelListList faceZoneMap(meshes.size());
+    labelListList cellZoneMap(meshes.size());
+    forAll(meshes, proci)
+    {
+        if (meshes.set(proci))
+        {
+            const polyMesh& mesh = meshes[proci];
+            patchMap[proci] = identity(mesh.boundaryMesh().nNonProcessor());
+            pointZoneMap[proci] = identity(mesh.pointZones().size());
+            faceZoneMap[proci] = identity(mesh.faceZones().size());
+            cellZoneMap[proci] = identity(mesh.cellZones().size());
+        }
+    }
+
+
+    // Swap myProci to 0th element
+    if (myProci != 0)
+    {
+        polyMesh* pm0 = meshes.get(0);
+        polyMesh* pmi = meshes.get(myProci);
+        meshes.set(0, pmi);
+        meshes.set(myProci, pm0);
+
+        fvMesh* fvm0 = fvMeshes.get(0);
+        fvMesh* fvmi = fvMeshes.get(myProci);
+        fvMeshes.set(0, fvmi);
+        fvMeshes.set(myProci, fvm0);
+
+        //Pout<< "swapped from 0 to " << myProci << endl;
+        //forAll(meshes, meshi)
+        //{
+        //    Pout<< "meshi:" << meshi << endl;
+        //    if (meshes.set(meshi))
+        //    {
+        //        Pout<< "    nCells:" << meshes[meshi].nCells() << endl;
+        //    }
+        //}
+
+
+        // Swap (and renumber) patch face information
+        labelListList& lbf = const_cast<labelListList&>(localBoundaryFace);
+        std::swap(lbf[0], lbf[myProci]);
+
+        labelListList& rfp = const_cast<labelListList&>(remoteFaceProc);
+        std::swap(rfp[0], rfp[myProci]);
+        forAll(rfp, proci)
+        {
+            for (label& proc : rfp[proci])
+            {
+                if (proc == 0) proc = myProci;
+                else if (proc == myProci) proc = 0;
+            }
+        }
+        labelListList& rbf = const_cast<labelListList&>(remoteBoundaryFace);
+        std::swap(rbf[0], rbf[myProci]);
+
+        labelListList& rfs = const_cast<labelListList&>(remoteFaceStart);
+        std::swap(rfs[0], rfs[myProci]);
+
+        // Swap optional renumbering maps
+        std::swap(patchMap[0], patchMap[myProci]);
+        std::swap(pointZoneMap[0], pointZoneMap[myProci]);
+        std::swap(faceZoneMap[0], faceZoneMap[myProci]);
+        std::swap(cellZoneMap[0], cellZoneMap[myProci]);
+    }
+
+    polyTopoChange meshMod(meshes[0].boundaryMesh().size(), true);
+    // Collect statistics for sizing
+    label nCells = 0;
+    label nFaces = 0;
+    label nPoints = 0;
+    forAll(meshes, proci)
+    {
+        if (meshes.set(proci))
+        {
+            const polyMesh& mesh = meshes[proci];
+            nCells += mesh.nCells();
+            nFaces += mesh.nFaces();
+            nPoints += mesh.nPoints();
+        }
+    }
+    meshMod.setCapacity(nPoints, nFaces, nCells);
+
+    // Add all cells in meshes' order
+    polyMeshAdder::add
+    (
+        meshes,
+        patchMap,
+
+        // Information on (one-to-one) boundary faces to stitch
+        localBoundaryFace,
+        remoteFaceProc,
+        remoteBoundaryFace,
+        remoteFaceStart,
+
+        pointZoneMap,
+        faceZoneMap,
+        cellZoneMap,
+
+        meshMod,
+
+        constructCellMap,
+        constructFaceMap,
+        constructPointMap
+    );
+
+    // Replace mesh
+    autoPtr<mapPolyMesh> mapPtr
+    (
+        meshMod.changeMesh
+        (
+            fvMeshes[0],        // note: still swapped to position 0
+            false               // no inflation
+        )
+    );
+
+    // Update fields. Note that this tries to interpolate from the mesh
+    // before adding all the remote meshes so is quite wrong.
+    //fvMeshes[0.updateMesh(mapPtr());
+    // Update polyMesh but not fvMesh to avoid mapping all the fields
+    fvMeshes[0].polyMesh::updateMesh(mapPtr());
+
+    // Now reverseFaceMap contains from order in which face was added
+    // to mesh face (after e.g. upper-triangular ordering)
+
+    // Renumber output of any mesh ordering done by changeMesh
+    for (labelList& cellMap : constructCellMap)
+    {
+        inplaceRenumber(mapPtr().reverseCellMap(), cellMap);
+    }
+    for (labelList& faceMap : constructFaceMap)
+    {
+        inplaceRenumber(mapPtr().reverseFaceMap(), faceMap);
+    }
+    for (labelList& pointMap : constructPointMap)
+    {
+        inplaceRenumber(mapPtr().reversePointMap(), pointMap);
+    }
+
+    // constructPatchMap is patchMap with -1 for the removed
+    // patches
+    forAll(meshes, meshi)
+    {
+        if (meshes.set(meshi))
+        {
+            constructPatchMap[meshi] = patchMap[meshi];
+        }
+    }
+
+    // Map all fields
+    fvMeshAdder::MapVolFields<scalar>
+    (
+        fvMeshes,
+        mapPtr().oldPatchStarts(),
+        mapPtr().oldPatchSizes(),
+        patchMap,
+        constructCellMap,
+        constructFaceMap,
+        constructPointMap
+    );
+    fvMeshAdder::MapVolFields<vector>
+    (
+        fvMeshes,
+        mapPtr().oldPatchStarts(),
+        mapPtr().oldPatchSizes(),
+        patchMap,
+        constructCellMap,
+        constructFaceMap,
+        constructPointMap
+    );
+    fvMeshAdder::MapVolFields<sphericalTensor>
+    (
+        fvMeshes,
+        mapPtr().oldPatchStarts(),
+        mapPtr().oldPatchSizes(),
+        patchMap,
+        constructCellMap,
+        constructFaceMap,
+        constructPointMap
+    );
+    fvMeshAdder::MapVolFields<symmTensor>
+    (
+        fvMeshes,
+        mapPtr().oldPatchStarts(),
+        mapPtr().oldPatchSizes(),
+        patchMap,
+        constructCellMap,
+        constructFaceMap,
+        constructPointMap
+    );
+    fvMeshAdder::MapVolFields<tensor>
+    (
+        fvMeshes,
+        mapPtr().oldPatchStarts(),
+        mapPtr().oldPatchSizes(),
+        patchMap,
+        constructCellMap,
+        constructFaceMap,
+        constructPointMap
+    );
+    fvMeshAdder::MapSurfaceFields<scalar>
+    (
+        fvMeshes,
+        oldFaceOwner,
+        mapPtr().oldPatchStarts(),
+        mapPtr().oldPatchSizes(),
+        patchMap,
+        constructCellMap,
+        constructFaceMap,
+        constructPointMap
+    );
+    fvMeshAdder::MapSurfaceFields<vector>
+    (
+        fvMeshes,
+        oldFaceOwner,
+        mapPtr().oldPatchStarts(),
+        mapPtr().oldPatchSizes(),
+        patchMap,
+        constructCellMap,
+        constructFaceMap,
+        constructPointMap
+    );
+    fvMeshAdder::MapSurfaceFields<sphericalTensor>
+    (
+        fvMeshes,
+        oldFaceOwner,
+        mapPtr().oldPatchStarts(),
+        mapPtr().oldPatchSizes(),
+        patchMap,
+        constructCellMap,
+        constructFaceMap,
+        constructPointMap
+    );
+    fvMeshAdder::MapSurfaceFields<symmTensor>
+    (
+        fvMeshes,
+        oldFaceOwner,
+        mapPtr().oldPatchStarts(),
+        mapPtr().oldPatchSizes(),
+        patchMap,
+        constructCellMap,
+        constructFaceMap,
+        constructPointMap
+    );
+    fvMeshAdder::MapSurfaceFields<tensor>
+    (
+        fvMeshes,
+        oldFaceOwner,
+        mapPtr().oldPatchStarts(),
+        mapPtr().oldPatchSizes(),
+        patchMap,
+        constructCellMap,
+        constructFaceMap,
+        constructPointMap
+    );
+
+
+    // Swap returned data back to processor order
+    if (myProci != 0)
+    {
+        fvMesh* fvm0 = fvMeshes.get(0);
+        fvMesh* fvmi = fvMeshes.get(myProci);
+        fvMeshes.set(0, fvmi);
+        fvMeshes.set(myProci, fvm0);
+
+        // Swap (and renumber) patch face information
+        labelListList& lbf = const_cast<labelListList&>(localBoundaryFace);
+        std::swap(lbf[0], lbf[myProci]);
+        labelListList& rfp = const_cast<labelListList&>(remoteFaceProc);
+        std::swap(rfp[0], rfp[myProci]);
+        forAll(rfp, proci)
+        {
+            for (label& proc : rfp[proci])
+            {
+                if (proc == 0) proc = myProci;
+                else if (proc == myProci) proc = 0;
+            }
+        }
+        labelListList& rbf = const_cast<labelListList&>(remoteBoundaryFace);
+        std::swap(rbf[0], rbf[myProci]);
+
+        std::swap(constructPatchMap[0], constructPatchMap[myProci]);
+        std::swap(constructCellMap[0], constructCellMap[myProci]);
+        std::swap(constructFaceMap[0], constructFaceMap[myProci]);
+        std::swap(constructPointMap[0], constructPointMap[myProci]);
+    }
+
+    return mapPtr;
+}
+
+
 // ************************************************************************* //
diff --git a/src/dynamicMesh/fvMeshAdder/fvMeshAdder.H b/src/dynamicMesh/fvMeshAdder/fvMeshAdder.H
index 0c10219b89..30c459466f 100644
--- a/src/dynamicMesh/fvMeshAdder/fvMeshAdder.H
+++ b/src/dynamicMesh/fvMeshAdder/fvMeshAdder.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -164,6 +165,109 @@ public:
             const fvMesh& mesh,
             const fvMesh& meshToAdd
         );
+
+
+        // Multi-mesh merging. Expects same patches (apart from processor
+        // patches). Adds to element 0 of flds, meshes
+
+            //- In-place add to fvMeshes[myProci]. Stitch boundary faces.
+            static autoPtr<mapPolyMesh> add
+            (
+                const label myProci,    // index of mesh to modify (== mesh_)
+                UPtrList<fvMesh>& fvMeshes,
+                const labelList& oldFaceOwner,  // face owner for myProci mesh
+
+                // Coupling info
+                const labelListList& localBoundaryFace,
+                const labelListList& remoteFaceProc,
+                const labelListList& remoteBoundaryFace,
+
+                labelListList& constructPatchMap,
+                labelListList& constructCellMap,
+                labelListList& constructFaceMap,
+                labelListList& constructPointMap
+            );
+
+            //- Update single dimensionedField.
+            template<class Type>
+            static void MapDimField
+            (
+                UPtrList<DimensionedField<Type, volMesh>>& flds,
+                const labelListList& cellProcAddressing,
+                const bool fullyMapped
+            );
+
+            //- Update single volField.
+            template<class Type>
+            static void MapVolField
+            (
+                UPtrList<GeometricField<Type, fvPatchField, volMesh>>& flds,
+                const labelList& oldPatchStarts0,
+                const labelList& oldPatchSizes0,
+                const labelListList& patchProcAddressing,
+                const labelListList& cellProcAddressing,
+                const labelListList& faceProcAddressing,
+                const bool fullyMapped
+            );
+
+            //- Update single surfaceField.
+            template<class Type>
+            static void MapSurfaceField
+            (
+                UPtrList<GeometricField<Type, fvsPatchField, surfaceMesh>>&,
+                const labelList& oldFaceOwner0,
+                const labelList& oldPatchStarts0,
+                const labelList& oldPatchSizes0,
+                const labelListList& patchProcAddressing,
+                const labelListList& cellProcAddressing,
+                const labelListList& faceProcAddressing,
+                const bool fullyMapped
+            );
+
+
+            //- Map all dimensionedField of Type.
+            //  Optionally override mapping detection
+            //  of unmapped values (e.g. used in fvMeshDistribute since it
+            //  fixes up mapping itself)
+            template<class Type>
+            static void MapDimFields
+            (
+                const UPtrList<fvMesh>& meshes,
+                const labelListList& cellProcAddressing,
+                const bool fullyMapped = false
+            );
+
+            //- Map all volFields of Type.
+            //  Optionally override mapping detection
+            //  of unmapped values (e.g. used in fvMeshDistribute since it
+            //  fixes up mapping itself)
+            template<class Type>
+            static void MapVolFields
+            (
+                const UPtrList<fvMesh>& meshes,
+                const labelList& oldPatchStarts0,
+                const labelList& oldPatchSizes0,
+                const labelListList& patchProcAddressing,
+                const labelListList& cellProcAddressing,
+                const labelListList& faceProcAddressing,
+                const labelListList& pointProcAddressing,
+                const bool fullyMapped = false
+            );
+
+            //- Map all surfaceFields of Type
+            template<class Type>
+            static void MapSurfaceFields
+            (
+                const UPtrList<fvMesh>& meshes,
+                const labelList& oldFaceOwner0,
+                const labelList& oldPatchStarts0,
+                const labelList& oldPatchSizes0,
+                const labelListList& patchProcAddressing,
+                const labelListList& cellProcAddressing,
+                const labelListList& faceProcAddressing,
+                const labelListList& pointProcAddressing,
+                const bool fullyMapped = false
+            );
 };
 
 
diff --git a/src/dynamicMesh/fvMeshAdder/fvMeshAdderTemplates.C b/src/dynamicMesh/fvMeshAdder/fvMeshAdderTemplates.C
index f7ad465343..45db30ae6d 100644
--- a/src/dynamicMesh/fvMeshAdder/fvMeshAdderTemplates.C
+++ b/src/dynamicMesh/fvMeshAdder/fvMeshAdderTemplates.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2019 OpenCFD Ltd.
+    Copyright (C) 2015-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -50,7 +50,7 @@ void Foam::fvMeshAdder::MapVolField
 
     {
         // Store old internal field
-        Field<Type> oldInternalField(fld.primitiveField());
+        const Field<Type> oldInternalField(fld.primitiveField());
 
         // Modify internal field
         Field<Type>& intFld = fld.primitiveFieldRef();
@@ -127,7 +127,7 @@ void Foam::fvMeshAdder::MapVolField
 
         forAll(oldPatchMap, patchi)
         {
-            label newPatchi = oldPatchMap[patchi];
+            const label newPatchi = oldPatchMap[patchi];
 
             if (newPatchi != -1)
             {
@@ -185,7 +185,7 @@ void Foam::fvMeshAdder::MapVolField
         // Add addedMesh patches
         forAll(addedPatchMap, patchi)
         {
-            label newPatchi = addedPatchMap[patchi];
+            const label newPatchi = addedPatchMap[patchi];
 
             if (newPatchi != -1)
             {
@@ -344,7 +344,7 @@ void Foam::fvMeshAdder::MapSurfaceField
 
     // Store old internal field
     {
-        Field<Type> oldField(fld);
+        const Field<Type> oldField(fld);
 
         // Modify internal field
         Field<Type>& intFld = fld.primitiveFieldRef();
@@ -390,7 +390,7 @@ void Foam::fvMeshAdder::MapSurfaceField
 
         forAll(oldPatchMap, patchi)
         {
-            label newPatchi = oldPatchMap[patchi];
+            const label newPatchi = oldPatchMap[patchi];
 
             if (newPatchi != -1)
             {
@@ -439,7 +439,7 @@ void Foam::fvMeshAdder::MapSurfaceField
 
         forAll(oldPatchMap, patchi)
         {
-            label newPatchi = oldPatchMap[patchi];
+            const label newPatchi = oldPatchMap[patchi];
 
             if (newPatchi != -1)
             {
@@ -497,7 +497,7 @@ void Foam::fvMeshAdder::MapSurfaceField
         // Add addedMesh patches
         forAll(addedPatchMap, patchi)
         {
-            label newPatchi = addedPatchMap[patchi];
+            const label newPatchi = addedPatchMap[patchi];
 
             if (newPatchi != -1)
             {
@@ -648,7 +648,7 @@ void Foam::fvMeshAdder::MapDimField
     const fvMesh& mesh = fld.mesh();
 
     // Store old field
-    Field<Type> oldField(fld);
+    const Field<Type> oldField(fld);
 
     fld.setSize(mesh.nCells());
 
@@ -704,4 +704,588 @@ void Foam::fvMeshAdder::MapDimFields
 }
 
 
+//- Multi-mesh mapping
+
+template<class Type>
+void Foam::fvMeshAdder::MapDimField
+(
+    UPtrList<DimensionedField<Type, volMesh>>& flds,
+    const labelListList& cellProcAddressing,
+    const bool fullyMapped
+)
+{
+    // Add fields to fields[0] after adding the meshes to meshes[0].
+    // Mesh[0] is the sum of all meshes. Fields are not yet mapped.
+
+    if (flds.size() == 0 || !flds.set(0))
+    {
+        FatalErrorInFunction << "Not valid field at element 0"
+            << " in field list of size " << flds.size() << exit(FatalError);
+    }
+
+
+    // Internal field
+    // ~~~~~~~~~~~~~~
+
+    {
+        // Store old internal field
+        const Field<Type> oldInternalField(flds[0].primitiveField());
+
+        // Modify internal field
+        Field<Type>& intFld = flds[0].primitiveFieldRef();
+
+        // Set to new mesh size
+        intFld.setSize(flds[0].mesh().nCells());
+        // Add fld0
+        intFld.rmap(oldInternalField, cellProcAddressing[0]);
+
+        for (label meshi = 1; meshi < flds.size(); meshi++)
+        {
+            if (flds.set(meshi))
+            {
+                const Field<Type>& addFld = flds[meshi].primitiveFieldRef();
+                intFld.rmap(addFld, cellProcAddressing[meshi]);
+            }
+        }
+    }
+}
+
+
+template<class Type>
+void Foam::fvMeshAdder::MapVolField
+(
+    UPtrList<GeometricField<Type, fvPatchField, volMesh>>& flds,
+    const labelList& oldPatchStarts0,
+    const labelList& oldPatchSizes0,
+    const labelListList& patchProcAddressing,
+    const labelListList& cellProcAddressing,
+    const labelListList& faceProcAddressing,
+    const bool fullyMapped
+)
+{
+    // Add fields to fields[0] after adding the meshes to meshes[0].
+    // Mesh[0] is the sum of all meshes. Fields are not yet mapped.
+
+    if (flds.size() == 0 || !flds.set(0))
+    {
+        FatalErrorInFunction << "Not valid field at element 0"
+            << " in field list of size " << flds.size() << exit(FatalError);
+    }
+
+
+    // Internal field
+    // ~~~~~~~~~~~~~~
+
+    {
+        // Store old internal field
+        const Field<Type> oldInternalField(flds[0].primitiveField());
+
+        // Modify internal field
+        Field<Type>& intFld = flds[0].primitiveFieldRef();
+
+        // Set to new mesh size
+        intFld.setSize(flds[0].mesh().nCells());
+        // Add fld0
+        intFld.rmap(oldInternalField, cellProcAddressing[0]);
+
+        for (label meshi = 1; meshi < flds.size(); meshi++)
+        {
+            if (flds.set(meshi))
+            {
+                const Field<Type>& addFld = flds[meshi].primitiveFieldRef();
+                intFld.rmap(addFld, cellProcAddressing[meshi]);
+            }
+        }
+    }
+
+
+    // Patch fields from old mesh
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    auto& bfld0 = flds[0].boundaryFieldRef();
+    //Pout<< "   Adding patchfields for field " << flds[0].name() << endl;
+    forAll(bfld0, patchi)
+    {
+        //Pout<< "    patch:" << patchi
+        //    << " patch:" << flds[0].mesh().boundaryMesh()[patchi].name()
+        //    << " size:" << flds[0].mesh().boundaryMesh()[patchi].size()
+        //    << endl;
+        //Pout<< "    patchField:" << bfld0[patchi].size()
+        //    << endl;
+
+        labelList newToOld
+        (
+            calcPatchMap
+            (
+                oldPatchStarts0[patchi],
+                oldPatchSizes0[patchi],
+                faceProcAddressing[0],
+                bfld0[patchi].patch().patch(),
+                -1              // unmapped value
+            )
+        );
+
+        directFvPatchFieldMapper patchMapper(newToOld);
+
+        // Override mapping (for use in e.g. fvMeshDistribute where
+        // it sorts mapping out itself)
+        if (fullyMapped)
+        {
+            patchMapper.hasUnmapped() = false;
+        }
+
+        bfld0[patchi].autoMap(patchMapper);
+    }
+
+    for (label meshi = 1; meshi < flds.size(); meshi++)
+    {
+        if (flds.set(meshi))
+        {
+            const auto& bfld = flds[meshi].boundaryFieldRef();
+
+            const labelList& patchMap = patchProcAddressing[meshi];
+
+            forAll(patchMap, oldPatchi)
+            {
+                const auto& fvp = bfld[oldPatchi].patch();
+                const label newPatchi = patchMap[oldPatchi];
+
+                //Pout<< "    oldPatch:" << oldPatchi
+                //    << " newPatch:" << newPatchi << endl;
+
+                if (newPatchi >= 0 && newPatchi < bfld0.size())
+                {
+                    const auto& fvp0 = bfld0[newPatchi].patch();
+                    labelList addedToNew(bfld[oldPatchi].size(), -1);
+                    forAll(addedToNew, i)
+                    {
+                        const label newFacei =
+                            faceProcAddressing[meshi][fvp.start()+i];
+                        const label patchFacei = newFacei-fvp0.start();
+                        if
+                        (
+                            patchFacei >= 0
+                         && patchFacei < fvp0.size()
+                        )
+                        {
+                            addedToNew[i] = patchFacei;
+                        }
+                    }
+
+                    bfld0[newPatchi].rmap(bfld[oldPatchi], addedToNew);
+                }
+                else
+                {
+                    WarningInFunction << "Ignoring old patch "
+                        << bfld[oldPatchi].patch().name() << " on field "
+                        << flds[meshi].name() << endl;  //exit(FatalError);
+                }
+            }
+        }
+    }
+}
+
+
+template<class Type>
+void Foam::fvMeshAdder::MapSurfaceField
+(
+    UPtrList<GeometricField<Type, fvsPatchField, surfaceMesh>>& flds,
+    const labelList& oldFaceOwner0,
+    const labelList& oldPatchStarts0,
+    const labelList& oldPatchSizes0,
+    const labelListList& patchProcAddressing,
+    const labelListList& cellProcAddressing,
+    const labelListList& faceProcAddressing,
+    const bool fullyMapped
+)
+{
+    // Add fields to fields[0] after adding the meshes to meshes[0].
+    // Mesh[0] is the sum of all meshes. Fields are not yet mapped.
+
+    if (flds.size() == 0 || !flds.set(0))
+    {
+        FatalErrorInFunction << "Not valid field at element 0"
+            << " in field list of size " << flds.size() << exit(FatalError);
+    }
+
+    const fvMesh& mesh0 = flds[0].mesh();
+
+
+    // Internal field
+    // ~~~~~~~~~~~~~~
+
+    {
+        // Store old internal field
+        const Field<Type> oldInternalField(flds[0].primitiveField());
+
+        // Modify internal field
+        Field<Type>& intFld = flds[0].primitiveFieldRef();
+
+        // Set to new mesh size
+        intFld.setSize(mesh0.nInternalFaces());
+
+        // Map
+        forAll(flds, meshi)
+        {
+            if (flds.set(meshi))
+            {
+                const labelList& faceMap = faceProcAddressing[meshi];
+                const auto& fld = flds[meshi];
+
+                // Map internal field
+                if (meshi == 0)
+                {
+                    intFld.rmap(oldInternalField, faceMap);
+                }
+                else
+                {
+                    intFld.rmap(fld.primitiveField(), faceMap);
+                }
+
+                // Map faces that were boundary faces but are not anymore.
+                // There will be two meshes that provide the same face. Use
+                // owner one.
+                const auto& bfld = flds[meshi].boundaryField();
+
+                forAll(bfld, oldPatchi)
+                {
+                    const fvsPatchField<Type>& pf = bfld[oldPatchi];
+                    //Pout<< "patch:" << pf.patch().name() << endl;
+                    forAll(pf, patchFacei)
+                    {
+                        // Get new face, mapped face owner
+                        label newFacei;
+                        label newOwn;
+                        if (meshi == 0)
+                        {
+                            // Do not access mesh information since in-place
+                            // modified
+                            const label oldFacei =
+                                oldPatchStarts0[oldPatchi]+patchFacei;
+                            newFacei = faceProcAddressing[meshi][oldFacei];
+                            const label oldOwn = oldFaceOwner0[oldFacei];
+                            newOwn = cellProcAddressing[meshi][oldOwn];
+
+                            //Pout<< "MESH0: pfi:" << patchFacei
+                            //    << " old face:" << oldFacei
+                            //    << " new face:" << newFacei
+                            //    << " at:" << mesh0.faceCentres()[newFacei]
+                            //    << " oldOwn:" << oldOwn
+                            //    << " newOwn:" << newOwn << endl;
+                        }
+                        else
+                        {
+                            const label oldFacei =
+                                pf.patch().start()+patchFacei;
+                            newFacei = faceProcAddressing[meshi][oldFacei];
+                            const label oldOwn =
+                                fld.mesh().faceOwner()[oldFacei];
+                            newOwn = cellProcAddressing[meshi][oldOwn];
+
+                            //Pout<< "MESH:" << meshi << " pfi:" << patchFacei
+                            //    << " old face:" << oldFacei
+                            //    << " new face:" << newFacei
+                            //    << " at:" << mesh0.faceCentres()[newFacei]
+                            //    << " oldOwn:" << oldOwn
+                            //    << " newOwn:" << newOwn << endl;
+                        }
+
+                        if
+                        (
+                            newFacei >= 0
+                         && newFacei < mesh0.nInternalFaces()
+                         && (newOwn == mesh0.faceOwner()[newFacei])
+                        )
+                        {
+                            intFld[newFacei] = pf[patchFacei];
+                        }
+                        //else
+                        //{
+                        //    Pout<< "    ignoring pfi:" << patchFacei
+                        //        << " value:" << pf[patchFacei]
+                        //        << " since newFacei:" << newFacei
+                        //        << " since newOwn:" << newOwn
+                        //        << " own:" << mesh0.faceOwner()[newFacei]
+                        //        << endl;
+                        //}
+                    }
+                }
+            }
+        }
+    }
+
+
+    // Patch fields from old mesh
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    auto& bfld0 = flds[0].boundaryFieldRef();
+    forAll(bfld0, patchi)
+    {
+        labelList newToOld
+        (
+            calcPatchMap
+            (
+                oldPatchStarts0[patchi],
+                oldPatchSizes0[patchi],
+                faceProcAddressing[0],
+                bfld0[patchi].patch().patch(),
+                -1              // unmapped value
+            )
+        );
+
+        directFvPatchFieldMapper patchMapper(newToOld);
+
+        // Override mapping (for use in e.g. fvMeshDistribute where
+        // it sorts mapping out itself)
+        if (fullyMapped)
+        {
+            patchMapper.hasUnmapped() = false;
+        }
+
+        bfld0[patchi].autoMap(patchMapper);
+    }
+
+    for (label meshi = 1; meshi < flds.size(); meshi++)
+    {
+        if (flds.set(meshi))
+        {
+            const auto& bfld = flds[meshi].boundaryFieldRef();
+
+            const labelList& patchMap = patchProcAddressing[meshi];
+
+            forAll(patchMap, oldPatchi)
+            {
+                const auto& fvp = bfld[oldPatchi].patch();
+                const label newPatchi = patchMap[oldPatchi];
+                if (newPatchi >= 0 && newPatchi < bfld0.size())
+                {
+                    const auto& fvp0 = bfld0[newPatchi].patch();
+                    labelList addedToNew(bfld[oldPatchi].size(), -1);
+                    forAll(addedToNew, i)
+                    {
+                        const label newFacei =
+                            faceProcAddressing[meshi][fvp.start()+i];
+                        const label patchFacei = newFacei-fvp0.start();
+                        if
+                        (
+                            patchFacei >= 0
+                         && patchFacei < fvp0.size()
+                        )
+                        {
+                            addedToNew[i] = patchFacei;
+                        }
+                    }
+
+                    bfld0[newPatchi].rmap(bfld[oldPatchi], addedToNew);
+                }
+                else
+                {
+                    WarningInFunction << "Ignoring old patch "
+                        << bfld[oldPatchi].patch().name() << " on field "
+                        << flds[meshi].name() << endl;  //exit(FatalError);
+                }
+            }
+        }
+    }
+}
+
+
+template<class Type>
+void Foam::fvMeshAdder::MapVolFields
+(
+    const UPtrList<fvMesh>& meshes,
+    const labelList& oldPatchStarts0,
+    const labelList& oldPatchSizes0,
+    const labelListList& patchProcAddressing,
+    const labelListList& cellProcAddressing,
+    const labelListList& faceProcAddressing,
+    const labelListList& pointProcAddressing,
+    const bool fullyMapped
+)
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> fldType;
+
+    if (meshes.size() == 0 || !meshes.set(0))
+    {
+        FatalErrorInFunction << "Not valid field at element 0"
+            << " in mesh list of size " << meshes.size() << exit(FatalError);
+    }
+    const fvMesh& mesh0 = meshes[0];
+
+    HashTable<const fldType*> fields
+    (
+        mesh0.objectRegistry::lookupClass<fldType>()
+    );
+
+
+    // It is necessary to enforce that all old-time fields are stored
+    // before the mapping is performed.  Otherwise, if the
+    // old-time-level field is mapped before the field itself, sizes
+    // will not match.
+
+    for (const auto& fld : fields)
+    {
+        DebugPout
+            << "MapVolFields : Storing old time for " << fld->name()
+            << endl;
+
+        const_cast<fldType&>(*fld).storeOldTimes();
+    }
+
+
+    for (const auto& fld : fields)
+    {
+        const word& name0 = fld->name();
+
+        DebugPout
+            << "MapVolFields : mapping " << name0 << endl;
+
+        UPtrList<fldType> meshToField(meshes.size());
+        forAll(meshes, meshi)
+        {
+            if (meshes.set(meshi))
+            {
+                auto& meshFld = meshes[meshi].
+                    objectRegistry::lookupObjectRef<fldType>(name0);
+                meshToField.set(meshi, &meshFld);
+            }
+        }
+
+        MapVolField
+        (
+            meshToField,
+            oldPatchStarts0,
+            oldPatchSizes0,
+            patchProcAddressing,
+            cellProcAddressing,
+            faceProcAddressing,
+            fullyMapped
+        );
+    }
+}
+
+
+template<class Type>
+void Foam::fvMeshAdder::MapDimFields
+(
+    const UPtrList<fvMesh>& meshes,
+    const labelListList& cellProcAddressing,
+    const bool fullyMapped
+)
+{
+    typedef DimensionedField<Type, volMesh> fldType;
+
+    if (meshes.size() == 0 || !meshes.set(0))
+    {
+        FatalErrorInFunction << "Not valid field at element 0"
+            << " in mesh list of size " << meshes.size() << exit(FatalError);
+    }
+    const fvMesh& mesh0 = meshes[0];
+
+    HashTable<const fldType*> fields
+    (
+        mesh0.objectRegistry::lookupClass<fldType>()
+    );
+
+
+    for (const auto& fld : fields)
+    {
+        const word& name0 = fld->name();
+
+        DebugPout
+            << "MapDimFields : mapping " << name0 << endl;
+
+        UPtrList<fldType> meshToField(meshes.size());
+        forAll(meshes, meshi)
+        {
+            if (meshes.set(meshi))
+            {
+                auto& meshFld = meshes[meshi].
+                    objectRegistry::lookupObjectRef<fldType>(name0);
+                meshToField.set(meshi, &meshFld);
+            }
+        }
+
+        MapDimField(meshToField, cellProcAddressing, fullyMapped);
+    }
+}
+
+
+template<class Type>
+void Foam::fvMeshAdder::MapSurfaceFields
+(
+    const UPtrList<fvMesh>& meshes,
+    const labelList& oldFaceOwner0,
+    const labelList& oldPatchStarts0,
+    const labelList& oldPatchSizes0,
+    const labelListList& patchProcAddressing,
+    const labelListList& cellProcAddressing,
+    const labelListList& faceProcAddressing,
+    const labelListList& pointProcAddressing,
+    const bool fullyMapped
+)
+{
+    typedef GeometricField<Type, fvsPatchField, surfaceMesh> fldType;
+
+    if (meshes.size() == 0 || !meshes.set(0))
+    {
+        FatalErrorInFunction << "Not valid field at element 0"
+            << " in mesh list of size " << meshes.size() << exit(FatalError);
+    }
+    const auto& mesh0 = meshes[0];
+
+    HashTable<const fldType*> fields
+    (
+        mesh0.objectRegistry::lookupClass<fldType>()
+    );
+
+
+    // It is necessary to enforce that all old-time fields are stored
+    // before the mapping is performed.  Otherwise, if the
+    // old-time-level field is mapped before the field itself, sizes
+    // will not match.
+
+    for (const auto& fld : fields)
+    {
+        DebugPout
+            << "MapSurfaceFields : Storing old time for " << fld->name()
+            << endl;
+
+        const_cast<fldType&>(*fld).storeOldTimes();
+    }
+
+
+    for (const auto& fld : fields)
+    {
+        const word& name0 = fld->name();
+
+        DebugPout
+            << "MapSurfaceFields : Mapping " << fld->name() << endl;
+
+        UPtrList<fldType> meshToField(meshes.size());
+        forAll(meshes, meshi)
+        {
+            if (meshes.set(meshi))
+            {
+                auto& meshFld = meshes[meshi].
+                    objectRegistry::lookupObjectRef<fldType>(name0);
+                meshToField.set(meshi, &meshFld);
+            }
+        }
+
+        MapSurfaceField
+        (
+            meshToField,
+            oldFaceOwner0,
+            oldPatchStarts0,
+            oldPatchSizes0,
+            patchProcAddressing,
+            cellProcAddressing,
+            faceProcAddressing,
+            fullyMapped
+        );
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/src/dynamicMesh/polyMeshAdder/polyMeshAdder.C b/src/dynamicMesh/polyMeshAdder/polyMeshAdder.C
index 5fc8eff9a9..0dfaa92d41 100644
--- a/src/dynamicMesh/polyMeshAdder/polyMeshAdder.C
+++ b/src/dynamicMesh/polyMeshAdder/polyMeshAdder.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019,2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -38,6 +38,7 @@ License
 #include "polyModifyFace.H"
 #include "polyRemovePoint.H"
 #include "polyTopoChange.H"
+#include "globalIndex.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -2121,8 +2122,7 @@ Foam::Map<Foam::label> Foam::polyMeshAdder::findSharedPoints
 
                             //Pout<< "Merging point " << pointi
                             //    << " at " << mesh.points()[pointi]
-                            //    << " into master point "
-                            //    << masterPointi
+                            //    << " into master point " << masterPointi
                             //    << " at " << mesh.points()[masterPointi]
                             //    << endl;
 
@@ -2302,4 +2302,747 @@ void Foam::polyMeshAdder::mergePoints
 }
 
 
+Foam::label Foam::polyMeshAdder::procPatchIndex
+(
+    const polyBoundaryMesh& pbm,
+    const label nbrProci,
+    const label n
+)
+{
+    // Find n'th processor patch going to nbrProci. Usually n=0 but in some
+    // cases (e.g. processorCyclic, implicit cht) there can be more than one
+    // processor patch between two processors.
+
+    label index = n;
+
+    for (label patchi = pbm.nNonProcessor(); patchi < pbm.size(); patchi++)
+    {
+        const processorPolyPatch& pp =
+            refCast<const processorPolyPatch>(pbm[patchi]);
+        if (pp.neighbProcNo() == nbrProci)
+        {
+            if (index == 0)
+            {
+                return patchi;
+            }
+            else
+            {
+                --index;
+            }
+        }
+    }
+
+    FatalErrorInFunction << "no patch found to processor " << nbrProci
+        << ". Current patches:" << pbm.names() << exit(FatalError);
+    return -1;
+}
+
+
+Foam::label Foam::polyMeshAdder::procPatchPairs
+(
+    const UPtrList<polyMesh>& meshes,
+    List<DynamicList<label>>& localPatch,
+    List<DynamicList<label>>& remoteProc,
+    List<DynamicList<label>>& remotePatch
+    //List<labelListList>& remotePatchFace,
+    //List<labelListList>& remotePatchFaceStart
+)
+{
+    // Determine pairs of processor patches:
+    // - remote processor
+    // - patch on remote processor
+    // - face on remote patch
+    // - starting index on remote face
+    localPatch.setSize(meshes.size());
+    remoteProc.setSize(meshes.size());
+    remotePatch.setSize(meshes.size());
+    //remotePatchFace.setSize(meshes.size());
+    //remotePatchFaceStart.setSize(meshes.size());
+
+
+    // Check that all processors have the same globalPatches
+    {
+        const polyBoundaryMesh& pbm = meshes[0].boundaryMesh();
+        const wordList names0(SubList<word>(pbm.names(), pbm.nNonProcessor()));
+        for (label proci = 1; proci < meshes.size(); proci++)
+        {
+            const polyBoundaryMesh& pbm = meshes[proci].boundaryMesh();
+            const wordList names(pbm.names());
+
+            if (SubList<word>(names, pbm.nNonProcessor()) != names0)
+            {
+                FatalErrorInFunction
+                    << "Patch names should be identical on all processors."
+                    << " Processor 0 has " << names0
+                    << ". Processor " << proci
+                    << " has " << names
+                    << exit(FatalError);
+            }
+        }
+    }
+
+
+    // Work array
+    labelList nNeighbourProcs(meshes.size());
+
+    forAll(meshes, proci)
+    {
+        const polyBoundaryMesh& pbm = meshes[proci].boundaryMesh();
+
+        // Running count of number of patches communicating with same nbr
+        // (usually 0)
+        nNeighbourProcs = 0;
+
+        for (label patchi = pbm.nNonProcessor(); patchi < pbm.size(); patchi++)
+        {
+            const processorPolyPatch& pp =
+                refCast<const processorPolyPatch>(pbm[patchi]);
+            if (pp.owner())
+            {
+                const label nbrProci = pp.neighbProcNo();
+                const label nbrPatchi = procPatchIndex
+                (
+                    meshes[nbrProci].boundaryMesh(),
+                    proci,
+                    nNeighbourProcs[nbrProci]
+                );
+
+                const auto& nbrPbm = meshes[nbrProci].boundaryMesh();
+                const auto& nbrPp = nbrPbm[nbrPatchi];
+                if (pp.size() != nbrPp.size())
+                {
+                    FatalErrorInFunction
+                        << "at proc:" << proci
+                        << " processor patch "
+                        << pp.name() << " is not same size " << pp.size()
+                        << " as coupled patch " << nbrPp.name()
+                        << " on proc:" << nbrProci
+                        << " size:" << nbrPp.size()
+                        << exit(FatalError);
+                }
+
+                localPatch[proci].append(patchi);
+                remoteProc[proci].append(nbrProci);
+                remotePatch[proci].append(nbrPatchi);
+
+                localPatch[nbrProci].append(nbrPatchi);
+                remoteProc[nbrProci].append(proci);
+                remotePatch[nbrProci].append(patchi);
+
+                nNeighbourProcs[nbrProci]++;
+            }
+        }
+    }
+
+    //// Fill in the patch face ordering. Assume correct ordering.
+    //forAll(meshes, proci)
+    //{
+    //    const polyBoundaryMesh& pbm = meshes[proci].boundaryMesh();
+    //
+    //    const DynamicList<label>& localPatches = localPatch[proci];
+    //
+    //    labelListList& rpf = remotePatchFace[proci];
+    //    rpf.setSize(localPatches.size());
+    //
+    //    labelListList& rps = remotePatchFaceStart[proci];
+    //    rps.setSize(localPatches.size());
+    //
+    //    forAll(localPatches, i)
+    //    {
+    //        const auto& pp = pbm[localPatches[i]];
+    //        rpf[i] = identity(pp.size());
+    //        rps[i].setSize(pp.size(), 0);
+    //    }
+    //}
+
+    return meshes[0].boundaryMesh().nNonProcessor();
+}
+
+
+void Foam::polyMeshAdder::patchFacePairs
+(
+    const UPtrList<polyMesh>& meshes,
+    const List<DynamicList<label>>& localPatch,
+    const List<DynamicList<label>>& remoteMesh,
+    const List<DynamicList<label>>& remotePatch,
+    labelListList& localBoundaryFace,
+    labelListList& remoteFaceMesh,
+    labelListList& remoteBoundaryFace
+)
+{
+    // Calculates pairs of matching (boundary) faces. Returns
+    //  localBoundaryFace[meshi]  : index of local face (in boundary face
+    //                              indexing!)
+    //  remoteMesh[meshi]         : index of remote mesh
+    //  remoteBoundaryFace[meshi] : index of remote face (in boundary face
+    //                              indexing!)
+    localBoundaryFace.setSize(meshes.size());
+    remoteFaceMesh.setSize(meshes.size());
+    remoteBoundaryFace.setSize(meshes.size());
+
+    forAll(meshes, meshi)
+    {
+        const auto& mesh = meshes[meshi];
+        const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+        const DynamicList<label>& procPatches = localPatch[meshi];
+        const DynamicList<label>& procNbrs = remoteMesh[meshi];
+        const DynamicList<label>& procNbrPatches = remotePatch[meshi];
+
+
+        // Count number of processor faces
+        label nFaces = 0;
+        for (const label patchi : procPatches)
+        {
+            nFaces += pbm[patchi].size();
+        }
+
+        labelList& procFaces = localBoundaryFace[meshi];
+        labelList& remoteMeshes = remoteFaceMesh[meshi];
+        labelList& remoteFaces = remoteBoundaryFace[meshi];
+
+        procFaces.setSize(nFaces);
+        remoteMeshes.setSize(nFaces);
+        remoteFaces.setSize(nFaces);
+
+        nFaces = 0;
+
+        forAll(procPatches, i)
+        {
+            const label patchi = procPatches[i];
+            const label nbrMeshi = procNbrs[i];
+            const label nbrPatchi = procNbrPatches[i];
+
+            const auto& pp = pbm[patchi];
+            const label offset = pp.start()-mesh.nInternalFaces();
+
+            const auto& nbrMesh = meshes[nbrMeshi];
+            const auto& nbrPp = nbrMesh.boundaryMesh()[nbrPatchi];
+            const label nbrOffset = nbrPp.start()-nbrMesh.nInternalFaces();
+
+            forAll(pp, patchFacei)
+            {
+                procFaces[nFaces] = offset+patchFacei;
+                remoteMeshes[nFaces] = nbrMeshi;
+                remoteFaces[nFaces] = nbrOffset+patchFacei;
+                nFaces++;
+            }
+        }
+    }
+}
+
+
+void Foam::polyMeshAdder::compactPoints
+(
+    const UPtrList<polyMesh>& meshes,
+    const labelListList& localBoundaryFace,
+    const labelListList& remoteFaceMesh,
+    const labelListList& remoteBoundaryFace,
+    const labelListList& remoteFaceStart,
+    const globalIndex& globalPoints,
+
+    labelListList& pointProcAddressing,
+    labelListList& localPoints
+)
+{
+    // Determines 'master' points for shared points. Returns
+    //  pointProcAddressing[meshi] : global point for every mesh point
+    //  localPoints[meshi]         : indices of local points that are unique
+
+    // Choose minimum point for shared ones. Iterate minimising global points
+    // until nothing changes
+    while (true)
+    {
+        label nChanged = 0;
+        forAll(meshes, meshi)
+        {
+            if (meshes.set(meshi))
+            {
+                const polyMesh& mesh = meshes[meshi];
+                const faceList& faces = mesh.faces();
+
+                labelList& pAddressing = pointProcAddressing[meshi];
+
+                const labelList& localBFaces = localBoundaryFace[meshi];
+                const labelList& procNbrs = remoteFaceMesh[meshi];
+                const labelList& procNbrBFaces = remoteBoundaryFace[meshi];
+                const labelList& procNbrIndex = remoteFaceStart[meshi];
+
+                forAll(localBFaces, i)
+                {
+                    const label bFacei = localBFaces[i];
+                    const label nbrMeshi = procNbrs[i];
+                    const label nbrBFacei = procNbrBFaces[i];
+
+                    // Local mesh face
+                    const label facei = mesh.nInternalFaces()+bFacei;
+                    const face& f = faces[facei];
+
+                    // Matched (remote) face
+                    const auto& nbrMesh = meshes[nbrMeshi];
+                    const label nbrFacei = nbrMesh.nInternalFaces()+nbrBFacei;
+                    const face& nbrF = nbrMesh.faces()[nbrFacei];
+
+
+                    labelList& nbrAddressing = pointProcAddressing[nbrMeshi];
+
+                    // Starting index is indexed wrt matched, not original
+                    // face. Check!
+                    label nbrFp = procNbrIndex[i];
+
+                    forAll(f, fp)
+                    {
+                        label& ppPointi = pAddressing[f[fp]];
+                        label& nbrPointi = nbrAddressing[nbrF[nbrFp]];
+
+                        //const point& pt = mesh.points()[f[fp]];
+                        //const point& nbrPt = nbrMesh.points()[nbrF[nbrFp]];
+                        //if (mag(pt-nbrPt) > SMALL)
+                        //{
+                        //    FatalErrorInFunction
+                        //        << "Merging differing points :"
+                        //        << " pt:" << pt
+                        //        << " nbrPt:" << nbrPt
+                        //        << exit(FatalError);
+                        //}
+
+                        if (ppPointi < nbrPointi)
+                        {
+                            //Pout<< "on proc:" << nbrMeshi
+                            //    << " point:" << nbrF[nbrFp]
+                            //    << " changing from:" << nbrPointi
+                            //    << " to:" << ppPointi << endl;
+                            nbrPointi = ppPointi;
+                            nChanged++;
+                        }
+                        else if (nbrPointi < ppPointi)
+                        {
+                            //Pout<< "on proc:" << meshi
+                            //    << " point:" << f[fp]
+                            //    << " changing from:" << ppPointi
+                            //    << " to:" << nbrPointi << endl;
+                            ppPointi = nbrPointi;
+                            nChanged++;
+                        }
+
+                        nbrFp = nbrF.rcIndex(nbrFp);
+                    }
+                }
+            }
+        }
+
+        if (nChanged == 0)
+        {
+            break;
+        }
+    }
+
+
+
+    // Compact out unused points
+    localPoints.setSize(meshes.size());
+    labelList globalToCompact(globalPoints.size(), -1);
+    label nGlobal = 0;
+    forAll(meshes, meshi)
+    {
+        if (meshes.set(meshi))
+        {
+            const polyMesh& mesh = meshes[meshi];
+
+            labelList& compactPoints = localPoints[meshi];
+            compactPoints.setSize(mesh.nPoints());
+            label nCompact = 0;
+
+            const labelList& pAddressing = pointProcAddressing[meshi];
+            forAll(pAddressing, pointi)
+            {
+                //Pout<< "proc:" << meshi
+                //    << " localpoint:" << pointi
+                //    << " globalpoint:"
+                //    << globalPoints.toGlobal(meshi, pointi)
+                //    << " merged:" << pAddressing[pointi]
+                //    << endl;
+
+                const label globali = globalPoints.toGlobal(meshi, pointi);
+
+                if (pAddressing[pointi] == globali)
+                {
+                    // Unchanged address
+                    globalToCompact[globali] = nGlobal++;
+                    compactPoints[nCompact++] = pointi;
+                }
+            }
+            compactPoints.setSize(nCompact);
+        }
+    }
+
+    forAll(meshes, meshi)
+    {
+        labelList& pAddressing = pointProcAddressing[meshi];
+        pAddressing = UIndirectList<label>(globalToCompact, pAddressing);
+    }
+}
+
+
+void Foam::polyMeshAdder::add
+(
+    const UPtrList<polyMesh>& meshes,
+    const UList<labelList>& patchMap,
+
+    const labelListList& localBoundaryFace,
+    const labelListList& remoteFaceMesh,
+    const labelListList& remoteBoundaryFace,
+    const labelListList& remoteFaceStart,
+
+    const UList<labelList>& pointZoneMap,
+    const UList<labelList>& faceZoneMap,
+    const UList<labelList>& cellZoneMap,
+
+    polyTopoChange& meshMod,
+    labelListList& cellProcAddressing,
+    labelListList& faceProcAddressing,
+    labelListList& pointProcAddressing
+    //labelListList& boundaryProcAddressing
+)
+{
+    cellProcAddressing.setSize(meshes.size());
+    faceProcAddressing.setSize(meshes.size());
+    pointProcAddressing.setSize(meshes.size());
+    //boundaryProcAddressing.setSize(meshes.size());
+    forAll(cellProcAddressing, meshi)
+    {
+        cellProcAddressing[meshi].clear();
+        faceProcAddressing[meshi].clear();
+        pointProcAddressing[meshi].clear();
+    }
+
+    // Start off with points allocated in processor order.
+    labelList offsets(meshes.size()+1);
+    offsets[0] = 0;
+    forAll(meshes, meshi)
+    {
+        label nPoints = 0;
+        if (meshes.set(meshi))
+        {
+            const polyMesh& mesh = meshes[meshi];
+            nPoints = mesh.nPoints();
+        }
+        pointProcAddressing[meshi] = identity(nPoints, offsets[meshi]);
+        offsets[meshi+1] = offsets[meshi]+nPoints;
+    }
+    const globalIndex globalPoints(std::move(offsets));
+
+    labelListList uniquePoints(meshes.size());
+    compactPoints
+    (
+        meshes,
+
+        localBoundaryFace,
+        remoteFaceMesh,
+        remoteBoundaryFace,
+        remoteFaceStart,
+
+        globalPoints,
+
+        pointProcAddressing,    // per proc, per point the compact,global point
+        uniquePoints            // per proc indices of local preserved points
+    );
+
+
+    // Add cells
+    // ~~~~~~~~~
+    // These are used by face adding later on
+
+    forAll(meshes, meshi)
+    {
+        if (meshes.set(meshi))
+        {
+            const polyMesh& mesh = meshes[meshi];
+            const cellZoneMesh& cellZones = mesh.cellZones();
+
+            // Precalc offset zones
+            labelList newZoneID(mesh.nCells(), -1);
+
+            forAll(cellZones, zonei)
+            {
+                const labelList& cellLabels = cellZones[zonei];
+
+                for (const label celli : cellLabels)
+                {
+                    if (newZoneID[celli] != -1)
+                    {
+                        WarningInFunction
+                            << "Cell:" << celli
+                            << " centre:" << mesh.cellCentres()[celli]
+                            << " is in two zones:"
+                            << cellZones[newZoneID[celli]].name()
+                            << " and " << cellZones[zonei].name() << endl
+                            << "    This is not supported."
+                            << " Continuing with first zone only." << endl;
+                    }
+                    else
+                    {
+                        newZoneID[celli] = cellZoneMap[meshi][zonei];
+                    }
+                }
+            }
+
+            // Add cells in mesh order
+            cellProcAddressing[meshi].setSize(mesh.nCells());
+            for (label celli = 0; celli < mesh.nCells(); celli++)
+            {
+                // Add cell from cell
+                cellProcAddressing[meshi][celli] = meshMod.addCell
+                (
+                    -1,
+                    -1,
+                    -1,
+                    celli,
+                    newZoneID[celli]
+                );
+            }
+        }
+    }
+
+
+    // Add points
+    // ~~~~~~~~~~
+    // These are used by face adding later on
+
+    forAll(meshes, meshi)
+    {
+        if (meshes.set(meshi))
+        {
+            const polyMesh& mesh = meshes[meshi];
+            const pointField& points = mesh.points();
+            const pointZoneMesh& pointZones = mesh.pointZones();
+
+            // Precalc offset zones
+            labelList newZoneID(points.size(), -1);
+
+            forAll(pointZones, zonei)
+            {
+                const labelList& pointLabels = pointZones[zonei];
+
+                for (const label pointi : pointLabels)
+                {
+                    newZoneID[pointi] = pointZoneMap[meshi][zonei];
+                }
+            }
+
+            // Add points in mesh order
+            const labelList& myUniquePoints = uniquePoints[meshi];
+            for (const label pointi : myUniquePoints)
+            {
+                // Rewrite the addressing (should already be compact) just in
+                // case the polyTopoChange does some special ordering
+                pointProcAddressing[meshi][pointi] = meshMod.addPoint
+                (
+                    points[pointi],
+                    pointi,
+                    newZoneID[pointi],
+                    true
+                );
+            }
+        }
+    }
+
+
+
+    // Add faces
+    // ~~~~~~~~~
+
+    forAll(meshes, meshi)
+    {
+        if (meshes.set(meshi))
+        {
+            const polyMesh& mesh = meshes[meshi];
+            faceProcAddressing[meshi].setSize(mesh.nFaces());
+            faceProcAddressing[meshi] = -1;
+            //boundaryProcAddressing[meshi].setSize(mesh.boundaryMesh().size());
+            //boundaryProcAddressing[meshi] = -1;
+        }
+    }
+
+    forAll(meshes, meshi)
+    {
+        if (meshes.set(meshi))
+        {
+            const polyMesh& mesh = meshes[meshi];
+            const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+            const faceList& faces = mesh.faces();
+            const labelList& faceOwner = mesh.faceOwner();
+            const labelList& faceNeighbour = mesh.faceNeighbour();
+            const faceZoneMesh& faceZones = mesh.faceZones();
+
+            // Precalc offset zones
+            labelList newZoneID(faces.size(), -1);
+            boolList zoneFlip(faces.size(), false);
+
+            forAll(faceZones, zonei)
+            {
+                const labelList& faceLabels = faceZones[zonei];
+                const boolList& flipMap = faceZones[zonei].flipMap();
+
+                forAll(faceLabels, facei)
+                {
+                    newZoneID[faceLabels[facei]] = faceZoneMap[meshi][zonei];
+                    zoneFlip[faceLabels[facei]] = flipMap[facei];
+                }
+            }
+
+            // Add faces in mesh order
+
+            // 1. Internal faces
+
+            face newFace;
+            for (label facei = 0; facei < mesh.nInternalFaces(); facei++)
+            {
+                const face& f = faces[facei];
+
+                label newOwn = cellProcAddressing[meshi][faceOwner[facei]];
+                label newNei = cellProcAddressing[meshi][faceNeighbour[facei]];
+
+                newFace.setSize(f.size());
+                forAll(f, fp)
+                {
+                    newFace[fp] = pointProcAddressing[meshi][f[fp]];
+                }
+                bool flipFaceFlux = false;
+                bool flip = zoneFlip[facei];
+                if (newNei < newOwn)
+                {
+                    Swap(newOwn, newNei);
+                    newFace = newFace.reverseFace();
+                    flipFaceFlux = !flipFaceFlux;
+                    flip = !flip;
+                }
+
+                const label newFacei = meshMod.addFace
+                (
+                    newFace,
+                    newOwn,
+                    newNei,
+                    -1,                         // masterPointID
+                    -1,                         // masterEdgeID
+                    facei,                      // masterFaceID
+                    flipFaceFlux,               // flipFaceFlux
+                    -1,                         // patchID
+                    newZoneID[facei],           // zoneID
+                    flip                        // zoneFlip
+                );
+
+                faceProcAddressing[meshi][facei] = newFacei;
+            }
+
+            // 1b. Owner side of coupled faces (since owner side cells are
+            //     smallest)
+
+            const labelList& localBFaces = localBoundaryFace[meshi];
+            const labelList& procNbrs = remoteFaceMesh[meshi];
+            const labelList& procNbrBFaces = remoteBoundaryFace[meshi];
+
+            forAll(localBFaces, i)
+            {
+                const label bFacei = localBFaces[i];
+                const label nbrMeshi = procNbrs[i];
+                const label nbrBFacei = procNbrBFaces[i];
+
+                // Local mesh face
+                const label facei = mesh.nInternalFaces()+bFacei;
+                const face& f = mesh.faces()[facei];
+                const label newOwn =
+                    cellProcAddressing[meshi][faceOwner[facei]];
+
+                // Neighbour mesh face
+                const auto& nbrMesh = meshes[nbrMeshi];
+                const label nbrFacei = nbrMesh.nInternalFaces()+nbrBFacei;
+                const label nbrOwn = nbrMesh.faceOwner()[nbrFacei];
+                const label newNei = cellProcAddressing[nbrMeshi][nbrOwn];
+
+                //Pout<< "** connection between face:" << facei
+                //    << " at:" << mesh.faceCentres()[facei]
+                //    << " and nbrMesh:" << nbrMeshi
+                //    << " nbrFacei:" << nbrMesh.faceCentres()[nbrFacei]
+                //    << endl;
+
+                if (newOwn < newNei)
+                {
+                    newFace.setSize(f.size());
+                    forAll(f, fp)
+                    {
+                        newFace[fp] = pointProcAddressing[meshi][f[fp]];
+                    }
+                    const bool flipFaceFlux = false;
+                    const bool flip = zoneFlip[facei];
+                    const label newFacei = meshMod.addFace
+                    (
+                        newFace,
+                        newOwn,
+                        newNei,                     // neighbour
+                        -1,                         // masterPointID
+                        -1,                         // masterEdgeID
+                        facei,                      // masterFaceID
+                        flipFaceFlux,               // flipFaceFlux
+                        -1,                         // patchID
+                        newZoneID[facei],           // zoneID
+                        flip                        // zoneFlip
+                    );
+
+                    faceProcAddressing[meshi][facei] = newFacei;
+                    faceProcAddressing[nbrMeshi][nbrFacei] = newFacei;
+                }
+            }
+
+
+            // 2. Remaining patch faces
+
+            forAll(pbm, patchi)
+            {
+                const auto& pp = pbm[patchi];
+
+                if (patchi < patchMap[meshi].size())
+                {
+                    const label newPatchi = patchMap[meshi][patchi];
+
+                    //boundaryProcAddressing[meshi][patchi] = newPatchi;
+
+                    forAll(pp, patchFacei)
+                    {
+                        const label facei = pp.start() + patchFacei;
+                        if (faceProcAddressing[meshi][facei] == -1)
+                        {
+                            const face& f = faces[facei];
+
+                            const label newOwn =
+                                cellProcAddressing[meshi][faceOwner[facei]];
+                            newFace.setSize(f.size());
+                            forAll(f, fp)
+                            {
+                                newFace[fp] = pointProcAddressing[meshi][f[fp]];
+                            }
+
+                            const label newFacei = meshMod.addFace
+                            (
+                                newFace,
+                                newOwn,
+                                -1,                         // neighbour
+                                -1,                         // masterPointID
+                                -1,                         // masterEdgeID
+                                facei,                      // masterFaceID
+                                false,                      // flipFaceFlux
+                                newPatchi,                  // patchID
+                                newZoneID[facei],           // zoneID
+                                zoneFlip[facei]             // zoneFlip
+                            );
+
+                            faceProcAddressing[meshi][facei] = newFacei;
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/src/dynamicMesh/polyMeshAdder/polyMeshAdder.H b/src/dynamicMesh/polyMeshAdder/polyMeshAdder.H
index f59adcb5f7..2df5b03e2a 100644
--- a/src/dynamicMesh/polyMeshAdder/polyMeshAdder.H
+++ b/src/dynamicMesh/polyMeshAdder/polyMeshAdder.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -246,6 +246,33 @@ class polyMeshAdder
         );
 
 
+        //- Multi-mesh addition
+
+            //- Find n'th processor patch going to nbrProci. Usually n=0 except
+            //- for e.g. processorCyclic
+            static label procPatchIndex
+            (
+                const polyBoundaryMesh& pbm,
+                const label nbrProci,
+                const label n
+            );
+
+            //- Determine shared points. Takes pairs of one-to-one patches
+            //0 that are converted into internal faces
+            static void compactPoints
+            (
+                const UPtrList<polyMesh>& meshes,
+                const labelListList& localBoundaryFace,
+                const labelListList& remoteFaceProc,
+                const labelListList& remoteBoundaryFace,
+                const labelListList& remoteFaceStart,
+                const globalIndex& globalPoints,
+
+                labelListList& pointProcAddressing,
+                labelListList& localPoints
+            );
+
+
 public:
 
     // Member Functions
@@ -295,6 +322,54 @@ public:
                 const Map<label>& pointToMaster,
                 polyTopoChange& meshMod
             );
+
+
+        // Multi-mesh merging. Expects same patches (apart from processor
+        // patches)
+
+            //- Helper: find pairs of processor patches. Return number of
+            //- non-processor patches
+            static label procPatchPairs
+            (
+                const UPtrList<polyMesh>& meshes,
+                List<DynamicList<label>>& localPatch,
+                List<DynamicList<label>>& remoteMesh,
+                List<DynamicList<label>>& remotePatch
+            );
+
+            //- Helper: expand list of coupled patches into pairs of coupled
+            //- faces
+            static void patchFacePairs
+            (
+                const UPtrList<polyMesh>& meshes,
+                const List<DynamicList<label>>& localPatch,
+                const List<DynamicList<label>>& remoteMesh,
+                const List<DynamicList<label>>& remotePatch,
+                labelListList& localBoundaryFace,
+                labelListList& remoteFaceMesh,
+                labelListList& remoteBoundaryFace
+            );
+
+            //- Add multiple meshes
+            static void add
+            (
+                const UPtrList<polyMesh>& meshes,
+                const UList<labelList>& patchMap,
+
+                const labelListList& localBoundaryFace,
+                const labelListList& remoteFaceMesh,
+                const labelListList& remoteBoundaryFace,
+                const labelListList& remoteFaceStart,
+
+                const UList<labelList>& pointZoneMap,
+                const UList<labelList>& faceZoneMap,
+                const UList<labelList>& cellZoneMap,
+
+                polyTopoChange& meshMod,
+                labelListList& cellProcAddressing,
+                labelListList& faceProcAddressing,
+                labelListList& pointProcAddressing
+            );
 };