diff --git a/applications/utilities/mesh/manipulation/createPatch/createPatch.C b/applications/utilities/mesh/manipulation/createPatch/createPatch.C
index 1f38a37854..0276c89bb3 100644
--- a/applications/utilities/mesh/manipulation/createPatch/createPatch.C
+++ b/applications/utilities/mesh/manipulation/createPatch/createPatch.C
@@ -98,11 +98,13 @@ void matchPatchFaces
     const labelList& patchesj,
 
     DynamicList<labelList>& interfaceMesh0,
+    DynamicList<label>& interfaceSource0,
     DynamicList<labelList>& interfacePatch0,
     DynamicList<wordList>& interfaceNames0,
     DynamicList<List<DynamicList<label>>>& interfaceFaces0,
 
     DynamicList<labelList>& interfaceMesh1,
+    DynamicList<label>& interfaceSource1,
     DynamicList<labelList>& interfacePatch1,
     DynamicList<wordList>& interfaceNames1,
     DynamicList<List<DynamicList<label>>>& interfaceFaces1
@@ -136,8 +138,10 @@ void matchPatchFaces
 
                 Info<< "Introducing interface " << inti << " between"
                     << " mesh " << meshes[meshi].name()
+                    //<< " source:" << sourcei
                     << " patch " << ppi.name()
                     << " and mesh " << meshes[meshj].name()
+                    //<< " source:" << sourcej
                     << " patch " << ppj.name()
                     << endl;
 
@@ -149,6 +153,8 @@ void matchPatchFaces
                 intMesh0.setSize(nSourcei, -1);
                 intMesh0[sourcei] = meshi;
 
+                interfaceSource0.append(sourcei);
+
                 interfacePatch0.append(labelList());
                 auto& intPatch0 = interfacePatch0.last();
                 intPatch0.setSize(nSourcei, -1);
@@ -157,15 +163,11 @@ void matchPatchFaces
                 interfaceNames0.append(wordList());
                 auto& intNames0 = interfaceNames0.last();
                 intNames0.setSize(nSourcei);
-                //intNames0[sourcei] =
-                //    meshes[meshi].name()
-                //  + "_to_"
-                //  + meshes[meshj].name();
                 intNames0[sourcei] =
                     patchName(entryName, meshes[meshi], meshes[meshj]);
 
 
-                // Mesh 0
+                // Mesh 1
                 //~~~~~~~
 
                 interfaceMesh1.append(labelList());
@@ -173,6 +175,8 @@ void matchPatchFaces
                 intMesh1.setSize(nSourcej, -1);
                 intMesh1[sourcej] = meshj;
 
+                interfaceSource1.append(sourcej);
+
                 interfacePatch1.append(labelList());
                 auto& intPatch1 = interfacePatch1.last();
                 intPatch1.setSize(nSourcej, -1);
@@ -181,10 +185,6 @@ void matchPatchFaces
                 interfaceNames1.append(wordList());
                 auto& intNames1 = interfaceNames1.last();
                 intNames1.setSize(nSourcej);
-                //intNames1[sourcej] =
-                //    meshes[meshj].name()
-                //  + "_to_"
-                //  + meshes[meshi].name();
                 intNames1[sourcej] =
                     patchName(entryName, meshes[meshj], meshes[meshi]);
 
@@ -283,11 +283,13 @@ void matchPatchFaces
     List<PtrList<dictionary>> patchInfoDicts,
 
     DynamicList<labelList>& interfaceMesh0,
+    DynamicList<label>& interfaceSource0,
     DynamicList<labelList>& interfacePatch0,
     DynamicList<List<DynamicList<label>>>& interfaceFaces0,
     DynamicList<wordList>& interfaceNames0,
 
     DynamicList<labelList>& interfaceMesh1,
+    DynamicList<label>& interfaceSource1,
     DynamicList<labelList>& interfacePatch1,
     DynamicList<List<DynamicList<label>>>& interfaceFaces1,
     DynamicList<wordList>& interfaceNames1
@@ -338,11 +340,13 @@ void matchPatchFaces
                         patchesj,
 
                         interfaceMesh0,
+                        interfaceSource0,
                         interfacePatch0,
                         interfaceNames0,
                         interfaceFaces0,
 
                         interfaceMesh1,
+                        interfaceSource1,
                         interfacePatch1,
                         interfaceNames1,
                         interfaceFaces1
@@ -413,10 +417,23 @@ void changePatchID
 
         if (!isRepatchedBoundary.set(facei-mesh.nInternalFaces()))
         {
-            FatalErrorInFunction
-                << "Face " << facei << " is already marked to be moved"
-                << " to patch " << meshMod.region()[facei]
-                << exit(FatalError);
+            static label nWarnings = 0;
+            if (nWarnings == 0)
+            {
+                const label newPatchi = meshMod.region()[facei];
+                //FatalErrorInFunction
+                WarningInFunction
+                    << "Face " << facei
+                    << " at " << mesh.faceCentres()[facei]
+                    << " marked for patch " << patchID
+                    << " name " << mesh.boundaryMesh()[patchID].name()
+                    << " is already marked for patch " << newPatchi
+                    << " name " << mesh.boundaryMesh()[newPatchi].name()
+                    << ". Suppressing further warnings"
+                    //<< exit(FatalError);
+                    << endl;
+            }
+            nWarnings++;
         }
 
         changePatchID(mesh, facei, patchID, meshMod);
@@ -721,7 +738,6 @@ int main(int argc, char *argv[])
     );
 
     #include "addOverwriteOption.H"
-    //#include "addRegionOption.H"
     #include "addAllRegionOptions.H"
 
     argList::addOption("dict", "file", "Alternative createPatchDict");
@@ -780,17 +796,25 @@ int main(int argc, char *argv[])
         forAll(patchSources, sourcei)
         {
             const auto& pDict = patchSources[sourcei];
-            patchNames[meshi][sourcei] = pDict.get<word>("name");
+            patchNames[meshi][sourcei] = pDict.getOrDefault<word>
+            (
+                "name",
+                word::null,
+                keyType::LITERAL
+            );
 
             patchInfoDicts[meshi].set
             (
                 sourcei,
                 new dictionary(pDict.subDict("patchInfo"))
             );
-            const dictionary& patchDict = patchInfoDicts[meshi][sourcei];
+            dictionary& patchDict = patchInfoDicts[meshi][sourcei];
             if (patchDict.found("AMIMethod"))
             {
                 matchMethods[meshi][sourcei] = patchDict.get<word>("AMIMethod");
+                // Disable full matching since we're trying to use AMIMethod to
+                // find out actual overlap
+                patchDict.add("requireMatch", false);
             }
 
             wordRes matchNames;
@@ -817,9 +841,11 @@ int main(int argc, char *argv[])
     // It matches all mesh against each other. Lower numbered mesh gets
     // postfix 0, higher numbered mesh postfix 1.
 
-    // Per interface, per mesh, per patchSource:
+    // Per interface, per patchSource:
     //      1. the lower numbered mesh
     DynamicList<labelList> interfaceMesh0;
+    //      1b. the source index (i.e. the patch dictionary)
+    DynamicList<label> interfaceSource0;
     //      2. the patch on the interfaceMesh0
     DynamicList<labelList> interfacePatch0;
     //      3. the facelabels on the interfaceMesh0
@@ -829,6 +855,7 @@ int main(int argc, char *argv[])
 
     // Same for the higher numbered mesh
     DynamicList<labelList> interfaceMesh1;
+    DynamicList<label> interfaceSource1;
     DynamicList<labelList> interfacePatch1;
     DynamicList<List<DynamicList<label>>> interfaceFaces1;
     DynamicList<wordList> interfaceNames1;
@@ -843,15 +870,17 @@ int main(int argc, char *argv[])
             interRegionSources,
             patchNames,
             matchPatchIDs,
-            matchMethods,   //faceAreaWeightAMI2D::typeName,
+            matchMethods,
             patchInfoDicts,
 
             interfaceMesh0,
+            interfaceSource0,
             interfacePatch0,
             interfaceFaces0,
             interfaceNames0,
 
             interfaceMesh1,
+            interfaceSource1,
             interfacePatch1,
             interfaceFaces1,
             interfaceNames1
@@ -859,6 +888,7 @@ int main(int argc, char *argv[])
     }
 
 
+
     // Read fields
     List<PtrList<volScalarField>> vsFlds(meshes.size());
     List<PtrList<volVectorField>> vvFlds(meshes.size());
@@ -914,6 +944,9 @@ int main(int argc, char *argv[])
     forAll(meshes, meshi)
     {
         fvMesh& mesh = meshes[meshi];
+
+        Info<< "\n\nAdding patches to mesh " << mesh.name() << nl << endl;
+
         const polyBoundaryMesh& patches = mesh.boundaryMesh();
         const dictionary& dict = dicts[meshi];
 
@@ -933,11 +966,20 @@ int main(int argc, char *argv[])
         {
             const dictionary& dict = patchSources[sourcei];
             const word sourceType(dict.get<word>("constructFrom"));
-            const word patchName(dict.get<word>("name"));
+            const word patchName
+            (
+                dict.getOrDefault<word>
+                (
+                    "name",
+                    word::null,
+                    keyType::LITERAL
+                )
+            );
 
-            dictionary patchDict(dict.subDict("patchInfo"));
+            dictionary patchDict(patchInfoDicts[meshi][sourcei]);
             patchDict.set("nFaces", 0);
-            patchDict.set("startFace", 0);  //startFacei);
+            patchDict.set("startFace", 0);    // Gets overwritten
+
 
             if (sourceType == "autoPatch")
             {
@@ -951,42 +993,45 @@ int main(int argc, char *argv[])
                     const labelList& allMeshes1 = interfaceMesh1[inti];
                     const wordList& allNames1 = interfaceNames1[inti];
 
-                    forAll(allMeshes0, sourcei)
+                    if
+                    (
+                        interfaceSource0[inti] == sourcei
+                     && allMeshes0[sourcei] == meshi
+                    )
                     {
-                        if (allMeshes0[sourcei] == meshi)
+                        // Current mesh is mesh0. mesh1 is the remote mesh.
+                        const label sourcej = interfaceSource1[inti];
+                        const word& patchName = allNames0[sourcei];
+                        if (patches.findPatchID(patchName) == -1)
                         {
-                            const auto& mesh1 = meshes[allMeshes1[sourcei]];
-                            const word& patchName = allNames0[sourcei];
-                            if (patches.findPatchID(patchName) == -1)
-                            {
-                                dictionary allDict(patchDict);
-                                allDict.set("sampleRegion", mesh1.name());
-                                const auto& destPatch = allNames1[sourcei];
-                                allDict.set("samplePatch", destPatch);
-                                allDict.set("neighbourPatch", destPatch);
+                            dictionary allDict(patchDict);
+                            const auto& mesh1 = meshes[allMeshes1[sourcej]];
+                            allDict.set("sampleRegion", mesh1.name());
+                            const auto& destPatch = allNames1[sourcej];
+                            allDict.set("samplePatch", destPatch);
+                            allDict.set("neighbourPatch", destPatch);
 
-                                Info<< "Adding new patch " << patchName
-                                    << " from " << allDict << endl;
+                            Info<< "Adding new patch " << patchName
+                                << " from " << allDict << endl;
 
-                                autoPtr<polyPatch> ppPtr
+                            autoPtr<polyPatch> ppPtr
+                            (
+                                polyPatch::New
                                 (
-                                    polyPatch::New
-                                    (
-                                        patchName,
-                                        allDict,
-                                        0,          // overwritten
-                                        patches
-                                    )
-                                );
-                                fvMeshTools::addPatch
-                                (
-                                    mesh,
-                                    ppPtr(),
-                                    patchDict.subOrEmptyDict("patchFields"),
-                                    calculatedFvPatchScalarField::typeName,
-                                    true
-                                );
-                            }
+                                    patchName,
+                                    allDict,
+                                    0,          // overwritten
+                                    patches
+                                )
+                            );
+                            fvMeshTools::addPatch
+                            (
+                                mesh,
+                                ppPtr(),
+                                patchDict.subOrEmptyDict("patchFields"),
+                                calculatedFvPatchScalarField::typeName,
+                                true
+                            );
                         }
                     }
                 }
@@ -998,42 +1043,46 @@ int main(int argc, char *argv[])
                     const labelList& allMeshes1 = interfaceMesh1[inti];
                     const wordList& allNames1 = interfaceNames1[inti];
 
-                    forAll(allMeshes1, sourcei)
+                    if
+                    (
+                        interfaceSource1[inti] == sourcei
+                     && allMeshes1[sourcei] == meshi
+                    )
                     {
-                        if (allMeshes1[sourcei] == meshi)
+                        // Current mesh is mesh1. mesh0 is the remote mesh.
+
+                        const label sourcej = interfaceSource0[inti];
+                        const word& patchName = allNames1[sourcei];
+                        if (patches.findPatchID(patchName) == -1)
                         {
-                            const auto& mesh0 = meshes[allMeshes0[sourcei]];
-                            const word& patchName = allNames1[sourcei];
-                            if (patches.findPatchID(patchName) == -1)
-                            {
-                                dictionary allDict(patchDict);
-                                const auto& destPatch = allNames0[sourcei];
-                                allDict.set("sampleRegion", mesh0.name());
-                                allDict.set("samplePatch", destPatch);
-                                allDict.set("neighbourPatch", destPatch);
+                            dictionary allDict(patchDict);
+                            const auto& destPatch = allNames0[sourcej];
+                            const auto& mesh0 = meshes[allMeshes0[sourcej]];
+                            allDict.set("sampleRegion", mesh0.name());
+                            allDict.set("samplePatch", destPatch);
+                            allDict.set("neighbourPatch", destPatch);
 
-                                Info<< "Adding new patch " << patchName
-                                    << " from " << allDict << endl;
+                            Info<< "Adding new patch " << patchName
+                                << " from " << allDict << endl;
 
-                                autoPtr<polyPatch> ppPtr
+                            autoPtr<polyPatch> ppPtr
+                            (
+                                polyPatch::New
                                 (
-                                    polyPatch::New
-                                    (
-                                        patchName,
-                                        allDict,
-                                        0,          // overwritten
-                                        patches
-                                    )
-                                );
-                                fvMeshTools::addPatch
-                                (
-                                    mesh,
-                                    ppPtr(),
-                                    patchDict.subOrEmptyDict("patchFields"),
-                                    calculatedFvPatchScalarField::typeName,
-                                    true
-                                );
-                            }
+                                    patchName,
+                                    allDict,
+                                    0,          // overwritten
+                                    patches
+                                )
+                            );
+                            fvMeshTools::addPatch
+                            (
+                                mesh,
+                                ppPtr(),
+                                patchDict.subOrEmptyDict("patchFields"),
+                                calculatedFvPatchScalarField::typeName,
+                                true
+                            );
                         }
                     }
                 }
