diff --git a/src/edgeMesh/Make/files b/src/edgeMesh/Make/files
index 70c5178cbe..ed4cec5d95 100644
--- a/src/edgeMesh/Make/files
+++ b/src/edgeMesh/Make/files
@@ -1,13 +1,20 @@
-edgeMesh.C
-edgeMeshIO.C
-edgeMeshNew.C
+em = .
 
-edgeMeshFormats = edgeMeshFormats
+$(em)/edgeMesh.C
+$(em)/edgeMeshIO.C
+$(em)/edgeMeshNew.C
+
+
+edgeMeshFormats = $(em)/edgeMeshFormats
 $(edgeMeshFormats)/edgeMeshFormatsCore.C
 
 $(edgeMeshFormats)/edgeMesh/edgeMeshFormat.C
 $(edgeMeshFormats)/edgeMesh/edgeMeshFormatRunTime.C
 
+$(edgeMeshFormats)/extendedFeatureEdgeMesh/extendedFeatureEdgeMeshFormat.C
+$(edgeMeshFormats)/extendedFeatureEdgeMesh/extendedFeatureEdgeMeshFormatRunTime.C
+
+
 $(edgeMeshFormats)/nas/NASedgeFormat.C
 $(edgeMeshFormats)/nas/NASedgeFormatRunTime.C
 
@@ -20,9 +27,22 @@ $(edgeMeshFormats)/starcd/STARCDedgeFormatRunTime.C
 $(edgeMeshFormats)/vtk/VTKedgeFormat.C
 $(edgeMeshFormats)/vtk/VTKedgeFormatRunTime.C
 
-extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C
 
-featureEdgeMesh/featureEdgeMesh.C
+$(em)/featureEdgeMesh/featureEdgeMesh.C
+
+eem = $(em)/extendedEdgeMesh
+
+$(eem)/extendedEdgeMesh.C
+$(eem)/extendedEdgeMeshNew.C
+
+$(eem)/extendedEdgeMeshFormats/obj/OBJextendedEdgeFormat.C
+$(eem)/extendedEdgeMeshFormats/obj/OBJextendedEdgeFormatRunTime.C
+$(eem)/extendedEdgeMeshFormats/extendedEdgeMeshFormat/extendedEdgeMeshFormat.C
+$(eem)/extendedEdgeMeshFormats/extendedEdgeMeshFormat/extendedEdgeMeshFormatRunTime.C
+
+efm = $(eem)/extendedFeatureEdgeMesh
+
+$(efm)/extendedFeatureEdgeMesh.C
 
 
 
diff --git a/src/edgeMesh/Make/options b/src/edgeMesh/Make/options
index 3c1b5258d3..61439e40e2 100644
--- a/src/edgeMesh/Make/options
+++ b/src/edgeMesh/Make/options
@@ -1,9 +1,11 @@
 EXE_INC = \
     -I$(LIB_SRC)/fileFormats/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/triSurface/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
 
 LIB_LIBS = \
     -ltriSurface \
     -lmeshTools \
-    -lfileFormats
+    -lfileFormats \
+    -lsurfMesh
diff --git a/src/edgeMesh/edgeMesh.H b/src/edgeMesh/edgeMesh.H
index c8abb55169..d3397a290b 100644
--- a/src/edgeMesh/edgeMesh.H
+++ b/src/edgeMesh/edgeMesh.H
@@ -27,6 +27,9 @@ Class
 Description
     Points connected by edges.
 
+    Can be read from fileName based on extension. Uses ::New factory method
+    to select the reader and transfer the result.
+
 SourceFiles
     edgeMeshI.H
     edgeMesh.C
@@ -253,7 +256,7 @@ public:
 
     // Write
 
-        void writeStats(Ostream&) const;
+        virtual void writeStats(Ostream&) const;
 
         //- Generic write routine. Chooses writer based on extension.
         virtual void write(const fileName& name) const
