From cd2f614a48932c7c567bbe497205bb0948ffe301 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Thu, 31 Mar 2022 23:46:50 +0100
Subject: [PATCH] Added new polyTopoChangeMap and polyDistributionMap
 directories

---
 .../polyDistributionMap/IOdistributionMap.C   |  158 ++
 .../polyDistributionMap/IOdistributionMap.H   |   98 ++
 .../polyDistributionMap/distributionMap.C     |  556 +++++++
 .../polyDistributionMap/distributionMap.H     |  672 ++++++++
 .../polyDistributionMap/distributionMapBase.C | 1272 +++++++++++++++
 .../polyDistributionMap/distributionMapBase.H |  495 ++++++
 .../distributionMapBaseTemplates.C            | 1420 +++++++++++++++++
 .../distributionMapTemplates.C                |  271 ++++
 .../lagrangianDistributionMap.H               |  122 ++
 .../polyDistributionMap/polyDistributionMap.C |  366 +++++
 .../polyDistributionMap/polyDistributionMap.H |  312 ++++
 .../polyTopoChangeMap/cellMapper/cellMapper.C |  463 ++++++
 .../polyTopoChangeMap/cellMapper/cellMapper.H |  165 ++
 .../polyTopoChangeMap/faceMapper/faceMapper.C |  401 +++++
 .../polyTopoChangeMap/faceMapper/faceMapper.H |  181 +++
 .../polyTopoChangeMap/mapAddedPolyMesh.C      |   81 +
 .../polyTopoChangeMap/mapAddedPolyMesh.H      |  245 +++
 .../mapPatchChange/mapPatchChange.H           |  152 ++
 .../mapSubsetMesh/mapSubsetMesh.H             |  202 +++
 .../polyTopoChangeMap/morphFieldMapper.H      |   87 +
 .../polyTopoChangeMap/objectMap/objectMap.H   |  121 ++
 .../polyTopoChangeMap/objectMap/objectMapI.H  |  142 ++
 .../polyTopoChangeMap/polyTopoChangeMap.C     |  210 +++
 .../polyTopoChangeMap/polyTopoChangeMap.H     |  662 ++++++++
 24 files changed, 8854 insertions(+)
 create mode 100644 src/OpenFOAM/meshes/polyMesh/polyDistributionMap/IOdistributionMap.C
 create mode 100644 src/OpenFOAM/meshes/polyMesh/polyDistributionMap/IOdistributionMap.H
 create mode 100644 src/OpenFOAM/meshes/polyMesh/polyDistributionMap/distributionMap.C
 create mode 100644 src/OpenFOAM/meshes/polyMesh/polyDistributionMap/distributionMap.H
 create mode 100644 src/OpenFOAM/meshes/polyMesh/polyDistributionMap/distributionMapBase.C
 create mode 100644 src/OpenFOAM/meshes/polyMesh/polyDistributionMap/distributionMapBase.H
 create mode 100644 src/OpenFOAM/meshes/polyMesh/polyDistributionMap/distributionMapBaseTemplates.C
 create mode 100644 src/OpenFOAM/meshes/polyMesh/polyDistributionMap/distributionMapTemplates.C
 create mode 100644 src/OpenFOAM/meshes/polyMesh/polyDistributionMap/lagrangianDistributionMap.H
 create mode 100644 src/OpenFOAM/meshes/polyMesh/polyDistributionMap/polyDistributionMap.C
 create mode 100644 src/OpenFOAM/meshes/polyMesh/polyDistributionMap/polyDistributionMap.H
 create mode 100644 src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/cellMapper/cellMapper.C
 create mode 100644 src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/cellMapper/cellMapper.H
 create mode 100644 src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/faceMapper/faceMapper.C
 create mode 100644 src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/faceMapper/faceMapper.H
 create mode 100644 src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/mapAddedPolyMesh.C
 create mode 100644 src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/mapAddedPolyMesh.H
 create mode 100644 src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/mapPatchChange/mapPatchChange.H
 create mode 100644 src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/mapSubsetMesh/mapSubsetMesh.H
 create mode 100644 src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/morphFieldMapper.H
 create mode 100644 src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/objectMap/objectMap.H
 create mode 100644 src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/objectMap/objectMapI.H
 create mode 100644 src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/polyTopoChangeMap.C
 create mode 100644 src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/polyTopoChangeMap.H

