diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options
index 2fcfff41267da303f74240ace88b730e728f87c0..c896f86d1654b9687670bc9339d4bb3d8c323deb 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options
@@ -8,7 +8,9 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basicSolidThermo/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude
+    -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \
+    -I$(LIB_SRC)/parallel/decompose/decompose/lnInclude \
+    -I$(LIB_SRC)/parallel/reconstruct/reconstruct/lnInclude
 
 EXE_LIBS = \
     -lbasicThermophysicalModels \
@@ -19,4 +21,6 @@ EXE_LIBS = \
     -lcompressibleLESModels \
     -lmeshTools \
     -lfiniteVolume \
-    -lradiationModels
+    -lradiationModels \
+    -ldecompose \
+    -lreconstruct
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/allhEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/allhEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..cfe2e84a0f77533303ec5bc569d42c18709fcab1
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/allhEqn.H
@@ -0,0 +1,84 @@
+
+// Get mapped alpha (surfaceScalarField)
+tmp<surfaceScalarField> tallAlpha = procToAllMapper().reconstructFvSurfaceField
+    (
+        IOobject
+        (
+            "alpha",
+            allMesh().time().timeName(),
+            allMesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        procAlpha
+    );
+
+// Get alpha from harmonic interpolation of vol quantities
+// (Note: really only needed at internal faces originating patches
+//  inbetween regions)
+tmp<surfaceScalarField> allHarmonicAlpha
+(
+    harmonic(allMesh()).interpolate(allVolAlpha())
+);
+
+// Loop over all fluid and solid regions to transfer
+// allHarmonicAlpha to allAlpha
+surfaceScalarField& allAlpha = tallAlpha();
+forAll(boundaryProcAddressing, procI)
+{
+    forAll(boundaryProcAddressing[procI], patchI)
+    {
+        if (boundaryProcAddressing[procI][patchI] == -1)
+        {
+            // Interface patch
+            const labelList::subList cp =
+                procMeshes[procI].boundary()[patchI].patchSlice
+                (
+                    faceProcAddressing[procI]
+                );
+
+            forAll(cp, faceI)
+            {
+                label curF = mag(cp[faceI])-1;
+                if (curF < allMesh().nInternalFaces())
+                {
+                    allAlpha[curF] = allHarmonicAlpha()[curF];
+                }
+            }
+        }
+    }
+}
+
+
+tmp<surfaceScalarField> allPhi
+(
+    procToAllMapper().reconstructFvSurfaceField
+    (
+        IOobject
+        (
+            "phi",
+            allMesh().time().timeName(),
+            allMesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        procPhi
+    )
+);
+
+// So we have nNonOrthCorr
+//#include "readSolidMultiRegionPIMPLEControls.H"
+//for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+{
+    fvScalarMatrix hEqn
+    (
+        fvm::ddt(allRho(), allh())
+      + fvm::div(allPhi(), allh())
+      - fvm::laplacian(allAlpha, allh())
+     ==
+        allSource()
+    );
+
+    hEqn.relax();
+    hEqn.solve(allMesh().solver(allh().select(finalIter)));
+}
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
index b65ecf5118e0908642ff045171af0f2dad0e58a9..927562c144dacf0a83f2ab5c26d9f76e9e1eaab2 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
@@ -39,6 +39,11 @@ Description
 #include "solidRegionDiffNo.H"
 #include "basicSolidThermo.H"
 #include "radiationModel.H"
+#include "fvFieldReconstructor.H"
+#include "mixedFvPatchFields.H"
+#include "fvFieldDecomposer.H"
+#include "harmonic.H"
+#include "rmap.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -49,12 +54,40 @@ int main(int argc, char *argv[])
 
     regionProperties rp(runTime);
 
+    const label nAllRegions =
+        rp.fluidRegionNames().size()
+      + rp.solidRegionNames().size();
+    PtrList<labelIOList> cellProcAddressing(nAllRegions);
+    PtrList<labelIOList> faceProcAddressing(nAllRegions);
+    PtrList<labelIOList> boundaryProcAddressing(nAllRegions);
+    PtrList<fvMesh> procMeshes(nAllRegions);
+
+    // Load meshes, fluid first
+    labelList fluidToProc(identity(rp.fluidRegionNames().size()));
+    labelList solidToProc(rp.solidRegionNames().size());
+    forAll(solidToProc, i)
+    {
+        solidToProc[i] = fluidToProc.size()+i;
+    }
+
+    // Get the coupled solution flag
+    #include "readPIMPLEControls.H"
+
+    if (temperatureCoupled)
+    {
+        Info<< "Solving single enthalpy for all equations" << nl << endl;
+    }
+
     #include "createFluidMeshes.H"
     #include "createSolidMeshes.H"
 
     #include "createFluidFields.H"
     #include "createSolidFields.H"
 
+    // Temperature solved on single mesh
+    #include "createAllMesh.H"
+    #include "createAllFields.H"
+
     #include "initContinuityErrs.H"
 
     #include "readTimeControls.H"
@@ -81,9 +114,10 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
+
         if (nOuterCorr != 1)
         {
-            forAll(fluidRegions, i)
+            forAll(fluidToProc, i)
             {
                 #include "setRegionFluidFields.H"
                 #include "storeOldFluidFields.H"
@@ -96,22 +130,133 @@ int main(int argc, char *argv[])
         {
             bool finalIter = oCorr == nOuterCorr-1;
 
-            forAll(fluidRegions, i)
+            if (finalIter)
+            {
+                forAll(procMeshes, procI)
+                {
+                    procMeshes[procI].data::add("finalIteration", true);
+                }
+            }
+
+
+            PtrList<surfaceScalarField> procPhi(nAllRegions);
+            PtrList<surfaceScalarField> procAlpha(nAllRegions);
+
+
+            // Solve (uncoupled) or set up (coupled) the temperature equation
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+            forAll(solidToProc, i)
+            {
+                label procI = solidToProc[i];
+
+                Info<< "\nSolving temperature for solid region "
+                    << procMeshes[procI].name() << endl;
+                #include "setRegionSolidFields.H"
+                #include "readSolidMultiRegionPIMPLEControls.H"
+
+                if (temperatureCoupled)
+                {
+                    // Map my properties to overall h equation
+                    #include "rmapSolid.H"
+                }
+                else
+                {
+                    #include "solveSolid.H"
+                }
+            }
+
+
+            forAll(fluidToProc, i)
+            {
+                label procI = fluidToProc[i];
+
+                Info<< "\nSolving temperature for fluid region "
+                    << procMeshes[procI].name() << endl;
+                #include "setRegionFluidFields.H"
+                #include "readFluidMultiRegionPIMPLEControls.H"
+
+                if (oCorr == 0)
+                {
+                    #include "rhoEqn.H"
+                }
+
+                if (temperatureCoupled)
+                {
+                    // Map my properties to overall h equation
+                    #include "rmapFluid.H"
+                }
+                else
+                {
+                    #include "hEqn.H"
+                }
+            }
+
+
+            // Solve combined h equation
+            // ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+            if (temperatureCoupled)
             {
-                Info<< "\nSolving for fluid region "
-                    << fluidRegions[i].name() << endl;
+                Info<< "\nSolving single enthalpy for all regions"
+                    << endl;
+
+                // Solve combined h
+                #include "allhEqn.H"
+
+                forAll(solidToProc, i)
+                {
+                    label procI = solidToProc[i];
+                    #include "setRegionSolidFields.H"
+                    #include "readSolidMultiRegionPIMPLEControls.H"
+
+                    #include "mapSolid.H"
+                }
+
+                forAll(fluidToProc, i)
+                {
+                    label procI = fluidToProc[i];
+                    #include "setRegionFluidFields.H"
+                    #include "readFluidMultiRegionPIMPLEControls.H"
+
+                    #include "mapFluid.H"
+                }
+            }
+
+
+            // Update thermos
+            // ~~~~~~~~~~~~~~
+
+            forAll(fluidToProc, i)
+            {
+                Info<< "\nUpdating thermo for fluid region "
+                    << procMeshes[fluidToProc[i]].name() << endl;
+
                 #include "setRegionFluidFields.H"
                 #include "readFluidMultiRegionPIMPLEControls.H"
-                #include "solveFluid.H"
+
+                thermo.correct();
+                rad.correct();
+                #include "solvePressureVelocityFluid.H"
             }
 
-            forAll(solidRegions, i)
+            forAll(solidToProc, i)
             {
-                Info<< "\nSolving for solid region "
-                    << solidRegions[i].name() << endl;
+                Info<< "\nUpdating thermo for solid region "
+                    << procMeshes[solidToProc[i]].name() << endl;
                 #include "setRegionSolidFields.H"
                 #include "readSolidMultiRegionPIMPLEControls.H"
-                #include "solveSolid.H"
+
+                thermo.correct();
+            }
+
+
+            if (finalIter)
+            {
+                forAll(procMeshes, procI)
+                {
+                    procMeshes[procI].data::remove("finalIteration");
+                }
             }
         }
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/createAllFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/createAllFields.H
new file mode 100644
index 0000000000000000000000000000000000000000..faa8eaac5872eded563c1d464978a13a2c874710
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/createAllFields.H
@@ -0,0 +1,64 @@
+    autoPtr<volScalarField> allRho;
+    autoPtr<volScalarField> allh;
+    autoPtr<volScalarField> allVolAlpha;
+    autoPtr<fvScalarMatrix> allSource;
+
+    if (temperatureCoupled)
+    {
+        allRho.reset
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    "rho",
+                    runTime.timeName(),
+                    allMesh(),
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                allMesh(),
+                dimensionedScalar("rho", dimDensity, 0.0)
+            )
+        );
+
+        allh.reset
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    "h",
+                    runTime.timeName(),
+                    allMesh(),
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                allMesh(),
+                dimensionedScalar("h", dimEnergy/dimMass, 0.0),
+                mixedFvPatchScalarField::typeName
+            )
+        );
+
+        allVolAlpha.reset
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    "volAlpha",
+                    runTime.timeName(),
+                    allMesh(),
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                allMesh(),
+                dimensionedScalar("volAlpha", dimMass/dimLength/dimTime, 0.0)
+            )
+        );
+
+        allSource.reset
+        (
+            new fvMatrix<scalar>(allh(), allh().dimensions()*dimMass/dimTime)
+        );
+    }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/createAllMesh.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/createAllMesh.H
new file mode 100644
index 0000000000000000000000000000000000000000..8bff40e48a02a4059f76a5e383a2c1316db55a2d
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/createAllMesh.H
@@ -0,0 +1,63 @@
+//
+// createAllMesh.H
+// ~~~~~~~~~~~~~~~
+
+    autoPtr<fvMesh> allMesh;
+    autoPtr<fvFieldReconstructor> procToAllMapper;
+    PtrList<fvFieldDecomposer> allToProcMappers;
+
+    if (temperatureCoupled)
+    {
+        Foam::Info
+            << "Create mesh for time = "
+            << runTime.timeName() << Foam::nl << Foam::endl;
+
+        allMesh.reset
+        (
+            new Foam::fvMesh
+            (
+                Foam::IOobject
+                (
+                    Foam::fvMesh::defaultRegion,
+                    runTime.timeName(),
+                    runTime,
+                    Foam::IOobject::MUST_READ
+                )
+            )
+        );
+
+        procToAllMapper.reset
+        (
+            new fvFieldReconstructor
+            (
+                allMesh(),
+                procMeshes,
+                faceProcAddressing,
+                cellProcAddressing,
+                boundaryProcAddressing
+            )
+        );
+
+
+        allToProcMappers.setSize
+        (
+            rp.fluidRegionNames().size()
+          + rp.solidRegionNames().size()
+        );
+
+        forAll(allToProcMappers, i)
+        {
+            allToProcMappers.set
+            (
+                i,
+                new fvFieldDecomposer
+                (
+                    allMesh(),
+                    procMeshes[i],
+                    faceProcAddressing[i],
+                    cellProcAddressing[i],
+                    boundaryProcAddressing[i]
+                )
+            );
+        }
+    }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H
index b5db5078f1ab81134f8efe86fbc6d67a8fab3984..61b7b1034089b01475bb2f00316c64012ae96ceb 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H
@@ -1,12 +1,12 @@
     scalar CoNum = -GREAT;
 
