diff --git a/applications/utilities/preProcessing/createExternalCoupledPatchGeometry/addRegionsOption.H b/applications/utilities/preProcessing/createExternalCoupledPatchGeometry/addRegionsOption.H
new file mode 100644
index 0000000000..9f57e35a40
--- /dev/null
+++ b/applications/utilities/preProcessing/createExternalCoupledPatchGeometry/addRegionsOption.H
@@ -0,0 +1,10 @@
+//
+// addRegionOption.H
+// ~~~~~~~~~~~~~~~~~
+
+    Foam::argList::addOption
+    (
+        "regions",
+        "(name1 .. nameN)",
+        "specify alternative mesh regions"
+    );
diff --git a/applications/utilities/preProcessing/createExternalCoupledPatchGeometry/createExternalCoupledPatchGeometry.C b/applications/utilities/preProcessing/createExternalCoupledPatchGeometry/createExternalCoupledPatchGeometry.C
index 51818b3787..6a8b78c683 100644
--- a/applications/utilities/preProcessing/createExternalCoupledPatchGeometry/createExternalCoupledPatchGeometry.C
+++ b/applications/utilities/preProcessing/createExternalCoupledPatchGeometry/createExternalCoupledPatchGeometry.C
@@ -38,6 +38,11 @@ Usage
     \param -region \<name\> \n
     Specify an alternative mesh region.
 
+    \param -regions (\<name1\> \<name2\> .. \<namen\>) \n
+    Specify alternative mesh regions. The region names will be sorted
+    alphabetically and a single composite name will be created
+        \<nameX\>_\<nameY\>.._\<nameZ\>
+
     On execution, the combined patch geometry (points and faces) are output
     to the communications directory.
 
