From 0a241800ce5ad902f5b9173a12d697a943a8b6c5 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Wed, 5 May 2021 10:52:15 +0200
Subject: [PATCH] ENH: collect and cleanup decompose and reconstruct methods
 (#2084)

- new faDecompose and faReconstruct libraries

- provide common readFields in the faDecompose library
---
 .../decomposePar/Make/files                   |   5 +-
 .../decomposePar/Make/options                 |   4 +-
 .../decomposePar/decomposePar.C               | 247 +++++++++---------
 .../decomposePar/domainDecomposition.C        |  57 ++--
 .../domainDecompositionDistribute.C           |   6 +-
 .../decomposePar/lagrangianFieldDecomposer.H  |   6 +-
 ...ds.C => lagrangianFieldDecomposerFields.C} |   0
 .../decomposePar/readFields.C                 |  19 +-
 .../reconstructPar/Make/files                 |   2 -
 .../reconstructPar/Make/options               |   6 +-
 src/parallel/decompose/Allwclean              |   1 +
 src/parallel/decompose/Allwmake               |   1 +
 src/parallel/decompose/decompose/Make/files   |   3 +
 .../decompose/decompose}/dimFieldDecomposer.C |  21 +-
 .../decompose/decompose}/dimFieldDecomposer.H |  48 ++--
 .../decompose/dimFieldDecomposerFields.C      |  11 +-
 .../decompose/decompose/fvFieldDecomposer.C   |  15 +-
 .../decompose/decompose/fvFieldDecomposer.H   |  41 +--
 ...poseFields.C => fvFieldDecomposerFields.C} |  31 +++
 .../decompose}/pointFieldDecomposer.C         |  40 +--
 .../decompose}/pointFieldDecomposer.H         |  24 +-
 .../decompose/pointFieldDecomposerFields.C    |  13 +-
 src/parallel/decompose/faDecompose/Make/files |   4 +
 .../decompose/faDecompose/Make/options        |   9 +
 .../faDecompose}/faFieldDecomposer.C          | 101 +++----
 .../faDecompose}/faFieldDecomposer.H          |  45 +++-
 .../faDecompose/faFieldDecomposerFields.C     |  58 ++--
 .../faDecompose/faFieldDecomposerReadFields.C |  99 +++++++
 .../faDecompose}/faMeshDecomposition.C        | 103 ++++----
 .../faDecompose}/faMeshDecomposition.H        |  66 +++--
 src/parallel/reconstruct/Allwmake             |   1 +
 .../reconstruct/faReconstruct/Make/files      |   4 +
 .../reconstruct/faReconstruct/Make/options    |   9 +
 .../faReconstruct}/faFieldReconstructor.C     |   0
 .../faReconstruct}/faFieldReconstructor.H     |   4 +-
 .../faFieldReconstructorFields.C              |   0
 .../faReconstruct}/processorFaMeshes.C        | 112 --------
 .../faReconstruct}/processorFaMeshes.H        |  18 +-
 38 files changed, 652 insertions(+), 582 deletions(-)
 rename applications/utilities/parallelProcessing/decomposePar/{lagrangianFieldDecomposerDecomposeFields.C => lagrangianFieldDecomposerFields.C} (100%)
 rename {applications/utilities/parallelProcessing/decomposePar => src/parallel/decompose/decompose}/dimFieldDecomposer.C (84%)
 rename {applications/utilities/parallelProcessing/decomposePar => src/parallel/decompose/decompose}/dimFieldDecomposer.H (77%)
 rename applications/utilities/parallelProcessing/decomposePar/dimFieldDecomposerDecomposeFields.C => src/parallel/decompose/decompose/dimFieldDecomposerFields.C (94%)
 rename src/parallel/decompose/decompose/{fvFieldDecomposerDecomposeFields.C => fvFieldDecomposerFields.C} (92%)
 rename {applications/utilities/parallelProcessing/decomposePar => src/parallel/decompose/decompose}/pointFieldDecomposer.C (80%)
 rename {applications/utilities/parallelProcessing/decomposePar => src/parallel/decompose/decompose}/pointFieldDecomposer.H (91%)
 rename applications/utilities/parallelProcessing/decomposePar/pointFieldDecomposerDecomposeFields.C => src/parallel/decompose/decompose/pointFieldDecomposerFields.C (93%)
 create mode 100644 src/parallel/decompose/faDecompose/Make/files
 create mode 100644 src/parallel/decompose/faDecompose/Make/options
 rename {applications/utilities/parallelProcessing/decomposePar => src/parallel/decompose/faDecompose}/faFieldDecomposer.C (69%)
 rename {applications/utilities/parallelProcessing/decomposePar => src/parallel/decompose/faDecompose}/faFieldDecomposer.H (87%)
 rename applications/utilities/parallelProcessing/decomposePar/faFieldDecomposerDecomposeFields.C => src/parallel/decompose/faDecompose/faFieldDecomposerFields.C (81%)
 create mode 100644 src/parallel/decompose/faDecompose/faFieldDecomposerReadFields.C
 rename {applications/utilities/parallelProcessing/decomposePar => src/parallel/decompose/faDecompose}/faMeshDecomposition.C (96%)
 rename {applications/utilities/parallelProcessing/decomposePar => src/parallel/decompose/faDecompose}/faMeshDecomposition.H (78%)
 create mode 100644 src/parallel/reconstruct/faReconstruct/Make/files
 create mode 100644 src/parallel/reconstruct/faReconstruct/Make/options
 rename {applications/utilities/parallelProcessing/reconstructPar => src/parallel/reconstruct/faReconstruct}/faFieldReconstructor.C (100%)
 rename {applications/utilities/parallelProcessing/reconstructPar => src/parallel/reconstruct/faReconstruct}/faFieldReconstructor.H (98%)
 rename applications/utilities/parallelProcessing/reconstructPar/faFieldReconstructorReconstructFields.C => src/parallel/reconstruct/faReconstruct/faFieldReconstructorFields.C (100%)
 rename {applications/utilities/parallelProcessing/reconstructPar => src/parallel/reconstruct/faReconstruct}/processorFaMeshes.C (57%)
 rename {applications/utilities/parallelProcessing/reconstructPar => src/parallel/reconstruct/faReconstruct}/processorFaMeshes.H (92%)