diff --git a/src/edgeMesh/edgeMeshFormats/extendedFeatureEdgeMesh/extendedFeatureEdgeMeshFormat.C b/src/edgeMesh/edgeMeshFormats/extendedFeatureEdgeMesh/extendedFeatureEdgeMeshFormat.C
new file mode 100644
index 0000000000..5f20b3016b
--- /dev/null
+++ b/src/edgeMesh/edgeMeshFormats/extendedFeatureEdgeMesh/extendedFeatureEdgeMeshFormat.C
@@ -0,0 +1,70 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "extendedFeatureEdgeMeshFormat.H"
+#include "edgeMeshFormat.H"
+#include "IFstream.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fileFormats::extendedFeatureEdgeMeshFormat::extendedFeatureEdgeMeshFormat
+(
+    const fileName& filename
+)
+{
+    read(filename);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::fileFormats::extendedFeatureEdgeMeshFormat::read
+(
+    const fileName& filename
+)
+{
+    clear();
+
+    IFstream is(filename);
+    if (!is.good())
+    {
+        FatalErrorIn
+        (
+            "fileFormats::extendedFeatureEdgeMeshFormat::read(const fileName&)"
+        )
+            << "Cannot read file " << filename
+            << exit(FatalError);
+    }
+
+    return fileFormats::edgeMeshFormat::read
+    (
+        is,
+        this->storedPoints(),
+        this->storedEdges()
+    );
+}
+
+
+// ************************************************************************* //
diff --git a/src/edgeMesh/edgeMeshFormats/extendedFeatureEdgeMesh/extendedFeatureEdgeMeshFormat.H b/src/edgeMesh/edgeMeshFormats/extendedFeatureEdgeMesh/extendedFeatureEdgeMeshFormat.H
new file mode 100644
index 0000000000..cabc9beec0
--- /dev/null
+++ b/src/edgeMesh/edgeMeshFormats/extendedFeatureEdgeMesh/extendedFeatureEdgeMeshFormat.H
@@ -0,0 +1,106 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::fileFormats::extendedFeatureEdgeMeshFormat
+
+Description
+    Provide a means of reading extendedFeatureEdgeMesh as featureEdgeMesh
+
+SourceFiles
+    extendedFeatureEdgeMeshFormat.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef extendedFeatureEdgeMeshFormat_H
+#define extendedFeatureEdgeMeshFormat_H
+
+#include "edgeMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fileFormats
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class extendedFeatureEdgeMeshFormat Declaration
+\*---------------------------------------------------------------------------*/
+
+class extendedFeatureEdgeMeshFormat
+:
+    public edgeMesh
+{
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        extendedFeatureEdgeMeshFormat(const extendedFeatureEdgeMeshFormat&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const extendedFeatureEdgeMeshFormat&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from file name
+        extendedFeatureEdgeMeshFormat(const fileName&);
+
+
+    // Selectors
+
+        //- Read file and return surface
+        static autoPtr<edgeMesh> New(const fileName& name)
+        {
+            return autoPtr<edgeMesh>
+            (
+                new extendedFeatureEdgeMeshFormat(name)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~extendedFeatureEdgeMeshFormat()
+    {}
+
+
+    // Member Functions
+
+        //- Read from file
+        virtual bool read(const fileName&);
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fileFormats
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/edgeMesh/edgeMeshFormats/extendedFeatureEdgeMesh/extendedFeatureEdgeMeshFormatRunTime.C b/src/edgeMesh/edgeMeshFormats/extendedFeatureEdgeMesh/extendedFeatureEdgeMeshFormatRunTime.C
new file mode 100644
index 0000000000..958f99facb
--- /dev/null
+++ b/src/edgeMesh/edgeMeshFormats/extendedFeatureEdgeMesh/extendedFeatureEdgeMeshFormatRunTime.C
@@ -0,0 +1,50 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "extendedFeatureEdgeMeshFormat.H"
+
+#include "addToRunTimeSelectionTable.H"
+#include "addToMemberFunctionSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fileFormats
+{
+
+// read extendedEdgeMesh
+addNamedToRunTimeSelectionTable
+(
+    edgeMesh,
+    extendedFeatureEdgeMeshFormat,
+    fileExtension,
+    featureEdgeMesh
+);
+
+}
+}
+
+// ************************************************************************* //
diff --git a/src/edgeMesh/edgeMeshIO.C b/src/edgeMesh/edgeMeshIO.C
index 08e1398a36..3887f9d679 100644
--- a/src/edgeMesh/edgeMeshIO.C
+++ b/src/edgeMesh/edgeMeshIO.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -122,9 +122,9 @@ void Foam::edgeMesh::write
 
 void Foam::edgeMesh::writeStats(Ostream& os) const
 {
-    os  << "points      : " << points().size() << nl;
-    os  << "edges       : " << edges().size() << nl;
-    os  << "boundingBox : " << boundBox(this->points()) << endl;
+    os  << indent << "points      : " << points().size() << nl;
+    os  << indent << "edges       : " << edges().size() << nl;
+    os  << indent << "boundingBox : " << boundBox(this->points()) << endl;
 }
 
 
diff --git a/src/edgeMesh/edgeMeshNew.C b/src/edgeMesh/edgeMeshNew.C
index da7091c73a..df52d7cf63 100644
--- a/src/edgeMesh/edgeMeshNew.C
+++ b/src/edgeMesh/edgeMeshNew.C
@@ -42,8 +42,9 @@ Foam::autoPtr<Foam::edgeMesh> Foam::edgeMesh::New
         (
             "edgeMesh<Face>::New(const fileName&, const word&) : "
             "constructing edgeMesh"
-        )   << "Unknown file extension " << ext << nl << nl
-            << "Valid types are :" << nl
+        )   << "Unknown file extension " << ext
+            << " for file " << name << nl << nl
+            << "Valid extensions are :" << nl
             << fileExtensionConstructorTablePtr_->sortedToc()
             << exit(FatalError);
     }
diff --git a/src/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C
new file mode 100644
index 0000000000..999d2f0651
--- /dev/null
+++ b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C
@@ -0,0 +1,1518 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "extendedEdgeMesh.H"
+#include "surfaceFeatures.H"
+#include "triSurface.H"
+#include "Random.H"
+#include "Time.H"
+#include "OBJstream.H"
+#include "DynamicField.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(extendedEdgeMesh, 0);
+
+    template<>
+    const char* Foam::NamedEnum
+    <
+        Foam::extendedEdgeMesh::pointStatus,
+        4
+    >::names[] =
+    {
+        "convex",
+        "concave",
+        "mixed",
+        "nonFeature"
+    };
+
+    template<>
+    const char* Foam::NamedEnum
+    <
+        Foam::extendedEdgeMesh::edgeStatus,
+        6
+    >::names[] =
+    {
+        "external",
+        "internal",
+        "flat",
+        "open",
+        "multiple",
+        "none"
+    };
+
+    template<>
+    const char* Foam::NamedEnum
+    <
+        Foam::extendedEdgeMesh::sideVolumeType,
+        4
+    >::names[] =
+    {
+        "inside",
+        "outside",
+        "both",
+        "neither"
+    };
+}
+
+const Foam::NamedEnum<Foam::extendedEdgeMesh::pointStatus, 4>
+    Foam::extendedEdgeMesh::pointStatusNames_;
+
+const Foam::NamedEnum<Foam::extendedEdgeMesh::edgeStatus, 6>
+    Foam::extendedEdgeMesh::edgeStatusNames_;
+
+const Foam::NamedEnum<Foam::extendedEdgeMesh::sideVolumeType, 4>
+    Foam::extendedEdgeMesh::sideVolumeTypeNames_;
+
+Foam::scalar Foam::extendedEdgeMesh::cosNormalAngleTol_ =
+    Foam::cos(degToRad(0.1));
+
+
+Foam::label Foam::extendedEdgeMesh::convexStart_ = 0;
+
+
+Foam::label Foam::extendedEdgeMesh::externalStart_ = 0;
+
+
+Foam::label Foam::extendedEdgeMesh::nPointTypes = 4;
+
+
+Foam::label Foam::extendedEdgeMesh::nEdgeTypes = 5;
+
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+Foam::extendedEdgeMesh::pointStatus
+Foam::extendedEdgeMesh::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::extendedEdgeMesh::edgeStatus
+Foam::extendedEdgeMesh::classifyEdge
+(
+    const List<vector>& norms,
+    const labelList& edNorms,
+    const vector& fC0tofC1
+)
+{
+    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;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::extendedEdgeMesh::extendedEdgeMesh()
+:
+    edgeMesh(pointField(0), edgeList(0)),
+    concaveStart_(0),
+    mixedStart_(0),
+    nonFeatureStart_(0),
+    internalStart_(0),
+    flatStart_(0),
+    openStart_(0),
+    multipleStart_(0),
+    normals_(0),
+    normalVolumeTypes_(0),
+    edgeDirections_(0),
+    normalDirections_(0),
+    edgeNormals_(0),
+    featurePointNormals_(0),
+    featurePointEdges_(0),
+    regionEdges_(0),
+    pointTree_(),
+    edgeTree_(),
+    edgeTreesByType_()
+{}
+
+
+Foam::extendedEdgeMesh::extendedEdgeMesh
+(
+    const extendedEdgeMesh& fem
+)
+:
+    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()),
+    normalVolumeTypes_(fem.normalVolumeTypes()),
+    edgeDirections_(fem.edgeDirections()),
+    normalDirections_(fem.normalDirections()),
+    edgeNormals_(fem.edgeNormals()),
+    featurePointNormals_(fem.featurePointNormals()),
+    featurePointEdges_(fem.featurePointEdges()),
+    regionEdges_(fem.regionEdges()),
+    pointTree_(),
+    edgeTree_(),
+    edgeTreesByType_()
+{}
+
+
+Foam::extendedEdgeMesh::extendedEdgeMesh(Istream& is)
+{
+    is >> *this;
+}
+
+
+Foam::extendedEdgeMesh::extendedEdgeMesh
+(
+    const Xfer<pointField>& pointLst,
+    const Xfer<edgeList>& edgeLst
+)
+:
+    edgeMesh(pointLst, edgeLst),
+    concaveStart_(0),
+    mixedStart_(0),
+    nonFeatureStart_(0),
+    internalStart_(0),
+    flatStart_(0),
+    openStart_(0),
+    multipleStart_(0),
+    normals_(0),
+    normalVolumeTypes_(0),
+    edgeDirections_(0),
+    normalDirections_(0),
+    edgeNormals_(0),
+    featurePointNormals_(0),
+    featurePointEdges_(0),
+    regionEdges_(0),
+    pointTree_(),
+    edgeTree_(),
+    edgeTreesByType_()
+{}
+
+
+Foam::extendedEdgeMesh::extendedEdgeMesh
+(
+    const surfaceFeatures& sFeat,
+    const boolList& surfBaffleRegions
+)
+:
+    edgeMesh(pointField(0), edgeList(0)),
+    concaveStart_(-1),
+    mixedStart_(-1),
+    nonFeatureStart_(-1),
+    internalStart_(-1),
+    flatStart_(-1),
+    openStart_(-1),
+    multipleStart_(-1),
+    normals_(0),
+    normalVolumeTypes_(0),
+    edgeDirections_(0),
+    normalDirections_(0),
+    edgeNormals_(0),
+    featurePointNormals_(0),
+    featurePointEdges_(0),
+    regionEdges_(0),
+    pointTree_(),
+    edgeTree_(),
+    edgeTreesByType_()
+{
+    // Extract and reorder the data from surfaceFeatures
+    const triSurface& surf = sFeat.surface();
+    const labelList& featureEdges = sFeat.featureEdges();
+    const labelList& featurePoints = sFeat.featurePoints();
+
+    // Get a labelList of all the featureEdges that are region edges
+    const labelList regionFeatureEdges(identity(sFeat.nRegionEdges()));
+
+    sortPointsAndEdges
+    (
+        surf,
+        featureEdges,
+        regionFeatureEdges,
+        featurePoints
+    );
+
+    const labelListList& edgeFaces = surf.edgeFaces();
+
+    normalVolumeTypes_.setSize(normals_.size());
+
+    // Noting when the normal of a face has been used so not to duplicate
+    labelList faceMap(surf.size(), -1);
+
+    label nAdded = 0;
+
+    forAll(featureEdges, i)
+    {
+        label sFEI = featureEdges[i];
+
+        // Pick up the faces adjacent to the feature edge
+        const labelList& eFaces = edgeFaces[sFEI];
+
+        forAll(eFaces, j)
+        {
+            label eFI = eFaces[j];
+
+            // Check to see if the points have been already used
+            if (faceMap[eFI] == -1)
+            {
+                normalVolumeTypes_[nAdded++] =
+                    (
+                        surfBaffleRegions[surf[eFI].region()]
+                      ? BOTH
+                      : INSIDE
+                    );
+
+                faceMap[eFI] = nAdded - 1;
+            }
+        }
+    }
+}
+
+
+Foam::extendedEdgeMesh::extendedEdgeMesh
+(
+    const PrimitivePatch<face, List, pointField, point>& surf,
+    const labelList& featureEdges,
+    const labelList& regionFeatureEdges,
+    const labelList& featurePoints
+)
+:
+    edgeMesh(pointField(0), edgeList(0)),
+    concaveStart_(-1),
+    mixedStart_(-1),
+    nonFeatureStart_(-1),
+    internalStart_(-1),
+    flatStart_(-1),
+    openStart_(-1),
+    multipleStart_(-1),
+    normals_(0),
+    normalVolumeTypes_(0),
+    edgeDirections_(0),
+    normalDirections_(0),
+    edgeNormals_(0),
+    featurePointNormals_(0),
+    featurePointEdges_(0),
+    regionEdges_(0),
+    pointTree_(),
+    edgeTree_(),
+    edgeTreesByType_()
+{
+    sortPointsAndEdges
+    (
+        surf,
+        featureEdges,
+        regionFeatureEdges,
+        featurePoints
+    );
+}
+
+
+Foam::extendedEdgeMesh::extendedEdgeMesh
+(
+    const pointField& pts,
+    const edgeList& eds,
+    label concaveStart,
+    label mixedStart,
+    label nonFeatureStart,
+    label internalStart,
+    label flatStart,
+    label openStart,
+    label multipleStart,
+    const vectorField& normals,
+    const List<sideVolumeType>& normalVolumeTypes,
+    const vectorField& edgeDirections,
+    const labelListList& normalDirections,
+    const labelListList& edgeNormals,
+    const labelListList& featurePointNormals,
+    const labelListList& featurePointEdges,
+    const labelList& regionEdges
+)
+:
+    edgeMesh(pts, eds),
+    concaveStart_(concaveStart),
+    mixedStart_(mixedStart),
+    nonFeatureStart_(nonFeatureStart),
+    internalStart_(internalStart),
+    flatStart_(flatStart),
+    openStart_(openStart),
+    multipleStart_(multipleStart),
+    normals_(normals),
+    normalVolumeTypes_(normalVolumeTypes),
+    edgeDirections_(edgeDirections),
+    normalDirections_(normalDirections),
+    edgeNormals_(edgeNormals),
+    featurePointNormals_(featurePointNormals),
+    featurePointEdges_(featurePointEdges),
+    regionEdges_(regionEdges),
+    pointTree_(),
+    edgeTree_(),
+    edgeTreesByType_()
+{}
+
+
+Foam::extendedEdgeMesh::extendedEdgeMesh
+(
+    const fileName& name,
+    const word& ext
+)
+:
+    edgeMesh(pointField(0), edgeList(0)),
+    concaveStart_(0),
+    mixedStart_(0),
+    nonFeatureStart_(0),
+    internalStart_(0),
+    flatStart_(0),
+    openStart_(0),
+    multipleStart_(0),
+    normals_(0),
+    normalVolumeTypes_(0),
+    edgeDirections_(0),
+    normalDirections_(0),
+    edgeNormals_(0),
+    featurePointNormals_(0),
+    featurePointEdges_(0),
+    regionEdges_(0),
+    pointTree_(),
+    edgeTree_(),
+    edgeTreesByType_()
+{
+    read(name, ext);
+}
+
+
+Foam::extendedEdgeMesh::extendedEdgeMesh(const fileName& name)
+:
+    edgeMesh(pointField(0), edgeList(0)),
+    concaveStart_(0),
+    mixedStart_(0),
+    nonFeatureStart_(0),
+    internalStart_(0),
+    flatStart_(0),
+    openStart_(0),
+    multipleStart_(0),
+    normals_(0),
+    normalVolumeTypes_(0),
+    edgeDirections_(0),
+    normalDirections_(0),
+    edgeNormals_(0),
+    featurePointNormals_(0),
+    featurePointEdges_(0),
+    regionEdges_(0),
+    pointTree_(),
+    edgeTree_(),
+    edgeTreesByType_()
+{
+    read(name);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::extendedEdgeMesh::~extendedEdgeMesh()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+bool Foam::extendedEdgeMesh::read(const fileName& name)
+{
+    word ext = name.ext();
+    if (ext == "gz")
+    {
+        fileName unzipName = name.lessExt();
+        return read(unzipName, unzipName.ext());
+    }
+    else
+    {
+        return read(name, ext);
+    }
+}
+
+
+// Read from file in given format
+bool Foam::extendedEdgeMesh::read
+(
+    const fileName& name,
+    const word& ext
+)
+{
+    // read via selector mechanism
+    transfer(New(name, ext)());
+    return true;
+}
+
+
+void Foam::extendedEdgeMesh::nearestFeaturePoint
+(
+    const point& sample,
+    scalar searchDistSqr,
+    pointIndexHit& info
+) const
+{
+    info = pointTree().findNearest
+    (
+        sample,
+        searchDistSqr
+    );
+}
+
+
+void Foam::extendedEdgeMesh::nearestFeatureEdge
+(
+    const point& sample,
+    scalar searchDistSqr,
+    pointIndexHit& info
+) const
+{
+    info = edgeTree().findNearest
+    (
+        sample,
+        searchDistSqr
+    );
+}
+
+
+void Foam::extendedEdgeMesh::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::extendedEdgeMesh::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]);
+    }
+}
+
+
+void Foam::extendedEdgeMesh::allNearestFeaturePoints
+(
+    const point& sample,
+    scalar searchRadiusSqr,
+    List<pointIndexHit>& info
+) const
+{
+    // Pick up all the feature points that intersect the search sphere
+    labelList elems = pointTree().findSphere
+    (
+        sample,
+        searchRadiusSqr
+    );
+
+    DynamicList<pointIndexHit> dynPointHit(elems.size());
+
+    forAll(elems, elemI)
+    {
+        label index = elems[elemI];
+        label ptI = pointTree().shapes().pointLabels()[index];
+        const point& pt = points()[ptI];
+
+        pointIndexHit nearHit(true, pt, index);
+
+        dynPointHit.append(nearHit);
+    }
+
+    info.transfer(dynPointHit);
+}
+
+
+void Foam::extendedEdgeMesh::allNearestFeatureEdges
+(
+    const point& sample,
+    const scalar searchRadiusSqr,
+    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_;
+
+    DynamicList<pointIndexHit> dynEdgeHit(edgeTrees.size()*3);
+
+    // Loop over all the feature edge types
+    forAll(edgeTrees, i)
+    {
+        // Pick up all the edges that intersect the search sphere
+        labelList elems = edgeTrees[i].findSphere
+        (
+            sample,
+            searchRadiusSqr
+        );
+
+        forAll(elems, elemI)
+        {
+            label index = elems[elemI];
+            label edgeI = edgeTrees[i].shapes().edgeLabels()[index];
+            const edge& e = edges()[edgeI];
+
+            pointHit hitPoint = e.line(points()).nearestDist(sample);
+
+            label hitIndex = index + sliceStarts[i];
+
+            pointIndexHit nearHit
+            (
+                hitPoint.hit(),
+                hitPoint.rawPoint(),
+                hitIndex
+            );
+
+            dynEdgeHit.append(nearHit);
+        }
+    }
+
+    info.transfer(dynEdgeHit);
+}
+
+
+const Foam::indexedOctree<Foam::treeDataPoint>&
+Foam::extendedEdgeMesh::pointTree() const
+{
+    if (pointTree_.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);
+
+        const labelList featurePointLabels = identity(nonFeatureStart_);
+
+        pointTree_.reset
+        (
+            new indexedOctree<treeDataPoint>
+            (
+                treeDataPoint
+                (
+                    points(),
+                    featurePointLabels
+                ),
+                bb,     // bb
+                8,      // maxLevel
+                10,     // leafsize
+                3.0     // duplicity
+            )
+        );
+    }
+
+    return pointTree_();
+}
+
+
+const Foam::indexedOctree<Foam::treeDataEdge>&
+Foam::extendedEdgeMesh::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::extendedEdgeMesh::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::extendedEdgeMesh::transfer(extendedEdgeMesh& mesh)
+{
+    edgeMesh::transfer(mesh);
+
+    concaveStart_ = mesh.concaveStart_;
+    mixedStart_ = mesh.mixedStart_;
+    nonFeatureStart_ = mesh.nonFeatureStart_;
+    internalStart_ = mesh.internalStart_;
+    flatStart_ = mesh.flatStart_;
+    openStart_ = mesh.openStart_;
+    multipleStart_ = mesh.multipleStart_;
+    normals_.transfer(mesh.normals_);
+    normalVolumeTypes_.transfer(mesh.normalVolumeTypes_);
+    edgeDirections_.transfer(mesh.edgeDirections_);
+    normalDirections_.transfer(mesh.normalDirections_);
+    edgeNormals_.transfer(mesh.edgeNormals_);
+    featurePointNormals_.transfer(mesh.featurePointNormals_);
+    featurePointEdges_.transfer(mesh.featurePointEdges_);
+    regionEdges_.transfer(mesh.regionEdges_);
+    pointTree_ = mesh.pointTree_;
+    edgeTree_ = mesh.edgeTree_;
+    edgeTreesByType_.transfer(mesh.edgeTreesByType_);
+}
+
+
+Foam::Xfer<Foam::extendedEdgeMesh> Foam::extendedEdgeMesh::xfer()
+{
+    return xferMoveTo<extendedEdgeMesh, extendedEdgeMesh>(*this);
+}
+
+
+void Foam::extendedEdgeMesh::clear()
+{
+    edgeMesh::clear();
+    concaveStart_ = 0;
+    mixedStart_ = 0;
+    nonFeatureStart_ = 0;
+    internalStart_ = 0;
+    flatStart_ = 0;
+    openStart_ = 0;
+    multipleStart_ = 0;
+    normals_.clear();
+    normalVolumeTypes_.clear();
+    edgeDirections_.clear();
+    normalDirections_.clear();
+    edgeNormals_.clear();
+    featurePointNormals_.clear();
+    featurePointEdges_.clear();
+    regionEdges_.clear();
+    pointTree_.clear();
+    edgeTree_.clear();
+    edgeTreesByType_.clear();
+}
+
+
+void Foam::extendedEdgeMesh::add(const extendedEdgeMesh& fem)
+{
+    // Points
+    // ~~~~~~
+
+    // From current points into combined points
+    labelList reversePointMap(points().size());
+    labelList reverseFemPointMap(fem.points().size());
+
+    label newPointI = 0;
+    for (label i = 0; i < concaveStart(); i++)
+    {
+        reversePointMap[i] = newPointI++;
+    }
+    for (label i = 0; i < fem.concaveStart(); i++)
+    {
+        reverseFemPointMap[i] = newPointI++;
+    }
+
+    // Concave
+    label newConcaveStart = newPointI;
+    for (label i = concaveStart(); i < mixedStart(); i++)
+    {
+        reversePointMap[i] = newPointI++;
+    }
+    for (label i = fem.concaveStart(); i < fem.mixedStart(); i++)
+    {
+        reverseFemPointMap[i] = newPointI++;
+    }
+
+    // Mixed
+    label newMixedStart = newPointI;
+    for (label i = mixedStart(); i < nonFeatureStart(); i++)
+    {
+        reversePointMap[i] = newPointI++;
+    }
+    for (label i = fem.mixedStart(); i < fem.nonFeatureStart(); i++)
+    {
+        reverseFemPointMap[i] = newPointI++;
+    }
+
+    // Non-feature
+    label newNonFeatureStart = newPointI;
+    for (label i = nonFeatureStart(); i < points().size(); i++)
+    {
+        reversePointMap[i] = newPointI++;
+    }
+    for (label i = fem.nonFeatureStart(); i < fem.points().size(); i++)
+    {
+        reverseFemPointMap[i] = newPointI++;
+    }
+
+    pointField newPoints(newPointI);
+    newPoints.rmap(points(), reversePointMap);
+    newPoints.rmap(fem.points(), reverseFemPointMap);
+
+
+    // Edges
+    // ~~~~~
+
+    // From current edges into combined edges
+    labelList reverseEdgeMap(edges().size());
+    labelList reverseFemEdgeMap(fem.edges().size());
+
+    // External
+    label newEdgeI = 0;
+    for (label i = 0; i < internalStart(); i++)
+    {
+        reverseEdgeMap[i] = newEdgeI++;
+    }
+    for (label i = 0; i < fem.internalStart(); i++)
+    {
+        reverseFemEdgeMap[i] = newEdgeI++;
+    }
+
+    // Internal
+    label newInternalStart = newEdgeI;
+    for (label i = internalStart(); i < flatStart(); i++)
+    {
+        reverseEdgeMap[i] = newEdgeI++;
+    }
+    for (label i = fem.internalStart(); i < fem.flatStart(); i++)
+    {
+        reverseFemEdgeMap[i] = newEdgeI++;
+    }
+
+    // Flat
+    label newFlatStart = newEdgeI;
+    for (label i = flatStart(); i < openStart(); i++)
+    {
+        reverseEdgeMap[i] = newEdgeI++;
+    }
+    for (label i = fem.flatStart(); i < fem.openStart(); i++)
+    {
+        reverseFemEdgeMap[i] = newEdgeI++;
+    }
+
+    // Open
+    label newOpenStart = newEdgeI;
+    for (label i = openStart(); i < multipleStart(); i++)
+    {
+        reverseEdgeMap[i] = newEdgeI++;
+    }
+    for (label i = fem.openStart(); i < fem.multipleStart(); i++)
+    {
+        reverseFemEdgeMap[i] = newEdgeI++;
+    }
+
+    // Multiple
+    label newMultipleStart = newEdgeI;
+    for (label i = multipleStart(); i < edges().size(); i++)
+    {
+        reverseEdgeMap[i] = newEdgeI++;
+    }
+    for (label i = fem.multipleStart(); i < fem.edges().size(); i++)
+    {
+        reverseFemEdgeMap[i] = newEdgeI++;
+    }
+
+    edgeList newEdges(newEdgeI);
+    forAll(edges(), i)
+    {
+        const edge& e = edges()[i];
+        newEdges[reverseEdgeMap[i]] = edge
+        (
+            reversePointMap[e[0]],
+            reversePointMap[e[1]]
+        );
+    }
+    forAll(fem.edges(), i)
+    {
+        const edge& e = fem.edges()[i];
+        newEdges[reverseFemEdgeMap[i]] = edge
+        (
+            reverseFemPointMap[e[0]],
+            reverseFemPointMap[e[1]]
+        );
+    }
+
+    pointField newEdgeDirections(newEdgeI);
+    newEdgeDirections.rmap(edgeDirections(), reverseEdgeMap);
+    newEdgeDirections.rmap(fem.edgeDirections(), reverseFemEdgeMap);
+
+
+
+
+    // Normals
+    // ~~~~~~~
+
+    // Combine normals
+    DynamicField<point> newNormals(normals().size()+fem.normals().size());
+    newNormals.append(normals());
+    newNormals.append(fem.normals());
+
+
+    // Combine and re-index into newNormals
+    labelListList newEdgeNormals(edgeNormals().size()+fem.edgeNormals().size());
+    UIndirectList<labelList>(newEdgeNormals, reverseEdgeMap) =
+        edgeNormals();
+    UIndirectList<labelList>(newEdgeNormals, reverseFemEdgeMap) =
+        fem.edgeNormals();
+    forAll(reverseFemEdgeMap, i)
+    {
+        label mapI = reverseFemEdgeMap[i];
+        labelList& en = newEdgeNormals[mapI];
+        forAll(en, j)
+        {
+            en[j] += normals().size();
+        }
+    }
+
+
+    // Combine and re-index into newFeaturePointNormals
+    labelListList newFeaturePointNormals
+    (
+       featurePointNormals().size()
+     + fem.featurePointNormals().size()
+    );
+
+    // Note: featurePointNormals only go up to nonFeatureStart
+    UIndirectList<labelList>
+    (
+        newFeaturePointNormals,
+        SubList<label>(reversePointMap, featurePointNormals().size())
+    ) = featurePointNormals();
+    UIndirectList<labelList>
+    (
+        newFeaturePointNormals,
+        SubList<label>(reverseFemPointMap, fem.featurePointNormals().size())
+    ) = fem.featurePointNormals();
+    forAll(fem.featurePointNormals(), i)
+    {
+        label mapI = reverseFemPointMap[i];
+        labelList& fn = newFeaturePointNormals[mapI];
+        forAll(fn, j)
+        {
+            fn[j] += normals().size();
+        }
+    }
+
+
+    // Combine regionEdges
+    DynamicList<label> newRegionEdges
+    (
+        regionEdges().size()
+      + fem.regionEdges().size()
+    );
+    forAll(regionEdges(), i)
+    {
+        newRegionEdges.append(reverseEdgeMap[regionEdges()[i]]);
+    }
+    forAll(fem.regionEdges(), i)
+    {
+        newRegionEdges.append(reverseFemEdgeMap[fem.regionEdges()[i]]);
+    }
+
+
+    // Assign
+    // ~~~~~~
+
+    // Transfer
+    concaveStart_ = newConcaveStart;
+    mixedStart_ = newMixedStart;
+    nonFeatureStart_ = newNonFeatureStart;
+
+    // Reset points and edges
+    reset(xferMove(newPoints), newEdges.xfer());
+
+    // Transfer
+    internalStart_ = newInternalStart;
+    flatStart_ = newFlatStart;
+    openStart_ = newOpenStart;
+    multipleStart_ = newMultipleStart;
+
+    edgeDirections_.transfer(newEdgeDirections);
+
+    normals_.transfer(newNormals);
+    edgeNormals_.transfer(newEdgeNormals);
+    featurePointNormals_.transfer(newFeaturePointNormals);
+
+    regionEdges_.transfer(newRegionEdges);
+
+    pointTree_.clear();
+    edgeTree_.clear();
+    edgeTreesByType_.clear();
+}
+
+
+void Foam::extendedEdgeMesh::flipNormals()
+{
+    // Points
+    // ~~~~~~
+
+    // From current points into new points
+    labelList reversePointMap(identity(points().size()));
+
+    // Flip convex and concave points
+
+    label newPointI = 0;
+    // Concave points become convex
+    for (label i = concaveStart(); i < mixedStart(); i++)
+    {
+        reversePointMap[i] = newPointI++;
+    }
+    // Convex points become concave
+    label newConcaveStart = newPointI;
+    for (label i = 0; i < concaveStart(); i++)
+    {
+        reversePointMap[i] = newPointI++;
+    }
+
+
+    // Edges
+    // ~~~~~~
+
+    // From current edges into new edges
+    labelList reverseEdgeMap(identity(edges().size()));
+
+    // Flip external and internal edges
+
+    label newEdgeI = 0;
+    // Internal become external
+    for (label i = internalStart(); i < flatStart(); i++)
+    {
+        reverseEdgeMap[i] = newEdgeI++;
+    }
+    // External become internal
+    label newInternalStart = newEdgeI;
+    for (label i = 0; i < internalStart(); i++)
+    {
+        reverseEdgeMap[i] = newEdgeI++;
+    }
+
+
+    pointField newPoints(points().size());
+    newPoints.rmap(points(), reversePointMap);
+
+    edgeList newEdges(edges().size());
+    forAll(edges(), i)
+    {
+        const edge& e = edges()[i];
+        newEdges[reverseEdgeMap[i]] = edge
+        (
+            reversePointMap[e[0]],
+            reversePointMap[e[1]]
+        );
+    }
+
+
+    // Normals are flipped
+    // ~~~~~~~~~~~~~~~~~~~
+
+    pointField newEdgeDirections(edges().size());
+    newEdgeDirections.rmap(-1.0*edgeDirections(), reverseEdgeMap);
+
+    pointField newNormals(-1.0*normals());
+
+    labelListList newEdgeNormals(edgeNormals().size());
+    UIndirectList<labelList>(newEdgeNormals, reverseEdgeMap) = edgeNormals();
+
+    labelListList newFeaturePointNormals(featurePointNormals().size());
+
+    // Note: featurePointNormals only go up to nonFeatureStart
+    UIndirectList<labelList>
+    (
+        newFeaturePointNormals,
+        SubList<label>(reversePointMap, featurePointNormals().size())
+    ) = featurePointNormals();
+
+    labelList newRegionEdges(regionEdges().size());
+    forAll(regionEdges(), i)
+    {
+        newRegionEdges[i] = reverseEdgeMap[regionEdges()[i]];
+    }
+
+    // Transfer
+    concaveStart_ = newConcaveStart;
+
+    // Reset points and edges
+    reset(xferMove(newPoints), newEdges.xfer());
+
+    // Transfer
+    internalStart_ = newInternalStart;
+
+    edgeDirections_.transfer(newEdgeDirections);
+    normals_.transfer(newNormals);
+    edgeNormals_.transfer(newEdgeNormals);
+    featurePointNormals_.transfer(newFeaturePointNormals);
+    regionEdges_.transfer(newRegionEdges);
+
+    pointTree_.clear();
+    edgeTree_.clear();
+    edgeTreesByType_.clear();
+}
+
+
+void Foam::extendedEdgeMesh::writeObj
+(
+    const fileName& prefix
+) const
+{
+    Info<< nl << "Writing extendedEdgeMesh components to " << prefix
+        << endl;
+
+    edgeMesh::write(prefix + "_edgeMesh.obj");
+
+    OBJstream convexFtPtStr(prefix + "_convexFeaturePts.obj");
+    Info<< "Writing convex feature points to " << convexFtPtStr.name() << endl;
+
+    for(label i = 0; i < concaveStart_; i++)
+    {
+        convexFtPtStr.write(points()[i]);
+    }
+
+    OBJstream concaveFtPtStr(prefix + "_concaveFeaturePts.obj");
+    Info<< "Writing concave feature points to "
+        << concaveFtPtStr.name() << endl;
+
+    for(label i = concaveStart_; i < mixedStart_; i++)
+    {
+        convexFtPtStr.write(points()[i]);
+    }
+
+    OBJstream mixedFtPtStr(prefix + "_mixedFeaturePts.obj");
+    Info<< "Writing mixed feature points to " << mixedFtPtStr.name() << endl;
+
+    for(label i = mixedStart_; i < nonFeatureStart_; i++)
+    {
+        mixedFtPtStr.write(points()[i]);
+    }
+
+    OBJstream mixedFtPtStructureStr(prefix + "_mixedFeaturePtsStructure.obj");
+    Info<< "Writing mixed feature point structure to "
+        << mixedFtPtStructureStr.name() << endl;
+
+    for(label i = mixedStart_; i < nonFeatureStart_; i++)
+    {
+        const labelList& ptEds = pointEdges()[i];
+
+        forAll(ptEds, j)
+        {
+            const edge& e = edges()[ptEds[j]];
+            mixedFtPtStructureStr.write
+            (
+                linePointRef(points()[e[0]],
+                points()[e[1]])
+            );
+        }
+    }
+
+    OBJstream externalStr(prefix + "_externalEdges.obj");
+    Info<< "Writing external edges to " << externalStr.name() << endl;
+
+    for (label i = externalStart_; i < internalStart_; i++)
+    {
+        const edge& e = edges()[i];
+        externalStr.write(linePointRef(points()[e[0]], points()[e[1]]));
+    }
+
+    OBJstream internalStr(prefix + "_internalEdges.obj");
+    Info<< "Writing internal edges to " << internalStr.name() << endl;
+
+    for (label i = internalStart_; i < flatStart_; i++)
+    {
+        const edge& e = edges()[i];
+        internalStr.write(linePointRef(points()[e[0]], points()[e[1]]));
+    }
+
+    OBJstream flatStr(prefix + "_flatEdges.obj");
+    Info<< "Writing flat edges to " << flatStr.name() << endl;
+
+    for (label i = flatStart_; i < openStart_; i++)
+    {
+        const edge& e = edges()[i];
+        flatStr.write(linePointRef(points()[e[0]], points()[e[1]]));
+    }
+
+    OBJstream openStr(prefix + "_openEdges.obj");
+    Info<< "Writing open edges to " << openStr.name() << endl;
+
+    for (label i = openStart_; i < multipleStart_; i++)
+    {
+        const edge& e = edges()[i];
+        openStr.write(linePointRef(points()[e[0]], points()[e[1]]));
+    }
+
+    OBJstream multipleStr(prefix + "_multipleEdges.obj");
+    Info<< "Writing multiple edges to " << multipleStr.name() << endl;
+
+    for (label i = multipleStart_; i < edges().size(); i++)
+    {
+        const edge& e = edges()[i];
+        multipleStr.write(linePointRef(points()[e[0]], points()[e[1]]));
+    }
+
+    OBJstream regionStr(prefix + "_regionEdges.obj");
+    Info<< "Writing region edges to " << regionStr.name() << endl;
+
+    forAll(regionEdges_, i)
+    {
+        const edge& e = edges()[regionEdges_[i]];
+        regionStr.write(linePointRef(points()[e[0]], points()[e[1]]));
+    }
+
+    OBJstream edgeDirsStr(prefix + "_edgeDirections.obj");
+    Info<< "Writing edge directions to " << edgeDirsStr.name() << endl;
+
+    forAll(edgeDirections_, i)
+    {
+        const vector& eVec = edgeDirections_[i];
+        const edge& e = edges()[i];
+
+        edgeDirsStr.write
+        (
+            linePointRef(points()[e.start()], eVec + points()[e.start()])
+        );
+    }
+}
+
+
+void Foam::extendedEdgeMesh::writeStats(Ostream& os) const
+{
+    edgeMesh::writeStats(os);
+
+    os  << indent << "point offsets :" << nl;
+    os  << incrIndent;
+    os  << indent << "convex feature points      : " << convexStart_ << nl;
+    os  << indent << "concave feature points     : " << concaveStart_ << nl;
+    os  << indent << "mixed feature points       : " << mixedStart_ << nl;
+    os  << indent << "other (non-feature) points : " << nonFeatureStart_ << nl;
+    os  << decrIndent;
+
+    os  << indent << "edge offsets :" << nl;
+    os  << incrIndent;
+    os  << indent << "external (convex angle) edges  : " << externalStart_
+        << nl;
+    os  << indent << "internal (concave angle) edges : " << internalStart_
+        << nl;
+    os  << indent << "flat region edges              : " << flatStart_ << nl;
+    os  << indent << "open edges                     : " << openStart_ << nl;
+    os  << indent << "multiply connected edges       : " << multipleStart_
+        << nl;
+    os  << decrIndent;
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Istream& Foam::operator>>
+(
+    Istream& is,
+    Foam::extendedEdgeMesh::sideVolumeType& vt
+)
+{
+    label type;
+    is  >> type;
+
+    vt = static_cast<Foam::extendedEdgeMesh::sideVolumeType>(type);
+
+    // Check state of Istream
+    is.check("operator>>(Istream&, sideVolumeType&)");
+
+    return is;
+}
+
+
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const Foam::extendedEdgeMesh::sideVolumeType& vt
+)
+{
+    os  << static_cast<label>(vt);
+
+    return os;
+}
+
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const extendedEdgeMesh& em)
+{
+    //fileFormats::extendedEdgeMeshFormat::write(os, em.points_, em.edges_);
+    os  << "// points" << nl
+        << em.points() << nl
+        << "// edges" << nl
+        << em.edges() << nl
+        << "// concaveStart mixedStart nonFeatureStart" << nl
+        << em.concaveStart_ << token::SPACE
+        << em.mixedStart_ << token::SPACE
+        << em.nonFeatureStart_ << nl
+        << "// internalStart flatStart openStart multipleStart" << nl
+        << em.internalStart_ << token::SPACE
+        << em.flatStart_ << token::SPACE
+        << em.openStart_ << token::SPACE
+        << em.multipleStart_ << nl
+        << "// normals" << nl
+        << em.normals_ << nl
+        << "// normal volume types" << nl
+        << em.normalVolumeTypes_ << nl
+        << "// normalDirections" << nl
+        << em.normalDirections_ << nl
+        << "// edgeNormals" << nl
+        << em.edgeNormals_ << nl
+        << "// featurePointNormals" << nl
+        << em.featurePointNormals_ << nl
+        << "// featurePointEdges" << nl
+        << em.featurePointEdges_ << nl
+        << "// regionEdges" << nl
+        << em.regionEdges_
+        << endl;
+
+    // Check state of Ostream
+    os.check("Ostream& operator<<(Ostream&, const extendedEdgeMesh&)");
+
+    return os;
+}
+
+
+Foam::Istream& Foam::operator>>(Istream& is, extendedEdgeMesh& em)
+{
+    //fileFormats::extendedEdgeMeshFormat::read(is, em.points_, em.edges_);
+    is  >> static_cast<edgeMesh&>(em)
+        >> em.concaveStart_
+        >> em.mixedStart_
+        >> em.nonFeatureStart_
+        >> em.internalStart_
+        >> em.flatStart_
+        >> em.openStart_
+        >> em.multipleStart_
+        >> em.normals_
+        >> em.normalVolumeTypes_
+        >> em.normalDirections_
+        >> em.edgeNormals_
+        >> em.featurePointNormals_
+        >> em.featurePointEdges_
+        >> em.regionEdges_;
+
+    // Check state of Istream
+    is.check("Istream& operator>>(Istream&, extendedEdgeMesh&)");
+
+    return is;
+}
+
+
+// ************************************************************************* //
diff --git a/src/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.H b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.H
new file mode 100644
index 0000000000..2648321888
--- /dev/null
+++ b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.H
@@ -0,0 +1,547 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::extendedEdgeMesh
+
+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_ .. mixedStart_-1      : concave points
+        mixedStart_ .. nonFeatureStart_-1   : 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
+    extendedEdgeMeshI.H
+    extendedEdgeMesh.C
+    extendedEdgeMeshNew.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef extendedEdgeMesh_H
+#define extendedEdgeMesh_H
+
+#include "edgeMesh.H"
+#include "indexedOctree.H"
+#include "treeDataEdge.H"
+#include "treeDataPoint.H"
+#include "PrimitivePatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class surfaceFeatures;
+class objectRegistry;
+
+// Forward declaration of friend functions and operators
+class extendedEdgeMesh;
+Istream& operator>>(Istream&, extendedEdgeMesh&);
+Ostream& operator<<(Ostream&, const extendedEdgeMesh&);
+
+/*---------------------------------------------------------------------------*\
+                       Class extendedEdgeMesh Declaration
+\*---------------------------------------------------------------------------*/
+
+class extendedEdgeMesh
+:
+    public edgeMesh
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("extendedEdgeMesh");
+
+    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
+    };
+
+    static const Foam::NamedEnum<pointStatus, 4> pointStatusNames_;
+
+    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)
+    };
+
+    static const Foam::NamedEnum<edgeStatus, 6> edgeStatusNames_;
+
+    //- Normals point to the outside
+    enum sideVolumeType
+    {
+        INSIDE  = 0,  // mesh inside
+        OUTSIDE = 1,  // mesh outside
+        BOTH    = 2,  // e.g. a baffle
+        NEITHER = 3   // not sure when this may be used
+    };
+
+    static const Foam::NamedEnum<sideVolumeType, 4> sideVolumeTypeNames_;
+
+    //- Angular closeness tolerance for treating normals as the same
+    static scalar cosNormalAngleTol_;
+
+
+protected:
+
+    // Static data
+
+        //- 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_;
+
+        //- Type per normal: which side of normal to mesh
+        List<sideVolumeType> normalVolumeTypes_;
+
+        //- Flat and open edges require the direction of the edge
+        vectorField edgeDirections_;
+
+        //- Starting directions for the edges.
+        //  This vector points to the half of the plane defined by the first
+        //  edge normal.
+        labelListList normalDirections_;
+
+        //- Indices of the normals that are adjacent to the feature edges
+        labelListList edgeNormals_;
+
+        //- Indices of the normals that are adjacent to the feature points
+        //  (only valid for 0..nonFeatureStart_-1)
+        labelListList featurePointNormals_;
+
+        //- Indices of feature edges attached to feature points. The edges are
+        //  ordered so that they can be circulated.
+        labelListList featurePointEdges_;
+
+        //- Feature edges which are on the boundary between regions
+        labelList regionEdges_;
+
+        //- Search tree for all feature points
+        mutable autoPtr<indexedOctree<treeDataPoint> > pointTree_;
+
+        //- 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;
+
+        template<class Patch>
+        void sortPointsAndEdges
+        (
+            const Patch&,
+            const labelList& featureEdges,
+            const labelList& regionFeatureEdges,
+            const labelList& feaurePoints
+        );
+
+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 null
+        extendedEdgeMesh();
+
+        //- Construct as copy
+        explicit extendedEdgeMesh(const extendedEdgeMesh&);
+
+        //- Construct from file name (uses extension to determine type)
+        extendedEdgeMesh(const fileName&);
+
+        //- Construct from file name (uses extension to determine type)
+        extendedEdgeMesh(const fileName&, const word& ext);
+
+        //- Construct from Istream
+        extendedEdgeMesh(Istream&);
+
+        //- Construct by transferring components (points, edges)
+        extendedEdgeMesh
+        (
+            const Xfer<pointField>&,
+            const Xfer<edgeList>&
+        );
+
+        //- Construct given a surface with selected edges,point
+        //  (surfaceFeatures), an objectRegistry and a
+        //  fileName to write to.
+        //  Extracts, classifies and reorders the data from surfaceFeatures.
+        extendedEdgeMesh
+        (
+            const surfaceFeatures& sFeat,
+            const boolList& surfBaffleRegions
+        );
+
+        //- Construct from PrimitivePatch
+        extendedEdgeMesh
+        (
+            const PrimitivePatch<face, List, pointField, point>& surf,
+            const labelList& featureEdges,
+            const labelList& regionFeatureEdges,
+            const labelList& featurePoints
+        );
+
+        //- Construct from all components
+        extendedEdgeMesh
+        (
+            const pointField& pts,
+            const edgeList& eds,
+            label concaveStart,
+            label mixedStart,
+            label nonFeatureStart,
+            label internalStart,
+            label flatStart,
+            label openStart,
+            label multipleStart,
+            const vectorField& normals,
+            const List<sideVolumeType>& normalVolumeTypes,
+            const vectorField& edgeDirections,
+            const labelListList& normalDirections,
+            const labelListList& edgeNormals,
+            const labelListList& featurePointNormals,
+            const labelListList& featurePointEdges,
+            const labelList& regionEdges
+        );
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            extendedEdgeMesh,
+            fileExtension,
+            (
+                const fileName& name
+            ),
+            (name)
+        );
+
+
+    // Selectors
+
+        //- Select constructed from filename (explicit extension)
+        static autoPtr<extendedEdgeMesh> New
+        (
+            const fileName&,
+            const word& ext
+        );
+
+        //- Select constructed from filename (implicit extension)
+        static autoPtr<extendedEdgeMesh> New(const fileName&);
+
+
+    //- Destructor
+    ~extendedEdgeMesh();
+
+
+    // Member Functions
+
+        // Find
+
+            //- Find nearest surface edge for the sample point.
+            void nearestFeaturePoint
+            (
+                const point& sample,
+                scalar searchDistSqr,
+                pointIndexHit& info
+            ) const;
+
+            //- 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;
+
+            //- Find all the feature points within searchDistSqr of sample
+            void allNearestFeaturePoints
+            (
+                const point& sample,
+                scalar searchRadiusSqr,
+                List<pointIndexHit>& info
+            ) const;
+
+            //- Find all the feature edges within searchDistSqr of sample
+            void allNearestFeatureEdges
+            (
+                const point& sample,
+                const scalar searchRadiusSqr,
+                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
+            inline const List<sideVolumeType>& normalVolumeTypes() const;
+
+            //- Return the edgeDirection vectors
+            inline const vectorField& edgeDirections() const;
+
+            //-
+            inline const labelListList& normalDirections() 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 edge labels for a given feature point. Edges are
+            //  ordered by the faces that they share. The edge labels
+            //  correspond to the entry in edges().
+            inline const labelListList& featurePointEdges() 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;
+
+            //- Return the baffle faces of a specified edge
+            inline PackedList<2> edgeBaffles(label edgeI) const;
+
+            //- Demand driven construction of octree for feature points
+            const indexedOctree<treeDataPoint>& pointTree() 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;
+
+
+        // Edit
+
+            //- Transfer the contents of the argument and annul the argument
+            void transfer(extendedEdgeMesh&);
+
+            //- Transfer contents to the Xfer container
+            Xfer<extendedEdgeMesh > xfer();
+
+            //- Clear all storage
+            virtual void clear();
+
+            //- Add extendedEdgeMesh. No filtering of duplicates.
+            void add(const extendedEdgeMesh&);
+
+            //- Flip normals. All concave become convex, all internal external
+            //  etc.
+            void flipNormals();
+
+
+        // Read
+
+            //- Read from file. Chooses reader based on explicit extension
+            bool read(const fileName&, const word& ext);
+
+            //- Read from file. Chooses reader based on detected extension
+            virtual bool read(const fileName&);
+
+
+        // Write
+
+            //- Write all components of the extendedEdgeMesh as obj files
+            void writeObj(const fileName& prefix) const;
+
+            //- Dump some information
+            virtual void writeStats(Ostream& os) const;
+
+            friend Istream& operator>>(Istream& is, sideVolumeType& vt);
+            friend Ostream& operator<<(Ostream& os, const sideVolumeType& vt);
+
+
+        //- Classify the type of feature edge.  Requires face centre 0 to face
+        //  centre 1 vector to distinguish internal from external
+        static edgeStatus classifyEdge
+        (
+            const List<vector>& norms,
+            const labelList& edNorms,
+            const vector& fC0tofC1
+        );
+
+
+        // Ostream Operator
+
+            friend Ostream& operator<<(Ostream&, const extendedEdgeMesh&);
+            friend Istream& operator>>(Istream&, extendedEdgeMesh&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "extendedEdgeMeshI.H"
+
+#ifdef NoRepository
+#   include "extendedEdgeMeshTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshFormats/extendedEdgeMeshFormat/extendedEdgeMeshFormat.C b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshFormats/extendedEdgeMeshFormat/extendedEdgeMeshFormat.C
new file mode 100644
index 0000000000..a8ee72c434
--- /dev/null
+++ b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshFormats/extendedEdgeMeshFormat/extendedEdgeMeshFormat.C
@@ -0,0 +1,103 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "extendedEdgeMeshFormat.H"
+#include "edgeMeshFormat.H"
+#include "IFstream.H"
+#include "Time.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fileFormats::extendedEdgeMeshFormat::extendedEdgeMeshFormat
+(
+    const fileName& filename
+)
+{
+    read(filename);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::fileFormats::extendedEdgeMeshFormat::read
+(
+    const fileName& filename
+)
+{
+    clear();
+
+    fileName dir = filename.path();
+    fileName caseName = dir.name();
+    fileName rootPath = dir.path();
+
+    // Construct dummy time to use as an objectRegistry
+    Time dummyTime
+    (
+        ".",        //rootPath,
+        ".",        //caseName,
+        "system",   //systemName,
+        "constant", //constantName,
+        false       //enableFunctionObjects
+    );
+
+    // Construct IOobject to re-use the headerOk & readHeader
+    // (so we can read ascii and binary)
+    IOobject io
+    (
+        filename,
+        dummyTime,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE,
+        false
+    );
+
+    if (!io.headerOk())
+    {
+        FatalErrorIn
+        ("fileFormats::extendedEdgeMeshFormat::read(const fileName&)")
+            << "Cannot read file " << filename
+            << exit(FatalError);
+    }
+
+    autoPtr<IFstream> isPtr(new IFstream(io.filePath()));
+    bool ok = false;
+    if (isPtr().good())
+    {
+        Istream& is = isPtr();
+        ok = io.readHeader(is);
+
+        if (ok)
+        {
+            // Use extendedEdgeMesh IO
+            is >> *this;
+            ok = is.good();
+        }
+    }
+
+    return ok;
+}
+
+
+// ************************************************************************* //
diff --git a/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshFormats/extendedEdgeMeshFormat/extendedEdgeMeshFormat.H b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshFormats/extendedEdgeMeshFormat/extendedEdgeMeshFormat.H
new file mode 100644
index 0000000000..0bd618e073
--- /dev/null
+++ b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshFormats/extendedEdgeMeshFormat/extendedEdgeMeshFormat.H
@@ -0,0 +1,95 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::fileFormats::extendedEdgeMeshFormat
+
+Description
+    Provide a means of reading/writing the single-file OpenFOAM
+    extendedEdgeMesh format
+
+SourceFiles
+    extendedEdgeMeshFormat.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef extendedEdgeMeshFormat_H
+#define extendedEdgeMeshFormat_H
+
+#include "extendedEdgeMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fileFormats
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class extendedEdgeMeshFormat Declaration
+\*---------------------------------------------------------------------------*/
+
+class extendedEdgeMeshFormat
+:
+    public extendedEdgeMesh
+{
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        extendedEdgeMeshFormat(const extendedEdgeMeshFormat&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const extendedEdgeMeshFormat&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from file name
+        extendedEdgeMeshFormat(const fileName&);
+
+
+    //- Destructor
+    virtual ~extendedEdgeMeshFormat()
+    {}
+
+
+    // Member Functions
+
+        //- Read from file
+        virtual bool read(const fileName&);
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fileFormats
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshFormats/extendedEdgeMeshFormat/extendedEdgeMeshFormatRunTime.C b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshFormats/extendedEdgeMeshFormat/extendedEdgeMeshFormatRunTime.C
new file mode 100644
index 0000000000..a98e8ecbe1
--- /dev/null
+++ b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshFormats/extendedEdgeMeshFormat/extendedEdgeMeshFormatRunTime.C
@@ -0,0 +1,51 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "extendedEdgeMeshFormat.H"
+#include "extendedEdgeMesh.H"
+
+#include "addToRunTimeSelectionTable.H"
+#include "addToMemberFunctionSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fileFormats
+{
+
+// read extendedEdgeMesh
+addNamedToRunTimeSelectionTable
+(
+    extendedEdgeMesh,
+    extendedEdgeMeshFormat,
+    fileExtension,
+    extendedFeatureEdgeMesh         // extension
+);
+
+}
+}
+
+// ************************************************************************* //
diff --git a/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshFormats/obj/OBJextendedEdgeFormat.C b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshFormats/obj/OBJextendedEdgeFormat.C
new file mode 100644
index 0000000000..d63ad27ed9
--- /dev/null
+++ b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshFormats/obj/OBJextendedEdgeFormat.C
@@ -0,0 +1,66 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "OBJextendedEdgeFormat.H"
+#include "OBJedgeFormat.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fileFormats::OBJextendedEdgeFormat::OBJextendedEdgeFormat
+(
+    const fileName& filename
+)
+{
+    read(filename);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::fileFormats::OBJextendedEdgeFormat::read(const fileName& filename)
+{
+    edgeMesh em(filename);
+
+    clear();
+
+    // Note: should transfer here instead
+    storedPoints() = em.points();
+    storedEdges() = em.edges();
+
+    return true;
+}
+
+
+void Foam::fileFormats::OBJextendedEdgeFormat::write
+(
+    const fileName& filename,
+    const extendedEdgeMesh& mesh
+)
+{
+    OBJedgeFormat::write(filename, mesh);
+}
+
+
+// ************************************************************************* //
diff --git a/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshFormats/obj/OBJextendedEdgeFormat.H b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshFormats/obj/OBJextendedEdgeFormat.H
new file mode 100644
index 0000000000..a65409f446
--- /dev/null
+++ b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshFormats/obj/OBJextendedEdgeFormat.H
@@ -0,0 +1,119 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::fileFormats::OBJextendedEdgeFormat
+
+Description
+    Provide a means of reading/writing Alias/Wavefront OBJ format.
+
+    Does not handle negative vertex indices.
+
+SourceFiles
+    OBJextendedEdgeFormat.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef OBJextendedEdgeFormat_H
+#define OBJextendedEdgeFormat_H
+
+#include "extendedEdgeMesh.H"
+//#include "IFstream.H"
+//#include "Ostream.H"
+#include "OFstream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fileFormats
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class OBJextendedEdgeFormat Declaration
+\*---------------------------------------------------------------------------*/
+
+class OBJextendedEdgeFormat
+:
+    public extendedEdgeMesh
+{
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        OBJextendedEdgeFormat(const OBJextendedEdgeFormat&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const OBJextendedEdgeFormat&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from file name
+        OBJextendedEdgeFormat(const fileName&);
+
+
+//    // Selectors
+//
+//        //- Read file and return surface
+//        static autoPtr<extendedEdgeMesh> New(const fileName& name)
+//        {
+//            return autoPtr<extendedEdgeMesh>
+//            (
+//                new OBJextendedEdgeFormat(name)
+//            );
+//        }
+//
+
+    //- Destructor
+    virtual ~OBJextendedEdgeFormat()
+    {}
+
+
+    // Member Functions
+
+        //- Write surface mesh components by proxy
+        static void write(const fileName&, const extendedEdgeMesh&);
+
+        //- Read from file
+        virtual bool read(const fileName&);
+
+        //- Write object file
+        virtual void write(const fileName& name) const
+        {
+            write(name, *this);
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fileFormats
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshFormats/obj/OBJextendedEdgeFormatRunTime.C b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshFormats/obj/OBJextendedEdgeFormatRunTime.C
new file mode 100644
index 0000000000..c88513a104
--- /dev/null
+++ b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshFormats/obj/OBJextendedEdgeFormatRunTime.C
@@ -0,0 +1,60 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "OBJextendedEdgeFormat.H"
+
+#include "addToRunTimeSelectionTable.H"
+#include "addToMemberFunctionSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fileFormats
+{
+
+// read extendedEdgeMesh
+addNamedToRunTimeSelectionTable
+(
+    extendedEdgeMesh,
+    OBJextendedEdgeFormat,
+    fileExtension,
+    obj
+);
+
+//// write extendedEdgeMesh
+//addNamedToMemberFunctionSelectionTable
+//(
+//    extendedEdgeMesh,
+//    OBJextendedEdgeFormat,
+//    write,
+//    fileExtension,
+//    obj
+//);
+
+}
+}
+
+// ************************************************************************* //
diff --git a/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshI.H b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshI.H
new file mode 100644
index 0000000000..8506fae9e9
--- /dev/null
+++ b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshI.H
@@ -0,0 +1,293 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline Foam::label Foam::extendedEdgeMesh::convexStart() const
+{
+    return convexStart_;
+}
+
+
+inline Foam::label Foam::extendedEdgeMesh::concaveStart() const
+{
+    return concaveStart_;
+}
+
+
+inline Foam::label Foam::extendedEdgeMesh::mixedStart() const
+{
+    return mixedStart_;
+}
+
+
+inline Foam::label Foam::extendedEdgeMesh::nonFeatureStart() const
+{
+    return nonFeatureStart_;
+}
+
+
+inline Foam::label Foam::extendedEdgeMesh::externalStart() const
+{
+    return externalStart_;
+}
+
+
+inline Foam::label Foam::extendedEdgeMesh::internalStart() const
+{
+    return internalStart_;
+}
+
+
+inline Foam::label Foam::extendedEdgeMesh::flatStart() const
+{
+    return flatStart_;
+}
+
+
+inline Foam::label Foam::extendedEdgeMesh::openStart() const
+{
+    return openStart_;
+}
+
+
+inline Foam::label Foam::extendedEdgeMesh::multipleStart() const
+{
+    return multipleStart_;
+}
+
+
+inline bool Foam::extendedEdgeMesh::featurePoint(label ptI) const
+{
+    return ptI < nonFeatureStart_;
+}
+
+
+inline const Foam::vectorField& Foam::extendedEdgeMesh::normals() const
+{
+    return normals_;
+}
+
+
+inline const Foam::List<Foam::extendedEdgeMesh::sideVolumeType>&
+Foam::extendedEdgeMesh::normalVolumeTypes() const
+{
+    return normalVolumeTypes_;
+}
+
+
+inline const Foam::vectorField& Foam::extendedEdgeMesh::edgeDirections()
+const
+{
+    return edgeDirections_;
+}
+
+
+inline const Foam::labelListList&
+Foam::extendedEdgeMesh::normalDirections() const
+{
+    return normalDirections_;
+}
+
+
+inline Foam::vector Foam::extendedEdgeMesh::edgeDirection
+(
+    label edgeI,
+    label ptI
+) const
+{
+    const edge& e = edges()[edgeI];
+
+    if (ptI == e.start())
+    {
+        return edgeDirections()[edgeI];
+    }
+    else if (ptI == e.end())
+    {
+        return -edgeDirections()[edgeI];
+    }
+    else
+    {
+        FatalErrorIn("Foam::extendedEdgeMesh::edgeDirection")
+            << "Requested ptI " << ptI << " is not a point on the requested "
+            << "edgeI " << edgeI << ". edgeI start and end: "
+            << e.start() << " " << e.end()
+            << exit(FatalError);
+
+        return vector::zero;
+    }
+}
+
+
+inline const Foam::labelListList& Foam::extendedEdgeMesh::edgeNormals()
+const
+{
+    return edgeNormals_;
+}
+
+
+inline Foam::vectorField Foam::extendedEdgeMesh::edgeNormals
+(
+    const labelList& edgeNormIs
+) const
+{
+    vectorField norms(edgeNormIs.size());
+
+    forAll(edgeNormIs, i)
+    {
+        norms[i] = normals_[edgeNormIs[i]];
+    }
+
+    return norms;
+}
+
+
+inline Foam::vectorField Foam::extendedEdgeMesh::edgeNormals(label edgeI)
+const
+{
+    return edgeNormals(edgeNormals_[edgeI]);
+}
+
+
+inline const Foam::labelListList&
+Foam::extendedEdgeMesh::featurePointNormals() const
+{
+    return featurePointNormals_;
+}
+
+
+inline Foam::vectorField Foam::extendedEdgeMesh::featurePointNormals
+(
+    label ptI
+) const
+{
+    if (!featurePoint(ptI))
+    {
+        WarningIn("vectorField extendedEdgeMesh::featurePointNormals")
+            << "Requesting the normals of a non-feature point. "
+            << "Returned zero length vectorField."
+            << endl;
+
+        return vectorField(0);
+    }
+
+    labelList featPtNormIs(featurePointNormals_[ptI]);
+
+    vectorField norms(featPtNormIs.size());
+
+    forAll(featPtNormIs, i)
+    {
+        norms[i] = normals_[featPtNormIs[i]];
+    }
+
+    return norms;
+}
+
+
+inline const Foam::labelListList&
+Foam::extendedEdgeMesh::featurePointEdges() const
+{
+    return featurePointEdges_;
+}
+
+
+inline const Foam::labelList& Foam::extendedEdgeMesh::regionEdges() const
+{
+    return regionEdges_;
+}
+
+
+inline Foam::extendedEdgeMesh::pointStatus
+Foam::extendedEdgeMesh::getPointStatus(label ptI) const
+{
+    if (ptI < concaveStart_)
+    {
+        return CONVEX;
+    }
+    else if (ptI < mixedStart_)
+    {
+        return CONCAVE;
+    }
+    else if (ptI < nonFeatureStart_)
+    {
+        return MIXED;
+    }
+    else
+    {
+        return NONFEATURE;
+    }
+}
+
+
+inline Foam::extendedEdgeMesh::edgeStatus
+Foam::extendedEdgeMesh::getEdgeStatus(label edgeI) const
+{
+    if (edgeI < internalStart_)
+    {
+        return EXTERNAL;
+    }
+    else if (edgeI < flatStart_)
+    {
+        return INTERNAL;
+    }
+    else if (edgeI < openStart_)
+    {
+        return FLAT;
+    }
+    else if (edgeI < multipleStart_)
+    {
+        return OPEN;
+    }
+    else
+    {
+        return MULTIPLE;
+    }
+}
+
+
+inline Foam::PackedList<2> Foam::extendedEdgeMesh::edgeBaffles
+(
+    label edgeI
+) const
+{
+    const labelList& eNormals = edgeNormals_[edgeI];
+
+    DynamicList<label> edgeBaffles(eNormals.size());
+
+    forAll(eNormals, enI)
+    {
+        const label normI = eNormals[enI];
+
+        if (normalVolumeTypes_[normI])
+        {
+            edgeBaffles.append(normI);
+        }
+    }
+
+    return PackedList<2>(edgeBaffles);
+}
+
+
+// ************************************************************************* //
diff --git a/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshNew.C b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshNew.C
new file mode 100644
index 0000000000..5fa8572b6b
--- /dev/null
+++ b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshNew.C
@@ -0,0 +1,79 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "extendedEdgeMesh.H"
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineRunTimeSelectionTable(extendedEdgeMesh, fileExtension);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::extendedEdgeMesh> Foam::extendedEdgeMesh::New
+(
+    const fileName& name,
+    const word& ext
+)
+{
+    fileExtensionConstructorTable::iterator cstrIter =
+        fileExtensionConstructorTablePtr_->find(ext);
+
+    if (cstrIter == fileExtensionConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "extendedEdgeMesh::New(const fileName&, const word&) : "
+            "constructing extendedEdgeMesh"
+        )   << "Unknown file extension " << ext
+            << " for file " << name << nl << nl
+            << "Valid extensions are :" << nl
+            << fileExtensionConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<extendedEdgeMesh>(cstrIter()(name));
+}
+
+
+Foam::autoPtr<Foam::extendedEdgeMesh> Foam::extendedEdgeMesh::New
+(
+    const fileName& name
+)
+{
+    word ext = name.ext();
+    if (ext == "gz")
+    {
+        ext = name.lessExt().ext();
+    }
+    return New(name, ext);
+}
+
+
+// ************************************************************************* //
diff --git a/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshTemplates.C b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshTemplates.C
new file mode 100644
index 0000000000..cdb19b441d
--- /dev/null
+++ b/src/edgeMesh/extendedEdgeMesh/extendedEdgeMeshTemplates.C
@@ -0,0 +1,432 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "extendedEdgeMesh.H"
+#include "ListListOps.H"
+#include "unitConversion.H"
+#include "PackedBoolList.H"
+#include "PatchTools.H"
+#include "searchableBox.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Patch>
+void Foam::extendedEdgeMesh::sortPointsAndEdges
+(
+    const Patch& surf,
+    const labelList& featureEdges,
+    const labelList& regionFeatureEdges,// subset of featureEdges: inter-region
+    const labelList& featurePoints
+)
+{
+    const pointField& sFeatLocalPts(surf.localPoints());
+    const edgeList& sFeatEds(surf.edges());
+    const labelListList edgeFaces = PatchTools::sortedEdgeFaces(surf);
+    const vectorField& faceNormals = surf.faceNormals();
+    const labelListList pointEdges = PatchTools::sortedPointEdges(surf);
+
+    // Extract and reorder the data from surfaceFeatures
+
+    // References to the surfaceFeatures data
+
+    // Filling the extendedEdgeMesh with the raw geometrical data.
+
+    label nFeatEds = featureEdges.size();
+    label nFeatPts = featurePoints.size();
+
+    DynamicList<point> tmpPts;
+    edgeList eds(nFeatEds);
+    DynamicList<vector> norms;
+    vectorField edgeDirections(nFeatEds);
+    labelListList edgeNormals(nFeatEds);
+    labelListList normalDirections(nFeatEds);
+    DynamicList<label> regionEdges;
+
+    // Keep track of the ordered feature point feature edges
+    labelListList featurePointFeatureEdges(nFeatPts);
+    forAll(featurePointFeatureEdges, pI)
+    {
+        featurePointFeatureEdges[pI] = pointEdges[featurePoints[pI]];
+    }
+
+    // Mapping between old and new indices, there is entry in the map for each
+    // of surf.localPoints, -1 means that this point hasn't been used (yet),
+    // >= 0 corresponds to the index
+    labelList pointMap(sFeatLocalPts.size(), -1);
+
+    // Mapping between surface edge index and its feature edge index. -1 if it
+    // is not a feature edge
+    labelList edgeMap(sFeatEds.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(featurePoints, i)
+    {
+        label sFPI = featurePoints[i];
+
+        tmpPts.append(sFeatLocalPts[sFPI]);
+
+        pointMap[sFPI] = tmpPts.size() - 1;
+    }
+
+    // All feature points have been added
+    nonFeatureStart_ = tmpPts.size();
+
+    PackedBoolList isRegionFeatureEdge(regionFeatureEdges);
+
+    forAll(featureEdges, i)
+    {
+        label sFEI = featureEdges[i];
+
+        edgeMap[sFEI] = i;
+
+        const edge& fE = sFeatEds[sFEI];
+
+        edgeDirections[i] = fE.vec(sFeatLocalPts);
+
+        // 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 = edgeFaces[sFEI];
+
+        edgeNormals[i].setSize(eFaces.size());
+        normalDirections[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(faceNormals[eFI]);
+
+                faceMap[eFI] = norms.size() - 1;
+            }
+
+            edgeNormals[i][j] = faceMap[eFI];
+
+            const vector cross = (faceNormals[eFI] ^ edgeDirections[i]);
+            const vector fC0tofE0 =
+                surf[eFI].centre(surf.points())
+              - sFeatLocalPts[fE.start()];
+
+            normalDirections[i][j] =
+                (
+                    (
+                        (cross/(mag(cross) + VSMALL))
+                      & (fC0tofE0/(mag(fC0tofE0)+ VSMALL))
+                    )
+                  > 0.0
+                    ? 1
+                    : -1
+                );
+        }
+
+        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);
+
+        if (isRegionFeatureEdge[i])
+        {
+            regionEdges.append(i);
+        }
+    }
+
+    // Populate feature point feature edges
+    DynamicList<label> newFeatureEdges;
+
+    forAll(featurePointFeatureEdges, pI)
+    {
+        const labelList& fpfe = featurePointFeatureEdges[pI];
+
+        newFeatureEdges.setCapacity(fpfe.size());
+
+        forAll(fpfe, eI)
+        {
+            const label oldEdgeIndex = fpfe[eI];
+
+            const label newFeatureEdgeIndex = edgeMap[oldEdgeIndex];
+
+            if (newFeatureEdgeIndex != -1)
+            {
+                newFeatureEdges.append(newFeatureEdgeIndex);
+            }
+        }
+
+        featurePointFeatureEdges[pI].transfer(newFeatureEdges);
+    }
+
+    // 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
+            (
+                ":extendedEdgeMesh::sortPointsAndEdges"
+                "("
+                "    const Patch&,"
+                "    const labelList& featureEdges,"
+                "    const labelList& feaurePoints"
+                ")"
+            )   << 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);
+    inplaceReorder(edMap, normalDirections);
+    inplaceRenumber(edMap, regionEdges);
+
+    forAll(featurePointFeatureEdges, pI)
+    {
+        inplaceRenumber(edMap, featurePointFeatureEdges[pI]);
+    }
+
+    pointField pts(tmpPts);
+
+    // Initialise the edgeMesh
+    edgeMesh::operator=(edgeMesh(pts, eds));
+
+    // Initialise sorted edge related data
+    edgeDirections_ = edgeDirections/(mag(edgeDirections) + VSMALL);
+    edgeNormals_ = edgeNormals;
+    normalDirections_ = normalDirections;
+    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
+            (
+                ":extendedEdgeMesh::sortPointsAndEdges"
+                "("
+                "    const Patch&,"
+                "    const labelList& featureEdges,"
+                "    const labelList& feaurePoints"
+                ")"
+            )   << 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);
+    inplaceReorder(ptMap, featurePointFeatureEdges);
+
+    forAll(eds, i)
+    {
+        inplaceRenumber(ptMap, eds[i]);
+    }
+
+    // Reinitialise the edgeMesh with sorted feature points and
+    // renumbered edges
+    reset(xferMove(pts), xferMove(eds));
+
+    // Generate the featurePointNormals
+
+    labelListList featurePointNormals(nonFeatureStart_);
+
+    for (label i = 0; i < nonFeatureStart_; i++)
+    {
+        DynamicList<label> tmpFtPtNorms;
+
+        const labelList& ptEds = edgeMesh::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;
+    featurePointEdges_ = featurePointFeatureEdges;
+}
+
+
+// ************************************************************************* //
diff --git a/src/edgeMesh/extendedEdgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C b/src/edgeMesh/extendedEdgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C
index 0d03912cee..126b9dcdf8 100644
--- a/src/edgeMesh/extendedEdgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C
+++ b/src/edgeMesh/extendedEdgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C
@@ -24,114 +24,21 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "extendedFeatureEdgeMesh.H"
-#include "surfaceFeatures.H"
-#include "triSurface.H"
-#include "Random.H"
 #include "Time.H"
-#include "meshTools.H"
-#include "ListListOps.H"
-#include "OFstream.H"
-#include "IFstream.H"
-#include "unitConversion.H"
-#include "DynamicField.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
     defineTypeNameAndDebug(extendedFeatureEdgeMesh, 0);
-
-    template<>
-    const char* Foam::NamedEnum
-    <
-        Foam::extendedFeatureEdgeMesh::pointStatus,
-        4
-    >::names[] =
-    {
-        "convex",
-        "concave",
-        "mixed",
-        "nonFeature"
-    };
-
-    template<>
-    const char* Foam::NamedEnum
-    <
-        Foam::extendedFeatureEdgeMesh::edgeStatus,
-        6
-    >::names[] =
-    {
-        "external",
-        "internal",
-        "flat",
-        "open",
-        "multiple",
-        "none"
-    };
-
-    template<>
-    const char* Foam::NamedEnum
-    <
-        Foam::extendedFeatureEdgeMesh::sideVolumeType,
-        4
-    >::names[] =
-    {
-        "inside",
-        "outside",
-        "both",
-        "neither"
-    };
 }
 
-const Foam::NamedEnum<Foam::extendedFeatureEdgeMesh::pointStatus, 4>
-    Foam::extendedFeatureEdgeMesh::pointStatusNames_;
-
-const Foam::NamedEnum<Foam::extendedFeatureEdgeMesh::edgeStatus, 6>
-    Foam::extendedFeatureEdgeMesh::edgeStatusNames_;
-
-const Foam::NamedEnum<Foam::extendedFeatureEdgeMesh::sideVolumeType, 4>
-    Foam::extendedFeatureEdgeMesh::sideVolumeTypeNames_;
-
-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),
-    normalVolumeTypes_(0),
-    edgeDirections_(0),
-    normalDirections_(0),
-    edgeNormals_(0),
-    featurePointNormals_(0),
-    featurePointEdges_(0),
-    regionEdges_(0),
-    pointTree_(),
-    edgeTree_(),
-    edgeTreesByType_()
+    extendedEdgeMesh()
 {
     if
     (
@@ -151,24 +58,7 @@ Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh(const IOobject& io)
                 << endl;
         }
 
-        Istream& is = readStream(typeName);
-
-        is  >> *this
-            >> concaveStart_
-            >> mixedStart_
-            >> nonFeatureStart_
-            >> internalStart_
-            >> flatStart_
-            >> openStart_
-            >> multipleStart_
-            >> normals_
-            >> normalVolumeTypes_
-            >> normalDirections_
-            >> edgeNormals_
-            >> featurePointNormals_
-            >> featurePointEdges_
-            >> regionEdges_;
-
+        readStream(typeName) >> *this;
         close();
 
         {
@@ -203,59 +93,11 @@ Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh(const IOobject& io)
 Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh
 (
     const IOobject& io,
-    const extendedFeatureEdgeMesh& fem
+    const extendedEdgeMesh& em
 )
 :
     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()),
-    normalVolumeTypes_(fem.normalVolumeTypes()),
-    edgeDirections_(fem.edgeDirections()),
-    normalDirections_(fem.normalDirections()),
-    edgeNormals_(fem.edgeNormals()),
-    featurePointNormals_(fem.featurePointNormals()),
-    featurePointEdges_(fem.featurePointEdges()),
-    regionEdges_(fem.regionEdges()),
-    pointTree_(),
-    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),
-    normalVolumeTypes_(0),
-    edgeDirections_(0),
-    normalDirections_(0),
-    edgeNormals_(0),
-    featurePointNormals_(0),
-    featurePointEdges_(0),
-    regionEdges_(0),
-    pointTree_(),
-    edgeTree_(),
-    edgeTreesByType_()
+    extendedEdgeMesh(em)
 {}
 
 
@@ -279,77 +121,8 @@ Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh
             IOobject::NO_WRITE
         )
     ),
-    edgeMesh(pointField(0), edgeList(0)),
-    concaveStart_(-1),
-    mixedStart_(-1),
-    nonFeatureStart_(-1),
-    internalStart_(-1),
-    flatStart_(-1),
-    openStart_(-1),
-    multipleStart_(-1),
-    normals_(0),
-    normalVolumeTypes_(0),
-    edgeDirections_(0),
-    normalDirections_(0),
-    edgeNormals_(0),
-    featurePointNormals_(0),
-    featurePointEdges_(0),
-    regionEdges_(0),
-    pointTree_(),
-    edgeTree_(),
-    edgeTreesByType_()
-{
-    // Extract and reorder the data from surfaceFeatures
-    const triSurface& surf = sFeat.surface();
-    const labelList& featureEdges = sFeat.featureEdges();
-    const labelList& featurePoints = sFeat.featurePoints();
-
-    // Get a labelList of all the featureEdges that are region edges
-    const labelList regionFeatureEdges(identity(sFeat.nRegionEdges()));
-
-    sortPointsAndEdges
-    (
-        surf,
-        featureEdges,
-        regionFeatureEdges,
-        featurePoints
-    );
-
-    const labelListList& edgeFaces = surf.edgeFaces();
-
-    normalVolumeTypes_.setSize(normals_.size());
-
-    // Noting when the normal of a face has been used so not to duplicate
-    labelList faceMap(surf.size(), -1);
-
-    label nAdded = 0;
-
-    forAll(featureEdges, i)
-    {
-        label sFEI = featureEdges[i];
-
-        // Pick up the faces adjacent to the feature edge
-        const labelList& eFaces = edgeFaces[sFEI];
-
-        forAll(eFaces, j)
-        {
-            label eFI = eFaces[j];
-
-            // Check to see if the points have been already used
-            if (faceMap[eFI] == -1)
-            {
-                normalVolumeTypes_[nAdded++] =
-                    (
-                        surfBaffleRegions[surf[eFI].region()]
-                      ? BOTH
-                      : INSIDE
-                    );
-
-                faceMap[eFI] = nAdded - 1;
-            }
-        }
-    }
-}
+    extendedEdgeMesh(sFeat, surfBaffleRegions)
+{}
 
 
 Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh
@@ -362,34 +135,8 @@ Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh
 )
 :
     regIOobject(io),
