Merge branch 'finiteArea-decompose' into 'develop'

support parallel creation of finiteArea meshes with on-the-fly decomposition of fields #2084

See merge request Development/openfoam!457

applications/utilities/finiteArea/makeFaMesh/Make/options

@@ -1,8 +1,13 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteArea/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/cfdTools/general/lnInclude
+    -I$(LIB_SRC)/fileFormats/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/parallel/decompose/faDecompose/lnInclude \
+    -I$(LIB_SRC)/parallel/reconstruct/faReconstruct/lnInclude
 
 EXE_LIBS = \
     -lfiniteArea \
-    -lfiniteVolume
+    -lfileFormats \
+    -lmeshTools \
+    -lfaDecompose \
+    -lfaReconstruct

applications/utilities/finiteArea/makeFaMesh/decomposeFaFields.H

@@ -0,0 +1,94 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+    Decompose area fields, when mesh was generated in parallel
+
+\*---------------------------------------------------------------------------*/
+
+if (Pstream::parRun())
+{
+    faMeshReconstructor reconstructor(areaMesh);
+    reconstructor.writeAddressing();
+
+    // Handle area fields
+    // ------------------
+
+    PtrList<areaScalarField> areaScalarFields;
+    PtrList<areaVectorField> areaVectorFields;
+    PtrList<areaSphericalTensorField> areaSphTensorFields;
+    PtrList<areaSymmTensorField> areaSymmTensorFields;
+    PtrList<areaTensorField> areaTensorFields;
+
+    const faMesh& fullMesh = reconstructor.mesh();
+
+    {
+        // Use uncollated (or master uncollated) file handler here.
+        // - each processor is reading in the identical serial fields.
+        // - nothing should be parallel-coordinated.
+
+        // Similarly, if we write the serial finite-area mesh, this is only
+        // done from one processor!
+
+        reconstructor.writeMesh();
+
+        const bool oldDistributed = fileHandler().distributed();
+        auto oldHandler = fileHandler(fileOperation::NewUncollated());
+        fileHandler().distributed(true);
+
+        IOobjectList objects(fullMesh.time(), runTime.timeName());
+
+        faFieldDecomposer::readFields(fullMesh, objects, areaScalarFields);
+        faFieldDecomposer::readFields(fullMesh, objects, areaVectorFields);
+        faFieldDecomposer::readFields(fullMesh, objects, areaSphTensorFields);
+        faFieldDecomposer::readFields(fullMesh, objects, areaSymmTensorFields);
+        faFieldDecomposer::readFields(fullMesh, objects, areaTensorFields);
+
+        // Restore old settings
+        if (oldHandler)
+        {
+            fileHandler(std::move(oldHandler));
+        }
+        fileHandler().distributed(oldDistributed);
+    }
+
+    const label nAreaFields =
+    (
+        areaScalarFields.size()
+      + areaVectorFields.size()
+      + areaSphTensorFields.size()
+      + areaSymmTensorFields.size()
+      + areaTensorFields.size()
+    );
+
+    if (nAreaFields)
+    {
+        Info<< "Decomposing " << nAreaFields << " area fields" << nl << endl;
+
+        faFieldDecomposer fieldDecomposer
+        (
+            fullMesh,
+            areaMesh,
+            reconstructor.edgeProcAddressing(),
+            reconstructor.faceProcAddressing(),
+            reconstructor.boundaryProcAddressing()
+        );
+
+        fieldDecomposer.decomposeFields(areaScalarFields);
+        fieldDecomposer.decomposeFields(areaVectorFields);
+        fieldDecomposer.decomposeFields(areaSphTensorFields);
+        fieldDecomposer.decomposeFields(areaSymmTensorFields);
+        fieldDecomposer.decomposeFields(areaTensorFields);
+    }
+}
+
+// ************************************************************************* //

applications/utilities/finiteArea/makeFaMesh/faMeshWriteEdgesOBJ.H

@@ -0,0 +1,43 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+    OBJ output of faMesh edges
+
+\*---------------------------------------------------------------------------*/
+
+{
+    Info<< "Writing edges in obj format" << endl;
+
+    word outputName("faMesh-edges.obj");
+
+    if (Pstream::parRun())
+    {
+        outputName = word
+        (
+            "faMesh-edges-" + Foam::name(Pstream::myProcNo()) + ".obj"
+        );
+    }
+
+    OBJstream os(runTime.globalPath()/outputName);
+
+    os.writeQuoted
+    (
+        ("# " + outputName + "\n"),
+        false
+    );
+
+    os.write(areaMesh.patch().edges(), areaMesh.patch().localPoints());
+}
+
+
+// ************************************************************************* //

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;
+
+
+// ************************************************************************* //

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.
@@ -28,6 +29,8 @@ Application
 
 Description
     A mesh generator for finiteArea mesh.
+    When called in parallel, it will also try to act like decomposePar,
+    create procAddressing and decompose serial finite-area fields.
 
 Author
     Zeljko Tukovic, FAMENA
@@ -35,36 +38,20 @@ Author
 
 \*---------------------------------------------------------------------------*/
 
-#include "objectRegistry.H"
 #include "Time.H"
 #include "argList.H"
 #include "OSspecific.H"
 #include "faMesh.H"
-#include "fvMesh.H"
+#include "IOdictionary.H"
+#include "IOobjectList.H"
+
+#include "areaFields.H"
+#include "faFieldDecomposer.H"
+#include "faMeshReconstructor.H"
+#include "OBJstream.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 +60,64 @@ 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");
+
+    argList::addBoolOption
+    (
+        "write-edges-obj",
+        "Write mesh edges as obj files and exit",
+        false  // could make an advanced option
+    );
 
     #include "addRegionOption.H"
-    argList::noParallel();
-
     #include "setRootCase.H"
     #include "createTime.H"
-    #include "createNamedMesh.H"
+    #include "createNamedPolyMesh.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);
+        meshDefDict.add("emptyPatch", patchName, true);
     }
 
-        const word neiName(curPatchDict.get<word>("neighbourPolyPatch"));
+    // Create
+    faMesh areaMesh(mesh, meshDefDict);
 
-        faPatches[patchI].ngbPolyPatchID_ =
-            mesh.boundaryMesh().findPatchID(neiName);
+    bool quickExit = false;
 
-        if (faPatches[patchI].ngbPolyPatchID_ < 0)
+    if (args.found("write-edges-obj"))
     {
-            FatalErrorIn("makeFaMesh:")
-                << "neighbourPolyPatch " << neiName << " does not exist"
-                << exit(FatalError);
-        }
+        quickExit = true;
+        #include "faMeshWriteEdgesOBJ.H"
     }
 
-    // Setting faceLabels list size
-    label size = 0;
-
-    labelList patchIDs(polyMeshPatches.size(), -1);
-
-    forAll(polyMeshPatches, patchI)
+    if (quickExit)
     {
-        patchIDs[patchI] =
-            mesh.boundaryMesh().findPatchID(polyMeshPatches[patchI]);
-
-        if (patchIDs[patchI] < 0)
-        {
-            FatalErrorIn("makeFaMesh:")
-                << "Patch " << polyMeshPatches[patchI] << " does not exist"
-                << exit(FatalError);
+        Info<< "\nEnd\n" << endl;
+        return 0;
     }
 
-        size += mesh.boundaryMesh()[patchIDs[patchI]].size();
-    }
+    // Set the precision of the points data to 10
+    IOstream::defaultPrecision(10);
 
-    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;
-
-    // Writing faMesh
-    Info << "Write finite area mesh ... ";
+    Info<< nl << "Write finite area mesh." << nl;
     areaMesh.write();
+    Info<< endl;
 
-    Info << "\nEnd" << endl;
+    #include "decomposeFaFields.H"
+
+    Info << "\nEnd\n" << endl;
 
     return 0;
 }

applications/utilities/parallelProcessing/decomposePar/Make/files

@@ -1,4 +1,5 @@
 decomposePar.C
+
 domainDecomposition.C
 domainDecompositionMesh.C
 domainDecompositionDistribute.C
@@ -7,10 +8,6 @@ domainDecompositionWrite.C
 domainDecompositionDryRun.C
 domainDecompositionDryRunWrite.C
 
-dimFieldDecomposer.C
-pointFieldDecomposer.C
-faMeshDecomposition.C
-faFieldDecomposer.C
 lagrangianFieldDecomposer.C
 
 EXE = $(FOAM_APPBIN)/decomposePar

applications/utilities/parallelProcessing/decomposePar/Make/options

@@ -6,7 +6,8 @@ EXE_INC = \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
-    -I$(LIB_SRC)/parallel/decompose/decompose/lnInclude
+    -I$(LIB_SRC)/parallel/decompose/decompose/lnInclude \
+    -I$(LIB_SRC)/parallel/decompose/faDecompose/lnInclude
 
 EXE_LIBS = \
     -lfiniteArea \
@@ -18,5 +19,6 @@ EXE_LIBS = \
     -lgenericPatchFields \
     -ldecompositionMethods \
     -ldecompose \
+    -lfaDecompose \
     -L$(FOAM_LIBBIN)/dummy \
     -lkahipDecomp -lmetisDecomp -lscotchDecomp

applications/utilities/parallelProcessing/decomposePar/decomposePar.C

@@ -162,7 +162,6 @@ Usage
 #include "regionProperties.H"
 
 #include "readFields.H"
-#include "dimFieldDecomposer.H"
 #include "fvFieldDecomposer.H"
 #include "pointFieldDecomposer.H"
 #include "lagrangianFieldDecomposer.H"
@@ -173,11 +172,52 @@ Usage
 #include "faMeshDecomposition.H"
 #include "faFieldDecomposer.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
+// Read proc addressing at specific instance.
+// Uses polyMesh/fvMesh meshSubDir by default
+autoPtr<labelIOList> procAddressing
+(
+    const fvMesh& procMesh,
+    const word& name,
+    const word& instance,
+    const word& local = fvMesh::meshSubDir
+)
+{
+    return autoPtr<labelIOList>::New
+    (
+        IOobject
+        (
+            name,
+            instance,
+            local,
+            procMesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false  // do not register
+        )
+    );
+}
+
+
+// Read proc addressing at specific instance.
+// Uses the finiteArea meshSubDir
+autoPtr<labelIOList> faProcAddressing
+(
+    const fvMesh& procMesh,
+    const word& name,
+    const word& instance,
+    const word& local = faMesh::meshSubDir
+)
+{
+    return procAddressing(procMesh, name, instance, local);
+}
+
+
+// Return cached or read proc addressing from facesInstance
 const labelIOList& procAddressing
 (
     const PtrList<fvMesh>& procMeshList,
@@ -193,19 +233,7 @@ const labelIOList& procAddressing
         procAddressingList.set
         (
             proci,
-            new labelIOList
-            (
-                IOobject
-                (
-                    name,
-                    procMesh.facesInstance(),
-                    procMesh.meshSubDir,
-                    procMesh,
-                    IOobject::MUST_READ,
-                    IOobject::NO_WRITE,
-                    false
-                )
-            )
+            procAddressing(procMesh, name, procMesh.facesInstance())
         );
     }
     return procAddressingList[proci];
@@ -275,7 +303,7 @@ void decomposeUniform
     }
 }
 
-}
+} // End namespace Foam
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -682,7 +710,6 @@ int main(int argc, char *argv[])
             PtrList<labelIOList> cellProcAddressingList(mesh.nProcs());
             PtrList<labelIOList> boundaryProcAddressingList(mesh.nProcs());
             PtrList<fvFieldDecomposer> fieldDecomposerList(mesh.nProcs());
-            PtrList<dimFieldDecomposer> dimFieldDecomposerList(mesh.nProcs());
             PtrList<labelIOList> pointProcAddressingList(mesh.nProcs());
             PtrList<pointFieldDecomposer> pointFieldDecomposerList
             (
@@ -1047,7 +1074,6 @@ int main(int argc, char *argv[])
                 {
                     Info<< "Processor " << proci << ": field transfer" << endl;
 
-
                     // open the database
                     if (!processorDbList.set(proci))
                     {
@@ -1112,7 +1138,7 @@ int main(int argc, char *argv[])
                     );
 
 
-                    // FV fields
+                    // FV fields: volume, surface, internal
                     {
                         if (!fieldDecomposerList.set(proci))
                         {
@@ -1132,6 +1158,7 @@ int main(int argc, char *argv[])
                         const fvFieldDecomposer& fieldDecomposer =
                             fieldDecomposerList[proci];
 
+                        // vol fields
                         fieldDecomposer.decomposeFields(volScalarFields);
                         fieldDecomposer.decomposeFields(volVectorFields);
                         fieldDecomposer.decomposeFields
@@ -1141,6 +1168,7 @@ int main(int argc, char *argv[])
                         fieldDecomposer.decomposeFields(volSymmTensorFields);
                         fieldDecomposer.decomposeFields(volTensorFields);
 
+                        // surface fields
                         fieldDecomposer.decomposeFields(surfaceScalarFields);
                         fieldDecomposer.decomposeFields(surfaceVectorFields);
                         fieldDecomposer.decomposeFields
@@ -1153,6 +1181,13 @@ int main(int argc, char *argv[])
                         );
                         fieldDecomposer.decomposeFields(surfaceTensorFields);
 
+                        // internal fields
+                        fieldDecomposer.decomposeFields(dimScalarFields);
+                        fieldDecomposer.decomposeFields(dimVectorFields);
+                        fieldDecomposer.decomposeFields(dimSphericalTensorFields);
+                        fieldDecomposer.decomposeFields(dimSymmTensorFields);
+                        fieldDecomposer.decomposeFields(dimTensorFields);
+
                         if (times.size() == 1)
                         {
                             // Clear cached decomposer
@@ -1160,37 +1195,6 @@ int main(int argc, char *argv[])
                         }
                     }
 
-                    // Dimensioned fields
-                    {
-                        if (!dimFieldDecomposerList.set(proci))
-                        {
-                            dimFieldDecomposerList.set
-                            (
-                                proci,
-                                new dimFieldDecomposer
-                                (
-                                    mesh,
-                                    procMesh,
-                                    faceProcAddressing,
-                                    cellProcAddressing
-                                )
-                            );
-                        }
-                        const dimFieldDecomposer& dimDecomposer =
-                            dimFieldDecomposerList[proci];
-
-                        dimDecomposer.decomposeFields(dimScalarFields);
-                        dimDecomposer.decomposeFields(dimVectorFields);
-                        dimDecomposer.decomposeFields(dimSphericalTensorFields);
-                        dimDecomposer.decomposeFields(dimSymmTensorFields);
-                        dimDecomposer.decomposeFields(dimTensorFields);
-
-                        if (times.size() == 1)
-                        {
-                            dimFieldDecomposerList.set(proci, nullptr);
-                        }
-                    }
-
 
                     // Point fields
                     if
@@ -1366,18 +1370,24 @@ int main(int argc, char *argv[])
                     faMesh::meshSubDir,
                     mesh,
                     IOobject::READ_IF_PRESENT,
-                    IOobject::NO_WRITE
+                    IOobject::NO_WRITE,
+                    false  // not registered
                 );
 
 
                 if (faMeshBoundaryIOobj.typeHeaderOk<faBoundaryMesh>(true))
                 {
-                    Info << "\nFinite area mesh decomposition" << endl;
+                    Info<< "\nFinite area mesh decomposition" << endl;
 
-                    faMeshDecomposition aMesh(mesh);
+                    // Always based on the volume decomposition!
+                    faMeshDecomposition aMesh
+                    (
+                        mesh,
+                        mesh.nProcs(),
+                        mesh.model()
+                    );
 
                     aMesh.decomposeMesh();
-
                     aMesh.writeDecomposition();
 
 
@@ -1404,13 +1414,29 @@ int main(int argc, char *argv[])
                     PtrList<edgeScalarField> edgeScalarFields;
                     readFields(aMesh, objects, edgeScalarFields);
 
-                    Info << endl;
+                    const label nAreaFields =
+                    (
+                        areaScalarFields.size()
+                      + areaVectorFields.size()
+                      + areaSphericalTensorFields.size()
+                      + areaSymmTensorFields.size()
+                      + areaTensorFields.size()
+                      + edgeScalarFields.size()
+                    );
+
+                    Info<< endl;
+                    Info<< "Finite area field transfer: "
+                        << nAreaFields << " fields" << endl;
 
                     // Split the fields over processors
-                    for (label procI = 0; procI < mesh.nProcs(); procI++)
+                    for
+                    (
+                        label proci = 0;
+                        nAreaFields && proci < mesh.nProcs();
+                        ++proci
+                    )
                     {
-                        Info<< "Processor " << procI
-                            << ": finite area field transfer" << endl;
+                        Info<< "    Processor " << proci << endl;
 
                         // open the database
                         Time processorDb
@@ -1418,7 +1444,7 @@ int main(int argc, char *argv[])
                             Time::controlDictName,
                             args.rootPath(),
                             args.caseName()
-                          / ("processor" + Foam::name(procI))
+                          / ("processor" + Foam::name(proci))
                         );
 
                         processorDb.setTime(runTime);
@@ -1441,7 +1467,7 @@ int main(int argc, char *argv[])
                         //     procAddressing
                         //     (
                         //         procMeshList,
-                        //         procI,
+                        //         proci,
                         //         "faceProcAddressing",
                         //         faceProcAddressingList
                         //     );
@@ -1450,60 +1476,39 @@ int main(int argc, char *argv[])
                         //     procAddressing
                         //     (
                         //         procMeshList,
-                        //         procI,
+                        //         proci,
                         //         "boundaryProcAddressing",
                         //         boundaryProcAddressingList
                         //     );
 
-                        labelIOList faceProcAddressing
-                        (
-                            IOobject
+                        // Addressing from faMesh (not polyMesh) meshSubDir
+
+                        autoPtr<labelIOList> tfaceProcAddr =
+                            faProcAddressing
                             (
+                                procFvMesh,
                                 "faceProcAddressing",
-                                "constant",
-                                procMesh.meshSubDir,
-                                procFvMesh,
-                                IOobject::MUST_READ,
-                                IOobject::NO_WRITE
-                            )
+                                runTime.constant()
                             );
+                        auto& faceProcAddressing = *tfaceProcAddr;
 
-                        labelIOList boundaryProcAddressing
-                        (
-                            IOobject
+                        autoPtr<labelIOList> tboundaryProcAddr =
+                            faProcAddressing
                             (
+                                procFvMesh,
                                 "boundaryProcAddressing",
-                                "constant",
-                                procMesh.meshSubDir,
-                                procFvMesh,
-                                IOobject::MUST_READ,
-                                IOobject::NO_WRITE
-                            )
+                                runTime.constant()
                             );
+                        auto& boundaryProcAddressing = *tboundaryProcAddr;
 
-                        // FA fields
-                        if
+                        autoPtr<labelIOList> tedgeProcAddr =
+                            faProcAddressing
                             (
-                            areaScalarFields.size()
-                         || areaVectorFields.size()
-                         || areaSphericalTensorFields.size()
-                         || areaSymmTensorFields.size()
-                         || areaTensorFields.size()
-                         || edgeScalarFields.size()
-                        )
-                        {
-                            labelIOList edgeProcAddressing
-                            (
-                                IOobject
-                                (
-                                    "edgeProcAddressing",
-                                    "constant",
-                                    procMesh.meshSubDir,
                                 procFvMesh,
-                                    IOobject::MUST_READ,
-                                    IOobject::NO_WRITE
-                                )
+                                "edgeProcAddressing",
+                                runTime.constant()
                             );
+                        const auto& edgeProcAddressing = *tedgeProcAddr;
 
                         faFieldDecomposer fieldDecomposer
                         (
@@ -1520,10 +1525,7 @@ int main(int argc, char *argv[])
                         (
                             areaSphericalTensorFields
                         );
-                            fieldDecomposer.decomposeFields
-                            (
-                                areaSymmTensorFields
-                            );
+                        fieldDecomposer.decomposeFields(areaSymmTensorFields);
                         fieldDecomposer.decomposeFields(areaTensorFields);
 
                         fieldDecomposer.decomposeFields(edgeScalarFields);
@@ -1532,7 +1534,6 @@ int main(int argc, char *argv[])
             }
         }
     }
-    }
 
     Info<< "\nEnd\n" << endl;
 

applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C