diff --git a/applications/utilities/parallelProcessing/decomposePar/Make/files b/applications/utilities/parallelProcessing/decomposePar/Make/files
index 03abb7247f2..d8bc47a3151 100644
--- a/applications/utilities/parallelProcessing/decomposePar/Make/files
+++ b/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
diff --git a/applications/utilities/parallelProcessing/decomposePar/Make/options b/applications/utilities/parallelProcessing/decomposePar/Make/options
index 104ca56d7c6..4acbf00e42e 100644
--- a/applications/utilities/parallelProcessing/decomposePar/Make/options
+++ b/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
diff --git a/applications/utilities/parallelProcessing/decomposePar/decomposePar.C b/applications/utilities/parallelProcessing/decomposePar/decomposePar.C
index 7c9784ce8db..292e694f71e 100644
--- a/applications/utilities/parallelProcessing/decomposePar/decomposePar.C
+++ b/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,84 +1476,59 @@ int main(int argc, char *argv[])
                         //     procAddressing
                         //     (
                         //         procMeshList,
-                        //         procI,
+                        //         proci,
                         //         "boundaryProcAddressing",
                         //         boundaryProcAddressingList
                         //     );
 
-                        labelIOList faceProcAddressing
-                        (
-                            IOobject
-                            (
-                                "faceProcAddressing",
-                                "constant",
-                                procMesh.meshSubDir,
-                                procFvMesh,
-                                IOobject::MUST_READ,
-                                IOobject::NO_WRITE
-                            )
-                        );
+                        // Addressing from faMesh (not polyMesh) meshSubDir
 
-                        labelIOList boundaryProcAddressing
-                        (
-                            IOobject
+                        autoPtr<labelIOList> tfaceProcAddr =
+                            faProcAddressing
                             (
-                                "boundaryProcAddressing",
-                                "constant",
-                                procMesh.meshSubDir,
                                 procFvMesh,
-                                IOobject::MUST_READ,
-                                IOobject::NO_WRITE
-                            )
-                        );
-
-                        // FA fields
-                        if
-                        (
-                            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
-                                )
+                                "faceProcAddressing",
+                                runTime.constant()
                             );
+                        auto& faceProcAddressing = *tfaceProcAddr;
 
-                            faFieldDecomposer fieldDecomposer
+                        autoPtr<labelIOList> tboundaryProcAddr =
+                            faProcAddressing
                             (
-                                aMesh,
-                                procMesh,
-                                edgeProcAddressing,
-                                faceProcAddressing,
-                                boundaryProcAddressing
+                                procFvMesh,
+                                "boundaryProcAddressing",
+                                runTime.constant()
                             );
+                        auto& boundaryProcAddressing = *tboundaryProcAddr;
 
-                            fieldDecomposer.decomposeFields(areaScalarFields);
-                            fieldDecomposer.decomposeFields(areaVectorFields);
-                            fieldDecomposer.decomposeFields
+                        autoPtr<labelIOList> tedgeProcAddr =
+                            faProcAddressing
                             (
-                                areaSphericalTensorFields
-                            );
-                            fieldDecomposer.decomposeFields
-                            (
-                                areaSymmTensorFields
+                                procFvMesh,
+                                "edgeProcAddressing",
+                                runTime.constant()
                             );
-                            fieldDecomposer.decomposeFields(areaTensorFields);
+                        const auto& edgeProcAddressing = *tedgeProcAddr;
 
