From 8bb44c7410f352a0052f7a7169b98c0779070153 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Wed, 19 Oct 2011 18:25:00 +0100
Subject: [PATCH] ENH: extrudeToRegionMesh: parallel operation

---
 .../extrude/extrudeToRegionMesh/Make/files    |    1 -
 .../extrudeToRegionMesh/createShellMesh.C     |  415 ++-
 .../extrudeToRegionMesh/createShellMesh.H     |   29 +-
 .../extrudeToRegionMesh/extrudeToRegionMesh.C | 2218 ++++++++++-------
 .../patchPointEdgeCirculator.C                |   67 -
 .../patchPointEdgeCirculator.H                |  205 --
 .../patchPointEdgeCirculatorI.H               |  308 ---
 7 files changed, 1733 insertions(+), 1510 deletions(-)
 delete mode 100644 applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/patchPointEdgeCirculator.C
 delete mode 100644 applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/patchPointEdgeCirculator.H
 delete mode 100644 applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/patchPointEdgeCirculatorI.H

diff --git a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/Make/files b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/Make/files
index be2c8d6975..28ac7de544 100644
--- a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/Make/files
+++ b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/Make/files
@@ -1,4 +1,3 @@
-patchPointEdgeCirculator.C
 createShellMesh.C
 extrudeToRegionMesh.C
 
diff --git a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh.C b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh.C
index 935b30a5bd..6bfe6ca62e 100644
--- a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh.C
+++ b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh.C
@@ -31,130 +31,381 @@ License
 #include "polyAddFace.H"
 #include "polyModifyFace.H"
 #include "polyAddCell.H"
-#include "patchPointEdgeCirculator.H"
+#include "labelPair.H"
+#include "indirectPrimitivePatch.H"
+#include "mapDistribute.H"
+#include "globalMeshData.H"
+#include "PatchTools.H"
+#include "globalIndex.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 defineTypeNameAndDebug(Foam::createShellMesh, 0);
 
+namespace Foam
+{
+template<>
+class minEqOp<labelPair>
+{
+public:
+    void operator()(labelPair& x, const labelPair& y) const
+    {
+        x[0] = min(x[0], y[0]);
+        x[1] = min(x[1], y[1]);
+    }
+};
+}
+
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::createShellMesh::calcPointRegions
+// Synchronise edges
+void Foam::createShellMesh::syncEdges
 (
-    const primitiveFacePatch& patch,
-    const PackedBoolList& nonManifoldEdge,
-    faceList& pointRegions,
-    labelList& regionPoints
+    const globalMeshData& globalData,
+
+    const labelList& patchEdges,
+    const labelList& coupledEdges,
+    const PackedBoolList& sameEdgeOrientation,
+
+    PackedBoolList& isChangedEdge,
+    DynamicList<label>& changedEdges,
+    labelPairList& allEdgeData
 )
 {
-    pointRegions.setSize(patch.size());
-    forAll(pointRegions, faceI)
+    const mapDistribute& map = globalData.globalEdgeSlavesMap();
+    const PackedBoolList& cppOrientation = globalData.globalEdgeOrientation();
+
+    // Convert patch-edge data into cpp-edge data
+    labelPairList cppEdgeData
+    (
+        map.constructSize(),
+        labelPair(labelMax, labelMax)
+    );
+
+    forAll(patchEdges, i)
     {
-        const face& f = patch.localFaces()[faceI];
-        pointRegions[faceI].setSize(f.size(), -1);
+        label patchEdgeI = patchEdges[i];
+        label coupledEdgeI = coupledEdges[i];
+
+        if (isChangedEdge[patchEdgeI])
+        {
+            const labelPair& data = allEdgeData[patchEdgeI];
+
+            // Patch-edge data needs to be converted into coupled-edge data
+            // (optionally flipped) and consistent in orientation with
+            // other coupled edge (optionally flipped)
+            if (sameEdgeOrientation[i] == cppOrientation[coupledEdgeI])
+            {
+                cppEdgeData[coupledEdgeI] = data;
+            }
+            else
+            {
+                cppEdgeData[coupledEdgeI] = labelPair(data[1], data[0]);
+            }
+        }
     }
 
-    label nRegions = 0;
+    // Synchronise
+    globalData.syncData
+    (
+        cppEdgeData,
+        globalData.globalEdgeSlaves(),
+        globalData.globalEdgeTransformedSlaves(),
+        map,
+        minEqOp<labelPair>()
+    );
 
-    forAll(pointRegions, faceI)
+    // Back from cpp-edge to patch-edge data
+    forAll(patchEdges, i)
+    {
+        label patchEdgeI = patchEdges[i];
+        label coupledEdgeI = coupledEdges[i];
+
+        if (cppEdgeData[coupledEdgeI] != labelPair(labelMax, labelMax))
+        {
+            const labelPair& data = cppEdgeData[coupledEdgeI];
+
+            if (sameEdgeOrientation[i] == cppOrientation[coupledEdgeI])
+            {
+                allEdgeData[patchEdgeI] = data;
+            }
+            else
+            {
+                allEdgeData[patchEdgeI] = labelPair(data[1], data[0]);
+            }
+
+            if (!isChangedEdge[patchEdgeI])
+            {
+                changedEdges.append(patchEdgeI);
+                isChangedEdge[patchEdgeI] = true;
+            }
+        }
+    }
+}
+
+
+void Foam::createShellMesh::calcPointRegions
+(
+    const globalMeshData& globalData,
+    const primitiveFacePatch& patch,
+    const PackedBoolList& nonManifoldEdge,
+
+    faceList& pointGlobalRegions,
+    faceList& pointLocalRegions,
+    labelList& localToGlobalRegion
+)
+{
+    const indirectPrimitivePatch& cpp = globalData.coupledPatch();
+
+    // Calculate correspondence between patch and globalData.coupledPatch.
+    labelList patchEdges;
+    labelList coupledEdges;
+    PackedBoolList sameEdgeOrientation;
+    PatchTools::matchEdges
+    (
+        cpp,
+        patch,
+
+        coupledEdges,
+        patchEdges,
+        sameEdgeOrientation
+    );
+
+
+    // Initial unique regions
+    // ~~~~~~~~~~~~~~~~~~~~~~
+    // These get merged later on across connected edges.
+
+    // 1. Count
+    label nMaxRegions = 0;
+    forAll(patch.localFaces(), faceI)
     {
         const face& f = patch.localFaces()[faceI];
+        nMaxRegions += f.size();
+    }
 
-        forAll(f, fp)
+    const globalIndex globalRegions(nMaxRegions);
+
+    // 2. Assign unique regions
+    label nRegions = 0;
+
+    pointGlobalRegions.setSize(patch.size());
+    forAll(pointGlobalRegions, faceI)
+    {
+        const face& f = patch.localFaces()[faceI];
+        labelList& pRegions = pointGlobalRegions[faceI];
+        pRegions.setSize(f.size());
+        forAll(pRegions, fp)
         {
-            if (pointRegions[faceI][fp] == -1)
+            pRegions[fp] = globalRegions.toGlobal(nRegions++);
+        }
+    }
+
+
+    DynamicList<label> changedEdges(patch.nEdges());
+    labelPairList allEdgeData(patch.nEdges(), labelPair(labelMax, labelMax));
+    PackedBoolList isChangedEdge(patch.nEdges());
+
+
+    // Fill initial seed
+    // ~~~~~~~~~~~~~~~~~
+
+    forAll(patch.edgeFaces(), edgeI)
+    {
+        if (!nonManifoldEdge[edgeI])
+        {
+            // Take over value from one face only.
+            const edge& e = patch.edges()[edgeI];
+            label faceI = patch.edgeFaces()[edgeI][0];
+            const face& f = patch.localFaces()[faceI];
+
+            label fp0 = findIndex(f, e[0]);
+            label fp1 = findIndex(f, e[1]);
+            allEdgeData[edgeI] = labelPair
+            (
+                pointGlobalRegions[faceI][fp0],
+                pointGlobalRegions[faceI][fp1]
+            );
+            if (!isChangedEdge[edgeI])
             {
-                // Found unassigned point. Distribute current region.
-                label pointI = f[fp];
-                label edgeI = patch.faceEdges()[faceI][fp];
+                changedEdges.append(edgeI);
+                isChangedEdge[edgeI] = true;
+            }
+        }
+    }
 
-                patchPointEdgeCirculator circ
-                (
-                    patch,
-                    nonManifoldEdge,
-                    edgeI,
-                    findIndex(patch.edgeFaces()[edgeI], faceI),
-                    pointI
-                );
 
-                for
-                (
-                    patchPointEdgeCirculator iter = circ.begin();
-                    iter != circ.end();
-                    ++iter
-                )
+    syncEdges
+    (
+        globalData,
+
+        patchEdges,
+        coupledEdges,
+        sameEdgeOrientation,
+
+        isChangedEdge,
+        changedEdges,
+        allEdgeData
+    );
+
+
+    // Edge-Face-Edge walk across patch
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Across edge minimum regions win
+
+    while (true)
+    {
+        // From edge to face
+        // ~~~~~~~~~~~~~~~~~
+
+        DynamicList<label> changedFaces(patch.size());
+        PackedBoolList isChangedFace(patch.size());
+
+        forAll(changedEdges, changedI)
+        {
+            label edgeI = changedEdges[changedI];
+            const labelPair& edgeData = allEdgeData[edgeI];
+
+            const edge& e = patch.edges()[edgeI];
+            const labelList& eFaces = patch.edgeFaces()[edgeI];
+
+            forAll(eFaces, i)
+            {
+                label faceI = eFaces[i];
+                const face& f = patch.localFaces()[faceI];
+
+                // Combine edgeData with face data
+                label fp0 = findIndex(f, e[0]);
+                if (pointGlobalRegions[faceI][fp0] > edgeData[0])
                 {
-                    label face2 = iter.faceID();
-
-                    if (face2 != -1)
+                    pointGlobalRegions[faceI][fp0] = edgeData[0];
+                    if (!isChangedFace[faceI])
                     {
-                        const face& f2 = patch.localFaces()[face2];
-                        label fp2 = findIndex(f2, pointI);
-                        label& region = pointRegions[face2][fp2];
-                        if (region != -1)
-                        {
-                            FatalErrorIn
-                            (
-                                "createShellMesh::calcPointRegions(..)"
-                            )   << "On point " << pointI
-                                << " at:" << patch.localPoints()[pointI]
-                                << " found region:" << region
-                                << abort(FatalError);
-                        }
-                        region = nRegions;
+                        isChangedFace[faceI] = true;
+                        changedFaces.append(faceI);
                     }
                 }
 
-                nRegions++;
+                label fp1 = findIndex(f, e[1]);
+                if (pointGlobalRegions[faceI][fp1] > edgeData[1])
+                {
+                    pointGlobalRegions[faceI][fp1] = edgeData[1];
+                    if (!isChangedFace[faceI])
+                    {
+                        isChangedFace[faceI] = true;
+                        changedFaces.append(faceI);
+                    }
+                }
             }
         }
+
+
+        label nChangedFaces = returnReduce(changedFaces.size(), sumOp<label>());
+        if (nChangedFaces == 0)
+        {
+            break;
+        }
+
+
+        // From face to edge
+        // ~~~~~~~~~~~~~~~~~
+
+        isChangedEdge = false;
+        changedEdges.clear();
+
+        forAll(changedFaces, i)
+        {
+            label faceI = changedFaces[i];
+            const face& f = patch.localFaces()[faceI];
+            const labelList& fEdges = patch.faceEdges()[faceI];
+
+            forAll(fEdges, fp)
+            {
+                label edgeI = fEdges[fp];
+
+                if (!nonManifoldEdge[edgeI])
+                {
+                    const edge& e = patch.edges()[edgeI];
+                    label fp0 = findIndex(f, e[0]);
+                    label region0 = pointGlobalRegions[faceI][fp0];
+                    label fp1 = findIndex(f, e[1]);
+                    label region1 = pointGlobalRegions[faceI][fp1];
+
+                    if
+                    (
+                        (allEdgeData[edgeI][0] > region0)
+                     || (allEdgeData[edgeI][1] > region1)
+                    )
+                    {
+                        allEdgeData[edgeI] = labelPair(region0, region1);
+                        if (!isChangedEdge[edgeI])
+                        {
+                            changedEdges.append(edgeI);
+                            isChangedEdge[edgeI] = true;
+                        }
+                    }
+                }
+            }
+        }
+
+        syncEdges
+        (
+            globalData,
+
+            patchEdges,
+            coupledEdges,
+            sameEdgeOrientation,
+
+            isChangedEdge,
+            changedEdges,
+            allEdgeData
+        );
+
+
+        label nChangedEdges = returnReduce(changedEdges.size(), sumOp<label>());
+        if (nChangedEdges == 0)
+        {
+            break;
+        }
     }
 
 
-    // From region back to originating point (many to one, a point might
-    // have multiple regions though)
-    regionPoints.setSize(nRegions);
-    forAll(pointRegions, faceI)
+
+    // Assign local regions
+    // ~~~~~~~~~~~~~~~~~~~~
+
+    // Calculate addressing from global region back to local region
+    pointLocalRegions.setSize(patch.size());
+    Map<label> globalToLocalRegion(globalRegions.localSize()/4);
+    DynamicList<label> dynLocalToGlobalRegion(globalToLocalRegion.size());
+    forAll(patch.localFaces(), faceI)
     {
         const face& f = patch.localFaces()[faceI];
-
+        face& pRegions = pointLocalRegions[faceI];
+        pRegions.setSize(f.size());
         forAll(f, fp)
         {
-            regionPoints[pointRegions[faceI][fp]] = f[fp];
-        }
-    }
+            label globalRegionI = pointGlobalRegions[faceI][fp];
 
+            Map<label>::iterator fnd = globalToLocalRegion.find(globalRegionI);
 
-    if (debug)
-    {
-        const labelListList& pointFaces = patch.pointFaces();
-        forAll(pointFaces, pointI)
-        {
-            label region = -1;
-            const labelList& pFaces = pointFaces[pointI];
-            forAll(pFaces, i)
+            if (fnd != globalToLocalRegion.end())
             {
-                label faceI = pFaces[i];
-                const face& f = patch.localFaces()[faceI];
-                label fp = findIndex(f, pointI);
-
-                if (region == -1)
-                {
-                    region = pointRegions[faceI][fp];
-                }
-                else if (region != pointRegions[faceI][fp])
-                {
-                    Pout<< "Non-manifold point:" << pointI
-                        << " at " << patch.localPoints()[pointI]
-                        << " region:" << region
-                        << " otherRegion:" << pointRegions[faceI][fp]
-                        << endl;
-
-                }
+                // Already encountered this global region. Assign same local one
+                pRegions[fp] = fnd();
+            }
+            else
+            {
+                // Region not yet seen. Create new one
+                label localRegionI = globalToLocalRegion.size();
+                pRegions[fp] = localRegionI;
+                globalToLocalRegion.insert(globalRegionI, localRegionI);
+                dynLocalToGlobalRegion.append(globalRegionI);
             }
         }
     }
+    localToGlobalRegion.transfer(dynLocalToGlobalRegion);
 }
 
 