@@ -120,16 +120,27 @@ Foam::domainDecomposition::domainDecomposition
     procProcessorPatchSubPatchIDs_(nProcs_),
     procProcessorPatchSubPatchStarts_(nProcs_)
 {
-    decompositionModel::New
-    (
-        *this,
-        decompDictFile
-    ).readIfPresent("distributed", distributed_);
+    updateParameters(this->model());
 }
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+const Foam::decompositionModel& Foam::domainDecomposition::model() const
+{
+    return decompositionModel::New(*this, decompDictFile_);
+}
+
+
+void Foam::domainDecomposition::updateParameters
+(
+    const dictionary& params
+)
+{
+    params.readIfPresent("distributed", distributed_);
+}
+
+
 bool Foam::domainDecomposition::writeDecomposition(const bool decomposeSets)
 {
     Info<< "\nConstructing processor meshes" << endl;
@@ -408,10 +419,9 @@ bool Foam::domainDecomposition::writeDecomposition(const bool decomposeSets)
             nInterProcPatches += curSubPatchIDs[procPatchi].size();
         }
 
-        List<polyPatch*> procPatches
+        PtrList<polyPatch> procPatches
         (
-            curPatchSizes.size() + nInterProcPatches,
-            reinterpret_cast<polyPatch*>(0)
+            curPatchSizes.size() + nInterProcPatches
         );
 
         label nPatches = 0;
@@ -434,13 +444,17 @@ bool Foam::domainDecomposition::writeDecomposition(const bool decomposeSets)
             );
 
             // Map existing patches
-            procPatches[nPatches] = meshPatch.clone
+            procPatches.set
+            (
+                nPatches,
+                meshPatch.clone
                 (
                     procMesh.boundaryMesh(),
                     nPatches,
                     patchMapper.directAddressing(),
                     curPatchStarts[patchi]
-            ).ptr();
+                )
+            );
 
             nPatches++;
         }
@@ -464,7 +478,9 @@ bool Foam::domainDecomposition::writeDecomposition(const bool decomposeSets)
                 if (subPatchID[i] == -1)
                 {
                     // From internal faces
-                    procPatches[nPatches] =
+                    procPatches.set
+                    (
+                        nPatches,
                         new processorPolyPatch
                         (
                             size,
@@ -473,6 +489,7 @@ bool Foam::domainDecomposition::writeDecomposition(const bool decomposeSets)
                             procMesh.boundaryMesh(),
                             proci,
                             curNeighbourProcessors[procPatchi]
+                        )
                     );
                 }
                 else
@@ -483,7 +500,9 @@ bool Foam::domainDecomposition::writeDecomposition(const bool decomposeSets)
                               boundaryMesh()[subPatchID[i]]
                           );
 
-                    procPatches[nPatches] =
+                    procPatches.set
+                    (
+                        nPatches,
                         new processorCyclicPolyPatch
                         (
                             size,
@@ -494,12 +513,12 @@ bool Foam::domainDecomposition::writeDecomposition(const bool decomposeSets)
                             curNeighbourProcessors[procPatchi],
                             pcPatch.name(),
                             pcPatch.transform()
+                        )
                     );
                 }
 
                 curStart += size;
-
-                nPatches++;
+                ++nPatches;
             }
         }
 

applications/utilities/parallelProcessing/decomposePar/domainDecompositionDistribute.C

@@ -26,12 +26,8 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "domainDecomposition.H"
-#include "decompositionMethod.H"
 #include "cpuTime.H"
-#include "cellSet.H"
-#include "regionSplit.H"
-#include "Tuple2.H"
-#include "faceSet.H"
+#include "decompositionMethod.H"
 #include "decompositionModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/utilities/parallelProcessing/decomposePar/faFieldDecomposer.C