-                            fieldDecomposer.decomposeFields(edgeScalarFields);
-                        }
+                        faFieldDecomposer fieldDecomposer
+                        (
+                            aMesh,
+                            procMesh,
+                            edgeProcAddressing,
+                            faceProcAddressing,
+                            boundaryProcAddressing
+                        );
+
+                        fieldDecomposer.decomposeFields(areaScalarFields);
+                        fieldDecomposer.decomposeFields(areaVectorFields);
+                        fieldDecomposer.decomposeFields
+                        (
+                            areaSphericalTensorFields
+                        );
+                        fieldDecomposer.decomposeFields(areaSymmTensorFields);
+                        fieldDecomposer.decomposeFields(areaTensorFields);
+
+                        fieldDecomposer.decomposeFields(edgeScalarFields);
                     }
                 }
             }
diff --git a/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C b/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C
index 40efd741514..f3e7c2dfc7f 100644
--- a/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C
+++ b/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
             (
-                procMesh.boundaryMesh(),
                 nPatches,
-                patchMapper.directAddressing(),
-                curPatchStarts[patchi]
-            ).ptr();
+                meshPatch.clone
+                (
+                    procMesh.boundaryMesh(),
+                    nPatches,
+                    patchMapper.directAddressing(),
+                    curPatchStarts[patchi]
+                )
+            );
 
             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,7 +489,8 @@ 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;
             }
         }
 
diff --git a/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDistribute.C b/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDistribute.C
index aedbbb34e59..a448982b1f6 100644
--- a/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDistribute.C
+++ b/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"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposer.H b/applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposer.H
index 31e209c871b..18a1746cb4e 100644
--- a/applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposer.H
+++ b/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
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposerDecomposeFields.C b/applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposerFields.C
similarity index 100%
rename from applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposerDecomposeFields.C
rename to applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposerFields.C
diff --git a/applications/utilities/parallelProcessing/decomposePar/readFields.C b/applications/utilities/parallelProcessing/decomposePar/readFields.C
index 231981bcee1..1b0c948c9c1 100644
--- a/applications/utilities/parallelProcessing/decomposePar/readFields.C
+++ b/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)
     {
diff --git a/applications/utilities/parallelProcessing/reconstructPar/Make/files b/applications/utilities/parallelProcessing/reconstructPar/Make/files
index 660641033b0..37acabe62db 100644
--- a/applications/utilities/parallelProcessing/reconstructPar/Make/files
+++ b/applications/utilities/parallelProcessing/reconstructPar/Make/files
@@ -1,5 +1,3 @@
-processorFaMeshes.C
-faFieldReconstructor.C
 reconstructPar.C
 
 EXE = $(FOAM_APPBIN)/reconstructPar
diff --git a/applications/utilities/parallelProcessing/reconstructPar/Make/options b/applications/utilities/parallelProcessing/reconstructPar/Make/options
index 75f62d17291..a526d94fbc4 100644
--- a/applications/utilities/parallelProcessing/reconstructPar/Make/options
+++ b/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
diff --git a/src/parallel/decompose/Allwclean b/src/parallel/decompose/Allwclean
index e6b78d10c68..a0a7d6335f8 100755
--- a/src/parallel/decompose/Allwclean
+++ b/src/parallel/decompose/Allwclean
@@ -8,6 +8,7 @@ wclean kahipDecomp
 wclean scotchDecomp
 wclean decompositionMethods
 wclean decompose
+wclean faDecompose
 
 ./Allwclean-mpi
 
diff --git a/src/parallel/decompose/Allwmake b/src/parallel/decompose/Allwmake
index c9bc6b610af..b1681aa0b7f 100755
--- a/src/parallel/decompose/Allwmake
+++ b/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
diff --git a/src/parallel/decompose/decompose/Make/files b/src/parallel/decompose/decompose/Make/files
index 225a27dfad8..e55de05d2e5 100644
--- a/src/parallel/decompose/decompose/Make/files
+++ b/src/parallel/decompose/decompose/Make/files
@@ -1,5 +1,8 @@
 decompositionInformation.C
 decompositionModel.C
+
+dimFieldDecomposer.C
 fvFieldDecomposer.C
+pointFieldDecomposer.C
 
 LIB = $(FOAM_LIBBIN)/libdecompose
diff --git a/applications/utilities/parallelProcessing/decomposePar/dimFieldDecomposer.C b/src/parallel/decompose/decompose/dimFieldDecomposer.C
similarity index 84%
rename from applications/utilities/parallelProcessing/decomposePar/dimFieldDecomposer.C
rename to src/parallel/decompose/decompose/dimFieldDecomposer.C
index 391d2bf5c85..09dc649ad34 100644
--- a/applications/utilities/parallelProcessing/decomposePar/dimFieldDecomposer.C
+++ b/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,31 @@ 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& completeMesh,
+    const fvMesh& procMesh,
+    const labelList& faceAddressing,
+    const labelList& cellAddressing
+)
+:
+    //UNUSED: completeMesh_(completeMesh),
+    procMesh_(procMesh),
+    //UNUSED: faceAddressing_(faceAddressing),
+    cellAddressing_(cellAddressing)
 {}
 
 