@@ -1077,6 +1126,10 @@ int main(int argc, char *argv[])
     forAll(meshes, meshi)
     {
         fvMesh& mesh = meshes[meshi];
+
+        Info<< "\n\nRepatching mesh " << mesh.name() << nl << endl;
+
+
         const polyBoundaryMesh& patches = mesh.boundaryMesh();
         const dictionary& dict = dicts[meshi];
 
@@ -1096,7 +1149,15 @@ int main(int argc, char *argv[])
         forAll(patchSources, sourcei)
         {
             const dictionary& dict = patchSources[sourcei];
-            const word patchName(dict.get<word>("name"));
+            const word patchName
+            (
+                dict.getOrDefault<word>
+                (
+                    "name",
+                    word::null,
+                    keyType::LITERAL
+                )
+            );
             const word sourceType(dict.get<word>("constructFrom"));
 
             if (sourceType == "autoPatch")
@@ -1111,49 +1172,57 @@ int main(int argc, char *argv[])
                     const wordList& allNames1 = interfaceNames1[inti];
                     const auto& allFaces1 = interfaceFaces1[inti];
 
-                    forAll(allMeshes0, sourcei)
+                    if
+                    (
+                        interfaceSource0[inti] == sourcei
+                     && allMeshes0[sourcei] == meshi
+                    )
                     {
-                        if (allMeshes0[sourcei] == meshi)
-                        {
-                            const label destPatchi =
-                                patches.findPatchID(allNames0[sourcei], false);
+                        // Current mesh is mesh0. mesh1 is the remote mesh.
 
-                            Pout<< "Matched mesh:" << mesh.name()
-                                << " to mesh:"
-                                << meshes[allMeshes1[sourcei]].name()
-                                << " through:" << allNames0[sourcei] << endl;
+                        const label destPatchi =
+                            patches.findPatchID(allNames0[sourcei], false);
 
-                            changePatchID
-                            (
-                                mesh,
-                                allFaces0[sourcei],
-                                destPatchi,
-                                isRepatchedBoundary,
-                                meshMod
-                            );
-                        }
+                        //const auto& mesh1 =
+                        //    meshes[allMeshes1[interfaceSource1[inti]]];
+                        //Pout<< "Matched mesh:" << mesh.name()
+                        //    << " to mesh:" << mesh1.name()
+                        //    << " through:" << allNames0[sourcei] << endl;
+
+                        changePatchID
+                        (
+                            mesh,
+                            allFaces0[sourcei],
+                            destPatchi,
+                            isRepatchedBoundary,
+                            meshMod
+                        );
                     }
-                    forAll(allMeshes1, sourcei)
+                    if
+                    (
+                        interfaceSource1[inti] == sourcei
+                     && allMeshes1[sourcei] == meshi
+                    )
                     {
-                        if (allMeshes1[sourcei] == meshi)
-                        {
-                            const label destPatchi =
-                                patches.findPatchID(allNames1[sourcei], false);
+                        // Current mesh is mesh1. mesh0 is the remote mesh.
 
-                            Pout<< "Matched mesh:" << mesh.name()
-                                << " to mesh:"
-                                << meshes[allMeshes0[sourcei]].name()
-                                << " through:" << allNames1[sourcei] << endl;
+                        const label destPatchi =
+                            patches.findPatchID(allNames1[sourcei], false);
 
-                            changePatchID
-                            (
-                                mesh,
-                                allFaces1[sourcei],
-                                destPatchi,
-                                isRepatchedBoundary,
-                                meshMod
-                            );
-                        }
+                        //const auto& mesh0 =
+                        //    meshes[allMeshes0[interfaceSource0[inti]]];
+                        //Pout<< "Matched mesh:" << mesh.name()
+                        //    << " to mesh:" << mesh0.name()
+                        //    << " through:" << allNames1[sourcei] << endl;
+
+                        changePatchID
+                        (
+                            mesh,
+                            allFaces1[sourcei],
+                            destPatchi,
+                            isRepatchedBoundary,
+                            meshMod
+                        );
                     }
                 }
             }
@@ -1213,7 +1282,6 @@ int main(int argc, char *argv[])
                     << exit(FatalError);
             }
         }
-        Info<< endl;
 
 
         // Change mesh, use inflation to reforce calculation of transformation
@@ -1398,10 +1466,6 @@ int main(int argc, char *argv[])
         {
             dumpCyclicMatch("final_", mesh);
         }
-
-
-        // Set the precision of the points data to 10
-        IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
     }
 
     if (!overwrite)
@@ -1418,11 +1482,14 @@ int main(int argc, char *argv[])
         }
     }
 
