From 71ccf51ba59a438410a71bfd456841891acceaae Mon Sep 17 00:00:00 2001
From: Will Bainbridge <http://cfd.direct>
Date: Thu, 20 Jul 2023 11:42:18 +0100
Subject: [PATCH] decomposePar, reconstructPar: Do all regions simultaneously

DecomposePar and reconstructPar now interleave the processing of
multiple regions. This means that with the -allRegions option, the
earlier times are completed in their entirety before later times are
considered. It also lets regions to access each other during
decomposition and reconstruction, which will be important for
non-conformal region interfaces.

To aid interpretation of the log, region prefixing is now used by both
utilities in the same way as is done by foamMultiRun.

DecomposePar has been overhauled so that it matches reconstructPar much
more closely, both in terms of output and of iteration sequence. All
meshes and addressing are loaded simultaneously and each field is
considered in turn. Previously, all the fields were loaded, and each
process and addressing set was considered in turn. This new strategy
optimises memory usage for cases with lots of fields.
---
 .../decomposePar/Make/files                   |    1 -
 .../decomposePar/decomposePar.C               | 1143 ++++++-----------
 .../decomposePar/dimFieldDecomposer.C         |   49 -
 .../decomposePar/dimFieldDecomposer.H         |  121 --
 .../dimFieldDecomposerDecomposeFields.C       |   75 --
 .../decomposePar/fvFieldDecomposer.C          |   76 +-
 .../decomposePar/fvFieldDecomposer.H          |   74 +-
 .../fvFieldDecomposerDecomposeFields.C        |  364 ------
 .../decomposePar/fvFieldDecomposerTemplates.C |  550 ++++++++
 .../decomposePar/lagrangianFieldDecomposer.C  |  200 ++-
 .../decomposePar/lagrangianFieldDecomposer.H  |  111 +-
 ...lagrangianFieldDecomposerDecomposeFields.C |  213 ---
 .../lagrangianFieldDecomposerTemplates.C      |  119 ++
 .../decomposePar/pointFieldDecomposer.C       |   62 +-
 .../decomposePar/pointFieldDecomposer.H       |   46 +-
 .../pointFieldDecomposerDecomposeFields.C     |  111 --
 .../pointFieldDecomposerTemplates.C           |  163 +++
 .../decomposePar/readFields.C                 |  107 --
 .../decomposePar/readFields.H                 |   75 --
 .../reconstructPar/Make/files                 |    2 +-
 .../reconstructPar/fvFieldReconstructor.C     |   27 +-
 .../reconstructPar/fvFieldReconstructor.H     |   87 +-
 ...elds.C => fvFieldReconstructorTemplates.C} |  269 ++--
 .../lagrangianFieldReconstructor.C            |  124 ++
 .../lagrangianFieldReconstructor.H            |  165 +++
 .../lagrangianFieldReconstructorTemplates.C   |  179 +++
 .../reconstructPar/pointFieldReconstructor.C  |   25 +-
 .../reconstructPar/pointFieldReconstructor.H  |   36 +-
 ...s.C => pointFieldReconstructorTemplates.C} |   48 +-
 .../reconstructPar/reconstructLagrangian.H    |  117 --
 .../reconstructLagrangianFields.C             |  273 ----
 .../reconstructLagrangianPositions.C          |   83 --
 .../reconstructPar/reconstructPar.C           |  788 +++++-------
 src/parallel/parallel/Make/files              |    3 +
 src/parallel/parallel/domainDecomposition.C   |  743 +----------
 src/parallel/parallel/domainDecomposition.H   |  104 +-
 .../parallel/domainDecompositionDecompose.C   |   57 +-
 .../domainDecompositionNonConformal.C         |  744 +++++++++++
 .../parallel/domainDecompositionReconstruct.C |   40 +-
 .../parallel/multiDomainDecomposition.C       |  180 +++
 .../parallel/multiDomainDecomposition.H       |  156 +++
 src/parallel/parallel/multiRegionPrefixer.C   |  123 ++
 src/parallel/parallel/multiRegionPrefixer.H   |  241 ++++
 43 files changed, 4246 insertions(+), 4028 deletions(-)
 delete mode 100644 applications/utilities/parallelProcessing/decomposePar/dimFieldDecomposer.C
 delete mode 100644 applications/utilities/parallelProcessing/decomposePar/dimFieldDecomposer.H
 delete mode 100644 applications/utilities/parallelProcessing/decomposePar/dimFieldDecomposerDecomposeFields.C
 delete mode 100644 applications/utilities/parallelProcessing/decomposePar/fvFieldDecomposerDecomposeFields.C
 create mode 100644 applications/utilities/parallelProcessing/decomposePar/fvFieldDecomposerTemplates.C
 delete mode 100644 applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposerDecomposeFields.C
 create mode 100644 applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposerTemplates.C
 delete mode 100644 applications/utilities/parallelProcessing/decomposePar/pointFieldDecomposerDecomposeFields.C
 create mode 100644 applications/utilities/parallelProcessing/decomposePar/pointFieldDecomposerTemplates.C
 delete mode 100644 applications/utilities/parallelProcessing/decomposePar/readFields.C
 delete mode 100644 applications/utilities/parallelProcessing/decomposePar/readFields.H
 rename applications/utilities/parallelProcessing/reconstructPar/{fvFieldReconstructorReconstructFields.C => fvFieldReconstructorTemplates.C} (84%)
 create mode 100644 applications/utilities/parallelProcessing/reconstructPar/lagrangianFieldReconstructor.C
 create mode 100644 applications/utilities/parallelProcessing/reconstructPar/lagrangianFieldReconstructor.H
 create mode 100644 applications/utilities/parallelProcessing/reconstructPar/lagrangianFieldReconstructorTemplates.C
 rename applications/utilities/parallelProcessing/reconstructPar/{pointFieldReconstructorReconstructFields.C => pointFieldReconstructorTemplates.C} (87%)
 delete mode 100644 applications/utilities/parallelProcessing/reconstructPar/reconstructLagrangian.H
 delete mode 100644 applications/utilities/parallelProcessing/reconstructPar/reconstructLagrangianFields.C
 delete mode 100644 applications/utilities/parallelProcessing/reconstructPar/reconstructLagrangianPositions.C
 create mode 100644 src/parallel/parallel/domainDecompositionNonConformal.C
 create mode 100644 src/parallel/parallel/multiDomainDecomposition.C
 create mode 100644 src/parallel/parallel/multiDomainDecomposition.H
 create mode 100644 src/parallel/parallel/multiRegionPrefixer.C
 create mode 100644 src/parallel/parallel/multiRegionPrefixer.H

diff --git a/applications/utilities/parallelProcessing/decomposePar/Make/files b/applications/utilities/parallelProcessing/decomposePar/Make/files
index e5f8107a96..cbc2c6bf7e 100644
--- a/applications/utilities/parallelProcessing/decomposePar/Make/files
+++ b/applications/utilities/parallelProcessing/decomposePar/Make/files
@@ -1,5 +1,4 @@
 decomposePar.C
-dimFieldDecomposer.C
 fvFieldDecomposer.C
 pointFieldDecomposer.C
 lagrangianFieldDecomposer.C
diff --git a/applications/utilities/parallelProcessing/decomposePar/decomposePar.C b/applications/utilities/parallelProcessing/decomposePar/decomposePar.C
index 608c07a8b6..19db79a679 100644
--- a/applications/utilities/parallelProcessing/decomposePar/decomposePar.C
+++ b/applications/utilities/parallelProcessing/decomposePar/decomposePar.C
@@ -50,11 +50,10 @@ Usage
         Copy any \a uniform directories too.
 
       - \par -constant
+        Decompose mesh and fields in the constant directory.
 
       - \par -time xxx:yyy \n
-        Override controlDict settings and decompose selected times. Does not
-        re-decompose the mesh i.e. does not handle moving mesh or changing
-        mesh cases.
+        Override controlDict settings and decompose selected times.
 
       - \par -fields \n
         Use existing geometry decomposition and convert fields only.
@@ -68,27 +67,12 @@ Usage
 
 \*---------------------------------------------------------------------------*/
 
-#include "processorRunTimes.H"
-#include "domainDecomposition.H"
-#include "decompositionMethod.H"
 #include "argList.H"
 #include "timeSelector.H"
-
-#include "labelIOField.H"
-#include "labelFieldIOField.H"
-#include "scalarIOField.H"
-#include "scalarFieldIOField.H"
-#include "vectorIOField.H"
-#include "vectorFieldIOField.H"
-#include "sphericalTensorIOField.H"
-#include "sphericalTensorFieldIOField.H"
-#include "symmTensorIOField.H"
-#include "symmTensorFieldIOField.H"
-#include "tensorIOField.H"
-#include "tensorFieldIOField.H"
-
-#include "readFields.H"
-#include "dimFieldDecomposer.H"
+#include "IOobjectList.H"
+#include "processorRunTimes.H"
+#include "multiDomainDecomposition.H"
+#include "decompositionMethod.H"
 #include "fvFieldDecomposer.H"
 #include "pointFieldDecomposer.H"
 #include "lagrangianFieldDecomposer.H"
@@ -100,40 +84,53 @@ using namespace Foam;
 namespace Foam
 {
 
+bool haveUniform
+(
+    const processorRunTimes& runTimes,
+    const word& regionDir = word::null
+)
+{
+    return
+        fileHandler().isDir
+        (
+            runTimes.completeTime().timePath()/regionDir/"uniform"
+        );
+}
+
+
 void decomposeUniform
 (
     const bool copyUniform,
-    const bool distributeUniform,
-    const Time& runTime,
-    const Time& procRunTime,
+    const processorRunTimes& runTimes,
     const word& regionDir = word::null
 )
 {
     const fileName uniformDir(regionDir/"uniform");
 
-    if (fileHandler().isDir(runTime.timePath()/uniformDir))
+    forAll(runTimes.procTimes(), proci)
     {
-        Info<< "Detected additional non-decomposed files in "
-            << runTime.timePath()/uniformDir
-            << endl;
+        const fileName procTimePath =
+            fileHandler().filePath(runTimes.procTimes()[proci].timePath());
 
-        const fileName timePath =
-            fileHandler().filePath(procRunTime.timePath());
-
-        if (copyUniform || distributeUniform)
+        if (!fileHandler().isDir(procTimePath))
         {
-            if (!fileHandler().exists(timePath/uniformDir))
+            fileHandler().mkDir(procTimePath);
+        }
+
+        if (copyUniform)
+        {
+            if (!fileHandler().exists(procTimePath/uniformDir))
             {
                 fileHandler().cp
                 (
-                    runTime.timePath()/uniformDir,
-                    timePath/uniformDir
+                    runTimes.completeTime().timePath()/uniformDir,
+                    procTimePath/uniformDir
                 );
             }
         }
         else
         {
-            // link with relative paths
+            // Link with relative paths
             string parentPath = string("..")/"..";
 
             if (regionDir != word::null)
@@ -142,16 +139,18 @@ void decomposeUniform
             }
 
             fileName currentDir(cwd());
-            chDir(timePath);
+
+            chDir(procTimePath);
 
             if (!fileHandler().exists(uniformDir))
             {
                 fileHandler().ln
                 (
-                    parentPath/runTime.name()/uniformDir,
+                    parentPath/runTimes.completeTime().name()/uniformDir,
                     uniformDir
                 );
             }
+
             chDir(currentDir);
         }
     }
@@ -179,14 +178,41 @@ void writeDecomposition(const domainDecomposition& meshes)
     cellProc.write();
 
     Info<< "Wrote decomposition as volScalarField::Internal to "
-        << cellProc.name() << " for use in postprocessing."
-        << nl << endl;
+        << cellProc.name() << " for use in postprocessing"
+        << endl;
 }
 
 
 }
 
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class delayedNewLine
+{
+    mutable bool first_;
+
+public:
+
+    delayedNewLine()
+    :
+        first_(true)
+    {}
+
+    friend Ostream& operator<<(Ostream& os, const delayedNewLine& dnl)
+    {
+        if (!dnl.first_) os << nl;
+        dnl.first_ = false;
+        return os;
+    }
+};
+
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
@@ -203,7 +229,7 @@ int main(int argc, char *argv[])
     (
         "cellProc",
         "write cell processor indices as a volScalarField::Internal for "
-        "post-processing."
+        "post-processing"
     );
     argList::addBoolOption
     (
@@ -236,19 +262,19 @@ int main(int argc, char *argv[])
         "remove existing processor*/ subdirs before decomposing the geometry"
     );
 
-    // Include explicit constant options, have zero from time range
+    // Include explicit constant option, execute from zero by default
     timeSelector::addOptions(true, false);
 
     #include "setRootCase.H"
 
-    bool region                  = args.optionFound("region");
-    bool writeCellProc           = args.optionFound("cellProc");
-    bool copyZero                = args.optionFound("copyZero");
-    bool copyUniform             = args.optionFound("copyUniform");
-    bool decomposeFieldsOnly     = args.optionFound("fields");
-    bool decomposeGeomOnly       = args.optionFound("noFields");
-    bool decomposeSets           = !args.optionFound("noSets");
-    bool forceOverwrite          = args.optionFound("force");
+    const bool region              = args.optionFound("region");
+    const bool writeCellProc       = args.optionFound("cellProc");
+    const bool copyZero            = args.optionFound("copyZero");
+    const bool copyUniform         = args.optionFound("copyUniform");
+    const bool decomposeFieldsOnly = args.optionFound("fields");
+    const bool decomposeGeomOnly   = args.optionFound("noFields");
+    const bool decomposeSets       = !args.optionFound("noSets");
+    const bool forceOverwrite      = args.optionFound("force");
 
     if (decomposeGeomOnly)
     {
@@ -264,14 +290,13 @@ int main(int argc, char *argv[])
     }
 
     // Set time from database
-    Info<< "Create time\n" << endl;
+    Info<< "Create time" << nl << endl;
     processorRunTimes runTimes(Foam::Time::controlDictName, args);
+    const Time& runTime = runTimes.completeTime();
 
     // Allow override of time
     const instantList times = runTimes.selectComplete(args);
 
-    const Time& runTime = runTimes.completeTime();
-
     #include "setRegionNames.H"
 
     // Remove existing processor directories if requested
@@ -288,7 +313,7 @@ int main(int argc, char *argv[])
             fileHandler().nProcs(runTimes.completeTime().path());
 
         Info<< "Removing " << nProcs0
-            << " existing processor directories" << endl;
+            << " existing processor directories" << nl << endl;
 
         // Remove existing processor directories
         const fileNameList dirs
@@ -356,17 +381,12 @@ int main(int argc, char *argv[])
     const dictionary decomposeParDict =
         decompositionMethod::decomposeParDict(runTimes.completeTime());
 
-    // Decompose all regions
+    // Check existing decomposition
     forAll(regionNames, regioni)
     {
         const word& regionName = regionNames[regioni];
-
-        const word& regionDir =
-            regionName == polyMesh::defaultRegion
-          ? word::null
-          : regionName;
-
-        Info<< "\n\nDecomposing mesh " << regionName << nl << endl;
+        const word regionDir =
+            regionName == polyMesh::defaultRegion ? word::null : regionName;
 
         // Determine the existing processor count directly
         const label nProcs =
@@ -388,690 +408,337 @@ int main(int argc, char *argv[])
                 << " domains as specified in decomposeParDict" << nl
                 << exit(FatalError);
         }
+    }
 
-        // Get flag to determine whether or not to distribute uniform data
-        const label distributeUniform =
-            decomposeParDict.lookupOrDefault<bool>("distributed", false);
+    // Create meshes
+    multiDomainDecomposition regionMeshes(runTimes, regionNames);
+    if
+    (
+        !(decomposeFieldsOnly && copyZero)
+     && regionMeshes.readDecompose(decomposeSets)
+    )
+    {
+        Info<< endl;
 
-        // Create meshes
-        Info<< "Create mesh" << endl;
-        domainDecomposition meshes(runTimes, regionName);
-        if (!decomposeFieldsOnly || !copyZero)
+        if (writeCellProc)
         {
-            if (meshes.readDecompose(decomposeSets) && writeCellProc)
+            forAll(regionNames, regioni)
             {
-                writeDecomposition(meshes);
-                fileHandler().flush();
-            }
-        }
-
-        // Field maps. These are preserved if decomposing multiple times.
-        PtrList<fvFieldDecomposer> fieldDecomposerList
-        (
-            meshes.nProcs()
-        );
-        PtrList<dimFieldDecomposer> dimFieldDecomposerList
-        (
-            meshes.nProcs()
-        );
-        PtrList<pointFieldDecomposer> pointFieldDecomposerList
-        (
-            meshes.nProcs()
-        );
-
-        // Loop over all times
-        forAll(times, timei)
-        {
-            // Set the time
-            runTimes.setTime(times[timei], timei);
-
-            Info<< "Time = " << runTimes.completeTime().userTimeName() << endl;
-
-            // Update the meshes, if necessary
-            fvMesh::readUpdateState state = fvMesh::UNCHANGED;
-            if (!decomposeFieldsOnly || !copyZero)
-            {
-                state = meshes.readUpdateDecompose();
-            }
-
-            // Write the mesh out, if necessary
-            if (decomposeFieldsOnly)
-            {
-                // Nothing to do
-            }
-            else if (state != fvMesh::UNCHANGED)
-            {
-                meshes.writeProcs(decomposeSets);
-            }
-
-            // Write the decomposition, if necessary
-            if
-            (
-                writeCellProc
-             && meshes.completeMesh().facesInstance()
-             == runTimes.completeTime().name()
-            )
-            {
-                writeDecomposition(meshes);
-                fileHandler().flush();
-            }
-
-            // Clear the field maps if there has been topology change
-            if (state >= fvMesh::TOPO_CHANGE)
-            {
-                for (label proci = 0; proci < meshes.nProcs(); proci++)
-                {
-                    fieldDecomposerList.set(proci, nullptr);
-                    dimFieldDecomposerList.set(proci, nullptr);
-                    pointFieldDecomposerList.set(proci, nullptr);
-                }
-            }
-
-            // Decompose the fields at this time
-            if (decomposeGeomOnly)
-            {
-                // Do nothing
-            }
-            else if (copyZero)
-            {
-                // Copy the field files from the <time> directory to the
-                // processor*/<time> directories without altering them
-                const fileName completeTimePath =
-                    runTimes.completeTime().timePath();
-
-                fileName prevProcTimePath;
-                for (label proci = 0; proci < runTimes.nProcs(); proci++)
-                {
-                    const Time& procRunTime = runTimes.procTimes()[proci];
-
-                    if (fileHandler().isDir(completeTimePath))
-                    {
-                        const fileName procTimePath
-                        (
-                            fileHandler().objectPath
-                            (
-                                IOobject
-                                (
-                                    "",
-                                    procRunTime.name(),
-                                    procRunTime
-                                ),
-                                word::null
-                            )
-                        );
-
-                        if (procTimePath != prevProcTimePath)
-                        {
-                            Info<< "Processor " << proci
-                                << ": copying " << completeTimePath << nl
-                                << " to " << procTimePath << endl;
-                            fileHandler().cp(completeTimePath, procTimePath);
-                            prevProcTimePath = procTimePath;
-                        }
-                    }
-                }
-            }
-            else
-            {
-                // Decompose the fields
-
-                // Search for list of objects for this time
-                IOobjectList objects
-                (
-                    meshes.completeMesh(),
-                    runTimes.completeTime().name()
-                );
-
-                // Construct the vol fields
-                PtrList<volScalarField> volScalarFields;
-                readFields(meshes.completeMesh(), objects, volScalarFields);
-                PtrList<volVectorField> volVectorFields;
-                readFields(meshes.completeMesh(), objects, volVectorFields);
-                PtrList<volSphericalTensorField> volSphericalTensorFields;
-                readFields
-                (
-                    meshes.completeMesh(),
-                    objects,
-                    volSphericalTensorFields
-                );
-                PtrList<volSymmTensorField> volSymmTensorFields;
-                readFields(meshes.completeMesh(), objects, volSymmTensorFields);
-                PtrList<volTensorField> volTensorFields;
-                readFields(meshes.completeMesh(), objects, volTensorFields);
-
-                // Construct the dimensioned fields
-                PtrList<DimensionedField<scalar, volMesh>> dimScalarFields;
-                readFields(meshes.completeMesh(), objects, dimScalarFields);
-                PtrList<DimensionedField<vector, volMesh>> dimVectorFields;
-                readFields(meshes.completeMesh(), objects, dimVectorFields);
-                PtrList<DimensionedField<sphericalTensor, volMesh>>
-                    dimSphericalTensorFields;
-                readFields
-                (
-                    meshes.completeMesh(),
-                    objects,
-                    dimSphericalTensorFields
-                );
-                PtrList<DimensionedField<symmTensor, volMesh>>
-                    dimSymmTensorFields;
-                readFields(meshes.completeMesh(), objects, dimSymmTensorFields);
-                PtrList<DimensionedField<tensor, volMesh>> dimTensorFields;
-                readFields(meshes.completeMesh(), objects, dimTensorFields);
-
-                // Construct the surface fields
-                PtrList<surfaceScalarField> surfaceScalarFields;
-                readFields(meshes.completeMesh(), objects, surfaceScalarFields);
-                PtrList<surfaceVectorField> surfaceVectorFields;
-                readFields(meshes.completeMesh(), objects, surfaceVectorFields);
-                PtrList<surfaceSphericalTensorField>
-                    surfaceSphericalTensorFields;
-                readFields
-                (
-                    meshes.completeMesh(),
-                    objects,
-                    surfaceSphericalTensorFields
-                );
-                PtrList<surfaceSymmTensorField> surfaceSymmTensorFields;
-                readFields
-                (
-                    meshes.completeMesh(),
-                    objects,
-                    surfaceSymmTensorFields
-                );
-                PtrList<surfaceTensorField> surfaceTensorFields;
-                readFields(meshes.completeMesh(), objects, surfaceTensorFields);
-
-                // Construct the point fields
-                const pointMesh& pMesh = pointMesh::New(meshes.completeMesh());
-                PtrList<pointScalarField> pointScalarFields;
-                readFields(pMesh, objects, pointScalarFields);
-                PtrList<pointVectorField> pointVectorFields;
-                readFields(pMesh, objects, pointVectorFields);
-                PtrList<pointSphericalTensorField> pointSphericalTensorFields;
-                readFields(pMesh, objects, pointSphericalTensorFields);
-                PtrList<pointSymmTensorField> pointSymmTensorFields;
-                readFields(pMesh, objects, pointSymmTensorFields);
-                PtrList<pointTensorField> pointTensorFields;
-                readFields(pMesh, objects, pointTensorFields);
-
-                // Construct the Lagrangian fields
-                fileNameList cloudDirs
-                (
-                    fileHandler().readDir
-                    (
-                        runTimes.completeTime().timePath()/cloud::prefix,
-                        fileType::directory
-                    )
-                );
-                PtrList<Cloud<indexedParticle>>
-                    lagrangianPositions(cloudDirs.size());
-                PtrList<List<SLList<indexedParticle*>*>>
-                    cellParticles(cloudDirs.size());
-                PtrList<PtrList<labelIOField>>
-                    lagrangianLabelFields(cloudDirs.size());
-                PtrList<PtrList<labelFieldCompactIOField>>
-                    lagrangianLabelFieldFields(cloudDirs.size());
-                PtrList<PtrList<scalarIOField>>
-                    lagrangianScalarFields(cloudDirs.size());
-                PtrList<PtrList<scalarFieldCompactIOField>>
-                    lagrangianScalarFieldFields(cloudDirs.size());
-                PtrList<PtrList<vectorIOField>>
-                    lagrangianVectorFields(cloudDirs.size());
-                PtrList<PtrList<vectorFieldCompactIOField>>
-                    lagrangianVectorFieldFields(cloudDirs.size());
-                PtrList<PtrList<sphericalTensorIOField>>
-                    lagrangianSphericalTensorFields(cloudDirs.size());
-                PtrList<PtrList<sphericalTensorFieldCompactIOField>>
-                    lagrangianSphericalTensorFieldFields(cloudDirs.size());
-                PtrList<PtrList<symmTensorIOField>>
-                    lagrangianSymmTensorFields(cloudDirs.size());
-                PtrList<PtrList<symmTensorFieldCompactIOField>>
-                    lagrangianSymmTensorFieldFields(cloudDirs.size());
-                PtrList<PtrList<tensorIOField>>
-                    lagrangianTensorFields(cloudDirs.size());
-                PtrList<PtrList<tensorFieldCompactIOField>>
-                    lagrangianTensorFieldFields(cloudDirs.size());
-
-                label cloudI = 0;
-
-                forAll(cloudDirs, i)
-                {
-                    IOobjectList sprayObjs
-                    (
-                        meshes.completeMesh(),
-                        runTimes.completeTime().name(),
-                        cloud::prefix/cloudDirs[i],
-                        IOobject::MUST_READ,
-                        IOobject::NO_WRITE,
-                        false
-                    );
-
-                    IOobject* positionsPtr = sprayObjs.lookup
-                    (
-                        word("positions")
-                    );
-
-                    if (positionsPtr)
-                    {
-                        // Read lagrangian particles
-                        Info<< "Identified lagrangian data set: "
-                            << cloudDirs[i] << endl;
-                        lagrangianPositions.set
-                        (
-                            cloudI,
-                            new Cloud<indexedParticle>
-                            (
-                                meshes.completeMesh(),
-                                cloudDirs[i],
-                                false
-                            )
-                        );
-
-                        // Sort particles per cell
-                        cellParticles.set
-                        (
-                            cloudI,
-                            new List<SLList<indexedParticle*>*>
-                            (
-                                meshes.completeMesh().nCells(),
-                                static_cast<SLList<indexedParticle*>*>(nullptr)
-                            )
-                        );
-
-                        // Populate the cloud
-                        label index = 0;
-                        forAllIter
-                        (
-                            Cloud<indexedParticle>,
-                            lagrangianPositions[cloudI],
-                            iter
-                        )
-                        {
-                            iter().index() = index ++;
-
-                            label celli = iter().cell();
-
-                            // Check
-                            if
-                            (
-                                celli < 0
-                             || celli >= meshes.completeMesh().nCells()
-                            )
-                            {
-                                FatalErrorInFunction
-                                    << "Illegal cell number " << celli
-                                    << " for particle with index "
-                                    << iter().index()
-                                    << " at position "
-                                    << iter().position(meshes.completeMesh())
-                                    << nl
-                                    << "Cell number should be between 0 and "
-                                    << meshes.completeMesh().nCells()-1 << nl
-                                    << "On this mesh the particle should"
-                                    << " be in cell "
-                                    << meshes.completeMesh().findCell
-                                       (
-                                           iter().position
-                                           (
-                                               meshes.completeMesh()
-                                           )
-                                       )
-                                    << exit(FatalError);
-                            }
-
-                            if (!cellParticles[cloudI][celli])
-                            {
-                                cellParticles[cloudI][celli] =
-                                    new SLList<indexedParticle*>();
-                            }
-
-                            cellParticles[cloudI][celli]->append(&iter());
-                        }
-
-                        // Read fields
-                        IOobjectList lagrangianObjects
-                        (
-                            meshes.completeMesh(),
-                            runTimes.completeTime().name(),
-                            cloud::prefix/cloudDirs[cloudI],
-                            IOobject::MUST_READ,
-                            IOobject::NO_WRITE,
-                            false
-                        );
-                        lagrangianFieldDecomposer::readFields
-                        (
-                            cloudI,
-                            lagrangianObjects,
-                            lagrangianLabelFields
-                        );
-                        lagrangianFieldDecomposer::readFieldFields
-                        (
-                            cloudI,
-                            lagrangianObjects,
-                            lagrangianLabelFieldFields
-                        );
-                        lagrangianFieldDecomposer::readFields
-                        (
-                            cloudI,
-                            lagrangianObjects,
-                            lagrangianScalarFields
-                        );
-                        lagrangianFieldDecomposer::readFieldFields
-                        (
-                            cloudI,
-                            lagrangianObjects,
-                            lagrangianScalarFieldFields
-                        );
-                        lagrangianFieldDecomposer::readFields
-                        (
-                            cloudI,
-                            lagrangianObjects,
-                            lagrangianVectorFields
-                        );
-                        lagrangianFieldDecomposer::readFieldFields
-                        (
-                            cloudI,
-                            lagrangianObjects,
-                            lagrangianVectorFieldFields
-                        );
-                        lagrangianFieldDecomposer::readFields
-                        (
-                            cloudI,
-                            lagrangianObjects,
-                            lagrangianSphericalTensorFields
-                        );
-                        lagrangianFieldDecomposer::readFieldFields
-                        (
-                            cloudI,
-                            lagrangianObjects,
-                            lagrangianSphericalTensorFieldFields
-                        );
-                        lagrangianFieldDecomposer::readFields
-                        (
-                            cloudI,
-                            lagrangianObjects,
-                            lagrangianSymmTensorFields
-                        );
-                        lagrangianFieldDecomposer::readFieldFields
-                        (
-                            cloudI,
-                            lagrangianObjects,
-                            lagrangianSymmTensorFieldFields
-                        );
-                        lagrangianFieldDecomposer::readFields
-                        (
-                            cloudI,
-                            lagrangianObjects,
-                            lagrangianTensorFields
-                        );
-                        lagrangianFieldDecomposer::readFieldFields
-                        (
-                            cloudI,
-                            lagrangianObjects,
-                            lagrangianTensorFieldFields
-                        );
-
-                        cloudI++;
-                    }
-                }
-
-                lagrangianPositions.setSize(cloudI);
-                cellParticles.setSize(cloudI);
-                lagrangianLabelFields.setSize(cloudI);
-                lagrangianLabelFieldFields.setSize(cloudI);
-                lagrangianScalarFields.setSize(cloudI);
-                lagrangianScalarFieldFields.setSize(cloudI);
-                lagrangianVectorFields.setSize(cloudI);
-                lagrangianVectorFieldFields.setSize(cloudI);
-                lagrangianSphericalTensorFields.setSize(cloudI);
-                lagrangianSphericalTensorFieldFields.setSize(cloudI);
-                lagrangianSymmTensorFields.setSize(cloudI);
-                lagrangianSymmTensorFieldFields.setSize(cloudI);
-                lagrangianTensorFields.setSize(cloudI);
-                lagrangianTensorFieldFields.setSize(cloudI);
-
+                writeDecomposition(regionMeshes.meshes(regioni)());
                 Info<< endl;
-
-                // split the fields over processors
-                for (label proci = 0; proci < meshes.nProcs(); proci++)
-                {
-                    Info<< "Processor " << proci << ": field transfer" << endl;
-
-                    // FV fields
-                    {
-                        if (!fieldDecomposerList.set(proci))
-                        {
-                            fieldDecomposerList.set
-                            (
-                                proci,
-                                new fvFieldDecomposer
-                                (
-                                    meshes.completeMesh(),
-                                    meshes.procMeshes()[proci],
-                                    meshes.procFaceAddressing()[proci],
-                                    meshes.procCellAddressing()[proci],
-                                    meshes.procFaceAddressingBf()[proci]
-                                )
-                            );
-                        }
-                        const fvFieldDecomposer& fieldDecomposer =
-                            fieldDecomposerList[proci];
-
-                        fieldDecomposer.decomposeFields(volScalarFields);
-                        fieldDecomposer.decomposeFields(volVectorFields);
-                        fieldDecomposer.decomposeFields
-                        (
-                            volSphericalTensorFields
-                        );
-                        fieldDecomposer.decomposeFields(volSymmTensorFields);
-                        fieldDecomposer.decomposeFields(volTensorFields);
-
-                        fieldDecomposer.decomposeFields(surfaceScalarFields);
-                        fieldDecomposer.decomposeFields(surfaceVectorFields);
-                        fieldDecomposer.decomposeFields
-                        (
-                            surfaceSphericalTensorFields
-                        );
-                        fieldDecomposer.decomposeFields
-                        (
-                            surfaceSymmTensorFields
-                        );
-                        fieldDecomposer.decomposeFields(surfaceTensorFields);
-
-                        if (times.size() == 1)
-                        {
-                            // Clear cached decomposer
-                            fieldDecomposerList.set(proci, nullptr);
-                        }
-                    }
-
-                    // Dimensioned fields
-                    {
-                        if (!dimFieldDecomposerList.set(proci))
-                        {
-                            dimFieldDecomposerList.set
-                            (
-                                proci,
-                                new dimFieldDecomposer
-                                (
-                                    meshes.completeMesh(),
-                                    meshes.procMeshes()[proci],
-                                    meshes.procFaceAddressing()[proci],
-                                    meshes.procCellAddressing()[proci]
-                                )
-                            );
-                        }
-                        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
-                    (
-                        pointScalarFields.size()
-                     || pointVectorFields.size()
-                     || pointSphericalTensorFields.size()
-                     || pointSymmTensorFields.size()
-                     || pointTensorFields.size()
-                    )
-                    {
-                        const pointMesh& procPMesh =
-                            pointMesh::New(meshes.procMeshes()[proci]);
-
-                        if (!pointFieldDecomposerList.set(proci))
-                        {
-                            pointFieldDecomposerList.set
-                            (
-                                proci,
-                                new pointFieldDecomposer
-                                (
-                                    pMesh,
-                                    procPMesh,
-                                    meshes.procPointAddressing()[proci]
-                                )
-                            );
-                        }
-                        const pointFieldDecomposer& pointDecomposer =
-                            pointFieldDecomposerList[proci];
-
-                        pointDecomposer.decomposeFields(pointScalarFields);
-                        pointDecomposer.decomposeFields(pointVectorFields);
-                        pointDecomposer.decomposeFields
-                        (
-                            pointSphericalTensorFields
-                        );
-                        pointDecomposer.decomposeFields(pointSymmTensorFields);
-                        pointDecomposer.decomposeFields(pointTensorFields);
-
-                        if (times.size() == 1)
-                        {
-                            pointFieldDecomposerList.set(proci, nullptr);
-                        }
-                    }
-
-                    // If there is lagrangian data write it out
-                    forAll(lagrangianPositions, cloudI)
-                    {
-                        if (lagrangianPositions[cloudI].size())
-                        {
-                            lagrangianFieldDecomposer fieldDecomposer
-                            (
-                                meshes.completeMesh(),
-                                meshes.procMeshes()[proci],
-                                meshes.procFaceAddressing()[proci],
-                                meshes.procCellAddressing()[proci],
-                                cloudDirs[cloudI],
-                                lagrangianPositions[cloudI],
-                                cellParticles[cloudI]
-                            );
-
-                            // Lagrangian fields
-                            {
-                                fieldDecomposer.decomposeFields
-                                (
-                                    cloudDirs[cloudI],
-                                    lagrangianLabelFields[cloudI]
-                                );
-                                fieldDecomposer.decomposeFieldFields
-                                (
-                                    cloudDirs[cloudI],
-                                    lagrangianLabelFieldFields[cloudI]
-                                );
-                                fieldDecomposer.decomposeFields
-                                (
-                                    cloudDirs[cloudI],
-                                    lagrangianScalarFields[cloudI]
-                                );
-                                fieldDecomposer.decomposeFieldFields
-                                (
-                                    cloudDirs[cloudI],
-                                    lagrangianScalarFieldFields[cloudI]
-                                );
-                                fieldDecomposer.decomposeFields
-                                (
-                                    cloudDirs[cloudI],
-                                    lagrangianVectorFields[cloudI]
-                                );
-                                fieldDecomposer.decomposeFieldFields
-                                (
-                                    cloudDirs[cloudI],
-                                    lagrangianVectorFieldFields[cloudI]
-                                );
-                                fieldDecomposer.decomposeFields
-                                (
-                                    cloudDirs[cloudI],
-                                    lagrangianSphericalTensorFields[cloudI]
-                                );
-                                fieldDecomposer.decomposeFieldFields
-                                (
-                                    cloudDirs[cloudI],
-                                    lagrangianSphericalTensorFieldFields[cloudI]
-                                );
-                                fieldDecomposer.decomposeFields
-                                (
-                                    cloudDirs[cloudI],
-                                    lagrangianSymmTensorFields[cloudI]
-                                );
-                                fieldDecomposer.decomposeFieldFields
-                                (
-                                    cloudDirs[cloudI],
-                                    lagrangianSymmTensorFieldFields[cloudI]
-                                );
-                                fieldDecomposer.decomposeFields
-                                (
-                                    cloudDirs[cloudI],
-                                    lagrangianTensorFields[cloudI]
-                                );
-                                fieldDecomposer.decomposeFieldFields
-                                (
-                                    cloudDirs[cloudI],
-                                    lagrangianTensorFieldFields[cloudI]
-                                );
-                            }
-                        }
-                    }
-
-                    // Decompose the "uniform" directory in the region time
-                    // directory
-                    decomposeUniform
-                    (
-                        copyUniform,
-                        distributeUniform,
-                        runTimes.completeTime(),
-                        meshes.procMeshes()[proci].time(),
-                        regionDir
-                    );
-
-                    // For the first region of a multi-region case additionally
-                    // decompose the "uniform" directory in the no-region time
-                    // directory
-                    if (regionNames.size() > 1 && regioni == 0)
-                    {
-                        decomposeUniform
-                        (
-                            copyUniform,
-                            distributeUniform,
-                            runTimes.completeTime(),
-                            meshes.procMeshes()[proci].time()
-                        );
-                    }
-                }
+                fileHandler().flush();
             }
         }
     }
 
-    Info<< "\nEnd\n" << endl;
+    // Get flag to determine whether or not to distribute uniform data
+    const bool distributed =
+        decomposeParDict.lookupOrDefault<bool>("distributed", false);
+
+    // Loop over all times
+    forAll(times, timei)
+    {
+        // Set the time
+        runTimes.setTime(times[timei], timei);
+
+        Info<< "Time = " << runTimes.completeTime().userTimeName()
+            << nl << endl;
+
+        // Update the meshes, if necessary
+        const fvMesh::readUpdateState stat =
+            !(decomposeFieldsOnly && copyZero)
+          ? regionMeshes.readUpdateDecompose()
+          : fvMesh::UNCHANGED;
+        if (stat >= fvMesh::TOPO_CHANGE) Info<< endl;
+
+        // Write the mesh out (if anything has changed), if necessary
+        if (!decomposeFieldsOnly)
+        {
+            regionMeshes.writeProcs(decomposeSets);
+        }
+
+        // Write the decomposition, if necessary
+        forAll(regionNames, regioni)
+        {
+            if (writeCellProc && stat >= fvMesh::TOPO_CHANGE)
+            {
+                writeDecomposition(regionMeshes.meshes(regioni)());
+                Info<< endl;
+                fileHandler().flush();
+            }
+        }
+
+        // If only decomposing geometry then there is no more to do
+        if (decomposeGeomOnly)
+        {
+            continue;
+        }
+
+        // If copying from  zero then just copy everything from the <time>
+        // directory to the processor*/<time> directories without altering them
+        if (copyZero)
+        {
+            const fileName completeTimePath =
+                runTimes.completeTime().timePath();
+
+            fileName prevProcTimePath;
+            for (label proci = 0; proci < runTimes.nProcs(); proci++)
+            {
+                const Time& procRunTime = runTimes.procTimes()[proci];
+
+                if (fileHandler().isDir(completeTimePath))
+                {
+                    const fileName procTimePath
+                    (
+                        fileHandler().objectPath
+                        (
+                            IOobject
+                            (
+                                "",
+                                procRunTime.name(),
+                                procRunTime
+                            ),
+                            word::null
+                        )
+                    );
+
+                    if (procTimePath != prevProcTimePath)
+                    {
+                        Info<< "Processor " << proci
+                            << ": copying " << completeTimePath << nl
+                            << " to " << procTimePath << endl;
+                        fileHandler().cp(completeTimePath, procTimePath);
+                        prevProcTimePath = procTimePath;
+                    }
+                }
+            }
+
+            Info<< endl;
+
+            continue;
+        }
+
+        // Otherwise we are doing full region-by-region decomposition of all
+        // the available fields
+        forAll(regionNames, regioni)
+        {
+            const word& regionName = regionNames[regioni];
+            const word regionDir =
+                regionName == polyMesh::defaultRegion ? word::null : regionName;
+
+            const delayedNewLine dnl;
+
+            // Prefixed scope
+            {
+                const RegionConstRef<domainDecomposition> meshes =
+                    regionMeshes.meshes(regioni);
+
+                // Search for objects at this time
+                IOobjectList objects
+                (
+                    meshes().completeMesh(),
+                    runTimes.completeTime().name()
+                );
+
+                {
+                    Info<< dnl << "Decomposing FV fields" << endl;
+
+                    if (fvFieldDecomposer::decomposes(objects))
+                    {
+                        fvFieldDecomposer fvDecomposer
+                        (
+                            meshes().completeMesh(),
+                            meshes().procMeshes(),
+                            meshes().procFaceAddressing(),
+                            meshes().procCellAddressing(),
+                            meshes().procFaceAddressingBf()
+                        );
+
+                        #define DO_FV_VOL_INTERNAL_FIELDS_TYPE(Type, nullArg)  \
+                            fvDecomposer.decomposeVolInternalFields<Type>      \
+                            (objects);
+                        FOR_ALL_FIELD_TYPES(DO_FV_VOL_INTERNAL_FIELDS_TYPE)
+                        #undef DO_FV_VOL_INTERNAL_FIELDS_TYPE
+
+                        #define DO_FV_VOL_FIELDS_TYPE(Type, nullArg)           \
+                            fvDecomposer.decomposeVolFields<Type>              \
+                            (objects);
+                        FOR_ALL_FIELD_TYPES(DO_FV_VOL_FIELDS_TYPE)
+                        #undef DO_FV_VOL_FIELDS_TYPE
+
+                        #define DO_FV_SURFACE_FIELDS_TYPE(Type, nullArg)       \
+                            fvDecomposer.decomposeFvSurfaceFields<Type>        \
+                            (objects);
+                        FOR_ALL_FIELD_TYPES(DO_FV_SURFACE_FIELDS_TYPE)
+                        #undef DO_FV_SURFACE_FIELDS_TYPE
+                    }
+                    else
+                    {
+                        Info<< dnl << "    (no FV fields)" << endl;
+                    }
+                }
+
+                {
+                    Info<< dnl << "Decomposing point fields" << endl;
+
+                    if (pointFieldDecomposer::decomposes(objects))
+                    {
+                        pointFieldDecomposer pointDecomposer
+                        (
+                            pointMesh::New(meshes().completeMesh()),
+                            meshes().procMeshes(),
+                            meshes().procPointAddressing()
+                        );
+
+                        #define DO_POINT_FIELDS_TYPE(Type, nullArg)            \
+                            pointDecomposer.decomposeFields<Type>              \
+                            (objects);
+                        FOR_ALL_FIELD_TYPES(DO_POINT_FIELDS_TYPE)
+                        #undef DO_POINT_FIELDS_TYPE
+                    }
+                    else
+                    {
+                        Info<< dnl << "    (no point fields)" << endl;
+                    }
+                }
+
+                {
+                    // Find cloud directories
+                    fileNameList cloudDirs
+                    (
+                        fileHandler().readDir
+                        (
+                            runTimes.completeTime().timePath()
+                           /regionDir
+                           /cloud::prefix,
+                            fileType::directory
+                        )
+                    );
+
+                    // Add objects in any found cloud directories
+                    HashTable<IOobjectList> cloudsObjects;
+                    forAll(cloudDirs, i)
+                    {
+                        // Do local scan for valid cloud objects
+                        IOobjectList cloudObjs
+                        (
+                            meshes().completeMesh(),
+                            runTimes.completeTime().name(),
+                            cloud::prefix/cloudDirs[i],
+                            IOobject::MUST_READ,
+                            IOobject::NO_WRITE,
+                            false
+                        );
+
+                        // If "positions" is present, then add to the table
+                        if (cloudObjs.lookup(word("positions")))
+                        {
+                            cloudsObjects.insert(cloudDirs[i], cloudObjs);
+                        }
+                    }
+
+                    // Decompose the objects found above
+                    if (cloudsObjects.size())
+                    {
+                        forAllConstIter
+                        (
+                            HashTable<IOobjectList>,
+                            cloudsObjects,
+                            iter
+                        )
+                        {
+                            const word cloudName =
+                                string::validate<word>(iter.key());
+
+                            const IOobjectList& cloudObjects = iter();
+
+                            Info<< dnl << "Decomposing lagrangian fields for "
+                                << "cloud " << cloudName << endl;
+
+                            if
+                            (
+                                lagrangianFieldDecomposer::decomposes
+                                (
+                                    cloudObjects
+                                )
+                            )
+                            {
+                                lagrangianFieldDecomposer
+                                    lagrangianDecomposer
+                                    (
+                                        meshes().completeMesh(),
+                                        meshes().procMeshes(),
+                                        meshes().procFaceAddressing(),
+                                        meshes().procCellAddressing(),
+                                        cloudName
+                                    );
+
+                                #define DO_CLOUD_FIELDS_TYPE(Type, nullArg)    \
+                                    lagrangianDecomposer.decomposeFields<Type> \
+                                    (cloudObjects);
+                                DO_CLOUD_FIELDS_TYPE(label, );
+                                FOR_ALL_FIELD_TYPES(DO_CLOUD_FIELDS_TYPE)
+                                #undef DO_CLOUD_FIELDS_TYPE
+                            }
+                            else
+                            {
+                                Info<< dnl << "    (no lagrangian fields)"
+                                    << endl;
+                            }
+                        }
+                    }
+                }
+            }
+
+            Info<< dnl;
+        }
+
+        // Distribute the uniform directory
+        if (haveUniform(runTimes))
+        {
+            Info<< "Distributing uniform files" << endl;
+
+            decomposeUniform
+            (
+                copyUniform || distributed,
+                runTimes
+            );
+
+            Info<< endl;
+        }
+
+        if (regionNames == wordList(1, polyMesh::defaultRegion)) continue;
+
+        // Distribute the region uniform directories
+        forAll(regionNames, regioni)
+        {
+            const word& regionName = regionNames[regioni];
+            const word regionDir =
+                regionName == polyMesh::defaultRegion ? word::null : regionName;
+
+            if (haveUniform(runTimes, regionDir))
+            {
+                // Prefixed scope
+                {
+                    const RegionConstRef<domainDecomposition> meshes =
+                        regionMeshes.meshes(regioni);
+
+                    Info<< "Distributing uniform files" << endl;
+
+                    decomposeUniform
+                    (
+                        copyUniform || distributed,
+                        runTimes,
+                        regionDir
+                    );
+                }
+
+                Info<< endl;
+            }
+        }
+    }
+
+    Info<< "End" << nl << endl;
 
     return 0;
 }