diff --git a/applications/utilities/parallelProcessing/decomposePar/dimFieldDecomposer.H b/src/parallel/decompose/decompose/dimFieldDecomposer.H
similarity index 77%
rename from applications/utilities/parallelProcessing/decomposePar/dimFieldDecomposer.H
rename to src/parallel/decompose/decompose/dimFieldDecomposer.H
index 26710ce5e8a..ec80c621cd9 100644
--- a/applications/utilities/parallelProcessing/decomposePar/dimFieldDecomposer.H
+++ b/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,69 +47,68 @@ 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
-
-        //- No copy construct
-        dimFieldDecomposer(const dimFieldDecomposer&) = delete;
+public:
 
-        //- No copy assignment
-        void operator=(const dimFieldDecomposer&) = delete;
+    //- No copy construct
+    dimFieldDecomposer(const dimFieldDecomposer&) = delete;
 
+    //- No copy assignment
+    void operator=(const dimFieldDecomposer&) = delete;
 
-public:
 
     // Constructors
 
-        //- Construct from components
+        //- Construct from minimal components
+        dimFieldDecomposer
+        (
+            const fvMesh& procMesh,
+            const labelList& cellAddressing
+        );
+
+        //- Construct from components with API as per fvFieldDecomposer
         dimFieldDecomposer
         (
-            const fvMesh& completeMesh,
+            const fvMesh& completeMesh,         //!< unused
             const fvMesh& procMesh,
-            const labelList& faceAddressing,
+            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
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/utilities/parallelProcessing/decomposePar/dimFieldDecomposerDecomposeFields.C b/src/parallel/decompose/decompose/dimFieldDecomposerFields.C
similarity index 94%
rename from applications/utilities/parallelProcessing/decomposePar/dimFieldDecomposerDecomposeFields.C
rename to src/parallel/decompose/decompose/dimFieldDecomposerFields.C
index 2d4c2bf9878..e4785f414b9 100644
--- a/applications/utilities/parallelProcessing/decomposePar/dimFieldDecomposerDecomposeFields.C
+++ b/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,9 +55,8 @@ Foam::dimFieldDecomposer::decomposeField
             ),
             procMesh_,
             field.dimensions(),
-            mappedField
-        )
-    );
+            std::move(mappedField)
+        );
 }
 
 