diff --git a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh.H b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh.H
index 627191c6a6..7517a956df 100644
--- a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh.H
+++ b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh.H
@@ -41,6 +41,7 @@ SourceFiles
 
 #include "primitiveFacePatch.H"
 #include "PackedBoolList.H"
+#include "labelPair.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -50,6 +51,7 @@ namespace Foam
 // Forward declaration of classes
 class mapPolyMesh;
 class polyTopoChange;
+class globalMeshData;
 
 /*---------------------------------------------------------------------------*\
                            Class createShellMesh Declaration
@@ -80,6 +82,18 @@ class createShellMesh
 
     // Private Member Functions
 
+        static void syncEdges
+        (
+            const globalMeshData&,
+            const labelList&,
+            const labelList&,
+            const PackedBoolList& sameEdgeOrientation,
+            PackedBoolList& isChangedEdge,
+            DynamicList<label>& changedEdges,
+            labelPairList& allEdgeData
+        );
+
+
         //- Disallow default bitwise copy construct
         createShellMesh(const createShellMesh&);
 
@@ -142,13 +156,22 @@ public:
 
         // Other
 
-            //- Helper: calculate point regions
+            //- Helper: calculate point regions. The point region is the
+            //  same on all faces connected to a point if they can be
+            //  reached through a face-edge-face walk without crossing
+            //  the nonManifoldEdge.
+            //  pointGlobalRegions : non-compact. Guaranteed to be the same
+            //  across processors.
+            //  pointLocalRegions : compact.
+            //  localToGlobalRegion : local to global region.
             static void calcPointRegions
             (
+                const globalMeshData& globalData,
                 const primitiveFacePatch& patch,
                 const PackedBoolList& nonManifoldEdge,
-                faceList& pointRegions,
-                labelList& regionPoints
+                faceList& pointGlobalRegions,
+                faceList& pointLocalRegions,
+                labelList& localToGlobalRegion
             );
 
             //- Play commands into polyTopoChange to create layer mesh.
diff --git a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C
index 0a085485f1..40e4021284 100644
--- a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C
+++ b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C
@@ -120,58 +120,93 @@ Usage
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
-#include "Time.H"
 #include "fvMesh.H"
 #include "polyTopoChange.H"
-#include "patchPointEdgeCirculator.H"
 #include "OFstream.H"
 #include "meshTools.H"
 #include "mappedWallPolyPatch.H"
 #include "createShellMesh.H"
-#include "volFields.H"
-#include "surfaceFields.H"
 #include "syncTools.H"
 #include "cyclicPolyPatch.H"
 #include "wedgePolyPatch.H"
 #include "nonuniformTransformCyclicPolyPatch.H"
 #include "extrudeModel.H"
+#include "globalIndex.H"
+#include "addPatchCellLayer.H"
+
+#include "volFields.H"
+#include "surfaceFields.H"
+#include "pointFields.H"
+//#include "ReadFields.H"
+
 
 using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-template<class GeoField>
-void addPatchFields(fvMesh& mesh, const word& patchFieldType)
-{
-    HashTable<const GeoField*> flds
-    (
-        mesh.objectRegistry::lookupClass<GeoField>()
-    );
-
-    forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
-    {
-        const GeoField& fld = *iter();
-
-        typename GeoField::GeometricBoundaryField& bfld =
-            const_cast<typename GeoField::GeometricBoundaryField&>
-            (
-                fld.boundaryField()
-            );
-
-        label sz = bfld.size();
-        bfld.setSize(sz+1);
-        bfld.set
-        (
-            sz,
-            GeoField::PatchFieldType::New
-            (
-                patchFieldType,
-                mesh.boundary()[sz],
-                fld.dimensionedInternalField()
-            )
-        );
-    }
-}
+//template<class GeoField>
+//void addPatchFields(const fvMesh& mesh, const word& patchFieldType)
+//{
+//    HashTable<const GeoField*> flds
+//    (
+//        mesh.objectRegistry::lookupClass<GeoField>()
+//    );
+//
+//    forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
+//    {
+//        const GeoField& fld = *iter();
+//
+//        typename GeoField::GeometricBoundaryField& bfld =
+//            const_cast<typename GeoField::GeometricBoundaryField&>
+//            (
+//                fld.boundaryField()
+//            );
+//
+//        label sz = bfld.size();
+//
+//        for (label i = 0; i < sz; i++)
+//        {
+//            bfld.set
+//            (
+//                i,
+//                bfld.clone(GeoField::PatchFieldType::New
+//                (
+//                    patchFieldType,
+//                    fld.mesh().boundary()[sz],
+//                    fld.dimensionedInternalField()
+//                )
+//            );
+//
+//
+//
+//        Pout<< "fld:" << fld.name() << " had " << sz << " patches." << endl;
+//        Pout<< "fld before:" << fld << endl;
+//        Pout<< "adding on patch:" << fld.mesh().boundary()[sz].name() << endl;
+//
+//        bfld.setSize(sz+1);
+//        bfld.set
+//        (
+//            sz,
+//            GeoField::PatchFieldType::New
+//            (
+//                patchFieldType,
+//                fld.mesh().boundary()[sz],
+//                fld.dimensionedInternalField()
+//            )
+//        );
+//
+//        bfld[sz].operator=(pTraits<typename GeoField::value_type>::zero);
+//
+//        Pout<< "fld:" << fld.name() << " now " << bfld.size() << " patches."
+//            << endl;
+//
+//        const typename GeoField::PatchFieldType& pfld = bfld[sz];
+//        Pout<< "pfld:" << pfld << endl;
+//
+//
+//        Pout<< "fld value:" << fld << endl;
+//    }
+//}
 
 
 // Remove last patch field
@@ -219,138 +254,448 @@ void reorderPatchFields(fvMesh& mesh, const labelList& oldToNew)
 }
 
 
-void addAllPatchFields(fvMesh& mesh, const label insertPatchI)
+//void addCalculatedPatchFields(const fvMesh& mesh)
+//{
+//    addPatchFields<volScalarField>
+//    (
+//        mesh,
+//        calculatedFvPatchField<scalar>::typeName
+//    );
+//    addPatchFields<volVectorField>
+//    (
+//        mesh,
+//        calculatedFvPatchField<vector>::typeName
+//    );
+//    addPatchFields<volSphericalTensorField>
+//    (
+//        mesh,
+//        calculatedFvPatchField<sphericalTensor>::typeName
+//    );
+//    addPatchFields<volSymmTensorField>
+//    (
+//        mesh,
+//        calculatedFvPatchField<symmTensor>::typeName
+//    );
+//    addPatchFields<volTensorField>
+//    (
+//        mesh,
+//        calculatedFvPatchField<tensor>::typeName
+//    );
+//
+//    // Surface fields
+//
+//    addPatchFields<surfaceScalarField>
+//    (
+//        mesh,
+//        calculatedFvPatchField<scalar>::typeName
+//    );
+//    addPatchFields<surfaceVectorField>
+//    (
+//        mesh,
+//        calculatedFvPatchField<vector>::typeName
+//    );
+//    addPatchFields<surfaceSphericalTensorField>
+//    (
+//        mesh,
+//        calculatedFvPatchField<sphericalTensor>::typeName
+//    );
+//    addPatchFields<surfaceSymmTensorField>
+//    (
+//        mesh,
+//        calculatedFvPatchField<symmTensor>::typeName
+//    );
+//    addPatchFields<surfaceTensorField>
+//    (
+//        mesh,
+//        calculatedFvPatchField<tensor>::typeName
+//    );
+//
+//    // Point fields
+//
+//    addPatchFields<pointScalarField>
+//    (
+//        mesh,
+//        calculatedFvPatchField<scalar>::typeName
+//    );
+//    addPatchFields<pointVectorField>
+//    (
+//        mesh,
+//        calculatedFvPatchField<vector>::typeName
+//    );
+//}
+//
+//
+//void addAllPatchFields(fvMesh& mesh, const label insertPatchI)
+//{
+//    polyBoundaryMesh& polyPatches =
+//        const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
+//    fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
+//
+//    label sz = polyPatches.size();
+//
+//    addPatchFields<volScalarField>
+//    (
+//        mesh,
+//        calculatedFvPatchField<scalar>::typeName
+//    );
+//    addPatchFields<volVectorField>
+//    (
+//        mesh,
+//        calculatedFvPatchField<vector>::typeName
+//    );
+//    addPatchFields<volSphericalTensorField>
+//    (
+//        mesh,
+//        calculatedFvPatchField<sphericalTensor>::typeName
+//    );
+//    addPatchFields<volSymmTensorField>
+//    (
+//        mesh,
+//        calculatedFvPatchField<symmTensor>::typeName
+//    );
+//    addPatchFields<volTensorField>
+//    (
+//        mesh,
+//        calculatedFvPatchField<tensor>::typeName
+//    );
+//
+//    // Surface fields
+//
+//    addPatchFields<surfaceScalarField>
+//    (
+//        mesh,
+//        calculatedFvPatchField<scalar>::typeName
+//    );
+//    addPatchFields<surfaceVectorField>
+//    (
+//        mesh,
+//        calculatedFvPatchField<vector>::typeName
+//    );
+//    addPatchFields<surfaceSphericalTensorField>
+//    (
+//        mesh,
+//        calculatedFvPatchField<sphericalTensor>::typeName
+//    );
+//    addPatchFields<surfaceSymmTensorField>
+//    (
+//        mesh,
+//        calculatedFvPatchField<symmTensor>::typeName
+//    );
+//    addPatchFields<surfaceTensorField>
+//    (
+//        mesh,
+//        calculatedFvPatchField<tensor>::typeName
+//    );
+//
+//    // Create reordering list
+//    // patches before insert position stay as is
+//    labelList oldToNew(sz);
+//    for (label i = 0; i < insertPatchI; i++)
+//    {
+//        oldToNew[i] = i;
+//    }
+//    // patches after insert position move one up
+//    for (label i = insertPatchI; i < sz-1; i++)
+//    {
+//        oldToNew[i] = i+1;
+//    }
+//    // appended patch gets moved to insert position
+//    oldToNew[sz-1] = insertPatchI;
+//
+//    // Shuffle into place
+//    polyPatches.reorder(oldToNew);
+//    fvPatches.reorder(oldToNew);
+//
+//    reorderPatchFields<volScalarField>(mesh, oldToNew);
+//    reorderPatchFields<volVectorField>(mesh, oldToNew);
+//    reorderPatchFields<volSphericalTensorField>(mesh, oldToNew);
+//    reorderPatchFields<volSymmTensorField>(mesh, oldToNew);
+//    reorderPatchFields<volTensorField>(mesh, oldToNew);
+//    reorderPatchFields<surfaceScalarField>(mesh, oldToNew);
+//    reorderPatchFields<surfaceVectorField>(mesh, oldToNew);
+//    reorderPatchFields<surfaceSphericalTensorField>(mesh, oldToNew);
+//    reorderPatchFields<surfaceSymmTensorField>(mesh, oldToNew);
+//    reorderPatchFields<surfaceTensorField>(mesh, oldToNew);
+//}
+
+
+//// Adds patch if not yet there. Returns patchID.
+//template<class PatchType>
+//label addPatch(fvMesh& mesh, const word& patchName, const dictionary& dict)
+//{
+//    polyBoundaryMesh& polyPatches =
+//        const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
+//
+//    label patchI = polyPatches.findPatchID(patchName);
+//    if (patchI != -1)
+//    {
+//        if (isA<PatchType>(polyPatches[patchI]))
+//        {
+//            // Already there
+//            return patchI;
+//        }
+//        else
+//        {
+//            FatalErrorIn("addPatch<PatchType>(fvMesh&, const word&)")
+//                << "Already have patch " << patchName
+//                << " but of type " << PatchType::typeName
+//                << exit(FatalError);
+//        }
+//    }
+//
+//
+//    label insertPatchI = polyPatches.size();
+//    label startFaceI = mesh.nFaces();
+//
+//    forAll(polyPatches, patchI)
+//    {
+//        const polyPatch& pp = polyPatches[patchI];
+//
+//        if (isA<processorPolyPatch>(pp))
+//        {
+//            insertPatchI = patchI;
+//            startFaceI = pp.start();
+//            break;
+//        }
+//    }
+//
+//    dictionary patchDict(dict);
+//    patchDict.set("type", PatchType::typeName);
+//    patchDict.set("nFaces", 0);
+//    patchDict.set("startFace", startFaceI);
+//
+//
+//    // Below is all quite a hack. Feel free to change once there is a better
+//    // mechanism to insert and reorder patches.
+//
+//    // Clear local fields and e.g. polyMesh parallelInfo.
+//    mesh.clearOut();
+//
+//    label sz = polyPatches.size();
+//
+//    fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
+//
+//    // Add polyPatch at the end
+//    polyPatches.setSize(sz+1);
+//    polyPatches.set
+//    (
+//        sz,
+//        polyPatch::New
+//        (
+//            patchName,
+//            patchDict,
+//            insertPatchI,
+//            polyPatches
+//        )
+//    );
+//    fvPatches.setSize(sz+1);
+//    fvPatches.set
+//    (
+//        sz,
+//        fvPatch::New
+//        (
+//            polyPatches[sz],  // point to newly added polyPatch
+//            mesh.boundary()
+//        )
+//    );
+//
+//    addAllPatchFields(mesh, insertPatchI);
+//
+//    return insertPatchI;
+//}
+//
+//
+//template<class PatchType>
+//label addPatch(fvMesh& mesh, const word& patchName)
+//{
+//Pout<< "addPatch:" << patchName << endl;
+//
+//    polyBoundaryMesh& polyPatches =
+//        const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
+//
+//    label patchI = polyPatches.findPatchID(patchName);
+//    if (patchI != -1)
+//    {
+//        if (isA<PatchType>(polyPatches[patchI]))
+//        {
+//            // Already there
+//            return patchI;
+//        }
+//        else
+//        {
+//            FatalErrorIn("addPatch<PatchType>(fvMesh&, const word&)")
+//                << "Already have patch " << patchName
+//                << " but of type " << PatchType::typeName
+//                << exit(FatalError);
+//        }
+//    }
+//
+//
+//    label insertPatchI = polyPatches.size();
+//    label startFaceI = mesh.nFaces();
+//
+//    forAll(polyPatches, patchI)
+//    {
+//        const polyPatch& pp = polyPatches[patchI];
+//
+//        if (isA<processorPolyPatch>(pp))
+//        {
+//            insertPatchI = patchI;
+//            startFaceI = pp.start();
+//            break;
+//        }
+//    }
+//
+//    // Below is all quite a hack. Feel free to change once there is a better
+//    // mechanism to insert and reorder patches.
+//
+//    // Clear local fields and e.g. polyMesh parallelInfo.
+//    mesh.clearOut();
+//
+//    label sz = polyPatches.size();
+//
+//    fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
+//
+//    // Add polyPatch at the end
+//    polyPatches.setSize(sz+1);
+//    polyPatches.set
+//    (
+//        sz,
+//        polyPatch::New
+//        (
+//            PatchType::typeName,
+//            patchName,
+//            0,              // size
+//            startFaceI,
+//            insertPatchI,
+//            polyPatches
+//        )
+//    );
+//    fvPatches.setSize(sz+1);
+//    fvPatches.set
+//    (
+//        sz,
+//        fvPatch::New
+//        (
+//            polyPatches[sz],  // point to newly added polyPatch
+//            mesh.boundary()
+//        )
+//    );
+//
+//    addAllPatchFields(mesh, insertPatchI);
+//
+//    return insertPatchI;
+//}
+
+
+label findPatchID(const List<polyPatch*>& newPatches, const word& name)
 {
-    polyBoundaryMesh& polyPatches =
-        const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
-    fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
-
-    label sz = polyPatches.size();
-
-    addPatchFields<volScalarField>
-    (
-        mesh,
-        calculatedFvPatchField<scalar>::typeName
-    );
-    addPatchFields<volVectorField>
-    (
-        mesh,
-        calculatedFvPatchField<vector>::typeName
-    );
-    addPatchFields<volSphericalTensorField>
-    (
-        mesh,
-        calculatedFvPatchField<sphericalTensor>::typeName
-    );
-    addPatchFields<volSymmTensorField>
-    (
-        mesh,
-        calculatedFvPatchField<symmTensor>::typeName
-    );
-    addPatchFields<volTensorField>
-    (
-        mesh,
-        calculatedFvPatchField<tensor>::typeName
-    );
-
-    // Surface fields
-
-    addPatchFields<surfaceScalarField>
-    (
-        mesh,
-        calculatedFvPatchField<scalar>::typeName
-    );
-    addPatchFields<surfaceVectorField>
-    (
-        mesh,
-        calculatedFvPatchField<vector>::typeName
-    );
-    addPatchFields<surfaceSphericalTensorField>
-    (
-        mesh,
-        calculatedFvPatchField<sphericalTensor>::typeName
-    );
-    addPatchFields<surfaceSymmTensorField>
-    (
-        mesh,
-        calculatedFvPatchField<symmTensor>::typeName
-    );
-    addPatchFields<surfaceTensorField>
-    (
-        mesh,
-        calculatedFvPatchField<tensor>::typeName
-    );
-
-    // Create reordering list
-    // patches before insert position stay as is
-    labelList oldToNew(sz);
-    for (label i = 0; i < insertPatchI; i++)
+    forAll(newPatches, i)
     {
-        oldToNew[i] = i;
+        if (newPatches[i]->name() == name)
+        {
+            return i;
+        }
     }
-    // patches after insert position move one up
-    for (label i = insertPatchI; i < sz-1; i++)
-    {
-        oldToNew[i] = i+1;
-    }
-    // appended patch gets moved to insert position
-    oldToNew[sz-1] = insertPatchI;
-
-    // Shuffle into place
-    polyPatches.reorder(oldToNew);
-    fvPatches.reorder(oldToNew);
-
-    reorderPatchFields<volScalarField>(mesh, oldToNew);
-    reorderPatchFields<volVectorField>(mesh, oldToNew);
-    reorderPatchFields<volSphericalTensorField>(mesh, oldToNew);
-    reorderPatchFields<volSymmTensorField>(mesh, oldToNew);
-    reorderPatchFields<volTensorField>(mesh, oldToNew);
-    reorderPatchFields<surfaceScalarField>(mesh, oldToNew);
-    reorderPatchFields<surfaceVectorField>(mesh, oldToNew);
-    reorderPatchFields<surfaceSphericalTensorField>(mesh, oldToNew);
-    reorderPatchFields<surfaceSymmTensorField>(mesh, oldToNew);
-    reorderPatchFields<surfaceTensorField>(mesh, oldToNew);
+    return -1;
 }
 
 
-// Adds patch if not yet there. Returns patchID.
 template<class PatchType>
-label addPatch(fvMesh& mesh, const word& patchName, const dictionary& dict)
+label addPatch
+(
+    const polyBoundaryMesh& patches,
+    const word& patchName,
+    DynamicList<polyPatch*>& newPatches
+)
 {
-    polyBoundaryMesh& polyPatches =
-        const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
+    label patchI = findPatchID(newPatches, patchName);
 
-    label patchI = polyPatches.findPatchID(patchName);
     if (patchI != -1)
     {
-        if (isA<PatchType>(polyPatches[patchI]))
+        if (isA<PatchType>(*newPatches[patchI]))
         {
             // Already there
             return patchI;
         }
         else
         {
-            FatalErrorIn("addPatch<PatchType>(fvMesh&, const word&)")
-                << "Already have patch " << patchName
-                << " but of type " << PatchType::typeName
+            FatalErrorIn
+            (
+                "addPatch<PatchType>(const polyBoundaryMesh&,"
+                " const word&, DynamicList<polyPatch*>)"
+            )   << "Already have patch " << patchName
+                << " but of type " << newPatches[patchI]->type()
                 << exit(FatalError);
         }
     }
 
 
-    label insertPatchI = polyPatches.size();
-    label startFaceI = mesh.nFaces();
+    patchI = newPatches.size();
 
-    forAll(polyPatches, patchI)
+    label startFaceI = 0;
+    if (patchI > 0)
     {
-        const polyPatch& pp = polyPatches[patchI];
+        const polyPatch& pp = *newPatches.last();
+        startFaceI = pp.start()+pp.size();
+    }
 
-        if (isA<processorPolyPatch>(pp))
+
+    newPatches.append
+    (
+        polyPatch::New
+        (
+            PatchType::typeName,
+            patchName,
+            0,                          // size
+            startFaceI,                 // nFaces
+            patchI,
+            patches
+        ).ptr()
+    );
+
+    return patchI;
+}
+
+
+template<class PatchType>
+label addPatch
+(
+    const polyBoundaryMesh& patches,
+    const word& patchName,
+    const dictionary& dict,
+    DynamicList<polyPatch*>& newPatches
+)
+{
+    label patchI = findPatchID(newPatches, patchName);
+
+    if (patchI != -1)
+    {
+        if (isA<PatchType>(*newPatches[patchI]))
         {
-            insertPatchI = patchI;
-            startFaceI = pp.start();
-            break;
+            // Already there
+            return patchI;
         }
+        else
+        {
+            FatalErrorIn
+            (
+                "addPatch<PatchType>(const polyBoundaryMesh&,"
+                " const word&, DynamicList<polyPatch*>)"
+            )   << "Already have patch " << patchName
+                << " but of type " << newPatches[patchI]->type()
+                << exit(FatalError);
+        }
+    }
+
+
+    patchI = newPatches.size();
+
+    label startFaceI = 0;
+    if (patchI > 0)
+    {
+        const polyPatch& pp = *newPatches.last();
+        startFaceI = pp.start()+pp.size();
     }
 
     dictionary patchDict(dict);
@@ -358,125 +703,18 @@ label addPatch(fvMesh& mesh, const word& patchName, const dictionary& dict)
     patchDict.set("nFaces", 0);
     patchDict.set("startFace", startFaceI);
 
-
-    // Below is all quite a hack. Feel free to change once there is a better
-    // mechanism to insert and reorder patches.
-
-    // Clear local fields and e.g. polyMesh parallelInfo.
-    mesh.clearOut();
-
-    label sz = polyPatches.size();
-
-    fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
-
-    // Add polyPatch at the end
-    polyPatches.setSize(sz+1);
-    polyPatches.set
+    newPatches.append
     (
-        sz,
         polyPatch::New
         (
             patchName,
             patchDict,
-            insertPatchI,
-            polyPatches
-        )
-    );
-    fvPatches.setSize(sz+1);
-    fvPatches.set
-    (
-        sz,
-        fvPatch::New
-        (
-            polyPatches[sz],  // point to newly added polyPatch
-            mesh.boundary()
-        )
+            patchI,
+            patches
+        ).ptr()
     );
 
-    addAllPatchFields(mesh, insertPatchI);
-
-    return insertPatchI;
-}
-
-
-template<class PatchType>
-label addPatch(fvMesh& mesh, const word& patchName)
-{
-    polyBoundaryMesh& polyPatches =
-        const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
-
-    label patchI = polyPatches.findPatchID(patchName);
-    if (patchI != -1)
-    {
-        if (isA<PatchType>(polyPatches[patchI]))
-        {
-            // Already there
-            return patchI;
-        }
-        else
-        {
-            FatalErrorIn("addPatch<PatchType>(fvMesh&, const word&)")
-                << "Already have patch " << patchName
-                << " but of type " << PatchType::typeName
-                << exit(FatalError);
-        }
-    }
-
-
-    label insertPatchI = polyPatches.size();
-    label startFaceI = mesh.nFaces();
-
-    forAll(polyPatches, patchI)
-    {
-        const polyPatch& pp = polyPatches[patchI];
-
-        if (isA<processorPolyPatch>(pp))
-        {
-            insertPatchI = patchI;
-            startFaceI = pp.start();
-            break;
-        }
-    }
-
-    // Below is all quite a hack. Feel free to change once there is a better
-    // mechanism to insert and reorder patches.
-
-    // Clear local fields and e.g. polyMesh parallelInfo.
-    mesh.clearOut();
-
-    label sz = polyPatches.size();
-
-    fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
-
-    // Add polyPatch at the end
-    polyPatches.setSize(sz+1);
-    polyPatches.set
-    (
-        sz,
-        polyPatch::New
-        (
-            PatchType::typeName,
-            patchName,
-            0,              // size
-            startFaceI,
-            insertPatchI,
-            polyPatches
-        )
-    );
-    fvPatches.setSize(sz+1);
-    fvPatches.set
-    (
-        sz,
-        fvPatch::New
-        (
-            polyPatches[sz],  // point to newly added polyPatch
-            mesh.boundary()
-        )
-    );
-
-    addAllPatchFields(mesh, insertPatchI);
-
-    return insertPatchI;
+    return patchI;
 }
 
 
@@ -506,6 +744,8 @@ void reorderPatches
     reorderPatchFields<surfaceSphericalTensorField>(mesh, oldToNew);
     reorderPatchFields<surfaceSymmTensorField>(mesh, oldToNew);
     reorderPatchFields<surfaceTensorField>(mesh, oldToNew);
+    reorderPatchFields<pointScalarField>(mesh, oldToNew);
+    reorderPatchFields<pointVectorField>(mesh, oldToNew);
 
     // Remove last.
     polyPatches.setSize(nNewPatches);
@@ -520,6 +760,8 @@ void reorderPatches
     trimPatchFields<surfaceSphericalTensorField>(mesh, nNewPatches);
     trimPatchFields<surfaceSymmTensorField>(mesh, nNewPatches);
     trimPatchFields<surfaceTensorField>(mesh, nNewPatches);
+    trimPatchFields<pointScalarField>(mesh, nNewPatches);
+    trimPatchFields<pointVectorField>(mesh, nNewPatches);
 }
 
 
@@ -531,15 +773,53 @@ void deleteEmptyPatches(fvMesh& mesh)
     labelList oldToNew(patches.size());
     label usedI = 0;
     label notUsedI = patches.size();
+
+    //Pout<< "deleteEmptyPatches:" << endl;
+    //forAll(patches, patchI)
+    //{
+    //    Pout<< "    patch:" << patchI << " name:" << patches[patchI].name()
+    //        << " start:" << patches[patchI].start()
+    //        << " nFaces:" << patches[patchI].size()
+    //        << " index:" << patches[patchI].index()
+    //        << endl;
+    //}
+    //Pout<< endl;
+
+
+    // Add all the non-empty, non-processor patches
     forAll(patches, patchI)
     {
-        if (returnReduce(patches[patchI].size(), sumOp<label>()) == 0)
+        if (isA<processorPolyPatch>(patches[patchI]))
         {
-            oldToNew[patchI] = --notUsedI;
+            // Processor patches are unique per processor so look at local
+            // size only
+            if (patches[patchI].size() == 0)
+            {
+                Pout<< "Deleting processor patch " << patchI
+                    << " name:" << patches[patchI].name()
+                    << endl;
+                oldToNew[patchI] = --notUsedI;
+            }
+            else
+            {
+                oldToNew[patchI] = usedI++;
+            }
         }
         else
         {
-            oldToNew[patchI] = usedI++;
+            // All non-processor patches are present everywhere to reduce
+            // size
+            if (returnReduce(patches[patchI].size(), sumOp<label>()) == 0)
+            {
+                Pout<< "Deleting patch " << patchI
+                    << " name:" << patches[patchI].name()
+                    << endl;
+                oldToNew[patchI] = --notUsedI;
+            }
+            else
+            {
+                oldToNew[patchI] = usedI++;
+            }
         }
     }
 
@@ -628,7 +908,11 @@ label findUncoveredPatchFace
 }
 
 
-// Count the number of faces in patches that need to be created
+// Count the number of faces in patches that need to be created. Calculates:
+//  zoneSidePatch[zoneI]         : the number of side faces to be created
+//  zoneZonePatch[zoneA,zoneB]   : the number of faces inbetween zoneA and B
+// Since this only counts we're not taking the processor patches into
+// account.
 void countExtrudePatches
 (
     const fvMesh& mesh,
@@ -691,200 +975,410 @@ void countExtrudePatches
 }
 
 
-// Constrain&sync normals on points that are on coupled patches to make sure
-// the face extruded from the edge has a valid normal with its coupled
-// equivalent.
-// Note that only points on cyclic edges need to be constrained and not
-// all points touching cyclics since only edges become faces.
-void constrainCoupledNormals
+void addCouplingPatches
+(
+    const fvMesh& mesh,
+    const word& regionName,
+    const word& shellRegionName,
+    const wordList& zoneNames,
+    const wordList& zoneShadowNames,
+    const boolList& isInternal,
+    DynamicList<polyPatch*>& newPatches,
+    labelList& interRegionTopPatch,
+    labelList& interRegionBottomPatch
+)
+{
+    Pout<< "Adding coupling patches:" << nl << nl
+        << "patchID\tpatch\ttype" << nl
+        << "-------\t-----\t----"
+        << endl;
+
+    interRegionTopPatch.setSize(zoneNames.size());
+    interRegionBottomPatch.setSize(zoneNames.size());
+
+    label nCoupled = 0;
+    forAll(zoneNames, i)
+    {
+        word interName(regionName+"_to_"+shellRegionName+'_'+zoneNames[i]);
+
+        if (isInternal[i])
+        {
+            interRegionTopPatch[i] = addPatch<mappedWallPolyPatch>
+            (
+                mesh.boundaryMesh(),
+                interName + "_top",
+                newPatches
+            );
+            nCoupled++;
+            Pout<< interRegionTopPatch[i]
+                << '\t' << newPatches[interRegionTopPatch[i]]->name()
+                << '\t' << newPatches[interRegionTopPatch[i]]->type()
+                << nl;
+
+            interRegionBottomPatch[i] = addPatch<mappedWallPolyPatch>
+            (
+                mesh.boundaryMesh(),
+                interName + "_bottom",
+                newPatches
+            );
+            nCoupled++;
+            Pout<< interRegionBottomPatch[i]
+                << '\t' << newPatches[interRegionBottomPatch[i]]->name()
+                << '\t' << newPatches[interRegionBottomPatch[i]]->type()
+                << nl;
+        }
+        else if (zoneShadowNames.size() == 0)
+        {
+            interRegionTopPatch[i] = addPatch<polyPatch>
+            (
+                mesh.boundaryMesh(),
+                zoneNames[i] + "_top",
+                newPatches
+            );
+            nCoupled++;
+            Pout<< interRegionTopPatch[i]
+                << '\t' << newPatches[interRegionTopPatch[i]]->name()
+                << '\t' << newPatches[interRegionTopPatch[i]]->type()
+                << nl;
+
+            interRegionBottomPatch[i] = addPatch<mappedWallPolyPatch>
+            (
+                mesh.boundaryMesh(),
+                interName,
+                newPatches
+            );
+            nCoupled++;
+            Pout<< interRegionBottomPatch[i]
+                << '\t' << newPatches[interRegionBottomPatch[i]]->name()
+                << '\t' << newPatches[interRegionBottomPatch[i]]->type()
+                << nl;
+        }
+        else    //patch using shadow face zones.
+        {
+            interRegionTopPatch[i] = addPatch<mappedWallPolyPatch>
+            (
+                mesh.boundaryMesh(),
+                zoneShadowNames[i] + "_top",
+                newPatches
+            );
+            nCoupled++;
+            Pout<< interRegionTopPatch[i]
+                << '\t' << newPatches[interRegionTopPatch[i]]->name()
+                << '\t' << newPatches[interRegionTopPatch[i]]->type()
+                << nl;
+
+            interRegionBottomPatch[i] = addPatch<mappedWallPolyPatch>
+            (
+                mesh.boundaryMesh(),
+                interName,
+                newPatches
+            );
+            nCoupled++;
+            Pout<< interRegionBottomPatch[i]
+                << '\t' << newPatches[interRegionBottomPatch[i]]->name()
+                << '\t' << newPatches[interRegionBottomPatch[i]]->type()
+                << nl;
+        }
+    }
+    Pout<< "Added " << nCoupled << " inter-region patches." << nl
+        << endl;
+}
+
+
+void addProcPatches
 (
     const fvMesh& mesh,
     const primitiveFacePatch& extrudePatch,
-    const labelList& meshEdges,
-    const faceList& pointRegions,   // per face, per index the region
+    const labelList& extrudeMeshFaces,
 
-    vectorField& regionNormals
+    labelList& sidePatchID,
+    DynamicList<polyPatch*>& newPatches
 )
 {
+    // Global face indices engine
+    const globalIndex globalFaces(mesh.nFaces());
+
+    indirectPrimitivePatch pp
+    (
+        IndirectList<face>(mesh.faces(), extrudeMeshFaces),
+        mesh.points()
+    );
+
+    forAll(extrudePatch.edges(), edgeI)
+    {
+        const edge& extrudeEdge = extrudePatch.edges()[edgeI];
+        const edge& ppEdge = pp.edges()[edgeI];
+
+        if (extrudeEdge != ppEdge)
+        {
+            FatalErrorIn("addProcPatches()")
+                << "Problem: extrudeEdge:" << extrudeEdge
+                << " ppEdge:" << ppEdge
+                << exit(FatalError);
+        }
+    }
+
+
+    // Determine extrudePatch.edgeFaces in global numbering (so across
+    // coupled patches).
+    labelListList edgeGlobalFaces
+    (
+        addPatchCellLayer::globalEdgeFaces
+        (
+            mesh,
+            globalFaces,
+            pp
+        )
+    );
+
+    // Calculate for every edge the patch to use. This will be an existing
+    // patch for uncoupled edge and a possibly new patch (in patchToNbrProc)
+    // for processor patches.
+    label nNewPatches;
+    Map<label> nbrProcToPatch;
+    Map<label> patchToNbrProc;
+
+    addPatchCellLayer::calcSidePatch
+    (
+        mesh,
+        globalFaces,
+        edgeGlobalFaces,
+        pp,
+
+        sidePatchID,
+        nNewPatches,
+        nbrProcToPatch,
+        patchToNbrProc
+    );
+
+
+    // All patchIDs calcSidePatch are in mesh.boundaryMesh() numbering.
+    // Redo in newPatches numbering.
+
+    Pout<< "Adding inter-processor patches:" << nl << nl
+        << "patchID\tpatch" << nl
+        << "-------\t-----"
+        << endl;
+
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
-    // Mark edges that are on boundary of extrusion.
-    Map<label> meshToExtrudEdge
-    (
-        2*(extrudePatch.nEdges()-extrudePatch.nInternalEdges())
-    );
-    for
-    (
-        label extrudeEdgeI = extrudePatch.nInternalEdges();
-        extrudeEdgeI < extrudePatch.nEdges();
-        extrudeEdgeI++
-    )
+    forAll(sidePatchID, edgeI)
     {
-        meshToExtrudEdge.insert(meshEdges[extrudeEdgeI], extrudeEdgeI);
-    }
-
-
-    // For owner: normal at first point of edge when walking through faces
-    // in order.
-    vectorField edgeNormals0(mesh.nEdges(), vector::zero);
-    vectorField edgeNormals1(mesh.nEdges(), vector::zero);
-
-    // Loop through all edges of patch. If they are to be extruded mark the
-    // point normals in order.
-    forAll(patches, patchI)
-    {
-        const polyPatch& pp = patches[patchI];
-
-        if (isA<cyclicPolyPatch>(pp))
+        label meshPatchI = sidePatchID[edgeI];
+        if (meshPatchI != -1)
         {
-            bool isOwner = refCast<const cyclicPolyPatch>(pp).owner();
-
-            forAll(pp.faceEdges(), faceI)
+            label newPatchI = -1;
+            if (meshPatchI < patches.size())
             {
-                const labelList& fEdges = pp.faceEdges()[faceI];
-                forAll(fEdges, fp)
-                {
-                    label meshEdgeI = pp.meshEdges()[fEdges[fp]];
-                    if (meshToExtrudEdge.found(meshEdgeI))
-                    {
-                        // Edge corresponds to a extrusion edge. Store extrusion
-                        // normals on edge so we can syncTools it.
-
-                        //const edge& ppE = pp.edges()[fEdges[fp]];
-                        //Pout<< "ppedge:" << pp.localPoints()[ppE[0]]
-                        //    << pp.localPoints()[ppE[1]]
-                        //    << endl;
-
-                        const face& f = pp.localFaces()[faceI];
-                        label fp0 = fp;
-                        label fp1 = f.fcIndex(fp0);
-                        label mp0 = pp[faceI][fp0];
-                        label mp1 = pp[faceI][fp1];
-
-                        // Find corresponding face and indices.
-                        vector regionN0;
-                        vector regionN1;
-                        {
-                            label exEdgeI = meshToExtrudEdge[meshEdgeI];
-                            const labelList& eFaces =
-                                extrudePatch.edgeFaces()[exEdgeI];
-                            // Use 0th face.
-                            label exFaceI = eFaces[0];
-                            const face& exF = extrudePatch[exFaceI];
-                            const face& exRegions = pointRegions[exFaceI];
-                            // Find points
-                            label r0 = exRegions[findIndex(exF, mp0)];
-                            regionN0 = regionNormals[r0];
-                            label r1 = exRegions[findIndex(exF, mp1)];
-                            regionN1 = regionNormals[r1];
-                        }
-
-                        vector& nA =
-                        (
-                            isOwner
-                          ? edgeNormals0[meshEdgeI]
-                          : edgeNormals1[meshEdgeI]
-                        );
-
-                        nA = regionN0;
-                        const vector& cyc0 = pp.pointNormals()[f[fp0]];
-                        nA -= (nA&cyc0)*cyc0;
-
-                        vector& nB =
-                        (
-                            isOwner
-                          ? edgeNormals1[meshEdgeI]
-                          : edgeNormals0[meshEdgeI]
-                        );
-
-                        nB = regionN1;
-                        const vector& cyc1 = pp.pointNormals()[f[fp1]];
-                        nB -= (nB&cyc1)*cyc1;
-                    }
-                }
+                // Existing mesh patch. See if I already have it.
+                newPatchI = findPatchID
+                (
+                    newPatches,
+                    patches[meshPatchI].name()
+                );
             }
+
+            if (newPatchI == -1)
+            {
+                // Newly added processor patch
+                label nbrProcI = patchToNbrProc[meshPatchI];
+                word name =
+                        "procBoundary"
+                      + Foam::name(Pstream::myProcNo())
+                      + "to"
+                      + Foam::name(nbrProcI);
+
+                dictionary patchDict;
+                patchDict.add("myProcNo", Pstream::myProcNo());
+                patchDict.add("neighbProcNo", nbrProcI);
+
+                newPatchI = addPatch<processorPolyPatch>
+                (
+                    mesh.boundaryMesh(),
+                    name,
+                    patchDict,
+                    newPatches
+                );
+
+                Pout<< newPatchI << '\t' << name
+                    << nl;
+            }
+
+            sidePatchID[edgeI] = newPatchI;
         }
     }
 
 
-    // Synchronise regionNormals
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Synchronise
-    syncTools::syncEdgeList
-    (
-        mesh,
-        edgeNormals0,
-        maxMagSqrEqOp<vector>(),
-        vector::zero            // nullValue
-    );
-    syncTools::syncEdgeList
-    (
-        mesh,
-        edgeNormals1,
-        maxMagSqrEqOp<vector>(),
-        vector::zero            // nullValue
-    );
-
-
-    // Re-work back into regionNormals
-    forAll(patches, patchI)
+    // Clear out sidePatchID for uncoupled edges. Just so we don't have
+    // to expose all the globalEdgeFaces info.
+    forAll(sidePatchID, edgeI)
     {
-        const polyPatch& pp = patches[patchI];
-
-        if (isA<cyclicPolyPatch>(pp))
+        if
+        (
+            edgeGlobalFaces[edgeI].size() == 2
+         && pp.edgeFaces()[edgeI].size() == 1
+        )
+        {}
+        else
         {
-            bool isOwner = refCast<const cyclicPolyPatch>(pp).owner();
+            sidePatchID[edgeI] = -1;
+        }
+    }
+}
 
-            forAll(pp.faceEdges(), faceI)
+
+void addZoneSidePatches
+(
+    const fvMesh& mesh,
+    const word& oneDPolyPatchType,
+    const wordList& zoneNames,
+
+    DynamicList<polyPatch*>& newPatches,
+    labelList& zoneSidePatch
+)
+{
+    Pout<< "Adding patches for sides on zones:" << nl << nl
+        << "patchID\tpatch" << nl
+        << "-------\t-----"
+        << endl;
+
+    label nSide = 0;
+
+    forAll(zoneNames, zoneI)
+    {
+        if (oneDPolyPatchType != word::null)
+        {
+            // Reuse single empty patch.
+            word patchName;
+            if (oneDPolyPatchType == "emptyPolyPatch")
             {
-                const labelList& fEdges = pp.faceEdges()[faceI];
-                forAll(fEdges, fp)
+                patchName = "oneDEmptyPatch";
+                zoneSidePatch[zoneI] = addPatch<emptyPolyPatch>
+                (
+                    mesh.boundaryMesh(),
+                    patchName,
+                    newPatches
+                );
+            }
+            else if (oneDPolyPatchType == "wedgePolyPatch")
+            {
+                patchName = "oneDWedgePatch";
+                zoneSidePatch[zoneI] = addPatch<wedgePolyPatch>
+                (
+                    mesh.boundaryMesh(),
+                    patchName,
+                    newPatches
+                );
+            }
+            else
+            {
+                FatalErrorIn("addZoneSidePatches()")
+                    << "Type " << oneDPolyPatchType << " does not exist "
+                    << exit(FatalError);
+            }
+
+            Pout<< zoneSidePatch[zoneI] << '\t' << patchName << nl;
+
+            nSide++;
+        }
+        else if (zoneSidePatch[zoneI] > 0)
+        {
+            word patchName = zoneNames[zoneI] + "_" + "side";
+
+            zoneSidePatch[zoneI] = addPatch<polyPatch>
+            (
+                mesh.boundaryMesh(),
+                patchName,
+                newPatches
+            );
+
+            Pout<< zoneSidePatch[zoneI] << '\t' << patchName << nl;
+
+            nSide++;
+        }
+    }
+    Pout<< "Added " << nSide << " zone-side patches." << nl
+        << endl;
+}
+
+
+void addInterZonePatches
+(
+    const fvMesh& mesh,
+    const wordList& zoneNames,
+    const bool oneD,
+
+    labelList& zoneZonePatch_min,
+    labelList& zoneZonePatch_max,
+    DynamicList<polyPatch*>& newPatches
+)
+{
+    Pout<< "Adding inter-zone patches:" << nl << nl
+        << "patchID\tpatch" << nl
+        << "-------\t-----"
+        << endl;
+
+    dictionary transformDict;
+    transformDict.add
+    (
+        "transform",
+        cyclicPolyPatch::transformTypeNames[cyclicPolyPatch::NOORDERING]
+    );
+
+    label nInter = 0;
+    if (!oneD)
+    {
+        forAll(zoneZonePatch_min, minZone)
+        {
+            for (label maxZone = minZone; maxZone < zoneNames.size(); maxZone++)
+            {
+                label index = minZone*zoneNames.size()+maxZone;
+
+                if (zoneZonePatch_min[index] > 0)
                 {
-                    label meshEdgeI = pp.meshEdges()[fEdges[fp]];
-                    if (meshToExtrudEdge.found(meshEdgeI))
+                    word minToMax =
+                        zoneNames[minZone]
+                      + "_to_"
+                      + zoneNames[maxZone];
+                    word maxToMin =
+                        zoneNames[maxZone]
+                      + "_to_"
+                      + zoneNames[minZone];
+
                     {
-                        const face& f = pp.localFaces()[faceI];
-                        label fp0 = fp;
-                        label fp1 = f.fcIndex(fp0);
-                        label mp0 = pp[faceI][fp0];
-                        label mp1 = pp[faceI][fp1];
-
-
-                        const vector& nA =
+                        transformDict.set("neighbourPatch", maxToMin);
+                        zoneZonePatch_min[index] =
+                        addPatch<nonuniformTransformCyclicPolyPatch>
                         (
-                            isOwner
-                          ? edgeNormals0[meshEdgeI]
-                          : edgeNormals1[meshEdgeI]
+                            mesh.boundaryMesh(),
+                            minToMax,
+                            transformDict,
+                            newPatches
                         );
-
-                        const vector& nB =
-                        (
-                            isOwner
-                          ? edgeNormals1[meshEdgeI]
-                          : edgeNormals0[meshEdgeI]
-                        );
-
-                        // Find corresponding face and indices.
-                        {
-                            label exEdgeI = meshToExtrudEdge[meshEdgeI];
-                            const labelList& eFaces =
-                                extrudePatch.edgeFaces()[exEdgeI];
-                            // Use 0th face.
-                            label exFaceI = eFaces[0];
-                            const face& exF = extrudePatch[exFaceI];
-                            const face& exRegions = pointRegions[exFaceI];
-                            // Find points
-                            label r0 = exRegions[findIndex(exF, mp0)];
-                            regionNormals[r0] = nA;
-                            label r1 = exRegions[findIndex(exF, mp1)];
-                            regionNormals[r1] = nB;
-                        }
+                        Pout<< zoneZonePatch_min[index] << '\t' << minToMax
+                            << nl;
+                        nInter++;
                     }
+                    {
+                        transformDict.set("neighbourPatch", minToMax);
+                        zoneZonePatch_max[index] =
+                        addPatch<nonuniformTransformCyclicPolyPatch>
+                        (
+                            mesh.boundaryMesh(),
+                            maxToMin,
+                            transformDict,
+                            newPatches
+                        );
+                        Pout<< zoneZonePatch_max[index] << '\t' << maxToMin
+                            << nl;
+                        nInter++;
+                    }
+
                 }
             }
         }
     }
+    Pout<< "Added " << nInter << " inter-zone patches." << nl
+        << endl;
 }
 
 
@@ -914,23 +1408,83 @@ tmp<pointField> calcOffset
 }
 
 
+void setCouplingInfo
+(
+    fvMesh& mesh,
+    const labelList& zoneToPatch,
+    const word& sampleRegion,
+    const mappedWallPolyPatch::sampleMode mode,
+    const List<pointField>& offsets
+)
+{
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    List<polyPatch*> newPatches(patches.size(), NULL);
+
+    forAll(zoneToPatch, zoneI)
+    {
+        label patchI = zoneToPatch[zoneI];
+
+        const polyPatch& pp = patches[patchI];
+
+        if (isA<mappedWallPolyPatch>(pp))
+        {
+            newPatches[patchI] = new mappedWallPolyPatch
+            (
+                pp.name(),
+                pp.size(),
+                pp.start(),
+                patchI,
+                sampleRegion,                           // sampleRegion
+                mode,                                   // sampleMode
+                pp.name(),                              // samplePatch
+                offsets[zoneI],                         // offset
+                patches
+            );
+        }
+    }
+
+    forAll(newPatches, patchI)
+    {
+        if (!newPatches[patchI])
+        {
+            newPatches[patchI] = patches[patchI].clone(patches).ptr();
+        }
+    }
+
+    mesh.removeFvBoundary();
+    mesh.addFvPatches(newPatches, true);
+}
+
 
 // Main program:
 
 int main(int argc, char *argv[])
 {
-    argList::noParallel();
     argList::addNote("Create region mesh by extruding a faceZone");
 
     #include "addRegionOption.H"
     #include "addOverwriteOption.H"
-    argList::addOption("dict", "name", "specify alternative dictionary");
-    argList::addBoolOption("AMI", "apply mapped AMI boundary type");
+    argList::addNote("Create region mesh by extruding a faceZone");
 
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createNamedMesh.H"
 
+    if (mesh.boundaryMesh().checkParallelSync(true))
+    {
+        List<wordList> allNames(Pstream::nProcs());
+        allNames[Pstream::myProcNo()] = mesh.boundaryMesh().names();
+        Pstream::gatherList(allNames);
+        Pstream::scatterList(allNames);
+
+        FatalErrorIn(args.executable())
+            << "Patches are not synchronised on all processors."
+            << " Per processor patches " << allNames
+            << exit(FatalError);
+    }
+
+
     const word oldInstance = mesh.pointsInstance();
     bool overwrite = args.optionFound("overwrite");
     const word dictName
@@ -966,7 +1520,7 @@ int main(int argc, char *argv[])
     const Switch oneD(dict.lookup("oneD"));
     const Switch adaptMesh(dict.lookup("adaptMesh"));
 
-    Info<< "Extruding zones " << zoneNames
+    Pout<< "Extruding zones " << zoneNames
         << " on mesh " << regionName
         << " into shell mesh " << shellRegionName
         << endl;
@@ -982,6 +1536,53 @@ int main(int argc, char *argv[])
 
 
 
+    //// Read objects in time directory
+    //IOobjectList objects(mesh, runTime.timeName());
+    //
+    //// Read vol fields.
+    //
+    //PtrList<volScalarField> vsFlds;
+    //ReadFields(mesh, objects, vsFlds);
+    //
+    //PtrList<volVectorField> vvFlds;
+    //ReadFields(mesh, objects, vvFlds);
+    //
+    //PtrList<volSphericalTensorField> vstFlds;
+    //ReadFields(mesh, objects, vstFlds);
+    //
+    //PtrList<volSymmTensorField> vsymtFlds;
+    //ReadFields(mesh, objects, vsymtFlds);
+    //
+    //PtrList<volTensorField> vtFlds;
+    //ReadFields(mesh, objects, vtFlds);
+    //
+    //// Read surface fields.
+    //
+    //PtrList<surfaceScalarField> ssFlds;
+    //ReadFields(mesh, objects, ssFlds);
+    //
+    //PtrList<surfaceVectorField> svFlds;
+    //ReadFields(mesh, objects, svFlds);
+    //
+    //PtrList<surfaceSphericalTensorField> sstFlds;
+    //ReadFields(mesh, objects, sstFlds);
+    //
+    //PtrList<surfaceSymmTensorField> ssymtFlds;
+    //ReadFields(mesh, objects, ssymtFlds);
+    //
+    //PtrList<surfaceTensorField> stFlds;
+    //ReadFields(mesh, objects, stFlds);
+    //
+    //// Read point fields.
+    //
+    //PtrList<pointScalarField> psFlds;
+    //ReadFields(pointMesh::New(mesh), objects, psFlds);
+    //
+    //PtrList<pointVectorField> pvFlds;
+    //ReadFields(pointMesh::New(mesh), objects, pvFlds);
+
+
+
     // Create dummy fv* files
     createDummyFvMeshFiles(mesh, shellRegionName);
 
@@ -996,7 +1597,7 @@ int main(int argc, char *argv[])
     {
         meshInstance = oldInstance;
     }
-    Info<< "Writing meshes to " << meshInstance << nl << endl;
+    Pout<< "Writing meshes to " << meshInstance << nl << endl;
 
 
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
@@ -1093,7 +1694,7 @@ int main(int argc, char *argv[])
     const labelListList& edgeFaces = extrudePatch.edgeFaces();
 
 
-    Info<< "extrudePatch :"
+    Pout<< "extrudePatch :"
         << " faces:" << extrudePatch.size()
         << " points:" << extrudePatch.nPoints()
         << " edges:" << extrudePatch.nEdges()
@@ -1149,112 +1750,92 @@ int main(int argc, char *argv[])
         const faceZone& fz = faceZones[zoneIDs[i]];
         if (isInternal[i])
         {
-            Info<< "FaceZone " << fz.name() << " has internal faces" << endl;
+            Pout<< "FaceZone " << fz.name() << " has internal faces" << endl;
         }
         else
         {
-            Info<< "FaceZone " << fz.name() << " has boundary faces" << endl;
+            Pout<< "FaceZone " << fz.name() << " has boundary faces" << endl;
         }
     }
-    Info<< endl;
+    Pout<< endl;
 
 
 
+    DynamicList<polyPatch*> regionPatches(patches.size());
+    // Copy all non-local patches since these are used on boundary edges of
+    // the extrusion
+    forAll(patches, patchI)
+    {
+        if (!isA<processorPolyPatch>(patches[patchI]))
+        {
+            label newPatchI = regionPatches.size();
+            regionPatches.append
+            (
+                patches[patchI].clone
+                (
+                    patches,
+                    newPatchI,
+                    0,              // size
+                    0               // start
+                ).ptr()
+            );
+        }
+    }
+
+
     // Add interface patches
     // ~~~~~~~~~~~~~~~~~~~~~
-    // Note that these actually get added to the original mesh
-    // so the shell mesh creation copies them. They then get removed
-    // from the original mesh.
 
-    Info<< "Adding coupling patches:" << nl << nl
-        << "patchID\tpatch\ttype" << nl
-        << "-------\t-----\t----"
-        << endl;
+    // From zone to interface patch (region side)
     labelList interRegionTopPatch(zoneNames.size());
     labelList interRegionBottomPatch(zoneNames.size());
-    label nCoupled = 0;
-    forAll(zoneIDs, i)
+
+    addCouplingPatches
+    (
+        mesh,
+        regionName,
+        shellRegionName,
+        zoneNames,
+        zoneShadowNames,
+        isInternal,
+
+        regionPatches,
+        interRegionTopPatch,
+        interRegionBottomPatch
+    );
+
+    // From zone to interface patch (mesh side)
+    labelList interMeshTopPatch(zoneNames.size());
+    labelList interMeshBottomPatch(zoneNames.size());
+
+    if (adaptMesh)
     {
-        word interName(regionName+"_to_"+shellRegionName+'_'+zoneNames[i]);
-
-        if (isInternal[i])
+        DynamicList<polyPatch*> newPatches(patches.size());
+        forAll(patches, patchI)
         {
-            interRegionTopPatch[i] = addPatch<mappedWallPolyPatch>
-            (
-                mesh,
-                interName + "_top"
-            );
-            nCoupled++;
-            Info<< interRegionTopPatch[i]
-                << '\t' << patches[interRegionTopPatch[i]].name()
-                << '\t' << patches[interRegionTopPatch[i]].type()
-                << nl;
-
-            interRegionBottomPatch[i] = addPatch<mappedWallPolyPatch>
-            (
-                mesh,
-                interName + "_bottom"
-            );
-            nCoupled++;
-            Info<< interRegionBottomPatch[i]
-                << '\t' << patches[interRegionBottomPatch[i]].name()
-                << '\t' << patches[interRegionBottomPatch[i]].type()
-                << nl;
-        }
-        else if (zoneShadowNames.size() == 0)
-        {
-            interRegionTopPatch[i] = addPatch<polyPatch>
-            (
-                mesh,
-                zoneNames[i] + "_top"
-            );
-            nCoupled++;
-            Info<< interRegionTopPatch[i]
-                << '\t' << patches[interRegionTopPatch[i]].name()
-                << '\t' << patches[interRegionTopPatch[i]].type()
-                << nl;
-
-            interRegionBottomPatch[i] = addPatch<mappedWallPolyPatch>
-            (
-                mesh,
-                interName
-            );
-            nCoupled++;
-            Info<< interRegionBottomPatch[i]
-                << '\t' << patches[interRegionBottomPatch[i]].name()
-                << '\t' << patches[interRegionBottomPatch[i]].type()
-                << nl;
-        }
-        else if (zoneShadowNames.size() > 0) //patch using shadow face zones.
-        {
-            interRegionTopPatch[i] = addPatch<mappedWallPolyPatch>
-            (
-                mesh,
-                zoneShadowNames[i] + "_top"
-            );
-            nCoupled++;
-            Info<< interRegionTopPatch[i]
-                << '\t' << patches[interRegionTopPatch[i]].name()
-                << '\t' << patches[interRegionTopPatch[i]].type()
-                << nl;
-
-            interRegionBottomPatch[i] = addPatch<mappedWallPolyPatch>
-            (
-                mesh,
-                interName
-            );
-            nCoupled++;
-            Info<< interRegionBottomPatch[i]
-                << '\t' << patches[interRegionBottomPatch[i]].name()
-                << '\t' << patches[interRegionBottomPatch[i]].type()
-                << nl;
+            newPatches.append(patches[patchI].clone(patches).ptr());
         }
 
+        addCouplingPatches
+        (
+            mesh,
+            regionName,
+            shellRegionName,
+            zoneNames,
+            zoneShadowNames,
+            isInternal,
+
+            newPatches,
+            interMeshTopPatch,
+            interMeshBottomPatch
+        );
+        mesh.clearOut();
+        mesh.removeFvBoundary();
+        mesh.addFvPatches(newPatches, true);
     }
-    Info<< "Added " << nCoupled << " inter-region patches." << nl
-        << endl;
 
 
+    // Patch per extruded face
     labelList extrudeTopPatchID(extrudePatch.size());
     labelList extrudeBottomPatchID(extrudePatch.size());
 
@@ -1294,136 +1875,79 @@ int main(int argc, char *argv[])
         zoneZonePatch_min   // reuse for counting
     );
 
-    // Now check which patches to add.
-    Info<< "Adding patches for edges on zones:" << nl << nl
-        << "patchID\tpatch" << nl
-        << "-------\t-----"
-        << endl;
-
-    label nSide = 0;
-
-    forAll(zoneSidePatch, zoneI)
-    {
-        if (oneD)
-        {
-            // Reuse single empty patch.
-            word patchType = dict.lookup("oneDPolyPatchType");
-            word patchName;
-            if (patchType == "emptyPolyPatch")
-            {
-                patchName = "oneDEmptyPatch";
-                zoneSidePatch[zoneI] = addPatch<emptyPolyPatch>
-                (
-                    mesh,
-                    patchName
-                );
-            }
-            else if (patchType == "wedgePolyPatch")
-            {
-                patchName = "oneDWedgePatch";
-                zoneSidePatch[zoneI] = addPatch<wedgePolyPatch>
-                (
-                    mesh,
-                    patchName
-                );
-            }
-            else
-            {
-                FatalErrorIn(args.executable())
-                << "Type " << patchType << " does not exist "
-                << exit(FatalError);
-            }
-
-            Info<< zoneSidePatch[zoneI] << '\t' << patchName << nl;
-
-            nSide++;
-        }
-        else if (zoneSidePatch[zoneI] > 0)
-        {
-            word patchName = faceZones[zoneI].name() + "_" + "side";
-
-            zoneSidePatch[zoneI] = addPatch<polyPatch>
-            (
-                mesh,
-                patchName
-            );
-
-            Info<< zoneSidePatch[zoneI] << '\t' << patchName << nl;
-
-            nSide++;
-        }
-    }
-    Info<< "Added " << nSide << " zone-edge patches." << nl
-        << endl;
+    // Now we'll have:
+    //  zoneSidePatch[zoneA] : number of faces needed on the side of zoneA
+    //  zoneZonePatch_min[zoneA,zoneB] : number of faces needed inbetween A,B
 
 
-
-    Info<< "Adding inter-zone patches:" << nl << nl
-        << "patchID\tpatch" << nl
-        << "-------\t-----"
-        << endl;
-
-    dictionary transformDict;
-    transformDict.add
+    // Add the zone-side patches.
+    addZoneSidePatches
     (
-        "transform",
-        cyclicPolyPatch::transformTypeNames[cyclicPolyPatch::NOORDERING]
+        mesh,
+        (oneD ? dict.lookup("oneDPolyPatchType") : word::null),
+        zoneNames,
+
+        regionPatches,
+        zoneSidePatch
     );
 
-    label nInter = 0;
-    if (!oneD)
-    {
-        forAll(zoneZonePatch_min, minZone)
-        {
-            for (label maxZone = minZone; maxZone < faceZones.size(); maxZone++)
-            {
-                label index = minZone*faceZones.size()+maxZone;
 
-                if (zoneZonePatch_min[index] > 0)
-                {
-                    word minToMax =
-                        faceZones[minZone].name()
-                      + "_to_"
-                      + faceZones[maxZone].name();
-                    word maxToMin =
-                        faceZones[maxZone].name()
-                      + "_to_"
-                      + faceZones[minZone].name();
+    // Add the patches inbetween zones
+    addInterZonePatches
+    (
+        mesh,
+        zoneNames,
+        oneD,
 
-                    {
-                        transformDict.set("neighbourPatch", maxToMin);
-                        zoneZonePatch_min[index] =
-                        addPatch<nonuniformTransformCyclicPolyPatch>
-                        (
-                            mesh,
-                            minToMax,
-                            transformDict
-                        );
-                        Info<< zoneZonePatch_min[index] << '\t' << minToMax
-                            << nl;
-                        nInter++;
-                    }
-                    {
-                        transformDict.set("neighbourPatch", minToMax);
-                        zoneZonePatch_max[index] =
-                        addPatch<nonuniformTransformCyclicPolyPatch>
-                        (
-                            mesh,
-                            maxToMin,
-                            transformDict
-                        );
-                        Info<< zoneZonePatch_max[index] << '\t' << maxToMin
-                            << nl;
-                        nInter++;
-                    }
+        zoneZonePatch_min,
+        zoneZonePatch_max,
+        regionPatches
+    );
 
-                }
-            }
-        }
-    }
-    Info<< "Added " << nInter << " inter-zone patches." << nl
-        << endl;
 
+    // Additionally check if any interprocessor patches need to be added.
+    // Reuses addPatchCellLayer functionality.
+    // Note: does not handle edges with > 2 faces?
+    labelList sidePatchID;
+    addProcPatches
+    (
+        mesh,
+        extrudePatch,
+        extrudeMeshFaces,
+
+        sidePatchID,
+        regionPatches
+    );
+
+
+//    // Add all the newPatches to the mesh and fields
+//    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+//    {
+//        forAll(newPatches, patchI)
+//        {
+//            Pout<< "Adding patch " << patchI
+//                << " name:" << newPatches[patchI]->name()
+//                << endl;
+//        }
+//        //label nOldPatches = mesh.boundary().size();
+//        mesh.clearOut();
+//        mesh.removeFvBoundary();
+//        mesh.addFvPatches(newPatches, true);
+//        //// Add calculated fvPatchFields for the added patches
+//        //for
+//        //(
+//        //    label patchI = nOldPatches;
+//        //    patchI < mesh.boundary().size();
+//        //    patchI++
+//        //)
+//        //{
+//        //    Pout<< "ADDing calculated to patch " << patchI
+//        //        << endl;
+//        //    addCalculatedPatchFields(mesh);
+//        //}
+//        //Pout<< "** Added " << mesh.boundary().size()-nOldPatches
+//        //    << " patches." << endl;
+//    }
 
 
     // Set patches to use for edges to be extruded into boundary faces
@@ -1433,7 +1957,7 @@ int main(int argc, char *argv[])
     // If empty size create an internal face.
     labelListList extrudeEdgePatches(extrudePatch.nEdges());
 
-    // Is edge an non-manifold edge
+    // Is edge a non-manifold edge
     PackedBoolList nonManifoldEdge(extrudePatch.nEdges());
 
     // Note: logic has to be same as in countExtrudePatches.
@@ -1483,6 +2007,15 @@ int main(int argc, char *argv[])
                 nonManifoldEdge[edgeI] = 1;
             }
         }
+        else if (sidePatchID[edgeI] != -1)
+        {
+            // Coupled extrusion
+            ePatches.setSize(eFaces.size());
+            forAll(eFaces, i)
+            {
+                ePatches[i] = sidePatchID[edgeI];
+            }
+        }
         else
         {
             label faceI = findUncoveredPatchFace
@@ -1494,8 +2027,12 @@ int main(int argc, char *argv[])
 
             if (faceI != -1)
             {
-                label patchI = mesh.boundaryMesh().whichPatch(faceI);
-                ePatches.setSize(eFaces.size(), patchI);
+                label newPatchI = findPatchID
+                (
+                    regionPatches,
+                    patches[patches.whichPatch(faceI)].name()
+                );
+                ePatches.setSize(eFaces.size(), newPatchI);
             }
             else
             {
@@ -1515,76 +2052,102 @@ int main(int argc, char *argv[])
     // ~~~~~~~~~~~~~~~~~~~~
 
     // Per face, per point the region number.
-    faceList pointRegions(extrudePatch.size());
-    // Per region the originating point
-    labelList regionPoints;
+    faceList pointGlobalRegions;
+    faceList pointLocalRegions;
+    labelList localToGlobalRegion;
 
     createShellMesh::calcPointRegions
     (
+        mesh.globalData(),
         extrudePatch,
         nonManifoldEdge,
 
-        pointRegions,
-        regionPoints
+        pointGlobalRegions,
+        pointLocalRegions,
+        localToGlobalRegion
     );
 
-
-    // Calculate a normal per region
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    vectorField regionNormals(regionPoints.size(), vector::zero);
-    vectorField regionCentres(regionPoints.size(), vector::zero);
-    labelList nRegionFaces(regionPoints.size(), 0);
-
-    forAll(pointRegions, faceI)
+    // Per local region an originating point
+    labelList localRegionPoints(localToGlobalRegion.size());
+    forAll(pointLocalRegions, faceI)
     {
-        const face& f = extrudeFaces[faceI];
-
-        forAll(f, fp)
+        const face& f = extrudePatch.localFaces()[faceI];
+        const face& pRegions = pointLocalRegions[faceI];
+        forAll(pRegions, fp)
         {
-            label region = pointRegions[faceI][fp];
-            regionNormals[region] += extrudePatch.faceNormals()[faceI];
-            regionCentres[region] += extrudePatch.faceCentres()[faceI];
-            nRegionFaces[region]++;
+            localRegionPoints[pRegions[fp]] = f[fp];
         }
     }
-    regionNormals /= mag(regionNormals);
-    forAll(regionCentres, regionI)
+
+    // Calculate region normals by reducing local region normals
+    pointField localRegionNormals(localToGlobalRegion.size());
     {
-        regionCentres[regionI] /= nRegionFaces[regionI];
+        pointField localSum(localToGlobalRegion.size(), vector::zero);
+        labelList localNum(localToGlobalRegion.size(), 0);
+
+        forAll(pointLocalRegions, faceI)
+        {
+            const face& pRegions = pointLocalRegions[faceI];
+            forAll(pRegions, fp)
+            {
+                label localRegionI = pRegions[fp];
+                localSum[localRegionI] += extrudePatch.faceNormals()[faceI];
+                localNum[localRegionI]++;
+            }
+        }
+
+        Map<point> globalSum(2*localToGlobalRegion.size());
+        Map<label> globalNum(2*localToGlobalRegion.size());
+
+        forAll(localSum, localRegionI)
+        {
+            label globalRegionI = localToGlobalRegion[localRegionI];
+            globalSum.insert(globalRegionI, localSum[localRegionI]);
+            globalNum.insert(globalRegionI, localNum[localRegionI]);
+        }
+
+        // Reduce
+        Pstream::mapCombineGather(globalSum, plusEqOp<point>());
+        Pstream::mapCombineScatter(globalSum);
+        Pstream::mapCombineGather(globalNum, plusEqOp<label>());
+        Pstream::mapCombineScatter(globalNum);
+
+        forAll(localToGlobalRegion, localRegionI)
+        {
+            label globalRegionI = localToGlobalRegion[localRegionI];
+            localRegionNormals[localRegionI] =
+                globalSum[globalRegionI]
+              / globalNum[globalRegionI];
+        }
+        localRegionNormals /= mag(localRegionNormals);
     }
 
 
-    // Constrain&sync normals on points that are on coupled patches.
-    constrainCoupledNormals
-    (
-        mesh,
-        extrudePatch,
-        extrudeMeshEdges,
-        pointRegions,
-
-        regionNormals
-    );
 
     // For debugging: dump hedgehog plot of normals
     if (false)
     {
-        OFstream str(runTime.path()/"regionNormals.obj");
+        OFstream str(runTime.path()/"localRegionNormals.obj");
         label vertI = 0;
 
         scalar thickness = model().sumThickness(1);
 
-        forAll(pointRegions, faceI)
+        forAll(pointLocalRegions, faceI)
         {
             const face& f = extrudeFaces[faceI];
 
             forAll(f, fp)
             {
-                label region = pointRegions[faceI][fp];
+                label region = pointLocalRegions[faceI][fp];
                 const point& pt = extrudePoints[f[fp]];
 
                 meshTools::writeOBJ(str, pt);
                 vertI++;
-                meshTools::writeOBJ(str, pt+thickness*regionNormals[region]);
+                meshTools::writeOBJ
+                (
+                    str,
+                    pt+thickness*localRegionNormals[region]
+                );
                 vertI++;
                 str << "l " << vertI-1 << ' ' << vertI << nl;
             }
@@ -1593,11 +2156,18 @@ int main(int argc, char *argv[])
 
 
     // Use model to create displacements of first layer
-    vectorField firstDisp(regionNormals.size());
+    vectorField firstDisp(localRegionNormals.size());
     forAll(firstDisp, regionI)
     {
-        const point& regionPt = regionCentres[regionI];
-        const vector& n = regionNormals[regionI];
+        //const point& regionPt = regionCentres[regionI];
+        const point& regionPt = extrudePatch.points()
+        [
+            extrudePatch.meshPoints()
+            [
+                localRegionPoints[regionI]
+            ]
+        ];
+        const vector& n = localRegionNormals[regionI];
         firstDisp[regionI] = model()(regionPt, n, 1) - regionPt;
     }
 
@@ -1606,13 +2176,47 @@ int main(int argc, char *argv[])
     // Create a new mesh
     // ~~~~~~~~~~~~~~~~~
 
-    createShellMesh extruder(extrudePatch, pointRegions, regionPoints);
+    createShellMesh extruder
+    (
+        extrudePatch,
+        pointLocalRegions,
+        localRegionPoints
+    );
 
 
-    autoPtr<fvMesh> regionMeshPtr;
     autoPtr<mapPolyMesh> shellMap;
+    fvMesh regionMesh
+    (
+        IOobject
+        (
+            shellRegionName,
+            meshInstance,
+            runTime,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE,
+            false
+        ),
+        xferCopy(pointField()),
+        xferCopy(faceList()),
+        xferCopy(labelList()),
+        xferCopy(labelList()),
+        false
+    );
+
+    // Add the new patches
+    forAll(regionPatches, patchI)
     {
-        polyTopoChange meshMod(mesh.boundaryMesh().size());
+        regionPatches[patchI] = regionPatches[patchI]->clone
+        (
+            regionMesh.boundaryMesh()
+        ).ptr();
+    }
+    regionMesh.clearOut();
+    regionMesh.removeFvBoundary();
+    regionMesh.addFvPatches(regionPatches, true);
+
+    {
+        polyTopoChange meshMod(regionPatches.size());
 
         extruder.setRefinement
         (
@@ -1625,22 +2229,12 @@ int main(int argc, char *argv[])
             meshMod
         );
 
-        shellMap = meshMod.makeMesh
+        shellMap = meshMod.changeMesh
         (
-            regionMeshPtr,     // mesh to create
-            IOobject
-            (
-                shellRegionName,
-                meshInstance,
-                runTime,
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE,
-                false
-            ),
-            mesh            // mesh to get original patch info from
+            regionMesh,     // mesh to change
+            false           // inflate
         );
     }
-    fvMesh& regionMesh = regionMeshPtr();
 
     // Necessary?
     regionMesh.setInstance(meshInstance);
@@ -1650,82 +2244,92 @@ int main(int argc, char *argv[])
     extruder.updateMesh(shellMap);
 
 
-    // Change top and bottom boundary conditions on regionMesh
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Calculate offsets from shell mesh back to original mesh
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    // Save offsets from shell mesh back to original mesh
     List<pointField> topOffsets(zoneIDs.size());
     List<pointField> bottomOffsets(zoneIDs.size());
 
+    forAll(regionMesh.boundaryMesh(), patchI)
     {
-        const polyBoundaryMesh& regionPatches = regionMesh.boundaryMesh();
-        List<polyPatch*> newPatches(regionPatches.size());
-        forAll(regionPatches, patchI)
+        const polyPatch& pp = regionMesh.boundaryMesh()[patchI];
+
+        if
+        (
+            isA<mappedWallPolyPatch>(pp)
+         && (findIndex(interRegionTopPatch, patchI) != -1)
+        )
         {
-            const polyPatch& pp = regionPatches[patchI];
-
-            if
-            (
-                isA<mappedWallPolyPatch>(pp)
-             && (findIndex(interRegionTopPatch, patchI) != -1)
-            )
-            {
-                label index = findIndex(interRegionTopPatch, patchI);
-
-                topOffsets[index] = calcOffset(extrudePatch, extruder, pp);
-
-                newPatches[patchI] = new mappedWallPolyPatch
-                (
-                    pp.name(),
-                    pp.size(),
-                    pp.start(),
-                    patchI,
-                    regionName,                             // sampleRegion
-                    sampleMode,                             // sampleMode
-                    pp.name(),                              // samplePatch
-                    topOffsets[index],                      // offset
-                    patches
-                );
-            }
-            else if
-            (
-                isA<mappedWallPolyPatch>(pp)
-             && (findIndex(interRegionBottomPatch, patchI) != -1)
-            )
-            {
-                label index = findIndex(interRegionBottomPatch, patchI);
-
-                bottomOffsets[index] = calcOffset(extrudePatch, extruder, pp);
-
-                newPatches[patchI] = new mappedWallPolyPatch
-                (
-                    pp.name(),
-                    pp.size(),
-                    pp.start(),
-                    patchI,
-                    regionName,                             // sampleRegion
-                    sampleMode,                             // sampleMode
-                    pp.name(),                              // samplePatch
-                    bottomOffsets[index],                   // offset
-                    patches
-                );
-            }
-            else
-            {
-                newPatches[patchI] = pp.clone
-                (
-                    regionPatches,
-                    patchI,
-                    pp.size(),
-                    pp.start()
-                ).ptr();
-            }
+            label zoneI = findIndex(interRegionTopPatch, patchI);
+            topOffsets[zoneI] = calcOffset(extrudePatch, extruder, pp);
         }
-        regionMesh.removeFvBoundary();
-        regionMesh.addFvPatches(newPatches, true);
+        else if
+        (
+            isA<mappedWallPolyPatch>(pp)
+         && (findIndex(interRegionBottomPatch, patchI) != -1)
+        )
+        {
+            label zoneI = findIndex(interRegionBottomPatch, patchI);
+
+            bottomOffsets[zoneI] = calcOffset(extrudePatch, extruder, pp);
+        }
+    }
+
+
+    // Change top and bottom boundary conditions on regionMesh
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    {
+        // Correct top patches for offset
+        setCouplingInfo
+        (
+            regionMesh,
+            interRegionTopPatch,
+            regionName,                 // name of main mesh
+            sampleMode,                 // sampleMode
+            topOffsets
+        );
+
+        // Correct bottom patches for offset
+        setCouplingInfo
+        (
+            regionMesh,
+            interRegionBottomPatch,
+            regionName,
+            sampleMode,                 // sampleMode
+            bottomOffsets
+        );
+
+        // Remove any unused patches
         deleteEmptyPatches(regionMesh);
     }
 
+    // Change top and bottom boundary conditions on main mesh
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    if (adaptMesh)
+    {
+        // Correct top patches for offset
+        setCouplingInfo
+        (
+            mesh,
+            interMeshTopPatch,
+            shellRegionName,                        // name of shell mesh
+            sampleMode,                             // sampleMode
+            -topOffsets
+        );
+
+        // Correct bottom patches for offset
+        setCouplingInfo
+        (
+            mesh,
+            interMeshBottomPatch,
+            shellRegionName,
+            sampleMode,
+            -bottomOffsets
+        );
+    }
+
 
 
     // Write addressing from region mesh back to originating patch
@@ -1799,7 +2403,7 @@ int main(int argc, char *argv[])
         "point to patch point addressing";
 
 
-    Info<< "Writing mesh " << regionMesh.name()
+    Pout<< "Writing mesh " << regionMesh.name()
         << " to " << regionMesh.facesInstance() << nl
         << endl;
 
@@ -1847,7 +2451,7 @@ int main(int argc, char *argv[])
                     mesh.faceOwner()[meshFaceI],// owner
                     -1,                         // neighbour
                     false,                      // face flip
-                    extrudeBottomPatchID[zoneFaceI],// patch for face
+                    interMeshBottomPatch[zoneI],// patch for face
                     zoneI,                      // zone for face
                     flip                        // face flip in zone
                 );
@@ -1861,7 +2465,7 @@ int main(int argc, char *argv[])
                     mesh.faceNeighbour()[meshFaceI],// owner
                     -1,                             // neighbour
                     true,                           // face flip
-                    extrudeBottomPatchID[zoneFaceI],// patch for face
+                    interMeshBottomPatch[zoneI],    // patch for face
                     zoneI,                          // zone for face
                     !flip                           // face flip in zone
                 );
@@ -1886,7 +2490,7 @@ int main(int argc, char *argv[])
                         mesh.faceOwner()[meshFaceI],// owner
                         -1,                         // neighbour
                         false,                      // face flip
-                        extrudeTopPatchID[zoneFaceI],// patch for face
+                        interMeshTopPatch[zoneI],   // patch for face
                         zoneI,                      // zone for face
                         flip                        // face flip in zone
                     );
@@ -1900,7 +2504,7 @@ int main(int argc, char *argv[])
                         mesh.faceNeighbour()[meshFaceI],// owner
                         -1,                             // neighbour
                         true,                           // face flip
-                        extrudeTopPatchID[zoneFaceI],   // patch for face
+                        interMeshTopPatch[zoneI],       // patch for face
                         zoneI,                          // zone for face
                         !flip                           // face flip in zone
                     );
@@ -1914,6 +2518,7 @@ int main(int argc, char *argv[])
             forAll(extrudeMeshFaces, zoneFaceI)
             {
                 label meshFaceI = extrudeMeshFaces[zoneFaceI];
+                label zoneI = zoneID[zoneFaceI];
                 bool flip = zoneFlipMap[zoneFaceI];
                 const face& f = mesh.faces()[meshFaceI];
 
@@ -1930,7 +2535,7 @@ int main(int argc, char *argv[])
                             -1,                             // master edge
                             meshFaceI,                      // master face
                             true,                           // flip flux
-                            extrudeTopPatchID[zoneFaceI],   // patch for face
+                            interMeshTopPatch[zoneI],       // patch for face
                             -1,                             // zone for face
                             false                           //face flip in zone
                         );
@@ -1947,7 +2552,7 @@ int main(int argc, char *argv[])
                         -1,                             // master edge
                         meshFaceI,                      // master face
                         false,                          // flip flux
-                        extrudeTopPatchID[zoneFaceI],   // patch for face
+                        interMeshTopPatch[zoneI],       // patch for face
                         -1,                             // zone for face
                         false                           // zone flip
                     );
@@ -1973,84 +2578,9 @@ int main(int argc, char *argv[])
         }
 
         mesh.setInstance(meshInstance);
-    }
 
 
-
-    // Change master and slave boundary conditions on originating mesh
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    if (adaptMesh)
-    {
-        const polyBoundaryMesh& patches = mesh.boundaryMesh();
-        List<polyPatch*> newPatches(patches.size());
-        forAll(patches, patchI)
-        {
-            const polyPatch& pp = patches[patchI];
-
-            if
-            (
-                isA<mappedWallPolyPatch>(pp)
-             && (findIndex(interRegionTopPatch, patchI) != -1)
-            )
-            {
-                label index = findIndex(interRegionTopPatch, patchI);
-                newPatches[patchI] = new mappedWallPolyPatch
-                (
-                    pp.name(),
-                    pp.size(),
-                    pp.start(),
-                    patchI,
-                    shellRegionName,                        // sampleRegion
-                    sampleMode,                             // sampleMode
-                    pp.name(),                              // samplePatch
-                    -topOffsets[index],                     // offset
-                    patches
-                );
-            }
-            else if
-            (
-                isA<mappedWallPolyPatch>(pp)
-             && (findIndex(interRegionBottomPatch, patchI) != -1)
-            )
-            {
-                label index = findIndex(interRegionBottomPatch, patchI);
-
-                newPatches[patchI] = new mappedWallPolyPatch
-                (
-                    pp.name(),
-                    pp.size(),
-                    pp.start(),
-                    patchI,
-                    shellRegionName,                        // sampleRegion
-                    sampleMode,                             // sampleMode
-                    pp.name(),                              // samplePatch
-                    -bottomOffsets[index],                  // offset
-                    patches
-                );
-            }
-            else
-            {
-                newPatches[patchI] = pp.clone
-                (
-                    patches,
-                    patchI,
-                    pp.size(),
-                    pp.start()
-                ).ptr();
-            }
-        }
-        mesh.removeFvBoundary();
-        mesh.addFvPatches(newPatches, true);
-
-        // Remove any unused patches
-        deleteEmptyPatches(mesh);
-    }
-
-
-    if (adaptMesh)
-    {
-        Info<< "Writing mesh " << mesh.name()
+        Pout<< "Writing mesh " << mesh.name()
             << " to " << mesh.facesInstance() << nl
             << endl;
 
diff --git a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/patchPointEdgeCirculator.C b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/patchPointEdgeCirculator.C
deleted file mode 100644
index ea8f509a29..0000000000
--- a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/patchPointEdgeCirculator.C
+++ /dev/null
@@ -1,67 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "patchPointEdgeCirculator.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-const Foam::patchPointEdgeCirculator
-Foam::patchPointEdgeCirculator::endConstIter
-(
-    *reinterpret_cast<primitiveFacePatch*>(0),  // primitiveFacePatch
-    *reinterpret_cast<PackedBoolList*>(0),      // PackedBoolList
-    -1,                                         // edgeID
-    -1,                                         // index
-    -1                                          // pointID
-);
-
-
-Foam::Ostream& Foam::operator<<
-(
-    Ostream& os,
-    const InfoProxy<patchPointEdgeCirculator>& ip
-)
-{
-    const patchPointEdgeCirculator& c = ip.t_;
-
-    const pointField& pts = c.patch_.localPoints();
-    const edge& e = c.patch_.edges()[c.edgeID_];
-    label faceI = c.faceID();
-
-    os  << "around point:" << c.pointID_
-        << " coord:" << pts[c.pointID_]
-        << " at edge:" << c.edgeID_
-        << " coord:" << pts[e.otherVertex(c.pointID_)]
-        << " in direction of face:" << faceI;
-
-    if (faceI != -1)
-    {
-        os  << " fc:" << c.patch_.localFaces()[faceI].centre(pts);
-    }
-    return os;
-}
-
-
-// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/patchPointEdgeCirculator.H b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/patchPointEdgeCirculator.H
deleted file mode 100644
index 4bc46bfb24..0000000000
--- a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/patchPointEdgeCirculator.H
+++ /dev/null
@@ -1,205 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::patchPointEdgeCirculator
-
-Description
-    Walks from starting edge/face around point on patch.
-
-    -# Use in-place: \n
-        \code
-            patchPointEdgeCirculator circ(..);
-
-            // Walk
-            do
-            {
-                Info<< "edge:" << circ.edgeID() << endl;
-                ++circ;
-            }
-            while (circ != circ.end());
-        \endcode
-
-    -# Use like STL iterator: \n
-        \code
-            patchPointEdgeCirculator circ(..);
-            for
-            (
-                patchPointEdgeCirculator iter = circ.begin();
-                iter != circ.end();
-                ++iter
-            )
-            {
-                Info<< "edge:" << iter.edgeID() << endl;
-            }
-        \endcode
-
-
-SourceFiles
-    patchPointEdgeCirculator.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef patchPointEdgeCirculator_H
-#define patchPointEdgeCirculator_H
-
-#include "face.H"
-#include "ListOps.H"
-#include "primitiveFacePatch.H"
-#include "PackedBoolList.H"
-#include "InfoProxy.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Forward declaration of classes
-
-/*---------------------------------------------------------------------------*\
-                           Class patchPointEdgeCirculator Declaration
-\*---------------------------------------------------------------------------*/
-
-class patchPointEdgeCirculator
-{
-    // Static data members
-
-        //- end iterator
-        static const patchPointEdgeCirculator endConstIter;
-
-
-    // Private data
-
-        //- Patch
-        const primitiveFacePatch& patch_;
-
-        const PackedBoolList& nonManifoldEdge_;
-
-        //- Current edge
-        label edgeID_;
-
-        //- Current direction (face, expressed as index into edgeFaces()[edgeID]
-        label index_;
-
-        //- Point
-        label pointID_;
-
-        //- Starting edge so we know when to stop.
-        label startEdgeID_;
-
-
-    // Private Member Functions
-
-        //- Set to end() iterator
-        inline void setEnd();
-
-        //- Set edgeID_ to be the other edge on the face that uses the
-        //  point.
-        inline void otherEdge(const label cellI);
-
-public:
-
-    // Constructors
-
-        //- Construct from components
-        inline patchPointEdgeCirculator
-        (
-            const primitiveFacePatch&,
-            const PackedBoolList& nonManifoldEdge,
-            const label edgeID,
-            const label index,
-            const label pointID
-        );
-
-        //- Construct as copy
-        inline patchPointEdgeCirculator(const patchPointEdgeCirculator&);
-
-
-    // Member Functions
-
-        inline label edgeID() const;
-
-        inline label index() const;
-
-        inline label pointID() const;
-
-        //- Helper: get face according to the index.
-        //  Returns -1 if on end.
-        inline label faceID() const;
-
-        //- Walk back until non-manifold edge (if any) or minimum edge.
-        inline void setCanonical();
-
-
-    // Member Operators
-
-        inline void operator=(const patchPointEdgeCirculator& iter);
-
-        inline bool operator==(const patchPointEdgeCirculator& iter) const;
-
-        inline bool operator!=(const patchPointEdgeCirculator& iter) const;
-
-        //- Step to next face.
-        inline patchPointEdgeCirculator& operator++();
-
-        //- iterator set to the beginning face. For internal edges this is
-        //  the current face. For boundary edges this is the first boundary face
-        //  reached from walking back (i.e. in opposite direction to ++)
-        inline patchPointEdgeCirculator begin() const;
-        inline patchPointEdgeCirculator cbegin() const;
-
-        //- iterator set to beyond the end of the walk.
-        inline const patchPointEdgeCirculator& end() const;
-        inline const patchPointEdgeCirculator& cend() const;
-
-    // Info
-
-        //- Return info proxy.
-        //  Used to print token information to a stream
-        InfoProxy<patchPointEdgeCirculator> info() const
-        {
-            return *this;
-        }
-
-        friend Ostream& operator<<
-        (
-            Ostream&,
-            const InfoProxy<patchPointEdgeCirculator>&
-        );
-};
-
-Ostream& operator<<(Ostream&, const InfoProxy<patchPointEdgeCirculator>&);
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "patchPointEdgeCirculatorI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/patchPointEdgeCirculatorI.H b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/patchPointEdgeCirculatorI.H
deleted file mode 100644
index 55844065a7..0000000000
--- a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/patchPointEdgeCirculatorI.H
+++ /dev/null
@@ -1,308 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-void Foam::patchPointEdgeCirculator::setEnd()
-{
-    edgeID_ = -1;
-    pointID_ = -1;
-}
-
-
-// Cross face to other edge on point
-void Foam::patchPointEdgeCirculator::otherEdge(const label faceI)
-{
-    const labelList& fEdges = patch_.faceEdges()[faceI];
-    const face& f = patch_.localFaces()[faceI];
-    label fp = findIndex(f, pointID_);
-
-    if (fEdges[fp] == edgeID_)
-    {
-        edgeID_ = fEdges[f.rcIndex(fp)];
-    }
-    else
-    {
-        // Check for now
-        if (fEdges[f.rcIndex(fp)] != edgeID_)
-        {
-            FatalErrorIn("patchPointEdgeCirculator::otherEdge(const label)")
-                << "face:" << faceI
-                << " verts:" << f
-                << " edges:" << fEdges
-                << " looking for other edge than " << edgeID_
-                << abort(FatalError);
-        }
-        edgeID_ = fEdges[fp];
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-//- Construct from components
-Foam::patchPointEdgeCirculator::patchPointEdgeCirculator
-(
-    const primitiveFacePatch& patch,
-    const PackedBoolList& nonManifoldEdge,
-    const label edgeID,
-    const label index,
-
-    const label pointID
-)
-:
-    patch_(patch),
-    nonManifoldEdge_(nonManifoldEdge),
-    edgeID_(edgeID),
-    index_(index),
-    pointID_(pointID),
-    startEdgeID_(edgeID_)
-{
-    if (edgeID_ != -1)
-    {
-        const edge& e = patch_.edges()[edgeID_];
-
-        if (pointID_ != e[0] && pointID_ != e[1])
-        {
-            FatalErrorIn
-            (
-                "patchPointEdgeCirculator::patchPointEdgeCirculator(..)"
-            )   << "edge " << edgeID_ << " verts " << e
-                << " does not contain point " << pointID_ << abort(FatalError);
-        }
-    }
-}
-
-
-//- Construct copy
-Foam::patchPointEdgeCirculator::patchPointEdgeCirculator
-(
-    const patchPointEdgeCirculator& circ
-)
-:
-    patch_(circ.patch_),
-    nonManifoldEdge_(circ.nonManifoldEdge_),
-    edgeID_(circ.edgeID_),
-    index_(circ.index_),
-    pointID_(circ.pointID_),
-    startEdgeID_(circ.startEdgeID_)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-Foam::label Foam::patchPointEdgeCirculator::edgeID() const
-{
-    return edgeID_;
-}
-
-
-Foam::label Foam::patchPointEdgeCirculator::index() const
-{
-    return index_;
-}
-
-
-Foam::label Foam::patchPointEdgeCirculator::pointID() const
-{
-    return pointID_;
-}
-
-
-Foam::label Foam::patchPointEdgeCirculator::faceID() const
-{
-    if (edgeID_ != -1 && index_ != -1)
-    {
-        return patch_.edgeFaces()[edgeID_][index_];
-    }
-    else
-    {
-        return -1;
-    }
-}
-
-
-void Foam::patchPointEdgeCirculator::operator=
-(
-    const patchPointEdgeCirculator& circ
-)
-{
-    edgeID_ = circ.edgeID_;
-    index_ = circ.index_;
-    pointID_ = circ.pointID_;
-    startEdgeID_ = circ.startEdgeID_;
-}
-
-
-bool Foam::patchPointEdgeCirculator::operator==
-(
-    const patchPointEdgeCirculator& circ
-) const
-{
-    // Do just enough to have setEnd() produce an iterator equal to end().
-    // Could include the direction(i.e. index_) to make sure that two
-    // circulators around same point but in different direction are not equal.
-    return edgeID_ == circ.edgeID_ && pointID_ == circ.pointID_;
-}
-
-
-bool Foam::patchPointEdgeCirculator::operator!=
-(
-    const patchPointEdgeCirculator& circ
-) const
-{
-    return !(*this == circ);
-}
-
-
-void Foam::patchPointEdgeCirculator::setCanonical()
-{
-    if (edgeID_ == -1)
-    {
-        FatalErrorIn("patchPointEdgeCirculator::setCanonical()")
-            << "Already reached end(). Cannot walk any further."
-            << abort(FatalError);
-    }
-
-    label minEdgeID = edgeID_;
-    label minIndex = index_;
-
-    while (true)
-    {
-        if (nonManifoldEdge_[edgeID_])
-        {
-            break;
-        }
-
-        // Step back
-        const labelList& eFaces = patch_.edgeFaces()[edgeID_];
-
-        if (eFaces.size() != 2)
-        {
-            FatalErrorIn("patchPointEdgeCirculator::setCanonical()")
-                << "problem eFaces:" << eFaces << abort(FatalError);
-        }
-
-        label faceI = (index_ == 0 ? eFaces[1] : eFaces[0]);
-
-        // Step to other edge on face
-        otherEdge(faceI);
-
-        // Update index
-        index_ = findIndex(patch_.edgeFaces()[edgeID_], faceI);
-
-        if (edgeID_ < minEdgeID)
-        {
-            minEdgeID = edgeID_;
-            minIndex = index_;
-        }
-
-        if (edgeID_ == startEdgeID_)
-        {
-            edgeID_ = minEdgeID;
-            index_ = minIndex;
-            break;
-        }
-    }
-
-    startEdgeID_ = edgeID_;
-}
-
-
-//- Step to next edge.
-Foam::patchPointEdgeCirculator&
-Foam::patchPointEdgeCirculator::operator++()
-{
-    if (index_ == -1)
-    {
-        setEnd();
-    }
-    else
-    {
-        // Step to other edge on face
-        label faceI = patch_.edgeFaces()[edgeID_][index_];
-        otherEdge(faceI);
-
-        if (edgeID_ == startEdgeID_)
-        {
-            setEnd();
-        }
-        else if (nonManifoldEdge_[edgeID_])
-        {
-            // Reached non-manifold edge. Cannot walk further.
-            // Mark so it gets set to end next time.
-            index_ = -1;
-        }
-        else
-        {
-            const labelList& eFaces = patch_.edgeFaces()[edgeID_];
-
-            if (eFaces.size() != 2)
-            {
-                FatalErrorIn("patchPointEdgeCirculator:::operator++()")
-                    << "problem eFaces:" << eFaces << abort(FatalError);
-            }
-            // Point to the face that is not faceI
-            index_ = (eFaces[0] != faceI ? 0 : 1);
-        }
-    }
-
-    return *this;
-}
-
-
-Foam::patchPointEdgeCirculator Foam::patchPointEdgeCirculator::begin() const
-{
-    patchPointEdgeCirculator iter(*this);
-    iter.setCanonical();
-
-    return iter;
-}
-
-
-Foam::patchPointEdgeCirculator Foam::patchPointEdgeCirculator::cbegin() const
-{
-    patchPointEdgeCirculator iter(*this);
-    iter.setCanonical();
-
-    return iter;
-}
-
-
-const Foam::patchPointEdgeCirculator& Foam::patchPointEdgeCirculator::end()
- const
-{
-    return endConstIter;
-}
-
-const Foam::patchPointEdgeCirculator& Foam::patchPointEdgeCirculator::cend()
- const
-{
-    return endConstIter;
-}
-
-
-// ************************************************************************* //