-    edgeMesh(pointField(0), edgeList(0)),
-    concaveStart_(-1),
-    mixedStart_(-1),
-    nonFeatureStart_(-1),
-    internalStart_(-1),
-    flatStart_(-1),
-    openStart_(-1),
-    multipleStart_(-1),
-    normals_(0),
-    normalVolumeTypes_(0),
-    edgeDirections_(0),
-    normalDirections_(0),
-    edgeNormals_(0),
-    featurePointNormals_(0),
-    featurePointEdges_(0),
-    regionEdges_(0),
-    pointTree_(),
-    edgeTree_(),
-    edgeTreesByType_()
-{
-    sortPointsAndEdges
-    (
-        surf,
-        featureEdges,
-        regionFeatureEdges,
-        featurePoints
-    );
-}
+    extendedEdgeMesh(surf, featureEdges, regionFeatureEdges, featurePoints)
+{}
 
 
 Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh
@@ -415,25 +162,26 @@ Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh
 )
 :
     regIOobject(io),
-    edgeMesh(pts, eds),
-    concaveStart_(concaveStart),
-    mixedStart_(mixedStart),
-    nonFeatureStart_(nonFeatureStart),
-    internalStart_(internalStart),
-    flatStart_(flatStart),
-    openStart_(openStart),
-    multipleStart_(multipleStart),
-    normals_(normals),
-    normalVolumeTypes_(normalVolumeTypes),
-    edgeDirections_(edgeDirections),
-    normalDirections_(normalDirections),
-    edgeNormals_(edgeNormals),
-    featurePointNormals_(featurePointNormals),
-    featurePointEdges_(featurePointEdges),
-    regionEdges_(regionEdges),
-    pointTree_(),
-    edgeTree_(),
-    edgeTreesByType_()
+    extendedEdgeMesh
+    (
+        pts,
+        eds,
+        concaveStart,
+        mixedStart,
+        nonFeatureStart,
+        internalStart,
+        flatStart,
+        openStart,
+        multipleStart,
+        normals,
+        normalVolumeTypes,
+        edgeDirections,
+        normalDirections,
+        edgeNormals,
+        featurePointNormals,
+        featurePointEdges,
+        regionEdges
+    )
 {}
 
 