diff --git a/src/parallel/decompose/decompose/fvFieldDecomposer.C b/src/parallel/decompose/decompose/fvFieldDecomposer.C
index 8cc5acdf683..077c177fc13 100644
--- a/src/parallel/decompose/decompose/fvFieldDecomposer.C
+++ b/src/parallel/decompose/decompose/fvFieldDecomposer.C
@@ -27,7 +27,6 @@ License
 
 #include "fvFieldDecomposer.H"
 
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::fvFieldDecomposer::patchFieldDecomposer::patchFieldDecomposer
@@ -135,11 +134,12 @@ Foam::fvFieldDecomposer::fvFieldDecomposer
 {
     forAll(boundaryAddressing_, patchi)
     {
+        const label oldPatchi = boundaryAddressing_[patchi];
         const fvPatch& fvp = procMesh_.boundary()[patchi];
 
         if
         (
-            boundaryAddressing_[patchi] >= 0
+            oldPatchi >= 0
         && !isA<processorLduInterface>(procMesh.boundary()[patchi])
         )
         {
@@ -149,10 +149,7 @@ Foam::fvFieldDecomposer::fvFieldDecomposer
                 new patchFieldDecomposer
                 (
                     fvp.patchSlice(faceAddressing_),
-                    completeMesh_.boundaryMesh()
-                    [
-                        boundaryAddressing_[patchi]
-                    ].start()
+                    completeMesh_.boundaryMesh()[oldPatchi].start()
                 )
             );
         }
@@ -202,10 +199,4 @@ Foam::fvFieldDecomposer::fvFieldDecomposer
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::fvFieldDecomposer::~fvFieldDecomposer()
-{}
-
-
 // ************************************************************************* //
diff --git a/src/parallel/decompose/decompose/fvFieldDecomposer.H b/src/parallel/decompose/decompose/fvFieldDecomposer.H
index 34a0904a580..d9358224e63 100644
--- a/src/parallel/decompose/decompose/fvFieldDecomposer.H
+++ b/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,7 +109,7 @@ public:
         :
             public fvPatchFieldMapper
         {
-            // Private data
+            // Private Data
 
                 labelList directAddressing_;
 
@@ -124,7 +123,7 @@ public:
             );
 
 
-            // Member functions
+            // Member Functions
 
                 label size() const
                 {
@@ -199,7 +198,7 @@ public:
 
 private:
 
-    // Private data
+    // Private Data
 
         //- Reference to complete mesh
         const fvMesh& completeMesh_;
@@ -225,20 +224,17 @@ private:
         PtrList<processorSurfacePatchFieldDecomposer>
             processorSurfacePatchFieldDecomposerPtrs_;
 
-
         PtrList<scalarField> faceSign_;
 
 
-    // Private Member Functions
-
-        //- No copy construct
-        fvFieldDecomposer(const fvFieldDecomposer&) = delete;
+public:
 
-        //- No copy assignment
-        void operator=(const fvFieldDecomposer&) = delete;
+    //- No copy construct
+    fvFieldDecomposer(const fvFieldDecomposer&) = delete;
 
+    //- No copy assignment
+    void operator=(const fvFieldDecomposer&) = delete;
 
-public:
 
     // Constructors
 
@@ -254,11 +250,19 @@ public:
 
 
     //- Destructor
-    ~fvFieldDecomposer();
+    ~fvFieldDecomposer() = default;
 
 
     // Member Functions
 
+        //- 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 +280,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 +293,7 @@ public:
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-    #include "fvFieldDecomposerDecomposeFields.C"
+    #include "fvFieldDecomposerFields.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/parallel/decompose/decompose/fvFieldDecomposerDecomposeFields.C b/src/parallel/decompose/decompose/fvFieldDecomposerFields.C
similarity index 92%
rename from src/parallel/decompose/decompose/fvFieldDecomposerDecomposeFields.C
rename to src/parallel/decompose/decompose/fvFieldDecomposerFields.C
index 2e9bcd9a5a5..c625c30196b 100644
--- a/src/parallel/decompose/decompose/fvFieldDecomposerDecomposeFields.C
+++ b/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
diff --git a/applications/utilities/parallelProcessing/decomposePar/pointFieldDecomposer.C b/src/parallel/decompose/decompose/pointFieldDecomposer.C
similarity index 80%
rename from applications/utilities/parallelProcessing/decomposePar/pointFieldDecomposer.C
rename to src/parallel/decompose/decompose/pointFieldDecomposer.C
index a4c40ee9806..aaa794e218c 100644
--- a/applications/utilities/parallelProcessing/decomposePar/pointFieldDecomposer.C
+++ b/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.
@@ -88,41 +89,28 @@ Foam::pointFieldDecomposer::pointFieldDecomposer
     procMesh_(procMesh),
     pointAddressing_(pointAddressing),
     boundaryAddressing_(boundaryAddressing),
-    patchFieldDecomposerPtrs_
-    (
-        procMesh_.boundary().size(),
-        static_cast<patchFieldDecomposer*>(nullptr)
-    )
+
+    patchFieldDecomposerPtrs_(procMesh_.boundary().size())
 {
     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]],
-                procMesh_.boundary()[patchi],
-                pointAddressing_
+                patchi,
+                new patchFieldDecomposer
+                (
+                    completeMesh_.boundary()[oldPatchi],
+                    procMesh_.boundary()[patchi],
+                    pointAddressing_
+                )
             );
         }
     }
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::pointFieldDecomposer::~pointFieldDecomposer()
-{
-    forAll(patchFieldDecomposerPtrs_, patchi)
-    {
-        if (patchFieldDecomposerPtrs_[patchi])
-        {
-            delete patchFieldDecomposerPtrs_[patchi];
-        }
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // ************************************************************************* //
diff --git a/applications/utilities/parallelProcessing/decomposePar/pointFieldDecomposer.H b/src/parallel/decompose/decompose/pointFieldDecomposer.H
similarity index 91%
rename from applications/utilities/parallelProcessing/decomposePar/pointFieldDecomposer.H
rename to src/parallel/decompose/decompose/pointFieldDecomposer.H
index f59b4f33b20..a7379026ad9 100644
--- a/applications/utilities/parallelProcessing/decomposePar/pointFieldDecomposer.H
+++ b/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,7 +107,7 @@ public:
 
 private:
 
-    // Private data
+    // Private Data
 
         //- Reference to complete mesh
         const pointMesh& completeMesh_;
@@ -122,19 +122,17 @@ private:
         const labelList& boundaryAddressing_;
 
         //- List of patch field decomposers
-        List<patchFieldDecomposer*> patchFieldDecomposerPtrs_;
-
+        PtrList<patchFieldDecomposer> patchFieldDecomposerPtrs_;
 
-    // Private Member Functions
 
-        //- No copy construct
-        pointFieldDecomposer(const pointFieldDecomposer&) = delete;
+public:
 
-        //- No copy assignment
-        void operator=(const pointFieldDecomposer&) = delete;
+    //- No copy construct
+    pointFieldDecomposer(const pointFieldDecomposer&) = delete;
 
+    //- No copy assignment
+    void operator=(const pointFieldDecomposer&) = delete;
 
-public:
 
     // Constructors
 
@@ -149,7 +147,7 @@ public:
 
 
     //- Destructor
-    ~pointFieldDecomposer();
+    ~pointFieldDecomposer() = default;
 
 
     // Member Functions
@@ -174,7 +172,7 @@ public:
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-    #include "pointFieldDecomposerDecomposeFields.C"
+    #include "pointFieldDecomposerFields.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/utilities/parallelProcessing/decomposePar/pointFieldDecomposerDecomposeFields.C b/src/parallel/decompose/decompose/pointFieldDecomposerFields.C
similarity index 93%
rename from applications/utilities/parallelProcessing/decomposePar/pointFieldDecomposerDecomposeFields.C
rename to src/parallel/decompose/decompose/pointFieldDecomposerFields.C
index bdbea0cb43c..e9c1b72a481 100644
--- a/applications/utilities/parallelProcessing/decomposePar/pointFieldDecomposerDecomposeFields.C
+++ b/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,8 +92,7 @@ Foam::pointFieldDecomposer::decomposeField
             field.dimensions(),
             internalField,
             patchFields
-        )
-    );
+        );
 }
 
 