-    forAll(fluidRegions, regionI)
+    forAll(fluidToProc, regionI)
     {
         CoNum = max
         (
             compressibleCourantNo
             (
-                fluidRegions[regionI],
+                procMeshes[fluidToProc[regionI]],
                 runTime,
                 rhoFluid[regionI],
                 phiFluid[regionI]
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
index f6788d26eb3ecccb2bea1484949a781e7f3010a3..667b39d64631e31c1b6d69854f42025ee20dd2fc 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
@@ -1,30 +1,36 @@
     // Initialise fluid field pointer lists
-    PtrList<basicRhoThermo> thermoFluid(fluidRegions.size());
-    PtrList<volScalarField> rhoFluid(fluidRegions.size());
-    PtrList<volScalarField> KFluid(fluidRegions.size());
-    PtrList<volVectorField> UFluid(fluidRegions.size());
-    PtrList<surfaceScalarField> phiFluid(fluidRegions.size());
-    PtrList<uniformDimensionedVectorField> gFluid(fluidRegions.size());
-    PtrList<compressible::turbulenceModel> turbulence(fluidRegions.size());
-    PtrList<volScalarField> p_rghFluid(fluidRegions.size());
-    PtrList<volScalarField> ghFluid(fluidRegions.size());
-    PtrList<surfaceScalarField> ghfFluid(fluidRegions.size());
-    PtrList<radiation::radiationModel> radiation(fluidRegions.size());
-    PtrList<volScalarField> DpDtFluid(fluidRegions.size());
-
-    List<scalar> initialMassFluid(fluidRegions.size());
+    PtrList<basicRhoThermo> thermoFluid(rp.fluidRegionNames().size());
+    PtrList<volScalarField> rhoFluid(rp.fluidRegionNames().size());
+    PtrList<volScalarField> KFluid(rp.fluidRegionNames().size());
+    PtrList<volVectorField> UFluid(rp.fluidRegionNames().size());
+    PtrList<surfaceScalarField> phiFluid(rp.fluidRegionNames().size());
+    PtrList<uniformDimensionedVectorField> gFluid(rp.fluidRegionNames().size());
+    PtrList<compressible::turbulenceModel> turbulence
+    (
+        rp.fluidRegionNames().size()
+    );
+    PtrList<volScalarField> p_rghFluid(rp.fluidRegionNames().size());
+    PtrList<volScalarField> ghFluid(rp.fluidRegionNames().size());
+    PtrList<surfaceScalarField> ghfFluid(rp.fluidRegionNames().size());
+    PtrList<radiation::radiationModel> radiation(rp.fluidRegionNames().size());
+    PtrList<volScalarField> DpDtFluid(rp.fluidRegionNames().size());
+
+    List<scalar> initialMassFluid(rp.fluidRegionNames().size());
 
     // Populate fluid field pointer lists
-    forAll(fluidRegions, i)
+    forAll(rp.fluidRegionNames(), i)
     {
         Info<< "*** Reading fluid mesh thermophysical properties for region "
-            << fluidRegions[i].name() << nl << endl;
+            << rp.fluidRegionNames()[i] << nl << endl;
+
+        label procI = fluidToProc[i];
+
 
         Info<< "    Adding to thermoFluid\n" << endl;
         thermoFluid.set
         (
             i,
-            basicRhoThermo::New(fluidRegions[i]).ptr()
+            basicRhoThermo::New(procMeshes[procI]).ptr()
         );
 
         Info<< "    Adding to rhoFluid\n" << endl;
@@ -37,7 +43,7 @@
                 (
                     "rho",
                     runTime.timeName(),
-                    fluidRegions[i],
+                    procMeshes[procI],
                     IOobject::NO_READ,
                     IOobject::AUTO_WRITE
                 ),
@@ -55,7 +61,7 @@
                 (
                     "K",
                     runTime.timeName(),
-                    fluidRegions[i],
+                    procMeshes[procI],
                     IOobject::NO_READ,
                     IOobject::NO_WRITE
                 ),
@@ -73,11 +79,11 @@
                 (
                     "U",
                     runTime.timeName(),
-                    fluidRegions[i],
+                    procMeshes[procI],
                     IOobject::MUST_READ,
                     IOobject::AUTO_WRITE
                 ),
-                fluidRegions[i]
+                procMeshes[procI]
             )
         );
 
@@ -91,12 +97,12 @@
                 (
                     "phi",
                     runTime.timeName(),
-                    fluidRegions[i],
+                    procMeshes[procI],
                     IOobject::READ_IF_PRESENT,
                     IOobject::AUTO_WRITE
                 ),
                 linearInterpolate(rhoFluid[i]*UFluid[i])
-                    & fluidRegions[i].Sf()
+                    & procMeshes[procI].Sf()
             )
         );
 