@@ -445,991 +193,88 @@ 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
-)
-{
-    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::nearestFeaturePoint
-(
-    const point& sample,
-    scalar searchDistSqr,
-    pointIndexHit& info
-) const
+bool Foam::extendedFeatureEdgeMesh::readData(Istream& is)
 {
-    info = pointTree().findNearest
-    (
-        sample,
-        searchDistSqr
-    );
-}
-
-
-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]);
-    }
-}
-
-
-void Foam::extendedFeatureEdgeMesh::allNearestFeaturePoints
-(
-    const point& sample,
-    scalar searchRadiusSqr,
-    List<pointIndexHit>& info
-) const
-{
-    // Pick up all the feature points that intersect the search sphere
-    labelList elems = pointTree().findSphere
-    (
-        sample,
-        searchRadiusSqr
-    );
-
-    DynamicList<pointIndexHit> dynPointHit(elems.size());
-
-    forAll(elems, elemI)
-    {
-        label index = elems[elemI];
-        label ptI = pointTree().shapes().pointLabels()[index];
-        const point& pt = points()[ptI];
-
-        pointIndexHit nearHit(true, pt, index);
-
-        dynPointHit.append(nearHit);
-    }
-
-    info.transfer(dynPointHit);
-}
-
-
-void Foam::extendedFeatureEdgeMesh::allNearestFeatureEdges
-(
-    const point& sample,
-    const scalar searchRadiusSqr,
-    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_;
-
-    DynamicList<pointIndexHit> dynEdgeHit(edgeTrees.size()*3);
-
-    // Loop over all the feature edge types
-    forAll(edgeTrees, i)
-    {
-        // Pick up all the edges that intersect the search sphere
-        labelList elems = edgeTrees[i].findSphere
-        (
-            sample,
-            searchRadiusSqr
-        );
-
-        forAll(elems, elemI)
-        {
-            label index = elems[elemI];
-            label edgeI = edgeTrees[i].shapes().edgeLabels()[index];
-            const edge& e = edges()[edgeI];
-
-            pointHit hitPoint = e.line(points()).nearestDist(sample);
-
-            label hitIndex = index + sliceStarts[i];
-
-            pointIndexHit nearHit
-            (
-                hitPoint.hit(),
-                hitPoint.rawPoint(),
-                hitIndex
-            );
-
-            dynEdgeHit.append(nearHit);
-        }
-    }
-
-    info.transfer(dynEdgeHit);
-}
-
-
-const Foam::indexedOctree<Foam::treeDataPoint>&
-Foam::extendedFeatureEdgeMesh::pointTree() const
-{
-    if (pointTree_.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);
-
-        const labelList featurePointLabels = identity(nonFeatureStart_);
-
-        pointTree_.reset
-        (
-            new indexedOctree<treeDataPoint>
-            (
-                treeDataPoint
-                (
-                    points(),
-                    featurePointLabels
-                ),
-                bb,     // bb
-                8,      // maxLevel
-                10,     // leafsize
-                3.0     // duplicity
-            )
-        );
-    }
-
-    return pointTree_();
-}
-
-
-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::add(const extendedFeatureEdgeMesh& fem)
-{
-    // Points
-    // ~~~~~~
-
-    // From current points into combined points
-    labelList reversePointMap(points().size());
-    labelList reverseFemPointMap(fem.points().size());
-
-    label newPointI = 0;
-    for (label i = 0; i < concaveStart(); i++)
-    {
-        reversePointMap[i] = newPointI++;
-    }
-    for (label i = 0; i < fem.concaveStart(); i++)
-    {
-        reverseFemPointMap[i] = newPointI++;
-    }
-
-    // Concave
-    label newConcaveStart = newPointI;
-    for (label i = concaveStart(); i < mixedStart(); i++)
-    {
-        reversePointMap[i] = newPointI++;
-    }
-    for (label i = fem.concaveStart(); i < fem.mixedStart(); i++)
-    {
-        reverseFemPointMap[i] = newPointI++;
-    }
-
-    // Mixed
-    label newMixedStart = newPointI;
-    for (label i = mixedStart(); i < nonFeatureStart(); i++)
-    {
-        reversePointMap[i] = newPointI++;
-    }
-    for (label i = fem.mixedStart(); i < fem.nonFeatureStart(); i++)
-    {
-        reverseFemPointMap[i] = newPointI++;
-    }
-
-    // Non-feature
-    label newNonFeatureStart = newPointI;
-    for (label i = nonFeatureStart(); i < points().size(); i++)
-    {
-        reversePointMap[i] = newPointI++;
-    }
-    for (label i = fem.nonFeatureStart(); i < fem.points().size(); i++)
-    {
-        reverseFemPointMap[i] = newPointI++;
-    }
-
-    pointField newPoints(newPointI);
-    newPoints.rmap(points(), reversePointMap);
-    newPoints.rmap(fem.points(), reverseFemPointMap);
-
-
-    // Edges
-    // ~~~~~
-
-    // From current edges into combined edges
-    labelList reverseEdgeMap(edges().size());
-    labelList reverseFemEdgeMap(fem.edges().size());
-
-    // External
-    label newEdgeI = 0;
-    for (label i = 0; i < internalStart(); i++)
-    {
-        reverseEdgeMap[i] = newEdgeI++;
-    }
-    for (label i = 0; i < fem.internalStart(); i++)
-    {
-        reverseFemEdgeMap[i] = newEdgeI++;
-    }
-
-    // Internal
-    label newInternalStart = newEdgeI;
-    for (label i = internalStart(); i < flatStart(); i++)
-    {
-        reverseEdgeMap[i] = newEdgeI++;
-    }
-    for (label i = fem.internalStart(); i < fem.flatStart(); i++)
-    {
-        reverseFemEdgeMap[i] = newEdgeI++;
-    }
-
-    // Flat
-    label newFlatStart = newEdgeI;
-    for (label i = flatStart(); i < openStart(); i++)
-    {
-        reverseEdgeMap[i] = newEdgeI++;
-    }
-    for (label i = fem.flatStart(); i < fem.openStart(); i++)
-    {
-        reverseFemEdgeMap[i] = newEdgeI++;
-    }
-
-    // Open
-    label newOpenStart = newEdgeI;
-    for (label i = openStart(); i < multipleStart(); i++)
-    {
-        reverseEdgeMap[i] = newEdgeI++;
-    }
-    for (label i = fem.openStart(); i < fem.multipleStart(); i++)
-    {
-        reverseFemEdgeMap[i] = newEdgeI++;
-    }
-
-    // Multiple
-    label newMultipleStart = newEdgeI;
-    for (label i = multipleStart(); i < edges().size(); i++)
-    {
-        reverseEdgeMap[i] = newEdgeI++;
-    }
-    for (label i = fem.multipleStart(); i < fem.edges().size(); i++)
-    {
-        reverseFemEdgeMap[i] = newEdgeI++;
-    }
-
-    edgeList newEdges(newEdgeI);
-    forAll(edges(), i)
-    {
-        const edge& e = edges()[i];
-        newEdges[reverseEdgeMap[i]] = edge
-        (
-            reversePointMap[e[0]],
-            reversePointMap[e[1]]
-        );
-    }
-    forAll(fem.edges(), i)
-    {
-        const edge& e = fem.edges()[i];
-        newEdges[reverseFemEdgeMap[i]] = edge
-        (
-            reverseFemPointMap[e[0]],
-            reverseFemPointMap[e[1]]
-        );
-    }
-
-    pointField newEdgeDirections(newEdgeI);
-    newEdgeDirections.rmap(edgeDirections(), reverseEdgeMap);
-    newEdgeDirections.rmap(fem.edgeDirections(), reverseFemEdgeMap);
-
-
-
-
-    // Normals
-    // ~~~~~~~
-
-    // Combine normals
-    DynamicField<point> newNormals(normals().size()+fem.normals().size());
-    newNormals.append(normals());
-    newNormals.append(fem.normals());
-
-
-    // Combine and re-index into newNormals
-    labelListList newEdgeNormals(edgeNormals().size()+fem.edgeNormals().size());
-    UIndirectList<labelList>(newEdgeNormals, reverseEdgeMap) =
-        edgeNormals();
-    UIndirectList<labelList>(newEdgeNormals, reverseFemEdgeMap) =
-        fem.edgeNormals();
-    forAll(reverseFemEdgeMap, i)
-    {
-        label mapI = reverseFemEdgeMap[i];
-        labelList& en = newEdgeNormals[mapI];
-        forAll(en, j)
-        {
-            en[j] += normals().size();
-        }
-    }
-
-
-    // Combine and re-index into newFeaturePointNormals
-    labelListList newFeaturePointNormals
-    (
-       featurePointNormals().size()
-     + fem.featurePointNormals().size()
-    );
-
-    // Note: featurePointNormals only go up to nonFeatureStart
-    UIndirectList<labelList>
-    (
-        newFeaturePointNormals,
-        SubList<label>(reversePointMap, featurePointNormals().size())
-    ) = featurePointNormals();
-    UIndirectList<labelList>
-    (
-        newFeaturePointNormals,
-        SubList<label>(reverseFemPointMap, fem.featurePointNormals().size())
-    ) = fem.featurePointNormals();
-    forAll(fem.featurePointNormals(), i)
-    {
-        label mapI = reverseFemPointMap[i];
-        labelList& fn = newFeaturePointNormals[mapI];
-        forAll(fn, j)
-        {
-            fn[j] += normals().size();
-        }
-    }
-
-
-    // Combine regionEdges
-    DynamicList<label> newRegionEdges
-    (
-        regionEdges().size()
-      + fem.regionEdges().size()
-    );
-    forAll(regionEdges(), i)
-    {
-        newRegionEdges.append(reverseEdgeMap[regionEdges()[i]]);
-    }
-    forAll(fem.regionEdges(), i)
-    {
-        newRegionEdges.append(reverseFemEdgeMap[fem.regionEdges()[i]]);
-    }
-
-
-    // Assign
-    // ~~~~~~
-
-    // Transfer
-    concaveStart_ = newConcaveStart;
-    mixedStart_ = newMixedStart;
-    nonFeatureStart_ = newNonFeatureStart;
-
-    // Reset points and edges
-    reset(xferMove(newPoints), newEdges.xfer());
-
-    // Transfer
-    internalStart_ = newInternalStart;
-    flatStart_ = newFlatStart;
-    openStart_ = newOpenStart;
-    multipleStart_ = newMultipleStart;
-
-    edgeDirections_.transfer(newEdgeDirections);
-
-    normals_.transfer(newNormals);
-    edgeNormals_.transfer(newEdgeNormals);
-    featurePointNormals_.transfer(newFeaturePointNormals);
-
-    regionEdges_.transfer(newRegionEdges);
-
-    pointTree_.clear();
-    edgeTree_.clear();
-    edgeTreesByType_.clear();
-}
-
-
-void Foam::extendedFeatureEdgeMesh::flipNormals()
-{
-    // Points
-    // ~~~~~~
-
-    // From current points into new points
-    labelList reversePointMap(identity(points().size()));
-
-    // Flip convex and concave points
-
-    label newPointI = 0;
-    // Concave points become convex
-    for (label i = concaveStart(); i < mixedStart(); i++)
-    {
-        reversePointMap[i] = newPointI++;
-    }
-    // Convex points become concave
-    label newConcaveStart = newPointI;
-    for (label i = 0; i < concaveStart(); i++)
-    {
-        reversePointMap[i] = newPointI++;
-    }
-
-
-    // Edges
-    // ~~~~~~
-
-    // From current edges into new edges
-    labelList reverseEdgeMap(identity(edges().size()));
-
-    // Flip external and internal edges
-
-    label newEdgeI = 0;
-    // Internal become external
-    for (label i = internalStart(); i < flatStart(); i++)
-    {
-        reverseEdgeMap[i] = newEdgeI++;
-    }
-    // External become internal
-    label newInternalStart = newEdgeI;
-    for (label i = 0; i < internalStart(); i++)
-    {
-        reverseEdgeMap[i] = newEdgeI++;
-    }
-
-
-    pointField newPoints(points().size());
-    newPoints.rmap(points(), reversePointMap);
-
-    edgeList newEdges(edges().size());
-    forAll(edges(), i)
-    {
-        const edge& e = edges()[i];
-        newEdges[reverseEdgeMap[i]] = edge
-        (
-            reversePointMap[e[0]],
-            reversePointMap[e[1]]
-        );
-    }
-
-
-    // Normals are flipped
-    // ~~~~~~~~~~~~~~~~~~~
-
-    pointField newEdgeDirections(edges().size());
-    newEdgeDirections.rmap(-1.0*edgeDirections(), reverseEdgeMap);
-
-    pointField newNormals(-1.0*normals());
-
-    labelListList newEdgeNormals(edgeNormals().size());
-    UIndirectList<labelList>(newEdgeNormals, reverseEdgeMap) = edgeNormals();
-
-    labelListList newFeaturePointNormals(featurePointNormals().size());
-
-    // Note: featurePointNormals only go up to nonFeatureStart
-    UIndirectList<labelList>
-    (
-        newFeaturePointNormals,
-        SubList<label>(reversePointMap, featurePointNormals().size())
-    ) = featurePointNormals();
-
-    labelList newRegionEdges(regionEdges().size());
-    forAll(regionEdges(), i)
-    {
-        newRegionEdges[i] = reverseEdgeMap[regionEdges()[i]];
-    }
-
-    // Transfer
-    concaveStart_ = newConcaveStart;
-
-    // Reset points and edges
-    reset(xferMove(newPoints), newEdges.xfer());
-
-    // Transfer
-    internalStart_ = newInternalStart;
-
-    edgeDirections_.transfer(newEdgeDirections);
-    normals_.transfer(newNormals);
-    edgeNormals_.transfer(newEdgeNormals);
-    featurePointNormals_.transfer(newFeaturePointNormals);
-    regionEdges_.transfer(newRegionEdges);
-
-    pointTree_.clear();
-    edgeTree_.clear();
-    edgeTreesByType_.clear();
-}
-//XXXXX
-
-void Foam::extendedFeatureEdgeMesh::writeObj
-(
-    const fileName& prefix
-) const
-{
-    Info<< nl << "Writing extendedFeatureEdgeMesh components to " << prefix
-        << endl;
-
-    label verti = 0;
-
-    edgeMesh::write(prefix + "_edgeMesh.obj");
-
-    OFstream convexFtPtStr(prefix + "_convexFeaturePts.obj");
-    Info<< "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");
-    Info<< "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");
-    Info<< "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");
-    Info<< "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");
-    Info<< "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");
-    Info<< "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");
-    Info<< "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");
-    Info<< "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");
-    Info<< "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");
-    Info<< "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;
-    }
-
-    OFstream edgeDirsStr(prefix + "_edgeDirections.obj");
-    Info<< "Writing edge directions to " << edgeDirsStr.name() << endl;
-
-    forAll(edgeDirections_, i)
-    {
-        const vector& eVec = edgeDirections_[i];
-        const edge& e = edges()[i];
-
-        meshTools::writeOBJ
-        (
-            edgeDirsStr,
-            points()[e.start()],
-            eVec + points()[e.start()]
-        );
-    }
+    is >> *this;
+    return !is.bad();
 }
 
 
 bool Foam::extendedFeatureEdgeMesh::writeData(Ostream& os) const
 {
-    os  << "// points" << nl
-        << points() << nl
-        << "// edges" << nl
-        << edges() << nl
-        << "// concaveStart mixedStart nonFeatureStart" << nl
-        << concaveStart_ << token::SPACE
-        << mixedStart_ << token::SPACE
-        << nonFeatureStart_ << nl
-        << "// internalStart flatStart openStart multipleStart" << nl
-        << internalStart_ << token::SPACE
-        << flatStart_ << token::SPACE
-        << openStart_ << token::SPACE
-        << multipleStart_ << nl
-        << "// normals" << nl
-        << normals_ << nl
-        << "// normal volume types" << nl
-        << normalVolumeTypes_ << nl
-        << "// normalDirections" << nl
-        << normalDirections_ << nl
-        << "// edgeNormals" << nl
-        << edgeNormals_ << nl
-        << "// featurePointNormals" << nl
-        << featurePointNormals_ << nl
-        << "// featurePointEdges" << nl
-        << featurePointEdges_ << nl
-        << "// regionEdges" << nl
-        << regionEdges_
-        << endl;
+    os << *this;
 
     return os.good();
 }
 
 