@@ -1,238 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2016-2017 Wikki 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 "faFieldDecomposer.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-faFieldDecomposer::patchFieldDecomposer::patchFieldDecomposer
-(
-    const label sizeBeforeMapping,
-    const labelUList& addressingSlice,
-    const label addressingOffset
-)
-:
-    sizeBeforeMapping_(sizeBeforeMapping),
-    directAddressing_(addressingSlice)
-{
-    forAll(directAddressing_, i)
-    {
-        // Subtract one to align addressing.
-        // directAddressing_[i] -= addressingOffset + 1;
-        // ZT, 12/Nov/2010
-        directAddressing_[i] -= addressingOffset;
-    }
-}
-
-
-faFieldDecomposer::processorAreaPatchFieldDecomposer::
-processorAreaPatchFieldDecomposer
-(
-    const faMesh& mesh,
-    const labelUList& addressingSlice
-)
-:
-    sizeBeforeMapping_(mesh.nFaces()),
-    addressing_(addressingSlice.size()),
-    weights_(addressingSlice.size())
-{
-    const scalarField& weights = mesh.weights().internalField();
-    const labelList& own = mesh.edgeOwner();
-    const labelList& neighb = mesh.edgeNeighbour();
-
-    forAll(addressing_, i)
-    {
-        // Subtract one to align addressing.
-        label ai = addressingSlice[i];
-//         label ai = mag(addressingSlice[i]) - 1;
-
-        if (ai < neighb.size())
-        {
-            // This is a regular edge. it has been an internal edge
-            // of the original mesh and now it has become a edge
-            // on the parallel boundary
-            addressing_[i].setSize(2);
-            weights_[i].setSize(2);
-
-            addressing_[i][0] = own[ai];
-            addressing_[i][1] = neighb[ai];
-
-            weights_[i][0] = weights[ai];
-            weights_[i][1] = 1.0 - weights[ai];
-        }
-        else
-        {
-            // This is a edge that used to be on a cyclic boundary
-            // but has now become a parallel patch edge. I cannot
-            // do the interpolation properly (I would need to look
-            // up the different (edge) list of data), so I will
-            // just grab the value from the owner face
-            //
-            addressing_[i].setSize(1);
-            weights_[i].setSize(1);
-
-            addressing_[i][0] = own[ai];
-
-            weights_[i][0] = 1.0;
-        }
-    }
-}
-
-
-faFieldDecomposer::processorEdgePatchFieldDecomposer::
-processorEdgePatchFieldDecomposer
-(
-    label sizeBeforeMapping,
-    const labelUList& addressingSlice
-)
-:
-    sizeBeforeMapping_(sizeBeforeMapping),
-    addressing_(addressingSlice.size()),
-    weights_(addressingSlice.size())
-{
-    forAll(addressing_, i)
-    {
-        addressing_[i].setSize(1);
-        weights_[i].setSize(1);
-
-        addressing_[i][0] = mag(addressingSlice[i]) - 1;
-        weights_[i][0] = sign(addressingSlice[i]);
-    }
-}
-
-
-faFieldDecomposer::faFieldDecomposer
-(
-    const faMesh& completeMesh,
-    const faMesh& procMesh,
-    const labelList& edgeAddressing,
-    const labelList& faceAddressing,
-    const labelList& boundaryAddressing
-)
-:
-    completeMesh_(completeMesh),
-    procMesh_(procMesh),
-    edgeAddressing_(edgeAddressing),
-    faceAddressing_(faceAddressing),
-    boundaryAddressing_(boundaryAddressing),
-    patchFieldDecomposerPtrs_
-    (
-        procMesh_.boundary().size(),
-        static_cast<patchFieldDecomposer*>(NULL)
-    ),
-    processorAreaPatchFieldDecomposerPtrs_
-    (
-        procMesh_.boundary().size(),
-        static_cast<processorAreaPatchFieldDecomposer*>(NULL)
-    ),
-    processorEdgePatchFieldDecomposerPtrs_
-    (
-        procMesh_.boundary().size(),
-        static_cast<processorEdgePatchFieldDecomposer*>(NULL)
-    )
-{
-    forAll(boundaryAddressing_, patchi)
-    {
-        if (boundaryAddressing_[patchi] >= 0)
-        {
-            patchFieldDecomposerPtrs_[patchi] = new patchFieldDecomposer
-            (
-                completeMesh_.boundary()[boundaryAddressing_[patchi]].size(),
-                procMesh_.boundary()[patchi].patchSlice(edgeAddressing_),
-//                 completeMesh_.boundaryMesh()
-                completeMesh_.boundary()
-                [
-                    boundaryAddressing_[patchi]
-                ].start()
-            );
-        }
-        else
-        {
-            processorAreaPatchFieldDecomposerPtrs_[patchi] =
-                new processorAreaPatchFieldDecomposer
-                (
-                    completeMesh_,
-                    procMesh_.boundary()[patchi].patchSlice(edgeAddressing_)
-                );
-
-            processorEdgePatchFieldDecomposerPtrs_[patchi] =
-                new processorEdgePatchFieldDecomposer
-                (
-                    procMesh_.boundary()[patchi].size(),
-                    static_cast<const labelUList&>
-                    (
-                        procMesh_.boundary()[patchi].patchSlice
-                        (
-                            edgeAddressing_
-                        )
-                    )
-                );
-        }
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-faFieldDecomposer::~faFieldDecomposer()
-{
-    forAll(patchFieldDecomposerPtrs_, patchi)
-    {
-        if (patchFieldDecomposerPtrs_[patchi])
-        {
-            delete patchFieldDecomposerPtrs_[patchi];
-        }
-    }
-
-    forAll(processorAreaPatchFieldDecomposerPtrs_, patchi)
-    {
-        if (processorAreaPatchFieldDecomposerPtrs_[patchi])
-        {
-            delete processorAreaPatchFieldDecomposerPtrs_[patchi];
-        }
-    }
-
-    forAll(processorEdgePatchFieldDecomposerPtrs_, patchi)
-    {
-        if (processorEdgePatchFieldDecomposerPtrs_[patchi])
-        {
-            delete processorEdgePatchFieldDecomposerPtrs_[patchi];
-        }
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //

applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposer.H

@@ -31,7 +31,7 @@ Description
 
 SourceFiles
     lagrangianFieldDecomposer.C
-    lagrangianFieldDecomposerDecomposeFields.C
+    lagrangianFieldDecomposerFields.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -84,7 +84,7 @@ public:
         //- Construct from components
         lagrangianFieldDecomposer
         (
-            const polyMesh& mesh,
+            const polyMesh& mesh,       //<! unused
             const polyMesh& procMesh,
             const labelList& faceProcAddressing,
             const labelList& cellProcAddressing,
@@ -156,7 +156,7 @@ public:
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-    #include "lagrangianFieldDecomposerDecomposeFields.C"
+    #include "lagrangianFieldDecomposerFields.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/utilities/parallelProcessing/decomposePar/readFields.C

@@ -42,7 +42,7 @@ void Foam::readFields
     typedef GeometricField<Type, PatchField, GeoMesh> GeoField;
 
     // Search list of objects for fields of type GeoField
-    IOobjectList fieldObjects(objects.lookupClass(GeoField::typeName));
+    IOobjectList fieldObjects(objects.lookupClass<GeoField>());
 
     // Remove the cellDist field
     auto iter = fieldObjects.find("cellDist");
@@ -51,12 +51,12 @@ void Foam::readFields
         fieldObjects.erase(iter);
     }
 
-    // Get sorted set of names (different processors might read objects in
-    // different order)
+    // Use sorted set of names
+    // (different processors might read objects in different order)
     const wordList masterNames(fieldObjects.sortedNames());
 
     // Construct the fields
-    fields.setSize(masterNames.size());
+    fields.resize(masterNames.size());
 
     forAll(masterNames, i)
     {
@@ -76,17 +76,14 @@ void Foam::readFields
 )
 {
     // Search list of objects for fields of type GeomField
-    IOobjectList fieldObjects(objects.lookupClass(GeoField::typeName));
+    IOobjectList fieldObjects(objects.lookupClass<GeoField>());
 
-    // Construct the fields
-    fields.setSize(fieldObjects.size());
-
-    // Get sorted set of names (different processors might read objects in
-    // different order)
+    // Use sorted set of names
+    // (different processors might read objects in different order)
     const wordList masterNames(fieldObjects.sortedNames());
 
     // Construct the fields
-    fields.setSize(masterNames.size());
+    fields.resize(masterNames.size());
 
     forAll(masterNames, i)
     {

applications/utilities/parallelProcessing/reconstructPar/Make/files

@@ -1,5 +1,3 @@
-processorFaMeshes.C
-faFieldReconstructor.C
 reconstructPar.C
 
 EXE = $(FOAM_APPBIN)/reconstructPar

applications/utilities/parallelProcessing/reconstructPar/Make/options

@@ -4,7 +4,8 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/parallel/reconstruct/reconstruct/lnInclude
+    -I$(LIB_SRC)/parallel/reconstruct/reconstruct/lnInclude \
+    -I$(LIB_SRC)/parallel/reconstruct/faReconstruct/lnInclude
 
 EXE_LIBS = \
     -lfiniteArea \
@@ -13,4 +14,5 @@ EXE_LIBS = \
     -llagrangian \
     -lgenericPatchFields \
     -ldynamicMesh \
-    -lreconstruct
+    -lreconstruct \
+    -lfaReconstruct

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

bin/tools/CleanFunctions

@@ -6,7 +6,7 @@
 #    \\/     M anipulation  |
 #------------------------------------------------------------------------------
 #     Copyright (C) 2011-2016 OpenFOAM Foundation
-#     Copyright (C) 2015-2020 OpenCFD Ltd.
+#     Copyright (C) 2015-2021 OpenCFD Ltd.
 #------------------------------------------------------------------------------
 # License
 #     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -138,15 +138,16 @@ cleanCase()
         then
             rm -f polyMesh/blockMeshDict
             echo
-            echo "Warning: not removing constant/polyMesh/ "
-            echo "   it contains a blockMeshDict, which should normally be under system/ instead"
+            echo "Warning: not removing constant/polyMesh/"
+            echo "    It contains a 'blockMeshDict.m4' file."
+            echo "    Relocate the file to system/ to avoid this warning"
             echo
         elif [ -e polyMesh/blockMeshDict ]
         then
             echo
-            echo "Warning: not removing constant/polyMesh/ "
-            echo "   it contains a blockMeshDict, which should normally be under system/ instead"
-            echo
+            echo "Warning: not removing constant/polyMesh/"
+            echo "    It contains a 'blockMeshDict' file."
+            echo "    Relocate the file to system/ to avoid this warning"
         else
             # Remove polyMesh entirely if there is no blockMeshDict
             rm -rf polyMesh
@@ -187,12 +188,26 @@ cleanUcomponents()
 }
 
 
-cleanFaMesh ()
+cleanFaMesh()
 {
-    rm -rf \
-        constant/faMesh/faceLabels* \
-        constant/faMesh/faBoundary* \
-    ;
+    (
+        cd constant 2>/dev/null || exit 0
+
+        # Old constant/polyMesh location for blockMeshDict still in use?
+        # - emit a gentle warning
+        if [ -e faMesh/faMeshDefinition ]
+        then
+            rm -f faMesh/faceLabels* faMesh/faBoundary*
+            echo
+            echo "Warning: not removing the constant/faMesh/ directory"
+            echo "    It contains a 'faMeshDefinition' file"
+            echo "    Relocate the file to system/ to avoid this warning"
+            echo
+        else
+            # Remove faMesh/ entirely if there is no faMeshDefinition
+            rm -rf faMesh
+        fi
+    )
 }
 
 

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,

src/faOptions/faOption/faOption.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -40,17 +40,6 @@ namespace Foam
 }
 
 
-// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
-
-void Foam::fa::option::constructMeshObjects()
-{
-    regionMeshPtr_.reset(new faMesh(mesh_));
-
-    vsmPtr_.reset(new volSurfaceMapping(regionMeshPtr_()));
-}
-
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::fa::option::option
@@ -67,16 +56,15 @@ Foam::fa::option::option
     patch_(patch),
     dict_(dict),
     coeffs_(dict.optionalSubDict(modelType + "Coeffs")),
-    active_(dict.getOrDefault<Switch>("active", true)),
     fieldNames_(),
     applied_(),
     regionName_(dict.get<word>("region")),
     regionMeshPtr_(nullptr),
-    vsmPtr_(nullptr)
+    vsmPtr_(nullptr),
+    active_(dict.getOrDefault("active", true)),
+    log(true)
 {
-    constructMeshObjects();
-
-    Info<< incrIndent << indent << "Source: " << name_ << endl << decrIndent;
+    Log << incrIndent << indent << "Source: " << name_ << endl << decrIndent;
 }
 
 

src/faOptions/faOption/faOption.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -75,13 +75,12 @@ SourceFiles
 #ifndef faOption_H
 #define faOption_H
 
-#include "faMatrices.H"
-#include "areaFields.H"
+#include "faMatricesFwd.H"
+#include "areaFieldsFwd.H"
 #include "dictionary.H"
-#include "Switch.H"
-#include "runTimeSelectionTables.H"
 #include "fvMesh.H"
 #include "volSurfaceMapping.H"
+#include "runTimeSelectionTables.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -124,27 +123,35 @@ protected:
         //- Dictionary containing source coefficients
         dictionary coeffs_;
 
-        //- Source active flag
-        Switch active_;
-
         //- Field names to apply source to - populated by derived models
         wordList fieldNames_;
 
         //- Applied flag list - corresponds to each fieldNames_ entry
         List<bool> applied_;
 
-         //- Region name
+        //- Region name (finite-area)
         word regionName_;
 
-        //- Pointer to the region mesh database
-        autoPtr<faMesh> regionMeshPtr_;
+private:
 
-        //-Volume-to surface mapping
-        autoPtr<volSurfaceMapping> vsmPtr_;
+    // Private Data
+
+        //- Demand-driven: pointer to region mesh database
+        mutable autoPtr<faMesh> regionMeshPtr_;
+
+        //- Demand-driven: volume-to-surface mapping
+        mutable autoPtr<volSurfaceMapping> vsmPtr_;
+
+        //- Source active flag
+        bool active_;
 
 
 public:
 
+    //- Switch write log to Info
+    bool log;
+
+
     //- Runtime type information
     TypeName("option");
 
@@ -188,7 +195,6 @@ public:
         //- on the freestore from an Istream
         class iNew
         {
-
             //- Reference to the patch
             const fvPatch& patch_;
 
@@ -239,37 +245,37 @@ public:
         // Access
 
             //- Return const access to the source name
-            inline const word& name() const;
+            inline const word& name() const noexcept;
 
             //- Return const access to the mesh database
-            inline const fvMesh& mesh() const;
+            inline const fvMesh& mesh() const noexcept;
 
             //- Return const access to fvPatch
-            inline const fvPatch& patch() const;
+            inline const fvPatch& patch() const noexcept;
 
             //- Return dictionary
-            inline const dictionary& coeffs() const;
+            inline const dictionary& coeffs() const noexcept;
 
             //- Return const access to the source active flag
-            inline bool active() const;
+            inline bool active() const noexcept;
 
             //- Set the applied flag to true for field index fieldi
             inline void setApplied(const label fieldi);
 
-            //- Return the region mesh database
+            //- The region name
+            inline const word& regionName() const noexcept;
+
+            //- Return the region mesh database (demand-driven)
             inline const faMesh& regionMesh() const;
 
-             //- Return volSurfaceMapping
+            //- Return volSurfaceMapping (demand-driven)
             inline const volSurfaceMapping& vsm() const;
 
-            //- Region name
-            inline const word& regionName() const;
-
 
         // Edit
 
-            //- Return access to the source active flag
-            inline Switch& active();
+            //- Change source active flag, return previous value
+            inline bool active(const bool on) noexcept;
 
 
         // Checks

src/faOptions/faOption/faOptionI.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 ------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,49 +27,51 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline const Foam::word& Foam::fa::option::name() const
+inline const Foam::word& Foam::fa::option::name() const noexcept
 {
     return name_;
 }
 
 
-inline const Foam::fvMesh& Foam::fa::option::mesh() const
+inline const Foam::fvMesh& Foam::fa::option::mesh() const noexcept
 {
     return mesh_;
 }
 
 
-inline const Foam::fvPatch& Foam::fa::option::patch() const
+inline const Foam::fvPatch& Foam::fa::option::patch() const noexcept
 {
     return patch_;
 }
 
 
-inline const Foam::dictionary& Foam::fa::option::coeffs() const
+inline const Foam::dictionary& Foam::fa::option::coeffs() const noexcept
 {
     return coeffs_;
 }
 
 
-inline bool Foam::fa::option::active() const
+inline bool Foam::fa::option::active() const noexcept
 {
     return active_;
 }
 
 
+inline bool Foam::fa::option::active(const bool on) noexcept
+{
+    bool old(active_);
+    active_ = on;
+    return old;
+}
+
+
 inline void Foam::fa::option::setApplied(const label fieldi)
 {
     applied_[fieldi] = true;
 }
 
 
-inline Foam::Switch& Foam::fa::option::active()
-{
-    return active_;
-}
-
-
-inline const Foam::word& Foam::fa::option::regionName() const
+inline const Foam::word& Foam::fa::option::regionName() const noexcept
 {
     return regionName_;
 }
@@ -77,31 +79,21 @@ inline const Foam::word& Foam::fa::option::regionName() const
 
 inline const Foam::faMesh& Foam::fa::option::regionMesh() const
 {
-    if (regionMeshPtr_.valid())
+    if (!regionMeshPtr_)
     {
-        return regionMeshPtr_();
+        regionMeshPtr_.reset(new faMesh(mesh_));
     }
-    else
-    {
-        FatalErrorInFunction
-            << "Region mesh not available" << abort(FatalError);
-    }
-    return *(new faMesh(mesh_));
+    return *regionMeshPtr_;
 }
 
 
 inline const Foam::volSurfaceMapping& Foam::fa::option::vsm() const
 {
-    if (vsmPtr_.valid())
+    if (!vsmPtr_)
     {
-        return vsmPtr_();
+        vsmPtr_.reset(new volSurfaceMapping(this->regionMesh()));
     }
-    else
-    {
-        FatalErrorInFunction
-            << "vsmPtr not available" << abort(FatalError);
-    }
-    return *(new volSurfaceMapping(regionMeshPtr_()));
+    return *vsmPtr_;
 }
 
 

src/faOptions/faOption/faOptionList.C

@@ -92,7 +92,7 @@ Foam::fa::optionList::optionList
     const dictionary& dict
 )
 :
-    PtrList<option>(),
+    PtrList<fa::option>(),
     mesh_(p.boundaryMesh().mesh()),
     patch_(p),
     checkTimeIndex_(mesh_.time().startTimeIndex() + 2)
@@ -103,7 +103,7 @@ Foam::fa::optionList::optionList
 
 Foam::fa::optionList::optionList(const fvPatch& p)
 :
-    PtrList<option>(),
+    PtrList<fa::option>(),
     mesh_(p.boundaryMesh().mesh()),
     patch_(p),
     checkTimeIndex_(mesh_.time().startTimeIndex() + 2)

src/faOptions/faOption/faOptionList.H

@@ -46,10 +46,9 @@ SourceFile
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+// Forward Declarations
 namespace Foam
 {
-
-// Forward declaration of friend functions and operators
 namespace fa
 {
     class optionList;
@@ -66,7 +65,7 @@ namespace fa
 
 class optionList
 :
-    public PtrList<option>
+    public PtrList<fa::option>
 {
 protected:
 
@@ -118,7 +117,7 @@ public:
 
     // Constructors
 
-        //- Construct null
+        //- Construct from patch
         optionList(const fvPatch& p);
 
         //- Construct from mesh and dictionary
@@ -126,8 +125,7 @@ public:
 
 
     //- Destructor
-    virtual ~optionList()
-    {}
+    virtual ~optionList() = default;
 
 
     // Member Functions

src/faOptions/faOption/faOptionListTemplates.C

@@ -26,6 +26,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "profiling.H"
+#include "areaFields.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 

src/faOptions/faceSetOption/faceSetOption.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -102,7 +102,7 @@ namespace fa
 
 class faceSetOption
 :
-    public option
+    public fa::option
 {
 public:
 

src/faOptions/sources/derived/contactHeatFluxSource/contactHeatFluxSource.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -98,7 +98,7 @@ namespace fa
 
 class contactHeatFluxSource
 :
-    public faceSetOption,
+    public fa::faceSetOption,
     public temperatureCoupledBase
 {
     // Private Data

src/faOptions/sources/derived/externalFileSource/externalFileSource.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -26,8 +26,8 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "externalFileSource.H"
-#include "faMatrices.H"
-#include "faCFD.H"
+#include "fam.H"
+#include "faScalarMatrix.H"
 #include "zeroGradientFaPatchFields.H"
 #include "addToRunTimeSelectionTable.H"
 

src/faOptions/sources/derived/externalFileSource/externalFileSource.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -89,7 +89,7 @@ namespace fa
 
 class externalFileSource
 :
-    public faceSetOption
+    public fa::faceSetOption
 {
     // Private Data
 

src/faOptions/sources/derived/externalHeatFluxSource/externalHeatFluxSource.C

@@ -26,7 +26,8 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "externalHeatFluxSource.H"
-#include "addToRunTimeSelectionTable.H"
+#include "fam.H"
+#include "faScalarMatrix.H"
 #include "physicoChemicalConstants.H"
 #include "zeroGradientFaPatchFields.H"
 #include "addToRunTimeSelectionTable.H"

src/faOptions/sources/derived/externalHeatFluxSource/externalHeatFluxSource.H

@@ -121,7 +121,6 @@ SourceFiles
 #include "Function1.H"
 #include "areaFields.H"
 #include "faceSetOption.H"
-#include "faCFD.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -136,7 +135,7 @@ namespace fa
 
 class externalHeatFluxSource
 :
-    public faceSetOption
+    public fa::faceSetOption
 {
 public:
 

src/faOptions/sources/derived/jouleHeatingSource/jouleHeatingSource.C

@@ -26,8 +26,8 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "jouleHeatingSource.H"
-#include "faMatrices.H"
-#include "faCFD.H"
+#include "fam.H"
+#include "faScalarMatrix.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //

src/faOptions/sources/derived/jouleHeatingSource/jouleHeatingSource.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -156,7 +156,7 @@ namespace fa
 
 class jouleHeatingSource
 :
-    public faceSetOption
+    public fa::faceSetOption
 {
     // Private Data
 

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

src/finiteArea/faMatrices/faMatricesFwd.H

@@ -0,0 +1,59 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Description
+    Forward declarations of standard faMatrix types/specializations.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef faMatricesFwd_H
+#define faMatricesFwd_H
+
+#include "fieldTypes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type> class faMatrix;
+
+typedef faMatrix<scalar> faScalarMatrix;
+typedef faMatrix<vector> faVectorMatrix;
+typedef faMatrix<sphericalTensor> faSphericalTensorMatrix;
+typedef faMatrix<symmTensor> faSymmTensorMatrix;
+typedef faMatrix<tensor> faTensorMatrix;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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);

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,25 +109,47 @@ public:
 
     // Member Functions
 
-        // Access
-
-            //- Calculate the geometry for the patches
-            //  (transformation tensors etc.)
-            void calcGeometry();
-
         //- Return the mesh reference
-            const faMesh& mesh() const;
+        const faMesh& mesh() const noexcept
+        {
+            return mesh_;
+        }
 
         //- Return a list of pointers for each patch
         //- with only those pointing to interfaces being set
         lduInterfacePtrsList interfaces() const;
 
+        //- The number of patches before the first processor patch.
+        label nNonProcessor() const;
+
         //- Return a list of patch names
         wordList names() const;
 
         //- Return a list of patch types
         wordList types() const;
 
+        //- Return a list of patch start indices
+        labelList patchStarts() const;
+
+        //- Return a list of patch sizes (number of edges in each patch)
+        labelList patchSizes() const;
+
+        //- Return a list of patch ranges
+        List<labelRange> patchRanges() const;
+
+        //- The start label of the edges in the faMesh edges list
+        //  Same as mesh.nInternalEdges()
+        label start() const;
+
+        //- The number of boundary edges for the underlying mesh
+        //  Same as mesh.nBoundaryEdges()
+        label nEdges() const;
+
+        //- The edge range for all boundary edges
+        //  Spans [nInternalEdges, nEdges) of the underlying mesh
+        labelRange range() const;
+
+
         //- Return patch indices for all matches.
         //  A no-op (returns -1) for an empty key
         //  \note Matching patchGroups currently not supported
@@ -156,17 +176,28 @@ public:
 
     // 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();
 
-            //- writeData member function required by regIOobject
-            bool writeData(Ostream&) const;
+        //- 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
+    // Ostream Operator
 
         friend Ostream& operator<<(Ostream&, const faBoundaryMesh&);
 

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
-        (
-            "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)
-{
-    DebugInFunction << "Creating faMesh from definition file" << endl;
-
-    // Reading faMeshDefinition dictionary
-    IOdictionary faMeshDefinition
-    (
-        IOobject
-        (
-            defFile,
-            mesh().time().constant(),
-            meshSubDir,
-            mesh(),
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    );
-
-    const wordList polyMeshPatches
+        pMesh,
+        selectPatchFaces
         (
+            pMesh.boundaryMesh(),
             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;
-            }
-        }
-    }
+{
+    DebugInFunction << "Creating from definition (dictionary)" << endl;
 
-
-    // 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
+    PtrList<faPatch> newPatches
     (
-            "ngbPolyPatchIndex",
-            faPatches[pI].ngbPolyPatchID_
+        createPatchList
+        (
+            faMeshDefinition.subDict("boundary"),
+
+            // Optional 'empty' patch
+            faMeshDefinition.getOrDefault<word>("emptyPatch", word::null),
+
+            // Optional specification for default patch
+            faMeshDefinition.findDict("defaultPatch")
+        )
     );
 
-        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_;

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,11 @@ Author
 namespace Foam
 {
 
-// Class forward declarations
+// Forward Declarations
 class faMeshLduAddressing;
 class faMeshMapper;
+class faPatchData;
+template<class T> class LabelledItem;
 
 /*---------------------------------------------------------------------------*\
                            Class faMesh Declaration
@@ -77,7 +79,6 @@ class faMeshMapper;
 
 class faMesh
 :
-    public GeoMesh<polyMesh>,
     public MeshObject<polyMesh, Foam::UpdateableMeshObject, faMesh>,
     public lduMesh,
     public edgeInterpolation,
@@ -85,7 +86,7 @@ class faMesh
     public faSolution,
     public data
 {
-    // Private data
+    // Private Data
 
         //- Face labels
         labelIOList faceLabels_;
@@ -130,7 +131,7 @@ class faMesh
     // Demand-driven data
 
         //- Primitive patch
-        mutable indirectPrimitivePatch* patchPtr_;
+        mutable uindirectPrimitivePatch* patchPtr_;
 
         //- Ldu addressing data
         mutable faMeshLduAddressing* lduPtr_;
@@ -187,7 +188,7 @@ class faMesh
         // Other mesh-related data
 
         //- Parallel info
-            mutable faGlobalMeshData* globalMeshDataPtr_;
+        mutable autoPtr<faGlobalMeshData> globalMeshDataPtr_;
 
 
     // Static Private Data
@@ -204,6 +205,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 +265,39 @@ class faMesh
             //- Clear demand-driven data
             void clearOut() const;
 
+
+    // Helpers
+
+        //- Get the polyPatch pairs for the boundary edges (natural order)
+        List<LabelledItem<edge>> getBoundaryEdgePatchPairs() const;
+
+        //- 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 List<LabelledItem<edge>>& bndEdgePatchPairs
+        ) const;
+
+
 public:
 
-    // Public typedefs
+    // Public Typedefs
 
         typedef faMesh Mesh;
         typedef faBoundaryMesh BoundaryMesh;
@@ -273,36 +306,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")
+        //- 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
         );
 
 
@@ -315,22 +350,27 @@ public:
     // Helpers
 
         //- Add boundary patches. Constructor helper
-            void addFaPatches(const List<faPatch*> &);
+        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 +387,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 +449,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 +480,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 +514,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 +573,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 +592,8 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#include "faMeshI.H"
+
 #ifdef NoRepository
     #include "faPatchFaMeshTemplates.C"
 #endif

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 "

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_;
+}
+
+
+// ************************************************************************* //

src/finiteArea/faMesh/faMeshPatches.C

@@ -0,0 +1,954 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "globalMeshData.H"
+#include "indirectPrimitivePatch.H"
+#include "edgeHashes.H"
+#include "LabelledItem.H"
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Manage patch pairs with a 'labelled' edge.
+// The edge first/second correspond to the owner/neighbour patches.
+// The index is a face index on the neighbour patch (FUTURE).
+
+// Local typedefs
+
+typedef LabelledItem<edge> patchPairInfo;
+typedef List<patchPairInfo> patchPairInfoList;
+typedef UIndirectList<patchPairInfo> patchPairInfoUIndList;
+
+
+// Handling of dangling coupled edges.
+// Tag values to "push" with special -(patchId+2)
+struct combineDanglingEdge
+{
+    const label upperLimit;
+
+    // Set dangling patchId from real patchId
+    static void setDangling(patchPairInfo& pairing, const label patchId)
+    {
+        pairing.first() = pairing.second() = -(patchId + 2);
+        pairing.setIndex(-1);  // Invalidate
+    }
+
+    // Convert dangling patchId to real patchId
+    static void correct(patchPairInfo& pairing)
+    {
+        if (pairing.first() < -1)
+        {
+            pairing.first() = -(pairing.first() + 2);
+        }
+        if (pairing.second() < -1)
+        {
+            pairing.second() = -(pairing.second() + 2);
+        }
+    }
+
+    //- Construct with upper limit (the number of non-processor patches)
+    explicit combineDanglingEdge(const label nNonProcessor)
+    :
+        upperLimit(nNonProcessor)
+    {}
+
+
+    // Combine operation: overwrite unused or processor patches with
+    // 'dangling' patch information only
+    void operator()(patchPairInfo& x, const patchPairInfo& y) const
+    {
+        if (y.first() < -1 && edge::compare(x, y) == 0)
+        {
+            if (x.first() == -1 || x.first() >= upperLimit)
+            {
+                x.first() = y.first();
+            }
+            if (x.second() == -1 || x.second() >= upperLimit)
+            {
+                x.second() = y.first();
+                x.index() = y.index();
+            }
+        }
+    }
+};
+
+
+// Populate patch pairings according to the boundary edges
+void findEdgePatchPairing
+(
+    const polyBoundaryMesh& pbm,
+    const EdgeMap<label>& edgeToIndex,
+    patchPairInfoList& patchPairs,
+    label nMissing = -1
+)
+{
+    // Count how many slots (both sides) to be filled
+    if (nMissing < 0)
+    {
+        nMissing = 0;
+        for (const patchPairInfo& pairing : patchPairs)
+        {
+            if (pairing.first() == -1) ++nMissing;
+            if (pairing.second() == -1) ++nMissing;
+        }
+    }
+
+    forAll(pbm, patchID)
+    {
+        if (!nMissing) break;  // Everything filled
+
+        const polyPatch& pp = pbm[patchID];
+
+        const bool isProcPatch = isA<processorPolyPatch>(pp);
+
+        // Examine neighbour boundary edges
+        for (label edgei = pp.nInternalEdges(); edgei < pp.nEdges(); ++edgei)
+        {
+            // Lookup global edge of this neighbour,
+            // find matching owner boundary edge
+
+            const label edgeIndex =
+                edgeToIndex.lookup(pp.meshEdge(edgei), -1);
+
+            if (edgeIndex != -1)
+            {
+                // Add patchId.
+                // - hash-like so will only insert once
+                // - also add in the attached patch face (there is only one),
+                //   saved in mesh face numbering for processor patches
+
+                patchPairInfo& pairing = patchPairs[edgeIndex];
+
+                if (pairing.insert(patchID))
+                {
+                    if (isProcPatch)
+                    {
+                        // Save the mesh face
+
+                        const label patchFacei = pp.edgeFaces()[edgei][0];
+
+                        pairing.index() = (pp.start() + patchFacei);
+                    }
+
+                    --nMissing;
+
+                    if (!nMissing) break;  // Early exit
+                }
+            }
+        }
+    }
+}
+
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::faMesh::reorderProcEdges
+(
+    faPatchData& patchDef,
+    const List<LabelledItem<edge>>& bndEdgePatchPairs
+) const
+{
+    if (!patchDef.coupled() || patchDef.edgeLabels_.empty())
+    {
+        return;
+    }
+
+    const polyBoundaryMesh& pbm = mesh().boundaryMesh();
+
+    const label procPatchID = patchDef.neighPolyPatchId_;
+
+    const auto* procPatch = isA<processorPolyPatch>(pbm[procPatchID]);
+
+    if (!procPatch)
+    {
+        FatalErrorInFunction
+            << "Internal addressing error. Patch " << procPatchID
+            << " is not a processor patch" << nl
+            << abort(FatalError);
+    }
+
+    // Reorder processor edges using order of neighbour processorPolyPatch
+
+    const label nProcEdges = patchDef.edgeLabels_.size();
+
+    labelList procFaces(nProcEdges, -1);
+
+    forAll(procFaces, edgei)
+    {
+        const label bndEdgei =
+            (patchDef.edgeLabels_[edgei] - patch().nInternalEdges());
+
+        procFaces[edgei] = bndEdgePatchPairs[bndEdgei].index();
+    }
+
+    // Ascending proc-face numbering
+    const labelList sortIndices(Foam::sortedOrder(procFaces));
+
+    if (procFaces.size() && procFaces[sortIndices[0]] < 0)
+    {
+        FatalErrorInFunction
+            << "Internal addressing error. Patch " << procPatchID
+            << " with negative face" << nl
+            << abort(FatalError);
+    }
+
+
+    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 meshFacei = procFaces[sortIndices[edgei]];
+
+        // Find all identical faces
+        label endEdgei = edgei + 1;  // one beyond
+        while
+        (
+            (endEdgei < nProcEdges)
+         && (meshFacei == 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 patch local edgeId
+            for (label i = edgei; i < endEdgei; ++i)
+            {
+                const label patchEdgei = oldEdgeLabels[sortIndices[i]];
+
+                // The edge in mesh numbering
+                multihit.insert
+                (
+                    patch().meshEdge(patchEdgei),
+                    patchEdgei
+                );
+            }
+
+            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()[meshFacei];
+
+            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: " << meshFacei
+                    << " ... indicates serious geometry issue" << nl
+                    << multihit << nl
+                    << abort(FatalError);
+            }
+            if (!multihit.empty())
+            {
+                FatalErrorInFunction
+                    << "Missed edges for face: " << meshFacei
+                    << " ... 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,
+        word::null,     // Name for empty patch placeholder
+        &onePatchDict   // Definitions for defaultPatch
+    );
+}
+
+
+Foam::List<Foam::LabelledItem<Foam::edge>>
+Foam::faMesh::getBoundaryEdgePatchPairs() const
+{
+    const polyBoundaryMesh& pbm = mesh().boundaryMesh();
+
+    const label nInternalEdges = patch().nInternalEdges();
+    const label nBoundaryEdges = patch().nBoundaryEdges();
+
+    // Map edges (mesh numbering) back to a boundary index
+    EdgeMap<label> edgeToBoundaryIndex(2*nBoundaryEdges);
+
+    // Use labelled 'edge' for accounting of patch pairs
+    patchPairInfoList bndEdgePatchPairs(nBoundaryEdges);
+
+
+    // Pass 1:
+    // - setup lookup (edge -> bnd index)
+    // - add owner patch for each boundary edge
+    {
+        const SubList<labelList> bndEdgeToPatchFace
+        (
+            patch().edgeFaces(),
+            patch().nBoundaryEdges(),
+            patch().nInternalEdges()
+        );
+
+        for (label bndEdgei = 0; bndEdgei < nBoundaryEdges; ++bndEdgei)
+        {
+            edgeToBoundaryIndex.insert
+            (
+                patch().meshEdge(bndEdgei + nInternalEdges),
+                edgeToBoundaryIndex.size()
+            );
+
+            // The attached patch face (there is only one):
+            const label patchFacei = bndEdgeToPatchFace[bndEdgei][0];
+
+            const label meshFacei = faceLabels_[patchFacei];
+
+            const label patchId = pbm.whichPatch(meshFacei);
+
+            bndEdgePatchPairs[bndEdgei].insert(patchId);
+        }
+    }
+
+    // Pass 2:
+    // - Add in first neighbour patch for the boundary edges
+    // - examine all possible connecting neighbours
+
+    findEdgePatchPairing
+    (
+        pbm,
+        edgeToBoundaryIndex,
+        bndEdgePatchPairs,
+        edgeToBoundaryIndex.size()  // Number of places to still fill
+    );
+
+
+    // Nothing dangling if running in serial - can return already
+    if (!Pstream::parRun())
+    {
+        return bndEdgePatchPairs;
+    }
+
+    // In parallel need to check for "dangling" edges, which are finiteArea
+    // boundary edges that only exist on one side of a proc boundary.
+    // Eg, proc boundary coincides with the outer extent of the finiteArea
+
+    const globalMeshData& globalData = mesh().globalData();
+    const indirectPrimitivePatch& cpp = globalData.coupledPatch();
+    const mapDistribute& map = globalData.globalEdgeSlavesMap();
+
+    // Construct coupled edge usage with all data
+    List<unsigned char> coupledEdgesUsed(map.constructSize(), 0u);
+
+    forAll(cpp.edges(), coupledEdgei)
+    {
+        const auto iter =
+            edgeToBoundaryIndex.cfind(cpp.meshEdge(coupledEdgei));
+
+        // Used from this side or other other side
+        coupledEdgesUsed[coupledEdgei] = (iter.found() ? 1 : 2);
+    }
+
+    // Save the original (pre-sync) coupling state
+    const List<unsigned char> coupledEdgesOrig(coupledEdgesUsed);
+
+    globalData.syncData
+    (
+        coupledEdgesUsed,
+        globalData.globalEdgeSlaves(),
+        globalData.globalEdgeTransformedSlaves(),  // probably not used
+        map,
+        bitOrEqOp<unsigned char>()
+    );
+
+
+    // Check for one-sided edge coupling (coupled value == 3)
+    // original == 1:
+    // - coupled to a real finiteArea edge.
+    // - receive a patch-pair value
+    //
+    // original == 2:
+    // - a "dangled" edge. Information required for other procs.
+    // - push a patch-pair value
+
+    // Map edges (mesh numbering) back to a coupled index.
+    // These are the edges to 'push' information for.
+
+    EdgeMap<label> edgeToCoupledIndex;
+
+    label nEdgesPull = 0;
+
+    forAll(coupledEdgesUsed, coupledEdgei)
+    {
+        if (coupledEdgesUsed[coupledEdgei] == 3)
+        {
+            if (coupledEdgesOrig[coupledEdgei] == 1)
+            {
+                // Coupled side with finiteArea
+                ++nEdgesPull;
+            }
+            else if (coupledEdgesOrig[coupledEdgei] == 2)
+            {
+                // Coupled side without finiteArea
+                edgeToCoupledIndex.insert
+                (
+                    cpp.meshEdge(coupledEdgei),
+                    coupledEdgei
+                );
+            }
+        }
+    }
+
+    // Nothing to do - can return already
+    if (returnReduce(edgeToCoupledIndex.empty(), andOp<bool>()))
+    {
+        return bndEdgePatchPairs;
+    }
+
+    // Data locations to pull
+    labelList patchEdgeLabels(nEdgesPull);
+    labelList coupledEdgeLabels(nEdgesPull);
+
+    // Populate the locations
+    {
+        nEdgesPull = 0;
+
+        forAll(cpp.edges(), coupledEdgei)
+        {
+            if
+            (
+                coupledEdgesUsed[coupledEdgei] == 3
+             && coupledEdgesOrig[coupledEdgei] == 1
+            )
+            {
+                // Pull this edge
+                const auto iter =
+                    edgeToBoundaryIndex.cfind(cpp.meshEdge(coupledEdgei));
+
+                if (iter.found())
+                {
+                    patchEdgeLabels[nEdgesPull] = iter.val();
+                    coupledEdgeLabels[nEdgesPull] = coupledEdgei;
+                    ++nEdgesPull;
+                }
+                else
+                {
+                    // Should be impossible to fail here
+                    FatalErrorInFunction
+                        << "Failed on second lookup of "
+                        << cpp.meshEdge(coupledEdgei) << nl
+                        << abort(FatalError);
+                }
+            }
+        }
+
+        if (nEdgesPull != coupledEdgeLabels.size())
+        {
+            FatalErrorInFunction
+                << "Failed lookup of some coupled edges" << nl
+                << abort(FatalError);
+        }
+    }
+
+    //- Construct edge sync with all data
+    patchPairInfoList cppEdgeData(map.constructSize());
+
+    // Fill in for 'push' locations. Only really interested in the owner
+    // (corresponds to non-proc connection), but grab everything
+    findEdgePatchPairing
+    (
+        pbm,
+        edgeToCoupledIndex,
+        cppEdgeData,
+        2*edgeToCoupledIndex.size()  // Accept both sides?
+    );
+
+
+    const label nNonProcessor = pbm.nNonProcessor();
+
+    // Adjust patch information to reflect dangling patch neighbour
+    // Tag with -(value+2)
+    forAllConstIters(edgeToCoupledIndex, iter)
+    {
+        const edge& e = iter.key();
+        const label coupledEdgei = iter.val();
+
+        patchPairInfo& pairing = cppEdgeData[coupledEdgei];
+        const label ownerPatchId = pairing.first();
+
+        // Some sanity checks
+        if (ownerPatchId < 0)
+        {
+            FatalErrorInFunction
+                << "Error finding dangling edge at "
+                << cpp.points()[e.first()] << ' '
+                << cpp.points()[e.second()] << nl
+                << abort(FatalError);
+        }
+        else if (ownerPatchId >= nNonProcessor)
+        {
+            FatalErrorInFunction
+                << "Cannot handle edge on processor-processor connection at "
+                << cpp.points()[e.first()] << ' '
+                << cpp.points()[e.second()] << nl
+                << abort(FatalError);
+        }
+
+        combineDanglingEdge::setDangling(pairing, ownerPatchId);
+
+        // TBD:
+        // may wish to remember the corresponding proc number,
+        // if we wish to bridge across 'fan-like' connections.
+        //
+        // pairing.setIndex(-(Pstream::myProcNo() + 2));
+    }
+
+
+    // Synchronize edge information
+
+    const combineDanglingEdge edgeCombineOp(nNonProcessor);
+
+    globalMeshData::syncData
+    (
+        cppEdgeData,
+        globalData.globalEdgeSlaves(),
+        globalData.globalEdgeTransformedSlaves(),  // probably not used
+        map,
+        edgeCombineOp
+    );
+
+
+    // Combine back from pushed cpp-edge data
+
+    forAll(patchEdgeLabels, i)
+    {
+        patchPairInfo& pairing = bndEdgePatchPairs[patchEdgeLabels[i]];
+        const patchPairInfo& other = cppEdgeData[coupledEdgeLabels[i]];
+
+        edgeCombineOp(pairing, other);
+
+        // Resolve special tagging
+        combineDanglingEdge::correct(pairing);
+    }
+
+    return bndEdgePatchPairs;
+}
+
+
+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_);
+        }
+    }
+
+    // ----------------------------------------------------------------------
+
+    const label nInternalEdges = patch().nInternalEdges();
+    const label nBoundaryEdges = patch().nBoundaryEdges();
+
+    patchPairInfoList bndEdgePatchPairs(this->getBoundaryEdgePatchPairs());
+
+    labelList bndEdgeFaPatchIDs(nBoundaryEdges, -1);
+
+    for (label bndEdgei = 0; bndEdgei < nBoundaryEdges; ++bndEdgei)
+    {
+        const patchPairInfo& patchPair = bndEdgePatchPairs[bndEdgei];
+
+        if (patchPair.valid())
+        {
+            // Non-negative, unique pairing
+            // - find corresponding definition
+
+            for (label patchi = 0; patchi < faPatchDefs.size(); ++patchi)
+            {
+                if (faPatchDefs[patchi].foundPatchPair(patchPair))
+                {
+                    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 patchEdgei = undefinedEdges[edgei];
+            const label bndEdgei = (patchEdgei - nInternalEdges);
+
+            edgeNbrPolyPatch[edgei] = bndEdgePatchPairs[bndEdgei].second();
+        }
+
+        // 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, bndEdgePatchPairs);
+            patchDef.neighPolyPatchId_ = -1; // No lookup of neighbour faces
+        }
+    }
+
+    // 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;
+}
+
+
+// ************************************************************************* //

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

src/finiteArea/faMesh/faPatches/basic/coupled/coupledFaPatch.C

@@ -117,10 +117,4 @@ void Foam::coupledFaPatch::calcTransformTensors
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::coupledFaPatch::~coupledFaPatch()
-{}
-
-
 // ************************************************************************* //

src/finiteArea/faMesh/faPatches/basic/coupled/coupledFaPatch.H

@@ -135,7 +135,7 @@ public:
 
 
     //- Destructor
-    virtual ~coupledFaPatch();
+    virtual ~coupledFaPatch() = default;
 
 
     // Member Functions

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
-
 // ************************************************************************* //

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,6 +66,15 @@ public:
 
     // Constructors
 
+        //- Minimal construct from components
+        emptyFaPatch
+        (
+            const word& name,
+            const label index,
+            const faBoundaryMesh& bm,
+            const label ngbPolyPatchIndex = -1
+        );
+
         //- Construct from components
         emptyFaPatch
         (
@@ -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;
-
 };
 
 

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

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 "

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>
+    const auto* wedgePtr = isA<wedgePolyPatch>
     (
         bm.mesh()().boundaryMesh()[ngbPolyPatchIndex()]
     );
 
-        wedgePolyPatchPtr_ = &wedge;
+    if (wedgePtr)
+    {
+        wedgePolyPatchPtr_ = wedgePtr;
     }
     else
     {

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);
 }
 
 

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

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()
+        )
+    );
+}
+
+
+// ************************************************************************* //

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
 {
+public:
+
+    // Data Members
+
         word name_;
         word type_;
-    dictionary dict_;
-    label ownPolyPatchID_;
-    label ngbPolyPatchID_;
+
+        label ownerPolyPatchId_;
+        label neighPolyPatchId_;
+
+        //- The owner/neighbour for processor patches
+        int ownerProcId_;
+        int neighProcId_;
+
+        // Storge (temporary or otherwise) for edge labels
         labelList edgeLabels_;
-    faPatchData()
-    :
-        name_(word::null),
-        type_(word::null),
-        ownPolyPatchID_(-1),
-        ngbPolyPatchID_(-1)
-    {}
+
+
+    // 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;
 };
 
 

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)
     );
 }
 

src/finiteVolume/cfdTools/general/fvOptions/fvOption.C

@@ -58,9 +58,9 @@ Foam::fv::option::option
     mesh_(mesh),
     dict_(dict),
     coeffs_(dict.optionalSubDict(modelType + "Coeffs")),
-    active_(dict_.getOrDefault<Switch>("active", true)),
     fieldNames_(),
     applied_(),
+    active_(dict_.getOrDefault("active", true)),
     log(true)
 {
     Log << incrIndent << indent << "Source: " << name_ << endl << decrIndent;
@@ -101,7 +101,7 @@ Foam::autoPtr<Foam::fv::option> Foam::fv::option::New
         ) << exit(FatalIOError);
     }
 
-    return autoPtr<option>(cstrIter()(name, modelType, coeffs, mesh));
+    return autoPtr<fv::option>(cstrIter()(name, modelType, coeffs, mesh));
 }
 
 

src/finiteVolume/cfdTools/general/fvOptions/fvOption.H

@@ -113,23 +113,25 @@ protected:
         //- Dictionary containing source coefficients
         dictionary coeffs_;
 
-        //- Source active flag
-        Switch active_;
-
         //- Field names to apply source to - populated by derived models
         wordList fieldNames_;
 
         //- Applied flag list - corresponds to each fieldNames_ entry
         List<bool> applied_;
 
+        //- Source active flag
+        bool active_;
+
 
 public:
 
+    //- Switch write log to Info
+    bool log;
+
+
     //- Runtime type information
     TypeName("option");
 
-    //- Switch write log to Info
-    bool log;
 
     // Declare run-time constructor selection table
 
@@ -219,16 +221,16 @@ public:
         // Access
 
             //- Return const access to the source name
-            inline const word& name() const;
+            inline const word& name() const noexcept;
 
             //- Return const access to the mesh database
-            inline const fvMesh& mesh() const;
+            inline const fvMesh& mesh() const noexcept;
 
             //- Return dictionary
-            inline const dictionary& coeffs() const;
+            inline const dictionary& coeffs() const noexcept;
 
-            //- Return const access to the source active flag
-            inline bool active() const;
+            //- True if source is active
+            inline bool active() const noexcept;
 
             //- Set the applied flag to true for field index fieldi
             inline void setApplied(const label fieldi);
@@ -236,8 +238,8 @@ public:
 
         // Edit
 
-            //- Return access to the source active flag
-            inline Switch& active();
+            //- Change source active flag, return previous value
+            inline bool active(const bool on) noexcept;
 
 
         // Checks

src/finiteVolume/cfdTools/general/fvOptions/fvOptionI.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,40 +28,42 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline const Foam::word& Foam::fv::option::name() const
+inline const Foam::word& Foam::fv::option::name() const noexcept
 {
     return name_;
 }
 
 
-inline const Foam::fvMesh& Foam::fv::option::mesh() const
+inline const Foam::fvMesh& Foam::fv::option::mesh() const noexcept
 {
     return mesh_;
 }
 
 
-inline const Foam::dictionary& Foam::fv::option::coeffs() const
+inline const Foam::dictionary& Foam::fv::option::coeffs() const noexcept
 {
     return coeffs_;
 }
 
 
-inline bool Foam::fv::option::active() const
+inline bool Foam::fv::option::active() const noexcept
 {
     return active_;
 }
 
 
+inline bool Foam::fv::option::active(const bool on) noexcept
+{
+    bool old(active_);
+    active_ = on;
+    return old;
+}
+
+
 inline void Foam::fv::option::setApplied(const label fieldi)
 {
     applied_[fieldi] = true;
 }
 
 
-inline Foam::Switch& Foam::fv::option::active()
-{
-    return active_;
-}
-
-
 // ************************************************************************* //

src/finiteVolume/cfdTools/general/fvOptions/fvOptionList.C

@@ -89,7 +89,7 @@ void Foam::fv::optionList::checkApplied() const
 
 Foam::fv::optionList::optionList(const fvMesh& mesh, const dictionary& dict)
 :
-    PtrList<option>(),
+    PtrList<fv::option>(),
     mesh_(mesh),
     checkTimeIndex_(mesh_.time().startTimeIndex() + 2)
 {
@@ -99,7 +99,7 @@ Foam::fv::optionList::optionList(const fvMesh& mesh, const dictionary& dict)
 
 Foam::fv::optionList::optionList(const fvMesh& mesh)
 :
-    PtrList<option>(),
+    PtrList<fv::option>(),
     mesh_(mesh),
     checkTimeIndex_(mesh_.time().startTimeIndex() + 2)
 {}

src/finiteVolume/cfdTools/general/fvOptions/fvOptionList.H

@@ -69,7 +69,7 @@ namespace fv
 
 class optionList
 :
-    public PtrList<option>
+    public PtrList<fv::option>
 {
 protected:
 

src/fvOptions/cellSetOption/cellSetOption.H

@@ -117,7 +117,7 @@ namespace fv
 
 class cellSetOption
 :
-    public option
+    public fv::option
 {
 public:
 

src/fvOptions/interRegionOption/interRegionOption.H

@@ -87,7 +87,7 @@ namespace fv
 
 class interRegionOption
 :
-    public option
+    public fv::option
 {
 protected:
 

src/fvOptions/sources/derived/buoyancyEnergy/buoyancyEnergy.H

@@ -118,7 +118,7 @@ namespace fv
 
 class buoyancyEnergy
 :
-    public option
+    public fv::option
 {
     // Private Data
 

src/fvOptions/sources/derived/buoyancyForce/buoyancyForce.H

@@ -110,7 +110,7 @@ namespace fv
 
 class buoyancyForce
 :
-    public option
+    public fv::option
 {
     // Private Data
 

src/fvOptions/sources/derived/jouleHeatingSource/jouleHeatingSource.H

@@ -175,7 +175,7 @@ namespace fv
 
 class jouleHeatingSource
 :
-    public option
+    public fv::option
 {
     // Private Data
 

src/fvOptions/sources/derived/multiphaseStabilizedTurbulence/multiphaseStabilizedTurbulence.H

@@ -133,7 +133,7 @@ namespace fv
 
 class multiphaseStabilizedTurbulence
 :
-    public option
+    public fv::option
 {
     // Private Data
 

src/fvOptions/sources/derived/phaseLimitStabilization/PhaseLimitStabilization.H

@@ -108,7 +108,7 @@ namespace fv
 template<class Type>
 class PhaseLimitStabilization
 :
-    public option
+    public fv::option
 {
     // Private Data
 

src/fvOptions/sources/derived/tabulatedAccelerationSource/tabulatedAccelerationSource.H

@@ -87,7 +87,7 @@ namespace fv
 
 class tabulatedAccelerationSource
 :
-    public option
+    public fv::option
 {
 protected:
 

src/fvOptions/sources/derived/viscousDissipation/viscousDissipation.H

@@ -111,7 +111,7 @@ namespace fv
 
 class viscousDissipation
 :
-    public option
+    public fv::option
 {
     // Private Data
 

src/parallel/decompose/Allwclean

@@ -8,6 +8,7 @@ wclean kahipDecomp
 wclean scotchDecomp
 wclean decompositionMethods
 wclean decompose
+wclean faDecompose
 
 ./Allwclean-mpi
 

src/parallel/decompose/Allwmake

@@ -12,6 +12,7 @@ export FOAM_EXT_LIBBIN
 
 wmake $targetType decompositionMethods
 wmake $targetType decompose
+wmake $targetType faDecompose
 
 if have_kahip
 then

src/parallel/decompose/decompose/Make/files

@@ -1,5 +1,8 @@
 decompositionInformation.C
 decompositionModel.C
+
+dimFieldDecomposer.C
 fvFieldDecomposer.C
+pointFieldDecomposer.C
 
 LIB = $(FOAM_LIBBIN)/libdecompose

src/parallel/decompose/decompose/dimFieldDecomposer.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,27 +28,30 @@ License
 
 #include "dimFieldDecomposer.H"
 
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::dimFieldDecomposer::dimFieldDecomposer
 (
-    const fvMesh& completeMesh,
     const fvMesh& procMesh,
-    const labelList& faceAddressing,
     const labelList& cellAddressing
 )
 :
-    completeMesh_(completeMesh),
     procMesh_(procMesh),
-    faceAddressing_(faceAddressing),
     cellAddressing_(cellAddressing)
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::dimFieldDecomposer::~dimFieldDecomposer()
+Foam::dimFieldDecomposer::dimFieldDecomposer
+(
+    const fvMesh& /* unused: completeMesh */,
+    const fvMesh& procMesh,
+    const labelList& /* unused: faceAddressing */,
+    const labelList& cellAddressing
+)
+:
+    procMesh_(procMesh),
+    //UNUSED: faceAddressing_(faceAddressing),
+    cellAddressing_(cellAddressing)
 {}
 
 

src/parallel/decompose/decompose/dimFieldDecomposer.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,7 +32,7 @@ Description
 
 SourceFiles
     dimFieldDecomposer.C
-    dimFieldDecomposerDecomposeFields.C
+    dimFieldDecomposerFields.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -46,32 +47,25 @@ SourceFiles
 namespace Foam
 {
 
-class IOobjectList;
-
 /*---------------------------------------------------------------------------*\
                     Class fvFieldDecomposer Declaration
 \*---------------------------------------------------------------------------*/
 
 class dimFieldDecomposer
 {
-private:
-
-    // Private data
-
-        //- Reference to complete mesh
-        const fvMesh& completeMesh_;
+    // Private Data
 
         //- Reference to processor mesh
         const fvMesh& procMesh_;
 
         //- Reference to face addressing
-        const labelList& faceAddressing_;
+        //UNUSED: const labelList& faceAddressing_;
 
         //- Reference to cell addressing
         const labelList& cellAddressing_;
 
 
-    // Private Member Functions
+public:
 
     //- No copy construct
     dimFieldDecomposer(const dimFieldDecomposer&) = delete;
@@ -80,35 +74,41 @@ private:
     void operator=(const dimFieldDecomposer&) = delete;
 
 
-public:
-
     // Constructors
 
-        //- Construct from components
+        //- Construct from minimal components
         dimFieldDecomposer
         (
-            const fvMesh& completeMesh,
             const fvMesh& procMesh,
-            const labelList& faceAddressing,
+            const labelList& cellAddressing
+        );
+
+        //- Construct from components with API as per fvFieldDecomposer
+        dimFieldDecomposer
+        (
+            const fvMesh& completeMesh,         //!< unused
+            const fvMesh& procMesh,
+            const labelList& faceAddressing,    //!< unused
             const labelList& cellAddressing
         );
 
 
     //- Destructor
-    ~dimFieldDecomposer();
+    ~dimFieldDecomposer() = default;
 
 
     // Member Functions
 
         //- Decompose field
         template<class Type>
-        tmp<DimensionedField<Type, volMesh>> decomposeField
+        tmp<DimensionedField<Type, volMesh>>
+        decomposeField
         (
             const DimensionedField<Type, volMesh>& field
         ) const;
 
 
-        //- Decompose llist of fields
+        //- Decompose list of fields
         template<class GeoField>
         void decomposeFields(const PtrList<GeoField>& fields) const;
 };
@@ -121,7 +121,7 @@ public:
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-    #include "dimFieldDecomposerDecomposeFields.C"
+    #include "dimFieldDecomposerFields.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/parallel/decompose/decompose/dimFieldDecomposerFields.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -40,9 +41,8 @@ Foam::dimFieldDecomposer::decomposeField
     Field<Type> mappedField(field, cellAddressing_);
 
     // Create the field for the processor
-    return tmp<DimensionedField<Type, volMesh>>
-    (
-        new DimensionedField<Type, volMesh>
+    return
+        tmp<DimensionedField<Type, volMesh>>::New
         (
             IOobject
             (
@@ -55,8 +55,7 @@ Foam::dimFieldDecomposer::decomposeField
             ),
             procMesh_,
             field.dimensions(),
-            mappedField
-        )
+            std::move(mappedField)
         );
 }
 
@@ -67,9 +66,9 @@ void Foam::dimFieldDecomposer::decomposeFields
     const PtrList<GeoField>& fields
 ) const
 {
-    forAll(fields, fieldi)
+    for (const auto& fld : fields)
     {
-        decomposeField(fields[fieldi])().write();
+        decomposeField(fld)().write();
     }
 }
 

src/parallel/decompose/decompose/fvFieldDecomposer.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,7 +28,6 @@ License
 
 #include "fvFieldDecomposer.H"
 
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::fvFieldDecomposer::patchFieldDecomposer::patchFieldDecomposer
@@ -49,21 +49,19 @@ Foam::fvFieldDecomposer::patchFieldDecomposer::patchFieldDecomposer
 Foam::fvFieldDecomposer::processorVolPatchFieldDecomposer::
 processorVolPatchFieldDecomposer
 (
-    const fvMesh& mesh,
+    const labelUList& owner,  // == mesh.faceOwner()
+    const labelUList& neigh,  // == mesh.faceNeighbour()
     const labelUList& addressingSlice
 )
 :
     directAddressing_(addressingSlice.size())
 {
-    const labelList& own = mesh.faceOwner();
-    const labelList& neighb = mesh.faceNeighbour();
-
     forAll(directAddressing_, i)
     {
         // Subtract one to align addressing.
         label ai = mag(addressingSlice[i]) - 1;
 
-        if (ai < neighb.size())
+        if (ai < neigh.size())
         {
             // This is a regular face. it has been an internal face
             // of the original mesh and now it has become a face
@@ -73,11 +71,11 @@ processorVolPatchFieldDecomposer
             if (addressingSlice[i] >= 0)
             {
                 // I have the owner so use the neighbour value
-                directAddressing_[i] = neighb[ai];
+                directAddressing_[i] = neigh[ai];
             }
             else
             {
-                directAddressing_[i] = own[ai];
+                directAddressing_[i] = owner[ai];
             }
         }
         else
@@ -88,12 +86,28 @@ processorVolPatchFieldDecomposer
             // up the different (face) list of data), so I will
             // just grab the value from the owner cell
 
-            directAddressing_[i] = own[ai];
+            directAddressing_[i] = owner[ai];
         }
     }
 }
 
 
+Foam::fvFieldDecomposer::processorVolPatchFieldDecomposer::
+processorVolPatchFieldDecomposer
+(
+    const fvMesh& mesh,
+    const labelUList& addressingSlice
+)
+:
+    processorVolPatchFieldDecomposer
+    (
+        mesh.faceOwner(),
+        mesh.faceNeighbour(),
+        addressingSlice
+    )
+{}
+
+
 Foam::fvFieldDecomposer::processorSurfacePatchFieldDecomposer::
 processorSurfacePatchFieldDecomposer
 (
@@ -105,8 +119,8 @@ processorSurfacePatchFieldDecomposer
 {
     forAll(addressing_, i)
     {
-        addressing_[i].setSize(1);
-        weights_[i].setSize(1);
+        addressing_[i].resize(1);
+        weights_[i].resize(1);
 
         addressing_[i][0] = mag(addressingSlice[i]) - 1;
         weights_[i][0] = 1;
@@ -114,6 +128,27 @@ processorSurfacePatchFieldDecomposer
 }
 
 
+Foam::fvFieldDecomposer::fvFieldDecomposer
+(
+    const Foam::zero,
+    const fvMesh& procMesh,
+    const labelList& faceAddressing,
+    const labelList& cellAddressing,
+    const labelList& boundaryAddressing
+)
+:
+    procMesh_(procMesh),
+    faceAddressing_(faceAddressing),
+    cellAddressing_(cellAddressing),
+    boundaryAddressing_(boundaryAddressing),
+    // Mappers
+    patchFieldDecomposerPtrs_(),
+    processorVolPatchFieldDecomposerPtrs_(),
+    processorSurfacePatchFieldDecomposerPtrs_(),
+    faceSign_()
+{}
+
+
 Foam::fvFieldDecomposer::fvFieldDecomposer
 (
     const fvMesh& completeMesh,
@@ -123,24 +158,85 @@ Foam::fvFieldDecomposer::fvFieldDecomposer
     const labelList& boundaryAddressing
 )
 :
-    completeMesh_(completeMesh),
-    procMesh_(procMesh),
-    faceAddressing_(faceAddressing),
-    cellAddressing_(cellAddressing),
-    boundaryAddressing_(boundaryAddressing),
-    patchFieldDecomposerPtrs_(procMesh_.boundary().size()),
-    processorVolPatchFieldDecomposerPtrs_(procMesh_.boundary().size()),
-    processorSurfacePatchFieldDecomposerPtrs_(procMesh_.boundary().size()),
-    faceSign_(procMesh_.boundary().size())
+    fvFieldDecomposer
+    (
+        zero{},
+        procMesh,
+        faceAddressing,
+        cellAddressing,
+        boundaryAddressing
+    )
 {
+    reset(completeMesh);
+}
+
+
+Foam::fvFieldDecomposer::fvFieldDecomposer
+(
+    const List<labelRange>& boundaryRanges,
+    const labelUList& faceOwner,
+    const labelUList& faceNeighbour,
+
+    const fvMesh& procMesh,
+    const labelList& faceAddressing,
+    const labelList& cellAddressing,
+    const labelList& boundaryAddressing
+)
+:
+    fvFieldDecomposer
+    (
+        zero{},
+        procMesh,
+        faceAddressing,
+        cellAddressing,
+        boundaryAddressing
+    )
+{
+    reset(boundaryRanges, faceOwner, faceNeighbour);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::fvFieldDecomposer::empty() const
+{
+    return patchFieldDecomposerPtrs_.empty();
+}
+
+
+void Foam::fvFieldDecomposer::clear()
+{
+    patchFieldDecomposerPtrs_.clear();
+    processorVolPatchFieldDecomposerPtrs_.clear();
+    processorSurfacePatchFieldDecomposerPtrs_.clear();
+    faceSign_.clear();
+}
+
+
+void Foam::fvFieldDecomposer::reset
+(
+    const List<labelRange>& boundaryRanges,
+    const labelUList& faceOwner,
+    const labelUList& faceNeighbour
+)
+{
+    clear();
+    const label nMappers = procMesh_.boundary().size();
+    patchFieldDecomposerPtrs_.resize(nMappers);
+    processorVolPatchFieldDecomposerPtrs_.resize(nMappers);
+    processorSurfacePatchFieldDecomposerPtrs_.resize(nMappers);
+    faceSign_.resize(nMappers);
+
     forAll(boundaryAddressing_, patchi)
     {
+        const label oldPatchi = boundaryAddressing_[patchi];
         const fvPatch& fvp = procMesh_.boundary()[patchi];
+        const labelSubList localPatchSlice(fvp.patchSlice(faceAddressing_));
 
         if
         (
-            boundaryAddressing_[patchi] >= 0
-        && !isA<processorLduInterface>(procMesh.boundary()[patchi])
+            oldPatchi >= 0
+        && !isA<processorLduInterface>(procMesh_.boundary()[patchi])
         )
         {
             patchFieldDecomposerPtrs_.set
@@ -148,11 +244,8 @@ Foam::fvFieldDecomposer::fvFieldDecomposer
                 patchi,
                 new patchFieldDecomposer
                 (
-                    fvp.patchSlice(faceAddressing_),
-                    completeMesh_.boundaryMesh()
-                    [
-                        boundaryAddressing_[patchi]
-                    ].start()
+                    localPatchSlice,
+                    boundaryRanges[oldPatchi].start()
                 )
             );
         }
@@ -163,8 +256,9 @@ Foam::fvFieldDecomposer::fvFieldDecomposer
                 patchi,
                 new processorVolPatchFieldDecomposer
                 (
-                    completeMesh_,
-                    fvp.patchSlice(faceAddressing_)
+                    faceOwner,
+                    faceNeighbour,
+                    localPatchSlice
                 )
             );
 
@@ -173,28 +267,21 @@ Foam::fvFieldDecomposer::fvFieldDecomposer
                 patchi,
                 new processorSurfacePatchFieldDecomposer
                 (
-                    static_cast<const labelUList&>
-                    (
-                        fvp.patchSlice
-                        (
-                            faceAddressing_
-                        )
-                    )
+                    static_cast<const labelUList&>(localPatchSlice)
                 )
             );
 
             faceSign_.set
             (
                 patchi,
-                new scalarField(fvp.patchSlice(faceAddressing_).size())
+                new scalarField(localPatchSlice.size())
             );
 
             {
-                const SubList<label> fa = fvp.patchSlice(faceAddressing_);
                 scalarField& s = faceSign_[patchi];
                 forAll(s, i)
                 {
-                    s[i] = sign(fa[i]);
+                    s[i] = sign(localPatchSlice[i]);
                 }
             }
         }
@@ -202,10 +289,74 @@ Foam::fvFieldDecomposer::fvFieldDecomposer
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+void Foam::fvFieldDecomposer::reset(const fvMesh& completeMesh)
+{
+    clear();
+    const label nMappers = procMesh_.boundary().size();
+    patchFieldDecomposerPtrs_.resize(nMappers);
+    processorVolPatchFieldDecomposerPtrs_.resize(nMappers);
+    processorSurfacePatchFieldDecomposerPtrs_.resize(nMappers);
+    faceSign_.resize(nMappers);
 
-Foam::fvFieldDecomposer::~fvFieldDecomposer()
-{}
+    forAll(boundaryAddressing_, patchi)
+    {
+        const label oldPatchi = boundaryAddressing_[patchi];
+        const fvPatch& fvp = procMesh_.boundary()[patchi];
+        const labelSubList localPatchSlice(fvp.patchSlice(faceAddressing_));
+
+        if
+        (
+            oldPatchi >= 0
+        && !isA<processorLduInterface>(procMesh_.boundary()[patchi])
+        )
+        {
+            patchFieldDecomposerPtrs_.set
+            (
+                patchi,
+                new patchFieldDecomposer
+                (
+                    localPatchSlice,
+                    completeMesh.boundaryMesh()[oldPatchi].start()
+                )
+            );
+        }
+        else
+        {
+            processorVolPatchFieldDecomposerPtrs_.set
+            (
+                patchi,
+                new processorVolPatchFieldDecomposer
+                (
+                    completeMesh,
+                    localPatchSlice
+                )
+            );
+
+            processorSurfacePatchFieldDecomposerPtrs_.set
+            (
+                patchi,
+                new processorSurfacePatchFieldDecomposer
+                (
+                    static_cast<const labelUList&>(localPatchSlice)
+                )
+            );
+
+            faceSign_.set
+            (
+                patchi,
+                new scalarField(localPatchSlice.size())
+            );
+
+            {
+                scalarField& s = faceSign_[patchi];
+                forAll(s, i)
+                {
+                    s[i] = sign(localPatchSlice[i]);
+                }
+            }
+        }
+    }
+}
 
 
 // ************************************************************************* //

src/parallel/decompose/decompose/fvFieldDecomposer.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,7 +32,7 @@ Description
 
 SourceFiles
     fvFieldDecomposer.C
-    fvFieldDecomposerDecomposeFields.C
+    fvFieldDecomposerFields.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -47,8 +48,6 @@ SourceFiles
 namespace Foam
 {
 
-class IOobjectList;
-
 /*---------------------------------------------------------------------------*\
                     Class fvFieldDecomposer Declaration
 \*---------------------------------------------------------------------------*/
@@ -62,7 +61,7 @@ public:
         :
             public fvPatchFieldMapper
         {
-            // Private data
+            // Private Data
 
                 labelList directAddressing_;
 
@@ -78,7 +77,7 @@ public:
                 );
 
 
-            // Member functions
+            // Member Functions
 
                 label size() const
                 {
@@ -110,13 +109,21 @@ public:
         :
             public fvPatchFieldMapper
         {
-            // Private data
+            // Private Data
 
                 labelList directAddressing_;
 
         public:
 
-            //- Construct given addressing
+            //- Construct addressing from details
+            processorVolPatchFieldDecomposer
+            (
+                const labelUList& faceOwner,
+                const labelUList& faceNeigbour,
+                const labelUList& addressingSlice
+            );
+
+            //- Construct given addressing from complete mesh
             processorVolPatchFieldDecomposer
             (
                 const fvMesh& mesh,
@@ -124,7 +131,7 @@ public:
             );
 
 
-            // Member functions
+            // Member Functions
 
                 label size() const
                 {
@@ -199,10 +206,7 @@ public:
 
 private:
 
-    // Private data
-
-        //- Reference to complete mesh
-        const fvMesh& completeMesh_;
+    // Private Data
 
         //- Reference to processor mesh
         const fvMesh& procMesh_;
@@ -225,11 +229,10 @@ private:
         PtrList<processorSurfacePatchFieldDecomposer>
             processorSurfacePatchFieldDecomposerPtrs_;
 
-
         PtrList<scalarField> faceSign_;
 
 
-    // Private Member Functions
+public:
 
     //- No copy construct
     fvFieldDecomposer(const fvFieldDecomposer&) = delete;
@@ -238,11 +241,19 @@ private:
     void operator=(const fvFieldDecomposer&) = delete;
 
 
-public:
-
     // Constructors
 
-        //- Construct from components
+        //- Construct without mappers, added later with reset()
+        fvFieldDecomposer
+        (
+            const Foam::zero,
+            const fvMesh& procMesh,
+            const labelList& faceAddressing,
+            const labelList& cellAddressing,
+            const labelList& boundaryAddressing
+        );
+
+        //- Construct from components using information from the complete mesh
         fvFieldDecomposer
         (
             const fvMesh& completeMesh,
@@ -252,13 +263,56 @@ public:
             const labelList& boundaryAddressing
         );
 
+        //- Construct from components without the complete mesh
+        fvFieldDecomposer
+        (
+            // Information about the complete mesh
+            const List<labelRange>& boundaryRanges,
+            const labelUList& faceOwner,
+            const labelUList& faceNeigbour,
+
+            // Addressing for processor mesh
+            const fvMesh& procMesh,
+            const labelList& faceAddressing,
+            const labelList& cellAddressing,
+            const labelList& boundaryAddressing
+        );
+
 
     //- Destructor
-    ~fvFieldDecomposer();
+    ~fvFieldDecomposer() = default;
 
 
     // Member Functions
 
+        //- True if no mappers have been allocated
+        bool empty() const;
+
+        //- Remove all mappers
+        void clear();
+
+        //- Reset mappers using information from the complete mesh
+        void reset(const fvMesh& completeMesh);
+
+        //- Reset mapper using information about the complete mesh
+        void reset
+        (
+            const List<labelRange>& boundaryRanges,
+            const labelUList& faceOwner,
+            const labelUList& faceNeigbour
+        );
+
+
+    // Mapping
+
+        //- Decompose internal field
+        template<class Type>
+        tmp<DimensionedField<Type, volMesh>>
+        decomposeField
+        (
+            const DimensionedField<Type, volMesh>& field
+        ) const;
+
         //- Decompose volume field
         template<class Type>
         tmp<GeometricField<Type, fvPatchField, volMesh>>
@@ -276,6 +330,7 @@ public:
             const GeometricField<Type, fvsPatchField, surfaceMesh>& field
         ) const;
 
+        //- Decompose list of fields
         template<class GeoField>
         void decomposeFields(const PtrList<GeoField>& fields) const;
 };
@@ -288,7 +343,7 @@ public:
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-    #include "fvFieldDecomposerDecomposeFields.C"
+    #include "fvFieldDecomposerFields.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/parallel/decompose/decompose/fvFieldDecomposerFields.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -34,6 +35,36 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class Type>
+Foam::tmp<Foam::DimensionedField<Type, Foam::volMesh>>
+Foam::fvFieldDecomposer::decomposeField
+(
+    const DimensionedField<Type, volMesh>& field
+) const
+{
+    // Create and map the internal field values
+    Field<Type> mappedField(field, cellAddressing_);
+
+    // Create the field for the processor
+    return
+        tmp<DimensionedField<Type, volMesh>>::New
+        (
+            IOobject
+            (
+                field.name(),
+                procMesh_.time().timeName(),
+                procMesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            procMesh_,
+            field.dimensions(),
+            std::move(mappedField)
+        );
+}
+
+
 template<class Type>
 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
 Foam::fvFieldDecomposer::decomposeField
@@ -338,9 +369,9 @@ void Foam::fvFieldDecomposer::decomposeFields
     const PtrList<GeoField>& fields
 ) const
 {
-    forAll(fields, fieldi)
+    for (const auto& fld : fields)
     {
-        decomposeField(fields[fieldi])().write();
+        decomposeField(fld)().write();
     }
 }
 

src/parallel/decompose/decompose/pointFieldDecomposer.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -76,6 +77,22 @@ Foam::pointFieldDecomposer::patchFieldDecomposer::patchFieldDecomposer
 }
 
 
+Foam::pointFieldDecomposer::pointFieldDecomposer
+(
+    const Foam::zero,
+    const pointMesh& procMesh,
+    const labelList& pointAddressing,
+    const labelList& boundaryAddressing
+)
+:
+    procMesh_(procMesh),
+    pointAddressing_(pointAddressing),
+    boundaryAddressing_(boundaryAddressing),
+    // Mappers
+    patchFieldDecomposerPtrs_()
+{}
+
+
 Foam::pointFieldDecomposer::pointFieldDecomposer
 (
     const pointMesh& completeMesh,
@@ -84,45 +101,60 @@ Foam::pointFieldDecomposer::pointFieldDecomposer
     const labelList& boundaryAddressing
 )
 :
-    completeMesh_(completeMesh),
-    procMesh_(procMesh),
-    pointAddressing_(pointAddressing),
-    boundaryAddressing_(boundaryAddressing),
-    patchFieldDecomposerPtrs_
+    pointFieldDecomposer
     (
-        procMesh_.boundary().size(),
-        static_cast<patchFieldDecomposer*>(nullptr)
+        zero{},
+        procMesh,
+        pointAddressing,
+        boundaryAddressing
     )
 {
+    reset(completeMesh);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::pointFieldDecomposer::empty() const
+{
+    return patchFieldDecomposerPtrs_.empty();
+}
+
+
+void Foam::pointFieldDecomposer::clear()
+{
+    patchFieldDecomposerPtrs_.clear();
+}
+
+
+void Foam::pointFieldDecomposer::reset
+(
+    const pointMesh& completeMesh
+)
+{
+    clear();
+    const label nMappers = procMesh_.boundary().size();
+    patchFieldDecomposerPtrs_.resize(nMappers);
+
     forAll(boundaryAddressing_, patchi)
     {
-        if (boundaryAddressing_[patchi] >= 0)
+        const label oldPatchi = boundaryAddressing_[patchi];
+
+        if (oldPatchi >= 0)
         {
-            patchFieldDecomposerPtrs_[patchi] = new patchFieldDecomposer
+            patchFieldDecomposerPtrs_.set
             (
-                completeMesh_.boundary()[boundaryAddressing_[patchi]],
+                patchi,
+                new patchFieldDecomposer
+                (
+                    completeMesh.boundary()[oldPatchi],
                     procMesh_.boundary()[patchi],
                     pointAddressing_
+                )
             );
         }
     }
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::pointFieldDecomposer::~pointFieldDecomposer()
-{
-    forAll(patchFieldDecomposerPtrs_, patchi)
-    {
-        if (patchFieldDecomposerPtrs_[patchi])
-        {
-            delete patchFieldDecomposerPtrs_[patchi];
-        }
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // ************************************************************************* //

src/parallel/decompose/decompose/pointFieldDecomposer.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,7 +32,7 @@ Description
 
 SourceFiles
     pointFieldDecomposer.C
-    pointFieldDecomposerDecomposeFields.C
+    pointFieldDecomposerFields.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -53,7 +54,6 @@ namespace Foam
 
 class pointFieldDecomposer
 {
-
 public:
 
         //- Point patch field decomposer class
@@ -107,10 +107,7 @@ public:
 
 private:
 
-    // Private data
-
-        //- Reference to complete mesh
-        const pointMesh& completeMesh_;
+    // Private Data
 
         //- Reference to processor mesh
         const pointMesh& procMesh_;
@@ -122,10 +119,10 @@ private:
         const labelList& boundaryAddressing_;
 
         //- List of patch field decomposers
-        List<patchFieldDecomposer*> patchFieldDecomposerPtrs_;
+        PtrList<patchFieldDecomposer> patchFieldDecomposerPtrs_;
 
 
-    // Private Member Functions
+public:
 
     //- No copy construct
     pointFieldDecomposer(const pointFieldDecomposer&) = delete;
@@ -134,10 +131,17 @@ private:
     void operator=(const pointFieldDecomposer&) = delete;
 
 
-public:
-
     // Constructors
 
+        //- Construct without mappers, added later with reset()
+        pointFieldDecomposer
+        (
+            const Foam::zero,
+            const pointMesh& procMesh,
+            const labelList& pointAddressing,
+            const labelList& boundaryAddressing
+        );
+
         //- Construct from components
         pointFieldDecomposer
         (
@@ -149,11 +153,23 @@ public:
 
 
     //- Destructor
-    ~pointFieldDecomposer();
+    ~pointFieldDecomposer() = default;
 
 
     // Member Functions
 
+        //- True if no mappers have been allocated
+        bool empty() const;
+
+        //- Remove all mappers
+        void clear();
+
+        //- Reset mappers using information from the complete mesh
+        void reset(const pointMesh& completeMesh);
+
+
+    // Mapping
+
         //- Decompose point field
         template<class Type>
         tmp<GeometricField<Type, pointPatchField, pointMesh>>
@@ -162,6 +178,7 @@ public:
             const GeometricField<Type, pointPatchField, pointMesh>&
         ) const;
 
+        //- Decompose list of fields
         template<class GeoField>
         void decomposeFields(const PtrList<GeoField>& fields) const;
 };
@@ -174,7 +191,7 @@ public:
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-    #include "pointFieldDecomposerDecomposeFields.C"
+    #include "pointFieldDecomposerFields.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/parallel/decompose/decompose/pointFieldDecomposerFields.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -46,7 +47,7 @@ Foam::pointFieldDecomposer::decomposeField
     // Create and map the patch field values
     forAll(boundaryAddressing_, patchi)
     {
-        if (patchFieldDecomposerPtrs_[patchi])
+        if (patchFieldDecomposerPtrs_.set(patchi))
         {
             patchFields.set
             (
@@ -56,7 +57,7 @@ Foam::pointFieldDecomposer::decomposeField
                     field.boundaryField()[boundaryAddressing_[patchi]],
                     procMesh_.boundary()[patchi],
                     DimensionedField<Type, pointMesh>::null(),
-                    *patchFieldDecomposerPtrs_[patchi]
+                    patchFieldDecomposerPtrs_[patchi]
                 )
             );
         }
@@ -75,9 +76,8 @@ Foam::pointFieldDecomposer::decomposeField
     }
 
     // Create the field for the processor
-    return tmp<GeometricField<Type, pointPatchField, pointMesh>>
-    (
-        new GeometricField<Type, pointPatchField, pointMesh>
+    return
+        tmp<GeometricField<Type, pointPatchField, pointMesh>>::New
         (
             IOobject
             (
@@ -92,7 +92,6 @@ Foam::pointFieldDecomposer::decomposeField
             field.dimensions(),
             internalField,
             patchFields
-        )
         );
 }
 
@@ -103,9 +102,9 @@ void Foam::pointFieldDecomposer::decomposeFields
     const PtrList<GeoField>& fields
 ) const
 {
-    forAll(fields, fieldi)
+    for (const auto& fld : fields)
     {
-        decomposeField(fields[fieldi])().write();
+        decomposeField(fld)().write();
     }
 }
 

src/parallel/decompose/faDecompose/Make/files

@@ -0,0 +1,4 @@
+faFieldDecomposer.C
+faMeshDecomposition.C
+
+LIB = $(FOAM_LIBBIN)/libfaDecompose

src/parallel/decompose/faDecompose/Make/options

@@ -0,0 +1,9 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+
+LIB_LIBS = \
+    -lfiniteVolume \
+    -lfiniteArea \
+    -lmeshTools

src/parallel/decompose/faDecompose/faFieldDecomposer.C

@@ -0,0 +1,357 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "faFieldDecomposer.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::faFieldDecomposer::patchFieldDecomposer::patchFieldDecomposer
+(
+    const label sizeBeforeMapping,
+    const labelUList& addressingSlice,
+    const label addressingOffset
+)
+:
+    sizeBeforeMapping_(sizeBeforeMapping),
+    directAddressing_(addressingSlice)
+{
+    forAll(directAddressing_, i)
+    {
+        // Subtract one to align addressing.
+        // directAddressing_[i] -= addressingOffset + 1;
+        // ZT, 12/Nov/2010
+        directAddressing_[i] -= addressingOffset;
+    }
+}
+
+
+Foam::faFieldDecomposer::processorAreaPatchFieldDecomposer::
+processorAreaPatchFieldDecomposer
+(
+    const label nTotalFaces,
+    const labelUList& owner,  // == mesh.edgeOwner()
+    const labelUList& neigh,  // == mesh.edgeNeighbour()
+    const labelUList& addressingSlice,
+    const scalarField& weights
+)
+:
+    sizeBeforeMapping_(nTotalFaces),
+    addressing_(addressingSlice.size()),
+    weights_(addressingSlice.size())
+{
+    forAll(addressing_, i)
+    {
+        // Subtract one to align addressing.
+        label ai = addressingSlice[i];
+//         label ai = mag(addressingSlice[i]) - 1;
+
+        if (ai < neigh.size())
+        {
+            // This is a regular edge. it has been an internal edge
+            // of the original mesh and now it has become a edge
+            // on the parallel boundary
+            addressing_[i].resize(2);
+            weights_[i].resize(2);
+
+            addressing_[i][0] = owner[ai];
+            addressing_[i][1] = neigh[ai];
+
+            if (ai < weights.size())
+            {
+                // Edge weights exist/are usable
+                weights_[i][0] = weights[ai];
+                weights_[i][1] = 1.0 - weights[ai];
+            }
+            else
+            {
+                // No edge weights. use equal weighting
+                weights_[i][0] = 0.5;
+                weights_[i][1] = 0.5;
+            }
+        }
+        else
+        {
+            // This is a edge that used to be on a cyclic boundary
+            // but has now become a parallel patch edge. I cannot
+            // do the interpolation properly (I would need to look
+            // up the different (edge) list of data), so I will
+            // just grab the value from the owner face
+
+            addressing_[i].resize(1);
+            weights_[i].resize(1);
+
+            addressing_[i][0] = owner[ai];
+
+            weights_[i][0] = 1.0;
+        }
+    }
+}
+
+
+Foam::faFieldDecomposer::processorAreaPatchFieldDecomposer::
+processorAreaPatchFieldDecomposer
+(
+    const faMesh& mesh,
+    const labelUList& addressingSlice
+)
+:
+    processorAreaPatchFieldDecomposer
+    (
+        mesh.nFaces(),
+        mesh.edgeOwner(),
+        mesh.edgeNeighbour(),
+        addressingSlice,
+        mesh.weights().internalField()
+    )
+{}
+
+
+Foam::faFieldDecomposer::processorEdgePatchFieldDecomposer::
+processorEdgePatchFieldDecomposer
+(
+    label sizeBeforeMapping,
+    const labelUList& addressingSlice
+)
+:
+    sizeBeforeMapping_(sizeBeforeMapping),
+    addressing_(addressingSlice.size()),
+    weights_(addressingSlice.size())
+{
+    forAll(addressing_, i)
+    {
+        addressing_[i].resize(1);
+        weights_[i].resize(1);
+
+        addressing_[i][0] = mag(addressingSlice[i]) - 1;
+        weights_[i][0] = sign(addressingSlice[i]);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::faFieldDecomposer::faFieldDecomposer
+(
+    const Foam::zero,
+    const faMesh& procMesh,
+    const labelList& edgeAddressing,
+    const labelList& faceAddressing,
+    const labelList& boundaryAddressing
+)
+:
+    procMesh_(procMesh),
+    edgeAddressing_(edgeAddressing),
+    faceAddressing_(faceAddressing),
+    boundaryAddressing_(boundaryAddressing),
+    // Mappers
+    patchFieldDecomposerPtrs_(),
+    processorAreaPatchFieldDecomposerPtrs_(),
+    processorEdgePatchFieldDecomposerPtrs_()
+{}
+
+
+Foam::faFieldDecomposer::faFieldDecomposer
+(
+    const faMesh& completeMesh,
+    const faMesh& procMesh,
+    const labelList& edgeAddressing,
+    const labelList& faceAddressing,
+    const labelList& boundaryAddressing
+)
+:
+    faFieldDecomposer
+    (
+        zero{},
+        procMesh,
+        edgeAddressing,
+        faceAddressing,
+        boundaryAddressing
+    )
+{
+    reset(completeMesh);
+}
+
+
+Foam::faFieldDecomposer::faFieldDecomposer
+(
+    const label nTotalFaces,
+    const List<labelRange>& boundaryRanges,
+    const labelUList& edgeOwner,
+    const labelUList& edgeNeigbour,
+
+    const faMesh& procMesh,
+    const labelList& edgeAddressing,
+    const labelList& faceAddressing,
+    const labelList& boundaryAddressing
+)
+:
+    faFieldDecomposer
+    (
+        zero{},
+        procMesh,
+        edgeAddressing,
+        faceAddressing,
+        boundaryAddressing
+    )
+{
+    reset(nTotalFaces, boundaryRanges, edgeOwner, edgeNeigbour);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::faFieldDecomposer::empty() const
+{
+    return patchFieldDecomposerPtrs_.empty();
+}
+
+
+void Foam::faFieldDecomposer::clear()
+{
+    patchFieldDecomposerPtrs_.clear();
+    processorAreaPatchFieldDecomposerPtrs_.clear();
+    processorEdgePatchFieldDecomposerPtrs_.clear();
+}
+
+
+void Foam::faFieldDecomposer::reset
+(
+    const label nTotalFaces,
+    const List<labelRange>& boundaryRanges,
+    const labelUList& edgeOwner,
+    const labelUList& edgeNeigbour
+)
+{
+    clear();
+    const label nMappers = procMesh_.boundary().size();
+
+    patchFieldDecomposerPtrs_.resize(nMappers);
+    processorAreaPatchFieldDecomposerPtrs_.resize(nMappers);
+    processorEdgePatchFieldDecomposerPtrs_.resize(nMappers);
+
+    forAll(boundaryAddressing_, patchi)
+    {
+        const label oldPatchi = boundaryAddressing_[patchi];
+        const faPatch& fap = procMesh_.boundary()[patchi];
+        const labelSubList localPatchSlice(fap.patchSlice(edgeAddressing_));
+
+        if (oldPatchi >= 0)
+        {
+            patchFieldDecomposerPtrs_.set
+            (
+                patchi,
+                new patchFieldDecomposer
+                (
+                    boundaryRanges[oldPatchi].size(),
+                    localPatchSlice,
+                    boundaryRanges[oldPatchi].start()
+                )
+            );
+        }
+        else
+        {
+            processorAreaPatchFieldDecomposerPtrs_.set
+            (
+                patchi,
+                new processorAreaPatchFieldDecomposer
+                (
+                    nTotalFaces,
+                    edgeOwner,
+                    edgeNeigbour,
+                    localPatchSlice
+                )
+            );
+
+            processorEdgePatchFieldDecomposerPtrs_.set
+            (
+                patchi,
+                new processorEdgePatchFieldDecomposer
+                (
+                    procMesh_.boundary()[patchi].size(),
+                    static_cast<const labelUList&>(localPatchSlice)
+                )
+            );
+        }
+    }
+}
+
+
+void Foam::faFieldDecomposer::reset(const faMesh& completeMesh)
+{
+    clear();
+    const label nMappers = procMesh_.boundary().size();
+    patchFieldDecomposerPtrs_.resize(nMappers);
+    processorAreaPatchFieldDecomposerPtrs_.resize(nMappers);
+    processorEdgePatchFieldDecomposerPtrs_.resize(nMappers);
+
+    forAll(boundaryAddressing_, patchi)
+    {
+        const label oldPatchi = boundaryAddressing_[patchi];
+        const faPatch& fap = procMesh_.boundary()[patchi];
+        const labelSubList localPatchSlice(fap.patchSlice(edgeAddressing_));
+
+        if (oldPatchi >= 0)
+        {
+            patchFieldDecomposerPtrs_.set
+            (
+                patchi,
+                new patchFieldDecomposer
+                (
+                    completeMesh.boundary()[oldPatchi].size(),
+                    localPatchSlice,
+                    completeMesh.boundary()[oldPatchi].start()
+                )
+            );
+        }
+        else
+        {
+            processorAreaPatchFieldDecomposerPtrs_.set
+            (
+                patchi,
+                new processorAreaPatchFieldDecomposer
+                (
+                    completeMesh,
+                    localPatchSlice
+                )
+            );
+
+            processorEdgePatchFieldDecomposerPtrs_.set
+            (
+                patchi,
+                new processorEdgePatchFieldDecomposer
+                (
+                    procMesh_.boundary()[patchi].size(),
+                    static_cast<const labelUList&>(localPatchSlice)
+                )
+            );
+        }
+    }
+}
+
+
+// ************************************************************************* //

src/parallel/decompose/faDecompose/faFieldDecomposer.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.
@@ -35,7 +36,7 @@ Author
 
 SourceFiles
     faFieldDecomposer.C
-    faFieldDecomposerDecomposeFields.C
+    faFieldDecomposerFields.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -51,6 +52,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward Declarations
 class IOobjectList;
 
 /*---------------------------------------------------------------------------*\
@@ -66,7 +68,7 @@ public:
         :
             public faPatchFieldMapper
         {
-            // Private data
+            // Private Data
 
                 label sizeBeforeMapping_;
                 labelList directAddressing_;
@@ -118,7 +120,7 @@ public:
         :
             public faPatchFieldMapper
         {
-            // Private data
+            // Private Data
 
                 label sizeBeforeMapping_;
                 labelListList addressing_;
@@ -126,7 +128,17 @@ public:
 
         public:
 
-            //- Construct given addressing
+            //- Construct addressing from details
+            processorAreaPatchFieldDecomposer
+            (
+                const label nTotalFaces,
+                const labelUList& edgeOwner,
+                const labelUList& edgeNeigbour,
+                const labelUList& addressingSlice,
+                const scalarField& edgeWeights = scalarField::null()
+            );
+
+            //- Construct given addressing from complete mesh
             processorAreaPatchFieldDecomposer
             (
                 const faMesh& mesh,
@@ -134,7 +146,7 @@ public:
             );
 
 
-            // Member functions
+            // Member Functions
 
                 label size() const
                 {
@@ -187,7 +199,7 @@ public:
             );
 
 
-            // Member functions
+            // Member Functions
 
                 label size() const
                 {
@@ -223,10 +235,7 @@ public:
 
 private:
 
-    // Private data
-
-        //- Reference to complete mesh
-        const faMesh& completeMesh_;
+    // Private Data
 
         //- Reference to processor mesh
         const faMesh& procMesh_;
@@ -241,12 +250,12 @@ private:
         const labelList& boundaryAddressing_;
 
         //- List of patch field decomposers
-        List<patchFieldDecomposer*> patchFieldDecomposerPtrs_;
+        PtrList<patchFieldDecomposer> patchFieldDecomposerPtrs_;
 
-        List<processorAreaPatchFieldDecomposer*>
+        PtrList<processorAreaPatchFieldDecomposer>
             processorAreaPatchFieldDecomposerPtrs_;
 
-        List<processorEdgePatchFieldDecomposer*>
+        PtrList<processorEdgePatchFieldDecomposer>
             processorEdgePatchFieldDecomposerPtrs_;
 
 
@@ -263,7 +272,17 @@ public:
 
     // Constructors
 
-        //- Construct from components
+        //- Construct without mappers, added later with reset()
+        faFieldDecomposer
+        (
+            const Foam::zero,
+            const faMesh& procMesh,
+            const labelList& edgeAddressing,
+            const labelList& faceAddressing,
+            const labelList& boundaryAddressing
+        );
+
+        //- Construct from components using information from the complete mesh
         faFieldDecomposer
         (
             const faMesh& completeMesh,
@@ -273,14 +292,50 @@ public:
             const labelList& boundaryAddressing
         );
 
+        //- Construct from components without the complete mesh
+        faFieldDecomposer
+        (
+            // Information about the complete mesh
+            const label nTotalFaces,
+            const List<labelRange>& boundaryRanges,
+            const labelUList& edgeOwner,
+            const labelUList& edgeNeigbour,
 
-    // Destructor
+            // Addressing for processor mesh
+            const faMesh& procMesh,
+            const labelList& edgeAddressing,
+            const labelList& faceAddressing,
+            const labelList& boundaryAddressing
+        );
 
-        ~faFieldDecomposer();
+
+    //- Destructor
+    ~faFieldDecomposer() = default;
 
 
     // Member Functions
 
+        //- True if no mappers have been allocated
+        bool empty() const;
+
+        //- Remove all mappers
+        void clear();
+
+        //- Reset mappers using information from the complete mesh
+        void reset(const faMesh& completeMesh);
+
+        //- Reset mapper using information about the complete mesh
+        void reset
+        (
+            const label nTotalFaces,
+            const List<labelRange>& boundaryRanges,
+            const labelUList& edgeOwner,
+            const labelUList& edgeNeigbour
+        );
+
+
+    // Mapping
+
         //- Decompose area field
         template<class Type>
         tmp<GeometricField<Type, faPatchField, areaMesh>>
@@ -299,6 +354,33 @@ public:
 
         template<class GeoField>
         void decomposeFields(const PtrList<GeoField>& fields) const;
+
+
+    // Reading helpers
+
+        //- Read the fields and hold on the pointer list
+        template
+        <
+            class Type,
+            template<class> class PatchField,
+            class GeoMesh
+        >
+        static void readFields
+        (
+            const typename GeoMesh::Mesh& mesh,
+            const IOobjectList& objects,
+            PtrList<GeometricField<Type, PatchField, GeoMesh>>& fields,
+            const bool readOldTime
+        );
+
+        //- Read fields and hold on the pointer list
+        template<class Mesh, class GeoField>
+        static void readFields
+        (
+            const Mesh& mesh,
+            const IOobjectList& objects,
+            PtrList<GeoField>& fields
+        );
 };
 
 
@@ -309,7 +391,8 @@ public:
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-#   include "faFieldDecomposerDecomposeFields.C"
+    #include "faFieldDecomposerFields.C"
+    #include "faFieldDecomposerReadFields.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/parallel/decompose/faDecompose/faFieldDecomposerFields.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.
@@ -29,16 +30,11 @@ License
 #include "processorFaPatchField.H"
 #include "processorFaePatchField.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
-tmp<GeometricField<Type, faPatchField, areaMesh>>
-faFieldDecomposer::decomposeField
+Foam::tmp<Foam::GeometricField<Type, Foam::faPatchField, Foam::areaMesh>>
+Foam::faFieldDecomposer::decomposeField
 (
     const GeometricField<Type, faPatchField, areaMesh>& field
 ) const
@@ -51,17 +47,19 @@ faFieldDecomposer::decomposeField
 
     forAll(boundaryAddressing_, patchi)
     {
-        if (boundaryAddressing_[patchi] >= 0)
+        const label oldPatchi = boundaryAddressing_[patchi];
+
+        if (oldPatchi >= 0)
         {
             patchFields.set
             (
                 patchi,
                 faPatchField<Type>::New
                 (
-                    field.boundaryField()[boundaryAddressing_[patchi]],
+                    field.boundaryField()[oldPatchi],
                     procMesh_.boundary()[patchi],
                     DimensionedField<Type, areaMesh>::null(),
-                    *patchFieldDecomposerPtrs_[patchi]
+                    patchFieldDecomposerPtrs_[patchi]
                 )
             );
         }
@@ -77,7 +75,7 @@ faFieldDecomposer::decomposeField
                     Field<Type>
                     (
                         field.internalField(),
-                        *processorAreaPatchFieldDecomposerPtrs_[patchi]
+                        processorAreaPatchFieldDecomposerPtrs_[patchi]
                     )
                 )
             );
@@ -85,9 +83,8 @@ faFieldDecomposer::decomposeField
     }
 
     // Create the field for the processor
-    return tmp<GeometricField<Type, faPatchField, areaMesh>>
-    (
-        new GeometricField<Type, faPatchField, areaMesh>
+    return
+        tmp<GeometricField<Type, faPatchField, areaMesh>>::New
         (
             IOobject
             (
@@ -101,14 +98,13 @@ faFieldDecomposer::decomposeField
             field.dimensions(),
             internalField,
             patchFields
-        )
         );
 }
 
 
 template<class Type>
-tmp<GeometricField<Type, faePatchField, edgeMesh>>
-faFieldDecomposer::decomposeField
+Foam::tmp<Foam::GeometricField<Type, Foam::faePatchField, Foam::edgeMesh>>
+Foam::faFieldDecomposer::decomposeField
 (
     const GeometricField<Type, faePatchField, edgeMesh>& field
 ) const
@@ -147,7 +143,7 @@ faFieldDecomposer::decomposeField
 
     forAll(field.boundaryField(), patchi)
     {
-        const Field<Type> & p = field.boundaryField()[patchi];
+        const Field<Type>& p = field.boundaryField()[patchi];
 
         const label patchStart = field.mesh().boundary()[patchi].start();
 
@@ -162,17 +158,19 @@ faFieldDecomposer::decomposeField
 
     forAll(boundaryAddressing_, patchi)
     {
-        if (boundaryAddressing_[patchi] >= 0)
+        const label oldPatchi = boundaryAddressing_[patchi];
+
+        if (oldPatchi >= 0)
         {
             patchFields.set
             (
                 patchi,
                 faePatchField<Type>::New
                 (
-                    field.boundaryField()[boundaryAddressing_[patchi]],
+                    field.boundaryField()[oldPatchi],
                     procMesh_.boundary()[patchi],
                     DimensionedField<Type, edgeMesh>::null(),
-                    *patchFieldDecomposerPtrs_[patchi]
+                    patchFieldDecomposerPtrs_[patchi]
                 )
             );
         }
@@ -188,7 +186,7 @@ faFieldDecomposer::decomposeField
                     Field<Type>
                     (
                         allEdgeField,
-                        *processorEdgePatchFieldDecomposerPtrs_[patchi]
+                        processorEdgePatchFieldDecomposerPtrs_[patchi]
                     )
                 )
             );
@@ -196,9 +194,8 @@ faFieldDecomposer::decomposeField
     }
 
     // Create the field for the processor
-    return tmp<GeometricField<Type, faePatchField, edgeMesh>>
-    (
-        new GeometricField<Type, faePatchField, edgeMesh>
+    return
+        tmp<GeometricField<Type, faePatchField, edgeMesh>>::New
         (
             IOobject
             (
@@ -212,13 +209,12 @@ faFieldDecomposer::decomposeField
             field.dimensions(),
             internalField,
             patchFields
-        )
         );
 }
 
 
 template<class GeoField>
-void faFieldDecomposer::decomposeFields
+void Foam::faFieldDecomposer::decomposeFields
 (
     const PtrList<GeoField>& fields
 ) const
@@ -230,8 +226,4 @@ void faFieldDecomposer::decomposeFields
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //

src/parallel/decompose/faDecompose/faFieldDecomposerReadFields.C

@@ -0,0 +1,99 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+    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 "faFieldDecomposer.H"
+#include "GeometricField.H"
+#include "IOobjectList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type, template<class> class PatchField, class GeoMesh>
+void Foam::faFieldDecomposer::readFields
+(
+    const typename GeoMesh::Mesh& mesh,
+    const IOobjectList& objects,
+    PtrList<GeometricField<Type, PatchField, GeoMesh>>& fields,
+    const bool readOldTime
+)
+{
+    typedef GeometricField<Type, PatchField, GeoMesh> GeoField;
+
+    // Search list of objects for fields of type GeoField
+    IOobjectList fieldObjects(objects.lookupClass<GeoField>());
+
+    /// // Remove the cellDist field
+    /// auto iter = fieldObjects.find("cellDist");
+    /// if (iter.found())
+    /// {
+    ///     fieldObjects.erase(iter);
+    /// }
+
+    // Use sorted set of names
+    // (different processors might read objects in different order)
+    const wordList masterNames(fieldObjects.sortedNames());
+
+    // Construct the fields
+    fields.resize(masterNames.size());
+
+    forAll(masterNames, i)
+    {
+        const IOobject& io = *fieldObjects[masterNames[i]];
+
+        fields.set(i, new GeoField(io, mesh, readOldTime));
+    }
+}
+
+
+template<class Mesh, class GeoField>
+void Foam::faFieldDecomposer::readFields
+(
+    const Mesh& mesh,
+    const IOobjectList& objects,
+    PtrList<GeoField>& fields
+)
+{
+    // Search list of objects for fields of type GeomField
+    IOobjectList fieldObjects(objects.lookupClass<GeoField>());
+
+    // Use sorted set of names
+    // (different processors might read objects in different order)
+    const wordList masterNames(fieldObjects.sortedNames());
+
+    // Construct the fields
+    fields.resize(masterNames.size());
+
+    forAll(masterNames, i)
+    {
+        const IOobject& io = *fieldObjects[masterNames[i]];
+
+        fields.set(i, new GeoField(io, mesh));
+    }
+}
+
+
+// ************************************************************************* //

src/parallel/decompose/faDecompose/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.
@@ -34,14 +34,16 @@ License
 #include "faMesh.H"
 #include "OSspecific.H"
 #include "Map.H"
+#include "SLList.H"
 #include "globalMeshData.H"
-#include "decompositionModel.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void Foam::faMeshDecomposition::distributeFaces()
 {
-    Info<< "\nCalculating distribution of faces" << endl;
+    const word& polyMeshRegionName = mesh().name();
+
+    Info<< "\nCalculating distribution of finiteArea faces" << endl;
 
     cpuTime decompositionTime;
 
@@ -58,7 +60,7 @@ void Foam::faMeshDecomposition::distributeFaces()
         (
             IOobject
             (
-                GeoMesh<polyMesh>::mesh_.name(),
+                polyMeshRegionName,
                 processorDb.timeName(),
                 processorDb
             )
@@ -160,22 +162,12 @@ void Foam::faMeshDecomposition::distributeFaces()
 Foam::faMeshDecomposition::faMeshDecomposition
 (
     const fvMesh& mesh,
-    const fileName& decompDictFile
+    const label nProcessors,
+    const dictionary& params
 )
 :
     faMesh(mesh),
-    decompDictFile_(decompDictFile),
-    nProcs_
-    (
-        decompositionMethod::nDomains
-        (
-            decompositionModel::New
-            (
-                mesh,
-                decompDictFile
-            )
-        )
-    ),
+    nProcs_(nProcessors),
     distributed_(false),
     hasGlobalFaceZones_(false),
     faceToProc_(nFaces()),
@@ -193,53 +185,73 @@ Foam::faMeshDecomposition::faMeshDecomposition
     procNeighbourProcessors_(nProcs_),
     procProcessorPatchSize_(nProcs_),
     procProcessorPatchStartIndex_(nProcs_),
-    globallySharedPoints_(0),
+    globallySharedPoints_(),
     cyclicParallel_(false)
 {
-    const decompositionModel& model = decompositionModel::New
-    (
-        mesh,
-        decompDictFile
-    );
-
-    model.readIfPresent("distributed", distributed_);
-    hasGlobalFaceZones_ = model.found("globalFaceZones");
+    updateParameters(params);
 }
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+void Foam::faMeshDecomposition::updateParameters
+(
+    const dictionary& params
+)
+{
+    params.readIfPresent("distributed", distributed_);
+    if (params.found("globalFaceZones"))
+    {
+        hasGlobalFaceZones_ = true;
+    }
+}
+
+
 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
-    {
-        List<SLList<label>> procFaceList(nProcs());
+    labelList nLocalFaces(nProcs_, Zero);
 
-        forAll(faceToProc_, faceI)
+    // Pass 1: determine local sizes, sanity check
+
+    forAll(faceToProc_, facei)
     {
-            if (faceToProc_[faceI] >= nProcs())
+        const label proci = faceToProc_[facei];
+
+        if (proci < 0 || proci >= nProcs_)
         {
-                FatalErrorIn("Finite area mesh decomposition")
-                    << "Impossible processor label " << faceToProc_[faceI]
-                    << "for face " << faceI
+            FatalErrorInFunction
+                << "Invalid processor label " << proci
+                << " for face " << facei << nl
                 << abort(FatalError);
         }
         else
         {
-                procFaceList[faceToProc_[faceI]].append(faceI);
+            ++nLocalFaces[proci];
         }
     }
 
-        // Convert linked lists into normal lists
-        forAll(procFaceList, procI)
+    // Adjust lengths
+    forAll(nLocalFaces, proci)
     {
-            procFaceAddressing_[procI] = procFaceList[procI];
+        procFaceAddressing_[proci].resize(nLocalFaces[proci]);
+        nLocalFaces[proci] = 0;  // restart list
     }
+
+    // Pass 2: fill in local lists
+    forAll(faceToProc_, facei)
+    {
+        const label proci = faceToProc_[facei];
+        const label localFacei = nLocalFaces[proci];
+        ++nLocalFaces[proci];
+
+        procFaceAddressing_[proci][localFacei] = facei;
     }
 
 
@@ -258,7 +270,7 @@ void Foam::faMeshDecomposition::decomposeMesh()
         (
             IOobject
             (
-                GeoMesh<polyMesh>::mesh_.name(),
+                polyMeshRegionName,
                 processorDb.timeName(),
                 processorDb
             )
@@ -314,7 +326,7 @@ void Foam::faMeshDecomposition::decomposeMesh()
                 fvFaceProcAddressingHash.find
                 (
                     faceLabels()[curProcFaceAddressing[faceI]] + 1
-                )();
+                ).val();
         }
 
         // create processor finite area mesh
@@ -324,38 +336,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 +379,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 +1063,7 @@ void Foam::faMeshDecomposition::decomposeMesh()
         (
             IOobject
             (
-                GeoMesh<polyMesh>::mesh_.name(),
+                polyMeshRegionName,
                 processorDb.timeName(),
                 processorDb
             )
@@ -1083,8 +1089,7 @@ void Foam::faMeshDecomposition::decomposeMesh()
             procProcessorPatchSize_[procI];
 
         labelListList& curPatchEdgeLabels = procPatchEdgeLabels_[procI];
-        curPatchEdgeLabels =
-            labelListList
+        curPatchEdgeLabels.resize
         (
             curPatchSize.size()
           + curProcessorPatchSize.size()
@@ -1137,8 +1142,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 +1181,7 @@ bool Foam::faMeshDecomposition::writeDecomposition()
         (
             IOobject
             (
-                GeoMesh<polyMesh>::mesh_.name(),
+                polyMeshRegionName,
                 processorDb.timeName(),
                 processorDb
             )
@@ -1195,7 +1201,7 @@ bool Foam::faMeshDecomposition::writeDecomposition()
         );
 
 
-        // create finite area mesh
+        // Create finite area mesh
         faMesh procMesh
         (
             procFvMesh,
@@ -1219,11 +1225,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 +1236,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 +1302,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
             {

src/parallel/decompose/faDecompose/faMeshDecomposition.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.
@@ -44,7 +45,6 @@ SourceFiles
 #include "fvMesh.H"
 #include "faMesh.H"
 #include "labelList.H"
-#include "SLList.H"
 #include "PtrList.H"
 #include "point.H"
 
@@ -59,10 +59,7 @@ class faMeshDecomposition
 :
     public faMesh
 {
-    // Private data
-
-        //- Optional non-standard file for decomposeParDict
-        const fileName decompDictFile_;
+    // Private Data
 
         //- Number of processors in decomposition
         label nProcs_;
@@ -86,7 +83,7 @@ class faMeshDecomposition
         labelList procNInternalEdges_;
 
         //- Edge labels for patches of processor meshes
-        List<List<List<label>>> procPatchEdgeLabels_;
+        List<labelListList> procPatchEdgeLabels_;
 
         //- Labels of points for each processor
         labelListList procPatchPointAddressing_;
@@ -131,18 +128,21 @@ class faMeshDecomposition
 
         void distributeFaces();
 
+
 public:
 
     // Constructors
 
         //- Construct from components.
+        //- Values will come from the volume decomposition
         //  \param mesh the fvMesh
-        //  \param decompDictFile optional non-standard location for the
-        //      decomposeParDict file
+        //  \param nProcessors the number of processors
+        //  \param params additional parameters, sent to updateParameters
         faMeshDecomposition
         (
             const fvMesh& mesh,
-            const fileName& decompDictFile = ""
+            const label nProcessors,
+            const dictionary& params = dictionary::null
         );
 
 
@@ -152,29 +152,63 @@ public:
 
     // Member Functions
 
+    // Settings
+
         //- Number of processor in decomposition
         label nProcs() const
         {
             return nProcs_;
         }
 
-        //- Is the decomposition data to be distributed for each processor
+        //- Is decomposition data to be distributed for each processor
         bool distributed() const
         {
             return distributed_;
         }
 
+        //- Change distributed flag
+        bool distributed(const bool on)
+        {
+            bool old(distributed_);
+            distributed_ = on;
+            return old;
+        }
+
+        //- Are global face zones used
+        bool useGlobalFaceZones() const
+        {
+            return distributed_;
+        }
+
+        //- Change global face zones flag
+        bool useGlobalFaceZones(const bool on)
+        {
+            bool old(hasGlobalFaceZones_);
+            hasGlobalFaceZones_ = on;
+            return old;
+        }
+
+        //- Update flags based on the decomposition model settings
+        //  Sets "distributed", detects presence of "globalFaceZones"
+        void updateParameters(const dictionary& params);
+
+
+    // Mappings
+
+        //- Face-processor decomposition labels
+        const labelList& faceToProc() const noexcept
+        {
+            return faceToProc_;
+        }
+
+
+    // Decompose
+
         //- Decompose mesh
         void decomposeMesh();
 
         //- Write decomposition
         bool writeDecomposition();
-
-        //- Cell-processor decomposition labels
-        const labelList& faceToProc() const
-        {
-            return faceToProc_;
-        }
 };
 
 

src/parallel/reconstruct/Allwmake

@@ -5,5 +5,6 @@ cd "${0%/*}" || exit                                # Run from this directory
 #------------------------------------------------------------------------------
 
 wmake $targetType reconstruct
+wmake $targetType faReconstruct
 
 #------------------------------------------------------------------------------

src/parallel/reconstruct/faReconstruct/Make/files

@@ -0,0 +1,5 @@
+processorFaMeshes.C
+faFieldReconstructor.C
+faMeshReconstructor.C
+
+LIB = $(FOAM_LIBBIN)/libfaReconstruct

src/parallel/reconstruct/faReconstruct/Make/options

@@ -0,0 +1,9 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -lfiniteArea \
+    -lfiniteVolume \
+    -lmeshTools

src/parallel/reconstruct/faReconstruct/faFieldReconstructor.H

@@ -35,7 +35,7 @@ Author
 
 SourceFiles
     faFieldReconstructor.C
-    faFieldReconstructorReconstructFields.C
+    faFieldReconstructorFields.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -194,7 +194,7 @@ public:
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-#   include "faFieldReconstructorReconstructFields.C"
+#   include "faFieldReconstructorFields.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/parallel/reconstruct/faReconstruct/faMeshReconstructor.C

@@ -0,0 +1,637 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "faMeshReconstructor.H"
+#include "globalIndex.H"
+#include "globalMeshData.H"
+#include "edgeHashes.H"
+#include "Time.H"
+#include "PstreamCombineReduceOps.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+int Foam::faMeshReconstructor::debug = 0;
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::faMeshReconstructor::calcAddressing
+(
+    const labelUList& fvFaceProcAddr
+)
+{
+    const globalIndex globalFaceNum(procMesh_.nFaces());
+
+    // ----------------------
+    // boundaryProcAddressing
+    //
+    // Trivial since non-processor boundary ordering is identical
+
+    const label nPatches = procMesh_.boundary().size();
+
+    faBoundaryProcAddr_ = identity(nPatches);
+
+    // Mark processor patches
+    for
+    (
+        label patchi = procMesh_.boundary().nNonProcessor();
+        patchi < nPatches;
+        ++patchi
+    )
+    {
+        faBoundaryProcAddr_[patchi] = -1;
+    }
+
+
+    // ------------------
+    // faceProcAddressing
+    //
+    // Transcribe/rewrite based on finiteVolume faceProcAddressing
+
+    faFaceProcAddr_ = procMesh_.faceLabels();
+
+    // From local faceLabels to proc values
+    for (label& facei : faFaceProcAddr_)
+    {
+        // Use finiteVolume info, ignoring face flips
+        facei = mag(fvFaceProcAddr[facei] - 1);
+    }
+
+
+    // Make global consistent.
+    // Starting at '0', without gaps in numbering etc.
+    // - the sorted order is the obvious solution
+    {
+        globalFaceNum.gather(faFaceProcAddr_, singlePatchFaceLabels_);
+
+        labelList order(Foam::sortedOrder(singlePatchFaceLabels_));
+
+        singlePatchFaceLabels_ = labelList(singlePatchFaceLabels_, order);
+
+        globalFaceNum.scatter(order, faFaceProcAddr_);
+    }
+
+    // Broadcast the same information everywhere
+    Pstream::scatter(singlePatchFaceLabels_);
+
+
+    // ------------------
+    // edgeProcAddressing
+    //
+    // This is more complicated.
+
+    // Create a single (serial) patch of the finiteArea mesh without a
+    // corresponding volume mesh, otherwise it would be the same as
+    // reconstructPar on the volume mesh (big, slow, etc).
+    //
+    // Do manual point-merge and relabeling to avoid dependency on
+    // finiteVolume pointProcAddressing
+
+    faPointProcAddr_.clear();  // Final size == procMesh_.nPoints()
+
+    // 1.
+    // - Topological point merge of the area meshes
+    // - use the local patch faces and points
+
+    // 2.
+    // - build a single (serial) patch of the finiteArea mesh only
+    //   without any point support from the volume mesh
+    // - it may be possible to skip this step, but not obvious how
+
+    // The collected faces are contiguous per processor, which gives us
+    // the possibility to easily identify their source by using the
+    // global face indexing (if needed).
+    // The ultimate face locations are handled with a separate ordering
+    // list.
+
+    const uindirectPrimitivePatch& procPatch = procMesh_.patch();
+
+
+    {
+        faceList singlePatchProcFaces;  // [proc0faces, proc1faces ...]
+        labelList uniqueMeshPointLabels;
+
+        // Local face points to compact merged point index
+        labelList pointToGlobal;
+
+        autoPtr<globalIndex> globalPointsPtr =
+            procMesh_.mesh().globalData().mergePoints
+            (
+                procPatch.meshPoints(),
+                procPatch.meshPointMap(),  // unused
+                pointToGlobal,
+                uniqueMeshPointLabels
+            );
+
+        // Gather faces, renumbered for the *merged* points
+        faceList tmpFaces(globalFaceNum.localSize());
+
+        forAll(tmpFaces, facei)
+        {
+            tmpFaces[facei] =
+                face(pointToGlobal, procPatch.localFaces()[facei]);
+        }
+
+        globalFaceNum.gather
+        (
+            tmpFaces,
+            singlePatchProcFaces,
+            UPstream::msgType(),
+            Pstream::commsTypes::scheduled
+        );
+
+        globalPointsPtr().gather
+        (
+            UIndirectList<point>
+            (
+                procPatch.points(),  // NB: mesh points (non-local)
+                uniqueMeshPointLabels
+            ),
+            singlePatchPoints_
+        );
+
+        // Transcribe into final assembled order
+        singlePatchFaces_.resize(singlePatchProcFaces.size());
+
+        // Target face ordering
+        labelList singlePatchProcAddr;
+        globalFaceNum.gather(faFaceProcAddr_, singlePatchProcAddr);
+
+        forAll(singlePatchProcAddr, facei)
+        {
+            const label targetFacei = singlePatchProcAddr[facei];
+            singlePatchFaces_[targetFacei] = singlePatchProcFaces[facei];
+        }
+
+        // Use initial equivalent "global" (serial) patch
+        // to establish the correct point and face walking order
+        //
+        // - only meaningful on master
+        const primitivePatch initialPatch
+        (
+            SubList<face>(singlePatchFaces_),
+            singlePatchPoints_
+        );
+
+        // Grab the faces/points in local point ordering
+        tmpFaces = initialPatch.localFaces();
+        pointField tmpPoints(initialPatch.localPoints());
+
+        // The meshPointMap is contiguous, so flatten as linear list
+        /// Map<label> mpm(initialPatch.meshPointMap());
+        labelList mpm(initialPatch.nPoints());
+        forAllConstIters(initialPatch.meshPointMap(), iter)
+        {
+            mpm[iter.key()] = iter.val();
+        }
+        Pstream::scatter(mpm);
+
+        // Rewrite pointToGlobal according to the correct point order
+        for (label& pointi : pointToGlobal)
+        {
+            pointi = mpm[pointi];
+        }
+
+        // Finalize. overwrite
+        faPointProcAddr_ = std::move(pointToGlobal);
+
+        singlePatchFaces_ = std::move(tmpFaces);
+        singlePatchPoints_ = std::move(tmpPoints);
+    }
+
+    // Broadcast the same information everywhere
+    Pstream::scatter(singlePatchFaces_);
+    Pstream::scatter(singlePatchPoints_);
+
+    // Now have enough global information to determine global edge mappings
+
+    // Equivalent "global" (serial) patch
+    const primitivePatch onePatch
+    (
+        SubList<face>(singlePatchFaces_),
+        singlePatchPoints_
+    );
+
+    faEdgeProcAddr_.clear();
+    faEdgeProcAddr_.resize(procMesh_.nEdges(), -1);
+
+    {
+        EdgeMap<label> globalEdgeMapping(2*onePatch.nEdges());
+
+        // Pass 1: edge-hash lookup with edges in "natural" patch order
+
+        // Can use local edges() directly without remap via meshPoints()
+        // since the previous sorting means that we are already working
+        // with faces that are in the local point order and even
+        // the points themselves are also in the local point order
+        for (const edge& e : onePatch.edges())
+        {
+            globalEdgeMapping.insert(e, globalEdgeMapping.size());
+        }
+
+        // Lookup proc local edges (in natural order) to global equivalent
+        for (label edgei = 0; edgei < procPatch.nEdges(); ++edgei)
+        {
+            const edge globalEdge(faPointProcAddr_, procPatch.edges()[edgei]);
+
+            const auto fnd = globalEdgeMapping.cfind(globalEdge);
+
+            if (fnd.found())
+            {
+                faEdgeProcAddr_[edgei] = fnd.val();
+            }
+            else
+            {
+                FatalErrorInFunction
+                    << "Failed to find edge " << globalEdge
+                    << " this indicates a programming error" << nl
+                    << exit(FatalError);
+            }
+        }
+    }
+
+    // Now have consistent global edge
+    // This of course would all be too easy.
+    // The finiteArea edge lists have been defined as their own sorted
+    // order with sublists etc.
+
+    // Gather edge ids for nonProcessor boundaries.
+    // These will also be in the serial geometry
+
+    Map<label> remapGlobal(2*onePatch.nEdges());
+    for (label edgei = 0; edgei < onePatch.nInternalEdges(); ++edgei)
+    {
+        remapGlobal.insert(edgei, remapGlobal.size());
+    }
+
+    //
+    singlePatchEdgeLabels_.clear();
+    singlePatchEdgeLabels_.resize(procMesh_.boundary().nNonProcessor());
+
+    forAll(singlePatchEdgeLabels_, patchi)
+    {
+        const faPatch& fap = procMesh_.boundary()[patchi];
+
+        labelList& patchEdgeLabels = singlePatchEdgeLabels_[patchi];
+        patchEdgeLabels = fap.edgeLabels();
+
+        // Renumber from local edges to global edges (natural order)
+        for (label& edgeId : patchEdgeLabels)
+        {
+            edgeId = faEdgeProcAddr_[edgeId];
+        }
+
+        // OR patchEdgeLabels =
+        // UIndirectList<label>(faEdgeProcAddr_, fap.edgeLabels());
+
+
+        // Collect from all processors
+        combineReduce
+        (
+            patchEdgeLabels,
+            ListOps::appendEqOp<label>()
+        );
+
+        // Sorted order will be the original non-decomposed order
+        Foam::sort(patchEdgeLabels);
+
+        for (const label sortedEdgei : patchEdgeLabels)
+        {
+            remapGlobal.insert(sortedEdgei, remapGlobal.size());
+        }
+    }
+
+    {
+        // Use the map to rewrite the local faEdgeProcAddr_
+
+        labelList newEdgeProcAddr(faEdgeProcAddr_);
+
+        label edgei = procMesh_.nInternalEdges();
+
+        for (const faPatch& fap : procMesh_.boundary())
+        {
+            for (const label patchEdgei : fap.edgeLabels())
+            {
+                const label globalEdgei = faEdgeProcAddr_[patchEdgei];
+
+                const auto fnd = remapGlobal.cfind(globalEdgei);
+                if (fnd.found())
+                {
+                    newEdgeProcAddr[edgei] = fnd.val();
+                    ++edgei;
+                }
+                else
+                {
+                    FatalErrorInFunction
+                        << "Failed to find edge " << globalEdgei
+                        << " this indicates a programming error" << nl
+                        << exit(FatalError);
+                }
+            }
+        }
+
+        faEdgeProcAddr_ = std::move(newEdgeProcAddr);
+    }
+}
+
+
+void Foam::faMeshReconstructor::initPatch() const
+{
+    serialPatchPtr_.reset
+    (
+        new primitivePatch
+        (
+            SubList<face>(singlePatchFaces_),
+            singlePatchPoints_
+        )
+    );
+}
+
+
+void Foam::faMeshReconstructor::createMesh()
+{
+    const Time& runTime = procMesh_.mesh().time();
+
+    // Time for non-parallel case (w/o functionObjects or libs)
+    serialRunTime_ = Time::New(runTime.globalPath().toAbsolute());
+
+
+    // Trivial polyMesh only containing points and faces.
+    // This is valid, provided we don't use it for anything complicated.
+
+    serialVolMesh_.reset
+    (
+        new polyMesh
+        (
+            IOobject
+            (
+                procMesh_.mesh().name(),  // Volume region name
+                procMesh_.mesh().facesInstance(),
+
+                *serialRunTime_,
+
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false  // not registered
+            ),
+            pointField(singlePatchPoints_),  // copy
+            faceList(singlePatchFaces_),     // copy
+            labelList(singlePatchFaces_.size(), Zero),  // faceOwner
+            labelList(),  // faceNeighbour
+            false  // no syncPar!
+        )
+    );
+
+    // A simple area-mesh with one-to-one mapping of faceLabels
+    serialAreaMesh_.reset
+    (
+        new faMesh
+        (
+            *serialVolMesh_,
+            identity(singlePatchFaces_.size())  // faceLabels
+        )
+    );
+
+    auto& completeMesh = *serialAreaMesh_;
+
+    // Add in non-processor boundary patches
+    PtrList<faPatch> completePatches(singlePatchEdgeLabels_.size());
+    forAll(completePatches, patchi)
+    {
+        const labelList& patchEdgeLabels = singlePatchEdgeLabels_[patchi];
+
+        const faPatch& fap = procMesh_.boundary()[patchi];
+
+        const label neiPolyPatchId = fap.ngbPolyPatchIndex();
+
+        completePatches.set
+        (
+            patchi,
+            fap.clone
+            (
+                completeMesh.boundary(),
+                patchEdgeLabels,
+                patchi,  // index
+                neiPolyPatchId
+            )
+        );
+    }
+
+    completeMesh.addFaPatches(completePatches);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::faMeshReconstructor::faMeshReconstructor
+(
+    const faMesh& procMesh
+)
+:
+    procMesh_(procMesh)
+{
+    if (!Pstream::parRun())
+    {
+        FatalErrorInFunction
+            << "Can only be called in parallel!!" << nl
+            << exit(FatalError);
+    }
+
+    // Require faceProcAddressing from finiteVolume decomposition
+    labelIOList fvFaceProcAddressing
+    (
+        IOobject
+        (
+            "faceProcAddressing",
+            procMesh_.mesh().facesInstance(),
+
+            // Or search?
+            // procMesh_.time().findInstance
+            // (
+            //     // Search for polyMesh face instance
+            //     // mesh.facesInstance()
+            //     procMesh_.mesh().meshDir(),
+            //     "faceProcAddressing"
+            // ),
+
+            polyMesh::meshSubDir,
+            procMesh_.mesh(),    // The polyMesh db
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false  // not registered
+        )
+    );
+
+    calcAddressing(fvFaceProcAddressing);
+}
+
+
+Foam::faMeshReconstructor::faMeshReconstructor
+(
+    const faMesh& procMesh,
+    const labelUList& fvFaceProcAddressing
+)
+:
+    procMesh_(procMesh)
+{
+    if (!Pstream::parRun())
+    {
+        FatalErrorInFunction
+            << "Can only be called in parallel!!" << nl
+            << exit(FatalError);
+    }
+
+    calcAddressing(fvFaceProcAddressing);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::faMeshReconstructor::~faMeshReconstructor()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::faMeshReconstructor::clearGeom()
+{
+    serialAreaMesh_.reset(nullptr);
+    serialVolMesh_.reset(nullptr);
+    serialRunTime_.reset(nullptr);
+}
+
+
+const Foam::primitivePatch& Foam::faMeshReconstructor::patch() const
+{
+    if (!serialPatchPtr_)
+    {
+        initPatch();
+    }
+
+    return *serialPatchPtr_;
+}
+
+
+Foam::primitivePatch& Foam::faMeshReconstructor::patch()
+{
+    if (!serialPatchPtr_)
+    {
+        initPatch();
+    }
+
+    return *serialPatchPtr_;
+}
+
+
+const Foam::faMesh& Foam::faMeshReconstructor::mesh() const
+{
+    if (!serialAreaMesh_)
+    {
+        const_cast<faMeshReconstructor&>(*this).createMesh();
+    }
+
+    return *serialAreaMesh_;
+}
+
+
+void Foam::faMeshReconstructor::writeAddressing() const
+{
+    writeAddressing(procMesh_.mesh().facesInstance());
+}
+
+
+void Foam::faMeshReconstructor::writeAddressing(const word& timeName) const
+{
+    // Write copies
+
+    IOobject ioAddr
+    (
+        "procAddressing",
+        timeName,
+        faMesh::meshSubDir,
+        procMesh_.mesh(),  // The polyMesh
+        IOobject::NO_READ,
+        IOobject::NO_WRITE,
+        false  // not registered
+    );
+
+    // boundaryProcAddressing
+    ioAddr.rename("boundaryProcAddressing");
+    labelIOList(ioAddr, faBoundaryProcAddr_).write();
+
+    // faceProcAddressing
+    ioAddr.rename("faceProcAddressing");
+    labelIOList(ioAddr, faFaceProcAddr_).write();
+
+    // pointProcAddressing
+    ioAddr.rename("pointProcAddressing");
+    labelIOList(ioAddr, faPointProcAddr_).write();
+
+    // edgeProcAddressing
+    ioAddr.rename("edgeProcAddressing");
+    labelIOList(ioAddr, faEdgeProcAddr_).write();
+}
+
+
+void Foam::faMeshReconstructor::writeMesh() const
+{
+    writeMesh(procMesh_.mesh().facesInstance());
+}
+
+
+void Foam::faMeshReconstructor::writeMesh(const word& timeName) const
+{
+    const faMesh& fullMesh = this->mesh();
+
+    const bool oldDistributed = fileHandler().distributed();
+    auto oldHandler = fileHandler(fileOperation::NewUncollated());
+    fileHandler().distributed(true);
+
+    if (Pstream::master())
+    {
+        const bool oldParRun = Pstream::parRun(false);
+
+        IOobject io(fullMesh.boundary());
+
+        io.rename("faceLabels");
+        labelIOList(io, singlePatchFaceLabels_).write();
+
+        fullMesh.boundary().write();
+
+        Pstream::parRun(oldParRun);
+    }
+
+    // Restore old settings
+    if (oldHandler)
+    {
+        fileHandler(std::move(oldHandler));
+    }
+    fileHandler().distributed(oldDistributed);
+}
+
+
+// ************************************************************************* //

src/parallel/reconstruct/faReconstruct/faMeshReconstructor.H

@@ -0,0 +1,222 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Class
+    Foam::faMeshReconstructor
+
+Description
+    A bare-bones reconstructor for finiteArea meshes when processor
+    meshes are available (in parallel) but an equivalent serial faMesh
+    is needed for reconstruction or decomposition.
+    In these situations, a serial version of the faMesh is needed,
+    but preferably without reconstructing the entire volume mesh.
+
+    It uses the finiteVolume faceProcAddressing in addition to
+    the geometric information available from the underlying polyMesh.
+
+    The resulting equivalent faMesh can be used for basic operations,
+    but caution should be exercised before attempting large operations.
+
+SourceFiles
+    faMeshReconstructor.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef faMeshReconstructor_H
+#define faMeshReconstructor_H
+
+#include "faMesh.H"
+#include "primitivePatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+class Time;
+
+/*---------------------------------------------------------------------------*\
+                     Class faMeshReconstructor Declaration
+\*---------------------------------------------------------------------------*/
+
+class faMeshReconstructor
+{
+    // Private Data
+
+        //- The processor-specific faMesh
+        const faMesh& procMesh_;
+
+
+    // Addressing
+
+        //- Processor face addressing, derived from finite volume information
+        labelList faFaceProcAddr_;
+
+        //- Processor boundary addressing
+        labelList faBoundaryProcAddr_;
+
+        //- Processor point addressing
+        labelList faPointProcAddr_;
+
+        //- Processor edge addressing
+        labelList faEdgeProcAddr_;
+
+
+    // Equivalent surface information
+
+        //- Faces labels for a single patch
+        labelList singlePatchFaceLabels_;
+
+        //- Faces for a single patch
+        faceList singlePatchFaces_;
+
+        //- Support points for a single patch
+        pointField singlePatchPoints_;
+
+        //- Lists of edge-labels (per edge patch) for the single patch
+        labelListList singlePatchEdgeLabels_;
+
+
+    // Demand-driven data
+
+        //- Primitive patch
+        mutable autoPtr<primitivePatch> serialPatchPtr_;
+
+        //- Time database for serial meshes
+        autoPtr<Time> serialRunTime_;
+
+        //- Dummy volume mesh, used for serial area mesh
+        autoPtr<polyMesh> serialVolMesh_;
+
+        //- Equivalent serial area mesh
+        autoPtr<faMesh> serialAreaMesh_;
+
+
+    // Private Member Functions
+
+        //- Calculate all addressing, using finiteVolume faceProcAddressing
+        void calcAddressing(const labelUList& fvFaceProcAddr);
+
+        //- Set primitive patch, removing any old one
+        void initPatch() const;
+
+        //- Create the serial geometry
+        void createMesh();
+
+        //- No copy construct
+        faMeshReconstructor(const faMeshReconstructor&) = delete;
+
+        //- No copy assignment
+        void operator=(const faMeshReconstructor&) = delete;
+
+
+public:
+
+    //- Debug flag
+    static int debug;
+
+
+    // Constructors
+
+        //- Construct from components
+        explicit faMeshReconstructor(const faMesh& procMesh);
+
+        //- Construct from components
+        faMeshReconstructor
+        (
+            const faMesh& procMesh,
+            const labelUList& fvFaceProcAddressing
+        );
+
+
+    //- Destructor
+    ~faMeshReconstructor();
+
+        void clearGeom();
+
+
+    // Member Functions
+
+        //- Processor point addressing
+        const labelList& pointProcAddressing() const noexcept
+        {
+            return faPointProcAddr_;
+        }
+
+        //- Processor edge addressing
+        const labelList& edgeProcAddressing() const noexcept
+        {
+            return faEdgeProcAddr_;
+        }
+
+        //- Processor face addressing
+        const labelList& faceProcAddressing() const noexcept
+        {
+            return faFaceProcAddr_;
+        }
+
+        //- Processor boundary addressing
+        const labelList& boundaryProcAddressing() const noexcept
+        {
+            return faBoundaryProcAddr_;
+        }
+
+
+        //- Serial equivalent patch
+        const primitivePatch& patch() const;
+
+        //- Serial equivalent patch
+        primitivePatch& patch();
+
+        //- Serial equivalent faMesh
+        const faMesh& mesh() const;
+
+
+    // Write
+
+        //- Write proc addressing at the polyMesh faceInstances time
+        void writeAddressing() const;
+
+        //- Write proc addressing at the given time
+        void writeAddressing(const word& timeName) const;
+
+        //- Write equivalent mesh information at the polyMesh faceInstances time
+        void writeMesh() const;
+
+        //- Write equivalent mesh information at the given time
+        void writeMesh(const word& timeName) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/parallel/reconstruct/faReconstruct/processorFaMeshes.C

@@ -145,116 +145,4 @@ Foam::processorFaMeshes::processorFaMeshes
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-// Foam::fvMesh::readUpdateState Foam::processorFaMeshes::readUpdate()
-// {
-//     fvMesh::readUpdateState stat = fvMesh::UNCHANGED;
-
-//     forAll(databases_, procI)
-//     {
-//         // Check if any new meshes need to be read.
-//         fvMesh::readUpdateState procStat = meshes_[procI].readUpdate();
-
-//         /*
-//         if (procStat != fvMesh::UNCHANGED)
-//         {
-//             Info<< "Processor " << procI
-//                 << " at time " << databases_[procI].timeName()
-//                 << " detected mesh change " << procStat
-//                 << endl;
-//         }
-//         */
-
-//         // Combine into overall mesh change status
-//         if (stat == fvMesh::UNCHANGED)
-//         {
-//             stat = procStat;
-//         }
-//         else
-//         {
-//             if (stat != procStat)
-//             {
-//                 FatalErrorIn("processorFaMeshes::readUpdate()")
-//                     << "Processor " << procI
-//                     << " has a different polyMesh at time "
-//                     << databases_[procI].timeName()
-//                     << " compared to any previous processors." << nl
-//                     << "Please check time " << databases_[procI].timeName()
-//                     << " directories on all processors for consistent"
-//                     << " mesh files."
-//                     << exit(FatalError);
-//             }
-//         }
-//     }
-
-//     if
-//     (
-//         stat == fvMesh::TOPO_CHANGE
-//      || stat == fvMesh::TOPO_PATCH_CHANGE
-//     )
-//     {
-//         // Reread all meshes and addressing
-//         read();
-//     }
-//     return stat;
-// }
-
-
-// void Foam::processorFaMeshes::reconstructPoints(fvMesh& mesh)
-// {
-//     // Read the field for all the processors
-//     PtrList<pointIOField> procsPoints(meshes_.size());
-
-//     forAll(meshes_, procI)
-//     {
-//         procsPoints.set
-//         (
-//             procI,
-//             new pointIOField
-//             (
-//                 IOobject
-//                 (
-//                     "points",
-//                     meshes_[procI].time().timeName(),
-//                     polyMesh::meshSubDir,
-//                     meshes_[procI],
-//                     IOobject::MUST_READ,
-//                     IOobject::NO_WRITE
-//                 )
-//             )
-//         );
-//     }
-
-//     // Create the new points
-//     vectorField newPoints(mesh.nPoints());
-
-//     forAll(meshes_, procI)
-//     {
-//         const vectorField& procPoints = procsPoints[procI];
-
-//         // Set the cell values in the reconstructed field
-
-//         const labelList& pointProcAddressingI = pointProcAddressing_[procI];
-
-//         if (pointProcAddressingI.size() != procPoints.size())
-//         {
-//             FatalErrorIn("processorFaMeshes")
-//                 << "problem :"
-//                 << " pointProcAddressingI:" << pointProcAddressingI.size()
-//                 << " procPoints:" << procPoints.size()
-//                 << abort(FatalError);
-//         }
-
-//         forAll(pointProcAddressingI, pointI)
-//         {
-//             newPoints[pointProcAddressingI[pointI]] = procPoints[pointI];
-//         }
-//     }
-
-//     mesh.movePoints(newPoints);
-//     mesh.write();
-// }
-
-
 // ************************************************************************* //

src/parallel/reconstruct/faReconstruct/processorFaMeshes.H

@@ -27,7 +27,7 @@ Class
     Foam::processorFaMeshes
 
 Description
-    Container for processor mesh addressing.
+    Container for finite-area processor mesh addressing.
 
 Author
     Zeljko Tukovic, FSB Zagreb
@@ -41,8 +41,8 @@ SourceFiles
 #define processorFaMeshes_H
 
 #include "PtrList.H"
-#include "fvMesh.H"
 #include "faMesh.H"
+#include "fvMesh.H"
 #include "IOobjectList.H"
 #include "labelIOList.H"
 
@@ -100,29 +100,31 @@ public:
 
     // Member Functions
 
-        //- Update the meshes based on the mesh files saved in
-        //  time directories
-//         fvMesh::readUpdateState readUpdate();
-
-        //- Reconstruct point position after motion in parallel
-//         void reconstructPoints(faMesh& mesh);
+        const PtrList<faMesh>& meshes() const
+        {
+            return meshes_;
+        }
 
         PtrList<faMesh>& meshes()
         {
             return meshes_;
         }
+
         const PtrList<labelIOList>& pointProcAddressing() const
         {
             return pointProcAddressing_;
         }
+
         PtrList<labelIOList>& edgeProcAddressing()
         {
             return edgeProcAddressing_;
         }
+
         const PtrList<labelIOList>& faceProcAddressing() const
         {
             return faceProcAddressing_;
         }
+
         const PtrList<labelIOList>& boundaryProcAddressing() const
         {
             return boundaryProcAddressing_;

src/thermophysicalModels/radiation/fvOptions/radiation/radiation.H

@@ -80,7 +80,7 @@ namespace fv
 
 class radiation
 :
-    public option
+    public fv::option
 {
     // Private Data