diff --git a/src/parallel/decompose/faDecompose/Make/files b/src/parallel/decompose/faDecompose/Make/files
new file mode 100644
index 00000000000..af2a87bf961
--- /dev/null
+++ b/src/parallel/decompose/faDecompose/Make/files
@@ -0,0 +1,4 @@
+faFieldDecomposer.C
+faMeshDecomposition.C
+
+LIB = $(FOAM_LIBBIN)/libfaDecompose
diff --git a/src/parallel/decompose/faDecompose/Make/options b/src/parallel/decompose/faDecompose/Make/options
new file mode 100644
index 00000000000..b5d6f595dc9
--- /dev/null
+++ b/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
diff --git a/applications/utilities/parallelProcessing/decomposePar/faFieldDecomposer.C b/src/parallel/decompose/faDecompose/faFieldDecomposer.C
similarity index 69%
rename from applications/utilities/parallelProcessing/decomposePar/faFieldDecomposer.C
rename to src/parallel/decompose/faDecompose/faFieldDecomposer.C
index 2ee76091d3e..122eaa5fb38 100644
--- a/applications/utilities/parallelProcessing/decomposePar/faFieldDecomposer.C
+++ b/src/parallel/decompose/faDecompose/faFieldDecomposer.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.
@@ -27,14 +28,9 @@ License
 
 #include "faFieldDecomposer.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-faFieldDecomposer::patchFieldDecomposer::patchFieldDecomposer
+Foam::faFieldDecomposer::patchFieldDecomposer::patchFieldDecomposer
 (
     const label sizeBeforeMapping,
     const labelUList& addressingSlice,
@@ -54,7 +50,7 @@ faFieldDecomposer::patchFieldDecomposer::patchFieldDecomposer
 }
 
 
-faFieldDecomposer::processorAreaPatchFieldDecomposer::
+Foam::faFieldDecomposer::processorAreaPatchFieldDecomposer::
 processorAreaPatchFieldDecomposer
 (
     const faMesh& mesh,
@@ -108,7 +104,7 @@ processorAreaPatchFieldDecomposer
 }
 
 
-faFieldDecomposer::processorEdgePatchFieldDecomposer::
+Foam::faFieldDecomposer::processorEdgePatchFieldDecomposer::
 processorEdgePatchFieldDecomposer
 (
     label sizeBeforeMapping,
@@ -130,7 +126,7 @@ processorEdgePatchFieldDecomposer
 }
 
 
-faFieldDecomposer::faFieldDecomposer
+Foam::faFieldDecomposer::faFieldDecomposer
 (
     const faMesh& completeMesh,
     const faMesh& procMesh,
@@ -144,47 +140,43 @@ faFieldDecomposer::faFieldDecomposer
     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)
-    )
+
+    patchFieldDecomposerPtrs_(procMesh_.boundary().size()),
+    processorAreaPatchFieldDecomposerPtrs_(procMesh_.boundary().size()),
+    processorEdgePatchFieldDecomposerPtrs_(procMesh_.boundary().size())
 {
     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]].size(),
-                procMesh_.boundary()[patchi].patchSlice(edgeAddressing_),
-//                 completeMesh_.boundaryMesh()
-                completeMesh_.boundary()
-                [
-                    boundaryAddressing_[patchi]
-                ].start()
+                patchi,
+                new patchFieldDecomposer
+                (
+                    completeMesh_.boundary()[oldPatchi].size(),
+                    procMesh_.boundary()[patchi].patchSlice(edgeAddressing_),
+                    completeMesh_.boundary()[oldPatchi].start()
+                )
             );
         }
         else
         {
-            processorAreaPatchFieldDecomposerPtrs_[patchi] =
+            processorAreaPatchFieldDecomposerPtrs_.set
+            (
+                patchi,
                 new processorAreaPatchFieldDecomposer
                 (
                     completeMesh_,
                     procMesh_.boundary()[patchi].patchSlice(edgeAddressing_)
-                );
+                )
+            );
 
-            processorEdgePatchFieldDecomposerPtrs_[patchi] =
+            processorEdgePatchFieldDecomposerPtrs_.set
+            (
+                patchi,
                 new processorEdgePatchFieldDecomposer
                 (
                     procMesh_.boundary()[patchi].size(),
@@ -195,44 +187,11 @@ faFieldDecomposer::faFieldDecomposer
                             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
-
 // ************************************************************************* //
diff --git a/applications/utilities/parallelProcessing/decomposePar/faFieldDecomposer.H b/src/parallel/decompose/faDecompose/faFieldDecomposer.H
similarity index 87%
rename from applications/utilities/parallelProcessing/decomposePar/faFieldDecomposer.H
rename to src/parallel/decompose/faDecompose/faFieldDecomposer.H
index 16928957bb5..18008a20773 100644
--- a/applications/utilities/parallelProcessing/decomposePar/faFieldDecomposer.H
+++ b/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_;
@@ -241,12 +243,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_;
 
 
@@ -275,8 +277,7 @@ public:
 
 
     // Destructor
-
-        ~faFieldDecomposer();
+    ~faFieldDecomposer() = default;
 
 
     // Member Functions
@@ -299,6 +300,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 +337,8 @@ public:
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-#   include "faFieldDecomposerDecomposeFields.C"
+    #include "faFieldDecomposerFields.C"
+    #include "faFieldDecomposerReadFields.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/utilities/parallelProcessing/decomposePar/faFieldDecomposerDecomposeFields.C b/src/parallel/decompose/faDecompose/faFieldDecomposerFields.C
similarity index 81%
rename from applications/utilities/parallelProcessing/decomposePar/faFieldDecomposerDecomposeFields.C
rename to src/parallel/decompose/faDecompose/faFieldDecomposerFields.C
index ad0bebd67d8..1b7eaea40f2 100644
--- a/applications/utilities/parallelProcessing/decomposePar/faFieldDecomposerDecomposeFields.C
+++ b/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
-
 // ************************************************************************* //
diff --git a/src/parallel/decompose/faDecompose/faFieldDecomposerReadFields.C b/src/parallel/decompose/faDecompose/faFieldDecomposerReadFields.C
new file mode 100644
index 00000000000..43529e3172d
--- /dev/null
+++ b/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));
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/parallelProcessing/decomposePar/faMeshDecomposition.C b/src/parallel/decompose/faDecompose/faMeshDecomposition.C
similarity index 96%
rename from applications/utilities/parallelProcessing/decomposePar/faMeshDecomposition.C
rename to src/parallel/decompose/faDecompose/faMeshDecomposition.C
index bdbb92dc79e..aecad1c65e8 100644
--- a/applications/utilities/parallelProcessing/decomposePar/faMeshDecomposition.C
+++ b/src/parallel/decompose/faDecompose/faMeshDecomposition.C
@@ -34,8 +34,8 @@ License
 #include "faMesh.H"
 #include "OSspecific.H"
 #include "Map.H"
+#include "SLList.H"
 #include "globalMeshData.H"
-#include "decompositionModel.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -43,7 +43,7 @@ void Foam::faMeshDecomposition::distributeFaces()
 {
     const word& polyMeshRegionName = mesh().name();
 
-    Info<< "\nCalculating distribution of faces" << endl;
+    Info<< "\nCalculating distribution of finiteArea faces" << endl;
 
     cpuTime decompositionTime;
 
@@ -162,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()),
@@ -195,22 +185,28 @@ 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
@@ -220,32 +216,44 @@ void Foam::faMeshDecomposition::decomposeMesh()
 
     Info<< "\nDistributing faces to processors" << endl;
 
-    // Memory management
+    labelList nLocalFaces(nProcs_, Zero);
+
+    // Pass 1: determine local sizes, sanity check
+
+    forAll(faceToProc_, facei)
     {
-        List<SLList<label>> procFaceList(nProcs());
+        const label proci = faceToProc_[facei];
 
-        forAll(faceToProc_, faceI)
+        if (proci < 0 || proci >= nProcs_)
         {
-            if (faceToProc_[faceI] >= nProcs())
-            {
-                FatalErrorInFunction
-                    << "Impossible processor label " << faceToProc_[faceI]
-                    << "for face " << faceI << nl
-                    << abort(FatalError);
-            }
-            else
-            {
-                procFaceList[faceToProc_[faceI]].append(faceI);
-            }
+            FatalErrorInFunction
+                << "Invalid processor label " << proci
+                << " for face " << facei << nl
+                << abort(FatalError);
         }
-
-        // Convert linked lists into normal lists
-        forAll(procFaceList, procI)
+        else
         {
-            procFaceAddressing_[procI] = procFaceList[procI];
+            ++nLocalFaces[proci];
         }
     }
 
+    // Adjust lengths
+    forAll(nLocalFaces, 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;
+    }
+
 
     // Find processor mesh faceLabels and ...
 
@@ -1081,12 +1089,11 @@ void Foam::faMeshDecomposition::decomposeMesh()
             procProcessorPatchSize_[procI];
 
         labelListList& curPatchEdgeLabels = procPatchEdgeLabels_[procI];