-Foam::Istream& Foam::operator>>
-(
-    Istream& is,
-    Foam::extendedFeatureEdgeMesh::sideVolumeType& vt
-)
-{
-    label type;
-    is  >> type;
 
-    vt = static_cast<Foam::extendedFeatureEdgeMesh::sideVolumeType>(type);
-
-    // Check state of Istream
-    is.check("operator>>(Istream&, sideVolumeType&)");
-
-    return is;
-}
-
-
-Foam::Ostream& Foam::operator<<
-(
-    Ostream& os,
-    const Foam::extendedFeatureEdgeMesh::sideVolumeType& vt
-)
-{
-    os  << static_cast<label>(vt);
-
-    return os;
-}
+//bool Foam::extendedFeatureEdgeMesh::writeData(Ostream& os) const
+//{
+//    os  << "// points" << nl
+//        << points() << nl
+//        << "// edges" << nl
+//        << edges() << nl
+//        << "// concaveStart mixedStart nonFeatureStart" << nl
+//        << concaveStart_ << token::SPACE
+//        << mixedStart_ << token::SPACE
+//        << nonFeatureStart_ << nl
+//        << "// internalStart flatStart openStart multipleStart" << nl
+//        << internalStart_ << token::SPACE
+//        << flatStart_ << token::SPACE
+//        << openStart_ << token::SPACE
+//        << multipleStart_ << nl
+//        << "// normals" << nl
+//        << normals_ << nl
+//        << "// normal volume types" << nl
+//        << normalVolumeTypes_ << nl
+//        << "// normalDirections" << nl
+//        << normalDirections_ << nl
+//        << "// edgeNormals" << nl
+//        << edgeNormals_ << nl
+//        << "// featurePointNormals" << nl
+//        << featurePointNormals_ << nl
+//        << "// featurePointEdges" << nl
+//        << featurePointEdges_ << nl
+//        << "// regionEdges" << nl
+//        << regionEdges_
+//        << endl;
+//
+//    return os.good();
+//}
 