+    // Set the precision of the points data to 10
+    IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
+
     // Write resulting mesh
     forAll(meshes, meshi)
     {
         fvMesh& mesh = meshes[meshi];
-        Info<< "Writing repatched mesh " << mesh.name()
+        Info<< "\n\nWriting repatched mesh " << mesh.name()
             << " to " << runTime.timeName() << nl << endl;
         mesh.clearOut();    // remove meshPhi
         mesh.write();
diff --git a/applications/utilities/mesh/manipulation/createPatch/createPatch.C.orig b/applications/utilities/mesh/manipulation/createPatch/createPatch.C.orig
new file mode 100644
index 0000000000..1f38a37854
--- /dev/null
+++ b/applications/utilities/mesh/manipulation/createPatch/createPatch.C.orig
@@ -0,0 +1,1439 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2016-2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    createPatch
+
+Group
+    grpMeshManipulationUtilities
+
+Description
+    Create patches out of selected boundary faces, which are either
+    from existing patches or from a faceSet.
+
+    More specifically it:
+    - creates new patches (from selected boundary faces).
+      Synchronise faces on coupled patches.
+    - synchronises points on coupled boundaries
+    - remove patches with 0 faces in them
+
+\*---------------------------------------------------------------------------*/
+
+#include "cyclicPolyPatch.H"
+#include "syncTools.H"
+#include "argList.H"
+#include "Time.H"
+#include "OFstream.H"
+#include "meshTools.H"
+#include "faceSet.H"
+#include "IOPtrList.H"
+#include "polyTopoChange.H"
+#include "polyModifyFace.H"
+#include "wordRes.H"
+#include "processorMeshes.H"
+#include "IOdictionary.H"
+#include "regionProperties.H"
+#include "faceAreaWeightAMI2D.H"
+#include "fvMeshTools.H"
+#include "ReadFields.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTemplateTypeNameAndDebug(IOPtrList<dictionary>, 0);
+}
+
+
+word patchName(const word& name, const fvMesh& mesh0, const fvMesh& mesh1)
+{
+    word pName(name != "none" ? name : word::null);
+    pName += mesh0.name();
+    pName += "_to_";
+    pName += mesh1.name();
+    return pName;
+}
+
+
+void matchPatchFaces
+(
+    const word& entryName,
+    const word& AMIMethod,
+    const dictionary& AMIDict,
+    const PtrList<fvMesh>& meshes,
+
+    const label meshi,
+    const label nSourcei,
+    const label sourcei,
+    const labelList& patchesi,
+
+    const label meshj,
+    const label nSourcej,
+    const label sourcej,
+    const labelList& patchesj,
+
+    DynamicList<labelList>& interfaceMesh0,
+    DynamicList<labelList>& interfacePatch0,
+    DynamicList<wordList>& interfaceNames0,
+    DynamicList<List<DynamicList<label>>>& interfaceFaces0,
+
+    DynamicList<labelList>& interfaceMesh1,
+    DynamicList<labelList>& interfacePatch1,
+    DynamicList<wordList>& interfaceNames1,
+    DynamicList<List<DynamicList<label>>>& interfaceFaces1
+)
+{
+    // Now we have:
+    // - meshi, sourcei, patchesi
+    // - meshj, sourcej, patchesj
+
+
+    // Attempt to match patches
+    forAll(patchesi, i)
+    {
+        const auto& ppi = meshes[meshi].boundaryMesh()[patchesi[i]];
+
+        forAll(patchesj, j)
+        {
+            const auto& ppj = meshes[meshj].boundaryMesh()[patchesj[j]];
+
+            // Use AMI to try and find matches
+            auto AMPtr(AMIInterpolation::New(AMIMethod, AMIDict));
+
+            AMPtr->calculate(ppi, ppj, nullptr);
+            if
+            (
+                gAverage(AMPtr->tgtWeightsSum()) > SMALL
+             || gAverage(AMPtr->srcWeightsSum()) > SMALL
+            )
+            {
+                const label inti = interfaceMesh0.size();
+
+                Info<< "Introducing interface " << inti << " between"
+                    << " mesh " << meshes[meshi].name()
+                    << " patch " << ppi.name()
+                    << " and mesh " << meshes[meshj].name()
+                    << " patch " << ppj.name()
+                    << endl;
+
+                // Mesh 0
+                //~~~~~~~
+
+                interfaceMesh0.append(labelList());
+                auto& intMesh0 = interfaceMesh0.last();
+                intMesh0.setSize(nSourcei, -1);
+                intMesh0[sourcei] = meshi;
+
+                interfacePatch0.append(labelList());
+                auto& intPatch0 = interfacePatch0.last();
+                intPatch0.setSize(nSourcei, -1);
+                intPatch0[sourcei] = ppi.index();
+
+                interfaceNames0.append(wordList());
+                auto& intNames0 = interfaceNames0.last();
+                intNames0.setSize(nSourcei);
+                //intNames0[sourcei] =
+                //    meshes[meshi].name()
+                //  + "_to_"
+                //  + meshes[meshj].name();
+                intNames0[sourcei] =
+                    patchName(entryName, meshes[meshi], meshes[meshj]);
+
+
+                // Mesh 0
+                //~~~~~~~
+
+                interfaceMesh1.append(labelList());
+                auto& intMesh1 = interfaceMesh1.last();
+                intMesh1.setSize(nSourcej, -1);
+                intMesh1[sourcej] = meshj;
+
+                interfacePatch1.append(labelList());
+                auto& intPatch1 = interfacePatch1.last();
+                intPatch1.setSize(nSourcej, -1);
+                intPatch1[sourcej] = ppj.index();
+
+                interfaceNames1.append(wordList());
+                auto& intNames1 = interfaceNames1.last();
+                intNames1.setSize(nSourcej);
+                //intNames1[sourcej] =
+                //    meshes[meshj].name()
+                //  + "_to_"
+                //  + meshes[meshi].name();
+                intNames1[sourcej] =
+                    patchName(entryName, meshes[meshj], meshes[meshi]);
+
+                interfaceFaces0.append(List<DynamicList<label>>());
+                auto& intFaces0 = interfaceFaces0.last();
+                intFaces0.setSize(nSourcei);
+                DynamicList<label>& faces0 = intFaces0[sourcei];
+                faces0.setCapacity(ppi.size());
+
+                interfaceFaces1.append(List<DynamicList<label>>());
+                auto& intFaces1 = interfaceFaces1.last();
+                intFaces1.setSize(nSourcej);
+                DynamicList<label>& faces1 = intFaces1[sourcej];
+                faces1.setCapacity(ppj.size());
+
+
+                // Mark any mesh faces:
+                // - that have a weight > 0.5.
+                // - contribute to these faces (even if < 0.5)
+
+                faces0.clear();
+                faces1.clear();
+
+                // Marked as donor of face with high weight
+                scalarField targetMask;
+                {
+                    const scalarField& weights = AMPtr->srcWeightsSum();
+                    scalarField mask(weights.size(), Zero);
+                    forAll(weights, facei)
+                    {
+                        const scalar sum = weights[facei];
+                        if (sum > 0.5)
+                        {
+                            mask[facei] = sum;
+                        }
+                    }
+
+                    // Push source field mask to target
+                    targetMask = AMPtr->interpolateToTarget(mask);
+                }
+
+                scalarField sourceMask;
+                {
+                    const scalarField& weights = AMPtr->tgtWeightsSum();
+                    scalarField mask(weights.size(), Zero);
+                    forAll(weights, facei)
+                    {
+                        const scalar sum = weights[facei];
+                        if (sum > 0.5)
+                        {
+                            mask[facei] = sum;
+                        }
+                    }
+
+                    // Push target mask back to source
+                    sourceMask = AMPtr->interpolateToSource(mask);
+                }
+
+                {
+                    const scalarField& weights = AMPtr->srcWeightsSum();
+                    forAll(weights, facei)
+                    {
+                        if (weights[facei] > 0.5 || sourceMask[facei] > SMALL)
+                        {
+                            faces0.append(ppi.start()+facei);
+                        }
+                    }
+                }
+                {
+                    const scalarField& weights = AMPtr->tgtWeightsSum();
+                    forAll(weights, facei)
+                    {
+                        if (weights[facei] > 0.5 || targetMask[facei] > SMALL)
+                        {
+                            faces1.append(ppj.start()+facei);
+                        }
+                    }
+                }
+
+                faces0.shrink();
+                faces1.shrink();
+            }
+        }
+    }
+}
+
+
+void matchPatchFaces
+(
+    const bool includeOwn,
+    const PtrList<fvMesh>& meshes,
+    const List<DynamicList<label>>& interRegionSources,
+    const List<wordList>& patchNames,
+    const labelListListList& matchPatchIDs,
+    const List<wordList>& AMIMethods,
+    List<PtrList<dictionary>> patchInfoDicts,
+
+    DynamicList<labelList>& interfaceMesh0,
+    DynamicList<labelList>& interfacePatch0,
+    DynamicList<List<DynamicList<label>>>& interfaceFaces0,
+    DynamicList<wordList>& interfaceNames0,
+
+    DynamicList<labelList>& interfaceMesh1,
+    DynamicList<labelList>& interfacePatch1,
+    DynamicList<List<DynamicList<label>>>& interfaceFaces1,
+    DynamicList<wordList>& interfaceNames1
+)
+{
+    // Add approximate matches
+    forAll(meshes, meshi)
+    {
+        const labelListList& allPatchesi = matchPatchIDs[meshi];
+
+        const label meshjStart(includeOwn ? meshi : meshi+1);
+
+        for (label meshj = meshjStart; meshj < meshes.size(); meshj++)
+        {
+            const labelListList& allPatchesj = matchPatchIDs[meshj];
+
+            for (const label sourcei : interRegionSources[meshi])
+            {
+                const labelList& patchesi = allPatchesi[sourcei];
+
+                for (const label sourcej : interRegionSources[meshj])
+                {
+                    const labelList& patchesj = allPatchesj[sourcej];
+
+                    // Now we have:
+                    // - meshi, sourcei, patchesi
+                    // - meshj, sourcej, patchesj
+
+                    matchPatchFaces
+                    (
+                        patchNames[meshi][sourcei],
+                        (
+                            AMIMethods[meshi][sourcei].size()
+                          ? AMIMethods[meshi][sourcei]
+                          : faceAreaWeightAMI2D::typeName
+                        ),
+                        patchInfoDicts[meshi][sourcei],
+                        meshes,
+
+                        meshi,
+                        allPatchesi.size(), // nSourcei,
+                        sourcei,
+                        patchesi,
+
+                        meshj,
+                        allPatchesj.size(), // nSourcej,
+                        sourcej,
+                        patchesj,
+
+                        interfaceMesh0,
+                        interfacePatch0,
+                        interfaceNames0,
+                        interfaceFaces0,
+
+                        interfaceMesh1,
+                        interfacePatch1,
+                        interfaceNames1,
+                        interfaceFaces1
+                    );
+                }
+            }
+        }
+    }
+}
+
+
+void changePatchID
+(
+    const fvMesh& mesh,
+    const label faceID,
+    const label patchID,
+    polyTopoChange& meshMod
+)
+{
+    const label zoneID = mesh.faceZones().whichZone(faceID);
+
+    bool zoneFlip = false;
+
+    if (zoneID >= 0)
+    {
+        const faceZone& fZone = mesh.faceZones()[zoneID];
+
+        zoneFlip = fZone.flipMap()[fZone.whichFace(faceID)];
+    }
+
+    meshMod.setAction
+    (
+        polyModifyFace
+        (
+            mesh.faces()[faceID],               // face
+            faceID,                             // face ID
+            mesh.faceOwner()[faceID],           // owner
+            -1,                                 // neighbour
+            false,                              // flip flux
+            patchID,                            // patch ID
+            false,                              // remove from zone
+            zoneID,                             // zone ID
+            zoneFlip                            // zone flip
+        )
+    );
+}
+
+
+void changePatchID
+(
+    const fvMesh& mesh,
+    const labelList& faceLabels,
+    const label patchID,
+    bitSet& isRepatchedBoundary,
+    polyTopoChange& meshMod
+)
+{
+    for (const label facei : faceLabels)
+    {
+        if (mesh.isInternalFace(facei))
+        {
+            FatalErrorInFunction
+                << "Face " << facei
+                << " is not an external face of the mesh." << endl
+                << "This application can only repatch"
+                << " existing boundary faces." << exit(FatalError);
+        }
+
+        if (!isRepatchedBoundary.set(facei-mesh.nInternalFaces()))
+        {
+            FatalErrorInFunction
+                << "Face " << facei << " is already marked to be moved"
+                << " to patch " << meshMod.region()[facei]
+                << exit(FatalError);
+        }
+
+        changePatchID(mesh, facei, patchID, meshMod);
+    }
+}
+
+
+// Dump for all patches the current match
+void dumpCyclicMatch(const fileName& prefix, const polyMesh& mesh)
+{
+    for (const polyPatch& pp : mesh.boundaryMesh())
+    {
+        const cyclicPolyPatch* cpp = isA<cyclicPolyPatch>(pp);
+
+        if (cpp && cpp->owner())
+        {
+            const auto& cycPatch = *cpp;
+            const auto& nbrPatch = cycPatch.neighbPatch();
+
+            // Dump patches
+            {
+                OFstream str(prefix+cycPatch.name()+".obj");
+                Pout<< "Dumping " << cycPatch.name()
+                    << " faces to " << str.name() << endl;
+                meshTools::writeOBJ
+                (
+                    str,
+                    cycPatch,
+                    cycPatch.points()
+                );
+            }
+
+            {
+                OFstream str(prefix+nbrPatch.name()+".obj");
+                Pout<< "Dumping " << nbrPatch.name()
+                    << " faces to " << str.name() << endl;
+                meshTools::writeOBJ
+                (
+                    str,
+                    nbrPatch,
+                    nbrPatch.points()
+                );
+            }
+
+
+            // Lines between corresponding face centres
+            OFstream str(prefix+cycPatch.name()+nbrPatch.name()+"_match.obj");
+            label vertI = 0;
+
+            Pout<< "Dumping cyclic match as lines between face centres to "
+                << str.name() << endl;
+
+            forAll(cycPatch, facei)
+            {
+                const point& fc0 = mesh.faceCentres()[cycPatch.start()+facei];
+                meshTools::writeOBJ(str, fc0);
+                vertI++;
+                const point& fc1 = mesh.faceCentres()[nbrPatch.start()+facei];
+                meshTools::writeOBJ(str, fc1);
+                vertI++;
+
+                str<< "l " << vertI-1 << ' ' << vertI << nl;
+            }
+        }
+    }
+}
+
+
+void separateList
+(
+    const vectorField& separation,
+    UList<vector>& field
+)
+{
+    if (separation.size() == 1)
+    {
+        // Single value for all.
+
+        forAll(field, i)
+        {
+            field[i] += separation[0];
+        }
+    }
+    else if (separation.size() == field.size())
+    {
+        forAll(field, i)
+        {
+            field[i] += separation[i];
+        }
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Sizes of field and transformation not equal. field:"
+            << field.size() << " transformation:" << separation.size()
+            << abort(FatalError);
+    }
+}
+
+
+// Synchronise points on both sides of coupled boundaries.
+template<class CombineOp>
+void syncPoints
+(
+    const polyMesh& mesh,
+    pointField& points,
+    const CombineOp& cop,
+    const point& nullValue
+)
+{
+    if (points.size() != mesh.nPoints())
+    {
+        FatalErrorInFunction
+            << "Number of values " << points.size()
+            << " is not equal to the number of points in the mesh "
+            << mesh.nPoints() << abort(FatalError);
+    }
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    // Is there any coupled patch with transformation?
+    bool hasTransformation = false;
+
+    if (Pstream::parRun())
+    {
+        const labelList& procPatches = mesh.globalData().processorPatches();
+
+        // Send
+        for (const label patchi : procPatches)
+        {
+            const polyPatch& pp = patches[patchi];
+            const auto& procPatch = refCast<const processorPolyPatch>(pp);
+
+            if (pp.nPoints() && procPatch.owner())
+            {
+                // Get data per patchPoint in neighbouring point numbers.
+                pointField patchInfo(procPatch.nPoints(), nullValue);
+
+                const labelList& meshPts = procPatch.meshPoints();
+                const labelList& nbrPts = procPatch.neighbPoints();
+
+                forAll(nbrPts, pointi)
+                {
+                    label nbrPointi = nbrPts[pointi];
+                    if (nbrPointi >= 0 && nbrPointi < patchInfo.size())
+                    {
+                        patchInfo[nbrPointi] = points[meshPts[pointi]];
+                    }
+                }
+
+                OPstream toNbr
+                (
+                    Pstream::commsTypes::blocking,
+                    procPatch.neighbProcNo()
+                );
+                toNbr << patchInfo;
+            }
+        }
+
+
+        // Receive and set.
+
+        for (const label patchi : procPatches)
+        {
+            const polyPatch& pp = patches[patchi];
+            const auto& procPatch = refCast<const processorPolyPatch>(pp);
+
+            if (pp.nPoints() && !procPatch.owner())
+            {
+                pointField nbrPatchInfo(procPatch.nPoints());
+                {
+                    // We do not know the number of points on the other side
+                    // so cannot use Pstream::read.
+                    IPstream fromNbr
+                    (
+                        Pstream::commsTypes::blocking,
+                        procPatch.neighbProcNo()
+                    );
+                    fromNbr >> nbrPatchInfo;
+                }
+                // Null any value which is not on neighbouring processor
+                nbrPatchInfo.setSize(procPatch.nPoints(), nullValue);
+
+                if (!procPatch.parallel())
+                {
+                    hasTransformation = true;
+                    transformList(procPatch.forwardT(), nbrPatchInfo);
+                }
+                else if (procPatch.separated())
+                {
+                    hasTransformation = true;
+                    separateList(-procPatch.separation(), nbrPatchInfo);
+                }
+
+                const labelList& meshPts = procPatch.meshPoints();
+
+                forAll(meshPts, pointi)
+                {
+                    label meshPointi = meshPts[pointi];
+                    points[meshPointi] = nbrPatchInfo[pointi];
+                }
+            }
+        }
+    }
+
+    // Do the cyclics.
+    for (const polyPatch& pp : patches)
+    {
+        const cyclicPolyPatch* cpp = isA<cyclicPolyPatch>(pp);
+
+        if (cpp && cpp->owner())
+        {
+            // Owner does all.
+
+            const auto& cycPatch = *cpp;
+            const auto& nbrPatch = cycPatch.neighbPatch();
+
+            const edgeList& coupledPoints = cycPatch.coupledPoints();
+            const labelList& meshPts = cycPatch.meshPoints();
+            const labelList& nbrMeshPts = nbrPatch.meshPoints();
+
+            pointField half0Values(coupledPoints.size());
+
+            forAll(coupledPoints, i)
+            {
+                const edge& e = coupledPoints[i];
+                label point0 = meshPts[e[0]];
+                half0Values[i] = points[point0];
+            }
+
+            if (!cycPatch.parallel())
+            {
+                hasTransformation = true;
+                transformList(cycPatch.reverseT(), half0Values);
+            }
+            else if (cycPatch.separated())
+            {
+                hasTransformation = true;
+                separateList(cycPatch.separation(), half0Values);
+            }
+
+            forAll(coupledPoints, i)
+            {
+                const edge& e = coupledPoints[i];
+                label point1 = nbrMeshPts[e[1]];
+                points[point1] = half0Values[i];
+            }
+        }
+    }
+
+    //- Note: hasTransformation is only used for warning messages so
+    //  reduction not strictly necessary.
+    //reduce(hasTransformation, orOp<bool>());
+
+    // Synchronize multiple shared points.
+    const globalMeshData& pd = mesh.globalData();
+
+    if (pd.nGlobalPoints() > 0)
+    {
+        if (hasTransformation)
+        {
+            WarningInFunction
+                << "There are decomposed cyclics in this mesh with"
+                << " transformations." << endl
+                << "This is not supported. The result will be incorrect"
+                << endl;
+        }
+
+
+        // Values on shared points.
+        pointField sharedPts(pd.nGlobalPoints(), nullValue);
+
+        forAll(pd.sharedPointLabels(), i)
+        {
+            label meshPointi = pd.sharedPointLabels()[i];
+            // Fill my entries in the shared points
+            sharedPts[pd.sharedPointAddr()[i]] = points[meshPointi];
+        }
+
+        // Combine on master.
+        Pstream::listCombineGather(sharedPts, cop);
+        Pstream::listCombineScatter(sharedPts);
+
+        // Now we will all have the same information. Merge it back with
+        // my local information.
+        forAll(pd.sharedPointLabels(), i)
+        {
+            label meshPointi = pd.sharedPointLabels()[i];
+            points[meshPointi] = sharedPts[pd.sharedPointAddr()[i]];
+        }
+    }
+}
+
+
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Create patches out of selected boundary faces, which are either"
+        " from existing patches or from a faceSet"
+    );
+
+    #include "addOverwriteOption.H"
+    //#include "addRegionOption.H"
+    #include "addAllRegionOptions.H"
+
+    argList::addOption("dict", "file", "Alternative createPatchDict");
+    argList::addBoolOption
+    (
+        "writeObj",
+        "Write obj files showing the cyclic matching process"
+    );
+
+    argList::noFunctionObjects();  // Never use function objects
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "getAllRegionOptions.H"
+
+    const bool overwrite = args.found("overwrite");
+
+    #include "createNamedMeshes.H"
+
+    const bool writeObj = args.found("writeObj");
+
+
+    // Read dictionaries and extract various
+    wordList oldInstances(meshes.size());
+    PtrList<dictionary> dicts(meshes.size());
+    List<wordList> patchNames(meshes.size());
+    List<labelListList> matchPatchIDs(meshes.size());
+    List<wordList> matchMethods(meshes.size());
+    List<PtrList<dictionary>> patchInfoDicts(meshes.size());
+
+    List<DynamicList<label>> interRegionSources(meshes.size());
+    forAll(meshes, meshi)
+    {
+        fvMesh& mesh = meshes[meshi];
+        const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+        // If running parallel check same patches everywhere
+        patches.checkParallelSync(true);
+
+        oldInstances[meshi] = mesh.pointsInstance();
+
+        const word dictName("createPatchDict");
+        #include "setSystemMeshDictionaryIO.H"
+        Info<< "Reading " << dictIO.instance()/dictIO.name() << nl << endl;
+
+        dicts.set(meshi, new IOdictionary(dictIO));
+
+        const auto& dict = dicts[meshi];
+        PtrList<dictionary> patchSources(dict.lookup("patches"));
+        patchNames[meshi].setSize(patchSources.size());
+        matchPatchIDs[meshi].setSize(patchSources.size());
+        matchMethods[meshi].setSize(patchSources.size());
+        patchInfoDicts[meshi].setSize(patchSources.size());
+
+        // Read patch construct info from dictionary
+        forAll(patchSources, sourcei)
+        {
+            const auto& pDict = patchSources[sourcei];
+            patchNames[meshi][sourcei] = pDict.get<word>("name");
+
+            patchInfoDicts[meshi].set
+            (
+                sourcei,
+                new dictionary(pDict.subDict("patchInfo"))
+            );
+            const dictionary& patchDict = patchInfoDicts[meshi][sourcei];
+            if (patchDict.found("AMIMethod"))
+            {
+                matchMethods[meshi][sourcei] = patchDict.get<word>("AMIMethod");
+            }
+
+            wordRes matchNames;
+            if (pDict.readIfPresent("patches", matchNames, keyType::LITERAL))
+            {
+                matchPatchIDs[meshi][sourcei] =
+                    patches.patchSet(matchNames).sortedToc();
+            }
+
+            if (pDict.get<word>("constructFrom") == "autoPatch")
+            {
+                interRegionSources[meshi].append(sourcei);
+            }
+        }
+    }
+
+
+
+    // Determine matching faces. This is an interface between two
+    // sets of faces:
+    //  - originating mesh
+    //  - originating patch
+    //  - originating (mesh)facelabels
+    // It matches all mesh against each other. Lower numbered mesh gets
+    // postfix 0, higher numbered mesh postfix 1.
+
+    // Per interface, per mesh, per patchSource:
+    //      1. the lower numbered mesh
+    DynamicList<labelList> interfaceMesh0;
+    //      2. the patch on the interfaceMesh0
+    DynamicList<labelList> interfacePatch0;
+    //      3. the facelabels on the interfaceMesh0
+    DynamicList<List<DynamicList<label>>> interfaceFaces0;
+    //      4. generated interface name
+    DynamicList<wordList> interfaceNames0;
+
+    // Same for the higher numbered mesh
+    DynamicList<labelList> interfaceMesh1;
+    DynamicList<labelList> interfacePatch1;
+    DynamicList<List<DynamicList<label>>> interfaceFaces1;
+    DynamicList<wordList> interfaceNames1;
+    {
+        // Whether to match to patches in own mesh
+        const bool includeOwn = (meshes.size() == 1);
+
+        matchPatchFaces
+        (
+            includeOwn,
+            meshes,
+            interRegionSources,
+            patchNames,
+            matchPatchIDs,
+            matchMethods,   //faceAreaWeightAMI2D::typeName,
+            patchInfoDicts,
+
+            interfaceMesh0,
+            interfacePatch0,
+            interfaceFaces0,
+            interfaceNames0,
+
+            interfaceMesh1,
+            interfacePatch1,
+            interfaceFaces1,
+            interfaceNames1
+        );
+    }
+
+
+    // Read fields
+    List<PtrList<volScalarField>> vsFlds(meshes.size());
+    List<PtrList<volVectorField>> vvFlds(meshes.size());
+    List<PtrList<volSphericalTensorField>> vstFlds(meshes.size());
+    List<PtrList<volSymmTensorField>> vsymtFlds(meshes.size());
+    List<PtrList<surfaceScalarField>> ssFlds(meshes.size());
+    List<PtrList<volTensorField>> vtFlds(meshes.size());
+    List<PtrList<surfaceVectorField>> svFlds(meshes.size());
+    List<PtrList<surfaceSphericalTensorField>> sstFlds(meshes.size());
+    List<PtrList<surfaceSymmTensorField>> ssymtFlds(meshes.size());
+    List<PtrList<surfaceTensorField>> stFlds(meshes.size());
+
+    {
+        forAll(meshes, meshi)
+        {
+            const fvMesh& mesh = meshes[meshi];
+
+            bool noFields = true;
+            for (const auto& d : patchInfoDicts[meshi])
+            {
+                if (d.found("patchFields"))
+                {
+                    noFields = false;
+                }
+            }
+
+            if (!noFields)
+            {
+                // Read objects in time directory
+                IOobjectList objects(mesh, runTime.timeName());
+
+                // Read volume fields.
+                ReadFields(mesh, objects, vsFlds[meshi]);
+                ReadFields(mesh, objects, vvFlds[meshi]);
+                ReadFields(mesh, objects, vstFlds[meshi]);
+                ReadFields(mesh, objects, vsymtFlds[meshi]);
+                ReadFields(mesh, objects, vtFlds[meshi]);
+
+                // Read surface fields.
+                ReadFields(mesh, objects, ssFlds[meshi]);
+                ReadFields(mesh, objects, svFlds[meshi]);
+                ReadFields(mesh, objects, sstFlds[meshi]);
+                ReadFields(mesh, objects, ssymtFlds[meshi]);
+                ReadFields(mesh, objects, stFlds[meshi]);
+            }
+        }
+    }
+
+
+
+    // Loop over all regions
+
+    forAll(meshes, meshi)
+    {
+        fvMesh& mesh = meshes[meshi];
+        const polyBoundaryMesh& patches = mesh.boundaryMesh();
+        const dictionary& dict = dicts[meshi];
+
+        if (writeObj)
+        {
+            dumpCyclicMatch("initial_", mesh);
+        }
+
+        // Read patch construct info from dictionary
+        PtrList<dictionary> patchSources(dict.lookup("patches"));
+
+
+        // 1. Add all new patches
+        // ~~~~~~~~~~~~~~~~~~~~~~
+
+        forAll(patchSources, sourcei)
+        {
+            const dictionary& dict = patchSources[sourcei];
+            const word sourceType(dict.get<word>("constructFrom"));
+            const word patchName(dict.get<word>("name"));
+
+            dictionary patchDict(dict.subDict("patchInfo"));
+            patchDict.set("nFaces", 0);
+            patchDict.set("startFace", 0);  //startFacei);
+
+            if (sourceType == "autoPatch")
+            {
+                // Special : automatically create the necessary inter-region
+                // coupling patches
+
+                forAll(interfaceMesh0, inti)
+                {
+                    const labelList& allMeshes0 = interfaceMesh0[inti];
+                    const wordList& allNames0 = interfaceNames0[inti];
+                    const labelList& allMeshes1 = interfaceMesh1[inti];
+                    const wordList& allNames1 = interfaceNames1[inti];
+
+                    forAll(allMeshes0, sourcei)
+                    {
+                        if (allMeshes0[sourcei] == meshi)
+                        {
+                            const auto& mesh1 = meshes[allMeshes1[sourcei]];
+                            const word& patchName = allNames0[sourcei];
+                            if (patches.findPatchID(patchName) == -1)
+                            {
+                                dictionary allDict(patchDict);
+                                allDict.set("sampleRegion", mesh1.name());
+                                const auto& destPatch = allNames1[sourcei];
+                                allDict.set("samplePatch", destPatch);
+                                allDict.set("neighbourPatch", destPatch);
+
+                                Info<< "Adding new patch " << patchName
+                                    << " from " << allDict << endl;
+
+                                autoPtr<polyPatch> ppPtr
+                                (
+                                    polyPatch::New
+                                    (
+                                        patchName,
+                                        allDict,
+                                        0,          // overwritten
+                                        patches
+                                    )
+                                );
+                                fvMeshTools::addPatch
+                                (
+                                    mesh,
+                                    ppPtr(),
+                                    patchDict.subOrEmptyDict("patchFields"),
+                                    calculatedFvPatchScalarField::typeName,
+                                    true
+                                );
+                            }
+                        }
+                    }
+                }
+
+                forAll(interfaceMesh1, inti)
+                {
+                    const labelList& allMeshes0 = interfaceMesh0[inti];
+                    const wordList& allNames0 = interfaceNames0[inti];
+                    const labelList& allMeshes1 = interfaceMesh1[inti];
+                    const wordList& allNames1 = interfaceNames1[inti];
+
+                    forAll(allMeshes1, sourcei)
+                    {
+                        if (allMeshes1[sourcei] == meshi)
+                        {
+                            const auto& mesh0 = meshes[allMeshes0[sourcei]];
+                            const word& patchName = allNames1[sourcei];
+                            if (patches.findPatchID(patchName) == -1)
+                            {
+                                dictionary allDict(patchDict);
+                                const auto& destPatch = allNames0[sourcei];
+                                allDict.set("sampleRegion", mesh0.name());
+                                allDict.set("samplePatch", destPatch);
+                                allDict.set("neighbourPatch", destPatch);
+
+                                Info<< "Adding new patch " << patchName
+                                    << " from " << allDict << endl;
+
+                                autoPtr<polyPatch> ppPtr
+                                (
+                                    polyPatch::New
+                                    (
+                                        patchName,
+                                        allDict,
+                                        0,          // overwritten
+                                        patches
+                                    )
+                                );
+                                fvMeshTools::addPatch
+                                (
+                                    mesh,
+                                    ppPtr(),
+                                    patchDict.subOrEmptyDict("patchFields"),
+                                    calculatedFvPatchScalarField::typeName,
+                                    true
+                                );
+                            }
+                        }
+                    }
+                }
+            }
+            else
+            {
+                if (patches.findPatchID(patchName) == -1)
+                {
+                    Info<< "Adding new patch " << patchName
+                        << " from " << patchDict << endl;
+
+                    autoPtr<polyPatch> ppPtr
+                    (
+                        polyPatch::New
+                        (
+                            patchName,
+                            patchDict,
+                            0,          // overwritten
+                            patches
+                        )
+                    );
+                    fvMeshTools::addPatch
+                    (
+                        mesh,
+                        ppPtr(),
+                        patchDict.subOrEmptyDict("patchFields"),
+                        calculatedFvPatchScalarField::typeName,
+                        true
+                    );
+                }
+            }
+        }
+
+        Info<< endl;
+    }
+
+
+    // 2. Repatch faces
+    // ~~~~~~~~~~~~~~~~
+
+    forAll(meshes, meshi)
+    {
+        fvMesh& mesh = meshes[meshi];
+        const polyBoundaryMesh& patches = mesh.boundaryMesh();
+        const dictionary& dict = dicts[meshi];
+
+        // Whether to synchronise points
+        const bool pointSync(dict.get<bool>("pointSync"));
+
+        // Read patch construct info from dictionary
+        PtrList<dictionary> patchSources(dict.lookup("patches"));
+
+
+        polyTopoChange meshMod(mesh);
+
+        // Mark all repatched faces. This makes sure that the faces to repatch
+        // do not overlap
+        bitSet isRepatchedBoundary(mesh.nBoundaryFaces());
+
+        forAll(patchSources, sourcei)
+        {
+            const dictionary& dict = patchSources[sourcei];
+            const word patchName(dict.get<word>("name"));
+            const word sourceType(dict.get<word>("constructFrom"));
+
+            if (sourceType == "autoPatch")
+            {
+                forAll(interfaceMesh0, inti)
+                {
+                    const labelList& allMeshes0 = interfaceMesh0[inti];
+                    const wordList& allNames0 = interfaceNames0[inti];
+                    const auto& allFaces0 = interfaceFaces0[inti];
+
+                    const labelList& allMeshes1 = interfaceMesh1[inti];
+                    const wordList& allNames1 = interfaceNames1[inti];
+                    const auto& allFaces1 = interfaceFaces1[inti];
+
+                    forAll(allMeshes0, sourcei)
+                    {
+                        if (allMeshes0[sourcei] == meshi)
+                        {
+                            const label destPatchi =
+                                patches.findPatchID(allNames0[sourcei], false);
+
+                            Pout<< "Matched mesh:" << mesh.name()
+                                << " to mesh:"
+                                << meshes[allMeshes1[sourcei]].name()
+                                << " through:" << allNames0[sourcei] << endl;
+
+                            changePatchID
+                            (
+                                mesh,
+                                allFaces0[sourcei],
+                                destPatchi,
+                                isRepatchedBoundary,
+                                meshMod
+                            );
+                        }
+                    }
+                    forAll(allMeshes1, sourcei)
+                    {
+                        if (allMeshes1[sourcei] == meshi)
+                        {
+                            const label destPatchi =
+                                patches.findPatchID(allNames1[sourcei], false);
+
+                            Pout<< "Matched mesh:" << mesh.name()
+                                << " to mesh:"
+                                << meshes[allMeshes0[sourcei]].name()
+                                << " through:" << allNames1[sourcei] << endl;
+
+                            changePatchID
+                            (
+                                mesh,
+                                allFaces1[sourcei],
+                                destPatchi,
+                                isRepatchedBoundary,
+                                meshMod
+                            );
+                        }
+                    }
+                }
+            }
+            else if (sourceType == "patches")
+            {
+                const label destPatchi = patches.findPatchID(patchName, false);
+                labelHashSet patchSources
+                (
+                    patches.patchSet(dict.get<wordRes>("patches"))
+                );
+
+                // Repatch faces of the patches.
+                for (const label patchi : patchSources)
+                {
+                    const polyPatch& pp = patches[patchi];
+
+                    Info<< "Moving faces from patch " << pp.name()
+                        << " to patch " << destPatchi << endl;
+
+                    changePatchID
+                    (
+                        mesh,
+                        pp.start()+identity(pp.size()),
+                        destPatchi,
+                        isRepatchedBoundary,
+                        meshMod
+                    );
+                }
+            }
+            else if (sourceType == "set")
+            {
+                const label destPatchi = patches.findPatchID(patchName, false);
+                const word setName(dict.get<word>("set"));
+
+                faceSet set(mesh, setName);
+
+                Info<< "Read " << returnReduce(set.size(), sumOp<label>())
+                    << " faces from faceSet " << set.name() << endl;
+
+                // Sort (since faceSet contains faces in arbitrary order)
+                labelList faceLabels(set.sortedToc());
+
+                changePatchID
+                (
+                    mesh,
+                    faceLabels,
+                    destPatchi,
+                    isRepatchedBoundary,
+                    meshMod
+                );
+            }
+            else
+            {
+                FatalErrorInFunction
+                    << "Invalid source type " << sourceType << endl
+                    << "Valid source types are 'patches' 'set'"
+                    << exit(FatalError);
+            }
+        }
+        Info<< endl;
+
+
+        // Change mesh, use inflation to reforce calculation of transformation
+        // tensors.
+        Info<< "Doing topology modification to order faces." << nl << endl;
+        autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, true);
+
+        if (map().hasMotionPoints())
+        {
+            mesh.movePoints(map().preMotionPoints());
+        }
+        else
+        {
+           // Force calculation of transformation tensors
+            mesh.movePoints(pointField(mesh.points()));
+        }
+
+
+        // Update fields
+        mesh.updateMesh(map());
+
+
+        // Update numbering pointing to meshi
+        const auto& oldToNew = map().reverseFaceMap();
+        forAll(interfaceMesh0, inti)
+        {
+            const labelList& allMeshes0 = interfaceMesh0[inti];
+            forAll(allMeshes0, sourcei)
+            {
+                if (allMeshes0[sourcei] == meshi)
+                {
+                    inplaceRenumber(oldToNew, interfaceFaces0[inti][sourcei]);
+                }
+            }
+        }
+        forAll(interfaceMesh1, inti)
+        {
+            const labelList& allMeshes1 = interfaceMesh1[inti];
+            forAll(allMeshes1, sourcei)
+            {
+                if (allMeshes1[sourcei] == meshi)
+                {
+                    inplaceRenumber(oldToNew, interfaceFaces1[inti][sourcei]);
+                }
+            }
+        }
+
+
+        if (writeObj)
+        {
+            dumpCyclicMatch("coupled_", mesh);
+        }
+
+        // Synchronise points.
+        if (!pointSync)
+        {
+            Info<< "Not synchronising points." << nl << endl;
+        }
+        else
+        {
+            Info<< "Synchronising points." << nl << endl;
+
+            // This is a bit tricky. Both normal and position might be out and
+            // current separation also includes the normal
+            // ( separation_ = (nf&(Cr - Cf))*nf ).
+
+            // For cyclic patches:
+            // - for separated ones use user specified offset vector
+
+            forAll(mesh.boundaryMesh(), patchi)
+            {
+                const polyPatch& pp = mesh.boundaryMesh()[patchi];
+
+                if (pp.size() && isA<coupledPolyPatch>(pp))
+                {
+                    const coupledPolyPatch& cpp =
+                        refCast<const coupledPolyPatch>(pp);
+
+                    if (cpp.separated())
+                    {
+                        Info<< "On coupled patch " << pp.name()
+                            << " separation[0] was "
+                            << cpp.separation()[0] << endl;
+
+                        if (isA<cyclicPolyPatch>(pp) && pp.size())
+                        {
+                            const auto& cycpp =
+                                refCast<const cyclicPolyPatch>(pp);
+
+                            if
+                            (
+                                cycpp.transform()
+                             == cyclicPolyPatch::TRANSLATIONAL
+                            )
+                            {
+                                // Force to wanted separation
+                                Info<< "On cyclic translation patch "
+                                    << pp.name()
+                                    << " forcing uniform separation of "
+                                    << cycpp.separationVector() << endl;
+                                const_cast<vectorField&>(cpp.separation()) =
+                                    pointField(1, cycpp.separationVector());
+                            }
+                            else
+                            {
+                                const auto& nbr = cycpp.neighbPatch();
+                                const_cast<vectorField&>(cpp.separation()) =
+                                    pointField
+                                    (
+                                        1,
+                                        nbr[0].centre(mesh.points())
+                                      - cycpp[0].centre(mesh.points())
+                                    );
+                            }
+                        }
+                        Info<< "On coupled patch " << pp.name()
+                            << " forcing uniform separation of "
+                            << cpp.separation() << endl;
+                    }
+                    else if (!cpp.parallel())
+                    {
+                        Info<< "On coupled patch " << pp.name()
+                            << " forcing uniform rotation of "
+                            << cpp.forwardT()[0] << endl;
+
+                        const_cast<tensorField&>
+                        (
+                            cpp.forwardT()
+                        ).setSize(1);
+                        const_cast<tensorField&>
+                        (
+                            cpp.reverseT()
+                        ).setSize(1);
+
+                        Info<< "On coupled patch " << pp.name()
+                            << " forcing uniform rotation of "
+                            << cpp.forwardT() << endl;
+                    }
+                }
+            }
+
+            Info<< "Synchronising points." << endl;
+
+            pointField newPoints(mesh.points());
+
+            syncPoints
+            (
+                mesh,
+                newPoints,
+                minMagSqrEqOp<vector>(),
+                point(GREAT, GREAT, GREAT)
+            );
+
+            scalarField diff(mag(newPoints-mesh.points()));
+            Info<< "Points changed by average:" << gAverage(diff)
+                << " max:" << gMax(diff) << nl << endl;
+
+            mesh.movePoints(newPoints);
+        }
+
+
+        // 3. Remove zeros-sized patches
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        Info<< "Removing patches with no faces in them." << nl << endl;
+        const wordList oldPatchNames(mesh.boundaryMesh().names());
+        const wordList oldPatchTypes(mesh.boundaryMesh().types());
+        fvMeshTools::removeEmptyPatches(mesh, true);
+        forAll(oldPatchNames, patchi)
+        {
+            const word& pName = oldPatchNames[patchi];
+            if (mesh.boundaryMesh().findPatchID(pName) == -1)
+            {
+                Info<< "Removed zero-sized patch " << pName
+                    << " type " << oldPatchTypes[patchi]
+                    << " at position " << patchi << endl;
+            }
+        }
+
+
+        if (writeObj)
+        {
+            dumpCyclicMatch("final_", mesh);
+        }
+
+
+        // Set the precision of the points data to 10
+        IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
+    }
+
+    if (!overwrite)
+    {
+        ++runTime;
+    }
+    else
+    {
+        forAll(meshes, meshi)
+        {
+            fvMesh& mesh = meshes[meshi];
+
+            mesh.setInstance(oldInstances[meshi]);
+        }
+    }
+
+    // Write resulting mesh
+    forAll(meshes, meshi)
+    {
+        fvMesh& mesh = meshes[meshi];
+        Info<< "Writing repatched mesh " << mesh.name()
+            << " to " << runTime.timeName() << nl << endl;
+        mesh.clearOut();    // remove meshPhi
+        mesh.write();
+        topoSet::removeFiles(mesh);
+        processorMeshes::removeFiles(mesh);
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/bottomSolid/T b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/bottomSolid/T
new file mode 100644
index 0000000000..0b3a1a6b2e
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/bottomSolid/T
@@ -0,0 +1,65 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2112                                  |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    arch        "LSB;label=32;scalar=64";
+    class       volScalarField;
+    location    "0/bottomWater";
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [ 0 0 0 1 0 0 0 ];
+
+internalField   uniform 300;
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    minX
+    {
+        type            zeroGradient;
+    }
+    maxX
+    {
+        type            zeroGradient;
+    }
+    minY
+    {
+        type            uniformFixedValue;
+        uniformValue    600;
+    }
+    minZ
+    {
+        type            zeroGradient;
+    }
+    maxZ
+    {
+        type            zeroGradient;
+    }
+    group_solid
+    {
+        type            compressible::turbulentTemperatureRadCoupledMixed;
+        value           uniform 300;
+        Tnbr            T;
+        kappaMethod     solidThermo;
+    }
+    group_fluid
+    {
+        type            compressible::turbulentTemperatureRadCoupledMixed;
+        value           uniform 300;
+        Tnbr            T;
+        kappaMethod     solidThermo;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/bottomSolid/p b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/bottomSolid/p
new file mode 100644
index 0000000000..d388059d75
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/bottomSolid/p
@@ -0,0 +1,35 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2112                                  |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    arch        "LSB;label=32;scalar=64";
+    class       volScalarField;
+    location    "0/bottomWater";
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [ 1 -1 -2 0 0 0 0 ];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    ".*"
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/leftSolid/T b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/leftSolid/T
new file mode 100644
index 0000000000..e4f0a27e05
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/leftSolid/T
@@ -0,0 +1,60 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2112                                  |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    arch        "LSB;label=32;scalar=64";
+    class       volScalarField;
+    location    "0/bottomWater";
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [ 0 0 0 1 0 0 0 ];
+
+internalField   uniform 300;
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    minX
+    {
+        type            zeroGradient;
+    }
+    maxX
+    {
+        type            zeroGradient;
+    }
+    minZ
+    {
+        type            zeroGradient;
+    }
+    maxZ
+    {
+        type            zeroGradient;
+    }
+    group_solid
+    {
+        type            compressible::turbulentTemperatureRadCoupledMixed;
+        value           uniform 300;
+        Tnbr            T;
+        kappaMethod     solidThermo;
+    }
+    group_fluid
+    {
+        type            compressible::turbulentTemperatureRadCoupledMixed;
+        value           uniform 300;
+        Tnbr            T;
+        kappaMethod     solidThermo;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/leftSolid/p b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/leftSolid/p
new file mode 100644
index 0000000000..d388059d75
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/leftSolid/p
@@ -0,0 +1,35 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2112                                  |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    arch        "LSB;label=32;scalar=64";
+    class       volScalarField;
+    location    "0/bottomWater";
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [ 1 -1 -2 0 0 0 0 ];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    ".*"
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/rightSolid/T b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/rightSolid/T
new file mode 100644
index 0000000000..e4f0a27e05
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/rightSolid/T
@@ -0,0 +1,60 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2112                                  |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    arch        "LSB;label=32;scalar=64";
+    class       volScalarField;
+    location    "0/bottomWater";
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [ 0 0 0 1 0 0 0 ];
+
+internalField   uniform 300;
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    minX
+    {
+        type            zeroGradient;
+    }
+    maxX
+    {
+        type            zeroGradient;
+    }
+    minZ
+    {
+        type            zeroGradient;
+    }
+    maxZ
+    {
+        type            zeroGradient;
+    }
+    group_solid
+    {
+        type            compressible::turbulentTemperatureRadCoupledMixed;
+        value           uniform 300;
+        Tnbr            T;
+        kappaMethod     solidThermo;
+    }
+    group_fluid
+    {
+        type            compressible::turbulentTemperatureRadCoupledMixed;
+        value           uniform 300;
+        Tnbr            T;
+        kappaMethod     solidThermo;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/rightSolid/p b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/rightSolid/p
new file mode 100644
index 0000000000..d388059d75
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/rightSolid/p
@@ -0,0 +1,35 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2112                                  |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    arch        "LSB;label=32;scalar=64";
+    class       volScalarField;
+    location    "0/bottomWater";
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [ 1 -1 -2 0 0 0 0 ];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    ".*"
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/topAir/T b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/topAir/T
new file mode 100644
index 0000000000..a17739d894
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/topAir/T
@@ -0,0 +1,67 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2112                                  |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    arch        "LSB;label=32;scalar=64";
+    class       volScalarField;
+    location    "0/topAir";
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [ 0 0 0 1 0 0 0 ];
+
+internalField   uniform 300;
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    maxY
+    {
+        type            zeroGradient;
+    }
+    minX
+    {
+        type            uniformFixedValue;
+        uniformValue    300;
+    }
+    maxX
+    {
+        type            inletOutlet;
+        value           uniform 300;
+        inletValue      uniform 300;
+    }
+    minZ
+    {
+        type            zeroGradient;
+    }
+    maxZ
+    {
+        type            zeroGradient;
+    }
+    group_solid
+    {
+        type            compressible::turbulentTemperatureRadCoupledMixed;
+        value           uniform 300;
+        Tnbr            T;
+        kappaMethod     fluidThermo;
+    }
+    group_fluid
+    {
+        type            compressible::turbulentTemperatureRadCoupledMixed;
+        value           uniform 300;
+        Tnbr            T;
+        kappaMethod     fluidThermo;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/topAir/U b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/topAir/U
new file mode 100644
index 0000000000..fc6207e5eb
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/topAir/U
@@ -0,0 +1,61 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2112                                  |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    arch        "LSB;label=32;scalar=64";
+    class       volVectorField;
+    location    "0/topAir";
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [ 0 1 -1 0 0 0 0 ];
+
+internalField   uniform ( 0.1 0 0 );
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    maxY
+    {
+        type            noSlip;
+    }
+    minX
+    {
+        type            uniformFixedValue;
+        uniformValue    ( 0.1 0 0 );
+    }
+    maxX
+    {
+        type            inletOutlet;
+        value           uniform ( 0.1 0 0 );
+        inletValue      uniform ( 0 0 0 );
+    }
+    minZ
+    {
+        type            noSlip;
+    }
+    maxZ
+    {
+        type            noSlip;
+    }
+    group_solid
+    {
+        type            noSlip;
+    }
+    group_fluid
+    {
+        type            noSlip;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/topAir/epsilon b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/topAir/epsilon
new file mode 100644
index 0000000000..dda3920826
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/topAir/epsilon
@@ -0,0 +1,66 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2112                                  |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    arch        "LSB;label=32;scalar=64";
+    class       volScalarField;
+    location    "0/topAir";
+    object      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [ 0 2 -3 0 0 0 0 ];
+
+internalField   uniform 0.01;
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    maxY
+    {
+        type            epsilonWallFunction;
+        value           uniform 0.01;
+    }
+    minX
+    {
+        type            uniformFixedValue;
+        uniformValue    0.01;
+    }
+    maxX
+    {
+        type            inletOutlet;
+        value           uniform 0.01;
+        inletValue      uniform 0.01;
+    }
+    minZ
+    {
+        type            epsilonWallFunction;
+        value           uniform 0.01;
+    }
+    maxZ
+    {
+        type            epsilonWallFunction;
+        value           uniform 0.01;
+    }
+    group_solid
+    {
+        type            epsilonWallFunction;
+        value           uniform 0.01;
+    }
+    group_fluid
+    {
+        type            epsilonWallFunction;
+        value           uniform 0.01;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/topAir/k b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/topAir/k
new file mode 100644
index 0000000000..d148990adf
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/topAir/k
@@ -0,0 +1,66 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2112                                  |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    arch        "LSB;label=32;scalar=64";
+    class       volScalarField;
+    location    "0/topAir";
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [ 0 2 -2 0 0 0 0 ];
+
+internalField   uniform 0.1;
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    maxY
+    {
+        type            kqRWallFunction;
+        value           uniform 0.1;
+    }
+    minX
+    {
+        type            uniformFixedValue;
+        uniformValue    0.1;
+    }
+    maxX
+    {
+        type            inletOutlet;
+        value           uniform 0.1;
+        inletValue      uniform 0.1;
+    }
+    minZ
+    {
+        type            kqRWallFunction;
+        value           uniform 0.1;
+    }
+    maxZ
+    {
+        type            kqRWallFunction;
+        value           uniform 0.1;
+    }
+    group_solid
+    {
+        type            kqRWallFunction;
+        value           uniform 0.1;
+    }
+    group_fluid
+    {
+        type            kqRWallFunction;
+        value           uniform 0.1;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/topAir/p b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/topAir/p
new file mode 100644
index 0000000000..6ea2201b22
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/topAir/p
@@ -0,0 +1,65 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2112                                  |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    arch        "LSB;label=32;scalar=64";
+    class       volScalarField;
+    location    "0/topAir";
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [ 1 -1 -2 0 0 0 0 ];
+
+internalField   uniform 100000;
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    maxY
+    {
+        type            calculated;
+        value           uniform 100000;
+    }
+    minX
+    {
+        type            calculated;
+        value           uniform 100000;
+    }
+    maxX
+    {
+        type            calculated;
+        value           uniform 100000;
+    }
+    minZ
+    {
+        type            calculated;
+        value           uniform 100000;
+    }
+    maxZ
+    {
+        type            calculated;
+        value           uniform 100000;
+    }
+    group_solid
+    {
+        type            calculated;
+        value           uniform 100000;
+    }
+    group_fluid
+    {
+        type            calculated;
+        value           uniform 100000;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/topAir/p_rgh b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/topAir/p_rgh
new file mode 100644
index 0000000000..b559a63780
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/0.orig/topAir/p_rgh
@@ -0,0 +1,65 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2112                                  |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    arch        "LSB;label=32;scalar=64";
+    class       volScalarField;
+    location    "0/topAir";
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [ 1 -1 -2 0 0 0 0 ];
+
+internalField   uniform 100000;
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    maxY
+    {
+        type            fixedFluxPressure;
+        value           uniform 100000;
+    }
+    minX
+    {
+        type            fixedFluxPressure;
+        value           uniform 100000;
+    }
+    maxX
+    {
+        type            uniformFixedValue;
+        uniformValue    100000;
+    }
+    minZ
+    {
+        type            fixedFluxPressure;
+        value           uniform 100000;
+    }
+    maxZ
+    {
+        type            fixedFluxPressure;
+        value           uniform 100000;
+    }
+    group_solid
+    {
+        type            fixedFluxPressure;
+        value           uniform 100000;
+    }
+    group_fluid
+    {
+        type            fixedFluxPressure;
+        value           uniform 100000;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/Allclean b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/Allclean
new file mode 100755
index 0000000000..67390b667b
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/Allclean
@@ -0,0 +1,12 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/CleanFunctions      # Tutorial clean functions
+#------------------------------------------------------------------------------
+
+cleanCase0
+rm -r constant/bottomSolid/polyMesh
+rm -r constant/leftSolid/polyMesh
+rm -r constant/rightSolid/polyMesh
+rm -r constant/topAir/polyMesh
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/Allrun b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/Allrun
new file mode 100755
index 0000000000..d5d83a85a5
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/Allrun
@@ -0,0 +1,20 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+runApplication ./Allrun.pre
+
+#-- Run on single processor
+#runApplication $(getApplication)
+
+# Decompose
+runApplication decomposePar -allRegions
+
+# Run
+runParallel $(getApplication)
+
+# Reconstruct
+runApplication reconstructPar -allRegions
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/Allrun.pre b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/Allrun.pre
new file mode 100755
index 0000000000..2c5eef5c8f
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/Allrun.pre
@@ -0,0 +1,25 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+for region in $(foamListRegions)
+do
+    runApplication -s $region blockMesh -region $region
+done
+
+#- createPatch creates patchFields out of dictionary only but only then
+#  moves faces to it. Hence most of its state (refValue, valueFraction etc)
+#  will not be sized. Make sure to not to fail construction so disable
+#  floating point error checking
+#unset FOAM_SIGFPE
+runApplication createPatch -allRegions -overwrite
+
+restore0Dir
+
+echo
+echo "Use paraFoam -touch-all to create files for paraview post-processing"
+echo
+echo "End"
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/README.txt b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/README.txt
new file mode 100644
index 0000000000..34427f116a
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/README.txt
@@ -0,0 +1,4 @@
+2022-02-24: Testcase for:
+- createPatch with autoPatch functionality
+- AMIMethod limiting the search distance
+- use of patch groups to ease set-up
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/bottomSolid/radiationProperties b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/bottomSolid/radiationProperties
new file mode 100644
index 0000000000..145878911f
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/bottomSolid/radiationProperties
@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2112                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      radiationProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+radiation off;
+
+radiationModel  none;
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/bottomSolid/thermophysicalProperties b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/bottomSolid/thermophysicalProperties
new file mode 100644
index 0000000000..eb1e7407c5
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/bottomSolid/thermophysicalProperties
@@ -0,0 +1,50 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2112                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heSolidThermo;
+    mixture         pureMixture;
+    transport       constIso;
+    thermo          hConst;
+    equationOfState rhoConst;
+    specie          specie;
+    energy          sensibleEnthalpy;
+}
+
+mixture
+{
+    specie
+    {
+        molWeight   50;
+    }
+    transport
+    {
+        kappa   80;
+    }
+    thermodynamics
+    {
+        Hf      0;
+        Cp      450;
+    }
+    equationOfState
+    {
+        rho     8000;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/g b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/g
new file mode 100644
index 0000000000..74e8297361
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/g
@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2112                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       uniformDimensionedVectorField;
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           (0 -9.81 0);
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/leftSolid/radiationProperties b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/leftSolid/radiationProperties
new file mode 120000
index 0000000000..e5151e317c
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/leftSolid/radiationProperties
@@ -0,0 +1 @@
+../bottomSolid/radiationProperties
\ No newline at end of file
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/leftSolid/thermophysicalProperties b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/leftSolid/thermophysicalProperties
new file mode 120000
index 0000000000..9dfe396a89
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/leftSolid/thermophysicalProperties
@@ -0,0 +1 @@
+../bottomSolid/thermophysicalProperties
\ No newline at end of file
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/regionProperties b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/regionProperties
new file mode 100644
index 0000000000..6c1281778e
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/regionProperties
@@ -0,0 +1,24 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2112                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      regionProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+regions
+(
+    fluid       (topAir)
+    solid       (bottomSolid leftSolid rightSolid)
+);
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/rightSolid/radiationProperties b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/rightSolid/radiationProperties
new file mode 120000
index 0000000000..e5151e317c
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/rightSolid/radiationProperties
@@ -0,0 +1 @@
+../bottomSolid/radiationProperties
\ No newline at end of file
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/rightSolid/thermophysicalProperties b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/rightSolid/thermophysicalProperties
new file mode 120000
index 0000000000..9dfe396a89
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/rightSolid/thermophysicalProperties
@@ -0,0 +1 @@
+../bottomSolid/thermophysicalProperties
\ No newline at end of file
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/topAir/radiationProperties b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/topAir/radiationProperties
new file mode 100644
index 0000000000..12d23ce2ac
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/topAir/radiationProperties
@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2112                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      radiationProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+radiation       off;
+
+radiationModel  none;
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/topAir/thermophysicalProperties b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/topAir/thermophysicalProperties
new file mode 100644
index 0000000000..9533066d7e
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/topAir/thermophysicalProperties
@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2112                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         pureMixture;
+    transport       const;
+    thermo          hConst;
+    equationOfState perfectGas;
+    specie          specie;
+    energy          sensibleEnthalpy;
+}
+
+mixture
+{
+    specie
+    {
+        molWeight       28.9;
+    }
+    thermodynamics
+    {
+        Cp              1000;
+        Hf              0;
+    }
+    transport
+    {
+        mu              1.8e-05;
+        Pr              0.7;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/topAir/turbulenceProperties b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/topAir/turbulenceProperties
new file mode 100644
index 0000000000..8e51de51a2
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/constant/topAir/turbulenceProperties
@@ -0,0 +1,20 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2112                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  laminar;
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/README b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/README
new file mode 100644
index 0000000000..f75a0ab534
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/README
@@ -0,0 +1,2 @@
+- fvSolution is used for outer correctors specification.
+- fvSchemes is only so that pre-processing activities can proceed
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/bottomSolid/blockMeshDict b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/bottomSolid/blockMeshDict
new file mode 100644
index 0000000000..aa89241a01
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/bottomSolid/blockMeshDict
@@ -0,0 +1,93 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2112                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scale   1;
+
+vertices
+(
+    (-0.1 -0.04  -0.05)
+    ( 0.1 -0.04  -0.05)
+    ( 0.1  0.00  -0.05)
+    (-0.1  0.00  -0.05)
+    (-0.1 -0.04   0.05)
+    ( 0.1 -0.04   0.05)
+    ( 0.1  0.00   0.05)
+    (-0.1  0.00   0.05)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (10 4 1) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    minX
+    {
+        type patch;
+        faces
+        (
+            (0 4 7 3)
+        );
+    }
+    maxX
+    {
+        type patch;
+        faces
+        (
+            (2 6 5 1)
+        );
+    }
+    minY
+    {
+        type patch;
+        faces
+        (
+            (1 5 4 0)
+        );
+    }
+    maxY
+    {
+        type patch;
+        inGroups    (group_solid);
+        faces
+        (
+            (3 7 6 2)
+        );
+    }
+    minZ
+    {
+        type patch;
+        faces
+        (
+            (0 3 2 1)
+        );
+    }
+    maxZ
+    {
+        type patch;
+        faces
+        (
+            (4 5 6 7)
+        );
+    }
+);
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/bottomSolid/createPatchDict b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/bottomSolid/createPatchDict
new file mode 100644
index 0000000000..71b1120950
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/bottomSolid/createPatchDict
@@ -0,0 +1,163 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2106                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      createPatchDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// This application/dictionary controls:
+// - optional: create new patches from boundary faces (either given as
+//   a set of patches or as a faceSet)
+// - always: order faces on coupled patches such that they are opposite. This
+//   is done for all coupled faces, not just for any patches created.
+// - optional: synchronise points on coupled patches.
+// - always: remove zero-sized (non-coupled) patches (that were not added)
+
+// 1. Create cyclic:
+// - specify where the faces should come from
+// - specify the type of cyclic. If a rotational specify the rotationAxis
+//   and centre to make matching easier
+// - always create both halves in one invocation with correct 'neighbourPatch'
+//   setting.
+// - optionally pointSync true to guarantee points to line up.
+
+// 2. Correct incorrect cyclic:
+// This will usually fail upon loading:
+//  "face 0 area does not match neighbour 2 by 0.0100005%"
+//  " -- possible face ordering problem."
+// - in polyMesh/boundary file:
+//      - loosen matchTolerance of all cyclics to get case to load
+//      - or change patch type from 'cyclic' to 'patch'
+//        and regenerate cyclic as above
+
+// Do a synchronisation of coupled points after creation of any patches.
+// Note: this does not work with points that are on multiple coupled patches
+//       with transformations (i.e. cyclics).
+pointSync false;
+
+// Patches to create.
+patches
+(
+    // Patches to couple to solids
+    {
+        // Dictionary to construct new patch from
+        patchInfo
+        {
+            type            mappedPatch;
+            inGroups        (group_solid);
+            sampleMode      nearestPatchFaceAMI;
+
+            AMIMethod       faceAreaWeightAMI;
+            maxDistance2    1e-6;   // max (square) distance
+            minCosAngle     0.5;    // min alignment
+
+            // Overwritten
+            //sampleRegion    otherRegion;
+            //samplePatch     otherPatch;
+
+            //- Optional override of added patchfields. If not specified
+            //  any added patchfields are of type calculated.
+            patchFields
+            {
+                ////- Problem is with patch fields that need additional
+                ////  state. Patches get created with 0 faces and
+                ////  then mapped to the correct size. The autoMap generally
+                ////  does not know value to give the new faces.
+                //T
+                //{
+                //    type    compressible::turbulentTemperatureRadCoupledMixed;
+                //    Tnbr            T;
+                //    kappaMethod     solidThermo;
+                //    value           uniform 300;
+                //    refValue        uniform 300;
+                //    valueFraction   uniform 1;
+                //    refGradient     uniform 0;
+                //}
+                //p
+                //{
+                //    type    calculated;
+                //    value   uniform 100000;
+                //}
+            }
+        }
+
+        // How to select the faces:
+        //  - set : specify faceSet in 'set'
+        //  - patches : specify names in 'patches'
+        //  - autoPatch : attempts automatic patching of the specified
+        //                candidates in 'patches'.
+        //      - single region : match in the region itself
+        //      - multi regions : match in between regions only
+        constructFrom autoPatch;
+
+        // If constructFrom = patches or autoPatch: names of patches.
+        // Wildcards&patchGroups allowed.
+        patches (group_solid);
+    }
+
+
+    // Patches to couple to fluids
+    {
+        // Dictionary to construct new patch from
+        patchInfo
+        {
+            type            mappedPatch;
+            inGroups        (group_fluid);
+            sampleMode      nearestPatchFaceAMI;
+
+            AMIMethod       faceAreaWeightAMI;
+            maxDistance2    1e-6;   // max (square) distance
+            minCosAngle     0.5;    // min alignment
+
+            // Overwritten
+            //sampleRegion    otherRegion;
+            //samplePatch     otherPatch;
+
+            //- Optional override of added patchfields. If not specified
+            //  any added patchfields are of type calculated.
+            patchFields
+            {
+                //- Problem is with patch fields that need additional
+                //  state. Patches get created with 0 faces and
+                //  then mapped to the correct size. The autoMap genrally
+                //  does not know value to give the new faces.
+                //T
+                //{
+                //    type    compressible::turbulentTemperatureRadCoupledMixed;
+                //    Tnbr    T;
+                //    kappaMethod solidThermo;
+                //    value   uniform 300;
+                //}
+                //p
+                //{
+                //    type    calculated;
+                //    value   uniform 100000;
+                //}
+            }
+        }
+
+        // How to select the faces:
+        //  - set : specify faceSet in 'set'
+        //  - patches : specify names in 'patches'
+        //  - autoPatch : attempts automatic patching of the specified
+        //                candidates in 'patches'.
+        //      - single region : match in the region itself
+        //      - multi regions : match in between regions only
+        constructFrom autoPatch;
+
+        // If constructFrom = patches or autoPatch: names of patches.
+        // Wildcards&patchGroups allowed.
+        patches (group_fluid);
+    }
+);
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/bottomSolid/fvSchemes b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/bottomSolid/fvSchemes
new file mode 100644
index 0000000000..083ce3c9b2
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/bottomSolid/fvSchemes
@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2112                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     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;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/bottomSolid/fvSolution b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/bottomSolid/fvSolution
new file mode 100644
index 0000000000..c792128442
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/bottomSolid/fvSolution
@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2112                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     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/mesh/createPatch/multiRegionHeater_autoPatch/system/controlDict b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/controlDict
new file mode 100644
index 0000000000..cb173aef7d
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/controlDict
@@ -0,0 +1,60 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2112                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     chtMultiRegionFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         100;
+
+deltaT          0.001;
+
+writeControl    adjustable;
+
+writeInterval   10;
+
+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
+{
+    #include "vtkWrite"
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/fvSchemes b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/fvSchemes
new file mode 100644
index 0000000000..ac3cd1a53d
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/fvSchemes
@@ -0,0 +1,42 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2112                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+}
+
+gradSchemes
+{
+}
+
+divSchemes
+{
+}
+
+laplacianSchemes
+{
+}
+
+interpolationSchemes
+{
+}
+
+snGradSchemes
+{
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/fvSolution b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/fvSolution
new file mode 100644
index 0000000000..4c9956b7cd
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/fvSolution
@@ -0,0 +1,23 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2112                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+PIMPLE
+{
+    nOuterCorrectors 1;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/leftSolid/blockMeshDict b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/leftSolid/blockMeshDict
new file mode 100644
index 0000000000..6a84cecab6
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/leftSolid/blockMeshDict
@@ -0,0 +1,95 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2112                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scale   1;
+
+vertices
+(
+    (-0.1  0.00  -0.05)
+    ( 0.0  0.00  -0.05)
+    ( 0.0  0.01  -0.05)
+    (-0.1  0.01  -0.05)
+    (-0.1  0.00   0.05)
+    ( 0.0  0.00   0.05)
+    ( 0.0  0.01   0.05)
+    (-0.1  0.01   0.05)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (10 4 1) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    minX
+    {
+        type patch;
+        faces
+        (
+            (0 4 7 3)
+        );
+    }
+    maxX
+    {
+        type patch;
+        inGroups    (group_solid);
+        faces
+        (
+            (2 6 5 1)
+        );
+    }
+    minY
+    {
+        type patch;
+        inGroups    (group_solid);
+        faces
+        (
+            (1 5 4 0)
+        );
+    }
+    maxY
+    {
+        type patch;
+        inGroups    (group_fluid);
+        faces
+        (
+            (3 7 6 2)
+        );
+    }
+    minZ
+    {
+        type patch;
+        faces
+        (
+            (0 3 2 1)
+        );
+    }
+    maxZ
+    {
+        type patch;
+        faces
+        (
+            (4 5 6 7)
+        );
+    }
+);
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/leftSolid/createPatchDict b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/leftSolid/createPatchDict
new file mode 120000
index 0000000000..4a18f41d0d
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/leftSolid/createPatchDict
@@ -0,0 +1 @@
+../bottomSolid/createPatchDict
\ No newline at end of file
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/leftSolid/fvSchemes b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/leftSolid/fvSchemes
new file mode 120000
index 0000000000..7630faab04
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/leftSolid/fvSchemes
@@ -0,0 +1 @@
+../bottomSolid/fvSchemes
\ No newline at end of file
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/leftSolid/fvSolution b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/leftSolid/fvSolution
new file mode 120000
index 0000000000..4661656087
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/leftSolid/fvSolution
@@ -0,0 +1 @@
+../bottomSolid/fvSolution
\ No newline at end of file
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/rightSolid/blockMeshDict b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/rightSolid/blockMeshDict
new file mode 100644
index 0000000000..87fe32a411
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/rightSolid/blockMeshDict
@@ -0,0 +1,95 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2112                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scale   1;
+
+vertices
+(
+    (-0.0  0.00  -0.05)
+    ( 0.1  0.00  -0.05)
+    ( 0.1  0.01  -0.05)
+    (-0.0  0.01  -0.05)
+    (-0.0  0.00   0.05)
+    ( 0.1  0.00   0.05)
+    ( 0.1  0.01   0.05)
+    (-0.0  0.01   0.05)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (10 4 1) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    minX
+    {
+        type patch;
+        inGroups    (group_solid);
+        faces
+        (
+            (0 4 7 3)
+        );
+    }
+    maxX
+    {
+        type patch;
+        faces
+        (
+            (2 6 5 1)
+        );
+    }
+    minY
+    {
+        type patch;
+        inGroups    (group_solid);
+        faces
+        (
+            (1 5 4 0)
+        );
+    }
+    maxY
+    {
+        type patch;
+        inGroups    (group_fluid);
+        faces
+        (
+            (3 7 6 2)
+        );
+    }
+    minZ
+    {
+        type patch;
+        faces
+        (
+            (0 3 2 1)
+        );
+    }
+    maxZ
+    {
+        type patch;
+        faces
+        (
+            (4 5 6 7)
+        );
+    }
+);
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/rightSolid/createPatchDict b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/rightSolid/createPatchDict
new file mode 120000
index 0000000000..4a18f41d0d
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/rightSolid/createPatchDict
@@ -0,0 +1 @@
+../bottomSolid/createPatchDict
\ No newline at end of file
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/rightSolid/fvSchemes b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/rightSolid/fvSchemes
new file mode 120000
index 0000000000..7630faab04
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/rightSolid/fvSchemes
@@ -0,0 +1 @@
+../bottomSolid/fvSchemes
\ No newline at end of file
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/rightSolid/fvSolution b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/rightSolid/fvSolution
new file mode 120000
index 0000000000..4661656087
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/rightSolid/fvSolution
@@ -0,0 +1 @@
+../bottomSolid/fvSolution
\ No newline at end of file
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/topAir/blockMeshDict b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/topAir/blockMeshDict
new file mode 100644
index 0000000000..a2afc7f431
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/topAir/blockMeshDict
@@ -0,0 +1,93 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2112                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scale   1;
+
+vertices
+(
+    (-0.1  0.01  -0.05)
+    ( 0.1  0.01  -0.05)
+    ( 0.1  0.04  -0.05)
+    (-0.1  0.04  -0.05)
+    (-0.1  0.01   0.05)
+    ( 0.1  0.01   0.05)
+    ( 0.1  0.04   0.05)
+    (-0.1  0.04   0.05)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (20 4 10) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    minX
+    {
+        type patch;
+        faces
+        (
+            (0 4 7 3)
+        );
+    }
+    maxX
+    {
+        type patch;
+        faces
+        (
+            (2 6 5 1)
+        );
+    }
+    minY
+    {
+        type patch;
+        inGroups    (group_fluid);
+        faces
+        (
+            (1 5 4 0)
+        );
+    }
+    maxY
+    {
+        type patch;
+        faces
+        (
+            (3 7 6 2)
+        );
+    }
+    minZ
+    {
+        type patch;
+        faces
+        (
+            (0 3 2 1)
+        );
+    }
+    maxZ
+    {
+        type patch;
+        faces
+        (
+            (4 5 6 7)
+        );
+    }
+);
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/topAir/createPatchDict b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/topAir/createPatchDict
new file mode 100644
index 0000000000..fbfb281279
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/topAir/createPatchDict
@@ -0,0 +1,99 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2106                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      createPatchDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// This application/dictionary controls:
+// - optional: create new patches from boundary faces (either given as
+//   a set of patches or as a faceSet)
+// - always: order faces on coupled patches such that they are opposite. This
+//   is done for all coupled faces, not just for any patches created.
+// - optional: synchronise points on coupled patches.
+// - always: remove zero-sized (non-coupled) patches (that were not added)
+
+// 1. Create cyclic:
+// - specify where the faces should come from
+// - specify the type of cyclic. If a rotational specify the rotationAxis
+//   and centre to make matching easier
+// - always create both halves in one invocation with correct 'neighbourPatch'
+//   setting.
+// - optionally pointSync true to guarantee points to line up.
+
+// 2. Correct incorrect cyclic:
+// This will usually fail upon loading:
+//  "face 0 area does not match neighbour 2 by 0.0100005%"
+//  " -- possible face ordering problem."
+// - in polyMesh/boundary file:
+//      - loosen matchTolerance of all cyclics to get case to load
+//      - or change patch type from 'cyclic' to 'patch'
+//        and regenerate cyclic as above
+
+// Do a synchronisation of coupled points after creation of any patches.
+// Note: this does not work with points that are on multiple coupled patches
+//       with transformations (i.e. cyclics).
+pointSync false;
+
+// Patches to create.
+patches
+(
+    {
+        // Dictionary to construct new patch from
+        patchInfo
+        {
+            type            mappedPatch;
+            inGroups        (group_fluid);
+            sampleMode      nearestPatchFaceAMI;
+
+            AMIMethod       faceAreaWeightAMI;
+            maxDistance2    1e-6;   // max (square) distance
+            minCosAngle     0.5;    // min alignment
+
+            // Overwritten
+            //sampleRegion    otherRegion;
+            //samplePatch     otherPatch;
+
+            //- Optional override of added patchfields. If not specified
+            //  any added patchfields are of type calculated.
+            patchFields
+            {
+                //- Problem is with patch fields that need additional
+                //  state. Patches get created with 0 faces and
+                //  then mapped to the correct size. The autoMap genrally
+                //  does not know value to give the new faces.
+                //T
+                //{
+                //    type    compressible::turbulentTemperatureRadCoupledMixed;
+                //    Tnbr    T;
+                //    kappaMethod solidThermo;
+                //    value   uniform 300;
+                //}
+            }
+        }
+
+        // How to select the faces:
+        //  - set : specify faceSet in 'set'
+        //  - patches : specify names in 'patches'
+        //  - autoPatch : attempts automatic patching of the specified
+        //                candidates in 'patches'.
+        //      - single region : match in the region itself
+        //      - multi regions : match in between regions only
+        constructFrom autoPatch;
+
+        // If constructFrom = patches or autoPatch: names of patches.
+        // Wildcards&patchGroups allowed.
+        patches (group_fluid);
+    }
+);
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/topAir/fvSchemes b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/topAir/fvSchemes
new file mode 100644
index 0000000000..4726d5cdac
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/topAir/fvSchemes
@@ -0,0 +1,61 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2112                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+
+    div(phi,U)      Gauss upwind;
+
+    div(phi,K)      Gauss linear;
+    div(phi,h)      Gauss upwind;
+
+    turbulence      Gauss upwind;
+    div(phi,k)      $turbulence;
+    div(phi,epsilon) $turbulence;
+    div(phi,R)      $turbulence;
+    div(R)          Gauss linear;
+
+    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default        Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/topAir/fvSolution b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/topAir/fvSolution
new file mode 100644
index 0000000000..61cf8131e2
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/topAir/fvSolution
@@ -0,0 +1,82 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2112                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    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         GaussSeidel;
+    }
+
+    p_rghFinal
+    {
+        $p_rgh;
+        tolerance        1e-7;
+        relTol           0;
+    }
+
+    "(U|h|k|epsilon|R)"
+    {
+        solver           PBiCGStab;
+        preconditioner   DILU;
+        tolerance        1e-7;
+        relTol           0.1;
+    }
+
+    "(U|h|k|epsilon|R)Final"
+    {
+        $U;
+        tolerance        1e-7;
+        relTol           0;
+    }
+}
+
+PIMPLE
+{
+    momentumPredictor   yes;
+    nCorrectors         2;
+    nNonOrthogonalCorrectors 0;
+}
+
+relaxationFactors
+{
+    equations
+    {
+        "h.*"           1;
+        "U.*"           1;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/vtkWrite b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/vtkWrite
new file mode 100644
index 0000000000..b9b5dd2301
--- /dev/null
+++ b/tutorials/mesh/createPatch/multiRegionHeater_autoPatch/system/vtkWrite
@@ -0,0 +1,61 @@
+// -*- C++ -*-
+// Use the vtkWrite function object
+
+vtkWrite
+{
+    type    vtkWrite;
+    libs    (utilityFunctionObjects);
+    log     true;
+
+    writeControl    writeTime;
+    writeInterval   1;
+
+    regions (".*");
+
+    internal    true;
+
+    boundary    true;
+
+    single      false;
+
+    interpolate true;
+
+    // Fields to output (words or regex)
+    fields  (".*");
+
+    //- Output format (ascii | binary) - Default=binary
+    // format  binary;
+
+    //- Use legacy output format - Default=false
+    // legacy  false;
+
+    //- Output directory name - Default="postProcessing/<name>"
+    // directory   "VTK";
+
+    //- Write cell ids as field - Default=true
+    writeIds    false;
+}
+
+
+// Solid walls only
+walls
+{
+    type    vtkWrite;
+    libs    (utilityFunctionObjects);
+    log     true;
+
+    writeControl    writeTime;
+    writeInterval   1;
+
+    internal    false;
+
+    // single      true;
+
+    regions     ( heater "(?i).*solid" );
+    patches     ( "(?i).*solid_to.*" "heater.*(Air|Water)" );
+
+    fields      (T);
+}
+
+
+// ************************************************************************* //