-        curPatchEdgeLabels =
-            labelListList
-            (
-                curPatchSize.size()
-              + curProcessorPatchSize.size()
-            );
+        curPatchEdgeLabels.resize
+        (
+            curPatchSize.size()
+          + curProcessorPatchSize.size()
+        );
 
         forAll(curPatchSize, patchI)
         {
diff --git a/applications/utilities/parallelProcessing/decomposePar/faMeshDecomposition.H b/src/parallel/decompose/faDecompose/faMeshDecomposition.H
similarity index 78%
rename from applications/utilities/parallelProcessing/decomposePar/faMeshDecomposition.H
rename to src/parallel/decompose/faDecompose/faMeshDecomposition.H
index 3c0526820c3..44a5323353a 100644
--- a/applications/utilities/parallelProcessing/decomposePar/faMeshDecomposition.H
+++ b/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_;
         }
 
-        //- Decompose mesh
-        void decomposeMesh();
+        //- Change distributed flag
+        bool distributed(const bool on)
+        {
+            bool old(distributed_);
+            distributed_ = on;
+            return old;
+        }
 
-        //- Write decomposition
-        bool writeDecomposition();
+        //- Are global face zones used
+        bool useGlobalFaceZones() const
+        {
+            return distributed_;
+        }
 
-        //- Cell-processor decomposition labels
-        const labelList& faceToProc() const
+        //- 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();
 };
 
 
diff --git a/src/parallel/reconstruct/Allwmake b/src/parallel/reconstruct/Allwmake
index 77b67cd1ecb..583901ce7c0 100755
--- a/src/parallel/reconstruct/Allwmake
+++ b/src/parallel/reconstruct/Allwmake
@@ -5,5 +5,6 @@ cd "${0%/*}" || exit                                # Run from this directory
 #------------------------------------------------------------------------------
 
 wmake $targetType reconstruct
+wmake $targetType faReconstruct
 
 #------------------------------------------------------------------------------
diff --git a/src/parallel/reconstruct/faReconstruct/Make/files b/src/parallel/reconstruct/faReconstruct/Make/files
new file mode 100644
index 00000000000..4d7028d13cb
--- /dev/null
+++ b/src/parallel/reconstruct/faReconstruct/Make/files
@@ -0,0 +1,4 @@
+processorFaMeshes.C
+faFieldReconstructor.C
+
+LIB = $(FOAM_LIBBIN)/libfaReconstruct
diff --git a/src/parallel/reconstruct/faReconstruct/Make/options b/src/parallel/reconstruct/faReconstruct/Make/options
new file mode 100644
index 00000000000..2526f4eb0e6
--- /dev/null
+++ b/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
diff --git a/applications/utilities/parallelProcessing/reconstructPar/faFieldReconstructor.C b/src/parallel/reconstruct/faReconstruct/faFieldReconstructor.C
similarity index 100%
rename from applications/utilities/parallelProcessing/reconstructPar/faFieldReconstructor.C
rename to src/parallel/reconstruct/faReconstruct/faFieldReconstructor.C
diff --git a/applications/utilities/parallelProcessing/reconstructPar/faFieldReconstructor.H b/src/parallel/reconstruct/faReconstruct/faFieldReconstructor.H
similarity index 98%
rename from applications/utilities/parallelProcessing/reconstructPar/faFieldReconstructor.H
rename to src/parallel/reconstruct/faReconstruct/faFieldReconstructor.H
index 77f312191da..0c1e0005732 100644
--- a/applications/utilities/parallelProcessing/reconstructPar/faFieldReconstructor.H
+++ b/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
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/utilities/parallelProcessing/reconstructPar/faFieldReconstructorReconstructFields.C b/src/parallel/reconstruct/faReconstruct/faFieldReconstructorFields.C
similarity index 100%
rename from applications/utilities/parallelProcessing/reconstructPar/faFieldReconstructorReconstructFields.C
rename to src/parallel/reconstruct/faReconstruct/faFieldReconstructorFields.C
diff --git a/applications/utilities/parallelProcessing/reconstructPar/processorFaMeshes.C b/src/parallel/reconstruct/faReconstruct/processorFaMeshes.C
similarity index 57%
rename from applications/utilities/parallelProcessing/reconstructPar/processorFaMeshes.C
rename to src/parallel/reconstruct/faReconstruct/processorFaMeshes.C
index 980d652ce41..645a166900b 100644
--- a/applications/utilities/parallelProcessing/reconstructPar/processorFaMeshes.C
+++ b/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();
-// }
-
-
 // ************************************************************************* //
diff --git a/applications/utilities/parallelProcessing/reconstructPar/processorFaMeshes.H b/src/parallel/reconstruct/faReconstruct/processorFaMeshes.H
similarity index 92%
rename from applications/utilities/parallelProcessing/reconstructPar/processorFaMeshes.H
rename to src/parallel/reconstruct/faReconstruct/processorFaMeshes.H
index b302a287d76..342beb83751 100644
--- a/applications/utilities/parallelProcessing/reconstructPar/processorFaMeshes.H
+++ b/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_;
-- 
GitLab