diff --git a/applications/utilities/parallelProcessing/decomposePar/dimFieldDecomposer.C b/applications/utilities/parallelProcessing/decomposePar/dimFieldDecomposer.C
deleted file mode 100644
index 63f7808c12..0000000000
--- a/applications/utilities/parallelProcessing/decomposePar/dimFieldDecomposer.C
+++ /dev/null
@@ -1,49 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "dimFieldDecomposer.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::dimFieldDecomposer::dimFieldDecomposer
-(
-    const fvMesh& completeMesh,
-    const fvMesh& procMesh,
-    const labelList& faceAddressing,
-    const labelList& cellAddressing
-)
-:
-    procMesh_(procMesh),
-    cellAddressing_(cellAddressing)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::dimFieldDecomposer::~dimFieldDecomposer()
-{}
-
-
-// ************************************************************************* //
diff --git a/applications/utilities/parallelProcessing/decomposePar/dimFieldDecomposer.H b/applications/utilities/parallelProcessing/decomposePar/dimFieldDecomposer.H
deleted file mode 100644
index 632eb6e87e..0000000000
--- a/applications/utilities/parallelProcessing/decomposePar/dimFieldDecomposer.H
+++ /dev/null
@@ -1,121 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::dimFieldDecomposer
-
-Description
-    Dimensioned field decomposer.
-
-SourceFiles
-    dimFieldDecomposer.C
-    dimFieldDecomposerDecomposeFields.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef dimFieldDecomposer_H
-#define dimFieldDecomposer_H
-
-#include "fvMesh.H"
-#include "surfaceFields.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-class IOobjectList;
-
-/*---------------------------------------------------------------------------*\
-                      Class fvFieldDecomposer Declaration
-\*---------------------------------------------------------------------------*/
-
-class dimFieldDecomposer
-{
-    // Private Data
-
-        //- Reference to processor mesh
-        const fvMesh& procMesh_;
-
-        //- Reference to cell addressing
-        const labelList& cellAddressing_;
-
-
-public:
-
-    // Constructors
-
-        //- Construct from components
-        dimFieldDecomposer
-        (
-            const fvMesh& completeMesh,
-            const fvMesh& procMesh,
-            const labelList& faceAddressing,
-            const labelList& cellAddressing
-        );
-
-        //- Disallow default bitwise copy construction
-        dimFieldDecomposer(const dimFieldDecomposer&) = delete;
-
-
-    //- Destructor
-    ~dimFieldDecomposer();
-
-
-    // Member Functions
-
-        //- Decompose field
-        template<class Type>
-        tmp<DimensionedField<Type, volMesh>> decomposeField
-        (
-            const DimensionedField<Type, volMesh>& field
-        ) const;
-
-
-        //- Decompose llist of fields
-        template<class GeoField>
-        void decomposeFields(const PtrList<GeoField>& fields) const;
-
-
-    // Member Operators
-
-        //- Disallow default bitwise assignment
-        void operator=(const dimFieldDecomposer&) = delete;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-    #include "dimFieldDecomposerDecomposeFields.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/applications/utilities/parallelProcessing/decomposePar/dimFieldDecomposerDecomposeFields.C b/applications/utilities/parallelProcessing/decomposePar/dimFieldDecomposerDecomposeFields.C
deleted file mode 100644
index 66b82b863d..0000000000
--- a/applications/utilities/parallelProcessing/decomposePar/dimFieldDecomposerDecomposeFields.C
+++ /dev/null
@@ -1,75 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "dimFieldDecomposer.H"
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class Type>
-Foam::tmp<Foam::DimensionedField<Type, Foam::volMesh>>
-Foam::dimFieldDecomposer::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 DimensionedField<Type, volMesh>
-        (
-            IOobject
-            (
-                field.name(),
-                procMesh_.time().name(),
-                procMesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            procMesh_,
-            field.dimensions(),
-            mappedField
-        )
-    );
-}
-
-
-template<class GeoField>
-void Foam::dimFieldDecomposer::decomposeFields
-(
-    const PtrList<GeoField>& fields
-) const
-{
-    forAll(fields, fieldi)
-    {
-        decomposeField(fields[fieldi])().write();
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/applications/utilities/parallelProcessing/decomposePar/fvFieldDecomposer.C b/applications/utilities/parallelProcessing/decomposePar/fvFieldDecomposer.C
index 61005c2fe7..0fcf49dfeb 100644
--- a/applications/utilities/parallelProcessing/decomposePar/fvFieldDecomposer.C
+++ b/applications/utilities/parallelProcessing/decomposePar/fvFieldDecomposer.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,10 +29,11 @@ License
 
 Foam::label Foam::fvFieldDecomposer::completePatchID
 (
+    const label proci,
     const label procPatchi
 ) const
 {
-    const fvPatch& procPatch = procMesh_.boundary()[procPatchi];
+    const fvPatch& procPatch = procMeshes_[proci].boundary()[procPatchi];
 
     if (procPatchi < completeMesh_.boundary().size())
     {
@@ -64,34 +65,45 @@ Foam::fvFieldDecomposer::patchFieldDecomposer::patchFieldDecomposer
 Foam::fvFieldDecomposer::fvFieldDecomposer
 (
     const fvMesh& completeMesh,
-    const fvMesh& procMesh,
-    const labelList& faceAddressing,
-    const labelList& cellAddressing,
-    const surfaceLabelField::Boundary& faceAddressingBf
+    const PtrList<fvMesh>& procMeshes,
+    const labelListList& faceProcAddressing,
+    const labelListList& cellProcAddressing,
+    const PtrList<surfaceLabelField::Boundary>& faceProcAddressingBf
 )
 :
     completeMesh_(completeMesh),
-    procMesh_(procMesh),
-    faceAddressing_(faceAddressing),
-    cellAddressing_(cellAddressing),
-    faceAddressingBf_(faceAddressingBf),
-    patchFieldDecomposers_(procMesh_.boundary().size())
+    procMeshes_(procMeshes),
+    faceProcAddressing_(faceProcAddressing),
+    cellProcAddressing_(cellProcAddressing),
+    faceProcAddressingBf_(faceProcAddressingBf),
+    patchFieldDecomposers_(procMeshes_.size())
 {
-    forAll(procMesh_.boundary(), procPatchi)
+    forAll(procMeshes_, proci)
     {
-        const label completePatchi = completePatchID(procPatchi);
-
-        // If there is a corresponding complete patch then create a patch mapper
-        if (completePatchi >= 0)
-        {
-            patchFieldDecomposers_.set
+        patchFieldDecomposers_.set
+        (
+            proci,
+            new PtrList<patchFieldDecomposer>
             (
-                procPatchi,
-                new patchFieldDecomposer
+                procMeshes_[proci].boundary().size()
+            )
+        );
+
+        forAll(procMeshes_[proci].boundary(), procPatchi)
+        {
+            const label completePatchi = completePatchID(proci, procPatchi);
+
+            if (completePatchi >= 0)
+            {
+                patchFieldDecomposers_[proci].set
                 (
-                    faceAddressingBf[completePatchi]
-                )
-            );
+                    procPatchi,
+                    new patchFieldDecomposer
+                    (
+                        faceProcAddressingBf[proci][completePatchi]
+                    )
+                );
+            }
         }
     }
 }
@@ -103,4 +115,22 @@ Foam::fvFieldDecomposer::~fvFieldDecomposer()
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::fvFieldDecomposer::decomposes(const IOobjectList& objects)
+{
+    bool result = false;
+
+    #define DO_FV_FIELDS_TYPE(Type, nullArg)                                   \
+        result = result                                                        \
+         || !objects.lookupClass(VolField<Type>::Internal::typeName).empty()   \
+         || !objects.lookupClass(VolField<Type>::typeName).empty()             \
+         || !objects.lookupClass(SurfaceField<Type>::typeName).empty();
+    FOR_ALL_FIELD_TYPES(DO_FV_FIELDS_TYPE)
+    #undef DO_FV_FIELDS_TYPE
+
+    return result;
+}
+
+
 // ************************************************************************* //
diff --git a/applications/utilities/parallelProcessing/decomposePar/fvFieldDecomposer.H b/applications/utilities/parallelProcessing/decomposePar/fvFieldDecomposer.H
index fe9dc25f39..3a4f050589 100644
--- a/applications/utilities/parallelProcessing/decomposePar/fvFieldDecomposer.H
+++ b/applications/utilities/parallelProcessing/decomposePar/fvFieldDecomposer.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,7 +29,7 @@ Description
 
 SourceFiles
     fvFieldDecomposer.C
-    fvFieldDecomposerDecomposeFields.C
+    fvFieldDecomposerTemplates.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -37,8 +37,10 @@ SourceFiles
 #define fvFieldDecomposer_H
 
 #include "fvMesh.H"
-#include "directFvPatchFieldMapper.H"
+#include "IOobjectList.H"
+#include "volFields.H"
 #include "surfaceFields.H"
+#include "directFvPatchFieldMapper.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -79,26 +81,26 @@ private:
         //- Reference to complete mesh
         const fvMesh& completeMesh_;
 
-        //- Reference to processor mesh
-        const fvMesh& procMesh_;
+        //- List of processor meshes
+        const PtrList<fvMesh>& procMeshes_;
 
         //- Reference to face addressing
-        const labelList& faceAddressing_;
+        const labelListList& faceProcAddressing_;
 
         //- Reference to cell addressing
-        const labelList& cellAddressing_;
+        const labelListList& cellProcAddressing_;
 
         //- Reference to face addressing boundary field
-        const surfaceLabelField::Boundary& faceAddressingBf_;
+        const PtrList<surfaceLabelField::Boundary>& faceProcAddressingBf_;
 
         //- List of patch field decomposers
-        PtrList<patchFieldDecomposer> patchFieldDecomposers_;
+        PtrList<PtrList<patchFieldDecomposer>> patchFieldDecomposers_;
 
 
     // Private Member Functions
 
         //- Convert a processor patch to the corresponding complete patch index
-        label completePatchID(const label procPatchi) const;
+        label completePatchID(const label proci, const label procPatchi) const;
 
         //- Map cell values to faces
         template<class Type>
@@ -119,6 +121,21 @@ private:
             const bool isFlux
         );
 
+        //- Decompose a volume internal field
+        template<class Type>
+        PtrList<typename VolField<Type>::Internal>
+        decomposeVolInternalField(const IOobject& fieldIoObject) const;
+
+        //- Decompose a volume field
+        template<class Type>
+        PtrList<VolField<Type>>
+        decomposeVolField(const IOobject& fieldIoObject) const;
+
+        //- Decompose a surface field
+        template<class Type>
+        PtrList<SurfaceField<Type>>
+        decomposeFvSurfaceField(const IOobject& fieldIoObject) const;
+
 
 public:
 
@@ -128,10 +145,10 @@ public:
         fvFieldDecomposer
         (
             const fvMesh& completeMesh,
-            const fvMesh& procMesh,
-            const labelList& faceAddressing,
-            const labelList& cellAddressing,
-            const surfaceLabelField::Boundary& faceAddressingBf
+            const PtrList<fvMesh>& procMeshes,
+            const labelListList& faceProcAddressing,
+            const labelListList& cellProcAddressing,
+            const PtrList<surfaceLabelField::Boundary>& faceProcAddressingBf
         );
 
         //- Disallow default bitwise copy construction
@@ -144,25 +161,20 @@ public:
 
     // Member Functions
 
-        //- Decompose volume field
-        template<class Type>
-        tmp<VolField<Type>>
-        decomposeField
-        (
-            const VolField<Type>& field
-        ) const;
+        //- Return whether anything in the object list gets decomposed
+        static bool decomposes(const IOobjectList& objects);
 
-        //- Decompose surface field
+        //- Read, decompose and write all volume internal fields
         template<class Type>
-        tmp<SurfaceField<Type>>
-        decomposeField
-        (
-            const SurfaceField<Type>& field
-        ) const;
+        void decomposeVolInternalFields(const IOobjectList& objects);
 
-        //- Decompose a list of fields
-        template<class GeoField>
-        void decomposeFields(const PtrList<GeoField>& fields) const;
+        //- Read, decompose and write all volume fields
+        template<class Type>
+        void decomposeVolFields(const IOobjectList& objects);
+
+        //- Read, decompose and write all surface fields
+        template<class Type>
+        void decomposeFvSurfaceFields(const IOobjectList& objects);
 
 
     // Member Operators
@@ -179,7 +191,7 @@ public:
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-    #include "fvFieldDecomposerDecomposeFields.C"
+    #include "fvFieldDecomposerTemplates.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/utilities/parallelProcessing/decomposePar/fvFieldDecomposerDecomposeFields.C b/applications/utilities/parallelProcessing/decomposePar/fvFieldDecomposerDecomposeFields.C
deleted file mode 100644
index 963ead583e..0000000000
--- a/applications/utilities/parallelProcessing/decomposePar/fvFieldDecomposerDecomposeFields.C
+++ /dev/null
@@ -1,364 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvFieldDecomposer.H"
-#include "processorFvPatchField.H"
-#include "processorFvsPatchField.H"
-#include "processorCyclicFvPatchField.H"
-#include "processorCyclicFvsPatchField.H"
-#include "emptyFvPatchFields.H"
-#include "stringOps.H"
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class Type>
-Foam::tmp<Foam::Field<Type>> Foam::fvFieldDecomposer::mapCellToFace
-(
-    const labelUList& owner,
-    const labelUList& neighbour,
-    const Field<Type>& field,
-    const labelUList& addressing
-)
-{
-    tmp<Field<Type>> tfld(new Field<Type>(addressing.size()));
-    Field<Type>& fld = tfld.ref();
-
-    forAll(addressing, i)
-    {
-        fld[i] =
-            field
-            [
-                addressing[i] > 0
-              ? neighbour[addressing[i] - 1]
-              : owner[- addressing[i] - 1]
-            ];
-    }
-
-    return tfld;
-}
-
-
-template<class Type>
-Foam::tmp<Foam::Field<Type>> Foam::fvFieldDecomposer::mapFaceToFace
-(
-    const Field<Type>& field,
-    const labelUList& addressing,
-    const bool isFlux
-)
-{
-    tmp<Field<Type>> tfld(new Field<Type>(addressing.size()));
-    Field<Type>& fld = tfld.ref();
-
-    forAll(addressing, i)
-    {
-        fld[i] =
-            (isFlux && addressing[i] < 0 ? -1 : +1)
-           *field[mag(addressing[i]) - 1];
-    }
-
-    return tfld;
-}
-
-
-template<class Type>
-Foam::tmp<Foam::VolField<Type>>
-Foam::fvFieldDecomposer::decomposeField
-(
-    const VolField<Type>& field
-) const
-{
-    // Create dummy patch fields
-    PtrList<fvPatchField<Type>> patchFields(procMesh_.boundary().size());
-    forAll(procMesh_.boundary(), procPatchi)
-    {
-        patchFields.set
-        (
-            procPatchi,
-            fvPatchField<Type>::New
-            (
-                calculatedFvPatchField<Type>::typeName,
-                procMesh_.boundary()[procPatchi],
-                DimensionedField<Type, volMesh>::null()
-            )
-        );
-    }
-
-    // Create the processor field with the dummy patch fields
-    tmp<VolField<Type>> tresF
-    (
-        new VolField<Type>
-        (
-            IOobject
-            (
-                field.name(),
-                procMesh_.time().name(),
-                procMesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            procMesh_,
-            field.dimensions(),
-            Field<Type>(field.primitiveField(), cellAddressing_),
-            patchFields
-        )
-    );
-    VolField<Type>& resF = tresF.ref();
-
-    // Change the patch fields to the correct type using a mapper constructor
-    // (with reference to the now correct internal field)
-    typename VolField<Type>::
-        Boundary& bf = resF.boundaryFieldRef();
-    forAll(bf, procPatchi)
-    {
-        const fvPatch& procPatch = procMesh_.boundary()[procPatchi];
-
-        const label completePatchi = completePatchID(procPatchi);
-
-        if (completePatchi == procPatchi)
-        {
-            bf.set
-            (
-                procPatchi,
-                fvPatchField<Type>::New
-                (
-                    field.boundaryField()[completePatchi],
-                    procPatch,
-                    resF(),
-                    patchFieldDecomposers_[procPatchi]
-                )
-            );
-        }
-        else if (isA<processorCyclicFvPatch>(procPatch))
-        {
-            if (field.boundaryField()[completePatchi].overridesConstraint())
-            {
-                OStringStream str;
-                str << "\nThe field \"" << field.name()
-                    << "\" on cyclic patch \""
-                    << field.boundaryField()[completePatchi].patch().name()
-                    << "\" cannot be decomposed as it is not a cyclic "
-                    << "patch field. A \"patchType cyclic;\" setting has "
-                    << "been used to override the cyclic patch type.\n\n"
-                    << "Cyclic patches like this with non-cyclic boundary "
-                    << "conditions should be confined to a single "
-                    << "processor using decomposition constraints.";
-                FatalErrorInFunction
-                    << stringOps::breakIntoIndentedLines(str.str()).c_str()
-                    << exit(FatalError);
-            }
-
-            const label nbrCompletePatchi =
-                refCast<const processorCyclicFvPatch>(procPatch)
-               .referPatch().nbrPatchID();
-
-            bf.set
-            (
-                procPatchi,
-                new processorCyclicFvPatchField<Type>
-                (
-                    procPatch,
-                    resF(),
-                    mapCellToFace
-                    (
-                        labelUList(),
-                        completeMesh_.lduAddr().patchAddr(nbrCompletePatchi),
-                        field.primitiveField(),
-                        faceAddressingBf_[procPatchi]
-                    )
-                )
-            );
-        }
-        else if (isA<processorFvPatch>(procPatch))
-        {
-            bf.set
-            (
-                procPatchi,
-                new processorFvPatchField<Type>
-                (
-                    procPatch,
-                    resF(),
-                    mapCellToFace
-                    (
-                        completeMesh_.owner(),
-                        completeMesh_.neighbour(),
-                        field.primitiveField(),
-                        faceAddressingBf_[procPatchi]
-                    )
-                )
-            );
-        }
-        else
-        {
-            FatalErrorInFunction
-                << "Unknown type." << abort(FatalError);
-        }
-    }
-
-    return tresF;
-}
-
-
-template<class Type>
-Foam::tmp<Foam::SurfaceField<Type>>
-Foam::fvFieldDecomposer::decomposeField
-(
-    const SurfaceField<Type>& field
-) const
-{
-    const SubList<label> faceAddressingIf
-    (
-        faceAddressing_,
-        procMesh_.nInternalFaces()
-    );
-
-    // Create dummy patch fields
-    PtrList<fvsPatchField<Type>> patchFields(procMesh_.boundary().size());
-    forAll(procMesh_.boundary(), procPatchi)
-    {
-        patchFields.set
-        (
-            procPatchi,
-            fvsPatchField<Type>::New
-            (
-                calculatedFvsPatchField<Type>::typeName,
-                procMesh_.boundary()[procPatchi],
-                DimensionedField<Type, surfaceMesh>::null()
-            )
-        );
-    }
-
-    // Create the processor field with the dummy patch fields
-    tmp<SurfaceField<Type>> tresF
-    (
-        new SurfaceField<Type>
-        (
-            IOobject
-            (
-                field.name(),
-                procMesh_.time().name(),
-                procMesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            procMesh_,
-            field.dimensions(),
-            mapFaceToFace
-            (
-                field,
-                faceAddressingIf,
-                isFlux(field)
-            ),
-            patchFields
-        )
-    );
-    SurfaceField<Type>& resF = tresF.ref();
-
-    // Change the patch fields to the correct type using a mapper constructor
-    // (with reference to the now correct internal field)
-    typename SurfaceField<Type>::
-        Boundary& bf = resF.boundaryFieldRef();
-    forAll(procMesh_.boundary(), procPatchi)
-    {
-        const fvPatch& procPatch = procMesh_.boundary()[procPatchi];
-
-        const label completePatchi = completePatchID(procPatchi);
-
-        if (completePatchi == procPatchi)
-        {
-            bf.set
-            (
-                procPatchi,
-                fvsPatchField<Type>::New
-                (
-                    field.boundaryField()[procPatchi],
-                    procPatch,
-                    resF(),
-                    patchFieldDecomposers_[procPatchi]
-                )
-            );
-        }
-        else if (isA<processorCyclicFvPatch>(procPatch))
-        {
-            bf.set
-            (
-                procPatchi,
-                new processorCyclicFvsPatchField<Type>
-                (
-                    procPatch,
-                    resF(),
-                    mapFaceToFace
-                    (
-                        field.boundaryField()[completePatchi],
-                        faceAddressingBf_[procPatchi],
-                        isFlux(field)
-                    )
-                )
-            );
-        }
-        else if (isA<processorFvPatch>(procPatch))
-        {
-            bf.set
-            (
-                procPatchi,
-                new processorFvsPatchField<Type>
-                (
-                    procPatch,
-                    resF(),
-                    mapFaceToFace
-                    (
-                        field.primitiveField(),
-                        faceAddressingBf_[procPatchi],
-                        isFlux(field)
-                    )
-                )
-            );
-        }
-        else
-        {
-            FatalErrorInFunction
-                << "Unknown type." << abort(FatalError);
-        }
-    }
-
-    return tresF;
-}
-
-
-template<class GeoField>
-void Foam::fvFieldDecomposer::decomposeFields
-(
-    const PtrList<GeoField>& fields
-) const
-{
-    forAll(fields, fieldi)
-    {
-        decomposeField(fields[fieldi])().write();
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/applications/utilities/parallelProcessing/decomposePar/fvFieldDecomposerTemplates.C b/applications/utilities/parallelProcessing/decomposePar/fvFieldDecomposerTemplates.C
new file mode 100644
index 0000000000..bc4f0b4491
--- /dev/null
+++ b/applications/utilities/parallelProcessing/decomposePar/fvFieldDecomposerTemplates.C
@@ -0,0 +1,550 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvFieldDecomposer.H"
+#include "processorFvPatchField.H"
+#include "processorFvsPatchField.H"
+#include "processorCyclicFvPatchField.H"
+#include "processorCyclicFvsPatchField.H"
+#include "emptyFvPatchFields.H"
+#include "stringOps.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>> Foam::fvFieldDecomposer::mapCellToFace
+(
+    const labelUList& owner,
+    const labelUList& neighbour,
+    const Field<Type>& field,
+    const labelUList& addressing
+)
+{
+    tmp<Field<Type>> tfld(new Field<Type>(addressing.size()));
+    Field<Type>& fld = tfld.ref();
+
+    forAll(addressing, i)
+    {
+        fld[i] =
+            field
+            [
+                addressing[i] > 0
+              ? neighbour[addressing[i] - 1]
+              : owner[- addressing[i] - 1]
+            ];
+    }
+
+    return tfld;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>> Foam::fvFieldDecomposer::mapFaceToFace
+(
+    const Field<Type>& field,
+    const labelUList& addressing,
+    const bool isFlux
+)
+{
+    tmp<Field<Type>> tfld(new Field<Type>(addressing.size()));
+    Field<Type>& fld = tfld.ref();
+
+    forAll(addressing, i)
+    {
+        fld[i] =
+            (isFlux && addressing[i] < 0 ? -1 : +1)
+           *field[mag(addressing[i]) - 1];
+    }
+
+    return tfld;
+}
+
+
+template<class Type>
+Foam::PtrList<typename Foam::VolField<Type>::Internal>
+Foam::fvFieldDecomposer::decomposeVolInternalField
+(
+    const IOobject& fieldIoObject
+) const
+{
+    // Read the field
+    const typename VolField<Type>::Internal field
+    (
+        IOobject
+        (
+            fieldIoObject.name(),
+            completeMesh_.time().name(),
+            completeMesh_,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        completeMesh_
+    );
+
+    // Construct the processor fields
+    PtrList<typename VolField<Type>::Internal> procFields(procMeshes_.size());
+    forAll(procMeshes_, proci)
+    {
+        // Create the processor field with the dummy patch fields
+        procFields.set
+        (
+            proci,
+            new typename VolField<Type>::Internal
+            (
+                IOobject
+                (
+                    fieldIoObject.name(),
+                    procMeshes_[proci].time().name(),
+                    procMeshes_[proci],
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                procMeshes_[proci],
+                field.dimensions(),
+                Field<Type>(field.field(), cellProcAddressing_[proci])
+            )
+        );
+    }
+
+    return procFields;
+}
+
+
+template<class Type>
+Foam::PtrList<Foam::VolField<Type>>
+Foam::fvFieldDecomposer::decomposeVolField
+(
+    const IOobject& fieldIoObject
+) const
+{
+    // Read the field
+    const VolField<Type> field
+    (
+        IOobject
+        (
+            fieldIoObject.name(),
+            completeMesh_.time().name(),
+            completeMesh_,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        completeMesh_
+    );
+
+    // Construct the processor fields
+    PtrList<VolField<Type>> procFields(procMeshes_.size());
+    forAll(procMeshes_, proci)
+    {
+        // Create dummy patch fields
+        PtrList<fvPatchField<Type>> patchFields
+        (
+            procMeshes_[proci].boundary().size()
+        );
+        forAll(procMeshes_[proci].boundary(), procPatchi)
+        {
+            patchFields.set
+            (
+                procPatchi,
+                fvPatchField<Type>::New
+                (
+                    calculatedFvPatchField<Type>::typeName,
+                    procMeshes_[proci].boundary()[procPatchi],
+                    DimensionedField<Type, volMesh>::null()
+                )
+            );
+        }
+
+        // Create the processor field with the dummy patch fields
+        procFields.set
+        (
+            proci,
+            new VolField<Type>
+            (
+                IOobject
+                (
+                    fieldIoObject.name(),
+                    procMeshes_[proci].time().name(),
+                    procMeshes_[proci],
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                procMeshes_[proci],
+                field.dimensions(),
+                Field<Type>(field.primitiveField(), cellProcAddressing_[proci]),
+                patchFields
+            )
+        );
+
+        // Alias the created proc field
+        VolField<Type>& vf = procFields[proci];
+
+        // Change the patch fields to the correct type using a mapper
+        // constructor (with reference to the now correct internal field)
+        typename VolField<Type>::Boundary& bf = vf.boundaryFieldRef();
+        forAll(bf, procPatchi)
+        {
+            const fvPatch& procPatch =
+                procMeshes_[proci].boundary()[procPatchi];
+
+            const label completePatchi = completePatchID(proci, procPatchi);
+
+            if (completePatchi == procPatchi)
+            {
+                bf.set
+                (
+                    procPatchi,
+                    fvPatchField<Type>::New
+                    (
+                        field.boundaryField()[completePatchi],
+                        procPatch,
+                        vf(),
+                        patchFieldDecomposers_[proci][procPatchi]
+                    )
+                );
+            }
+            else if (isA<processorCyclicFvPatch>(procPatch))
+            {
+                if (field.boundaryField()[completePatchi].overridesConstraint())
+                {
+                    OStringStream str;
+                    str << "\nThe field \"" << field.name()
+                        << "\" on cyclic patch \""
+                        << field.boundaryField()[completePatchi].patch().name()
+                        << "\" cannot be decomposed as it is not a cyclic "
+                        << "patch field. A \"patchType cyclic;\" setting has "
+                        << "been used to override the cyclic patch type.\n\n"
+                        << "Cyclic patches like this with non-cyclic boundary "
+                        << "conditions should be confined to a single "
+                        << "processor using decomposition constraints.";
+                    FatalErrorInFunction
+                        << stringOps::breakIntoIndentedLines(str.str()).c_str()
+                        << exit(FatalError);
+                }
+
+                const label nbrCompletePatchi =
+                    refCast<const processorCyclicFvPatch>(procPatch)
+                   .referPatch().nbrPatchID();
+
+                bf.set
+                (
+                    procPatchi,
+                    new processorCyclicFvPatchField<Type>
+                    (
+                        procPatch,
+                        vf(),
+                        mapCellToFace
+                        (
+                            labelUList(),
+                            completeMesh_.lduAddr().patchAddr
+                            (
+                                nbrCompletePatchi
+                            ),
+                            field.primitiveField(),
+                            faceProcAddressingBf_[proci][procPatchi]
+                        )
+                    )
+                );
+            }
+            else if (isA<processorFvPatch>(procPatch))
+            {
+                bf.set
+                (
+                    procPatchi,
+                    new processorFvPatchField<Type>
+                    (
+                        procPatch,
+                        vf(),
+                        mapCellToFace
+                        (
+                            completeMesh_.owner(),
+                            completeMesh_.neighbour(),
+                            field.primitiveField(),
+                            faceProcAddressingBf_[proci][procPatchi]
+                        )
+                    )
+                );
+            }
+            else
+            {
+                FatalErrorInFunction
+                    << "Unknown type." << abort(FatalError);
+            }
+        }
+    }
+
+    return procFields;
+}
+
+
+template<class Type>
+Foam::PtrList<Foam::SurfaceField<Type>>
+Foam::fvFieldDecomposer::decomposeFvSurfaceField
+(
+    const IOobject& fieldIoObject
+) const
+{
+    // Read the field
+    const SurfaceField<Type> field
+    (
+        IOobject
+        (
+            fieldIoObject.name(),
+            completeMesh_.time().name(),
+            completeMesh_,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        completeMesh_
+    );
+
+    // Construct the processor fields
+    PtrList<SurfaceField<Type>> procFields(procMeshes_.size());
+    forAll(procMeshes_, proci)
+    {
+        const SubList<label> faceAddressingIf
+        (
+            faceProcAddressing_[proci],
+            procMeshes_[proci].nInternalFaces()
+        );
+
+        // Create dummy patch fields
+        PtrList<fvsPatchField<Type>> patchFields
+        (
+            procMeshes_[proci].boundary().size()
+        );
+        forAll(procMeshes_[proci].boundary(), procPatchi)
+        {
+            patchFields.set
+            (
+                procPatchi,
+                fvsPatchField<Type>::New
+                (
+                    calculatedFvsPatchField<Type>::typeName,
+                    procMeshes_[proci].boundary()[procPatchi],
+                    DimensionedField<Type, surfaceMesh>::null()
+                )
+            );
+        }
+
+        // Create the processor field with the dummy patch fields
+        procFields.set
+        (
+            proci,
+            new SurfaceField<Type>
+            (
+                IOobject
+                (
+                    field.name(),
+                    procMeshes_[proci].time().name(),
+                    procMeshes_[proci],
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                procMeshes_[proci],
+                field.dimensions(),
+                mapFaceToFace
+                (
+                    field,
+                    faceAddressingIf,
+                    isFlux(field)
+                ),
+                patchFields
+            )
+        );
+
+        // Alias the created proc field
+        SurfaceField<Type>& sf = procFields[proci];
+
+        // Change the patch fields to the correct type using a mapper
+        // constructor (with reference to the now correct internal field)
+        typename SurfaceField<Type>::Boundary& bf = sf.boundaryFieldRef();
+        forAll(procMeshes_[proci].boundary(), procPatchi)
+        {
+            const fvPatch& procPatch =
+                procMeshes_[proci].boundary()[procPatchi];
+
+            const label completePatchi = completePatchID(proci, procPatchi);
+
+            if (completePatchi == procPatchi)
+            {
+                bf.set
+                (
+                    procPatchi,
+                    fvsPatchField<Type>::New
+                    (
+                        field.boundaryField()[procPatchi],
+                        procPatch,
+                        sf(),
+                        patchFieldDecomposers_[proci][procPatchi]
+                    )
+                );
+            }
+            else if (isA<processorCyclicFvPatch>(procPatch))
+            {
+                bf.set
+                (
+                    procPatchi,
+                    new processorCyclicFvsPatchField<Type>
+                    (
+                        procPatch,
+                        sf(),
+                        mapFaceToFace
+                        (
+                            field.boundaryField()[completePatchi],
+                            faceProcAddressingBf_[proci][procPatchi],
+                            isFlux(field)
+                        )
+                    )
+                );
+            }
+            else if (isA<processorFvPatch>(procPatch))
+            {
+                bf.set
+                (
+                    procPatchi,
+                    new processorFvsPatchField<Type>
+                    (
+                        procPatch,
+                        sf(),
+                        mapFaceToFace
+                        (
+                            field.primitiveField(),
+                            faceProcAddressingBf_[proci][procPatchi],
+                            isFlux(field)
+                        )
+                    )
+                );
+            }
+            else
+            {
+                FatalErrorInFunction
+                    << "Unknown type." << abort(FatalError);
+            }
+        }
+    }
+
+    return procFields;
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::fvFieldDecomposer::decomposeVolInternalFields
+(
+    const IOobjectList& objects
+)
+{
+    const word& fieldClassName = VolField<Type>::Internal::typeName;
+
+    IOobjectList fields = objects.lookupClass(fieldClassName);
+
+    if (fields.size())
+    {
+        Info<< nl << "    Decomposing " << fieldClassName << "s" << nl << endl;
+
+        forAllConstIter(IOobjectList, fields, fieldIter)
+        {
+            Info<< "        " << fieldIter()->name() << endl;
+
+            PtrList<typename VolField<Type>::Internal> procFields =
+                decomposeVolInternalField<Type>(*fieldIter());
+
+            forAll(procFields, proci)
+            {
+                procFields[proci].write();
+            }
+        }
+    }
+}
+
+
+template<class Type>
+void Foam::fvFieldDecomposer::decomposeVolFields
+(
+    const IOobjectList& objects
+)
+{
+    const word& fieldClassName = VolField<Type>::typeName;
+
+    IOobjectList fields = objects.lookupClass(fieldClassName);
+
+    if (fields.size())
+    {
+        Info<< nl << "    Decomposing " << fieldClassName << "s" << nl << endl;
+
+        forAllConstIter(IOobjectList, fields, fieldIter)
+        {
+            Info<< "        " << fieldIter()->name() << endl;
+
+            PtrList<VolField<Type>> procFields =
+                decomposeVolField<Type>(*fieldIter());
+
+            forAll(procFields, proci)
+            {
+                procFields[proci].write();
+            }
+        }
+    }
+}
+
+
+template<class Type>
+void Foam::fvFieldDecomposer::decomposeFvSurfaceFields
+(
+    const IOobjectList& objects
+)
+{
+    const word& fieldClassName = SurfaceField<Type>::typeName;
+
+    IOobjectList fields = objects.lookupClass(fieldClassName);
+
+    if (fields.size())
+    {
+        Info<< nl << "    Decomposing " << fieldClassName << "s" << nl << endl;
+
+        forAllConstIter(IOobjectList, fields, fieldIter)
+        {
+            Info<< "        " << fieldIter()->name() << endl;
+
+            PtrList<SurfaceField<Type>> procFields =
+                decomposeFvSurfaceField<Type>(*fieldIter());
+
+            forAll(procFields, proci)
+            {
+                procFields[proci].write();
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposer.C b/applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposer.C
index 97ee475d74..386f1b6dcb 100644
--- a/applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposer.C
+++ b/applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposer.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -21,83 +21,161 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-Description
-    Lagrangian field decomposer.
-
 \*---------------------------------------------------------------------------*/
 
 #include "lagrangianFieldDecomposer.H"
+#include "passiveParticleCloud.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::lagrangianFieldDecomposer::lagrangianFieldDecomposer
 (
-    const polyMesh& mesh,
-    const polyMesh& procMesh,
-    const labelList& faceProcAddressing,
-    const labelList& cellProcAddressing,
-    const word& cloudName,
-    const Cloud<indexedParticle>& lagrangianPositions,
-    const List<SLList<indexedParticle*>*>& cellParticles
+    const fvMesh& completeMesh,
+    const PtrList<fvMesh>& procMeshes,
+    const labelListList& faceProcAddressing,
+    const labelListList& cellProcAddressing,
+    const word& cloudName
 )
 :
-    procMesh_(procMesh),
-    positions_(procMesh, cloudName, IDLList<passiveParticle>()),
-    particleIndices_(lagrangianPositions.size())
+    completeMesh_(completeMesh),
+    procMeshes_(procMeshes),
+    particleProcAddressing_(procMeshes_.size()),
+    cloudName_(cloudName)
 {
-    label pi = 0;
-
-    labelList decodedProcFaceAddressing(faceProcAddressing.size());
-
-    forAll(faceProcAddressing, i)
+    // Create reverse cell addressing
+    List<remote> cellProcCell(completeMesh_.nCells());
+    forAll(cellProcAddressing, proci)
     {
-        decodedProcFaceAddressing[i] = mag(faceProcAddressing[i]) - 1;
-    }
-
-    forAll(cellProcAddressing, procCelli)
-    {
-        label celli = cellProcAddressing[procCelli];
-
-        if (cellParticles[celli])
+        forAll(cellProcAddressing[proci], procCelli)
         {
-            SLList<indexedParticle*>& particlePtrs = *cellParticles[celli];
-
-            forAllConstIter(SLList<indexedParticle*>, particlePtrs, iter)
-            {
-                const indexedParticle& ppi = *iter();
-                particleIndices_[pi++] = ppi.index();
-
-                label mappedTetFace = findIndex
-                (
-                    decodedProcFaceAddressing,
-                    ppi.tetFace()
-                );
-
-                if (mappedTetFace == -1)
-                {
-                    FatalErrorInFunction
-                        << "Face lookup failure." << nl
-                        << abort(FatalError);
-                }
-
-                positions_.append
-                (
-                    new passiveParticle
-                    (
-                        procMesh,
-                        ppi.coordinates(),
-                        procCelli,
-                        mappedTetFace,
-                        ppi.procTetPt(mesh, procMesh, procCelli, mappedTetFace)
-                    )
-                );
-            }
+            cellProcCell[cellProcAddressing[proci][procCelli]] =
+                remote(proci, procCelli);
         }
     }
 
-    particleIndices_.setSize(pi);
+    // Create reverse face addressing
+    List<remote> faceOwnerProcFace(completeMesh_.nFaces());
+    List<remote> faceNeighbourProcFace(completeMesh_.nFaces());
+    forAll(faceProcAddressing, proci)
+    {
+        forAll(faceProcAddressing[proci], procFacei)
+        {
+            const bool owner = faceProcAddressing[proci][procFacei] > 0;
+            const label facei = mag(faceProcAddressing[proci][procFacei]) - 1;
 
-    IOPosition<Cloud<passiveParticle>>(positions_).write();
+            (owner ? faceOwnerProcFace : faceNeighbourProcFace)[facei] =
+                remote(proci, procFacei);
+        }
+    }
+
+    // Read the complete positions
+    const passiveParticleCloud completePositions
+    (
+        completeMesh_,
+        cloudName_,
+        false
+    );
+
+    // Construct empty clouds for processor positions
+    PtrList<passiveParticleCloud> procPositions(procMeshes_.size());
+    forAll(procMeshes_, proci)
+    {
+        procPositions.set
+        (
+            proci,
+            new passiveParticleCloud
+            (
+                procMeshes_[proci],
+                cloudName_,
+                IDLList<passiveParticle>()
+            )
+        );
+    }
+
+    // Count the number of particles on each processor
+    labelList procNParticles(procMeshes_.size(), 0);
+    forAllConstIter(passiveParticleCloud, completePositions, iter)
+    {
+        const passiveParticle& p = iter();
+        const label proci = cellProcCell[p.cell()].proci;
+
+        procNParticles[proci] ++;
+    }
+
+    // Resize the addressing
+    forAll(procMeshes_, proci)
+    {
+        particleProcAddressing_[proci].resize(procNParticles[proci], -1);
+    }
+
+    // Distribute positions to the processor meshes
+    label completeParticlei = 0;
+    labelList procParticlei(procMeshes_.size(), 0);
+    forAllConstIter(passiveParticleCloud, completePositions, iter)
+    {
+        const passiveParticle& p = iter();
+        const label proci = cellProcCell[p.cell()].proci;
+        const label procCelli = cellProcCell[p.cell()].elementi;
+        const label procFacei =
+            faceOwnerProcFace[p.tetFace()].proci == proci
+          ? faceOwnerProcFace[p.tetFace()].elementi
+          : faceNeighbourProcFace[p.tetFace()].elementi;
+
+        particleProcAddressing_[proci][procParticlei[proci]] =
+            completeParticlei;
+
+        procPositions[proci].append
+        (
+            new passiveParticle
+            (
+                procMeshes_[proci],
+                p.coordinates(),
+                procCelli,
+                procFacei,
+                p.procTetPt
+                (
+                    completeMesh_,
+                    procMeshes_[proci],
+                    procCelli,
+                    procFacei
+                )
+            )
+        );
+
+        completeParticlei ++;
+        procParticlei[proci] ++;
+    }
+
+    // Write
+    forAll(procPositions, proci)
+    {
+        IOPosition<passiveParticleCloud>(procPositions[proci]).write();
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::lagrangianFieldDecomposer::~lagrangianFieldDecomposer()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::lagrangianFieldDecomposer::decomposes(const IOobjectList& objects)
+{
+    bool result = false;
+
+    #define DO_LAGRANGIAN_FIELDS_TYPE(Type, nullArg)                           \
+        result = result                                                        \
+         || !objects.lookupClass(IOField<Type>::typeName).empty()              \
+         || !objects.lookupClass(IOField<Field<Type>>::typeName).empty()       \
+         || !objects.lookupClass(CompactIOField<Field<Type>>::typeName).empty();
+    DO_LAGRANGIAN_FIELDS_TYPE(label, )
+    FOR_ALL_FIELD_TYPES(DO_LAGRANGIAN_FIELDS_TYPE)
+    #undef DO_LAGRANGIAN_FIELDS_TYPE
+
+    return result;
 }
 
 
diff --git a/applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposer.H b/applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposer.H
index 09660a8f75..27a5f35a71 100644
--- a/applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposer.H
+++ b/applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposer.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,17 +29,15 @@ Description
 
 SourceFiles
     lagrangianFieldDecomposer.C
-    lagrangianFieldDecomposerDecomposeFields.C
+    lagrangianFieldDecomposerTemplates.C
 
 \*---------------------------------------------------------------------------*/
 
 #ifndef lagrangianFieldDecomposer_H
 #define lagrangianFieldDecomposer_H
 
-#include "Cloud.H"
-#include "CompactIOField.H"
-#include "indexedParticle.H"
-#include "passiveParticle.H"
+#include "cloud.H"
+#include "fvMesh.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -56,14 +54,34 @@ class lagrangianFieldDecomposer
 {
     // Private Data
 
-        //- Reference to processor mesh
-        const polyMesh& procMesh_;
+        //- Reference to complete mesh
+        const fvMesh& completeMesh_;
 
-        //- Lagrangian positions for this processor
-        Cloud<passiveParticle> positions_;
+        //- List of processor meshes
+        const PtrList<fvMesh>& procMeshes_;
 
-        //- The indices of the particles on this processor
-        labelList particleIndices_;
+        //- For each processor, for each particle, the global particle index
+        labelListList particleProcAddressing_;
+
+        //- The name of the cloud
+        const word cloudName_;
+
+
+    // Private Member Functions
+
+        //- Decompose a field
+        template<class Type, template<class> class IOContainer>
+        PtrList<IOContainer<Type>>
+        decomposeField(const IOobject& fieldIoObject) const;
+
+        //- Read, decompose and write all fields
+        template
+        <
+            class Type,
+            template<class> class IOContainer,
+            template<class> class IOContainerType
+        >
+        void decomposeFields(const IOobjectList& objects);
 
 
 public:
@@ -73,68 +91,29 @@ public:
         //- Construct from components
         lagrangianFieldDecomposer
         (
-            const polyMesh& mesh,
-            const polyMesh& procMesh,
-            const labelList& faceProcAddressing,
-            const labelList& cellProcAddressing,
-            const word& cloudName,
-            const Cloud<indexedParticle>& lagrangianPositions,
-            const List<SLList<indexedParticle*>*>& cellParticles
+            const fvMesh& completeMesh,
+            const PtrList<fvMesh>& procMeshes,
+            const labelListList& faceProcAddressing,
+            const labelListList& cellProcAddressing,
+            const word& cloudName
         );
 
         //- Disallow default bitwise copy construction
         lagrangianFieldDecomposer(const lagrangianFieldDecomposer&) = delete;
 
 
+    //- Destructor
+    ~lagrangianFieldDecomposer();
+
+
     // Member Functions
 
-        // Read the fields and hold on the pointer list
+        //- Return whether anything in the object list gets decomposed
+        static bool decomposes(const IOobjectList& objects);
+
+        //- Read, decompose and write all fields
         template<class Type>
-        static void readFields
-        (
-            const label cloudI,
-            const IOobjectList& lagrangianObjects,
-            PtrList<PtrList<IOField<Type>>>& lagrangianFields
-        );
-
-        template<class Type>
-        static void readFieldFields
-        (
-            const label cloudI,
-            const IOobjectList& lagrangianObjects,
-            PtrList<PtrList<CompactIOField<Field<Type>>>>& lagrangianFields
-        );
-
-
-        //- Decompose volume field
-        template<class Type>
-        tmp<IOField<Type>> decomposeField
-        (
-            const word& cloudName,
-            const IOField<Type>& field
-        ) const;
-
-        template<class Type>
-        tmp<CompactIOField<Field<Type>>> decomposeFieldField
-        (
-            const word& cloudName,
-            const CompactIOField<Field<Type>>& field
-        ) const;
-
-
-        template<class GeoField>
-        void decomposeFields
-        (
-            const word& cloudName,
-            const PtrList<GeoField>& fields
-        ) const;
-
-        template<class GeoField>
-        void decomposeFieldFields
-        (
-            const word& cloudName,
-            const PtrList<GeoField>& fields
-        ) const;
+        void decomposeFields(const IOobjectList& objects);
 
 
     // Member Operators
@@ -151,7 +130,7 @@ public:
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-    #include "lagrangianFieldDecomposerDecomposeFields.C"
+    #include "lagrangianFieldDecomposerTemplates.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposerDecomposeFields.C b/applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposerDecomposeFields.C
deleted file mode 100644
index afed75e1d3..0000000000
--- a/applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposerDecomposeFields.C
+++ /dev/null
@@ -1,213 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "lagrangianFieldDecomposer.H"
-#include "IOobjectList.H"
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class Type>
-void Foam::lagrangianFieldDecomposer::readFields
-(
-    const label cloudI,
-    const IOobjectList& lagrangianObjects,
-    PtrList<PtrList<IOField<Type>>>& lagrangianFields
-)
-{
-    // Search list of objects for lagrangian fields
-    IOobjectList lagrangianTypeObjects
-    (
-        lagrangianObjects.lookupClass(IOField<Type>::typeName)
-    );
-
-    lagrangianFields.set
-    (
-        cloudI,
-        new PtrList<IOField<Type>>
-        (
-            lagrangianTypeObjects.size()
-        )
-    );
-
-    label lagrangianFieldi = 0;
-    forAllIter(IOobjectList, lagrangianTypeObjects, iter)
-    {
-        lagrangianFields[cloudI].set
-        (
-            lagrangianFieldi++,
-            new IOField<Type>(*iter())
-        );
-    }
-}
-
-
-template<class Type>
-void Foam::lagrangianFieldDecomposer::readFieldFields
-(
-    const label cloudI,
-    const IOobjectList& lagrangianObjects,
-    PtrList<PtrList<CompactIOField<Field<Type>>>>& lagrangianFields
-)
-{
-    // Search list of objects for lagrangian fields
-    IOobjectList lagrangianTypeObjectsA
-    (
-        lagrangianObjects.lookupClass(IOField<Field<Type>>::typeName)
-    );
-
-    IOobjectList lagrangianTypeObjectsB
-    (
-        lagrangianObjects.lookupClass
-        (
-            CompactIOField<Field<Type>>::typeName
-        )
-    );
-
-    lagrangianFields.set
-    (
-        cloudI,
-        new PtrList<CompactIOField<Field<Type>>>
-        (
-            lagrangianTypeObjectsA.size() + lagrangianTypeObjectsB.size()
-        )
-    );
-
-    label lagrangianFieldi = 0;
-
-    forAllIter(IOobjectList, lagrangianTypeObjectsA, iter)
-    {
-        lagrangianFields[cloudI].set
-        (
-            lagrangianFieldi++,
-            new CompactIOField<Field<Type>>(*iter())
-        );
-    }
-
-    forAllIter(IOobjectList, lagrangianTypeObjectsB, iter)
-    {
-        lagrangianFields[cloudI].set
-        (
-            lagrangianFieldi++,
-            new CompactIOField<Field<Type>>(*iter())
-        );
-    }
-}
-
-
-template<class Type>
-Foam::tmp<Foam::IOField<Type>>
-Foam::lagrangianFieldDecomposer::decomposeField
-(
-    const word& cloudName,
-    const IOField<Type>& field
-) const
-{
-    // Create and map the internal field values
-    Field<Type> procField(field, particleIndices_);
-
-    // Create the field for the processor
-    return tmp<IOField<Type>>
-    (
-        new IOField<Type>
-        (
-            IOobject
-            (
-                field.name(),
-                procMesh_.time().name(),
-                cloud::prefix/cloudName,
-                procMesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            procField
-        )
-    );
-}
-
-
-template<class Type>
-Foam::tmp<Foam::CompactIOField<Foam::Field<Type>>>
-Foam::lagrangianFieldDecomposer::decomposeFieldField
-(
-    const word& cloudName,
-    const CompactIOField<Field<Type>>& field
-) const
-{
-    // Create and map the internal field values
-    Field<Field<Type>> procField(field, particleIndices_);
-
-    // Create the field for the processor
-    return tmp<CompactIOField<Field<Type>>>
-    (
-        new CompactIOField<Field<Type>>
-        (
-            IOobject
-            (
-                field.name(),
-                procMesh_.time().name(),
-                cloud::prefix/cloudName,
-                procMesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            procField
-        )
-    );
-}
-
-
-template<class GeoField>
-void Foam::lagrangianFieldDecomposer::decomposeFields
-(
-    const word& cloudName,
-    const PtrList<GeoField>& fields
-) const
-{
-    const bool write = particleIndices_.size() > 0;
-    forAll(fields, fieldi)
-    {
-        decomposeField(cloudName, fields[fieldi])().write(write);
-    }
-}
-
-
-template<class GeoField>
-void Foam::lagrangianFieldDecomposer::decomposeFieldFields
-(
-    const word& cloudName,
-    const PtrList<GeoField>& fields
-) const
-{
-    const bool write = particleIndices_.size() > 0;
-    forAll(fields, fieldi)
-    {
-        decomposeFieldField(cloudName, fields[fieldi])().write(write);
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposerTemplates.C b/applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposerTemplates.C
new file mode 100644
index 0000000000..7e06cc5327
--- /dev/null
+++ b/applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposerTemplates.C
@@ -0,0 +1,119 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "lagrangianFieldDecomposer.H"
+#include "IOobjectList.H"
+#include "CompactIOField.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type, template<class> class IOContainer>
+Foam::PtrList<IOContainer<Type>>
+Foam::lagrangianFieldDecomposer::decomposeField
+(
+    const IOobject& fieldIoObject
+) const
+{
+    // Read the complete field
+    const IOContainer<Type> field(fieldIoObject);
+
+    // Construct the processor fields
+    PtrList<IOContainer<Type>> procFields(procMeshes_.size());
+    forAll(procMeshes_, proci)
+    {
+        procFields.set
+        (
+            proci,
+            new IOContainer<Type>
+            (
+                IOobject
+                (
+                    field.name(),
+                    procMeshes_[proci].time().name(),
+                    cloud::prefix/cloudName_,
+                    procMeshes_[proci],
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                Field<Type>(field, particleProcAddressing_[proci])
+            )
+        );
+    }
+
+    return procFields;
+}
+
+
+template
+<
+    class Type,
+    template<class> class IOContainer,
+    template<class> class IOContainerType
+>
+void Foam::lagrangianFieldDecomposer::decomposeFields
+(
+    const IOobjectList& objects
+)
+{
+    const word& fieldClassName = IOContainerType<Type>::typeName;
+
+    IOobjectList fields = objects.lookupClass(fieldClassName);
+
+    if (fields.size())
+    {
+        Info<< nl << "    Decomposing " << fieldClassName << "s" << nl << endl;
+
+        forAllConstIter(IOobjectList, fields, fieldIter)
+        {
+            Info<< "        " << fieldIter()->name() << endl;
+
+            PtrList<IOContainer<Type>> procFields =
+                decomposeField<Type, IOContainer>(*fieldIter());
+
+            forAll(procFields, proci)
+            {
+                procFields[proci].write();
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::lagrangianFieldDecomposer::decomposeFields
+(
+    const IOobjectList& objects
+)
+{
+    decomposeFields<Type, IOField, IOField>(objects);
+    decomposeFields<Field<Type>, CompactIOField, IOField>(objects);
+    decomposeFields<Field<Type>, CompactIOField, CompactIOField>(objects);
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/parallelProcessing/decomposePar/pointFieldDecomposer.C b/applications/utilities/parallelProcessing/decomposePar/pointFieldDecomposer.C
index 5f2b9d7260..5e13703208 100644
--- a/applications/utilities/parallelProcessing/decomposePar/pointFieldDecomposer.C
+++ b/applications/utilities/parallelProcessing/decomposePar/pointFieldDecomposer.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "pointFieldDecomposer.H"
+#include "fvMesh.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -77,29 +78,40 @@ Foam::pointFieldDecomposer::patchFieldDecomposer::patchFieldDecomposer
 Foam::pointFieldDecomposer::pointFieldDecomposer
 (
     const pointMesh& completeMesh,
-    const pointMesh& procMesh,
-    const labelList& pointAddressing
+    const PtrList<fvMesh>& procMeshes,
+    const labelListList& pointProcAddressing
 )
 :
     completeMesh_(completeMesh),
-    procMesh_(procMesh),
-    pointAddressing_(pointAddressing),
-    patchFieldDecomposers_(procMesh_.boundary().size())
+    procMeshes_(procMeshes),
+    pointProcAddressing_(pointProcAddressing),
+    patchFieldDecomposers_(procMeshes_.size())
 {
-    forAll(procMesh_.boundary(), patchi)
+    forAll(procMeshes_, proci)
     {
-        if (patchi < completeMesh_.boundary().size())
+        const pointMesh& procMesh = pointMesh::New(procMeshes_[proci]);
+
+        patchFieldDecomposers_.set
+        (
+            proci,
+            new PtrList<patchFieldDecomposer>(procMesh.boundary().size())
+        );
+
+        forAll(procMesh.boundary(), procPatchi)
         {
-            patchFieldDecomposers_.set
-            (
-                patchi,
-                new patchFieldDecomposer
+            if (procPatchi < completeMesh_.boundary().size())
+            {
+                patchFieldDecomposers_[proci].set
                 (
-                    completeMesh_.boundary()[patchi],
-                    procMesh_.boundary()[patchi],
-                    pointAddressing_
-                )
-            );
+                    procPatchi,
+                    new patchFieldDecomposer
+                    (
+                        completeMesh_.boundary()[procPatchi],
+                        procMesh.boundary()[procPatchi],
+                        pointProcAddressing_[proci]
+                    )
+                );
+            }
         }
     }
 }
@@ -111,4 +123,20 @@ Foam::pointFieldDecomposer::~pointFieldDecomposer()
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::pointFieldDecomposer::decomposes(const IOobjectList& objects)
+{
+    bool result = false;
+
+    #define DO_POINT_FIELDS_TYPE(Type, nullArg)                                \
+        result = result                                                        \
+         || !objects.lookupClass(PointField<Type>::typeName).empty();
+    FOR_ALL_FIELD_TYPES(DO_POINT_FIELDS_TYPE)
+    #undef DO_POINT_FIELDS_TYPE
+
+    return result;
+}
+
+
 // ************************************************************************* //
diff --git a/applications/utilities/parallelProcessing/decomposePar/pointFieldDecomposer.H b/applications/utilities/parallelProcessing/decomposePar/pointFieldDecomposer.H
index b6e559b9e2..079f4c5b76 100644
--- a/applications/utilities/parallelProcessing/decomposePar/pointFieldDecomposer.H
+++ b/applications/utilities/parallelProcessing/decomposePar/pointFieldDecomposer.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,7 +29,7 @@ Description
 
 SourceFiles
     pointFieldDecomposer.C
-    pointFieldDecomposerDecomposeFields.C
+    pointFieldDecomposerTemplates.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -37,14 +37,17 @@ SourceFiles
 #define pointFieldDecomposer_H
 
 #include "pointMesh.H"
-#include "pointPatchFieldMapperPatchRef.H"
 #include "pointFields.H"
+#include "pointPatchFieldMapperPatchRef.H"
+#include "IOobjectList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
+class fvMesh;
+
 /*---------------------------------------------------------------------------*\
                     Class pointFieldDecomposer Declaration
 \*---------------------------------------------------------------------------*/
@@ -111,14 +114,22 @@ private:
         //- Reference to complete mesh
         const pointMesh& completeMesh_;
 
-        //- Reference to processor mesh
-        const pointMesh& procMesh_;
+        //- Reference to processor meshes
+        const PtrList<fvMesh>& procMeshes_;
 
         //- Reference to point addressing
-        const labelList& pointAddressing_;
+        const labelListList& pointProcAddressing_;
 
         //- List of patch field decomposers
-        PtrList<patchFieldDecomposer> patchFieldDecomposers_;
+        PtrList<PtrList<patchFieldDecomposer>> patchFieldDecomposers_;
+
+
+    // Private Member Functions
+
+        //- Decompose a field
+        template<class Type>
+        PtrList<PointField<Type>>
+        decomposeField(const IOobject& fieldIoObject) const;
 
 
 public:
@@ -129,8 +140,8 @@ public:
         pointFieldDecomposer
         (
             const pointMesh& completeMesh,
-            const pointMesh& procMesh,
-            const labelList& pointAddressing
+            const PtrList<fvMesh>& procMeshes,
+            const labelListList& pointAddressing
         );
 
         //- Disallow default bitwise copy construction
@@ -143,17 +154,12 @@ public:
 
     // Member Functions
 
-        //- Decompose point field
-        template<class Type>
-        tmp<PointField<Type>>
-        decomposeField
-        (
-            const PointField<Type>&
-        ) const;
+        //- Return whether anything in the object list gets decomposed
+        static bool decomposes(const IOobjectList& objects);
 
-        //- Decompose a list of fields
-        template<class GeoField>
-        void decomposeFields(const PtrList<GeoField>& fields) const;
+        //- Read, decompose and write all fields
+        template<class Type>
+        void decomposeFields(const IOobjectList& objects);
 
 
     // Member Operators
@@ -170,7 +176,7 @@ public:
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-    #include "pointFieldDecomposerDecomposeFields.C"
+    #include "pointFieldDecomposerTemplates.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/utilities/parallelProcessing/decomposePar/pointFieldDecomposerDecomposeFields.C b/applications/utilities/parallelProcessing/decomposePar/pointFieldDecomposerDecomposeFields.C
deleted file mode 100644
index 33b5c7afd1..0000000000
--- a/applications/utilities/parallelProcessing/decomposePar/pointFieldDecomposerDecomposeFields.C
+++ /dev/null
@@ -1,111 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "pointFieldDecomposer.H"
-#include "processorPointPatchFields.H"
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class Type>
-Foam::tmp<Foam::PointField<Type>>
-Foam::pointFieldDecomposer::decomposeField
-(
-    const PointField<Type>& field
-) const
-{
-    // Create and map the internal field values
-    Field<Type> internalField(field.primitiveField(), pointAddressing_);
-
-    // Create a list of pointers for the patchFields
-    PtrList<pointPatchField<Type>> patchFields(procMesh_.boundary().size());
-
-    // Create and map the patch field values
-    forAll(procMesh_.boundary(), patchi)
-    {
-        if (patchi < completeMesh_.boundary().size())
-        {
-            patchFields.set
-            (
-                patchi,
-                pointPatchField<Type>::New
-                (
-                    field.boundaryField()[patchi],
-                    procMesh_.boundary()[patchi],
-                    DimensionedField<Type, pointMesh>::null(),
-                    patchFieldDecomposers_[patchi]
-                )
-            );
-        }
-        else
-        {
-            patchFields.set
-            (
-                patchi,
-                new processorPointPatchField<Type>
-                (
-                    procMesh_.boundary()[patchi],
-                    DimensionedField<Type, pointMesh>::null()
-                )
-            );
-        }
-    }
-
-    // Create the field for the processor
-    return tmp<PointField<Type>>
-    (
-        new PointField<Type>
-        (
-            IOobject
-            (
-                field.name(),
-                procMesh_().time().name(),
-                procMesh_(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            procMesh_,
-            field.dimensions(),
-            internalField,
-            patchFields
-        )
-    );
-}
-
-
-template<class GeoField>
-void Foam::pointFieldDecomposer::decomposeFields
-(
-    const PtrList<GeoField>& fields
-) const
-{
-    forAll(fields, fieldi)
-    {
-        decomposeField(fields[fieldi])().write();
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/applications/utilities/parallelProcessing/decomposePar/pointFieldDecomposerTemplates.C b/applications/utilities/parallelProcessing/decomposePar/pointFieldDecomposerTemplates.C
new file mode 100644
index 0000000000..233bf76d7c
--- /dev/null
+++ b/applications/utilities/parallelProcessing/decomposePar/pointFieldDecomposerTemplates.C
@@ -0,0 +1,163 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "pointFieldDecomposer.H"
+#include "fvMesh.H"
+#include "processorPointPatchFields.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+Foam::PtrList<Foam::PointField<Type>>
+Foam::pointFieldDecomposer::decomposeField
+(
+    const IOobject& fieldIoObject
+) const
+{
+    // Read the field
+    const PointField<Type> field
+    (
+        IOobject
+        (
+            fieldIoObject.name(),
+            completeMesh_.db().time().name(),
+            completeMesh_.db(),
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        completeMesh_
+    );
+
+    // Construct the processor fields
+    PtrList<PointField<Type>> procFields(procMeshes_.size());
+    forAll(procMeshes_, proci)
+    {
+        const pointMesh& procMesh = pointMesh::New(procMeshes_[proci]);
+
+        // Create and map the internal field values
+        Field<Type> internalField
+        (
+            field.primitiveField(),
+            pointProcAddressing_[proci]
+        );
+
+        // Create a list of pointers for the patchFields
+        PtrList<pointPatchField<Type>> patchFields
+        (
+            procMesh.boundary().size()
+        );
+
+        // Create and map the patch field values
+        forAll(procMesh.boundary(), patchi)
+        {
+            if (patchi < completeMesh_.boundary().size())
+            {
+                patchFields.set
+                (
+                    patchi,
+                    pointPatchField<Type>::New
+                    (
+                        field.boundaryField()[patchi],
+                        procMesh.boundary()[patchi],
+                        DimensionedField<Type, pointMesh>::null(),
+                        patchFieldDecomposers_[proci][patchi]
+                    )
+                );
+            }
+            else
+            {
+                patchFields.set
+                (
+                    patchi,
+                    new processorPointPatchField<Type>
+                    (
+                        procMesh.boundary()[patchi],
+                        DimensionedField<Type, pointMesh>::null()
+                    )
+                );
+            }
+        }
+
+        // Create the field for the processor
+        procFields.set
+        (
+            proci,
+            new PointField<Type>
+            (
+                IOobject
+                (
+                    field.name(),
+                    procMesh().time().name(),
+                    procMesh(),
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                procMesh,
+                field.dimensions(),
+                internalField,
+                patchFields
+            )
+        );
+    }
+
+    return procFields;
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::pointFieldDecomposer::decomposeFields
+(
+    const IOobjectList& objects
+)
+{
+    const word& fieldClassName = PointField<Type>::typeName;
+
+    IOobjectList fields = objects.lookupClass(fieldClassName);
+
+    if (fields.size())
+    {
+        Info<< nl << "    Decomposing " << fieldClassName << "s" << nl << endl;
+
+        forAllConstIter(IOobjectList, fields, fieldIter)
+        {
+            Info<< "        " << fieldIter()->name() << endl;
+
+            PtrList<PointField<Type>> procFields =
+                decomposeField<Type>(*fieldIter());
+
+            forAll(procFields, proci)
+            {
+                procFields[proci].write();
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/parallelProcessing/decomposePar/readFields.C b/applications/utilities/parallelProcessing/decomposePar/readFields.C
deleted file mode 100644
index d8fa341279..0000000000
--- a/applications/utilities/parallelProcessing/decomposePar/readFields.C
+++ /dev/null
@@ -1,107 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "readFields.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-template<class Type, template<class> class PatchField, class GeoMesh>
-void Foam::readFields
-(
-    const typename GeoMesh::Mesh& mesh,
-    const IOobjectList& objects,
-    PtrList<GeometricField<Type, PatchField, GeoMesh>>& fields,
-    const bool readOldTime
-)
-{
-    // Search list of objects for fields
-    IOobjectList fieldObjects(objects.lookupClass
-    (
-        GeometricField<Type, PatchField, GeoMesh>::typeName)
-    );
-
-    // Remove the cellProc field
-    IOobjectList::iterator cellProcIter = fieldObjects.find("cellProc");
-    if (cellProcIter != fieldObjects.end())
-    {
-        fieldObjects.erase(cellProcIter);
-    }
-
-    // Get sorted set of names (different processors might read objects in
-    // different order)
-    const wordList masterNames(fieldObjects.sortedNames());
-
-    // Construct the fields
-    fields.setSize(masterNames.size());
-
-    forAll(masterNames, i)
-    {
-        const IOobject& io = *fieldObjects[masterNames[i]];
-
-        fields.set
-        (
-            i,
-            new GeometricField<Type, PatchField, GeoMesh>
-            (
-                io,
-                mesh,
-                readOldTime
-            )
-        );
-    }
-}
-
-
-template<class Mesh, class GeoField>
-void Foam::readFields
-(
-    const Mesh& mesh,
-    const IOobjectList& objects,
-    PtrList<GeoField>& fields
-)
-{
-    // Search list of objects for fields of type GeomField
-    IOobjectList fieldObjects(objects.lookupClass(GeoField::typeName));
-
-    // Construct the fields
-    fields.setSize(fieldObjects.size());
-
-    // Get sorted set of names (different processors might read objects in
-    // different order)
-    const wordList masterNames(fieldObjects.sortedNames());
-
-    // Construct the fields
-    fields.setSize(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/readFields.H b/applications/utilities/parallelProcessing/decomposePar/readFields.H
deleted file mode 100644
index 8f373d5554..0000000000
--- a/applications/utilities/parallelProcessing/decomposePar/readFields.H
+++ /dev/null
@@ -1,75 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Global
-    readFields
-
-Description
-
-SourceFiles
-    readFields.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef readFields_H
-#define readFields_H
-
-#include "IOobjectList.H"
-#include "PtrList.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    // Read the fields and hold on the pointer list
-    template<class Type, template<class> class PatchField, class GeoMesh>
-    void readFields
-    (
-        const typename GeoMesh::Mesh& mesh,
-        const IOobjectList& objects,
-        PtrList<GeometricField<Type, PatchField, GeoMesh>>& fields,
-        const bool readOldTime
-    );
-
-    // Read the fields and hold on the pointer list
-    template<class Mesh, class GeoField>
-    void readFields
-    (
-        const Mesh& mesh,
-        const IOobjectList& objects,
-        PtrList<GeoField>& fields
-    );
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-    #include "readFields.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/applications/utilities/parallelProcessing/reconstructPar/Make/files b/applications/utilities/parallelProcessing/reconstructPar/Make/files
index b0c0918a4b..f10f5fdd90 100644
--- a/applications/utilities/parallelProcessing/reconstructPar/Make/files
+++ b/applications/utilities/parallelProcessing/reconstructPar/Make/files
@@ -1,6 +1,6 @@
 reconstructPar.C
 fvFieldReconstructor.C
 pointFieldReconstructor.C
-reconstructLagrangianPositions.C
+lagrangianFieldReconstructor.C
 
 EXE = $(FOAM_APPBIN)/reconstructPar
diff --git a/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructor.C b/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructor.C
index fb713410f5..13c5f78bc6 100644
--- a/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructor.C
+++ b/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructor.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -66,8 +66,7 @@ Foam::fvFieldReconstructor::fvFieldReconstructor
     procMeshes_(procMeshes),
     faceProcAddressing_(faceProcAddressing),
     cellProcAddressing_(cellProcAddressing),
-    faceProcAddressingBf_(faceProcAddressingBf),
-    nReconstructed_(0)
+    faceProcAddressingBf_(faceProcAddressingBf)
 {
     forAll(procMeshes_, proci)
     {
@@ -93,4 +92,26 @@ Foam::fvFieldReconstructor::fvFieldReconstructor
 }
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::fvFieldReconstructor::reconstructs
+(
+    const IOobjectList& objects,
+    const HashSet<word>& selectedFields
+)
+{
+    bool result = false;
+
+    #define DO_FV_FIELDS_TYPE(Type, nullArg)                                   \
+        result = result                                                        \
+         || reconstructs<VolField<Type>::Internal>(objects, selectedFields)    \
+         || reconstructs<VolField<Type>>(objects, selectedFields)              \
+         || reconstructs<SurfaceField<Type>>(objects, selectedFields);
+    FOR_ALL_FIELD_TYPES(DO_FV_FIELDS_TYPE)
+    #undef DO_FV_FIELDS_TYPE
+
+    return result;
+}
+
+
 // ************************************************************************* //
diff --git a/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructor.H b/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructor.H
index 7fb34138fc..1577116cc8 100644
--- a/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructor.H
+++ b/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructor.H
@@ -29,7 +29,7 @@ Description
 
 SourceFiles
     fvFieldReconstructor.C
-    fvFieldReconstructorReconstructFields.C
+    fvFieldReconstructorTemplates.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -71,12 +71,17 @@ class fvFieldReconstructor
         //- Boundary field of face addressing
         const PtrList<surfaceLabelField::Boundary>& faceProcAddressingBf_;
 
-        //- Number of fields reconstructed
-        label nReconstructed_;
-
 
     // Private Member Functions
 
+        //- Return whether anything in the object list gets reconstructed
+        template<class FieldType>
+        static bool reconstructs
+        (
+            const IOobjectList& objects,
+            const HashSet<word>& selectedFields
+        );
+
         //- Convert a processor patch to the corresponding complete patch index
         label completePatchID(const label proci, const label procPatchi) const;
 
@@ -90,6 +95,21 @@ class fvFieldReconstructor
             const bool isFlux
         );
 
+        //- Read and reconstruct a volume internal field
+        template<class Type>
+        tmp<DimensionedField<Type, volMesh>>
+        reconstructVolInternalField(const IOobject& fieldIoObject) const;
+
+        //- Read and reconstruct a volume field
+        template<class Type>
+        tmp<VolField<Type>>
+        reconstructVolField(const IOobject& fieldIoObject) const;
+
+        //- Read and reconstruct a surface field
+        template<class Type>
+        tmp<SurfaceField<Type>>
+        reconstructFvSurfaceField(const IOobject& fieldIoObject) const;
+
 
 public:
 
@@ -111,59 +131,16 @@ public:
 
     // Member Functions
 
-        //- Return number of fields reconstructed
-        label nReconstructed() const
-        {
-            return nReconstructed_;
-        }
-
-        //- Reconstruct volume internal field
-        template<class Type>
-        tmp<DimensionedField<Type, volMesh>>
-        reconstructFvVolumeInternalField
+        //- Return whether anything in the object list gets reconstructed
+        static bool reconstructs
         (
-            const IOobject& fieldIoObject,
-            const PtrList<DimensionedField<Type, volMesh>>& procFields
-        ) const;
-
-        //- Read and reconstruct volume internal field
-        template<class Type>
-        tmp<DimensionedField<Type, volMesh>>
-        reconstructFvVolumeInternalField(const IOobject& fieldIoObject) const;
-
-
-        //- Reconstruct volume field
-        template<class Type>
-        tmp<VolField<Type>>
-        reconstructFvVolumeField
-        (
-            const IOobject& fieldIoObject,
-            const PtrList<VolField<Type>>&
-        ) const;
-
-        //- Read and reconstruct volume field
-        template<class Type>
-        tmp<VolField<Type>>
-        reconstructFvVolumeField(const IOobject& fieldIoObject) const;
-
-
-        //- Reconstruct surface field
-        template<class Type>
-        tmp<SurfaceField<Type>>
-        reconstructFvSurfaceField
-        (
-            const IOobject& fieldIoObject,
-            const PtrList<SurfaceField<Type>>&
-        ) const;
-
-        //- Read and reconstruct surface field
-        template<class Type>
-        tmp<SurfaceField<Type>>
-        reconstructFvSurfaceField(const IOobject& fieldIoObject) const;
+            const IOobjectList& objects,
+            const HashSet<word>& selectedFields
+        );
 
         //- Read, reconstruct and write all/selected volume internal fields
         template<class Type>
-        void reconstructFvVolumeInternalFields
+        void reconstructVolInternalFields
         (
             const IOobjectList& objects,
             const HashSet<word>& selectedFields
@@ -171,7 +148,7 @@ public:
 
         //- Read, reconstruct and write all/selected volume fields
         template<class Type>
-        void reconstructFvVolumeFields
+        void reconstructVolFields
         (
             const IOobjectList& objects,
             const HashSet<word>& selectedFields
@@ -200,7 +177,7 @@ public:
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-    #include "fvFieldReconstructorReconstructFields.C"
+    #include "fvFieldReconstructorTemplates.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructorReconstructFields.C b/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructorTemplates.C
similarity index 84%
rename from applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructorReconstructFields.C
rename to applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructorTemplates.C
index f55765dbf3..2412a6955b 100644
--- a/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructorReconstructFields.C
+++ b/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructorTemplates.C
@@ -34,7 +34,33 @@ License
 #include "reverseFvPatchFieldMapper.H"
 #include "stringOps.H"
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class FieldType>
+bool Foam::fvFieldReconstructor::reconstructs
+(
+    const IOobjectList& objects,
+    const HashSet<word>& selectedFields
+)
+{
+    IOobjectList fields = objects.lookupClass(FieldType::typeName);
+
+    if (fields.size() && selectedFields.empty())
+    {
+        return true;
+    }
+
+    forAllConstIter(IOobjectList, fields, fieldIter)
+    {
+        if (selectedFields.found(fieldIter()->name()))
+        {
+            return true;
+        }
+    }
+
+    return false;
+}
+
 
 template<class Type>
 void Foam::fvFieldReconstructor::rmapFaceToFace
@@ -53,54 +79,15 @@ void Foam::fvFieldReconstructor::rmapFaceToFace
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
 template<class Type>
 Foam::tmp<Foam::DimensionedField<Type, Foam::volMesh>>
-Foam::fvFieldReconstructor::reconstructFvVolumeInternalField
-(
-    const IOobject& fieldIoObject,
-    const PtrList<DimensionedField<Type, volMesh>>& procFields
-) const
-{
-    // Create the internalField
-    Field<Type> internalField(completeMesh_.nCells());
-
-    forAll(procMeshes_, proci)
-    {
-        const DimensionedField<Type, volMesh>& procField = procFields[proci];
-
-        // Set the cell values in the reconstructed field
-        internalField.rmap
-        (
-            procField.field(),
-            cellProcAddressing_[proci]
-        );
-    }
-
-    return tmp<DimensionedField<Type, volMesh>>
-    (
-        new DimensionedField<Type, volMesh>
-        (
-            fieldIoObject,
-            completeMesh_,
-            procFields[0].dimensions(),
-            internalField
-        )
-    );
-}
-
-
-template<class Type>
-Foam::tmp<Foam::DimensionedField<Type, Foam::volMesh>>
-Foam::fvFieldReconstructor::reconstructFvVolumeInternalField
+Foam::fvFieldReconstructor::reconstructVolInternalField
 (
     const IOobject& fieldIoObject
 ) const
 {
-    PtrList<DimensionedField<Type, volMesh>>
-        procFields(procMeshes_.size());
-
+    // Read the field for all the processors
+    PtrList<DimensionedField<Type, volMesh>> procFields(procMeshes_.size());
     forAll(procMeshes_, proci)
     {
         procFields.set
@@ -122,30 +109,72 @@ Foam::fvFieldReconstructor::reconstructFvVolumeInternalField
         );
     }
 
-    return reconstructFvVolumeInternalField
-    (
-        IOobject
+    // Create the internalField
+    Field<Type> internalField(completeMesh_.nCells());
+
+    forAll(procMeshes_, proci)
+    {
+        const DimensionedField<Type, volMesh>& procField = procFields[proci];
+
+        // Set the cell values in the reconstructed field
+        internalField.rmap
         (
-            fieldIoObject.name(),
-            completeMesh_.time().name(),
+            procField.field(),
+            cellProcAddressing_[proci]
+        );
+    }
+
+    return tmp<DimensionedField<Type, volMesh>>
+    (
+        new DimensionedField<Type, volMesh>
+        (
+            IOobject
+            (
+                fieldIoObject.name(),
+                completeMesh_.time().name(),
+                completeMesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
             completeMesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        procFields
+            procFields[0].dimensions(),
+            internalField
+        )
     );
 }
 
 
 template<class Type>
 Foam::tmp<Foam::VolField<Type>>
-Foam::fvFieldReconstructor::reconstructFvVolumeField
+Foam::fvFieldReconstructor::reconstructVolField
 (
-    const IOobject& fieldIoObject,
-    const PtrList<VolField<Type>>& procFields
+    const IOobject& fieldIoObject
 ) const
 {
+    // Read the field for all the processors
+    PtrList<VolField<Type>> procFields(procMeshes_.size());
+    forAll(procMeshes_, proci)
+    {
+        procFields.set
+        (
+            proci,
+            new VolField<Type>
+            (
+                IOobject
+                (
+                    fieldIoObject.name(),
+                    procMeshes_[proci].time().name(),
+                    procMeshes_[proci],
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                procMeshes_[proci]
+            )
+        );
+    }
+
     // Create the internalField
     Field<Type> internalField(completeMesh_.nCells());
 
@@ -251,7 +280,15 @@ Foam::fvFieldReconstructor::reconstructFvVolumeField
     (
         new VolField<Type>
         (
-            fieldIoObject,
+            IOobject
+            (
+                fieldIoObject.name(),
+                completeMesh_.time().name(),
+                completeMesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
             completeMesh_,
             procFields[0].dimensions(),
             internalField,
@@ -262,21 +299,20 @@ Foam::fvFieldReconstructor::reconstructFvVolumeField
 
 
 template<class Type>
-Foam::tmp<Foam::VolField<Type>>
-Foam::fvFieldReconstructor::reconstructFvVolumeField
+Foam::tmp<Foam::SurfaceField<Type>>
+Foam::fvFieldReconstructor::reconstructFvSurfaceField
 (
     const IOobject& fieldIoObject
 ) const
 {
-    PtrList<VolField<Type>>
-        procFields(procMeshes_.size());
-
+    // Read the field for all the processors
+    PtrList<SurfaceField<Type>> procFields(procMeshes_.size());
     forAll(procMeshes_, proci)
     {
         procFields.set
         (
             proci,
-            new VolField<Type>
+            new SurfaceField<Type>
             (
                 IOobject
                 (
@@ -292,30 +328,6 @@ Foam::fvFieldReconstructor::reconstructFvVolumeField
         );
     }
 
-    return reconstructFvVolumeField
-    (
-        IOobject
-        (
-            fieldIoObject.name(),
-            completeMesh_.time().name(),
-            completeMesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        procFields
-    );
-}
-
-
-template<class Type>
-Foam::tmp<Foam::SurfaceField<Type>>
-Foam::fvFieldReconstructor::reconstructFvSurfaceField
-(
-    const IOobject& fieldIoObject,
-    const PtrList<SurfaceField<Type>>& procFields
-) const
-{
     // Create the internalField
     Field<Type> internalField(completeMesh_.nInternalFaces());
 
@@ -420,7 +432,15 @@ Foam::fvFieldReconstructor::reconstructFvSurfaceField
     (
         new SurfaceField<Type>
         (
-            fieldIoObject,
+            IOobject
+            (
+                fieldIoObject.name(),
+                completeMesh_.time().name(),
+                completeMesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
             completeMesh_,
             procFields[0].dimensions(),
             internalField,
@@ -430,55 +450,10 @@ Foam::fvFieldReconstructor::reconstructFvSurfaceField
 }
 
 
-template<class Type>
-Foam::tmp<Foam::SurfaceField<Type>>
-Foam::fvFieldReconstructor::reconstructFvSurfaceField
-(
-    const IOobject& fieldIoObject
-) const
-{
-    PtrList<SurfaceField<Type>>
-        procFields(procMeshes_.size());
-
-    forAll(procMeshes_, proci)
-    {
-        procFields.set
-        (
-            proci,
-            new SurfaceField<Type>
-            (
-                IOobject
-                (
-                    fieldIoObject.name(),
-                    procMeshes_[proci].time().name(),
-                    procMeshes_[proci],
-                    IOobject::MUST_READ,
-                    IOobject::NO_WRITE,
-                    false
-                ),
-                procMeshes_[proci]
-            )
-        );
-    }
-
-    return reconstructFvSurfaceField
-    (
-        IOobject
-        (
-            fieldIoObject.name(),
-            completeMesh_.time().name(),
-            completeMesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        procFields
-    );
-}
-
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
-void Foam::fvFieldReconstructor::reconstructFvVolumeInternalFields
+void Foam::fvFieldReconstructor::reconstructVolInternalFields
 (
     const IOobjectList& objects,
     const HashSet<word>& selectedFields
@@ -490,7 +465,8 @@ void Foam::fvFieldReconstructor::reconstructFvVolumeInternalFields
 
     if (fields.size())
     {
-        Info<< "    Reconstructing " << fieldClassName << "s\n" << endl;
+        Info<< nl << "    Reconstructing " << fieldClassName << "s"
+            << nl << endl;
 
         forAllConstIter(IOobjectList, fields, fieldIter)
         {
@@ -502,18 +478,15 @@ void Foam::fvFieldReconstructor::reconstructFvVolumeInternalFields
             {
                 Info<< "        " << fieldIter()->name() << endl;
 
-                reconstructFvVolumeInternalField<Type>(*fieldIter())().write();
-
-                nReconstructed_++;
+                reconstructVolInternalField<Type>(*fieldIter())().write();
             }
         }
-        Info<< endl;
     }
 }
 
 
 template<class Type>
-void Foam::fvFieldReconstructor::reconstructFvVolumeFields
+void Foam::fvFieldReconstructor::reconstructVolFields
 (
     const IOobjectList& objects,
     const HashSet<word>& selectedFields
@@ -526,7 +499,8 @@ void Foam::fvFieldReconstructor::reconstructFvVolumeFields
 
     if (fields.size())
     {
-        Info<< "    Reconstructing " << fieldClassName << "s\n" << endl;
+        Info<< nl << "    Reconstructing " << fieldClassName << "s"
+            << nl << endl;
 
         forAllConstIter(IOobjectList, fields, fieldIter)
         {
@@ -538,12 +512,9 @@ void Foam::fvFieldReconstructor::reconstructFvVolumeFields
             {
                 Info<< "        " << fieldIter()->name() << endl;
 
-                reconstructFvVolumeField<Type>(*fieldIter())().write();
-
-                nReconstructed_++;
+                reconstructVolField<Type>(*fieldIter())().write();
             }
         }
-        Info<< endl;
     }
 }
 
@@ -562,7 +533,8 @@ void Foam::fvFieldReconstructor::reconstructFvSurfaceFields
 
     if (fields.size())
     {
-        Info<< "    Reconstructing " << fieldClassName << "s\n" << endl;
+        Info<< nl << "    Reconstructing " << fieldClassName << "s"
+            << nl << endl;
 
         forAllConstIter(IOobjectList, fields, fieldIter)
         {
@@ -575,11 +547,8 @@ void Foam::fvFieldReconstructor::reconstructFvSurfaceFields
                 Info<< "        " << fieldIter()->name() << endl;
 
                 reconstructFvSurfaceField<Type>(*fieldIter())().write();
-
-                nReconstructed_++;
             }
         }
-        Info<< endl;
     }
 }
 
diff --git a/applications/utilities/parallelProcessing/reconstructPar/lagrangianFieldReconstructor.C b/applications/utilities/parallelProcessing/reconstructPar/lagrangianFieldReconstructor.C
new file mode 100644
index 0000000000..d41d7bdde4
--- /dev/null
+++ b/applications/utilities/parallelProcessing/reconstructPar/lagrangianFieldReconstructor.C
@@ -0,0 +1,124 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "lagrangianFieldReconstructor.H"
+#include "passiveParticleCloud.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::lagrangianFieldReconstructor::lagrangianFieldReconstructor
+(
+    const fvMesh& completeMesh,
+    const PtrList<fvMesh>& procMeshes,
+    const labelListList& faceProcAddressing,
+    const labelListList& cellProcAddressing,
+    const word& cloudName
+)
+:
+    completeMesh_(completeMesh),
+    procMeshes_(procMeshes),
+    cloudName_(cloudName)
+{
+    // Construct and empty cloud for the complete positions
+    passiveParticleCloud completePositions
+    (
+        completeMesh_,
+        cloudName_,
+        IDLList<passiveParticle>()
+    );
+
+    forAll(procMeshes_, proci)
+    {
+        // Read the processor positions
+        Cloud<passiveParticle> procPositions
+        (
+            procMeshes_[proci],
+            cloudName_,
+            false
+        );
+
+        // Combine the processor's positions into the complete cloud
+        forAllConstIter(Cloud<passiveParticle>, procPositions, iter)
+        {
+            const passiveParticle& p = iter();
+            const label completeCelli = cellProcAddressing[proci][p.cell()];
+            const label completeFacei =
+                mag(faceProcAddressing[proci][p.tetFace()]) - 1;
+
+            completePositions.append
+            (
+                new passiveParticle
+                (
+                    completeMesh_,
+                    p.coordinates(),
+                    completeCelli,
+                    completeFacei,
+                    p.procTetPt
+                    (
+                        procMeshes_[proci],
+                        completeMesh_,
+                        completeCelli,
+                        completeFacei
+                    )
+                )
+            );
+        }
+    }
+
+    // Write
+    IOPosition<Cloud<passiveParticle>>(completePositions).write();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::lagrangianFieldReconstructor::~lagrangianFieldReconstructor()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::lagrangianFieldReconstructor::reconstructs
+(
+    const IOobjectList& objects,
+    const HashSet<word>& selectedFields
+)
+{
+    bool result = false;
+
+    #define DO_LAGRANGIAN_FIELDS_TYPE(Type, nullArg)                           \
+        result = result                                                        \
+         || reconstructs<IOField<Type>>(objects, selectedFields)               \
+         || reconstructs<IOField<Field<Type>>>(objects, selectedFields)        \
+         || reconstructs<CompactIOField<Field<Type>>>(objects, selectedFields);
+    DO_LAGRANGIAN_FIELDS_TYPE(label, )
+    FOR_ALL_FIELD_TYPES(DO_LAGRANGIAN_FIELDS_TYPE)
+    #undef DO_LAGRANGIAN_FIELDS_TYPE
+
+    return result;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/parallelProcessing/reconstructPar/lagrangianFieldReconstructor.H b/applications/utilities/parallelProcessing/reconstructPar/lagrangianFieldReconstructor.H
new file mode 100644
index 0000000000..21d74e6268
--- /dev/null
+++ b/applications/utilities/parallelProcessing/reconstructPar/lagrangianFieldReconstructor.H
@@ -0,0 +1,165 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::lagrangianFieldReconstructor
+
+Description
+    Lagrangian field reconstructor.
+
+SourceFiles
+    lagrangianFieldReconstructor.C
+    lagrangianFieldReconstructorTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef lagrangianFieldReconstructor_H
+#define lagrangianFieldReconstructor_H
+
+#include "cloud.H"
+#include "fvMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class IOobjectList;
+
+/*---------------------------------------------------------------------------*\
+                  Class lagrangianFieldReconstructor Declaration
+\*---------------------------------------------------------------------------*/
+
+class lagrangianFieldReconstructor
+{
+    // Private Data
+
+        //- Reference to complete mesh
+        const fvMesh& completeMesh_;
+
+        //- List of processor meshes
+        const PtrList<fvMesh>& procMeshes_;
+
+        //- The name of the cloud
+        const word cloudName_;
+
+
+    // Private Member Functions
+
+        //- Return whether anything in the object list gets reconstructed
+        template<class FieldType>
+        static bool reconstructs
+        (
+            const IOobjectList& objects,
+            const HashSet<word>& selectedFields
+        );
+
+        //- Reconstruct a field
+        template
+        <
+            class Type,
+            template<class> class IOContainer,
+            template<class> class IOContainerType
+        >
+        tmp<IOContainer<Type>>
+        reconstructField(const IOobject& fieldIoObject) const;
+
+        //- Read, reconstruct and write all fields
+        template
+        <
+            class Type,
+            template<class> class IOContainer,
+            template<class> class IOContainerType
+        >
+        void reconstructFields
+        (
+            const IOobjectList& objects,
+            const HashSet<word>& selectedFields
+        );
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        lagrangianFieldReconstructor
+        (
+            const fvMesh& completeMesh,
+            const PtrList<fvMesh>& procMeshes,
+            const labelListList& faceProcAddressing,
+            const labelListList& cellProcAddressing,
+            const word& cloudName
+        );
+
+        //- Disallow default bitwise copy construction
+        lagrangianFieldReconstructor
+        (
+            const lagrangianFieldReconstructor&
+        ) = delete;
+
+
+    //- Destructor
+    ~lagrangianFieldReconstructor();
+
+
+    // Member Functions
+
+        //- Return whether anything in the object list gets reconstructed
+        static bool reconstructs
+        (
+            const IOobjectList& objects,
+            const HashSet<word>& selectedFields
+        );
+
+        //- Read, reconstruct and write all fields
+        template<class Type>
+        void reconstructFields
+        (
+            const IOobjectList& objects,
+            const HashSet<word>& selectedFields
+        );
+
+
+    // Member Operators
+
+        //- Disallow default bitwise assignment
+        void operator=(const lagrangianFieldReconstructor&) = delete;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "lagrangianFieldReconstructorTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/parallelProcessing/reconstructPar/lagrangianFieldReconstructorTemplates.C b/applications/utilities/parallelProcessing/reconstructPar/lagrangianFieldReconstructorTemplates.C
new file mode 100644
index 0000000000..5d485f493a
--- /dev/null
+++ b/applications/utilities/parallelProcessing/reconstructPar/lagrangianFieldReconstructorTemplates.C
@@ -0,0 +1,179 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "lagrangianFieldReconstructor.H"
+#include "IOobjectList.H"
+#include "CompactIOField.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class FieldType>
+bool Foam::lagrangianFieldReconstructor::reconstructs
+(
+    const IOobjectList& objects,
+    const HashSet<word>& selectedFields
+)
+{
+    IOobjectList fields = objects.lookupClass(FieldType::typeName);
+
+    if (fields.size() && selectedFields.empty())
+    {
+        return true;
+    }
+
+    forAllConstIter(IOobjectList, fields, fieldIter)
+    {
+        if (selectedFields.found(fieldIter()->name()))
+        {
+            return true;
+        }
+    }
+
+    return false;
+}
+
+
+template
+<
+    class Type,
+    template<class> class IOContainer,
+    template<class> class IOContainerType
+>
+Foam::tmp<IOContainer<Type>>
+Foam::lagrangianFieldReconstructor::reconstructField
+(
+    const IOobject& fieldIoObject
+) const
+{
+    // Construct the complete field
+    tmp<IOContainer<Type>> tfield
+    (
+        new IOContainer<Type>
+        (
+            IOobject
+            (
+                fieldIoObject.name(),
+                completeMesh_.time().name(),
+                cloud::prefix/cloudName_,
+                completeMesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            Field<Type>(0)
+        )
+    );
+    Field<Type>& field = tfield.ref();
+
+    // Combine the processor fields into the complete field
+    forAll(procMeshes_, proci)
+    {
+        typeIOobject<IOContainerType<Type>> localIOobject
+        (
+            fieldIoObject.name(),
+            procMeshes_[proci].time().name(),
+            cloud::prefix/cloudName_,
+            procMeshes_[proci],
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        );
+
+        if (localIOobject.headerOk())
+        {
+            IOContainer<Type> fieldi(localIOobject);
+
+            label offset = field.size();
+            field.setSize(offset + fieldi.size());
+
+            forAll(fieldi, j)
+            {
+                field[offset + j] = fieldi[j];
+            }
+        }
+    }
+
+    return tfield;
+}
+
+
+template
+<
+    class Type,
+    template<class> class IOContainer,
+    template<class> class IOContainerType
+>
+void Foam::lagrangianFieldReconstructor::reconstructFields
+(
+    const IOobjectList& objects,
+    const HashSet<word>& selectedFields
+)
+{
+    const word& fieldClassName = IOContainerType<Type>::typeName;
+
+    IOobjectList fields = objects.lookupClass(fieldClassName);
+
+    if (fields.size())
+    {
+        Info<< nl << "    Reconstructing lagrangian "
+            << fieldClassName << "s" << nl << endl;
+
+        forAllConstIter(IOobjectList, fields, fieldIter)
+        {
+            if
+            (
+                selectedFields.empty()
+             || selectedFields.found(fieldIter()->name())
+            )
+            {
+                Info<< "        " << fieldIter()->name() << endl;
+
+                reconstructField<Type, IOContainer, IOContainerType>
+                (
+                    *fieldIter()
+                )().write();
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::lagrangianFieldReconstructor::reconstructFields
+(
+    const IOobjectList& objects,
+    const HashSet<word>& selectedFields
+)
+{
+    reconstructFields<Type, IOField, IOField>
+    (objects, selectedFields);
+    reconstructFields<Field<Type>, CompactIOField, IOField>
+    (objects, selectedFields);
+    reconstructFields<Field<Type>, CompactIOField, CompactIOField>
+    (objects, selectedFields);
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/parallelProcessing/reconstructPar/pointFieldReconstructor.C b/applications/utilities/parallelProcessing/reconstructPar/pointFieldReconstructor.C
index f81b562c3e..b66f7ed9a3 100644
--- a/applications/utilities/parallelProcessing/reconstructPar/pointFieldReconstructor.C
+++ b/applications/utilities/parallelProcessing/reconstructPar/pointFieldReconstructor.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -38,8 +38,7 @@ Foam::pointFieldReconstructor::pointFieldReconstructor
     completeMesh_(completeMesh),
     procMeshes_(procMeshes),
     pointProcAddressing_(pointProcAddressing),
-    patchPointAddressing_(procMeshes.size()),
-    nReconstructed_(0)
+    patchPointAddressing_(procMeshes.size())
 {
     // Inverse-addressing of the patch point labels.
     labelList pointMap(completeMesh_.size(), -1);
@@ -91,4 +90,24 @@ Foam::pointFieldReconstructor::pointFieldReconstructor
 }
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::pointFieldReconstructor::reconstructs
+(
+    const IOobjectList& objects,
+    const HashSet<word>& selectedFields
+)
+{
+    bool result = false;
+
+    #define DO_POINT_FIELDS_TYPE(Type, nullArg)                                \
+        result = result                                                        \
+         || reconstructs<PointField<Type>>(objects, selectedFields);
+    FOR_ALL_FIELD_TYPES(DO_POINT_FIELDS_TYPE)
+    #undef DO_POINT_FIELDS_TYPE
+
+    return result;
+}
+
+
 // ************************************************************************* //
diff --git a/applications/utilities/parallelProcessing/reconstructPar/pointFieldReconstructor.H b/applications/utilities/parallelProcessing/reconstructPar/pointFieldReconstructor.H
index 6f3d282282..479a59afa4 100644
--- a/applications/utilities/parallelProcessing/reconstructPar/pointFieldReconstructor.H
+++ b/applications/utilities/parallelProcessing/reconstructPar/pointFieldReconstructor.H
@@ -29,6 +29,7 @@ Description
 
 SourceFiles
     pointFieldReconstructor.C
+    pointFieldReconstructorTemplates.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -68,8 +69,21 @@ class pointFieldReconstructor
         //- Point patch addressing
         labelListListList patchPointAddressing_;
 
-        //- Number of fields reconstructed
-        label nReconstructed_;
+
+    // Private Member Functions
+
+        //- Return whether anything in the object list gets reconstructed
+        template<class FieldType>
+        static bool reconstructs
+        (
+            const IOobjectList& objects,
+            const HashSet<word>& selectedFields
+        );
+
+        //- Reconstruct field
+        template<class Type>
+        tmp<PointField<Type>>
+        reconstructField(const IOobject& fieldIoObject);
 
 
 public:
@@ -90,16 +104,12 @@ public:
 
     // Member Functions
 
-        //- Return number of fields reconstructed
-        label nReconstructed() const
-        {
-            return nReconstructed_;
-        }
-
-        //- Reconstruct field
-        template<class Type>
-        tmp<PointField<Type>>
-        reconstructField(const IOobject& fieldIoObject);
+        //- Return whether anything in the object list gets reconstructed
+        static bool reconstructs
+        (
+            const IOobjectList& objects,
+            const HashSet<word>& selectedFields
+        );
 
         //- Reconstruct and write all fields
         template<class Type>
@@ -124,7 +134,7 @@ public:
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-    #include "pointFieldReconstructorReconstructFields.C"
+    #include "pointFieldReconstructorTemplates.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/applications/utilities/parallelProcessing/reconstructPar/pointFieldReconstructorReconstructFields.C b/applications/utilities/parallelProcessing/reconstructPar/pointFieldReconstructorTemplates.C
similarity index 87%
rename from applications/utilities/parallelProcessing/reconstructPar/pointFieldReconstructorReconstructFields.C
rename to applications/utilities/parallelProcessing/reconstructPar/pointFieldReconstructorTemplates.C
index 7d5bd3a87e..e5e1c583f2 100644
--- a/applications/utilities/parallelProcessing/reconstructPar/pointFieldReconstructorReconstructFields.C
+++ b/applications/utilities/parallelProcessing/reconstructPar/pointFieldReconstructorTemplates.C
@@ -27,17 +27,40 @@ License
 #include "fvMesh.H"
 #include "reversePointPatchFieldMapper.H"
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class FieldType>
+bool Foam::pointFieldReconstructor::reconstructs
+(
+    const IOobjectList& objects,
+    const HashSet<word>& selectedFields
+)
+{
+    IOobjectList fields = objects.lookupClass(FieldType::typeName);
+
+    if (fields.size() && selectedFields.empty())
+    {
+        return true;
+    }
+
+    forAllConstIter(IOobjectList, fields, fieldIter)
+    {
+        if (selectedFields.found(fieldIter()->name()))
+        {
+            return true;
+        }
+    }
+
+    return false;
+}
+
 
 template<class Type>
 Foam::tmp<Foam::PointField<Type>>
 Foam::pointFieldReconstructor::reconstructField(const IOobject& fieldIoObject)
 {
     // Read the field for all the processors
-    PtrList<PointField<Type>> procFields
-    (
-        procMeshes_.size()
-    );
+    PtrList<PointField<Type>> procFields(procMeshes_.size());
 
     forAll(procMeshes_, proci)
     {
@@ -59,14 +82,12 @@ Foam::pointFieldReconstructor::reconstructField(const IOobject& fieldIoObject)
         );
     }
 
-
     // Create the internalField
     Field<Type> internalField(completeMesh_.size());
 
     // Create the patch fields
     PtrList<pointPatchField<Type>> patchFields(completeMesh_.boundary().size());
 
-
     forAll(procMeshes_, proci)
     {
         const PointField<Type>&
@@ -122,8 +143,7 @@ Foam::pointFieldReconstructor::reconstructField(const IOobject& fieldIoObject)
         }
     }
 
-    // Construct and write the field
-    // setting the internalField and patchFields
+    // Construct and return the field
     return tmp<PointField<Type>>
     (
         new PointField<Type>
@@ -145,7 +165,8 @@ Foam::pointFieldReconstructor::reconstructField(const IOobject& fieldIoObject)
 }
 
 
-// Reconstruct and write all point fields
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
 template<class Type>
 void Foam::pointFieldReconstructor::reconstructFields
 (
@@ -162,7 +183,8 @@ void Foam::pointFieldReconstructor::reconstructFields
 
     if (fields.size())
     {
-        Info<< "    Reconstructing " << fieldClassName << "s\n" << endl;
+        Info<< nl << "    Reconstructing " << fieldClassName << "s"
+            << nl << endl;
 
         forAllConstIter(IOobjectList, fields, fieldIter)
         {
@@ -175,12 +197,8 @@ void Foam::pointFieldReconstructor::reconstructFields
                 Info<< "        " << fieldIter()->name() << endl;
 
                 reconstructField<Type>(*fieldIter())().write();
-
-                nReconstructed_++;
             }
         }
-
-        Info<< endl;
     }
 }
 
diff --git a/applications/utilities/parallelProcessing/reconstructPar/reconstructLagrangian.H b/applications/utilities/parallelProcessing/reconstructPar/reconstructLagrangian.H
deleted file mode 100644
index 4072dce2f7..0000000000
--- a/applications/utilities/parallelProcessing/reconstructPar/reconstructLagrangian.H
+++ /dev/null
@@ -1,117 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-InClass
-    Foam::reconstructLagrangian
-
-Description
-
-SourceFiles
-    reconstructLagrangianPositions.C
-    reconstructLagrangianFields.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef reconstructLagrangian_H
-#define reconstructLagrangian_H
-
-#include "cloud.H"
-#include "polyMesh.H"
-#include "IOobjectList.H"
-#include "CompactIOField.H"
-#include "fvMesh.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-void reconstructLagrangianPositions
-(
-    const polyMesh& mesh,
-    const word& cloudName,
-    const PtrList<fvMesh>& meshes,
-    const labelListList& faceProcAddressing,
-    const labelListList& cellProcAddressing
-);
-
-
-template<class Type>
-tmp<IOField<Type>> reconstructLagrangianField
-(
-    const word& cloudName,
-    const polyMesh& mesh,
-    const PtrList<fvMesh>& meshes,
-    const word& fieldName
-);
-
-
-template<class Type>
-tmp<CompactIOField<Field<Type>>> reconstructLagrangianFieldField
-(
-    const word& cloudName,
-    const polyMesh& mesh,
-    const PtrList<fvMesh>& meshes,
-    const word& fieldName
-);
-
-
-template<class Type>
-void reconstructLagrangianFields
-(
-    const word& cloudName,
-    const polyMesh& mesh,
-    const PtrList<fvMesh>& meshes,
-    const IOobjectList& objects,
-    const HashSet<word>& selectedFields
-);
-
-
-template<class Type>
-void reconstructLagrangianFieldFields
-(
-    const word& cloudName,
-    const polyMesh& mesh,
-    const PtrList<fvMesh>& meshes,
-    const IOobjectList& objects,
-    const HashSet<word>& selectedFields
-);
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-    #include "reconstructLagrangianFields.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/applications/utilities/parallelProcessing/reconstructPar/reconstructLagrangianFields.C b/applications/utilities/parallelProcessing/reconstructPar/reconstructLagrangianFields.C
deleted file mode 100644
index b50ff39442..0000000000
--- a/applications/utilities/parallelProcessing/reconstructPar/reconstructLagrangianFields.C
+++ /dev/null
@@ -1,273 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "IOField.H"
-#include "CompactIOField.H"
-#include "Time.H"
-
-// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
-
-template<class Type>
-Foam::tmp<Foam::IOField<Type>> Foam::reconstructLagrangianField
-(
-    const word& cloudName,
-    const polyMesh& mesh,
-    const PtrList<fvMesh>& meshes,
-    const word& fieldName
-)
-{
-    // Construct empty field on mesh
-    tmp<IOField<Type>> tfield
-    (
-        new IOField<Type>
-        (
-            IOobject
-            (
-                fieldName,
-                mesh.time().name(),
-                cloud::prefix/cloudName,
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            Field<Type>(0)
-        )
-    );
-    Field<Type>& field = tfield.ref();
-
-    forAll(meshes, i)
-    {
-        // Check object on local mesh
-        typeIOobject<IOField<Type>> localIOobject
-        (
-            fieldName,
-            meshes[i].time().name(),
-            cloud::prefix/cloudName,
-            meshes[i],
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        );
-
-        if (localIOobject.headerOk())
-        {
-            IOField<Type> fieldi(localIOobject);
-
-            label offset = field.size();
-            field.setSize(offset + fieldi.size());
-
-            forAll(fieldi, j)
-            {
-                field[offset + j] = fieldi[j];
-            }
-        }
-    }
-
-    return tfield;
-}
-
-
-template<class Type>
-Foam::tmp<Foam::CompactIOField<Foam::Field<Type>>>
-Foam::reconstructLagrangianFieldField
-(
-    const word& cloudName,
-    const polyMesh& mesh,
-    const PtrList<fvMesh>& meshes,
-    const word& fieldName
-)
-{
-    // Construct empty field on mesh
-    tmp<CompactIOField<Field<Type>>> tfield
-    (
-        new CompactIOField<Field<Type>>
-        (
-            IOobject
-            (
-                fieldName,
-                mesh.time().name(),
-                cloud::prefix/cloudName,
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            Field<Field<Type>>(0)
-        )
-    );
-    Field<Field<Type>>& field = tfield.ref();
-
-    forAll(meshes, i)
-    {
-        // No type checking is done to handle CompactIOField and IOField
-        IOobject localIOobject
-        (
-            fieldName,
-            meshes[i].time().name(),
-            cloud::prefix/cloudName,
-            meshes[i],
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        );
-
-        if (localIOobject.headerOk())
-        {
-            CompactIOField<Field<Type>> fieldi(localIOobject);
-
-            label offset = field.size();
-            field.setSize(offset + fieldi.size());
-
-            forAll(fieldi, j)
-            {
-                field[offset + j] = fieldi[j];
-            }
-        }
-    }
-
-    return tfield;
-}
-
-
-
-template<class Type>
-void Foam::reconstructLagrangianFields
-(
-    const word& cloudName,
-    const polyMesh& mesh,
-    const PtrList<fvMesh>& meshes,
-    const IOobjectList& objects,
-    const HashSet<word>& selectedFields
-)
-{
-    const word fieldClassName(IOField<Type>::typeName);
-
-    IOobjectList fields = objects.lookupClass(fieldClassName);
-
-    if (fields.size())
-    {
-        Info<< "    Reconstructing lagrangian "
-            << fieldClassName << "s\n" << endl;
-
-        forAllConstIter(IOobjectList, fields, fieldIter)
-        {
-            if
-            (
-                selectedFields.empty()
-             || selectedFields.found(fieldIter()->name())
-            )
-            {
-                Info<< "        " << fieldIter()->name() << endl;
-                reconstructLagrangianField<Type>
-                (
-                    cloudName,
-                    mesh,
-                    meshes,
-                    fieldIter()->name()
-                )().write();
-            }
-        }
-
-        Info<< endl;
-    }
-}
-
-
-template<class Type>
-void Foam::reconstructLagrangianFieldFields
-(
-    const word& cloudName,
-    const polyMesh& mesh,
-    const PtrList<fvMesh>& meshes,
-    const IOobjectList& objects,
-    const HashSet<word>& selectedFields
-)
-{
-    {
-        const word fieldClassName(CompactIOField<Field<Type>>::typeName);
-
-        IOobjectList fields = objects.lookupClass(fieldClassName);
-
-        if (fields.size())
-        {
-            Info<< "    Reconstructing lagrangian "
-                << fieldClassName << "s\n" << endl;
-
-            forAllConstIter(IOobjectList, fields, fieldIter)
-            {
-                if
-                (
-                    selectedFields.empty()
-                 || selectedFields.found(fieldIter()->name())
-                )
-                {
-                    Info<< "        " << fieldIter()->name() << endl;
-                    reconstructLagrangianFieldField<Type>
-                    (
-                        cloudName,
-                        mesh,
-                        meshes,
-                        fieldIter()->name()
-                    )().write();
-                }
-            }
-
-            Info<< endl;
-        }
-    }
-
-    {
-        const word fieldClassName(IOField<Field<Type>>::typeName);
-
-        IOobjectList fields = objects.lookupClass(fieldClassName);
-
-        if (fields.size())
-        {
-            Info<< "    Reconstructing lagrangian "
-                << fieldClassName << "s\n" << endl;
-
-            forAllConstIter(IOobjectList, fields, fieldIter)
-            {
-                if
-                (
-                    selectedFields.empty()
-                 || selectedFields.found(fieldIter()->name())
-                )
-                {
-                    Info<< "        " << fieldIter()->name() << endl;
-                    reconstructLagrangianFieldField<Type>
-                    (
-                        cloudName,
-                        mesh,
-                        meshes,
-                        fieldIter()->name()
-                    )().write();
-                }
-            }
-
-            Info<< endl;
-        }
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/applications/utilities/parallelProcessing/reconstructPar/reconstructLagrangianPositions.C b/applications/utilities/parallelProcessing/reconstructPar/reconstructLagrangianPositions.C
deleted file mode 100644
index 283aba65e4..0000000000
--- a/applications/utilities/parallelProcessing/reconstructPar/reconstructLagrangianPositions.C
+++ /dev/null
@@ -1,83 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "reconstructLagrangian.H"
-#include "labelIOList.H"
-#include "passiveParticleCloud.H"
-
-// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
-
-void Foam::reconstructLagrangianPositions
-(
-    const polyMesh& mesh,
-    const word& cloudName,
-    const PtrList<fvMesh>& meshes,
-    const labelListList& faceProcAddressing,
-    const labelListList& cellProcAddressing
-)
-{
-    passiveParticleCloud lagrangianPositions
-    (
-        mesh,
-        cloudName,
-        IDLList<passiveParticle>()
-    );
-
-    forAll(meshes, i)
-    {
-        const labelList& cellMap = cellProcAddressing[i];
-        const labelList& faceMap = faceProcAddressing[i];
-
-        Cloud<passiveParticle> lpi(meshes[i], cloudName, false);
-
-        forAllConstIter(Cloud<passiveParticle>, lpi, iter)
-        {
-            const passiveParticle& ppi = iter();
-
-            const label mappedCell = cellMap[ppi.cell()];
-
-            // Inverting sign if necessary and subtracting 1 from
-            // faceProcAddressing
-            label mappedTetFace = mag(faceMap[ppi.tetFace()]) - 1;
-
-            lagrangianPositions.append
-            (
-                new passiveParticle
-                (
-                    mesh,
-                    ppi.coordinates(),
-                    mappedCell,
-                    mappedTetFace,
-                    ppi.procTetPt(meshes[i], mesh, mappedCell, mappedTetFace)
-                )
-            );
-        }
-    }
-
-    IOPosition<Cloud<passiveParticle>>(lagrangianPositions).write();
-}
-
-
-// ************************************************************************* //
diff --git a/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C b/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C
index 903602522e..2b46fd0005 100644
--- a/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C
+++ b/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C
@@ -34,10 +34,10 @@ Description
 #include "timeSelector.H"
 #include "IOobjectList.H"
 #include "processorRunTimes.H"
-#include "domainDecomposition.H"
+#include "multiDomainDecomposition.H"
 #include "fvFieldReconstructor.H"
 #include "pointFieldReconstructor.H"
-#include "reconstructLagrangian.H"
+#include "lagrangianFieldReconstructor.H"
 
 using namespace Foam;
 
@@ -46,21 +46,37 @@ using namespace Foam;
 namespace Foam
 {
 
-bool haveAllTimes
+bool haveUniform
 (
-    const HashSet<word>& masterTimeDirSet,
-    const instantList& timeDirs
+    const processorRunTimes& runTimes,
+    const word& regionDir = word::null
 )
 {
-    // Loop over all times
-    forAll(timeDirs, timei)
-    {
-        if (!masterTimeDirSet.found(timeDirs[timei].name()))
-        {
-            return false;
-        }
-    }
-    return true;
+    return
+        fileHandler().isDir
+        (
+            fileHandler().filePath
+            (
+                runTimes.procTimes()[0].timePath()/regionDir/"uniform"
+            )
+        );
+}
+
+
+void reconstructUniform
+(
+    const processorRunTimes& runTimes,
+    const word& regionDir = word::null
+)
+{
+    fileHandler().cp
+    (
+        fileHandler().filePath
+        (
+            runTimes.procTimes()[0].timePath()/regionDir/"uniform"
+        ),
+        runTimes.completeTime().timePath()/regionDir
+    );
 }
 
 
@@ -85,14 +101,41 @@ void writeDecomposition(const domainDecomposition& meshes)
     cellProc.write();
 
     Info<< "Wrote decomposition as volScalarField::Internal to "
-        << cellProc.name() << " for use in postprocessing."
-        << nl << endl;
+        << cellProc.name() << " for use in postprocessing"
+        << endl;
 }
 
 
 }
 
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class delayedNewLine
+{
+    mutable bool first_;
+
+public:
+
+    delayedNewLine()
+    :
+        first_(true)
+    {}
+
+    friend Ostream& operator<<(Ostream& os, const delayedNewLine& dnl)
+    {
+        if (!dnl.first_) os << nl;
+        dnl.first_ = false;
+        return os;
+    }
+};
+
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
@@ -102,9 +145,6 @@ int main(int argc, char *argv[])
         "Reconstruct fields of a parallel case"
     );
 
-    // Enable -constant ... if someone really wants it
-    // Enable -withZero to prevent accidentally trashing the initial fields
-    timeSelector::addOptions(true, true);
     argList::noParallel();
     #include "addRegionOption.H"
     #include "addAllRegionsOption.H"
@@ -112,7 +152,7 @@ int main(int argc, char *argv[])
     (
         "cellProc",
         "write cell processor indices as a volScalarField::Internal for "
-        "post-processing."
+        "post-processing"
     );
     argList::addOption
     (
@@ -132,7 +172,7 @@ int main(int argc, char *argv[])
         "list",
         "specify a list of lagrangian fields to be reconstructed. Eg, '(U d)' -"
         "regular expressions not currently supported, "
-        "positions always included."
+        "positions always included"
     );
     argList::addBoolOption
     (
@@ -150,6 +190,10 @@ int main(int argc, char *argv[])
         "only reconstruct new times (i.e. that do not exist already)"
     );
 
+    // Include explicit constant options, and explicit zero option (to prevent
+    // the user accidentally trashing the initial fields)
+    timeSelector::addOptions(true, true);
+
     #include "setRootCase.H"
 
     const bool writeCellProc = args.optionFound("cellProc");
@@ -164,8 +208,7 @@ int main(int argc, char *argv[])
 
     if (noFields)
     {
-        Info<< "Skipping reconstructing fields"
-            << nl << endl;
+        Info<< "Skipping reconstructing fields" << nl << endl;
     }
 
     const bool noLagrangian = args.optionFound("noLagrangian");
@@ -191,7 +234,7 @@ int main(int argc, char *argv[])
         {
             FatalErrorInFunction
                 << "Cannot specify noLagrangian and lagrangianFields "
-                << "options together."
+                << "options together"
                 << exit(FatalError);
         }
 
@@ -199,11 +242,11 @@ int main(int argc, char *argv[])
     }
 
     // Set time from database
-    Info<< "Create time\n" << endl;
+    Info<< "Create time" << nl << endl;
     processorRunTimes runTimes(Foam::Time::controlDictName, args);
 
-    // Allow override of time
-    const instantList times = runTimes.selectProc(args);
+    // Get the times to reconstruct
+    instantList times = runTimes.selectProc(args);
 
     const Time& runTime = runTimes.procTimes()[0];
 
@@ -228,484 +271,335 @@ int main(int argc, char *argv[])
     // Warn fileHandler of number of processors
     const_cast<fileOperation&>(fileHandler()).setNProcs(nProcs);
 
-    // Note that we do not set the runTime time so it is still the
-    // one set through the controlDict. The -time option
-    // only affects the selected set of times from processor0.
-    // - can be illogical
-    // + any point motion handled through mesh.readUpdate
+    // Quit if no times
     if (times.empty())
     {
-        WarningInFunction << "No times selected" << endl;
+        WarningInFunction << "No times selected" << nl << endl;
         exit(1);
     }
 
-    // Get current times if -newTimes
-    const bool newTimes   = args.optionFound("newTimes");
-    instantList masterTimeDirs;
-    if (newTimes)
+    // If only reconstructing new times then filter out existing times
+    if (args.optionFound("newTimes"))
     {
-        masterTimeDirs = runTimes.completeTime().times();
+        // Get all existing times
+        const instantList existingTimes = runTimes.completeTime().times();
+
+        // Put into a set
+        HashSet<word> existingTimesSet;
+        existingTimesSet.resize(2*existingTimes.size());
+        forAll(existingTimes, i)
+        {
+            existingTimesSet.insert(existingTimes[i].name());
+        }
+
+        // Remove times from the existing time set by shuffling up
+        label timei = 0;
+        forAll(times, timej)
+        {
+            if (!existingTimesSet.found(times[timej].name()))
+            {
+                times[timei ++] = times[timej];
+            }
+        }
+        times.resize(timei);
     }
-    HashSet<word> masterTimeDirSet(2*masterTimeDirs.size());
-    forAll(masterTimeDirs, i)
+
+    // Quit if no times
+    if (times.empty())
     {
-        masterTimeDirSet.insert(masterTimeDirs[i].name());
-    }
-    if
-    (
-        newTimes
-     && regionNames.size() == 1
-     && regionNames[0] == fvMesh::defaultRegion
-     && haveAllTimes(masterTimeDirSet, times)
-    )
-    {
-        Info<< "All times already reconstructed.\n\nEnd\n" << endl;
+        Info<< "All times already reconstructed" << nl << nl
+            << "End" << nl << endl;
         return 0;
     }
 
-    // Reconstruct all regions
-    forAll(regionNames, regioni)
+    // Create meshes
+    multiDomainDecomposition regionMeshes(runTimes, regionNames);
+    if (regionMeshes.readReconstruct(!noReconstructSets))
     {
-        const word& regionName = regionNames[regioni];
+        Info<< endl;
 
-        const word& regionDir =
-            regionName == polyMesh::defaultRegion
-          ? word::null
-          : regionName;
-
-        // Create meshes
-        Info<< "\n\nReconstructing mesh " << regionName << nl << endl;
-        domainDecomposition meshes(runTimes, regionName);
-        if (meshes.readReconstruct(!noReconstructSets) && writeCellProc)
+        if (writeCellProc)
         {
-            writeDecomposition(meshes);
-            fileHandler().flush();
-        }
-
-        // Loop over all times
-        forAll(times, timei)
-        {
-            if (newTimes && masterTimeDirSet.found(times[timei].name()))
+            forAll(regionNames, regioni)
             {
-                Info<< "Skipping time " << times[timei].name()
-                    << endl << endl;
-                continue;
-            }
-
-            // Set the time
-            runTimes.setTime(times[timei], timei);
-
-            Info<< "Time = " << runTimes.completeTime().userTimeName()
-                << nl << endl;
-
-            // Update the meshes
-            const fvMesh::readUpdateState state =
-                meshes.readUpdateReconstruct();
-
-            // Write the mesh out, if necessary
-            if (state != fvMesh::UNCHANGED)
-            {
-                meshes.writeComplete(!noReconstructSets);
-            }
-
-            // Write the decomposition, if necessary
-            if
-            (
-                writeCellProc
-             && meshes.completeMesh().facesInstance()
-             == runTimes.completeTime().name()
-            )
-            {
-                writeDecomposition(meshes);
+                writeDecomposition(regionMeshes.meshes(regioni)());
+                Info<< endl;
                 fileHandler().flush();
             }
+        }
+    }
 
-            // Get list of objects from processor0 database
-            IOobjectList objects
-            (
-                meshes.procMeshes()[0],
-                runTimes.procTimes()[0].name()
-            );
+    // Loop over all times
+    forAll(times, timei)
+    {
+        // Set the time
+        runTimes.setTime(times[timei], timei);
 
-            if (!noFields)
+        Info<< "Time = " << runTimes.completeTime().userTimeName()
+            << nl << endl;
+
+        // Update the meshes
+        const fvMesh::readUpdateState stat =
+            regionMeshes.readUpdateReconstruct();
+        if (stat >= fvMesh::TOPO_CHANGE) Info<< endl;
+
+        // Write the mesh out (if anything has changed)
+        regionMeshes.writeComplete(!noReconstructSets);
+
+        // Write the decomposition, if necessary
+        forAll(regionNames, regioni)
+        {
+            if (writeCellProc && stat >= fvMesh::TOPO_CHANGE)
             {
-                // If there are any FV fields, reconstruct them
-                Info<< "Reconstructing FV fields" << nl << endl;
-
-                fvFieldReconstructor fvReconstructor
-                (
-                    meshes.completeMesh(),
-                    meshes.procMeshes(),
-                    meshes.procFaceAddressing(),
-                    meshes.procCellAddressing(),
-                    meshes.procFaceAddressingBf()
-                );
-
-                fvReconstructor.reconstructFvVolumeInternalFields<scalar>
-                (
-                    objects,
-                    selectedFields
-                );
-                fvReconstructor.reconstructFvVolumeInternalFields<vector>
-                (
-                    objects,
-                    selectedFields
-                );
-                fvReconstructor.reconstructFvVolumeInternalFields
-                <sphericalTensor>
-                (
-                    objects,
-                    selectedFields
-                );
-                fvReconstructor.reconstructFvVolumeInternalFields<symmTensor>
-                (
-                    objects,
-                    selectedFields
-                );
-                fvReconstructor.reconstructFvVolumeInternalFields<tensor>
-                (
-                    objects,
-                    selectedFields
-                );
-
-                fvReconstructor.reconstructFvVolumeFields<scalar>
-                (
-                    objects,
-                    selectedFields
-                );
-                fvReconstructor.reconstructFvVolumeFields<vector>
-                (
-                    objects,
-                    selectedFields
-                );
-                fvReconstructor.reconstructFvVolumeFields<sphericalTensor>
-                (
-                    objects,
-                    selectedFields
-                );
-                fvReconstructor.reconstructFvVolumeFields<symmTensor>
-                (
-                    objects,
-                    selectedFields
-                );
-                fvReconstructor.reconstructFvVolumeFields<tensor>
-                (
-                    objects,
-                    selectedFields
-                );
-
-                fvReconstructor.reconstructFvSurfaceFields<scalar>
-                (
-                    objects,
-                    selectedFields
-                );
-                fvReconstructor.reconstructFvSurfaceFields<vector>
-                (
-                    objects,
-                    selectedFields
-                );
-                fvReconstructor.reconstructFvSurfaceFields<sphericalTensor>
-                (
-                    objects,
-                    selectedFields
-                );
-                fvReconstructor.reconstructFvSurfaceFields<symmTensor>
-                (
-                    objects,
-                    selectedFields
-                );
-                fvReconstructor.reconstructFvSurfaceFields<tensor>
-                (
-                    objects,
-                    selectedFields
-                );
-
-                if (fvReconstructor.nReconstructed() == 0)
-                {
-                    Info<< "No FV fields" << nl << endl;
-                }
+                writeDecomposition(regionMeshes.meshes(regioni)());
+                Info<< endl;
+                fileHandler().flush();
             }
+        }
 
-            if (!noFields)
+        // Do a region-by-region reconstruction of all the available fields
+        forAll(regionNames, regioni)
+        {
+            const word& regionName = regionNames[regioni];
+            const word regionDir =
+                regionName == polyMesh::defaultRegion ? word::null : regionName;
+
+            const delayedNewLine dnl;
+
+            // Prefixed scope
             {
-                Info<< "Reconstructing point fields" << nl << endl;
+                const RegionConstRef<domainDecomposition> meshes =
+                    regionMeshes.meshes(regioni);
 
-                const pointMesh& completePMesh =
-                    pointMesh::New(meshes.completeMesh());
-
-                pointFieldReconstructor pointReconstructor
+                // Search for objects at this time
+                IOobjectList objects
                 (
-                    completePMesh,
-                    meshes.procMeshes(),
-                    meshes.procPointAddressing()
+                    meshes().procMeshes()[0],
+                    runTimes.procTimes()[0].name()
                 );
 
-                pointReconstructor.reconstructFields<scalar>
-                (
-                    objects,
-                    selectedFields
-                );
-                pointReconstructor.reconstructFields<vector>
-                (
-                    objects,
-                    selectedFields
-                );
-                pointReconstructor.reconstructFields<sphericalTensor>
-                (
-                    objects,
-                    selectedFields
-                );
-                pointReconstructor.reconstructFields<symmTensor>
-                (
-                    objects,
-                    selectedFields
-                );
-                pointReconstructor.reconstructFields<tensor>
-                (
-                    objects,
-                    selectedFields
-                );
-
-                if (pointReconstructor.nReconstructed() == 0)
+                if (!noFields)
                 {
-                    Info<< "No point fields" << nl << endl;
-                }
-            }
+                    Info<< dnl << "Reconstructing FV fields" << endl;
 
-
-            // If there are any clouds, reconstruct them.
-            // The problem is that a cloud of size zero will not get written so
-            // in pass 1 we determine the cloud names and per cloud name the
-            // fields. Note that the fields are stored as IOobjectList from
-            // the first processor that has them. They are in pass2 only used
-            // for name and type (scalar, vector etc).
-
-            if (!noLagrangian)
-            {
-                HashTable<IOobjectList> cloudObjects;
-
-                forAll(runTimes.procTimes(), proci)
-                {
-                    fileName lagrangianDir
+                    if
                     (
-                        fileHandler().filePath
+                        fvFieldReconstructor::reconstructs
                         (
-                            runTimes.procTimes()[proci].timePath()
-                           /regionDir
-                           /cloud::prefix
+                            objects,
+                            selectedFields
                         )
-                    );
-
-                    fileNameList cloudDirs;
-                    if (!lagrangianDir.empty())
+                    )
                     {
-                        cloudDirs = fileHandler().readDir
+                        fvFieldReconstructor fvReconstructor
                         (
-                            lagrangianDir,
-                            fileType::directory
+                            meshes().completeMesh(),
+                            meshes().procMeshes(),
+                            meshes().procFaceAddressing(),
+                            meshes().procCellAddressing(),
+                            meshes().procFaceAddressingBf()
                         );
+
+                        #define DO_FV_VOL_INTERNAL_FIELDS_TYPE(Type, nullArg)  \
+                            fvReconstructor.reconstructVolInternalFields<Type> \
+                            (objects, selectedFields);
+                        FOR_ALL_FIELD_TYPES(DO_FV_VOL_INTERNAL_FIELDS_TYPE)
+                        #undef DO_FV_VOL_INTERNAL_FIELDS_TYPE
+
+                        #define DO_FV_VOL_FIELDS_TYPE(Type, nullArg)           \
+                            fvReconstructor.reconstructVolFields<Type>         \
+                            (objects, selectedFields);
+                        FOR_ALL_FIELD_TYPES(DO_FV_VOL_FIELDS_TYPE)
+                        #undef DO_FV_VOL_FIELDS_TYPE
+
+                        #define DO_FV_SURFACE_FIELDS_TYPE(Type, nullArg)       \
+                            fvReconstructor.reconstructFvSurfaceFields<Type>   \
+                            (objects, selectedFields);
+                        FOR_ALL_FIELD_TYPES(DO_FV_SURFACE_FIELDS_TYPE)
+                        #undef DO_FV_SURFACE_FIELDS_TYPE
                     }
-
-                    forAll(cloudDirs, i)
+                    else
                     {
-                        // Check if we already have cloud objects for this
-                        // cloudname
-                        HashTable<IOobjectList>::const_iterator iter =
-                            cloudObjects.find(cloudDirs[i]);
+                        Info<< dnl << "    (no FV fields)" << endl;
+                    }
+                }
 
-                        if (iter == cloudObjects.end())
-                        {
-                            // Do local scan for valid cloud objects
-                            IOobjectList sprayObjs
+                if (!noFields)
+                {
+                    Info<< dnl << "Reconstructing point fields" << endl;
+
+                    if
+                    (
+                        pointFieldReconstructor::reconstructs
+                        (
+                            objects,
+                            selectedFields
+                        )
+                    )
+                    {
+                        pointFieldReconstructor pointReconstructor
+                        (
+                            pointMesh::New(meshes().completeMesh()),
+                            meshes().procMeshes(),
+                            meshes().procPointAddressing()
+                        );
+
+                        #define DO_POINT_FIELDS_TYPE(Type, nullArg)            \
+                            pointReconstructor.reconstructFields<Type>         \
+                            (objects, selectedFields);
+                        FOR_ALL_FIELD_TYPES(DO_POINT_FIELDS_TYPE)
+                        #undef DO_POINT_FIELDS_TYPE
+                    }
+                    else
+                    {
+                        Info<< dnl << "    (no point fields)" << endl;
+                    }
+                }
+
+                if (!noLagrangian)
+                {
+                    // Search for clouds that exist on any processor and add
+                    // them into this table of cloud objects
+                    HashTable<IOobjectList> cloudsObjects;
+                    forAll(runTimes.procTimes(), proci)
+                    {
+                        // Find cloud directories
+                        fileNameList cloudDirs
+                        (
+                            fileHandler().readDir
                             (
-                                meshes.procMeshes()[proci],
+                                fileHandler().filePath
+                                (
+                                    runTimes.procTimes()[proci].timePath()
+                                   /regionDir
+                                   /cloud::prefix
+                                ),
+                                fileType::directory
+                            )
+                        );
+
+                        // Add objects in any found cloud directories
+                        forAll(cloudDirs, i)
+                        {
+                            // Pass if we already have an objects for this name
+                            HashTable<IOobjectList>::const_iterator iter =
+                                cloudsObjects.find(cloudDirs[i]);
+                            if (iter != cloudsObjects.end()) continue;
+
+                            // Do local scan for valid cloud objects
+                            IOobjectList cloudObjs
+                            (
+                                meshes().procMeshes()[proci],
                                 runTimes.procTimes()[proci].name(),
-                                cloud::prefix/cloudDirs[i]
+                                cloud::prefix/cloudDirs[i],
+                                IOobject::MUST_READ,
+                                IOobject::NO_WRITE,
+                                false
                             );
 
-                            IOobject* positionsPtr =
-                                sprayObjs.lookup(word("positions"));
-
-                            if (positionsPtr)
+                            // If "positions" is present, then add to the table
+                            if (cloudObjs.lookup(word("positions")))
                             {
-                                cloudObjects.insert(cloudDirs[i], sprayObjs);
+                                cloudsObjects.insert(cloudDirs[i], cloudObjs);
+                            }
+                        }
+                    }
+
+                    // Reconstruct the objects found above
+                    if (cloudsObjects.size())
+                    {
+                        forAllConstIter
+                        (
+                            HashTable<IOobjectList>,
+                            cloudsObjects,
+                            iter
+                        )
+                        {
+                            const word cloudName =
+                                string::validate<word>(iter.key());
+
+                            const IOobjectList& cloudObjects = iter();
+
+                            Info<< dnl << "Reconstructing lagrangian fields "
+                                << "for cloud " << cloudName << endl;
+
+                            if
+                            (
+                                lagrangianFieldReconstructor::reconstructs
+                                (
+                                    cloudObjects,
+                                    selectedLagrangianFields
+                                )
+                            )
+                            {
+                                lagrangianFieldReconstructor
+                                    lagrangianReconstructor
+                                    (
+                                        meshes().completeMesh(),
+                                        meshes().procMeshes(),
+                                        meshes().procFaceAddressing(),
+                                        meshes().procCellAddressing(),
+                                        cloudName
+                                    );
+
+                                #define DO_CLOUD_FIELDS_TYPE(Type, nullArg)    \
+                                    lagrangianReconstructor                    \
+                                   .reconstructFields<Type>                    \
+                                    (cloudObjects, selectedLagrangianFields);
+                                DO_CLOUD_FIELDS_TYPE(label, );
+                                FOR_ALL_FIELD_TYPES(DO_CLOUD_FIELDS_TYPE)
+                                #undef DO_CLOUD_FIELDS_TYPE
+                            }
+                            else
+                            {
+                                Info<< dnl << "    (no lagrangian fields)"
+                                    << endl;
                             }
                         }
                     }
                 }
-
-                if (cloudObjects.size())
-                {
-                    // Pass2: reconstruct the cloud
-                    forAllConstIter(HashTable<IOobjectList>, cloudObjects, iter)
-                    {
-                        const word cloudName =
-                            string::validate<word>(iter.key());
-
-                        // Objects (on arbitrary processor)
-                        const IOobjectList& sprayObjs = iter();
-
-                        Info<< "Reconstructing lagrangian fields for cloud "
-                            << cloudName << nl << endl;
-
-                        reconstructLagrangianPositions
-                        (
-                            meshes.completeMesh(),
-                            cloudName,
-                            meshes.procMeshes(),
-                            meshes.procFaceAddressing(),
-                            meshes.procCellAddressing()
-                        );
-                        reconstructLagrangianFields<label>
-                        (
-                            cloudName,
-                            meshes.completeMesh(),
-                            meshes.procMeshes(),
-                            sprayObjs,
-                            selectedLagrangianFields
-                        );
-                        reconstructLagrangianFieldFields<label>
-                        (
-                            cloudName,
-                            meshes.completeMesh(),
-                            meshes.procMeshes(),
-                            sprayObjs,
-                            selectedLagrangianFields
-                        );
-                        reconstructLagrangianFields<scalar>
-                        (
-                            cloudName,
-                            meshes.completeMesh(),
-                            meshes.procMeshes(),
-                            sprayObjs,
-                            selectedLagrangianFields
-                        );
-                        reconstructLagrangianFieldFields<scalar>
-                        (
-                            cloudName,
-                            meshes.completeMesh(),
-                            meshes.procMeshes(),
-                            sprayObjs,
-                            selectedLagrangianFields
-                        );
-                        reconstructLagrangianFields<vector>
-                        (
-                            cloudName,
-                            meshes.completeMesh(),
-                            meshes.procMeshes(),
-                            sprayObjs,
-                            selectedLagrangianFields
-                        );
-                        reconstructLagrangianFieldFields<vector>
-                        (
-                            cloudName,
-                            meshes.completeMesh(),
-                            meshes.procMeshes(),
-                            sprayObjs,
-                            selectedLagrangianFields
-                        );
-                        reconstructLagrangianFields<sphericalTensor>
-                        (
-                            cloudName,
-                            meshes.completeMesh(),
-                            meshes.procMeshes(),
-                            sprayObjs,
-                            selectedLagrangianFields
-                        );
-                        reconstructLagrangianFieldFields<sphericalTensor>
-                        (
-                            cloudName,
-                            meshes.completeMesh(),
-                            meshes.procMeshes(),
-                            sprayObjs,
-                            selectedLagrangianFields
-                        );
-                        reconstructLagrangianFields<symmTensor>
-                        (
-                            cloudName,
-                            meshes.completeMesh(),
-                            meshes.procMeshes(),
-                            sprayObjs,
-                            selectedLagrangianFields
-                        );
-                        reconstructLagrangianFieldFields<symmTensor>
-                        (
-                            cloudName,
-                            meshes.completeMesh(),
-                            meshes.procMeshes(),
-                            sprayObjs,
-                            selectedLagrangianFields
-                        );
-                        reconstructLagrangianFields<tensor>
-                        (
-                            cloudName,
-                            meshes.completeMesh(),
-                            meshes.procMeshes(),
-                            sprayObjs,
-                            selectedLagrangianFields
-                        );
-                        reconstructLagrangianFieldFields<tensor>
-                        (
-                            cloudName,
-                            meshes.completeMesh(),
-                            meshes.procMeshes(),
-                            sprayObjs,
-                            selectedLagrangianFields
-                        );
-                    }
-                }
-                else
-                {
-                    Info<< "No lagrangian fields" << nl << endl;
-                }
             }
 
-            // If there is a "uniform" directory in the time region
-            // directory copy from the master processor
+            Info<< dnl;
+        }
+
+        // Collect the uniform directory
+        if (haveUniform(runTimes))
+        {
+            Info<< "Collecting uniform files" << endl;
+
+            reconstructUniform(runTimes);
+
+            Info<< endl;
+        }
+
+        if (regionNames == wordList(1, polyMesh::defaultRegion)) continue;
+
+        // Collect the region uniform directories
+        forAll(regionNames, regioni)
+        {
+            const word& regionName = regionNames[regioni];
+            const word regionDir =
+                regionName == polyMesh::defaultRegion ? word::null : regionName;
+
+            if (haveUniform(runTimes, regionDir))
             {
-                fileName uniformDir0
-                (
-                    fileHandler().filePath
-                    (
-                        runTimes.procTimes()[0].timePath()/regionDir/"uniform"
-                    )
-                );
-
-                if (!uniformDir0.empty() && fileHandler().isDir(uniformDir0))
+                // Prefixed scope
                 {
-                    fileHandler().cp
-                    (
-                        uniformDir0,
-                        runTimes.completeTime().timePath()/regionDir
-                    );
-                }
-            }
+                    const RegionConstRef<domainDecomposition> meshes =
+                        regionMeshes.meshes(regioni);
 
-            // For the first region of a multi-region case additionally
-            // copy the "uniform" directory in the time directory
-            if (regioni == 0 && regionDir != word::null)
-            {
-                fileName uniformDir0
-                (
-                    fileHandler().filePath
-                    (
-                        runTimes.procTimes()[0].timePath()/"uniform"
-                    )
-                );
+                    Info<< "Collecting uniform files" << endl;
 
-                if (!uniformDir0.empty() && fileHandler().isDir(uniformDir0))
-                {
-                    fileHandler().cp
-                    (
-                        uniformDir0,
-                        runTimes.completeTime().timePath()
-                    );
+                    reconstructUniform(runTimes, regionDir);
                 }
+
+                Info<< endl;
             }
         }
     }
 
-    Info<< "\nEnd\n" << endl;
+    Info<< "End" << nl << endl;
 
     return 0;
 }
diff --git a/src/parallel/parallel/Make/files b/src/parallel/parallel/Make/files
index aa49671549..1b3416e703 100644
--- a/src/parallel/parallel/Make/files
+++ b/src/parallel/parallel/Make/files
@@ -2,5 +2,8 @@ processorRunTimes.C
 domainDecomposition.C
 domainDecompositionDecompose.C
 domainDecompositionReconstruct.C
+domainDecompositionNonConformal.C
+multiRegionPrefixer.C
+multiDomainDecomposition.C
 
 LIB = $(FOAM_LIBBIN)/libparallel
diff --git a/src/parallel/parallel/domainDecomposition.C b/src/parallel/parallel/domainDecomposition.C
index 3015fafe4f..16f6e30270 100644
--- a/src/parallel/parallel/domainDecomposition.C
+++ b/src/parallel/parallel/domainDecomposition.C
@@ -32,9 +32,7 @@ License
 #include "hexRef8Data.H"
 #include "cyclicFvPatch.H"
 #include "processorCyclicFvPatch.H"
-#include "nonConformalCyclicFvPatch.H"
-#include "nonConformalProcessorCyclicFvPatch.H"
-#include "nonConformalErrorFvPatch.H"
+#include "nonConformalFvPatch.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -46,471 +44,6 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::domainDecomposition::decomposePoints()
-{
-    for (label proci = 0; proci < nProcs(); proci++)
-    {
-        fvMesh& procMesh = procMeshes_[proci];
-
-        const label pointsCompare =
-            compareInstances
-            (
-                completeMesh().pointsInstance(),
-                procMeshes_[proci].pointsInstance()
-            );
-
-        if (pointsCompare == -1)
-        {
-            procMesh.setPoints
-            (
-                pointField
-                (
-                    completeMesh().points(),
-                    procPointAddressing_[proci]
-                )
-            );
-        }
-    }
-}
-
-
-void Foam::domainDecomposition::reconstructPoints()
-{
-    const label pointsCompare =
-        compareInstances
-        (
-            completeMesh().pointsInstance(),
-            procMeshes_[0].pointsInstance()
-        );
-
-    if (pointsCompare == 1)
-    {
-        pointField completePoints(completeMesh().nPoints());
-
-        for (label proci = 0; proci < nProcs(); proci++)
-        {
-            const fvMesh& procMesh = procMeshes_[proci];
-
-            completePoints.rmap(procMesh.points(), procPointAddressing_[proci]);
-        }
-
-        completeMesh_->setPoints(completePoints);
-    }
-}
-
-
-bool Foam::domainDecomposition::completeConformal() const
-{
-    return completeMesh().conformal();
-}
-
-
-bool Foam::domainDecomposition::procsConformal() const
-{
-    forAll(procMeshes_, proci)
-    {
-        if (!procMeshes_[proci].conformal())
-        {
-            return false;
-        }
-    }
-
-    return true;
-}
-
-
-void Foam::domainDecomposition::unconform()
-{
-    // ...
-    auto checkNonConformalCoupledPatchOrdering = []
-    (
-        const labelPair& procs,
-        const fvPatch& fvp,
-        const fvPatch& nbrFvp,
-        const labelList& polyFacesPf,
-        const labelList& nbrPolyFacesPf
-    )
-    {
-        if (fvp.size() != nbrFvp.size())
-        {
-            FatalErrorInFunction
-                << "Coupled patches " << fvp.name() << " and "
-                << nbrFvp.name() << " are not the same size"
-                << exit(FatalError);
-        }
-
-        if (fvp.size())
-        {
-            for (label i = 1; i < fvp.size(); ++ i)
-            {
-                if
-                (
-                    polyFacesPf[i - 1] > polyFacesPf[i]
-                  ? true
-                  : polyFacesPf[i - 1] == polyFacesPf[i]
-                  ? nbrPolyFacesPf[i - 1] >= nbrPolyFacesPf[i]
-                  : false
-                )
-                {
-                    FatalErrorInFunction
-                        << "Coupled patches " << fvp.name()
-                        << " and " << nbrFvp.name()
-                        << " are not in order";
-
-                    if (procs[0] == procs[1])
-                    {
-                        FatalErrorInFunction
-                            << " on process #" << procs[0];
-                    }
-
-                    FatalErrorInFunction
-                        << exit(FatalError);
-                }
-            }
-        }
-
-        InfoInFunction
-            << "Coupled patches " << fvp.name()
-            << " and " << nbrFvp.name()
-            << " are in order" << endl;
-    };
-
-    // ...
-    auto checkNonConformalErrorPatchOrdering = []
-    (
-        const label& proci,
-        const fvPatch& fvp,
-        const labelList& polyFacesPf
-    )
-    {
-        if (fvp.size())
-        {
-            for (label i = 1; i < fvp.size(); ++ i)
-            {
-                if (polyFacesPf[i - 1] > polyFacesPf[i])
-                {
-                    FatalErrorInFunction
-                        << "Error patch " << fvp.name()
-                        << " is not in order";
-
-                    if (proci > 0)
-                    {
-                        FatalErrorInFunction
-                            << " on process #" << proci;
-                    }
-
-                    FatalErrorInFunction
-                        << exit(FatalError);
-                }
-            }
-        }
-
-        InfoInFunction
-            << "Error patch " << fvp.name()
-            << " is in order" << endl;
-    };
-
-    if (completeConformal() && procsConformal())
-    {
-        // Nothing to do
-    }
-    else if (!completeConformal() && procsConformal())
-    {
-        // Construct the reverse of proc-face-face addressing
-        labelList faceProcFace(completeMesh().nFaces(), -labelMax);
-        forAll(procMeshes_, proci)
-        {
-            forAll(procFaceAddressing()[proci], procFacei)
-            {
-                const label facei =
-                    mag(procFaceAddressing()[proci][procFacei]) - 1;
-
-                faceProcFace[facei] =
-                    faceProcFace[facei] == -labelMax ? procFacei : -1;
-            }
-        }
-
-        // Map the polyFacesBf from the complete to the processor meshes and
-        // unconform the processor meshes. Set dummy data for the face
-        // geometry. This should not be used during decomposition.
-        forAll(procMeshes_, proci)
-        {
-            fvMesh& procMesh = procMeshes_[proci];
-
-            const surfaceLabelField::Boundary& faceAddressingBf =
-                procFaceAddressingBf()[proci];
-
-            surfaceLabelField::Boundary polyFacesBf
-            (
-                surfaceLabelField::null(),
-                procMesh.polyFacesBf()
-            );
-            surfaceVectorField Sf(procMesh.Sf().cloneUnSliced());
-            surfaceVectorField Cf(procMesh.Cf().cloneUnSliced());
-
-            forAll(procMesh.boundary(), procNccPatchi)
-            {
-                const fvPatch& fvp = procMesh.boundary()[procNccPatchi];
-
-                if (isA<nonConformalFvPatch>(fvp))
-                {
-                    const label completeNccPatchi =
-                        isA<processorCyclicFvPatch>(fvp)
-                      ? refCast<const processorCyclicFvPatch>(fvp)
-                       .referPatchID()
-                      : procNccPatchi;
-
-                    polyFacesBf[procNccPatchi] =
-                        labelField
-                        (
-                            faceProcFace,
-                            labelField
-                            (
-                                completeMesh().polyFacesBf()[completeNccPatchi],
-                                mag(faceAddressingBf[procNccPatchi]) - 1
-                            )
-                        );
-
-                    const label size = polyFacesBf[procNccPatchi].size();
-
-                    Sf.boundaryFieldRef()[procNccPatchi].resize(size, Zero);
-                    Cf.boundaryFieldRef()[procNccPatchi].resize(size, Zero);
-                }
-            }
-
-            procMesh.unconform
-            (
-                polyFacesBf,
-                Sf,
-                Cf,
-                NullObjectRef<surfaceScalarField>(),
-                false
-            );
-        }
-
-        // Check ordering
-        if (debug)
-        {
-            forAll(procMeshes_, proci)
-            {
-                forAll(procMeshes_[proci].boundary(), patchi)
-                {
-                    const fvPatch& fvp =
-                        procMeshes_[proci].boundary()[patchi];
-
-                    // Coupled patches
-                    if
-                    (
-                        isA<nonConformalCoupledFvPatch>(fvp)
-                     && refCast<const nonConformalCoupledFvPatch>(fvp).owner()
-                    )
-                    {
-                        const label nccPatchi = patchi;
-
-                        label nbrProci = -1, nbrNccPatchi = -1;
-                        if (isA<cyclicFvPatch>(fvp))
-                        {
-                            nbrProci = proci;
-                            nbrNccPatchi =
-                                refCast<const cyclicFvPatch>(fvp).nbrPatchID();
-                        }
-                        else if (isA<processorCyclicFvPatch>(fvp))
-                        {
-                            typedef processorCyclicFvPatch PcFvp;
-
-                            const PcFvp& pcFvp = refCast<const PcFvp>(fvp);
-
-                            nbrProci = pcFvp.neighbProcNo();
-
-                            const fvBoundaryMesh& nbrFvPatches =
-                                procMeshes_[nbrProci].boundary();
-
-                            forAll(nbrFvPatches, nbrNccPatchj)
-                            {
-                                const fvPatch& nbrFvp =
-                                    nbrFvPatches[nbrNccPatchj];
-
-                                if (isA<PcFvp>(nbrFvp))
-                                {
-                                    const PcFvp& nbrPcFvp =
-                                        refCast<const PcFvp>(nbrFvp);
-
-                                    if
-                                    (
-                                        nbrPcFvp.neighbProcNo()
-                                     == proci
-                                     && nbrPcFvp.referPatchID()
-                                     == pcFvp.referPatch().nbrPatchID()
-                                    )
-                                    {
-                                        nbrNccPatchi = nbrNccPatchj;
-                                        break;
-                                    }
-                                }
-                            }
-
-                            if (nbrNccPatchi == -1)
-                            {
-                                FatalErrorInFunction
-                                    << "Opposite processor patch not found for "
-                                    << "patch " << fvp.name() << " on proc #"
-                                    << proci << exit(FatalError);
-                            }
-                        }
-                        else
-                        {
-                            FatalErrorInFunction
-                                << "Non-conformal-coupled type not recognised "
-                                << "for patch " << fvp.name() << " on proc #"
-                                << proci << exit(FatalError);
-                        }
-
-                        checkNonConformalCoupledPatchOrdering
-                        (
-                            {proci, nbrProci},
-                            procMeshes_[proci].boundary()[nccPatchi],
-                            procMeshes_[nbrProci].boundary()[nbrNccPatchi],
-                            procMeshes_[proci].polyFacesBf()[nccPatchi],
-                            procMeshes_[nbrProci].polyFacesBf()[nbrNccPatchi]
-                        );
-                    }
-
-                    // Error patches
-                    if (isA<nonConformalErrorFvPatch>(fvp))
-                    {
-                        const label ncePatchi = patchi;
-
-                        checkNonConformalErrorPatchOrdering
-                        (
-                            proci,
-                            procMeshes_[proci].boundary()[ncePatchi],
-                            procMeshes_[proci].polyFacesBf()[ncePatchi]
-                        );
-                    }
-                }
-            }
-        }
-    }
-    else if (completeConformal() && !procsConformal())
-    {
-        // Map the polyFacesBf from the processor to the complete meshes and
-        // unconform the complete mesh. Set dummy data for the face
-        // geometry. This should not be used during decomposition.
-        surfaceLabelField::Boundary polyFacesBf
-        (
-            surfaceLabelField::null(),
-            completeMesh().polyFacesBf()
-        );
-        surfaceVectorField Sf(completeMesh().Sf().cloneUnSliced());
-        surfaceVectorField Cf(completeMesh().Cf().cloneUnSliced());
-
-        forAll(procMeshes_, proci)
-        {
-            const fvMesh& procMesh = procMeshes_[proci];
-
-            const surfaceLabelField::Boundary& faceAddressingBf =
-                procFaceAddressingBf()[proci];
-
-            forAll(procMesh.boundary(), procNccPatchi)
-            {
-                const fvPatch& fvp = procMesh.boundary()[procNccPatchi];
-
-                if (isA<nonConformalFvPatch>(fvp))
-                {
-                    const label completeNccPatchi =
-                        isA<processorCyclicFvPatch>(fvp)
-                      ? refCast<const processorCyclicFvPatch>(fvp)
-                       .referPatchID()
-                      : procNccPatchi;
-
-                    const label size =
-                        max
-                        (
-                            max(mag(faceAddressingBf[procNccPatchi])),
-                            polyFacesBf[completeNccPatchi].size()
-                        );
-
-                    polyFacesBf[completeNccPatchi].resize(size, -1);
-                    Sf.boundaryFieldRef()[completeNccPatchi].resize(size, Zero);
-                    Cf.boundaryFieldRef()[completeNccPatchi].resize(size, Zero);
-
-                    polyFacesBf[completeNccPatchi].labelField::rmap
-                    (
-                        mag
-                        (
-                            labelField
-                            (
-                                procFaceAddressing_[proci],
-                                procMesh.polyFacesBf()[procNccPatchi]
-                            )
-                        ) - 1,
-                        mag(faceAddressingBf[procNccPatchi]) - 1
-                    );
-                }
-            }
-        }
-
-        completeMesh_->unconform
-        (
-            polyFacesBf,
-            Sf,
-            Cf,
-            NullObjectRef<surfaceScalarField>(),
-            false
-        );
-
-        // Check ordering
-        if (debug)
-        {
-            forAll(completeMesh().boundary(), patchi)
-            {
-                const fvPatch& fvp = completeMesh().boundary()[patchi];
-
-                // Coupled patches
-                if
-                (
-                    isA<nonConformalCyclicFvPatch>(fvp)
-                 && refCast<const nonConformalCoupledFvPatch>(fvp).owner()
-                )
-                {
-                    const label nccPatchi = patchi;
-
-                    const label nbrNccPatchi =
-                        refCast<const nonConformalCyclicFvPatch>(fvp)
-                       .nbrPatchID();
-
-                    checkNonConformalCoupledPatchOrdering
-                    (
-                        {-labelMax, labelMax},
-                        completeMesh().boundary()[nccPatchi],
-                        completeMesh().boundary()[nbrNccPatchi],
-                        completeMesh().polyFacesBf()[nccPatchi],
-                        completeMesh().polyFacesBf()[nbrNccPatchi]
-                    );
-                }
-
-                // Error patches
-                if (isA<nonConformalErrorFvPatch>(fvp))
-                {
-                    const label ncePatchi = patchi;
-
-                    checkNonConformalErrorPatchOrdering
-                    (
-                        -labelMax,
-                        completeMesh().boundary()[ncePatchi],
-                        completeMesh().polyFacesBf()[ncePatchi]
-                    );
-                }
-            }
-        }
-    }
-    else // if (!completeConformal() && !procsConformal())
-    {
-        // Assume everything is consistent and do nothing
-    }
-}
-
-
 Foam::label Foam::domainDecomposition::compareInstances
 (
     const fileName& a,
@@ -1284,277 +817,9 @@ Foam::domainDecomposition::procFaceAddressingBf() const
 
     if (procFaceAddressingBf_.empty())
     {
-        // Map from reference patch and processors to the interface patch
-        typedef HashTable<label, labelPair, Hash<labelPair>> labelPairTable;
-        List<labelPairTable> refPatchProcPatchTable
-        (
-            completeMesh().boundary().size()
-        );
-        forAll(procMeshes_, proci)
-        {
-            const fvMesh& procMesh = procMeshes_[proci];
-
-            forAll(procMesh.boundary(), procPatchi)
-            {
-                const fvPatch& fvp = procMesh.boundary()[procPatchi];
-
-                if (isA<cyclicFvPatch>(fvp))
-                {
-                    refPatchProcPatchTable[procPatchi].insert
-                    (
-                        labelPair(proci, proci),
-                        procPatchi
-                    );
-                }
-                else if (isA<processorCyclicFvPatch>(fvp))
-                {
-                    const processorCyclicFvPatch& pcFvp =
-                        refCast<const processorCyclicFvPatch>(fvp);
-
-                    refPatchProcPatchTable[pcFvp.referPatchID()].insert
-                    (
-                        labelPair(proci, pcFvp.neighbProcNo()),
-                        procPatchi
-                    );
-                }
-            }
-        }
-
-        // Build non-conformal finite volume face addressing for each processor
-        List<List<DynamicList<label>>>
-            nonConformalProcFaceAddressingBf(nProcs());
-        forAll(nonConformalProcFaceAddressingBf, proci)
-        {
-            nonConformalProcFaceAddressingBf[proci].resize
-            (
-                procMeshes_[proci].boundary().size()
-            );
-        }
-        if (completeConformal() && procsConformal())
-        {
-            // Nothing to do
-        }
-        else if (!completeConformal())
-        {
-            // Decompose non-conformal addressing
-
-            const surfaceLabelField::Boundary& polyFacesBf =
-                completeMesh().polyFacesBf();
-
-            // Cyclic patches
-            forAll(completeMesh().boundary(), nccPatchi)
-            {
-                const fvPatch& fvp = completeMesh().boundary()[nccPatchi];
-
-                if (!isA<nonConformalCyclicFvPatch>(fvp)) continue;
-
-                const nonConformalCyclicFvPatch& nccFvp =
-                    refCast<const nonConformalCyclicFvPatch>(fvp);
-
-                if (!nccFvp.owner()) continue;
-
-                const label nccNbrPatchi = nccFvp.nbrPatchID();
-
-                forAll(polyFacesBf[nccPatchi], nccPatchFacei)
-                {
-                    const label facei = polyFacesBf[nccPatchi][nccPatchFacei];
-                    const label celli = completeMesh().faceOwner()[facei];
-                    const label proci = cellProc_[celli];
-
-                    const label nbrFacei =
-                        polyFacesBf[nccNbrPatchi][nccPatchFacei];
-                    const label nbrCelli =
-                        completeMesh().faceOwner()[nbrFacei];
-                    const label nbrProci = cellProc_[nbrCelli];
-
-                    const label procNccPatchi =
-                        refPatchProcPatchTable
-                        [nccPatchi][labelPair(proci, nbrProci)];
-                    const label nbrProcNccPatchi =
-                        refPatchProcPatchTable
-                        [nccNbrPatchi][labelPair(nbrProci, proci)];
-
-                    nonConformalProcFaceAddressingBf[proci][procNccPatchi]
-                       .append(nccPatchFacei + 1);
-                    nonConformalProcFaceAddressingBf[nbrProci][nbrProcNccPatchi]
-                       .append(nccPatchFacei + 1);
-                }
-            }
-
-            // Error patches
-            forAll(completeMesh().boundary(), ncePatchi)
-            {
-                const fvPatch& fvp = completeMesh().boundary()[ncePatchi];
-
-                if (!isA<nonConformalErrorFvPatch>(fvp)) continue;
-
-                forAll(polyFacesBf[ncePatchi], ncePatchFacei)
-                {
-                    const label facei = polyFacesBf[ncePatchi][ncePatchFacei];
-                    const label celli = completeMesh().faceOwner()[facei];
-                    const label proci = cellProc_[celli];
-
-                    nonConformalProcFaceAddressingBf[proci][ncePatchi]
-                       .append(ncePatchFacei + 1);
-                }
-            }
-        }
-        else // if (!procsConformal())
-        {
-            // Reconstruct non-conformal addressing
-
-            // Cyclic patches
-            forAll(completeMesh().boundary(), nccPatchi)
-            {
-                const fvPatch& fvp = completeMesh().boundary()[nccPatchi];
-
-                if (!isA<nonConformalCyclicFvPatch>(fvp)) continue;
-
-                const nonConformalCyclicFvPatch& nccFvp =
-                    refCast<const nonConformalCyclicFvPatch>(fvp);
-
-                if (!nccFvp.owner()) continue;
-
-                const label nccNbrPatchi = nccFvp.nbrPatchID();
-
-                label count = 0;
-                labelPairTable procCounts;
-                forAllConstIter
-                (
-                    labelPairTable,
-                    refPatchProcPatchTable[nccPatchi],
-                    iter
-                )
-                {
-                    procCounts.insert(iter.key(), 0);
-                }
-
-                while (true)
-                {
-                    labelPair procNbrProc(labelMax, labelMax);
-                    labelPair faceNbrFace(labelMax, labelMax);
-
-                    forAllConstIter(labelPairTable, procCounts, iter)
-                    {
-                        const label proci = iter.key().first();
-                        const label nbrProci = iter.key().second();
-
-                        const labelPair procNbrProcStar(proci, nbrProci);
-                        const labelPair nbrProcProcStar(nbrProci, proci);
-
-                        const label procNccPatchi =
-                            refPatchProcPatchTable
-                            [nccPatchi][procNbrProcStar];
-                        const label nbrProcNccPatchi =
-                            refPatchProcPatchTable
-                            [nccNbrPatchi][nbrProcProcStar];
-
-                        const label size =
-                            procMeshes_[proci]
-                           .polyFacesBf()[procNccPatchi]
-                           .size();
-
-                        if (iter() >= size) continue;
-
-                        const label procFacei =
-                            procMeshes_[proci].polyFacesBf()
-                            [procNccPatchi][iter()];
-                        const label nbrProcFacei =
-                            procMeshes_[nbrProci].polyFacesBf()
-                            [nbrProcNccPatchi][iter()];
-
-                        const labelPair faceNbrFaceStar
-                        (
-                            procFaceAddressing_[proci][procFacei] - 1,
-                            procFaceAddressing_[nbrProci][nbrProcFacei] - 1
-                        );
-
-                        if (faceNbrFace > faceNbrFaceStar)
-                        {
-                            procNbrProc = procNbrProcStar;
-                            faceNbrFace = faceNbrFaceStar;
-                        }
-                    }
-
-                    if (faceNbrFace == labelPair(labelMax, labelMax))
-                    {
-                        break;
-                    }
-                    else
-                    {
-                        const label proci = procNbrProc.first();
-                        const label nbrProci = procNbrProc.second();
-
-                        const labelPair nbrProcProc(nbrProci, proci);
-
-                        const label procNccPatchi =
-                            refPatchProcPatchTable[nccPatchi][procNbrProc];
-                        const label nbrProcNccPatchi =
-                            refPatchProcPatchTable[nccNbrPatchi][nbrProcProc];
-
-                        nonConformalProcFaceAddressingBf
-                            [proci][procNccPatchi].append(count + 1);
-                        nonConformalProcFaceAddressingBf
-                            [nbrProci][nbrProcNccPatchi].append(count + 1);
-
-                        count ++;
-                        procCounts[procNbrProc] ++;
-                    }
-                }
-            }
-
-            // Error patches
-            forAll(completeMesh().boundary(), ncePatchi)
-            {
-                const fvPatch& fvp = completeMesh().boundary()[ncePatchi];
-
-                if (!isA<nonConformalErrorFvPatch>(fvp)) continue;
-
-                label count = 0;
-                labelList procCounts(nProcs(), 0);
-
-                while (true)
-                {
-                    label facei = labelMax, proci = labelMax;
-
-                    forAll(procCounts, procStari)
-                    {
-                        const label size =
-                            procMeshes_[procStari]
-                           .polyFacesBf()[ncePatchi]
-                           .size();
-
-                        if (procCounts[procStari] >= size) continue;
-
-                        const label procFacei =
-                            procMeshes_[procStari].polyFacesBf()
-                            [ncePatchi][procCounts[procStari]];
-
-                        const label faceStari =
-                            procFaceAddressing_[procStari][procFacei] - 1;
-
-                        if (facei > faceStari)
-                        {
-                            facei = faceStari;
-                            proci = procStari;
-                        }
-                    }
-
-                    if (facei == labelMax)
-                    {
-                        break;
-                    }
-                    else
-                    {
-                        nonConformalProcFaceAddressingBf
-                            [proci][ncePatchi].append(count + 1);
-
-                        count ++;
-                        procCounts[proci] ++;
-                    }
-                }
-            }
-        }
+        // Get any non-conformal proc-face addressing
+        List<List<DynamicList<label>>> nonConformalProcFaceAddressingBf =
+            this->nonConformalProcFaceAddressingBf();
 
         // Build finite volume face addressing boundary fields
         procFaceAddressingBf_.resize(nProcs());
diff --git a/src/parallel/parallel/domainDecomposition.H b/src/parallel/parallel/domainDecomposition.H
index 1ff1192fa8..2b53c3cbfb 100644
--- a/src/parallel/parallel/domainDecomposition.H
+++ b/src/parallel/parallel/domainDecomposition.H
@@ -31,6 +31,7 @@ SourceFiles
     domainDecomposition.C
     domainDecompositionDecompose.C
     domainDecompositionReconstruct.C
+    domainDecompositionReconstruct.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -47,6 +48,8 @@ SourceFiles
 namespace Foam
 {
 
+class faceCoupleInfo;
+
 /*---------------------------------------------------------------------------*\
                      Class domainDecomposition Declaration
 \*---------------------------------------------------------------------------*/
@@ -100,8 +103,6 @@ class domainDecomposition
             mutable PtrList<surfaceLabelField::Boundary> procFaceAddressingBf_;
 
 
-
-
     // Private Member Functions
 
         // Decomposition
@@ -145,25 +146,62 @@ class domainDecomposition
                 BinaryOp bop
             ) const;
 
+            //- Decompose the complete mesh to create the processor meshes and
+            //  populate the addressing
+            void decompose();
 
-        //- Decompose the complete mesh's points and apply the result to the
-        //  processor meshes
-        void decomposePoints();
+            //- Decompose the complete mesh's points and apply the result to
+            //  the processor meshes
+            void decomposePoints();
 
-        //- Reconstruct the processor meshes' points and apply the result to
-        //  the complete mesh
-        void reconstructPoints();
 
-        //- Is the complete mesh conformal?
-        bool completeConformal() const;
+        // Reconstruction
 
-        //- Are the processor meshes conformal?
-        bool procsConformal() const;
+            //- Find shared points and faces between two meshes that are to be
+            //  added together
+            static autoPtr<faceCoupleInfo> determineCoupledFaces
+            (
+                const label masterMeshProcStart,
+                const label masterMeshProcEnd,
+                const polyMesh& masterMesh,
+                const label meshToAddProcStart,
+                const label meshToAddProcEnd,
+                const polyMesh& meshToAdd
+            );
+
+            //- Reconstruct the processor meshes to create the complete mesh
+            //  and populate the addressing
+            void reconstruct();
+
+            //- Reconstruct the processor meshes' points and apply the result
+            //  to the complete mesh
+            void reconstructPoints();
+
+
+        // Non conformal
+
+            //- Is the complete mesh conformal?
+            bool completeConformal() const;
+
+            //- Are the processor meshes conformal?
+            bool procsConformal() const;
+
+            //- ...
+            List<List<DynamicList<label>>>
+                nonConformalProcFaceAddressingBf() const;
+
+            // Map the polyFacesBf from the processor to the complete meshes and
+            // unconform the complete mesh
+            void unconformComplete();
+
+            // Map the polyFacesBf from the complete to the processor meshes and
+            // unconform the processor meshes
+            void unconformProcs();
+
+            //- If the complete mesh is non-conformal and the processor meshes
+            //  are not, then unconform the processor meshes. And vice versa.
+            void unconform();
 
-        //- If the complete mesh is non-conformal and the processor meshes are
-        //  not, then map the non-conformal addressing across and unconform the
-        //  processor meshes. And vice versa.
-        void unconform();
 
         //- Compare two instances. A return value of -1 means a is newer than b
         //  (i.e., the arguments are in reverse order), 0 means a is equal to
@@ -194,14 +232,6 @@ class domainDecomposition
         //- Read-update the complete and processor meshes for a change in time
         fvMesh::readUpdateState readUpdate();
 
-        //- Decompose the complete mesh to create the processor meshes and
-        //  populate the addressing
-        void decompose();
-
-        //- Reconstruct the processor meshes to create the complete mesh and
-        //  populate the addressing
-        void reconstruct();
-
         //- Write the decomposition addressing
         void writeCompleteAddressing() const;
 
@@ -244,22 +274,34 @@ public:
 
     // Member Functions
 
+        //- Access the run times
+        inline const processorRunTimes& procRunTimes() const
+        {
+            return runTimes_;
+        }
+
+        //- Access the region name
+        inline const word& regionName() const
+        {
+            return regionName_;
+        }
+
         //- Access the global mesh
-        const fvMesh& completeMesh() const
+        inline const fvMesh& completeMesh() const
         {
             validateComplete();
             return completeMesh_();
         }
 
         //- Access the processor meshes
-        const PtrList<fvMesh>& procMeshes() const
+        inline const PtrList<fvMesh>& procMeshes() const
         {
             validateProcs();
             return procMeshes_;
         }
 
         //- Return the number of processors in the decomposition
-        label nProcs() const
+        inline label nProcs() const
         {
             return runTimes_.nProcs();
         }
@@ -281,27 +323,27 @@ public:
         fvMesh::readUpdateState readUpdateReconstruct();
 
         //- Return the distribution as an FV field for writing
-        const labelList& cellProc() const
+        inline const labelList& cellProc() const
         {
             return cellProc_;
         }
 
         //- Access the labels of points for each processor
-        const labelListList& procPointAddressing() const
+        inline const labelListList& procPointAddressing() const
         {
             validateProcs();
             return procPointAddressing_;
         }
 
         //- Access the labels of faces for each processor (see notes above)
-        const labelListList& procFaceAddressing() const
+        inline const labelListList& procFaceAddressing() const
         {
             validateProcs();
             return procFaceAddressing_;
         }
 
         //- Access the labels of cells for each processor
-        const labelListList& procCellAddressing() const
+        inline const labelListList& procCellAddressing() const
         {
             validateProcs();
             return procCellAddressing_;
diff --git a/src/parallel/parallel/domainDecompositionDecompose.C b/src/parallel/parallel/domainDecompositionDecompose.C
index 3d4e1eb5a1..45cddbb7e4 100644
--- a/src/parallel/parallel/domainDecompositionDecompose.C
+++ b/src/parallel/parallel/domainDecompositionDecompose.C
@@ -133,7 +133,7 @@ void Foam::domainDecomposition::addInterProcFace
 
 Foam::labelList Foam::domainDecomposition::distributeCells()
 {
-    Info<< "\nCalculating distribution of cells" << endl;
+    Info<< "Calculating distribution of cells" << nl << endl;
 
     cpuTime decompositionTime;
 
@@ -167,7 +167,7 @@ Foam::labelList Foam::domainDecomposition::distributeCells()
             cellWeights
         );
 
-    Info<< "\nFinished decomposition in "
+    Info<< nl << "Finished decomposition in "
         << decompositionTime.elapsedCpuTime()
         << " s" << endl;
 
@@ -269,17 +269,16 @@ inline void Foam::domainDecomposition::processInterCyclics
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
 void Foam::domainDecomposition::decompose()
 {
     // Decide which cell goes to which processor
     cellProc_ = distributeCells();
+    Info<< nl;
 
     // Distribute the cells according to the given processor label
 
     // calculate the addressing information for the original mesh
-    Info<< "\nCalculating original mesh data" << endl;
+    Info<< "Calculating original mesh data" << nl << endl;
 
     // set references to the original mesh
     const polyBoundaryMesh& patches = completeMesh().boundaryMesh();
@@ -290,12 +289,12 @@ void Foam::domainDecomposition::decompose()
     // loop through the list of processor labels for the cell and add the
     // cell shape to the list of cells for the appropriate processor
 
-    Info<< "\nDistributing cells to processors" << endl;
+    Info<< "Distributing cells to processors" << nl << endl;
 
     // Cells per processor
     procCellAddressing_ = invertOneToMany(nProcs(), cellProc_);
 
-    Info<< "\nDistributing faces to processors" << endl;
+    Info<< "Distributing faces to processors" << nl << endl;
 
     // Loop through all internal faces and decide which processor they belong to
     // First visit all internal faces. If cells at both sides belong to the
@@ -606,7 +605,7 @@ void Foam::domainDecomposition::decompose()
         }
     }
 
-    Info<< "\nDistributing points to processors" << endl;
+    Info<< "Distributing points to processors" << nl << endl;
 
     // For every processor, loop through the list of faces for the processor.
     // For every face, loop through the list of points and mark the point as
@@ -653,7 +652,7 @@ void Foam::domainDecomposition::decompose()
         procPointLabels.setSize(nUsedPoints);
     }
 
-    Info<< "\nConstructing processor meshes" << endl;
+    Info<< "Constructing processor meshes" << nl << endl;
 
     // Mark point/faces/cells that are in zones.
     // -1   : not in zone
@@ -1150,8 +1149,7 @@ void Foam::domainDecomposition::decompose()
 
         // Report processor and update global statistics
         {
-            Info<< endl
-                << "Processor " << proci << nl
+            Info<< "Processor " << proci << nl
                 << "    Number of cells = " << procMesh.nCells()
                 << endl;
 
@@ -1185,7 +1183,8 @@ void Foam::domainDecomposition::decompose()
 
             Info<< "    Number of processor patches = " << nProcPatches << nl
                 << "    Number of processor faces = " << nProcFaces << nl
-                << "    Number of boundary faces = " << nBoundaryFaces << endl;
+                << "    Number of boundary faces = " << nBoundaryFaces << nl
+                << endl;
 
             totProcFaces += nProcFaces;
             totProcPatches += nProcPatches;
@@ -1209,8 +1208,7 @@ void Foam::domainDecomposition::decompose()
         avgProcFaces = 1;
     }
 
-    Info<< nl
-        << "Number of processor faces = " << totProcFaces/2 << nl
+    Info<< "Number of processor faces = " << totProcFaces/2 << nl
         << "Max number of cells = " << maxProcCells
         << " (" << 100.0*(maxProcCells-avgProcCells)/avgProcCells
         << "% above average " << avgProcCells << ")" << nl
@@ -1219,8 +1217,7 @@ void Foam::domainDecomposition::decompose()
         << "% above average " << avgProcPatches << ")" << nl
         << "Max number of faces between processors = " << maxProcFaces
         << " (" << 100.0*(maxProcFaces-avgProcFaces)/avgProcFaces
-        << "% above average " << avgProcFaces << ")" << nl
-        << endl;
+        << "% above average " << avgProcFaces << ")" << endl;
 
     // Clear (and thus trigger re-generation) of finite volume face addressing
     procFaceAddressingBf_.clear();
@@ -1230,4 +1227,32 @@ void Foam::domainDecomposition::decompose()
 }
 
 
+void Foam::domainDecomposition::decomposePoints()
+{
+    for (label proci = 0; proci < nProcs(); proci++)
+    {
+        fvMesh& procMesh = procMeshes_[proci];
+
+        const label pointsCompare =
+            compareInstances
+            (
+                completeMesh().pointsInstance(),
+                procMeshes_[proci].pointsInstance()
+            );
+
+        if (pointsCompare == -1)
+        {
+            procMesh.setPoints
+            (
+                pointField
+                (
+                    completeMesh().points(),
+                    procPointAddressing_[proci]
+                )
+            );
+        }
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/src/parallel/parallel/domainDecompositionNonConformal.C b/src/parallel/parallel/domainDecompositionNonConformal.C
new file mode 100644
index 0000000000..37fa479865
--- /dev/null
+++ b/src/parallel/parallel/domainDecompositionNonConformal.C
@@ -0,0 +1,744 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "domainDecomposition.H"
+#include "cyclicFvPatch.H"
+#include "processorCyclicFvPatch.H"
+#include "nonConformalCyclicFvPatch.H"
+#include "nonConformalProcessorCyclicFvPatch.H"
+#include "nonConformalErrorFvPatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+void checkNonConformalCoupledPatchOrdering
+(
+    const labelPair& procs,
+    const fvPatch& fvp,
+    const fvPatch& nbrFvp,
+    const labelList& polyFacesPf,
+    const labelList& nbrPolyFacesPf
+)
+{
+    if (fvp.size() != nbrFvp.size())
+    {
+        FatalErrorInFunction
+            << "Coupled patches " << fvp.name() << " and "
+            << nbrFvp.name() << " are not the same size"
+            << exit(FatalError);
+    }
+
+    if (fvp.size())
+    {
+        for (label i = 1; i < fvp.size(); ++ i)
+        {
+            if
+            (
+                polyFacesPf[i - 1] > polyFacesPf[i]
+              ? true
+              : polyFacesPf[i - 1] == polyFacesPf[i]
+              ? nbrPolyFacesPf[i - 1] >= nbrPolyFacesPf[i]
+              : false
+            )
+            {
+                FatalErrorInFunction
+                    << "Coupled patches " << fvp.name()
+                    << " and " << nbrFvp.name()
+                    << " are not in order";
+
+                if (procs[0] == procs[1])
+                {
+                    FatalErrorInFunction
+                        << " on process #" << procs[0];
+                }
+
+                FatalErrorInFunction
+                    << exit(FatalError);
+            }
+        }
+    }
+
+    InfoInFunction
+        << "Coupled patches " << fvp.name()
+        << " and " << nbrFvp.name()
+        << " are in order" << endl;
+}
+
+
+void checkNonConformalErrorPatchOrdering
+(
+    const label& proci,
+    const fvPatch& fvp,
+    const labelList& polyFacesPf
+)
+{
+    if (fvp.size())
+    {
+        for (label i = 1; i < fvp.size(); ++ i)
+        {
+            if (polyFacesPf[i - 1] > polyFacesPf[i])
+            {
+                FatalErrorInFunction
+                    << "Error patch " << fvp.name()
+                    << " is not in order";
+
+                if (proci > 0)
+                {
+                    FatalErrorInFunction
+                        << " on process #" << proci;
+                }
+
+                FatalErrorInFunction
+                    << exit(FatalError);
+            }
+        }
+    }
+
+    InfoInFunction
+        << "Error patch " << fvp.name()
+        << " is in order" << endl;
+}
+
+
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+bool Foam::domainDecomposition::completeConformal() const
+{
+    return completeMesh().conformal();
+}
+
+
+bool Foam::domainDecomposition::procsConformal() const
+{
+    forAll(procMeshes_, proci)
+    {
+        if (!procMeshes_[proci].conformal())
+        {
+            return false;
+        }
+    }
+
+    return true;
+}
+
+
+Foam::List<Foam::List<Foam::DynamicList<Foam::label>>>
+Foam::domainDecomposition::nonConformalProcFaceAddressingBf() const
+{
+    validateComplete();
+    validateProcs();
+
+    // Map from reference patch and processors to the interface patch
+    typedef HashTable<label, labelPair, Hash<labelPair>> labelPairTable;
+    List<labelPairTable> refPatchProcPatchTable
+    (
+        completeMesh().boundary().size()
+    );
+    forAll(procMeshes_, proci)
+    {
+        const fvMesh& procMesh = procMeshes_[proci];
+
+        forAll(procMesh.boundary(), procPatchi)
+        {
+            const fvPatch& fvp = procMesh.boundary()[procPatchi];
+
+            if (isA<cyclicFvPatch>(fvp))
+            {
+                refPatchProcPatchTable[procPatchi].insert
+                (
+                    labelPair(proci, proci),
+                    procPatchi
+                );
+            }
+            else if (isA<processorCyclicFvPatch>(fvp))
+            {
+                const processorCyclicFvPatch& pcFvp =
+                    refCast<const processorCyclicFvPatch>(fvp);
+
+                refPatchProcPatchTable[pcFvp.referPatchID()].insert
+                (
+                    labelPair(proci, pcFvp.neighbProcNo()),
+                    procPatchi
+                );
+            }
+        }
+    }
+
+    // Build non-conformal finite volume face addressing for each processor
+    List<List<DynamicList<label>>> result(nProcs());
+    forAll(result, proci)
+    {
+        result[proci].resize
+        (
+            procMeshes_[proci].boundary().size()
+        );
+    }
+
+    if (completeConformal() && procsConformal())
+    {
+        // Nothing to do
+    }
+    else if (!completeConformal())
+    {
+        // Decompose non-conformal addressing
+
+        const surfaceLabelField::Boundary& polyFacesBf =
+            completeMesh().polyFacesBf();
+
+        // Cyclic patches
+        forAll(completeMesh().boundary(), nccPatchi)
+        {
+            const fvPatch& fvp = completeMesh().boundary()[nccPatchi];
+
+            if (!isA<nonConformalCyclicFvPatch>(fvp)) continue;
+
+            const nonConformalCyclicFvPatch& nccFvp =
+                refCast<const nonConformalCyclicFvPatch>(fvp);
+
+            if (!nccFvp.owner()) continue;
+
+            const label nccNbrPatchi = nccFvp.nbrPatchID();
+
+            forAll(polyFacesBf[nccPatchi], nccPatchFacei)
+            {
+                const label facei = polyFacesBf[nccPatchi][nccPatchFacei];
+                const label celli = completeMesh().faceOwner()[facei];
+                const label proci = cellProc_[celli];
+
+                const label nbrFacei =
+                    polyFacesBf[nccNbrPatchi][nccPatchFacei];
+                const label nbrCelli =
+                    completeMesh().faceOwner()[nbrFacei];
+                const label nbrProci = cellProc_[nbrCelli];
+
+                const label procNccPatchi =
+                    refPatchProcPatchTable
+                    [nccPatchi][labelPair(proci, nbrProci)];
+                const label nbrProcNccPatchi =
+                    refPatchProcPatchTable
+                    [nccNbrPatchi][labelPair(nbrProci, proci)];
+
+                result[proci][procNccPatchi].append(nccPatchFacei + 1);
+                result[nbrProci][nbrProcNccPatchi].append(nccPatchFacei + 1);
+            }
+        }
+
+        // Error patches
+        forAll(completeMesh().boundary(), ncePatchi)
+        {
+            const fvPatch& fvp = completeMesh().boundary()[ncePatchi];
+
+            if (!isA<nonConformalErrorFvPatch>(fvp)) continue;
+
+            forAll(polyFacesBf[ncePatchi], ncePatchFacei)
+            {
+                const label facei = polyFacesBf[ncePatchi][ncePatchFacei];
+                const label celli = completeMesh().faceOwner()[facei];
+                const label proci = cellProc_[celli];
+
+                result[proci][ncePatchi].append(ncePatchFacei + 1);
+            }
+        }
+    }
+    else // if (!procsConformal())
+    {
+        // Reconstruct non-conformal addressing
+
+        // Cyclic patches
+        forAll(completeMesh().boundary(), nccPatchi)
+        {
+            const fvPatch& fvp = completeMesh().boundary()[nccPatchi];
+
+            if (!isA<nonConformalCyclicFvPatch>(fvp)) continue;
+
+            const nonConformalCyclicFvPatch& nccFvp =
+                refCast<const nonConformalCyclicFvPatch>(fvp);
+
+            if (!nccFvp.owner()) continue;
+
+            const label nccNbrPatchi = nccFvp.nbrPatchID();
+
+            label nccPatchFacei = 0;
+            labelPairTable procNccPatchFaceis;
+            forAllConstIter
+            (
+                labelPairTable,
+                refPatchProcPatchTable[nccPatchi],
+                iter
+            )
+            {
+                procNccPatchFaceis.insert(iter.key(), 0);
+            }
+
+            while (true)
+            {
+                labelPair procNbrProc(labelMax, labelMax);
+                labelPair faceNbrFace(labelMax, labelMax);
+
+                forAllConstIter(labelPairTable, procNccPatchFaceis, iter)
+                {
+                    const label proci = iter.key().first();
+                    const label nbrProci = iter.key().second();
+
+                    const labelPair procNbrProcStar(proci, nbrProci);
+                    const labelPair nbrProcProcStar(nbrProci, proci);
+
+                    const label procNccPatchi =
+                        refPatchProcPatchTable
+                        [nccPatchi][procNbrProcStar];
+                    const label nbrProcNccPatchi =
+                        refPatchProcPatchTable
+                        [nccNbrPatchi][nbrProcProcStar];
+
+                    const label size =
+                        procMeshes_[proci]
+                       .polyFacesBf()[procNccPatchi]
+                       .size();
+
+                    if (iter() >= size) continue;
+
+                    const label procFacei =
+                        procMeshes_[proci].polyFacesBf()
+                        [procNccPatchi][iter()];
+                    const label nbrProcFacei =
+                        procMeshes_[nbrProci].polyFacesBf()
+                        [nbrProcNccPatchi][iter()];
+
+                    const labelPair faceNbrFaceStar
+                    (
+                        procFaceAddressing_[proci][procFacei] - 1,
+                        procFaceAddressing_[nbrProci][nbrProcFacei] - 1
+                    );
+
+                    if (faceNbrFace > faceNbrFaceStar)
+                    {
+                        procNbrProc = procNbrProcStar;
+                        faceNbrFace = faceNbrFaceStar;
+                    }
+                }
+
+                if (faceNbrFace == labelPair(labelMax, labelMax))
+                {
+                    break;
+                }
+                else
+                {
+                    const label proci = procNbrProc.first();
+                    const label nbrProci = procNbrProc.second();
+
+                    const labelPair nbrProcProc(nbrProci, proci);
+
+                    const label procNccPatchi =
+                        refPatchProcPatchTable[nccPatchi][procNbrProc];
+                    const label nbrProcNccPatchi =
+                        refPatchProcPatchTable[nccNbrPatchi][nbrProcProc];
+
+                    result[proci][procNccPatchi]
+                        .append(nccPatchFacei + 1);
+                    result[nbrProci][nbrProcNccPatchi]
+                        .append(nccPatchFacei + 1);
+
+                    nccPatchFacei ++;
+                    procNccPatchFaceis[procNbrProc] ++;
+                }
+            }
+        }
+
+        // Error patches
+        forAll(completeMesh().boundary(), ncePatchi)
+        {
+            const fvPatch& fvp = completeMesh().boundary()[ncePatchi];
+
+            if (!isA<nonConformalErrorFvPatch>(fvp)) continue;
+
+            label ncePatchFacei = 0;
+            labelList procNcePatchFaceis(nProcs(), 0);
+
+            while (true)
+            {
+                label facei = labelMax, proci = labelMax;
+
+                forAll(procNcePatchFaceis, procStari)
+                {
+                    const label size =
+                        procMeshes_[procStari]
+                       .polyFacesBf()[ncePatchi]
+                       .size();
+
+                    if (procNcePatchFaceis[procStari] >= size) continue;
+
+                    const label procFacei =
+                        procMeshes_[procStari].polyFacesBf()
+                        [ncePatchi][procNcePatchFaceis[procStari]];
+
+                    const label faceStari =
+                        procFaceAddressing_[procStari][procFacei] - 1;
+
+                    if (facei > faceStari)
+                    {
+                        facei = faceStari;
+                        proci = procStari;
+                    }
+                }
+
+                if (facei == labelMax)
+                {
+                    break;
+                }
+                else
+                {
+                    result[proci][ncePatchi].append(ncePatchFacei + 1);
+
+                    ncePatchFacei ++;
+                    procNcePatchFaceis[proci] ++;
+                }
+            }
+        }
+    }
+
+    return result;
+}
+
+
+void Foam::domainDecomposition::unconformComplete()
+{
+    surfaceLabelField::Boundary polyFacesBf
+    (
+        surfaceLabelField::null(),
+        completeMesh().polyFacesBf()
+    );
+    surfaceVectorField Sf(completeMesh().Sf().cloneUnSliced());
+    surfaceVectorField Cf(completeMesh().Cf().cloneUnSliced());
+
+    forAll(procMeshes_, proci)
+    {
+        const fvMesh& procMesh = procMeshes_[proci];
+
+        const surfaceLabelField::Boundary& faceAddressingBf =
+            procFaceAddressingBf()[proci];
+
+        forAll(procMesh.boundary(), procNccPatchi)
+        {
+            const fvPatch& fvp = procMesh.boundary()[procNccPatchi];
+
+            if (isA<nonConformalFvPatch>(fvp))
+            {
+                const label completeNccPatchi =
+                    isA<processorCyclicFvPatch>(fvp)
+                  ? refCast<const processorCyclicFvPatch>(fvp)
+                   .referPatchID()
+                  : procNccPatchi;
+
+                const label size =
+                    max
+                    (
+                        max(mag(faceAddressingBf[procNccPatchi])),
+                        polyFacesBf[completeNccPatchi].size()
+                    );
+
+                polyFacesBf[completeNccPatchi].resize(size, -1);
+                polyFacesBf[completeNccPatchi].labelField::rmap
+                (
+                    mag
+                    (
+                        labelField
+                        (
+                            procFaceAddressing_[proci],
+                            procMesh.polyFacesBf()[procNccPatchi]
+                        )
+                    ) - 1,
+                    mag(faceAddressingBf[procNccPatchi]) - 1
+                );
+
+                // Set dummy data for the face geometry. This should not be
+                // used during decomposition.
+                Sf.boundaryFieldRef()[completeNccPatchi].resize(size, Zero);
+                Cf.boundaryFieldRef()[completeNccPatchi].resize(size, Zero);
+            }
+        }
+    }
+
+    completeMesh_->unconform
+    (
+        polyFacesBf,
+        Sf,
+        Cf,
+        NullObjectRef<surfaceScalarField>(),
+        false
+    );
+
+    // Check ordering
+    if (debug)
+    {
+        forAll(completeMesh().boundary(), patchi)
+        {
+            const fvPatch& fvp = completeMesh().boundary()[patchi];
+
+            // Coupled patches
+            if
+            (
+                isA<nonConformalCyclicFvPatch>(fvp)
+             && refCast<const nonConformalCoupledFvPatch>(fvp).owner()
+            )
+            {
+                const label nccPatchi = patchi;
+
+                const label nbrNccPatchi =
+                    refCast<const nonConformalCyclicFvPatch>(fvp)
+                   .nbrPatchID();
+
+                checkNonConformalCoupledPatchOrdering
+                (
+                    {-labelMax, labelMax},
+                    completeMesh().boundary()[nccPatchi],
+                    completeMesh().boundary()[nbrNccPatchi],
+                    completeMesh().polyFacesBf()[nccPatchi],
+                    completeMesh().polyFacesBf()[nbrNccPatchi]
+                );
+            }
+
+            // Error patches
+            if (isA<nonConformalErrorFvPatch>(fvp))
+            {
+                const label ncePatchi = patchi;
+
+                checkNonConformalErrorPatchOrdering
+                (
+                    -labelMax,
+                    completeMesh().boundary()[ncePatchi],
+                    completeMesh().polyFacesBf()[ncePatchi]
+                );
+            }
+        }
+    }
+}
+
+
+void Foam::domainDecomposition::unconformProcs()
+{
+    // Construct the reverse of proc-face-face addressing. -1 indicates a face
+    // that is on a (conformal) processor boundary and hence has multiple
+    // associated proc-face indices.
+    labelList faceProcFace(completeMesh().nFaces(), -labelMax);
+    forAll(procMeshes_, proci)
+    {
+        forAll(procFaceAddressing()[proci], procFacei)
+        {
+            const label facei =
+                mag(procFaceAddressing()[proci][procFacei]) - 1;
+
+            faceProcFace[facei] =
+                faceProcFace[facei] == -labelMax ? procFacei : -1;
+        }
+    }
+
+    forAll(procMeshes_, proci)
+    {
+        fvMesh& procMesh = procMeshes_[proci];
+
+        const surfaceLabelField::Boundary& faceAddressingBf =
+            procFaceAddressingBf()[proci];
+
+        surfaceLabelField::Boundary polyFacesBf
+        (
+            surfaceLabelField::null(),
+            procMesh.polyFacesBf()
+        );
+        surfaceVectorField Sf(procMesh.Sf().cloneUnSliced());
+        surfaceVectorField Cf(procMesh.Cf().cloneUnSliced());
+
+        forAll(procMesh.boundary(), procNccPatchi)
+        {
+            const fvPatch& fvp = procMesh.boundary()[procNccPatchi];
+
+            if (isA<nonConformalFvPatch>(fvp))
+            {
+                const label completeNccPatchi =
+                    isA<processorCyclicFvPatch>(fvp)
+                  ? refCast<const processorCyclicFvPatch>(fvp)
+                   .referPatchID()
+                  : procNccPatchi;
+
+                polyFacesBf[procNccPatchi] =
+                    labelField
+                    (
+                        faceProcFace,
+                        labelField
+                        (
+                            completeMesh().polyFacesBf()[completeNccPatchi],
+                            mag(faceAddressingBf[procNccPatchi]) - 1
+                        )
+                    );
+
+                const label size = polyFacesBf[procNccPatchi].size();
+
+                // Set dummy data for the face geometry. This should not be
+                // used during decomposition.
+                Sf.boundaryFieldRef()[procNccPatchi].resize(size, Zero);
+                Cf.boundaryFieldRef()[procNccPatchi].resize(size, Zero);
+            }
+        }
+
+        procMesh.unconform
+        (
+            polyFacesBf,
+            Sf,
+            Cf,
+            NullObjectRef<surfaceScalarField>(),
+            false
+        );
+    }
+
+    if (debug)
+    {
+        forAll(procMeshes_, proci)
+        {
+            forAll(procMeshes_[proci].boundary(), patchi)
+            {
+                const fvPatch& fvp =
+                    procMeshes_[proci].boundary()[patchi];
+
+                // Coupled patches
+                if
+                (
+                    isA<nonConformalCoupledFvPatch>(fvp)
+                 && refCast<const nonConformalCoupledFvPatch>(fvp).owner()
+                )
+                {
+                    const label nccPatchi = patchi;
+
+                    label nbrProci = -1, nbrNccPatchi = -1;
+                    if (isA<cyclicFvPatch>(fvp))
+                    {
+                        nbrProci = proci;
+                        nbrNccPatchi =
+                            refCast<const cyclicFvPatch>(fvp).nbrPatchID();
+                    }
+                    else if (isA<processorCyclicFvPatch>(fvp))
+                    {
+                        typedef processorCyclicFvPatch PcFvp;
+
+                        const PcFvp& pcFvp = refCast<const PcFvp>(fvp);
+
+                        nbrProci = pcFvp.neighbProcNo();
+
+                        const fvBoundaryMesh& nbrFvPatches =
+                            procMeshes_[nbrProci].boundary();
+
+                        forAll(nbrFvPatches, nbrNccPatchj)
+                        {
+                            const fvPatch& nbrFvp =
+                                nbrFvPatches[nbrNccPatchj];
+
+                            if (isA<PcFvp>(nbrFvp))
+                            {
+                                const PcFvp& nbrPcFvp =
+                                    refCast<const PcFvp>(nbrFvp);
+
+                                if
+                                (
+                                    nbrPcFvp.neighbProcNo()
+                                 == proci
+                                 && nbrPcFvp.referPatchID()
+                                 == pcFvp.referPatch().nbrPatchID()
+                                )
+                                {
+                                    nbrNccPatchi = nbrNccPatchj;
+                                    break;
+                                }
+                            }
+                        }
+
+                        if (nbrNccPatchi == -1)
+                        {
+                            FatalErrorInFunction
+                                << "Opposite processor patch not found for "
+                                << "patch " << fvp.name() << " on proc #"
+                                << proci << exit(FatalError);
+                        }
+                    }
+                    else
+                    {
+                        FatalErrorInFunction
+                            << "Non-conformal-coupled type not recognised "
+                            << "for patch " << fvp.name() << " on proc #"
+                            << proci << exit(FatalError);
+                    }
+
+                    checkNonConformalCoupledPatchOrdering
+                    (
+                        {proci, nbrProci},
+                        procMeshes_[proci].boundary()[nccPatchi],
+                        procMeshes_[nbrProci].boundary()[nbrNccPatchi],
+                        procMeshes_[proci].polyFacesBf()[nccPatchi],
+                        procMeshes_[nbrProci].polyFacesBf()[nbrNccPatchi]
+                    );
+                }
+
+                // Error patches
+                if (isA<nonConformalErrorFvPatch>(fvp))
+                {
+                    const label ncePatchi = patchi;
+
+                    checkNonConformalErrorPatchOrdering
+                    (
+                        proci,
+                        procMeshes_[proci].boundary()[ncePatchi],
+                        procMeshes_[proci].polyFacesBf()[ncePatchi]
+                    );
+                }
+            }
+        }
+    }
+}
+
+
+void Foam::domainDecomposition::unconform()
+{
+    if (completeConformal() && procsConformal())
+    {
+        // Nothing to do
+    }
+    else if (!completeConformal() && procsConformal())
+    {
+        unconformProcs();
+    }
+    else if (completeConformal() && !procsConformal())
+    {
+        unconformComplete();
+    }
+    else // if (!completeConformal() && !procsConformal())
+    {
+        // Assume everything is consistent and do nothing
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/parallel/parallel/domainDecompositionReconstruct.C b/src/parallel/parallel/domainDecompositionReconstruct.C
index 2a478e9a01..a36df6d11f 100644
--- a/src/parallel/parallel/domainDecompositionReconstruct.C
+++ b/src/parallel/parallel/domainDecompositionReconstruct.C
@@ -28,12 +28,10 @@ License
 #include "processorPolyPatch.H"
 #include "processorCyclicPolyPatch.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-namespace Foam
-{
-
-autoPtr<faceCoupleInfo> determineCoupledFaces
+Foam::autoPtr<Foam::faceCoupleInfo>
+Foam::domainDecomposition::determineCoupledFaces
 (
     const label masterMeshProcStart,
     const label masterMeshProcEnd,
@@ -136,9 +134,6 @@ autoPtr<faceCoupleInfo> determineCoupledFaces
     );
 }
 
-}
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void Foam::domainDecomposition::reconstruct()
 {
@@ -225,7 +220,7 @@ void Foam::domainDecomposition::reconstruct()
             Info<< indent << "Merging mesh " << proci
                 << " with " << procj << endl;
 
-            // Find shared points/faces
+            // Find shared points and faces
             autoPtr<faceCoupleInfo> couples = determineCoupledFaces
             (
                 proci,
@@ -521,7 +516,32 @@ void Foam::domainDecomposition::reconstruct()
     completeMesh_->setInstance(procMeshes()[0].facesInstance());
     completeMesh_->setPointsInstance(procMeshes()[0].pointsInstance());
 
-    Info<< decrIndent << endl;
+    Info<< decrIndent;
+}
+
+
+void Foam::domainDecomposition::reconstructPoints()
+{
+    const label pointsCompare =
+        compareInstances
+        (
+            completeMesh().pointsInstance(),
+            procMeshes_[0].pointsInstance()
+        );
+
+    if (pointsCompare == 1)
+    {
+        pointField completePoints(completeMesh().nPoints());
+
+        for (label proci = 0; proci < nProcs(); proci++)
+        {
+            const fvMesh& procMesh = procMeshes_[proci];
+
+            completePoints.rmap(procMesh.points(), procPointAddressing_[proci]);
+        }
+
+        completeMesh_->setPoints(completePoints);
+    }
 }
 
 
diff --git a/src/parallel/parallel/multiDomainDecomposition.C b/src/parallel/parallel/multiDomainDecomposition.C
new file mode 100644
index 0000000000..41ac00ca71
--- /dev/null
+++ b/src/parallel/parallel/multiDomainDecomposition.C
@@ -0,0 +1,180 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2023 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "multiDomainDecomposition.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(multiDomainDecomposition, 0);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::multiDomainDecomposition::multiDomainDecomposition
+(
+    const processorRunTimes& runTimes,
+    const wordList& regionNames
+)
+:
+    multiRegionPrefixer(false, regionNames),
+    runTimes_(runTimes),
+    regionNames_(regionNames),
+    regionMeshes_(regionNames.size())
+{
+    forAll(regionMeshes_, regioni)
+    {
+        regionMeshes_.set
+        (
+            regioni,
+            new domainDecomposition(runTimes, regionNames[regioni])
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::multiDomainDecomposition::~multiDomainDecomposition()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::multiDomainDecomposition::readDecompose(const bool doSets)
+{
+    bool result = false;
+
+    forAll(regionMeshes_, regioni)
+    {
+        if (meshes(regioni)().readDecompose(doSets))
+        {
+            result = true;
+
+            if (regioni != regionMeshes_.size() - 1)
+            {
+                Info<< endl;
+            }
+        }
+    }
+
+    return result;
+}
+
+
+bool Foam::multiDomainDecomposition::readReconstruct(const bool doSets)
+{
+    bool result = false;
+
+    forAll(regionMeshes_, regioni)
+    {
+        if (meshes(regioni)().readReconstruct(doSets))
+        {
+            result = true;
+
+            if (regioni != regionMeshes_.size() - 1)
+            {
+                Info<< endl;
+            }
+        }
+    }
+
+    return result;
+}
+
+
+Foam::fvMesh::readUpdateState
+Foam::multiDomainDecomposition::readUpdateDecompose()
+{
+    fvMesh::readUpdateState result = fvMesh::UNCHANGED;
+
+    forAll(regionMeshes_, regioni)
+    {
+        const fvMesh::readUpdateState regionResult =
+            meshes(regioni)().readUpdateDecompose();
+
+        if
+        (
+            regioni != regionMeshes_.size() - 1
+         && regionResult >= fvMesh::TOPO_CHANGE
+        )
+        {
+            Info<< endl;
+        }
+
+        result = result > regionResult ? result : regionResult;
+    }
+
+    return result;
+}
+
+
+Foam::fvMesh::readUpdateState
+Foam::multiDomainDecomposition::readUpdateReconstruct()
+{
+    fvMesh::readUpdateState result = fvMesh::UNCHANGED;
+
+    forAll(regionMeshes_, regioni)
+    {
+        const fvMesh::readUpdateState regionResult =
+            meshes(regioni)().readUpdateReconstruct();
+
+        if
+        (
+            regioni != regionMeshes_.size() - 1
+         && regionResult >= fvMesh::TOPO_CHANGE
+        )
+        {
+            Info<< endl;
+        }
+
+        result = result > regionResult ? result : regionResult;
+    }
+
+    return result;
+}
+
+
+void Foam::multiDomainDecomposition::writeComplete(const bool doSets) const
+{
+    forAll(regionMeshes_, regioni)
+    {
+        meshes(regioni)().writeComplete(doSets);
+    }
+}
+
+
+void Foam::multiDomainDecomposition::writeProcs(const bool doSets) const
+{
+    forAll(regionMeshes_, regioni)
+    {
+        meshes(regioni)().writeProcs(doSets);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/parallel/parallel/multiDomainDecomposition.H b/src/parallel/parallel/multiDomainDecomposition.H
new file mode 100644
index 0000000000..ba58a51786
--- /dev/null
+++ b/src/parallel/parallel/multiDomainDecomposition.H
@@ -0,0 +1,156 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2023 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::multiDomainDecomposition
+
+Description
+    ...
+
+SourceFiles
+    multiDomainDecomposition.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef multiDomainDecomposition_H
+#define multiDomainDecomposition_H
+
+#include "domainDecomposition.H"
+#include "multiRegionPrefixer.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class multiDomainDecomposition Declaration
+\*---------------------------------------------------------------------------*/
+
+class multiDomainDecomposition
+:
+    private multiRegionPrefixer
+{
+    // Private Data
+
+        //- Run times
+        const processorRunTimes& runTimes_;
+
+        //- Region names
+        const wordList& regionNames_;
+
+        //- The region complete and processor meshes
+        PtrList<domainDecomposition> regionMeshes_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("multiDomainDecomposition");
+
+
+    // Constructors
+
+        //- Construct from processor run times and region names
+        multiDomainDecomposition
+        (
+            const processorRunTimes& procRunTimes,
+            const wordList& regionNames
+        );
+
+
+    //- Destructor
+    virtual ~multiDomainDecomposition();
+
+
+    // Member Functions
+
+        //- Access the meshes for a region
+        inline RegionConstRef<domainDecomposition> meshes
+        (
+            const label regioni
+        ) const
+        {
+            return
+                RegionConstRef<domainDecomposition>
+                (
+                    *this,
+                    regioni,
+                    regionMeshes_[regioni]
+                );
+        }
+
+        //- Access the meshes for a region
+        inline RegionRef<domainDecomposition> meshes
+        (
+            const label regioni
+        )
+        {
+            return
+                RegionRef<domainDecomposition>
+                (
+                    *this,
+                    regioni,
+                    regionMeshes_[regioni]
+                );
+        }
+
+        //- Return the number of processors in the decomposition
+        inline label nProcs() const
+        {
+            return runTimes_.nProcs();
+        }
+
+        //- Read in the complete mesh. Read the processor meshes if they are
+        //  available and up to date relative to the complete mesh, or
+        //  decompose if not. Return whether or not decomposition happened.
+        bool readDecompose(const bool doSets);
+
+        //- Read in the processor meshes. Read the complete mesh if it is
+        //  available and up to date relative to the processor meshes, or
+        //  reconstruct if not. Return whether or not reconstruction happened.
+        bool readReconstruct(const bool doSets);
+
+        //- Read-update for decomposition
+        fvMesh::readUpdateState readUpdateDecompose();
+
+        //- Read-update for reconstruction
+        fvMesh::readUpdateState readUpdateReconstruct();
+
+        //- Write the decomposed meshes and associated data
+        void writeComplete(const bool doSets) const;
+
+        //- Write the decomposed meshes and associated data
+        void writeProcs(const bool doSets) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/parallel/parallel/multiRegionPrefixer.C b/src/parallel/parallel/multiRegionPrefixer.C
new file mode 100644
index 0000000000..a48dc1f105
--- /dev/null
+++ b/src/parallel/parallel/multiRegionPrefixer.C
@@ -0,0 +1,123 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2023 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "multiRegionPrefixer.H"
+#include "IOstreams.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+bool Foam::multiRegionPrefixer::prefixes() const
+{
+    return prefixSingleRegion_ || regionNames_.size() > 1;
+}
+
+
+Foam::string::size_type Foam::multiRegionPrefixer::prefixWidth() const
+{
+    string::size_type n = 0;
+
+    if (prefixes())
+    {
+        forAll(regionNames_, regionj)
+        {
+            n = max(n, regionNames_[regionj].size() + 1);
+        }
+    }
+
+    return n;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::multiRegionPrefixer::regionPrefixer::regionPrefixer
+(
+    const multiRegionPrefixer& mrp,
+    const label regioni
+)
+:
+    mrp_(mrp),
+    regioni_(regioni)
+{
+    if (mrp_.prefixes() && regioni_ != -1)
+    {
+        const string::size_type dn =
+            mrp_.prefixWidth()
+          - mrp_.regionNames_[regioni].size();
+
+        Sout.prefix() =
+            mrp_.regionNames_[regioni] + string(dn, ' ');
+    }
+}
+
+
+Foam::multiRegionPrefixer::regionPrefixer::regionPrefixer
+(
+    regionPrefixer&& rp
+)
+:
+    regionPrefixer(rp.mrp_, rp.regioni_)
+{
+    rp.regioni_ = -1;
+}
+
+
+Foam::multiRegionPrefixer::multiRegionPrefixer
+(
+    const bool prefixSingleRegion,
+    const wordList& regionNames
+)
+:
+    prefixSingleRegion_(prefixSingleRegion),
+    regionNames_(regionNames)
+{
+    if (prefixes())
+    {
+        Sout.prefix() = string(prefixWidth(), ' ');
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
+
+Foam::multiRegionPrefixer::regionPrefixer::~regionPrefixer()
+{
+    if (mrp_.prefixes() && regioni_ != -1)
+    {
+        Sout.prefix() = string(mrp_.prefixWidth(), ' ');
+    }
+}
+
+
+Foam::multiRegionPrefixer::~multiRegionPrefixer()
+{
+    if (prefixes())
+    {
+        Sout.prefix() = string::null;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/parallel/parallel/multiRegionPrefixer.H b/src/parallel/parallel/multiRegionPrefixer.H
new file mode 100644
index 0000000000..8ee2cf0c17
--- /dev/null
+++ b/src/parallel/parallel/multiRegionPrefixer.H
@@ -0,0 +1,241 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2023 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::multiRegionPrefixer
+
+SourceFiles
+    multiRegionPrefixer.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef multiRegionPrefixer_H
+#define multiRegionPrefixer_H
+
+#include "wordList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class multiRegionPrefixer Declaration
+\*---------------------------------------------------------------------------*/
+
+class multiRegionPrefixer
+{
+private:
+
+    // Private Data
+
+        //- Should we prefix a single region?
+        const bool prefixSingleRegion_;
+
+        //- Region names
+        const wordList& regionNames_;
+
+
+    // Private Member Functions
+
+        //- Are we prefixing?
+        bool prefixes() const;
+
+        //- Width of the write prefix
+        string::size_type prefixWidth() const;
+
+
+public:
+
+    // Public classes
+
+        //- Prefixer for a single region
+        class regionPrefixer
+        {
+        private:
+
+            // Private data
+
+                //- The multi-region prefixing engine
+                const multiRegionPrefixer& mrp_;
+
+                //- The region index
+                label regioni_;
+
+
+        public:
+
+            // Constructors
+
+                //- Construct from components
+                regionPrefixer
+                (
+                    const multiRegionPrefixer& mrp,
+                    const label regioni
+                );
+
+                //- Disallow copy construct
+                regionPrefixer(const regionPrefixer& rp) = delete;
+
+                //- Move construct
+                regionPrefixer(regionPrefixer&& rp);
+
+
+            //- Destructor
+            ~regionPrefixer();
+
+
+            // Member operators
+
+                //- Disallow copy assign
+                regionPrefixer& operator=(const regionPrefixer&) = delete;
+
+                //- Disallow move assign
+                regionPrefixer& operator=(regionPrefixer&&) = delete;
+        };
+
+
+    // Constructors
+
+        //- Construct from components
+        multiRegionPrefixer
+        (
+            const bool prefixSingleRegion,
+            const wordList& regionNames
+        );
+
+
+    //- Destructor
+    ~multiRegionPrefixer();
+};
+
+
+/*---------------------------------------------------------------------------*\
+                         Class RegionConstRef Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Region>
+class RegionConstRef
+:
+    private multiRegionPrefixer::regionPrefixer
+{
+private:
+
+    // Private Data
+
+        //- Reference to the region object
+        const Region& r_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct form components
+        inline RegionConstRef
+        (
+            const multiRegionPrefixer& mrp,
+            const label regioni,
+            const Region& r
+        )
+        :
+            multiRegionPrefixer::regionPrefixer(mrp, regioni),
+            r_(r)
+        {}
+
+
+    // Member operators
+
+        //- Cast to reference
+        inline operator const Region&() const
+        {
+            return r_;
+        }
+
+        //- Obtain the reference
+        inline const Region& operator()() const
+        {
+            return r_;
+        }
+};
+
+
+/*---------------------------------------------------------------------------*\
+                         Class RegionRef Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Region>
+class RegionRef
+:
+    public RegionConstRef<Region>
+{
+private:
+
+    // Private Data
+
+        //- Reference to the region object
+        Region& r_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct form components
+        inline RegionRef
+        (
+            const multiRegionPrefixer& mrp,
+            const label regioni,
+            Region& r
+        )
+        :
+            RegionConstRef<Region>(mrp, regioni, r),
+            r_(r)
+        {}
+
+
+    // Member operators
+
+        //- Cast to reference
+        inline operator Region&() const
+        {
+            return r_;
+        }
+
+        //- Obtain the reference
+        inline Region& operator()()
+        {
+            return r_;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //