diff --git a/ReleaseNotes-dev b/ReleaseNotes-dev
index 227041536f1268daf4b3ab5b1356ea9d9fe131d0..230a1168263fc7d0fe9ad64a916a5b3ef9557263 100644
--- a/ReleaseNotes-dev
+++ b/ReleaseNotes-dev
@@ -104,6 +104,11 @@
       taking film into account
     + Parallel aware
 *** *New* ptscotch decomposition method
+*** *Updated* scotch decomposition method to run in parallel by doing
+    decomposition on the master. Unfortunately scotch and ptscotch cannot
+    be linked in to the same executable.
+*** *Updated* simple decomposition method to run in parallel by doing
+    decomposition on the master.
 *** *Updated* decomposePar maps polyPatches instead of recreating them so
     polyPatches holding data can map the data.
 *** *Updated* particle tracking algorithm
@@ -183,6 +188,9 @@
       e.g. pitzDailyDirectMapped tutorial.
     + =setSet=: allows time range (e.g. 0:100) in combination with -batch argument
     to execute the commands for multiple times.
+    + =extrudeMesh=:
+            - option to add extrusion to existing mesh.
+            - works in parallel
 * Post-processing
   + =paraFoam=, =foamToVTK=: full support for polyhedral cell type in recent
      Paraview versions.
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H
index cac8e52a8efb7d6a1b9115afc595e4b4617f88eb..8f1ad9df8e3387677a5e4d1b9f61e88236dd3524 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H
@@ -2,8 +2,8 @@
     word scheme("div(phi,alpha)");
     word schemer("div(phir,alpha)");
 
-    surfaceScalarField phic(phi);
-    surfaceScalarField phir(phia - phib);
+    surfaceScalarField phic("phic", phi);
+    surfaceScalarField phir("phir", phia - phib);
 
     if (g0.value() > 0.0)
     {
@@ -41,7 +41,7 @@
         alphaEqn.relax();
         alphaEqn.solve();
 
-#       include "packingLimiter.H"
+        #include "packingLimiter.H"
 
         beta = scalar(1) - alpha;
 
diff --git a/applications/test/syncTools/Test-syncTools.C b/applications/test/syncTools/Test-syncTools.C
index 6f9d5f408a802e5e2c08b79b264952138073b0af..2638f0fb67ef3aeb939afdc118feb61b84a3f2f4 100644
--- a/applications/test/syncTools/Test-syncTools.C
+++ b/applications/test/syncTools/Test-syncTools.C
@@ -210,7 +210,7 @@ void testSparseData(const polyMesh& mesh, Random& rndGen)
     {
         // Create some data. Use slightly perturbed positions.
         Map<point> sparseData;
-        pointField fullData(mesh.nPoints(), point::max);
+        pointField fullData(mesh.nPoints(), point(GREAT, GREAT, GREAT));
 
         forAll(localPoints, i)
         {
@@ -236,7 +236,7 @@ void testSparseData(const polyMesh& mesh, Random& rndGen)
             mesh,
             fullData,
             minMagSqrEqOp<point>(),
-            point::max
+            point(GREAT, GREAT, GREAT)
             // true                    // apply separation
         );
 
@@ -246,7 +246,7 @@ void testSparseData(const polyMesh& mesh, Random& rndGen)
         {
             const point& fullPt = fullData[meshPointI];
 
-            if (fullPt != point::max)
+            if (fullPt != point(GREAT, GREAT, GREAT))
             {
                 const point& sparsePt = sparseData[meshPointI];
 
@@ -286,7 +286,7 @@ void testSparseData(const polyMesh& mesh, Random& rndGen)
     {
         // Create some data. Use slightly perturbed positions.
         EdgeMap<point> sparseData;
-        pointField fullData(mesh.nEdges(), point::max);
+        pointField fullData(mesh.nEdges(), point(GREAT, GREAT, GREAT));
 
         const edgeList& edges = allBoundary.edges();
         const labelList meshEdges = allBoundary.meshEdges
@@ -320,7 +320,7 @@ void testSparseData(const polyMesh& mesh, Random& rndGen)
             mesh,
             fullData,
             minMagSqrEqOp<point>(),
-            point::max
+            point(GREAT, GREAT, GREAT)
         );
 
         // Compare.
@@ -329,7 +329,7 @@ void testSparseData(const polyMesh& mesh, Random& rndGen)
         {
             const point& fullPt = fullData[meshEdgeI];
 
-            if (fullPt != point::max)
+            if (fullPt != point(GREAT, GREAT, GREAT))
             {
                 const point& sparsePt = sparseData[mesh.edges()[meshEdgeI]];
 
@@ -385,7 +385,7 @@ void testPointSync(const polyMesh& mesh, Random& rndGen)
             mesh,
             syncedPoints,
             minMagSqrEqOp<point>(),
-            point::max
+            point(GREAT, GREAT, GREAT)
         );
 
         forAll(syncedPoints, pointI)
@@ -461,7 +461,7 @@ void testEdgeSync(const polyMesh& mesh, Random& rndGen)
             mesh,
             syncedMids,
             minMagSqrEqOp<point>(),
-            point::max
+            point(GREAT, GREAT, GREAT)
         );
 
         forAll(syncedMids, edgeI)
diff --git a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C
index 8a7296eca00a1627cdc50abb9e53a27c9ad0ec64..5804d962c81949135978a2a61af2ec67b2c68221 100644
--- a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C
+++ b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C
@@ -43,6 +43,7 @@ Description
 #include "addPatchCellLayer.H"
 #include "fvMesh.H"
 #include "MeshedSurfaces.H"
+#include "globalIndex.H"
 
 #include "extrudedMesh.H"
 #include "extrudeModel.H"
@@ -261,6 +262,12 @@ int main(int argc, char *argv[])
         sourceCasePath.expand();
         fileName sourceRootDir = sourceCasePath.path();
         fileName sourceCaseDir = sourceCasePath.name();
+        if (Pstream::parRun())
+        {
+            sourceCaseDir =
+                sourceCaseDir
+               /"processor" + Foam::name(Pstream::myProcNo());
+        }
         wordList sourcePatches;
         dict.lookup("sourcePatches") >> sourcePatches;
 
@@ -279,54 +286,173 @@ int main(int argc, char *argv[])
             sourceRootDir,
             sourceCaseDir
         );
+
         #include "createMesh.H"
 
         const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
 
-        // Topo change container. Either copy an existing mesh or start
-        // with empty storage (number of patches only needed for checking)
-        autoPtr<polyTopoChange> meshMod
-        (
-            (
-                mode == MESH
-              ? new polyTopoChange(mesh)
-              : new polyTopoChange(patches.size())
-            )
-        );
-
         // Extrusion engine. Either adding to existing mesh or
         // creating separate mesh.
         addPatchCellLayer layerExtrude(mesh, (mode == MESH));
 
+        const labelList meshFaces(patchFaces(patches, sourcePatches));
         indirectPrimitivePatch extrudePatch
         (
             IndirectList<face>
             (
                 mesh.faces(),
-                patchFaces(patches, sourcePatches)
+                meshFaces
             ),
             mesh.points()
         );
 
+        // Determine extrudePatch normal
+        pointField extrudePatchPointNormals
+        (
+            PatchTools::pointNormals    //calcNormals
+            (
+                mesh,
+                extrudePatch,
+                meshFaces
+            )
+        );
+
+
+        // Precalculate mesh edges for pp.edges.
+        const labelList meshEdges
+        (
+            extrudePatch.meshEdges
+            (
+                mesh.edges(),
+                mesh.pointEdges()
+            )
+        );
+
+        // Global face indices engine
+        const globalIndex globalFaces(mesh.nFaces());
+
+        // Determine extrudePatch.edgeFaces in global numbering (so across
+        // coupled patches)
+        labelListList edgeGlobalFaces
+        (
+            addPatchCellLayer::globalEdgeFaces
+            (
+                mesh,
+                globalFaces,
+                extrudePatch
+            )
+        );
+
+
+        // Determine what patches boundary edges need to get extruded into.
+        // This might actually cause edge-connected processors to become
+        // face-connected so might need to introduce new processor boundaries.
+        // Calculates:
+        //  - per pp.edge the patch to extrude into
+        //  - any additional processor boundaries (patchToNbrProc = map from
+        //    new patchID to neighbour processor)
+        //  - number of new patches (nPatches)
+
+        labelList sidePatchID;
+        label nPatches;
+        Map<label> nbrProcToPatch;
+        Map<label> patchToNbrProc;
+        addPatchCellLayer::calcSidePatch
+        (
+            mesh,
+            globalFaces,
+            edgeGlobalFaces,
+            extrudePatch,
+
+            sidePatchID,
+            nPatches,
+            nbrProcToPatch,
+            patchToNbrProc
+        );
+
+
+        // Add any patches.
+
+        label nAdded = nPatches - mesh.boundaryMesh().size();
+        reduce(nAdded, sumOp<label>());
+
+        Info<< "Adding overall " << nAdded << " processor patches." << endl;
+
+        if (nAdded > 0)
+        {
+            DynamicList<polyPatch*> newPatches(nPatches);
+            forAll(mesh.boundaryMesh(), patchI)
+            {
+                newPatches.append
+                (
+                    mesh.boundaryMesh()[patchI].clone
+                    (
+                        mesh.boundaryMesh()
+                    ).ptr()
+                );
+            }
+            for
+            (
+                label patchI = mesh.boundaryMesh().size();
+                patchI < nPatches;
+                patchI++
+            )
+            {
+                label nbrProcI = patchToNbrProc[patchI];
+                word name =
+                        "procBoundary"
+                      + Foam::name(Pstream::myProcNo())
+                      + "to"
+                      + Foam::name(nbrProcI);
+
+                Pout<< "Adding patch " << patchI
+                    << " name:" << name
+                    << " between " << Pstream::myProcNo()
+                    << " and " << nbrProcI
+                    << endl;
+
+
+                newPatches.append
+                (
+                    new processorPolyPatch
+                    (
+                        name,
+                        0,                  // size
+                        mesh.nFaces(),      // start
+                        patchI,             // index
+                        mesh.boundaryMesh(),// polyBoundaryMesh
+                        Pstream::myProcNo(),// myProcNo
+                        nbrProcI            // neighbProcNo
+                    )
+                );
+            }
+
+            // Add patches. Do no parallel updates.
+            mesh.removeFvBoundary();
+            mesh.addFvPatches(newPatches, true);
+        }
+
+
+
         // Only used for addPatchCellLayer into new mesh
-        labelList exposedPatchIDs;
+        labelList exposedPatchID;
         if (mode == PATCH)
         {
             dict.lookup("exposedPatchName") >> backPatchName;
-            exposedPatchIDs.setSize
+            exposedPatchID.setSize
             (
                 extrudePatch.size(),
                 findPatchID(patches, backPatchName)
             );
         }
 
-
+        // Determine points and extrusion
         pointField layer0Points(extrudePatch.nPoints());
         pointField displacement(extrudePatch.nPoints());
         forAll(displacement, pointI)
         {
-            const vector& patchNormal = extrudePatch.pointNormals()[pointI];
+            const vector& patchNormal = extrudePatchPointNormals[pointI];
 
             // layer0 point
             layer0Points[pointI] = model()
@@ -373,15 +499,31 @@ int main(int argc, char *argv[])
         // Layers per point
         labelList nPointLayers(extrudePatch.nPoints(), model().nLayers());
         // Displacement for first layer
-        vectorField firstLayerDisp(displacement*model().sumThickness(1));
+        vectorField firstLayerDisp = displacement*model().sumThickness(1);
+
         // Expansion ratio not used.
         scalarField ratio(extrudePatch.nPoints(), 1.0);
 
+        // Topo change container. Either copy an existing mesh or start
+        // with empty storage (number of patches only needed for checking)
+        autoPtr<polyTopoChange> meshMod
+        (
+            (
+                mode == MESH
+              ? new polyTopoChange(mesh)
+              : new polyTopoChange(patches.size())
+            )
+        );
+
         layerExtrude.setRefinement
         (
+            globalFaces,
+            edgeGlobalFaces,
+
             ratio,              // expansion ratio
             extrudePatch,       // patch faces to extrude
-            exposedPatchIDs,    // if new mesh: patches for exposed faces
+            sidePatchID,        // if boundary edge: patch to use
+            exposedPatchID,     // if new mesh: patches for exposed faces
             nFaceLayers,
             nPointLayers,
             firstLayerDisp,
@@ -401,7 +543,7 @@ int main(int argc, char *argv[])
                     model()
                     (
                         extrudePatch.localPoints()[pointI],
-                        extrudePatch.pointNormals()[pointI],
+                        extrudePatchPointNormals[pointI],
                         pPointI+1       // layer
                     )
                 );
diff --git a/applications/utilities/miscellaneous/foamDebugSwitches/Make/options b/applications/utilities/miscellaneous/foamDebugSwitches/Make/options
index 826bc4c144f01bf27c3e9f1ac301e3c612c58611..11303fd0e875b1c46ef306abeb0af14782ccb3fb 100644
--- a/applications/utilities/miscellaneous/foamDebugSwitches/Make/options
+++ b/applications/utilities/miscellaneous/foamDebugSwitches/Make/options
@@ -1,25 +1,32 @@
 EXE_LIBS = \
     -lautoMesh \
     -lbarotropicCompressibilityModel \
+    -lbasicSolidThermo \
     -lbasicThermophysicalModels \
     -lblockMesh \
     -lchemistryModel \
     -lcoalCombustion \
+    -lcombustionModels \
     -lcompressibleLESModels \
     -lcompressibleRASModels \
     -lcompressibleTurbulenceModel \
     -lconversion \
     -ldecompositionMethods \
     -ldieselSpray \
+    -ldistributed \
+    -ldistributionModels \
     -ldsmc \
     -ldynamicFvMesh \
     -ldynamicMesh \
     -ledgeMesh \
     -lengine \
     -lerrorEstimation \
+    -lEulerianInterfacialModels \
+    -lextrudeModel \
     -lfieldFunctionObjects \
     -lfileFormats \
     -lfiniteVolume \
+    -lfoamCalcFunctions \
     -lforces \
     -lfvMotionSolvers \
     -lgenericPatchFields \
@@ -29,6 +36,8 @@ EXE_LIBS = \
     -lincompressibleTurbulenceModel \
     -linterfaceProperties \
     -lIOFunctionObjects \
+    -ljobControl \
+    -lkineticTheoryModel \
     -llagrangian \
     -llagrangianIntermediate \
     -llaminarFlameSpeedModels \
@@ -37,22 +46,29 @@ EXE_LIBS = \
     -lliquidMixtureProperties \
     -lliquidProperties \
     -lmeshTools \
+    -lMGridGenGAMGAgglomeration \
     -lmolecularMeasurements \
     -lmolecule \
+    -lmultiphaseInterFoam \
     -lODE \
     -lOpenFOAM \
-    -ldistributionModels \
+    -lphaseModel \
     -lpotential \
     -lradiationModels \
     -lrandomProcesses \
     -lreactionThermophysicalModels \
+    -lreconstruct \
     -lsampling \
+    -lSLGThermo \
     -lsolidMixtureProperties \
     -lsolidParticle \
     -lsolidProperties \
+    -lsolid \
     -lspecie \
+    -lsurfaceFilmModels \
     -lsurfMesh \
     -lsystemCall \
+    -ltabulatedWallFunctions \
     -lthermalPorousZone \
     -lthermophysicalFunctions \
     -ltopoChangerFvMesh \
diff --git a/applications/utilities/parallelProcessing/redistributeMeshPar/redistributeMeshPar.C b/applications/utilities/parallelProcessing/redistributeMeshPar/redistributeMeshPar.C
index f429d63f571b97ff373e48d25d8da5fede364357..11922c3c380dc55c9df1c77d9ad102929f8a7c70 100644
--- a/applications/utilities/parallelProcessing/redistributeMeshPar/redistributeMeshPar.C
+++ b/applications/utilities/parallelProcessing/redistributeMeshPar/redistributeMeshPar.C
@@ -73,8 +73,8 @@ autoPtr<fvMesh> createMesh
     const bool haveMesh
 )
 {
-    Pout<< "Create mesh for time = "
-        << runTime.timeName() << nl << endl;
+    //Pout<< "Create mesh for time = "
+    //    << runTime.timeName() << nl << endl;
 
     IOobject io
     (
@@ -98,12 +98,12 @@ autoPtr<fvMesh> createMesh
             xferCopy(labelList()),
             false
         );
-        Pout<< "Writing dummy mesh to " << dummyMesh.polyMesh::objectPath()
-            << endl;
+        //Pout<< "Writing dummy mesh to " << dummyMesh.polyMesh::objectPath()
+        //    << endl;
         dummyMesh.write();
     }
 
-    Pout<< "Reading mesh from " << io.objectPath() << endl;
+    //Pout<< "Reading mesh from " << io.objectPath() << endl;
     autoPtr<fvMesh> meshPtr(new fvMesh(io));
     fvMesh& mesh = meshPtr();
 
@@ -201,9 +201,8 @@ autoPtr<fvMesh> createMesh
                     break;
                 }
 
-                Pout<< "Adding patch:" << nPatches
-                    << " name:" << name
-                    << " type:" << type << endl;
+                //Pout<< "Adding patch:" << nPatches
+                //    << " name:" << name << " type:" << type << endl;
 
                 dictionary patchDict(e.dict());
                 patchDict.remove("nFaces");
@@ -282,8 +281,7 @@ autoPtr<fvMesh> createMesh
     if (!haveMesh)
     {
         // We created a dummy mesh file above. Delete it.
-        Pout<< "Removing dummy mesh " << io.objectPath()
-            << endl;
+        //Pout<< "Removing dummy mesh " << io.objectPath() << endl;
         rmDir(io.objectPath());
     }
 
@@ -351,17 +349,19 @@ scalar getMergeDistance
 }
 
 
-void printMeshData(Ostream& os, const polyMesh& mesh)
-{
-    os  << "Number of points:           " << mesh.points().size() << nl
-        << "          faces:            " << mesh.faces().size() << nl
-        << "          internal faces:   " << mesh.faceNeighbour().size() << nl
-        << "          cells:            " << mesh.cells().size() << nl
-        << "          boundary patches: " << mesh.boundaryMesh().size() << nl
-        << "          point zones:      " << mesh.pointZones().size() << nl
-        << "          face zones:       " << mesh.faceZones().size() << nl
-        << "          cell zones:       " << mesh.cellZones().size() << nl;
-}
+//void printMeshData(Ostream& os, const polyMesh& mesh)
+//{
+//    os  << "Number of points:           " << mesh.points().size() << nl
+//        << "          faces:            " << mesh.faces().size() << nl
+//        << "          internal faces:   " << mesh.faceNeighbour().size() << nl
+//        << "          cells:            " << mesh.cells().size() << nl
+//        << "          boundary patches: " << mesh.boundaryMesh().size() << nl
+//        << "          point zones:      " << mesh.pointZones().size() << nl
+//        << "          face zones:       " << mesh.faceZones().size() << nl
+//        << "          cell zones:       " << mesh.cellZones().size() << nl;
+//}
+
+
 void printMeshData(const polyMesh& mesh)
 {
     // Collect all data on master
@@ -525,6 +525,7 @@ void readFields
 
             // Receive field
             IPstream fromMaster(Pstream::blocking, Pstream::masterNo());
+            dictionary fieldDict(fromMaster);
 
             fields.set
             (
@@ -540,7 +541,7 @@ void readFields
                         IOobject::AUTO_WRITE
                     ),
                     mesh,
-                    fromMaster
+                    fieldDict
                 )
             );
 
@@ -916,8 +917,8 @@ int main(int argc, char *argv[])
     // Mesh distribution engine
     fvMeshDistribute distributor(mesh, tolDim);
 