@@ -110,7 +116,7 @@
                 (
                     "g",
                     runTime.constant(),
-                    fluidRegions[i],
+                    procMeshes[procI],
                     IOobject::MUST_READ,
                     IOobject::NO_WRITE
                 )
@@ -137,14 +143,14 @@
         ghFluid.set
         (
             i,
-            new volScalarField("gh", gFluid[i] & fluidRegions[i].C())
+            new volScalarField("gh", gFluid[i] & procMeshes[procI].C())
         );
 
         Info<< "    Adding to ghfFluid\n" << endl;
         ghfFluid.set
         (
             i,
-            new surfaceScalarField("ghf", gFluid[i] & fluidRegions[i].Cf())
+            new surfaceScalarField("ghf", gFluid[i] & procMeshes[procI].Cf())
         );
 
         p_rghFluid.set
@@ -156,11 +162,11 @@
                 (
                     "p_rgh",
                     runTime.timeName(),
-                    fluidRegions[i],
+                    procMeshes[procI],
                     IOobject::MUST_READ,
                     IOobject::AUTO_WRITE
                 ),
-                fluidRegions[i]
+                procMeshes[procI]
             )
         );
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidMeshes.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidMeshes.H
index 30a2e1089f8875cf507ac0d4d76492830cf9647c..73b533b0d91136a1a16d1395752373e034c43bb7 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidMeshes.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidMeshes.H
@@ -1,13 +1,13 @@
-    PtrList<fvMesh> fluidRegions(rp.fluidRegionNames().size());
-
     forAll(rp.fluidRegionNames(), i)
     {
         Info<< "Create fluid mesh for region " << rp.fluidRegionNames()[i]
             << " for time = " << runTime.timeName() << nl << endl;
 
-        fluidRegions.set
+        label procI = fluidToProc[i];
+
+        procMeshes.set
         (
-            i,
+            procI,
             new fvMesh
             (
                 IOobject
@@ -19,4 +19,58 @@
                 )
             )
         );
+
+        if (temperatureCoupled)
+        {
+            cellProcAddressing.set
+            (
+                procI,
+                new labelIOList
+                (
+                    IOobject
+                    (
+                        "cellRegionAddressing",
+                        procMeshes[procI].facesInstance(),
+                        procMeshes[procI].meshSubDir,
+                        procMeshes[procI],
+                        IOobject::MUST_READ,
+                        IOobject::NO_WRITE
+                    )
+                )
+            );
+
+            faceProcAddressing.set
+            (
+                procI,
+                new labelIOList
+                (
+                    IOobject
+                    (
+                        "faceRegionAddressing",
+                        procMeshes[procI].facesInstance(),
+                        procMeshes[procI].meshSubDir,
+                        procMeshes[procI],
+                        IOobject::MUST_READ,
+                        IOobject::NO_WRITE
+                    )
+                )
+            );
+
+            boundaryProcAddressing.set
+            (
+                procI,
+                new labelIOList
+                (
+                    IOobject
+                    (
+                        "boundaryRegionAddressing",
+                        procMeshes[procI].facesInstance(),
+                        procMeshes[procI].meshSubDir,
+                        procMeshes[procI],
+                        IOobject::MUST_READ,
+                        IOobject::NO_WRITE
+                    )
+                )
+            );
+        }
     }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/initContinuityErrs.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/initContinuityErrs.H
index 1a7f5a3262c79506f0459fd1bedbbbe1bd086e72..8cc47b1ca01a5fb88563b8ec14c205a097a68e19 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/initContinuityErrs.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/initContinuityErrs.H
@@ -1 +1 @@
-List<scalar> cumulativeContErr(fluidRegions.size(), 0.0);
+List<scalar> cumulativeContErr(fluidToProc.size(), 0.0);
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/mapFluid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/mapFluid.H
new file mode 100644
index 0000000000000000000000000000000000000000..75d83afd09ffe23a1f85014b71f747b4e35920b2
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/mapFluid.H
@@ -0,0 +1,3 @@
+h = allToProcMappers[procI].decomposeField(allh(), true);
+h.oldTime().timeIndex() = allh().oldTime().timeIndex();
+h.correctBoundaryConditions();
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/rmapFluid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/rmapFluid.H
new file mode 100644
index 0000000000000000000000000000000000000000..576d883f0129b711b2d50ba317aa42d8dd2b0536
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/rmapFluid.H
@@ -0,0 +1,56 @@
+// Note:Map rho and rho.oldTime() since fluid rho assigned to at
+// end of iteration.
+rmap
+(
+    allRho(),
+    rho,
+    faceProcAddressing[procI],
+    cellProcAddressing[procI],
+    boundaryProcAddressing[procI]
+);
+rmap
+(
+    allRho().oldTime(),
+    rho.oldTime(),
+    faceProcAddressing[procI],
+    cellProcAddressing[procI],
+    boundaryProcAddressing[procI]
+);
+
+// Necessary? Probably only for boundary values since bcs on
+// h are not the same as those on allh
+
+rmap
+(
+    allh(),
+    h,
+    faceProcAddressing[procI],
+    cellProcAddressing[procI],
+    boundaryProcAddressing[procI]
+);
+
+procAlpha.set(procI, fvc::interpolate(turb.alphaEff()));
+
+rmap
+(
+    allVolAlpha(),
+    turb.alphaEff()(),
+    faceProcAddressing[procI],
+    cellProcAddressing[procI],
+    boundaryProcAddressing[procI]
+);
+
+rmap
+(
+    allSource(),
+    (DpDt + rad.Sh(thermo))(),
+    faceProcAddressing[procI],
+    cellProcAddressing[procI],
+    boundaryProcAddressing[procI]
+);
+
+procPhi.set
+(
+    procI,
+    new surfaceScalarField(phi)
+);
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
index 81c6d25bb0ccc68597a0b154e0b7c3baa98552ba..5bc0fffa4048d57a290dfd5296fbc0beea62823a 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
@@ -1,4 +1,4 @@
-    fvMesh& mesh = fluidRegions[i];
+    fvMesh& mesh = procMeshes[fluidToProc[i]];
 
     basicRhoThermo& thermo = thermoFluid[i];
     volScalarField& rho = rhoFluid[i];
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solvePressureVelocityFluid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solvePressureVelocityFluid.H
new file mode 100644
index 0000000000000000000000000000000000000000..d04a301112c86483cdd81493fd5985f3caca2f8b
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solvePressureVelocityFluid.H
@@ -0,0 +1,11 @@
+#include "UEqn.H"
+
+// --- PISO loop
+for (int corr=0; corr<nCorr; corr++)
+{
+    #include "pEqn.H"
+}
+
+turb.correct();
+
+rho = thermo.rho();
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/readPIMPLEControls.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/readPIMPLEControls.H
index 7405dee8d5ae5a8c9eb1caa8f867526a6eac4823..6d34d4388fd5f41d9134bb9c9b3de23bd07bca69 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/readPIMPLEControls.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/readPIMPLEControls.H
@@ -6,3 +6,6 @@
 
     const int nOuterCorr =
         pimple.lookupOrDefault<int>("nOuterCorrectors", 1);