+//
+//Foam::Istream& Foam::operator>>
+//(
+//    Istream& is,
+//    Foam::extendedFeatureEdgeMesh::sideVolumeType& vt
+//)
+//{
+//    label type;
+//    is  >> type;
+//
+//    vt = static_cast<Foam::extendedFeatureEdgeMesh::sideVolumeType>(type);
+//
+//    // Check state of Istream
+//    is.check("operator>>(Istream&, sideVolumeType&)");
+//
+//    return is;
+//}
+//
+//
+//Foam::Ostream& Foam::operator<<
+//(
+//    Ostream& os,
+//    const Foam::extendedFeatureEdgeMesh::sideVolumeType& vt
+//)
+//{
+//    os  << static_cast<label>(vt);
+//
+//    return os;
+//}
+//
 
 // ************************************************************************* //
diff --git a/src/edgeMesh/extendedEdgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.H b/src/edgeMesh/extendedEdgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.H
index fdeeeb8dd0..3026f15b4e 100644
--- a/src/edgeMesh/extendedEdgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.H
+++ b/src/edgeMesh/extendedEdgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.H
@@ -26,28 +26,9 @@ Class
 
 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_ .. mixedStart_-1      : concave points
-        mixedStart_ .. nonFeatureStart_-1   : 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.
+    extendedEdgeMesh + IO.
 
 SourceFiles