-    Pout<< "Wanted distribution:"
-        << distributor.countCells(finalDecomp) << nl << endl;
+    //Pout<< "Wanted distribution:"
+    //    << distributor.countCells(finalDecomp) << nl << endl;
 
     // Do actual sending/receiving of mesh
     autoPtr<mapDistributePolyMesh> map = distributor.distribute(finalDecomp);
@@ -987,7 +988,7 @@ int main(int argc, char *argv[])
     Info<< endl;
 
 
-    Pout<< "End\n" << endl;
+    Info<< "End\n" << endl;
 
     return 0;
 }
diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
index b0d3d83ec02638c43d00d11a491d3a1a8462b73d..6d83d1836324e3ee2516c9059dadd826f83c52be 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,6 +36,7 @@ Description
 #include "Time.H"
 #include "surfaceFeatures.H"
 #include "featureEdgeMesh.H"
+#include "extendedFeatureEdgeMesh.H"
 #include "treeBoundBox.H"
 #include "meshTools.H"
 #include "OFstream.H"
@@ -298,20 +299,46 @@ int main(int argc, char *argv[])
 
     // Extracting and writing a featureEdgeMesh
 
-    Pout<< nl << "Writing featureEdgeMesh to constant/featureEdgeMesh."
-        << endl;
-
-    featureEdgeMesh feMesh
+    extendedFeatureEdgeMesh feMesh
     (
         newSet,
         runTime,
         surfFileName.lessExt().name() + ".featureEdgeMesh"
     );
 
+    Info<< nl << "Writing extendedFeatureEdgeMesh to " << feMesh.objectPath()
+        << endl;
+
+
     feMesh.writeObj(surfFileName.lessExt().name());
 
     feMesh.write();
 