diff --git a/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/IOdistributionMap.C b/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/IOdistributionMap.C
new file mode 100644
index 0000000000..9efc758d7c
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/IOdistributionMap.C
@@ -0,0 +1,158 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2014-2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "IOdistributionMap.H"
+
+/* * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * * */
+
+namespace Foam
+{
+    defineTypeNameAndDebug(IOdistributionMap, 0);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::IOdistributionMap::IOdistributionMap(const IOobject& io)
+:
+    regIOobject(io)
+{
+    // Temporary warning
+    if (io.readOpt() == IOobject::MUST_READ_IF_MODIFIED)
+    {
+        WarningInFunction
+            << "IOdistributionMap " << name()
+            << " constructed with IOobject::MUST_READ_IF_MODIFIED"
+            " but IOdistributionMap does not support automatic rereading."
+            << endl;
+    }
+
+    if
+    (
+        (
+            io.readOpt() == IOobject::MUST_READ
+         || io.readOpt() == IOobject::MUST_READ_IF_MODIFIED
+        )
+     || (io.readOpt() == IOobject::READ_IF_PRESENT && headerOk())
+    )
+    {
+        readStream(typeName) >> *this;
+        close();
+    }
+}
+
+
+Foam::IOdistributionMap::IOdistributionMap
+(
+    const IOobject& io,
+    const distributionMap& map
+)
+:
+    regIOobject(io)
+{
+     // Temporary warning
+    if (io.readOpt() == IOobject::MUST_READ_IF_MODIFIED)
+    {
+        WarningInFunction
+            << "IOdistributionMap " << name()
+            << " constructed with IOobject::MUST_READ_IF_MODIFIED"
+            " but IOdistributionMap does not support automatic rereading."
+            << endl;
+    }
+
+    if
+    (
+        (
+            io.readOpt() == IOobject::MUST_READ
+         || io.readOpt() == IOobject::MUST_READ_IF_MODIFIED
+        )
+     || (io.readOpt() == IOobject::READ_IF_PRESENT && headerOk())
+    )
+    {
+        readStream(typeName) >> *this;
+        close();
+    }
+    else
+    {
+        distributionMap::operator=(map);
+    }
+}
+
+
+Foam::IOdistributionMap::IOdistributionMap
+(
+    const IOobject& io,
+    distributionMap&& map
+)
+:
+    regIOobject(io),
+    distributionMap(move(map))
+{
+    // Temporary warning
+    if (io.readOpt() == IOobject::MUST_READ_IF_MODIFIED)
+    {
+        WarningInFunction
+            << "IOdistributionMap " << name()
+            << " constructed with IOobject::MUST_READ_IF_MODIFIED"
+            " but IOdistributionMap does not support automatic rereading."
+            << endl;
+    }
+
+    if
+    (
+        (
+            io.readOpt() == IOobject::MUST_READ
+         || io.readOpt() == IOobject::MUST_READ_IF_MODIFIED
+        )
+     || (io.readOpt() == IOobject::READ_IF_PRESENT && headerOk())
+    )
+    {
+        readStream(typeName) >> *this;
+        close();
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * * //
+
+Foam::IOdistributionMap::~IOdistributionMap()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::IOdistributionMap::readData(Istream& is)
+{
+    return (is >> *this).good();
+}
+
+
+bool Foam::IOdistributionMap::writeData(Ostream& os) const
+{
+    return (os << *this).good();
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/IOdistributionMap.H b/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/IOdistributionMap.H
new file mode 100644
index 0000000000..e78a31307e
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/IOdistributionMap.H
@@ -0,0 +1,98 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2014-2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::IOdistributionMap
+
+Description
+    IOdistributionMap is derived from distributionMap and
+    IOobject to give the distributionMap
+    automatic IO functionality via the objectRegistry.
+
+SourceFiles
+    IOdistributionMap.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef IOdistributionMap_H
+#define IOdistributionMap_H
+
+#include "distributionMap.H"
+#include "regIOobject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class IOdistributionMap Declaration
+\*---------------------------------------------------------------------------*/
+
+class IOdistributionMap
+:
+    public regIOobject,
+    public distributionMap
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("distributionMap");
+
+    // Constructors
+
+        //- Construct given an IOobject
+        IOdistributionMap(const IOobject&);
+
+        //- Construct given an IOobject and distributionMap
+        IOdistributionMap(const IOobject&, const distributionMap&);
+
+        //- Move constructor transferring the distributionMap contents
+        IOdistributionMap(const IOobject&, distributionMap&&);
+
+
+    //- Destructor
+    virtual ~IOdistributionMap();
+
+
+    // Member Functions
+
+        //- ReadData function required for regIOobject read operation
+        virtual bool readData(Istream&);
+
+        //- WriteData function required for regIOobject write operation
+        virtual bool writeData(Ostream&) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/distributionMap.C b/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/distributionMap.C
new file mode 100644
index 0000000000..21626cbaeb
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/distributionMap.C
@@ -0,0 +1,556 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "distributionMap.H"
+#include "globalIndexAndTransform.H"
+#include "transformField.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(distributionMap, 0);
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<>
+void Foam::distributionMap::transform::operator()
+(
+    const transformer&,
+    const bool,
+    List<label>&
+) const
+{}
+template<>
+void Foam::distributionMap::transform::operator()
+(
+    const coupledPolyPatch&,
+    UList<label>&
+) const
+{}
+template<>
+void Foam::distributionMap::transform::operator()
+(
+    const coupledPolyPatch&,
+    Map<label>&
+) const
+{}
+template<>
+void Foam::distributionMap::transform::operator()
+(
+    const coupledPolyPatch&,
+    EdgeMap<label>&
+) const
+{}
+
+
+template<>
+void Foam::distributionMap::transform::operator()
+(
+    const transformer&,
+    const bool,
+    List<scalar>&
+) const
+{}
+template<>
+void Foam::distributionMap::transform::operator()
+(
+    const coupledPolyPatch&,
+    UList<scalar>&
+) const
+{}
+template<>
+void Foam::distributionMap::transform::operator()
+(
+    const coupledPolyPatch&,
+    Map<scalar>&
+) const
+{}
+template<>
+void Foam::distributionMap::transform::operator()
+(
+    const coupledPolyPatch&,
+    EdgeMap<scalar>&
+) const
+{}
+
+
+template<>
+void Foam::distributionMap::transform::operator()
+(
+    const transformer&,
+    const bool,
+    List<bool>&
+) const
+{}
+template<>
+void Foam::distributionMap::transform::operator()
+(
+    const coupledPolyPatch&,
+    UList<bool>&
+) const
+{}
+template<>
+void Foam::distributionMap::transform::operator()
+(
+    const coupledPolyPatch&,
+    Map<bool>&
+) const
+{}
+template<>
+void Foam::distributionMap::transform::operator()
+(
+    const coupledPolyPatch&,
+    EdgeMap<bool>&
+) const
+{}
+
+
+void Foam::distributionMap::printLayout(Ostream& os) const
+{
+    distributionMapBase::printLayout(os);
+
+    forAll(transformElements_, trafoI)
+    {
+        if (transformElements_[trafoI].size() > 0)
+        {
+            os  << "transform " << trafoI << ':' << endl
+                << "    start : " << transformStart_[trafoI] << endl
+                << "    size  : " << transformElements_[trafoI].size() << endl;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::distributionMap::distributionMap()
+:
+    distributionMapBase()
+{}
+
+
+Foam::distributionMap::distributionMap
+(
+    const label constructSize,
+    labelListList&& subMap,
+    labelListList&& constructMap,
+    const bool subHasFlip,
+    const bool constructHasFlip
+)
+:
+    distributionMapBase
+    (
+        constructSize,
+        move(subMap),
+        move(constructMap),
+        subHasFlip,
+        constructHasFlip
+    )
+{}
+
+
+Foam::distributionMap::distributionMap
+(
+    const label constructSize,
+    labelListList&& subMap,
+    labelListList&& constructMap,
+    labelListList&& transformElements,
+    labelList&& transformStart,
+    const bool subHasFlip,
+    const bool constructHasFlip
+)
+:
+    distributionMapBase
+    (
+        constructSize,
+        move(subMap),
+        move(constructMap),
+        subHasFlip,
+        constructHasFlip
+    ),
+    transformElements_(move(transformElements)),
+    transformStart_(move(transformStart))
+{}
+
+
+Foam::distributionMap::distributionMap
+(
+    const labelList& sendProcs,
+    const labelList& recvProcs
+)
+:
+    distributionMapBase(sendProcs, recvProcs)
+{}
+
+
+Foam::distributionMap::distributionMap
+(
+    const globalIndex& globalNumbering,
+    labelList& elements,
+    List<Map<label>>& compactMap,
+    const int tag
+)
+:
+    distributionMapBase
+    (
+        globalNumbering,
+        elements,
+        compactMap,
+        tag
+    )
+{}
+
+
+Foam::distributionMap::distributionMap
+(
+    const globalIndex& globalNumbering,
+    labelListList& cellCells,
+    List<Map<label>>& compactMap,
+    const int tag
+)
+:
+    distributionMapBase
+    (
+        globalNumbering,
+        cellCells,
+        compactMap,
+        tag
+    )
+{}
+
+
+Foam::distributionMap::distributionMap
+(
+    const globalIndex& globalNumbering,
+    labelList& elements,
+    const globalIndexAndTransform& globalTransforms,
+    const labelPairList& transformedElements,
+    labelList& transformedIndices,
+    List<Map<label>>& compactMap,
+    const int tag
+)
+:
+    distributionMapBase()
+{
+    // Construct per processor compact addressing of the global elements
+    // needed. The ones from the local processor are not included since
+    // these are always all needed.
+    calcCompactAddressing
+    (
+        globalNumbering,
+        elements,
+        compactMap
+    );
+
+    // Add all (non-local) transformed elements needed.
+    forAll(transformedElements, i)
+    {
+        labelPair elem = transformedElements[i];
+        label proci = globalTransforms.processor(elem);
+        if (proci != Pstream::myProcNo())
+        {
+            label index = globalTransforms.index(elem);
+            label nCompact = compactMap[proci].size();
+            compactMap[proci].insert(index, nCompact);
+        }
+    }
+
+
+    // Exchange what I need with processor that supplies it. Renumber elements
+    // into compact numbering
+    labelList compactStart;
+    exchangeAddressing
+    (
+        tag,
+        globalNumbering,
+        elements,
+        compactMap,
+        compactStart
+    );
+
+
+    // Renumber the transformed elements
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Count per transformIndex
+    label nTrafo = globalTransforms.transformPermutations().size();
+    labelList nPerTransform(nTrafo, 0);
+    forAll(transformedElements, i)
+    {
+        labelPair elem = transformedElements[i];
+        label trafoI = globalTransforms.transformIndex(elem);
+        nPerTransform[trafoI]++;
+    }
+    // Offset per transformIndex
+    transformStart_.setSize(nTrafo);
+    transformElements_.setSize(nTrafo);
+    forAll(transformStart_, trafoI)
+    {
+        transformStart_[trafoI] = constructSize_;
+        constructSize_ += nPerTransform[trafoI];
+        transformElements_[trafoI].setSize(nPerTransform[trafoI]);
+    }
+
+    // Sort transformed elements into their new slot.
+    nPerTransform = 0;
+
+    transformedIndices.setSize(transformedElements.size());
+    forAll(transformedElements, i)
+    {
+        labelPair elem = transformedElements[i];
+        label proci = globalTransforms.processor(elem);
+        label index = globalTransforms.index(elem);
+        label trafoI = globalTransforms.transformIndex(elem);
+
+        // Get compact index for untransformed element
+        label rawElemI =
+        (
+            proci == Pstream::myProcNo()
+          ? index
+          : compactMap[proci][index]
+        );
+
+        label& n = nPerTransform[trafoI];
+        // index of element to transform
+        transformElements_[trafoI][n] = rawElemI;
+        // destination of transformed element
+        transformedIndices[i] = transformStart_[trafoI]+n;
+        n++;
+    }
+
+    if (debug)
+    {
+        printLayout(Pout);
+    }
+}
+
+
+Foam::distributionMap::distributionMap
+(
+    const globalIndex& globalNumbering,
+    labelListList& cellCells,
+    const globalIndexAndTransform& globalTransforms,
+    const List<labelPairList>& transformedElements,
+    labelListList& transformedIndices,
+    List<Map<label>>& compactMap,
+    const int tag
+)
+:
+    distributionMapBase()
+{
+    // Construct per processor compact addressing of the global elements
+    // needed. The ones from the local processor are not included since
+    // these are always all needed.
+    calcCompactAddressing
+    (
+        globalNumbering,
+        cellCells,
+        compactMap
+    );
+
+    // Add all (non-local) transformed elements needed.
+    forAll(transformedElements, celli)
+    {
+        const labelPairList& elems = transformedElements[celli];
+
+        forAll(elems, i)
+        {
+            label proci = globalTransforms.processor(elems[i]);
+            if (proci != Pstream::myProcNo())
+            {
+                label index = globalTransforms.index(elems[i]);
+                label nCompact = compactMap[proci].size();
+                compactMap[proci].insert(index, nCompact);
+            }
+        }
+    }
+
+
+    // Exchange what I need with processor that supplies it. Renumber elements
+    // into compact numbering
+    labelList compactStart;
+    exchangeAddressing
+    (
+        tag,
+        globalNumbering,
+        cellCells,
+        compactMap,
+        compactStart
+    );
+
+
+    // Renumber the transformed elements
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Count per transformIndex
+    label nTrafo = globalTransforms.transformPermutations().size();
+    labelList nPerTransform(nTrafo, 0);
+    forAll(transformedElements, celli)
+    {
+        const labelPairList& elems = transformedElements[celli];
+
+        forAll(elems, i)
+        {
+            label trafoI = globalTransforms.transformIndex(elems[i]);
+            nPerTransform[trafoI]++;
+        }
+    }
+    // Offset per transformIndex
+    transformStart_.setSize(nTrafo);
+    transformElements_.setSize(nTrafo);
+    forAll(transformStart_, trafoI)
+    {
+        transformStart_[trafoI] = constructSize_;
+        constructSize_ += nPerTransform[trafoI];
+        transformElements_[trafoI].setSize(nPerTransform[trafoI]);
+    }
+
+    // Sort transformed elements into their new slot.
+    nPerTransform = 0;
+
+    transformedIndices.setSize(transformedElements.size());
+    forAll(transformedElements, celli)
+    {
+        const labelPairList& elems = transformedElements[celli];
+        transformedIndices[celli].setSize(elems.size());
+
+        forAll(elems, i)
+        {
+            label proci = globalTransforms.processor(elems[i]);
+            label index = globalTransforms.index(elems[i]);
+            label trafoI = globalTransforms.transformIndex(elems[i]);
+
+            // Get compact index for untransformed element
+            label rawElemI =
+            (
+                proci == Pstream::myProcNo()
+              ? index
+              : compactMap[proci][index]
+            );
+
+            label& n = nPerTransform[trafoI];
+            // index of element to transform
+            transformElements_[trafoI][n] = rawElemI;
+            // destination of transformed element
+            transformedIndices[celli][i] = transformStart_[trafoI]+n;
+            n++;
+        }
+    }
+
+    if (debug)
+    {
+        printLayout(Pout);
+    }
+}
+
+
+Foam::distributionMap::distributionMap(const distributionMap& map)
+:
+    distributionMapBase(map),
+    transformElements_(map.transformElements_),
+    transformStart_(map.transformStart_)
+{}
+
+
+Foam::distributionMap::distributionMap(distributionMap&& map)
+:
+    distributionMapBase(move(map)),
+    transformElements_(move(map.transformElements_)),
+    transformStart_(move(map.transformStart_))
+{}
+
+
+Foam::distributionMap::distributionMap(Istream& is)
+{
+    is  >> *this;
+}
+
+
+Foam::autoPtr<Foam::distributionMap> Foam::distributionMap::clone() const
+{
+    return autoPtr<distributionMap>(new distributionMap(*this));
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::label Foam::distributionMap::whichTransform(const label index)
+const
+{
+    return findLower(transformStart_, index+1);
+}
+
+
+void Foam::distributionMap::transfer(distributionMap& rhs)
+{
+    distributionMapBase::transfer(rhs);
+    transformElements_.transfer(rhs.transformElements_);
+    transformStart_.transfer(rhs.transformStart_);
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+void Foam::distributionMap::operator=(const distributionMap& rhs)
+{
+    // Check for assignment to self
+    if (this == &rhs)
+    {
+        FatalErrorInFunction
+            << "Attempted assignment to self"
+            << abort(FatalError);
+    }
+    distributionMapBase::operator=(rhs);
+    transformElements_ = rhs.transformElements_;
+    transformStart_ = rhs.transformStart_;
+}
+
+
+// * * * * * * * * * * * * * * Istream Operator  * * * * * * * * * * * * * * //
+
+Foam::Istream& Foam::operator>>(Istream& is, distributionMap& map)
+{
+    is.fatalCheck("operator>>(Istream&, distributionMap&)");
+
+    is  >> static_cast<distributionMapBase&>(map)
+        >> map.transformElements_ >> map.transformStart_;
+
+    return is;
+}
+
+
+// * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const distributionMap& map)
+{
+    os  << static_cast<const distributionMapBase&>(map) << token::NL
+        << map.transformElements_ << token::NL
+        << map.transformStart_;
+
+    return os;
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/distributionMap.H b/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/distributionMap.H
new file mode 100644
index 0000000000..95eab374f0
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/distributionMap.H
@@ -0,0 +1,672 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::distributionMap
+
+Description
+    Class containing processor-to-processor mapping information.
+
+    We store mapping from the bits-to-send to the complete starting list
+    (subXXXMap) and from the received bits to their location in the new
+    list (constructXXXMap).
+
+Note:
+    Schedule is a list of processor pairs (one send, one receive. One of
+    them will be myself) which forms a scheduled (i.e. non-buffered) exchange.
+    See distribute on how to use it.
+    Note2: number of items sent on one processor have to equal the number
+    of items received on the other processor.
+
+    To aid constructing these maps there are the constructors from global
+    numbering, either with or without transforms.
+
+    - without transforms:
+    Constructors using compact numbering: layout is
+    - all my own elements first (whether used or not)
+    - followed by used-only remote elements sorted by remote processor.
+    So e.g 4 procs and on proc 1 the compact
+    table will first have all globalIndex.localSize() elements from proc1
+    followed by used-only elements of proc0, proc2, proc3.
+    The constructed distributionMap sends the local elements from and
+    receives the remote elements into their compact position.
+    compactMap[proci] is the position of elements from proci in the compact
+    map. compactMap[myProcNo()] is empty since trivial addressing.
+
+    It rewrites the input global indices into indices into the constructed
+    data.
+
+
+    - with transforms:
+    This requires the precalculated set of possible transforms
+    (globalIndexAndTransform). These are given as permutations (+, -, or none)
+    of up to 3 independent transforms.
+    The layout of the data is
+    - all my own elements first (whether used or not)
+    - followed by used-only remote elements sorted by remote processor.
+    - followed by - for each transformation index - the set of local or
+    remote elements with that transformation.
+    The inputs for the constructor are
+    - the set of untransformed local or remote indices in globalIndex
+    numbering. These get rewritten to be indices into the layout of the data.
+    - the set of transformed local or remote indices in globalIndexAndTransform
+    encoding. These are labelPairs.
+
+    Any distribute with transforms is now done as:
+    1. exchange data with other processors and receive these into the
+    slots for that processor
+    2. for all transformations transform a subset of the data according
+    to transformElements_[transformI] and store this starting from
+    transformStart_[transformI]
+
+    In the same way a reverse distribute will
+    1. apply the inverse transform to the data starting at
+    transformStart_[transformI] and copy the result back into the
+    transformElements_[transformI]. These might be local or remote slots.
+    2. the data in the remote slots will now be sent back to the correct
+    location in the originating processor.
+
+    E.g. a map to handle
+    - mesh points on a mesh with
+    - 1 cyclic so 3 permutations (+,-,none) will have layout
+    - on e.g. processor 1 out of 2:
+
+        +------+ <- transformStart[2]
+        |      |
+        |      | <- transform2 applied to data in local or remote slots
+        |      |
+        +------+ <- transformStart[1]
+        |      |
+        |      | <- transform1 applied to data in local or remote slots
+        |      |
+        +------+ <- transformStart[1]
+        |      |
+        |      | <- transform0 applied to data in local or remote slots
+        |      |
+        +------+ <- transformStart[0]
+        |      |
+        |      | <- data from proc2
+        |      |
+        +------+
+        |      |
+        |      | <- data from proc0
+        |      |
+        +------+ <- mesh.nPoints()
+        |      |
+        |      |
+        |      |
+        +------+ 0
+
+
+    When constructing from components optionally a 'flip' on
+    the maps can be specified. This will interpret the map
+    values as index+flip, similar to e.g. faceProcAddressing. The flip
+    will only be applied to fieldTypes (scalar, vector, .. triad)
+
+
+SourceFiles
+    distributionMap.C
+    distributionMapTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef distributionMap_H
+#define distributionMap_H
+
+#include "distributionMapBase.H"
+#include "transformer.H"
+#include "coupledPolyPatch.H"
+#include "EdgeMap.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class globalIndexAndTransform;
+
+
+// Forward declaration of friend functions and operators
+
+class distributionMap;
+
+Istream& operator>>(Istream&, distributionMap&);
+Ostream& operator<<(Ostream&, const distributionMap&);
+
+
+/*---------------------------------------------------------------------------*\
+                        Class distributionMap Declaration
+\*---------------------------------------------------------------------------*/
+
+class distributionMap
+:
+    public distributionMapBase
+{
+    // Private Data
+
+        //- For every globalIndexAndTransform::transformPermutations
+        //  gives the elements that need to be transformed
+        labelListList transformElements_;
+
+        //- Destination in constructMap for transformed elements
+        labelList transformStart_;
+
+
+    // Private Member Functions
+
+        //- Helper function: copy transformElements without transformation
+        template<class T>
+        void applyDummyTransforms(List<T>& field) const;
+
+        template<class T, class TransformOp>
+        void applyTransforms
+        (
+            const globalIndexAndTransform& globalTransforms,
+            List<T>& field,
+            const TransformOp& top
+        ) const;
+
+        //- Helper function: copy transformElements without transformation
+        template<class T>
+        void applyDummyInverseTransforms(List<T>& field) const;
+
+        template<class T, class TransformOp>
+        void applyInverseTransforms
+        (
+            const globalIndexAndTransform& globalTransforms,
+            List<T>& field,
+            const TransformOp& top
+        ) const;
+
+
+public:
+
+    // Public classes
+
+        //- Default transformation behaviour
+        class transform
+        {
+        public:
+
+            template<class Type>
+            void operator()
+            (
+                const transformer& vt,
+                const bool forward,
+                List<Type>& fld
+            ) const
+            {
+                if (forward)
+                {
+                    vt.transformList(fld);
+                }
+                else
+                {
+                    vt.invTransformList(fld);
+                }
+            }
+
+            template<class Type>
+            void operator()
+            (
+                const transformer& vt,
+                const bool forward,
+                List<List<Type>>& flds
+            ) const
+            {
+                forAll(flds, i)
+                {
+                    operator()(vt, forward, flds[i]);
+                }
+            }
+
+            //- Transform patch-based field
+            template<class Type>
+            void operator()(const coupledPolyPatch& cpp, UList<Type>& fld) const
+            {
+                if (cpp.transform().transforms())
+                {
+                    cpp.transform().transformList(fld);
+                }
+            }
+
+            //- Transform sparse field
+            template<class Type, template<class> class Container>
+            void operator()(const coupledPolyPatch& cpp, Container<Type>& map)
+            const
+            {
+                if (cpp.transform().transforms())
+                {
+                    cpp.transform().transformList(map);
+                }
+            }
+        };
+
+        //- Default transformation behaviour for position
+        class transformPosition
+        {
+        public:
+
+            void operator()
+            (
+                const transformer& vt,
+                const bool forward,
+                List<point>& fld
+            ) const
+            {
+                pointField pfld(move(fld));
+                if (forward)
+                {
+                    fld = vt.transformPosition(pfld);
+                }
+                else
+                {
+                    fld = vt.invTransformPosition(pfld);
+                }
+            }
+
+            void operator()
+            (
+                const transformer& vt,
+                const bool forward,
+                List<List<point>>& flds
+            ) const
+            {
+                forAll(flds, i)
+                {
+                    operator()(vt, forward, flds[i]);
+                }
+            }
+
+            //- Transform patch-based field
+            void operator()(const coupledPolyPatch& cpp, pointField& fld) const
+            {
+                cpp.transform().transformPosition(fld, fld);
+            }
+
+            template<template<class> class Container>
+            void operator()(const coupledPolyPatch& cpp, Container<point>& map)
+            const
+            {
+                Field<point> fld(map.size());
+
+                label i = 0;
+                forAllConstIter(typename Container<point>, map, iter)
+                {
+                    fld[i++] = iter();
+                }
+
+                cpp.transform().transformPosition(fld, fld);
+
+                i = 0;
+                forAllIter(typename Container<point>, map, iter)
+                {
+                    iter() = fld[i++];
+                }
+            }
+        };
+
+
+    // Declare name of the class and its debug switch
+    ClassName("distributionMap");
+
+
+    // Constructors
+
+        //- Construct null
+        distributionMap();
+
+        //- Construct from components
+        distributionMap
+        (
+            const label constructSize,
+            labelListList&& subMap,
+            labelListList&& constructMap,
+            const bool subHasFlip = false,
+            const bool constructHasFlip = false
+        );
+
+        //- Construct from components
+        distributionMap
+        (
+            const label constructSize,
+            labelListList&& subMap,
+            labelListList&& constructMap,
+            labelListList&& transformElements,
+            labelList&& transformStart,
+            const bool subHasFlip = false,
+            const bool constructHasFlip = false
+        );
+
+        //- Construct from reverse addressing: per data item the send
+        //  processor and the receive processor. (note: data is not stored
+        //  sorted per processor so cannot use printLayout).
+        distributionMap
+        (
+            const labelList& sendProcs,
+            const labelList& recvProcs
+        );
+
+        //- Construct from list of (possibly) remote elements in globalIndex
+        //  numbering (or -1). Determines compact numbering (see above) and
+        //  distribute map to get data into this ordering and renumbers the
+        //  elements to be in compact numbering.
+        distributionMap
+        (
+            const globalIndex&,
+            labelList& elements,
+            List<Map<label>>& compactMap,
+            const int tag = Pstream::msgType()
+        );
+
+        //- Special variant that works with the info sorted into bins
+        //  according to local indices. E.g. think cellCells where
+        //  cellCells[localCelli] is a list of global cells
+        distributionMap
+        (
+            const globalIndex&,
+            labelListList& cellCells,
+            List<Map<label>>& compactMap,
+            const int tag = Pstream::msgType()
+        );
+
+        //- Construct from list of (possibly remote) untransformed elements
+        //  in globalIndex numbering (or -1) and (possibly remote)
+        //  transformded elements in globalIndexAndTransform numbering.
+        //  Determines compact numbering (see above) and
+        //  distribute map to get data into this ordering and renumbers the
+        //  elements to be in compact numbering.
+        distributionMap
+        (
+            const globalIndex&,
+            labelList& untransformedElements,
+            const globalIndexAndTransform&,
+            const labelPairList& transformedElements,
+            labelList& transformedIndices,
+            List<Map<label>>& compactMap,
+            const int tag = Pstream::msgType()
+        );
+
+        //- As above but with ListLists.
+        distributionMap
+        (
+            const globalIndex&,
+            labelListList& cellCells,
+            const globalIndexAndTransform&,
+            const List<labelPairList>& transformedElements,
+            labelListList& transformedIndices,
+            List<Map<label>>& compactMap,
+            const int tag = Pstream::msgType()
+        );
+
+        //- Copy constructor
+        distributionMap(const distributionMap&);
+
+        //- Move constructor
+        distributionMap(distributionMap&&);
+
+        //- Construct from Istream
+        distributionMap(Istream&);
+
+        //- Clone
+        autoPtr<distributionMap> clone() const;
+
+
+    //- Destructor
+    virtual ~distributionMap()
+    {}
+
+
+    // Member Functions
+
+        // Access
+
+            //- For every globalIndexAndTransform::transformPermutations
+            //  gives the elements that need to be transformed
+            const labelListList& transformElements() const
+            {
+                return transformElements_;
+            }
+
+            //- Destination in constructMap for transformed elements
+            const labelList& transformStart() const
+            {
+                return transformStart_;
+            }
+
+            //- Find transform from transformElements
+            label whichTransform(const label index) const;
+
+
+        // Other
+
+            //- Transfer the contents of the argument and annul the argument.
+            void transfer(distributionMap&);
+
+            //- Distribute data using default commsType.
+            template<class T>
+            void distribute
+            (
+                List<T>& fld,
+                const bool dummyTransform = true,
+                const int tag = UPstream::msgType()
+            ) const;
+
+            //- Distribute data using default commsType.
+            template<class T, class negateOp>
+            void distribute
+            (
+                List<T>& fld,
+                const negateOp& negOp,
+                const bool dummyTransform = true,
+                const int tag = UPstream::msgType()
+            ) const;
+
+            //- Distribute data using default commsType.
+            template<class T>
+            void distribute
+            (
+                DynamicList<T>& fld,
+                const bool dummyTransform = true,
+                const int tag = UPstream::msgType()
+            ) const;
+
+            //- Reverse distribute data using default commsType.
+            template<class T>
+            void reverseDistribute
+            (
+                const label constructSize,
+                List<T>&,
+                const bool dummyTransform = true,
+                const int tag = UPstream::msgType()
+            ) const;
+
+            //- Reverse distribute data using default commsType.
+            //  Since constructSize might be larger than supplied size supply
+            //  a nullValue
+            template<class T>
+            void reverseDistribute
+            (
+                const label constructSize,
+                const T& nullValue,
+                List<T>& fld,
+                const bool dummyTransform = true,
+                const int tag = UPstream::msgType()
+            ) const;
+
+            //- Distribute with transforms
+            template<class T, class TransformOp>
+            void distribute
+            (
+                const globalIndexAndTransform&,
+                List<T>& fld,
+                const TransformOp& top,
+                const int tag = UPstream::msgType()
+            ) const;
+
+            //- Reverse distribute with transforms
+            template<class T, class TransformOp>
+            void reverseDistribute
+            (
+                const globalIndexAndTransform&,
+                const label constructSize,
+                List<T>& fld,
+                const TransformOp& top,
+                const int tag = UPstream::msgType()
+            ) const;
+
+            //- Reverse distribute with transforms
+            template<class T, class TransformOp>
+            void reverseDistribute
+            (
+                const globalIndexAndTransform&,
+                const label constructSize,
+                const T& nullValue,
+                List<T>& fld,
+                const TransformOp& top,
+                const int tag = UPstream::msgType()
+            ) const;
+
+            //- Debug: print layout. Can only be used on maps with sorted
+            //  storage (local data first, then non-local data)
+            void printLayout(Ostream& os) const;
+
+            //- Correct for topo change.
+            void updateMesh(const polyTopoChangeMap&)
+            {
+                notImplemented
+                (
+                    "distributionMap::updateMesh(const polyTopoChangeMap&)"
+                );
+            }
+
+    // Member Operators
+
+        void operator=(const distributionMap&);
+
+    // IOstream Operators
+
+        //- Read dictionary from Istream
+        friend Istream& operator>>(Istream&, distributionMap&);
+
+        //- Write dictionary to Ostream
+        friend Ostream& operator<<(Ostream&, const distributionMap&);
+
+};
+
+
+// Template specialisation for primitives that do not need transform
+template<>
+void distributionMap::transform::operator()
+(
+    const transformer&,
+    const bool,
+    List<label>&
+) const;
+template<>
+void distributionMap::transform::operator()
+(
+    const coupledPolyPatch&,
+    UList<label>&
+) const;
+template<>
+void distributionMap::transform::operator()
+(
+    const coupledPolyPatch&,
+    Map<label>&
+) const;
+template<>
+void distributionMap::transform::operator()
+(
+    const coupledPolyPatch&,
+    EdgeMap<label>&
+) const;
+
+template<>
+void distributionMap::transform::operator()
+(
+    const coupledPolyPatch&,
+    UList<scalar>&
+) const;
+template<>
+void distributionMap::transform::operator()
+(
+    const transformer&,
+    const bool,
+    List<scalar>&
+) const;
+template<>
+void distributionMap::transform::operator()
+(
+    const coupledPolyPatch&,
+    Map<scalar>&
+) const;
+template<>
+void distributionMap::transform::operator()
+(
+    const coupledPolyPatch&,
+    EdgeMap<scalar>&
+) const;
+
+template<>
+void distributionMap::transform::operator()
+(
+    const coupledPolyPatch& cpp,
+    UList<bool>& fld
+) const;
+template<>
+void distributionMap::transform::operator()
+(
+    const transformer&,
+    const bool,
+    List<bool>&
+) const;
+template<>
+void distributionMap::transform::operator()
+(
+    const coupledPolyPatch&,
+    Map<bool>&
+) const;
+template<>
+void distributionMap::transform::operator()
+(
+    const coupledPolyPatch&,
+    EdgeMap<bool>&
+) const;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "distributionMapTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/distributionMapBase.C b/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/distributionMapBase.C
new file mode 100644
index 0000000000..511b5e66b2
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/distributionMapBase.C
@@ -0,0 +1,1272 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2015-2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "distributionMapBase.H"
+#include "commSchedule.H"
+#include "HashSet.H"
+#include "globalIndex.H"
+#include "ListOps.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(distributionMapBase, 0);
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+Foam::List<Foam::labelPair> Foam::distributionMapBase::schedule
+(
+    const labelListList& subMap,
+    const labelListList& constructMap,
+    const int tag
+)
+{
+    // Communications: send and receive processor
+    List<labelPair> allComms;
+
+    {
+        HashSet<labelPair, labelPair::Hash<>> commsSet(Pstream::nProcs());
+
+        // Find what communication is required
+        forAll(subMap, proci)
+        {
+            if (proci != Pstream::myProcNo())
+            {
+                if (subMap[proci].size())
+                {
+                    // I need to send to proci
+                    commsSet.insert(labelPair(Pstream::myProcNo(), proci));
+                }
+                if (constructMap[proci].size())
+                {
+                    // I need to receive from proci
+                    commsSet.insert(labelPair(proci, Pstream::myProcNo()));
+                }
+            }
+        }
+        allComms = commsSet.toc();
+    }
+
+
+    // Reduce
+    if (Pstream::master())
+    {
+        // Receive and merge
+        for
+        (
+            int slave=Pstream::firstSlave();
+            slave<=Pstream::lastSlave();
+            slave++
+        )
+        {
+            IPstream fromSlave(Pstream::commsTypes::scheduled, slave, 0, tag);
+            List<labelPair> nbrData(fromSlave);
+
+            forAll(nbrData, i)
+            {
+                if (findIndex(allComms, nbrData[i]) == -1)
+                {
+                    label sz = allComms.size();
+                    allComms.setSize(sz+1);
+                    allComms[sz] = nbrData[i];
+                }
+            }
+        }
+        // Send back
+        for
+        (
+            int slave=Pstream::firstSlave();
+            slave<=Pstream::lastSlave();
+            slave++
+        )
+        {
+            OPstream toSlave(Pstream::commsTypes::scheduled, slave, 0, tag);
+            toSlave << allComms;
+        }
+    }
+    else
+    {
+        {
+            OPstream toMaster
+            (
+                Pstream::commsTypes::scheduled,
+                Pstream::masterNo(),
+                0,
+                tag
+            );
+            toMaster << allComms;
+        }
+        {
+            IPstream fromMaster
+            (
+                Pstream::commsTypes::scheduled,
+                Pstream::masterNo(),
+                0,
+                tag
+            );
+            fromMaster >> allComms;
+        }
+    }
+
+
+    // Determine my schedule.
+    labelList mySchedule
+    (
+        commSchedule
+        (
+            Pstream::nProcs(),
+            allComms
+        ).procSchedule()[Pstream::myProcNo()]
+    );
+
+    // Processors involved in my schedule
+    return List<labelPair>(UIndirectList<labelPair>(allComms, mySchedule));
+
+
+    // if (debug)
+    //{
+    //    Pout<< "I need to:" << endl;
+    //    const List<labelPair>& comms = schedule();
+    //    forAll(comms, i)
+    //    {
+    //        const labelPair& twoProcs = comms[i];
+    //        label sendProc = twoProcs[0];
+    //        label recvProc = twoProcs[1];
+    //
+    //        if (recvProc == Pstream::myProcNo())
+    //        {
+    //            Pout<< "    receive from " << sendProc << endl;
+    //        }
+    //        else
+    //        {
+    //            Pout<< "    send to " << recvProc << endl;
+    //        }
+    //    }
+    //}
+}
+
+
+const Foam::List<Foam::labelPair>& Foam::distributionMapBase::schedule() const
+{
+    if (schedulePtr_.empty())
+    {
+        schedulePtr_.reset
+        (
+            new List<labelPair>
+            (
+                schedule(subMap_, constructMap_, Pstream::msgType())
+            )
+        );
+    }
+    return schedulePtr_();
+}
+
+
+void Foam::distributionMapBase::checkReceivedSize
+(
+    const label proci,
+    const label expectedSize,
+    const label receivedSize
+)
+{
+    if (receivedSize != expectedSize)
+    {
+        FatalErrorInFunction
+            << "Expected from processor " << proci
+            << " " << expectedSize << " but received "
+            << receivedSize << " elements."
+            << abort(FatalError);
+    }
+}
+
+
+void Foam::distributionMapBase::printLayout(Ostream& os) const
+{
+    // Determine offsets of remote data.
+    labelList minIndex(Pstream::nProcs(), labelMax);
+    labelList maxIndex(Pstream::nProcs(), labelMin);
+    forAll(constructMap_, proci)
+    {
+        const labelList& construct = constructMap_[proci];
+        if (constructHasFlip_)
+        {
+            forAll(construct, i)
+            {
+                label index = mag(construct[i])-1;
+                minIndex[proci] = min(minIndex[proci], index);
+                maxIndex[proci] = max(maxIndex[proci], index);
+            }
+        }
+        else
+        {
+            forAll(construct, i)
+            {
+                label index = construct[i];
+                minIndex[proci] = min(minIndex[proci], index);
+                maxIndex[proci] = max(maxIndex[proci], index);
+            }
+        }
+    }
+
+    label localSize;
+    if (maxIndex[Pstream::myProcNo()] == labelMin)
+    {
+        localSize = 0;
+    }
+    else
+    {
+        localSize = maxIndex[Pstream::myProcNo()]+1;
+    }
+
+    os  << "Layout: (constructSize:" << constructSize_
+        << " subHasFlip:" << subHasFlip_
+        << " constructHasFlip:" << constructHasFlip_
+        << ")" << endl
+        << "local (processor " << Pstream::myProcNo() << "):" << endl
+        << "    start : 0" << endl
+        << "    size  : " << localSize << endl;
+
+    label offset = localSize;
+    forAll(minIndex, proci)
+    {
+        if (proci != Pstream::myProcNo())
+        {
+            if (constructMap_[proci].size() > 0)
+            {
+                if (minIndex[proci] != offset)
+                {
+                    FatalErrorInFunction
+                        << "offset:" << offset
+                        << " proci:" << proci
+                        << " minIndex:" << minIndex[proci]
+                        << abort(FatalError);
+                }
+
+                label size = maxIndex[proci]-minIndex[proci]+1;
+                os  << "processor " << proci << ':' << endl
+                    << "    start : " << offset << endl
+                    << "    size  : " << size << endl;
+
+                offset += size;
+            }
+        }
+    }
+}
+
+
+void Foam::distributionMapBase::calcCompactAddressing
+(
+    const globalIndex& globalNumbering,
+    const labelList& elements,
+    List<Map<label>>& compactMap
+) const
+{
+    compactMap.setSize(Pstream::nProcs());
+
+    // Count all (non-local) elements needed. Just for presizing map.
+    labelList nNonLocal(Pstream::nProcs(), 0);
+
+    forAll(elements, i)
+    {
+        label globalIndex = elements[i];
+
+        if (globalIndex != -1 && !globalNumbering.isLocal(globalIndex))
+        {
+            label proci = globalNumbering.whichProcID(globalIndex);
+            nNonLocal[proci]++;
+        }
+    }
+
+    forAll(compactMap, proci)
+    {
+        compactMap[proci].clear();
+        if (proci != Pstream::myProcNo())
+        {
+            compactMap[proci].resize(2*nNonLocal[proci]);
+        }
+    }
+
+
+    // Collect all (non-local) elements needed.
+    forAll(elements, i)
+    {
+        label globalIndex = elements[i];
+
+        if (globalIndex != -1 && !globalNumbering.isLocal(globalIndex))
+        {
+            label proci = globalNumbering.whichProcID(globalIndex);
+            label index = globalNumbering.toLocal(proci, globalIndex);
+            label nCompact = compactMap[proci].size();
+            compactMap[proci].insert(index, nCompact);
+        }
+    }
+}
+
+
+void Foam::distributionMapBase::calcCompactAddressing
+(
+    const globalIndex& globalNumbering,
+    const labelListList& cellCells,
+    List<Map<label>>& compactMap
+) const
+{
+    compactMap.setSize(Pstream::nProcs());
+
+    // Count all (non-local) elements needed. Just for presizing map.
+    labelList nNonLocal(Pstream::nProcs(), 0);
+
+    forAll(cellCells, cellI)
+    {
+        const labelList& cCells = cellCells[cellI];
+
+        forAll(cCells, i)
+        {
+            label globalIndex = cCells[i];
+
+            if (globalIndex != -1 && !globalNumbering.isLocal(globalIndex))
+            {
+                label proci = globalNumbering.whichProcID(globalIndex);
+                nNonLocal[proci]++;
+            }
+        }
+    }
+
+    forAll(compactMap, proci)
+    {
+        compactMap[proci].clear();
+        if (proci != Pstream::myProcNo())
+        {
+            compactMap[proci].resize(2*nNonLocal[proci]);
+        }
+    }
+
+
+    // Collect all (non-local) elements needed.
+    forAll(cellCells, cellI)
+    {
+        const labelList& cCells = cellCells[cellI];
+
+        forAll(cCells, i)
+        {
+            label globalIndex = cCells[i];
+
+            if (globalIndex != -1 && !globalNumbering.isLocal(globalIndex))
+            {
+                label proci = globalNumbering.whichProcID(globalIndex);
+                label index = globalNumbering.toLocal(proci, globalIndex);
+                label nCompact = compactMap[proci].size();
+                compactMap[proci].insert(index, nCompact);
+            }
+        }
+    }
+}
+
+
+void Foam::distributionMapBase::exchangeAddressing
+(
+    const int tag,
+    const globalIndex& globalNumbering,
+    labelList& elements,
+    List<Map<label>>& compactMap,
+    labelList& compactStart
+)
+{
+    // The overall compact addressing is
+    // - myProcNo data first (uncompacted)
+    // - all other processors consecutively
+
+    compactStart.setSize(Pstream::nProcs());
+    compactStart[Pstream::myProcNo()] = 0;
+    constructSize_ = globalNumbering.localSize();
+    forAll(compactStart, proci)
+    {
+        if (proci != Pstream::myProcNo())
+        {
+            compactStart[proci] = constructSize_;
+            constructSize_ += compactMap[proci].size();
+        }
+    }
+
+
+
+    // Find out what to receive/send in compact addressing.
+
+    // What I want to receive is what others have to send
+    labelListList wantedRemoteElements(Pstream::nProcs());
+    // Compact addressing for received data
+    constructMap_.setSize(Pstream::nProcs());
+    forAll(compactMap, proci)
+    {
+        if (proci == Pstream::myProcNo())
+        {
+            // All my own elements are used
+            label nLocal = globalNumbering.localSize();
+            wantedRemoteElements[proci] = identity(nLocal);
+            constructMap_[proci] = identity(nLocal);
+        }
+        else
+        {
+            // Remote elements wanted from processor proci
+            labelList& remoteElem = wantedRemoteElements[proci];
+            labelList& localElem = constructMap_[proci];
+            remoteElem.setSize(compactMap[proci].size());
+            localElem.setSize(compactMap[proci].size());
+            label i = 0;
+            forAllIter(Map<label>, compactMap[proci], iter)
+            {
+                const label compactI = compactStart[proci] + iter();
+                remoteElem[i] = iter.key();
+                localElem[i]  = compactI;
+                iter() = compactI;
+                i++;
+            }
+        }
+    }
+
+    subMap_.setSize(Pstream::nProcs());
+    Pstream::exchange<labelList, label>
+    (
+        wantedRemoteElements,
+        subMap_,
+        tag,
+        Pstream::worldComm  // TBD
+    );
+
+    // Renumber elements
+    forAll(elements, i)
+    {
+        elements[i] = renumber(globalNumbering, compactMap, elements[i]);
+    }
+}
+
+
+void Foam::distributionMapBase::exchangeAddressing
+(
+    const int tag,
+    const globalIndex& globalNumbering,
+    labelListList& cellCells,
+    List<Map<label>>& compactMap,
+    labelList& compactStart
+)
+{
+    // The overall compact addressing is
+    // - myProcNo data first (uncompacted)
+    // - all other processors consecutively
+
+    compactStart.setSize(Pstream::nProcs());
+    compactStart[Pstream::myProcNo()] = 0;
+    constructSize_ = globalNumbering.localSize();
+    forAll(compactStart, proci)
+    {
+        if (proci != Pstream::myProcNo())
+        {
+            compactStart[proci] = constructSize_;
+            constructSize_ += compactMap[proci].size();
+        }
+    }
+
+
+
+    // Find out what to receive/send in compact addressing.
+
+    // What I want to receive is what others have to send
+    labelListList wantedRemoteElements(Pstream::nProcs());
+    // Compact addressing for received data
+    constructMap_.setSize(Pstream::nProcs());
+    forAll(compactMap, proci)
+    {
+        if (proci == Pstream::myProcNo())
+        {
+            // All my own elements are used
+            label nLocal = globalNumbering.localSize();
+            wantedRemoteElements[proci] = identity(nLocal);
+            constructMap_[proci] = identity(nLocal);
+        }
+        else
+        {
+            // Remote elements wanted from processor proci
+            labelList& remoteElem = wantedRemoteElements[proci];
+            labelList& localElem = constructMap_[proci];
+            remoteElem.setSize(compactMap[proci].size());
+            localElem.setSize(compactMap[proci].size());
+            label i = 0;
+            forAllIter(Map<label>, compactMap[proci], iter)
+            {
+                const label compactI = compactStart[proci] + iter();
+                remoteElem[i] = iter.key();
+                localElem[i]  = compactI;
+                iter() = compactI;
+                i++;
+            }
+        }
+    }
+
+    subMap_.setSize(Pstream::nProcs());
+    Pstream::exchange<labelList, label>
+    (
+        wantedRemoteElements,
+        subMap_,
+        tag,
+        Pstream::worldComm      // TBD
+    );
+
+    // Renumber elements
+    forAll(cellCells, cellI)
+    {
+        labelList& cCells = cellCells[cellI];
+
+        forAll(cCells, i)
+        {
+            cCells[i] = renumber(globalNumbering, compactMap, cCells[i]);
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::distributionMapBase::distributionMapBase()
+:
+    constructSize_(0),
+    subHasFlip_(false),
+    constructHasFlip_(false),
+    schedulePtr_()
+{}
+
+
+Foam::distributionMapBase::distributionMapBase
+(
+    const label constructSize,
+    const labelListList&& subMap,
+    const labelListList&& constructMap,
+    const bool subHasFlip,
+    const bool constructHasFlip
+)
+:
+    constructSize_(constructSize),
+    subMap_(move(subMap)),
+    constructMap_(move(constructMap)),
+    subHasFlip_(subHasFlip),
+    constructHasFlip_(constructHasFlip),
+    schedulePtr_()
+{}
+
+
+Foam::distributionMapBase::distributionMapBase
+(
+    const labelList& sendProcs,
+    const labelList& recvProcs
+)
+:
+    constructSize_(0),
+    subHasFlip_(false),
+    constructHasFlip_(false),
+    schedulePtr_()
+{
+    if (sendProcs.size() != recvProcs.size())
+    {
+        FatalErrorInFunction
+            << "The send and receive data is not the same length. sendProcs:"
+            << sendProcs.size() << " recvProcs:" << recvProcs.size()
+            << abort(FatalError);
+    }
+
+    // Per processor the number of samples we have to send/receive.
+    labelList nSend(Pstream::nProcs(), 0);
+    labelList nRecv(Pstream::nProcs(), 0);
+
+    forAll(sendProcs, sampleI)
+    {
+        label sendProc = sendProcs[sampleI];
+        label recvProc = recvProcs[sampleI];
+
+        // Note that also need to include local communication (both
+        // RecvProc and sendProc on local processor)
+
+        if (Pstream::myProcNo() == sendProc)
+        {
+            // I am the sender. Count destination processor.
+            nSend[recvProc]++;
+        }
+        if (Pstream::myProcNo() == recvProc)
+        {
+            // I am the receiver.
+            nRecv[sendProc]++;
+        }
+    }
+
+    subMap_.setSize(Pstream::nProcs());
+    constructMap_.setSize(Pstream::nProcs());
+    forAll(nSend, proci)
+    {
+        subMap_[proci].setSize(nSend[proci]);
+        constructMap_[proci].setSize(nRecv[proci]);
+    }
+    nSend = 0;
+    nRecv = 0;
+
+    forAll(sendProcs, sampleI)
+    {
+        label sendProc = sendProcs[sampleI];
+        label recvProc = recvProcs[sampleI];
+
+        if (Pstream::myProcNo() == sendProc)
+        {
+            // I am the sender. Store index I need to send.
+            subMap_[recvProc][nSend[recvProc]++] = sampleI;
+        }
+        if (Pstream::myProcNo() == recvProc)
+        {
+            // I am the receiver.
+            constructMap_[sendProc][nRecv[sendProc]++] = sampleI;
+            // Largest entry inside constructMap
+            constructSize_ = sampleI+1;
+        }
+    }
+}
+
+
+Foam::distributionMapBase::distributionMapBase
+(
+    const globalIndex& globalNumbering,
+    labelList& elements,
+    List<Map<label>>& compactMap,
+    const int tag
+)
+:
+    constructSize_(0),
+    subHasFlip_(false),
+    constructHasFlip_(false),
+    schedulePtr_()
+{
+    // Construct per processor compact addressing of the global elements
+    // needed. The ones from the local processor are not included since
+    // these are always all needed.
+    calcCompactAddressing
+    (
+        globalNumbering,
+        elements,
+        compactMap
+    );
+
+    //// Sort remote elements needed (not really necessary)
+    // forAll(compactMap, proci)
+    //{
+    //    if (proci != Pstream::myProcNo())
+    //    {
+    //        Map<label>& globalMap = compactMap[proci];
+    //
+    //        SortableList<label> sorted(move(globalMap.toc()));
+    //
+    //        forAll(sorted, i)
+    //        {
+    //            Map<label>::iterator iter = globalMap.find(sorted[i]);
+    //            iter() = i;
+    //        }
+    //    }
+    //}
+
+
+    // Exchange what I need with processor that supplies it. Renumber elements
+    // into compact numbering
+    labelList compactStart;
+    exchangeAddressing
+    (
+        tag,
+        globalNumbering,
+        elements,
+        compactMap,
+        compactStart
+    );
+
+    if (debug)
+    {
+        printLayout(Pout);
+    }
+}
+
+
+Foam::distributionMapBase::distributionMapBase
+(
+    const globalIndex& globalNumbering,
+    labelListList& cellCells,
+    List<Map<label>>& compactMap,
+    const int tag
+)
+:
+    constructSize_(0),
+    subHasFlip_(false),
+    constructHasFlip_(false),
+    schedulePtr_()
+{
+    // Construct per processor compact addressing of the global elements
+    // needed. The ones from the local processor are not included since
+    // these are always all needed.
+    calcCompactAddressing
+    (
+        globalNumbering,
+        cellCells,
+        compactMap
+    );
+
+    //// Sort remote elements needed (not really necessary)
+    // forAll(compactMap, proci)
+    //{
+    //    if (proci != Pstream::myProcNo())
+    //    {
+    //        Map<label>& globalMap = compactMap[proci];
+    //
+    //        SortableList<label> sorted(move(globalMap.toc()));
+    //
+    //        forAll(sorted, i)
+    //        {
+    //            Map<label>::iterator iter = globalMap.find(sorted[i]);
+    //            iter() = i;
+    //        }
+    //    }
+    //}
+
+
+    // Exchange what I need with processor that supplies it. Renumber elements
+    // into compact numbering
+    labelList compactStart;
+    exchangeAddressing
+    (
+        tag,
+        globalNumbering,
+        cellCells,
+        compactMap,
+        compactStart
+    );
+
+    if (debug)
+    {
+        printLayout(Pout);
+    }
+}
+
+
+Foam::distributionMapBase::distributionMapBase(const distributionMapBase& map)
+:
+    constructSize_(map.constructSize_),
+    subMap_(map.subMap_),
+    constructMap_(map.constructMap_),
+    subHasFlip_(map.subHasFlip_),
+    constructHasFlip_(map.constructHasFlip_),
+    schedulePtr_()
+{}
+
+
+Foam::distributionMapBase::distributionMapBase(distributionMapBase&& map)
+:
+    constructSize_(map.constructSize_),
+    subMap_(move(map.subMap_)),
+    constructMap_(move(map.constructMap_)),
+    subHasFlip_(map.subHasFlip_),
+    constructHasFlip_(map.constructHasFlip_),
+    schedulePtr_()
+{}
+
+
+Foam::distributionMapBase::distributionMapBase(Istream& is)
+{
+    is >> *this;
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::distributionMapBase::transfer(distributionMapBase& rhs)
+{
+    constructSize_ = rhs.constructSize_;
+    subMap_.transfer(rhs.subMap_);
+    constructMap_.transfer(rhs.constructMap_);
+    subHasFlip_ = rhs.subHasFlip_;
+    constructHasFlip_ = rhs.constructHasFlip_;
+    schedulePtr_.clear();
+}
+
+
+Foam::label Foam::distributionMapBase::renumber
+(
+    const globalIndex& globalNumbering,
+    const List<Map<label>>& compactMap,
+    const label globalI
+)
+{
+    if (globalI == -1)
+    {
+        return globalI;
+    }
+    if (globalNumbering.isLocal(globalI))
+    {
+        return globalNumbering.toLocal(globalI);
+    }
+    else
+    {
+        label proci = globalNumbering.whichProcID(globalI);
+        label index = globalNumbering.toLocal(proci, globalI);
+        return compactMap[proci][index];
+    }
+}
+
+
+void Foam::distributionMapBase::compact
+(
+    const boolList& elemIsUsed,
+    const int tag
+)
+{
+    // 1. send back to sender. Have sender delete the corresponding element
+    //    from the submap and do the same to the constructMap locally
+    //    (and in same order).
+
+    // Send elemIsUsed field to neighbour. Use nonblocking code from
+    // distributionMapBase but in reverse order.
+    if (Pstream::parRun())
+    {
+        label startOfRequests = Pstream::nRequests();
+
+        // Set up receives from neighbours
+
+        List<boolList> recvFields(Pstream::nProcs());
+
+        for (label domain = 0; domain < Pstream::nProcs(); domain++)
+        {
+            const labelList& map = subMap_[domain];
+
+            if (domain != Pstream::myProcNo() && map.size())
+            {
+                recvFields[domain].setSize(map.size());
+                IPstream::read
+                (
+                    Pstream::commsTypes::nonBlocking,
+                    domain,
+                    reinterpret_cast<char*>(recvFields[domain].begin()),
+                    recvFields[domain].size()*sizeof(bool),
+                    tag
+                );
+            }
+        }
+
+
+        List<boolList> sendFields(Pstream::nProcs());
+
+        for (label domain = 0; domain < Pstream::nProcs(); domain++)
+        {
+            const labelList& map = constructMap_[domain];
+
+            if (domain != Pstream::myProcNo() && map.size())
+            {
+                boolList& subField = sendFields[domain];
+                subField.setSize(map.size());
+                forAll(map, i)
+                {
+                    subField[i] = accessAndFlip
+                    (
+                        elemIsUsed,
+                        map[i],
+                        constructHasFlip_,
+                        noOp()          // do not flip elemIsUsed value
+                    );
+                }
+
+                OPstream::write
+                (
+                    Pstream::commsTypes::nonBlocking,
+                    domain,
+                    reinterpret_cast<const char*>(subField.begin()),
+                    subField.size()*sizeof(bool),
+                    tag
+                );
+            }
+        }
+
+
+
+        // Set up 'send' to myself - write directly into recvFields
+
+        {
+            const labelList& map = constructMap_[Pstream::myProcNo()];
+
+            recvFields[Pstream::myProcNo()].setSize(map.size());
+            forAll(map, i)
+            {
+                recvFields[Pstream::myProcNo()][i] = accessAndFlip
+                (
+                    elemIsUsed,
+                    map[i],
+                    constructHasFlip_,
+                    noOp()              // do not flip elemIsUsed value
+                );
+            }
+        }
+
+
+        // Wait for all to finish
+
+        Pstream::waitRequests(startOfRequests);
+
+
+        // Compact out all submap entries that are referring to unused elements
+        for (label domain = 0; domain < Pstream::nProcs(); domain++)
+        {
+            const labelList& map = subMap_[domain];
+
+            labelList newMap(map.size());
+            label newI = 0;
+
+            forAll(map, i)
+            {
+                if (recvFields[domain][i])
+                {
+                    // So element is used on destination side
+                    newMap[newI++] = map[i];
+                }
+            }
+            if (newI < map.size())
+            {
+                newMap.setSize(newI);
+                subMap_[domain].transfer(newMap);
+            }
+        }
+    }
+
+
+    // 2. remove from construct map - since end-result (element in elemIsUsed)
+    //    not used.
+
+    label maxConstructIndex = -1;
+
+    for (label domain = 0; domain < Pstream::nProcs(); domain++)
+    {
+        const labelList& map = constructMap_[domain];
+
+        labelList newMap(map.size());
+        label newI = 0;
+
+        forAll(map, i)
+        {
+            label destinationI = map[i];
+            if (constructHasFlip_)
+            {
+                destinationI = mag(destinationI)-1;
+            }
+
+            // Is element is used on destination side
+            if (elemIsUsed[destinationI])
+            {
+                maxConstructIndex = max(maxConstructIndex, destinationI);
+
+                newMap[newI++] = map[i];
+            }
+        }
+        if (newI < map.size())
+        {
+            newMap.setSize(newI);
+            constructMap_[domain].transfer(newMap);
+        }
+    }
+
+    constructSize_ = maxConstructIndex+1;
+
+    // Clear the schedule (note:not necessary if nothing changed)
+    schedulePtr_.clear();
+}
+
+
+void Foam::distributionMapBase::compact
+(
+    const boolList& elemIsUsed,
+    const label localSize,            // max index for subMap
+    labelList& oldToNewSub,
+    labelList& oldToNewConstruct,
+    const int tag
+)
+{
+    // 1. send back to sender. Have sender delete the corresponding element
+    //    from the submap and do the same to the constructMap locally
+    //    (and in same order).
+
+    // Send elemIsUsed field to neighbour. Use nonblocking code from
+    // distributionMapBase but in reverse order.
+    if (Pstream::parRun())
+    {
+        label startOfRequests = Pstream::nRequests();
+
+        // Set up receives from neighbours
+
+        List<boolList> recvFields(Pstream::nProcs());
+
+        for (label domain = 0; domain < Pstream::nProcs(); domain++)
+        {
+            const labelList& map = subMap_[domain];
+
+            if (domain != Pstream::myProcNo() && map.size())
+            {
+                recvFields[domain].setSize(map.size());
+                IPstream::read
+                (
+                    Pstream::commsTypes::nonBlocking,
+                    domain,
+                    reinterpret_cast<char*>(recvFields[domain].begin()),
+                    recvFields[domain].size()*sizeof(bool),
+                    tag
+                );
+            }
+        }
+
+
+        List<boolList> sendFields(Pstream::nProcs());
+
+        for (label domain = 0; domain < Pstream::nProcs(); domain++)
+        {
+            const labelList& map = constructMap_[domain];
+
+            if (domain != Pstream::myProcNo() && map.size())
+            {
+                boolList& subField = sendFields[domain];
+                subField.setSize(map.size());
+                forAll(map, i)
+                {
+                    label index = map[i];
+                    if (constructHasFlip_)
+                    {
+                        index = mag(index)-1;
+                    }
+                    subField[i] = elemIsUsed[index];
+                }
+
+                OPstream::write
+                (
+                    Pstream::commsTypes::nonBlocking,
+                    domain,
+                    reinterpret_cast<const char*>(subField.begin()),
+                    subField.size()*sizeof(bool),
+                    tag
+                );
+            }
+        }
+
+
+
+        // Set up 'send' to myself - write directly into recvFields
+
+        {
+            const labelList& map = constructMap_[Pstream::myProcNo()];
+
+            recvFields[Pstream::myProcNo()].setSize(map.size());
+            forAll(map, i)
+            {
+                label index = map[i];
+                if (constructHasFlip_)
+                {
+                    index = mag(index)-1;
+                }
+                recvFields[Pstream::myProcNo()][i] = elemIsUsed[index];
+            }
+        }
+
+
+        // Wait for all to finish
+
+        Pstream::waitRequests(startOfRequests);
+
+
+
+
+        // Work out which elements on the sending side are needed
+        {
+            oldToNewSub.setSize(localSize, -1);
+
+            boolList sendElemIsUsed(localSize, false);
+
+            for (label domain = 0; domain < Pstream::nProcs(); domain++)
+            {
+                const labelList& map = subMap_[domain];
+                forAll(map, i)
+                {
+                    if (recvFields[domain][i])
+                    {
+                        label index = map[i];
+                        if (subHasFlip_)
+                        {
+                            index = mag(index)-1;
+                        }
+                        sendElemIsUsed[index] = true;
+                    }
+                }
+            }
+
+            label newI = 0;
+            forAll(sendElemIsUsed, i)
+            {
+                if (sendElemIsUsed[i])
+                {
+                    oldToNewSub[i] = newI++;
+                }
+            }
+        }
+
+
+        // Compact out all submap entries that are referring to unused elements
+        for (label domain = 0; domain < Pstream::nProcs(); domain++)
+        {
+            const labelList& map = subMap_[domain];
+
+            labelList newMap(map.size());
+            label newI = 0;
+
+            forAll(map, i)
+            {
+                if (recvFields[domain][i])
+                {
+                    // So element is used on destination side
+                    label index = map[i];
+                    label sign = 1;
+                    if (subHasFlip_)
+                    {
+                        if (index < 0)
+                        {
+                            sign = -1;
+                        }
+                        index = mag(index)-1;
+                    }
+                    label newIndex = oldToNewSub[index];
+                    if (subHasFlip_)
+                    {
+                        newIndex = sign*(newIndex+1);
+                    }
+                    newMap[newI++] = newIndex;
+                }
+            }
+            newMap.setSize(newI);
+            subMap_[domain].transfer(newMap);
+        }
+    }
+
+
+    // 2. remove from construct map - since end-result (element in elemIsUsed)
+    //    not used.
+
+
+    oldToNewConstruct.setSize(elemIsUsed.size(), -1);
+    constructSize_ = 0;
+    forAll(elemIsUsed, i)
+    {
+        if (elemIsUsed[i])
+        {
+            oldToNewConstruct[i] = constructSize_++;
+        }
+    }
+
+    for (label domain = 0; domain < Pstream::nProcs(); domain++)
+    {
+        const labelList& map = constructMap_[domain];
+
+        labelList newMap(map.size());
+        label newI = 0;
+
+        forAll(map, i)
+        {
+            label destinationI = map[i];
+            label sign = 1;
+            if (constructHasFlip_)
+            {
+                if (destinationI < 0)
+                {
+                    sign = -1;
+                }
+                destinationI = mag(destinationI)-1;
+            }
+
+            // Is element is used on destination side
+            if (elemIsUsed[destinationI])
+            {
+                label newIndex = oldToNewConstruct[destinationI];
+                if (constructHasFlip_)
+                {
+                    newIndex = sign*(newIndex+1);
+                }
+                newMap[newI++] = newIndex;
+            }
+        }
+        newMap.setSize(newI);
+        constructMap_[domain].transfer(newMap);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+void Foam::distributionMapBase::operator=(const distributionMapBase& rhs)
+{
+    // Check for assignment to self
+    if (this == &rhs)
+    {
+        FatalErrorInFunction
+            << "Attempted assignment to self"
+            << abort(FatalError);
+    }
+    constructSize_ = rhs.constructSize_;
+    subMap_ = rhs.subMap_;
+    constructMap_ = rhs.constructMap_;
+    subHasFlip_ = rhs.subHasFlip_;
+    constructHasFlip_ = rhs.constructHasFlip_;
+    schedulePtr_.clear();
+}
+
+
+// * * * * * * * * * * * * * * Istream Operator  * * * * * * * * * * * * * * //
+
+Foam::Istream& Foam::operator>>(Istream& is, distributionMapBase& map)
+{
+    is.fatalCheck("operator>>(Istream&, distributionMapBase&)");
+
+    is  >> map.constructSize_ >> map.subMap_ >> map.constructMap_
+        >> map.subHasFlip_ >> map.constructHasFlip_;
+
+    return is;
+}
+
+
+// * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const distributionMapBase& map)
+{
+    os  << map.constructSize_ << token::NL
+        << map.subMap_ << token::NL
+        << map.constructMap_ << token::NL
+        << map.subHasFlip_ << token::SPACE << map.constructHasFlip_
+        << token::NL;
+
+    return os;
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/distributionMapBase.H b/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/distributionMapBase.H
new file mode 100644
index 0000000000..95ea2af1bf
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/distributionMapBase.H
@@ -0,0 +1,495 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2015-2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::distributionMapBase
+
+Description
+    Class containing processor-to-processor mapping information.
+
+    We store mapping from the bits-to-send to the complete starting list
+    (subXXXMap) and from the received bits to their location in the new
+    list (constructXXXMap).
+
+Note:
+    Schedule is a list of processor pairs (one send, one receive. One of
+    them will be myself) which forms a scheduled (i.e. non-buffered) exchange.
+    See distribute on how to use it.
+    Note2: number of items sent on one processor have to equal the number
+    of items received on the other processor.
+
+    To aid constructing these maps there are the constructors from global
+    numbering, either with or without transforms.
+
+    Constructors using compact numbering: layout is
+    - all my own elements first (whether used or not)
+    - followed by used-only remote elements sorted by remote processor.
+    So e.g 4 procs and on proc 1 the compact
+    table will first have all globalIndex.localSize() elements from proc1
+    followed by used-only elements of proc0, proc2, proc3.
+    The constructed distributionMapBase sends the local elements from and
+    receives the remote elements into their compact position.
+    compactMap[proci] is the position of elements from proci in the compact
+    map. compactMap[myProcNo()] is empty since trivial addressing.
+
+    It rewrites the input global indices into indices into the constructed
+    data.
+
+    When constructing from components optionally a 'flip' on
+    the maps can be specified. This will interpret the map
+    values as index+flip, similar to e.g. faceProcAddressing. The flip
+    will only be applied to fieldTypes (scalar, vector, .. triad)
+
+
+SourceFiles
+    distributionMapBase.C
+    distributionMapBaseTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef distributionMapBase_H
+#define distributionMapBase_H
+
+#include "labelList.H"
+#include "labelPair.H"
+#include "Pstream.H"
+#include "boolList.H"
+#include "Map.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class polyTopoChangeMap;
+class globalIndex;
+class PstreamBuffers;
+
+
+// Forward declaration of friend functions and operators
+
+class distributionMapBase;
+
+Istream& operator>>(Istream&, distributionMapBase&);
+Ostream& operator<<(Ostream&, const distributionMapBase&);
+
+
+/*---------------------------------------------------------------------------*\
+                      Class distributionMapBase Declaration
+\*---------------------------------------------------------------------------*/
+
+class distributionMapBase
+{
+protected:
+
+    // Protected data
+
+        //- Size of reconstructed data
+        label constructSize_;
+
+        //- Maps from subsetted data back to original data
+        labelListList subMap_;
+
+        //- Maps from subsetted data to new reconstructed data
+        labelListList constructMap_;
+
+        //- Whether subMap includes flip or not
+        bool subHasFlip_;
+
+        //- Whether constructMap includes flip or not
+        bool constructHasFlip_;
+
+
+        //- Schedule
+        mutable autoPtr<List<labelPair>> schedulePtr_;
+
+
+    // Private Member Functions
+
+        static void checkReceivedSize
+        (
+            const label proci,
+            const label expectedSize,
+            const label receivedSize
+        );
+
+        //- Construct per processor compact addressing of the global elements
+        //  needed. The ones from the local processor are not included since
+        //  these are always all needed.
+        void calcCompactAddressing
+        (
+            const globalIndex& globalNumbering,
+            const labelList& elements,
+            List<Map<label>>& compactMap
+        ) const;
+
+        void calcCompactAddressing
+        (
+            const globalIndex& globalNumbering,
+            const labelListList& elements,
+            List<Map<label>>& compactMap
+        ) const;
+
+        void exchangeAddressing
+        (
+            const int tag,
+            const globalIndex& globalNumbering,
+            labelList& elements,
+            List<Map<label>>& compactMap,
+            labelList& compactStart
+        );
+        void exchangeAddressing
+        (
+            const int tag,
+            const globalIndex& globalNumbering,
+            labelListList& elements,
+            List<Map<label>>& compactMap,
+            labelList& compactStart
+        );
+
+        template<class T, class CombineOp, class negateOp>
+        static void flipAndCombine
+        (
+            const UList<label>& map,
+            const bool hasFlip,
+            const UList<T>& rhs,
+            const CombineOp& cop,
+            const negateOp& negOp,
+            List<T>& lhs
+        );
+
+        template<class T, class negateOp>
+        static T accessAndFlip
+        (
+            const UList<T>& fld,
+            const label index,
+            const bool hasFlip,
+            const negateOp& negOp
+        );
+
+public:
+
+    // Declare name of the class and its debug switch
+    ClassName("distributionMapBase");
+
+
+    // Constructors
+
+        //- Construct null
+        distributionMapBase();
+
+        //- Move construct from components
+        distributionMapBase
+        (
+            const label constructSize,
+            const labelListList&& subMap,
+            const labelListList&& constructMap,
+            const bool subHasFlip = false,
+            const bool constructHasFlip = false
+        );
+
+        //- Construct from reverse addressing: per data item the send
+        //  processor and the receive processor. (note: data is not stored
+        //  sorted per processor so cannot use printLayout).
+        distributionMapBase
+        (
+            const labelList& sendProcs,
+            const labelList& recvProcs
+        );
+
+        //- Construct from list of (possibly) remote elements in globalIndex
+        //  numbering (or -1). Determines compact numbering (see above) and
+        //  distribute map to get data into this ordering and renumbers the
+        //  elements to be in compact numbering.
+        distributionMapBase
+        (
+            const globalIndex&,
+            labelList& elements,
+            List<Map<label>>& compactMap,
+            const int tag = Pstream::msgType()
+        );
+
+        //- Special variant that works with the info sorted into bins
+        //  according to local indices. E.g. think cellCells where
+        //  cellCells[localCellI] is a list of global cells
+        distributionMapBase
+        (
+            const globalIndex&,
+            labelListList& cellCells,
+            List<Map<label>>& compactMap,
+            const int tag = Pstream::msgType()
+        );
+
+        //- Construct copy
+        distributionMapBase(const distributionMapBase&);
+
+        //- Move constructor
+        distributionMapBase(distributionMapBase&&);
+
+        //- Construct from Istream
+        distributionMapBase(Istream&);
+
+
+    // Member Functions
+
+        // Access
+
+            //- Constructed data size
+            label constructSize() const
+            {
+                return constructSize_;
+            }
+
+            //- Constructed data size
+            label& constructSize()
+            {
+                return constructSize_;
+            }
+
+            //- From subsetted data back to original data
+            const labelListList& subMap() const
+            {
+                return subMap_;
+            }
+
+            //- From subsetted data back to original data
+            labelListList& subMap()
+            {
+                return subMap_;
+            }
+
+            //- From subsetted data to new reconstructed data
+            const labelListList& constructMap() const
+            {
+                return constructMap_;
+            }
+
+            //- From subsetted data to new reconstructed data
+            labelListList& constructMap()
+            {
+                return constructMap_;
+            }
+
+            //- Does subMap include a sign
+            bool subHasFlip() const
+            {
+                return subHasFlip_;
+            }
+
+            //- Does subMap include a sign
+            bool& subHasFlip()
+            {
+                return subHasFlip_;
+            }
+
+            //- Does constructMap include a sign
+            bool constructHasFlip() const
+            {
+                return constructHasFlip_;
+            }
+
+            //- Does constructMap include a sign
+            bool& constructHasFlip()
+            {
+                return constructHasFlip_;
+            }
+
+            //- Calculate a schedule. See above.
+            static List<labelPair> schedule
+            (
+                const labelListList& subMap,
+                const labelListList& constructMap,
+                const int tag
+            );
+
+            //- Return a schedule. Demand driven. See above.
+            const List<labelPair>& schedule() const;
+
+
+        // Other
+
+            //- Transfer the contents of the argument and annul the argument.
+            void transfer(distributionMapBase&);
+
+            //- Helper for construct from globalIndex. Renumbers element
+            //  (in globalIndex numbering) into compact indices.
+            static label renumber
+            (
+                const globalIndex&,
+                const List<Map<label>>& compactMap,
+                const label globalElement
+            );
+
+            //- Compact maps. Gets per field a bool whether it is used (locally)
+            //  and works out itself what this side and sender side can remove
+            //  from maps. Only compacts non-local elements (i.e. the stuff
+            //  that gets sent over), does not change the local layout
+            void compact
+            (
+                const boolList& elemIsUsed,
+                const int tag = UPstream::msgType()
+            );
+
+            //- Compact all maps and layout. Returns compaction maps for
+            //  subMap and constructMap
+            void compact
+            (
+                const boolList& elemIsUsed,
+                const label localSize,            // max index for subMap
+                labelList& oldToNewSub,
+                labelList& oldToNewConstruct,
+                const int tag = UPstream::msgType()
+            );
+
+            //- Distribute data. Note:schedule only used for
+            //  Pstream::commsTypes::scheduled for now, all others just use
+            //  send-to-all, receive-from-all.
+            template<class T, class negateOp>
+            static void distribute
+            (
+                const Pstream::commsTypes commsType,
+                const List<labelPair>& schedule,
+                const label constructSize,
+                const labelListList& subMap,
+                const bool subHasFlip,
+                const labelListList& constructMap,
+                const bool constructHasFlip,
+                List<T>&,
+                const negateOp& negOp,
+                const int tag = UPstream::msgType()
+            );
+
+            //- Distribute data. If multiple processors writing to same
+            //  position adds contributions using cop.
+            template<class T, class CombineOp, class negateOp>
+            static void distribute
+            (
+                const Pstream::commsTypes commsType,
+                const List<labelPair>& schedule,
+                const label constructSize,
+                const labelListList& subMap,
+                const bool subHasFlip,
+                const labelListList& constructMap,
+                const bool constructHasFlip,
+                List<T>&,
+                const CombineOp& cop,
+                const negateOp& negOp,
+                const T& nullValue,
+                const int tag = UPstream::msgType()
+            );
+
+            //- Distribute data using default commsType.
+            template<class T>
+            void distribute
+            (
+                List<T>& fld,
+                const int tag = UPstream::msgType()
+            ) const;
+
+            //- Distribute data using default commsType.
+            template<class T, class negateOp>
+            void distribute
+            (
+                List<T>& fld,
+                const negateOp& negOp,
+                const int tag = UPstream::msgType()
+            ) const;
+
+            //- Distribute data using default commsType.
+            template<class T>
+            void distribute
+            (
+                DynamicList<T>& fld,
+                const int tag = UPstream::msgType()
+            ) const;
+
+            //- Reverse distribute data using default commsType.
+            template<class T>
+            void reverseDistribute
+            (
+                const label constructSize,
+                List<T>&,
+                const int tag = UPstream::msgType()
+            ) const;
+
+            //- Reverse distribute data using default commsType.
+            //  Since constructSize might be larger than supplied size supply
+            //  a nullValue
+            template<class T>
+            void reverseDistribute
+            (
+                const label constructSize,
+                const T& nullValue,
+                List<T>& fld,
+                const int tag = UPstream::msgType()
+            ) const;
+
+            //- Do all sends using PstreamBuffers
+            template<class T>
+            void send(PstreamBuffers&, const List<T>&) const;
+            //- Do all receives using PstreamBuffers
+            template<class T>
+            void receive(PstreamBuffers&, List<T>&) const;
+
+            //- Debug: print layout. Can only be used on maps with sorted
+            //  storage (local data first, then non-local data)
+            void printLayout(Ostream& os) const;
+
+            //- Correct for topo change.
+            void updateMesh(const polyTopoChangeMap&)
+            {
+                NotImplemented;
+            }
+
+
+    // Member Operators
+
+        void operator=(const distributionMapBase&);
+        void operator=(distributionMapBase&&);
+
+
+    // IOstream Operators
+
+        //- Read dictionary from Istream
+        friend Istream& operator>>(Istream&, distributionMapBase&);
+
+        //- Write dictionary to Ostream
+        friend Ostream& operator<<(Ostream&, const distributionMapBase&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "distributionMapBaseTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/distributionMapBaseTemplates.C b/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/distributionMapBaseTemplates.C
new file mode 100644
index 0000000000..a0c263079c
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/distributionMapBaseTemplates.C
@@ -0,0 +1,1420 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2015-2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "Pstream.H"
+#include "PstreamBuffers.H"
+#include "PstreamCombineReduceOps.H"
+#include "flipOp.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class T, class CombineOp, class negateOp>
+void Foam::distributionMapBase::flipAndCombine
+(
+    const UList<label>& map,
+    const bool hasFlip,
+    const UList<T>& rhs,
+    const CombineOp& cop,
+    const negateOp& negOp,
+    List<T>& lhs
+)
+{
+    if (hasFlip)
+    {
+        forAll(map, i)
+        {
+            if (map[i] > 0)
+            {
+                label index = map[i]-1;
+                cop(lhs[index], rhs[i]);
+            }
+            else if (map[i] < 0)
+            {
+                label index = -map[i]-1;
+                cop(lhs[index], negOp(rhs[i]));
+            }
+            else
+            {
+                FatalErrorInFunction
+                    << "At index " << i << " out of " << map.size()
+                    << " have illegal index " << map[i]
+                    << " for field " << rhs.size() << " with flipMap"
+                    << exit(FatalError);
+            }
+        }
+    }
+    else
+    {
+        forAll(map, i)
+        {
+            cop(lhs[map[i]], rhs[i]);
+        }
+    }
+}
+
+
+template<class T, class negateOp>
+T Foam::distributionMapBase::accessAndFlip
+(
+    const UList<T>& fld,
+    const label index,
+    const bool hasFlip,
+    const negateOp& negOp
+)
+{
+    T t;
+    if (hasFlip)
+    {
+        if (index > 0)
+        {
+            t = fld[index-1];
+        }
+        else if (index < 0)
+        {
+            t = negOp(fld[-index-1]);
+        }
+        else
+        {
+            FatalErrorInFunction
+                << "Illegal index " << index
+                << " into field of size " << fld.size()
+                << " with face-flipping"
+                << exit(FatalError);
+            t = fld[index];
+        }
+    }
+    else
+    {
+        t = fld[index];
+    }
+    return t;
+}
+
+
+// Distribute list.
+template<class T, class negateOp>
+void Foam::distributionMapBase::distribute
+(
+    const Pstream::commsTypes commsType,
+    const List<labelPair>& schedule,
+    const label constructSize,
+    const labelListList& subMap,
+    const bool subHasFlip,
+    const labelListList& constructMap,
+    const bool constructHasFlip,
+    List<T>& field,
+    const negateOp& negOp,
+    const int tag
+)
+{
+    if (!Pstream::parRun())
+    {
+        // Do only me to me.
+
+        const labelList& mySubMap = subMap[Pstream::myProcNo()];
+
+        List<T> subField(mySubMap.size());
+        forAll(mySubMap, i)
+        {
+            subField[i] = accessAndFlip(field, mySubMap[i], subHasFlip, negOp);
+        }
+
+        // Receive sub field from myself (subField)
+        const labelList& map = constructMap[Pstream::myProcNo()];
+
+        field.setSize(constructSize);
+
+        flipAndCombine
+        (
+            map,
+            constructHasFlip,
+            subField,
+            eqOp<T>(),
+            negOp,
+            field
+        );
+
+        return;
+    }
+
+    if (commsType == Pstream::commsTypes::blocking)
+    {
+        // Since buffered sending can reuse the field to collect the
+        // received data.
+
+        // Send sub field to neighbour
+        for (label domain = 0; domain < Pstream::nProcs(); domain++)
+        {
+            const labelList& map = subMap[domain];
+
+            if (domain != Pstream::myProcNo() && map.size())
+            {
+                OPstream toNbr(Pstream::commsTypes::blocking, domain, 0, tag);
+
+                List<T> subField(map.size());
+                forAll(subField, i)
+                {
+                    subField[i] = accessAndFlip
+                    (
+                        field,
+                        map[i],
+                        subHasFlip,
+                        negOp
+                    );
+                }
+                toNbr << subField;
+            }
+        }
+
+        // Subset myself
+        const labelList& mySubMap = subMap[Pstream::myProcNo()];
+
+        List<T> subField(mySubMap.size());
+        forAll(mySubMap, i)
+        {
+            subField[i] = accessAndFlip(field, mySubMap[i], subHasFlip, negOp);
+        }
+
+        // Receive sub field from myself (subField)
+        const labelList& map = constructMap[Pstream::myProcNo()];
+
+        field.setSize(constructSize);
+
+        flipAndCombine
+        (
+            map,
+            constructHasFlip,
+            subField,
+            eqOp<T>(),
+            negOp,
+            field
+        );
+
+        // Receive sub field from neighbour
+        for (label domain = 0; domain < Pstream::nProcs(); domain++)
+        {
+            const labelList& map = constructMap[domain];
+
+            if (domain != Pstream::myProcNo() && map.size())
+            {
+                IPstream fromNbr(Pstream::commsTypes::blocking, domain, 0, tag);
+                List<T> subField(fromNbr);
+
+                checkReceivedSize(domain, map.size(), subField.size());
+
+                flipAndCombine
+                (
+                    map,
+                    constructHasFlip,
+                    subField,
+                    eqOp<T>(),
+                    negOp,
+                    field
+                );
+            }
+        }
+    }
+    else if (commsType == Pstream::commsTypes::scheduled)
+    {
+        // Need to make sure I don't overwrite field with received data
+        // since the data might need to be sent to another processor. So
+        // allocate a new field for the results.
+        List<T> newField(constructSize);
+
+        // Receive sub field from myself
+        {
+            const labelList& mySubMap = subMap[Pstream::myProcNo()];
+
+            List<T> subField(mySubMap.size());
+            forAll(subField, i)
+            {
+                subField[i] = accessAndFlip
+                (
+                    field,
+                    mySubMap[i],
+                    subHasFlip,
+                    negOp
+                );
+            }
+
+            // Receive sub field from myself (subField)
+            flipAndCombine
+            (
+                constructMap[Pstream::myProcNo()],
+                constructHasFlip,
+                subField,
+                eqOp<T>(),
+                negOp,
+                newField
+            );
+        }
+
+        // Schedule will already have pruned 0-sized comms
+        forAll(schedule, i)
+        {
+            const labelPair& twoProcs = schedule[i];
+            // twoProcs is a swap pair of processors. The first one is the
+            // one that needs to send first and then receive.
+
+            label sendProc = twoProcs[0];
+            label recvProc = twoProcs[1];
+
+            if (Pstream::myProcNo() == sendProc)
+            {
+                // I am send first, receive next
+                {
+                    OPstream toNbr
+                    (
+                        Pstream::commsTypes::scheduled,
+                        recvProc,
+                        0,
+                        tag
+                    );
+
+                    const labelList& map = subMap[recvProc];
+                    List<T> subField(map.size());
+                    forAll(subField, i)
+                    {
+                        subField[i] = accessAndFlip
+                        (
+                            field,
+                            map[i],
+                            subHasFlip,
+                            negOp
+                        );
+                    }
+                    toNbr << subField;
+                }
+                {
+                    IPstream fromNbr
+                    (
+                        Pstream::commsTypes::scheduled,
+                        recvProc,
+                        0,
+                        tag
+                    );
+                    List<T> subField(fromNbr);
+
+                    const labelList& map = constructMap[recvProc];
+
+                    checkReceivedSize(recvProc, map.size(), subField.size());
+
+                    flipAndCombine
+                    (
+                        map,
+                        constructHasFlip,
+                        subField,
+                        eqOp<T>(),
+                        negOp,
+                        newField
+                    );
+                }
+            }
+            else
+            {
+                // I am receive first, send next
+                {
+                    IPstream fromNbr
+                    (
+                        Pstream::commsTypes::scheduled,
+                        sendProc,
+                        0,
+                        tag
+                    );
+                    List<T> subField(fromNbr);
+
+                    const labelList& map = constructMap[sendProc];
+
+                    checkReceivedSize(sendProc, map.size(), subField.size());
+
+                    flipAndCombine
+                    (
+                        map,
+                        constructHasFlip,
+                        subField,
+                        eqOp<T>(),
+                        negOp,
+                        newField
+                    );
+                }
+                {
+                    OPstream toNbr
+                    (
+                        Pstream::commsTypes::scheduled,
+                        sendProc,
+                        0,
+                        tag
+                    );
+
+                    const labelList& map = subMap[sendProc];
+                    List<T> subField(map.size());
+                    forAll(subField, i)
+                    {
+                        subField[i] = accessAndFlip
+                        (
+                            field,
+                            map[i],
+                            subHasFlip,
+                            negOp
+                        );
+                    }
+                    toNbr << subField;
+                }
+            }
+        }
+        field.transfer(newField);
+    }
+    else if (commsType == Pstream::commsTypes::nonBlocking)
+    {
+        label nOutstanding = Pstream::nRequests();
+
+        if (!contiguous<T>())
+        {
+            PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking, tag);
+
+            // Stream data into buffer
+            for (label domain = 0; domain < Pstream::nProcs(); domain++)
+            {
+                const labelList& map = subMap[domain];
+
+                if (domain != Pstream::myProcNo() && map.size())
+                {
+                    // Put data into send buffer
+                    UOPstream toDomain(domain, pBufs);
+
+                    List<T> subField(map.size());
+                    forAll(subField, i)
+                    {
+                        subField[i] = accessAndFlip
+                        (
+                            field,
+                            map[i],
+                            subHasFlip,
+                            negOp
+                        );
+                    }
+                    toDomain << subField;
+                }
+            }
+
+            // Start receiving. Do not block.
+            pBufs.finishedSends(false);
+
+            {
+                // Set up 'send' to myself
+                const labelList& mySub = subMap[Pstream::myProcNo()];
+                List<T> mySubField(mySub.size());
+                forAll(mySub, i)
+                {
+                    mySubField[i] = accessAndFlip
+                    (
+                        field,
+                        mySub[i],
+                        subHasFlip,
+                        negOp
+                    );
+                }
+                // Combine bits. Note that can reuse field storage
+                field.setSize(constructSize);
+                // Receive sub field from myself
+                {
+                    const labelList& map = constructMap[Pstream::myProcNo()];
+
+                    flipAndCombine
+                    (
+                        map,
+                        constructHasFlip,
+                        mySubField,
+                        eqOp<T>(),
+                        negOp,
+                        field
+                    );
+                }
+            }
+
+            // Block ourselves, waiting only for the current comms
+            Pstream::waitRequests(nOutstanding);
+
+            // Consume
+            for (label domain = 0; domain < Pstream::nProcs(); domain++)
+            {
+                const labelList& map = constructMap[domain];
+
+                if (domain != Pstream::myProcNo() && map.size())
+                {
+                    UIPstream str(domain, pBufs);
+                    List<T> recvField(str);
+
+                    checkReceivedSize(domain, map.size(), recvField.size());
+
+                    flipAndCombine
+                    (
+                        map,
+                        constructHasFlip,
+                        recvField,
+                        eqOp<T>(),
+                        negOp,
+                        field
+                    );
+                }
+            }
+        }
+        else
+        {
+            // Set up sends to neighbours
+
+            List<List<T>> sendFields(Pstream::nProcs());
+
+            for (label domain = 0; domain < Pstream::nProcs(); domain++)
+            {
+                const labelList& map = subMap[domain];
+
+                if (domain != Pstream::myProcNo() && map.size())
+                {
+                    List<T>& subField = sendFields[domain];
+                    subField.setSize(map.size());
+                    forAll(map, i)
+                    {
+                        subField[i] = accessAndFlip
+                        (
+                            field,
+                            map[i],
+                            subHasFlip,
+                            negOp
+                        );
+                    }
+
+                    OPstream::write
+                    (
+                        Pstream::commsTypes::nonBlocking,
+                        domain,
+                        reinterpret_cast<const char*>(subField.begin()),
+                        subField.byteSize(),
+                        tag
+                    );
+                }
+            }
+
+            // Set up receives from neighbours
+
+            List<List<T>> recvFields(Pstream::nProcs());
+
+            for (label domain = 0; domain < Pstream::nProcs(); domain++)
+            {
+                const labelList& map = constructMap[domain];
+
+                if (domain != Pstream::myProcNo() && map.size())
+                {
+                    recvFields[domain].setSize(map.size());
+                    IPstream::read
+                    (
+                        Pstream::commsTypes::nonBlocking,
+                        domain,
+                        reinterpret_cast<char*>(recvFields[domain].begin()),
+                        recvFields[domain].byteSize(),
+                        tag
+                    );
+                }
+            }
+
+
+            // Set up 'send' to myself
+
+            {
+                const labelList& map = subMap[Pstream::myProcNo()];
+
+                List<T>& subField = sendFields[Pstream::myProcNo()];
+                subField.setSize(map.size());
+                forAll(map, i)
+                {
+                    subField[i] = accessAndFlip
+                    (
+                        field,
+                        map[i],
+                        subHasFlip,
+                        negOp
+                    );
+                }
+            }
+
+
+            // Combine bits. Note that can reuse field storage
+
+            field.setSize(constructSize);
+
+
+            // Receive sub field from myself (sendFields[Pstream::myProcNo()])
+            {
+                const labelList& map = constructMap[Pstream::myProcNo()];
+                const List<T>& subField = sendFields[Pstream::myProcNo()];
+
+                flipAndCombine
+                (
+                    map,
+                    constructHasFlip,
+                    subField,
+                    eqOp<T>(),
+                    negOp,
+                    field
+                );
+            }
+
+
+            // Wait for all to finish
+
+            Pstream::waitRequests(nOutstanding);
+
+
+            // Collect neighbour fields
+
+            for (label domain = 0; domain < Pstream::nProcs(); domain++)
+            {
+                const labelList& map = constructMap[domain];
+
+                if (domain != Pstream::myProcNo() && map.size())
+                {
+                    const List<T>& subField = recvFields[domain];
+
+                    checkReceivedSize(domain, map.size(), subField.size());
+
+                    flipAndCombine
+                    (
+                        map,
+                        constructHasFlip,
+                        subField,
+                        eqOp<T>(),
+                        negOp,
+                        field
+                    );
+                }
+            }
+        }
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Unknown communication schedule " << int(commsType)
+            << abort(FatalError);
+    }
+}
+
+
+// Distribute list.
+template<class T, class CombineOp, class negateOp>
+void Foam::distributionMapBase::distribute
+(
+    const Pstream::commsTypes commsType,
+    const List<labelPair>& schedule,
+    const label constructSize,
+    const labelListList& subMap,
+    const bool subHasFlip,
+    const labelListList& constructMap,
+    const bool constructHasFlip,
+    List<T>& field,
+    const CombineOp& cop,
+    const negateOp& negOp,
+    const T& nullValue,
+    const int tag
+)
+{
+    if (!Pstream::parRun())
+    {
+        // Do only me to me.
+
+        const labelList& mySubMap = subMap[Pstream::myProcNo()];
+
+        List<T> subField(mySubMap.size());
+        forAll(mySubMap, i)
+        {
+            subField[i] = accessAndFlip(field, mySubMap[i], subHasFlip, negOp);
+        }
+
+        // Receive sub field from myself (subField)
+        const labelList& map = constructMap[Pstream::myProcNo()];
+
+        field.setSize(constructSize);
+        field = nullValue;
+
+        flipAndCombine(map, constructHasFlip, subField, cop, negOp, field);
+
+        return;
+    }
+
+    if (commsType == Pstream::commsTypes::blocking)
+    {
+        // Since buffered sending can reuse the field to collect the
+        // received data.
+
+        // Send sub field to neighbour
+        for (label domain = 0; domain < Pstream::nProcs(); domain++)
+        {
+            const labelList& map = subMap[domain];
+
+            if (domain != Pstream::myProcNo() && map.size())
+            {
+                OPstream toNbr(Pstream::commsTypes::blocking, domain, 0, tag);
+                List<T> subField(map.size());
+                forAll(subField, i)
+                {
+                    subField[i] = accessAndFlip
+                    (
+                        field,
+                        map[i],
+                        subHasFlip,
+                        negOp
+                    );
+                }
+                toNbr << subField;
+            }
+        }
+
+        // Subset myself
+        const labelList& mySubMap = subMap[Pstream::myProcNo()];
+
+        List<T> subField(mySubMap.size());
+        forAll(mySubMap, i)
+        {
+            subField[i] = accessAndFlip(field, mySubMap[i], subHasFlip, negOp);
+        }
+
+        // Receive sub field from myself (subField)
+        const labelList& map = constructMap[Pstream::myProcNo()];
+
+        field.setSize(constructSize);
+        field = nullValue;
+
+        flipAndCombine(map, constructHasFlip, subField, cop, negOp, field);
+
+        // Receive sub field from neighbour
+        for (label domain = 0; domain < Pstream::nProcs(); domain++)
+        {
+            const labelList& map = constructMap[domain];
+
+            if (domain != Pstream::myProcNo() && map.size())
+            {
+                IPstream fromNbr(Pstream::commsTypes::blocking, domain, 0, tag);
+                List<T> subField(fromNbr);
+
+                checkReceivedSize(domain, map.size(), subField.size());
+
+                flipAndCombine
+                (
+                    map,
+                    constructHasFlip,
+                    subField,
+                    cop,
+                    negOp,
+                    field
+                );
+            }
+        }
+    }
+    else if (commsType == Pstream::commsTypes::scheduled)
+    {
+        // Need to make sure I don't overwrite field with received data
+        // since the data might need to be sent to another processor. So
+        // allocate a new field for the results.
+        List<T> newField(constructSize, nullValue);
+
+        {
+            const labelList& mySubMap = subMap[Pstream::myProcNo()];
+
+            // Subset myself
+            List<T> subField(mySubMap.size());
+            forAll(subField, i)
+            {
+                subField[i] = accessAndFlip
+                (
+                    field,
+                    mySubMap[i],
+                    subHasFlip,
+                    negOp
+                );
+            }
+
+            // Receive sub field from myself (subField)
+            const labelList& map = constructMap[Pstream::myProcNo()];
+
+            flipAndCombine
+            (
+                map,
+                constructHasFlip,
+                subField,
+                cop,
+                negOp,
+                newField
+            );
+        }
+
+
+        // Schedule will already have pruned 0-sized comms
+        forAll(schedule, i)
+        {
+            const labelPair& twoProcs = schedule[i];
+            // twoProcs is a swap pair of processors. The first one is the
+            // one that needs to send first and then receive.
+
+            label sendProc = twoProcs[0];
+            label recvProc = twoProcs[1];
+
+            if (Pstream::myProcNo() == sendProc)
+            {
+                // I am send first, receive next
+                {
+                    OPstream toNbr
+                    (
+                        Pstream::commsTypes::scheduled,
+                        recvProc,
+                        0,
+                        tag
+                    );
+
+                    const labelList& map = subMap[recvProc];
+
+                    List<T> subField(map.size());
+                    forAll(subField, i)
+                    {
+                        subField[i] = accessAndFlip
+                        (
+                            field,
+                            map[i],
+                            subHasFlip,
+                            negOp
+                        );
+                    }
+                    toNbr << subField;
+                }
+                {
+                    IPstream fromNbr
+                    (
+                        Pstream::commsTypes::scheduled,
+                        recvProc,
+                        0,
+                        tag
+                    );
+                    List<T> subField(fromNbr);
+                    const labelList& map = constructMap[recvProc];
+
+                    checkReceivedSize(recvProc, map.size(), subField.size());
+
+                    flipAndCombine
+                    (
+                        map,
+                        constructHasFlip,
+                        subField,
+                        cop,
+                        negOp,
+                        newField
+                    );
+                }
+            }
+            else
+            {
+                // I am receive first, send next
+                {
+                    IPstream fromNbr
+                    (
+                        Pstream::commsTypes::scheduled,
+                        sendProc,
+                        0,
+                        tag
+                    );
+                    List<T> subField(fromNbr);
+                    const labelList& map = constructMap[sendProc];
+
+                    checkReceivedSize(sendProc, map.size(), subField.size());
+
+                    flipAndCombine
+                    (
+                        map,
+                        constructHasFlip,
+                        subField,
+                        cop,
+                        negOp,
+                        newField
+                    );
+                }
+                {
+                    OPstream toNbr
+                    (
+                        Pstream::commsTypes::scheduled,
+                        sendProc,
+                        0,
+                        tag
+                    );
+
+                    const labelList& map = subMap[sendProc];
+
+                    List<T> subField(map.size());
+                    forAll(subField, i)
+                    {
+                        subField[i] = accessAndFlip
+                        (
+                            field,
+                            map[i],
+                            subHasFlip,
+                            negOp
+                        );
+                    }
+                    toNbr << subField;
+                }
+            }
+        }
+        field.transfer(newField);
+    }
+    else if (commsType == Pstream::commsTypes::nonBlocking)
+    {
+        label nOutstanding = Pstream::nRequests();
+
+        if (!contiguous<T>())
+        {
+            PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking, tag);
+
+            // Stream data into buffer
+            for (label domain = 0; domain < Pstream::nProcs(); domain++)
+            {
+                const labelList& map = subMap[domain];
+
+                if (domain != Pstream::myProcNo() && map.size())
+                {
+                    // Put data into send buffer
+                    UOPstream toDomain(domain, pBufs);
+
+                    List<T> subField(map.size());
+                    forAll(subField, i)
+                    {
+                        subField[i] = accessAndFlip
+                        (
+                            field,
+                            map[i],
+                            subHasFlip,
+                            negOp
+                        );
+                    }
+                    toDomain << subField;
+                }
+            }
+
+            // Start receiving. Do not block.
+            pBufs.finishedSends(false);
+
+            {
+                // Set up 'send' to myself
+                const labelList& myMap = subMap[Pstream::myProcNo()];
+
+                List<T> mySubField(myMap.size());
+                forAll(myMap, i)
+                {
+                    mySubField[i] = accessAndFlip
+                    (
+                        field,
+                        myMap[i],
+                        subHasFlip,
+                        negOp
+                    );
+                }
+
+                // Combine bits. Note that can reuse field storage
+                field.setSize(constructSize);
+                field = nullValue;
+                // Receive sub field from myself
+                {
+                    const labelList& map = constructMap[Pstream::myProcNo()];
+
+                    flipAndCombine
+                    (
+                        map,
+                        constructHasFlip,
+                        mySubField,
+                        cop,
+                        negOp,
+                        field
+                    );
+                }
+            }
+
+            // Block ourselves, waiting only for the current comms
+            Pstream::waitRequests(nOutstanding);
+
+            // Consume
+            for (label domain = 0; domain < Pstream::nProcs(); domain++)
+            {
+                const labelList& map = constructMap[domain];
+
+                if (domain != Pstream::myProcNo() && map.size())
+                {
+                    UIPstream str(domain, pBufs);
+                    List<T> recvField(str);
+
+                    checkReceivedSize(domain, map.size(), recvField.size());
+
+                    flipAndCombine
+                    (
+                        map,
+                        constructHasFlip,
+                        recvField,
+                        cop,
+                        negOp,
+                        field
+                    );
+                }
+            }
+        }
+        else
+        {
+            // Set up sends to neighbours
+
+            List<List<T>> sendFields(Pstream::nProcs());
+
+            for (label domain = 0; domain < Pstream::nProcs(); domain++)
+            {
+                const labelList& map = subMap[domain];
+
+                if (domain != Pstream::myProcNo() && map.size())
+                {
+                    List<T>& subField = sendFields[domain];
+                    subField.setSize(map.size());
+                    forAll(map, i)
+                    {
+                        subField[i] = accessAndFlip
+                        (
+                            field,
+                            map[i],
+                            subHasFlip,
+                            negOp
+                        );
+                    }
+
+                    OPstream::write
+                    (
+                        Pstream::commsTypes::nonBlocking,
+                        domain,
+                        reinterpret_cast<const char*>(subField.begin()),
+                        subField.size()*sizeof(T),
+                        tag
+                    );
+                }
+            }
+
+            // Set up receives from neighbours
+
+            List<List<T>> recvFields(Pstream::nProcs());
+
+            for (label domain = 0; domain < Pstream::nProcs(); domain++)
+            {
+                const labelList& map = constructMap[domain];
+
+                if (domain != Pstream::myProcNo() && map.size())
+                {
+                    recvFields[domain].setSize(map.size());
+                    UIPstream::read
+                    (
+                        Pstream::commsTypes::nonBlocking,
+                        domain,
+                        reinterpret_cast<char*>(recvFields[domain].begin()),
+                        recvFields[domain].size()*sizeof(T),
+                        tag
+                    );
+                }
+            }
+
+            // Set up 'send' to myself
+
+            {
+                const labelList& map = subMap[Pstream::myProcNo()];
+
+                List<T>& subField = sendFields[Pstream::myProcNo()];
+                subField.setSize(map.size());
+                forAll(map, i)
+                {
+                    subField[i] = accessAndFlip
+                    (
+                        field,
+                        map[i],
+                        subHasFlip,
+                        negOp
+                    );
+                }
+            }
+
+
+            // Combine bits. Note that can reuse field storage
+
+            field.setSize(constructSize);
+            field = nullValue;
+
+            // Receive sub field from myself (subField)
+            {
+                const labelList& map = constructMap[Pstream::myProcNo()];
+                const List<T>& subField = sendFields[Pstream::myProcNo()];
+
+                flipAndCombine
+                (
+                    map,
+                    constructHasFlip,
+                    subField,
+                    cop,
+                    negOp,
+                    field
+                );
+            }
+
+
+            // Wait for all to finish
+
+            Pstream::waitRequests(nOutstanding);
+
+
+            // Collect neighbour fields
+
+            for (label domain = 0; domain < Pstream::nProcs(); domain++)
+            {
+                const labelList& map = constructMap[domain];
+
+                if (domain != Pstream::myProcNo() && map.size())
+                {
+                    const List<T>& subField = recvFields[domain];
+
+                    checkReceivedSize(domain, map.size(), subField.size());
+
+                    flipAndCombine
+                    (
+                        map,
+                        constructHasFlip,
+                        subField,
+                        cop,
+                        negOp,
+                        field
+                    );
+                }
+            }
+        }
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Unknown communication schedule " << int(commsType)
+            << abort(FatalError);
+    }
+}
+
+
+template<class T>
+void Foam::distributionMapBase::send
+(
+    PstreamBuffers& pBufs,
+    const List<T>& field
+)
+const
+{
+    // Stream data into buffer
+    for (label domain = 0; domain < Pstream::nProcs(); domain++)
+    {
+        const labelList& map = subMap_[domain];
+
+        if (map.size())
+        {
+            // Put data into send buffer
+            UOPstream toDomain(domain, pBufs);
+
+            List<T> subField(map.size());
+            forAll(subField, i)
+            {
+                subField[i] = accessAndFlip
+                (
+                    field,
+                    map[i],
+                    subHasFlip_,
+                    flipOp()
+                );
+            }
+            toDomain << subField;
+        }
+    }
+
+    // Start sending and receiving but do not block.
+    pBufs.finishedSends(false);
+}
+
+
+template<class T>
+void Foam::distributionMapBase::receive(PstreamBuffers& pBufs, List<T>& field)
+const
+{
+    // Consume
+    field.setSize(constructSize_);
+
+    for (label domain = 0; domain < Pstream::nProcs(); domain++)
+    {
+        const labelList& map = constructMap_[domain];
+
+        if (map.size())
+        {
+            UIPstream str(domain, pBufs);
+            List<T> recvField(str);
+
+            if (recvField.size() != map.size())
+            {
+                FatalErrorInFunction
+                    << "Expected from processor " << domain
+                    << " " << map.size() << " but received "
+                    << recvField.size() << " elements."
+                    << abort(FatalError);
+            }
+
+            flipAndCombine
+            (
+                map,
+                constructHasFlip_,
+                recvField,
+                eqOp<T>(),
+                flipOp(),
+                field
+            );
+        }
+    }
+}
+
+
+//- Distribute data using default commsType.
+template<class T, class negateOp>
+void Foam::distributionMapBase::distribute
+(
+    List<T>& fld,
+    const negateOp& negOp,
+    const int tag
+) const
+{
+    if (Pstream::defaultCommsType == Pstream::commsTypes::nonBlocking)
+    {
+        distribute
+        (
+            Pstream::commsTypes::nonBlocking,
+            List<labelPair>(),
+            constructSize_,
+            subMap_,
+            subHasFlip_,
+            constructMap_,
+            constructHasFlip_,
+            fld,
+            negOp,
+            tag
+        );
+    }
+    else if (Pstream::defaultCommsType == Pstream::commsTypes::scheduled)
+    {
+        distribute
+        (
+            Pstream::commsTypes::scheduled,
+            schedule(),
+            constructSize_,
+            subMap_,
+            subHasFlip_,
+            constructMap_,
+            constructHasFlip_,
+            fld,
+            negOp,
+            tag
+        );
+    }
+    else
+    {
+        distribute
+        (
+            Pstream::commsTypes::blocking,
+            List<labelPair>(),
+            constructSize_,
+            subMap_,
+            subHasFlip_,
+            constructMap_,
+            constructHasFlip_,
+            fld,
+            negOp,
+            tag
+        );
+    }
+}
+
+
+//- Distribute data using default commsType.
+template<class T>
+void Foam::distributionMapBase::distribute
+(
+    List<T>& fld,
+    const int tag
+) const
+{
+    distribute(fld, flipOp(), tag);
+}
+
+
+//- Distribute data using default commsType.
+template<class T>
+void Foam::distributionMapBase::distribute
+(
+    DynamicList<T>& fld,
+    const int tag
+) const
+{
+    fld.shrink();
+
+    List<T>& fldList = static_cast<List<T>& >(fld);
+
+    distribute(fldList, tag);
+
+    fld.setCapacity(fldList.size());
+}
+
+
+//- Reverse distribute data using default commsType.
+template<class T>
+void Foam::distributionMapBase::reverseDistribute
+(
+    const label constructSize,
+    List<T>& fld,
+    const int tag
+) const
+{
+    if (Pstream::defaultCommsType == Pstream::commsTypes::nonBlocking)
+    {
+        distribute
+        (
+            Pstream::commsTypes::nonBlocking,
+            List<labelPair>(),
+            constructSize,
+            constructMap_,
+            constructHasFlip_,
+            subMap_,
+            subHasFlip_,
+            fld,
+            flipOp(),
+            tag
+        );
+    }
+    else if (Pstream::defaultCommsType == Pstream::commsTypes::scheduled)
+    {
+        distribute
+        (
+            Pstream::commsTypes::scheduled,
+            schedule(),
+            constructSize,
+            constructMap_,
+            constructHasFlip_,
+            subMap_,
+            subHasFlip_,
+            fld,
+            flipOp(),
+            tag
+        );
+    }
+    else
+    {
+        distribute
+        (
+            Pstream::commsTypes::blocking,
+            List<labelPair>(),
+            constructSize,
+            constructMap_,
+            constructHasFlip_,
+            subMap_,
+            subHasFlip_,
+            fld,
+            flipOp(),
+            tag
+        );
+    }
+}
+
+
+//- Reverse distribute data using default commsType.
+//  Since constructSize might be larger than supplied size supply
+//  a nullValue
+template<class T>
+void Foam::distributionMapBase::reverseDistribute
+(
+    const label constructSize,
+    const T& nullValue,
+    List<T>& fld,
+    const int tag
+) const
+{
+    if (Pstream::defaultCommsType == Pstream::commsTypes::nonBlocking)
+    {
+        distribute
+        (
+            Pstream::commsTypes::nonBlocking,
+            List<labelPair>(),
+            constructSize,
+            constructMap_,
+            constructHasFlip_,
+            subMap_,
+            subHasFlip_,
+            fld,
+            eqOp<T>(),
+            flipOp(),
+            nullValue,
+            tag
+        );
+    }
+    else if (Pstream::defaultCommsType == Pstream::commsTypes::scheduled)
+    {
+        distribute
+        (
+            Pstream::commsTypes::scheduled,
+            schedule(),
+            constructSize,
+            constructMap_,
+            constructHasFlip_,
+            subMap_,
+            subHasFlip_,
+            fld,
+            eqOp<T>(),
+            flipOp(),
+            nullValue,
+            tag
+        );
+    }
+    else
+    {
+        distribute
+        (
+            Pstream::commsTypes::blocking,
+            List<labelPair>(),
+            constructSize,
+            constructMap_,
+            constructHasFlip_,
+            subMap_,
+            subHasFlip_,
+            fld,
+            eqOp<T>(),
+            flipOp(),
+            nullValue,
+            tag
+        );
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/distributionMapTemplates.C b/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/distributionMapTemplates.C
new file mode 100644
index 0000000000..a6bbc6cde7
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/distributionMapTemplates.C
@@ -0,0 +1,271 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "Pstream.H"
+#include "PstreamBuffers.H"
+#include "PstreamCombineReduceOps.H"
+#include "globalIndexAndTransform.H"
+#include "transformField.H"
+#include "flipOp.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class T>
+void Foam::distributionMap::applyDummyTransforms(List<T>& field) const
+{
+    forAll(transformElements_, trafoI)
+    {
+        const labelList& elems = transformElements_[trafoI];
+
+        label n = transformStart_[trafoI];
+
+        forAll(elems, i)
+        {
+            field[n++] = field[elems[i]];
+        }
+    }
+}
+
+
+template<class T>
+void Foam::distributionMap::applyDummyInverseTransforms(List<T>& field) const
+{
+    forAll(transformElements_, trafoI)
+    {
+        const labelList& elems = transformElements_[trafoI];
+        label n = transformStart_[trafoI];
+
+        forAll(elems, i)
+        {
+            field[elems[i]] = field[n++];
+        }
+    }
+}
+
+
+template<class T, class TransformOp>   //, class CombineOp>
+void Foam::distributionMap::applyTransforms
+(
+    const globalIndexAndTransform& globalTransforms,
+    List<T>& field,
+    const TransformOp& top
+) const
+{
+    const List<transformer>& totalTransform =
+        globalTransforms.transformPermutations();
+
+    forAll(totalTransform, trafoI)
+    {
+        const transformer& vt = totalTransform[trafoI];
+        const labelList& elems = transformElements_[trafoI];
+        label n = transformStart_[trafoI];
+
+        // Could be optimised to avoid memory allocations
+        List<T> transformFld(UIndirectList<T>(field, elems));
+        top(vt, true, transformFld);
+
+        forAll(transformFld, i)
+        {
+            // cop(field[n++], transformFld[i]);
+            field[n++] = transformFld[i];
+        }
+    }
+}
+
+
+template<class T, class TransformOp>   //, class CombineOp>
+void Foam::distributionMap::applyInverseTransforms
+(
+    const globalIndexAndTransform& globalTransforms,
+    List<T>& field,
+    const TransformOp& top
+) const
+{
+    const List<transformer>& totalTransform =
+        globalTransforms.transformPermutations();
+
+    forAll(totalTransform, trafoI)
+    {
+        const transformer& vt = totalTransform[trafoI];
+        const labelList& elems = transformElements_[trafoI];
+        label n = transformStart_[trafoI];
+
+        // Could be optimised to avoid memory allocations
+        List<T> transformFld(SubList<T>(field, elems.size(), n));
+        top(vt, false, transformFld);
+
+        forAll(transformFld, i)
+        {
+            // cop(field[elems[i]], transformFld[i]);
+            field[elems[i]] = transformFld[i];
+        }
+    }
+}
+
+
+template<class T, class negateOp>
+void Foam::distributionMap::distribute
+(
+    List<T>& fld,
+    const negateOp& negOp,
+    const bool dummyTransform,
+    const int tag
+) const
+{
+    distributionMapBase::distribute(fld, negOp, tag);
+
+    //- Fill in transformed slots with copies
+    if (dummyTransform)
+    {
+        applyDummyTransforms(fld);
+    }
+}
+
+
+template<class T>
+void Foam::distributionMap::distribute
+(
+    List<T>& fld,
+    const bool dummyTransform,
+    const int tag
+) const
+{
+    distribute(fld, flipOp(), dummyTransform, tag);
+}
+
+
+template<class T>
+void Foam::distributionMap::distribute
+(
+    DynamicList<T>& fld,
+    const bool dummyTransform,
+    const int tag
+) const
+{
+    fld.shrink();
+
+    List<T>& fldList = static_cast<List<T>& >(fld);
+
+    distribute(fldList, dummyTransform, tag);
+
+    fld.setCapacity(fldList.size());
+}
+
+
+template<class T>
+void Foam::distributionMap::reverseDistribute
+(
+    const label constructSize,
+    List<T>& fld,
+    const bool dummyTransform,
+    const int tag
+) const
+{
+    if (dummyTransform)
+    {
+        applyDummyInverseTransforms(fld);
+    }
+
+    distributionMapBase::reverseDistribute(constructSize, fld, tag);
+}
+
+
+template<class T>
+void Foam::distributionMap::reverseDistribute
+(
+    const label constructSize,
+    const T& nullValue,
+    List<T>& fld,
+    const bool dummyTransform,
+    const int tag
+) const
+{
+    if (dummyTransform)
+    {
+        applyDummyInverseTransforms(fld);
+    }
+
+    distributionMapBase::reverseDistribute(constructSize, nullValue, fld, tag);
+}
+
+
+template<class T, class TransformOp>
+void Foam::distributionMap::distribute
+(
+    const globalIndexAndTransform& git,
+    List<T>& fld,
+    const TransformOp& top,
+    const int tag
+) const
+{
+    // Distribute. Leave out dummy transforms since we're doing them ourselves
+    distribute(fld, false, tag);
+    // Do transforms
+    applyTransforms(git, fld, top);
+}
+
+
+template<class T, class TransformOp>
+void Foam::distributionMap::reverseDistribute
+(
+    const globalIndexAndTransform& git,
+    const label constructSize,
+    List<T>& fld,
+    const TransformOp& top,
+    const int tag
+) const
+{
+    // Fill slots with reverse-transformed data. Note that it also copies
+    // back into the non-remote part of fld even though these values are not
+    // used.
+    applyInverseTransforms(git, fld, top);
+
+    // And send back (the remote slots). Disable dummy transformations.
+    reverseDistribute(constructSize, fld, false, tag);
+}
+
+
+template<class T, class TransformOp>
+void Foam::distributionMap::reverseDistribute
+(
+    const globalIndexAndTransform& git,
+    const label constructSize,
+    const T& nullValue,
+    List<T>& fld,
+    const TransformOp& top,
+    const int tag
+) const
+{
+    // Fill slots with reverse-transformed data Note that it also copies
+    // back into the non-remote part of fld even though these values are not
+    // used.
+    applyInverseTransforms(git, fld, top);   //, eqOp<T>());
+
+    // And send back (the remote slots) Disable dummy transformations.
+    reverseDistribute(constructSize, nullValue, fld, false, tag);
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/lagrangianDistributionMap.H b/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/lagrangianDistributionMap.H
new file mode 100644
index 0000000000..c611fc4561
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/lagrangianDistributionMap.H
@@ -0,0 +1,122 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::lagrangianDistributionMap
+
+Description
+    Class containing mesh-to-mesh mapping information for particles
+
+SourceFiles
+    lagrangianDistributionMap.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef lagrangianDistributionMap_H
+#define lagrangianDistributionMap_H
+
+#include "distributionMap.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class polyTopoChangeMap;
+
+/*---------------------------------------------------------------------------*\
+                   Class lagrangianDistributionMap Declaration
+\*---------------------------------------------------------------------------*/
+
+class lagrangianDistributionMap
+{
+    // Private Data
+
+        //- Map to distribute particles
+        const distributionMap particleMap_;
+
+        //- Per element in subsetted mesh the cell label
+        const labelListList constructCellLabels_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        lagrangianDistributionMap
+        (
+            const label nNewParticles,
+            labelListList&& subParticleMap,
+            labelListList&& constructParticleMap,
+            labelListList&& constructCellLabels
+        )
+        :
+            particleMap_(nNewParticles, subParticleMap, constructParticleMap),
+            constructCellLabels_(constructCellLabels)
+        {}
+
+
+    // Member Functions
+
+        // Access
+
+            //- Distribution map
+            const distributionMap& particleMap() const
+            {
+                return particleMap_;
+            }
+
+            //- Per received particle the destination cell label
+            const labelListList& constructCellLabels() const
+            {
+                return constructCellLabels_;
+            }
+
+
+        // Edit
+
+            //- Distribute list of lagrangian data
+            template<class T>
+            void distributeLagrangianData(List<T>& lst) const
+            {
+                particleMap_.distribute(lst);
+            }
+
+            //- Correct for topo change.
+            void updateMesh(const polyTopoChangeMap&)
+            {
+                NotImplemented;
+            }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/polyDistributionMap.C b/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/polyDistributionMap.C
new file mode 100644
index 0000000000..5c23682468
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/polyDistributionMap.C
@@ -0,0 +1,366 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "polyDistributionMap.H"
+#include "polyMesh.H"
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::polyDistributionMap::calcPatchSizes()
+{
+    oldPatchSizes_.setSize(oldPatchStarts_.size());
+
+    if (oldPatchStarts_.size())
+    {
+        // Calculate old patch sizes
+        for (label patchi = 0; patchi < oldPatchStarts_.size() - 1; patchi++)
+        {
+            oldPatchSizes_[patchi] =
+                oldPatchStarts_[patchi + 1] - oldPatchStarts_[patchi];
+        }
+
+        // Set the last one by hand
+        const label lastPatchID = oldPatchStarts_.size() - 1;
+
+        oldPatchSizes_[lastPatchID] = nOldFaces_ - oldPatchStarts_[lastPatchID];
+
+        if (min(oldPatchSizes_) < 0)
+        {
+            FatalErrorInFunction
+                << "Calculated negative old patch size:" << oldPatchSizes_ << nl
+                << "Error in mapping data" << abort(FatalError);
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::polyDistributionMap::polyDistributionMap()
+:
+    nOldPoints_(0),
+    nOldFaces_(0),
+    nOldCells_(0),
+    oldPatchSizes_(0),
+    oldPatchStarts_(0),
+    oldPatchNMeshPoints_(0),
+    pointMap_(),
+    faceMap_(),
+    cellMap_(),
+    patchMap_()
+{}
+
+
+Foam::polyDistributionMap::polyDistributionMap
+(
+    const polyMesh& mesh,
+
+    // mesh before changes
+    const label nOldPoints,
+    const label nOldFaces,
+    const label nOldCells,
+    labelList&& oldPatchStarts,
+    labelList&& oldPatchNMeshPoints,
+
+    // how to subset pieces of mesh to send across
+    labelListList&& subPointMap,
+    labelListList&& subFaceMap,
+    labelListList&& subCellMap,
+    labelListList&& subPatchMap,
+
+    // how to reconstruct received mesh
+    labelListList&& constructPointMap,
+    labelListList&& constructFaceMap,
+    labelListList&& constructCellMap,
+    labelListList&& constructPatchMap,
+
+    const bool subFaceHasFlip,
+    const bool constructFaceHasFlip
+)
+:
+    nOldPoints_(nOldPoints),
+    nOldFaces_(nOldFaces),
+    nOldCells_(nOldCells),
+    oldPatchSizes_(oldPatchStarts.size()),
+    oldPatchStarts_(move(oldPatchStarts)),
+    oldPatchNMeshPoints_(move(oldPatchNMeshPoints)),
+    pointMap_(mesh.nPoints(), move(subPointMap), move(constructPointMap)),
+    faceMap_
+    (
+        mesh.nFaces(),
+        move(subFaceMap),
+        move(constructFaceMap),
+        move(subFaceHasFlip),
+        constructFaceHasFlip
+    ),
+    cellMap_(mesh.nCells(), move(subCellMap), move(constructCellMap)),
+    patchMap_
+    (
+        mesh.boundaryMesh().size(),
+        move(subPatchMap),
+        move(constructPatchMap)
+    )
+{
+    calcPatchSizes();
+}
+
+
+Foam::polyDistributionMap::polyDistributionMap
+(
+    // mesh before changes
+    const label nOldPoints,
+    const label nOldFaces,
+    const label nOldCells,
+    labelList&& oldPatchStarts,
+    labelList&& oldPatchNMeshPoints,
+
+    // how to transfer pieces of mesh
+    distributionMap&& pointMap,
+    distributionMap&& faceMap,
+    distributionMap&& cellMap,
+    distributionMap&& patchMap
+)
+:
+    nOldPoints_(nOldPoints),
+    nOldFaces_(nOldFaces),
+    nOldCells_(nOldCells),
+    oldPatchSizes_(oldPatchStarts.size()),
+    oldPatchStarts_(move(oldPatchStarts)),
+    oldPatchNMeshPoints_(move(oldPatchNMeshPoints)),
+    pointMap_(move(pointMap)),
+    faceMap_(move(faceMap)),
+    cellMap_(move(cellMap)),
+    patchMap_(move(patchMap))
+{
+    calcPatchSizes();
+}
+
+
+Foam::polyDistributionMap::polyDistributionMap
+(
+    polyDistributionMap&& map
+)
+:
+    nOldPoints_(map.nOldPoints_),
+    nOldFaces_(map.nOldFaces_),
+    nOldCells_(map.nOldCells_),
+    oldPatchSizes_(move(map.oldPatchSizes_)),
+    oldPatchStarts_(move(map.oldPatchStarts_)),
+    oldPatchNMeshPoints_(move(map.oldPatchNMeshPoints_)),
+    pointMap_(move(map.pointMap_)),
+    faceMap_(move(map.faceMap_)),
+    cellMap_(move(map.cellMap_)),
+    patchMap_(move(map.patchMap_))
+{}
+
+
+Foam::polyDistributionMap::polyDistributionMap(Istream& is)
+{
+    is  >> *this;
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::polyDistributionMap::transfer(polyDistributionMap& rhs)
+{
+    nOldPoints_ = rhs.nOldPoints_;
+    nOldFaces_ = rhs.nOldFaces_;
+    nOldCells_ = rhs.nOldCells_;
+    oldPatchSizes_.transfer(rhs.oldPatchSizes_);
+    oldPatchStarts_.transfer(rhs.oldPatchStarts_);
+    oldPatchNMeshPoints_.transfer(rhs.oldPatchNMeshPoints_);
+    pointMap_.transfer(rhs.pointMap_);
+    faceMap_.transfer(rhs.faceMap_);
+    cellMap_.transfer(rhs.cellMap_);
+    patchMap_.transfer(rhs.patchMap_);
+}
+
+
+void Foam::polyDistributionMap::distributePointIndices(labelList& lst) const
+{
+    // Construct boolList from selected elements
+    boolList isSelected
+    (
+        createWithValues<boolList>
+        (
+            nOldPoints(),
+            false,
+            lst,
+            true
+        )
+    );
+
+    // Distribute
+    distributePointData(isSelected);
+
+    // Collect selected elements
+    lst = findIndices(isSelected, true);
+}
+
+
+void Foam::polyDistributionMap::distributeFaceIndices(labelList& lst) const
+{
+    // Construct boolList from selected elements
+    boolList isSelected
+    (
+        createWithValues<boolList>
+        (
+            nOldFaces(),
+            false,
+            lst,
+            true
+        )
+    );
+
+    // Distribute
+    distributeFaceData(isSelected);
+
+    // Collect selected elements
+    lst = findIndices(isSelected, true);
+}
+
+
+void Foam::polyDistributionMap::distributeCellIndices(labelList& lst) const
+{
+    // Construct boolList from selected elements
+    boolList isSelected
+    (
+        createWithValues<boolList>
+        (
+            nOldCells(),
+            false,
+            lst,
+            true
+        )
+    );
+
+    // Distribute
+    distributeCellData(isSelected);
+
+    // Collect selected elements
+    lst = findIndices(isSelected, true);
+}
+
+
+void Foam::polyDistributionMap::distributePatchIndices(labelList& lst) const
+{
+    // Construct boolList from selected elements
+    boolList isSelected
+    (
+        createWithValues<boolList>
+        (
+            oldPatchStarts().size(),    // nOldPatches
+            false,
+            lst,
+            true
+        )
+    );
+
+    // Distribute
+    distributePatchData(isSelected);
+
+    // Collect selected elements
+    lst = findIndices(isSelected, true);
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+void Foam::polyDistributionMap::operator=
+(
+    const polyDistributionMap& rhs
+)
+{
+    nOldPoints_ = rhs.nOldPoints_;
+    nOldFaces_ = rhs.nOldFaces_;
+    nOldCells_ = rhs.nOldCells_;
+    oldPatchSizes_ = rhs.oldPatchSizes_;
+    oldPatchStarts_ = rhs.oldPatchStarts_;
+    oldPatchNMeshPoints_ = rhs.oldPatchNMeshPoints_;
+    pointMap_ = rhs.pointMap_;
+    faceMap_ = rhs.faceMap_;
+    cellMap_ = rhs.cellMap_;
+    patchMap_ = rhs.patchMap_;
+}
+
+
+void Foam::polyDistributionMap::operator=(polyDistributionMap&& rhs)
+{
+    nOldPoints_ = rhs.nOldPoints_;
+    nOldFaces_ = rhs.nOldFaces_;
+    nOldCells_ = rhs.nOldCells_;
+    oldPatchSizes_ = move(rhs.oldPatchSizes_);
+    oldPatchStarts_ = move(rhs.oldPatchStarts_);
+    oldPatchNMeshPoints_ = move(rhs.oldPatchNMeshPoints_);
+    pointMap_ = move(rhs.pointMap_);
+    faceMap_ = move(rhs.faceMap_);
+    cellMap_ = move(rhs.cellMap_);
+    patchMap_ = move(rhs.patchMap_);
+}
+
+
+// * * * * * * * * * * * * * * Istream Operator  * * * * * * * * * * * * * * //
+
+Foam::Istream& Foam::operator>>(Istream& is, polyDistributionMap& map)
+{
+    is.fatalCheck("operator>>(Istream&, polyDistributionMap&)");
+
+    is  >> map.nOldPoints_
+        >> map.nOldFaces_
+        >> map.nOldCells_
+        >> map.oldPatchSizes_
+        >> map.oldPatchStarts_
+        >> map.oldPatchNMeshPoints_
+        >> map.pointMap_
+        >> map.faceMap_
+        >> map.cellMap_
+        >> map.patchMap_;
+
+    return is;
+}
+
+
+// * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const polyDistributionMap& map)
+{
+    os  << map.nOldPoints_
+        << token::SPACE << map.nOldFaces_
+        << token::SPACE << map.nOldCells_ << token::NL
+        << map.oldPatchSizes_ << token::NL
+        << map.oldPatchStarts_ << token::NL
+        << map.oldPatchNMeshPoints_ << token::NL
+        << map.pointMap_ << token::NL
+        << map.faceMap_ << token::NL
+        << map.cellMap_ << token::NL
+        << map.patchMap_;
+
+    return os;
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/polyDistributionMap.H b/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/polyDistributionMap.H
new file mode 100644
index 0000000000..bfdc481cb8
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyDistributionMap/polyDistributionMap.H
@@ -0,0 +1,312 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::polyDistributionMap
+
+Description
+    Class containing mesh-to-mesh mapping information after a mesh distribution
+    where we send parts of meshes (using subsetting) to other processors
+    and receive and reconstruct mesh.
+
+    We store mapping from the bits-to-send to the complete starting mesh
+    (subXXXMap) and from the received bits to their location in the new
+    mesh (constructXXXMap).
+
+SourceFiles
+    polyDistributionMap.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef polyDistributionMap_H
+#define polyDistributionMap_H
+
+#include "distributionMap.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class polyTopoChangeMap;
+class polyMesh;
+
+
+// Forward declaration of friend functions and operators
+
+class polyDistributionMap;
+
+Istream& operator>>(Istream&, polyDistributionMap&);
+Ostream& operator<<(Ostream&, const polyDistributionMap&);
+
+
+/*---------------------------------------------------------------------------*\
+                    Class polyDistributionMap Declaration
+\*---------------------------------------------------------------------------*/
+
+class polyDistributionMap
+{
+    // Private Data
+
+        //- Number of old live points
+        label nOldPoints_;
+
+        //- Number of old live faces
+        label nOldFaces_;
+
+        //- Number of old live cells
+        label nOldCells_;
+
+        //- List of the old patch sizes
+        labelList oldPatchSizes_;
+
+        //- List of the old patch start labels
+        labelList oldPatchStarts_;
+
+        //- List of numbers of mesh points per old patch
+        labelList oldPatchNMeshPoints_;
+
+
+        //- Point distribute map
+        distributionMap pointMap_;
+
+        //- Face distribute map
+        distributionMap faceMap_;
+
+        //- Cell distribute map
+        distributionMap cellMap_;
+
+        //- Patch distribute map
+        distributionMap patchMap_;
+
+
+    // Private Member Functions
+
+        void calcPatchSizes();
+
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        polyDistributionMap();
+
+        //- Move constructor from components.
+        //  Note that mesh has to be changed already
+        //  since uses mesh.nPoints etc as the new size.
+        polyDistributionMap
+        (
+            const polyMesh& mesh,
+
+            // mesh before changes
+            const label nOldPoints,
+            const label nOldFaces,
+            const label nOldCells,
+            labelList&& oldPatchStarts,
+            labelList&& oldPatchNMeshPoints,
+
+            // how to subset pieces of mesh to send across
+            labelListList&& subPointMap,
+            labelListList&& subFaceMap,
+            labelListList&& subCellMap,
+            labelListList&& subPatchMap,
+
+            // how to reconstruct received mesh
+            labelListList&& constructPointMap,
+            labelListList&& constructFaceMap,
+            labelListList&& constructCellMap,
+            labelListList&& constructPatchMap,
+
+            const bool subFaceHasFlip = false,
+            const bool constructFaceHasFlip = false
+        );
+
+        //- Move constructor from components
+        polyDistributionMap
+        (
+            // mesh before changes
+            const label nOldPoints,
+            const label nOldFaces,
+            const label nOldCells,
+            labelList&& oldPatchStarts,
+            labelList&& oldPatchNMeshPoints,
+
+            // how to subset pieces of mesh to send across
+            distributionMap&& pointMap,
+            distributionMap&& faceMap,
+            distributionMap&& cellMap,
+            distributionMap&& patchMap
+        );
+
+        //- Move constructor
+        polyDistributionMap(polyDistributionMap&&);
+
+        //- Construct from Istream
+        polyDistributionMap(Istream&);
+
+        //- Disallow default bitwise copy construction
+        polyDistributionMap(const polyDistributionMap&) = delete;
+
+
+    // Member Functions
+
+        // Access
+
+            //- Number of points in mesh before distribution
+            label nOldPoints() const
+            {
+                return nOldPoints_;
+            }
+
+            //- Number of faces in mesh before distribution
+            label nOldFaces() const
+            {
+                return nOldFaces_;
+            }
+
+            //- Number of cells in mesh before distribution
+            label nOldCells() const
+            {
+                return nOldCells_;
+            }
+
+            //- List of the old patch sizes
+            const labelList& oldPatchSizes() const
+            {
+                return oldPatchSizes_;
+            }
+
+            //- List of the old patch start labels
+            const labelList& oldPatchStarts() const
+            {
+                return oldPatchStarts_;
+            }
+
+            //- List of numbers of mesh points per old patch
+            const labelList& oldPatchNMeshPoints() const
+            {
+                return oldPatchNMeshPoints_;
+            }
+
+            //- Point distribute map
+            const distributionMap& pointMap() const
+            {
+                return pointMap_;
+            }
+
+            //- Face distribute map
+            const distributionMap& faceMap() const
+            {
+                return faceMap_;
+            }
+
+            //- Cell distribute map
+            const distributionMap& cellMap() const
+            {
+                return cellMap_;
+            }
+
+            //- Patch distribute map
+            const distributionMap& patchMap() const
+            {
+                return patchMap_;
+            }
+
+
+        // Other
+
+            //- Transfer the contents of the argument and annul the argument.
+            void transfer(polyDistributionMap&);
+
+            //- Distribute list of point data
+            template<class T>
+            void distributePointData(List<T>& lst) const
+            {
+                pointMap_.distribute(lst);
+            }
+
+            //- Distribute list of face data
+            template<class T>
+            void distributeFaceData(List<T>& lst) const
+            {
+                faceMap_.distribute(lst);
+            }
+
+            //- Distribute list of cell data
+            template<class T>
+            void distributeCellData(List<T>& lst) const
+            {
+                cellMap_.distribute(lst);
+            }
+
+            //- Distribute list of patch data
+            template<class T>
+            void distributePatchData(List<T>& lst) const
+            {
+                patchMap_.distribute(lst);
+            }
+
+
+            //- Distribute list of point/face/cell/patch indices.
+            //  (Converts to boolList, distributes boolList and reconstructs)
+            void distributePointIndices(labelList& pointIDs) const;
+
+            void distributeFaceIndices(labelList& faceIDs) const;
+            void distributeCellIndices(labelList& cellIDs) const;
+            void distributePatchIndices(labelList& patchIDs) const;
+
+
+            //- Correct for topo change.
+            void updateMesh(const polyTopoChangeMap&)
+            {
+                NotImplemented;
+            }
+
+
+    // Member Operators
+
+        void operator=(const polyDistributionMap&);
+
+        void operator=(polyDistributionMap&&);
+
+
+    // IOstream Operators
+
+        //- Read dictionary from Istream
+        friend Istream& operator>>(Istream&, polyDistributionMap&);
+
+        //- Write dictionary to Ostream
+        friend Ostream& operator<<(Ostream&, const polyDistributionMap&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/cellMapper/cellMapper.C b/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/cellMapper/cellMapper.C
new file mode 100644
index 0000000000..efb381241e
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/cellMapper/cellMapper.C
@@ -0,0 +1,463 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "cellMapper.H"
+#include "demandDrivenData.H"
+#include "polyMesh.H"
+#include "polyTopoChangeMap.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::cellMapper::calcAddressing() const
+{
+    if
+    (
+        directAddrPtr_
+     || interpolationAddrPtr_
+     || weightsPtr_
+     || insertedCellLabelsPtr_
+    )
+    {
+        FatalErrorInFunction
+            << "Addressing already calculated."
+            << abort(FatalError);
+    }
+
+    if (direct())
+    {
+        // Direct addressing, no weights
+
+        directAddrPtr_ = new labelList(mpm_.cellMap());
+        labelList& directAddr = *directAddrPtr_;
+
+        // Not necessary to resize the list as there are no retired cells
+        // directAddr.setSize(mesh_.nCells());
+
+        insertedCellLabelsPtr_ = new labelList(mesh_.nCells());
+        labelList& insertedCells = *insertedCellLabelsPtr_;
+
+        label nInsertedCells = 0;
+
+        forAll(directAddr, celli)
+        {
+            if (directAddr[celli] < 0)
+            {
+                // Found inserted cell
+                directAddr[celli] = 0;
+                insertedCells[nInsertedCells] = celli;
+                nInsertedCells++;
+            }
+        }
+
+        insertedCells.setSize(nInsertedCells);
+    }
+    else
+    {
+        // Interpolative addressing
+
+        interpolationAddrPtr_ = new labelListList(mesh_.nCells());
+        labelListList& addr = *interpolationAddrPtr_;
+
+        weightsPtr_ = new scalarListList(mesh_.nCells());
+        scalarListList& w = *weightsPtr_;
+
+        const List<objectMap>& cfp = mpm_.cellsFromPointsMap();
+
+        forAll(cfp, cfpI)
+        {
+            // Get addressing
+            const labelList& mo = cfp[cfpI].masterObjects();
+
+            label celli = cfp[cfpI].index();
+
+            if (addr[celli].size())
+            {
+                FatalErrorInFunction
+                    << "Master cell " << celli
+                    << " mapped from point cells " << mo
+                    << " already destination of mapping." << abort(FatalError);
+            }
+
+            // Map from masters, uniform weights
+            addr[celli] = mo;
+            w[celli] = scalarList(mo.size(), 1.0/mo.size());
+        }
+
+        const List<objectMap>& cfe = mpm_.cellsFromEdgesMap();
+
+        forAll(cfe, cfeI)
+        {
+            // Get addressing
+            const labelList& mo = cfe[cfeI].masterObjects();
+
+            label celli = cfe[cfeI].index();
+
+            if (addr[celli].size())
+            {
+                FatalErrorInFunction
+                    << "Master cell " << celli
+                    << " mapped from edge cells " << mo
+                    << " already destination of mapping." << abort(FatalError);
+            }
+
+            // Map from masters, uniform weights
+            addr[celli] = mo;
+            w[celli] = scalarList(mo.size(), 1.0/mo.size());
+        }
+
+        const List<objectMap>& cff = mpm_.cellsFromFacesMap();
+
+        forAll(cff, cffI)
+        {
+            // Get addressing
+            const labelList& mo = cff[cffI].masterObjects();
+
+            label celli = cff[cffI].index();
+
+            if (addr[celli].size())
+            {
+                FatalErrorInFunction
+                    << "Master cell " << celli
+                    << " mapped from face cells " << mo
+                    << " already destination of mapping." << abort(FatalError);
+            }
+
+            // Map from masters, uniform weights
+            addr[celli] = mo;
+            w[celli] = scalarList(mo.size(), 1.0/mo.size());
+        }
+
+        // Volume conservative mapping if possible
+
+        const List<objectMap>& cfc = mpm_.cellsFromCellsMap();
+
+        forAll(cfc, cfcI)
+        {
+            // Get addressing
+            const labelList& mo = cfc[cfcI].masterObjects();
+
+            label celli = cfc[cfcI].index();
+
+            if (addr[celli].size())
+            {
+                FatalErrorInFunction
+                    << "Master cell " << celli
+                    << " mapped from cell cells " << mo
+                    << " already destination of mapping."
+                    << abort(FatalError);
+            }
+
+            // Map from masters
+            addr[celli] = mo;
+        }
+
+        if (mpm_.hasOldCellVolumes())
+        {
+            // Volume weighted
+
+            const scalarField& V = mpm_.oldCellVolumes();
+
+            if (V.size() != sizeBeforeMapping())
+            {
+                FatalErrorInFunction
+                    << "cellVolumes size " << V.size()
+                    << " is not the old number of cells " << sizeBeforeMapping()
+                    << ". Are your cellVolumes already mapped?"
+                    << " (new number of cells " << mpm_.cellMap().size() << ")"
+                    << abort(FatalError);
+            }
+
+            forAll(cfc, cfcI)
+            {
+                const labelList& mo = cfc[cfcI].masterObjects();
+
+                label celli = cfc[cfcI].index();
+
+                w[celli].setSize(mo.size());
+
+                if (mo.size())
+                {
+                    scalar sumV = 0;
+                    forAll(mo, ci)
+                    {
+                        w[celli][ci] = V[mo[ci]];
+                        sumV += V[mo[ci]];
+                    }
+                    if (sumV > vSmall)
+                    {
+                        forAll(mo, ci)
+                        {
+                            w[celli][ci] /= sumV;
+                        }
+                    }
+                    else
+                    {
+                        // Exception: zero volume. Use uniform mapping
+                        w[celli] = scalarList(mo.size(), 1.0/mo.size());
+                    }
+                }
+            }
+        }
+        else
+        {
+            // Uniform weighted
+
+            forAll(cfc, cfcI)
+            {
+                const labelList& mo = cfc[cfcI].masterObjects();
+
+                label celli = cfc[cfcI].index();
+
+                w[celli] = scalarList(mo.size(), 1.0/mo.size());
+            }
+        }
+
+
+        // Do mapped faces. Note that can already be set from cellsFromCells
+        // so check if addressing size still zero.
+
+        const labelList& cm = mpm_.cellMap();
+
+        forAll(cm, celli)
+        {
+            if (cm[celli] > -1 && addr[celli].empty())
+            {
+                // Mapped from a single cell
+                addr[celli] = labelList(1, cm[celli]);
+                w[celli] = scalarList(1, 1.0);
+            }
+        }
+
+        // Grab inserted points (for them the size of addressing is still zero)
+
+        insertedCellLabelsPtr_ = new labelList(mesh_.nCells());
+        labelList& insertedCells = *insertedCellLabelsPtr_;
+
+        label nInsertedCells = 0;
+
+        forAll(addr, celli)
+        {
+            if (addr[celli].empty())
+            {
+                // Mapped from a dummy cell
+                addr[celli] = labelList(1, label(0));
+                w[celli] = scalarList(1, 1.0);
+
+                insertedCells[nInsertedCells] = celli;
+                nInsertedCells++;
+            }
+        }
+
+        insertedCells.setSize(nInsertedCells);
+    }
+}
+
+
+void Foam::cellMapper::clearOut()
+{
+    deleteDemandDrivenData(directAddrPtr_);
+    deleteDemandDrivenData(interpolationAddrPtr_);
+    deleteDemandDrivenData(weightsPtr_);
+    deleteDemandDrivenData(insertedCellLabelsPtr_);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::cellMapper::cellMapper(const polyTopoChangeMap& mpm)
+:
+    mesh_(mpm.mesh()),
+    mpm_(mpm),
+    insertedCells_(true),
+    direct_(false),
+    directAddrPtr_(nullptr),
+    interpolationAddrPtr_(nullptr),
+    weightsPtr_(nullptr),
+    insertedCellLabelsPtr_(nullptr)
+{
+    // Check for possibility of direct mapping
+    if
+    (
+        mpm_.cellsFromPointsMap().empty()
+     && mpm_.cellsFromEdgesMap().empty()
+     && mpm_.cellsFromFacesMap().empty()
+     && mpm_.cellsFromCellsMap().empty()
+    )
+    {
+        direct_ = true;
+    }
+    else
+    {
+        direct_ = false;
+    }
+
+    // Check for inserted cells
+    if (direct_ && (mpm_.cellMap().empty() || min(mpm_.cellMap()) > -1))
+    {
+        insertedCells_ = false;
+    }
+    else
+    {
+        // Need to check all 3 lists to see if there are inserted cells
+        // with no owner
+
+        // Make a copy of the cell map, add the entried for cells from points,
+        // cells from edges and cells from faces and check for left-overs
+        labelList cm(mesh_.nCells(), -1);
+
+        const List<objectMap>& cfp = mpm_.cellsFromPointsMap();
+
+        forAll(cfp, cfpI)
+        {
+            cm[cfp[cfpI].index()] = 0;
+        }
+
+        const List<objectMap>& cfe = mpm_.cellsFromEdgesMap();
+
+        forAll(cfe, cfeI)
+        {
+            cm[cfe[cfeI].index()] = 0;
+        }
+
+        const List<objectMap>& cff = mpm_.cellsFromFacesMap();
+
+        forAll(cff, cffI)
+        {
+            cm[cff[cffI].index()] = 0;
+        }
+
+        const List<objectMap>& cfc = mpm_.cellsFromCellsMap();
+
+        forAll(cfc, cfcI)
+        {
+            cm[cfc[cfcI].index()] = 0;
+        }
+
+        if (min(cm) < 0)
+        {
+            insertedCells_ = true;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::cellMapper::~cellMapper()
+{
+    clearOut();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::label Foam::cellMapper::sizeBeforeMapping() const
+{
+    return mpm_.nOldCells();
+}
+
+
+const Foam::labelUList& Foam::cellMapper::directAddressing() const
+{
+    if (!direct())
+    {
+        FatalErrorInFunction
+            << "Requested direct addressing for an interpolative mapper."
+            << abort(FatalError);
+    }
+
+    if (!insertedObjects())
+    {
+        // No inserted cells.  Re-use cellMap
+        return mpm_.cellMap();
+    }
+    else
+    {
+        if (!directAddrPtr_)
+        {
+            calcAddressing();
+        }
+
+        return *directAddrPtr_;
+    }
+}
+
+
+const Foam::labelListList& Foam::cellMapper::addressing() const
+{
+    if (direct())
+    {
+        FatalErrorInFunction
+            << "Requested interpolative addressing for a direct mapper."
+            << abort(FatalError);
+    }
+
+    if (!interpolationAddrPtr_)
+    {
+        calcAddressing();
+    }
+
+    return *interpolationAddrPtr_;
+}
+
+
+const Foam::scalarListList& Foam::cellMapper::weights() const
+{
+    if (direct())
+    {
+        FatalErrorInFunction
+            << "Requested interpolative weights for a direct mapper."
+            << abort(FatalError);
+    }
+
+    if (!weightsPtr_)
+    {
+        calcAddressing();
+    }
+
+    return *weightsPtr_;
+}
+
+
+const Foam::labelList& Foam::cellMapper::insertedObjectLabels() const
+{
+    if (!insertedCellLabelsPtr_)
+    {
+        if (!insertedObjects())
+        {
+            // There are no inserted cells
+            insertedCellLabelsPtr_ = new labelList(0);
+        }
+        else
+        {
+            calcAddressing();
+        }
+    }
+
+    return *insertedCellLabelsPtr_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/cellMapper/cellMapper.H b/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/cellMapper/cellMapper.H
new file mode 100644
index 0000000000..ad57aa6a89
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/cellMapper/cellMapper.H
@@ -0,0 +1,165 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::cellMapper
+
+Description
+    This object provides mapping and fill-in information for cell data
+    between the two meshes after the topological change.  It is
+    constructed from polyTopoChangeMap.
+
+SourceFiles
+    cellMapper.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef cellMapper_H
+#define cellMapper_H
+
+#include "morphFieldMapper.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class polyMesh;
+class polyTopoChangeMap;
+
+/*---------------------------------------------------------------------------*\
+                         Class cellMapper Declaration
+\*---------------------------------------------------------------------------*/
+
+class cellMapper
+:
+    public morphFieldMapper
+{
+    // Private Data
+
+        //- Reference to polyMesh
+        const polyMesh& mesh_;
+
+        //- Reference to polyTopoChangeMap
+        const polyTopoChangeMap& mpm_;
+
+        //- Are there any inserted (unmapped) cells
+        bool insertedCells_;
+
+        //- Is the mapping direct
+        bool direct_;
+
+
+    // Demand-driven private data
+
+        //- Direct addressing (only one for of addressing is used)
+        mutable labelList* directAddrPtr_;
+
+        //- Interpolated addressing (only one for of addressing is used)
+        mutable labelListList* interpolationAddrPtr_;
+
+        //- Interpolation weights
+        mutable scalarListList* weightsPtr_;
+
+        //- Inserted cells
+        mutable labelList* insertedCellLabelsPtr_;
+
+
+    // Private Member Functions
+
+        //- Calculate addressing for mapping with inserted cells
+        void calcAddressing() const;
+
+        //- Clear out local storage
+        void clearOut();
+
+
+public:
+
+    // Static Data Members
+
+    // Constructors
+
+        //- Construct from polyTopoChangeMap
+        cellMapper(const polyTopoChangeMap& mpm);
+
+        //- Disallow default bitwise copy construction
+        cellMapper(const cellMapper&) = delete;
+
+
+    //- Destructor
+    virtual ~cellMapper();
+
+
+    // Member Functions
+
+        //- Return size before mapping
+        virtual label sizeBeforeMapping() const;
+
+        //- Is the mapping direct
+        virtual bool direct() const
+        {
+            return direct_;
+        }
+
+        virtual bool hasUnmapped() const
+        {
+            return insertedObjects();
+        }
+
+        //- Return direct addressing
+        virtual const labelUList& directAddressing() const;
+
+        //- Return interpolated addressing
+        virtual const labelListList& addressing() const;
+
+        //- Return interpolation weights
+        virtual const scalarListList& weights() const;
+
+        //- Are there any inserted cells
+        virtual bool insertedObjects() const
+        {
+            return insertedCells_;
+        }
+
+        //- Return list of inserted cells
+        const virtual labelList& insertedObjectLabels() const;
+
+
+    // Member Operators
+
+        //- Disallow default bitwise assignment
+        void operator=(const cellMapper&) = delete;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/faceMapper/faceMapper.C b/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/faceMapper/faceMapper.C
new file mode 100644
index 0000000000..06cce90a05
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/faceMapper/faceMapper.C
@@ -0,0 +1,401 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "faceMapper.H"
+#include "fieldTypes.H"
+#include "demandDrivenData.H"
+#include "polyMesh.H"
+#include "polyTopoChangeMap.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::faceMapper::calcAddressing() const
+{
+    if
+    (
+        directAddrPtr_
+     || interpolationAddrPtr_
+     || weightsPtr_
+     || insertedFaceLabelsPtr_
+    )
+    {
+        FatalErrorInFunction
+            << "Addressing already calculated."
+            << abort(FatalError);
+    }
+
+    if (direct())
+    {
+        // Direct addressing, no weights
+
+        directAddrPtr_ = new labelList(mpm_.faceMap());
+        labelList& directAddr = *directAddrPtr_;
+
+        // Reset the size of addressing list to contain only live faces
+        directAddr.setSize(mesh_.nFaces());
+
+        insertedFaceLabelsPtr_ = new labelList(mesh_.nFaces());
+        labelList& insertedFaces = *insertedFaceLabelsPtr_;
+
+        label nInsertedFaces = 0;
+
+        forAll(directAddr, facei)
+        {
+            if (directAddr[facei] < 0)
+            {
+                // Found inserted face
+                directAddr[facei] = 0;
+                insertedFaces[nInsertedFaces] = facei;
+                nInsertedFaces++;
+            }
+        }
+
+        insertedFaces.setSize(nInsertedFaces);
+    }
+    else
+    {
+        // Interpolative addressing
+
+        interpolationAddrPtr_ = new labelListList(mesh_.nFaces());
+        labelListList& addr = *interpolationAddrPtr_;
+
+        weightsPtr_ = new scalarListList(mesh_.nFaces());
+        scalarListList& w = *weightsPtr_;
+
+        const List<objectMap>& ffp = mpm_.facesFromPointsMap();
+
+        forAll(ffp, ffpI)
+        {
+            // Get addressing
+            const labelList& mo = ffp[ffpI].masterObjects();
+
+            label facei = ffp[ffpI].index();
+
+            if (addr[facei].size())
+            {
+                FatalErrorInFunction
+                    << "Master face " << facei
+                    << " mapped from point faces " << mo
+                    << " already destination of mapping." << abort(FatalError);
+            }
+
+            // Map from masters, uniform weights
+            addr[facei] = mo;
+            w[facei] = scalarList(mo.size(), 1.0/mo.size());
+        }
+
+        const List<objectMap>& ffe = mpm_.facesFromEdgesMap();
+
+        forAll(ffe, ffeI)
+        {
+            // Get addressing
+            const labelList& mo = ffe[ffeI].masterObjects();
+
+            label facei = ffe[ffeI].index();
+
+            if (addr[facei].size())
+            {
+                FatalErrorInFunction
+                    << "Master face " << facei
+                    << " mapped from edge faces " << mo
+                    << " already destination of mapping." << abort(FatalError);
+            }
+
+            // Map from masters, uniform weights
+            addr[facei] = mo;
+            w[facei] = scalarList(mo.size(), 1.0/mo.size());
+        }
+
+        const List<objectMap>& fff = mpm_.facesFromFacesMap();
+
+        forAll(fff, fffI)
+        {
+            // Get addressing
+            const labelList& mo = fff[fffI].masterObjects();
+
+            label facei = fff[fffI].index();
+
+            if (addr[facei].size())
+            {
+                FatalErrorInFunction
+                    << "Master face " << facei
+                    << " mapped from face faces " << mo
+                    << " already destination of mapping." << abort(FatalError);
+            }
+
+            // Map from masters, uniform weights
+            addr[facei] = mo;
+            w[facei] = scalarList(mo.size(), 1.0/mo.size());
+        }
+
+
+        // Do mapped faces. Note that can already be set from facesFromFaces
+        // so check if addressing size still zero.
+        const labelList& fm = mpm_.faceMap();
+
+        forAll(fm, facei)
+        {
+            if (fm[facei] > -1 && addr[facei].empty())
+            {
+                // Mapped from a single face
+                addr[facei] = labelList(1, fm[facei]);
+                w[facei] = scalarList(1, 1.0);
+            }
+        }
+
+
+        // Grab inserted faces (for them the size of addressing is still zero)
+
+        insertedFaceLabelsPtr_ = new labelList(mesh_.nFaces());
+        labelList& insertedFaces = *insertedFaceLabelsPtr_;
+
+        label nInsertedFaces = 0;
+
+        forAll(addr, facei)
+        {
+            if (addr[facei].empty())
+            {
+                // Mapped from a dummy face
+                addr[facei] = labelList(1, label(0));
+                w[facei] = scalarList(1, 1.0);
+
+                insertedFaces[nInsertedFaces] = facei;
+                nInsertedFaces++;
+            }
+        }
+
+        insertedFaces.setSize(nInsertedFaces);
+    }
+}
+
+
+void Foam::faceMapper::clearOut()
+{
+    deleteDemandDrivenData(directAddrPtr_);
+    deleteDemandDrivenData(interpolationAddrPtr_);
+    deleteDemandDrivenData(weightsPtr_);
+    deleteDemandDrivenData(insertedFaceLabelsPtr_);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::faceMapper::faceMapper(const polyTopoChangeMap& mpm)
+:
+    mesh_(mpm.mesh()),
+    mpm_(mpm),
+    insertedFaces_(true),
+    direct_(false),
+    directAddrPtr_(nullptr),
+    interpolationAddrPtr_(nullptr),
+    weightsPtr_(nullptr),
+    insertedFaceLabelsPtr_(nullptr)
+{
+    // Check for possibility of direct mapping
+    if
+    (
+        mpm_.facesFromPointsMap().empty()
+     && mpm_.facesFromEdgesMap().empty()
+     && mpm_.facesFromFacesMap().empty()
+    )
+    {
+        direct_ = true;
+    }
+    else
+    {
+        direct_ = false;
+    }
+
+    // Check for inserted faces
+    if (direct_ && (mpm_.faceMap().empty() || min(mpm_.faceMap()) > -1))
+    {
+        insertedFaces_ = false;
+    }
+    else
+    {
+        // Need to check all 3 lists to see if there are inserted faces
+        // with no owner
+
+        // Make a copy of the face map, add the entries for faces from points
+        // and faces from edges and check for left-overs
+        labelList fm(mesh_.nFaces(), -1);
+
+        const List<objectMap>& ffp = mpm_.facesFromPointsMap();
+
+        forAll(ffp, ffpI)
+        {
+            fm[ffp[ffpI].index()] = 0;
+        }
+
+        const List<objectMap>& ffe = mpm_.facesFromEdgesMap();
+
+        forAll(ffe, ffeI)
+        {
+            fm[ffe[ffeI].index()] = 0;
+        }
+
+        const List<objectMap>& fff = mpm_.facesFromFacesMap();
+
+        forAll(fff, fffI)
+        {
+            fm[fff[fffI].index()] = 0;
+        }
+
+        if (min(fm) < 0)
+        {
+            insertedFaces_ = true;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::faceMapper::~faceMapper()
+{
+    clearOut();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::label Foam::faceMapper::sizeBeforeMapping() const
+{
+    return mpm_.nOldFaces();
+}
+
+
+Foam::label Foam::faceMapper::internalSizeBeforeMapping() const
+{
+    return mpm_.nOldInternalFaces();
+}
+
+
+const Foam::labelUList& Foam::faceMapper::directAddressing() const
+{
+    if (!direct())
+    {
+        FatalErrorInFunction
+            << "Requested direct addressing for an interpolative mapper."
+            << abort(FatalError);
+    }
+
+    if (!insertedObjects())
+    {
+        // No inserted faces.  Re-use faceMap
+        return mpm_.faceMap();
+    }
+    else
+    {
+        if (!directAddrPtr_)
+        {
+            calcAddressing();
+        }
+
+        return *directAddrPtr_;
+    }
+}
+
+
+const Foam::labelListList& Foam::faceMapper::addressing() const
+{
+    if (direct())
+    {
+        FatalErrorInFunction
+            << "Requested interpolative addressing for a direct mapper."
+            << abort(FatalError);
+    }
+
+    if (!interpolationAddrPtr_)
+    {
+        calcAddressing();
+    }
+
+    return *interpolationAddrPtr_;
+}
+
+
+const Foam::scalarListList& Foam::faceMapper::weights() const
+{
+    if (direct())
+    {
+        FatalErrorInFunction
+            << "Requested interpolative weights for a direct mapper."
+            << abort(FatalError);
+    }
+
+    if (!weightsPtr_)
+    {
+        calcAddressing();
+    }
+
+    return *weightsPtr_;
+}
+
+
+const Foam::labelList& Foam::faceMapper::insertedObjectLabels() const
+{
+    if (!insertedFaceLabelsPtr_)
+    {
+        if (!insertedObjects())
+        {
+            // There are no inserted faces
+            insertedFaceLabelsPtr_ = new labelList(0);
+        }
+        else
+        {
+            calcAddressing();
+        }
+    }
+
+    return *insertedFaceLabelsPtr_;
+}
+
+
+const Foam::labelHashSet& Foam::faceMapper::flipFaceFlux() const
+{
+    return mpm_.flipFaceFlux();
+}
+
+
+Foam::label Foam::faceMapper::nOldInternalFaces() const
+{
+    return mpm_.nOldInternalFaces();
+}
+
+
+const Foam::labelList& Foam::faceMapper::oldPatchStarts() const
+{
+    return mpm_.oldPatchStarts();
+}
+
+
+const Foam::labelList& Foam::faceMapper::oldPatchSizes() const
+{
+    return mpm_.oldPatchSizes();
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/faceMapper/faceMapper.H b/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/faceMapper/faceMapper.H
new file mode 100644
index 0000000000..7ae2ce7877
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/faceMapper/faceMapper.H
@@ -0,0 +1,181 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::faceMapper
+
+Description
+    This object provides mapping and fill-in information for face data
+    between the two meshes after the topological change.  It is
+    constructed from polyTopoChangeMap.
+
+SourceFiles
+    faceMapper.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef faceMapper_H
+#define faceMapper_H
+
+#include "morphFieldMapper.H"
+#include "HashSet.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class polyMesh;
+class polyTopoChangeMap;
+
+/*---------------------------------------------------------------------------*\
+                         Class faceMapper Declaration
+\*---------------------------------------------------------------------------*/
+
+class faceMapper
+:
+    public morphFieldMapper
+{
+    // Private Data
+
+        //- Reference to polyMesh
+        const polyMesh& mesh_;
+
+        //- Reference to polyTopoChangeMap
+        const polyTopoChangeMap& mpm_;
+
+        //- Are there any inserted (unmapped) faces
+        bool insertedFaces_;
+
+        //- Is the mapping direct
+        bool direct_;
+
+
+    // Demand-driven private data
+
+        //- Direct addressing (only one for of addressing is used)
+        mutable labelList* directAddrPtr_;
+
+        //- Interpolated addressing (only one for of addressing is used)
+        mutable labelListList* interpolationAddrPtr_;
+
+        //- Interpolation weights
+        mutable scalarListList* weightsPtr_;
+
+        //- Inserted faces
+        mutable labelList* insertedFaceLabelsPtr_;
+
+
+    // Private Member Functions
+
+        //- Calculate addressing for mapping with inserted faces
+        void calcAddressing() const;
+
+        //- Clear out local storage
+        void clearOut();
+
+
+public:
+
+    // Static Data Members
+
+    // Constructors
+
+        //- Construct from polyTopoChangeMap
+        faceMapper(const polyTopoChangeMap& mpm);
+
+        //- Disallow default bitwise copy construction
+        faceMapper(const faceMapper&) = delete;
+
+
+    //- Destructor
+    virtual ~faceMapper();
+
+
+    // Member Functions
+
+        //- Return size of field before mapping
+        virtual label sizeBeforeMapping() const;
+
+        //- Return number of internal faces before mapping
+        virtual label internalSizeBeforeMapping() const;
+
+        //- Is the mapping direct
+        virtual bool direct() const
+        {
+            return direct_;
+        }
+
+        virtual bool hasUnmapped() const
+        {
+            return insertedObjects();
+        }
+
+        //- Return direct addressing
+        virtual const labelUList& directAddressing() const;
+
+        //- Return interpolated addressing
+        virtual const labelListList& addressing() const;
+
+        //- Return interpolation weights
+        virtual const scalarListList& weights() const;
+
+        //- Return flux flip map
+        virtual const labelHashSet& flipFaceFlux() const;
+
+        //- Return number of old internalFaces
+        virtual label nOldInternalFaces() const;
+
+        //- Return old patch starts
+        virtual const labelList& oldPatchStarts() const;
+
+        //- Return old patch sizes
+        virtual const labelList& oldPatchSizes() const;
+
+        //- Are there any inserted faces
+        virtual bool insertedObjects() const
+        {
+            return insertedFaces_;
+        }
+
+        //- Return list of inserted faces
+        virtual const labelList& insertedObjectLabels() const;
+
+
+    // Member Operators
+
+        //- Disallow default bitwise assignment
+        void operator=(const faceMapper&) = delete;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/mapAddedPolyMesh.C b/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/mapAddedPolyMesh.C
new file mode 100644
index 0000000000..2566669a5e
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/mapAddedPolyMesh.C
@@ -0,0 +1,81 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "mapAddedPolyMesh.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::mapAddedPolyMesh::mapAddedPolyMesh
+(
+    const label nOldPoints,
+    const label nOldFaces,
+    const label nOldCells,
+    const label nAddedPoints,
+    const label nAddedFaces,
+    const label nAddedCells,
+    const labelList& oldPointMap,
+    const labelList& oldFaceMap,
+    const labelList& oldCellMap,
+
+    const labelList& addedPointMap,
+    const labelList& addedFaceMap,
+    const labelList& addedCellMap,
+
+    const labelList& oldPatchMap,
+    const labelList& addedPatchMap,
+    const labelList& oldPatchSizes,
+    const labelList& oldPatchStarts
+)
+:
+    nOldPoints_(nOldPoints),
+    nOldFaces_(nOldFaces),
+    nOldCells_(nOldCells),
+    nAddedPoints_(nAddedPoints),
+    nAddedFaces_(nAddedFaces),
+    nAddedCells_(nAddedCells),
+
+    oldPointMap_(oldPointMap),
+    oldFaceMap_(oldFaceMap),
+    oldCellMap_(oldCellMap),
+
+    addedPointMap_(addedPointMap),
+    addedFaceMap_(addedFaceMap),
+    addedCellMap_(addedCellMap),
+
+    oldPatchMap_(oldPatchMap),
+    addedPatchMap_(addedPatchMap),
+
+    oldPatchSizes_(oldPatchSizes),
+    oldPatchStarts_(oldPatchStarts)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/mapAddedPolyMesh.H b/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/mapAddedPolyMesh.H
new file mode 100644
index 0000000000..5785fb6df4
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/mapAddedPolyMesh.H
@@ -0,0 +1,245 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::mapAddedPolyMesh
+
+Description
+    Class containing mesh-to-mesh mapping information after a mesh addition
+    where we add a mesh ('added mesh') to an old mesh, creating a new mesh.
+
+    We store mapping from the old to the new mesh and from the added mesh
+    to the new mesh.
+
+    Note: Might need some more access functions or maybe some zone maps?
+
+SourceFiles
+    mapAddedPolyMesh.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef mapAddedPolyMesh_H
+#define mapAddedPolyMesh_H
+
+#include "labelList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class polyTopoChangeMap;
+
+/*---------------------------------------------------------------------------*\
+                      Class mapAddedPolyMesh Declaration
+\*---------------------------------------------------------------------------*/
+
+class mapAddedPolyMesh
+{
+    // Private Data
+
+        //- Old mesh points/face/cells
+        label nOldPoints_;
+        label nOldFaces_;
+        label nOldCells_;
+
+        //- Added mesh points/faces/cells
+        label nAddedPoints_;
+        label nAddedFaces_;
+        label nAddedCells_;
+
+
+        //- From old mesh points to new points
+        labelList oldPointMap_;
+        //- From old mesh faces to new faces
+        labelList oldFaceMap_;
+        //- From old mesh cells to new cells
+        labelList oldCellMap_;
+
+        //- From added mesh points to new points
+        labelList addedPointMap_;
+        //- From added mesh faces to new faces
+        labelList addedFaceMap_;
+        //- From added mesh cells to new cells
+        labelList addedCellMap_;
+
+        //- Original mesh to new mesh patch map. -1 for deleted patches.
+        labelList oldPatchMap_;
+
+        //- Added mesh to new mesh patch map. -1 for deleted patches.
+        labelList addedPatchMap_;
+
+        //- Original patch sizes on old mesh
+        labelList oldPatchSizes_;
+
+        //- Original patch starts
+        labelList oldPatchStarts_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        mapAddedPolyMesh
+        (
+            const label nOldPoints,
+            const label nOldFaces,
+            const label nOldCells,
+            const label nAddedPoints,
+            const label nAddedFaces,
+            const label nAddedCells,
+            const labelList& oldPointMap,
+            const labelList& oldFaceMap,
+            const labelList& oldCellMap,
+
+            const labelList& addedPointMap,
+            const labelList& addedFaceMap,
+            const labelList& addedCellMap,
+
+            const labelList& oldPatchMap,
+            const labelList& addedPatchMap,
+            const labelList& oldPatchSizes,
+            const labelList& oldPatchStarts
+        );
+
+
+    // Member Functions
+
+        // Access
+
+            // Old mesh data
+
+                label nOldPoints() const
+                {
+                    return nOldPoints_;
+                }
+
+                label nOldFaces() const
+                {
+                    return nOldFaces_;
+                }
+
+                label nOldCells() const
+                {
+                    return nOldCells_;
+                }
+
+
+                //- From old mesh point/face/cell to new mesh point/face/cell.
+                const labelList& oldPointMap() const
+                {
+                    return oldPointMap_;
+                }
+                const labelList& oldFaceMap() const
+                {
+                    return oldFaceMap_;
+                }
+                const labelList& oldCellMap() const
+                {
+                    return oldCellMap_;
+                }
+
+                //- From old patch index to new patch index or -1 if patch
+                //  not present (since 0 size)
+                const labelList& oldPatchMap() const
+                {
+                    return oldPatchMap_;
+                }
+
+                //- Return list of the old patch sizes
+                const labelList& oldPatchSizes() const
+                {
+                    return oldPatchSizes_;
+                }
+
+                //- Return list of the old patch start labels
+                const labelList& oldPatchStarts() const
+                {
+                    return oldPatchStarts_;
+                }
+
+                //- Number of old internal faces
+                label nOldInternalFaces() const
+                {
+                    return oldPatchStarts_[0];
+                }
+
+
+            // Added mesh data
+
+                label nAddedPoints() const
+                {
+                    return nAddedPoints_;
+                }
+
+                label nAddedFaces() const
+                {
+                    return nAddedFaces_;
+                }
+
+                label nAddedCells() const
+                {
+                    return nAddedCells_;
+                }
+
+                //- From added mesh point/face/cell to new mesh point/face/cell.
+                const labelList& addedPointMap() const
+                {
+                    return addedPointMap_;
+                }
+                const labelList& addedFaceMap() const
+                {
+                    return addedFaceMap_;
+                }
+                const labelList& addedCellMap() const
+                {
+                    return addedCellMap_;
+                }
+
+                //- From added mesh patch index to new patch index or -1 if
+                //  patch not present (since 0 size)
+                const labelList& addedPatchMap() const
+                {
+                    return addedPatchMap_;
+                }
+
+
+        // Edit
+
+            void updateMesh(const polyTopoChangeMap&)
+            {
+                NotImplemented;
+            }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/mapPatchChange/mapPatchChange.H b/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/mapPatchChange/mapPatchChange.H
new file mode 100644
index 0000000000..7c38d6ef47
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/mapPatchChange/mapPatchChange.H
@@ -0,0 +1,152 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::mapPatchChange
+
+Description
+    Class containing mesh-to-mesh mapping information after a patch change
+    operation.
+
+SourceFiles
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef mapPatchChange_H
+#define mapPatchChange_H
+
+#include "labelList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class mapPatchChange Declaration
+\*---------------------------------------------------------------------------*/
+
+class mapPatchChange
+{
+    // Private Data
+
+        //- Old patches
+        const label nOldPatches_;
+
+        //- Patch mapping array
+        const labelList patchMap_;
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        mapPatchChange(const label nOldPatches, const labelList& patchMap)
+        :
+            nOldPatches_(nOldPatches),
+            patchMap_(patchMap)
+        {}
+
+
+    // Member Functions
+
+        // Access
+
+            //- Number of old patches
+            label nOldPatches() const
+            {
+                return nOldPatches_;
+            }
+
+            //- Patch map. Size of current patches.
+            //  -1  : patch was added
+            //  >=0 : old position of patch
+            //  any original patch which is not in the list has been deleted
+            const labelList& patchMap() const
+            {
+                return patchMap_;
+            }
+
+
+        // Utility functions
+
+            //- Labels of added patches
+            labelList addedPatches() const
+            {
+                labelList added(patchMap_.size());
+
+                label addedI = 0;
+
+                forAll(patchMap_, patchi)
+                {
+                    if (patchMap_[patchi] == -1)
+                    {
+                        added[addedI++] = patchi;
+                    }
+                }
+                added.setSize(addedI);
+                return added;
+            }
+
+            //- Labels (on old mesh) of deleted patches
+            labelList deletedPatches() const
+            {
+                labelList oldToNew(nOldPatches_, -1);
+
+                // Mark all preserved patches
+                forAll(patchMap_, patchi)
+                {
+                    if (patchMap_[patchi] != -1)
+                    {
+                        oldToNew[patchMap_[patchi]] = patchi;
+                    }
+                }
+
+                // Extract -1 elements from oldToNew. These are the deleted
+                // patches.
+                label deletedI = 0;
+
+                forAll(oldToNew, oldPatchi)
+                {
+                    if (oldToNew[oldPatchi] == -1)
+                    {
+                        oldToNew[deletedI++] = oldPatchi;
+                    }
+                }
+
+                oldToNew.setSize(deletedI);
+
+                return oldToNew;
+            }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/mapSubsetMesh/mapSubsetMesh.H b/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/mapSubsetMesh/mapSubsetMesh.H
new file mode 100644
index 0000000000..808ff34724
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/mapSubsetMesh/mapSubsetMesh.H
@@ -0,0 +1,202 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::mapSubsetMesh
+
+Description
+    Class containing mesh-to-mesh mapping information after a subset operation
+
+SourceFiles
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef mapSubsetMesh_H
+#define mapSubsetMesh_H
+
+#include "labelList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class mapSubsetMesh Declaration
+\*---------------------------------------------------------------------------*/
+
+class mapSubsetMesh
+{
+    // Private Data
+
+        //- Number of old live points
+        label nOldPoints_;
+
+        //- Number of old live faces
+        label nOldFaces_;
+
+        //- Number of old live cells
+        label nOldCells_;
+
+        //- Point mapping array
+        const labelList pointMap_;
+
+        //- Face mapping array
+        const labelList faceMap_;
+
+        //- Cell mapping array
+        const labelList cellMap_;
+
+        //- Patch for exposed faces
+        const label exposedPatchID_;
+
+        //- List of the old patch sizes
+        labelList oldPatchSizes_;
+
+        //- List of the old patch start labels
+        const labelList oldPatchStarts_;
+
+        //- List of numbers of mesh points per old patch
+        const labelList oldPatchNMeshPoints_;
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        mapSubsetMesh
+        (
+            const label nOldPoints,
+            const label nOldFaces,
+            const label nOldCells,
+            const labelList& pointMap,
+            const labelList& faceMap,
+            const labelList& cellMap,
+            const label exposedPatchID,
+            const labelList& oldPatchStarts,
+            const labelList& oldPatchNMeshPoints
+        );
+
+        //- Construct from components and optionally reuse storage
+        mapSubsetMesh
+        (
+            const label nOldPoints,
+            const label nOldFaces,
+            const label nOldCells,
+            labelList& pointMap,
+            labelList& faceMap,
+            labelList& cellMap,
+            const label exposedPatchID,
+            labelList& oldPatchStarts,
+            labelList& oldPatchNMeshPoints,
+            const bool reuse
+        );
+
+
+    // Member Functions
+
+        // Access
+
+            //- Number of old points
+            label nOldPoints() const
+            {
+                return nOldPoints_;
+            }
+
+            //- Number of old internal faces
+            label nOldInternalFaces() const
+            {
+                return oldPatchStarts_[0];
+            }
+
+            //- Number of old faces
+            label nOldFaces() const
+            {
+                return nOldFaces_;
+            }
+
+            //- Number of old cells
+            label nOldCells() const
+            {
+                return nOldCells_;
+            }
+
+            //- Patch that exposed faces were put into
+            label exposedPatchID() const
+            {
+                return exposedPatchID_;
+            }
+
+            //- Old point map.
+            //  Contains the old point label for all points of the subsetted
+            //  mesh
+            const labelList& pointMap() const
+            {
+                return pointMap_;
+            }
+
+            //- Old face map.
+            //  Contains the old point label for all faces of the subsetted
+            //  mesh
+            const labelList& faceMap() const
+            {
+                return faceMap_;
+            }
+
+            //- Old cell map.
+            //  Contains the old point label for all cells of the subsetted
+            //  mesh
+            const labelList& cellMap() const
+            {
+                return cellMap_;
+            }
+
+            //- Return list of the old patch sizes
+            const labelList& oldPatchSizes() const
+            {
+                return oldPatchSizes_;
+            }
+
+            //- Return list of the old patch start labels
+            const labelList& oldPatchStarts() const
+            {
+                return oldPatchStarts_;
+            }
+
+            //- Return numbers of mesh points per old patch
+            const labelList& oldPatchNMeshPoints() const
+            {
+                return oldPatchNMeshPoints_;
+            }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/morphFieldMapper.H b/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/morphFieldMapper.H
new file mode 100644
index 0000000000..4510841617
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/morphFieldMapper.H
@@ -0,0 +1,87 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::morphFieldMapper
+
+Description
+    Abstract base class to hold the Field mapping for mesh morphs.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef morphFieldMapper_H
+#define morphFieldMapper_H
+
+#include "generalFieldMapper.H"
+#include "Map.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class morphFieldMapper Declaration
+\*---------------------------------------------------------------------------*/
+
+class morphFieldMapper
+:
+    public generalFieldMapper
+{
+
+public:
+
+    // Constructors
+
+        //- Null constructor
+        morphFieldMapper()
+        {}
+
+
+    //- Destructor
+    virtual ~morphFieldMapper()
+    {}
+
+
+    // Member Functions
+
+        //- Return size of field before mapping
+        virtual label sizeBeforeMapping() const = 0;
+
+        //- Are there any inserted objects
+        virtual bool insertedObjects() const = 0;
+
+        //- Return list of inserted objects
+        virtual const labelList& insertedObjectLabels() const = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/objectMap/objectMap.H b/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/objectMap/objectMap.H
new file mode 100644
index 0000000000..1d265517a5
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/objectMap/objectMap.H
@@ -0,0 +1,121 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::objectMap
+
+Description
+    An objectMap is a pair of labels defining the mapping of an object from
+    another object, e.g. a cell mapped from a point.
+
+SourceFiles
+    objectMapI.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef objectMap_H
+#define objectMap_H
+
+#include "labelList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of friend functions and operators
+
+class objectMap;
+inline bool operator==(const objectMap& a, const objectMap& b);
+inline bool operator!=(const objectMap& a, const objectMap& b);
+inline Ostream& operator<<(Ostream&, const objectMap&);
+inline Istream& operator>>(Istream&, objectMap&);
+
+
+/*---------------------------------------------------------------------------*\
+                          Class objectMap Declaration
+\*---------------------------------------------------------------------------*/
+
+class objectMap
+{
+    // Private Data
+
+        //- Object index
+        label index_;
+
+        //- Master object index
+        labelList masterObjects_;
+
+
+public:
+
+    // Constructors
+
+        //- Null constructor for lists
+        inline objectMap();
+
+        //- Construct from components
+        inline objectMap(const label index, const labelList& master);
+
+        //- Construct from Istream
+        inline objectMap(Istream&);
+
+
+    // Member Functions
+
+        //- Return object index
+        inline label& index();
+        inline label index() const;
+
+        //- Return master object index
+        inline labelList& masterObjects();
+        inline const labelList& masterObjects() const;
+
+
+    // Friend Operators
+
+        friend bool operator==(const objectMap& a, const objectMap& b);
+        friend bool operator!=(const objectMap& a, const objectMap& b);
+
+
+    // IOstream Operators
+
+        friend Ostream& operator<<(Ostream&, const objectMap&);
+
+        friend Istream& operator>>(Istream&, objectMap&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "objectMapI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/objectMap/objectMapI.H b/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/objectMap/objectMapI.H
new file mode 100644
index 0000000000..8b55ccb627
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/objectMap/objectMapI.H
@@ -0,0 +1,142 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Description
+
+\*---------------------------------------------------------------------------*/
+
+#include "IOstreams.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+inline objectMap::objectMap()
+:
+    index_(-1),
+    masterObjects_(0)
+{}
+
+
+inline objectMap::objectMap(const label index, const labelList& master)
+:
+    index_(index),
+    masterObjects_(master)
+{}
+
+
+inline objectMap::objectMap(Istream& is)
+{
+    // Read beginning of objectMap
+    is.readBegin("objectMap");
+
+    is >> index_ >> static_cast<labelList&>(masterObjects_);
+
+    // Read master of objectMap
+    is.readEnd("objectMap");
+
+    // Check state of Istream
+    is.check("objectMap::objectMap(Istream&)");
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+label& objectMap::index()
+{
+    return index_;
+}
+
+
+inline label objectMap::index() const
+{
+    return index_;
+}
+
+
+inline labelList& objectMap::masterObjects()
+{
+    return masterObjects_;
+}
+
+
+inline const labelList& objectMap::masterObjects() const
+{
+    return masterObjects_;
+}
+
+
+// * * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * //
+
+inline bool operator==(const objectMap& a, const objectMap& b)
+{
+    return
+    (
+        (a.index_ == b.index_) && (a.masterObjects_ == b.masterObjects_)
+    );
+}
+
+
+inline bool operator!=(const objectMap& a, const objectMap& b)
+{
+    return (!(a == b));
+}
+
+
+// * * * * * * * * * * * * * * * Ostream Operator *  * * * * * * * * * * * * //
+
+inline Ostream& operator<<(Ostream& os, const objectMap& a)
+{
+    os  << token::BEGIN_LIST
+        << a.index_ << token::SPACE
+        << a.masterObjects_
+        << token::END_LIST;
+
+    // Check state of Ostream
+    os.check("Ostream& operator<<(Ostream&, const objectMap&)");
+
+    return os;
+}
+
+
+inline Istream& operator>>(Istream& is, objectMap& a)
+{
+    is.readBegin("objectMap");
+    is  >> a.index_ >> a.masterObjects_;
+    is.readEnd("objectMap");
+
+    // Check state of Istream
+    is.check("Istream& operator>>(Istream&, objectMap&)");
+
+    return is;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // Master namespace Foam
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/polyTopoChangeMap.C b/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/polyTopoChangeMap.C
new file mode 100644
index 0000000000..0dc081520d
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/polyTopoChangeMap.C
@@ -0,0 +1,210 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "polyTopoChangeMap.H"
+#include "polyMesh.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::polyTopoChangeMap::polyTopoChangeMap
+(
+    const polyMesh& mesh,
+    const label nOldPoints,
+    const label nOldFaces,
+    const label nOldCells,
+    const labelList& pointMap,
+    const List<objectMap>& pointsFromPoints,
+    const labelList& faceMap,
+    const List<objectMap>& facesFromPoints,
+    const List<objectMap>& facesFromEdges,
+    const List<objectMap>& facesFromFaces,
+    const labelList& cellMap,
+    const List<objectMap>& cellsFromPoints,
+    const List<objectMap>& cellsFromEdges,
+    const List<objectMap>& cellsFromFaces,
+    const List<objectMap>& cellsFromCells,
+    const labelList& reversePointMap,
+    const labelList& reverseFaceMap,
+    const labelList& reverseCellMap,
+    const labelHashSet& flipFaceFlux,
+    const labelListList& patchPointMap,
+    const labelListList& pointZoneMap,
+    const labelListList& faceZonePointMap,
+    const labelListList& faceZoneFaceMap,
+    const labelListList& cellZoneMap,
+    const pointField& preMotionPoints,
+    const labelList& oldPatchStarts,
+    const labelList& oldPatchNMeshPoints,
+    const autoPtr<scalarField>& oldCellVolumesPtr
+)
+:
+    mesh_(mesh),
+    nOldPoints_(nOldPoints),
+    nOldFaces_(nOldFaces),
+    nOldCells_(nOldCells),
+    pointMap_(pointMap),
+    pointsFromPointsMap_(pointsFromPoints),
+    faceMap_(faceMap),
+    facesFromPointsMap_(facesFromPoints),
+    facesFromEdgesMap_(facesFromEdges),
+    facesFromFacesMap_(facesFromFaces),
+    cellMap_(cellMap),
+    cellsFromPointsMap_(cellsFromPoints),
+    cellsFromEdgesMap_(cellsFromEdges),
+    cellsFromFacesMap_(cellsFromFaces),
+    cellsFromCellsMap_(cellsFromCells),
+    reversePointMap_(reversePointMap),
+    reverseFaceMap_(reverseFaceMap),
+    reverseCellMap_(reverseCellMap),
+    flipFaceFlux_(flipFaceFlux),
+    patchPointMap_(patchPointMap),
+    pointZoneMap_(pointZoneMap),
+    faceZonePointMap_(faceZonePointMap),
+    faceZoneFaceMap_(faceZoneFaceMap),
+    cellZoneMap_(cellZoneMap),
+    preMotionPoints_(preMotionPoints),
+    oldPatchSizes_(oldPatchStarts.size()),
+    oldPatchStarts_(oldPatchStarts),
+    oldPatchNMeshPoints_(oldPatchNMeshPoints),
+    oldCellVolumesPtr_(oldCellVolumesPtr)
+{
+    if (oldPatchStarts_.size())
+    {
+        // Calculate old patch sizes
+        for (label patchi = 0; patchi < oldPatchStarts_.size() - 1; patchi++)
+        {
+            oldPatchSizes_[patchi] =
+                oldPatchStarts_[patchi + 1] - oldPatchStarts_[patchi];
+        }
+
+        // Set the last one by hand
+        const label lastPatchID = oldPatchStarts_.size() - 1;
+
+        oldPatchSizes_[lastPatchID] = nOldFaces_ - oldPatchStarts_[lastPatchID];
+
+        if (polyMesh::debug)
+        {
+            if (min(oldPatchSizes_) < 0)
+            {
+                FatalErrorInFunction
+                    << abort(FatalError);
+            }
+        }
+    }
+}
+
+
+Foam::polyTopoChangeMap::polyTopoChangeMap
+(
+    const polyMesh& mesh,
+    const label nOldPoints,
+    const label nOldFaces,
+    const label nOldCells,
+    labelList& pointMap,
+    List<objectMap>& pointsFromPoints,
+    labelList& faceMap,
+    List<objectMap>& facesFromPoints,
+    List<objectMap>& facesFromEdges,
+    List<objectMap>& facesFromFaces,
+    labelList& cellMap,
+    List<objectMap>& cellsFromPoints,
+    List<objectMap>& cellsFromEdges,
+    List<objectMap>& cellsFromFaces,
+    List<objectMap>& cellsFromCells,
+    labelList& reversePointMap,
+    labelList& reverseFaceMap,
+    labelList& reverseCellMap,
+    labelHashSet& flipFaceFlux,
+    labelListList& patchPointMap,
+    labelListList& pointZoneMap,
+    labelListList& faceZonePointMap,
+    labelListList& faceZoneFaceMap,
+    labelListList& cellZoneMap,
+    pointField& preMotionPoints,
+    labelList& oldPatchStarts,
+    labelList& oldPatchNMeshPoints,
+    autoPtr<scalarField>& oldCellVolumesPtr,
+    const bool reuse
+)
+:
+    mesh_(mesh),
+    nOldPoints_(nOldPoints),
+    nOldFaces_(nOldFaces),
+    nOldCells_(nOldCells),
+    pointMap_(pointMap, reuse),
+    pointsFromPointsMap_(pointsFromPoints, reuse),
+    faceMap_(faceMap, reuse),
+    facesFromPointsMap_(facesFromPoints, reuse),
+    facesFromEdgesMap_(facesFromEdges, reuse),
+    facesFromFacesMap_(facesFromFaces, reuse),
+    cellMap_(cellMap, reuse),
+    cellsFromPointsMap_(cellsFromPoints, reuse),
+    cellsFromEdgesMap_(cellsFromEdges, reuse),
+    cellsFromFacesMap_(cellsFromFaces, reuse),
+    cellsFromCellsMap_(cellsFromCells, reuse),
+    reversePointMap_(reversePointMap, reuse),
+    reverseFaceMap_(reverseFaceMap, reuse),
+    reverseCellMap_(reverseCellMap, reuse),
+    flipFaceFlux_(flipFaceFlux),
+    patchPointMap_(patchPointMap, reuse),
+    pointZoneMap_(pointZoneMap, reuse),
+    faceZonePointMap_(faceZonePointMap, reuse),
+    faceZoneFaceMap_(faceZoneFaceMap, reuse),
+    cellZoneMap_(cellZoneMap, reuse),
+    preMotionPoints_(preMotionPoints, reuse),
+    oldPatchSizes_(oldPatchStarts.size()),
+    oldPatchStarts_(oldPatchStarts, reuse),
+    oldPatchNMeshPoints_(oldPatchNMeshPoints, reuse),
+    oldCellVolumesPtr_(oldCellVolumesPtr, reuse)
+{
+    if (oldPatchStarts_.size() > 0)
+    {
+        // Calculate old patch sizes
+        for (label patchi = 0; patchi < oldPatchStarts_.size() - 1; patchi++)
+        {
+            oldPatchSizes_[patchi] =
+                oldPatchStarts_[patchi + 1] - oldPatchStarts_[patchi];
+        }
+
+        // Set the last one by hand
+        const label lastPatchID = oldPatchStarts_.size() - 1;
+
+        oldPatchSizes_[lastPatchID] = nOldFaces_ - oldPatchStarts_[lastPatchID];
+
+        if (polyMesh::debug)
+        {
+            if (min(oldPatchSizes_) < 0)
+            {
+                FatalErrorInFunction
+                    << "Calculated negative old patch size."
+                    << "  Error in mapping data"
+                    << abort(FatalError);
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/polyTopoChangeMap.H b/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/polyTopoChangeMap.H
new file mode 100644
index 0000000000..e98cd0b9af
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyTopoChangeMap/polyTopoChangeMap.H
@@ -0,0 +1,662 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::polyTopoChangeMap
+
+Description
+    Class containing mesh-to-mesh mapping information after a change
+    in polyMesh topology.
+
+    General:
+        - pointMap/faceMap/cellMap: \n
+          from current mesh back to previous mesh.
+          (so to 'pull' the information onto the current mesh)
+        - reversePointMap/faceMap/cellMap: \n
+          from previous mesh to current. (so to 'push' information)
+
+    In the topology change points/faces/cells
+    - can be unchanged. (faces might be renumbered though)
+    - can be removed (into nothing)
+    - can be removed into/merged with existing same entity
+      (so point merged with other point, face with other face, cell with
+       other cell. Note that probably only cell with cell is relevant)
+    - can be added from existing same 'master' entity
+      (so point from point, face from face and cell from cell)
+    - can be inflated: face out of edge or point,
+        cell out of face, edge or point.
+    - can be appended: added 'out of nothing'.
+
+    All this information is necessary to correctly map fields.
+
+    \par points
+
+    - unchanged:
+        - pointMap[pointi] contains old point label
+        - reversePointMap[oldPointi] contains new point label
+    - removed:
+        - reversePointMap[oldPointi] contains -1
+    - merged into point:
+        - reversePointMap[oldPointi] contains <-1 : -newPointi-2
+        - pointMap[pointi] contains the old master point label
+        - pointsFromPoints gives for pointi all the old point labels
+          (including the old master point!)
+    - added-from-same:
+        - pointMap[pointi] contains the old master point label
+    - appended:
+        - pointMap[pointi] contains -1
+
+    \par faces
+
+    - unchanged:
+        - faceMap[facei] contains old face label
+        - reverseFaceMap[oldFacei] contains new face label
+    - removed:
+        - reverseFaceMap[oldFacei] contains -1
+    - merged into face:
+        - reverseFaceMap[oldFacei] contains <-1 : -newFacei-2
+        - faceMap[facei] contains the old master face label
+        - facesFromFaces gives for facei all the old face labels
+          (including the old master face!)
+    - added-from-same:
+        - faceMap[facei] contains the old master face label
+    - inflated-from-edge:
+        - faceMap[facei] contains -1
+        - facesFromEdges contains an entry with
+            - facei
+            - list of faces(*) on old mesh that connected to the old edge
+    - inflated-from-point:
+        - faceMap[facei] contains -1
+        - facesFromPoints contains an entry with
+            - facei
+            - list of faces(*) on old mesh that connected to the old point
+    - appended:
+        - faceMap[facei] contains -1
+
+    Note (*) \n
+       if the newly inflated face is a boundary face the list of faces will
+       only be boundary faces; if the new face is an internal face they
+       will only be internal faces.
+
+    \par cells
+
+    - unchanged:
+        - cellMap[celli] contains old cell label
+        - reverseCellMap[oldCelli] contains new cell label
+    - removed:
+        - reverseCellMap[oldCelli] contains -1
+    - merged into cell:
+        - reverseCellMap[oldCelli] contains <-1 : -newCelli-2
+        - cellMap[celli] contains the old master cell label
+        - cellsFromCells gives for celli all the old cell labels
+          (including the old master cell!)
+    - added-from-same:
+        - cellMap[celli] contains the old master cell label
+    - inflated-from-face:
+        - cellMap[celli] contains -1
+        - cellsFromFaces contains an entry with
+            - celli
+            - list of cells on old mesh that connected to the old face
+    - inflated-from-edge:
+        - cellMap[celli] contains -1
+        - cellsFromEdges contains an entry with
+            - celli
+            - list of cells on old mesh that connected to the old edge
+    - inflated-from-point:
+        - cellMap[celli] contains -1
+        - cellsFromPoints contains an entry with
+            - celli
+            - list of cells on old mesh that connected to the old point
+    - appended:
+        - cellMap[celli] contains -1
+
+
+SourceFiles
+    polyTopoChangeMap.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef polyTopoChangeMap_H
+#define polyTopoChangeMap_H
+
+#include "labelList.H"
+#include "objectMap.H"
+#include "pointField.H"
+#include "HashSet.H"
+#include "Map.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class polyMesh;
+
+/*---------------------------------------------------------------------------*\
+                         Class polyTopoChangeMap Declaration
+\*---------------------------------------------------------------------------*/
+
+class polyTopoChangeMap
+{
+    // Private Data
+
+        //- Reference to polyMesh
+        const polyMesh& mesh_;
+
+        //- Number of old live points
+        const label nOldPoints_;
+
+        //- Number of old live faces
+        const label nOldFaces_;
+
+        //- Number of old live cells
+        const label nOldCells_;
+
+        //- Old point map.
+        //  Contains the old point label for all new points.
+        //  - for preserved points this is the old point label.
+        //  - for added points this is the master point ID
+        //  - for points added with no master, this is -1
+        //  Size of the list equals the size of new points
+        const labelList pointMap_;
+
+        //- Points resulting from merging points
+        const List<objectMap> pointsFromPointsMap_;
+
+        //- Old face map.
+        //  Contains a list of old face labels for every new face.
+        //  Size of the list equals the number of new faces
+        //  - for preserved faces this is the old face label.
+        //  - for faces added from faces this is the master face ID
+        //  - for faces added with no master, this is -1
+        //  - for faces added from points or edges, this is -1
+        const labelList faceMap_;
+
+        //- Faces inflated from points
+        const List<objectMap> facesFromPointsMap_;
+
+        //- Faces inflated from edges
+        const List<objectMap> facesFromEdgesMap_;
+
+        //- Faces resulting from merging faces
+        const List<objectMap> facesFromFacesMap_;
+
+        //- Old cell map.
+        //  Contains old cell label for all preserved cells.
+        //  Size of the list equals the number or preserved cells
+        const labelList cellMap_;
+
+        //- Cells inflated from points
+        const List<objectMap> cellsFromPointsMap_;
+
+        //- Cells inflated from edges
+        const List<objectMap> cellsFromEdgesMap_;
+
+        //- Cells inflated from faces
+        const List<objectMap> cellsFromFacesMap_;
+
+        //- Cells resulting from merging cells
+        const List<objectMap> cellsFromCellsMap_;
+
+        //- Reverse point map
+        const labelList reversePointMap_;
+
+        //- Reverse face map
+        const labelList reverseFaceMap_;
+
+        //- Reverse cell map
+        const labelList reverseCellMap_;
+
+        //- Map of flipped face flux faces
+        const labelHashSet flipFaceFlux_;
+
+        //- Patch mesh point renumbering
+        const labelListList patchPointMap_;
+
+        //- Point zone renumbering
+        //  For every preserved point in zone give the old position.
+        //  For added points, the index is set to -1
+        const labelListList pointZoneMap_;
+
+        //- Face zone point renumbering
+        //  For every preserved point in zone give the old position.
+        //  For added points, the index is set to -1
+        const labelListList faceZonePointMap_;
+
+        //- Face zone face renumbering
+        //  For every preserved face in zone give the old position.
+        //  For added faces, the index is set to -1
+        const labelListList faceZoneFaceMap_;
+
+        //- Cell zone renumbering
+        //  For every preserved cell in zone give the old position.
+        //  For added cells, the index is set to -1
+        const labelListList cellZoneMap_;
+
+        //- Pre-motion point positions.
+        //  This specifies the correct way of blowing up zero-volume objects
+        const pointField preMotionPoints_;
+
+        //- List of the old patch sizes
+        labelList oldPatchSizes_;
+
+        //- List of the old patch start labels
+        const labelList oldPatchStarts_;
+
+        //- List of numbers of mesh points per old patch
+        const labelList oldPatchNMeshPoints_;
+
+        //- Optional old cell volumes (for mapping)
+        autoPtr<scalarField> oldCellVolumesPtr_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components. Copy (except for oldCellVolumes).
+        polyTopoChangeMap
+        (
+            const polyMesh& mesh,
+            const label nOldPoints,
+            const label nOldFaces,
+            const label nOldCells,
+            const labelList& pointMap,
+            const List<objectMap>& pointsFromPoints,
+            const labelList& faceMap,
+            const List<objectMap>& facesFromPoints,
+            const List<objectMap>& facesFromEdges,
+            const List<objectMap>& facesFromFaces,
+            const labelList& cellMap,
+            const List<objectMap>& cellsFromPoints,
+            const List<objectMap>& cellsFromEdges,
+            const List<objectMap>& cellsFromFaces,
+            const List<objectMap>& cellsFromCells,
+            const labelList& reversePointMap,
+            const labelList& reverseFaceMap,
+            const labelList& reverseCellMap,
+            const labelHashSet& flipFaceFlux,
+            const labelListList& patchPointMap,
+            const labelListList& pointZoneMap,
+            const labelListList& faceZonePointMap,
+            const labelListList& faceZoneFaceMap,
+            const labelListList& cellZoneMap,
+            const pointField& preMotionPoints,
+            const labelList& oldPatchStarts,
+            const labelList& oldPatchNMeshPoints,
+            const autoPtr<scalarField>& oldCellVolumesPtr
+        );
+
+        //- Construct from components and optionally reuse storage
+        polyTopoChangeMap
+        (
+            const polyMesh& mesh,
+            const label nOldPoints,
+            const label nOldFaces,
+            const label nOldCells,
+            labelList& pointMap,
+            List<objectMap>& pointsFromPoints,
+            labelList& faceMap,
+            List<objectMap>& facesFromPoints,
+            List<objectMap>& facesFromEdges,
+            List<objectMap>& facesFromFaces,
+            labelList& cellMap,
+            List<objectMap>& cellsFromPoints,
+            List<objectMap>& cellsFromEdges,
+            List<objectMap>& cellsFromFaces,
+            List<objectMap>& cellsFromCells,
+            labelList& reversePointMap,
+            labelList& reverseFaceMap,
+            labelList& reverseCellMap,
+            labelHashSet& flipFaceFlux,
+            labelListList& patchPointMap,
+            labelListList& pointZoneMap,
+            labelListList& faceZonePointMap,
+            labelListList& faceZoneFaceMap,
+            labelListList& cellZoneMap,
+            pointField& preMotionPoints,
+            labelList& oldPatchStarts,
+            labelList& oldPatchNMeshPoints,
+            autoPtr<scalarField>& oldCellVolumesPtr,
+            const bool reuse
+        );
+
+        //- Disallow default bitwise copy construction
+        polyTopoChangeMap(const polyTopoChangeMap&) = delete;
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return polyMesh
+            const polyMesh& mesh() const
+            {
+                return mesh_;
+            }
+
+            //- Number of old points
+            label nOldPoints() const
+            {
+                return nOldPoints_;
+            }
+
+            //- Number of old internal faces
+            label nOldInternalFaces() const
+            {
+                return oldPatchStarts_[0];
+            }
+
+            //- Number of old faces
+            label nOldFaces() const
+            {
+                return nOldFaces_;
+            }
+
+            //- Number of old cells
+            label nOldCells() const
+            {
+                return nOldCells_;
+            }
+
+            //- Old point map.
+            //  Contains the old point label for all new points.
+            //  For preserved points this is the old point label.
+            //  For added points this is the master point ID
+            const labelList& pointMap() const
+            {
+                return pointMap_;
+            }
+
+            //- Points originating from points
+            const List<objectMap>& pointsFromPointsMap() const
+            {
+                return pointsFromPointsMap_;
+            }
+
+            //- Old face map.
+            //  Contains a list of old face labels for every new face.
+            //  Warning: this map contains invalid entries for new faces
+            const labelList& faceMap() const
+            {
+                return faceMap_;
+            }
+
+            //- Faces inflated from points
+            const List<objectMap>& facesFromPointsMap() const
+            {
+                return facesFromPointsMap_;
+            }
+
+            //- Faces inflated from edges
+            const List<objectMap>& facesFromEdgesMap() const
+            {
+                return facesFromEdgesMap_;
+            }
+
+            //- Faces originating from faces
+            const List<objectMap>& facesFromFacesMap() const
+            {
+                return facesFromFacesMap_;
+            }
+
+            //- Old cell map.
+            //  Contains old cell label for all preserved cells.
+            const labelList& cellMap() const
+            {
+                return cellMap_;
+            }
+
+            //- Cells inflated from points
+            const List<objectMap>& cellsFromPointsMap() const
+            {
+                return cellsFromPointsMap_;
+            }
+
+            //- Cells inflated from edges
+            const List<objectMap>& cellsFromEdgesMap() const
+            {
+                return cellsFromEdgesMap_;
+            }
+
+            //- Cells inflated from faces
+            const List<objectMap>& cellsFromFacesMap() const
+            {
+                return cellsFromFacesMap_;
+            }
+
+            //- Cells originating from cells
+            const List<objectMap>& cellsFromCellsMap() const
+            {
+                return cellsFromCellsMap_;
+            }
+
+
+            // Reverse maps
+
+                //- Reverse point map
+                //  Contains new point label for all old and added points
+                const labelList& reversePointMap() const
+                {
+                    return reversePointMap_;
+                }
+
+                //- If point is removed return point (on new mesh) it merged
+                //  into
+                label mergedPoint(const label oldPointi) const
+                {
+                    label i = reversePointMap_[oldPointi];
+
+                    if (i == -1)
+                    {
+                        return i;
+                    }
+                    else if (i < -1)
+                    {
+                        return -i-2;
+                    }
+                    else
+                    {
+                        FatalErrorInFunction
+                            << "old point label " << oldPointi
+                            << " has reverseMap " << i << endl
+                            << "Only call mergedPoint for removed points."
+                            << abort(FatalError);
+                        return -1;
+                    }
+                }
+
+                //- Reverse face map
+                //  Contains new face label for all old and added faces
+                const labelList& reverseFaceMap() const
+                {
+                    return reverseFaceMap_;
+                }
+
+                //- If face is removed return face (on new mesh) it merged into
+                label mergedFace(const label oldFacei) const
+                {
+                    label i = reverseFaceMap_[oldFacei];
+
+                    if (i == -1)
+                    {
+                        return i;
+                    }
+                    else if (i < -1)
+                    {
+                        return -i-2;
+                    }
+                    else
+                    {
+                        FatalErrorInFunction
+                            << "old face label " << oldFacei
+                            << " has reverseMap " << i << endl
+                            << "Only call mergedFace for removed faces."
+                            << abort(FatalError);
+                        return -1;
+                    }
+                }
+
+                //- Reverse cell map
+                //  Contains new cell label for all old and added cells
+                const labelList& reverseCellMap() const
+                {
+                    return reverseCellMap_;
+                }
+
+                //- If cell is removed return cell (on new mesh) it merged into
+                label mergedCell(const label oldCelli) const
+                {
+                    label i = reverseCellMap_[oldCelli];
+
+                    if (i == -1)
+                    {
+                        return i;
+                    }
+                    else if (i < -1)
+                    {
+                        return -i-2;
+                    }
+                    else
+                    {
+                        FatalErrorInFunction
+                            << "old cell label " << oldCelli
+                            << " has reverseMap " << i << endl
+                            << "Only call mergedCell for removed cells."
+                            << abort(FatalError);
+                        return -1;
+                    }
+                }
+
+                //- Map of flipped face flux faces
+                const labelHashSet& flipFaceFlux() const
+                {
+                    return flipFaceFlux_;
+                }
+
+                //- Patch point renumbering
+                //  For every preserved point on a patch give the old position.
+                //  For added points, the index is set to -1
+                const labelListList& patchPointMap() const
+                {
+                    return patchPointMap_;
+                }
+
+
+            // Zone mapping
+
+                //- Point zone renumbering
+                //  For every preserved point in zone give the old position.
+                //  For added points, the index is set to -1
+                const labelListList& pointZoneMap() const
+                {
+                    return pointZoneMap_;
+                }
+
+                //- Face zone point renumbering
+                //  For every preserved point in zone give the old position.
+                //  For added points, the index is set to -1
+                const labelListList& faceZonePointMap() const
+                {
+                    return faceZonePointMap_;
+                }
+
+                //- Face zone face renumbering
+                //  For every preserved face in zone give the old position.
+                //  For added faces, the index is set to -1
+                const labelListList& faceZoneFaceMap() const
+                {
+                    return faceZoneFaceMap_;
+                }
+
+                //- Cell zone renumbering
+                //  For every preserved cell in zone give the old position.
+                //  For added cells, the index is set to -1
+                const labelListList& cellZoneMap() const
+                {
+                    return cellZoneMap_;
+                }
+
+            //- Pre-motion point positions.
+            //  This specifies the correct way of blowing up
+            //  zero-volume objects
+            const pointField& preMotionPoints() const
+            {
+                return preMotionPoints_;
+            }
+
+            //- Has valid preMotionPoints?
+            bool hasMotionPoints() const
+            {
+                return preMotionPoints_.size() > 0;
+            }
+
+
+            //- Return list of the old patch sizes
+            const labelList& oldPatchSizes() const
+            {
+                return oldPatchSizes_;
+            }
+
+            //- Return list of the old patch start labels
+            const labelList& oldPatchStarts() const
+            {
+                return oldPatchStarts_;
+            }
+
+            //- Return numbers of mesh points per old patch
+            const labelList& oldPatchNMeshPoints() const
+            {
+                return oldPatchNMeshPoints_;
+            }
+
+
+            // Geometric mapping data
+
+                bool hasOldCellVolumes() const
+                {
+                    return oldCellVolumesPtr_.valid();
+                }
+
+                const scalarField& oldCellVolumes() const
+                {
+                    return oldCellVolumesPtr_();
+                }
+
+
+    // Member Operators
+
+        //- Disallow default bitwise assignment
+        void operator=(const polyTopoChangeMap&) = delete;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //