From 06d2dc9e4cb6960c5ed4c8ec144f715baf661c2c Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Wed, 17 Mar 2010 11:17:17 +0000
Subject: [PATCH] ENH: application to extrude a faceZone into a separate region
 mesh

Extrudes either internal faces to inflate a baffle or extrudes boundary
faces to create separate mesh. Inserts directMappedWall between the two
meshes.
---
 .../generation/extrudeToRegionMesh/Make/files |    7 +
 .../extrudeToRegionMesh/Make/options          |   14 +
 .../extrudeToRegionMesh/createShellMesh.C     |  592 ++++++++
 .../extrudeToRegionMesh/createShellMesh.H     |  177 +++
 .../extrudeToRegionMesh/extrudeToRegionMesh.C | 1334 +++++++++++++++++
 .../patchPointEdgeCirculator.C                |   68 +
 .../patchPointEdgeCirculator.H                |  206 +++
 .../patchPointEdgeCirculatorI.H               |  309 ++++
 8 files changed, 2707 insertions(+)
 create mode 100644 applications/utilities/mesh/generation/extrudeToRegionMesh/Make/files
 create mode 100644 applications/utilities/mesh/generation/extrudeToRegionMesh/Make/options
 create mode 100644 applications/utilities/mesh/generation/extrudeToRegionMesh/createShellMesh.C
 create mode 100644 applications/utilities/mesh/generation/extrudeToRegionMesh/createShellMesh.H
 create mode 100644 applications/utilities/mesh/generation/extrudeToRegionMesh/extrudeToRegionMesh.C
 create mode 100644 applications/utilities/mesh/generation/extrudeToRegionMesh/patchPointEdgeCirculator.C
 create mode 100644 applications/utilities/mesh/generation/extrudeToRegionMesh/patchPointEdgeCirculator.H
 create mode 100644 applications/utilities/mesh/generation/extrudeToRegionMesh/patchPointEdgeCirculatorI.H

diff --git a/applications/utilities/mesh/generation/extrudeToRegionMesh/Make/files b/applications/utilities/mesh/generation/extrudeToRegionMesh/Make/files
new file mode 100644
index 0000000000..62431e23d4
--- /dev/null
+++ b/applications/utilities/mesh/generation/extrudeToRegionMesh/Make/files
@@ -0,0 +1,7 @@
+patchPointEdgeCirculator.C
+createShellMesh.C
+extrudeToRegionMesh.C
+
+EXE = $(FOAM_APPBIN)/extrudeToRegionMesh
+
+
diff --git a/applications/utilities/mesh/generation/extrudeToRegionMesh/Make/options b/applications/utilities/mesh/generation/extrudeToRegionMesh/Make/options
new file mode 100644
index 0000000000..dbf50a9b09
--- /dev/null
+++ b/applications/utilities/mesh/generation/extrudeToRegionMesh/Make/options
@@ -0,0 +1,14 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    /* -I$(LIB_SRC)/surfMesh/lnInclude */ \
+    /* -I$(LIB_SRC)/lagrangian/basic/lnInclude */ \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    /* -lsurfMesh */ \
+    /* -llagrangian */ \
+    -lmeshTools \
+    -ldynamicMesh
+
diff --git a/applications/utilities/mesh/generation/extrudeToRegionMesh/createShellMesh.C b/applications/utilities/mesh/generation/extrudeToRegionMesh/createShellMesh.C
new file mode 100644
index 0000000000..da1fc43066
--- /dev/null
+++ b/applications/utilities/mesh/generation/extrudeToRegionMesh/createShellMesh.C
@@ -0,0 +1,592 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "createShellMesh.H"
+#include "polyTopoChange.H"
+#include "meshTools.H"
+#include "mapPolyMesh.H"
+#include "polyAddPoint.H"
+#include "polyAddFace.H"
+#include "polyModifyFace.H"
+#include "polyAddCell.H"
+#include "patchPointEdgeCirculator.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(Foam::createShellMesh, 0);
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::createShellMesh::calcPointRegions
+(
+    const primitiveFacePatch& patch,
+    const PackedBoolList& nonManifoldEdge,
+    faceList& pointRegions,
+    labelList& regionPoints
+)
+{
+    pointRegions.setSize(patch.size());
+    forAll(pointRegions, faceI)
+    {
+        const face& f = patch.localFaces()[faceI];
+        pointRegions[faceI].setSize(f.size(), -1);
+    }
+
+    label nRegions = 0;
+
+    forAll(pointRegions, faceI)
+    {
+        const face& f = patch.localFaces()[faceI];
+
+        forAll(f, fp)
+        {
+            if (pointRegions[faceI][fp] == -1)
+            {
+                // Found unassigned point. Distribute current region.
+                label pointI = f[fp];
+                label edgeI = patch.faceEdges()[faceI][fp];
+
+                patchPointEdgeCirculator circ
+                (
+                    patch,
+                    nonManifoldEdge,
+                    edgeI,
+                    findIndex(patch.edgeFaces()[edgeI], faceI),
+                    pointI
+                );
+
+                for
+                (
+                    patchPointEdgeCirculator iter = circ.begin();
+                    iter != circ.end();
+                    ++iter
+                )
+                {
+                    label face2 = iter.faceID();
+
+                    if (face2 != -1)
+                    {
+                        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;
+                    }
+                }
+
+                nRegions++;
+            }
+        }
+    }
+
+
+    // From region back to originating point (many to one, a point might
+    // have multiple regions though)
+    regionPoints.setSize(nRegions);
+    forAll(pointRegions, faceI)
+    {
+        const face& f = patch.localFaces()[faceI];
+
+        forAll(f, fp)
+        {
+            regionPoints[pointRegions[faceI][fp]] = f[fp];
+        }
+    }
+
+
+    if (debug)
+    {
+        const labelListList& pointFaces = patch.pointFaces();
+        forAll(pointFaces, pointI)
+        {
+            label region = -1;
+            const labelList& pFaces = pointFaces[pointI];
+            forAll(pFaces, i)
+            {
+                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;
+
+                }
+            }
+        }
+    }
+}
+
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::createShellMesh::createShellMesh
+(
+    const primitiveFacePatch& patch,
+    const faceList& pointRegions,
+    const labelList& regionPoints
+)
+:
+    patch_(patch),
+    pointRegions_(pointRegions),
+    regionPoints_(regionPoints)
+{
+    if (pointRegions_.size() != patch_.size())
+    {
+        FatalErrorIn("createShellMesh::createShellMesh(..)")
+            << "nFaces:" << patch_.size()
+            << " pointRegions:" << pointRegions.size()
+            << exit(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+void Foam::createShellMesh::setRefinement
+(
+    const pointField& thickness,
+    const labelList& extrudeMasterPatchID,
+    const labelList& extrudeSlavePatchID,
+    const labelListList& extrudeEdgePatches,
+    polyTopoChange& meshMod
+)
+{
+    if (thickness.size() != regionPoints_.size())
+    {
+        FatalErrorIn("createShellMesh::setRefinement(..)")
+            << "nRegions:" << regionPoints_.size()
+            << " thickness:" << thickness.size()
+            << exit(FatalError);
+    }
+
+    if
+    (
+        extrudeMasterPatchID.size() != patch_.size()
+     && extrudeSlavePatchID.size() != patch_.size()
+    )
+    {
+        FatalErrorIn("createShellMesh::setRefinement(..)")
+            << "nFaces:" << patch_.size()
+            << " extrudeMasterPatchID:" << extrudeMasterPatchID.size()
+            << " extrudeSlavePatchID:" << extrudeSlavePatchID.size()
+            << exit(FatalError);
+    }
+
+    if (extrudeEdgePatches.size() != patch_.nEdges())
+    {
+        FatalErrorIn("createShellMesh::setRefinement(..)")
+            << "nEdges:" << patch_.nEdges()
+            << " extrudeEdgePatches:" << extrudeEdgePatches.size()
+            << exit(FatalError);
+    }
+
+
+
+    // From cell to patch (trivial)
+    DynamicList<label> cellToFaceMap(patch_.size());
+    // From face to patch+turning index
+    DynamicList<label> faceToFaceMap(2*patch_.size()+patch_.nEdges());
+    // From face to patch edge index
+    DynamicList<label> faceToEdgeMap(patch_.nEdges()+patch_.nEdges());
+    // From point to patch point index
+    DynamicList<label> pointToPointMap(2*patch_.nPoints());
+
+
+    // Introduce new cell for every face
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    labelList addedCells(patch_.size());
+    forAll(patch_, faceI)
+    {
+        addedCells[faceI] = meshMod.addCell
+        (
+            -1,                     // masterPointID
+            -1,                     // masterEdgeID
+            -1,                     // masterFaceID
+            cellToFaceMap.size(),   // masterCellID
+            -1                      // zoneID
+        );
+        cellToFaceMap.append(faceI);
+    }
+
+
+    // Introduce original points
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Original point numbers in local point ordering so no need to store.
+    forAll(patch_.localPoints(), pointI)
+    {
+        //label addedPointI =
+        meshMod.addPoint
+        (
+            patch_.localPoints()[pointI],   // point
+            pointToPointMap.size(),         // masterPointID 
+            -1,                             // zoneID
+            true                            // inCell
+        );
+        pointToPointMap.append(pointI);
+
+        //Pout<< "Added bottom point " << addedPointI
+        //    << " at " << patch_.localPoints()[pointI]
+        //    << "  from point " << pointI
+        //    << endl;
+    }
+
+
+    // Introduce new points (one for every region)
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    labelList addedPoints(regionPoints_.size());
+    forAll(regionPoints_, regionI)
+    {
+        label pointI = regionPoints_[regionI];
+        point extrudedPt = patch_.localPoints()[pointI] + thickness[regionI];
+
+        addedPoints[regionI] = meshMod.addPoint
+        (
+            extrudedPt,             // point
+            pointToPointMap.size(), // masterPointID - used only addressing
+            -1,                     // zoneID
+            true                    // inCell
+        );
+        pointToPointMap.append(pointI);
+
+        //Pout<< "Added top point " << addedPoints[regionI]
+        //    << " at " << extrudedPt
+        //    << "  from point " << pointI
+        //    << endl;
+    }
+
+
+    // Add face on patch' master side
+    //labelList masterFaces(patch_.localFaces().size());
+    forAll(patch_.localFaces(), faceI)
+    {
+        //masterFaces[faceI] =
+        meshMod.addFace
+        (
+            patch_.localFaces()[faceI].reverseFace(),// vertices
+            addedCells[faceI],          // own
+            -1,                         // nei
+            -1,                         // masterPointID
+            -1,                         // masterEdgeID
+            faceToFaceMap.size(),       // masterFaceID : current faceI
+            true,                       // flipFaceFlux
+            extrudeMasterPatchID[faceI],// patchID
+            -1,                         // zoneID
+            false                       // zoneFlip
+        );
+        faceToFaceMap.append(faceI+1);  // points to unflipped original face
+        faceToEdgeMap.append(-1);
+
+        //Pout<< "Added master face "
+        //    << patch_.localFaces()[faceI].reverseFace()
+        //    << " own " << addedCells[faceI]
+        //    << "  at " << patch_.faceCentres()[faceI]
+        //    << endl;
+
+    }
+
+    // Add face on patch' slave side
+    //labelList slaveFaces(patch_.localFaces().size());
+    forAll(patch_.localFaces(), faceI)
+    {
+        // Get face in original ordering
+        const face& f = patch_.localFaces()[faceI];
+
+        // Pick up point based on region
+        face newF(f.size());
+        forAll(f, fp)
+        {
+            label region = pointRegions_[faceI][fp];
+            newF[fp] = addedPoints[region];
+        }
+
+        //slaveFaces[faceI] =
+        meshMod.addFace
+        (
+            newF,                       // vertices
+            addedCells[faceI],          // own
+            -1,                         // nei
+            -1,                         // masterPointID
+            -1,                         // masterEdgeID
+            faceToFaceMap.size(),       // masterFaceID : current faceI
+            false,                      // flipFaceFlux
+            extrudeSlavePatchID[faceI], // patchID
+            -1,                         // zoneID
+            false                       // zoneFlip
+        );
+        faceToFaceMap.append(-faceI-1);
+        faceToEdgeMap.append(-1);
+
+        //Pout<< "Added slave face " << newF
+        //    << " own " << addedCells[faceI]
+        //    << "  at " << patch_.faceCentres()[faceI]
+        //    << endl;
+    }
+
+
+    // Add side faces
+    // ~~~~~~~~~~~~~~
+    // Note that we loop over edges multiple times so for edges with
+    // two cyclic faces they get added in two passes (for correct ordering)
+
+    // Pass1. Internal edges and first face of other edges
+    forAll(extrudeEdgePatches, edgeI)
+    {
+        const labelList& eFaces = patch_.edgeFaces()[edgeI];
+        const labelList& ePatches = extrudeEdgePatches[edgeI];
+
+        if (ePatches.size() == 0)
+        {
+            // internal face.
+
+            if (eFaces.size() != 2)
+            {
+                FatalErrorIn("createShellMesh::setRefinement(..)")
+                    << "edge:" << edgeI
+                    << " not internal but does not have side-patches defined."
+                    << exit(FatalError);
+            }
+
+            // Extrude
+
+            // Make face pointing in to eFaces[0] so out of new master face
+            const face& f = patch_.localFaces()[eFaces[0]];
+            const edge& e = patch_.edges()[edgeI];
+
+            label fp0 = findIndex(f, e[0]);
+            label fp1 = f.fcIndex(fp0);
+
+            if (f[fp1] != e[1])
+            {
+                fp1 = fp0;
+                fp0 = f.rcIndex(fp1);
+            }
+
+            face newF(4);
+            newF[0] = f[fp0];
+            newF[1] = f[fp1];
+            newF[2] = addedPoints[pointRegions_[eFaces[0]][fp1]];
+            newF[3] = addedPoints[pointRegions_[eFaces[0]][fp0]];
+
+            label minCellI = addedCells[eFaces[0]];
+            label maxCellI = addedCells[eFaces[1]];
+
+            if (minCellI > maxCellI)
+            {
+                // Swap
+                Swap(minCellI, maxCellI);
+                newF = newF.reverseFace();
+            }
+
+            //Pout<< "for internal edge:" << e
+            //    << " at:" << patch_.localPoints()[e[0]]
+            //    << patch_.localPoints()[e[1]]
+            //    << " adding face:" << newF
+            //    << " from f:" << f
+            //    << " inbetween " << minCellI << " and " << maxCellI << endl;
+
+            // newF already outwards pointing.
+            meshMod.addFace
+            (
+                newF,       // vertices
+                minCellI,   // own
+                maxCellI,   // nei
+                -1,         // masterPointID
+                -1,         // masterEdgeID
+                faceToFaceMap.size(),   // masterFaceID
+                false,      // flipFaceFlux
+                -1,         // patchID
+                -1,         // zoneID
+                false       // zoneFlip
+            );
+            faceToFaceMap.append(0);
+            faceToEdgeMap.append(edgeI);
+        }
+        else
+        {
+            if (eFaces.size() != ePatches.size())
+            {
+                FatalErrorIn("createShellMesh::setRefinement(..)")
+                    << "external/feature edge:" << edgeI
+                    << " has " << eFaces.size() << " connected extruded faces "
+                    << " but only " << ePatches.size()
+                    << " boundary faces defined." << exit(FatalError);
+            }
+
+
+            // Extrude eFaces[0]
+            label minFaceI = eFaces[0];
+
+            // Make face pointing in to eFaces[0] so out of new master face
+            const face& f = patch_.localFaces()[minFaceI];
+
+            const edge& e = patch_.edges()[edgeI];
+            label fp0 = findIndex(f, e[0]);
+            label fp1 = f.fcIndex(fp0);
+
+            if (f[fp1] != e[1])
+            {
+                fp1 = fp0;
+                fp0 = f.rcIndex(fp1);
+            }
+
+            face newF(4);
+            newF[0] = f[fp0];
+            newF[1] = f[fp1];
+            newF[2] = addedPoints[pointRegions_[minFaceI][fp1]];
+            newF[3] = addedPoints[pointRegions_[minFaceI][fp0]];
+
+
+            //Pout<< "for external edge:" << e
+            //    << " at:" << patch_.localPoints()[e[0]]
+            //    << patch_.localPoints()[e[1]]
+            //    << " adding first patch face:" << newF
+            //    << " from:" << f
+            //    << " into patch:" << ePatches[0]
+            //    << " own:" << addedCells[minFaceI]
+            //    << endl;
+
+
+            // newF already outwards pointing.
+            meshMod.addFace
+            (
+                newF,                   // vertices
+                addedCells[minFaceI],   // own
+                -1,                     // nei
+                -1,                     // masterPointID
+                -1,                     // masterEdgeID
+                faceToFaceMap.size(),   // masterFaceID
+                false,                  // flipFaceFlux
+                ePatches[0],            // patchID
+                -1,                     // zoneID
+                false                   // zoneFlip
+            );
+            faceToFaceMap.append(0);
+            faceToEdgeMap.append(edgeI);
+        }
+    }
+
+    // Pass2. Other faces of boundary edges
+    forAll(extrudeEdgePatches, edgeI)
+    {
+        const labelList& eFaces = patch_.edgeFaces()[edgeI];
+        const labelList& ePatches = extrudeEdgePatches[edgeI];
+
+        if (ePatches.size() >= 2)
+        {
+            for (label i = 1; i < ePatches.size(); i++)
+            {
+                // Extrude eFaces[i]
+                label minFaceI = eFaces[i];
+
+                // Make face pointing in to eFaces[0] so out of new master face
+                const face& f = patch_.localFaces()[minFaceI];
+
+                const edge& e = patch_.edges()[edgeI];
+                label fp0 = findIndex(f, e[0]);
+                label fp1 = f.fcIndex(fp0);
+
+                if (f[fp1] != e[1])
+                {
+                    fp1 = fp0;
+                    fp0 = f.rcIndex(fp1);
+                }
+
+                face newF(4);
+                newF[0] = f[fp0];
+                newF[1] = f[fp1];
+                newF[2] = addedPoints[pointRegions_[minFaceI][fp1]];
+                newF[3] = addedPoints[pointRegions_[minFaceI][fp0]];
+
+                //Pout<< "for external edge:" << e
+                //    << " at:" << patch_.localPoints()[e[0]]
+                //    << patch_.localPoints()[e[1]]
+                //    << " adding patch face:" << newF
+                //    << " from:" << f
+                //    << " into patch:" << ePatches[i]
+                //    << endl;
+
+                // newF already outwards pointing.
+                meshMod.addFace
+                (
+                    newF,                   // vertices
+                    addedCells[minFaceI],   // own
+                    -1,                     // nei
+                    -1,                     // masterPointID
+                    -1,                     // masterEdgeID
+                    faceToFaceMap.size(),   // masterFaceID
+                    false,                  // flipFaceFlux
+                    ePatches[i],            // patchID
+                    -1,                     // zoneID
+                    false                   // zoneFlip
+                );
+                faceToFaceMap.append(0);
+                faceToEdgeMap.append(edgeI);
+            }
+        }
+    }
+
+
+    cellToFaceMap_.transfer(cellToFaceMap);
+    faceToFaceMap_.transfer(faceToFaceMap);
+    faceToEdgeMap_.transfer(faceToEdgeMap);
+    pointToPointMap_.transfer(pointToPointMap);
+}
+
+
+void Foam::createShellMesh::updateMesh(const mapPolyMesh& map)
+{
+    inplaceReorder(map.reverseCellMap(), cellToFaceMap_);
+    inplaceReorder(map.reverseFaceMap(), faceToFaceMap_);
+    inplaceReorder(map.reverseFaceMap(), faceToEdgeMap_);
+    inplaceReorder(map.reversePointMap(), pointToPointMap_);
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/extrudeToRegionMesh/createShellMesh.H b/applications/utilities/mesh/generation/extrudeToRegionMesh/createShellMesh.H
new file mode 100644
index 0000000000..7823f11459
--- /dev/null
+++ b/applications/utilities/mesh/generation/extrudeToRegionMesh/createShellMesh.H
@@ -0,0 +1,177 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::createShellMesh
+
+Description
+    Creates mesh by extruding a patch.
+
+SourceFiles
+    createShellMesh.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef createShellMesh_H
+#define createShellMesh_H
+
+#include "primitiveFacePatch.H"
+#include "PackedBoolList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class mapPolyMesh;
+class polyTopoChange;
+
+/*---------------------------------------------------------------------------*\
+                           Class createShellMesh Declaration
+\*---------------------------------------------------------------------------*/
+
+class createShellMesh
+{
+    // Private data
+
+        //- Reference to patch to extrude
+        const primitiveFacePatch& patch_;
+
+//        //- Non-manifold edges on patch
+//        const PackedBoolList& nonManifoldEdge_;
+
+        //- Per point on face the region
+        const faceList& pointRegions_;
+
+        //- From region to original patch point
+        const labelList& regionPoints_;
+
+
+        labelList cellToFaceMap_;
+
+        labelList faceToFaceMap_;
+
+        labelList faceToEdgeMap_;
+
+        labelList pointToPointMap_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        createShellMesh(const createShellMesh&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const createShellMesh&);
+
+
+public:
+
+    //- Runtime type information
+    ClassName("createShellMesh");
+
+
+    // Constructors
+
+        //- Construct from mesh.
+        createShellMesh
+        (
+            const primitiveFacePatch& patch,
+//            const PackedBoolList& nonManifoldEdge,
+            const faceList& pointRegions,
+            const labelList& regionPoints
+        );
+
+
+    // Member Functions
+
+
+        // Access
+
+            //- From region cell to patch face
+            const labelList& cellToFaceMap() const
+            {
+                return cellToFaceMap_;
+            }
+
+            //- From region face to patch face. Contains turning index:
+            //      > 0 : face originating from master side
+            //            (owner side for unflipped face)
+            //      < 0 : (bottom) face originating from slave side
+            //            (neighbour side for unflipped face)
+            //      = 0 : for all side faces
+            const labelList& faceToFaceMap() const
+            {
+                return faceToFaceMap_;
+            }
+
+            //- From region side-face to patch edge. -1 for non-edge faces.
+            const labelList& faceToEdgeMap() const
+            {
+                return faceToEdgeMap_;
+            }
+
+            //- From region point to patch point.
+            const labelList& pointToPointMap() const
+            {
+                return pointToPointMap_;
+            }
+
+
+        // Other
+
+            //- Helper: calculate point regions
+            static void calcPointRegions
+            (
+                const primitiveFacePatch& patch,
+                const PackedBoolList& nonManifoldEdge,
+                faceList& pointRegions,
+                labelList& regionPoints
+            );
+
+            //- Play commands into polyTopoChange to create layer mesh.
+            void setRefinement
+            (
+                const pointField& thickness,
+                const labelList& extrudeMasterPatchID,
+                const labelList& extrudeSlavePatchID,
+                const labelListList& extrudeEdgePatches,
+                polyTopoChange& meshMod
+            );
+
+            //- Update any locally stored mesh information.
+            void updateMesh(const mapPolyMesh&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/extrudeToRegionMesh/extrudeToRegionMesh.C b/applications/utilities/mesh/generation/extrudeToRegionMesh/extrudeToRegionMesh.C
new file mode 100644
index 0000000000..1c81c07650
--- /dev/null
+++ b/applications/utilities/mesh/generation/extrudeToRegionMesh/extrudeToRegionMesh.C
@@ -0,0 +1,1334 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    Extrude faceZones into separate mesh (as a different region).
+
+    - used to e.g. extrude baffles (extrude internal faces) or create
+    liquid film regions.
+    - if extruding internal faces:
+        - create baffles in original mesh with directMappedWall patches
+    - if extruding boundary faces:
+        - convert boundary faces to directMappedWall patches
+    - extrude edges of faceZone as a <zone>_sidePatch
+    - extrude edges inbetween different faceZones as a cyclic <zoneA>_<zoneB>
+    - extrudes into master direction (i.e. away from the owner cell
+      if flipMap is false)
+    - not parallel
+
+Usage
+
+    - extrudeToRegionMesh <regionName> <faceZones> <thickness>
+
+    @param \<regionName\> \n
+    Name of mesh to create.
+
+    @param \<faceZones\> \n
+    List of faceZones to extrude
+
+    @param \<thickness\> \n
+    Thickness of extruded mesh.
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "fvMesh.H"
+#include "polyTopoChange.H"
+#include "patchPointEdgeCirculator.H"
+#include "OFstream.H"
+#include "meshTools.H"
+#include "directMappedWallPolyPatch.H"
+#include "createShellMesh.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+#include "cyclicPolyPatch.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class GeoField>
+void addPatchFields(fvMesh& mesh, const word& patchFieldType)
+{
+    HashTable<const GeoField*> flds
+    (
+        mesh.objectRegistry::lookupClass<GeoField>()
+    );
+
+    for
+    (
+        typename HashTable<const GeoField*>::const_iterator iter = flds.begin();
+        iter != flds.end();
+        ++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()
+            )
+        );
+    }
+}
+
+
+// Remove last patch field
+template<class GeoField>
+void trimPatchFields(fvMesh& mesh, const label nPatches)
+{
+    HashTable<const GeoField*> flds
+    (
+        mesh.objectRegistry::lookupClass<GeoField>()
+    );
+
+    for
+    (
+        typename HashTable<const GeoField*>::const_iterator iter = flds.begin();
+        iter != flds.end();
+        ++iter
+    )
+    {
+        const GeoField& fld = *iter();
+
+        const_cast<typename GeoField::GeometricBoundaryField&>
+        (
+            fld.boundaryField()
+        ).setSize(nPatches);
+    }
+}
+
+
+// Reorder patch field
+template<class GeoField>
+void reorderPatchFields(fvMesh& mesh, const labelList& oldToNew)
+{
+    HashTable<const GeoField*> flds
+    (
+        mesh.objectRegistry::lookupClass<GeoField>()
+    );
+
+    for
+    (
+        typename HashTable<const GeoField*>::const_iterator iter = flds.begin();
+        iter != flds.end();
+        ++iter
+    )
+    {
+        const GeoField& fld = *iter();
+
+        typename GeoField::GeometricBoundaryField& bfld =
+            const_cast<typename GeoField::GeometricBoundaryField&>
+            (
+                fld.boundaryField()
+            );
+
+        bfld.reorder(oldToNew);
+    }
+}
+
+
+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)
+{
+    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;
+}
+
+
+// Reorder and delete patches.
+void reorderPatches
+(
+    fvMesh& mesh,
+    const labelList& oldToNew,
+    const label nNewPatches
+)
+{
+    polyBoundaryMesh& polyPatches =
+        const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
+    fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
+
+    // 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);
+
+    // Remove last.
+    polyPatches.setSize(nNewPatches);
+    fvPatches.setSize(nNewPatches);
+    trimPatchFields<volScalarField>(mesh, nNewPatches);
+    trimPatchFields<volVectorField>(mesh, nNewPatches);
+    trimPatchFields<volSphericalTensorField>(mesh, nNewPatches);
+    trimPatchFields<volSymmTensorField>(mesh, nNewPatches);
+    trimPatchFields<volTensorField>(mesh, nNewPatches);
+    trimPatchFields<surfaceScalarField>(mesh, nNewPatches);
+    trimPatchFields<surfaceVectorField>(mesh, nNewPatches);
+    trimPatchFields<surfaceSphericalTensorField>(mesh, nNewPatches);
+    trimPatchFields<surfaceSymmTensorField>(mesh, nNewPatches);
+    trimPatchFields<surfaceTensorField>(mesh, nNewPatches);
+}
+
+
+void createDummyFvMeshFiles(const polyMesh& mesh, const word& regionName)
+{
+    // Create dummy system/fv*
+    {
+        IOobject io
+        (
+            "fvSchemes",
+            mesh.time().system(),
+            regionName,
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        );
+
+        Info<< "Testing:" << io.objectPath() << endl;
+
+        if (!io.headerOk())
+        {
+            Info<< "Writing dummy " << regionName/io.name() << endl;
+            dictionary dummyDict;
+            dictionary divDict;
+            dummyDict.add("divSchemes", divDict);
+            dictionary gradDict;
+            dummyDict.add("gradSchemes", gradDict);
+            dictionary laplDict;
+            dummyDict.add("laplacianSchemes", laplDict);
+
+            IOdictionary(io, dummyDict).regIOobject::write();
+        }
+    }
+    {
+        IOobject io
+        (
+            "fvSolution",
+            mesh.time().system(),
+            regionName,
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        );
+
+        if (!io.headerOk())
+        {
+            Info<< "Writing dummy " << regionName/io.name() << endl;
+            dictionary dummyDict;
+            IOdictionary(io, dummyDict).regIOobject::write();
+        }
+    }
+}
+
+
+tmp<pointField> calcOffset
+(
+    const primitiveFacePatch& extrudePatch,
+    const createShellMesh& extruder,
+    const polyPatch& pp
+)
+{
+    vectorField::subField fc = pp.faceCentres();
+
+    tmp<pointField> toffsets(new pointField(fc.size()));
+    pointField& offsets = toffsets();
+
+    forAll(fc, i)
+    {
+        label meshFaceI = pp.start()+i;
+        label patchFaceI = mag(extruder.faceToFaceMap()[meshFaceI])-1;
+        point patchFc = extrudePatch[patchFaceI].centre
+        (
+            extrudePatch.points()
+        );
+        offsets[i] = patchFc - fc[i];
+    }
+    return toffsets;
+}
+
+
+
+// Main program:
+
+int main(int argc, char *argv[])
+{
+    argList::noParallel();
+    argList::validArgs.append("shellRegion");
+    argList::validArgs.append("faceZones");
+    argList::validArgs.append("thickness");
+
+    #include "addRegionOption.H"
+    argList::validOptions.insert("overwrite", "");
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createNamedMesh.H"
+    const word oldInstance = mesh.pointsInstance();
+
+    word shellRegionName = args.additionalArgs()[0];
+    const wordList zoneNames(IStringStream(args.additionalArgs()[1])());
+    scalar thickness = readScalar(IStringStream(args.additionalArgs()[2])());
+    bool overwrite = args.optionFound("overwrite");
+
+
+    Info<< "Extruding zones " << zoneNames
+        << " on mesh " << regionName
+        << " into shell mesh " << shellRegionName
+        << " of thickness " << thickness << nl
+        << endl;
+
+    if (shellRegionName == regionName)
+    {
+        FatalErrorIn(args.executable())
+            << "Cannot extrude into same region as mesh." << endl
+            << "Mesh region : " << regionName << endl
+            << "Shell region : " << shellRegionName
+            << exit(FatalError);
+    }
+
+
+    // Add interface patches
+    // Note: should be one per faceZone?
+    word masterPatchName = regionName + "_to_" + shellRegionName + "_master";
+    label masterPatchID = addPatch<directMappedWallPolyPatch>
+    (
+        mesh,
+        masterPatchName
+    );
+    word slavePatchName = regionName + "_to_" + shellRegionName + "_slave";
+    label slavePatchID = addPatch<directMappedWallPolyPatch>
+    (
+        mesh,
+        slavePatchName
+    );
+
+
+
+
+    // Create dummy fv* files
+    createDummyFvMeshFiles(mesh, shellRegionName);
+
+
+    word meshInstance;
+    if (!overwrite)
+    {
+        runTime++;
+        meshInstance = runTime.timeName();
+    }
+    else
+    {
+        meshInstance = oldInstance;
+    }
+    Info<< "Writing meshes to " << meshInstance << nl << endl;
+
+
+
+    // Collect faces to extrude
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+    const faceZoneMesh& faceZones = mesh.faceZones();
+
+    label nExtrudeFaces = 0;
+    forAll(zoneNames, i)
+    {
+        label zoneI = faceZones.findZoneID(zoneNames[i]);
+        if (zoneI == -1)
+        {
+            FatalErrorIn(args.executable())
+                << "Cannot find zone " << zoneNames[i] << endl
+                << "Valid zones are " << faceZones.names()
+                << exit(FatalError);
+        }
+        nExtrudeFaces += faceZones[zoneI].size();
+    }
+    labelList extrudeMeshFaces(nExtrudeFaces);
+    faceList zoneFaces(nExtrudeFaces);
+    labelList zoneID(nExtrudeFaces);
+    boolList zoneFlipMap(nExtrudeFaces);
+    nExtrudeFaces = 0;
+    forAll(zoneNames, i)
+    {
+        label zoneI = faceZones.findZoneID(zoneNames[i]);
+        const faceZone& fz = faceZones[zoneI];
+        const primitiveFacePatch& fzp = fz();
+        forAll(fz, j)
+        {
+            extrudeMeshFaces[nExtrudeFaces] = fz[j];
+            zoneFaces[nExtrudeFaces] = fzp[j];
+            zoneID[nExtrudeFaces] = zoneI;
+            zoneFlipMap[nExtrudeFaces] = fz.flipMap()[j];
+            nExtrudeFaces++;
+        }
+    }
+    primitiveFacePatch extrudePatch(zoneFaces.xfer(), mesh.points());
+    const pointField& extrudePoints = extrudePatch.localPoints();
+    const faceList& extrudeFaces = extrudePatch.localFaces();
+    const labelListList& edgeFaces = extrudePatch.edgeFaces();
+
+    labelList extrudeSlavePatchID(extrudePatch.size());
+    labelList extrudeMasterPatchID(extrudePatch.size());
+    forAll(extrudeMeshFaces, i)
+    {
+        extrudeSlavePatchID[i] = slavePatchID;
+        extrudeMasterPatchID[i] = masterPatchID;
+    }
+
+    Info<< "extrudePatch :"
+        << " faces:" << extrudePatch.size()
+        << " points:" << extrudePatch.nPoints()
+        << " edges:" << extrudePatch.nEdges()
+        << nl
+        << endl;
+
+
+    // Count how many patches on special edges of extrudePatch are necessary
+    // - zoneXXX_sides
+    // - zoneXXX_zoneYYY
+    labelList zoneSidePatch(faceZones.size(), 0);
+    labelList zoneZonePatch(faceZones.size()*faceZones.size(), 0);
+
+    forAll(edgeFaces, edgeI)
+    {
+        const labelList& eFaces = edgeFaces[edgeI];
+        if (eFaces.size() == 2)
+        {
+            label zone0 = zoneID[eFaces[0]];
+            label zone1 = zoneID[eFaces[1]];
+
+            if (zone0 != zone1)
+            {
+                label minZone = min(zone0,zone1);
+                label maxZone = max(zone0,zone1);
+                zoneZonePatch[minZone*faceZones.size()+maxZone]++;
+            }
+        }
+        else
+        {
+            // Determine the min zone of all connected zones.
+            label minZone = zoneID[eFaces[0]];
+            for (label i = 1; i < eFaces.size(); i++)
+            {
+                minZone = min(minZone, zoneID[eFaces[i]]);
+            }
+            zoneSidePatch[minZone]++;
+        }
+    }
+    Pstream::listCombineGather(zoneSidePatch, plusEqOp<label>());
+    Pstream::listCombineScatter(zoneSidePatch);
+    Pstream::listCombineGather(zoneZonePatch, plusEqOp<label>());
+    Pstream::listCombineScatter(zoneZonePatch);
+
+    // Now check which patches to add.
+    Info<< "Adding patches for edges on zones:" << nl
+        << "\tpatchID\tpatch" << nl
+        << "\t-------\t-----"
+        << endl;
+
+    label nSide = 0;
+    forAll(zoneSidePatch, zoneI)
+    {
+        if (zoneSidePatch[zoneI] > 0)
+        {
+            word patchName = faceZones[zoneI].name() + "_" + "sidePatch";
+            zoneSidePatch[zoneI] = addPatch<polyPatch>
+            (
+                mesh,
+                patchName
+            );
+
+            Info<< '\t' << zoneSidePatch[zoneI]
+                << '\t' << patchName << nl;
+
+            nSide++;
+        }
+    }
+    Info<< "Added " << nSide << " zone-edge patches." << nl
+        << endl;
+
+
+
+    Info<< "Adding inter-zone patches:" << nl
+        << "\tpatchID\tpatch" << nl
+        << "\t-------\t-----"
+        << endl;
+
+    dictionary transformDict;
+    transformDict.add
+    (
+        "transform",
+        cyclicPolyPatch::transformTypeNames[cyclicPolyPatch::NOORDERING]
+    );
+
+    label nInter = 0;
+    forAll(zoneZonePatch, minZone)
+    {
+        for (label maxZone = minZone; maxZone < faceZones.size(); maxZone++)
+        {
+            label index = minZone*faceZones.size()+maxZone;
+
+            if (zoneZonePatch[index] > 0)
+            {
+                word patchName =
+                    faceZones[minZone].name()
+                  + "_"
+                  + faceZones[maxZone].name();
+
+                zoneZonePatch[index] = addPatch<cyclicPolyPatch>
+                (
+                    mesh,
+                    patchName,
+                    transformDict
+                );
+
+                Info<< '\t' << zoneZonePatch[index]
+                    << '\t' << patchName << nl;
+
+                nInter++;
+            }
+        }
+    }
+    Info<< "Added " << nInter << " inter-zone patches." << nl
+        << endl;
+
+
+
+    // Set patches to use for edges to be extruded into boundary faces
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // In order of edgeFaces: per edge, per originating face the
+    // patch to use for the side face (from the extruded edge).
+    // If empty size create an internal face.
+    labelListList extrudeEdgePatches(extrudePatch.nEdges());
+
+    // Is edge an non-manifold edge
+    PackedBoolList nonManifoldEdge(extrudePatch.nEdges());
+
+    forAll(edgeFaces, edgeI)
+    {
+        const labelList& eFaces = edgeFaces[edgeI];
+
+        labelList& ePatches = extrudeEdgePatches[edgeI];
+
+        if (eFaces.size() == 2)
+        {
+            label zone0 = zoneID[eFaces[0]];
+            label zone1 = zoneID[eFaces[1]];
+
+            if (zone0 != zone1) // || (cos(angle) > blabla))
+            {
+                label minZone = min(zone0,zone1);
+                label maxZone = max(zone0,zone1);
+                label patchI = zoneZonePatch[minZone*faceZones.size()+maxZone];
+                ePatches.setSize(eFaces.size());
+                ePatches = patchI;
+
+                nonManifoldEdge[edgeI] = 1;
+            }
+        }
+        else
+        {
+            ePatches.setSize(eFaces.size());
+            forAll(eFaces, i)
+            {
+                ePatches[i] = zoneSidePatch[zoneID[eFaces[i]]];
+            }
+            nonManifoldEdge[edgeI] = 1;
+        }
+    }
+
+    // Override constraint types
+    {
+        const edgeList& extrudeEdges = extrudePatch.edges();
+        const labelList& extrudeMeshPts = extrudePatch.meshPoints();
+
+        // Map from mesh edge to local patch edge index
+        EdgeMap<label> extrudeMeshEdges(extrudePatch.nEdges());
+
+        forAll(extrudeEdges, edgeI)
+        {
+            if (edgeFaces[edgeI].size() == 1)
+            {
+                const edge& e = extrudeEdges[edgeI];
+                const edge meshE(extrudeMeshPts[e[0]], extrudeMeshPts[e[1]]);
+                extrudeMeshEdges.insert(meshE, edgeI);
+            }
+        }
+        forAll(patches, patchI)
+        {
+            const polyPatch& pp = patches[patchI];
+
+            if (polyPatch::constraintType(pp.type()))
+            {
+                const edgeList& edges = pp.edges();
+
+                forAll(edges, ppEdgeI)
+                {
+                    const edge& e = edges[ppEdgeI];
+                    edge meshE(pp.meshPoints()[e[0]], pp.meshPoints()[e[1]]);
+                    EdgeMap<label>::const_iterator iter = extrudeMeshEdges.find
+                    (
+                        meshE
+                    );
+                    if (iter != extrudeMeshEdges.end())
+                    {
+                        label extrudeEdgeI = iter();
+                        extrudeEdgePatches[extrudeEdgeI] = labelList
+                        (
+                            edgeFaces[extrudeEdgeI].size(),
+                            patchI
+                        );
+                    }
+                }
+            }
+        }
+    }
+
+
+
+
+    // Assign point regions
+    // ~~~~~~~~~~~~~~~~~~~~
+
+    // Per face, per point the region number.
+    faceList pointRegions(extrudePatch.size());
+    // Per region the originating point
+    labelList regionPoints;
+
+    createShellMesh::calcPointRegions
+    (
+        extrudePatch,
+        nonManifoldEdge,
+        pointRegions,
+        regionPoints
+    );
+
+
+    // Calculate a normal per region
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    vectorField regionNormals(regionPoints.size(), vector::zero);
+
+    forAll(pointRegions, faceI)
+    {
+        const face& f = extrudeFaces[faceI];
+
+        forAll(f, fp)
+        {
+            label region = pointRegions[faceI][fp];
+            regionNormals[region] += extrudePatch.faceNormals()[faceI];
+        }
+    }
+    regionNormals /= mag(regionNormals);
+
+    // For debugging: dump hedgehog plot of normals
+    {
+        OFstream str(runTime.path()/"regionNormals.obj");
+        label vertI = 0;
+
+        forAll(pointRegions, faceI)
+        {
+            const face& f = extrudeFaces[faceI];
+
+            forAll(f, fp)
+            {
+                label region = pointRegions[faceI][fp];
+                const point& pt = extrudePoints[f[fp]];
+
+                meshTools::writeOBJ(str, pt);
+                vertI++;
+                meshTools::writeOBJ(str, pt+0.01*regionNormals[region]);
+                vertI++;
+                str << "l " << vertI-1 << ' ' << vertI << nl;
+            }
+        }
+    }
+
+
+
+    // Create a new mesh
+    // ~~~~~~~~~~~~~~~~~
+
+    createShellMesh extruder
+    (
+        extrudePatch,
+        pointRegions,
+        regionPoints
+    );
+
+
+    autoPtr<fvMesh> regionMeshPtr;
+    autoPtr<mapPolyMesh> shellMap;
+    {
+        polyTopoChange meshMod(mesh.boundaryMesh().size());
+
+        extruder.setRefinement
+        (
+            thickness*regionNormals,
+            extrudeMasterPatchID,
+            extrudeSlavePatchID,
+            extrudeEdgePatches,
+            meshMod
+        );
+
+        shellMap = meshMod.makeMesh
+        (
+            regionMeshPtr,     // mesh to create
+            IOobject
+            (
+                shellRegionName,
+                meshInstance,
+                runTime,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            mesh            // mesh to get original patch info from
+        );
+    }
+    fvMesh& regionMesh = regionMeshPtr();
+
+    // Necessary?
+    regionMesh.setInstance(meshInstance);
+
+
+    // Update numbering on extruder.
+    extruder.updateMesh(shellMap);
+
+
+    // Change master and slave boundary conditions on regionMesh
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Save offsets from shell mesh back to original mesh
+    pointField masterOffsets, slaveOffsets;
+
+    {
+        const polyBoundaryMesh& regionPatches = regionMesh.boundaryMesh();
+        List<polyPatch*> newPatches(regionPatches.size());
+        forAll(regionPatches, patchI)
+        {
+            const polyPatch& pp = regionPatches[patchI];
+            if (pp.name() == masterPatchName)
+            {
+                masterOffsets = calcOffset(extrudePatch, extruder, pp);
+
+                newPatches[patchI] = new directMappedWallPolyPatch
+                (
+                    pp.name(),
+                    pp.size(),
+                    pp.start(),
+                    patchI,
+                    regionName,                                 // sampleRegion
+                    directMappedPatchBase::NEARESTPATCHFACE,    // sampleMode
+                    pp.name(),                                  // samplePatch
+                    masterOffsets,                              // offset
+                    patches
+                );
+            }
+            else if (pp.name() == slavePatchName)
+            {
+                slaveOffsets = calcOffset(extrudePatch, extruder, pp);
+
+                newPatches[patchI] = new directMappedWallPolyPatch
+                (
+                    pp.name(),
+                    pp.size(),
+                    pp.start(),
+                    patchI,
+                    regionName,                                 // sampleRegion
+                    directMappedPatchBase::NEARESTPATCHFACE,    // sampleMode
+                    pp.name(),                                  // samplePatch
+                    slaveOffsets,                               // offset
+                    patches
+                );
+            }
+            else
+            {
+                newPatches[patchI] = pp.clone
+                (
+                    regionPatches,
+                    patchI,
+                    pp.size(),
+                    pp.start()
+                ).ptr();
+            }
+        }
+        regionMesh.removeFvBoundary();
+        regionMesh.addFvPatches(newPatches, true);
+    }
+
+
+
+    // Write addressing from region mesh back to originating patch
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    labelIOList cellToPatchFaceAddressing
+    (
+        IOobject
+        (
+            "cellToPatchFaceAddressing",
+            regionMesh.facesInstance(),
+            regionMesh.meshSubDir,
+            regionMesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        extruder.cellToFaceMap()
+    );
+    cellToPatchFaceAddressing.note() = "cell to patch face addressing";
+
+    labelIOList faceToPatchFaceAddressing
+    (
+        IOobject
+        (
+            "faceToPatchFaceAddressing",
+            regionMesh.facesInstance(),
+            regionMesh.meshSubDir,
+            regionMesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        extruder.faceToFaceMap()
+    );
+    faceToPatchFaceAddressing.note() =
+        "front/back face + turning index to patch face addressing";
+
+    labelIOList faceToPatchEdgeAddressing
+    (
+        IOobject
+        (
+            "faceToPatchEdgeAddressing",
+            regionMesh.facesInstance(),
+            regionMesh.meshSubDir,
+            regionMesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        extruder.faceToEdgeMap()
+    );
+    faceToPatchEdgeAddressing.note() =
+        "side face to patch edge addressing";
+
+    labelIOList pointToPatchPointAddressing
+    (
+        IOobject
+        (
+            "pointToPatchPointAddressing",
+            regionMesh.facesInstance(),
+            regionMesh.meshSubDir,
+            regionMesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        extruder.pointToPointMap()
+    );
+    pointToPatchPointAddressing.note() =
+        "point to patch point addressing";
+
+
+    Info<< "Writing mesh " << regionMesh.name()
+        << " to " << regionMesh.facesInstance() << nl
+        << endl;
+
+    bool ok =
+        regionMesh.write()
+     && cellToPatchFaceAddressing.write()
+     && faceToPatchFaceAddressing.write()
+     && faceToPatchEdgeAddressing.write()
+     && pointToPatchPointAddressing.write();
+
+    if (!ok)
+    {
+        FatalErrorIn(args.executable())
+            << "Failed writing mesh " << regionMesh.name()
+            << " at location " << regionMesh.facesInstance()
+            << exit(FatalError);
+    }
+
+
+
+
+    // Insert baffles into original mesh
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    autoPtr<mapPolyMesh> addBafflesMap;
+    {
+        polyTopoChange meshMod(mesh);
+
+        // Modify faces to be in master patch
+        forAll(extrudeMeshFaces, zoneFaceI)
+        {
+            label meshFaceI = extrudeMeshFaces[zoneFaceI];
+            label zoneI = zoneID[zoneFaceI];
+            bool flip = zoneFlipMap[zoneFaceI];
+            const face& f = mesh.faces()[meshFaceI];
+
+            if (!flip)
+            {
+                meshMod.modifyFace
+                (
+                    f,                          // modified face
+                    meshFaceI,                  // label of face being modified
+                    mesh.faceOwner()[meshFaceI],// owner
+                    -1,                         // neighbour
+                    false,                      // face flip
+                    masterPatchID,              // patch for face
+                    zoneI,                      // zone for face
+                    flip                        // face flip in zone
+                );
+            }
+            else if (mesh.isInternalFace(meshFaceI))
+            {
+                meshMod.modifyFace
+                (
+                    f.reverseFace(),                // modified face
+                    meshFaceI,                      // label of modified face
+                    mesh.faceNeighbour()[meshFaceI],// owner
+                    -1,                             // neighbour
+                    true,                           // face flip
+                    masterPatchID,                  // patch for face
+                    zoneI,                          // zone for face
+                    !flip                           // face flip in zone
+                );
+            }
+        }
+
+        // Add faces (using same points) to be in slave patch
+        forAll(extrudeMeshFaces, zoneFaceI)
+        {
+            label meshFaceI = extrudeMeshFaces[zoneFaceI];
+            bool flip = zoneFlipMap[zoneFaceI];
+            const face& f = mesh.faces()[meshFaceI];
+
+            if (!flip)
+            {
+                if (mesh.isInternalFace(meshFaceI))
+                {
+                    meshMod.addFace
+                    (
+                        f.reverseFace(),                // modified face
+                        mesh.faceNeighbour()[meshFaceI],// owner
+                        -1,                             // neighbour
+                        -1,                             // master point
+                        -1,                             // master edge
+                        meshFaceI,                      // master face
+                        true,                           // flip flux
+                        slavePatchID,                   // patch for face
+                        -1,                             // zone for face
+                        false                           // face flip in zone
+                    );
+                }
+            }
+            else
+            {
+                meshMod.addFace
+                (
+                    f,                              // face
+                    mesh.faceOwner()[meshFaceI],    // owner
+                    -1,                             // neighbour
+                    -1,                             // master point
+                    -1,                             // master edge
+                    meshFaceI,                      // master face
+                    false,                          // flip flux
+                    slavePatchID,                   // patch for face
+                    -1,                             // zone for face
+                    false                           // zone flip
+                );
+            }
+        }
+
+        // Change the mesh. Change points directly (no inflation).
+        addBafflesMap = meshMod.changeMesh(mesh, false);
+
+        // Update fields
+        mesh.updateMesh(addBafflesMap);
+
+
+//XXXXXX
+// Update maps! e.g. faceToPatchFaceAddressing
+//XXXXXX
+
+        // Move mesh (since morphing might not do this)
+        if (addBafflesMap().hasMotionPoints())
+        {
+            mesh.movePoints(addBafflesMap().preMotionPoints());
+        }
+
+        mesh.setInstance(meshInstance);
+    }
+    
+
+
+    // Change master and slave boundary conditions on originating mesh
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    {
+        const polyBoundaryMesh& patches = mesh.boundaryMesh();
+        List<polyPatch*> newPatches(patches.size());
+        forAll(patches, patchI)
+        {
+            const polyPatch& pp = patches[patchI];
+            if (pp.name() == masterPatchName)
+            {
+                newPatches[patchI] = new directMappedWallPolyPatch
+                (
+                    pp.name(),
+                    pp.size(),
+                    pp.start(),
+                    patchI,
+                    shellRegionName,                            // sampleRegion
+                    directMappedPatchBase::NEARESTPATCHFACE,    // sampleMode
+                    pp.name(),                                  // samplePatch
+                    -masterOffsets,                             // offset
+                    patches
+                );
+            }
+            else if (pp.name() == slavePatchName)
+            {
+                newPatches[patchI] = new directMappedWallPolyPatch
+                (
+                    pp.name(),
+                    pp.size(),
+                    pp.start(),
+                    patchI,
+                    shellRegionName,                            // sampleRegion
+                    directMappedPatchBase::NEARESTPATCHFACE,    // sampleMode
+                    pp.name(),                                  // samplePatch
+                    -slaveOffsets,                              // offset
+                    patches
+                );
+            }
+            else
+            {
+                newPatches[patchI] = pp.clone
+                (
+                    patches,
+                    patchI,
+                    pp.size(),
+                    pp.start()
+                ).ptr();
+            }
+        }
+        mesh.removeFvBoundary();
+        mesh.addFvPatches(newPatches, true);
+    }
+
+
+    Info<< "Writing mesh " << mesh.name()
+        << " to " << mesh.facesInstance() << nl
+        << endl;
+
+    if (!mesh.write())
+    {
+        FatalErrorIn(args.executable())
+            << "Failed writing mesh " << mesh.name()
+            << " at location " << mesh.facesInstance()
+            << exit(FatalError);
+    }
+
+
+    Info << "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/extrudeToRegionMesh/patchPointEdgeCirculator.C b/applications/utilities/mesh/generation/extrudeToRegionMesh/patchPointEdgeCirculator.C
new file mode 100644
index 0000000000..12d5744a91
--- /dev/null
+++ b/applications/utilities/mesh/generation/extrudeToRegionMesh/patchPointEdgeCirculator.C
@@ -0,0 +1,68 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#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/extrudeToRegionMesh/patchPointEdgeCirculator.H b/applications/utilities/mesh/generation/extrudeToRegionMesh/patchPointEdgeCirculator.H
new file mode 100644
index 0000000000..9cfcdbda5e
--- /dev/null
+++ b/applications/utilities/mesh/generation/extrudeToRegionMesh/patchPointEdgeCirculator.H
@@ -0,0 +1,206 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+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/extrudeToRegionMesh/patchPointEdgeCirculatorI.H b/applications/utilities/mesh/generation/extrudeToRegionMesh/patchPointEdgeCirculatorI.H
new file mode 100644
index 0000000000..7ea7abc1a9
--- /dev/null
+++ b/applications/utilities/mesh/generation/extrudeToRegionMesh/patchPointEdgeCirculatorI.H
@@ -0,0 +1,309 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+
+// * * * * * * * * * * * * * 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::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;
+}
+
+
+// ************************************************************************* //