+
+    const Switch temperatureCoupled =
+        pimple.lookupOrDefault<Switch>("temperatureCoupled", false);
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/rmap.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/rmap.H
new file mode 100644
index 0000000000000000000000000000000000000000..d9c8e350600585212ec7c5bca78bf0bc4a211e9a
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/rmap.H
@@ -0,0 +1,96 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 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
+    rmap
+
+Description
+    map field on subset of mesh onto overall field. Rewrites boundary
+    conditions to be mixed.
+
+    The source fields can have different patch types for the same destination
+    patch. To work around this it attempts to convert all patch fields into
+    mixed type since this can accomodate anything from fixedValue to
+    fixedGradient.
+
+SourceFiles
+    rmap.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef rmap_H
+#define rmap_H
+
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//- Map patchField
+template<class Type>
+static void rmap
+(
+    fvPatchField<Type>& destBC,
+    const labelList& reverseAddressing,
+    const fvPatchField<Type>& sourceBC
+);
+
+//- Map volField
+template<class Type>
+static void rmap
+(
+    GeometricField<Type, fvPatchField, volMesh>& dest,
+    const GeometricField<Type, fvPatchField, volMesh>& source,
+    const labelList& faceProcAddressing,
+    const labelList& cellProcAddressing,
+    const labelList& boundaryProcAddressing
+);
+
+//- Map fvMatrix
+template<class Type>
+static void rmap
+(
+    fvMatrix<Type>& dest,
+    const fvMatrix<Type>& source,
+    const labelList& faceProcAddressing,
+    const labelList& cellProcAddressing,
+    const labelList& boundaryProcAddressing
+);
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "rmapTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/rmapTemplates.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/rmapTemplates.C
new file mode 100644
index 0000000000000000000000000000000000000000..6dbf50a81c201d52570cdea50287fa609938bc23
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/rmapTemplates.C
@@ -0,0 +1,309 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 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 "rmap.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::rmap
+(
+    fvPatchField<Type>& destBC,
+    const labelList& reverseAddressing,
+    const fvPatchField<Type>& sourceBC
+)
+{
+    // Assign value
+    destBC.Field<Type>::rmap(sourceBC, reverseAddressing);
+
+    // Assign other properties
+    if (isA<mixedFvPatchField<Type> >(destBC))
+    {
+        mixedFvPatchField<Type>& mp =
+            refCast<mixedFvPatchField<Type> >(destBC);
+
+        if (isA<mixedFvPatchField<Type> >(sourceBC))
+        {
+            const mixedFvPatchField<Type>& Tp =
+                refCast<const mixedFvPatchField<Type> >(sourceBC);
+
+            mp.refValue().rmap(Tp.refValue(), reverseAddressing);
+            mp.refGrad().rmap(Tp.refGrad(), reverseAddressing);
+            mp.valueFraction().rmap(Tp.valueFraction(), reverseAddressing);
+        }
+        else if (isA<fixedGradientFvPatchField<Type> >(sourceBC))
+        {
+            const fixedGradientFvPatchField<Type>& Tp =
+                refCast<const fixedGradientFvPatchField<Type> >
+                (
+                    sourceBC
+                );
+            // Make pure fixedGradient
+            mp.refValue().rmap(Tp, reverseAddressing); // unused
+            mp.refGrad().rmap(Tp.gradient(), reverseAddressing);
+            mp.valueFraction().rmap
+            (
+                Field<Type>(reverseAddressing.size(), 0.0),
+                reverseAddressing
+            );
+        }
+        else if (isA<zeroGradientFvPatchField<Type> >(sourceBC))
+        {
+            // Make pure fixedGradient with gradient = 0
+            mp.refValue().rmap(sourceBC, reverseAddressing); // unused
+            mp.refGrad().rmap
+            (
+                Field<Type>(reverseAddressing.size(), 0.0),
+                reverseAddressing
+            );
+            mp.valueFraction().rmap
+            (
+                Field<Type>(reverseAddressing.size(), 0.0),
+                reverseAddressing
+            );
+        }
+        else if (isA<fixedValueFvPatchField<Type> >(sourceBC))
+        {
+            // Make pure fixedValue
+            mp.refValue().rmap(sourceBC, reverseAddressing);
+            mp.refGrad().rmap
+            (
+                Field<Type>(reverseAddressing.size(), 0.0),
+                reverseAddressing
+            ); // unused
+            mp.valueFraction().rmap
+            (
+                Field<Type>(reverseAddressing.size(), 1.0),
+                reverseAddressing
+            );
+        }
+        else if (isA<calculatedFvPatchField<Type> >(sourceBC))
+        {
+            // Make pure fixedValue
+            mp.refValue().rmap(sourceBC, reverseAddressing);
+            mp.refGrad().rmap
+            (
+                Field<Type>(reverseAddressing.size(), 0.0),
+                reverseAddressing
+            ); // unused
+            mp.valueFraction().rmap
+            (
+                Field<Type>(reverseAddressing.size(), 1.0),
+                reverseAddressing
+            );
+        }
+        else
+        {
+            FatalErrorIn("rmap(..)")
+                << "Don't know how to map source bc "
+                << sourceBC.type()
+                << " into a mixed boundary condition at "
+                << destBC.patch().name()
+                << exit(FatalError);
+        }
+    }
+    else if (isA<fixedGradientFvPatchField<Type> >(destBC))
+    {
+        fixedGradientFvPatchField<Type>& mp =
+            refCast<fixedGradientFvPatchField<Type> >(destBC);
+
+        if (isA<fixedGradientFvPatchField<Type> >(sourceBC))
+        {
+            const fixedGradientFvPatchField<Type>& Tp =
+                refCast<const fixedGradientFvPatchField<Type> >
+                (
+                    sourceBC
+                );
+            mp.gradient().rmap(Tp.gradient(), reverseAddressing);
+        }
+        else if (isA<mixedFvPatchField<Type> >(sourceBC))
+        {
+            const mixedFvPatchField<Type>& Tp =
+                refCast<const mixedFvPatchField<Type> >(sourceBC);
+            mp.gradient().rmap(Tp.snGrad(), reverseAddressing);
+        }
+        else if (isA<zeroGradientFvPatchField<Type> >(sourceBC))
+        {
+            mp.gradient().rmap
+            (
+                Field<Type>(reverseAddressing.size(), 0.0),
+                reverseAddressing
+            );
+        }
+        else
+        {
+            FatalErrorIn("rmap(..)")
+                << "Don't know how to map source bc "
+                << sourceBC.type()
+                << " into a fixedGradient boundary condition at "
+                << destBC.patch().name()
+                << exit(FatalError);
+        }
+    }
+}
+
+
+template<class Type>
+void Foam::rmap
+(
+    GeometricField<Type, fvPatchField, volMesh>& dest,
+    const GeometricField<Type, fvPatchField, volMesh>& source,
+    const labelList& faceProcAddressing,
+    const labelList& cellProcAddressing,
+    const labelList& boundaryProcAddressing
+)
+{
+    if (dest.dimensions() != source.dimensions())
+    {
+        FatalErrorIn("rmap(..)")
+            << "Different dimensions for = for fields " << dest.name()
+            << " and " << source.name() << endl
+            << "     dimensions : " << dest.dimensions()
+            << " = " << source.dimensions() << endl
+            << exit(FatalError);
+    }
+
+    // Copy internal field
+    dest.internalField().rmap(source.internalField(), cellProcAddressing);
+
+    // Copy boundary properties as mixed
+    forAll(source.boundaryField(), patchI)
+    {
+        label curBPatch = boundaryProcAddressing[patchI];
+
+        if (curBPatch == -1)
+        {
+            // Unknown patch. Do not change any values.
+        }
+        else
+        {
+            // Get addressing slice for this patch
+            const labelList::subList cp =
+                source.mesh().boundary()[patchI].patchSlice
+                (
+                    faceProcAddressing
+                );
+
+            const label curPatchStart =
+                dest.mesh().boundaryMesh()[curBPatch].start();
+
+            labelList reverseAddressing(cp.size());
+
+            forAll(cp, faceI)
+            {
+                // Subtract one to take into account offsets for
+                // face direction.
+                if (cp[faceI] <= 0)
+                {
+                    FatalErrorIn("rmap(..)")
+                        << "Problem:"
+                        << " patch:" << source.mesh().boundary()[patchI].name()
+                        << " field:" << source.name()
+                        << " local face:" << faceI
+                        << " mapped to:" << cp[faceI] << exit(FatalError);
+                }
+
+                reverseAddressing[faceI] = cp[faceI] - 1 - curPatchStart;
+            }
+
+            // Map curBPatch from source patch. Is like rmap but also
+            // copies non-value properties from alike patchFields.
+            rmap
+            (
+                dest.boundaryField()[curBPatch],
+                reverseAddressing,
+                source.boundaryField()[patchI]
+            );
+        }
+    }
+
+    // Copy timeIndex
+    dest.timeIndex() = source.timeIndex();
+}
+
+
+template<class Type>
+void Foam::rmap
+(
+    fvMatrix<Type>& dest,
+    const fvMatrix<Type>& source,
+    const labelList& faceProcAddressing,
+    const labelList& cellProcAddressing,
+    const labelList& boundaryProcAddressing
+)
+{
+    dest.source().rmap(source.source(), cellProcAddressing);
+
+    FieldField<Field, Type>& sourceInternal =
+        const_cast<fvMatrix<Type>&>(source).internalCoeffs();
+    FieldField<Field, Type>& sourceBoundary =
+        const_cast<fvMatrix<Type>&>(source).boundaryCoeffs();
+
+    forAll(sourceInternal, patchI)
+    {
+        label curBPatch = boundaryProcAddressing[patchI];
+
+        if (curBPatch == -1)
+        {
+            // Unknown patch. Do not change any values.
+        }
+        else
+        {
+            // Get addressing slice for this patch
+            const fvMesh& sourceMesh = source.psi().mesh();
+
+            const labelList::subList cp =
+                sourceMesh.boundary()[patchI].patchSlice
+                (
+                    faceProcAddressing
+                );
+
+            const label curPatchStart =
+                dest.psi().mesh().boundaryMesh()[curBPatch].start();
+
+            labelList reverseAddressing(cp.size());
+
+            forAll(cp, faceI)
+            {
+                // Subtract one to take into account offsets for
+                // face direction.
+                reverseAddressing[faceI] = cp[faceI] - 1 - curPatchStart;
+            }
+            dest.internalCoeffs()[curBPatch].rmap
+            (
+                sourceInternal[patchI],
+                reverseAddressing
+            );
+            dest.boundaryCoeffs()[curBPatch].rmap
+            (
+                sourceBoundary[patchI],
+                reverseAddressing
+            );
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidFields.H
index 837305659e6c088fa7048bb0e174a6da57fc418c..ab58752d80a0573d4fb8f51f2609375feaa4623c 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidFields.H
@@ -1,12 +1,36 @@
     // Initialise solid field pointer lists
-    PtrList<basicSolidThermo> thermos(solidRegions.size());
+    PtrList<basicSolidThermo> thermoSolid(rp.solidRegionNames().size());
+    PtrList<volScalarField> hSolid(rp.solidRegionNames().size());
 
     // Populate solid field pointer lists
-    forAll(solidRegions, i)
+    forAll(rp.solidRegionNames(), i)
     {
         Info<< "*** Reading solid mesh thermophysical properties for region "
-            << solidRegions[i].name() << nl << endl;
+            << rp.solidRegionNames()[i] << nl << endl;
 
-        Info<< "    Adding to thermos\n" << endl;
-        thermos.set(i, basicSolidThermo::New(solidRegions[i]));
+        label procI = solidToProc[i];
+
+        Info<< "    Adding to thermoSolid\n" << endl;
+        thermoSolid.set(i, basicSolidThermo::New(procMeshes[procI]));
+
+        if (temperatureCoupled)
+        {
+            Info<< "    Adding to hSolid\n" << endl;
+            hSolid.set
+            (
+                i,
+                new volScalarField
+                (
+                    IOobject
+                    (
+                        "h",
+                        runTime.timeName(),
+                        procMeshes[procI],
+                        IOobject::MUST_READ,
+                        IOobject::AUTO_WRITE
+                    ),
+                    procMeshes[procI]
+                )
+            );
+        }
     }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidMeshes.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidMeshes.H
index eb50be23808e6ffdf31c17fca398bb5366f24c7e..10b5f4160cee55cc88d1e4fa77f63f1fe229f3e7 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidMeshes.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidMeshes.H
@@ -1,13 +1,13 @@
-    PtrList<fvMesh> solidRegions(rp.solidRegionNames().size());
-
     forAll(rp.solidRegionNames(), i)
     {
         Info<< "Create solid mesh for region " << rp.solidRegionNames()[i]
             << " for time = " << runTime.timeName() << nl << endl;
 
-        solidRegions.set
+        label procI = solidToProc[i];
+
+        procMeshes.set
         (
-            i,
+            procI,
             new fvMesh
             (
                 IOobject
@@ -20,8 +20,57 @@
             )
         );
 
-        // Force calculation of geometric properties to prevent it being done
-        // later in e.g. some boundary evaluation
-        //(void)solidRegions[i].weights();
-        //(void)solidRegions[i].deltaCoeffs();
+        if (temperatureCoupled)
+        {
+            cellProcAddressing.set
+            (
+                procI,
+                new labelIOList
+                (
+                    IOobject
+                    (
+                        "cellRegionAddressing",
+                        procMeshes[procI].facesInstance(),
+                        procMeshes[procI].meshSubDir,
+                        procMeshes[procI],
+                        IOobject::MUST_READ,
+                        IOobject::NO_WRITE
+                    )
+                )
+            );
+
+            faceProcAddressing.set
+            (
+                procI,
+                new labelIOList
+                (
+                    IOobject
+                    (
+                        "faceRegionAddressing",
+                        procMeshes[procI].facesInstance(),
+                        procMeshes[procI].meshSubDir,
+                        procMeshes[procI],
+                        IOobject::MUST_READ,
+                        IOobject::NO_WRITE
+                    )
+                )
+            );
+
+            boundaryProcAddressing.set
+            (
+                procI,
+                new labelIOList
+                (
+                    IOobject
+                    (
+                        "boundaryRegionAddressing",
+                        procMeshes[procI].facesInstance(),
+                        procMeshes[procI].meshSubDir,
+                        procMeshes[procI],
+                        IOobject::MUST_READ,
+                        IOobject::NO_WRITE
+                    )
+                )
+            );
+        }
     }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/mapSolid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/mapSolid.H
new file mode 100644
index 0000000000000000000000000000000000000000..ca6e506347d5d76fba7b19299e545fabea583ca2
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/mapSolid.H
@@ -0,0 +1,15 @@
+{
+    volScalarField& h = hSolid[i];
+
+    h = allToProcMappers[procI].decomposeField(allh(), true);
+    h.oldTime().timeIndex() = allh().oldTime().timeIndex();
+
+    T += (h-h.oldTime())/cp;
+    // Correct T boundary conditions and update h boundary
+    // conditions accordingly.
+    volScalarField::GeometricBoundaryField Told = T.boundaryField();
+    T.correctBoundaryConditions();
+    h.boundaryField() +=
+        cp.boundaryField()
+      * (T.boundaryField()-Told);
+}
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/rmapSolid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/rmapSolid.H
new file mode 100644
index 0000000000000000000000000000000000000000..b36474a531cb66108eba65be7e6c8dbb61704767
--- /dev/null
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/rmapSolid.H
@@ -0,0 +1,90 @@
+{
+    volScalarField& h = hSolid[i];
+
+    procPhi.setSize(nAllRegions);
+    procAlpha.setSize(nAllRegions);
+
+    rmap
+    (
+        allRho(),
+        rho,
+        faceProcAddressing[procI],
+        cellProcAddressing[procI],
+        boundaryProcAddressing[procI]
+    );
+
+    // Necessary? Probably only for boundary values since bcs on
+    // h are not the same as those on allh
+    rmap
+    (
+        allh(),
+        h,
+        faceProcAddressing[procI],
+        cellProcAddressing[procI],
+        boundaryProcAddressing[procI]
+    );
+
+
+    tmp<volScalarField> Kcp(K/cp);
+
+    rmap
+    (
+        allVolAlpha(),
+        Kcp(),
+        faceProcAddressing[procI],
+        cellProcAddressing[procI],
+        boundaryProcAddressing[procI]
+    );
+
+    procAlpha.set(procI, fvc::interpolate(Kcp));
+
+    // allSource is initialised to zero already
+    //rmap
+    //(
+    //    allSource(),
+    //    volScalarField
+    //    (
+    //        IOobject
+    //        (
+    //            "procSource",
+    //            runTime.timeName(),
+    //            mesh,
+    //            IOobject::NO_READ,
+    //            IOobject::AUTO_WRITE
+    //        ),
+    //        mesh,
+    //        dimensionedScalar
+    //        (
+    //            "procSource",
+    //            allh().dimensions()*dimDensity/dimTime,
+    //            0.0
+    //        )
+    //    ),
+    //    faceProcAddressing[procI],
+    //    cellProcAddressing[procI],
+    //    boundaryProcAddressing[procI]
+    //);
+
+    procPhi.set
+    (
+        procI,
+        new surfaceScalarField
+        (
+            IOobject
+            (
+                "phi",
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            mesh,
+            dimensionedScalar
+            (
+                "phi",
+                dimDensity*dimVelocity*dimArea,
+                0.0
+            )
+        )
+    );
+}
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H
index a843ed8bd75f562bbba0dfb157356b36aaaaf5fd..a1dbe88c0abeb1a08bd53f4270cee3d830d385ff 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H
@@ -1,5 +1,5 @@
-    fvMesh& mesh = solidRegions[i];
-    basicSolidThermo& thermo = thermos[i];
+    fvMesh& mesh = procMeshes[solidToProc[i]];
+    basicSolidThermo& thermo = thermoSolid[i];
 
     tmp<volScalarField> trho = thermo.rho();
     const volScalarField& rho = trho();
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H
index 77dc6f04bf44ba748359c2f8bf3267307416b650..5564d3bd3a8ada7dd6d3f39e755ae8a51e903015 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H
@@ -1,6 +1,6 @@
     scalar DiNum = -GREAT;
 
-    forAll(solidRegions, i)
+    forAll(solidToProc, i)
     {
 #       include "setRegionSolidFields.H"
 
@@ -8,7 +8,7 @@
         (
             solidRegionDiffNo
             (
-                solidRegions[i],
+                procMeshes[solidToProc[i]],
                 runTime,
                 rho*cp,
                 K
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H
index d8aa03283b413d632634933e396f09ac4ab64e9b..96cd32ddb37a3e7ef2aebd1f14602692177c4b43 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H
@@ -1,8 +1,3 @@
-if (finalIter)
-{
-    mesh.data::add("finalIteration", true);
-}
-
 {
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
@@ -17,10 +12,3 @@ if (finalIter)
 
     Info<< "Min/max T:" << min(T) << ' ' << max(T) << endl;
 }
-
-thermo.correct();
-
-if (finalIter)
-{
-    mesh.data::remove("finalIteration");
-}
diff --git a/applications/utilities/parallelProcessing/decomposePar/Make/files b/applications/utilities/parallelProcessing/decomposePar/Make/files
index 125f6b81e70155299ceff3857b31f7c0326d7090..cad5ac1ce3fba213b585ad2d7e3dff96b3bb9054 100644
--- a/applications/utilities/parallelProcessing/decomposePar/Make/files
+++ b/applications/utilities/parallelProcessing/decomposePar/Make/files
@@ -3,7 +3,6 @@ domainDecomposition.C
 domainDecompositionMesh.C
 domainDecompositionDistribute.C
 dimFieldDecomposer.C
-fvFieldDecomposer.C
 pointFieldDecomposer.C
 lagrangianFieldDecomposer.C
 
diff --git a/applications/utilities/parallelProcessing/decomposePar/Make/options b/applications/utilities/parallelProcessing/decomposePar/Make/options
index 7a9f1df3f592ddf510c6a99c6754d3dc5cd20ecb..7194eb1bce4743f0b8a0e7961ef91c49d903a339 100644
--- a/applications/utilities/parallelProcessing/decomposePar/Make/options
+++ b/applications/utilities/parallelProcessing/decomposePar/Make/options
@@ -1,4 +1,5 @@
 EXE_INC = \
+    -I$(LIB_SRC)/parallel/decompose/decompose/lnInclude \
     -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
@@ -6,6 +7,7 @@ EXE_INC = \
 
 EXE_LIBS = \
     -lfiniteVolume \
+    -ldecompose \
     -lgenericPatchFields \
     -ldecompositionMethods -L$(FOAM_LIBBIN)/dummy -lmetisDecomp -lscotchDecomp \
     -llagrangian \
diff --git a/src/parallel/decompose/Allwmake b/src/parallel/decompose/Allwmake
index b620b27791e388a86f88b3b59474cc60769a7cf4..9d0acbd618a5b0fa9b7de804494e790f3d59f195 100755
--- a/src/parallel/decompose/Allwmake
+++ b/src/parallel/decompose/Allwmake
@@ -53,4 +53,6 @@ fi
 
 wmake $makeType decompositionMethods
 
+wmake $makeType decompose
+
 # ----------------------------------------------------------------- end-of-file
diff --git a/src/parallel/decompose/decompose/Make/files b/src/parallel/decompose/decompose/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..2d0f959eaed2fa913a61f4cdbaf873f4d465ece5
--- /dev/null
+++ b/src/parallel/decompose/decompose/Make/files
@@ -0,0 +1,3 @@
+fvFieldDecomposer.C
+
+LIB = $(FOAM_LIBBIN)/libdecompose
diff --git a/src/parallel/decompose/decompose/Make/options b/src/parallel/decompose/decompose/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..7a728f9dd7cf75800c9b44943a2964c781376576
--- /dev/null
+++ b/src/parallel/decompose/decompose/Make/options
@@ -0,0 +1,9 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude
+
+LIB_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -llagrangian
diff --git a/applications/utilities/parallelProcessing/decomposePar/fvFieldDecomposer.C b/src/parallel/decompose/decompose/fvFieldDecomposer.C
similarity index 100%
rename from applications/utilities/parallelProcessing/decomposePar/fvFieldDecomposer.C
rename to src/parallel/decompose/decompose/fvFieldDecomposer.C
diff --git a/applications/utilities/parallelProcessing/decomposePar/fvFieldDecomposer.H b/src/parallel/decompose/decompose/fvFieldDecomposer.H
similarity index 99%
rename from applications/utilities/parallelProcessing/decomposePar/fvFieldDecomposer.H
rename to src/parallel/decompose/decompose/fvFieldDecomposer.H
index ee87db59f41bf3e00d3ffbe17f5dac53656a4cdf..899b45b4b7550369f031279310f0456174cd9edc 100644
--- a/applications/utilities/parallelProcessing/decomposePar/fvFieldDecomposer.H
+++ b/src/parallel/decompose/decompose/fvFieldDecomposer.H
@@ -241,7 +241,8 @@ public:
         tmp<GeometricField<Type, fvPatchField, volMesh> >
         decomposeField
         (
-            const GeometricField<Type, fvPatchField, volMesh>& field
+            const GeometricField<Type, fvPatchField, volMesh>& field,
+            const bool allowUnknownPatchFields = false
         ) const;
 
         //- Decompose surface field
diff --git a/applications/utilities/parallelProcessing/decomposePar/fvFieldDecomposerDecomposeFields.C b/src/parallel/decompose/decompose/fvFieldDecomposerDecomposeFields.C
similarity index 94%
rename from applications/utilities/parallelProcessing/decomposePar/fvFieldDecomposerDecomposeFields.C
rename to src/parallel/decompose/decompose/fvFieldDecomposerDecomposeFields.C
index 0bce0023cc8b3c23e94e2213a79464351921ed08..ca9d0e31b5eaa8f2c361864c93e32ee5e6fd452c 100644
--- a/applications/utilities/parallelProcessing/decomposePar/fvFieldDecomposerDecomposeFields.C
+++ b/src/parallel/decompose/decompose/fvFieldDecomposerDecomposeFields.C
@@ -28,6 +28,7 @@ License
 #include "processorFvsPatchField.H"
 #include "processorCyclicFvPatchField.H"
 #include "processorCyclicFvsPatchField.H"
+#include "emptyFvPatchFields.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
@@ -35,7 +36,8 @@ template<class Type>
 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
 Foam::fvFieldDecomposer::decomposeField
 (
-    const GeometricField<Type, fvPatchField, volMesh>& field
+    const GeometricField<Type, fvPatchField, volMesh>& field,
+    const bool allowUnknownPatchFields
 ) const
 {
     // Create and map the internal field values
@@ -94,6 +96,18 @@ Foam::fvFieldDecomposer::decomposeField
                 )
             );
         }
+        else if (allowUnknownPatchFields)
+        {
+            patchFields.set
+            (
+                patchi,
+                new emptyFvPatchField<Type>
+                (
+                    procMesh_.boundary()[patchi],
+                    DimensionedField<Type, volMesh>::null()
+                )
+            );
+        }
         else
         {
             FatalErrorIn("fvFieldDecomposer::decomposeField()")
diff --git a/src/parallel/reconstruct/reconstruct/fvFieldReconstructor.C b/src/parallel/reconstruct/reconstruct/fvFieldReconstructor.C
index 93fbb8bb995831528b07f879b15d8d77efe122e0..3c85d65c6eb68e97191d0a3e7210e2d42c9ae75f 100644
--- a/src/parallel/reconstruct/reconstruct/fvFieldReconstructor.C
+++ b/src/parallel/reconstruct/reconstruct/fvFieldReconstructor.C
@@ -42,7 +42,40 @@ Foam::fvFieldReconstructor::fvFieldReconstructor
     cellProcAddressing_(cellProcAddressing),
     boundaryProcAddressing_(boundaryProcAddressing),
     nReconstructed_(0)
-{}
+{
+    forAll(procMeshes_, procI)
+    {
+        const fvMesh& procMesh = procMeshes_[procI];
+        if
+        (
+            faceProcAddressing[procI].size() != procMesh.nFaces()
+         || cellProcAddressing[procI].size() != procMesh.nCells()
+         || boundaryProcAddressing[procI].size() != procMesh.boundary().size()
+        )
+        {
+            FatalErrorIn
+            (
+                "fvFieldReconstructor::fvFieldReconstructor\n"
+                "(\n"
+                "   fvMesh&,\n"
+                "   const PtrList<fvMesh>&,\n"
+                "   const PtrList<labelIOList>&,\n"
+                "   const PtrList<labelIOList>&,\n"
+                "   const PtrList<labelIOList>&\n"
+                ")"
+            )   << "Size of maps does not correspond to size of mesh"
+                << " for processor " << procI << endl
+                << "faceProcAddressing : " << faceProcAddressing[procI].size()
+                << " nFaces : " << procMesh.nFaces() << endl
+                << "cellProcAddressing : " << cellProcAddressing[procI].size()
+                << " nCell : " << procMesh.nCells() << endl
+                << "boundaryProcAddressing : "
+                << boundaryProcAddressing[procI].size()
+                << " nFaces : " << procMesh.boundary().size()
+                << exit(FatalError);
+        }
+    }
+}
 
 
 // ************************************************************************* //
diff --git a/src/parallel/reconstruct/reconstruct/fvFieldReconstructor.H b/src/parallel/reconstruct/reconstruct/fvFieldReconstructor.H
index 2aa3011195c8d23d20e8d63ea51ebf0c161af818..43469d24b4f0b6899b08833516ee84e44b3bd709 100644
--- a/src/parallel/reconstruct/reconstruct/fvFieldReconstructor.H
+++ b/src/parallel/reconstruct/reconstruct/fvFieldReconstructor.H
@@ -86,6 +86,7 @@ class fvFieldReconstructor
 
 public:
 
+        //- Mapper for sizing only - does not do any actual mapping.
         class fvPatchFieldReconstructor
         :
             public fvPatchFieldMapper
@@ -146,19 +147,48 @@ public:
         //- Reconstruct volume internal field
         template<class Type>
         tmp<DimensionedField<Type, volMesh> >
-        reconstructFvVolumeInternalField(const IOobject& fieldIoObject);
+        reconstructFvVolumeInternalField
+        (
+            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<GeometricField<Type, fvPatchField, volMesh> >
-        reconstructFvVolumeField(const IOobject& fieldIoObject);
+        reconstructFvVolumeField
+        (
+            const IOobject& fieldIoObject,
+            const PtrList<GeometricField<Type, fvPatchField, volMesh> >&
+        ) const;
+
+        //- Read and reconstruct volume field
+        template<class Type>
+        tmp<GeometricField<Type, fvPatchField, volMesh> >
+        reconstructFvVolumeField(const IOobject& fieldIoObject) const;
+
 
         //- Reconstruct surface field
         template<class Type>
         tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
-        reconstructFvSurfaceField(const IOobject& fieldIoObject);
+        reconstructFvSurfaceField
+        (
+            const IOobject& fieldIoObject,
+            const PtrList<GeometricField<Type, fvsPatchField, surfaceMesh> >&
+        ) const;
+
+        //- Read and reconstruct surface field
+        template<class Type>
+        tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
+        reconstructFvSurfaceField(const IOobject& fieldIoObject) const;
 
-        //- Reconstruct and write all/selected volume internal fields
+        //- Read, reconstruct and write all/selected volume internal fields
         template<class Type>
         void reconstructFvVolumeInternalFields
         (
@@ -166,7 +196,7 @@ public:
             const HashSet<word>& selectedFields
         );
 
-        //- Reconstruct and write all/selected volume fields
+        //- Read, reconstruct and write all/selected volume fields
         template<class Type>
         void reconstructFvVolumeFields
         (
@@ -174,7 +204,7 @@ public:
             const HashSet<word>& selectedFields
         );
 
-        //- Reconstruct and write all/selected surface fields
+        //- Read, reconstruct and write all/selected surface fields
         template<class Type>
         void reconstructFvSurfaceFields
         (
diff --git a/src/parallel/reconstruct/reconstruct/fvFieldReconstructorReconstructFields.C b/src/parallel/reconstruct/reconstruct/fvFieldReconstructorReconstructFields.C
index 9a0dfb4b3cf3d519b45ea4dadb38b2fb5702189d..33105892ce685436b15d5364686cb0746558fac5 100644
--- a/src/parallel/reconstruct/reconstruct/fvFieldReconstructorReconstructFields.C
+++ b/src/parallel/reconstruct/reconstruct/fvFieldReconstructorReconstructFields.C
@@ -37,36 +37,10 @@ template<class Type>
 Foam::tmp<Foam::DimensionedField<Type, Foam::volMesh> >
 Foam::fvFieldReconstructor::reconstructFvVolumeInternalField
 (
-    const IOobject& fieldIoObject
-)
+    const IOobject& fieldIoObject,
+    const PtrList<DimensionedField<Type, volMesh> >& procFields
+) const
 {
-    // Read the field for all the processors
-    PtrList<DimensionedField<Type, volMesh> > procFields
-    (
-        procMeshes_.size()
-    );
-
-    forAll(procMeshes_, procI)
-    {
-        procFields.set
-        (
-            procI,
-            new DimensionedField<Type, volMesh>
-            (
-                IOobject
-                (
-                    fieldIoObject.name(),
-                    procMeshes_[procI].time().timeName(),
-                    procMeshes_[procI],
-                    IOobject::MUST_READ,
-                    IOobject::NO_WRITE
-                ),
-                procMeshes_[procI]
-            )
-        );
-    }
-
-
     // Create the internalField
     Field<Type> internalField(mesh_.nCells());
 
@@ -86,14 +60,7 @@ Foam::fvFieldReconstructor::reconstructFvVolumeInternalField
     (
         new DimensionedField<Type, volMesh>
         (
-            IOobject
-            (
-                fieldIoObject.name(),
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
+            fieldIoObject,
             mesh_,
             procFields[0].dimensions(),
             internalField
@@ -103,14 +70,14 @@ Foam::fvFieldReconstructor::reconstructFvVolumeInternalField
 
 
 template<class Type>
-Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
-Foam::fvFieldReconstructor::reconstructFvVolumeField
+Foam::tmp<Foam::DimensionedField<Type, Foam::volMesh> >
+Foam::fvFieldReconstructor::reconstructFvVolumeInternalField
 (
     const IOobject& fieldIoObject
-)
+) const
 {
     // Read the field for all the processors
-    PtrList<GeometricField<Type, fvPatchField, volMesh> > procFields
+    PtrList<DimensionedField<Type, volMesh> > procFields
     (
         procMeshes_.size()
     );
@@ -120,7 +87,7 @@ Foam::fvFieldReconstructor::reconstructFvVolumeField
         procFields.set
         (
             procI,
-            new GeometricField<Type, fvPatchField, volMesh>
+            new DimensionedField<Type, volMesh>
             (
                 IOobject
                 (
@@ -136,14 +103,36 @@ Foam::fvFieldReconstructor::reconstructFvVolumeField
     }
 
 
+    return reconstructFvVolumeInternalField
+    (
+        IOobject
+        (
+            fieldIoObject.name(),
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        procFields
+    );
+}
+
+
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
+Foam::fvFieldReconstructor::reconstructFvVolumeField
+(
+    const IOobject& fieldIoObject,
+    const PtrList<GeometricField<Type, fvPatchField, volMesh> >& procFields
+) const
+{
     // Create the internalField
     Field<Type> internalField(mesh_.nCells());
 
     // Create the patch fields
     PtrList<fvPatchField<Type> > patchFields(mesh_.boundary().size());
 
-
-    forAll(procMeshes_, procI)
+    forAll(procFields, procI)
     {
         const GeometricField<Type, fvPatchField, volMesh>& procField =
             procFields[procI];
@@ -163,7 +152,7 @@ Foam::fvFieldReconstructor::reconstructFvVolumeField
 
             // Get addressing slice for this patch
             const labelList::subList cp =
-                procMeshes_[procI].boundary()[patchI].patchSlice
+                procField.mesh().boundary()[patchI].patchSlice
                 (
                     faceProcAddressing_[procI]
                 );
@@ -198,11 +187,32 @@ Foam::fvFieldReconstructor::reconstructFvVolumeField
 
                 forAll(cp, faceI)
                 {
+                    // Check
+                    if (cp[faceI] <= 0)
+                    {
+                        FatalErrorIn
+                        (
+                            "fvFieldReconstructor::reconstructFvVolumeField\n"
+                            "(\n"
+                            "    const IOobject&,\n"
+                            "    const PtrList<GeometricField<Type,"
+                            " fvPatchField, volMesh> >&\n"
+                            ") const\n"
+                        )   << "Processor " << procI
+                            << " patch "
+                            << procField.mesh().boundary()[patchI].name()
+                            << " face " << faceI
+                            << " originates from reversed face since "
+                            << cp[faceI]
+                            << exit(FatalError);
+                    }
+
                     // Subtract one to take into account offsets for
                     // face direction.
                     reverseAddressing[faceI] = cp[faceI] - 1 - curPatchStart;
                 }
 
+
                 patchFields[curBPatch].rmap
                 (
                     procField.boundaryField()[patchI],
@@ -283,14 +293,7 @@ Foam::fvFieldReconstructor::reconstructFvVolumeField
     (
         new GeometricField<Type, fvPatchField, volMesh>
         (
-            IOobject
-            (
-                fieldIoObject.name(),
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
+            fieldIoObject,
             mesh_,
             procFields[0].dimensions(),
             internalField,
@@ -301,14 +304,14 @@ Foam::fvFieldReconstructor::reconstructFvVolumeField
 
 
 template<class Type>
-Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> >
-Foam::fvFieldReconstructor::reconstructFvSurfaceField
+Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
+Foam::fvFieldReconstructor::reconstructFvVolumeField
 (
     const IOobject& fieldIoObject
-)
+) const
 {
     // Read the field for all the processors
-    PtrList<GeometricField<Type, fvsPatchField, surfaceMesh> > procFields
+    PtrList<GeometricField<Type, fvPatchField, volMesh> > procFields
     (
         procMeshes_.size()
     );
@@ -318,7 +321,7 @@ Foam::fvFieldReconstructor::reconstructFvSurfaceField
         procFields.set
         (
             procI,
-            new GeometricField<Type, fvsPatchField, surfaceMesh>
+            new GeometricField<Type, fvPatchField, volMesh>
             (
                 IOobject
                 (
@@ -333,7 +336,29 @@ Foam::fvFieldReconstructor::reconstructFvSurfaceField
         );
     }
 
+    return reconstructFvVolumeField
+    (
+        IOobject
+        (
+            fieldIoObject.name(),
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        procFields
+    );
+}
+
 
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> >
+Foam::fvFieldReconstructor::reconstructFvSurfaceField
+(
+    const IOobject& fieldIoObject,
+    const PtrList<GeometricField<Type, fvsPatchField, surfaceMesh> >& procFields
+) const
+{
     // Create the internalField
     Field<Type> internalField(mesh_.nInternalFaces());
 
@@ -503,14 +528,7 @@ Foam::fvFieldReconstructor::reconstructFvSurfaceField
     (
         new GeometricField<Type, fvsPatchField, surfaceMesh>
         (
-            IOobject
-            (
-                fieldIoObject.name(),
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
+            fieldIoObject,
             mesh_,
             procFields[0].dimensions(),
             internalField,
@@ -520,6 +538,54 @@ Foam::fvFieldReconstructor::reconstructFvSurfaceField
 }
 
 
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> >
+Foam::fvFieldReconstructor::reconstructFvSurfaceField
+(
+    const IOobject& fieldIoObject
+) const
+{
+    // Read the field for all the processors
+    PtrList<GeometricField<Type, fvsPatchField, surfaceMesh> > procFields
+    (
+        procMeshes_.size()
+    );
+
+    forAll(procMeshes_, procI)
+    {
+        procFields.set
+        (
+            procI,
+            new GeometricField<Type, fvsPatchField, surfaceMesh>
+            (
+                IOobject
+                (
+                    fieldIoObject.name(),
+                    procMeshes_[procI].time().timeName(),
+                    procMeshes_[procI],
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE
+                ),
+                procMeshes_[procI]
+            )
+        );
+    }
+
+    return reconstructFvSurfaceField
+    (
+        IOobject
+        (
+            fieldIoObject.name(),
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        procFields
+    );
+}
+
+
 template<class Type>
 void Foam::fvFieldReconstructor::reconstructFvVolumeInternalFields
 (