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