+
+    // Write a featureEdgeMesh for backwards compatibility
+    {
+        featureEdgeMesh bfeMesh
+        (
+            IOobject
+            (
+                surfFileName.lessExt().name() + ".eMesh",   // name
+                runTime.constant(), // instance
+                "triSurface",
+                runTime,            // registry
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            feMesh.points(),
+            feMesh.edges()
+        );
+
+        Info<< nl << "Writing featureEdgeMesh to "
+            << bfeMesh.objectPath() << endl;
+
+        bfeMesh.regIOobject::write();
+    }
+
     Info<< "End\n" << endl;
 
     return 0;
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C
index 0a906aa7bcfd22bdd707f8bca9afbf5c16f3b452..f8f7ab41b0052a056d4bf2956dc0a7ab93452e59 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C
@@ -43,75 +43,6 @@ defineTypeNameAndDebug(Foam::addPatchCellLayer, 0);
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-// Calculate global faces per pp edge.
-Foam::labelListList Foam::addPatchCellLayer::calcGlobalEdgeFaces
-(
-    const polyMesh& mesh,
-    const globalIndex& globalFaces,
-    const indirectPrimitivePatch& pp,
-    const labelList& meshEdges
-)
-{
-    //// Determine coupled edges just so we don't have to have storage
-    //// for all non-coupled edges.
-    //
-    //PackedBoolList isCoupledEdge(mesh.nEdges());
-    //
-    //const polyBoundaryMesh& patches = mesh.boundaryMesh();
-    //
-    //forAll(patches, patchI)
-    //{
-    //    const polyPatch& pp = patches[patchI];
-    //
-    //    if (pp.coupled())
-    //    {
-    //        const labelList& meshEdges = pp.meshEdges();
-    //
-    //        forAll(meshEdges, i)
-    //        {
-    //            isCoupledEdge.set(meshEdges[i], 1);
-    //        }
-    //    }
-    //}
-
-    // From mesh edge to global face labels. Sized only for pp edges.
-    labelListList globalEdgeFaces(mesh.nEdges());
-
-    const labelListList& edgeFaces = pp.edgeFaces();
-
-    forAll(edgeFaces, edgeI)
-    {
-        label meshEdgeI = meshEdges[edgeI];
-
-        //if (isCoupledEdge.get(meshEdgeI) == 1)
-        {
-            const labelList& eFaces = edgeFaces[edgeI];
-
-            // Store face and processor as unique tag.
-            labelList& globalEFaces = globalEdgeFaces[meshEdgeI];
-            globalEFaces.setSize(eFaces.size());
-            forAll(eFaces, i)
-            {
-                globalEFaces[i] =
-                    globalFaces.toGlobal(pp.addressing()[eFaces[i]]);
-            }
-        }
-    }
-
-    // Synchronise across coupled edges.
-    syncTools::syncEdgeList
-    (
-        mesh,
-        globalEdgeFaces,
-        uniqueEqOp(),
-        labelList()     // null value
-    );
-
-    // Extract pp part
-    return labelListList(UIndirectList<labelList>(globalEdgeFaces, meshEdges));
-}
-
-
 Foam::label Foam::addPatchCellLayer::nbrFace
 (
     const labelListList& edgeFaces,
@@ -316,12 +247,12 @@ Foam::labelPair Foam::addPatchCellLayer::getEdgeString
 Foam::label Foam::addPatchCellLayer::addSideFace
 (
     const indirectPrimitivePatch& pp,
-    const labelList& patchID,           // prestored patch per pp face
     const labelListList& addedCells,    // per pp face the new extruded cell
     const face& newFace,
+    const label newPatchID,
+
     const label ownFaceI,               // pp face that provides owner
     const label nbrFaceI,
-    const label patchEdgeI,             // edge to add to
     const label meshEdgeI,              // corresponding mesh edge
     const label layerI,                 // layer
     const label numEdgeFaces,           // number of layers for edge
@@ -329,8 +260,9 @@ Foam::label Foam::addPatchCellLayer::addSideFace
     polyTopoChange& meshMod
 ) const
 {
-    // Edge to 'inflate' from
+    // Face or edge to 'inflate' from
     label inflateEdgeI = -1;
+    label inflateFaceI = -1;
 
     // Check mesh faces using edge
     if (addToMesh_)
@@ -346,8 +278,6 @@ Foam::label Foam::addPatchCellLayer::addSideFace
         }
     }
 
-    // Get my mesh face and its zone.
-    label meshFaceI = pp.addressing()[ownFaceI];
     // Zone info comes from any side patch face. Otherwise -1 since we
     // don't know what to put it in - inherit from the extruded faces?
     label zoneI = -1;   //mesh_.faceZones().whichZone(meshFaceI);
@@ -358,14 +288,15 @@ Foam::label Foam::addPatchCellLayer::addSideFace
     // Is patch edge external edge of indirectPrimitivePatch?
     if (nbrFaceI == -1)
     {
-        // External edge so external face. Patch id is obtained from
-        // any other patch connected to edge.
+        // External edge so external face.
 
         const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
         // Loop over all faces connected to edge to inflate and
-        // see if any boundary face (but not meshFaceI)
-        label otherPatchID = patchID[ownFaceI];
+        // see if we can find a face that is otherPatchID
+
+        // Get my mesh face and its zone.
+        label meshFaceI = pp.addressing()[ownFaceI];
 
         forAll(meshFaces, k)
         {
@@ -373,11 +304,14 @@ Foam::label Foam::addPatchCellLayer::addSideFace
 
             if
             (
-                faceI != meshFaceI
-             && !mesh_.isInternalFace(faceI)
+                (faceI != meshFaceI)
+             && (patches.whichPatch(faceI) == newPatchID)
             )
             {
-                otherPatchID = patches.whichPatch(faceI);
+                // Found the patch face. Use it to inflate from
+                inflateEdgeI = -1;
+                inflateFaceI = faceI;
+
                 zoneI = mesh_.faceZones().whichZone(faceI);
                 if (zoneI != -1)
                 {
@@ -414,7 +348,7 @@ Foam::label Foam::addPatchCellLayer::addSideFace
 
         //Pout<< "Added boundary face:" << newFace
         //    << " own:" << addedCells[ownFaceI][layerOwn]
-        //    << " patch:" << otherPatchID
+        //    << " patch:" << newPatchID
         //    << endl;
 
         addedFaceI = meshMod.setAction
@@ -426,9 +360,9 @@ Foam::label Foam::addPatchCellLayer::addSideFace
                 -1,                         // neighbour
                 -1,                         // master point
                 inflateEdgeI,               // master edge
-                -1,                         // master face
+                inflateFaceI,               // master face
                 false,                      // flux flip
-                otherPatchID,               // patch for face
+                newPatchID,                 // patch for face
                 zoneI,                      // zone for face
                 flip                        // face zone flip
             )
@@ -510,6 +444,51 @@ Foam::label Foam::addPatchCellLayer::addSideFace
 }
 
 
+Foam::label Foam::addPatchCellLayer::findProcPatch
+(
+    const polyMesh& mesh,
+    const label nbrProcID
+)
+{
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    forAll(mesh.globalData().processorPatches(), i)
+    {
+        label patchI = mesh.globalData().processorPatches()[i];
+
+        if
+        (
+            refCast<const processorPolyPatch>(patches[patchI]).neighbProcNo()
+         == nbrProcID
+        )
+        {
+            return patchI;
+        }
+    }
+    return -1;
+}
+
+
+void Foam::addPatchCellLayer::setFaceProps
+(
+    const polyMesh& mesh,
+    const label faceI,
+
+    label& patchI,
+    label& zoneI,
+    bool& zoneFlip
+)
+{
+    patchI = mesh.boundaryMesh().whichPatch(faceI);
+    zoneI = mesh.faceZones().whichZone(faceI);
+    if (zoneI != -1)
+    {
+        label index = mesh.faceZones()[zoneI].whichFace(faceI);
+        zoneFlip = mesh.faceZones()[zoneI].flipMap()[index];
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from mesh
@@ -561,10 +540,251 @@ Foam::labelListList Foam::addPatchCellLayer::addedCells() const
 }
 
 
+// Calculate global faces per pp edge.
+Foam::labelListList Foam::addPatchCellLayer::globalEdgeFaces
+(
+    const polyMesh& mesh,
+    const globalIndex& globalFaces,
+    const indirectPrimitivePatch& pp
+)
+{
+    // Precalculate mesh edges for pp.edges.
+    const labelList meshEdges(pp.meshEdges(mesh.edges(), mesh.pointEdges()));
+
+    // From mesh edge to global face labels. Non-empty sublists only for
+    // pp edges.
+    labelListList globalEdgeFaces(mesh.nEdges());
+
+    const labelListList& edgeFaces = pp.edgeFaces();
+
+    forAll(edgeFaces, edgeI)
+    {
+        label meshEdgeI = meshEdges[edgeI];
+
+        const labelList& eFaces = edgeFaces[edgeI];
+
+        // Store face and processor as unique tag.
+        labelList& globalEFaces = globalEdgeFaces[meshEdgeI];
+        globalEFaces.setSize(eFaces.size());
+        forAll(eFaces, i)
+        {
+            globalEFaces[i] = globalFaces.toGlobal(pp.addressing()[eFaces[i]]);
+        }
+    }
+
+    // Synchronise across coupled edges.
+    syncTools::syncEdgeList
+    (
+        mesh,
+        globalEdgeFaces,
+        uniqueEqOp(),
+        labelList()             // null value
+    );
+
+    // Extract pp part
+    return labelListList(UIndirectList<labelList>(globalEdgeFaces, meshEdges));
+}
+
+
+void Foam::addPatchCellLayer::calcSidePatch
+(
+    const polyMesh& mesh,
+    const globalIndex& globalFaces,
+    const labelListList& globalEdgeFaces,
+    const indirectPrimitivePatch& pp,
+
+    labelList& sidePatchID,
+    label& nPatches,
+    Map<label>& nbrProcToPatch,
+    Map<label>& patchToNbrProc
+)
+{
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    // Precalculate mesh edges for pp.edges.
+    const labelList meshEdges(pp.meshEdges(mesh.edges(), mesh.pointEdges()));
+
+    sidePatchID.setSize(pp.nEdges());
+    sidePatchID = -1;
+
+    // These also get determined but not (yet) exported:
+    // - whether face is created from other face or edge
+    // - what zone&orientation face should have
+
+    labelList inflateEdgeI(pp.nEdges(), -1);
+    labelList inflateFaceI(pp.nEdges(), -1);
+    labelList sideZoneID(pp.nEdges(), -1);
+    boolList sideFlip(pp.nEdges(), false);
+
+    nPatches = patches.size();
+
+    forAll(globalEdgeFaces, edgeI)
+    {
+        const labelList& eGlobalFaces = globalEdgeFaces[edgeI];
+        if
+        (
+            eGlobalFaces.size() == 2
+         && pp.edgeFaces()[edgeI].size() == 1
+        )
+        {
+            // Locally but not globally a boundary edge. Hence a coupled
+            // edge. Find the patch to use if on different
+            // processors.
+
+            label f0 = eGlobalFaces[0];
+            label f1 = eGlobalFaces[1];
+
+            label otherProcI = -1;
+            if (globalFaces.isLocal(f0) && !globalFaces.isLocal(f1))
+            {
+                otherProcI = globalFaces.whichProcID(f1);
+            }
+            else if (!globalFaces.isLocal(f0) && globalFaces.isLocal(f1))
+            {
+                otherProcI = globalFaces.whichProcID(f0);
+            }
+
+
+            if (otherProcI != -1)
+            {
+                sidePatchID[edgeI] = findProcPatch(mesh, otherProcI);
+                if (sidePatchID[edgeI] == -1)
+                {
+                    // Cannot find a patch to processor. See if already
+                    // marked for addition
+                    if (nbrProcToPatch.found(otherProcI))
+                    {
+                        sidePatchID[edgeI] = nbrProcToPatch[otherProcI];
+                    }
+                    else
+                    {
+                        sidePatchID[edgeI] = nPatches;
+                        nbrProcToPatch.insert(otherProcI, nPatches);
+                        patchToNbrProc.insert(nPatches, otherProcI);
+                        nPatches++;
+                    }
+                }
+            }
+        }
+    }
+
+
+
+    // Determine face properties for all other boundary edges
+    // ------------------------------------------------------
+
+    const labelListList& edgeFaces = pp.edgeFaces();
+    forAll(edgeFaces, edgeI)
+    {
+        if (edgeFaces[edgeI].size() == 1 && sidePatchID[edgeI] == -1)
+        {
+            // Proper, uncoupled patch edge.
+
+            label myFaceI = pp.addressing()[edgeFaces[edgeI][0]];
+
+            // Pick up any boundary face on this edge and use its properties
+            label meshEdgeI = meshEdges[edgeI];
+            const labelList& meshFaces = mesh.edgeFaces()[meshEdgeI];
+
+            forAll(meshFaces, k)
+            {
+                label faceI = meshFaces[k];
+
+                if (faceI != myFaceI && !mesh.isInternalFace(faceI))
+                {
+                    setFaceProps
+                    (
+                        mesh,
+                        faceI,
+
+                        sidePatchID[edgeI],
+                        sideZoneID[edgeI],
+                        sideFlip[edgeI]
+                    );
+                    inflateFaceI[edgeI] = faceI;
+                    inflateEdgeI[edgeI] = -1;
+
+                    break;
+                }
+            }
+        }
+    }
+
+
+
+    // Now hopefully every boundary edge has a side patch. Check
+    forAll(edgeFaces, edgeI)
+    {
+        if (edgeFaces[edgeI].size() == 1 && sidePatchID[edgeI] == -1)
+        {
+            const edge& e = pp.edges()[edgeI];
+            FatalErrorIn("addPatchCellLayer::calcSidePatch(..)")
+                << "Have no sidePatchID for edge " << edgeI << " points "
+                << pp.points()[pp.meshPoints()[e[0]]]
+                << pp.points()[pp.meshPoints()[e[1]]]
+                << abort(FatalError);
+        }
+    }
+
+
+
+    // Now we have sidepatch see if we have patchface or edge to inflate
+    // from.
+    forAll(edgeFaces, edgeI)
+    {
+        if (edgeFaces[edgeI].size() == 1 && inflateFaceI[edgeI] == -1)
+        {
+            // 1. Do we have a boundary face to inflate from
+
+            label myFaceI = pp.addressing()[edgeFaces[edgeI][0]];
+
+            // Pick up any boundary face on this edge and use its properties
+            label meshEdgeI = meshEdges[edgeI];
+            const labelList& meshFaces = mesh.edgeFaces()[meshEdgeI];
+
+            forAll(meshFaces, k)
+            {
+                label faceI = meshFaces[k];
+
+                if (faceI != myFaceI)
+                {
+                    if (mesh.isInternalFace(faceI))
+                    {
+                        inflateEdgeI[edgeI] = meshEdgeI;
+                    }
+                    else
+                    {
+                        if (patches.whichPatch(faceI) == sidePatchID[edgeI])
+                        {
+                            setFaceProps
+                            (
+                                mesh,
+                                faceI,
+
+                                sidePatchID[edgeI],
+                                sideZoneID[edgeI],
+                                sideFlip[edgeI]
+                            );
+                            inflateFaceI[edgeI] = faceI;
+                            inflateEdgeI[edgeI] = -1;
+
+                            break;
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
+
+
 void Foam::addPatchCellLayer::setRefinement
 (
+    const globalIndex& globalFaces,
+    const labelListList& globalEdgeFaces,
     const scalarField& expansionRatio,
     const indirectPrimitivePatch& pp,
+    const labelList& sidePatchID,
     const labelList& exposedPatchID,
     const labelList& nFaceLayers,
     const labelList& nPointLayers,
@@ -575,7 +795,7 @@ void Foam::addPatchCellLayer::setRefinement
     if (debug)
     {
         Pout<< "addPatchCellLayer::setRefinement : Adding up to "
-            << max(nPointLayers)
+            << gMax(nPointLayers)
             << " layers of cells to indirectPrimitivePatch with "
             << pp.nPoints() << " points" << endl;
     }
@@ -788,8 +1008,6 @@ void Foam::addPatchCellLayer::setRefinement
 
                 label meshEdgeI = meshEdges[edgeI];
 
-                // Mesh faces using edge
-
                 // Mesh faces using edge
                 const labelList& meshFaces = mesh_.edgeFaces(meshEdgeI, ef);
 
@@ -1213,22 +1431,6 @@ void Foam::addPatchCellLayer::setRefinement
     }
 
 
-    // Global indices engine
-    const globalIndex globalFaces(mesh_.nFaces());
-
-    // Get for all pp edgeFaces a unique faceID
-    labelListList globalEdgeFaces
-    (
-         calcGlobalEdgeFaces
-         (
-            mesh_,
-            globalFaces,
-            pp,
-            meshEdges
-        )
-    );
-
-
     // Mark off which edges have been extruded
     boolList doneEdge(pp.nEdges(), false);
 
@@ -1474,16 +1676,17 @@ void Foam::addPatchCellLayer::setRefinement
                         addSideFace
                         (
                             pp,
-                            patchID,
                             addedCells,
-                            newFace,
+
+                            newFace,                // vertices of new face
+                            sidePatchID[startEdgeI],// -1 or patch for face
+
                             patchFaceI,
                             nbrFaceI,
-                            startEdgeI,     // edge to inflate from
-                            meshEdgeI,      // corresponding mesh edge
-                            i,
-                            numEdgeSideFaces,
-                            meshFaces,
+                            meshEdgeI,          // (mesh) edge to inflate
+                            i,                  // layer
+                            numEdgeSideFaces,   // num layers
+                            meshFaces,          // edgeFaces
                             meshMod
                         );
                     }
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H b/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H
index 4835ac68e39ef79ba61246ea7893e81f10ee3b24..2458de9fa45324aabf7252fa6fd9e3f17396263f 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H
@@ -181,16 +181,6 @@ class addPatchCellLayer
 
     // Private Member Functions
 
-        //- Per patch edge the pp faces (in global indices) using it. Uses
-        //  uniqueEqOp() to remove duplicates.
-        labelListList calcGlobalEdgeFaces
-        (
-            const polyMesh& mesh,
-            const globalIndex& globalFaces,
-            const indirectPrimitivePatch& pp,
-            const labelList& meshEdges
-        );
-
         //- Get the face on the other side of the edge.
         static label nbrFace
         (
@@ -226,12 +216,13 @@ class addPatchCellLayer
         label addSideFace
         (
             const indirectPrimitivePatch&,
-            const labelList& patchID,
             const labelListList& addedCells,
+
             const face& newFace,
+            const label newPatchID,
+
             const label ownFaceI,
             const label nbrFaceI,
-            const label patchEdgeI,
             const label meshEdgeI,
             const label layerI,
             const label numEdgeFaces,
@@ -239,6 +230,18 @@ class addPatchCellLayer
             polyTopoChange&
         ) const;
 
+        //- Find patch to neighbouring processor
+        static label findProcPatch(const polyMesh&, const label nbrProcID);
+
+        //- Extract properties from mesh face
+        static void setFaceProps
+        (
+            const polyMesh&,
+            const label,
+            label&,
+            label&,
+            bool&
+        );
 
         //- Disallow default bitwise copy construct
         addPatchCellLayer(const addPatchCellLayer&);
@@ -256,7 +259,7 @@ public:
     // Constructors
 
         //- Construct from mesh.
-        addPatchCellLayer(const polyMesh& mesh, const bool addToMesh = true);
+        addPatchCellLayer(const polyMesh&, const bool addToMesh = true);
 
 
     // Member Functions
@@ -291,6 +294,33 @@ public:
 
         // Edit
 
+            //- Per patch edge the pp faces (in global indices) using it. Uses
+            //  uniqueEqOp() to remove duplicates.
+            static labelListList globalEdgeFaces
+            (
+                const polyMesh&,
+                const globalIndex& globalFaces,
+                const indirectPrimitivePatch& pp
+            );
+
+            //- Boundary edges get extruded into boundary faces. Determine patch
+            //  for these faces. This might be a to-be-created processor patch
+            //  (patchI >= mesh.boundaryMesh().size()) in which case the
+            //  nbrProcToPatch, patchToNbrProc give the correspondence. nPatches
+            //  is the new number of patches.
+            static void calcSidePatch
+            (
+                const polyMesh&,
+                const globalIndex& globalFaces,
+                const labelListList& globalEdgeFaces,
+                const indirectPrimitivePatch& pp,
+
+                labelList& sidePatchID,
+                label& nPatches,
+                Map<label>& nbrProcToPatch,
+                Map<label>& patchToNbrProc
+            );
+
             //- Play commands into polyTopoChange to create layers on top
             //  of indirectPrimitivePatch (have to be outside faces).
             //  Gets displacement per patch point.
@@ -313,8 +343,11 @@ public:
             //        (instead of e.g. from patch faces)
             void setRefinement
             (
+                const globalIndex& globalFaces,
+                const labelListList& globalEdgeFaces,
                 const scalarField& expansionRatio,
                 const indirectPrimitivePatch& pp,
+                const labelList& sidePatchID,
                 const labelList& exposedPatchID,
                 const labelList& nFaceLayers,
                 const labelList& nPointLayers,
@@ -326,20 +359,26 @@ public:
             //- Add with constant expansion ratio and same nLayers everywhere
             void setRefinement
             (
+                const globalIndex& globalFaces,
+                const labelListList& globalEdgeFaces,
                 const label nLayers,
                 const indirectPrimitivePatch& pp,
+                const labelList& sidePatchID,
                 const vectorField& overallDisplacement,
                 polyTopoChange& meshMod
             )
             {
                 setRefinement
                 (
-                    scalarField(pp.nPoints(), 1.0), // expansion ration
+                    globalFaces,
+                    globalEdgeFaces,
+                    scalarField(pp.nPoints(), 1.0),     // expansion ration
                     pp,
+                    sidePatchID,
                     labelList(0),
-                    labelList(pp.size(), nLayers),
-                    labelList(pp.nPoints(), nLayers),
-                    overallDisplacement / nLayers,
+                    labelList(pp.size(), nLayers),      // nFaceLayers
+                    labelList(pp.nPoints(), nLayers),   // nPointLayers
+                    overallDisplacement / nLayers,      // firstLayerDisp
                     meshMod
                 );
             }
diff --git a/src/edgeMesh/Make/files b/src/edgeMesh/Make/files
index cdad2399991b81b82fda8df34d4d62a796a1363b..52f985de8cd7745f26ff798972dd1ba0fcc7ce5b 100644
--- a/src/edgeMesh/Make/files
+++ b/src/edgeMesh/Make/files
@@ -20,7 +20,10 @@ $(edgeFormats)/starcd/STARCDedgeFormatRunTime.C
 $(edgeFormats)/vtk/VTKedgeFormat.C
 $(edgeFormats)/vtk/VTKedgeFormatRunTime.C
 
+extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C
+
 featureEdgeMesh/featureEdgeMesh.C
 
 
+
 LIB = $(FOAM_LIBBIN)/libedgeMesh
diff --git a/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C
new file mode 100644
index 0000000000000000000000000000000000000000..c4f536e2092d725db61f84b03a368b1c5ac68a1c
--- /dev/null
+++ b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C
@@ -0,0 +1,1028 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2004-2011 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 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 "extendedFeatureEdgeMesh.H"
+#include "triSurface.H"
+#include "Random.H"
+#include "Time.H"
+#include "meshTools.H"
+#include "linePointRef.H"
+#include "ListListOps.H"
+#include "OFstream.H"
+#include "IFstream.H"
+#include "unitConversion.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(Foam::extendedFeatureEdgeMesh, 0);
+
+Foam::scalar Foam::extendedFeatureEdgeMesh::cosNormalAngleTol_ =
+    Foam::cos(degToRad(0.1));
+
+
+Foam::label Foam::extendedFeatureEdgeMesh::convexStart_ = 0;
+
+
+Foam::label Foam::extendedFeatureEdgeMesh::externalStart_ = 0;
+
+
+Foam::label Foam::extendedFeatureEdgeMesh::nPointTypes = 4;
+
+
+Foam::label Foam::extendedFeatureEdgeMesh::nEdgeTypes = 5;
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh(const IOobject& io)
+:
+    regIOobject(io),
+    edgeMesh(pointField(0), edgeList(0)),
+    concaveStart_(0),
+    mixedStart_(0),
+    nonFeatureStart_(0),
+    internalStart_(0),
+    flatStart_(0),
+    openStart_(0),
+    multipleStart_(0),
+    normals_(0),
+    edgeDirections_(0),
+    edgeNormals_(0),
+    featurePointNormals_(0),
+    regionEdges_(0),
+    edgeTree_(),
+    edgeTreesByType_()
+{
+    if
+    (
+        io.readOpt() == IOobject::MUST_READ
+     || io.readOpt() == IOobject::MUST_READ_IF_MODIFIED
+     || (io.readOpt() == IOobject::READ_IF_PRESENT && headerOk())
+    )
+    {
+        if (readOpt() == IOobject::MUST_READ_IF_MODIFIED)
+        {
+            WarningIn
+            (
+                "extendedFeatureEdgeMesh::extendedFeatureEdgeMesh"
+                "(const IOobject&)"
+            )   << "Specified IOobject::MUST_READ_IF_MODIFIED but class"
+                << " does not support automatic rereading."
+                << endl;
+        }
+
+        Istream& is = readStream(typeName);
+
+        is  >> *this
+            >> concaveStart_
+            >> mixedStart_
+            >> nonFeatureStart_
+            >> internalStart_
+            >> flatStart_
+            >> openStart_
+            >> multipleStart_
+            >> normals_
+            >> edgeNormals_
+            >> featurePointNormals_
+            >> regionEdges_;
+
+        close();
+
+        {
+            // Calculate edgeDirections
+
+            const edgeList& eds(edges());
+
+            const pointField& pts(points());
+
+            edgeDirections_.setSize(eds.size());
+
+            forAll(eds, eI)
+            {
+                edgeDirections_[eI] = eds[eI].vec(pts);
+            }
+
+            edgeDirections_ /= mag(edgeDirections_);
+        }
+    }
+
+    if (debug)
+    {
+        Pout<< "extendedFeatureEdgeMesh::extendedFeatureEdgeMesh :"
+            << " constructed from IOobject :"
+            << " points:" << points().size()
+            << " edges:" << edges().size()
+            << endl;
+    }
+}
+
+
+Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh
+(
+    const IOobject& io,
+    const extendedFeatureEdgeMesh& fem
+)
+:
+    regIOobject(io),
+    edgeMesh(fem),
+    concaveStart_(fem.concaveStart()),
+    mixedStart_(fem.mixedStart()),
+    nonFeatureStart_(fem.nonFeatureStart()),
+    internalStart_(fem.internalStart()),
+    flatStart_(fem.flatStart()),
+    openStart_(fem.openStart()),
+    multipleStart_(fem.multipleStart()),
+    normals_(fem.normals()),
+    edgeDirections_(fem.edgeDirections()),
+    edgeNormals_(fem.edgeNormals()),
+    featurePointNormals_(fem.featurePointNormals()),
+    regionEdges_(fem.regionEdges()),
+    edgeTree_(),
+    edgeTreesByType_()
+{}
+
+
+Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh
+(
+    const IOobject& io,
+    const Xfer<pointField>& pointLst,
+    const Xfer<edgeList>& edgeLst
+)
+:
+    regIOobject(io),
+    edgeMesh(pointLst, edgeLst),
+    concaveStart_(0),
+    mixedStart_(0),
+    nonFeatureStart_(0),
+    internalStart_(0),
+    flatStart_(0),
+    openStart_(0),
+    multipleStart_(0),
+    normals_(0),
+    edgeDirections_(0),
+    edgeNormals_(0),
+    featurePointNormals_(0),
+    regionEdges_(0),
+    edgeTree_(),
+    edgeTreesByType_()
+{}
+
+
+Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh
+(
+    const surfaceFeatures& sFeat,
+    const objectRegistry& obr,
+    const fileName& sFeatFileName
+)
+:
+    regIOobject
+    (
+        IOobject
+        (
+            sFeatFileName,
+            obr.time().constant(),
+            "extendedFeatureEdgeMesh",
+            obr,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        )
+    ),
+    edgeMesh(pointField(0), edgeList(0)),
+    concaveStart_(-1),
+    mixedStart_(-1),
+    nonFeatureStart_(-1),
+    internalStart_(-1),
+    flatStart_(-1),
+    openStart_(-1),
+    multipleStart_(-1),
+    normals_(0),
+    edgeDirections_(0),
+    edgeNormals_(0),
+    featurePointNormals_(0),
+    regionEdges_(0),
+    edgeTree_(),
+    edgeTreesByType_()
+{
+    // Extract and reorder the data from surfaceFeatures
+
+    // References to the surfaceFeatures data
+    const triSurface& surf(sFeat.surface());
+    const pointField& sFeatLocalPts(surf.localPoints());
+    const edgeList& sFeatEds(surf.edges());
+
+    // Filling the extendedFeatureEdgeMesh with the raw geometrical data.
+
+    label nFeatEds = sFeat.featureEdges().size();
+
+    DynamicList<point> tmpPts;
+    edgeList eds(nFeatEds);
+    DynamicList<vector> norms;
+    vectorField edgeDirections(nFeatEds);
+    labelListList edgeNormals(nFeatEds);
+    DynamicList<label> regionEdges;
+
+    // Mapping between old and new indices, there is entry in the map for each
+    // of sFeatLocalPts, -1 means that this point hasn't been used (yet), >= 0
+    // corresponds to the index
+    labelList pointMap(sFeatLocalPts.size(), -1);
+
+    // Noting when the normal of a face has been used so not to duplicate
+    labelList faceMap(surf.size(), -1);
+
+    // Collecting the status of edge for subsequent sorting
+    List<edgeStatus> edStatus(nFeatEds, NONE);
+
+    forAll(sFeat.featurePoints(), i)
+    {
+        label sFPI = sFeat.featurePoints()[i];
+
+        tmpPts.append(sFeatLocalPts[sFPI]);
+
+        pointMap[sFPI] = tmpPts.size() - 1;
+    }
+
+    // All feature points have been added
+    nonFeatureStart_ = tmpPts.size();
+
+    forAll(sFeat.featureEdges(), i)
+    {
+        label sFEI = sFeat.featureEdges()[i];
+
+        const edge& fE(sFeatEds[sFEI]);
+
+        // Check to see if the points have been already used
+        if (pointMap[fE.start()] == -1)
+        {
+             tmpPts.append(sFeatLocalPts[fE.start()]);
+
+             pointMap[fE.start()] = tmpPts.size() - 1;
+        }
+
+        eds[i].start() = pointMap[fE.start()];
+
+        if (pointMap[fE.end()] == -1)
+        {
+            tmpPts.append(sFeatLocalPts[fE.end()]);
+
+            pointMap[fE.end()] = tmpPts.size() - 1;
+        }
+
+        eds[i].end() = pointMap[fE.end()];
+
+        // Pick up the faces adjacent to the feature edge
+        const labelList& eFaces = surf.edgeFaces()[sFEI];
+
+        edgeNormals[i].setSize(eFaces.size());
+
+        forAll(eFaces, j)
+        {
+            label eFI = eFaces[j];
+
+            // Check to see if the points have been already used
+            if (faceMap[eFI] == -1)
+            {
+                norms.append(surf.faceNormals()[eFI]);
+
+                faceMap[eFI] = norms.size() - 1;
+            }
+
+            edgeNormals[i][j] = faceMap[eFI];
+        }
+
+        vector fC0tofC1(vector::zero);
+
+        if (eFaces.size() == 2)
+        {
+            fC0tofC1 =
+                surf[eFaces[1]].centre(surf.points())
+              - surf[eFaces[0]].centre(surf.points());
+        }
+
+        edStatus[i] = classifyEdge(norms, edgeNormals[i], fC0tofC1);
+
+        edgeDirections[i] = fE.vec(sFeatLocalPts);
+
+        if (i < sFeat.nRegionEdges())
+        {
+            regionEdges.append(i);
+        }
+    }
+
+    // Reorder the edges by classification
+
+    List<DynamicList<label> > allEds(nEdgeTypes);
+
+    DynamicList<label>& externalEds(allEds[0]);
+    DynamicList<label>& internalEds(allEds[1]);
+    DynamicList<label>& flatEds(allEds[2]);
+    DynamicList<label>& openEds(allEds[3]);
+    DynamicList<label>& multipleEds(allEds[4]);
+
+    forAll(eds, i)
+    {
+        edgeStatus eStat = edStatus[i];
+
+        if (eStat == EXTERNAL)
+        {
+            externalEds.append(i);
+        }
+        else if (eStat == INTERNAL)
+        {
+            internalEds.append(i);
+        }
+        else if (eStat == FLAT)
+        {
+            flatEds.append(i);
+        }
+        else if (eStat == OPEN)
+        {
+            openEds.append(i);
+        }
+        else if (eStat == MULTIPLE)
+        {
+            multipleEds.append(i);
+        }
+        else if (eStat == NONE)
+        {
+            FatalErrorIn
+            (
+                "Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh"
+                "("
+                "    const surfaceFeatures& sFeat,"
+                "    const objectRegistry& obr,"
+                "    const fileName& sFeatFileName"
+                ")"
+            )
+                << nl << "classifyEdge returned NONE on edge "
+                << eds[i]
+                << ". There is a problem with definition of this edge."
+                << nl << abort(FatalError);
+        }
+    }
+
+    internalStart_ = externalEds.size();
+    flatStart_ = internalStart_ + internalEds.size();
+    openStart_ = flatStart_ + flatEds.size();
+    multipleStart_ = openStart_ + openEds.size();
+
+    labelList edMap
+    (
+        ListListOps::combine<labelList>
+        (
+            allEds,
+            accessOp<labelList>()
+        )
+    );
+
+    edMap = invert(edMap.size(), edMap);
+
+    inplaceReorder(edMap, eds);
+    inplaceReorder(edMap, edStatus);
+    inplaceReorder(edMap, edgeDirections);
+    inplaceReorder(edMap, edgeNormals);
+    inplaceRenumber(edMap, regionEdges);
+
+    pointField pts(tmpPts);
+
+    // Initialise the edgeMesh
+    edgeMesh::operator=(edgeMesh(pts, eds));
+
+    // Initialise sorted edge related data
+    edgeDirections_ = edgeDirections/mag(edgeDirections);
+    edgeNormals_ = edgeNormals;
+    regionEdges_ = regionEdges;
+
+    // Normals are all now found and indirectly addressed, can also be stored
+    normals_ = vectorField(norms);
+
+    // Reorder the feature points by classification
+
+    List<DynamicList<label> > allPts(3);
+
+    DynamicList<label>& convexPts(allPts[0]);
+    DynamicList<label>& concavePts(allPts[1]);
+    DynamicList<label>& mixedPts(allPts[2]);
+
+    for (label i = 0; i < nonFeatureStart_; i++)
+    {
+        pointStatus ptStatus = classifyFeaturePoint(i);
+
+        if (ptStatus == CONVEX)
+        {
+            convexPts.append(i);
+        }
+        else if (ptStatus == CONCAVE)
+        {
+            concavePts.append(i);
+        }
+        else if (ptStatus == MIXED)
+        {
+            mixedPts.append(i);
+        }
+        else if (ptStatus == NONFEATURE)
+        {
+            FatalErrorIn
+            (
+                "Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh"
+                "("
+                "    const surfaceFeatures& sFeat,"
+                "    const objectRegistry& obr,"
+                "    const fileName& sFeatFileName"
+                ")"
+            )
+                << nl << "classifyFeaturePoint returned NONFEATURE on point at "
+                << points()[i]
+                << ". There is a problem with definition of this feature point."
+                << nl << abort(FatalError);
+        }
+    }
+
+    concaveStart_ = convexPts.size();
+    mixedStart_ = concaveStart_ + concavePts.size();
+
+    labelList ftPtMap
+    (
+        ListListOps::combine<labelList>
+        (
+            allPts,
+            accessOp<labelList>()
+        )
+    );
+
+    ftPtMap = invert(ftPtMap.size(), ftPtMap);
+
+    // Creating the ptMap from the ftPtMap with identity values up to the size
+    // of pts to create an oldToNew map for inplaceReorder
+
+    labelList ptMap(identity(pts.size()));
+
+    forAll(ftPtMap, i)
+    {
+        ptMap[i] = ftPtMap[i];
+    }
+
+    inplaceReorder(ptMap, pts);
+
+    forAll(eds, i)
+    {
+        inplaceRenumber(ptMap, eds[i]);
+    }
+
+    // Reinitialise the edgeMesh with sorted feature points and
+    // renumbered edges
+    edgeMesh::operator=(edgeMesh(pts, eds));
+
+    // Generate the featurePointNormals
+
+    labelListList featurePointNormals(nonFeatureStart_);
+
+    for (label i = 0; i < nonFeatureStart_; i++)
+    {
+        DynamicList<label> tmpFtPtNorms;
+
+        const labelList& ptEds = pointEdges()[i];
+
+        forAll(ptEds, j)
+        {
+            const labelList& ptEdNorms(edgeNormals[ptEds[j]]);
+
+            forAll(ptEdNorms, k)
+            {
+                if (findIndex(tmpFtPtNorms, ptEdNorms[k]) == -1)
+                {
+                    bool addNormal = true;
+
+                    // Check that the normal direction is unique at this feature
+                    forAll(tmpFtPtNorms, q)
+                    {
+                        if
+                        (
+                            (normals_[ptEdNorms[k]] & normals_[tmpFtPtNorms[q]])
+                          > cosNormalAngleTol_
+                        )
+                        {
+                            // Parallel to an existing normal, do not add
+                            addNormal = false;
+
+                            break;
+                        }
+                    }
+
+                    if (addNormal)
+                    {
+                        tmpFtPtNorms.append(ptEdNorms[k]);
+                    }
+                }
+            }
+        }
+
+        featurePointNormals[i] = tmpFtPtNorms;
+    }
+
+    featurePointNormals_ = featurePointNormals;
+}
+
+
+Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh
+(
+    const IOobject& io,
+    const pointField& pts,
+    const edgeList& eds,
+    label concaveStart,
+    label mixedStart,
+    label nonFeatureStart,
+    label internalStart,
+    label flatStart,
+    label openStart,
+    label multipleStart,
+    const vectorField& normals,
+    const vectorField& edgeDirections,
+    const labelListList& edgeNormals,
+    const labelListList& featurePointNormals,
+    const labelList& regionEdges
+)
+:
+    regIOobject(io),
+    edgeMesh(pts, eds),
+    concaveStart_(concaveStart),
+    mixedStart_(mixedStart),
+    nonFeatureStart_(nonFeatureStart),
+    internalStart_(internalStart),
+    flatStart_(flatStart),
+    openStart_(openStart),
+    multipleStart_(multipleStart),
+    normals_(normals),
+    edgeDirections_(edgeDirections),
+    edgeNormals_(edgeNormals),
+    featurePointNormals_(featurePointNormals),
+    regionEdges_(regionEdges),
+    edgeTree_(),
+    edgeTreesByType_()
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::extendedFeatureEdgeMesh::~extendedFeatureEdgeMesh()
+{}
+
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+Foam::extendedFeatureEdgeMesh::pointStatus
+Foam::extendedFeatureEdgeMesh::classifyFeaturePoint
+(
+    label ptI
+) const
+{
+    labelList ptEds(pointEdges()[ptI]);
+
+    label nPtEds = ptEds.size();
+    label nExternal = 0;
+    label nInternal = 0;
+
+    if (nPtEds == 0)
+    {
+        // There are no edges attached to the point, this is a problem
+        return NONFEATURE;
+    }
+
+    forAll(ptEds, i)
+    {
+        edgeStatus edStat = getEdgeStatus(ptEds[i]);
+
+        if (edStat == EXTERNAL)
+        {
+            nExternal++;
+        }
+        else if (edStat == INTERNAL)
+        {
+            nInternal++;
+        }
+    }
+
+    if (nExternal == nPtEds)
+    {
+        return CONVEX;
+    }
+    else if (nInternal == nPtEds)
+    {
+        return CONCAVE;
+    }
+    else
+    {
+        return MIXED;
+    }
+}
+
+
+Foam::extendedFeatureEdgeMesh::edgeStatus
+Foam::extendedFeatureEdgeMesh::classifyEdge
+(
+    const List<vector>& norms,
+    const labelList& edNorms,
+    const vector& fC0tofC1
+) const
+{
+    label nEdNorms = edNorms.size();
+
+    if (nEdNorms == 1)
+    {
+        return OPEN;
+    }
+    else if (nEdNorms == 2)
+    {
+        const vector n0(norms[edNorms[0]]);
+        const vector n1(norms[edNorms[1]]);
+
+        if ((n0 & n1) > cosNormalAngleTol_)
+        {
+            return FLAT;
+        }
+        else if ((fC0tofC1 & n0) > 0.0)
+        {
+            return INTERNAL;
+        }
+        else
+        {
+            return EXTERNAL;
+        }
+    }
+    else if (nEdNorms > 2)
+    {
+        return MULTIPLE;
+    }
+    else
+    {
+        // There is a problem - the edge has no normals
+        return NONE;
+    }
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::extendedFeatureEdgeMesh::nearestFeatureEdge
+(
+    const point& sample,
+    scalar searchDistSqr,
+    pointIndexHit& info
+) const
+{
+    info = edgeTree().findNearest
+    (
+        sample,
+        searchDistSqr
+    );
+}
+
+
+void Foam::extendedFeatureEdgeMesh::nearestFeatureEdge
+(
+    const pointField& samples,
+    const scalarField& searchDistSqr,
+    List<pointIndexHit>& info
+) const
+{
+    info.setSize(samples.size());
+
+    forAll(samples, i)
+    {
+        nearestFeatureEdge
+        (
+            samples[i],
+            searchDistSqr[i],
+            info[i]
+        );
+    }
+}
+
+
+void Foam::extendedFeatureEdgeMesh::nearestFeatureEdgeByType
+(
+    const point& sample,
+    const scalarField& searchDistSqr,
+    List<pointIndexHit>& info
+) const
+{
+    const PtrList<indexedOctree<treeDataEdge> >& edgeTrees = edgeTreesByType();
+
+    info.setSize(edgeTrees.size());
+
+    labelList sliceStarts(edgeTrees.size());
+
+    sliceStarts[0] = externalStart_;
+    sliceStarts[1] = internalStart_;
+    sliceStarts[2] = flatStart_;
+    sliceStarts[3] = openStart_;
+    sliceStarts[4] = multipleStart_;
+
+    forAll(edgeTrees, i)
+    {
+        info[i] = edgeTrees[i].findNearest
+        (
+            sample,
+            searchDistSqr[i]
+        );
+
+        // The index returned by the indexedOctree is local to the slice of
+        // edges it was supplied with, return the index to the value in the
+        // complete edge list
+
+        info[i].setIndex(info[i].index() + sliceStarts[i]);
+    }
+}
+
+
+const Foam::indexedOctree<Foam::treeDataEdge>&
+Foam::extendedFeatureEdgeMesh::edgeTree() const
+{
+    if (edgeTree_.empty())
+    {
+        Random rndGen(17301893);
+
+        // Slightly extended bb. Slightly off-centred just so on symmetric
+        // geometry there are less face/edge aligned items.
+        treeBoundBox bb
+        (
+            treeBoundBox(points()).extend(rndGen, 1E-4)
+        );
+
+        bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+        bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+
+        labelList allEdges(identity(edges().size()));
+
+        edgeTree_.reset
+        (
+            new indexedOctree<treeDataEdge>
+            (
+                treeDataEdge
+                (
+                    false,          // cachebb
+                    edges(),        // edges
+                    points(),       // points
+                    allEdges        // selected edges
+                ),
+                bb,     // bb
+                8,      // maxLevel
+                10,     // leafsize
+                3.0     // duplicity
+            )
+        );
+    }
+
+    return edgeTree_();
+}
+
+
+const Foam::PtrList<Foam::indexedOctree<Foam::treeDataEdge> >&
+Foam::extendedFeatureEdgeMesh::edgeTreesByType() const
+{
+    if (edgeTreesByType_.size() == 0)
+    {
+        edgeTreesByType_.setSize(nEdgeTypes);
+
+        Random rndGen(872141);
+
+        // Slightly extended bb. Slightly off-centred just so on symmetric
+        // geometry there are less face/edge aligned items.
+        treeBoundBox bb
+        (
+            treeBoundBox(points()).extend(rndGen, 1E-4)
+        );
+
+        bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+        bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+
+        labelListList sliceEdges(nEdgeTypes);
+
+        // External edges
+        sliceEdges[0] =
+            identity(internalStart_ - externalStart_) + externalStart_;
+
+        // Internal edges
+        sliceEdges[1] = identity(flatStart_ - internalStart_) + internalStart_;
+
+        // Flat edges
+        sliceEdges[2] = identity(openStart_ - flatStart_) + flatStart_;
+
+        // Open edges
+        sliceEdges[3] = identity(multipleStart_ - openStart_) + openStart_;
+
+        // Multiple edges
+        sliceEdges[4] =
+            identity(edges().size() - multipleStart_) + multipleStart_;
+
+        forAll(edgeTreesByType_, i)
+        {
+            edgeTreesByType_.set
+            (
+                i,
+                new indexedOctree<treeDataEdge>
+                (
+                    treeDataEdge
+                    (
+                        false,          // cachebb
+                        edges(),        // edges
+                        points(),       // points
+                        sliceEdges[i]   // selected edges
+                    ),
+                    bb,     // bb
+                    8,      // maxLevel
+                    10,     // leafsize
+                    3.0     // duplicity
+                )
+            );
+        }
+    }
+
+    return edgeTreesByType_;
+}
+
+
+void Foam::extendedFeatureEdgeMesh::writeObj
+(
+    const fileName& prefix
+) const
+{
+    Pout<< nl << "Writing extendedFeatureEdgeMesh components to " << prefix
+        << endl;
+
+    label verti = 0;
+
+    edgeMesh::write(prefix + "_edgeMesh.obj");
+
+    OFstream convexFtPtStr(prefix + "_convexFeaturePts.obj");
+    Pout<< "Writing convex feature points to " << convexFtPtStr.name() << endl;
+
+    for(label i = 0; i < concaveStart_; i++)
+    {
+        meshTools::writeOBJ(convexFtPtStr, points()[i]);
+    }
+
+    OFstream concaveFtPtStr(prefix + "_concaveFeaturePts.obj");
+    Pout<< "Writing concave feature points to "
+        << concaveFtPtStr.name() << endl;
+
+    for(label i = concaveStart_; i < mixedStart_; i++)
+    {
+        meshTools::writeOBJ(concaveFtPtStr, points()[i]);
+    }
+
+    OFstream mixedFtPtStr(prefix + "_mixedFeaturePts.obj");
+    Pout<< "Writing mixed feature points to " << mixedFtPtStr.name() << endl;
+
+    for(label i = mixedStart_; i < nonFeatureStart_; i++)
+    {
+        meshTools::writeOBJ(mixedFtPtStr, points()[i]);
+    }
+
+    OFstream mixedFtPtStructureStr(prefix + "_mixedFeaturePtsStructure.obj");
+    Pout<< "Writing mixed feature point structure to "
+        << mixedFtPtStructureStr.name() << endl;
+
+    verti = 0;
+    for(label i = mixedStart_; i < nonFeatureStart_; i++)
+    {
+        const labelList& ptEds = pointEdges()[i];
+
+        forAll(ptEds, j)
+        {
+            const edge& e = edges()[ptEds[j]];
+
+            meshTools::writeOBJ(mixedFtPtStructureStr, points()[e[0]]); verti++;
+            meshTools::writeOBJ(mixedFtPtStructureStr, points()[e[1]]); verti++;
+            mixedFtPtStructureStr << "l " << verti-1 << ' ' << verti << endl;
+        }
+    }
+
+    OFstream externalStr(prefix + "_externalEdges.obj");
+    Pout<< "Writing external edges to " << externalStr.name() << endl;
+
+    verti = 0;
+    for (label i = externalStart_; i < internalStart_; i++)
+    {
+        const edge& e = edges()[i];
+
+        meshTools::writeOBJ(externalStr, points()[e[0]]); verti++;
+        meshTools::writeOBJ(externalStr, points()[e[1]]); verti++;
+        externalStr << "l " << verti-1 << ' ' << verti << endl;
+    }
+
+    OFstream internalStr(prefix + "_internalEdges.obj");
+    Pout<< "Writing internal edges to " << internalStr.name() << endl;
+
+    verti = 0;
+    for (label i = internalStart_; i < flatStart_; i++)
+    {
+        const edge& e = edges()[i];
+
+        meshTools::writeOBJ(internalStr, points()[e[0]]); verti++;
+        meshTools::writeOBJ(internalStr, points()[e[1]]); verti++;
+        internalStr << "l " << verti-1 << ' ' << verti << endl;
+    }
+
+    OFstream flatStr(prefix + "_flatEdges.obj");
+    Pout<< "Writing flat edges to " << flatStr.name() << endl;
+
+    verti = 0;
+    for (label i = flatStart_; i < openStart_; i++)
+    {
+        const edge& e = edges()[i];
+
+        meshTools::writeOBJ(flatStr, points()[e[0]]); verti++;
+        meshTools::writeOBJ(flatStr, points()[e[1]]); verti++;
+        flatStr << "l " << verti-1 << ' ' << verti << endl;
+    }
+
+    OFstream openStr(prefix + "_openEdges.obj");
+    Pout<< "Writing open edges to " << openStr.name() << endl;
+
+    verti = 0;
+    for (label i = openStart_; i < multipleStart_; i++)
+    {
+        const edge& e = edges()[i];
+
+        meshTools::writeOBJ(openStr, points()[e[0]]); verti++;
+        meshTools::writeOBJ(openStr, points()[e[1]]); verti++;
+        openStr << "l " << verti-1 << ' ' << verti << endl;
+    }
+
+    OFstream multipleStr(prefix + "_multipleEdges.obj");
+    Pout<< "Writing multiple edges to " << multipleStr.name() << endl;
+
+    verti = 0;
+    for (label i = multipleStart_; i < edges().size(); i++)
+    {
+        const edge& e = edges()[i];
+
+        meshTools::writeOBJ(multipleStr, points()[e[0]]); verti++;
+        meshTools::writeOBJ(multipleStr, points()[e[1]]); verti++;
+        multipleStr << "l " << verti-1 << ' ' << verti << endl;
+    }
+
+    OFstream regionStr(prefix + "_regionEdges.obj");
+    Pout<< "Writing region edges to " << regionStr.name() << endl;
+
+    verti = 0;
+    forAll(regionEdges_, i)
+    {
+        const edge& e = edges()[regionEdges_[i]];
+
+        meshTools::writeOBJ(regionStr, points()[e[0]]); verti++;
+        meshTools::writeOBJ(regionStr, points()[e[1]]); verti++;
+        regionStr << "l " << verti-1 << ' ' << verti << endl;
+    }
+}
+
+
+bool Foam::extendedFeatureEdgeMesh::writeData(Ostream& os) const
+{
+    os  << "// points, edges, concaveStart, mixedStart, nonFeatureStart, " << nl
+        << "// internalStart, flatStart, openStart, multipleStart, " << nl
+        << "// normals, edgeNormals, featurePointNormals, regionEdges" << nl
+        << endl;
+
+    os  << points() << nl
+        << edges() << nl
+        << concaveStart_ << token::SPACE
+        << mixedStart_ << token::SPACE
+        << nonFeatureStart_ << token::SPACE
+        << internalStart_ << token::SPACE
+        << flatStart_ << token::SPACE
+        << openStart_ << token::SPACE
+        << multipleStart_ << nl
+        << normals_ << nl
+        << edgeNormals_ << nl
+        << featurePointNormals_ << nl
+        << regionEdges_
+        << endl;
+
+    return os.good();
+}
+
+
+// ************************************************************************* //
diff --git a/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.H b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.H
new file mode 100644
index 0000000000000000000000000000000000000000..f1e85867d4530412fe49df70e4ebb519796a5e81
--- /dev/null
+++ b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.H
@@ -0,0 +1,377 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2004-2011 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 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::extendedFeatureEdgeMesh
+
+Description
+
+    Description of feature edges and points.
+
+    Feature points are a sorted subset at the start of the overall points list:
+        0 .. concaveStart_-1                : convex points (w.r.t normals)
+        concaveStart_-1 .. mixedStart_-1    : concave points
+        mixedStart_ .. nonFeatureStart_     : mixed internal/external points
+        nonFeatureStart_ .. size-1          : non-feature points
+
+    Feature edges are the edgeList of the edgeMesh and are sorted:
+        0 .. internalStart_-1           : external edges (convex w.r.t normals)
+        internalStart_ .. flatStart_-1  : internal edges (concave)
+        flatStart_ .. openStart_-1      : flat edges (neither concave or convex)
+                                          can arise from region interfaces on
+                                          flat surfaces
+        openStart_ .. multipleStart_-1  : open edges (e.g. from baffle surfaces)
+        multipleStart_ .. size-1        : multiply connected edges
+
+    The edge direction and feature edge and feature point adjacent normals
+    are stored.
+
+SourceFiles
+    extendedFeatureEdgeMeshI.H
+    extendedFeatureEdgeMesh.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef extendedFeatureEdgeMesh_H
+#define extendedFeatureEdgeMesh_H
+
+#include "edgeMesh.H"
+#include "surfaceFeatures.H"
+#include "objectRegistry.H"
+#include "IOdictionary.H"
+#include "indexedOctree.H"
+#include "treeDataEdge.H"
+#include "pointIndexHit.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class extendedFeatureEdgeMesh Declaration
+\*---------------------------------------------------------------------------*/
+
+class extendedFeatureEdgeMesh
+:
+    public regIOobject,
+    public edgeMesh
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("extendedFeatureEdgeMesh");
+
+    enum pointStatus
+    {
+        CONVEX,     // Fully convex point (w.r.t normals)
+        CONCAVE,    // Fully concave point
+        MIXED,      // A point surrounded by both convex and concave edges
+        NONFEATURE  // Not a feature point
+    };
+
+    enum edgeStatus
+    {
+        EXTERNAL,   // "Convex" edge
+        INTERNAL,   // "Concave" edge
+        FLAT,       // Neither concave or convex, on a flat surface
+        OPEN,       // i.e. only connected to one face
+        MULTIPLE,   // Multiply connected (connected to more than two faces)
+        NONE        // Not a classified feature edge (consistency with
+                    // surfaceFeatures)
+    };
+
+private:
+
+    // Static data
+
+        //- Angular closeness tolerance for treating normals as the same
+        static scalar cosNormalAngleTol_;
+
+        //- Index of the start of the convex feature points - static as 0
+        static label convexStart_;
+
+        //- Index of the start of the external feature edges - static as 0
+        static label externalStart_;
+
+
+    // Private data
+
+        //- Index of the start of the concave feature points
+        label concaveStart_;
+
+        //- Index of the start of the mixed type feature points
+        label mixedStart_;
+
+        //- Index of the start of the non-feature points
+        label nonFeatureStart_;
+
+        //- Index of the start of the internal feature edges
+        label internalStart_;
+
+        //- Index of the start of the flat feature edges
+        label flatStart_;
+
+        //- Index of the start of the open feature edges
+        label openStart_;
+
+        //- Index of the start of the multiply-connected feature edges
+        label multipleStart_;
+
+        //- Normals of the features, to be referred to by index by both feature
+        //  points and edges, unsorted
+        vectorField normals_;
+
+        //- Flat and open edges require the direction of the edge
+        vectorField edgeDirections_;
+
+        //- Indices of the normals that are adjacent to the feature edges
+        labelListList edgeNormals_;
+
+        //- Indices of the normals that are adjacent to the feature points
+        labelListList featurePointNormals_;
+
+        //- Feature edges which are on the boundary between regions
+        labelList regionEdges_;
+
+        //- Search tree for all edges
+        mutable autoPtr<indexedOctree<treeDataEdge> > edgeTree_;
+
+        //- Individual search trees for each type of edge
+        mutable PtrList<indexedOctree<treeDataEdge> > edgeTreesByType_;
+
+
+    // Private Member Functions
+
+        //- Classify the type of feature point.  Requires valid stored member
+        //  data for edges and normals.
+        pointStatus classifyFeaturePoint(label ptI) const;
+
+        //- Classify the type of feature edge.  Requires face centre 0 to face
+        //  centre 1 vector to distinguish internal from external
+        edgeStatus classifyEdge
+        (
+            const List<vector>& norms,
+            const labelList& edNorms,
+            const vector& fC0tofC1
+        ) const;
+
+
+public:
+
+    // Static data
+
+        //- Number of possible point types (i.e. number of slices)
+        static label nPointTypes;
+
+        //- Number of possible feature edge types (i.e. number of slices)
+        static label nEdgeTypes;
+
+    // Constructors
+
+        //- Construct (read) given an IOobject
+        extendedFeatureEdgeMesh(const IOobject&);
+
+        //- Construct as copy
+        explicit extendedFeatureEdgeMesh
+        (
+            const IOobject&,
+            const extendedFeatureEdgeMesh&
+        );
+
+        //- Construct by transferring components (points, edges)
+        extendedFeatureEdgeMesh
+        (
+            const IOobject&,
+            const Xfer<pointField>&,
+            const Xfer<edgeList>&
+        );
+
+        //- Construct (read) given surfaceFeatures, an objectRegistry and a
+        //  fileName to write to.  Extracts, classifies and reorders the data
+        //  from surfaceFeatures.
+        extendedFeatureEdgeMesh
+        (
+            const surfaceFeatures& sFeat,
+            const objectRegistry& obr,
+            const fileName& sFeatFileName
+        );
+
+        //- Construct from all components
+        extendedFeatureEdgeMesh
+        (
+            const IOobject& io,
+            const pointField& pts,
+            const edgeList& eds,
+            label concaveStart,
+            label mixedStart,
+            label nonFeatureStart,
+            label internalStart,
+            label flatStart,
+            label openStart,
+            label multipleStart,
+            const vectorField& normals,
+            const vectorField& edgeDirections,
+            const labelListList& edgeNormals,
+            const labelListList& featurePointNormals,
+            const labelList& regionEdges
+        );
+
+
+    //- Destructor
+    ~extendedFeatureEdgeMesh();
+
+
+    // Member Functions
+
+        // Find
+
+            //- Find nearest surface edge for the sample point.
+            void nearestFeatureEdge
+            (
+                const point& sample,
+                scalar searchDistSqr,
+                pointIndexHit& info
+            ) const;
+
+            //- Find nearest surface edge for each sample point.
+            void nearestFeatureEdge
+            (
+                const pointField& samples,
+                const scalarField& searchDistSqr,
+                List<pointIndexHit>& info
+            ) const;
+
+            //- Find the nearest point on each type of feature edge
+            void nearestFeatureEdgeByType
+            (
+                const point& sample,
+                const scalarField& searchDistSqr,
+                List<pointIndexHit>& info
+            ) const;
+
+        // Access
+
+            //- Return the index of the start of the convex feature points
+            inline label convexStart() const;
+
+            //- Return the index of the start of the concave feature points
+            inline label concaveStart() const;
+
+            //- Return the index of the start of the mixed type feature points
+            inline label mixedStart() const;
+
+            //- Return the index of the start of the non-feature points
+            inline label nonFeatureStart() const;
+
+            //- Return the index of the start of the external feature edges
+            inline label externalStart() const;
+
+            //- Return the index of the start of the internal feature edges
+            inline label internalStart() const;
+
+            //- Return the index of the start of the flat feature edges
+            inline label flatStart() const;
+
+            //- Return the index of the start of the open feature edges
+            inline label openStart() const;
+
+            //- Return the index of the start of the multiply-connected feature
+            //  edges
+            inline label multipleStart() const;
+
+            //- Return whether or not the point index is a feature point
+            inline bool featurePoint(label ptI) const;
+
+            //- Return the normals of the surfaces adjacent to the feature edges
+            //  and points
+            inline const vectorField& normals() const;
+
+            //- Return the edgeDirection vectors
+            inline const vectorField& edgeDirections() const;
+
+            //- Return the direction of edgeI, pointing away from ptI
+            inline vector edgeDirection(label edgeI, label ptI) const;
+
+            //- Return the indices of the normals that are adjacent to the
+            //  feature edges
+            inline const labelListList& edgeNormals() const;
+
+            //- Return the normal vectors for a given set of normal indices
+            inline vectorField edgeNormals(const labelList& edgeNormIs) const;
+
+            //- Return the normal vectors for a given edge
+            inline vectorField edgeNormals(label edgeI) const;
+
+            //- Return the indices of the normals that are adjacent to the
+            //  feature points
+            inline const labelListList& featurePointNormals() const;
+
+            //- Return the normal vectors for a given feature point
+            inline vectorField featurePointNormals(label ptI) const;
+
+            //- Return the feature edges which are on the boundary between
+            //  regions
+            inline const labelList& regionEdges() const;
+
+            //- Return the pointStatus of a specified point
+            inline pointStatus getPointStatus(label ptI) const;
+
+            //- Return the edgeStatus of a specified edge
+            inline edgeStatus getEdgeStatus(label edgeI) const;
+
+            //- Demand driven construction of octree for boundary edges
+            const indexedOctree<treeDataEdge>& edgeTree() const;
+
+            //- Demand driven construction of octree for boundary edges by type
+            const PtrList<indexedOctree<treeDataEdge> >&
+            edgeTreesByType() const;
+
+
+        // Write
+
+            //- Write all components of the extendedFeatureEdgeMesh as obj files
+            void writeObj(const fileName& prefix) const;
+
+            //- Give precedence to the regIOobject write
+            using regIOobject::write;
+
+            //- WriteData function required for regIOobject write operation
+            virtual bool writeData(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "extendedFeatureEdgeMeshI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/edgeMesh/featureEdgeMesh/featureEdgeMeshI.H b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMeshI.H
similarity index 66%
rename from src/edgeMesh/featureEdgeMesh/featureEdgeMeshI.H
rename to src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMeshI.H
index 6aac568b03cffc3dd8ffd2ca4bf3c634ed14a1b7..4feb72e6867bad53f7867a7d238d8dffe7b06b7f 100644
--- a/src/edgeMesh/featureEdgeMesh/featureEdgeMeshI.H
+++ b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMeshI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2009-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2009-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,78 +25,79 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline Foam::label Foam::featureEdgeMesh::convexStart() const
+inline Foam::label Foam::extendedFeatureEdgeMesh::convexStart() const
 {
     return convexStart_;
 }
 
 
-inline Foam::label Foam::featureEdgeMesh::concaveStart() const
+inline Foam::label Foam::extendedFeatureEdgeMesh::concaveStart() const
 {
     return concaveStart_;
 }
 
 
-inline Foam::label Foam::featureEdgeMesh::mixedStart() const
+inline Foam::label Foam::extendedFeatureEdgeMesh::mixedStart() const
 {
     return mixedStart_;
 }
 
 
-inline Foam::label Foam::featureEdgeMesh::nonFeatureStart() const
+inline Foam::label Foam::extendedFeatureEdgeMesh::nonFeatureStart() const
 {
     return nonFeatureStart_;
 }
 
 
-inline Foam::label Foam::featureEdgeMesh::externalStart() const
+inline Foam::label Foam::extendedFeatureEdgeMesh::externalStart() const
 {
     return externalStart_;
 }
 
 
-inline Foam::label Foam::featureEdgeMesh::internalStart() const
+inline Foam::label Foam::extendedFeatureEdgeMesh::internalStart() const
 {
     return internalStart_;
 }
 
 
-inline Foam::label Foam::featureEdgeMesh::flatStart() const
+inline Foam::label Foam::extendedFeatureEdgeMesh::flatStart() const
 {
     return flatStart_;
 }
 
 
-inline Foam::label Foam::featureEdgeMesh::openStart() const
+inline Foam::label Foam::extendedFeatureEdgeMesh::openStart() const
 {
     return openStart_;
 }
 
 
-inline Foam::label Foam::featureEdgeMesh::multipleStart() const
+inline Foam::label Foam::extendedFeatureEdgeMesh::multipleStart() const
 {
     return multipleStart_;
 }
 
 
-inline bool Foam::featureEdgeMesh::featurePoint(label ptI) const
+inline bool Foam::extendedFeatureEdgeMesh::featurePoint(label ptI) const
 {
     return ptI < nonFeatureStart_;
 }
 
 
-inline const Foam::vectorField& Foam::featureEdgeMesh::normals() const
+inline const Foam::vectorField& Foam::extendedFeatureEdgeMesh::normals() const
 {
     return normals_;
 }
 
-inline const Foam::vectorField& Foam::featureEdgeMesh::edgeDirections() const
+inline const Foam::vectorField& Foam::extendedFeatureEdgeMesh::edgeDirections()
+const
 {
     return edgeDirections_;
 }
 
 
-inline Foam::vector Foam::featureEdgeMesh::edgeDirection
+inline Foam::vector Foam::extendedFeatureEdgeMesh::edgeDirection
 (
     label edgeI,
     label ptI
@@ -114,7 +115,7 @@ inline Foam::vector Foam::featureEdgeMesh::edgeDirection
     }
     else
     {
-        FatalErrorIn("Foam::featureEdgeMesh::edgedirection")
+        FatalErrorIn("Foam::extendedFeatureEdgeMesh::edgedirection")
             << "Requested ptI " << ptI << " is not a point on the requested "
             << "edgeI " << edgeI << ". edgeI start and end: "
             << e.start() << " " << e.end()
@@ -125,13 +126,14 @@ inline Foam::vector Foam::featureEdgeMesh::edgeDirection
 }
 
 
-inline const Foam::labelListList& Foam::featureEdgeMesh::edgeNormals() const
+inline const Foam::labelListList& Foam::extendedFeatureEdgeMesh::edgeNormals()
+const
 {
     return edgeNormals_;
 }
 
 
-inline Foam::vectorField Foam::featureEdgeMesh::edgeNormals
+inline Foam::vectorField Foam::extendedFeatureEdgeMesh::edgeNormals
 (
     const labelList& edgeNormIs
 ) const
@@ -147,27 +149,28 @@ inline Foam::vectorField Foam::featureEdgeMesh::edgeNormals
 }
 
 
-inline Foam::vectorField Foam::featureEdgeMesh::edgeNormals(label edgeI) const
+inline Foam::vectorField Foam::extendedFeatureEdgeMesh::edgeNormals(label edgeI)
+const
 {
     return edgeNormals(edgeNormals_[edgeI]);
 }
 
 
 inline const Foam::labelListList&
-Foam::featureEdgeMesh::featurePointNormals() const
+Foam::extendedFeatureEdgeMesh::featurePointNormals() const
 {
     return featurePointNormals_;
 }
 
 
-inline Foam::vectorField Foam::featureEdgeMesh::featurePointNormals
+inline Foam::vectorField Foam::extendedFeatureEdgeMesh::featurePointNormals
 (
     label ptI
 ) const
 {
     if (!featurePoint(ptI))
     {
-        WarningIn("vectorField featureEdgeMesh::featurePointNormals")
+        WarningIn("vectorField extendedFeatureEdgeMesh::featurePointNormals")
             << "Requesting the normals of a non-feature point. "
             << "Returned zero length vectorField."
             << endl;
@@ -188,13 +191,14 @@ inline Foam::vectorField Foam::featureEdgeMesh::featurePointNormals
 }
 
 
-inline const Foam::labelList& Foam::featureEdgeMesh::regionEdges() const
+inline const Foam::labelList& Foam::extendedFeatureEdgeMesh::regionEdges() const
 {
     return regionEdges_;
 }
 
 
-inline Foam::featureEdgeMesh::pointStatus Foam::featureEdgeMesh::getPointStatus
+inline Foam::extendedFeatureEdgeMesh::pointStatus
+Foam::extendedFeatureEdgeMesh::getPointStatus
 (
     label ptI
 ) const
@@ -218,7 +222,8 @@ inline Foam::featureEdgeMesh::pointStatus Foam::featureEdgeMesh::getPointStatus
 }
 
 
-inline Foam::featureEdgeMesh::edgeStatus Foam::featureEdgeMesh::getEdgeStatus
+inline Foam::extendedFeatureEdgeMesh::edgeStatus
+Foam::extendedFeatureEdgeMesh::getEdgeStatus
 (
     label edgeI
 ) const
diff --git a/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.C b/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.C
index 78ee016cfd7be0c4be6beec6177d12b5660f346a..5b3cea13f43a93b1376578454820c011f06d1ff9 100644
--- a/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.C
+++ b/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,34 +24,15 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "featureEdgeMesh.H"
-#include "triSurface.H"
-#include "Random.H"
-#include "Time.H"
-#include "meshTools.H"
-#include "linePointRef.H"
-#include "ListListOps.H"
-#include "OFstream.H"
-#include "IFstream.H"
-#include "unitConversion.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-defineTypeNameAndDebug(Foam::featureEdgeMesh, 0);
-
-Foam::scalar Foam::featureEdgeMesh::cosNormalAngleTol_ =
-    Foam::cos(degToRad(0.1));
-
-
-Foam::label Foam::featureEdgeMesh::convexStart_ = 0;
-
-
-Foam::label Foam::featureEdgeMesh::externalStart_ = 0;
-
-
-Foam::label Foam::featureEdgeMesh::nPointTypes = 4;
+namespace Foam
+{
 
+defineTypeNameAndDebug(featureEdgeMesh, 0);
 
-Foam::label Foam::featureEdgeMesh::nEdgeTypes = 5;
+}
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -59,21 +40,7 @@ Foam::label Foam::featureEdgeMesh::nEdgeTypes = 5;
 Foam::featureEdgeMesh::featureEdgeMesh(const IOobject& io)
 :
     regIOobject(io),
-    edgeMesh(pointField(0), edgeList(0)),
-    concaveStart_(0),
-    mixedStart_(0),
-    nonFeatureStart_(0),
-    internalStart_(0),
-    flatStart_(0),
-    openStart_(0),
-    multipleStart_(0),
-    normals_(0),
-    edgeDirections_(0),
-    edgeNormals_(0),
-    featurePointNormals_(0),
-    regionEdges_(0),
-    edgeTree_(),
-    edgeTreesByType_()
+    edgeMesh(pointField(0), edgeList(0))
 {
     if
     (
@@ -82,47 +49,8 @@ Foam::featureEdgeMesh::featureEdgeMesh(const IOobject& io)
      || (io.readOpt() == IOobject::READ_IF_PRESENT && headerOk())
     )
     {
-        if (readOpt() == IOobject::MUST_READ_IF_MODIFIED)
-        {
-            WarningIn("featureEdgeMesh::featureEdgeMesh(const IOobject&)")
-                << "Specified IOobject::MUST_READ_IF_MODIFIED but class"
-                << " does not support automatic rereading."
-                << endl;
-        }
-
-        Istream& is = readStream(typeName);
-
-        is  >> *this
-            >> concaveStart_
-            >> mixedStart_
-            >> nonFeatureStart_
-            >> internalStart_
-            >> flatStart_
-            >> openStart_
-            >> multipleStart_
-            >> normals_
-            >> edgeNormals_
-            >> featurePointNormals_
-            >> regionEdges_;
-
+        readStream(typeName) >> *this;
         close();
-
-        {
-            // Calculate edgeDirections
-
-            const edgeList& eds(edges());
-
-            const pointField& pts(points());
-
-            edgeDirections_.setSize(eds.size());
-
-            forAll(eds, eI)
-            {
-                edgeDirections_[eI] = eds[eI].vec(pts);
-            }
-
-            edgeDirections_ /= mag(edgeDirections_);
-        }
     }
 
     if (debug)
@@ -136,884 +64,43 @@ Foam::featureEdgeMesh::featureEdgeMesh(const IOobject& io)
 }
 
 
+//- Construct from components
 Foam::featureEdgeMesh::featureEdgeMesh
 (
     const IOobject& io,
-    const featureEdgeMesh& fem
+    const pointField& points,
+    const edgeList& edges
 )
 :
     regIOobject(io),
-    edgeMesh(fem),
-    concaveStart_(fem.concaveStart()),
-    mixedStart_(fem.mixedStart()),
-    nonFeatureStart_(fem.nonFeatureStart()),
-    internalStart_(fem.internalStart()),
-    flatStart_(fem.flatStart()),
-    openStart_(fem.openStart()),
-    multipleStart_(fem.multipleStart()),
-    normals_(fem.normals()),
-    edgeDirections_(fem.edgeDirections()),
-    edgeNormals_(fem.edgeNormals()),
-    featurePointNormals_(fem.featurePointNormals()),
-    regionEdges_(fem.regionEdges()),
-    edgeTree_(),
-    edgeTreesByType_()
+    edgeMesh(points, edges)
 {}
 
 
+// Construct as copy
 Foam::featureEdgeMesh::featureEdgeMesh
 (
     const IOobject& io,
-    const Xfer<pointField>& pointLst,
-    const Xfer<edgeList>& edgeLst
-)
-:
-    regIOobject(io),
-    edgeMesh(pointLst, edgeLst),
-    concaveStart_(0),
-    mixedStart_(0),
-    nonFeatureStart_(0),
-    internalStart_(0),
-    flatStart_(0),
-    openStart_(0),
-    multipleStart_(0),
-    normals_(0),
-    edgeDirections_(0),
-    edgeNormals_(0),
-    featurePointNormals_(0),
-    regionEdges_(0),
-    edgeTree_(),
-    edgeTreesByType_()
-{}
-
-
-Foam::featureEdgeMesh::featureEdgeMesh
-(
-    const surfaceFeatures& sFeat,
-    const objectRegistry& obr,
-    const fileName& sFeatFileName
-)
-:
-    regIOobject
-    (
-        IOobject
-        (
-            sFeatFileName,
-            obr.time().constant(),
-            "featureEdgeMesh",
-            obr,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        )
-    ),
-    edgeMesh(pointField(0), edgeList(0)),
-    concaveStart_(-1),
-    mixedStart_(-1),
-    nonFeatureStart_(-1),
-    internalStart_(-1),
-    flatStart_(-1),
-    openStart_(-1),
-    multipleStart_(-1),
-    normals_(0),
-    edgeDirections_(0),
-    edgeNormals_(0),
-    featurePointNormals_(0),
-    regionEdges_(0),
-    edgeTree_(),
-    edgeTreesByType_()
-{
-    // Extract and reorder the data from surfaceFeatures
-
-    // References to the surfaceFeatures data
-    const triSurface& surf(sFeat.surface());
-    const pointField& sFeatLocalPts(surf.localPoints());
-    const edgeList& sFeatEds(surf.edges());
-
-    // Filling the featureEdgeMesh with the raw geometrical data.
-
-    label nFeatEds = sFeat.featureEdges().size();
-
-    DynamicList<point> tmpPts;
-    edgeList eds(nFeatEds);
-    DynamicList<vector> norms;
-    vectorField edgeDirections(nFeatEds);
-    labelListList edgeNormals(nFeatEds);
-    DynamicList<label> regionEdges;
-
-    // Mapping between old and new indices, there is entry in the map for each
-    // of sFeatLocalPts, -1 means that this point hasn't been used (yet), >= 0
-    // corresponds to the index
-    labelList pointMap(sFeatLocalPts.size(), -1);
-
-    // Noting when the normal of a face has been used so not to duplicate
-    labelList faceMap(surf.size(), -1);
-
-    // Collecting the status of edge for subsequent sorting
-    List<edgeStatus> edStatus(nFeatEds, NONE);
-
-    forAll(sFeat.featurePoints(), i)
-    {
-        label sFPI = sFeat.featurePoints()[i];
-
-        tmpPts.append(sFeatLocalPts[sFPI]);
-
-        pointMap[sFPI] = tmpPts.size() - 1;
-    }
-
-    // All feature points have been added
-    nonFeatureStart_ = tmpPts.size();
-
-    forAll(sFeat.featureEdges(), i)
-    {
-        label sFEI = sFeat.featureEdges()[i];
-
-        const edge& fE(sFeatEds[sFEI]);
-
-        // Check to see if the points have been already used
-        if (pointMap[fE.start()] == -1)
-        {
-             tmpPts.append(sFeatLocalPts[fE.start()]);
-
-             pointMap[fE.start()] = tmpPts.size() - 1;
-        }
-
-        eds[i].start() = pointMap[fE.start()];
-
-        if (pointMap[fE.end()] == -1)
-        {
-            tmpPts.append(sFeatLocalPts[fE.end()]);
-
-            pointMap[fE.end()] = tmpPts.size() - 1;
-        }
-
-        eds[i].end() = pointMap[fE.end()];
-
-        // Pick up the faces adjacent to the feature edge
-        const labelList& eFaces = surf.edgeFaces()[sFEI];
-
-        edgeNormals[i].setSize(eFaces.size());
-
-        forAll(eFaces, j)
-        {
-            label eFI = eFaces[j];
-
-            // Check to see if the points have been already used
-            if (faceMap[eFI] == -1)
-            {
-                norms.append(surf.faceNormals()[eFI]);
-
-                faceMap[eFI] = norms.size() - 1;
-            }
-
-            edgeNormals[i][j] = faceMap[eFI];
-        }
-
-        vector fC0tofC1(vector::zero);
-
-        if (eFaces.size() == 2)
-        {
-            fC0tofC1 =
-                surf[eFaces[1]].centre(surf.points())
-              - surf[eFaces[0]].centre(surf.points());
-        }
-
-        edStatus[i] = classifyEdge(norms, edgeNormals[i], fC0tofC1);
-
-        edgeDirections[i] = fE.vec(sFeatLocalPts);
-
-        if (i < sFeat.nRegionEdges())
-        {
-            regionEdges.append(i);
-        }
-    }
-
-    // Reorder the edges by classification
-
-    List<DynamicList<label> > allEds(nEdgeTypes);
-
-    DynamicList<label>& externalEds(allEds[0]);
-    DynamicList<label>& internalEds(allEds[1]);
-    DynamicList<label>& flatEds(allEds[2]);
-    DynamicList<label>& openEds(allEds[3]);
-    DynamicList<label>& multipleEds(allEds[4]);
-
-    forAll(eds, i)
-    {
-        edgeStatus eStat = edStatus[i];
-
-        if (eStat == EXTERNAL)
-        {
-            externalEds.append(i);
-        }
-        else if (eStat == INTERNAL)
-        {
-            internalEds.append(i);
-        }
-        else if (eStat == FLAT)
-        {
-            flatEds.append(i);
-        }
-        else if (eStat == OPEN)
-        {
-            openEds.append(i);
-        }
-        else if (eStat == MULTIPLE)
-        {
-            multipleEds.append(i);
-        }
-        else if (eStat == NONE)
-        {
-            FatalErrorIn
-            (
-                "Foam::featureEdgeMesh::featureEdgeMesh"
-                "("
-                "    const surfaceFeatures& sFeat,"
-                "    const objectRegistry& obr,"
-                "    const fileName& sFeatFileName"
-                ")"
-            )
-                << nl << "classifyEdge returned NONE on edge "
-                << eds[i]
-                << ". There is a problem with definition of this edge."
-                << nl << abort(FatalError);
-        }
-    }
-
-    internalStart_ = externalEds.size();
-    flatStart_ = internalStart_ + internalEds.size();
-    openStart_ = flatStart_ + flatEds.size();
-    multipleStart_ = openStart_ + openEds.size();
-
-    labelList edMap
-    (
-        ListListOps::combine<labelList>
-        (
-            allEds,
-            accessOp<labelList>()
-        )
-    );
-
-    edMap = invert(edMap.size(), edMap);
-
-    inplaceReorder(edMap, eds);
-    inplaceReorder(edMap, edStatus);
-    inplaceReorder(edMap, edgeDirections);
-    inplaceReorder(edMap, edgeNormals);
-    inplaceRenumber(edMap, regionEdges);
-
-    pointField pts(tmpPts);
-
-    // Initialise the edgeMesh
-    edgeMesh::operator=(edgeMesh(pts, eds));
-
-    // Initialise sorted edge related data
-    edgeDirections_ = edgeDirections/mag(edgeDirections);
-    edgeNormals_ = edgeNormals;
-    regionEdges_ = regionEdges;
-
-    // Normals are all now found and indirectly addressed, can also be stored
-    normals_ = vectorField(norms);
-
-    // Reorder the feature points by classification
-
-    List<DynamicList<label> > allPts(3);
-
-    DynamicList<label>& convexPts(allPts[0]);
-    DynamicList<label>& concavePts(allPts[1]);
-    DynamicList<label>& mixedPts(allPts[2]);
-
-    for (label i = 0; i < nonFeatureStart_; i++)
-    {
-        pointStatus ptStatus = classifyFeaturePoint(i);
-
-        if (ptStatus == CONVEX)
-        {
-            convexPts.append(i);
-        }
-        else if (ptStatus == CONCAVE)
-        {
-            concavePts.append(i);
-        }
-        else if (ptStatus == MIXED)
-        {
-            mixedPts.append(i);
-        }
-        else if (ptStatus == NONFEATURE)
-        {
-            FatalErrorIn
-            (
-                "Foam::featureEdgeMesh::featureEdgeMesh"
-                "("
-                "    const surfaceFeatures& sFeat,"
-                "    const objectRegistry& obr,"
-                "    const fileName& sFeatFileName"
-                ")"
-            )
-                << nl << "classifyFeaturePoint returned NONFEATURE on point at "
-                << points()[i]
-                << ". There is a problem with definition of this feature point."
-                << nl << abort(FatalError);
-        }
-    }
-
-    concaveStart_ = convexPts.size();
-    mixedStart_ = concaveStart_ + concavePts.size();
-
-    labelList ftPtMap
-    (
-        ListListOps::combine<labelList>
-        (
-            allPts,
-            accessOp<labelList>()
-        )
-    );
-
-    ftPtMap = invert(ftPtMap.size(), ftPtMap);
-
-    // Creating the ptMap from the ftPtMap with identity values up to the size
-    // of pts to create an oldToNew map for inplaceReorder
-
-    labelList ptMap(identity(pts.size()));
-
-    forAll(ftPtMap, i)
-    {
-        ptMap[i] = ftPtMap[i];
-    }
-
-    inplaceReorder(ptMap, pts);
-
-    forAll(eds, i)
-    {
-        inplaceRenumber(ptMap, eds[i]);
-    }
-
-    // Reinitialise the edgeMesh with sorted feature points and
-    // renumbered edges
-    edgeMesh::operator=(edgeMesh(pts, eds));
-
-    // Generate the featurePointNormals
-
-    labelListList featurePointNormals(nonFeatureStart_);
-
-    for (label i = 0; i < nonFeatureStart_; i++)
-    {
-        DynamicList<label> tmpFtPtNorms;
-
-        const labelList& ptEds = pointEdges()[i];
-
-        forAll(ptEds, j)
-        {
-            const labelList& ptEdNorms(edgeNormals[ptEds[j]]);
-
-            forAll(ptEdNorms, k)
-            {
-                if (findIndex(tmpFtPtNorms, ptEdNorms[k]) == -1)
-                {
-                    bool addNormal = true;
-
-                    // Check that the normal direction is unique at this feature
-                    forAll(tmpFtPtNorms, q)
-                    {
-                        if
-                        (
-                            (normals_[ptEdNorms[k]] & normals_[tmpFtPtNorms[q]])
-                          > cosNormalAngleTol_
-                        )
-                        {
-                            // Parallel to an existing normal, do not add
-                            addNormal = false;
-
-                            break;
-                        }
-                    }
-
-                    if (addNormal)
-                    {
-                        tmpFtPtNorms.append(ptEdNorms[k]);
-                    }
-                }
-            }
-        }
-
-        featurePointNormals[i] = tmpFtPtNorms;
-    }
-
-    featurePointNormals_ = featurePointNormals;
-}
-
-
-Foam::featureEdgeMesh::featureEdgeMesh
-(
-    const IOobject& io,
-    const pointField& pts,
-    const edgeList& eds,
-    label concaveStart,
-    label mixedStart,
-    label nonFeatureStart,
-    label internalStart,
-    label flatStart,
-    label openStart,
-    label multipleStart,
-    const vectorField& normals,
-    const vectorField& edgeDirections,
-    const labelListList& edgeNormals,
-    const labelListList& featurePointNormals,
-    const labelList& regionEdges
+    const featureEdgeMesh& em
 )
 :
     regIOobject(io),
-    edgeMesh(pts, eds),
-    concaveStart_(concaveStart),
-    mixedStart_(mixedStart),
-    nonFeatureStart_(nonFeatureStart),
-    internalStart_(internalStart),
-    flatStart_(flatStart),
-    openStart_(openStart),
-    multipleStart_(multipleStart),
-    normals_(normals),
-    edgeDirections_(edgeDirections),
-    edgeNormals_(edgeNormals),
-    featurePointNormals_(featurePointNormals),
-    regionEdges_(regionEdges),
-    edgeTree_(),
-    edgeTreesByType_()
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::featureEdgeMesh::~featureEdgeMesh()
+    edgeMesh(em)
 {}
 
 
-// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::featureEdgeMesh::pointStatus Foam::featureEdgeMesh::classifyFeaturePoint
-(
-    label ptI
-) const
+bool Foam::featureEdgeMesh::readData(Istream& is)
 {
-    labelList ptEds(pointEdges()[ptI]);
-
-    label nPtEds = ptEds.size();
-    label nExternal = 0;
-    label nInternal = 0;
-
-    if (nPtEds == 0)
-    {
-        // There are no edges attached to the point, this is a problem
-        return NONFEATURE;
-    }
-
-    forAll(ptEds, i)
-    {
-        edgeStatus edStat = getEdgeStatus(ptEds[i]);
-
-        if (edStat == EXTERNAL)
-        {
-            nExternal++;
-        }
-        else if (edStat == INTERNAL)
-        {
-            nInternal++;
-        }
-    }
-
-    if (nExternal == nPtEds)
-    {
-        return CONVEX;
-    }
-    else if (nInternal == nPtEds)
-    {
-        return CONCAVE;
-    }
-    else
-    {
-        return MIXED;
-    }
-}
-
-
-Foam::featureEdgeMesh::edgeStatus Foam::featureEdgeMesh::classifyEdge
-(
-    const List<vector>& norms,
-    const labelList& edNorms,
-    const vector& fC0tofC1
-) const
-{
-    label nEdNorms = edNorms.size();
-
-    if (nEdNorms == 1)
-    {
-        return OPEN;
-    }
-    else if (nEdNorms == 2)
-    {
-        const vector n0(norms[edNorms[0]]);
-        const vector n1(norms[edNorms[1]]);
-
-        if ((n0 & n1) > cosNormalAngleTol_)
-        {
-            return FLAT;
-        }
-        else if ((fC0tofC1 & n0) > 0.0)
-        {
-            return INTERNAL;
-        }
-        else
-        {
-            return EXTERNAL;
-        }
-    }
-    else if (nEdNorms > 2)
-    {
-        return MULTIPLE;
-    }
-    else
-    {
-        // There is a problem - the edge has no normals
-        return NONE;
-    }
-}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-void Foam::featureEdgeMesh::nearestFeatureEdge
-(
-    const point& sample,
-    scalar searchDistSqr,
-    pointIndexHit& info
-) const
-{
-    info = edgeTree().findNearest
-    (
-        sample,
-        searchDistSqr
-    );
-}
-
-
-void Foam::featureEdgeMesh::nearestFeatureEdge
-(
-    const pointField& samples,
-    const scalarField& searchDistSqr,
-    List<pointIndexHit>& info
-) const
-{
-    info.setSize(samples.size());
-
-    forAll(samples, i)
-    {
-        nearestFeatureEdge
-        (
-            samples[i],
-            searchDistSqr[i],
-            info[i]
-        );
-    }
-}
-
-
-void Foam::featureEdgeMesh::nearestFeatureEdgeByType
-(
-    const point& sample,
-    const scalarField& searchDistSqr,
-    List<pointIndexHit>& info
-) const
-{
-    const PtrList<indexedOctree<treeDataEdge> >& edgeTrees = edgeTreesByType();
-
-    info.setSize(edgeTrees.size());
-
-    labelList sliceStarts(edgeTrees.size());
-
-    sliceStarts[0] = externalStart_;
-    sliceStarts[1] = internalStart_;
-    sliceStarts[2] = flatStart_;
-    sliceStarts[3] = openStart_;
-    sliceStarts[4] = multipleStart_;
-
-    forAll(edgeTrees, i)
-    {
-        info[i] = edgeTrees[i].findNearest
-        (
-            sample,
-            searchDistSqr[i]
-        );
-
-        // The index returned by the indexedOctree is local to the slice of
-        // edges it was supplied with, return the index to the value in the
-        // complete edge list
-
-        info[i].setIndex(info[i].index() + sliceStarts[i]);
-    }
-}
-
-
-const Foam::indexedOctree<Foam::treeDataEdge>&
-Foam::featureEdgeMesh::edgeTree() const
-{
-    if (edgeTree_.empty())
-    {
-        Random rndGen(17301893);
-
-        // Slightly extended bb. Slightly off-centred just so on symmetric
-        // geometry there are less face/edge aligned items.
-        treeBoundBox bb
-        (
-            treeBoundBox(points()).extend(rndGen, 1E-4)
-        );
-
-        bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-        bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-
-        labelList allEdges(identity(edges().size()));
-
-        edgeTree_.reset
-        (
-            new indexedOctree<treeDataEdge>
-            (
-                treeDataEdge
-                (
-                    false,          // cachebb
-                    edges(),        // edges
-                    points(),       // points
-                    allEdges        // selected edges
-                ),
-                bb,     // bb
-                8,      // maxLevel
-                10,     // leafsize
-                3.0     // duplicity
-            )
-        );
-    }
-
-    return edgeTree_();
-}
-
-
-const Foam::PtrList<Foam::indexedOctree<Foam::treeDataEdge> >&
-Foam::featureEdgeMesh::edgeTreesByType() const
-{
-    if (edgeTreesByType_.size() == 0)
-    {
-        edgeTreesByType_.setSize(nEdgeTypes);
-
-        Random rndGen(872141);
-
-        // Slightly extended bb. Slightly off-centred just so on symmetric
-        // geometry there are less face/edge aligned items.
-        treeBoundBox bb
-        (
-            treeBoundBox(points()).extend(rndGen, 1E-4)
-        );
-
-        bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-        bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-
-        labelListList sliceEdges(nEdgeTypes);
-
-        // External edges
-        sliceEdges[0] =
-            identity(internalStart_ - externalStart_) + externalStart_;
-
-        // Internal edges
-        sliceEdges[1] = identity(flatStart_ - internalStart_) + internalStart_;
-
-        // Flat edges
-        sliceEdges[2] = identity(openStart_ - flatStart_) + flatStart_;
-
-        // Open edges
-        sliceEdges[3] = identity(multipleStart_ - openStart_) + openStart_;
-
-        // Multiple edges
-        sliceEdges[4] =
-            identity(edges().size() - multipleStart_) + multipleStart_;
-
-        forAll(edgeTreesByType_, i)
-        {
-            edgeTreesByType_.set
-            (
-                i,
-                new indexedOctree<treeDataEdge>
-                (
-                    treeDataEdge
-                    (
-                        false,          // cachebb
-                        edges(),        // edges
-                        points(),       // points
-                        sliceEdges[i]   // selected edges
-                    ),
-                    bb,     // bb
-                    8,      // maxLevel
-                    10,     // leafsize
-                    3.0     // duplicity
-                )
-            );
-        }
-    }
-
-    return edgeTreesByType_;
-}
-
-
-void Foam::featureEdgeMesh::writeObj
-(
-    const fileName& prefix
-) const
-{
-    Pout<< nl << "Writing featureEdgeMesh components to " << prefix << endl;
-
-    label verti = 0;
-
-    edgeMesh::write(prefix + "_edgeMesh.obj");
-
-    OFstream convexFtPtStr(prefix + "_convexFeaturePts.obj");
-    Pout<< "Writing convex feature points to " << convexFtPtStr.name() << endl;
-
-    for(label i = 0; i < concaveStart_; i++)
-    {
-        meshTools::writeOBJ(convexFtPtStr, points()[i]);
-    }
-
-    OFstream concaveFtPtStr(prefix + "_concaveFeaturePts.obj");
-    Pout<< "Writing concave feature points to "
-        << concaveFtPtStr.name() << endl;
-
-    for(label i = concaveStart_; i < mixedStart_; i++)
-    {
-        meshTools::writeOBJ(concaveFtPtStr, points()[i]);
-    }
-
-    OFstream mixedFtPtStr(prefix + "_mixedFeaturePts.obj");
-    Pout<< "Writing mixed feature points to " << mixedFtPtStr.name() << endl;
-
-    for(label i = mixedStart_; i < nonFeatureStart_; i++)
-    {
-        meshTools::writeOBJ(mixedFtPtStr, points()[i]);
-    }
-
-    OFstream mixedFtPtStructureStr(prefix + "_mixedFeaturePtsStructure.obj");
-    Pout<< "Writing mixed feature point structure to "
-        << mixedFtPtStructureStr.name() << endl;
-
-    verti = 0;
-    for(label i = mixedStart_; i < nonFeatureStart_; i++)
-    {
-        const labelList& ptEds = pointEdges()[i];
-
-        forAll(ptEds, j)
-        {
-            const edge& e = edges()[ptEds[j]];
-
-            meshTools::writeOBJ(mixedFtPtStructureStr, points()[e[0]]); verti++;
-            meshTools::writeOBJ(mixedFtPtStructureStr, points()[e[1]]); verti++;
-            mixedFtPtStructureStr << "l " << verti-1 << ' ' << verti << endl;
-        }
-    }
-
-    OFstream externalStr(prefix + "_externalEdges.obj");
-    Pout<< "Writing external edges to " << externalStr.name() << endl;
-
-    verti = 0;
-    for (label i = externalStart_; i < internalStart_; i++)
-    {
-        const edge& e = edges()[i];
-
-        meshTools::writeOBJ(externalStr, points()[e[0]]); verti++;
-        meshTools::writeOBJ(externalStr, points()[e[1]]); verti++;
-        externalStr << "l " << verti-1 << ' ' << verti << endl;
-    }
-
-    OFstream internalStr(prefix + "_internalEdges.obj");
-    Pout<< "Writing internal edges to " << internalStr.name() << endl;
-
-    verti = 0;
-    for (label i = internalStart_; i < flatStart_; i++)
-    {
-        const edge& e = edges()[i];
-
-        meshTools::writeOBJ(internalStr, points()[e[0]]); verti++;
-        meshTools::writeOBJ(internalStr, points()[e[1]]); verti++;
-        internalStr << "l " << verti-1 << ' ' << verti << endl;
-    }
-
-    OFstream flatStr(prefix + "_flatEdges.obj");
-    Pout<< "Writing flat edges to " << flatStr.name() << endl;
-
-    verti = 0;
-    for (label i = flatStart_; i < openStart_; i++)
-    {
-        const edge& e = edges()[i];
-
-        meshTools::writeOBJ(flatStr, points()[e[0]]); verti++;
-        meshTools::writeOBJ(flatStr, points()[e[1]]); verti++;
-        flatStr << "l " << verti-1 << ' ' << verti << endl;
-    }
-
-    OFstream openStr(prefix + "_openEdges.obj");
-    Pout<< "Writing open edges to " << openStr.name() << endl;
-
-    verti = 0;
-    for (label i = openStart_; i < multipleStart_; i++)
-    {
-        const edge& e = edges()[i];
-
-        meshTools::writeOBJ(openStr, points()[e[0]]); verti++;
-        meshTools::writeOBJ(openStr, points()[e[1]]); verti++;
-        openStr << "l " << verti-1 << ' ' << verti << endl;
-    }
-
-    OFstream multipleStr(prefix + "_multipleEdges.obj");
-    Pout<< "Writing multiple edges to " << multipleStr.name() << endl;
-
-    verti = 0;
-    for (label i = multipleStart_; i < edges().size(); i++)
-    {
-        const edge& e = edges()[i];
-
-        meshTools::writeOBJ(multipleStr, points()[e[0]]); verti++;
-        meshTools::writeOBJ(multipleStr, points()[e[1]]); verti++;
-        multipleStr << "l " << verti-1 << ' ' << verti << endl;
-    }
-
-    OFstream regionStr(prefix + "_regionEdges.obj");
-    Pout<< "Writing region edges to " << regionStr.name() << endl;
-
-    verti = 0;
-    forAll(regionEdges_, i)
-    {
-        const edge& e = edges()[regionEdges_[i]];
-
-        meshTools::writeOBJ(regionStr, points()[e[0]]); verti++;
-        meshTools::writeOBJ(regionStr, points()[e[1]]); verti++;
-        regionStr << "l " << verti-1 << ' ' << verti << endl;
-    }
+    is >> *this;
+    return !is.bad();
 }
 
 
 bool Foam::featureEdgeMesh::writeData(Ostream& os) const
 {
-    os  << "// points, edges, concaveStart, mixedStart, nonFeatureStart, " << nl
-        << "// internalStart, flatStart, openStart, multipleStart, " << nl
-        << "// normals, edgeNormals, featurePointNormals, regionEdges" << nl
-        << endl;
-
-    os  << points() << nl
-        << edges() << nl
-        << concaveStart_ << token::SPACE
-        << mixedStart_ << token::SPACE
-        << nonFeatureStart_ << token::SPACE
-        << internalStart_ << token::SPACE
-        << flatStart_ << token::SPACE
-        << openStart_ << token::SPACE
-        << multipleStart_ << nl
-        << normals_ << nl
-        << edgeNormals_ << nl
-        << featurePointNormals_ << nl
-        << regionEdges_
-        << endl;
+    os << *this;
 
     return os.good();
 }
diff --git a/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.H b/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.H
index 4e1e360fa5f9c127940f73ba2e9c031f6cb3e62a..9e7982a7c6f8cd320f1fe56e160afa0c38d820b1 100644
--- a/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.H
+++ b/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,27 +25,12 @@ Class
     Foam::featureEdgeMesh
 
 Description
+    edgeMesh + IO.
 
-    Feature points are a sorted subset at the start of the overall points list:
-        0 .. concaveStart_-1                : convex points (w.r.t normals)
-        concaveStart_-1 .. mixedStart_-1    : concave points
-        mixedStart_ .. nonFeatureStart_     : mixed internal/external points
-        nonFeatureStart_ .. size-1          : non-feature points
-
-    Feature edges are the edgeList of the edgeMesh and are sorted:
-        0 .. internalStart_-1           : external edges (convex w.r.t normals)
-        internalStart_ .. flatStart_-1  : internal edges (concave)
-        flatStart_ .. openStart_-1      : flat edges (neither concave or convex)
-                                          can arise from region interfaces on
-                                          flat surfaces
-        openStart_ .. multipleStart_-1  : open edges (e.g. from baffle surfaces)
-        multipleStart_ .. size-1        : multiply connected edges
-
-    The edge direction and feature edge and feature point adjacent normals
-    are stored.
+    See also extendedFeatureEdgeMesh type which stores additional classification
+    of features.
 
 SourceFiles
-    featureEdgeMeshI.H
     featureEdgeMesh.C
 
 \*---------------------------------------------------------------------------*/
@@ -54,12 +39,7 @@ SourceFiles
 #define featureEdgeMesh_H
 
 #include "edgeMesh.H"
-#include "surfaceFeatures.H"
-#include "objectRegistry.H"
-#include "IOdictionary.H"
-#include "indexedOctree.H"
-#include "treeDataEdge.H"
-#include "pointIndexHit.H"
+#include "regIOobject.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -67,7 +47,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                       Class featureEdgeMesh Declaration
+                           Class featureEdgeMesh Declaration
 \*---------------------------------------------------------------------------*/
 
 class featureEdgeMesh
@@ -78,281 +58,31 @@ class featureEdgeMesh
 
 public:
 
-    //- Runtime type information
     TypeName("featureEdgeMesh");
 
-    enum pointStatus
-    {
-        CONVEX,     // Fully convex point (w.r.t normals)
-        CONCAVE,    // Fully concave point
-        MIXED,      // A point surrounded by both convex and concave edges
-        NONFEATURE  // Not a feature point
-    };
-
-    enum edgeStatus
-    {
-        EXTERNAL,   // "Convex" edge
-        INTERNAL,   // "Concave" edge
-        FLAT,       // Neither concave or convex, on a flat surface
-        OPEN,       // i.e. only connected to one face
-        MULTIPLE,   // Multiply connected (connected to more than two faces)
-        NONE        // Not a classified feature edge (consistency with
-                    // surfaceFeatures)
-    };
-
-private:
-
-    // Static data
-
-        //- Angular closeness tolerance for treating normals as the same
-        static scalar cosNormalAngleTol_;
-
-        //- Index of the start of the convex feature points - static as 0
-        static label convexStart_;
-
-        //- Index of the start of the external feature edges - static as 0
-        static label externalStart_;
-
-
-    // Private data
-
-        //- Index of the start of the concave feature points
-        label concaveStart_;
-
-        //- Index of the start of the mixed type feature points
-        label mixedStart_;
-
-        //- Index of the start of the non-feature points
-        label nonFeatureStart_;
-
-        //- Index of the start of the internal feature edges
-        label internalStart_;
-
-        //- Index of the start of the flat feature edges
-        label flatStart_;
-
-        //- Index of the start of the open feature edges
-        label openStart_;
-
-        //- Index of the start of the multiply-connected feature edges
-        label multipleStart_;
-
-        //- Normals of the features, to be referred to by index by both feature
-        //  points and edges, unsorted
-        vectorField normals_;
-
-        //- Flat and open edges require the direction of the edge
-        vectorField edgeDirections_;
-
-        //- Indices of the normals that are adjacent to the feature edges
-        labelListList edgeNormals_;
-
-        //- Indices of the normals that are adjacent to the feature points
-        labelListList featurePointNormals_;
-
-        //- Feature edges which are on the boundary between regions
-        labelList regionEdges_;
-
-        //- Search tree for all edges
-        mutable autoPtr<indexedOctree<treeDataEdge> > edgeTree_;
-
-        //- Individual search trees for each type of edge
-        mutable PtrList<indexedOctree<treeDataEdge> > edgeTreesByType_;
-
-
-    // Private Member Functions
-
-        //- Classify the type of feature point.  Requires valid stored member
-        //  data for edges and normals.
-        pointStatus classifyFeaturePoint(label ptI) const;
-
-        //- Classify the type of feature edge.  Requires face centre 0 to face
-        //  centre 1 vector to distinguish internal from external
-        edgeStatus classifyEdge
-        (
-            const List<vector>& norms,
-            const labelList& edNorms,
-            const vector& fC0tofC1
-        ) const;
-
-
-public:
-
-    // Static data
-
-        //- Number of possible point types (i.e. number of slices)
-        static label nPointTypes;
-
-        //- Number of possible feature edge types (i.e. number of slices)
-        static label nEdgeTypes;
 
     // Constructors
 
         //- Construct (read) given an IOobject
         featureEdgeMesh(const IOobject&);
 
-        //- Construct as copy
-        explicit featureEdgeMesh(const IOobject&, const featureEdgeMesh&);
-
-        //- Construct by transferring components (points, edges)
+        //- Construct from featureEdgeMesh data
         featureEdgeMesh
         (
             const IOobject&,
-            const Xfer<pointField>&,
-            const Xfer<edgeList>&
-        );
-
-        //- Construct (read) given surfaceFeatures, an objectRegistry and a
-        //  fileName to write to.  Extracts, classifies and reorders the data
-        //  from surfaceFeatures.
-        featureEdgeMesh
-        (
-            const surfaceFeatures& sFeat,
-            const objectRegistry& obr,
-            const fileName& sFeatFileName
+            const pointField&,
+            const edgeList&
         );
 
-        //- Construct from all components
-        featureEdgeMesh
-        (
-            const IOobject& io,
-            const pointField& pts,
-            const edgeList& eds,
-            label concaveStart,
-            label mixedStart,
-            label nonFeatureStart,
-            label internalStart,
-            label flatStart,
-            label openStart,
-            label multipleStart,
-            const vectorField& normals,
-            const vectorField& edgeDirections,
-            const labelListList& edgeNormals,
-            const labelListList& featurePointNormals,
-            const labelList& regionEdges
-        );
-
-
-    //- Destructor
-    ~featureEdgeMesh();
-
-
-    // Member Functions
-
-        // Find
-
-            //- Find nearest surface edge for the sample point.
-            void nearestFeatureEdge
-            (
-                const point& sample,
-                scalar searchDistSqr,
-                pointIndexHit& info
-            ) const;
-
-            //- Find nearest surface edge for each sample point.
-            void nearestFeatureEdge
-            (
-                const pointField& samples,
-                const scalarField& searchDistSqr,
-                List<pointIndexHit>& info
-            ) const;
-
-            //- Find the nearest point on each type of feature edge
-            void nearestFeatureEdgeByType
-            (
-                const point& sample,
-                const scalarField& searchDistSqr,
-                List<pointIndexHit>& info
-            ) const;
-
-        // Access
-
-            //- Return the index of the start of the convex feature points
-            inline label convexStart() const;
-
-            //- Return the index of the start of the concave feature points
-            inline label concaveStart() const;
-
-            //- Return the index of the start of the mixed type feature points
-            inline label mixedStart() const;
-
-            //- Return the index of the start of the non-feature points
-            inline label nonFeatureStart() const;
-
-            //- Return the index of the start of the external feature edges
-            inline label externalStart() const;
-
-            //- Return the index of the start of the internal feature edges
-            inline label internalStart() const;
-
-            //- Return the index of the start of the flat feature edges
-            inline label flatStart() const;
-
-            //- Return the index of the start of the open feature edges
-            inline label openStart() const;
-
-            //- Return the index of the start of the multiply-connected feature
-            //  edges
-            inline label multipleStart() const;
-
-            //- Return whether or not the point index is a feature point
-            inline bool featurePoint(label ptI) const;
-
-            //- Return the normals of the surfaces adjacent to the feature edges
-            //  and points
-            inline const vectorField& normals() const;
-
-            //- Return the edgeDirection vectors
-            inline const vectorField& edgeDirections() const;
-
-            //- Return the direction of edgeI, pointing away from ptI
-            inline vector edgeDirection(label edgeI, label ptI) const;
-
-            //- Return the indices of the normals that are adjacent to the
-            //  feature edges
-            inline const labelListList& edgeNormals() const;
-
-            //- Return the normal vectors for a given set of normal indices
-            inline vectorField edgeNormals(const labelList& edgeNormIs) const;
-
-            //- Return the normal vectors for a given edge
-            inline vectorField edgeNormals(label edgeI) const;
-
-            //- Return the indices of the normals that are adjacent to the
-            //  feature points
-            inline const labelListList& featurePointNormals() const;
-
-            //- Return the normal vectors for a given feature point
-            inline vectorField featurePointNormals(label ptI) const;
-
-            //- Return the feature edges which are on the boundary between
-            //  regions
-            inline const labelList& regionEdges() const;
-
-            //- Return the pointStatus of a specified point
-            inline pointStatus getPointStatus(label ptI) const;
-
-            //- Return the edgeStatus of a specified edge
-            inline edgeStatus getEdgeStatus(label edgeI) const;
-
-            //- Demand driven construction of octree for boundary edges
-            const indexedOctree<treeDataEdge>& edgeTree() const;
-
-            //- Demand driven construction of octree for boundary edges by type
-            const PtrList<indexedOctree<treeDataEdge> >&
-            edgeTreesByType() const;
-
-
-        // Write
+        //- Construct as copy
+        featureEdgeMesh(const IOobject&, const featureEdgeMesh&);
 
-            //- Write all components of the featureEdgeMesh as obj files
-            void writeObj(const fileName& prefix) const;
 
-            //- Give precedence to the regIOobject write
-            using regIOobject::write;
+        //- ReadData function required for regIOobject read operation
+        virtual bool readData(Istream&);
 
-            //- WriteData function required for regIOobject write operation
-            virtual bool writeData(Ostream&) const;
+        //- WriteData function required for regIOobject write operation
+        virtual bool writeData(Ostream&) const;
 };
 
 
@@ -362,10 +92,6 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#include "featureEdgeMeshI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #endif
 
 // ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/adjustPhi/adjustPhi.C b/src/finiteVolume/cfdTools/general/adjustPhi/adjustPhi.C
index 1bfeda1018580efb4c4b1e7565f41899f6a9b3f0..93a23c94cb13c12ae35183bed6a7c92fc7b213df 100644
--- a/src/finiteVolume/cfdTools/general/adjustPhi/adjustPhi.C
+++ b/src/finiteVolume/cfdTools/general/adjustPhi/adjustPhi.C
@@ -40,7 +40,8 @@ bool Foam::adjustPhi
 {
     if (p.needReference())
     {
-        p.boundaryField().updateCoeffs();
+        // p coefficients should not be updated here
+        // p.boundaryField().updateCoeffs();
 
         scalar massIn = 0.0;
         scalar fixedMassOut = 0.0;
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
index 46d4ae134313aeea7bca212527f51394aa7dca98..164da09c5ba883db5d31735b984da0db95173908 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -45,6 +45,7 @@ Description
 #include "layerParameters.H"
 #include "combineFaces.H"
 #include "IOmanip.H"
+#include "globalIndex.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -888,68 +889,6 @@ void Foam::autoLayerDriver::handleNonManifolds
     Info<< "Outside of local patch is multiply connected across edges or"
         << " points at " << nNonManif << " points." << endl;
 
-
-    const labelList& meshPoints = pp.meshPoints();
-
-
-    // 2. Parallel check
-    // For now disable extrusion at any shared edges.
-    const labelHashSet sharedEdgeSet(mesh.globalData().sharedEdgeLabels());
-
-    forAll(pp.edges(), edgeI)
-    {
-        if (sharedEdgeSet.found(meshEdges[edgeI]))
-        {
-            const edge& e = mesh.edges()[meshEdges[edgeI]];
-
-            Pout<< "Disabling extrusion at edge "
-                << mesh.points()[e[0]] << mesh.points()[e[1]]
-                << " since it is non-manifold across coupled patches."
-                << endl;
-
-            nonManifoldPoints.insert(e[0]);
-            nonManifoldPoints.insert(e[1]);
-        }
-    }
-
-    // 3b. extrusion can produce multiple faces between 2 cells
-    // across processor boundary
-    // This occurs when a coupled face shares more than 1 edge with a
-    // non-coupled boundary face.
-    // This is now correctly handled by addPatchCellLayer in that it
-    // extrudes a single face from the stringed up edges.
-
-
-    nNonManif = returnReduce(nonManifoldPoints.size(), sumOp<label>());
-
-    if (nNonManif > 0)
-    {
-        Info<< "Outside of patches is multiply connected across edges or"
-            << " points at " << nNonManif << " points." << nl
-            << "Writing " << nNonManif
-            << " points where this happens to pointSet "
-            << nonManifoldPoints.name() << nl
-            << "and setting layer thickness to zero on these points."
-            << endl;
-
-        nonManifoldPoints.instance() = mesh.time().timeName();
-        nonManifoldPoints.write();
-
-        forAll(meshPoints, patchPointI)
-        {
-            if (nonManifoldPoints.found(meshPoints[patchPointI]))
-            {
-                unmarkExtrusion
-                (
-                    patchPointI,
-                    patchDisp,
-                    patchNLayers,
-                    extrudeStatus
-                );
-            }
-        }
-    }
-
     Info<< "Set displacement to zero for all " << nNonManif
         << " non-manifold points" << endl;
 }
@@ -1348,7 +1287,6 @@ void Foam::autoLayerDriver::setNumLayers
 }
 
 
-// Grow no-extrusion layer
 void Foam::autoLayerDriver::growNoExtrusion
 (
     const indirectPrimitivePatch& pp,
@@ -1439,6 +1377,103 @@ void Foam::autoLayerDriver::growNoExtrusion
 }
 
 
+void Foam::autoLayerDriver::determineSidePatches
+(
+    const globalIndex& globalFaces,
+    const labelListList& edgeGlobalFaces,
+    const indirectPrimitivePatch& pp,
+
+    labelList& sidePatchID
+)
+{
+    // Sometimes edges-to-be-extruded are on more than 2 processors.
+    // Work out which 2 hold the faces to be extruded and thus which procpatch
+    // the side-face should be in. As an additional complication this might
+    // mean that 2 procesors that were only edge-connected now suddenly need
+    // to become face-connected i.e. have a processor patch between them.
+
+    fvMesh& mesh = meshRefiner_.mesh();
+
+    // Determine sidePatchID. Any additional processor boundary gets added to
+    // patchToNbrProc,nbrProcToPatch and nPatches gets set to the new number
+    // of patches.
+    label nPatches;
+    Map<label> nbrProcToPatch;
+    Map<label> patchToNbrProc;
+    addPatchCellLayer::calcSidePatch
+    (
+        mesh,
+        globalFaces,
+        edgeGlobalFaces,
+        pp,
+
+        sidePatchID,
+        nPatches,
+        nbrProcToPatch,
+        patchToNbrProc
+    );
+
+    label nOldPatches = mesh.boundaryMesh().size();
+    label nAdded = returnReduce(nPatches-nOldPatches, sumOp<label>());
+    Info<< nl << "Adding in total " << nAdded/2 << " inter-processor patches to"
+        << " handle extrusion of non-manifold processor boundaries."
+        << endl;
+
+    if (nAdded > 0)
+    {
+        // We might not add patches in same order as in patchToNbrProc
+        // so prepare to renumber sidePatchID
+        Map<label> wantedToAddedPatch;
+
+        for (label patchI = nOldPatches; patchI < nPatches; patchI++)
+        {
+            label nbrProcI = patchToNbrProc[patchI];
+            word name =
+                    "procBoundary"
+                  + Foam::name(Pstream::myProcNo())
+                  + "to"
+                  + Foam::name(nbrProcI);
+
+            dictionary patchDict;
+            patchDict.add("type", processorPolyPatch::typeName);
+            patchDict.add("myProcNo", Pstream::myProcNo());
+            patchDict.add("neighbProcNo", nbrProcI);
+            patchDict.add("nFaces", 0);
+            patchDict.add("startFace", mesh.nFaces());
+
+            Pout<< "Adding patch " << patchI
+                << " name:" << name
+                << " between " << Pstream::myProcNo()
+                << " and " << nbrProcI
+                << endl;
+
+            label procPatchI = meshRefiner_.appendPatch
+            (
+                mesh,
+                mesh.boundaryMesh().size(), // new patch index
+                name,
+                patchDict
+            );
+            wantedToAddedPatch.insert(patchI, procPatchI);
+        }
+
+        // Renumber sidePatchID
+        forAll(sidePatchID, i)
+        {
+            label patchI = sidePatchID[i];
+            Map<label>::const_iterator fnd = wantedToAddedPatch.find(patchI);
+            if (fnd != wantedToAddedPatch.end())
+            {
+                sidePatchID[i] = fnd();
+            }
+        }
+
+        mesh.clearOut();
+        const_cast<polyBoundaryMesh&>(mesh.boundaryMesh()).updateMesh();
+    }
+}
+
+
 void Foam::autoLayerDriver::calculateLayerThickness
 (
     const indirectPrimitivePatch& pp,
@@ -2553,6 +2588,37 @@ void Foam::autoLayerDriver::addLayers
         )
     );
 
+
+    // Global face indices engine
+    const globalIndex globalFaces(mesh.nFaces());
+
+    // Determine extrudePatch.edgeFaces in global numbering (so across
+    // coupled patches). This is used only to string up edges between coupled
+    // faces (all edges between same (global)face indices get extruded).
+    labelListList edgeGlobalFaces
+    (
+        addPatchCellLayer::globalEdgeFaces
+        (
+            mesh,
+            globalFaces,
+            pp
+        )
+    );
+
+    // Determine patches for extruded boundary edges. Calculates if any
+    // additional processor patches need to be constructed.
+
+    labelList sidePatchID;
+    determineSidePatches
+    (
+        globalFaces,
+        edgeGlobalFaces,
+        pp,
+
+        sidePatchID
+    );
+
+
     // Construct iterative mesh mover.
     Info<< "Constructing mesh displacer ..." << endl;
 
@@ -3038,8 +3104,12 @@ void Foam::autoLayerDriver::addLayers
         // Not add layer if patchDisp is zero.
         addLayer.setRefinement
         (
+            globalFaces,
+            edgeGlobalFaces,
+
             invExpansionRatio,
             pp(),
+            sidePatchID,        // boundary patch for extruded boundary edges
             labelList(0),       // exposed patchIDs, not used for adding layers
             nPatchFaceLayers,   // layers per face
             nPatchPointLayers,  // layers per point
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
index c96e304983ead0e5d07a31a06141a386e85bbe4f..fd5a1f3aa881b53130ad551a6ae93dc1ee1756e5 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -243,6 +243,16 @@ class autoLayerDriver
                     List<extrudeMode>& extrudeStatus
                 ) const;
 
+                //- See what patches boundaryedges should be extruded into
+                void determineSidePatches
+                (
+                    const globalIndex& globalFaces,
+                    const labelListList& edgeGlobalFaces,
+                    const indirectPrimitivePatch& pp,
+
+                    labelList& sidePatchID
+                );
+
                 //- Calculate pointwise wanted and minimum thickness.
                 //  thickness: wanted thickness
                 //  minthickness: when to give up and not extrude
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
index 8a749cd60b6369d7a298c6d6df8f1d0ab779a18d..de0ba875879649100d882eb19b51a17b16695b34 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
@@ -1549,60 +1549,28 @@ void Foam::meshRefinement::checkCoupledFaceZones(const polyMesh& mesh)
 }
 
 
-// Adds patch if not yet there. Returns patchID.
-Foam::label Foam::meshRefinement::addPatch
+Foam::label Foam::meshRefinement::appendPatch
 (
     fvMesh& mesh,
+    const label insertPatchI,
     const word& patchName,
-    const dictionary& patchInfo
+    const dictionary& patchDict
 )
 {
-    polyBoundaryMesh& polyPatches =
-        const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
-
-    const label patchI = polyPatches.findPatchID(patchName);
-    if (patchI != -1)
-    {
-        // Already there
-        return patchI;
-    }
-
-
-    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();
-
+    polyBoundaryMesh& polyPatches =
+        const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
     fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
 
-    dictionary patchDict(patchInfo);
-    patchDict.set("nFaces", 0);
-    patchDict.set("startFace", startFaceI);
+    label patchI = polyPatches.size();
 
     // Add polyPatch at the end
-    polyPatches.setSize(sz+1);
+    polyPatches.setSize(patchI+1);
     polyPatches.set
     (
-        sz,
+        patchI,
         polyPatch::New
         (
             patchName,
@@ -1611,13 +1579,13 @@ Foam::label Foam::meshRefinement::addPatch
             polyPatches
         )
     );
-    fvPatches.setSize(sz+1);
+    fvPatches.setSize(patchI+1);
     fvPatches.set
     (
-        sz,
+        patchI,
         fvPatch::New
         (
-            polyPatches[sz],  // point to newly added polyPatch
+            polyPatches[patchI],  // point to newly added polyPatch
             mesh.boundary()
         )
     );
@@ -1675,21 +1643,69 @@ Foam::label Foam::meshRefinement::addPatch
         mesh,
         calculatedFvPatchField<tensor>::typeName
     );
+    return patchI;
+}
+
+
+// Adds patch if not yet there. Returns patchID.
+Foam::label Foam::meshRefinement::addPatch
+(
+    fvMesh& mesh,
+    const word& patchName,
+    const dictionary& patchInfo
+)
+{
+    polyBoundaryMesh& polyPatches =
+        const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
+    fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
+
+    const label patchI = polyPatches.findPatchID(patchName);
+    if (patchI != -1)
+    {
+        // Already there
+        return patchI;
+    }
+
+
+    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(patchInfo);
+    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.
+
+    label addedPatchI = appendPatch(mesh, insertPatchI, patchName, patchDict);
+
 
     // Create reordering list
     // patches before insert position stay as is
-    labelList oldToNew(sz+1);
+    labelList oldToNew(addedPatchI+1);
     for (label i = 0; i < insertPatchI; i++)
     {
         oldToNew[i] = i;
     }
     // patches after insert position move one up
-    for (label i = insertPatchI; i < sz; i++)
+    for (label i = insertPatchI; i < addedPatchI; i++)
     {
         oldToNew[i] = i+1;
     }
     // appended patch gets moved to insert position
-    oldToNew[sz] = insertPatchI;
+    oldToNew[addedPatchI] = insertPatchI;
 
     // Shuffle into place
     polyPatches.reorder(oldToNew);
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index 7337161bb750dbe941bef9bf87aa6bbe208b1d0b..5a1d8ff821194ff719021b7ce5d32834ab1250ca 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -747,8 +747,17 @@ public:
 
         // Other topo changes
 
+            //- Helper:append patch to end of mesh.
+            static label appendPatch
+            (
+                fvMesh&,
+                const label insertPatchI,
+                const word&,
+                const dictionary&
+            );
+
             //- Helper:add patch to mesh. Update all registered fields.
-            //  Use addMeshedPatch to add patches originating from surfaces.
+            //  Used by addMeshedPatch to add patches originating from surfaces.
             static label addPatch(fvMesh&, const word& name, const dictionary&);
 
             //- Add patch originating from meshing. Update meshedPatches_.
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
index 18e5bf01b46c5fcf98957f628394df0523b3b6c6..b9e2d4d3e077e0eec1eccf78296ca3b06dd24a31 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
+++ b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -1007,6 +1007,82 @@ void Foam::refinementSurfaces::findNearestRegion
 }
 
 
+void Foam::refinementSurfaces::findNearestRegion
+(
+    const labelList& surfacesToTest,
+    const pointField& samples,
+    const scalarField& nearestDistSqr,
+    labelList& hitSurface,
+    List<pointIndexHit>& hitInfo,
+    labelList& hitRegion,
+    vectorField& hitNormal
+) const
+{
+    labelList geometries(UIndirectList<label>(surfaces_, surfacesToTest));
+
+    // Do the tests. Note that findNearest returns index in geometries.
+    searchableSurfacesQueries::findNearest
+    (
+        allGeometry_,
+        geometries,
+        samples,
+        nearestDistSqr,
+        hitSurface,
+        hitInfo
+    );
+
+    // Rework the hitSurface to be surface (i.e. index into surfaces_)
+    forAll(hitSurface, pointI)
+    {
+        if (hitSurface[pointI] != -1)
+        {
+            hitSurface[pointI] = surfacesToTest[hitSurface[pointI]];
+        }
+    }
+
+    // Collect the region
+    hitRegion.setSize(hitSurface.size());
+    hitRegion = -1;
+    hitNormal.setSize(hitSurface.size());
+    hitNormal = vector::zero;
+
+    forAll(surfacesToTest, i)
+    {
+        label surfI = surfacesToTest[i];
+
+        // Collect hits for surfI
+        const labelList localIndices(findIndices(hitSurface, surfI));
+
+        List<pointIndexHit> localHits
+        (
+            UIndirectList<pointIndexHit>
+            (
+                hitInfo,
+                localIndices
+            )
+        );
+
+        // Region
+        labelList localRegion;
+        allGeometry_[surfaces_[surfI]].getRegion(localHits, localRegion);
+
+        forAll(localIndices, i)
+        {
+            hitRegion[localIndices[i]] = localRegion[i];
+        }
+
+        // Normal
+        vectorField localNormal;
+        allGeometry_[surfaces_[surfI]].getNormal(localHits, localNormal);
+
+        forAll(localIndices, i)
+        {
+            hitNormal[localIndices[i]] = localNormal[i];
+        }
+    }
+}
+
+
 //// Find intersection with max of edge. Return -1 or the surface
 //// with the highest maxLevel above currentLevel
 //Foam::label Foam::refinementSurfaces::findHighestIntersection
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
index dc243f316fb0ece4db4443cc624c18a5d72cac66..a7728b90824eca7117541159f8e7f1c5587df046 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
+++ b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -327,6 +327,19 @@ public:
                 labelList& hitRegion
             ) const;
 
+            //- Find nearest point on surfaces. Return surface, region and
+            //  normal on surface (so not global surface)
+            void findNearestRegion
+            (
+                const labelList& surfacesToTest,
+                const pointField& samples,
+                const scalarField& nearestDistSqr,
+                labelList& hitSurface,
+                List<pointIndexHit>& hitInfo,
+                labelList& hitRegion,
+                vectorField& hitNormal
+            ) const;
+
             //- Detect if a point is 'inside' (closed) surfaces.
             //  Returns -1 if not, returns first surface it is.
             void findInside
diff --git a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
index 291530963909ba870f0c58f2a18de3714abf3773..1cfb22887b7c795b97956e4b1336a75eaab9a84c 100644
--- a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
+++ b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
@@ -398,7 +398,15 @@ void Foam::decompositionMethod::calcCellCells
     // Create global cell numbers
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    label nAgglom = max(agglom)+1;
+    label nAgglom;
+    if (agglom.size())
+    {
+        nAgglom = max(agglom)+1;
+    }
+    else
+    {
+        nAgglom = 0;
+    }
     globalIndex globalAgglom(nAgglom);
 
     const labelList& faceOwner = mesh.faceOwner();
diff --git a/src/parallel/decompose/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.C b/src/parallel/decompose/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.C
index 9f5ff83436a0337cbc865a7eb8f499a77b33865d..4ce8706dea087ca9abea6eea014f425def63ce77 100644
--- a/src/parallel/decompose/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.C
+++ b/src/parallel/decompose/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,6 +26,7 @@ License
 #include "simpleGeomDecomp.H"
 #include "addToRunTimeSelectionTable.H"
 #include "SortableList.H"
+#include "globalIndex.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -53,7 +54,7 @@ void Foam::simpleGeomDecomp::assignToProcessorGroup
 (
     labelList& processorGroup,
     const label nProcGroup
-)
+) const
 {
     label jump = processorGroup.size()/nProcGroup;
     label jumpb = jump + 1;
@@ -90,7 +91,7 @@ void Foam::simpleGeomDecomp::assignToProcessorGroup
     const labelList& indices,
     const scalarField& weights,
     const scalar summedWeights
-)
+) const
 {
     // This routine gets the sorted points.
     // Easiest to explain with an example.
@@ -126,7 +127,10 @@ void Foam::simpleGeomDecomp::assignToProcessorGroup
 }
 
 
-Foam::labelList Foam::simpleGeomDecomp::decompose(const pointField& points)
+Foam::labelList Foam::simpleGeomDecomp::decomposeOneProc
+(
+    const pointField& points
+) const
 {
     // construct a list for the final result
     labelList finalDecomp(points.size());
@@ -195,11 +199,11 @@ Foam::labelList Foam::simpleGeomDecomp::decompose(const pointField& points)
 }
 
 
-Foam::labelList Foam::simpleGeomDecomp::decompose
+Foam::labelList Foam::simpleGeomDecomp::decomposeOneProc
 (
     const pointField& points,
     const scalarField& weights
-)
+) const
 {
     // construct a list for the final result
     labelList finalDecomp(points.size());
@@ -300,4 +304,162 @@ Foam::simpleGeomDecomp::simpleGeomDecomp(const dictionary& decompositionDict)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+Foam::labelList Foam::simpleGeomDecomp::decompose
+(
+    const pointField& points
+)
+{
+    if (!Pstream::parRun())
+    {
+        return decomposeOneProc(points);
+    }
+    else
+    {
+        globalIndex globalNumbers(points.size());
+
+        // Collect all points on master
+        if (Pstream::master())
+        {
+            pointField allPoints(globalNumbers.size());
+
+            label nTotalPoints = 0;
+            // Master first
+            SubField<point>(allPoints, points.size()).assign(points);
+            nTotalPoints += points.size();
+
+            // Add slaves
+            for (int slave=1; slave<Pstream::nProcs(); slave++)
+            {
+                IPstream fromSlave(Pstream::scheduled, slave);
+                pointField nbrPoints(fromSlave);
+                SubField<point>
+                (
+                    allPoints,
+                    nbrPoints.size(),
+                    nTotalPoints
+                ).assign(nbrPoints);
+                nTotalPoints += nbrPoints.size();
+            }
+
+            // Decompose
+            labelList finalDecomp(decomposeOneProc(allPoints));
+
+            // Send back
+            for (int slave=1; slave<Pstream::nProcs(); slave++)
+            {
+                OPstream toSlave(Pstream::scheduled, slave);
+                toSlave << SubField<label>
+                (
+                    finalDecomp,
+                    globalNumbers.localSize(slave),
+                    globalNumbers.offset(slave)
+                );
+            }
+            // Get my own part
+            finalDecomp.setSize(points.size());
+
+            return finalDecomp;
+        }
+        else
+        {
+            // Send my points
+            {
+                OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
+                toMaster<< points;
+            }
+
+            // Receive back decomposition
+            IPstream fromMaster(Pstream::scheduled, Pstream::masterNo());
+            labelList finalDecomp(fromMaster);
+
+            return finalDecomp;
+        }
+    }
+}
+
+
+Foam::labelList Foam::simpleGeomDecomp::decompose
+(
+    const pointField& points,
+    const scalarField& weights
+)
+{
+    if (!Pstream::parRun())
+    {
+        return decomposeOneProc(points, weights);
+    }
+    else
+    {
+        globalIndex globalNumbers(points.size());
+
+        // Collect all points on master
+        if (Pstream::master())
+        {
+            pointField allPoints(globalNumbers.size());
+            scalarField allWeights(allPoints.size());
+
+            label nTotalPoints = 0;
+            // Master first
+            SubField<point>(allPoints, points.size()).assign(points);
+            SubField<scalar>(allWeights, points.size()).assign(weights);
+            nTotalPoints += points.size();
+
+            // Add slaves
+            for (int slave=1; slave<Pstream::nProcs(); slave++)
+            {
+                IPstream fromSlave(Pstream::scheduled, slave);
+                pointField nbrPoints(fromSlave);
+                scalarField nbrWeights(fromSlave);
+                SubField<point>
+                (
+                    allPoints,
+                    nbrPoints.size(),
+                    nTotalPoints
+                ).assign(nbrPoints);
+                SubField<scalar>
+                (
+                    allWeights,
+                    nbrWeights.size(),
+                    nTotalPoints
+                ).assign(nbrWeights);
+                nTotalPoints += nbrPoints.size();
+            }
+
+            // Decompose
+            labelList finalDecomp(decomposeOneProc(allPoints, allWeights));
+
+            // Send back
+            for (int slave=1; slave<Pstream::nProcs(); slave++)
+            {
+                OPstream toSlave(Pstream::scheduled, slave);
+                toSlave << SubField<label>
+                (
+                    finalDecomp,
+                    globalNumbers.localSize(slave),
+                    globalNumbers.offset(slave)
+                );
+            }
+            // Get my own part
+            finalDecomp.setSize(points.size());
+
+            return finalDecomp;
+        }
+        else
+        {
+            // Send my points
+            {
+                OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
+                toMaster<< points << weights;
+            }
+
+            // Receive back decomposition
+            IPstream fromMaster(Pstream::scheduled, Pstream::masterNo());
+            labelList finalDecomp(fromMaster);
+
+            return finalDecomp;
+        }
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/src/parallel/decompose/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.H b/src/parallel/decompose/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.H
index d14b05abdee2be99964cd1995f4bb7b3a0ca63af..c829f4f9bf2ab71687ee0e754c392afd000909d8 100644
--- a/src/parallel/decompose/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.H
+++ b/src/parallel/decompose/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -49,7 +49,7 @@ class simpleGeomDecomp
 {
     // Private Member Functions
 
-        void assignToProcessorGroup(labelList& processorGroup, const label);
+        void assignToProcessorGroup(labelList&, const label) const;
 
         void assignToProcessorGroup
         (
@@ -58,7 +58,15 @@ class simpleGeomDecomp
             const labelList& indices,
             const scalarField& weights,
             const scalar summedWeights
-        );
+        ) const;
+
+        labelList decomposeOneProc(const pointField& points) const;
+
+        labelList decomposeOneProc
+        (
+            const pointField& points,
+            const scalarField& weights
+        ) const;
 
         //- Disallow default bitwise copy construct and assignment
         void operator=(const simpleGeomDecomp&);
@@ -86,9 +94,9 @@ public:
 
         virtual bool parallelAware() const
         {
-            // simpleDecomp does not attempt to do anything across proc
-            // boundaries
-            return false;
+            // simpleDecomp sends all points to the master which does
+            // the decomposition.
+            return true;
         }
 
         virtual labelList decompose(const pointField&);
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/Allrun b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/Allrun
index cbb64d4d0752a37d5c46112a800b437d8e681d8e..5929021b719cc35149559c82cf346271c22c6128 100755
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/Allrun
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/Allrun
@@ -16,11 +16,16 @@ setSet -batch system/sets.setSet > log.setSet1 2>&1
 setsToZones -noFlipMap > log.setsToZones 2>&1
 
 # create the first cyclic - lhs of porous zone
-unset FOAM_SETNAN
+# Note that we don't know what value to give these patches-out-of-nothing so
+# - use binary writing to avoid 'nan'
+# - use setFields to set values
+unset FOAM_SIGFPE
 createBaffles cycLeft '(cycLeft_half0 cycLeft_half1)' -overwrite > log.createBaffles1 2>&1
 
 # create the second cyclic - rhs of porous zone
 createBaffles cycRight '(cycRight_half0 cycRight_half1)' -overwrite > log.createBaffles2 2>&1
+# Initialise newly created patchFields to 0
+changeDictionary
 
 runApplication $application
 
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/changeDictionaryDict b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/changeDictionaryDict
new file mode 100644
index 0000000000000000000000000000000000000000..0f9c594877873d95345bf97aa2cdf6baac37f583
--- /dev/null
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/changeDictionaryDict
@@ -0,0 +1,197 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.7.1                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      changeDictionaryDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dictionaryReplacement
+{
+    alphat
+    {
+        boundaryField
+        {
+            cycLeft
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+            cycRight
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+        }
+    }
+    epsilon
+    {
+        boundaryField
+        {
+            cycLeft
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+            cycRight
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+        }
+    }
+    G
+    {
+        boundaryField
+        {
+            cycLeft
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+            cycRight
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+        }
+    }
+    H2O
+    {
+        boundaryField
+        {
+            cycLeft
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+            cycRight
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+        }
+    }
+    k
+    {
+        boundaryField
+        {
+            cycLeft
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+            cycRight
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+        }
+    }
+    mut
+    {
+        boundaryField
+        {
+            cycLeft
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+            cycRight
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+        }
+    }
+    N2
+    {
+        boundaryField
+        {
+            cycLeft
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+            cycRight
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+        }
+    }
+    O2
+    {
+        boundaryField
+        {
+            cycLeft
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+            cycRight
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+        }
+    }
+    p
+    {
+        boundaryField
+        {
+            cycLeft
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+            cycRight
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+        }
+    }
+    T
+    {
+        boundaryField
+        {
+            cycLeft
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+            cycRight
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+        }
+    }
+    U
+    {
+        boundaryField
+        {
+            cycLeft
+            {
+                type            cyclic;
+                value           uniform (0 0 0);
+            }
+            cycRight
+            {
+                type            cyclic;
+                value           uniform (0 0 0);
+            }
+        }
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/controlDict b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/controlDict
index 6190c66a0722fedd86b6b711d2df393480c11013..a1620b8346b19abc53af89248b5e818586d893a0 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/controlDict
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/controlDict
@@ -33,7 +33,7 @@ writeInterval   0.1;
 
 purgeWrite      0;
 
-writeFormat     ascii;
+writeFormat     binary;
 
 writePrecision  10;
 
diff --git a/wmake/wmakePrintBuild b/wmake/wmakePrintBuild
index 103e65cb07ee03d2b761f47570f36b96c2dbcd88..71ff327fe377f04da9bec8294b38c75aee682682 100755
--- a/wmake/wmakePrintBuild
+++ b/wmake/wmakePrintBuild
@@ -98,14 +98,20 @@ done
 #
 build="$WM_PROJECT_DIR/.build"
 
+unset oldVersion oldPackage
+
+#
 # retrieve old values from the $WM_PROJECT_DIR/.build cache, stored as
 #     version [packager]
-set -- $(tail -1 $build 2>/dev/null)
-
-oldVersion="$1"; shift
-oldPackage="$@"
-[ "${oldPackage:-none}" = none ] && unset oldPackage
-
+#
+getOldValues()
+{
+    set -- $(tail -1 $build 2>/dev/null)
+    oldVersion="$1"
+    [ "$#" -gt 0 ] && shift
+    oldPackage="$@"
+    [ "${oldPackage:-none}" = none ] && unset oldPackage
+}
 
 #
 # printTag - output the build tag
@@ -148,6 +154,9 @@ else
 fi
 
 
+# retrieve old values
+getOldValues
+
 #
 # update persistent build tag if possible
 #