-    extendedFeatureEdgeMeshI.H
     extendedFeatureEdgeMesh.C
 
 \*---------------------------------------------------------------------------*/
@@ -55,19 +36,14 @@ SourceFiles
 #ifndef extendedFeatureEdgeMesh_H
 #define extendedFeatureEdgeMesh_H
 
-#include "edgeMesh.H"
+#include "extendedEdgeMesh.H"
 #include "regIOobject.H"
-#include "indexedOctree.H"
-#include "treeDataEdge.H"
-#include "treeDataPoint.H"
-#include "PrimitivePatch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-class surfaceFeatures;
 class objectRegistry;
 
 /*---------------------------------------------------------------------------*\
@@ -77,7 +53,7 @@ class objectRegistry;
 class extendedFeatureEdgeMesh
 :
     public regIOobject,
-    public edgeMesh
+    public extendedEdgeMesh
 {
 
 public:
@@ -85,142 +61,6 @@ 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
-    };
-
-    static const Foam::NamedEnum<pointStatus, 4> pointStatusNames_;
-
-    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)
-    };
-
-    static const Foam::NamedEnum<edgeStatus, 6> edgeStatusNames_;
-
-    //- Normals point to the outside
-    enum sideVolumeType
-    {
-        INSIDE  = 0,  // mesh inside
-        OUTSIDE = 1,  // mesh outside
-        BOTH    = 2,  // e.g. a baffle
-        NEITHER = 3   // not sure when this may be used
-    };
-
-    static const Foam::NamedEnum<sideVolumeType, 4> sideVolumeTypeNames_;
-
-    //- Angular closeness tolerance for treating normals as the same
-    static scalar cosNormalAngleTol_;
-
-
-private:
-
-    // Static data
-
-        //- 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_;
-
-        //- Type per normal: which side of normal to mesh
-        List<sideVolumeType> normalVolumeTypes_;
-
-        //- Flat and open edges require the direction of the edge
-        vectorField edgeDirections_;
-
-        //- Starting directions for the edges.
-        //  This vector points to the half of the plane defined by the first
-        //  edge normal.
-        labelListList normalDirections_;
-
-        //- Indices of the normals that are adjacent to the feature edges
-        labelListList edgeNormals_;
-
-        //- Indices of the normals that are adjacent to the feature points
-        //  (only valid for 0..nonFeatureStart_-1)
-        labelListList featurePointNormals_;
-
-        //- Indices of feature edges attached to feature points. The edges are
-        //  ordered so that they can be circulated.
-        labelListList featurePointEdges_;
-
-        //- Feature edges which are on the boundary between regions
-        labelList regionEdges_;
-
-        //- Search tree for all feature points
-        mutable autoPtr<indexedOctree<treeDataPoint> > pointTree_;
-
-        //- 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;
-
-        template<class Patch>
-        void sortPointsAndEdges
-        (
-            const Patch&,
-            const labelList& featureEdges,
-            const labelList& regionFeatureEdges,
-            const labelList& feaurePoints
-        );
-
-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
 
@@ -228,18 +68,10 @@ public:
         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>&
+            const extendedEdgeMesh&
         );
 
         //- Construct given a surface with selected edges,point
@@ -289,190 +121,20 @@ public:
 
 
     //- Destructor
-    ~extendedFeatureEdgeMesh();
+    virtual ~extendedFeatureEdgeMesh();
 
 
-    // Member Functions
+    // IO
 
-        // Find
+        //- Give precedence to the regIOobject write
+        using regIOobject::write;
 
-            //- Find nearest surface edge for the sample point.
-            void nearestFeaturePoint
-            (
-                const point& sample,
-                scalar searchDistSqr,
-                pointIndexHit& info
-            ) const;
+        //- ReadData function required for regIOobject read operation
+        virtual bool readData(Istream&);
 
-            //- Find nearest surface edge for the sample point.
-            void nearestFeatureEdge
-            (
-                const point& sample,
-                scalar searchDistSqr,
-                pointIndexHit& info
-            ) const;
+        //- WriteData function required for regIOobject write operation
+        virtual bool writeData(Ostream&) 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;
-
-            //- Find all the feature points within searchDistSqr of sample
-            void allNearestFeaturePoints
-            (
-                const point& sample,
-                scalar searchRadiusSqr,
-                List<pointIndexHit>& info
-            ) const;
-
-            //- Find all the feature edges within searchDistSqr of sample
-            void allNearestFeatureEdges
-            (
-                const point& sample,
-                const scalar searchRadiusSqr,
-                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
-            inline const List<sideVolumeType>& normalVolumeTypes() const;
-
-            //- Return the edgeDirection vectors
-            inline const vectorField& edgeDirections() const;
-
-            //-
-            inline const labelListList& normalDirections() 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 edge labels for a given feature point. Edges are
-            //  ordered by the faces that they share. The edge labels
-            //  correspond to the entry in edges().
-            inline const labelListList& featurePointEdges() 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;
-
-            //- Return the baffle faces of a specified edge
-            inline PackedList<2> edgeBaffles(label edgeI) const;
-
-            //- Demand driven construction of octree for feature points
-            const indexedOctree<treeDataPoint>& pointTree() 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;
-
-
-        // Edit
-
-            //- Add extendedFeatureEdgeMesh. No filtering of duplicates.
-            void add(const extendedFeatureEdgeMesh&);
-
-            //- Flip normals. All concave become convex, all internal external
-            //  etc.
-            void flipNormals();
-
-
-        // 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;
-
-            friend Istream& operator>>(Istream& is, sideVolumeType& vt);
-            friend Ostream& operator<<(Ostream& os, const sideVolumeType& vt);
-
-
-        //- Classify the type of feature edge.  Requires face centre 0 to face
-        //  centre 1 vector to distinguish internal from external
-        static edgeStatus classifyEdge
-        (
-            const List<vector>& norms,
-            const labelList& edNorms,
-            const vector& fC0tofC1
-        );
 };
 
 
@@ -482,14 +144,6 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#include "extendedFeatureEdgeMeshI.H"
-
-#ifdef NoRepository
-#   include "extendedFeatureEdgeMeshTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #endif
 
 // ************************************************************************* //