@@ -59,6 +64,7 @@ SeeAlso
 int main(int argc, char *argv[])
 {
     #include "addRegionOption.H"
+    #include "addRegionsOption.H"
     argList::validArgs.append("patchGroup");
     argList::addOption
     (
@@ -68,16 +74,52 @@ int main(int argc, char *argv[])
     );
     #include "setRootCase.H"
     #include "createTime.H"
-    #include "createNamedMesh.H"
 
-    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    wordList regionNames(1, fvMesh::defaultRegion);
+    if (!args.optionReadIfPresent("region", regionNames[0]))
+    {
+        args.optionReadIfPresent("regions", regionNames);
+    }
 
-    const wordRe patchGroup(args[1]);
+    const wordRe patchGroup(args.argRead<wordRe>(1));
 
     fileName commsDir(runTime.path()/"comms");
     args.optionReadIfPresent("commsDir", commsDir);
 
-    externalCoupledFunctionObject::writeGeometry(mesh, commsDir, patchGroup);
+
+    // Make sure region names are in canonical order
+    stableSort(regionNames);
+
+
+    PtrList<const fvMesh> meshes(regionNames.size());
+    forAll(regionNames, i)
+    {
+        Info<< "Create mesh " << regionNames[i] << " for time = "
+            << runTime.timeName() << nl << endl;
+
+        meshes.set
+        (
+            i,
+            new fvMesh
+            (
+                Foam::IOobject
+                (
+                    regionNames[i],
+                    runTime.timeName(),
+                    runTime,
+                    Foam::IOobject::MUST_READ
+                )
+            )
+        );
+    }
+
+
+    externalCoupledFunctionObject::writeGeometry
+    (
+        UPtrList<const fvMesh>(meshes),
+        commsDir,
+        patchGroup
+    );
 
     Info<< "\nEnd\n" << endl;
 
diff --git a/src/postProcessing/functionObjects/jobControl/externalCoupled/externalCoupledFunctionObject.C b/src/postProcessing/functionObjects/jobControl/externalCoupled/externalCoupledFunctionObject.C
index ca7e2d3951..ca9905c037 100644
--- a/src/postProcessing/functionObjects/jobControl/externalCoupled/externalCoupledFunctionObject.C
+++ b/src/postProcessing/functionObjects/jobControl/externalCoupled/externalCoupledFunctionObject.C
@@ -31,6 +31,7 @@ License
 #include "volFields.H"
 #include "globalIndex.H"
 #include "fvMesh.H"
+#include "DynamicField.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -65,11 +66,16 @@ Foam::fileName Foam::externalCoupledFunctionObject::baseDir() const
 Foam::fileName Foam::externalCoupledFunctionObject::groupDir
 (
     const fileName& commsDir,
-    const word& regionName,
+    const word& regionGroupName,
     const wordRe& groupName
 )
 {
-    fileName result(commsDir/regionName/string::validate<fileName>(groupName));
+    fileName result
+    (
+        commsDir
+       /regionGroupName
+       /string::validate<fileName>(groupName)
+    );
     result.clean();
 
     return result;
@@ -126,11 +132,11 @@ void Foam::externalCoupledFunctionObject::removeReadFiles() const
 
     if (log_) Info<< type() << ": removing all read files" << endl;
 
-    forAll(regionNames_, regionI)
+    forAll(regionGroupNames_, regionI)
     {
-        const word& regionName = regionNames_[regionI];
-        const fvMesh& mesh = time_.lookupObject<fvMesh>(regionName);
-        const labelList& groups = regionToGroups_[regionName];
+        const word& compName = regionGroupNames_[regionI];
+
+        const labelList& groups = regionToGroups_[compName];
         forAll(groups, i)
         {
             label groupI = groups[i];
@@ -141,7 +147,7 @@ void Foam::externalCoupledFunctionObject::removeReadFiles() const
                 const word& fieldName = groupReadFields_[groupI][fieldI];
                 rm
                 (
-                    groupDir(commsDir_, mesh.dbDir(), groupName)
+                    groupDir(commsDir_, compName, groupName)
                   / fieldName + ".in"
                 );
             }
@@ -159,22 +165,22 @@ void Foam::externalCoupledFunctionObject::removeWriteFiles() const
 
     if (log_) Info<< type() << ": removing all write files" << endl;
 
-    forAll(regionNames_, regionI)
+    forAll(regionGroupNames_, regionI)
     {
-        const word& regionName = regionNames_[regionI];
-        const fvMesh& mesh = time_.lookupObject<fvMesh>(regionName);
-        const labelList& groups = regionToGroups_[regionName];
+        const word& compName = regionGroupNames_[regionI];
+
+        const labelList& groups = regionToGroups_[compName];
         forAll(groups, i)
         {
             label groupI = groups[i];
             const wordRe& groupName = groupNames_[groupI];
 
-            forAll(groupWriteFields_[groupI], fieldI)
+            forAll(groupReadFields_[groupI], fieldI)
             {
-                const word& fieldName = groupWriteFields_[groupI][fieldI];
+                const word& fieldName = groupReadFields_[groupI][fieldI];
                 rm
                 (
-                    groupDir(commsDir_, mesh.dbDir(), groupName)
+                    groupDir(commsDir_, compName, groupName)
                   / fieldName + ".out"
                 );
             }
@@ -376,12 +382,21 @@ void Foam::externalCoupledFunctionObject::readLines
 
 void Foam::externalCoupledFunctionObject::writeGeometry
 (
-    const fvMesh& mesh,
+    const UPtrList<const fvMesh>& meshes,
     const fileName& commsDir,
     const wordRe& groupName
 )
 {
-    fileName dir(groupDir(commsDir, mesh.dbDir(), groupName));
+    wordList regionNames(meshes.size());
+    forAll(meshes, i)
+    {
+        regionNames[i] = meshes[i].dbDir();
+    }
+
+    // Make sure meshes are provided in sorted order
+    checkOrder(regionNames);
+
+    fileName dir(groupDir(commsDir, compositeName(regionNames), groupName));
 
     //if (log_)
     {
@@ -397,99 +412,210 @@ void Foam::externalCoupledFunctionObject::writeGeometry
         osFacesPtr.reset(new OFstream(dir/"patchFaces"));
     }
 
-    const labelList patchIDs
-    (
-        mesh.boundaryMesh().patchSet
-        (
-            List<wordRe>(1, groupName)
-        ).sortedToc()
-    );
 
-    forAll(patchIDs, i)
+    DynamicList<face> allMeshesFaces;
+    DynamicField<point> allMeshesPoints;
+
+    forAll(meshes, meshI)
     {
-        label patchI = patchIDs[i];
+        const fvMesh& mesh = meshes[meshI];
 
-        const polyPatch& p = mesh.boundaryMesh()[patchI];
+        const labelList patchIDs
+        (
+            mesh.boundaryMesh().patchSet
+            (
+                List<wordRe>(1, groupName)
+            ).sortedToc()
+        );
+
+        // Count faces
+        label nFaces = 0;
+        forAll(patchIDs, i)
+        {
+            nFaces += mesh.boundaryMesh()[patchIDs[i]].size();
+        }
+
+        // Collect faces
+        DynamicList<label> allFaceIDs(nFaces);
+        forAll(patchIDs, i)
+        {
+            const polyPatch& p = mesh.boundaryMesh()[patchIDs[i]];
+
+            forAll(p, pI)
+            {
+                allFaceIDs.append(p.start()+pI);
+            }
+        }
+
+        // Construct overall patch
+        indirectPrimitivePatch allPatch
+        (
+            IndirectList<face>(mesh.faces(), allFaceIDs),
+            mesh.points()
+        );
 
         labelList pointToGlobal;
         labelList uniquePointIDs;
         mesh.globalData().mergePoints
         (
-            p.meshPoints(),
-            p.meshPointMap(),
+            allPatch.meshPoints(),
+            allPatch.meshPointMap(),
             pointToGlobal,
             uniquePointIDs
         );
 
         label procI = Pstream::myProcNo();
 
-        List<pointField> allPoints(Pstream::nProcs());
-        allPoints[procI] = pointField(mesh.points(), uniquePointIDs);
-        Pstream::gatherList(allPoints);
+        List<pointField> collectedPoints(Pstream::nProcs());
+        collectedPoints[procI] = pointField(mesh.points(), uniquePointIDs);
+        Pstream::gatherList(collectedPoints);
 
-        List<faceList> allFaces(Pstream::nProcs());
-        faceList& patchFaces = allFaces[procI];
-        patchFaces = p.localFaces();
+        List<faceList> collectedFaces(Pstream::nProcs());
+        faceList& patchFaces = collectedFaces[procI];
+        patchFaces = allPatch.localFaces();
         forAll(patchFaces, faceI)
         {
             inplaceRenumber(pointToGlobal, patchFaces[faceI]);
         }
-        Pstream::gatherList(allFaces);
+        Pstream::gatherList(collectedFaces);
 
         if (Pstream::master())
         {
-            pointField pts
-            (
-                ListListOps::combine<pointField>
-                (
-                    allPoints,
-                    accessOp<pointField>()
-                )
-            );
+            // Append and renumber
+            label nPoints = allMeshesPoints.size();
 
-            //if (log_)
+            forAll(collectedPoints, procI)
             {
-                Info<< typeName << ": for patch " << p.name()
-                    << " writing " << pts.size() << " points to "
-                    << osPointsPtr().name() << endl;
+                allMeshesPoints.append(collectedPoints[procI]);
+
             }
-
-            // Write points
-            osPointsPtr() << patchKey.c_str() << p.name() << pts << endl;
-
-            faceList fcs
-            (
-                ListListOps::combine<faceList>(allFaces, accessOp<faceList>())
-            );
-
-            //if (log_)
+            face newFace;
+            forAll(collectedFaces, procI)
             {
-                Info<< typeName << ": for patch " << p.name()
-                    << " writing " << fcs.size() << " faces to "
-                    << osFacesPtr().name() << endl;
-            }
+                const faceList& procFaces = collectedFaces[procI];
 
-            // Write faces
-            osFacesPtr() << patchKey.c_str() << p.name() << fcs << endl;
+                forAll(procFaces, faceI)
+                {
+                    const face& f = procFaces[faceI];
+
+                    newFace.setSize(f.size());
+                    forAll(f, fp)
+                    {
+                        newFace[fp] = f[fp]+nPoints;
+                    }
+                    allMeshesFaces.append(newFace);
+                }
+
+                nPoints += collectedPoints[procI].size();
+            }
         }
+
+        //if (log_)
+        {
+            Info<< typeName << ": for mesh " << mesh.name()
+                << " writing " << allMeshesPoints.size() << " points to "
+                << osPointsPtr().name() << endl;
+            Info<< typeName << ": for mesh " << mesh.name()
+                << " writing " << allMeshesFaces.size() << " faces to "
+                << osFacesPtr().name() << endl;
+        }
+    }
+
+    // Write points
+    if (osPointsPtr.valid())
+    {
+        osPointsPtr() << allMeshesPoints << endl;
+    }
+
+    // Write faces
+    if (osFacesPtr.valid())
+    {
+        osFacesPtr() << allMeshesFaces << endl;
+    }
+}
+
+
+Foam::word Foam::externalCoupledFunctionObject::compositeName
+(
+    const wordList& regionNames
+)
+{
+    if (regionNames.size() == 0)
+    {
+        FatalErrorIn
+        (
+            "externalCoupledFunctionObject::compositeName(const wordList&)"
+        )   << "Empty regionNames" << abort(FatalError);
+        return word::null;
+    }
+    else if (regionNames.size() == 1)
+    {
+        if (regionNames[0] == polyMesh::defaultRegion)
+        {
+            // For compatibility with single region cases suppress single
+            // region name
+            return word("");
+        }
+        else
+        {
+            return regionNames[0];
+        }
+    }
+    else
+    {
+        // Enforce lexical ordering
+        checkOrder(regionNames);
+
+        word composite(regionNames[0]);
+        for (label i = 1; i < regionNames.size(); i++)
+        {
+            composite += "_" + regionNames[i];
+        }
+
+        return composite;
+    }
+}
+
+
+void Foam::externalCoupledFunctionObject::checkOrder
+(
+    const wordList& regionNames
+)
+{
+    labelList order;
+    sortedOrder(regionNames, order);
+    if (order != identity(regionNames.size()))
+    {
+        FatalErrorIn
+        (
+            "externalCoupledFunctionObject::checkOrder(const wordList&)"
+        )   << "regionNames " << regionNames << " not in alphabetical order :"
+            << order << exit(FatalError);
     }
 }
 
 
 void Foam::externalCoupledFunctionObject::readData()
 {
-    forAll(regionNames_, regionI)
+    forAll(regionGroupNames_, regionI)
     {
-        const word& regionName = regionNames_[regionI];
-        const labelList& groups = regionToGroups_[regionName];
+        const word& compName = regionGroupNames_[regionI];
+        const wordList& regionNames = regionGroupRegions_[regionI];
 
-        const fvMesh& mesh = time_.lookupObject<fvMesh>(regionName);
+        // Get the meshes for the region-group
+        UPtrList<const fvMesh> meshes(regionNames.size());
+        forAll(regionNames, j)
+        {
+            const word& regionName = regionNames[j];
+            meshes.set(j, &time_.lookupObject<fvMesh>(regionName));
+        }
+
+        const labelList& groups = regionToGroups_[compName];
 
         forAll(groups, i)
         {
             label groupI = groups[i];
             const wordRe& groupName = groupNames_[groupI];
-            const labelList& patchIDs = groupPatchIDs_[groupI];
             const wordList& fieldNames = groupReadFields_[groupI];
 
             forAll(fieldNames, fieldI)
@@ -498,37 +624,32 @@ void Foam::externalCoupledFunctionObject::readData()
 
                 bool ok = readData<scalar>
                 (
-                    mesh,
+                    meshes,
                     groupName,
-                    patchIDs,
                     fieldName
                 );
                 ok = ok || readData<vector>
                 (
-                    mesh,
+                    meshes,
                     groupName,
-                    patchIDs,
                     fieldName
                 );
                 ok = ok || readData<sphericalTensor>
                 (
-                    mesh,
+                    meshes,
                     groupName,
-                    patchIDs,
                     fieldName
                 );
                 ok = ok || readData<symmTensor>
                 (
-                    mesh,
+                    meshes,
                     groupName,
-                    patchIDs,
                     fieldName
                 );
                 ok = ok || readData<tensor>
                 (
-                    mesh,
+                    meshes,
                     groupName,
-                    patchIDs,
                     fieldName
                 );
 
@@ -538,7 +659,7 @@ void Foam::externalCoupledFunctionObject::readData()
                     (
                         "void Foam::externalCoupledFunctionObject::readData()"
                     )
-                        << "Field " << fieldName << " in region " << mesh.name()
+                        << "Field " << fieldName << " in regions " << compName
                         << " was not found." << endl;
                 }
             }
@@ -549,56 +670,59 @@ void Foam::externalCoupledFunctionObject::readData()
 
 void Foam::externalCoupledFunctionObject::writeData() const
 {
-    forAll(regionNames_, regionI)
+    forAll(regionGroupNames_, regionI)
     {
-        const word& regionName = regionNames_[regionI];
-        const labelList& groups = regionToGroups_[regionName];
+        const word& compName = regionGroupNames_[regionI];
+        const wordList& regionNames = regionGroupRegions_[regionI];
 
-        const fvMesh& mesh = time_.lookupObject<fvMesh>(regionName);
+        // Get the meshes for the region-group
+        UPtrList<const fvMesh> meshes(regionNames.size());
+        forAll(regionNames, j)
+        {
+            const word& regionName = regionNames[j];
+            meshes.set(j, &time_.lookupObject<fvMesh>(regionName));
+        }
+
+        const labelList& groups = regionToGroups_[compName];
 
         forAll(groups, i)
         {
             label groupI = groups[i];
             const wordRe& groupName = groupNames_[groupI];
-            const labelList& patchIDs = groupPatchIDs_[groupI];
             const wordList& fieldNames = groupWriteFields_[groupI];
 
             forAll(fieldNames, fieldI)
             {
                 const word& fieldName = fieldNames[fieldI];
+
                 bool ok = writeData<scalar>
                 (
-                    mesh,
+                    meshes,
                     groupName,
-                    patchIDs,
                     fieldName
                 );
                 ok = ok || writeData<vector>
                 (
-                    mesh,
+                    meshes,
                     groupName,
-                    patchIDs,
                     fieldName
                 );
                 ok = ok || writeData<sphericalTensor>
                 (
-                    mesh,
+                    meshes,
                     groupName,
-                    patchIDs,
                     fieldName
                 );
                 ok = ok || writeData<symmTensor>
                 (
-                    mesh,
+                    meshes,
                     groupName,
-                    patchIDs,
                     fieldName
                 );
                 ok = ok || writeData<tensor>
                 (
-                    mesh,
+                    meshes,
                     groupName,
-                    patchIDs,
                     fieldName
                 );
 
@@ -608,7 +732,7 @@ void Foam::externalCoupledFunctionObject::writeData() const
                     (
                         "void Foam::externalCoupledFunctionObject::writeData()"
                     )
-                        << "Field " << fieldName << " in region " << mesh.name()
+                        << "Field " << fieldName << " in regions " << compName
                         << " was not found." << endl;
                 }
             }
@@ -625,12 +749,20 @@ void Foam::externalCoupledFunctionObject::initialise()
     }
 
     // Write the geometry if not already there
-    forAll(regionNames_, regionI)
+    forAll(regionGroupRegions_, i)
     {
-        const word& regionName = regionNames_[regionI];
-        const labelList& groups = regionToGroups_[regionName];
+        const word& compName = regionGroupNames_[i];
+        const wordList& regionNames = regionGroupRegions_[i];
 
-        const fvMesh& mesh = time_.lookupObject<fvMesh>(regionName);
+        // Get the meshes for the region-group
+        UPtrList<const fvMesh> meshes(regionNames.size());
+        forAll(regionNames, j)
+        {
+            const word& regionName = regionNames[j];
+            meshes.set(j, &time_.lookupObject<fvMesh>(regionName));
+        }
+
+        const labelList& groups = regionToGroups_[compName];
 
         forAll(groups, i)
         {
@@ -640,14 +772,16 @@ void Foam::externalCoupledFunctionObject::initialise()
             bool exists = false;
             if (Pstream::master())
             {
-                fileName dir(groupDir(commsDir_, mesh.dbDir(), groupName));
+                fileName dir(groupDir(commsDir_, compName, groupName));
 
-                exists = isFile(dir/"patchPoints") || isFile(dir/"patchFaces");
+                exists =
+                    isFile(dir/"patchPoints")
+                 || isFile(dir/"patchFaces");
             }
 
             if (!returnReduce(exists, orOp<bool>()))
             {
-                writeGeometry(mesh, commsDir_, groupName);
+                writeGeometry(meshes, commsDir_, groupName);
             }
         }
     }
@@ -802,6 +936,12 @@ bool Foam::externalCoupledFunctionObject::read(const dictionary& dict)
     initByExternal_ = readBool(dict.lookup("initByExternal"));
     log_ = dict.lookupOrDefault("log", false);
 
+
+    // Get names of all fvMeshes (and derived types)
+    wordList allRegionNames(time_.lookupClass<fvMesh>().sortedToc());
+
+
+
     const dictionary& allRegionsDict = dict.subDict("regions");
 
     forAllConstIter(dictionary, allRegionsDict, iter)
@@ -818,9 +958,16 @@ bool Foam::externalCoupledFunctionObject::read(const dictionary& dict)
                 << exit(FatalIOError);
         }
 
-        const word& regionName = iter().keyword();
+        const wordRe regionGroupName(iter().keyword());
         const dictionary& regionDict = iter().dict();
-        regionNames_.append(regionName);
+
+        labelList regionIDs = findStrings(regionGroupName, allRegionNames);
+
+        const wordList regionNames(allRegionNames, regionIDs);
+
+        regionGroupNames_.append(compositeName(regionNames));
+        regionGroupRegions_.append(regionNames);
+
 
         forAllConstIter(dictionary, regionDict, regionIter)
         {
@@ -844,7 +991,7 @@ bool Foam::externalCoupledFunctionObject::read(const dictionary& dict)
 
             HashTable<labelList>::iterator fnd = regionToGroups_.find
             (
-                regionName
+                regionGroupNames_.last()
             );
             if (fnd != regionToGroups_.end())
             {
@@ -852,21 +999,15 @@ bool Foam::externalCoupledFunctionObject::read(const dictionary& dict)
             }
             else
             {
-                regionToGroups_.insert(regionName, labelList(1, nGroups));
+                regionToGroups_.insert
+                (
+                    regionGroupNames_.last(),
+                    labelList(1, nGroups)
+                );
             }
             groupNames_.append(groupName);
             groupReadFields_.append(readFields);
             groupWriteFields_.append(writeFields);
-
-            // Pre-calculate the patchIDs
-            const fvMesh& mesh = time_.lookupObject<fvMesh>(regionName);
-            groupPatchIDs_.append
-            (
-                mesh.boundaryMesh().patchSet
-                (
-                    List<wordRe>(1, groupName)
-                ).sortedToc()
-            );
         }
     }
 
@@ -875,25 +1016,26 @@ bool Foam::externalCoupledFunctionObject::read(const dictionary& dict)
     if (log_)
     {
         Info<< type() << ": Communicating with regions:" << endl;
-        forAll(regionNames_, regionI)
+        forAll(regionGroupNames_, rgI)
         {
-            const word& regionName = regionNames_[regionI];
-            const fvMesh& mesh = time_.lookupObject<fvMesh>(regionName);
+            //const wordList& regionNames = regionGroupRegions_[rgI];
+            const word& compName = regionGroupNames_[rgI];
 
-            Info<< "Region: " << mesh.name() << endl << incrIndent;
-            const labelList& groups = regionToGroups_[regionName];
+            Info<< "Region: " << compName << endl << incrIndent;
+            const labelList& groups = regionToGroups_[compName];
             forAll(groups, i)
             {
                 label groupI = groups[i];
                 const wordRe& groupName = groupNames_[groupI];
-                const labelList& patchIDs = groupPatchIDs_[groupI];
 
-                Info<< indent << "Group: " << groupName << "\t"
-                    << " patches: " << patchIDs << endl
-                    << incrIndent
-                    << indent << "Reading fields: " << groupReadFields_[groupI]
+                Info<< indent << "patchGroup: " << groupName << "\t"
                     << endl
-                    << indent << "Writing fields: " << groupWriteFields_[groupI]
+                    << incrIndent
+                    << indent << "Reading fields: "
+                    << groupReadFields_[groupI]
+                    << endl
+                    << indent << "Writing fields: "
+                    << groupWriteFields_[groupI]
                     << endl
                     << decrIndent;
             }
@@ -907,25 +1049,24 @@ bool Foam::externalCoupledFunctionObject::read(const dictionary& dict)
     //       should already be written - but just make sure
     if (Pstream::master())
     {
-        forAll(regionNames_, regionI)
+        forAll(regionGroupNames_, rgI)
         {
-            const word& regionName = regionNames_[regionI];
-            const fvMesh& mesh = time_.lookupObject<fvMesh>(regionName);
-            const labelList& groups = regionToGroups_[regionName];
+            const word& compName = regionGroupNames_[rgI];
+
+            const labelList& groups = regionToGroups_[compName];
             forAll(groups, i)
             {
                 label groupI = groups[i];
                 const wordRe& groupName = groupNames_[groupI];
 
-                fileName dir(groupDir(commsDir_, mesh.dbDir(), groupName));
+                fileName dir(groupDir(commsDir_, compName, groupName));
                 if (!isDir(dir))
                 {
                     if (log_)
                     {
                         Info<< type() << ": creating communications directory "
-                        << dir << endl;
+                            << dir << endl;
                     }
-
                     mkDir(dir);
                 }
             }
diff --git a/src/postProcessing/functionObjects/jobControl/externalCoupled/externalCoupledFunctionObject.H b/src/postProcessing/functionObjects/jobControl/externalCoupled/externalCoupledFunctionObject.H
index 7806e64d64..5985f1f587 100644
--- a/src/postProcessing/functionObjects/jobControl/externalCoupled/externalCoupledFunctionObject.H
+++ b/src/postProcessing/functionObjects/jobControl/externalCoupled/externalCoupledFunctionObject.H
@@ -48,7 +48,10 @@ Description
     which gets read/written on the master processor only. In the
     communications directory the structure will be
 
-        <regionName>/<patchGroup>/<fieldName>.[in|out]
+        <regionsName>/<patchGroup>/<fieldName>.[in|out]
+
+    (where regionsName is either the name of a single region or a composite
+     of multiple region names)
 
     At start-up, the boundary creates a lock file, i.e..
 
@@ -58,13 +61,13 @@ Description
     execution the boundary values are written to files (one per region,
     per patch(group), per field), e.g.
 
-        <regionName>/<patchGroup>/<fieldName>.out
+        <regionsName>/<patchGroup>/<fieldName>.out
 
     The lock file is then removed, instructing the external source to take
     control of the program execution. When ready, the external program
     should create the return values, e.g. to files
 
-        <regionName>/<patchGroup>/<fieldName>.in
+        <regionsName>/<patchGroup>/<fieldName>.in
 
     ... and then re-instate the lock file. The functionObject will then
     read these values, apply them to the boundary conditions and pass
@@ -82,7 +85,7 @@ Description
 
         regions
         {
-            region0
+            "(region1|region0)"         // Name of region(s)
             {
                 TPatchGroup             // Name of patch(group)
                 {
@@ -95,7 +98,7 @@ Description
     \endverbatim
 
     This reads/writes (on the master processor) the directory:
-        comms/region0/TPatchGroup/
+        comms/region0_region1/TPatchGroup/
     with contents:
         patchPoints     (collected points)
         patchFaces      (collected faces)
@@ -120,6 +123,7 @@ SourceFiles
 #include "wordReList.H"
 #include "scalarField.H"
 #include "Switch.H"
+#include "UPtrList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -167,18 +171,18 @@ class externalCoupledFunctionObject
         //- Log flag
         bool log_;
 
-        //- Names of regions
-        DynamicList<word> regionNames_;
+        //- Names of (composite) regions
+        DynamicList<word> regionGroupNames_;
 
-        // Per region the indices of the group information
+        // Per (composite) region the names of the regions
+        DynamicList<wordList> regionGroupRegions_;
+
+        // Per (composite) region the indices of the group information
         HashTable<labelList> regionToGroups_;
 
         // Per group the names of the patches/patchGroups
         DynamicList<wordRe> groupNames_;
 
-        // Per group the indices of the patches
-        DynamicList<labelList> groupPatchIDs_;
-
         // Per group the names of the fields to read
         DynamicList<wordList> groupReadFields_;
 
@@ -195,7 +199,7 @@ class externalCoupledFunctionObject
         static fileName groupDir
         (
             const fileName& commsDir,
-            const word& regionName,
+            const word& regionsName,
             const wordRe& groupName
         );
 
@@ -225,9 +229,8 @@ class externalCoupledFunctionObject
         template<class Type>
         bool readData
         (
-            const fvMesh& mesh,
+            const UPtrList<const fvMesh>& meshes,
             const wordRe& groupName,
-            const labelList& patchIDs,
             const word& fieldName
         );
         //- Read data for all regions, all fields
@@ -237,9 +240,8 @@ class externalCoupledFunctionObject
         template<class Type>
         bool writeData
         (
-            const fvMesh& mesh,
+            const UPtrList<const fvMesh>& meshes,
             const wordRe& groupName,
-            const labelList& patchIDs,
             const word& fieldName
         ) const;
 
@@ -275,6 +277,7 @@ class externalCoupledFunctionObject
         template<class Type>
         static tmp<Field<Type> > gatherAndCombine(const Field<Type>& fld);
 
+        static void checkOrder(const wordList&);
 
         //- Disallow default bitwise copy construc
         externalCoupledFunctionObject(const externalCoupledFunctionObject&);
@@ -356,10 +359,14 @@ public:
 
         // Other
 
+            //- Create single name by appending words (in sorted order),
+            //  separated by '_'
+            static word compositeName(const wordList&);
+
             //- Write geometry for the group/patch
             static void writeGeometry
             (
-                const fvMesh& mesh,
+                const UPtrList<const fvMesh>& meshes,
                 const fileName& commsDir,
                 const wordRe& groupName
             );
diff --git a/src/postProcessing/functionObjects/jobControl/externalCoupled/externalCoupledFunctionObjectTemplates.C b/src/postProcessing/functionObjects/jobControl/externalCoupled/externalCoupledFunctionObjectTemplates.C
index 88477c5bb6..f483ecb71e 100644
--- a/src/postProcessing/functionObjects/jobControl/externalCoupled/externalCoupledFunctionObjectTemplates.C
+++ b/src/postProcessing/functionObjects/jobControl/externalCoupled/externalCoupledFunctionObjectTemplates.C
@@ -40,25 +40,20 @@ License
 template<class Type>
 bool Foam::externalCoupledFunctionObject::readData
 (
-    const fvMesh& mesh,
+    const UPtrList<const fvMesh>& meshes,
     const wordRe& groupName,
-    const labelList& patchIDs,
     const word& fieldName
 )
 {
     typedef GeometricField<Type, fvPatchField, volMesh> volFieldType;
     typedef externalCoupledMixedFvPatchField<Type> patchFieldType;
 
-    if (!mesh.foundObject<volFieldType>(fieldName))
+    wordList regionNames(meshes.size());
+    forAll(meshes, i)
     {
-        return false;
+        regionNames[i] = meshes[i].dbDir();
     }
 
-    const volFieldType& cvf = mesh.lookupObject<volFieldType>(fieldName);
-    const typename volFieldType::GeometricBoundaryField& bf =
-        cvf.boundaryField();
-
-
     // File only opened on master; contains data for all processors, for all
     // patchIDs.
     autoPtr<IFstream> masterFilePtr;
@@ -66,7 +61,7 @@ bool Foam::externalCoupledFunctionObject::readData
     {
         const fileName transferFile
         (
-            groupDir(commsDir_, mesh.dbDir(), groupName)
+            groupDir(commsDir_, compositeName(regionNames), groupName)
           / fieldName + ".in"
         );
 
@@ -82,181 +77,234 @@ bool Foam::externalCoupledFunctionObject::readData
             (
                 "void externalCoupledFunctionObject::readData"
                 "("
-                    "const fvMesh&, "
+                    "const UPtrList<const fvMesh>&, "
                     "const wordRe&, "
-                    "const labelList&, "
                     "const word&"
                ")",
                 masterFilePtr()
-            )   << "Cannot open file for region " << mesh.name()
-                << ", field " << fieldName << ", patches " << patchIDs
+            )   << "Cannot open file for region " << compositeName(regionNames)
+                << ", field " << fieldName
                 << exit(FatalIOError);
         }
     }
 
-    // Handle column-wise reading of patch data. Supports most easy types
-    forAll(patchIDs, i)
+
+    label nFound = 0;
+
+
+    forAll(meshes, i)
     {
-        label patchI = patchIDs[i];
+        const fvMesh& mesh = meshes[i];
 
-        if (isA<patchFieldType>(bf[patchI]))
+        if (!mesh.foundObject<volFieldType>(fieldName))
         {
-            // Explicit handling of externalCoupledObjectMixed bcs - they have
-            // specialised reading routines.
-
-            patchFieldType& pf = const_cast<patchFieldType&>
-            (
-                refCast<const patchFieldType>
-                (
-                    bf[patchI]
-                )
-            );
-
-            // Read from master into local stream
-            OStringStream os;
-            readLines
-            (
-                bf[patchI].size(),      // number of lines to read
-                masterFilePtr,
-                os
-            );
-
-            // Pass responsability for all reading over to bc
-            pf.readData(IStringStream(os.str())());
-
-            // Update the value from the read coefficicient. Bypass any
-            // additional processing by derived type.
-            pf.patchFieldType::evaluate();
+            continue;
         }
-        else if (isA<mixedFvPatchField<Type> >(bf[patchI]))
+
+        nFound++;
+
+        const volFieldType& cvf = mesh.lookupObject<volFieldType>(fieldName);
+        const typename volFieldType::GeometricBoundaryField& bf =
+            cvf.boundaryField();
+
+
+        // Get the patches
+        const labelList patchIDs
+        (
+            mesh.boundaryMesh().patchSet
+            (
+                List<wordRe>(1, groupName)
+            ).sortedToc()
+        );
+
+        // Handle column-wise reading of patch data. Supports most easy types
+        forAll(patchIDs, i)
         {
-            // Read columns from file for
-            // value, snGrad, refValue, refGrad, valueFraction
-            List<scalarField> data;
-            readColumns
-            (
-                bf[patchI].size(),              // number of lines to read
-                4*pTraits<Type>::nComponents+1, // nColumns: 4*Type + 1*scalar
-                masterFilePtr,
-                data
-            );
+            label patchI = patchIDs[i];
 
-            mixedFvPatchField<Type>& pf = const_cast<mixedFvPatchField<Type>&>
-            (
-                refCast<const mixedFvPatchField<Type> >
+            if (isA<patchFieldType>(bf[patchI]))
+            {
+                // Explicit handling of externalCoupledMixed bcs - they
+                // have specialised reading routines.
+
+                patchFieldType& pf = const_cast<patchFieldType&>
                 (
-                    bf[patchI]
-                )
-            );
+                    refCast<const patchFieldType>
+                    (
+                        bf[patchI]
+                    )
+                );
 
-            // Transfer read data to bc.
-            // Skip value, snGrad
-            direction columnI = 2*pTraits<Type>::nComponents;
-
-            Field<Type>& refValue = pf.refValue();
-            for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
-            {
-                refValue.replace(cmpt, data[columnI++]);
-            }
-            Field<Type>& refGrad = pf.refGrad();
-            for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
-            {
-                refGrad.replace(cmpt, data[columnI++]);
-            }
-            pf.valueFraction() = data[columnI];
-
-            // Update the value from the read coefficicient. Bypass any
-            // additional processing by derived type.
-            pf.mixedFvPatchField<Type>::evaluate();
-        }
-        else if (isA<fixedGradientFvPatchField<Type> >(bf[patchI]))
-        {
-            // Read columns for value and gradient
-            List<scalarField> data;
-            readColumns
-            (
-                bf[patchI].size(),              // number of lines to read
-                2*pTraits<Type>::nComponents,   // nColumns: Type
-                masterFilePtr,
-                data
-            );
-
-            fixedGradientFvPatchField<Type>& pf =
-            const_cast<fixedGradientFvPatchField<Type>&>
-            (
-                refCast<const fixedGradientFvPatchField<Type> >
+                // Read from master into local stream
+                OStringStream os;
+                readLines
                 (
-                    bf[patchI]
-                )
-            );
+                    bf[patchI].size(),      // number of lines to read
+                    masterFilePtr,
+                    os
+                );
 
-            // Transfer gradient to bc
-            Field<Type>& gradient = pf.gradient();
-            for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
+                // Pass responsability for all reading over to bc
+                pf.readData(IStringStream(os.str())());
+
+                // Update the value from the read coefficicient. Bypass any
+                // additional processing by derived type.
+                pf.patchFieldType::evaluate();
+            }
+            else if (isA<mixedFvPatchField<Type> >(bf[patchI]))
             {
-                gradient.replace(cmpt, data[pTraits<Type>::nComponents+cmpt]);
+                // Read columns from file for
+                // value, snGrad, refValue, refGrad, valueFraction
+                List<scalarField> data;
+                readColumns
+                (
+                    bf[patchI].size(),              // number of lines to read
+                    4*pTraits<Type>::nComponents+1, // nColumns: 4*Type+1*scalar
+                    masterFilePtr,
+                    data
+                );
+
+                mixedFvPatchField<Type>& pf =
+                const_cast<mixedFvPatchField<Type>&>
+                (
+                    refCast<const mixedFvPatchField<Type> >
+                    (
+                        bf[patchI]
+                    )
+                );
+
+                // Transfer read data to bc.
+                // Skip value, snGrad
+                direction columnI = 2*pTraits<Type>::nComponents;
+
+                Field<Type>& refValue = pf.refValue();
+                for
+                (
+                    direction cmpt = 0;
+                    cmpt < pTraits<Type>::nComponents;
+                    cmpt++
+                )
+                {
+                    refValue.replace(cmpt, data[columnI++]);
+                }
+                Field<Type>& refGrad = pf.refGrad();
+                for
+                (
+                    direction cmpt = 0;
+                    cmpt < pTraits<Type>::nComponents;
+                    cmpt++
+                )
+                {
+                    refGrad.replace(cmpt, data[columnI++]);
+                }
+                pf.valueFraction() = data[columnI];
+
+                // Update the value from the read coefficicient. Bypass any
+                // additional processing by derived type.
+                pf.mixedFvPatchField<Type>::evaluate();
+            }
+            else if (isA<fixedGradientFvPatchField<Type> >(bf[patchI]))
+            {
+                // Read columns for value and gradient
+                List<scalarField> data;
+                readColumns
+                (
+                    bf[patchI].size(),              // number of lines to read
+                    2*pTraits<Type>::nComponents,   // nColumns: Type
+                    masterFilePtr,
+                    data
+                );
+
+                fixedGradientFvPatchField<Type>& pf =
+                const_cast<fixedGradientFvPatchField<Type>&>
+                (
+                    refCast<const fixedGradientFvPatchField<Type> >
+                    (
+                        bf[patchI]
+                    )
+                );
+
+                // Transfer gradient to bc
+                Field<Type>& gradient = pf.gradient();
+                for
+                (
+                    direction cmpt = 0;
+                    cmpt < pTraits<Type>::nComponents;
+                    cmpt++
+                )
+                {
+                    gradient.replace
+                    (
+                        cmpt,
+                        data[pTraits<Type>::nComponents+cmpt]
+                    );
+                }
+
+                // Update the value from the read coefficicient. Bypass any
+                // additional processing by derived type.
+                pf.fixedGradientFvPatchField<Type>::evaluate();
+            }
+            else if (isA<fixedValueFvPatchField<Type> >(bf[patchI]))
+            {
+                // Read columns for value only
+                List<scalarField> data;
+                readColumns
+                (
+                    bf[patchI].size(),              // number of lines to read
+                    pTraits<Type>::nComponents,     // number of columns to read
+                    masterFilePtr,
+                    data
+                );
+
+                // Transfer read value to bc
+                Field<Type> value(bf[patchI].size());
+                for
+                (
+                    direction cmpt = 0;
+                    cmpt < pTraits<Type>::nComponents;
+                    cmpt++
+                )
+                {
+                    value.replace(cmpt, data[cmpt]);
+                }
+
+                fixedValueFvPatchField<Type>& pf =
+                const_cast<fixedValueFvPatchField<Type>&>
+                (
+                    refCast<const fixedValueFvPatchField<Type> >
+                    (
+                        bf[patchI]
+                    )
+                );
+
+                pf == value;
+
+                // Update the value from the read coefficicient. Bypass any
+                // additional processing by derived type.
+                pf.fixedValueFvPatchField<Type>::evaluate();
+            }
+            else
+            {
+                FatalErrorIn
+                (
+                    "void externalCoupledFunctionObject::readData"
+                    "("
+                        "const UPtrList<const fvMesh>&, "
+                        "const wordRe&, "
+                        "const word&"
+                   ")"
+                )
+                    << "Unsupported boundary condition " << bf[patchI].type()
+                    << " for patch " << bf[patchI].patch().name()
+                    << " in region " << mesh.name()
+                    << exit(FatalError);
             }
 
-            // Update the value from the read coefficicient. Bypass any
-            // additional processing by derived type.
-            pf.fixedGradientFvPatchField<Type>::evaluate();
+            initialised_ = true;
         }
-        else if (isA<fixedValueFvPatchField<Type> >(bf[patchI]))
-        {
-            // Read columns for value only
-            List<scalarField> data;
-            readColumns
-            (
-                bf[patchI].size(),              // number of lines to read
-                pTraits<Type>::nComponents,     // number of columns to read
-                masterFilePtr,
-                data
-            );
-
-            // Transfer read value to bc
-            Field<Type> value(bf[patchI].size());
-            for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
-            {
-                value.replace(cmpt, data[cmpt]);
-            }
-
-            fixedValueFvPatchField<Type>& pf =
-            const_cast<fixedValueFvPatchField<Type>&>
-            (
-                refCast<const fixedValueFvPatchField<Type> >
-                (
-                    bf[patchI]
-                )
-            );
-
-            pf == value;
-
-            // Update the value from the read coefficicient. Bypass any
-            // additional processing by derived type.
-            pf.fixedValueFvPatchField<Type>::evaluate();
-        }
-        else
-        {
-            FatalErrorIn
-            (
-                "void externalCoupledFunctionObject::readData"
-                "("
-                    "const fvMesh&, "
-                    "const wordRe&, "
-                    "const labelList&, "
-                    "const word&"
-               ")"
-            )
-                << "Unsupported boundary condition " << bf[patchI].type()
-                << " for patch " << bf[patchI].patch().name()
-                << " in region " << mesh.name()
-                << exit(FatalError);
-        }
-
-        initialised_ = true;
     }
 
-    return true;
+    return nFound > 0;
 }
 
 
@@ -308,25 +356,20 @@ Foam::externalCoupledFunctionObject::gatherAndCombine
 template<class Type>
 bool Foam::externalCoupledFunctionObject::writeData
 (
-    const fvMesh& mesh,
+    const UPtrList<const fvMesh>& meshes,
     const wordRe& groupName,
-    const labelList& patchIDs,
     const word& fieldName
 ) const
 {
     typedef GeometricField<Type, fvPatchField, volMesh> volFieldType;
     typedef externalCoupledMixedFvPatchField<Type> patchFieldType;
 
-    if (!mesh.foundObject<volFieldType>(fieldName))
+    wordList regionNames(meshes.size());
+    forAll(meshes, i)
     {
-        return false;
+        regionNames[i] = meshes[i].dbDir();
     }
 
-    const volFieldType& cvf = mesh.lookupObject<volFieldType>(fieldName);
-    const typename volFieldType::GeometricBoundaryField& bf =
-        cvf.boundaryField();
-
-
     // File only opened on master; contains data for all processors, for all
     // patchIDs
     autoPtr<OFstream> masterFilePtr;
@@ -334,7 +377,7 @@ bool Foam::externalCoupledFunctionObject::writeData
     {
         const fileName transferFile
         (
-            groupDir(commsDir_, mesh.dbDir(), groupName)
+            groupDir(commsDir_, compositeName(regionNames), groupName)
           / fieldName + ".out"
         );
 
@@ -350,14 +393,13 @@ bool Foam::externalCoupledFunctionObject::writeData
             (
                 "externalCoupledFunctionObject::writeData"
                 "("
-                    "const fvMesh&, "
+                    "const UPtrList<const fvMesh>&, "
                     "const wordRe&, "
-                    "const labelList&, "
                     "const word&"
                 ") const",
                 masterFilePtr()
-            )   << "Cannot open file for region " << mesh.name()
-                << ", field " << fieldName << ", patches " << patchIDs
+            )   << "Cannot open file for region " << compositeName(regionNames)
+                << ", field " << fieldName
                 << exit(FatalIOError);
         }
     }
@@ -365,94 +407,122 @@ bool Foam::externalCoupledFunctionObject::writeData
 
     bool headerDone = false;
 
-    // Handle column-wise writing of patch data. Supports most easy types
-    forAll(patchIDs, i)
+    label nFound = 0;
+
+    forAll(meshes, i)
     {
-        label patchI = patchIDs[i];
+        const fvMesh& mesh = meshes[i];
 
-        const globalIndex globalFaces(bf[patchI].size());
-
-        if (isA<patchFieldType>(bf[patchI]))
+        if (!mesh.foundObject<volFieldType>(fieldName))
         {
-            // Explicit handling of externalCoupledObjectMixed bcs - they have
-            // specialised writing routines
+            continue;
+        }
 
-            const patchFieldType& pf = refCast<const patchFieldType>
+        nFound++;
+
+        const volFieldType& cvf = mesh.lookupObject<volFieldType>(fieldName);
+        const typename volFieldType::GeometricBoundaryField& bf =
+            cvf.boundaryField();
+
+
+        // Get the patches
+        const labelList patchIDs
+        (
+            mesh.boundaryMesh().patchSet
             (
-                bf[patchI]
-            );
-            OStringStream os;
+                List<wordRe>(1, groupName)
+            ).sortedToc()
+        );
 
-            // Pass responsibility for all writing over to bc
-            pf.writeData(os);
+        // Handle column-wise writing of patch data. Supports most easy types
+        forAll(patchIDs, i)
+        {
+            label patchI = patchIDs[i];
 
-            // Collect contributions from all processors and output them on
-            // master
-            if (Pstream::master())
+            const globalIndex globalFaces(bf[patchI].size());
+
+            if (isA<patchFieldType>(bf[patchI]))
             {
-                // Output master data first
-                if (!headerDone)
-                {
-                    pf.writeHeader(masterFilePtr());
-                    headerDone = true;
-                }
-                masterFilePtr() << os.str().c_str();
+                // Explicit handling of externalCoupledMixed bcs - they
+                // have specialised writing routines
 
-                for (label procI = 1; procI < Pstream::nProcs(); procI++)
+                const patchFieldType& pf = refCast<const patchFieldType>
+                (
+                    bf[patchI]
+                );
+                OStringStream os;
+
+                // Pass responsibility for all writing over to bc
+                pf.writeData(os);
+
+                // Collect contributions from all processors and output them on
+                // master
+                if (Pstream::master())
                 {
-                    IPstream fromSlave(Pstream::scheduled, procI);
-                    string str(fromSlave);
-                    masterFilePtr() << str.c_str();
+                    // Output master data first
+                    if (!headerDone)
+                    {
+                        pf.writeHeader(masterFilePtr());
+                        headerDone = true;
+                    }
+                    masterFilePtr() << os.str().c_str();
+
+                    for (label procI = 1; procI < Pstream::nProcs(); procI++)
+                    {
+                        IPstream fromSlave(Pstream::scheduled, procI);
+                        string str(fromSlave);
+                        masterFilePtr() << str.c_str();
+                    }
+                }
+                else
+                {
+                    OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
+                    toMaster << os.str();
+                }
+            }
+            else if (isA<mixedFvPatchField<Type> >(bf[patchI]))
+            {
+                const mixedFvPatchField<Type>& pf =
+                    refCast<const mixedFvPatchField<Type> >(bf[patchI]);
+
+                Field<Type> value(gatherAndCombine(pf));
+                Field<Type> snGrad(gatherAndCombine(pf.snGrad()()));
+                Field<Type> refValue(gatherAndCombine(pf.refValue()));
+                Field<Type> refGrad(gatherAndCombine(pf.refGrad()));
+                scalarField valueFraction(gatherAndCombine(pf.valueFraction()));
+
+                if (Pstream::master())
+                {
+                    forAll(refValue, faceI)
+                    {
+                        masterFilePtr()
+                            << value[faceI] << token::SPACE
+                            << snGrad[faceI] << token::SPACE
+                            << refValue[faceI] << token::SPACE
+                            << refGrad[faceI] << token::SPACE
+                            << valueFraction[faceI] << nl;
+                    }
                 }
             }
             else
             {
-                OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
-                toMaster << os.str();
-            }
-        }
-        else if (isA<mixedFvPatchField<Type> >(bf[patchI]))
-        {
-            const mixedFvPatchField<Type>& pf =
-                refCast<const mixedFvPatchField<Type> >(bf[patchI]);
-
-            Field<Type> value(gatherAndCombine(pf));
-            Field<Type> snGrad(gatherAndCombine(pf.snGrad()()));
-            Field<Type> refValue(gatherAndCombine(pf.refValue()));
-            Field<Type> refGrad(gatherAndCombine(pf.refGrad()));
-            scalarField valueFraction(gatherAndCombine(pf.valueFraction()));
-
-            if (Pstream::master())
-            {
-                forAll(refValue, faceI)
+                // Output the value and snGrad
+                Field<Type> value(gatherAndCombine(bf[patchI]));
+                Field<Type> snGrad(gatherAndCombine(bf[patchI].snGrad()()));
+                if (Pstream::master())
                 {
-                    masterFilePtr()
-                        << value[faceI] << token::SPACE
-                        << snGrad[faceI] << token::SPACE
-                        << refValue[faceI] << token::SPACE
-                        << refGrad[faceI] << token::SPACE
-                        << valueFraction[faceI] << nl;
-                }
-            }
-        }
-        else
-        {
-            // Output the value and snGrad
-            Field<Type> value(gatherAndCombine(bf[patchI]));
-            Field<Type> snGrad(gatherAndCombine(bf[patchI].snGrad()()));
-            if (Pstream::master())
-            {
-                forAll(value, faceI)
-                {
-                    masterFilePtr()
-                        << value[faceI] << token::SPACE
-                        << snGrad[faceI] << nl;
+                    forAll(value, faceI)
+                    {
+                        masterFilePtr()
+                            << value[faceI] << token::SPACE
+                            << snGrad[faceI] << nl;
+                    }
                 }
             }
         }
     }
 
-    return true;
+    return nFound > 0;
 }
 
 
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/0/T b/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/0/T
deleted file mode 100644
index 3bf09662bb..0000000000
--- a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/0/T
+++ /dev/null
@@ -1,56 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    location    "0";
-    object      T;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 0 0 1 0 0 0];
-
-internalField   uniform 293;
-
-boundaryField
-{
-    frontAndBack
-    {
-        type            zeroGradient;
-    }
-
-    topAndBottom
-    {
-        type            zeroGradient;
-    }
-
-    hot
-    {
-        type            externalCoupledTemperature;
-        commsDir        "${FOAM_CASE}/comms";
-        fileName        "data";
-        initByExternal  yes;
-        log             true;
-        value           uniform 307.75; // 34.6 degC
-    }
-
-    cold
-    {
-        type            externalCoupledTemperature;
-        commsDir        "${FOAM_CASE}/comms";
-        fileName        "data";
-        initByExternal  yes;
-        log             true;
-        value           uniform 288.15; // 15 degC
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/0/alphat b/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/0/alphat
deleted file mode 100644
index 9e486ba215..0000000000
--- a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/0/alphat
+++ /dev/null
@@ -1,51 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    location    "0";
-    object      alphat;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [1 -1 -1 0 0 0 0];
-
-internalField   uniform 0;
-
-boundaryField
-{
-    frontAndBack
-    {
-        type            compressible::alphatWallFunction;
-        Prt             0.85;
-        value           uniform 0;
-    }
-    topAndBottom
-    {
-        type            compressible::alphatWallFunction;
-        Prt             0.85;
-        value           uniform 0;
-    }
-    hot
-    {
-        type            compressible::alphatWallFunction;
-        Prt             0.85;
-        value           uniform 0;
-    }
-    cold
-    {
-        type            compressible::alphatWallFunction;
-        Prt             0.85;
-        value           uniform 0;
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/0/k b/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/0/k
deleted file mode 100644
index ccdab2054e..0000000000
--- a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/0/k
+++ /dev/null
@@ -1,47 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    location    "0";
-    object      k;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 2 -2 0 0 0 0];
-
-internalField   uniform 3.75e-04;
-
-boundaryField
-{
-    frontAndBack
-    {
-        type            kqRWallFunction;
-        value           uniform 3.75e-04;
-    }
-    topAndBottom
-    {
-        type            kqRWallFunction;
-        value           uniform 3.75e-04;
-    }
-    hot
-    {
-        type            kqRWallFunction;
-        value           uniform 3.75e-04;
-    }
-    cold
-    {
-        type            kqRWallFunction;
-        value           uniform 3.75e-04;
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/0/omega b/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/0/omega
deleted file mode 100644
index 481f4e64d7..0000000000
--- a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/0/omega
+++ /dev/null
@@ -1,47 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    location    "0";
-    object      omega;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 0 -1 0 0 0 0];
-
-internalField   uniform 0.12;
-
-boundaryField
-{
-    frontAndBack
-    {
-        type            omegaWallFunction;
-        value           uniform 0.12;
-    }
-    topAndBottom
-    {
-        type            omegaWallFunction;
-        value           uniform 0.12;
-    }
-    hot
-    {
-        type            omegaWallFunction;
-        value           uniform 0.12;
-    }
-    cold
-    {
-        type            omegaWallFunction;
-        value           uniform 0.12;
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/Allclean b/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/Allclean
deleted file mode 100755
index 871192082a..0000000000
--- a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/Allclean
+++ /dev/null
@@ -1,14 +0,0 @@
-#!/bin/sh
-cd ${0%/*} || exit 1    # Run from this directory
-
-# Source tutorial run functions
-. $WM_PROJECT_DIR/bin/tools/CleanFunctions
-
-cleanCase
-
-rm -rf comms
-
-killall externalSolver > /dev/null 2>&1
-
-
-# ----------------------------------------------------------------- end-of-file
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/Allrun b/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/Allrun
deleted file mode 100755
index ba2d933543..0000000000
--- a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/Allrun
+++ /dev/null
@@ -1,14 +0,0 @@
-#!/bin/sh
-cd ${0%/*} || exit 1    # Run from this directory
-
-# Source tutorial run functions
-. $WM_PROJECT_DIR/bin/tools/RunFunctions
-
-./Allrun.pre
-
-runApplication $(getApplication) &
-
-./externalSolver
-
-
-# ----------------------------------------------------------------- end-of-file
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/Allrun-parallel b/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/Allrun-parallel
deleted file mode 100755
index 8c451a124e..0000000000
--- a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/Allrun-parallel
+++ /dev/null
@@ -1,16 +0,0 @@
-#!/bin/sh
-cd ${0%/*} || exit 1    # Run from this directory
-
-# Source tutorial run functions
-. $WM_PROJECT_DIR/bin/tools/RunFunctions
-
-./Allrun.pre
-
-runApplication decomposePar
-
-runParallel $(getApplication) 4 &
-
-./externalSolver
-
-
-# ----------------------------------------------------------------- end-of-file
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/Allrun.pre b/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/Allrun.pre
deleted file mode 100755
index 0fca0f1d79..0000000000
--- a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/Allrun.pre
+++ /dev/null
@@ -1,11 +0,0 @@
-#!/bin/sh
-cd ${0%/*} || exit 1    # Run from this directory
-
-# Source tutorial run functions
-. $WM_PROJECT_DIR/bin/tools/RunFunctions
-
-runApplication blockMesh
-runApplication createExternalCoupledPatchGeometry T
-
-
-# ----------------------------------------------------------------- end-of-file
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/README b/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/README
deleted file mode 100644
index 10d7392eb4..0000000000
--- a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/README
+++ /dev/null
@@ -1,5 +0,0 @@
-Example of an explicit coupling between OpenFOAM and an external application
-using the externalCoupled boundary conditions.
-
-The case is based on the buoyantCavity tutorial case, whereby on each iteration
-the 'hot' and 'cold' patch temperatures are incremented by 1K.
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/0/T b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/0/T
new file mode 100644
index 0000000000..a2e5fd849a
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/0/T
@@ -0,0 +1,30 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 300;
+
+boundaryField
+{
+    ".*"
+    {
+        type            calculated;
+        value           uniform 300;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/0/U b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/0/U
similarity index 66%
rename from tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/0/U
rename to tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/0/U
index 2d917d036f..5c6a42bb1d 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/0/U
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/0/U
@@ -10,41 +10,21 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volVectorField;
-    location    "0";
     object      U;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 dimensions      [0 1 -1 0 0 0 0];
 
-internalField   uniform (0 0 0);
+internalField   uniform (0.01 0 0);
 
 boundaryField
 {
-    frontAndBack
+    ".*"
     {
-        type            fixedValue;
-        value           uniform (0 0 0);
-    }
-
-    topAndBottom
-    {
-        type            fixedValue;
-        value           uniform (0 0 0);
-    }
-
-    hot
-    {
-        type            fixedValue;
-        value           uniform (0 0 0);
-    }
-
-    cold
-    {
-        type            fixedValue;
-        value           uniform (0 0 0);
+        type            calculated;
+        value           uniform (0.01 0 0);
     }
 }
 
-
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/0/epsilon b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/0/epsilon
new file mode 100644
index 0000000000..e20d661451
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/0/epsilon
@@ -0,0 +1,31 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -3 0 0 0 0];
+
+internalField   uniform 0.01;
+
+boundaryField
+{
+    ".*"
+    {
+        type            calculated;
+        value           uniform 0.01;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/0/k b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/0/k
new file mode 100644
index 0000000000..351dfc5606
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/0/k
@@ -0,0 +1,31 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0.1;
+
+boundaryField
+{
+    ".*"
+    {
+        type            calculated;
+        value           uniform 0.1;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/0/p b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/0/p
similarity index 71%
rename from tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/0/p
rename to tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/0/p
index 9b1e74247b..3e97215963 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/0/p
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/0/p
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "0";
     object      p;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -21,30 +20,11 @@ internalField   uniform 1e5;
 
 boundaryField
 {
-    frontAndBack
+    ".*"
     {
         type            calculated;
-        value           $internalField;
-    }
-
-    topAndBottom
-    {
-        type            calculated;
-        value           $internalField;
-    }
-
-    hot
-    {
-        type            calculated;
-        value           $internalField;
-    }
-
-    cold
-    {
-        type            calculated;
-        value           $internalField;
+        value           uniform 1e5;
     }
 }
 
-
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/0/p_rgh b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/0/p_rgh
similarity index 70%
rename from tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/0/p_rgh
rename to tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/0/p_rgh
index 15979d020e..8673f6e9bc 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/0/p_rgh
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/0/p_rgh
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "0";
     object      p_rgh;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -21,30 +20,11 @@ internalField   uniform 1e5;
 
 boundaryField
 {
-    frontAndBack
+    ".*"
     {
-        type            fixedFluxPressure;
-        value           uniform 1e5;
-    }
-
-    topAndBottom
-    {
-        type            fixedFluxPressure;
-        value           uniform 1e5;
-    }
-
-    hot
-    {
-        type            fixedFluxPressure;
-        value           uniform 1e5;
-    }
-
-    cold
-    {
-        type            fixedFluxPressure;
+        type            calculated;
         value           uniform 1e5;
     }
 }
 
-
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/Allclean b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/Allclean
new file mode 100755
index 0000000000..82e538bb43
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/Allclean
@@ -0,0 +1,23 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+cleanCase
+rm -rf comms
+rm -rf VTK
+rm -rf constant/cellToRegion constant/polyMesh/sets
+rm -rf 0/bottomWater
+rm -rf 0/topAir
+rm -rf 0/heater
+rm -rf 0/leftSolid
+rm -rf 0/rightSolid
+rm -f 0/cellToRegion
+rm -rf constant/bottomWater/polyMesh
+rm -rf constant/topAir/polyMesh
+rm -rf constant/heater/polyMesh
+rm -rf constant/leftSolid/polyMesh
+rm -rf constant/rightSolid/polyMesh
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/Allrun b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/Allrun
new file mode 100755
index 0000000000..bfb99ee4d1
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/Allrun
@@ -0,0 +1,28 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+./Allrun.pre
+
+#-- Run on single processor
+#runApplication `getApplication` &
+# Simulated external solver
+#runApplication ./externalSolver
+
+# Decompose
+runApplication decomposePar -allRegions
+
+# Run OpenFOAM
+runParallel `getApplication` 4 &
+
+# Simulated external solver
+runApplication ./externalSolver
+
+# Reconstruct
+runApplication reconstructPar -allRegions
+
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/Allrun.pre b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/Allrun.pre
new file mode 100755
index 0000000000..d841e3cfe0
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/Allrun.pre
@@ -0,0 +1,33 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+runApplication blockMesh
+runApplication topoSet
+runApplication splitMeshRegions -cellZones -overwrite
+
+# remove fluid fields from solid regions (important for post-processing)
+for i in heater leftSolid rightSolid
+do
+   rm -f 0*/$i/{nut,alphat,epsilon,k,U,p_rgh}
+done
+
+
+for i in bottomWater topAir heater leftSolid rightSolid
+do
+   changeDictionary -region $i > log.changeDictionary.$i 2>&1
+done
+
+# Create coupling geometry
+runApplication createExternalCoupledPatchGeometry \
+    -regions '(topAir heater)' coupleGroup
+
+echo
+echo "creating files for paraview post-processing"
+echo
+paraFoam -touchAll
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/README.txt b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/README.txt
new file mode 100644
index 0000000000..b0ca34a873
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/README.txt
@@ -0,0 +1,3 @@
+Modification of the heatTransfer chtMultiRegionFoam tutorial that demonstrates
+the externalCoupled functionObject in combination with the ./externalSolver
+script to simulate coupling to an external code.
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/constant/g b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/bottomWater/g
similarity index 96%
rename from tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/constant/g
rename to tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/bottomWater/g
index 0cc222ca34..6f32e33835 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/constant/g
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/bottomWater/g
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       uniformDimensionedVectorField;
-    location    "constant";
     object      g;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -18,5 +17,4 @@ FoamFile
 dimensions      [0 1 -2 0 0 0 0];
 value           (0 -9.81 0);
 
-
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/bottomWater/radiationProperties b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/bottomWater/radiationProperties
new file mode 100644
index 0000000000..42defceaf4
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/bottomWater/radiationProperties
@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      radiationProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+radiation       off;
+
+radiationModel  none;
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/0/epsilon b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/bottomWater/thermophysicalProperties
similarity index 59%
rename from tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/0/epsilon
rename to tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/bottomWater/thermophysicalProperties
index f1ac119104..828c78d208 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/0/epsilon
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/bottomWater/thermophysicalProperties
@@ -9,39 +9,43 @@ FoamFile
 {
     version     2.0;
     format      ascii;
-    class       volScalarField;
-    location    "0";
-    object      epsilon;
+    class       dictionary;
+    object      thermophysicalProperties;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [0 2 -3 0 0 0 0];
-
-internalField   uniform 4e-06;
-
-boundaryField
+thermoType
 {
-    frontAndBack
+    type            heRhoThermo;
+    mixture         pureMixture;
+    transport       const;
+    thermo          hConst;
+    equationOfState rhoConst;
+    specie          specie;
+    energy          sensibleEnthalpy;
+}
+
+mixture
+{
+    specie
     {
-        type            epsilonWallFunction;
-        value           uniform 4e-06;
+        nMoles          1;
+        molWeight       18;
     }
-    topAndBottom
+    equationOfState
     {
-        type            epsilonWallFunction;
-        value           uniform 4e-06;
+        rho             1000;
     }
-    hot
+    thermodynamics
     {
-        type            epsilonWallFunction;
-        value           uniform 4e-06;
+        Cp              4181;
+        Hf              0;
     }
-    cold
+    transport
     {
-        type            epsilonWallFunction;
-        value           uniform 4e-06;
+        mu              959e-6;
+        Pr              6.62;
     }
 }
 
-
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/bottomWater/turbulenceProperties b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/bottomWater/turbulenceProperties
new file mode 100644
index 0000000000..e63bbc5081
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/bottomWater/turbulenceProperties
@@ -0,0 +1,19 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  laminar;
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/heater/radiationProperties b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/heater/radiationProperties
new file mode 100644
index 0000000000..ef3aba9e4b
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/heater/radiationProperties
@@ -0,0 +1,23 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      radiationProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+radiation off;
+
+radiationModel  none;
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/0/nut b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/heater/thermophysicalProperties
similarity index 62%
rename from tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/0/nut
rename to tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/heater/thermophysicalProperties
index bc7b910408..edb01db8b8 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/0/nut
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/heater/thermophysicalProperties
@@ -9,39 +9,45 @@ FoamFile
 {
     version     2.0;
     format      ascii;
-    class       volScalarField;
-    location    "0";
-    object      nut;
+    class       dictionary;
+    object      thermophysicalProperties;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [0 2 -1 0 0 0 0];
-
-internalField   uniform 0;
-
-boundaryField
+thermoType
 {
-    frontAndBack
+    type            heSolidThermo;
+    mixture         pureMixture;
+    transport       constIso;
+    thermo          hConst;
+    equationOfState rhoConst;
+    specie          specie;
+    energy          sensibleEnthalpy;
+}
+
+mixture
+{
+    specie
     {
-        type            nutUWallFunction;
-        value           uniform 0;
+        nMoles      1;
+        molWeight   50;
     }
-    topAndBottom
+
+    transport
     {
-        type            nutUWallFunction;
-        value           uniform 0;
+        kappa   80;
     }
-    hot
+
+    thermodynamics
     {
-        type            nutUWallFunction;
-        value           uniform 0;
+        Hf      0;
+        Cp      450;
     }
-    cold
+
+    equationOfState
     {
-        type            nutUWallFunction;
-        value           uniform 0;
+        rho     8000;
     }
 }
 
-
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/leftSolid/radiationProperties b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/leftSolid/radiationProperties
new file mode 120000
index 0000000000..08087c37b4
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/leftSolid/radiationProperties
@@ -0,0 +1 @@
+../heater/radiationProperties
\ No newline at end of file
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/leftSolid/thermophysicalProperties b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/leftSolid/thermophysicalProperties
new file mode 120000
index 0000000000..dc4d3a18ee
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/leftSolid/thermophysicalProperties
@@ -0,0 +1 @@
+../heater/thermophysicalProperties
\ No newline at end of file
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/system/blockMeshDict b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/polyMesh/blockMeshDict
similarity index 63%
rename from tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/system/blockMeshDict
rename to tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/polyMesh/blockMeshDict
index 5f2f1e8fd2..2a69ccd09c 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/system/blockMeshDict
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/polyMesh/blockMeshDict
@@ -14,66 +14,77 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-convertToMeters 0.001;
+convertToMeters 1;
 
 vertices
 (
-    ( 0     0  -260)
-    (76     0  -260)
-    (76  2180  -260)
-    ( 0  2180  -260)
-    ( 0     0   260)
-    (76     0   260)
-    (76  2180   260)
-    ( 0  2180   260)
+    (-0.1 -0.04  -0.05)
+    ( 0.1 -0.04  -0.05)
+    ( 0.1  0.04  -0.05)
+    (-0.1  0.04  -0.05)
+    (-0.1 -0.04   0.05)
+    ( 0.1 -0.04   0.05)
+    ( 0.1  0.04   0.05)
+    (-0.1  0.04   0.05)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (30 10 10) simpleGrading (1 1 1)
 );
 
 edges
 (
 );
 
-blocks
-(
-    hex (0 1 2 3 4 5 6 7) (35 150 15) simpleGrading (1 1 1)
-);
-
 boundary
 (
-    frontAndBack
+    maxY
     {
         type wall;
         faces
         (
-            (0 1 5 4)
-            (2 3 7 6)
+            (3 7 6 2)
         );
     }
-
-    topAndBottom
+    minX
+    {
+        type patch;
+        faces
+        (
+            (0 4 7 3)
+        );
+    }
+    maxX
+    {
+        type patch;
+        faces
+        (
+            (2 6 5 1)
+        );
+    }
+    minY
+    {
+        type wall;
+        faces
+        (
+            (1 5 4 0)
+        );
+    }
+    minZ
+    {
+        type wall;
+        faces
+        (
+            (0 3 2 1)
+        );
+    }
+    maxZ
     {
         type wall;
         faces
         (
             (4 5 6 7)
-            (3 2 1 0)
-        );
-    }
-
-    hot
-    {
-        type wall;
-        faces
-        (
-            (6 5 1 2)
-        );
-    }
-
-    cold
-    {
-        type wall;
-        faces
-        (
-            (4 7 3 0)
         );
     }
 );
@@ -81,3 +92,5 @@ boundary
 mergePatchPairs
 (
 );
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/constant/polyMesh/boundary b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/polyMesh/boundary
similarity index 65%
rename from tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/constant/polyMesh/boundary
rename to tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/polyMesh/boundary
index 0e07c61c65..b692f177a0 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/constant/polyMesh/boundary
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/polyMesh/boundary
@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
+|  \\    /   O peration     | Version:  dev-OpenCFD.feature-externalCoupled   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
@@ -15,35 +15,47 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-4
+6
 (
-    frontAndBack
+    maxY
     {
         type            wall;
         inGroups        1(wall);
-        nFaces          1050;
-        startFace       228225;
+        nFaces          300;
+        startFace       8300;
     }
-    topAndBottom
+    minX
     {
-        type            wall;
-        inGroups        1(wall);
-        nFaces          10500;
-        startFace       229275;
+        type            patch;
+        nFaces          100;
+        startFace       8600;
     }
-    hot
+    maxX
     {
-        type            wall;
-        inGroups        1(wall);
-        nFaces          2250;
-        startFace       239775;
+        type            patch;
+        nFaces          100;
+        startFace       8700;
     }
-    cold
+    minY
     {
         type            wall;
         inGroups        1(wall);
-        nFaces          2250;
-        startFace       242025;
+        nFaces          300;
+        startFace       8800;
+    }
+    minZ
+    {
+        type            wall;
+        inGroups        1(wall);
+        nFaces          300;
+        startFace       9100;
+    }
+    maxZ
+    {
+        type            wall;
+        inGroups        1(wall);
+        nFaces          300;
+        startFace       9400;
     }
 )
 
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/regionProperties b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/regionProperties
new file mode 100644
index 0000000000..0666069c2c
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/regionProperties
@@ -0,0 +1,24 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      regionProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+regions
+(
+    fluid       (bottomWater topAir)
+    solid       (heater leftSolid rightSolid)
+);
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/rightSolid/radiationProperties b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/rightSolid/radiationProperties
new file mode 120000
index 0000000000..08087c37b4
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/rightSolid/radiationProperties
@@ -0,0 +1 @@
+../heater/radiationProperties
\ No newline at end of file
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/rightSolid/thermophysicalProperties b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/rightSolid/thermophysicalProperties
new file mode 120000
index 0000000000..dc4d3a18ee
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/rightSolid/thermophysicalProperties
@@ -0,0 +1 @@
+../heater/thermophysicalProperties
\ No newline at end of file
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/topAir/g b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/topAir/g
new file mode 120000
index 0000000000..fea3757006
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/topAir/g
@@ -0,0 +1 @@
+../bottomWater/g
\ No newline at end of file
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/topAir/radiationProperties b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/topAir/radiationProperties
new file mode 120000
index 0000000000..583b06cb0b
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/topAir/radiationProperties
@@ -0,0 +1 @@
+../bottomWater/radiationProperties
\ No newline at end of file
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/constant/thermophysicalProperties b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/topAir/thermophysicalProperties
similarity index 88%
rename from tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/constant/thermophysicalProperties
rename to tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/topAir/thermophysicalProperties
index 95579e34cb..551767d574 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/constant/thermophysicalProperties
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/topAir/thermophysicalProperties
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "constant";
+    location    "constant/bottomWater";
     object      thermophysicalProperties;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -31,17 +31,17 @@ mixture
     specie
     {
         nMoles          1;
-        molWeight       28.96;
+        molWeight       28.9;
     }
     thermodynamics
     {
-        Cp              1004.4;
+        Cp              1000;
         Hf              0;
     }
     transport
     {
-        mu              1.831e-05;
-        Pr              0.705;
+        mu              1.8e-05;
+        Pr              0.7;
     }
 }
 
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/topAir/turbulenceProperties b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/topAir/turbulenceProperties
new file mode 120000
index 0000000000..ec52cbd592
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/constant/topAir/turbulenceProperties
@@ -0,0 +1 @@
+../bottomWater/turbulenceProperties
\ No newline at end of file
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/externalSolver b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/externalSolver
similarity index 53%
rename from tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/externalSolver
rename to tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/externalSolver
index 03c25c644a..cf8caf9c6b 100755
--- a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/externalSolver
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/externalSolver
@@ -1,20 +1,29 @@
 #!/bin/sh
 #
 # Dummy external solver to communicate with OpenFOAM via externalCoupled
-# boundary conditions
+# functionObject
 #
 # Functionality is hard-coded for this particular test case
 # - patch temperatures increased by 1K on each step
 #
-cd ${0%/*} || exit 1    # run from this directory
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Check for unassigned variables
+set -u
 
 echo "Executing dummy external solver"
 
 commsDir="comms"
+regionGroupName="heater_topAir"
+patchGroupName="coupleGroup"
+fieldName="T"
+
 lockFile="${commsDir}/OpenFOAM.lock"
-dataFile="${commsDir}/data"
+dataFile="${commsDir}/${regionGroupName}/${patchGroupName}/${fieldName}"
 waitSec=1
 timeOut=10
+nSteps=200  # maximum number of time steps. Note: should be more than
+            # number of iterations on the OpenFOAM side
 refGrad=0
 valueFraction=1
 
@@ -27,16 +36,21 @@ init()
 {
     log "initialisation: creating ${dataFile}.in"
 
-    # Hard-coded for 2 patches of size 2250
-    n=2250
-    refCold=283
-    refHot=303
+    # Hard-coded for patch of size 8 (heater/minY)
+    n1=8
+    refValue1=500
     touch "${dataFile}.in"
-    for i in $(seq 1 $n); do
-        echo "$refHot $refGrad $valueFraction" >> "${dataFile}.in"
+    log "initialisation: adding $n1 data elements with refValue $refValue1"
+    for i in $(seq 1 $n1); do
+        echo "$refValue1 $refGrad $valueFraction" >> "${dataFile}.in"
     done
-    for i in $(seq 1 $n); do
-        echo "$refCold $refGrad $valueFraction" >> "${dataFile}.in"
+
+    # Hard-coded for patch of size 40 (topAir/minX)
+    n2=40
+    refValue2=300
+    log "initialisation: adding $n2 data elements with refValue $refValue2"
+    for i in $(seq 1 $n2); do
+        echo "$refValue2 $refGrad $valueFraction" >> "${dataFile}.in"
     done
 
     # create lock file to pass control to OF
@@ -44,6 +58,10 @@ init()
 }
 
 
+# create the comms directory
+mkdir -p ${commsDir}/${regionGroupName}/${patchGroupName}
+
+
 # tutorial case employs the 'initByExternalOption', so we need to provide
 # the initial values
 init
@@ -51,7 +69,7 @@ init
 
 totalWait=0
 step=0
-while [ 1 ]; do
+while [ $step -lt $nSteps ]; do
     if [ -f $lockFile ]; then
         log "found lock file ${lockFile} - waiting"
         totalWait=$(expr $totalWait + $waitSec)
@@ -70,9 +88,10 @@ while [ 1 ]; do
         log "sleeping for $waitSec secs to simulate external process"
         sleep $waitSec
 
-        log "creating ${dataFile}.in"
+        log "updating ${dataFile}.in from ${dataFile}.out"
 
-        awk '{if( $1 != "#" ){print $2+1 " 0 1"}}' ${dataFile}.out > ${dataFile}.in
+        awk '{if( $1 != "#" ){print $1+1 " 0 1"}}' \
+            ${dataFile}.out | tee ${dataFile}.in
 
         log "creating lock file ${lockFile}"
         touch ${lockFile}
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/README b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/README
new file mode 100644
index 0000000000..5a81b9a570
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/README
@@ -0,0 +1,3 @@
+fvSolution is used for outer correctors specification.
+fvSchemes is only so that pre-processing activities can proceed
+
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/bottomWater/changeDictionaryDict b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/bottomWater/changeDictionaryDict
new file mode 100644
index 0000000000..63a08c87a1
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/bottomWater/changeDictionaryDict
@@ -0,0 +1,173 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      changeDictionaryDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dictionaryReplacement
+{
+    U
+    {
+        internalField   uniform (0.001 0 0);
+
+        boundaryField
+        {
+            minX
+            {
+                type            fixedValue;
+                value           uniform (0.001 0 0);
+            }
+
+            maxX
+            {
+                type            inletOutlet;
+                inletValue      uniform (0 0 0);
+            }
+
+            ".*"
+            {
+                type            fixedValue;
+                value           uniform (0 0 0);
+            }
+        }
+    }
+
+    T
+    {
+        internalField   uniform 300;
+
+        boundaryField
+        {
+            minX
+            {
+                type            fixedValue;
+                value           uniform 300;
+            }
+
+            maxX
+            {
+                type            inletOutlet;
+                inletValue      uniform 300;
+            }
+
+            ".*"
+            {
+                type            zeroGradient;
+                value           uniform 300;
+            }
+
+            "bottomWater_to_.*"
+            {
+                type            compressible::turbulentTemperatureCoupledBaffleMixed;
+                Tnbr            T;
+                kappa           fluidThermo;
+                kappaName       none;
+                value           uniform 300;
+            }
+        }
+    }
+
+    epsilon
+    {
+        internalField   uniform 0.01;
+
+        boundaryField
+        {
+            minX
+            {
+                type            fixedValue;
+                value           uniform 0.01;
+            }
+
+            maxX
+            {
+                type            inletOutlet;
+                inletValue      uniform 0.01;
+            }
+
+            ".*"
+            {
+                type            epsilonWallFunction;
+                value           uniform 0.01;
+            }
+        }
+    }
+
+    k
+    {
+        internalField   uniform 0.1;
+
+        boundaryField
+        {
+            minX
+            {
+                type            inletOutlet;
+                inletValue      uniform 0.1;
+            }
+
+            maxX
+            {
+                type            zeroGradient;
+                value           uniform 0.1;
+            }
+
+            ".*"
+            {
+                type            kqRWallFunction;
+                value           uniform 0.1;
+            }
+        }
+    }
+
+    p_rgh
+    {
+        internalField   uniform 0;
+
+        boundaryField
+        {
+            minX
+            {
+                type            zeroGradient;
+                value           uniform 0;
+            }
+
+            maxX
+            {
+                type            fixedValue;
+                value           uniform 0;
+            }
+
+            ".*"
+            {
+                type            fixedFluxPressure;
+                value           uniform 0;
+            }
+        }
+    }
+
+    p
+    {
+        internalField   uniform 0;
+
+        boundaryField
+        {
+            ".*"
+            {
+                type            calculated;
+                value           uniform 0;
+            }
+        }
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/bottomWater/decomposeParDict b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/bottomWater/decomposeParDict
new file mode 100644
index 0000000000..9bdbb31cd2
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/bottomWater/decomposeParDict
@@ -0,0 +1,72 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    note        "mesh decomposition control dictionary";
+    location    "system";
+    object      decomposeParDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains  4;
+
+//- Keep owner and neighbour on same processor for faces in zones:
+// preserveFaceZones (heater solid1 solid3);
+
+method          scotch;
+// method          hierarchical;
+// method          simple;
+// method          manual;
+
+simpleCoeffs
+{
+    n           (2 2 1);
+    delta       0.001;
+}
+
+hierarchicalCoeffs
+{
+    n           (2 2 1);
+    delta       0.001;
+    order       xyz;
+}
+
+scotchCoeffs
+{
+    //processorWeights
+    //(
+    //    1
+    //    1
+    //    1
+    //    1
+    //);
+    //writeGraph  true;
+    //strategy "b";
+}
+
+manualCoeffs
+{
+    dataFile    "decompositionData";
+}
+
+
+//// Is the case distributed
+//distributed     yes;
+//// Per slave (so nProcs-1 entries) the directory above the case.
+//roots
+//(
+//    "/tmp"
+//    "/tmp"
+//);
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/system/fvSchemes b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/bottomWater/fvSchemes
similarity index 70%
rename from tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/system/fvSchemes
rename to tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/bottomWater/fvSchemes
index 0b96d46074..1278f2856a 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/system/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/bottomWater/fvSchemes
@@ -16,7 +16,7 @@ FoamFile
 
 ddtSchemes
 {
-    default steadyState;
+    default Euler;
 }
 
 gradSchemes
@@ -28,18 +28,19 @@ divSchemes
 {
     default         none;
 
-    div(phi,U)      bounded Gauss limitedLinear 0.2;
-    div(phi,K)      bounded Gauss limitedLinear 0.2;
-    div(phi,h)      bounded Gauss limitedLinear 0.2;
-    div(phi,k)      bounded Gauss limitedLinear 0.2;
-    div(phi,epsilon) bounded Gauss limitedLinear 0.2;
-    div(phi,omega) bounded Gauss limitedLinear 0.2;
+    div(phi,U)      Gauss upwind;
+    div(phi,K)      Gauss linear;
+    div(phi,h)      Gauss upwind;
+    div(phi,k)      Gauss upwind;
+    div(phi,epsilon) Gauss upwind;
+    div(phi,R)      Gauss upwind;
+    div(R)          Gauss linear;
     div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
 {
-    default         Gauss linear orthogonal;
+    default        Gauss linear corrected;
 }
 
 interpolationSchemes
@@ -49,13 +50,13 @@ interpolationSchemes
 
 snGradSchemes
 {
-    default         orthogonal;
+    default         corrected;
 }
 
-wallDist
+fluxRequired
 {
-    method meshWave;
+    default         no;
+    p_rgh;
 }
 
-
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/system/fvSolution b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/bottomWater/fvSolution
similarity index 62%
rename from tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/system/fvSolution
rename to tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/bottomWater/fvSolution
index d30cdf20b4..102345b38d 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/system/fvSolution
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/bottomWater/fvSolution
@@ -10,68 +10,78 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "system";
     object      fvSolution;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 solvers
 {
+    rho
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-7;
+        relTol          0.1;
+    }
+
+    rhoFinal
+    {
+        $rho;
+        tolerance       1e-7;
+        relTol          0;
+    }
+
     p_rgh
     {
         solver           GAMG;
         tolerance        1e-7;
         relTol           0.01;
 
-        smoother         DICGaussSeidel;
+        smoother         GaussSeidel;
 
-        cacheAgglomeration true;
+        cacheAgglomeration  true;
         nCellsInCoarsestLevel 10;
         agglomerator     faceAreaPair;
         mergeLevels      1;
     }
 
-    "(U|h|k|epsilon|omega)"
+    p_rghFinal
     {
-        solver          PBiCG;
-        preconditioner  DILU;
-        tolerance       1e-8;
-        relTol          0.1;
+        $p_rgh;
+        tolerance        1e-7;
+        relTol           0;
+    }
+
+    "(U|h|k|epsilon|R)"
+    {
+        solver           PBiCG;
+        preconditioner   DILU;
+        tolerance        1e-7;
+        relTol           0.1;
+    }
+
+    "(U|h|k|epsilon|R)Final"
+    {
+        $U;
+        tolerance        1e-7;
+        relTol           0;
     }
 }
 
-SIMPLE
+PIMPLE
 {
-    momentumPredictor yes;
+    momentumPredictor   on;
+    nCorrectors         2;
     nNonOrthogonalCorrectors 0;
-    pRefCell        0;
-    pRefValue       0;
-
-    residualControl
-    {
-        p_rgh           1e-2;
-        U               1e-3;
-        h               1e-3;
-
-        // possibly check turbulence fields
-        "(k|epsilon|omega)" 1e-3;
-    }
 }
 
 relaxationFactors
 {
-    fields
-    {
-        rho             1.0;
-        p_rgh           0.7;
-    }
     equations
     {
-        U               0.3;
-        h               0.3;
-        "(k|epsilon|omega)" 0.7;
+        "h.*"           1;
+        "U.*"           1;
     }
 }
 
-
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/controlDict b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/controlDict
new file mode 100644
index 0000000000..817b7f284a
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/controlDict
@@ -0,0 +1,96 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Library defines new boundary conditions
+libs            ("libOpenFOAM.so" "libjobControl.so");
+
+application     chtMultiRegionFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         1;
+
+deltaT          0.001;
+
+writeControl    adjustableRunTime;
+writeInterval   0.1;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  8;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+maxCo           0.6;
+
+// Maximum diffusion number
+maxDi           10.0;
+
+adjustTimeStep  yes;
+
+functions
+{
+    externalCoupled
+    {
+        // Where to load it from (if not already in solver)
+        functionObjectLibs ("libjobControl.so");
+
+        type            externalCoupled;
+
+        // Directory to use for communication
+        commsDir        "${FOAM_CASE}/comms";
+
+        // Does external process start first
+        initByExternal  true;
+
+        // Additional output
+        log             true;
+
+        regions
+        {
+            // Region name (wildcards allowed)
+            "(topAir|heater)"
+            {
+                // In topAir adjust the minX patch (fixedValue)
+
+                // Patch or patchGroup
+                coupleGroup
+                {
+                    // Fields to output in commsDir
+                    writeFields (T);
+                    // Fields to read from commsDir
+                    readFields  (T);
+                }
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/decomposeParDict b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/decomposeParDict
new file mode 100644
index 0000000000..3eb54ba02f
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/decomposeParDict
@@ -0,0 +1,60 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    note        "mesh decomposition control dictionary";
+    location    "system";
+    object      decomposeParDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains  4;
+
+//- Keep owner and neighbour on same processor for faces in zones:
+// preserveFaceZones (heater solid1 solid3);
+
+method          scotch;
+// method          hierarchical;
+// method          simple;
+// method          manual;
+
+simpleCoeffs
+{
+    n           (2 2 1);
+    delta       0.001;
+}
+
+hierarchicalCoeffs
+{
+    n           (2 2 1);
+    delta       0.001;
+    order       xyz;
+}
+
+
+manualCoeffs
+{
+    dataFile    "decompositionData";
+}
+
+
+//// Is the case distributed
+//distributed     yes;
+//// Per slave (so nProcs-1 entries) the directory above the case.
+//roots
+//(
+//    "/tmp"
+//    "/tmp"
+//);
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/system/decomposeParDict b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/fvSchemes
similarity index 82%
rename from tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/system/decomposeParDict
rename to tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/fvSchemes
index 70c6c85c26..e7d321e959 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/system/decomposeParDict
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/fvSchemes
@@ -10,19 +10,36 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    location    "system";
-    object      decomposeParDict;
+    object      fvSchemes;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-numberOfSubdomains 4;
-
-method          simple;
-
-simpleCoeffs
+ddtSchemes
+{
+}
+
+gradSchemes
+{
+}
+
+divSchemes
+{
+}
+
+laplacianSchemes
+{
+}
+
+interpolationSchemes
+{
+}
+
+snGradSchemes
+{
+}
+
+fluxRequired
 {
-    n               (2 2 1);
-    delta           0.001;
 }
 
 
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/constant/turbulenceProperties b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/fvSolution
similarity index 84%
rename from tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/constant/turbulenceProperties
rename to tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/fvSolution
index 691c02bb16..cde9e49189 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/constant/turbulenceProperties
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/fvSolution
@@ -10,20 +10,13 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    object      RASProperties;
+    object      fvSolution;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-simulationType RAS;
-
-RAS
+PIMPLE
 {
-    RASModel            kOmegaSST;
-
-    turbulence          on;
-
-    printCoeffs         on;
+    nOuterCorrectors 1;
 }
 
-
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/heater/changeDictionaryDict b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/heater/changeDictionaryDict
new file mode 100644
index 0000000000..d8d26b8376
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/heater/changeDictionaryDict
@@ -0,0 +1,77 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      changeDictionaryDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dictionaryReplacement
+{
+    boundary
+    {
+        minY
+        {
+            type            patch;
+            inGroups        (coupleGroup);
+        }
+        minZ
+        {
+            type            patch;
+        }
+        maxZ
+        {
+            type            patch;
+        }
+    }
+
+
+    T
+    {
+        internalField   uniform 300;
+
+        boundaryField
+        {
+            ".*"
+            {
+                type            zeroGradient;
+                value           uniform 300;
+            }
+            "heater_to_.*"
+            {
+                type            compressible::turbulentTemperatureCoupledBaffleMixed;
+                Tnbr            T;
+                kappa           solidThermo;
+                kappaName       none;
+                value           uniform 300;
+            }
+
+            heater_to_leftSolid
+            {
+                type            compressible::turbulentTemperatureCoupledBaffleMixed;
+                Tnbr            T;
+                kappa           solidThermo;
+                kappaName       none;
+                thicknessLayers (1e-3);
+                kappaLayers     (5e-4);
+                value           uniform 300;
+            }
+
+            minY
+            {
+                type            fixedValue;
+                value           uniform 500;
+            }
+        }
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/system/controlDict b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/heater/decomposeParDict
similarity index 68%
rename from tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/system/controlDict
rename to tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/heater/decomposeParDict
index 2df20103c1..e405fc009e 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/externalCoupledCavity/system/controlDict
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/heater/decomposeParDict
@@ -10,39 +10,35 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    object      controlDict;
+    location    "system";
+    object      decomposeParDict;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-application     buoyantSimpleFoam;
+numberOfSubdomains  4;
 
-startFrom       startTime;
+method          scotch;
 
-startTime       0;
+simpleCoeffs
+{
+    n           (2 2 1);
+    delta       0.001;
+}
 
-stopAt          endTime;
+hierarchicalCoeffs
+{
+    n           (2 2 1);
+    delta       0.001;
+    order       xyz;
+}
 
-endTime         100;
-
-deltaT          1;
-
-writeControl    timeStep;
-
-writeInterval   10;
-
-purgeWrite      0;
-
-writeFormat     ascii;
-
-writePrecision  6;
-
-writeCompression off;
-
-timeFormat      general;
-
-timePrecision   6;
-
-runTimeModifiable true;
+scotchCoeffs
+{
+}
 
+manualCoeffs
+{
+    dataFile    "decompositionData";
+}
 
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/heater/fvSchemes b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/heater/fvSchemes
new file mode 100644
index 0000000000..d1227f9fad
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/heater/fvSchemes
@@ -0,0 +1,53 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+}
+
+laplacianSchemes
+{
+    default             none;
+    laplacian(alpha,h)  Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+fluxRequired
+{
+    default         no;
+}
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/heater/fvSolution b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/heater/fvSolution
new file mode 100644
index 0000000000..ccb0d3183f
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/heater/fvSolution
@@ -0,0 +1,40 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    h
+    {
+        solver           PCG;
+        preconditioner   DIC;
+        tolerance        1e-06;
+        relTol           0.1;
+    }
+
+    hFinal
+    {
+        $h;
+        tolerance        1e-06;
+        relTol           0;
+    }
+}
+
+PIMPLE
+{
+    nNonOrthogonalCorrectors 0;
+}
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/leftSolid/changeDictionaryDict b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/leftSolid/changeDictionaryDict
new file mode 100644
index 0000000000..1709662b51
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/leftSolid/changeDictionaryDict
@@ -0,0 +1,65 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      changeDictionaryDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dictionaryReplacement
+{
+    boundary
+    {
+        minZ
+        {
+            type            patch;
+        }
+        maxZ
+        {
+            type            patch;
+        }
+    }
+
+    T
+    {
+        internalField   uniform 300;
+
+        boundaryField
+        {
+            ".*"
+            {
+                type            zeroGradient;
+                value           uniform 300;
+            }
+            "leftSolid_to_.*"
+            {
+                type            compressible::turbulentTemperatureCoupledBaffleMixed;
+                Tnbr            T;
+                kappa           solidThermo;
+                kappaName       none;
+                value           uniform 300;
+            }
+
+            leftSolid_to_heater
+            {
+                type            compressible::turbulentTemperatureCoupledBaffleMixed;
+                Tnbr            T;
+                kappa           solidThermo;
+                kappaName       none;
+                thicknessLayers (1e-3);
+                kappaLayers     (5e-4);
+                value           uniform 300;
+            }
+        }
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/leftSolid/decomposeParDict b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/leftSolid/decomposeParDict
new file mode 100644
index 0000000000..e405fc009e
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/leftSolid/decomposeParDict
@@ -0,0 +1,44 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains  4;
+
+method          scotch;
+
+simpleCoeffs
+{
+    n           (2 2 1);
+    delta       0.001;
+}
+
+hierarchicalCoeffs
+{
+    n           (2 2 1);
+    delta       0.001;
+    order       xyz;
+}
+
+scotchCoeffs
+{
+}
+
+manualCoeffs
+{
+    dataFile    "decompositionData";
+}
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/leftSolid/fvSchemes b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/leftSolid/fvSchemes
new file mode 120000
index 0000000000..63236f302c
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/leftSolid/fvSchemes
@@ -0,0 +1 @@
+../heater/fvSchemes
\ No newline at end of file
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/leftSolid/fvSolution b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/leftSolid/fvSolution
new file mode 120000
index 0000000000..0bde0fc62f
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/leftSolid/fvSolution
@@ -0,0 +1 @@
+../heater/fvSolution
\ No newline at end of file
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/rightSolid/changeDictionaryDict b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/rightSolid/changeDictionaryDict
new file mode 100644
index 0000000000..ba99b13705
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/rightSolid/changeDictionaryDict
@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      changeDictionaryDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dictionaryReplacement
+{
+    boundary
+    {
+        minZ
+        {
+            type            patch;
+        }
+        maxZ
+        {
+            type            patch;
+        }
+    }
+
+    T
+    {
+        internalField   uniform 300;
+
+        boundaryField
+        {
+            ".*"
+            {
+                type            zeroGradient;
+                value           uniform 300;
+            }
+            "rightSolid_to_.*"
+            {
+                type            compressible::turbulentTemperatureCoupledBaffleMixed;
+                Tnbr            T;
+                kappa           solidThermo;
+                kappaName       none;
+                value           uniform 300;
+            }
+        }
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/rightSolid/decomposeParDict b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/rightSolid/decomposeParDict
new file mode 100644
index 0000000000..e405fc009e
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/rightSolid/decomposeParDict
@@ -0,0 +1,44 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains  4;
+
+method          scotch;
+
+simpleCoeffs
+{
+    n           (2 2 1);
+    delta       0.001;
+}
+
+hierarchicalCoeffs
+{
+    n           (2 2 1);
+    delta       0.001;
+    order       xyz;
+}
+
+scotchCoeffs
+{
+}
+
+manualCoeffs
+{
+    dataFile    "decompositionData";
+}
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/rightSolid/fvSchemes b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/rightSolid/fvSchemes
new file mode 120000
index 0000000000..63236f302c
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/rightSolid/fvSchemes
@@ -0,0 +1 @@
+../heater/fvSchemes
\ No newline at end of file
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/rightSolid/fvSolution b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/rightSolid/fvSolution
new file mode 120000
index 0000000000..0bde0fc62f
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/rightSolid/fvSolution
@@ -0,0 +1 @@
+../heater/fvSolution
\ No newline at end of file
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/topAir/changeDictionaryDict b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/topAir/changeDictionaryDict
new file mode 100644
index 0000000000..d31389b049
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/topAir/changeDictionaryDict
@@ -0,0 +1,179 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      changeDictionaryDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dictionaryReplacement
+{
+    boundary
+    {
+        minX
+        {
+            inGroups        (coupleGroup);
+        }
+    }
+
+    U
+    {
+        internalField   uniform (0.1 0 0);
+
+        boundaryField
+        {
+            ".*"
+            {
+                type            fixedValue;
+                value           uniform (0 0 0);
+            }
+            minX
+            {
+                type            fixedValue;
+                value           uniform ( 0.1 0 0 );
+            }
+            maxX
+            {
+                type            inletOutlet;
+                inletValue      uniform ( 0 0 0 );
+                value           uniform ( 0.1 0 0 );
+            }
+        }
+    }
+
+    T
+    {
+        internalField   uniform 300;
+
+        boundaryField
+        {
+            ".*"
+            {
+                type            zeroGradient;
+            }
+
+            minX
+            {
+                type            fixedValue;
+                value           uniform 300;
+            }
+            maxX
+            {
+                type            inletOutlet;
+                inletValue      uniform 300;
+                value           uniform 300;
+            }
+
+            "topAir_to_.*"
+            {
+                type            compressible::turbulentTemperatureCoupledBaffleMixed;
+                Tnbr            T;
+                kappa           fluidThermo;
+                kappaName       none;
+                value           uniform 300;
+            }
+        }
+    }
+
+    epsilon
+    {
+        internalField   uniform 0.01;
+
+        boundaryField
+        {
+            ".*"
+            {
+                type            epsilonWallFunction;
+                value           uniform 0.01;
+            }
+
+            minX
+            {
+                type            fixedValue;
+                value           uniform 0.01;
+            }
+            maxX
+            {
+                type            inletOutlet;
+                inletValue      uniform 0.01;
+                value           uniform 0.01;
+            }
+        }
+    }
+
+    k
+    {
+        internalField   uniform 0.1;
+
+        boundaryField
+        {
+            ".*"
+            {
+                type            kqRWallFunction;
+                value           uniform 0.1;
+            }
+
+            minX
+            {
+                type            fixedValue;
+                value           uniform 0.1;
+            }
+            maxX
+            {
+                type            inletOutlet;
+                inletValue      uniform 0.1;
+                value           uniform 0.1;
+            }
+        }
+    }
+
+    p_rgh
+    {
+        internalField   uniform 1e5;
+
+        boundaryField
+        {
+            ".*"
+            {
+                type            fixedFluxPressure;
+                value           uniform 1e5;
+            }
+
+            maxX
+            {
+                type            fixedValue;
+                value           uniform 1e5;
+            }
+        }
+    }
+
+    p
+    {
+        internalField   uniform 1e5;
+
+        boundaryField
+        {
+            ".*"
+            {
+                type            calculated;
+                value           uniform 1e5;
+            }
+
+            maxX
+            {
+                type            calculated;
+                value           uniform 1e5;
+            }
+        }
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/topAir/decomposeParDict b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/topAir/decomposeParDict
new file mode 100644
index 0000000000..9bdbb31cd2
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/topAir/decomposeParDict
@@ -0,0 +1,72 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    note        "mesh decomposition control dictionary";
+    location    "system";
+    object      decomposeParDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains  4;
+
+//- Keep owner and neighbour on same processor for faces in zones:
+// preserveFaceZones (heater solid1 solid3);
+
+method          scotch;
+// method          hierarchical;
+// method          simple;
+// method          manual;
+
+simpleCoeffs
+{
+    n           (2 2 1);
+    delta       0.001;
+}
+
+hierarchicalCoeffs
+{
+    n           (2 2 1);
+    delta       0.001;
+    order       xyz;
+}
+
+scotchCoeffs
+{
+    //processorWeights
+    //(
+    //    1
+    //    1
+    //    1
+    //    1
+    //);
+    //writeGraph  true;
+    //strategy "b";
+}
+
+manualCoeffs
+{
+    dataFile    "decompositionData";
+}
+
+
+//// Is the case distributed
+//distributed     yes;
+//// Per slave (so nProcs-1 entries) the directory above the case.
+//roots
+//(
+//    "/tmp"
+//    "/tmp"
+//);
+
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/topAir/fvSchemes b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/topAir/fvSchemes
new file mode 120000
index 0000000000..323c0787e2
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/topAir/fvSchemes
@@ -0,0 +1 @@
+../bottomWater/fvSchemes
\ No newline at end of file
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/topAir/fvSolution b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/topAir/fvSolution
new file mode 120000
index 0000000000..90d9c92347
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/topAir/fvSolution
@@ -0,0 +1 @@
+../bottomWater/fvSolution
\ No newline at end of file
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/topoSetDict b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/topoSetDict
new file mode 100644
index 0000000000..feb3d898bc
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/externalCoupledMultiRegionHeater/system/topoSetDict
@@ -0,0 +1,178 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      topoSetDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+actions
+(
+    // Heater
+    {
+        name    heaterCellSet;
+        type    cellSet;
+        action  new;
+        source  boxToCell;
+        sourceInfo
+        {
+            box (-0.01001    0 -100 )(0.01001 0.00999 100);
+        }
+    }
+    {
+        name    heaterCellSet;
+        type    cellSet;
+        action  add;
+        source  boxToCell;
+        sourceInfo
+        {
+            box (-0.01001 -100 -0.01001)(0.01001 0.00999 0.01001);
+        }
+    }
+    {
+        name    heater;
+        type    cellZoneSet;
+        action  new;
+        source  setToCellZone;
+        sourceInfo
+        {
+            set heaterCellSet;
+        }
+    }
+
+    // leftSolid
+    {
+        name    leftSolidCellSet;
+        type    cellSet;
+        action  new;
+        source  boxToCell;
+        sourceInfo
+        {
+            box (-100 0 -100 )(-0.01001 0.00999 100);
+        }
+    }
+    {
+        name    leftSolid;
+        type    cellZoneSet;
+        action  new;
+        source  setToCellZone;
+        sourceInfo
+        {
+            set leftSolidCellSet;
+        }
+    }
+
+    // rightSolid
+    {
+        name    rightSolidCellSet;
+        type    cellSet;
+        action  new;
+        source  boxToCell;
+        sourceInfo
+        {
+            box (0.01001 0 -100 )(100 0.00999 100);
+        }
+    }
+    {
+        name    rightSolid;
+        type    cellZoneSet;
+        action  new;
+        source  setToCellZone;
+        sourceInfo
+        {
+            set rightSolidCellSet;
+        }
+    }
+
+    // topAir
+    {
+        name    topAirCellSet;
+        type    cellSet;
+        action  new;
+        source  boxToCell;
+        sourceInfo
+        {
+            box (-100 0.00999 -100 )(100 100 100);
+        }
+    }
+    {
+        name    topAir;
+        type    cellZoneSet;
+        action  new;
+        source  setToCellZone;
+        sourceInfo
+        {
+            set topAirCellSet;
+        }
+    }
+
+
+    // bottomWater is all the other cells
+    {
+        name    bottomWaterCellSet;
+        type    cellSet;
+        action  new;
+        source  cellToCell;
+        sourceInfo
+        {
+            set heaterCellSet;
+        }
+    }
+    {
+        name    bottomWaterCellSet;
+        type    cellSet;
+        action  add;
+        source  cellToCell;
+        sourceInfo
+        {
+            set leftSolidCellSet;
+        }
+    }
+    {
+        name    bottomWaterCellSet;
+        type    cellSet;
+        action  add;
+        source  cellToCell;
+        sourceInfo
+        {
+            set rightSolidCellSet;
+        }
+    }
+    {
+        name    bottomWaterCellSet;
+        type    cellSet;
+        action  add;
+        source  cellToCell;
+        sourceInfo
+        {
+            set topAirCellSet;
+        }
+    }
+    {
+        name    bottomWaterCellSet;
+        type    cellSet;
+        action  invert;
+    }
+    {
+        name    bottomWater;
+        type    cellZoneSet;
+        action  new;
+        source  setToCellZone;
+        sourceInfo
+        {
+            set bottomWaterCellSet;
+        }
+    }
+);
+
+
+// ************************************************************************* //