diff --git a/applications/utilities/finiteArea/makeFaMesh/Make/options b/applications/utilities/finiteArea/makeFaMesh/Make/options
index 98ba4f7225..3300bd23e7 100644
--- a/applications/utilities/finiteArea/makeFaMesh/Make/options
+++ b/applications/utilities/finiteArea/makeFaMesh/Make/options
@@ -1,8 +1,12 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteArea/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fileFormats/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/cfdTools/general/lnInclude
 
 EXE_LIBS = \
     -lfiniteArea \
-    -lfiniteVolume
+    -lfiniteVolume \
+    -lfileFormats \
+    -lmeshTools
diff --git a/applications/utilities/finiteArea/makeFaMesh/decomposeFaFields.H b/applications/utilities/finiteArea/makeFaMesh/decomposeFaFields.H
new file mode 100644
index 0000000000..23b0c2c1e6
--- /dev/null
+++ b/applications/utilities/finiteArea/makeFaMesh/decomposeFaFields.H
@@ -0,0 +1,23 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
+
+Description
+    placeholder for decomposing fields
+
+\*---------------------------------------------------------------------------*/
+
+if (Pstream::parRun())
+{
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/finiteArea/makeFaMesh/findMeshDefinitionDict.H b/applications/utilities/finiteArea/makeFaMesh/findMeshDefinitionDict.H
new file mode 100644
index 0000000000..8a678fdcc8
--- /dev/null
+++ b/applications/utilities/finiteArea/makeFaMesh/findMeshDefinitionDict.H
@@ -0,0 +1,97 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
+
+Description
+    Search for the appropriate faMeshDefinition dictionary....
+
+\*---------------------------------------------------------------------------*/
+
+const word dictName("faMeshDefinition");
+
+autoPtr<IOdictionary> meshDictPtr;
+
+{
+    fileName dictPath;
+
+    const word& regionDir =
+        (regionName == polyMesh::defaultRegion ? word::null : regionName);
+
+    if (args.readIfPresent("dict", dictPath))
+    {
+        // Dictionary specified on the command-line ...
+
+        if (isDir(dictPath))
+        {
+            dictPath /= dictName;
+        }
+    }
+    else if
+    (
+        // Check global location
+        exists
+        (
+            runTime.path()/runTime.caseConstant()
+            /regionDir/faMesh::meshSubDir/dictName
+        )
+    )
+    {
+        // Dictionary present in constant faMesh directory (old-style)
+
+        dictPath =
+            runTime.constant()
+            /regionDir/faMesh::meshSubDir/dictName;
+
+        // Warn that constant/faMesh/faMeshDefinition was used
+        // instead of system/faMeshDefinition
+        #if 0
+        WarningIn(args.executable())
+            << "Using the old faMeshDefinition location: "
+            << dictPath << nl
+            << "    instead of default location: "
+            << runTime.system()/regionDir/dictName << nl
+            << endl;
+        #endif
+    }
+    else
+    {
+        // Assume dictionary is in the system directory
+
+        dictPath = runTime.system()/regionDir/dictName;
+    }
+
+    IOobject meshDictIO
+    (
+        dictPath,
+        runTime,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE,
+        false,  // no registerObject
+        true    // is globalObject
+    );
+
+    if (!meshDictIO.typeHeaderOk<IOdictionary>(true))
+    {
+        FatalErrorInFunction
+            << meshDictIO.objectPath() << nl
+            << exit(FatalError);
+    }
+
+    Info<< "Creating faMesh from definition: "
+        << runTime.relativePath(meshDictIO.objectPath()) << endl;
+
+    meshDictPtr = autoPtr<IOdictionary>::New(meshDictIO);
+}
+
+IOdictionary& meshDefDict = *meshDictPtr;
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/finiteArea/makeFaMesh/makeFaMesh.C b/applications/utilities/finiteArea/makeFaMesh/makeFaMesh.C
index 0010be3dab..2751ca3154 100644
--- a/applications/utilities/finiteArea/makeFaMesh/makeFaMesh.C
+++ b/applications/utilities/finiteArea/makeFaMesh/makeFaMesh.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -41,30 +42,12 @@ Author
 #include "OSspecific.H"
 #include "faMesh.H"
 #include "fvMesh.H"
+#include "IOdictionary.H"
+#include "globalIndex.H"
+#include "globalMeshData.H"
 
 using namespace Foam;
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-class faPatchData
-{
-public:
-    word name_;
-    word type_;
-    dictionary dict_;
-    label ownPolyPatchID_;
-    label ngbPolyPatchID_;
-    labelList edgeLabels_;
-    faPatchData()
-    :
-        name_(word::null),
-        type_(word::null),
-        ownPolyPatchID_(-1),
-        ngbPolyPatchID_(-1)
-    {}
-};
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
@@ -73,286 +56,40 @@ int main(int argc, char *argv[])
     (
         "A mesh generator for finiteArea mesh"
     );
+    argList::addOption
+    (
+        "empty-patch",
+        "name",
+        "Specify name for a default empty patch",
+        false  // An advanced option, but not enough to worry about that
+    );
+    argList::addOption("dict", "file", "Alternative faMeshDefinition");
 
     #include "addRegionOption.H"
-    argList::noParallel();
-
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createNamedMesh.H"
 
     // Reading faMeshDefinition dictionary
-    IOdictionary faMeshDefinition
-    (
-        IOobject
-        (
-            "faMeshDefinition",
-            runTime.constant(),
-            "faMesh",
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    );
+    #include "findMeshDefinitionDict.H"
 
-    wordList polyMeshPatches
-    (
-        faMeshDefinition.get<wordList>("polyMeshPatches")
-    );
-
-    const dictionary& bndDict = faMeshDefinition.subDict("boundary");
-
-    const wordList faPatchNames(bndDict.toc());
-
-    List<faPatchData> faPatches(faPatchNames.size()+1);
-
-    forAll(faPatchNames, patchI)
+    // Inject/overwrite name for optional 'empty' patch
+    word patchName;
+    if (args.readIfPresent("empty-patch", patchName))
     {
-        const dictionary& curPatchDict = bndDict.subDict(faPatchNames[patchI]);
-
-        faPatches[patchI].name_ = faPatchNames[patchI];
-
-        faPatches[patchI].type_ = curPatchDict.get<word>("type");
-
-        const word ownName(curPatchDict.get<word>("ownerPolyPatch"));
-
-        faPatches[patchI].ownPolyPatchID_ =
-            mesh.boundaryMesh().findPatchID(ownName);
-
-        if (faPatches[patchI].ownPolyPatchID_ < 0)
-        {
-            FatalErrorIn("makeFaMesh:")
-                << "neighbourPolyPatch " << ownName << " does not exist"
-                << exit(FatalError);
-        }
-
-        const word neiName(curPatchDict.get<word>("neighbourPolyPatch"));
-
-        faPatches[patchI].ngbPolyPatchID_ =
-            mesh.boundaryMesh().findPatchID(neiName);
-
-        if (faPatches[patchI].ngbPolyPatchID_ < 0)
-        {
-            FatalErrorIn("makeFaMesh:")
-                << "neighbourPolyPatch " << neiName << " does not exist"
-                << exit(FatalError);
-        }
+        meshDefDict.add("emptyPatch", patchName, true);
     }
 
-    // Setting faceLabels list size
-    label size = 0;
-
-    labelList patchIDs(polyMeshPatches.size(), -1);
-
-    forAll(polyMeshPatches, patchI)
-    {
-        patchIDs[patchI] =
-            mesh.boundaryMesh().findPatchID(polyMeshPatches[patchI]);
-
-        if (patchIDs[patchI] < 0)
-        {
-            FatalErrorIn("makeFaMesh:")
-                << "Patch " << polyMeshPatches[patchI] << " does not exist"
-                << exit(FatalError);
-        }
-
-        size += mesh.boundaryMesh()[patchIDs[patchI]].size();
-    }
-
-    labelList faceLabels(size, -1);
-
-    sort(patchIDs);
-
-
-    // Filling of faceLabels list
-    label faceI = -1;
-
-    forAll(polyMeshPatches, patchI)
-    {
-        label start = mesh.boundaryMesh()[patchIDs[patchI]].start();
-
-        label size  = mesh.boundaryMesh()[patchIDs[patchI]].size();
-
-        for (label i = 0; i < size; ++i)
-        {
-            faceLabels[++faceI] = start + i;
-        }
-    }
-
-    // Creating faMesh
-    Info << "Create faMesh ... ";
-
-    faMesh areaMesh
-    (
-        mesh,
-        faceLabels
-    );
-    Info << "Done" << endl;
-
-
-    // Determination of faPatch ID for each boundary edge.
-    // Result is in the bndEdgeFaPatchIDs list
-    const indirectPrimitivePatch& patch = areaMesh.patch();
-
-    labelList faceCells(faceLabels.size(), -1);
-
-    forAll(faceCells, faceI)
-    {
-        label faceID = faceLabels[faceI];
-
-        faceCells[faceI] = mesh.faceOwner()[faceID];
-    }
-
-    labelList meshEdges =
-        patch.meshEdges
-        (
-            mesh.edges(),
-            mesh.cellEdges(),
-            faceCells
-        );
-
-    const labelListList& edgeFaces = mesh.edgeFaces();
-
-    const label nTotalEdges = patch.nEdges();
-    const label nInternalEdges = patch.nInternalEdges();
-
-    labelList bndEdgeFaPatchIDs(nTotalEdges - nInternalEdges, -1);
-
-    for (label edgeI = nInternalEdges; edgeI < nTotalEdges; ++edgeI)
-    {
-        label curMeshEdge = meshEdges[edgeI];
-
-        labelList curEdgePatchIDs(2, label(-1));
-
-        label patchI = -1;
-
-        forAll(edgeFaces[curMeshEdge], faceI)
-        {
-            label curFace = edgeFaces[curMeshEdge][faceI];
-
-            label curPatchID = mesh.boundaryMesh().whichPatch(curFace);
-
-            if (curPatchID != -1)
-            {
-                curEdgePatchIDs[++patchI] = curPatchID;
-            }
-        }
-
-        for (label pI = 0; pI < faPatches.size() - 1; ++pI)
-        {
-            if
-            (
-                (
-                    curEdgePatchIDs[0] == faPatches[pI].ownPolyPatchID_
-                 && curEdgePatchIDs[1] == faPatches[pI].ngbPolyPatchID_
-                )
-                ||
-                (
-                    curEdgePatchIDs[1] == faPatches[pI].ownPolyPatchID_
-                 && curEdgePatchIDs[0] == faPatches[pI].ngbPolyPatchID_
-                )
-            )
-            {
-                bndEdgeFaPatchIDs[edgeI - nInternalEdges] = pI;
-                break;
-            }
-        }
-    }
-
-
-    // Set edgeLabels for each faPatch
-    for (label pI=0; pI<(faPatches.size()-1); ++pI)
-    {
-        SLList<label> tmpList;
-
-        forAll(bndEdgeFaPatchIDs, eI)
-        {
-            if (bndEdgeFaPatchIDs[eI] == pI)
-            {
-                tmpList.append(nInternalEdges + eI);
-            }
-        }
-
-        faPatches[pI].edgeLabels_ = tmpList;
-    }
-
-    // Check for undefined edges
-    SLList<label> tmpList;
-
-    forAll(bndEdgeFaPatchIDs, eI)
-    {
-        if (bndEdgeFaPatchIDs[eI] == -1)
-        {
-            tmpList.append(nInternalEdges + eI);
-        }
-    }
-
-    if (tmpList.size() > 0)
-    {
-        label pI = faPatches.size()-1;
-
-        faPatches[pI].name_ = "undefined";
-        faPatches[pI].type_ = "patch";
-        faPatches[pI].edgeLabels_ = tmpList;
-    }
-
-    // Add good patches to faMesh
-    SLList<faPatch*> faPatchLst;
-
-    for (label pI = 0; pI < faPatches.size(); ++pI)
-    {
-        faPatches[pI].dict_.add("type", faPatches[pI].type_);
-        faPatches[pI].dict_.add("edgeLabels", faPatches[pI].edgeLabels_);
-        faPatches[pI].dict_.add
-        (
-            "ngbPolyPatchIndex",
-            faPatches[pI].ngbPolyPatchID_
-        );
-
-        if(faPatches[pI].edgeLabels_.size() > 0)
-        {
-            faPatchLst.append
-            (
-                faPatch::New
-                (
-                    faPatches[pI].name_,
-                    faPatches[pI].dict_,
-                    pI,
-                    areaMesh.boundary()
-                ).ptr()
-            );
-        }
-    }
-
-    word emptyPatchName;
-    if (args.readIfPresent("addEmptyPatch", emptyPatchName))
-    {
-        dictionary emptyPatchDict;
-        emptyPatchDict.add("type", "empty");
-        emptyPatchDict.add("edgeLabels", labelList());
-        emptyPatchDict.add("ngbPolyPatchIndex", -1);
-
-        faPatchLst.append
-        (
-            faPatch::New
-            (
-                emptyPatchName,
-                emptyPatchDict,
-                faPatchLst.size(),
-                areaMesh.boundary()
-            ).ptr()
-        );
-    }
-
-    Info << "Add faPatches ... ";
-    areaMesh.addFaPatches(List<faPatch*>(faPatchLst));
-    Info << "Done" << endl;
+    // Creation
+    faMesh areaMesh(mesh, meshDefDict);
 
     // Writing faMesh
     Info << "Write finite area mesh ... ";
     areaMesh.write();
 
-    Info << "\nEnd" << endl;
+    #include "decomposeFaFields.H"
+
+    Info << "\nEnd\n" << endl;
 
     return 0;
 }
diff --git a/applications/utilities/parallelProcessing/decomposePar/faMeshDecomposition.C b/applications/utilities/parallelProcessing/decomposePar/faMeshDecomposition.C
index c08daf4b48..bdbb92dc79 100644
--- a/applications/utilities/parallelProcessing/decomposePar/faMeshDecomposition.C
+++ b/applications/utilities/parallelProcessing/decomposePar/faMeshDecomposition.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2018-2019 OpenCFD Ltd.
+    Copyright (C) 2018-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -41,6 +41,8 @@ License
 
 void Foam::faMeshDecomposition::distributeFaces()
 {
+    const word& polyMeshRegionName = mesh().name();
+
     Info<< "\nCalculating distribution of faces" << endl;
 
     cpuTime decompositionTime;
@@ -58,7 +60,7 @@ void Foam::faMeshDecomposition::distributeFaces()
         (
             IOobject
             (
-                GeoMesh<polyMesh>::mesh_.name(),
+                polyMeshRegionName,
                 processorDb.timeName(),
                 processorDb
             )
@@ -214,6 +216,8 @@ void Foam::faMeshDecomposition::decomposeMesh()
     // Decide which cell goes to which processor
     distributeFaces();
 
+    const word& polyMeshRegionName = mesh().name();
+
     Info<< "\nDistributing faces to processors" << endl;
 
     // Memory management
@@ -224,9 +228,9 @@ void Foam::faMeshDecomposition::decomposeMesh()
         {
             if (faceToProc_[faceI] >= nProcs())
             {
-                FatalErrorIn("Finite area mesh decomposition")
+                FatalErrorInFunction
                     << "Impossible processor label " << faceToProc_[faceI]
-                    << "for face " << faceI
+                    << "for face " << faceI << nl
                     << abort(FatalError);
             }
             else
@@ -258,7 +262,7 @@ void Foam::faMeshDecomposition::decomposeMesh()
         (
             IOobject
             (
-                GeoMesh<polyMesh>::mesh_.name(),
+                polyMeshRegionName,
                 processorDb.timeName(),
                 processorDb
             )
@@ -314,7 +318,7 @@ void Foam::faMeshDecomposition::decomposeMesh()
                 fvFaceProcAddressingHash.find
                 (
                     faceLabels()[curProcFaceAddressing[faceI]] + 1
-                )();
+                ).val();
         }
 
         // create processor finite area mesh
@@ -324,38 +328,35 @@ void Foam::faMeshDecomposition::decomposeMesh()
             procFaceLabels_[procI]
         );
 
-        const indirectPrimitivePatch& patch = this->patch();
+        const uindirectPrimitivePatch& patch = this->patch();
         const Map<label>& map = patch.meshPointMap();
 
         EdgeMap<label> edgesHash;
 
-        label edgeI = -1;
-
         const label nIntEdges = patch.nInternalEdges();
 
-        for (label curEdge = 0; curEdge < nIntEdges; curEdge++)
+        for (label edgei = 0; edgei < nIntEdges; ++edgei)
         {
-            edgesHash.insert(patch.edges()[curEdge], ++edgeI);
+            edgesHash.insert(patch.edges()[edgei], edgesHash.size());
         }
 
-        forAll(boundary(), patchI)
+        forAll(boundary(), patchi)
         {
             // Include emptyFaPatch
+            const label size = boundary()[patchi].labelList::size();
 
-            const label size = boundary()[patchI].labelList::size();
-
-            for(int eI=0; eI<size; eI++)
+            for (label edgei=0; edgei < size; ++edgei)
             {
                 edgesHash.insert
                 (
-                    patch.edges()[boundary()[patchI][eI]],
-                    ++edgeI
+                    patch.edges()[boundary()[patchi][edgei]],
+                    edgesHash.size()
                 );
             }
         }
 
 
-        const indirectPrimitivePatch& procPatch = procMesh.patch();
+        const uindirectPrimitivePatch& procPatch = procMesh.patch();
         const vectorField& procPoints = procPatch.localPoints();
         const labelList& procMeshPoints = procPatch.meshPoints();
         const edgeList& procEdges = procPatch.edges();
@@ -370,21 +371,18 @@ void Foam::faMeshDecomposition::decomposeMesh()
         }
 
         labelList& curPatchEdgeAddressing = procPatchEdgeAddressing_[procI];
-        curPatchEdgeAddressing.setSize(procEdges.size(), -1);
+        curPatchEdgeAddressing.resize(procEdges.size(), -1);
+
+        Map<label>& curMap = procMeshEdgesMap_[procI];
+        curMap.clear();
+        curMap.resize(2*procEdges.size());
 
         forAll(procEdges, edgeI)
         {
-            edge curGlobalEdge = procEdges[edgeI];
-            curGlobalEdge[0] = curPatchPointAddressing[curGlobalEdge[0]];
-            curGlobalEdge[1] = curPatchPointAddressing[curGlobalEdge[1]];
-
-            curPatchEdgeAddressing[edgeI] = edgesHash.find(curGlobalEdge)();
+            edge curGlobalEdge(curPatchPointAddressing, procEdges[edgeI]);
+            curPatchEdgeAddressing[edgeI] = edgesHash.find(curGlobalEdge).val();
         }
 
-        Map<label>& curMap = procMeshEdgesMap_[procI];
-
-        curMap = Map<label>(2*procEdges.size());
-
         forAll(curPatchEdgeAddressing, edgeI)
         {
             curMap.insert(curPatchEdgeAddressing[edgeI], edgeI);
@@ -1057,7 +1055,7 @@ void Foam::faMeshDecomposition::decomposeMesh()
         (
             IOobject
             (
-                GeoMesh<polyMesh>::mesh_.name(),
+                polyMeshRegionName,
                 processorDb.timeName(),
                 processorDb
             )
@@ -1137,8 +1135,9 @@ void Foam::faMeshDecomposition::decomposeMesh()
 
 bool Foam::faMeshDecomposition::writeDecomposition()
 {
-    Info<< "\nConstructing processor FA meshes" << endl;
+    const word& polyMeshRegionName = mesh().name();
 
+    Info<< "\nConstructing processor FA meshes" << endl;
 
     // Make a lookup map for globally shared points
     Map<label> sharedPointLookup(2*globallySharedPoints_.size());
@@ -1175,7 +1174,7 @@ bool Foam::faMeshDecomposition::writeDecomposition()
         (
             IOobject
             (
-                GeoMesh<polyMesh>::mesh_.name(),
+                polyMeshRegionName,
                 processorDb.timeName(),
                 processorDb
             )
@@ -1195,7 +1194,7 @@ bool Foam::faMeshDecomposition::writeDecomposition()
         );
 
 
-        // create finite area mesh
+        // Create finite area mesh
         faMesh procMesh
         (
             procFvMesh,
@@ -1219,11 +1218,9 @@ bool Foam::faMeshDecomposition::writeDecomposition()
 
         const faPatchList& meshPatches = boundary();
 
-        List<faPatch*> procPatches
+        PtrList<faPatch> procPatches
         (
-            curPatchSizes.size()
-          + curProcessorPatchSizes.size(),
-            reinterpret_cast<faPatch*>(NULL)
+            curPatchSizes.size() + curProcessorPatchSizes.size()
         );
 
         label nPatches = 0;
@@ -1232,44 +1229,51 @@ bool Foam::faMeshDecomposition::writeDecomposition()
         {
             const labelList& curEdgeLabels = curPatchEdgeLabels[nPatches];
 
-            const label ngbPolyPatchIndex =
+            const label neiPolyPatchId =
                 fvBoundaryProcAddressing.find
                 (
                     meshPatches[curBoundaryAddressing[patchi]]
                     .ngbPolyPatchIndex()
                 );
 
-            procPatches[nPatches] =
+            procPatches.set
+            (
+                nPatches,
                 meshPatches[curBoundaryAddressing[patchi]].clone
                 (
                     procMesh.boundary(),
                     curEdgeLabels,
                     nPatches,
-                    ngbPolyPatchIndex
-                ).ptr();
-
-            nPatches++;
+                    neiPolyPatchId
+                )
+            );
+            ++nPatches;
         }
 
         forAll(curProcessorPatchSizes, procPatchI)
         {
             const labelList& curEdgeLabels = curPatchEdgeLabels[nPatches];
 
-            procPatches[nPatches] =
+            procPatches.set
+            (
+                nPatches,
                 new processorFaPatch
                 (
-                    word("procBoundary") + Foam::name(procI)
-                  + word("to")
-                  + Foam::name(curNeighbourProcessors[procPatchI]),
+                    processorPolyPatch::newName
+                    (
+                        procI,
+                        curNeighbourProcessors[procPatchI]
+                    ),
                     curEdgeLabels,
                     nPatches,
                     procMesh.boundary(),
                     -1,
                     procI,
                     curNeighbourProcessors[procPatchI]
-                );
+                )
+            );
 
-            nPatches++;
+            ++nPatches;
         }
 
         // Add boundary patches
@@ -1291,23 +1295,19 @@ bool Foam::faMeshDecomposition::writeDecomposition()
 
         forAll(procMesh.boundary(), patchi)
         {
-            if
-            (
-                procMesh.boundary()[patchi].type()
-             == processorFaPatch::typeName
-            )
+            const auto* ppp =
+                isA<processorFaPatch>(procMesh.boundary()[patchi]);
+
+            if (ppp)
             {
-                const processorFaPatch& ppp =
-                refCast<const processorFaPatch>
-                (
-                    procMesh.boundary()[patchi]
-                );
+                const auto& procPatch = *ppp;
 
                 Info<< "    Number of edges shared with processor "
-                    << ppp.neighbProcNo() << " = " << ppp.size() << endl;
+                    << procPatch.neighbProcNo() << " = "
+                    << procPatch.size() << endl;
 
-                nProcPatches++;
-                nProcEdges += ppp.size();
+                nProcEdges += procPatch.size();
+                ++nProcPatches;
             }
             else
             {
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/writeAreaFields.H b/applications/utilities/postProcessing/dataConversion/foamToVTK/writeAreaFields.H
index 959881e65a..9bd2cb8c37 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/writeAreaFields.H
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/writeAreaFields.H
@@ -39,7 +39,7 @@ SourceFiles
 #define writeAreaFields_H
 
 #include "readFields.H"
-#include "foamVtkIndPatchGeoFieldsWriter.H"
+#include "foamVtkUIndPatchGeoFieldsWriter.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -47,7 +47,7 @@ namespace Foam
 {
 
 // Writer type for finite-area mesh + fields
-typedef vtk::indirectPatchGeoFieldsWriter vtkWriterType_areaMesh;
+typedef vtk::uindirectPatchGeoFieldsWriter vtkWriterType_areaMesh;
 
 template<class GeoField>
 bool writeAreaField
diff --git a/src/dynamicFaMesh/interfaceTrackingFvMesh/interfaceTrackingFvMesh.C b/src/dynamicFaMesh/interfaceTrackingFvMesh/interfaceTrackingFvMesh.C
index bf90e5220a..5e2cb8a3ca 100644
--- a/src/dynamicFaMesh/interfaceTrackingFvMesh/interfaceTrackingFvMesh.C
+++ b/src/dynamicFaMesh/interfaceTrackingFvMesh/interfaceTrackingFvMesh.C
@@ -48,7 +48,7 @@ License
 #include "turbulentTransportModel.H"
 #include "demandDrivenData.H"
 #include "unitConversion.H"
-#include "foamVtkIndPatchWriter.H"
+#include "foamVtkUIndPatchWriter.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -2270,7 +2270,7 @@ bool Foam::interfaceTrackingFvMesh::update()
 
 void Foam::interfaceTrackingFvMesh::writeVTK() const
 {
-    vtk::indirectPatchWriter writer
+    vtk::uindirectPatchWriter writer
     (
         aMesh().patch(),
         vtk::formatType::LEGACY_ASCII,
diff --git a/src/finiteArea/Make/files b/src/finiteArea/Make/files
index e218791232..7a96d8c4a1 100644
--- a/src/finiteArea/Make/files
+++ b/src/finiteArea/Make/files
@@ -1,11 +1,13 @@
 faMesh/faGlobalMeshData/faGlobalMeshData.C
 faMesh/faMesh.C
 faMesh/faMeshDemandDrivenData.C
+faMesh/faMeshPatches.C
 faMesh/faMeshUpdate.C
 faMesh/faBoundaryMesh/faBoundaryMesh.C
 
 faPatches = faMesh/faPatches
 $(faPatches)/faPatch/faPatch.C
+$(faPatches)/faPatch/faPatchData.C
 $(faPatches)/faPatch/faPatchNew.C
 $(faPatches)/basic/coupled/coupledFaPatch.C
 $(faPatches)/constraint/empty/emptyFaPatch.C
diff --git a/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.C b/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.C
index b436ecf996..da7962e51b 100644
--- a/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.C
+++ b/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.C
@@ -54,7 +54,7 @@ Foam::faBoundaryMesh::faBoundaryMesh
     {
         faPatchList& patches = *this;
 
-        // Read polyPatchList
+        // Read faPatch list
         Istream& is = readStream(typeName);
 
         PtrList<entry> patchEntries(is);
@@ -117,12 +117,6 @@ void Foam::faBoundaryMesh::calcGeometry()
 }
 
 
-const Foam::faMesh& Foam::faBoundaryMesh::mesh() const
-{
-    return mesh_;
-}
-
-
 Foam::lduInterfacePtrsList Foam::faBoundaryMesh::interfaces() const
 {
     lduInterfacePtrsList interfaces(size());
@@ -143,6 +137,26 @@ Foam::lduInterfacePtrsList Foam::faBoundaryMesh::interfaces() const
 }
 
 
+Foam::label Foam::faBoundaryMesh::nNonProcessor() const
+{
+    const faPatchList& patches = *this;
+
+    label nonProc = 0;
+
+    for (const faPatch& p : patches)
+    {
+        if (isA<processorFaPatch>(p))
+        {
+            break;
+        }
+
+        ++nonProc;
+    }
+
+    return nonProc;
+}
+
+
 Foam::wordList Foam::faBoundaryMesh::names() const
 {
     return PtrListOps::get<word>(*this, nameOp<faPatch>());
@@ -155,6 +169,71 @@ Foam::wordList Foam::faBoundaryMesh::types() const
 }
 
 
+Foam::labelList Foam::faBoundaryMesh::patchStarts() const
+{
+    // Manually: faPatch does not have independent start() information
+
+    const faPatchList& patches = *this;
+
+    labelList list(patches.size());
+
+    label beg = mesh_.nInternalEdges();
+    forAll(patches, patchi)
+    {
+        const label len = patches[patchi].nEdges();
+        list[patchi] = beg;
+        beg += len;
+    }
+    return list;
+}
+
+
+Foam::labelList Foam::faBoundaryMesh::patchSizes() const
+{
+    return
+        PtrListOps::get<label>
+        (
+            *this,
+            [](const faPatch& p) { return p.nEdges(); }  // avoid virtual
+        );
+}
+
+
+Foam::List<Foam::labelRange> Foam::faBoundaryMesh::patchRanges() const
+{
+    const faPatchList& patches = *this;
+
+    List<labelRange> list(patches.size());
+
+    label beg = mesh_.nInternalEdges();
+    forAll(patches, patchi)
+    {
+        const label len = patches[patchi].nEdges();
+        list[patchi].reset(beg, len);
+        beg += len;
+    }
+    return list;
+}
+
+
+Foam::label Foam::faBoundaryMesh::start() const
+{
+    return mesh_.nInternalEdges();
+}
+
+
+Foam::label Foam::faBoundaryMesh::nEdges() const
+{
+    return mesh_.nBoundaryEdges();
+}
+
+
+Foam::labelRange Foam::faBoundaryMesh::range() const
+{
+    return labelRange(mesh_.nInternalEdges(), mesh_.nBoundaryEdges());
+}
+
+
 Foam::labelList Foam::faBoundaryMesh::indices
 (
     const wordRe& matcher,
@@ -355,6 +434,22 @@ bool Foam::faBoundaryMesh::writeData(Ostream& os) const
 }
 
 
+bool Foam::faBoundaryMesh::writeObject
+(
+    IOstreamOption streamOpt,
+    const bool valid
+) const
+{
+    // Allow/disallow compression?
+    // 1. keep readable
+    // 2. save some space
+    // ??? streamOpt.compression(IOstreamOption::UNCOMPRESSED);
+    return regIOobject::writeObject(streamOpt, valid);
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
 Foam::Ostream& Foam::operator<<(Ostream& os, const faBoundaryMesh& bm)
 {
     bm.writeData(os);
diff --git a/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.H b/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.H
index 4c2c69ead5..397719aedf 100644
--- a/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.H
+++ b/src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.H
@@ -53,8 +53,6 @@ Author
 namespace Foam
 {
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // Forward Declarations
 class faMesh;
 class faBoundaryMesh;
@@ -89,14 +87,14 @@ public:
 
     // Constructors
 
-        //- Construct from dictionary
+        //- Construct from faMesh
         faBoundaryMesh
         (
             const IOobject& io,
             const faMesh& fam
         );
 
-        //- Construct given size
+        //- Construct from faMesh and given size
         faBoundaryMesh
         (
             const IOobject& io,
@@ -111,62 +109,95 @@ public:
 
     // Member Functions
 
-        // Access
+        //- Return the mesh reference
+        const faMesh& mesh() const noexcept
+        {
+            return mesh_;
+        }
 
-            //- Calculate the geometry for the patches
-            //  (transformation tensors etc.)
-            void calcGeometry();
+        //- Return a list of pointers for each patch
+        //- with only those pointing to interfaces being set
+        lduInterfacePtrsList interfaces() const;
 
-            //- Return the mesh reference
-            const faMesh& mesh() const;
+        //- The number of patches before the first processor patch.
+        label nNonProcessor() const;
 
-            //- Return a list of pointers for each patch
-            //- with only those pointing to interfaces being set
-            lduInterfacePtrsList interfaces() const;
+        //- Return a list of patch names
+        wordList names() const;
 
-            //- Return a list of patch names
-            wordList names() const;
+        //- Return a list of patch types
+        wordList types() const;
 
-            //- Return a list of patch types
-            wordList types() const;
+        //- Return a list of patch start indices
+        labelList patchStarts() const;
 
-            //- Return patch indices for all matches.
-            //  A no-op (returns -1) for an empty key
-            //  \note Matching patchGroups currently not supported
-            labelList indices
-            (
-                const wordRe& matcher,
-                const bool useGroups = false  /* ignored */
-            ) const;
+        //- Return a list of patch sizes (number of edges in each patch)
+        labelList patchSizes() const;
 
-            //- Return patch index for the first match, return -1 if not found
-            //  A no-op (returns -1) for an empty key
-            label findIndex(const wordRe& key) const;
+        //- Return a list of patch ranges
+        List<labelRange> patchRanges() const;
 
-            //- Find patch index given a name, return -1 if not found
-            //  A no-op (returns -1) for an empty name
-            label findPatchID(const word& patchName) const;
+        //- The start label of the edges in the faMesh edges list
+        //  Same as mesh.nInternalEdges()
+        label start() const;
 
-            //- Return patch index for a given edge label
-            label whichPatch(const label edgeIndex) const;
+        //- The number of boundary edges for the underlying mesh
+        //  Same as mesh.nBoundaryEdges()
+        label nEdges() const;
 
-            //- Check boundary definition
-            bool checkDefinition(const bool report = false) const;
+        //- The edge range for all boundary edges
+        //  Spans [nInternalEdges, nEdges) of the underlying mesh
+        labelRange range() const;
 
 
-        // Edit
+        //- Return patch indices for all matches.
+        //  A no-op (returns -1) for an empty key
+        //  \note Matching patchGroups currently not supported
+        labelList indices
+        (
+            const wordRe& matcher,
+            const bool useGroups = false  /* ignored */
+        ) const;
 
-            //- Correct faBoundaryMesh after moving points
-            void movePoints(const pointField&);
+        //- Return patch index for the first match, return -1 if not found
+        //  A no-op (returns -1) for an empty key
+        label findIndex(const wordRe& key) const;
 
-            //- Correct faBoundaryMesh after topology update
-            void updateMesh();
+        //- Find patch index given a name, return -1 if not found
+        //  A no-op (returns -1) for an empty name
+        label findPatchID(const word& patchName) const;
 
-            //- writeData member function required by regIOobject
-            bool writeData(Ostream&) const;
+        //- Return patch index for a given edge label
+        label whichPatch(const label edgeIndex) const;
+
+        //- Check boundary definition
+        bool checkDefinition(const bool report = false) const;
 
 
-    // Ostream operator
+    // Edit
+
+        //- Calculate the geometry for the patches
+        //  (transformation tensors etc.)
+        void calcGeometry();
+
+        //- Correct faBoundaryMesh after moving points
+        void movePoints(const pointField&);
+
+        //- Correct faBoundaryMesh after topology update
+        void updateMesh();
+
+        //- The writeData member function required by regIOobject
+        bool writeData(Ostream& os) const;
+
+        //- Write using stream options
+        virtual bool writeObject
+        (
+            IOstreamOption streamOpt,
+            const bool valid
+        ) const;
+
+
+    // Ostream Operator
 
         friend Ostream& operator<<(Ostream&, const faBoundaryMesh&);
 
diff --git a/src/finiteArea/faMesh/faMesh.C b/src/finiteArea/faMesh/faMesh.C
index 1e06c4cc03..d53be60604 100644
--- a/src/finiteArea/faMesh/faMesh.C
+++ b/src/finiteArea/faMesh/faMesh.C
@@ -47,75 +47,122 @@ namespace Foam
     defineTypeNameAndDebug(faMesh, 0);
 }
 
+
+const Foam::word Foam::faMesh::prefix("finite-area");
+
 Foam::word Foam::faMesh::meshSubDir = "faMesh";
 
-const int Foam::faMesh::quadricsFit_ = 0;
+const int Foam::faMesh::quadricsFit_ = 0;  // Tuning
+
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Convert patch names to face labels. Preserve patch order
+static labelList selectPatchFaces
+(
+    const polyBoundaryMesh& pbm,
+    const wordList& polyPatchNames
+)
+{
+    labelHashSet patchIDs;
+
+    label nFaceLabels = 0;
+    for (const word& patchName : polyPatchNames)
+    {
+        const label polyPatchi = pbm.findPatchID(patchName);
+
+        if (polyPatchi < 0)
+        {
+            FatalErrorInFunction
+                << "Patch " << patchName << " not found"
+                << exit(FatalError);
+        }
+
+        if (patchIDs.insert(polyPatchi))
+        {
+            nFaceLabels += pbm[polyPatchi].size();
+        }
+    }
+
+    labelList faceLabels(nFaceLabels);
+
+    nFaceLabels = 0;
+    for (const label polyPatchi : patchIDs.sortedToc())
+    {
+        for (const label facei : pbm[polyPatchi].range())
+        {
+            faceLabels[nFaceLabels] = facei;
+            ++nFaceLabels;
+        }
+    }
+
+    return faceLabels;
+}
+
+} // End namespace Foam
 
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+void Foam::faMesh::initPatch() const
+{
+    if (patchPtr_)
+    {
+        delete patchPtr_;
+    }
+    patchPtr_ = new uindirectPrimitivePatch
+    (
+        UIndirectList<face>(mesh().faces(), faceLabels_),
+        mesh().points()
+    );
+}
+
+
 void Foam::faMesh::setPrimitiveMeshData()
 {
     DebugInFunction << "Setting primitive data" << endl;
 
-    const indirectPrimitivePatch& bp = patch();
+    const uindirectPrimitivePatch& bp = patch();
+    const labelListList& edgeFaces = bp.edgeFaces();
 
-    // Set faMesh edges
-    edges_.setSize(bp.nEdges());
+    // Dimensions
 
-    label edgeI = -1;
-
-    label nIntEdges = bp.nInternalEdges();
-
-    for (label curEdge = 0; curEdge < nIntEdges; ++curEdge)
-    {
-        edges_[++edgeI] = bp.edges()[curEdge];
-    }
-
-    forAll(boundary(), patchI)
-    {
-        const labelList& curP = boundary()[patchI];
-
-        forAll(curP, eI)
-        {
-            edges_[++edgeI] = bp.edges()[curP[eI]];
-        }
-    }
-
-    nEdges_ = edges_.size();
-    nInternalEdges_ = nIntEdges;
-
-
-    // Set edge owner and neighbour
-    edgeOwner_.setSize(nEdges());
-    edgeNeighbour_.setSize(nInternalEdges());
-
-    edgeI = -1;
-
-    const labelListList& bpef = bp.edgeFaces();
-
-    for (label curEdge = 0; curEdge < nIntEdges; ++curEdge)
-    {
-        edgeOwner_[++edgeI] = bpef[curEdge][0];
-
-        edgeNeighbour_[edgeI] = bpef[curEdge][1];
-    }
-
-    forAll(boundary(), patchI)
-    {
-        const labelList& curP = boundary()[patchI];
-
-        forAll(curP, eI)
-        {
-            edgeOwner_[++edgeI] = bpef[curP[eI]][0];
-        }
-    }
-
-    // Set number of faces
+    nEdges_ = bp.nEdges();
+    nInternalEdges_ = bp.nInternalEdges();
     nFaces_ = bp.size();
-
-    // Set number of points
     nPoints_ = bp.nPoints();
+
+    edges_.resize(nEdges_);
+    edgeOwner_.resize(nEdges_);
+    edgeNeighbour_.resize(nInternalEdges_);
+
+    // Internal edges
+    for (label edgei = 0; edgei < nInternalEdges_; ++edgei)
+    {
+        edges_[edgei] = bp.edges()[edgei];
+
+        edgeOwner_[edgei] = edgeFaces[edgei][0];
+
+        edgeNeighbour_[edgei] = edgeFaces[edgei][1];
+    }
+
+    // Continue with boundary edges
+    label edgei = nInternalEdges_;
+
+    for (const faPatch& p : boundary())
+    {
+        for (const label patchEdgei : p.edgeLabels())
+        {
+            edges_[edgei] = bp.edges()[patchEdgei];
+
+            edgeOwner_[edgei] = edgeFaces[patchEdgei][0];
+
+            ++edgei;
+        }
+    }
 }
 
 
@@ -161,15 +208,20 @@ void Foam::faMesh::clearOut() const
 {
     clearGeom();
     clearAddressing();
-    deleteDemandDrivenData(globalMeshDataPtr_);
+    globalMeshDataPtr_.reset(nullptr);
 }
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+Foam::faMesh::faMesh(const polyMesh& pMesh, const zero)
+:
+    faMesh(pMesh, labelList())
+{}
+
+
 Foam::faMesh::faMesh(const polyMesh& pMesh)
 :
-    GeoMesh<polyMesh>(pMesh),
     MeshObject<polyMesh, Foam::UpdateableMeshObject, faMesh>(pMesh),
     edgeInterpolation(*this),
     faSchemes(mesh()),
@@ -181,7 +233,7 @@ Foam::faMesh::faMesh(const polyMesh& pMesh)
         (
             "faceLabels",
             time().findInstance(meshDir(), "faceLabels"),
-            meshSubDir,
+            faMesh::meshSubDir,
             mesh(),
             IOobject::MUST_READ,
             IOobject::NO_WRITE
@@ -193,7 +245,7 @@ Foam::faMesh::faMesh(const polyMesh& pMesh)
         (
             "faBoundary",
             time().findInstance(meshDir(), "faBoundary"),
-            meshSubDir,
+            faMesh::meshSubDir,
             mesh(),
             IOobject::MUST_READ,
             IOobject::NO_WRITE
@@ -220,7 +272,7 @@ Foam::faMesh::faMesh(const polyMesh& pMesh)
     correctPatchPointNormalsPtr_(nullptr),
     globalMeshDataPtr_(nullptr)
 {
-    DebugInFunction << "Creating faMesh from IOobject" << endl;
+    DebugInFunction << "Creating from IOobject" << endl;
 
     setPrimitiveMeshData();
 
@@ -244,7 +296,7 @@ Foam::faMesh::faMesh(const polyMesh& pMesh)
             (
                 "S0",
                 time().timeName(),
-                meshSubDir,
+                faMesh::meshSubDir,
                 mesh(),
                 IOobject::MUST_READ,
                 IOobject::AUTO_WRITE
@@ -258,10 +310,9 @@ Foam::faMesh::faMesh(const polyMesh& pMesh)
 Foam::faMesh::faMesh
 (
     const polyMesh& pMesh,
-    const labelList& faceLabels
+    const UList<label>& faceLabels
 )
 :
-    GeoMesh<polyMesh>(pMesh),
     MeshObject<polyMesh, Foam::UpdateableMeshObject, faMesh>(pMesh),
     edgeInterpolation(*this),
     faSchemes(mesh()),
@@ -273,7 +324,7 @@ Foam::faMesh::faMesh
         (
             "faceLabels",
             mesh().facesInstance(),
-            meshSubDir,
+            faMesh::meshSubDir,
             mesh(),
             IOobject::NO_READ,
             IOobject::NO_WRITE
@@ -286,13 +337,13 @@ Foam::faMesh::faMesh
         (
             "faBoundary",
             mesh().facesInstance(),
-            meshSubDir,
+            faMesh::meshSubDir,
             mesh(),
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
         *this,
-        0
+        label(0)
     ),
     comm_(Pstream::worldComm),
     patchPtr_(nullptr),
@@ -313,434 +364,77 @@ Foam::faMesh::faMesh
     edgeTransformTensorsPtr_(nullptr),
     correctPatchPointNormalsPtr_(nullptr),
     globalMeshDataPtr_(nullptr)
+{}
+
+
+Foam::faMesh::faMesh(const polyPatch& pp)
+:
+    faMesh
+    (
+        pp.boundaryMesh().mesh(),
+        identity(pp.range())
+    )
 {
-    DebugInFunction << "Creating faMesh from components" << endl;
+    DebugInFunction << "Creating from polyPatch:" << pp.name() << endl;
+
+    // Add single faPatch "default", but with processor connections
+    PtrList<faPatch> newPatches
+    (
+        createOnePatch("default")
+    );
+
+    addFaPatches(newPatches);
+
+    setPrimitiveMeshData();
+
+    // Create global mesh data
+    if (Pstream::parRun())
+    {
+        globalData();
+    }
+
+    // Calculate topology for the patches (processor-processor comms etc.)
+    boundary_.updateMesh();
+
+    // Calculate the geometry for the patches (transformation tensors etc.)
+    boundary_.calcGeometry();
 }
 
 
 Foam::faMesh::faMesh
 (
     const polyMesh& pMesh,
-    const fileName& defFile
+    const dictionary& faMeshDefinition
 )
 :
-    GeoMesh<polyMesh>(pMesh),
-    MeshObject<polyMesh, Foam::UpdateableMeshObject, faMesh>(pMesh),
-    edgeInterpolation(*this),
-    faSchemes(mesh()),
-    faSolution(mesh()),
-    data(mesh()),
-    faceLabels_
+    faMesh
     (
-        IOobject
+        pMesh,
+        selectPatchFaces
         (
-            "faceLabels",
-            mesh().facesInstance(),
-            meshSubDir,
-            mesh(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        List<label>(0)
-    ),
-    boundary_
-    (
-        IOobject
-        (
-            "faBoundary",
-            mesh().facesInstance(),
-            meshSubDir,
-            mesh(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        *this,
-        0
-    ),
-    comm_(Pstream::worldComm),
-    patchPtr_(nullptr),
-    lduPtr_(nullptr),
-    curTimeIndex_(time().timeIndex()),
-    SPtr_(nullptr),
-    S0Ptr_(nullptr),
-    S00Ptr_(nullptr),
-    patchStartsPtr_(nullptr),
-    LePtr_(nullptr),
-    magLePtr_(nullptr),
-    centresPtr_(nullptr),
-    edgeCentresPtr_(nullptr),
-    faceAreaNormalsPtr_(nullptr),
-    edgeAreaNormalsPtr_(nullptr),
-    pointAreaNormalsPtr_(nullptr),
-    faceCurvaturesPtr_(nullptr),
-    edgeTransformTensorsPtr_(nullptr),
-    correctPatchPointNormalsPtr_(nullptr),
-    globalMeshDataPtr_(nullptr)
+            pMesh.boundaryMesh(),
+            faMeshDefinition.get<wordList>("polyMeshPatches")
+        )
+    )
 {
-    DebugInFunction << "Creating faMesh from definition file" << endl;
+    DebugInFunction << "Creating from definition (dictionary)" << endl;
 
-    // Reading faMeshDefinition dictionary
-    IOdictionary faMeshDefinition
+    PtrList<faPatch> newPatches
     (
-        IOobject
+        createPatchList
         (
-            defFile,
-            mesh().time().constant(),
-            meshSubDir,
-            mesh(),
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
+            faMeshDefinition.subDict("boundary"),
+
+            // Optional 'empty' patch
+            faMeshDefinition.getOrDefault<word>("emptyPatch", word::null),
+
+            // Optional specification for default patch
+            faMeshDefinition.findDict("defaultPatch")
         )
     );
 
-    const wordList polyMeshPatches
-    (
-        faMeshDefinition.get<wordList>("polyMeshPatches")
-    );
 
-    const dictionary& bndDict = faMeshDefinition.subDict("boundary");
-
-    const wordList faPatchNames(bndDict.toc());
-
-    List<faPatchData> faPatches(faPatchNames.size() + 1);
-
-    const polyBoundaryMesh& pbm = pMesh.boundaryMesh();
-
-    forAll(faPatchNames, patchI)
-    {
-        dictionary curPatchDict = bndDict.subDict(faPatchNames[patchI]);
-
-        faPatches[patchI].name_ = faPatchNames[patchI];
-
-        faPatches[patchI].type_ = curPatchDict.get<word>("type");
-
-        faPatches[patchI].ownPolyPatchID_ =
-            pbm.findPatchID(curPatchDict.get<word>("ownerPolyPatch"));
-
-        faPatches[patchI].ngbPolyPatchID_ =
-            pbm.findPatchID(curPatchDict.get<word>("neighbourPolyPatch"));
-    }
-
-
-    // Setting faceLabels list size
-    label size = 0;
-
-    labelList patchIDs(polyMeshPatches.size(), -1);
-
-    forAll(polyMeshPatches, patchI)
-    {
-        patchIDs[patchI] = pbm.findPatchID(polyMeshPatches[patchI]);
-
-        size += pbm[patchIDs[patchI]].size();
-    }
-
-    faceLabels_ = labelList(size, -1);
-
-
-    // Filling of faceLabels list
-    label faceI = -1;
-
-    sort(patchIDs);
-
-    forAll(polyMeshPatches, patchI)
-    {
-        label start = pbm[patchIDs[patchI]].start();
-        label size  = pbm[patchIDs[patchI]].size();
-
-        for (label i = 0; i < size; ++i)
-        {
-            faceLabels_[++faceI] = start + i;
-        }
-    }
-
-
-    // Determination of faPatch ID for each boundary edge.
-    // Result is in the bndEdgeFaPatchIDs list
-    labelList faceCells(faceLabels_.size(), -1);
-
-    forAll(faceCells, faceI)
-    {
-        label faceID = faceLabels_[faceI];
-
-        faceCells[faceI] = mesh().faceOwner()[faceID];
-    }
-
-    labelList meshEdges =
-        patch().meshEdges
-        (
-            mesh().edges(),
-            mesh().cellEdges(),
-            faceCells
-        );
-
-    const labelListList& edgeFaces = mesh().edgeFaces();
-
-    const label nTotalEdges = patch().nEdges();
-    const label nInternalEdges = patch().nInternalEdges();
-
-    labelList bndEdgeFaPatchIDs(nTotalEdges - nInternalEdges, -1);
-
-    for (label edgeI = nInternalEdges; edgeI < nTotalEdges; ++edgeI)
-    {
-        label curMeshEdge = meshEdges[edgeI];
-
-        labelList curEdgePatchIDs(2, label(-1));
-
-        label patchI = -1;
-
-        forAll(edgeFaces[curMeshEdge], faceI)
-        {
-            label curFace = edgeFaces[curMeshEdge][faceI];
-
-            label curPatchID = pbm.whichPatch(curFace);
-
-            if (curPatchID != -1)
-            {
-                curEdgePatchIDs[++patchI] = curPatchID;
-            }
-        }
-
-        for (label pI = 0; pI < faPatches.size() - 1; ++pI)
-        {
-            if
-            (
-                (
-                    curEdgePatchIDs[0] == faPatches[pI].ownPolyPatchID_
-                 && curEdgePatchIDs[1] == faPatches[pI].ngbPolyPatchID_
-                )
-             ||
-                (
-                    curEdgePatchIDs[1] == faPatches[pI].ownPolyPatchID_
-                 && curEdgePatchIDs[0] == faPatches[pI].ngbPolyPatchID_
-                )
-            )
-            {
-                bndEdgeFaPatchIDs[edgeI - nInternalEdges] = pI;
-                break;
-            }
-        }
-    }
-
-
-    // Set edgeLabels for each faPatch
-    for (label pI = 0; pI < (faPatches.size() - 1); ++pI)
-    {
-        SLList<label> tmpList;
-
-        forAll(bndEdgeFaPatchIDs, eI)
-        {
-            if (bndEdgeFaPatchIDs[eI] == pI)
-            {
-                tmpList.append(nInternalEdges + eI);
-            }
-        }
-
-        faPatches[pI].edgeLabels_ = tmpList;
-    }
-
-    // Check for undefined edges
-    SLList<label> tmpList;
-
-    forAll(bndEdgeFaPatchIDs, eI)
-    {
-        if (bndEdgeFaPatchIDs[eI] == -1)
-        {
-            tmpList.append(nInternalEdges + eI);
-        }
-    }
-
-    if (tmpList.size() > 0)
-    {
-        // Check for processor edges
-        labelList allUndefEdges(tmpList);
-        labelList ngbPolyPatch(allUndefEdges.size(), -1);
-        forAll(ngbPolyPatch, edgeI)
-        {
-            label curEdge = allUndefEdges[edgeI];
-
-            label curPMeshEdge = meshEdges[curEdge];
-
-            forAll(edgeFaces[curPMeshEdge], faceI)
-            {
-                label curFace = edgeFaces[curPMeshEdge][faceI];
-
-                if (faceLabels_.found(curFace))
-                {
-                    label polyPatchID = pbm.whichPatch(curFace);
-
-                    if (polyPatchID != -1)
-                    {
-                        ngbPolyPatch[edgeI] = polyPatchID;
-                    }
-                }
-            }
-        }
-
-        // Count ngb processorPolyPatch-es
-        labelHashSet processorPatchSet;
-        forAll(ngbPolyPatch, edgeI)
-        {
-            if (ngbPolyPatch[edgeI] != -1)
-            {
-                if (isA<processorPolyPatch>(pbm[ngbPolyPatch[edgeI]]))
-                {
-                    processorPatchSet.insert(ngbPolyPatch[edgeI]);
-                }
-            }
-        }
-        labelList processorPatches(processorPatchSet.toc());
-        faPatches.setSize(faPatches.size() + processorPatches.size());
-
-        for (label i = 0; i < processorPatches.size(); ++i)
-        {
-            SLList<label> tmpLst;
-
-            forAll(ngbPolyPatch, eI)
-            {
-                if (ngbPolyPatch[eI] == processorPatches[i])
-                {
-                    tmpLst.append(allUndefEdges[eI]);
-                }
-            }
-
-            faPatches[faPatchNames.size() + i].edgeLabels_ = tmpLst;
-
-            faPatches[faPatchNames.size() + i].name_ =
-                pbm[processorPatches[i]].name();
-
-            faPatches[faPatchNames.size() + i].type_ =
-                processorFaPatch::typeName;
-
-            faPatches[faPatchNames.size() + i].ngbPolyPatchID_ =
-                processorPatches[i];
-        }
-
-        // Remaining undefined edges
-        SLList<label> undefEdges;
-        forAll(ngbPolyPatch, eI)
-        {
-            if (ngbPolyPatch[eI] == -1)
-            {
-                undefEdges.append(allUndefEdges[eI]);
-            }
-            else if (!isA<processorPolyPatch>(pbm[ngbPolyPatch[eI]]))
-            {
-                undefEdges.append(allUndefEdges[eI]);
-            }
-        }
-
-        if (undefEdges.size() > 0)
-        {
-            label pI = faPatches.size()-1;
-
-            faPatches[pI].name_ = "undefined";
-            faPatches[pI].type_ = "patch";
-            faPatches[pI].edgeLabels_ = undefEdges;
-        }
-        else
-        {
-            faPatches.setSize(faPatches.size() - 1);
-        }
-    }
-    else
-    {
-        faPatches.setSize(faPatches.size() - 1);
-    }
-
-
-    // Reorder processorFaPatch using
-    // ordering of ngb processorPolyPatch
-    forAll(faPatches, patchI)
-    {
-        if (faPatches[patchI].type_ == processorFaPatch::typeName)
-        {
-            labelList ngbFaces(faPatches[patchI].edgeLabels_.size(), -1);
-
-            forAll(ngbFaces, edgeI)
-            {
-                label curEdge = faPatches[patchI].edgeLabels_[edgeI];
-
-                label curPMeshEdge = meshEdges[curEdge];
-
-                forAll(edgeFaces[curPMeshEdge], faceI)
-                {
-                    label curFace = edgeFaces[curPMeshEdge][faceI];
-
-                    label curPatchID = pbm.whichPatch(curFace);
-
-                    if (curPatchID == faPatches[patchI].ngbPolyPatchID_)
-                    {
-                        ngbFaces[edgeI] = curFace;
-                    }
-                }
-            }
-
-            SortableList<label> sortedNgbFaces(ngbFaces);
-            labelList reorderedEdgeLabels(ngbFaces.size(), -1);
-            for (label i = 0; i < reorderedEdgeLabels.size(); ++i)
-            {
-                reorderedEdgeLabels[i] =
-                    faPatches[patchI].edgeLabels_
-                    [
-                        sortedNgbFaces.indices()[i]
-                    ];
-            }
-
-            faPatches[patchI].edgeLabels_ = reorderedEdgeLabels;
-        }
-    }
-
-
-    // Add good patches to faMesh
-    SLList<faPatch*> faPatchLst;
-
-    for (label pI = 0; pI < faPatches.size(); ++pI)
-    {
-        faPatches[pI].dict_.add("type", faPatches[pI].type_);
-        faPatches[pI].dict_.add("edgeLabels", faPatches[pI].edgeLabels_);
-        faPatches[pI].dict_.add
-        (
-            "ngbPolyPatchIndex",
-            faPatches[pI].ngbPolyPatchID_
-        );
-
-        if (faPatches[pI].type_ == processorFaPatch::typeName)
-        {
-            if (faPatches[pI].ngbPolyPatchID_ == -1)
-            {
-                FatalErrorInFunction
-                    << "ngbPolyPatch is not defined for processorFaPatch: "
-                    << faPatches[pI].name_
-                    << abort(FatalError);
-            }
-
-            const processorPolyPatch& procPolyPatch =
-                refCast<const processorPolyPatch>
-                (
-                    pbm[faPatches[pI].ngbPolyPatchID_]
-                );
-
-            faPatches[pI].dict_.add("myProcNo", procPolyPatch.myProcNo());
-            faPatches[pI].dict_.add
-            (
-                "neighbProcNo",
-                procPolyPatch.neighbProcNo()
-            );
-        }
-
-        faPatchLst.append
-        (
-            faPatch::New
-            (
-                faPatches[pI].name_,
-                faPatches[pI].dict_,
-                pI,
-                boundary()
-            ).ptr()
-        );
-    }
-
-    addFaPatches(List<faPatch*>(faPatchLst));
+    addFaPatches(newPatches);
 
     // Create global mesh data
     if (Pstream::parRun())
@@ -762,7 +456,7 @@ Foam::faMesh::faMesh
             (
                 "S0",
                 time().timeName(),
-                meshSubDir,
+                faMesh::meshSubDir,
                 mesh(),
                 IOobject::MUST_READ,
                 IOobject::AUTO_WRITE
@@ -773,111 +467,6 @@ Foam::faMesh::faMesh
 }
 
 
-Foam::faMesh::faMesh
-(
-    const polyMesh& pMesh,
-    const label polyPatchID
-)
-:
-    GeoMesh<polyMesh>(pMesh),
-    MeshObject<polyMesh, Foam::UpdateableMeshObject, faMesh>(pMesh),
-    edgeInterpolation(*this),
-    faSchemes(mesh()),
-    faSolution(mesh()),
-    data(mesh()),
-    faceLabels_
-    (
-        IOobject
-        (
-            "faceLabels",
-            mesh().facesInstance(),
-            meshSubDir,
-            mesh(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        labelList(pMesh.boundaryMesh()[polyPatchID].size(), -1)
-    ),
-    boundary_
-    (
-        IOobject
-        (
-            "faBoundary",
-            mesh().facesInstance(),
-            meshSubDir,
-            mesh(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        *this,
-        0
-    ),
-    comm_(Pstream::worldComm),
-    patchPtr_(nullptr),
-    lduPtr_(nullptr),
-    curTimeIndex_(time().timeIndex()),
-    SPtr_(nullptr),
-    S0Ptr_(nullptr),
-    S00Ptr_(nullptr),
-    patchStartsPtr_(nullptr),
-    LePtr_(nullptr),
-    magLePtr_(nullptr),
-    centresPtr_(nullptr),
-    edgeCentresPtr_(nullptr),
-    faceAreaNormalsPtr_(nullptr),
-    edgeAreaNormalsPtr_(nullptr),
-    pointAreaNormalsPtr_(nullptr),
-    faceCurvaturesPtr_(nullptr),
-    edgeTransformTensorsPtr_(nullptr),
-    correctPatchPointNormalsPtr_(nullptr),
-    globalMeshDataPtr_(nullptr)
-{
-    DebugInFunction << "Creating faMesh from polyPatch" << endl;
-
-    const polyBoundaryMesh& pbm = pMesh.boundaryMesh();
-
-    // Set faceLabels
-    forAll(faceLabels_, faceI)
-    {
-        faceLabels_[faceI] = pbm[polyPatchID].start() + faceI;
-    }
-
-    // Add one faPatch
-    const indirectPrimitivePatch& bp = patch();
-
-    const label nTotalEdges = bp.nEdges();
-
-    const label nInternalEdges = bp.nInternalEdges();
-
-    labelList edgeLabels(nTotalEdges - nInternalEdges, -1);
-
-    forAll(edgeLabels, edgeI)
-    {
-        edgeLabels[edgeI] = nInternalEdges + edgeI;
-    }
-
-    dictionary patchDict;
-
-    patchDict.add("type", "patch");
-    patchDict.add("edgeLabels", edgeLabels);
-    patchDict.add("ngbPolyPatchIndex", -1);
-
-    List<faPatch*> faPatchLst(1);
-
-    faPatchLst[0] = faPatch::New("default", patchDict, 0, boundary()).ptr();
-
-    addFaPatches(faPatchLst);
-
-    setPrimitiveMeshData();
-
-    // Calculate topology for the patches (processor-processor comms etc.)
-    boundary_.updateMesh();
-
-    // Calculate the geometry for the patches (transformation tensors etc.)
-    boundary_.calcGeometry();
-}
-
-
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::faMesh::~faMesh()
@@ -890,7 +479,7 @@ Foam::faMesh::~faMesh()
 
 Foam::fileName Foam::faMesh::meshDir() const
 {
-    return mesh().dbDir()/meshSubDir;
+    return mesh().dbDir()/faMesh::meshSubDir;
 }
 
 
@@ -912,98 +501,6 @@ const Foam::fileName& Foam::faMesh::facesInstance() const
 }
 
 
-const Foam::indirectPrimitivePatch& Foam::faMesh::patch() const
-{
-    if (!patchPtr_)
-    {
-        patchPtr_ = new indirectPrimitivePatch
-        (
-            IndirectList<face>
-            (
-                mesh().faces(),
-                faceLabels_
-            ),
-            mesh().points()
-        );
-    }
-
-    return *patchPtr_;
-}
-
-
-Foam::indirectPrimitivePatch& Foam::faMesh::patch()
-{
-    if (!patchPtr_)
-    {
-        patchPtr_ = new indirectPrimitivePatch
-        (
-            IndirectList<face>
-            (
-                mesh().faces(),
-                faceLabels_
-            ),
-            mesh().points()
-        );
-    }
-
-    return *patchPtr_;
-}
-
-
-const Foam::pointField& Foam::faMesh::points() const
-{
-    return patch().localPoints();
-}
-
-
-const Foam::edgeList& Foam::faMesh::edges() const
-{
-    return edges_;
-}
-
-
-const Foam::faceList& Foam::faMesh::faces() const
-{
-    return patch().localFaces();
-}
-
-
-void Foam::faMesh::addFaPatches(const List<faPatch*>& p)
-{
-    DebugInFunction << "Adding patches to faMesh" << endl;
-
-    if (boundary().size() > 0)
-    {
-        FatalErrorInFunction
-            << "boundary already exists"
-            << abort(FatalError);
-    }
-
-    boundary_.setSize(p.size());
-
-    forAll(p, patchI)
-    {
-        boundary_.set(patchI, p[patchI]);
-    }
-
-    setPrimitiveMeshData();
-
-    boundary_.checkDefinition();
-}
-
-
-Foam::label Foam::faMesh::comm() const
-{
-    return comm_;
-}
-
-
-Foam::label& Foam::faMesh::comm()
-{
-    return comm_;
-}
-
-
 bool Foam::faMesh::hasDb() const
 {
     return true;
@@ -1016,12 +513,6 @@ const Foam::objectRegistry& Foam::faMesh::thisDb() const
 }
 
 
-const Foam::faBoundaryMesh& Foam::faMesh::boundary() const
-{
-    return boundary_;
-}
-
-
 const Foam::labelList& Foam::faMesh::patchStarts() const
 {
     if (!patchStartsPtr_)
@@ -1193,7 +684,7 @@ const Foam::faGlobalMeshData& Foam::faMesh::globalData() const
 {
     if (!globalMeshDataPtr_)
     {
-        globalMeshDataPtr_ = new faGlobalMeshData(*this);
+        globalMeshDataPtr_.reset(new faGlobalMeshData(*this));
     }
 
     return *globalMeshDataPtr_;
diff --git a/src/finiteArea/faMesh/faMesh.H b/src/finiteArea/faMesh/faMesh.H
index 8715a22f01..b6515b510f 100644
--- a/src/finiteArea/faMesh/faMesh.H
+++ b/src/finiteArea/faMesh/faMesh.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -42,7 +43,6 @@ Author
 #ifndef faMesh_H
 #define faMesh_H
 
-#include "GeoMesh.H"
 #include "MeshObject.H"
 #include "polyMesh.H"
 #include "lduMesh.H"
@@ -53,7 +53,7 @@ Author
 #include "DimensionedField.H"
 #include "areaFieldsFwd.H"
 #include "edgeFieldsFwd.H"
-#include "indirectPrimitivePatch.H"
+#include "uindirectPrimitivePatch.H"
 #include "edgeInterpolation.H"
 #include "labelIOList.H"
 #include "FieldFields.H"
@@ -67,9 +67,10 @@ Author
 namespace Foam
 {
 
-// Class forward declarations
+// Forward Declarations
 class faMeshLduAddressing;
 class faMeshMapper;
+class faPatchData;
 
 /*---------------------------------------------------------------------------*\
                            Class faMesh Declaration
@@ -77,7 +78,6 @@ class faMeshMapper;
 
 class faMesh
 :
-    public GeoMesh<polyMesh>,
     public MeshObject<polyMesh, Foam::UpdateableMeshObject, faMesh>,
     public lduMesh,
     public edgeInterpolation,
@@ -85,7 +85,7 @@ class faMesh
     public faSolution,
     public data
 {
-    // Private data
+    // Private Data
 
         //- Face labels
         labelIOList faceLabels_;
@@ -130,7 +130,7 @@ class faMesh
     // Demand-driven data
 
         //- Primitive patch
-        mutable indirectPrimitivePatch* patchPtr_;
+        mutable uindirectPrimitivePatch* patchPtr_;
 
         //- Ldu addressing data
         mutable faMeshLduAddressing* lduPtr_;
@@ -186,8 +186,8 @@ class faMesh
 
         // Other mesh-related data
 
-            //- Parallel info
-            mutable faGlobalMeshData* globalMeshDataPtr_;
+        //- Parallel info
+        mutable autoPtr<faGlobalMeshData> globalMeshDataPtr_;
 
 
     // Static Private Data
@@ -204,6 +204,8 @@ class faMesh
         //- No copy assignment
         void operator=(const faMesh&) = delete;
 
+        //- Set indirect patch, removing any old one
+        void initPatch() const;
 
         //- Set primitive mesh data
         void setPrimitiveMeshData();
@@ -262,9 +264,36 @@ class faMesh
             //- Clear demand-driven data
             void clearOut() const;
 
+
+    // Helpers
+
+        //- Create a single patch
+        PtrList<faPatch> createOnePatch
+        (
+            const word& patchName,
+            const word& patchType = ""
+        ) const;
+
+        //- Create list of patches from boundary definition
+        PtrList<faPatch> createPatchList
+        (
+            const dictionary& bndDict,
+            const word& emptyPatchName = "",
+            const dictionary* defaultPatchDefinition = nullptr
+        ) const;
+
+        //- Reorder processor edges using order of the
+        //- neighbour processorPolyPatch
+        void reorderProcEdges
+        (
+            faPatchData& patchDef,
+            const labelUList& meshEdges
+        ) const;
+
+
 public:
 
-    // Public typedefs
+    // Public Typedefs
 
         typedef faMesh Mesh;
         typedef faBoundaryMesh BoundaryMesh;
@@ -273,36 +302,38 @@ public:
     //- Runtime type information
     TypeName("faMesh");
 
+        //- The prefix to local: %finite-area
+        static const word prefix;
 
-    //- Return the mesh sub-directory name (usually "faMesh")
-    static word meshSubDir;
+        //- The mesh sub-directory name (usually "faMesh")
+        static word meshSubDir;
 
 
     // Constructors
 
-        //- Construct from polyMesh
-        explicit faMesh(const polyMesh& m);
+        //- Construct zero-sized from polyMesh
+        //  Boundary is added using addFaPatches() member function
+        faMesh(const polyMesh& pMesh, const Foam::zero);
 
-        //- Construct from components without boundary.
+        //- Construct from polyMesh
+        explicit faMesh(const polyMesh& pMesh);
+
+        //- Construct for specified face labels without boundary.
         //  Boundary is added using addFaPatches() member function
         faMesh
         (
-            const polyMesh& m,
-            const labelList& faceLabels
+            const polyMesh& pMesh,
+            const UList<label>& faceLabels
         );
 
-        //- Construct from finite area mesh definition file
-        faMesh
-        (
-            const polyMesh& m,
-            const fileName& defFile
-        );
+        //- Construct from single polyPatch
+        explicit faMesh(const polyPatch& pp);
 
-        //- Construct from polyPatch
+        //- Construct from definition
         faMesh
         (
-            const polyMesh& m,
-            const label polyPatchID
+            const polyMesh& pMesh,
+            const dictionary& faMeshDefinition
         );
 
 
@@ -312,25 +343,30 @@ public:
 
     // Member Functions
 
-        // Helpers
+    // Helpers
 
-            //- Add boundary patches. Constructor helper
-            void addFaPatches(const List<faPatch*> &);
+        //- Add boundary patches. Constructor helper
+        void addFaPatches
+        (
+            PtrList<faPatch>& plist,
+            const bool validBoundary = true
+        );
+
+        //- Add boundary patches. Constructor helper
+        void addFaPatches
+        (
+            const List<faPatch*>& p,
+            const bool validBoundary = true
+        );
 
 
         // Database
 
             //- Return access to polyMesh
-            const polyMesh& mesh() const
-            {
-                return
-                    MeshObject
-                    <
-                        polyMesh,
-                        Foam::UpdateableMeshObject,
-                        faMesh
-                    >::mesh();
-            }
+            inline const polyMesh& mesh() const;
+
+            //- Interface to referenced polyMesh (similar to GeoMesh)
+            const polyMesh& operator()() const { return mesh(); }
 
             //- Return the local mesh directory (dbDir()/meshSubDir)
             fileName meshDir() const;
@@ -347,61 +383,51 @@ public:
             const fileName& facesInstance() const;
 
 
+        // Communication support
+
+            //- Return communicator used for parallel communication
+            inline label comm() const noexcept;
+
+            //- Return communicator used for parallel communication
+            inline label& comm() noexcept;
+
+
             // Mesh size parameters
 
-                inline label nPoints() const
-                {
-                    return nPoints_;
-                }
+                //- Number of local mesh points
+                inline label nPoints() const noexcept;
 
-                inline label nEdges() const
-                {
-                    return nEdges_;
-                }
+                //- Number of local mesh edges
+                inline label nEdges() const noexcept;
 
-                inline label nInternalEdges() const
-                {
-                    return nInternalEdges_;
-                }
+                //- Number of internal faces
+                inline label nInternalEdges() const noexcept;
 
-                inline label nFaces() const
-                {
-                    return nFaces_;
-                }
+                //- Number of boundary edges (== nEdges - nInternalEdges)
+                inline label nBoundaryEdges() const noexcept;
+
+                //- Number of patch faces
+                inline label nFaces() const noexcept;
 
 
             // Primitive mesh data
 
-                //- Return mesh points
-                const pointField& points() const;
+                //- Return local patch points
+                inline const pointField& points() const;
 
-                //- Return edges
-                const edgeList& edges() const;
+                //- Return local patch edges with reordered boundary
+                inline const edgeList& edges() const noexcept;
 
-                //- Return faces
-                const faceList& faces() const;
+                //- Return local patch faces
+                inline const faceList& faces() const;
 
                 //- Edge owner addressing
-                inline const labelList& edgeOwner() const
-                {
-                    return edgeOwner_;
-                }
+                inline const labelList& edgeOwner() const noexcept;
 
                 //- Edge neighbour addressing
-                inline const labelList& edgeNeighbour() const
-                {
-                    return edgeNeighbour_;
-                }
+                inline const labelList& edgeNeighbour() const noexcept;
 
 
-            // Communication support
-
-                //- Return communicator used for parallel communication
-                label comm() const;
-
-                //- Return communicator used for parallel communication
-                label& comm();
-
 
         // Access
 
@@ -419,13 +445,10 @@ public:
             }
 
             //- Return constant reference to boundary mesh
-            const faBoundaryMesh& boundary() const;
+            inline const faBoundaryMesh& boundary() const noexcept;
 
             //- Return faMesh face labels
-            const labelList& faceLabels() const
-            {
-                return faceLabels_;
-            }
+            inline const labelList& faceLabels() const noexcept;
 
 
             //- Return parallel info
@@ -453,10 +476,10 @@ public:
                 return lduAddr().upperAddr();
             }
 
-            //- Return true if given edge label is internal to the mesh
-            inline bool isInternalEdge(const label edgeIndex) const
+            //- True if given edge label is internal to the mesh
+            bool isInternalEdge(const label edgeIndex) const
             {
-                return edgeIndex < nInternalEdges();
+                return (edgeIndex < nInternalEdges_);
             }
 
 
@@ -487,10 +510,10 @@ public:
         // Demand-driven data
 
             //- Return constant reference to primitive patch
-            const indirectPrimitivePatch& patch() const;
+            inline const uindirectPrimitivePatch& patch() const;
 
             //- Return reference to primitive patch
-            indirectPrimitivePatch& patch();
+            inline uindirectPrimitivePatch& patch();
 
             //- Return patch starts
             const labelList& patchStarts() const;
@@ -546,6 +569,7 @@ public:
             //- Set whether point normals should be corrected for a patch
             boolList& correctPatchPointNormals() const;
 
+
         //- Write mesh
         virtual bool write(const bool valid = true) const;
 
@@ -564,6 +588,8 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#include "faMeshI.H"
+
 #ifdef NoRepository
     #include "faPatchFaMeshTemplates.C"
 #endif
diff --git a/src/finiteArea/faMesh/faMeshDemandDrivenData.C b/src/finiteArea/faMesh/faMeshDemandDrivenData.C
index 7a363fd2ab..63ec9815bd 100644
--- a/src/finiteArea/faMesh/faMeshDemandDrivenData.C
+++ b/src/finiteArea/faMesh/faMeshDemandDrivenData.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -40,6 +40,30 @@ License
 #include "processorFaPatchFields.H"
 #include "emptyFaPatchFields.H"
 
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// A bitSet (size patch nPoints()) with boundary points marked
+static Foam::bitSet markupBoundaryPoints(const uindirectPrimitivePatch& p)
+{
+    // Initially all unmarked
+    bitSet markPoints(p.nPoints());
+    for (const edge& e : p.boundaryEdges())
+    {
+        // Mark boundary points
+        markPoints.set(e.first());
+        markPoints.set(e.second());
+    }
+
+    return markPoints;
+}
+
+} // End namespace Foam
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void Foam::faMesh::calcLduAddressing() const
@@ -70,16 +94,7 @@ void Foam::faMesh::calcPatchStarts() const
             << abort(FatalError);
     }
 
-    patchStartsPtr_ = new labelList(boundary().size(), -1);
-    labelList& patchStarts = *patchStartsPtr_;
-
-    patchStarts[0] = nInternalEdges();
-
-    for (label i = 1; i < boundary().size(); ++i)
-    {
-        patchStarts[i] =
-            patchStarts[i - 1] + boundary()[i - 1].faPatch::size();
-    }
+    patchStartsPtr_ = new labelList(boundary().patchStarts());
 }
 
 
@@ -211,13 +226,11 @@ void Foam::faMesh::calcMagLe() const
 
     const pointField& localPoints = points();
 
-    const edgeList::subList internalEdges =
-        edgeList::subList(edges(), nInternalEdges());
-
-
-    forAll(internalEdges, edgeI)
+    label edgei = 0;
+    for (const edge& e : patch().internalEdges())
     {
-        magLe.ref()[edgeI] = internalEdges[edgeI].mag(localPoints);
+        magLe.ref()[edgei] = e.mag(localPoints);
+        ++edgei;
     }
 
 
@@ -331,13 +344,12 @@ void Foam::faMesh::calcEdgeCentres() const
 
     const pointField& localPoints = points();
 
-    const edgeList::subList internalEdges =
-        edgeList::subList(edges(), nInternalEdges());
 
-
-    forAll(internalEdges, edgeI)
+    label edgei = 0;
+    for (const edge& e : patch().internalEdges())
     {
-        edgeCentres.ref()[edgeI] = internalEdges[edgeI].centre(localPoints);
+        edgeCentres.ref()[edgei] = e.centre(localPoints);
+        ++edgei;
     }
 
 
@@ -850,31 +862,10 @@ Foam::labelList Foam::faMesh::internalPoints() const
     DebugInFunction
         << "Calculating internal points" << endl;
 
-    const edgeList& edges = patch().edges();
-    label nIntEdges = patch().nInternalEdges();
+    bitSet markPoints(markupBoundaryPoints(this->patch()));
+    markPoints.flip();
 
-    List<bool> internal(nPoints(), true);
-
-    for (label curEdge = nIntEdges; curEdge < edges.size(); ++curEdge)
-    {
-        internal[edges[curEdge].start()] = false;
-
-        internal[edges[curEdge].end()] = false;
-    }
-
-    SLList<label> internalPoints;
-
-    forAll(internal, pointI)
-    {
-        if (internal[pointI])
-        {
-            internalPoints.append(pointI);
-        }
-    }
-
-    labelList result(internalPoints);
-
-    return result;
+    return markPoints.sortedToc();
 }
 
 
@@ -883,31 +874,9 @@ Foam::labelList Foam::faMesh::boundaryPoints() const
     DebugInFunction
         << "Calculating boundary points" << endl;
 
-    const edgeList& edges = patch().edges();
-    label nIntEdges = patch().nInternalEdges();
+    bitSet markPoints(markupBoundaryPoints(this->patch()));
 
-    List<bool> internal(nPoints(), true);
-
-    for (label curEdge = nIntEdges; curEdge < edges.size(); ++curEdge)
-    {
-        internal[edges[curEdge].start()] = false;
-
-        internal[edges[curEdge].end()] = false;
-    }
-
-    SLList<label> boundaryPoints;
-
-    forAll(internal, pointI)
-    {
-        if (!internal[pointI])
-        {
-            boundaryPoints.append(pointI);
-        }
-    }
-
-    labelList result(boundaryPoints);
-
-    return result;
+    return markPoints.sortedToc();
 }
 
 
@@ -1168,7 +1137,7 @@ void Foam::faMesh::calcPointAreaNormals() const
 
         if (correctPatchPointNormals(patchI) && !fap.coupled())
         {
-            if (fap.ngbPolyPatchIndex() == -1)
+            if (fap.ngbPolyPatchIndex() < 0)
             {
                 FatalErrorInFunction
                     << "Neighbour polyPatch index is not defined "
diff --git a/src/finiteArea/faMesh/faMeshI.H b/src/finiteArea/faMesh/faMeshI.H
new file mode 100644
index 0000000000..a3d7932860
--- /dev/null
+++ b/src/finiteArea/faMesh/faMeshI.H
@@ -0,0 +1,142 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2016-2017 Wikki Ltd
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const Foam::polyMesh& Foam::faMesh::mesh() const
+{
+    return
+        MeshObject<polyMesh, Foam::UpdateableMeshObject, faMesh>::mesh();
+}
+
+
+inline const Foam::faBoundaryMesh& Foam::faMesh::boundary() const noexcept
+{
+    return boundary_;
+}
+
+
+inline Foam::label Foam::faMesh::comm() const noexcept
+{
+    return comm_;
+}
+
+
+inline Foam::label& Foam::faMesh::comm() noexcept
+{
+    return comm_;
+}
+
+
+inline Foam::label Foam::faMesh::nPoints() const noexcept
+{
+    return nPoints_;
+}
+
+
+inline Foam::label Foam::faMesh::nEdges() const noexcept
+{
+    return nEdges_;
+}
+
+
+inline Foam::label Foam::faMesh::nInternalEdges() const noexcept
+{
+    return nInternalEdges_;
+}
+
+
+inline Foam::label Foam::faMesh::nBoundaryEdges() const noexcept
+{
+    return nEdges_ - nInternalEdges_;
+}
+
+
+inline Foam::label Foam::faMesh::nFaces() const noexcept
+{
+    return nFaces_;
+}
+
+
+inline const Foam::pointField& Foam::faMesh::points() const
+{
+    return patch().localPoints();
+}
+
+
+inline const Foam::edgeList& Foam::faMesh::edges() const noexcept
+{
+    return edges_;
+}
+
+
+inline const Foam::faceList& Foam::faMesh::faces() const
+{
+    return patch().localFaces();
+}
+
+
+inline const Foam::labelList& Foam::faMesh::edgeOwner() const noexcept
+{
+    return edgeOwner_;
+}
+
+
+inline const Foam::labelList& Foam::faMesh::edgeNeighbour() const noexcept
+{
+    return edgeNeighbour_;
+}
+
+
+inline const Foam::labelList& Foam::faMesh::faceLabels() const noexcept
+{
+    return faceLabels_;
+}
+
+
+inline const Foam::uindirectPrimitivePatch& Foam::faMesh::patch() const
+{
+    if (!patchPtr_)
+    {
+        initPatch();
+    }
+    return *patchPtr_;
+}
+
+
+inline Foam::uindirectPrimitivePatch& Foam::faMesh::patch()
+{
+    if (!patchPtr_)
+    {
+        initPatch();
+    }
+    return *patchPtr_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteArea/faMesh/faMeshPatches.C b/src/finiteArea/faMesh/faMeshPatches.C
new file mode 100644
index 0000000000..765252c670
--- /dev/null
+++ b/src/finiteArea/faMesh/faMeshPatches.C
@@ -0,0 +1,572 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2017 Wikki Ltd
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "faMesh.H"
+#include "IndirectList.H"
+#include "faPatchData.H"
+#include "processorPolyPatch.H"
+#include "processorFaPatch.H"
+#include "edgeHashes.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::faMesh::reorderProcEdges
+(
+    faPatchData& patchDef,
+    const labelUList& meshEdges
+) const
+{
+    if (!patchDef.coupled() || patchDef.edgeLabels_.empty())
+    {
+        return;
+    }
+
+    const polyBoundaryMesh& pbm = mesh().boundaryMesh();
+    const labelListList& edgeFaces = mesh().edgeFaces();
+
+    // Reorder processor edges using order of neighbour processorPolyPatch
+
+    const label procPatchID = patchDef.neighPolyPatchId_;
+    const label nProcEdges = patchDef.edgeLabels_.size();
+
+    labelList procFaces(nProcEdges, -1);
+
+    forAll(procFaces, edgei)
+    {
+        const label localEdgei = patchDef.edgeLabels_[edgei];
+        const label meshEdgei = meshEdges[localEdgei];
+
+        for (const label meshFacei : edgeFaces[meshEdgei])
+        {
+            if
+            (
+                !faceLabels_.found(meshFacei)
+             && (procPatchID == pbm.whichPatch(meshFacei))
+            )
+            {
+                // The edge's proc-face
+                procFaces[edgei] = meshFacei;
+                break;
+            }
+        }
+    }
+
+    // Ascending proc-face numbering
+    const labelList sortIndices(Foam::sortedOrder(procFaces));
+
+    const labelList& oldEdgeLabels = patchDef.edgeLabels_;
+    labelList newEdgeLabels(oldEdgeLabels.size());
+
+    // Most of the time, an individual proc-face will only be singly
+    // attached to the finite-area patch. In rarer case, there could
+    // multiple connections. For these cases, need to walk the face
+    // edges - the direction depends on owner vs neighbour side.
+
+    EdgeMap<label> multihit;
+
+    for (label edgei = 0; edgei < nProcEdges; /*nil*/)
+    {
+        const label procFacei = procFaces[sortIndices[edgei]];
+
+        // Find all identical faces
+        label endEdgei = edgei + 1;  // one beyond
+        while
+        (
+            (endEdgei < nProcEdges)
+         && (procFacei == procFaces[sortIndices[endEdgei]])
+        )
+        {
+            ++endEdgei;
+        }
+
+        if (edgei + 1 == endEdgei)
+        {
+            // Simplest case - a single connection
+
+            newEdgeLabels[edgei] = oldEdgeLabels[sortIndices[edgei]];
+        }
+        else
+        {
+            multihit.clear();
+
+            // Map from global edge to local edgeId
+            for (label i = edgei; i < endEdgei; ++i)
+            {
+                label localEdgei = oldEdgeLabels[sortIndices[i]];
+                label meshEdgei = meshEdges[localEdgei];
+
+                multihit.insert
+                (
+                    mesh().edges()[meshEdgei],
+                    localEdgei
+                );
+            }
+
+            if (multihit.size() != (endEdgei - edgei))
+            {
+                FatalErrorInFunction
+                    << "Could only hash " << multihit.size()
+                    << " edges from " << (endEdgei - edgei)
+                    << " ... indicates a non-manifold connection" << nl
+                    << multihit << nl
+                    << abort(FatalError);
+            }
+
+            const face& f = mesh().faces()[procFacei];
+
+            forAll(f, fedgei)  // Note size() == nEdges()
+            {
+                edge e =
+                (
+                    patchDef.owner()
+                  ? f.edge(fedgei)      // Forward walk
+                  : f.rcEdge(fedgei)    // Reverse walk
+                );
+
+                auto iter = multihit.find(e);
+                if (iter.found())
+                {
+                    newEdgeLabels[edgei++] = iter.val();
+                    multihit.erase(iter);
+                    if (multihit.empty())
+                    {
+                        break;
+                    }
+                }
+            }
+
+            if (edgei != endEdgei)
+            {
+                FatalErrorInFunction
+                    << "Missed " << (edgei < endEdgei)
+                    << " edges for face: " << procFacei
+                    << " ... indicates serious geometry issue" << nl
+                    << multihit << nl
+                    << abort(FatalError);
+            }
+            if (!multihit.empty())
+            {
+                FatalErrorInFunction
+                    << "Missed edges for face: " << procFacei
+                    << " ... indicates serious geometry issue" << nl
+                    << multihit << nl
+                    << abort(FatalError);
+            }
+        }
+        edgei = endEdgei;
+    }
+
+    patchDef.edgeLabels_.transfer(newEdgeLabels);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::faMesh::addFaPatches
+(
+    PtrList<faPatch>& plist,
+    const bool validBoundary
+)
+{
+    if (!boundary().empty())
+    {
+        FatalErrorInFunction
+            << "boundary already exists"
+            << abort(FatalError);
+    }
+
+    globalMeshDataPtr_.reset(nullptr);
+
+    boundary_.transfer(plist);
+
+    setPrimitiveMeshData();
+
+    if (validBoundary)
+    {
+        boundary_.checkDefinition();
+    }
+}
+
+
+void Foam::faMesh::addFaPatches
+(
+    const List<faPatch*>& p,
+    const bool validBoundary
+)
+{
+    // Acquire ownership of the pointers
+    PtrList<faPatch> plist(const_cast<List<faPatch*>&>(p));
+
+    addFaPatches(plist, validBoundary);
+}
+
+
+Foam::PtrList<Foam::faPatch> Foam::faMesh::createOnePatch
+(
+    const word& patchName,
+    const word& patchType
+) const
+{
+    dictionary onePatchDict;
+    if (!patchName.empty())
+    {
+        onePatchDict.add("name", patchName);
+    }
+    if (!patchType.empty())
+    {
+        onePatchDict.add("type", patchType);
+    }
+
+    return createPatchList
+    (
+        dictionary::null,
+        "",             // Name for empty patch placeholder
+        &onePatchDict   // Definitions for defaultPatch
+    );
+}
+
+
+Foam::PtrList<Foam::faPatch> Foam::faMesh::createPatchList
+(
+    const dictionary& bndDict,
+    const word& emptyPatchName,
+    const dictionary* defaultPatchDefinition
+) const
+{
+    const polyBoundaryMesh& pbm = mesh().boundaryMesh();
+
+    // Transcribe into patch definitions
+    DynamicList<faPatchData> faPatchDefs(bndDict.size() + 4);
+    for (const entry& dEntry : bndDict)
+    {
+        if (!dEntry.isDict())
+        {
+            WarningInFunction
+                << "Not a dictionary entry: " << dEntry.name() << nl;
+            continue;
+        }
+        const dictionary& patchDict = dEntry.dict();
+
+        // Add entry
+        faPatchDefs.append(faPatchData());
+
+        auto& patchDef = faPatchDefs.last();
+        patchDef.name_ = dEntry.keyword();
+        patchDef.type_ = patchDict.get<word>("type");
+
+        const word ownName(patchDict.get<word>("ownerPolyPatch"));
+        const word neiName(patchDict.get<word>("neighbourPolyPatch"));
+
+        patchDef.ownerPolyPatchId_ = pbm.findPatchID(ownName);
+        patchDef.neighPolyPatchId_ = pbm.findPatchID(neiName);
+
+        if (patchDef.ownerPolyPatchId_ < 0)
+        {
+            FatalErrorInFunction
+                << "ownerPolyPatch " << ownName << " not found"
+                << exit(FatalError);
+        }
+        if (patchDef.neighPolyPatchId_ < 0)
+        {
+            FatalErrorInFunction
+                << "neighbourPolyPatch " << neiName << " not found"
+                << exit(FatalError);
+        }
+    }
+
+    // Additional empty placeholder patch?
+    if (!emptyPatchName.empty())
+    {
+        faPatchDefs.append(faPatchData());
+
+        auto& patchDef = faPatchDefs.last();
+        patchDef.name_ = emptyPatchName;
+        patchDef.type_ = "empty";
+    }
+
+    // Placeholder for any undefined edges
+    const label undefPatchId = faPatchDefs.size();
+    {
+        faPatchDefs.append(faPatchData());
+
+        auto& patchDef = faPatchDefs.last();
+        patchDef.name_ = "undefined";
+        patchDef.type_ = "patch";
+
+        if (defaultPatchDefinition)
+        {
+            (*defaultPatchDefinition).readIfPresent("name", patchDef.name_);
+            (*defaultPatchDefinition).readIfPresent("type", patchDef.type_);
+        }
+    }
+
+    // ----------------------------------------------------------------------
+
+    // Determine faPatch ID for each boundary edge.
+    // Result is in the bndEdgeFaPatchIDs list
+
+    const labelList meshEdges
+    (
+        patch().meshEdges(mesh().edges(), mesh().pointEdges())
+    );
+
+    const labelListList& edgeFaces = mesh().edgeFaces();
+
+    const label nInternalEdges = patch().nInternalEdges();
+    const label nBoundaryEdges = patch().nBoundaryEdges();
+
+    labelList bndEdgeFaPatchIDs(nBoundaryEdges, -1);
+
+    for (label bndEdgei = 0; bndEdgei < nBoundaryEdges; ++bndEdgei)
+    {
+        const label patchEdgei = meshEdges[bndEdgei + nInternalEdges];
+
+        // Use 'edge' for accounting
+        edge curEdgePatchPair;
+
+        for (const label meshFacei : edgeFaces[patchEdgei])
+        {
+            const label polyPatchID = pbm.whichPatch(meshFacei);
+
+            if (polyPatchID != -1)
+            {
+                curEdgePatchPair.insert(polyPatchID);
+            }
+        }
+
+
+        if (curEdgePatchPair.valid())
+        {
+            // Non-negative, unique pairing
+            // - find corresponding definition
+
+            for (label patchi = 0; patchi < faPatchDefs.size(); ++patchi)
+            {
+                if (faPatchDefs[patchi].foundPatchPair(curEdgePatchPair))
+                {
+                    bndEdgeFaPatchIDs[bndEdgei] = patchi;
+                    break;
+                }
+            }
+        }
+    }
+
+
+    // Extract which edges map to which patch
+    // and set edgeLabels for each faPatch
+
+    DynamicList<label> selectEdges(bndEdgeFaPatchIDs.size());
+
+    for (label patchi = 0; patchi < faPatchDefs.size(); ++patchi)
+    {
+        auto& patchDef = faPatchDefs[patchi];
+
+        selectEdges.clear();
+
+        forAll(bndEdgeFaPatchIDs, bndEdgei)
+        {
+            if (bndEdgeFaPatchIDs[bndEdgei] == patchi)
+            {
+                selectEdges.append(bndEdgei + nInternalEdges);
+            }
+        }
+
+        patchDef.edgeLabels_ = selectEdges;
+    }
+
+    // Check for undefined edges
+    selectEdges.clear();
+
+    forAll(bndEdgeFaPatchIDs, bndEdgei)
+    {
+        if (bndEdgeFaPatchIDs[bndEdgei] == -1)
+        {
+            selectEdges.append(bndEdgei + nInternalEdges);
+        }
+    }
+
+    // Save the information
+    faPatchDefs[undefPatchId].edgeLabels_ = selectEdges;
+
+    bool hasUndefined = returnReduce(!selectEdges.empty(), orOp<bool>());
+
+    if (hasUndefined)
+    {
+        // The initial edges to consider
+        const labelList& undefinedEdges =
+            faPatchDefs[undefPatchId].edgeLabels_;
+
+        // Check for edges that butt against a processor (or other) patch
+
+        labelList edgeNbrPolyPatch(undefinedEdges.size(), -1);
+        forAll(edgeNbrPolyPatch, edgei)
+        {
+            const label localEdgei = undefinedEdges[edgei];
+            const label meshEdgei = meshEdges[localEdgei];
+            label polyPatchID;
+
+            for (const label meshFacei : edgeFaces[meshEdgei])
+            {
+                if
+                (
+                    !faceLabels_.found(meshFacei)
+                 && (polyPatchID = pbm.whichPatch(meshFacei)) != -1
+                )
+                {
+                    // Found the edge's off-patch neighbour face
+                    edgeNbrPolyPatch[edgei] = polyPatchID;
+                    break;
+                }
+            }
+        }
+
+        // Categorize as processor/non-processor associations
+        labelHashSet procPatchIDs;
+        labelHashSet nonProcPatchIDs;
+
+        for (const label polyPatchID : edgeNbrPolyPatch)
+        {
+            if (polyPatchID == -1)
+            {
+                nonProcPatchIDs.insert(polyPatchID);
+            }
+            else if
+            (
+                !nonProcPatchIDs.found(polyPatchID)
+             && !procPatchIDs.found(polyPatchID)
+            )
+            {
+                if (isA<processorPolyPatch>(pbm[polyPatchID]))
+                {
+                    procPatchIDs.insert(polyPatchID);
+                }
+                else
+                {
+                    nonProcPatchIDs.insert(polyPatchID);
+                }
+            }
+        }
+
+        // Select by processor association
+        for (const label polyPatchID : procPatchIDs.sortedToc())
+        {
+            selectEdges.clear();
+
+            forAll(edgeNbrPolyPatch, edgei)
+            {
+                if (edgeNbrPolyPatch[edgei] == polyPatchID)
+                {
+                    selectEdges.append(undefinedEdges[edgei]);
+                }
+            }
+
+            faPatchDefs.append(faPatchData());
+
+            auto& patchDef = faPatchDefs.last();
+            patchDef.name_ = pbm[polyPatchID].name();
+            patchDef.type_ = processorFaPatch::typeName;
+            patchDef.neighPolyPatchId_ = polyPatchID;  // Needed for reorder
+
+            const auto* ppp = isA<processorPolyPatch>(pbm[polyPatchID]);
+            if (ppp)
+            {
+                patchDef.ownerProcId_ = ppp->myProcNo();
+                patchDef.neighProcId_ = ppp->neighbProcNo();
+            }
+
+            patchDef.edgeLabels_ = selectEdges;
+        }
+
+
+        // Check for any remaining undefined edges
+        selectEdges.clear();
+
+        // Simply grab any/all (don't worry about which patch)
+        if (!nonProcPatchIDs.empty())
+        {
+            forAll(edgeNbrPolyPatch, edgei)
+            {
+                const label polyPatchID = edgeNbrPolyPatch[edgei];
+                if (nonProcPatchIDs.found(polyPatchID))
+                {
+                    selectEdges.append(undefinedEdges[edgei]);
+                }
+            }
+        }
+
+        // Complete the information
+        faPatchDefs[undefPatchId].edgeLabels_ = selectEdges;
+
+        hasUndefined = returnReduce(!selectEdges.empty(), orOp<bool>());
+    }
+
+    // Remove unnecessary entry
+    if (!hasUndefined)
+    {
+        faPatchDefs.remove(undefPatchId);
+    }
+
+    for (auto& patchDef : faPatchDefs)
+    {
+        if (patchDef.coupled())
+        {
+            reorderProcEdges(patchDef, meshEdges);
+            patchDef.neighPolyPatchId_ = -1; // No longer required + confusing
+        }
+    }
+
+
+    // Now convert list of definitions to list of patches
+
+    label nPatches = 0;
+    PtrList<faPatch> newPatches(faPatchDefs.size());
+
+    for (faPatchData& patchDef : faPatchDefs)
+    {
+        newPatches.set
+        (
+            nPatches,
+            faPatch::New
+            (
+                patchDef.name(),        // name
+                patchDef.dict(false),   // withEdgeLabels == false
+                nPatches,               // index
+                boundary()
+            )
+        );
+
+        // Transfer edge labels
+        newPatches[nPatches].resetEdges(std::move(patchDef.edgeLabels_));
+        ++nPatches;
+    }
+
+    return newPatches;
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteArea/faMesh/faMeshUpdate.C b/src/finiteArea/faMesh/faMeshUpdate.C
index cb518d77ef..1e8593b261 100644
--- a/src/finiteArea/faMesh/faMeshUpdate.C
+++ b/src/finiteArea/faMesh/faMeshUpdate.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -61,11 +61,12 @@ void Foam::faMesh::updateMesh(const mapPolyMesh& mpm)
     // Set new labels
     m.faceLabels_ = mapper.areaMap().newFaceLabels();
 
-    const indirectPrimitivePatch& bp = patch();
+    const uindirectPrimitivePatch& bp = patch();
 
     // Collect patch data
     const label nTotalEdges = bp.nEdges();
     const label nInternalEdges = bp.nInternalEdges();
+    const label nBoundaryEdges = bp.nBoundaryEdges();
     const labelListList& edgeFaces = bp.edgeFaces();
 
     labelListList patchEdges(boundary_.size());
@@ -73,7 +74,7 @@ void Foam::faMesh::updateMesh(const mapPolyMesh& mpm)
     // Special handling required for faces that have more than one edge
     // Each patch will be visited separately
 
-    labelList edgeToPatch(nTotalEdges - nInternalEdges, -1);
+    labelList edgeToPatch(nBoundaryEdges, -1);
     const labelList& newFaceLabelsMap = mapper.areaMap().newFaceLabelsMap();
 
     const labelListList& oldPatchEdgeFaces = mapper.oldPatchEdgeFaces();
@@ -81,7 +82,7 @@ void Foam::faMesh::updateMesh(const mapPolyMesh& mpm)
     forAll(oldPatchEdgeFaces, patchI)
     {
         labelList& curPatchEdges = patchEdges[patchI];
-        curPatchEdges.setSize(nTotalEdges - nInternalEdges);
+        curPatchEdges.resize(nBoundaryEdges);
         label nCurPatchEdges = 0;
 
         // Note: it is possible to pick up the old-to-new boundary patch
diff --git a/src/finiteArea/faMesh/faPatches/basic/coupled/coupledFaPatch.C b/src/finiteArea/faMesh/faPatches/basic/coupled/coupledFaPatch.C
index bbc073bfc2..8ca96bbdaa 100644
--- a/src/finiteArea/faMesh/faPatches/basic/coupled/coupledFaPatch.C
+++ b/src/finiteArea/faMesh/faPatches/basic/coupled/coupledFaPatch.C
@@ -117,10 +117,4 @@ void Foam::coupledFaPatch::calcTransformTensors
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::coupledFaPatch::~coupledFaPatch()
-{}
-
-
 // ************************************************************************* //
diff --git a/src/finiteArea/faMesh/faPatches/basic/coupled/coupledFaPatch.H b/src/finiteArea/faMesh/faPatches/basic/coupled/coupledFaPatch.H
index d5424a42a4..7f5ddc489b 100644
--- a/src/finiteArea/faMesh/faPatches/basic/coupled/coupledFaPatch.H
+++ b/src/finiteArea/faMesh/faPatches/basic/coupled/coupledFaPatch.H
@@ -135,7 +135,7 @@ public:
 
 
     //- Destructor
-    virtual ~coupledFaPatch();
+    virtual ~coupledFaPatch() = default;
 
 
     // Member Functions
diff --git a/src/finiteArea/faMesh/faPatches/constraint/empty/emptyFaPatch.C b/src/finiteArea/faMesh/faPatches/constraint/empty/emptyFaPatch.C
index 5b4c4bb971..93c69622e8 100644
--- a/src/finiteArea/faMesh/faPatches/constraint/empty/emptyFaPatch.C
+++ b/src/finiteArea/faMesh/faPatches/constraint/empty/emptyFaPatch.C
@@ -28,19 +28,50 @@ License
 #include "emptyFaPatch.H"
 #include "addToRunTimeSelectionTable.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // Patch name
 defineTypeNameAndDebug(emptyFaPatch, 0);
 
 // Add the patch constructor functions to the hash tables
 addToRunTimeSelectionTable(faPatch, emptyFaPatch, dictionary);
 
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::emptyFaPatch::emptyFaPatch
+(
+    const word& name,
+    const label index,
+    const faBoundaryMesh& bm,
+    const label ngbPolyPatchIndex
+)
+:
+    emptyFaPatch(name, labelList(), index, bm, ngbPolyPatchIndex)
+{}
+
+
+Foam::emptyFaPatch::emptyFaPatch
+(
+    const word& name,
+    const labelList& edgeLabels,
+    const label index,
+    const faBoundaryMesh& bm,
+    const label ngbPolyPatchIndex
+)
+:
+    faPatch(name, edgeLabels, index, bm, ngbPolyPatchIndex)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
 // Over-riding the face normals return from the underlying patch
 // This is the only piece of info used out of the underlying primitivePatch
 // I choose to store it there because it is used in primitive patch operations
@@ -54,8 +85,4 @@ addToRunTimeSelectionTable(faPatch, emptyFaPatch, dictionary);
 // }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //
diff --git a/src/finiteArea/faMesh/faPatches/constraint/empty/emptyFaPatch.H b/src/finiteArea/faMesh/faPatches/constraint/empty/emptyFaPatch.H
index 0ff674f566..e9ecd37363 100644
--- a/src/finiteArea/faMesh/faPatches/constraint/empty/emptyFaPatch.H
+++ b/src/finiteArea/faMesh/faPatches/constraint/empty/emptyFaPatch.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -57,7 +58,6 @@ class emptyFaPatch
 :
     public faPatch
 {
-
 public:
 
     //- Runtime type information
@@ -66,7 +66,16 @@ public:
 
     // Constructors
 
-         //- Construct from components
+        //- Minimal construct from components
+        emptyFaPatch
+        (
+            const word& name,
+            const label index,
+            const faBoundaryMesh& bm,
+            const label ngbPolyPatchIndex = -1
+        );
+
+        //- Construct from components
         emptyFaPatch
         (
             const word& name,
@@ -74,10 +83,7 @@ public:
             const label index,
             const faBoundaryMesh& bm,
             const label ngbPolyPatchIndex
-        )
-        :
-            faPatch(name, edgeLabels, index, bm, ngbPolyPatchIndex)
-        {}
+        );
 
         //- Construct from dictionary
         emptyFaPatch
@@ -114,6 +120,7 @@ public:
             );
         }
 
+
     // Member Functions
 
         virtual label size() const
@@ -124,7 +131,6 @@ public:
         //- Return face normals. Over-riding base class return to get zero size
         //
 //         virtual const vectorField& edgeNormals() const;
-
 };
 
 
diff --git a/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.H b/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.H
index 4bf3f9d683..c48962843c 100644
--- a/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.H
+++ b/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.H
@@ -122,7 +122,6 @@ protected:
 public:
 
     //- Runtime type information
-//     TypeName(processorPolyPatch::typeName_());
     TypeName("processor");
 
 
@@ -174,7 +173,7 @@ public:
     virtual ~processorFaPatch();
 
 
-    // Member functions
+    // Member Functions
 
         //- Return interface size
         virtual label interfaceSize() const
diff --git a/src/finiteArea/faMesh/faPatches/constraint/symmetry/symmetryFaPatch.C b/src/finiteArea/faMesh/faPatches/constraint/symmetry/symmetryFaPatch.C
index 437b03474b..bdf45fac31 100644
--- a/src/finiteArea/faMesh/faPatches/constraint/symmetry/symmetryFaPatch.C
+++ b/src/finiteArea/faMesh/faPatches/constraint/symmetry/symmetryFaPatch.C
@@ -70,7 +70,7 @@ Foam::symmetryFaPatch::symmetryFaPatch
 :
     faPatch(name, dict, index, bm)
 {
-    if (ngbPolyPatchIndex() == -1)
+    if (ngbPolyPatchIndex() < 0)
     {
         FatalErrorInFunction
             << "Neighbour polyPatch index is not specified for faPatch "
diff --git a/src/finiteArea/faMesh/faPatches/constraint/wedge/wedgeFaPatch.C b/src/finiteArea/faMesh/faPatches/constraint/wedge/wedgeFaPatch.C
index 8cca67583a..c841637f45 100644
--- a/src/finiteArea/faMesh/faPatches/constraint/wedge/wedgeFaPatch.C
+++ b/src/finiteArea/faMesh/faPatches/constraint/wedge/wedgeFaPatch.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -87,22 +88,21 @@ Foam::wedgeFaPatch::wedgeFaPatch
     axisPoint_(-1),
     axisPointChecked_(false)
 {
-    if (ngbPolyPatchIndex() == -1)
+    if (ngbPolyPatchIndex() < 0)
     {
         FatalErrorInFunction
             << "Neighbour polyPatch index is not specified for faPatch "
             << this->name() << exit(FatalError);
     }
 
-    if (isA<wedgePolyPatch>(bm.mesh()().boundaryMesh()[ngbPolyPatchIndex()]))
-    {
-        const wedgePolyPatch& wedge =
-            refCast<const wedgePolyPatch>
-            (
-                bm.mesh()().boundaryMesh()[ngbPolyPatchIndex()]
-            );
+    const auto* wedgePtr = isA<wedgePolyPatch>
+    (
+        bm.mesh()().boundaryMesh()[ngbPolyPatchIndex()]
+    );
 
-        wedgePolyPatchPtr_ = &wedge;
+    if (wedgePtr)
+    {
+        wedgePolyPatchPtr_ = wedgePtr;
     }
     else
     {
diff --git a/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C b/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C
index 52f7a3ba62..a050e497a8 100644
--- a/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C
+++ b/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -21,7 +21,7 @@ License
     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
+    You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 \*---------------------------------------------------------------------------*/
@@ -68,7 +68,7 @@ Foam::faPatch::faPatch
 :
     labelList(edgeLabels),
     patchIdentifier(name, index),
-    ngbPolyPatchIndex_(ngbPolyPatchIndex),
+    nbrPolyPatchId_(ngbPolyPatchIndex),
     boundaryMesh_(bm),
     edgeFacesPtr_(nullptr),
     pointLabelsPtr_(nullptr),
@@ -86,7 +86,7 @@ Foam::faPatch::faPatch
 :
     labelList(dict.get<labelList>("edgeLabels")),
     patchIdentifier(name, dict, index),
-    ngbPolyPatchIndex_(dict.get<label>("ngbPolyPatchIndex")),
+    nbrPolyPatchId_(dict.get<label>("ngbPolyPatchIndex")),
     boundaryMesh_(bm),
     edgeFacesPtr_(nullptr),
     pointLabelsPtr_(nullptr),
@@ -98,7 +98,7 @@ Foam::faPatch::faPatch(const faPatch& p, const faBoundaryMesh& bm)
 :
     labelList(p),
     patchIdentifier(p, p.index()),
-    ngbPolyPatchIndex_(p.ngbPolyPatchIndex_),
+    nbrPolyPatchId_(p.nbrPolyPatchId_),
     boundaryMesh_(bm),
     edgeFacesPtr_(nullptr),
     pointLabelsPtr_(nullptr),
@@ -116,13 +116,13 @@ Foam::faPatch::~faPatch()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::label Foam::faPatch::ngbPolyPatchIndex() const
+Foam::label Foam::faPatch::ngbPolyPatchIndex() const noexcept
 {
-    return ngbPolyPatchIndex_;
+    return nbrPolyPatchId_;
 }
 
 
-const Foam::faBoundaryMesh& Foam::faPatch::boundaryMesh() const
+const Foam::faBoundaryMesh& Foam::faPatch::boundaryMesh() const noexcept
 {
     return boundaryMesh_;
 }
@@ -240,38 +240,22 @@ const Foam::labelListList& Foam::faPatch::pointEdges() const
 
 Foam::labelList Foam::faPatch::ngbPolyPatchFaces() const
 {
-    labelList ngbFaces;
-
-    if (ngbPolyPatchIndex() == -1)
+    if (nbrPolyPatchId_ < 0)
     {
-        return ngbFaces;
+        return labelList();
     }
 
-    ngbFaces.setSize(faPatch::size());
+    labelList ngbFaces(faPatch::size());
 
     const faMesh& aMesh = boundaryMesh().mesh();
-    const polyMesh& pMesh = aMesh();
-    const indirectPrimitivePatch& patch = aMesh.patch();
+    const polyMesh& pMesh = aMesh.mesh();
+    const auto& patch = aMesh.patch();
 
     const labelListList& edgeFaces = pMesh.edgeFaces();
 
-    labelList faceCells(patch.size(), -1);
-
-    forAll(faceCells, faceI)
-    {
-        label faceID = aMesh.faceLabels()[faceI];
-
-        faceCells[faceI] = pMesh.faceOwner()[faceID];
-    }
-
-    labelList meshEdges
+    const labelList meshEdges
     (
-        patch.meshEdges
-        (
-            pMesh.edges(),
-            pMesh.cellEdges(),
-            faceCells
-        )
+        patch.meshEdges(pMesh.edges(), pMesh.pointEdges())
     );
 
     forAll(ngbFaces, edgeI)
@@ -288,7 +272,7 @@ Foam::labelList Foam::faPatch::ngbPolyPatchFaces() const
 
             label curPatchID = pMesh.boundaryMesh().whichPatch(curFace);
 
-            if (curPatchID == ngbPolyPatchIndex())
+            if (curPatchID == nbrPolyPatchId_)
             {
                 ngbFaces[edgeI] = curFace;
             }
@@ -308,14 +292,14 @@ Foam::labelList Foam::faPatch::ngbPolyPatchFaces() const
 
 Foam::tmp<Foam::vectorField> Foam::faPatch::ngbPolyPatchFaceNormals() const
 {
+    if (nbrPolyPatchId_ < 0)
+    {
+        return tmp<vectorField>::New();
+    }
+
     auto tfN = tmp<vectorField>::New();
     auto& fN = tfN.ref();
 
-    if (ngbPolyPatchIndex() == -1)
-    {
-        return tfN;
-    }
-
     fN.setSize(faPatch::size());
 
     labelList ngbFaces = ngbPolyPatchFaces();
@@ -336,7 +320,7 @@ Foam::tmp<Foam::vectorField> Foam::faPatch::ngbPolyPatchFaceNormals() const
 
 Foam::tmp<Foam::vectorField> Foam::faPatch::ngbPolyPatchPointNormals() const
 {
-    if (ngbPolyPatchIndex() == -1)
+    if (nbrPolyPatchId_ < 0)
     {
         return tmp<vectorField>::New();
     }
@@ -468,12 +452,17 @@ void Foam::faPatch::movePoints(const pointField& points)
 {}
 
 
-void Foam::faPatch::resetEdges(const labelList& newEdges)
+void Foam::faPatch::resetEdges(const UList<label>& newEdges)
 {
-    Info<< "Resetting patch edges" << endl;
-    labelList::operator=(newEdges);
-
     clearOut();
+    static_cast<labelList&>(*this) = newEdges;
+}
+
+
+void Foam::faPatch::resetEdges(labelList&& newEdges)
+{
+    clearOut();
+    static_cast<labelList&>(*this) = std::move(newEdges);
 }
 
 
@@ -483,9 +472,8 @@ void Foam::faPatch::write(Ostream& os) const
 
     patchIdentifier::write(os);
 
-    const labelList& edgeLabels = *this;
-    edgeLabels.writeEntry("edgeLabels", os);
-    os.writeEntry("ngbPolyPatchIndex", ngbPolyPatchIndex_);
+    os.writeEntry("ngbPolyPatchIndex", nbrPolyPatchId_);
+    static_cast<const labelList&>(*this).writeEntry("edgeLabels", os);
 }
 
 
diff --git a/src/finiteArea/faMesh/faPatches/faPatch/faPatch.H b/src/finiteArea/faMesh/faPatches/faPatch/faPatch.H
index 52fb33e613..98bb9efe76 100644
--- a/src/finiteArea/faMesh/faPatches/faPatch/faPatch.H
+++ b/src/finiteArea/faMesh/faPatches/faPatch/faPatch.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -58,6 +58,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward Declarations
 class faBoundaryMesh;
 class faPatch;
 Ostream& operator<<(Ostream&, const faPatch&);
@@ -74,7 +75,7 @@ class faPatch
     // Private Data
 
         //- Neighbour polyPatch index
-        const label ngbPolyPatchIndex_;
+        const label nbrPolyPatchId_;
 
         //- Reference to boundary mesh
         const faBoundaryMesh& boundaryMesh_;
@@ -100,7 +101,6 @@ class faPatch
         //- No copy assignment
         void operator=(const faPatch&) = delete;
 
-
         //- Clear out topological patch data
         void clearOut();
 
@@ -225,17 +225,31 @@ public:
 
     // Member Functions
 
-        //- Return number of patch points
+        //- Return the list of edges
+        const labelList& edgeLabels() const noexcept
+        {
+            return static_cast<const labelList&>(*this);
+        }
+
+        void edgeLabels(const UList<label>& newEdgeLabels);
+
+        //- Number of patch points
         label nPoints() const
         {
             return pointLabels().size();
         }
 
+        //- Number of edge labels (boundary edges) addressed by this patch
+        label nEdges() const noexcept
+        {
+            return labelList::size();
+        }
+
         //- Return neighbour polyPatch index
-        label ngbPolyPatchIndex() const;
+        label ngbPolyPatchIndex() const noexcept;
 
         //- Return boundaryMesh reference
-        const faBoundaryMesh& boundaryMesh() const;
+        const faBoundaryMesh& boundaryMesh() const noexcept;
 
         //- Return true if this patch is coupled
         virtual bool coupled() const
@@ -246,7 +260,7 @@ public:
         //- Patch start in edge list
         label start() const;
 
-        //- Patch size
+        //- Patch size is the number of edge labels
         virtual label size() const
         {
             return labelList::size();
@@ -330,8 +344,11 @@ public:
 
         // Topological changes
 
-            //- Reset edge list
-            void resetEdges(const labelList&);
+            //- Reset the list of edges (use with caution)
+            void resetEdges(const UList<label>& newEdges);
+
+            //- Reset the list of edges (use with caution)
+            void resetEdges(labelList&& newEdges);
 
 
         // Evaluation
diff --git a/src/finiteArea/faMesh/faPatches/faPatch/faPatchData.C b/src/finiteArea/faMesh/faPatches/faPatch/faPatchData.C
new file mode 100644
index 0000000000..7d712c99ef
--- /dev/null
+++ b/src/finiteArea/faMesh/faPatches/faPatch/faPatchData.C
@@ -0,0 +1,130 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "faPatchData.H"
+#include "edge.H"
+#include "dictionary.H"
+#include "processorFaPatch.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::faPatchData::faPatchData()
+:
+    ownerPolyPatchId_(-1),
+    neighPolyPatchId_(-1),
+    ownerProcId_(-1),
+    neighProcId_(-1)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::word& Foam::faPatchData::name() const noexcept
+{
+    return name_;
+}
+
+
+Foam::dictionary Foam::faPatchData::dict(const bool withEdgeLabels) const
+{
+    dictionary patchDict;
+    patchDict.add("type", type_);
+
+    if (withEdgeLabels)
+    {
+        patchDict.add("edgeLabels", edgeLabels_);
+    }
+    else
+    {
+        patchDict.add("edgeLabels", labelList());
+    }
+    patchDict.add("ngbPolyPatchIndex", neighPolyPatchId_);
+
+    if (coupled())
+    {
+        patchDict.add("myProcNo", ownerProcId_);
+        patchDict.add("neighbProcNo", neighProcId_);
+    }
+
+    return patchDict;
+}
+
+
+void Foam::faPatchData::clear()
+{
+    name_.clear();
+    type_.clear();
+
+    ownerPolyPatchId_ = -1;
+    neighPolyPatchId_ = -1;
+
+    ownerProcId_ = -1;
+    neighProcId_ = -1;
+
+    edgeLabels_.clear();
+}
+
+
+void Foam::faPatchData::assign(const faPatch& fap)
+{
+    clear();
+
+    // Copy information
+    name_ = fap.name();
+    type_ = fap.type();
+
+    neighPolyPatchId_ = fap.ngbPolyPatchIndex();
+    edgeLabels_ = fap.edgeLabels();
+
+    const auto* fapp = isA<processorFaPatch>(fap);
+    if (fapp)
+    {
+        ownerProcId_ = fapp->myProcNo();
+        neighProcId_ = fapp->neighbProcNo();
+    }
+}
+
+
+bool Foam::faPatchData::foundPatchPair(const edge& patchPair) const
+{
+    // Same as edge::compare
+    return
+    (
+        (
+            ownerPolyPatchId_ == patchPair.first()
+         && neighPolyPatchId_ == patchPair.second()
+        )
+     ||
+        (
+            ownerPolyPatchId_ == patchPair.second()
+         && neighPolyPatchId_ == patchPair.first()
+        )
+    );
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteArea/faMesh/faPatches/faPatch/faPatchData.H b/src/finiteArea/faMesh/faPatches/faPatch/faPatchData.H
index c8d4669791..4f785b3515 100644
--- a/src/finiteArea/faMesh/faPatches/faPatch/faPatchData.H
+++ b/src/finiteArea/faMesh/faPatches/faPatch/faPatchData.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,7 +28,8 @@ Class
     Foam::faPatchData
 
 Description
-    Class which holds data needed for faPatch construction
+    Helper class for holding data during faPatch construction.
+    Most data members are exposed at the moment.
 
 \*---------------------------------------------------------------------------*/
 
@@ -41,25 +43,81 @@ Description
 namespace Foam
 {
 
+// Forward Declarations
+class edge;
+class faPatch;
+class dictionary;
+
 /*---------------------------------------------------------------------------*\
                          Class faPatchData Declaration
 \*---------------------------------------------------------------------------*/
 
-struct faPatchData
+class faPatchData
 {
-    word name_;
-    word type_;
-    dictionary dict_;
-    label ownPolyPatchID_;
-    label ngbPolyPatchID_;
-    labelList edgeLabels_;
-    faPatchData()
-    :
-        name_(word::null),
-        type_(word::null),
-        ownPolyPatchID_(-1),
-        ngbPolyPatchID_(-1)
-    {}
+public:
+
+    // Data Members
+
+        word name_;
+        word type_;
+
+        label ownerPolyPatchId_;
+        label neighPolyPatchId_;
+
+        //- The owner/neighbour for processor patches
+        int ownerProcId_;
+        int neighProcId_;
+
+        // Storge (temporary or otherwise) for edge labels
+        labelList edgeLabels_;
+
+
+    // Constructors
+
+        //- Default construct
+        faPatchData();
+
+
+    // Member Functions
+
+    // Opaque read-only access
+
+        //- Return the name
+        const word& name() const noexcept;
+
+        //- Contents transcribed into a patch dictionary,
+        //- usually including the edge labels.
+        dictionary dict(const bool withEdgeLabels = true) const;
+
+
+    // Other Functions
+
+        //- Reset data
+        void clear();
+
+        //- Clear and populate with values from finiteArea patch
+        void assign(const faPatch& fap);
+
+        //- True if owner/neighbour processor ids are non-equal
+        bool coupled() const noexcept
+        {
+            return (ownerProcId_ != neighProcId_);
+        }
+
+        //- Does this side own the patch? Also true for non-coupled patches
+        bool owner() const noexcept
+        {
+            return (ownerProcId_ <= neighProcId_);
+        }
+
+        //- Does the other side own the patch?
+        bool neighbour() const noexcept
+        {
+            return !owner();
+        }
+
+        //- True it matches the owner/neighbour patch pair (any order)
+        bool foundPatchPair(const edge& patchPair) const;
 };
 
 
diff --git a/src/finiteArea/faMesh/faPatches/faPatch/faPatchFaMeshTemplates.C b/src/finiteArea/faMesh/faPatches/faPatch/faPatchFaMeshTemplates.C
index 6e3fab5ee8..ce4a007e11 100644
--- a/src/finiteArea/faMesh/faPatches/faPatch/faPatchFaMeshTemplates.C
+++ b/src/finiteArea/faMesh/faPatches/faPatch/faPatchFaMeshTemplates.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 Wikki Ltd
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -39,10 +40,8 @@ const typename GeometricField::Patch& Foam::faPatch::lookupPatchField
 {
     return patchField<GeometricField, Type>
     (
-        boundaryMesh().mesh()().objectRegistry::lookupObject<GeometricField>
-        (
-            name
-        )
+        boundaryMesh().mesh().mesh().objectRegistry::template
+            lookupObject<GeometricField>(name)
     );
 }