diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
index b0d3d83ec0..6d83d18363 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,6 +36,7 @@ Description
 #include "Time.H"
 #include "surfaceFeatures.H"
 #include "featureEdgeMesh.H"
+#include "extendedFeatureEdgeMesh.H"
 #include "treeBoundBox.H"
 #include "meshTools.H"
 #include "OFstream.H"
@@ -298,20 +299,46 @@ int main(int argc, char *argv[])
 
     // Extracting and writing a featureEdgeMesh
 
-    Pout<< nl << "Writing featureEdgeMesh to constant/featureEdgeMesh."
-        << endl;
-
-    featureEdgeMesh feMesh
+    extendedFeatureEdgeMesh feMesh
     (
         newSet,
         runTime,
         surfFileName.lessExt().name() + ".featureEdgeMesh"
     );
 
+    Info<< nl << "Writing extendedFeatureEdgeMesh to " << feMesh.objectPath()
+        << endl;
+
+
     feMesh.writeObj(surfFileName.lessExt().name());
 
     feMesh.write();
 
+
+    // Write a featureEdgeMesh for backwards compatibility
+    {
+        featureEdgeMesh bfeMesh
+        (
+            IOobject
+            (
+                surfFileName.lessExt().name() + ".eMesh",   // name
+                runTime.constant(), // instance
+                "triSurface",
+                runTime,            // registry
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            feMesh.points(),
+            feMesh.edges()
+        );
+
+        Info<< nl << "Writing featureEdgeMesh to "
+            << bfeMesh.objectPath() << endl;
+
+        bfeMesh.regIOobject::write();
+    }
+
     Info<< "End\n" << endl;
 
     return 0;
diff --git a/src/edgeMesh/Make/files b/src/edgeMesh/Make/files
index cdad239999..52f985de8c 100644
--- a/src/edgeMesh/Make/files
+++ b/src/edgeMesh/Make/files
@@ -20,7 +20,10 @@ $(edgeFormats)/starcd/STARCDedgeFormatRunTime.C
 $(edgeFormats)/vtk/VTKedgeFormat.C
 $(edgeFormats)/vtk/VTKedgeFormatRunTime.C
 
+extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C
+
 featureEdgeMesh/featureEdgeMesh.C
 
 
+
 LIB = $(FOAM_LIBBIN)/libedgeMesh
diff --git a/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C
new file mode 100644
index 0000000000..c4f536e209
--- /dev/null
+++ b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C
@@ -0,0 +1,1028 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "extendedFeatureEdgeMesh.H"
+#include "triSurface.H"
+#include "Random.H"
+#include "Time.H"
+#include "meshTools.H"
+#include "linePointRef.H"
+#include "ListListOps.H"
+#include "OFstream.H"
+#include "IFstream.H"
+#include "unitConversion.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(Foam::extendedFeatureEdgeMesh, 0);
+
+Foam::scalar Foam::extendedFeatureEdgeMesh::cosNormalAngleTol_ =
+    Foam::cos(degToRad(0.1));
+
+
+Foam::label Foam::extendedFeatureEdgeMesh::convexStart_ = 0;
+
+
+Foam::label Foam::extendedFeatureEdgeMesh::externalStart_ = 0;
+
+
+Foam::label Foam::extendedFeatureEdgeMesh::nPointTypes = 4;
+
+
+Foam::label Foam::extendedFeatureEdgeMesh::nEdgeTypes = 5;
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh(const IOobject& io)
+:
+    regIOobject(io),
+    edgeMesh(pointField(0), edgeList(0)),
+    concaveStart_(0),
+    mixedStart_(0),
+    nonFeatureStart_(0),
+    internalStart_(0),
+    flatStart_(0),
+    openStart_(0),
+    multipleStart_(0),
+    normals_(0),
+    edgeDirections_(0),
+    edgeNormals_(0),
+    featurePointNormals_(0),
+    regionEdges_(0),
+    edgeTree_(),
+    edgeTreesByType_()
+{
+    if
+    (
+        io.readOpt() == IOobject::MUST_READ
+     || io.readOpt() == IOobject::MUST_READ_IF_MODIFIED
+     || (io.readOpt() == IOobject::READ_IF_PRESENT && headerOk())
+    )
+    {
+        if (readOpt() == IOobject::MUST_READ_IF_MODIFIED)
+        {
+            WarningIn
+            (
+                "extendedFeatureEdgeMesh::extendedFeatureEdgeMesh"
+                "(const IOobject&)"
+            )   << "Specified IOobject::MUST_READ_IF_MODIFIED but class"
+                << " does not support automatic rereading."
+                << endl;
+        }
+
+        Istream& is = readStream(typeName);
+
+        is  >> *this
+            >> concaveStart_
+            >> mixedStart_
+            >> nonFeatureStart_
+            >> internalStart_
+            >> flatStart_
+            >> openStart_
+            >> multipleStart_
+            >> normals_
+            >> edgeNormals_
+            >> featurePointNormals_
+            >> regionEdges_;
+
+        close();
+
+        {
+            // Calculate edgeDirections
+
+            const edgeList& eds(edges());
+
+            const pointField& pts(points());
+
+            edgeDirections_.setSize(eds.size());
+
+            forAll(eds, eI)
+            {
+                edgeDirections_[eI] = eds[eI].vec(pts);
+            }
+
+            edgeDirections_ /= mag(edgeDirections_);
+        }
+    }
+
+    if (debug)
+    {
+        Pout<< "extendedFeatureEdgeMesh::extendedFeatureEdgeMesh :"
+            << " constructed from IOobject :"
+            << " points:" << points().size()
+            << " edges:" << edges().size()
+            << endl;
+    }
+}
+
+
+Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh
+(
+    const IOobject& io,
+    const extendedFeatureEdgeMesh& fem
+)
+:
+    regIOobject(io),
+    edgeMesh(fem),
+    concaveStart_(fem.concaveStart()),
+    mixedStart_(fem.mixedStart()),
+    nonFeatureStart_(fem.nonFeatureStart()),
+    internalStart_(fem.internalStart()),
+    flatStart_(fem.flatStart()),
+    openStart_(fem.openStart()),
+    multipleStart_(fem.multipleStart()),
+    normals_(fem.normals()),
+    edgeDirections_(fem.edgeDirections()),
+    edgeNormals_(fem.edgeNormals()),
+    featurePointNormals_(fem.featurePointNormals()),
+    regionEdges_(fem.regionEdges()),
+    edgeTree_(),
+    edgeTreesByType_()
+{}
+
+
+Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh
+(
+    const IOobject& io,
+    const Xfer<pointField>& pointLst,
+    const Xfer<edgeList>& edgeLst
+)
+:
+    regIOobject(io),
+    edgeMesh(pointLst, edgeLst),
+    concaveStart_(0),
+    mixedStart_(0),
+    nonFeatureStart_(0),
+    internalStart_(0),
+    flatStart_(0),
+    openStart_(0),
+    multipleStart_(0),
+    normals_(0),
+    edgeDirections_(0),
+    edgeNormals_(0),
+    featurePointNormals_(0),
+    regionEdges_(0),
+    edgeTree_(),
+    edgeTreesByType_()
+{}
+
+
+Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh
+(
+    const surfaceFeatures& sFeat,
+    const objectRegistry& obr,
+    const fileName& sFeatFileName
+)
+:
+    regIOobject
+    (
+        IOobject
+        (
+            sFeatFileName,
+            obr.time().constant(),
+            "extendedFeatureEdgeMesh",
+            obr,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        )
+    ),
+    edgeMesh(pointField(0), edgeList(0)),
+    concaveStart_(-1),
+    mixedStart_(-1),
+    nonFeatureStart_(-1),
+    internalStart_(-1),
+    flatStart_(-1),
+    openStart_(-1),
+    multipleStart_(-1),
+    normals_(0),
+    edgeDirections_(0),
+    edgeNormals_(0),
+    featurePointNormals_(0),
+    regionEdges_(0),
+    edgeTree_(),
+    edgeTreesByType_()
+{
+    // Extract and reorder the data from surfaceFeatures
+
+    // References to the surfaceFeatures data
+    const triSurface& surf(sFeat.surface());
+    const pointField& sFeatLocalPts(surf.localPoints());
+    const edgeList& sFeatEds(surf.edges());
+
+    // Filling the extendedFeatureEdgeMesh with the raw geometrical data.
+
+    label nFeatEds = sFeat.featureEdges().size();
+
+    DynamicList<point> tmpPts;
+    edgeList eds(nFeatEds);
+    DynamicList<vector> norms;
+    vectorField edgeDirections(nFeatEds);
+    labelListList edgeNormals(nFeatEds);
+    DynamicList<label> regionEdges;
+
+    // Mapping between old and new indices, there is entry in the map for each
+    // of sFeatLocalPts, -1 means that this point hasn't been used (yet), >= 0
+    // corresponds to the index
+    labelList pointMap(sFeatLocalPts.size(), -1);
+
+    // Noting when the normal of a face has been used so not to duplicate
+    labelList faceMap(surf.size(), -1);
+
+    // Collecting the status of edge for subsequent sorting
+    List<edgeStatus> edStatus(nFeatEds, NONE);
+
+    forAll(sFeat.featurePoints(), i)
+    {
+        label sFPI = sFeat.featurePoints()[i];
+
+        tmpPts.append(sFeatLocalPts[sFPI]);
+
+        pointMap[sFPI] = tmpPts.size() - 1;
+    }
+
+    // All feature points have been added
+    nonFeatureStart_ = tmpPts.size();
+
+    forAll(sFeat.featureEdges(), i)
+    {
+        label sFEI = sFeat.featureEdges()[i];
+
+        const edge& fE(sFeatEds[sFEI]);
+
+        // Check to see if the points have been already used
+        if (pointMap[fE.start()] == -1)
+        {
+             tmpPts.append(sFeatLocalPts[fE.start()]);
+
+             pointMap[fE.start()] = tmpPts.size() - 1;
+        }
+
+        eds[i].start() = pointMap[fE.start()];
+
+        if (pointMap[fE.end()] == -1)
+        {
+            tmpPts.append(sFeatLocalPts[fE.end()]);
+
+            pointMap[fE.end()] = tmpPts.size() - 1;
+        }
+
+        eds[i].end() = pointMap[fE.end()];
+
+        // Pick up the faces adjacent to the feature edge
+        const labelList& eFaces = surf.edgeFaces()[sFEI];
+
+        edgeNormals[i].setSize(eFaces.size());
+
+        forAll(eFaces, j)
+        {
+            label eFI = eFaces[j];
+
+            // Check to see if the points have been already used
+            if (faceMap[eFI] == -1)
+            {
+                norms.append(surf.faceNormals()[eFI]);
+
+                faceMap[eFI] = norms.size() - 1;
+            }
+
+            edgeNormals[i][j] = faceMap[eFI];
+        }
+
+        vector fC0tofC1(vector::zero);
+
+        if (eFaces.size() == 2)
+        {
+            fC0tofC1 =
+                surf[eFaces[1]].centre(surf.points())
+              - surf[eFaces[0]].centre(surf.points());
+        }
+
+        edStatus[i] = classifyEdge(norms, edgeNormals[i], fC0tofC1);
+
+        edgeDirections[i] = fE.vec(sFeatLocalPts);
+
+        if (i < sFeat.nRegionEdges())
+        {
+            regionEdges.append(i);
+        }
+    }
+
+    // Reorder the edges by classification
+
+    List<DynamicList<label> > allEds(nEdgeTypes);
+
+    DynamicList<label>& externalEds(allEds[0]);
+    DynamicList<label>& internalEds(allEds[1]);
+    DynamicList<label>& flatEds(allEds[2]);
+    DynamicList<label>& openEds(allEds[3]);
+    DynamicList<label>& multipleEds(allEds[4]);
+
+    forAll(eds, i)
+    {
+        edgeStatus eStat = edStatus[i];
+
+        if (eStat == EXTERNAL)
+        {
+            externalEds.append(i);
+        }
+        else if (eStat == INTERNAL)
+        {
+            internalEds.append(i);
+        }
+        else if (eStat == FLAT)
+        {
+            flatEds.append(i);
+        }
+        else if (eStat == OPEN)
+        {
+            openEds.append(i);
+        }
+        else if (eStat == MULTIPLE)
+        {
+            multipleEds.append(i);
+        }
+        else if (eStat == NONE)
+        {
+            FatalErrorIn
+            (
+                "Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh"
+                "("
+                "    const surfaceFeatures& sFeat,"
+                "    const objectRegistry& obr,"
+                "    const fileName& sFeatFileName"
+                ")"
+            )
+                << nl << "classifyEdge returned NONE on edge "
+                << eds[i]
+                << ". There is a problem with definition of this edge."
+                << nl << abort(FatalError);
+        }
+    }
+
+    internalStart_ = externalEds.size();
+    flatStart_ = internalStart_ + internalEds.size();
+    openStart_ = flatStart_ + flatEds.size();
+    multipleStart_ = openStart_ + openEds.size();
+
+    labelList edMap
+    (
+        ListListOps::combine<labelList>
+        (
+            allEds,
+            accessOp<labelList>()
+        )
+    );
+
+    edMap = invert(edMap.size(), edMap);
+
+    inplaceReorder(edMap, eds);
+    inplaceReorder(edMap, edStatus);
+    inplaceReorder(edMap, edgeDirections);
+    inplaceReorder(edMap, edgeNormals);
+    inplaceRenumber(edMap, regionEdges);
+
+    pointField pts(tmpPts);
+
+    // Initialise the edgeMesh
+    edgeMesh::operator=(edgeMesh(pts, eds));
+
+    // Initialise sorted edge related data
+    edgeDirections_ = edgeDirections/mag(edgeDirections);
+    edgeNormals_ = edgeNormals;
+    regionEdges_ = regionEdges;
+
+    // Normals are all now found and indirectly addressed, can also be stored
+    normals_ = vectorField(norms);
+
+    // Reorder the feature points by classification
+
+    List<DynamicList<label> > allPts(3);
+
+    DynamicList<label>& convexPts(allPts[0]);
+    DynamicList<label>& concavePts(allPts[1]);
+    DynamicList<label>& mixedPts(allPts[2]);
+
+    for (label i = 0; i < nonFeatureStart_; i++)
+    {
+        pointStatus ptStatus = classifyFeaturePoint(i);
+
+        if (ptStatus == CONVEX)
+        {
+            convexPts.append(i);
+        }
+        else if (ptStatus == CONCAVE)
+        {
+            concavePts.append(i);
+        }
+        else if (ptStatus == MIXED)
+        {
+            mixedPts.append(i);
+        }
+        else if (ptStatus == NONFEATURE)
+        {
+            FatalErrorIn
+            (
+                "Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh"
+                "("
+                "    const surfaceFeatures& sFeat,"
+                "    const objectRegistry& obr,"
+                "    const fileName& sFeatFileName"
+                ")"
+            )
+                << nl << "classifyFeaturePoint returned NONFEATURE on point at "
+                << points()[i]
+                << ". There is a problem with definition of this feature point."
+                << nl << abort(FatalError);
+        }
+    }
+
+    concaveStart_ = convexPts.size();
+    mixedStart_ = concaveStart_ + concavePts.size();
+
+    labelList ftPtMap
+    (
+        ListListOps::combine<labelList>
+        (
+            allPts,
+            accessOp<labelList>()
+        )
+    );
+
+    ftPtMap = invert(ftPtMap.size(), ftPtMap);
+
+    // Creating the ptMap from the ftPtMap with identity values up to the size
+    // of pts to create an oldToNew map for inplaceReorder
+
+    labelList ptMap(identity(pts.size()));
+
+    forAll(ftPtMap, i)
+    {
+        ptMap[i] = ftPtMap[i];
+    }
+
+    inplaceReorder(ptMap, pts);
+
+    forAll(eds, i)
+    {
+        inplaceRenumber(ptMap, eds[i]);
+    }
+
+    // Reinitialise the edgeMesh with sorted feature points and
+    // renumbered edges
+    edgeMesh::operator=(edgeMesh(pts, eds));
+
+    // Generate the featurePointNormals
+
+    labelListList featurePointNormals(nonFeatureStart_);
+
+    for (label i = 0; i < nonFeatureStart_; i++)
+    {
+        DynamicList<label> tmpFtPtNorms;
+
+        const labelList& ptEds = pointEdges()[i];
+
+        forAll(ptEds, j)
+        {
+            const labelList& ptEdNorms(edgeNormals[ptEds[j]]);
+
+            forAll(ptEdNorms, k)
+            {
+                if (findIndex(tmpFtPtNorms, ptEdNorms[k]) == -1)
+                {
+                    bool addNormal = true;
+
+                    // Check that the normal direction is unique at this feature
+                    forAll(tmpFtPtNorms, q)
+                    {
+                        if
+                        (
+                            (normals_[ptEdNorms[k]] & normals_[tmpFtPtNorms[q]])
+                          > cosNormalAngleTol_
+                        )
+                        {
+                            // Parallel to an existing normal, do not add
+                            addNormal = false;
+
+                            break;
+                        }
+                    }
+
+                    if (addNormal)
+                    {
+                        tmpFtPtNorms.append(ptEdNorms[k]);
+                    }
+                }
+            }
+        }
+
+        featurePointNormals[i] = tmpFtPtNorms;
+    }
+
+    featurePointNormals_ = featurePointNormals;
+}
+
+
+Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh
+(
+    const IOobject& io,
+    const pointField& pts,
+    const edgeList& eds,
+    label concaveStart,
+    label mixedStart,
+    label nonFeatureStart,
+    label internalStart,
+    label flatStart,
+    label openStart,
+    label multipleStart,
+    const vectorField& normals,
+    const vectorField& edgeDirections,
+    const labelListList& edgeNormals,
+    const labelListList& featurePointNormals,
+    const labelList& regionEdges
+)
+:
+    regIOobject(io),
+    edgeMesh(pts, eds),
+    concaveStart_(concaveStart),
+    mixedStart_(mixedStart),
+    nonFeatureStart_(nonFeatureStart),
+    internalStart_(internalStart),
+    flatStart_(flatStart),
+    openStart_(openStart),
+    multipleStart_(multipleStart),
+    normals_(normals),
+    edgeDirections_(edgeDirections),
+    edgeNormals_(edgeNormals),
+    featurePointNormals_(featurePointNormals),
+    regionEdges_(regionEdges),
+    edgeTree_(),
+    edgeTreesByType_()
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::extendedFeatureEdgeMesh::~extendedFeatureEdgeMesh()
+{}
+
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+Foam::extendedFeatureEdgeMesh::pointStatus
+Foam::extendedFeatureEdgeMesh::classifyFeaturePoint
+(
+    label ptI
+) const
+{
+    labelList ptEds(pointEdges()[ptI]);
+
+    label nPtEds = ptEds.size();
+    label nExternal = 0;
+    label nInternal = 0;
+
+    if (nPtEds == 0)
+    {
+        // There are no edges attached to the point, this is a problem
+        return NONFEATURE;
+    }
+
+    forAll(ptEds, i)
+    {
+        edgeStatus edStat = getEdgeStatus(ptEds[i]);
+
+        if (edStat == EXTERNAL)
+        {
+            nExternal++;
+        }
+        else if (edStat == INTERNAL)
+        {
+            nInternal++;
+        }
+    }
+
+    if (nExternal == nPtEds)
+    {
+        return CONVEX;
+    }
+    else if (nInternal == nPtEds)
+    {
+        return CONCAVE;
+    }
+    else
+    {
+        return MIXED;
+    }
+}
+
+
+Foam::extendedFeatureEdgeMesh::edgeStatus
+Foam::extendedFeatureEdgeMesh::classifyEdge
+(
+    const List<vector>& norms,
+    const labelList& edNorms,
+    const vector& fC0tofC1
+) const
+{
+    label nEdNorms = edNorms.size();
+
+    if (nEdNorms == 1)
+    {
+        return OPEN;
+    }
+    else if (nEdNorms == 2)
+    {
+        const vector n0(norms[edNorms[0]]);
+        const vector n1(norms[edNorms[1]]);
+
+        if ((n0 & n1) > cosNormalAngleTol_)
+        {
+            return FLAT;
+        }
+        else if ((fC0tofC1 & n0) > 0.0)
+        {
+            return INTERNAL;
+        }
+        else
+        {
+            return EXTERNAL;
+        }
+    }
+    else if (nEdNorms > 2)
+    {
+        return MULTIPLE;
+    }
+    else
+    {
+        // There is a problem - the edge has no normals
+        return NONE;
+    }
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::extendedFeatureEdgeMesh::nearestFeatureEdge
+(
+    const point& sample,
+    scalar searchDistSqr,
+    pointIndexHit& info
+) const
+{
+    info = edgeTree().findNearest
+    (
+        sample,
+        searchDistSqr
+    );
+}
+
+
+void Foam::extendedFeatureEdgeMesh::nearestFeatureEdge
+(
+    const pointField& samples,
+    const scalarField& searchDistSqr,
+    List<pointIndexHit>& info
+) const
+{
+    info.setSize(samples.size());
+
+    forAll(samples, i)
+    {
+        nearestFeatureEdge
+        (
+            samples[i],
+            searchDistSqr[i],
+            info[i]
+        );
+    }
+}
+
+
+void Foam::extendedFeatureEdgeMesh::nearestFeatureEdgeByType
+(
+    const point& sample,
+    const scalarField& searchDistSqr,
+    List<pointIndexHit>& info
+) const
+{
+    const PtrList<indexedOctree<treeDataEdge> >& edgeTrees = edgeTreesByType();
+
+    info.setSize(edgeTrees.size());
+
+    labelList sliceStarts(edgeTrees.size());
+
+    sliceStarts[0] = externalStart_;
+    sliceStarts[1] = internalStart_;
+    sliceStarts[2] = flatStart_;
+    sliceStarts[3] = openStart_;
+    sliceStarts[4] = multipleStart_;
+
+    forAll(edgeTrees, i)
+    {
+        info[i] = edgeTrees[i].findNearest
+        (
+            sample,
+            searchDistSqr[i]
+        );
+
+        // The index returned by the indexedOctree is local to the slice of
+        // edges it was supplied with, return the index to the value in the
+        // complete edge list
+
+        info[i].setIndex(info[i].index() + sliceStarts[i]);
+    }
+}
+
+
+const Foam::indexedOctree<Foam::treeDataEdge>&
+Foam::extendedFeatureEdgeMesh::edgeTree() const
+{
+    if (edgeTree_.empty())
+    {
+        Random rndGen(17301893);
+
+        // Slightly extended bb. Slightly off-centred just so on symmetric
+        // geometry there are less face/edge aligned items.
+        treeBoundBox bb
+        (
+            treeBoundBox(points()).extend(rndGen, 1E-4)
+        );
+
+        bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+        bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+
+        labelList allEdges(identity(edges().size()));
+
+        edgeTree_.reset
+        (
+            new indexedOctree<treeDataEdge>
+            (
+                treeDataEdge
+                (
+                    false,          // cachebb
+                    edges(),        // edges
+                    points(),       // points
+                    allEdges        // selected edges
+                ),
+                bb,     // bb
+                8,      // maxLevel
+                10,     // leafsize
+                3.0     // duplicity
+            )
+        );
+    }
+
+    return edgeTree_();
+}
+
+
+const Foam::PtrList<Foam::indexedOctree<Foam::treeDataEdge> >&
+Foam::extendedFeatureEdgeMesh::edgeTreesByType() const
+{
+    if (edgeTreesByType_.size() == 0)
+    {
+        edgeTreesByType_.setSize(nEdgeTypes);
+
+        Random rndGen(872141);
+
+        // Slightly extended bb. Slightly off-centred just so on symmetric
+        // geometry there are less face/edge aligned items.
+        treeBoundBox bb
+        (
+            treeBoundBox(points()).extend(rndGen, 1E-4)
+        );
+
+        bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+        bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+
+        labelListList sliceEdges(nEdgeTypes);
+
+        // External edges
+        sliceEdges[0] =
+            identity(internalStart_ - externalStart_) + externalStart_;
+
+        // Internal edges
+        sliceEdges[1] = identity(flatStart_ - internalStart_) + internalStart_;
+
+        // Flat edges
+        sliceEdges[2] = identity(openStart_ - flatStart_) + flatStart_;
+
+        // Open edges
+        sliceEdges[3] = identity(multipleStart_ - openStart_) + openStart_;
+
+        // Multiple edges
+        sliceEdges[4] =
+            identity(edges().size() - multipleStart_) + multipleStart_;
+
+        forAll(edgeTreesByType_, i)
+        {
+            edgeTreesByType_.set
+            (
+                i,
+                new indexedOctree<treeDataEdge>
+                (
+                    treeDataEdge
+                    (
+                        false,          // cachebb
+                        edges(),        // edges
+                        points(),       // points
+                        sliceEdges[i]   // selected edges
+                    ),
+                    bb,     // bb
+                    8,      // maxLevel
+                    10,     // leafsize
+                    3.0     // duplicity
+                )
+            );
+        }
+    }
+
+    return edgeTreesByType_;
+}
+
+
+void Foam::extendedFeatureEdgeMesh::writeObj
+(
+    const fileName& prefix
+) const
+{
+    Pout<< nl << "Writing extendedFeatureEdgeMesh components to " << prefix
+        << endl;
+
+    label verti = 0;
+
+    edgeMesh::write(prefix + "_edgeMesh.obj");
+
+    OFstream convexFtPtStr(prefix + "_convexFeaturePts.obj");
+    Pout<< "Writing convex feature points to " << convexFtPtStr.name() << endl;
+
+    for(label i = 0; i < concaveStart_; i++)
+    {
+        meshTools::writeOBJ(convexFtPtStr, points()[i]);
+    }
+
+    OFstream concaveFtPtStr(prefix + "_concaveFeaturePts.obj");
+    Pout<< "Writing concave feature points to "
+        << concaveFtPtStr.name() << endl;
+
+    for(label i = concaveStart_; i < mixedStart_; i++)
+    {
+        meshTools::writeOBJ(concaveFtPtStr, points()[i]);
+    }
+
+    OFstream mixedFtPtStr(prefix + "_mixedFeaturePts.obj");
+    Pout<< "Writing mixed feature points to " << mixedFtPtStr.name() << endl;
+
+    for(label i = mixedStart_; i < nonFeatureStart_; i++)
+    {
+        meshTools::writeOBJ(mixedFtPtStr, points()[i]);
+    }
+
+    OFstream mixedFtPtStructureStr(prefix + "_mixedFeaturePtsStructure.obj");
+    Pout<< "Writing mixed feature point structure to "
+        << mixedFtPtStructureStr.name() << endl;
+
+    verti = 0;
+    for(label i = mixedStart_; i < nonFeatureStart_; i++)
+    {
+        const labelList& ptEds = pointEdges()[i];
+
+        forAll(ptEds, j)
+        {
+            const edge& e = edges()[ptEds[j]];
+
+            meshTools::writeOBJ(mixedFtPtStructureStr, points()[e[0]]); verti++;
+            meshTools::writeOBJ(mixedFtPtStructureStr, points()[e[1]]); verti++;
+            mixedFtPtStructureStr << "l " << verti-1 << ' ' << verti << endl;
+        }
+    }
+
+    OFstream externalStr(prefix + "_externalEdges.obj");
+    Pout<< "Writing external edges to " << externalStr.name() << endl;
+
+    verti = 0;
+    for (label i = externalStart_; i < internalStart_; i++)
+    {
+        const edge& e = edges()[i];
+
+        meshTools::writeOBJ(externalStr, points()[e[0]]); verti++;
+        meshTools::writeOBJ(externalStr, points()[e[1]]); verti++;
+        externalStr << "l " << verti-1 << ' ' << verti << endl;
+    }
+
+    OFstream internalStr(prefix + "_internalEdges.obj");
+    Pout<< "Writing internal edges to " << internalStr.name() << endl;
+
+    verti = 0;
+    for (label i = internalStart_; i < flatStart_; i++)
+    {
+        const edge& e = edges()[i];
+
+        meshTools::writeOBJ(internalStr, points()[e[0]]); verti++;
+        meshTools::writeOBJ(internalStr, points()[e[1]]); verti++;
+        internalStr << "l " << verti-1 << ' ' << verti << endl;
+    }
+
+    OFstream flatStr(prefix + "_flatEdges.obj");
+    Pout<< "Writing flat edges to " << flatStr.name() << endl;
+
+    verti = 0;
+    for (label i = flatStart_; i < openStart_; i++)
+    {
+        const edge& e = edges()[i];
+
+        meshTools::writeOBJ(flatStr, points()[e[0]]); verti++;
+        meshTools::writeOBJ(flatStr, points()[e[1]]); verti++;
+        flatStr << "l " << verti-1 << ' ' << verti << endl;
+    }
+
+    OFstream openStr(prefix + "_openEdges.obj");
+    Pout<< "Writing open edges to " << openStr.name() << endl;
+
+    verti = 0;
+    for (label i = openStart_; i < multipleStart_; i++)
+    {
+        const edge& e = edges()[i];
+
+        meshTools::writeOBJ(openStr, points()[e[0]]); verti++;
+        meshTools::writeOBJ(openStr, points()[e[1]]); verti++;
+        openStr << "l " << verti-1 << ' ' << verti << endl;
+    }
+
+    OFstream multipleStr(prefix + "_multipleEdges.obj");
+    Pout<< "Writing multiple edges to " << multipleStr.name() << endl;
+
+    verti = 0;
+    for (label i = multipleStart_; i < edges().size(); i++)
+    {
+        const edge& e = edges()[i];
+
+        meshTools::writeOBJ(multipleStr, points()[e[0]]); verti++;
+        meshTools::writeOBJ(multipleStr, points()[e[1]]); verti++;
+        multipleStr << "l " << verti-1 << ' ' << verti << endl;
+    }
+
+    OFstream regionStr(prefix + "_regionEdges.obj");
+    Pout<< "Writing region edges to " << regionStr.name() << endl;
+
+    verti = 0;
+    forAll(regionEdges_, i)
+    {
+        const edge& e = edges()[regionEdges_[i]];
+
+        meshTools::writeOBJ(regionStr, points()[e[0]]); verti++;
+        meshTools::writeOBJ(regionStr, points()[e[1]]); verti++;
+        regionStr << "l " << verti-1 << ' ' << verti << endl;
+    }
+}
+
+
+bool Foam::extendedFeatureEdgeMesh::writeData(Ostream& os) const
+{
+    os  << "// points, edges, concaveStart, mixedStart, nonFeatureStart, " << nl
+        << "// internalStart, flatStart, openStart, multipleStart, " << nl
+        << "// normals, edgeNormals, featurePointNormals, regionEdges" << nl
+        << endl;
+
+    os  << points() << nl
+        << edges() << nl
+        << concaveStart_ << token::SPACE
+        << mixedStart_ << token::SPACE
+        << nonFeatureStart_ << token::SPACE
+        << internalStart_ << token::SPACE
+        << flatStart_ << token::SPACE
+        << openStart_ << token::SPACE
+        << multipleStart_ << nl
+        << normals_ << nl
+        << edgeNormals_ << nl
+        << featurePointNormals_ << nl
+        << regionEdges_
+        << endl;
+
+    return os.good();
+}
+
+
+// ************************************************************************* //
diff --git a/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.H b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.H
new file mode 100644
index 0000000000..f1e85867d4
--- /dev/null
+++ b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.H
@@ -0,0 +1,377 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::extendedFeatureEdgeMesh
+
+Description
+
+    Description of feature edges and points.
+
+    Feature points are a sorted subset at the start of the overall points list:
+        0 .. concaveStart_-1                : convex points (w.r.t normals)
+        concaveStart_-1 .. mixedStart_-1    : concave points
+        mixedStart_ .. nonFeatureStart_     : mixed internal/external points
+        nonFeatureStart_ .. size-1          : non-feature points
+
+    Feature edges are the edgeList of the edgeMesh and are sorted:
+        0 .. internalStart_-1           : external edges (convex w.r.t normals)
+        internalStart_ .. flatStart_-1  : internal edges (concave)
+        flatStart_ .. openStart_-1      : flat edges (neither concave or convex)
+                                          can arise from region interfaces on
+                                          flat surfaces
+        openStart_ .. multipleStart_-1  : open edges (e.g. from baffle surfaces)
+        multipleStart_ .. size-1        : multiply connected edges
+
+    The edge direction and feature edge and feature point adjacent normals
+    are stored.
+
+SourceFiles
+    extendedFeatureEdgeMeshI.H
+    extendedFeatureEdgeMesh.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef extendedFeatureEdgeMesh_H
+#define extendedFeatureEdgeMesh_H
+
+#include "edgeMesh.H"
+#include "surfaceFeatures.H"
+#include "objectRegistry.H"
+#include "IOdictionary.H"
+#include "indexedOctree.H"
+#include "treeDataEdge.H"
+#include "pointIndexHit.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class extendedFeatureEdgeMesh Declaration
+\*---------------------------------------------------------------------------*/
+
+class extendedFeatureEdgeMesh
+:
+    public regIOobject,
+    public edgeMesh
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("extendedFeatureEdgeMesh");
+
+    enum pointStatus
+    {
+        CONVEX,     // Fully convex point (w.r.t normals)
+        CONCAVE,    // Fully concave point
+        MIXED,      // A point surrounded by both convex and concave edges
+        NONFEATURE  // Not a feature point
+    };
+
+    enum edgeStatus
+    {
+        EXTERNAL,   // "Convex" edge
+        INTERNAL,   // "Concave" edge
+        FLAT,       // Neither concave or convex, on a flat surface
+        OPEN,       // i.e. only connected to one face
+        MULTIPLE,   // Multiply connected (connected to more than two faces)
+        NONE        // Not a classified feature edge (consistency with
+                    // surfaceFeatures)
+    };
+
+private:
+
+    // Static data
+
+        //- Angular closeness tolerance for treating normals as the same
+        static scalar cosNormalAngleTol_;
+
+        //- Index of the start of the convex feature points - static as 0
+        static label convexStart_;
+
+        //- Index of the start of the external feature edges - static as 0
+        static label externalStart_;
+
+
+    // Private data
+
+        //- Index of the start of the concave feature points
+        label concaveStart_;
+
+        //- Index of the start of the mixed type feature points
+        label mixedStart_;
+
+        //- Index of the start of the non-feature points
+        label nonFeatureStart_;
+
+        //- Index of the start of the internal feature edges
+        label internalStart_;
+
+        //- Index of the start of the flat feature edges
+        label flatStart_;
+
+        //- Index of the start of the open feature edges
+        label openStart_;
+
+        //- Index of the start of the multiply-connected feature edges
+        label multipleStart_;
+
+        //- Normals of the features, to be referred to by index by both feature
+        //  points and edges, unsorted
+        vectorField normals_;
+
+        //- Flat and open edges require the direction of the edge
+        vectorField edgeDirections_;
+
+        //- Indices of the normals that are adjacent to the feature edges
+        labelListList edgeNormals_;
+
+        //- Indices of the normals that are adjacent to the feature points
+        labelListList featurePointNormals_;
+
+        //- Feature edges which are on the boundary between regions
+        labelList regionEdges_;
+
+        //- Search tree for all edges
+        mutable autoPtr<indexedOctree<treeDataEdge> > edgeTree_;
+
+        //- Individual search trees for each type of edge
+        mutable PtrList<indexedOctree<treeDataEdge> > edgeTreesByType_;
+
+
+    // Private Member Functions
+
+        //- Classify the type of feature point.  Requires valid stored member
+        //  data for edges and normals.
+        pointStatus classifyFeaturePoint(label ptI) const;
+
+        //- Classify the type of feature edge.  Requires face centre 0 to face
+        //  centre 1 vector to distinguish internal from external
+        edgeStatus classifyEdge
+        (
+            const List<vector>& norms,
+            const labelList& edNorms,
+            const vector& fC0tofC1
+        ) const;
+
+
+public:
+
+    // Static data
+
+        //- Number of possible point types (i.e. number of slices)
+        static label nPointTypes;
+
+        //- Number of possible feature edge types (i.e. number of slices)
+        static label nEdgeTypes;
+
+    // Constructors
+
+        //- Construct (read) given an IOobject
+        extendedFeatureEdgeMesh(const IOobject&);
+
+        //- Construct as copy
+        explicit extendedFeatureEdgeMesh
+        (
+            const IOobject&,
+            const extendedFeatureEdgeMesh&
+        );
+
+        //- Construct by transferring components (points, edges)
+        extendedFeatureEdgeMesh
+        (
+            const IOobject&,
+            const Xfer<pointField>&,
+            const Xfer<edgeList>&
+        );
+
+        //- Construct (read) given surfaceFeatures, an objectRegistry and a
+        //  fileName to write to.  Extracts, classifies and reorders the data
+        //  from surfaceFeatures.
+        extendedFeatureEdgeMesh
+        (
+            const surfaceFeatures& sFeat,
+            const objectRegistry& obr,
+            const fileName& sFeatFileName
+        );
+
+        //- Construct from all components
+        extendedFeatureEdgeMesh
+        (
+            const IOobject& io,
+            const pointField& pts,
+            const edgeList& eds,
+            label concaveStart,
+            label mixedStart,
+            label nonFeatureStart,
+            label internalStart,
+            label flatStart,
+            label openStart,
+            label multipleStart,
+            const vectorField& normals,
+            const vectorField& edgeDirections,
+            const labelListList& edgeNormals,
+            const labelListList& featurePointNormals,
+            const labelList& regionEdges
+        );
+
+
+    //- Destructor
+    ~extendedFeatureEdgeMesh();
+
+
+    // Member Functions
+
+        // Find
+
+            //- Find nearest surface edge for the sample point.
+            void nearestFeatureEdge
+            (
+                const point& sample,
+                scalar searchDistSqr,
+                pointIndexHit& info
+            ) const;
+
+            //- Find nearest surface edge for each sample point.
+            void nearestFeatureEdge
+            (
+                const pointField& samples,
+                const scalarField& searchDistSqr,
+                List<pointIndexHit>& info
+            ) const;
+
+            //- Find the nearest point on each type of feature edge
+            void nearestFeatureEdgeByType
+            (
+                const point& sample,
+                const scalarField& searchDistSqr,
+                List<pointIndexHit>& info
+            ) const;
+
+        // Access
+
+            //- Return the index of the start of the convex feature points
+            inline label convexStart() const;
+
+            //- Return the index of the start of the concave feature points
+            inline label concaveStart() const;
+
+            //- Return the index of the start of the mixed type feature points
+            inline label mixedStart() const;
+
+            //- Return the index of the start of the non-feature points
+            inline label nonFeatureStart() const;
+
+            //- Return the index of the start of the external feature edges
+            inline label externalStart() const;
+
+            //- Return the index of the start of the internal feature edges
+            inline label internalStart() const;
+
+            //- Return the index of the start of the flat feature edges
+            inline label flatStart() const;
+
+            //- Return the index of the start of the open feature edges
+            inline label openStart() const;
+
+            //- Return the index of the start of the multiply-connected feature
+            //  edges
+            inline label multipleStart() const;
+
+            //- Return whether or not the point index is a feature point
+            inline bool featurePoint(label ptI) const;
+
+            //- Return the normals of the surfaces adjacent to the feature edges
+            //  and points
+            inline const vectorField& normals() const;
+
+            //- Return the edgeDirection vectors
+            inline const vectorField& edgeDirections() const;
+
+            //- Return the direction of edgeI, pointing away from ptI
+            inline vector edgeDirection(label edgeI, label ptI) const;
+
+            //- Return the indices of the normals that are adjacent to the
+            //  feature edges
+            inline const labelListList& edgeNormals() const;
+
+            //- Return the normal vectors for a given set of normal indices
+            inline vectorField edgeNormals(const labelList& edgeNormIs) const;
+
+            //- Return the normal vectors for a given edge
+            inline vectorField edgeNormals(label edgeI) const;
+
+            //- Return the indices of the normals that are adjacent to the
+            //  feature points
+            inline const labelListList& featurePointNormals() const;
+
+            //- Return the normal vectors for a given feature point
+            inline vectorField featurePointNormals(label ptI) const;
+
+            //- Return the feature edges which are on the boundary between
+            //  regions
+            inline const labelList& regionEdges() const;
+
+            //- Return the pointStatus of a specified point
+            inline pointStatus getPointStatus(label ptI) const;
+
+            //- Return the edgeStatus of a specified edge
+            inline edgeStatus getEdgeStatus(label edgeI) const;
+
+            //- Demand driven construction of octree for boundary edges
+            const indexedOctree<treeDataEdge>& edgeTree() const;
+
+            //- Demand driven construction of octree for boundary edges by type
+            const PtrList<indexedOctree<treeDataEdge> >&
+            edgeTreesByType() const;
+
+
+        // Write
+
+            //- Write all components of the extendedFeatureEdgeMesh as obj files
+            void writeObj(const fileName& prefix) const;
+
+            //- Give precedence to the regIOobject write
+            using regIOobject::write;
+
+            //- WriteData function required for regIOobject write operation
+            virtual bool writeData(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "extendedFeatureEdgeMeshI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/edgeMesh/featureEdgeMesh/featureEdgeMeshI.H b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMeshI.H
similarity index 66%
rename from src/edgeMesh/featureEdgeMesh/featureEdgeMeshI.H
rename to src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMeshI.H
index 6aac568b03..4feb72e686 100644
--- a/src/edgeMesh/featureEdgeMesh/featureEdgeMeshI.H
+++ b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMeshI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2009-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2009-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,78 +25,79 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline Foam::label Foam::featureEdgeMesh::convexStart() const
+inline Foam::label Foam::extendedFeatureEdgeMesh::convexStart() const
 {
     return convexStart_;
 }
 
 
-inline Foam::label Foam::featureEdgeMesh::concaveStart() const
+inline Foam::label Foam::extendedFeatureEdgeMesh::concaveStart() const
 {
     return concaveStart_;
 }
 
 
-inline Foam::label Foam::featureEdgeMesh::mixedStart() const
+inline Foam::label Foam::extendedFeatureEdgeMesh::mixedStart() const
 {
     return mixedStart_;
 }
 
 
-inline Foam::label Foam::featureEdgeMesh::nonFeatureStart() const
+inline Foam::label Foam::extendedFeatureEdgeMesh::nonFeatureStart() const
 {
     return nonFeatureStart_;
 }
 
 
-inline Foam::label Foam::featureEdgeMesh::externalStart() const
+inline Foam::label Foam::extendedFeatureEdgeMesh::externalStart() const
 {
     return externalStart_;
 }
 
 
-inline Foam::label Foam::featureEdgeMesh::internalStart() const
+inline Foam::label Foam::extendedFeatureEdgeMesh::internalStart() const
 {
     return internalStart_;
 }
 
 
-inline Foam::label Foam::featureEdgeMesh::flatStart() const
+inline Foam::label Foam::extendedFeatureEdgeMesh::flatStart() const
 {
     return flatStart_;
 }
 
 
-inline Foam::label Foam::featureEdgeMesh::openStart() const
+inline Foam::label Foam::extendedFeatureEdgeMesh::openStart() const
 {
     return openStart_;
 }
 
 
-inline Foam::label Foam::featureEdgeMesh::multipleStart() const
+inline Foam::label Foam::extendedFeatureEdgeMesh::multipleStart() const
 {
     return multipleStart_;
 }
 
 
-inline bool Foam::featureEdgeMesh::featurePoint(label ptI) const
+inline bool Foam::extendedFeatureEdgeMesh::featurePoint(label ptI) const
 {
     return ptI < nonFeatureStart_;
 }
 
 
-inline const Foam::vectorField& Foam::featureEdgeMesh::normals() const
+inline const Foam::vectorField& Foam::extendedFeatureEdgeMesh::normals() const
 {
     return normals_;
 }
 
-inline const Foam::vectorField& Foam::featureEdgeMesh::edgeDirections() const
+inline const Foam::vectorField& Foam::extendedFeatureEdgeMesh::edgeDirections()
+const
 {
     return edgeDirections_;
 }
 
 
-inline Foam::vector Foam::featureEdgeMesh::edgeDirection
+inline Foam::vector Foam::extendedFeatureEdgeMesh::edgeDirection
 (
     label edgeI,
     label ptI
@@ -114,7 +115,7 @@ inline Foam::vector Foam::featureEdgeMesh::edgeDirection
     }
     else
     {
-        FatalErrorIn("Foam::featureEdgeMesh::edgedirection")
+        FatalErrorIn("Foam::extendedFeatureEdgeMesh::edgedirection")
             << "Requested ptI " << ptI << " is not a point on the requested "
             << "edgeI " << edgeI << ". edgeI start and end: "
             << e.start() << " " << e.end()
@@ -125,13 +126,14 @@ inline Foam::vector Foam::featureEdgeMesh::edgeDirection
 }
 
 
-inline const Foam::labelListList& Foam::featureEdgeMesh::edgeNormals() const
+inline const Foam::labelListList& Foam::extendedFeatureEdgeMesh::edgeNormals()
+const
 {
     return edgeNormals_;
 }
 
 
-inline Foam::vectorField Foam::featureEdgeMesh::edgeNormals
+inline Foam::vectorField Foam::extendedFeatureEdgeMesh::edgeNormals
 (
     const labelList& edgeNormIs
 ) const
@@ -147,27 +149,28 @@ inline Foam::vectorField Foam::featureEdgeMesh::edgeNormals
 }
 
 
-inline Foam::vectorField Foam::featureEdgeMesh::edgeNormals(label edgeI) const
+inline Foam::vectorField Foam::extendedFeatureEdgeMesh::edgeNormals(label edgeI)
+const
 {
     return edgeNormals(edgeNormals_[edgeI]);
 }
 
 
 inline const Foam::labelListList&
-Foam::featureEdgeMesh::featurePointNormals() const
+Foam::extendedFeatureEdgeMesh::featurePointNormals() const
 {
     return featurePointNormals_;
 }
 
 
-inline Foam::vectorField Foam::featureEdgeMesh::featurePointNormals
+inline Foam::vectorField Foam::extendedFeatureEdgeMesh::featurePointNormals
 (
     label ptI
 ) const
 {
     if (!featurePoint(ptI))
     {
-        WarningIn("vectorField featureEdgeMesh::featurePointNormals")
+        WarningIn("vectorField extendedFeatureEdgeMesh::featurePointNormals")
             << "Requesting the normals of a non-feature point. "
             << "Returned zero length vectorField."
             << endl;
@@ -188,13 +191,14 @@ inline Foam::vectorField Foam::featureEdgeMesh::featurePointNormals
 }
 
 
-inline const Foam::labelList& Foam::featureEdgeMesh::regionEdges() const
+inline const Foam::labelList& Foam::extendedFeatureEdgeMesh::regionEdges() const
 {
     return regionEdges_;
 }
 
 
-inline Foam::featureEdgeMesh::pointStatus Foam::featureEdgeMesh::getPointStatus
+inline Foam::extendedFeatureEdgeMesh::pointStatus
+Foam::extendedFeatureEdgeMesh::getPointStatus
 (
     label ptI
 ) const
@@ -218,7 +222,8 @@ inline Foam::featureEdgeMesh::pointStatus Foam::featureEdgeMesh::getPointStatus
 }
 
 
-inline Foam::featureEdgeMesh::edgeStatus Foam::featureEdgeMesh::getEdgeStatus
+inline Foam::extendedFeatureEdgeMesh::edgeStatus
+Foam::extendedFeatureEdgeMesh::getEdgeStatus
 (
     label edgeI
 ) const
diff --git a/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.C b/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.C
index 78ee016cfd..5b3cea13f4 100644
--- a/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.C
+++ b/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,34 +24,15 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "featureEdgeMesh.H"
-#include "triSurface.H"
-#include "Random.H"
-#include "Time.H"
-#include "meshTools.H"
-#include "linePointRef.H"
-#include "ListListOps.H"
-#include "OFstream.H"
-#include "IFstream.H"
-#include "unitConversion.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-defineTypeNameAndDebug(Foam::featureEdgeMesh, 0);
+namespace Foam
+{
 
-Foam::scalar Foam::featureEdgeMesh::cosNormalAngleTol_ =
-    Foam::cos(degToRad(0.1));
+defineTypeNameAndDebug(featureEdgeMesh, 0);
 
-
-Foam::label Foam::featureEdgeMesh::convexStart_ = 0;
-
-
-Foam::label Foam::featureEdgeMesh::externalStart_ = 0;
-
-
-Foam::label Foam::featureEdgeMesh::nPointTypes = 4;
-
-
-Foam::label Foam::featureEdgeMesh::nEdgeTypes = 5;
+}
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -59,21 +40,7 @@ Foam::label Foam::featureEdgeMesh::nEdgeTypes = 5;
 Foam::featureEdgeMesh::featureEdgeMesh(const IOobject& io)
 :
     regIOobject(io),
-    edgeMesh(pointField(0), edgeList(0)),
-    concaveStart_(0),
-    mixedStart_(0),
-    nonFeatureStart_(0),
-    internalStart_(0),
-    flatStart_(0),
-    openStart_(0),
-    multipleStart_(0),
-    normals_(0),
-    edgeDirections_(0),
-    edgeNormals_(0),
-    featurePointNormals_(0),
-    regionEdges_(0),
-    edgeTree_(),
-    edgeTreesByType_()
+    edgeMesh(pointField(0), edgeList(0))
 {
     if
     (
@@ -82,47 +49,8 @@ Foam::featureEdgeMesh::featureEdgeMesh(const IOobject& io)
      || (io.readOpt() == IOobject::READ_IF_PRESENT && headerOk())
     )
     {
-        if (readOpt() == IOobject::MUST_READ_IF_MODIFIED)
-        {
-            WarningIn("featureEdgeMesh::featureEdgeMesh(const IOobject&)")
-                << "Specified IOobject::MUST_READ_IF_MODIFIED but class"
-                << " does not support automatic rereading."
-                << endl;
-        }
-
-        Istream& is = readStream(typeName);
-
-        is  >> *this
-            >> concaveStart_
-            >> mixedStart_
-            >> nonFeatureStart_
-            >> internalStart_
-            >> flatStart_
-            >> openStart_
-            >> multipleStart_
-            >> normals_
-            >> edgeNormals_
-            >> featurePointNormals_
-            >> regionEdges_;
-
+        readStream(typeName) >> *this;
         close();
-
-        {
-            // Calculate edgeDirections
-
-            const edgeList& eds(edges());
-
-            const pointField& pts(points());
-
-            edgeDirections_.setSize(eds.size());
-
-            forAll(eds, eI)
-            {
-                edgeDirections_[eI] = eds[eI].vec(pts);
-            }
-
-            edgeDirections_ /= mag(edgeDirections_);
-        }
     }
 
     if (debug)
@@ -136,884 +64,43 @@ Foam::featureEdgeMesh::featureEdgeMesh(const IOobject& io)
 }
 
 
+//- Construct from components
 Foam::featureEdgeMesh::featureEdgeMesh
 (
     const IOobject& io,
-    const featureEdgeMesh& fem
+    const pointField& points,
+    const edgeList& edges
 )
 :
     regIOobject(io),
-    edgeMesh(fem),
-    concaveStart_(fem.concaveStart()),
-    mixedStart_(fem.mixedStart()),
-    nonFeatureStart_(fem.nonFeatureStart()),
-    internalStart_(fem.internalStart()),
-    flatStart_(fem.flatStart()),
-    openStart_(fem.openStart()),
-    multipleStart_(fem.multipleStart()),
-    normals_(fem.normals()),
-    edgeDirections_(fem.edgeDirections()),
-    edgeNormals_(fem.edgeNormals()),
-    featurePointNormals_(fem.featurePointNormals()),
-    regionEdges_(fem.regionEdges()),
-    edgeTree_(),
-    edgeTreesByType_()
+    edgeMesh(points, edges)
 {}
 
 
+// Construct as copy
 Foam::featureEdgeMesh::featureEdgeMesh
 (
     const IOobject& io,
-    const Xfer<pointField>& pointLst,
-    const Xfer<edgeList>& edgeLst
+    const featureEdgeMesh& em
 )
 :
     regIOobject(io),
-    edgeMesh(pointLst, edgeLst),
-    concaveStart_(0),
-    mixedStart_(0),
-    nonFeatureStart_(0),
-    internalStart_(0),
-    flatStart_(0),
-    openStart_(0),
-    multipleStart_(0),
-    normals_(0),
-    edgeDirections_(0),
-    edgeNormals_(0),
-    featurePointNormals_(0),
-    regionEdges_(0),
-    edgeTree_(),
-    edgeTreesByType_()
+    edgeMesh(em)
 {}
 
 
-Foam::featureEdgeMesh::featureEdgeMesh
-(
-    const surfaceFeatures& sFeat,
-    const objectRegistry& obr,
-    const fileName& sFeatFileName
-)
-:
-    regIOobject
-    (
-        IOobject
-        (
-            sFeatFileName,
-            obr.time().constant(),
-            "featureEdgeMesh",
-            obr,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        )
-    ),
-    edgeMesh(pointField(0), edgeList(0)),
-    concaveStart_(-1),
-    mixedStart_(-1),
-    nonFeatureStart_(-1),
-    internalStart_(-1),
-    flatStart_(-1),
-    openStart_(-1),
-    multipleStart_(-1),
-    normals_(0),
-    edgeDirections_(0),
-    edgeNormals_(0),
-    featurePointNormals_(0),
-    regionEdges_(0),
-    edgeTree_(),
-    edgeTreesByType_()
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::featureEdgeMesh::readData(Istream& is)
 {
-    // Extract and reorder the data from surfaceFeatures
-
-    // References to the surfaceFeatures data
-    const triSurface& surf(sFeat.surface());
-    const pointField& sFeatLocalPts(surf.localPoints());
-    const edgeList& sFeatEds(surf.edges());
-
-    // Filling the featureEdgeMesh with the raw geometrical data.
-
-    label nFeatEds = sFeat.featureEdges().size();
-
-    DynamicList<point> tmpPts;
-    edgeList eds(nFeatEds);
-    DynamicList<vector> norms;
-    vectorField edgeDirections(nFeatEds);
-    labelListList edgeNormals(nFeatEds);
-    DynamicList<label> regionEdges;
-
-    // Mapping between old and new indices, there is entry in the map for each
-    // of sFeatLocalPts, -1 means that this point hasn't been used (yet), >= 0
-    // corresponds to the index
-    labelList pointMap(sFeatLocalPts.size(), -1);
-
-    // Noting when the normal of a face has been used so not to duplicate
-    labelList faceMap(surf.size(), -1);
-
-    // Collecting the status of edge for subsequent sorting
-    List<edgeStatus> edStatus(nFeatEds, NONE);
-
-    forAll(sFeat.featurePoints(), i)
-    {
-        label sFPI = sFeat.featurePoints()[i];
-
-        tmpPts.append(sFeatLocalPts[sFPI]);
-
-        pointMap[sFPI] = tmpPts.size() - 1;
-    }
-
-    // All feature points have been added
-    nonFeatureStart_ = tmpPts.size();
-
-    forAll(sFeat.featureEdges(), i)
-    {
-        label sFEI = sFeat.featureEdges()[i];
-
-        const edge& fE(sFeatEds[sFEI]);
-
-        // Check to see if the points have been already used
-        if (pointMap[fE.start()] == -1)
-        {
-             tmpPts.append(sFeatLocalPts[fE.start()]);
-
-             pointMap[fE.start()] = tmpPts.size() - 1;
-        }
-
-        eds[i].start() = pointMap[fE.start()];
-
-        if (pointMap[fE.end()] == -1)
-        {
-            tmpPts.append(sFeatLocalPts[fE.end()]);
-
-            pointMap[fE.end()] = tmpPts.size() - 1;
-        }
-
-        eds[i].end() = pointMap[fE.end()];
-
-        // Pick up the faces adjacent to the feature edge
-        const labelList& eFaces = surf.edgeFaces()[sFEI];
-
-        edgeNormals[i].setSize(eFaces.size());
-
-        forAll(eFaces, j)
-        {
-            label eFI = eFaces[j];
-
-            // Check to see if the points have been already used
-            if (faceMap[eFI] == -1)
-            {
-                norms.append(surf.faceNormals()[eFI]);
-
-                faceMap[eFI] = norms.size() - 1;
-            }
-
-            edgeNormals[i][j] = faceMap[eFI];
-        }
-
-        vector fC0tofC1(vector::zero);
-
-        if (eFaces.size() == 2)
-        {
-            fC0tofC1 =
-                surf[eFaces[1]].centre(surf.points())
-              - surf[eFaces[0]].centre(surf.points());
-        }
-
-        edStatus[i] = classifyEdge(norms, edgeNormals[i], fC0tofC1);
-
-        edgeDirections[i] = fE.vec(sFeatLocalPts);
-
-        if (i < sFeat.nRegionEdges())
-        {
-            regionEdges.append(i);
-        }
-    }
-
-    // Reorder the edges by classification
-
-    List<DynamicList<label> > allEds(nEdgeTypes);
-
-    DynamicList<label>& externalEds(allEds[0]);
-    DynamicList<label>& internalEds(allEds[1]);
-    DynamicList<label>& flatEds(allEds[2]);
-    DynamicList<label>& openEds(allEds[3]);
-    DynamicList<label>& multipleEds(allEds[4]);
-
-    forAll(eds, i)
-    {
-        edgeStatus eStat = edStatus[i];
-
-        if (eStat == EXTERNAL)
-        {
-            externalEds.append(i);
-        }
-        else if (eStat == INTERNAL)
-        {
-            internalEds.append(i);
-        }
-        else if (eStat == FLAT)
-        {
-            flatEds.append(i);
-        }
-        else if (eStat == OPEN)
-        {
-            openEds.append(i);
-        }
-        else if (eStat == MULTIPLE)
-        {
-            multipleEds.append(i);
-        }
-        else if (eStat == NONE)
-        {
-            FatalErrorIn
-            (
-                "Foam::featureEdgeMesh::featureEdgeMesh"
-                "("
-                "    const surfaceFeatures& sFeat,"
-                "    const objectRegistry& obr,"
-                "    const fileName& sFeatFileName"
-                ")"
-            )
-                << nl << "classifyEdge returned NONE on edge "
-                << eds[i]
-                << ". There is a problem with definition of this edge."
-                << nl << abort(FatalError);
-        }
-    }
-
-    internalStart_ = externalEds.size();
-    flatStart_ = internalStart_ + internalEds.size();
-    openStart_ = flatStart_ + flatEds.size();
-    multipleStart_ = openStart_ + openEds.size();
-
-    labelList edMap
-    (
-        ListListOps::combine<labelList>
-        (
-            allEds,
-            accessOp<labelList>()
-        )
-    );
-
-    edMap = invert(edMap.size(), edMap);
-
-    inplaceReorder(edMap, eds);
-    inplaceReorder(edMap, edStatus);
-    inplaceReorder(edMap, edgeDirections);
-    inplaceReorder(edMap, edgeNormals);
-    inplaceRenumber(edMap, regionEdges);
-
-    pointField pts(tmpPts);
-
-    // Initialise the edgeMesh
-    edgeMesh::operator=(edgeMesh(pts, eds));
-
-    // Initialise sorted edge related data
-    edgeDirections_ = edgeDirections/mag(edgeDirections);
-    edgeNormals_ = edgeNormals;
-    regionEdges_ = regionEdges;
-
-    // Normals are all now found and indirectly addressed, can also be stored
-    normals_ = vectorField(norms);
-
-    // Reorder the feature points by classification
-
-    List<DynamicList<label> > allPts(3);
-
-    DynamicList<label>& convexPts(allPts[0]);
-    DynamicList<label>& concavePts(allPts[1]);
-    DynamicList<label>& mixedPts(allPts[2]);
-
-    for (label i = 0; i < nonFeatureStart_; i++)
-    {
-        pointStatus ptStatus = classifyFeaturePoint(i);
-
-        if (ptStatus == CONVEX)
-        {
-            convexPts.append(i);
-        }
-        else if (ptStatus == CONCAVE)
-        {
-            concavePts.append(i);
-        }
-        else if (ptStatus == MIXED)
-        {
-            mixedPts.append(i);
-        }
-        else if (ptStatus == NONFEATURE)
-        {
-            FatalErrorIn
-            (
-                "Foam::featureEdgeMesh::featureEdgeMesh"
-                "("
-                "    const surfaceFeatures& sFeat,"
-                "    const objectRegistry& obr,"
-                "    const fileName& sFeatFileName"
-                ")"
-            )
-                << nl << "classifyFeaturePoint returned NONFEATURE on point at "
-                << points()[i]
-                << ". There is a problem with definition of this feature point."
-                << nl << abort(FatalError);
-        }
-    }
-
-    concaveStart_ = convexPts.size();
-    mixedStart_ = concaveStart_ + concavePts.size();
-
-    labelList ftPtMap
-    (
-        ListListOps::combine<labelList>
-        (
-            allPts,
-            accessOp<labelList>()
-        )
-    );
-
-    ftPtMap = invert(ftPtMap.size(), ftPtMap);
-
-    // Creating the ptMap from the ftPtMap with identity values up to the size
-    // of pts to create an oldToNew map for inplaceReorder
-
-    labelList ptMap(identity(pts.size()));
-
-    forAll(ftPtMap, i)
-    {
-        ptMap[i] = ftPtMap[i];
-    }
-
-    inplaceReorder(ptMap, pts);
-
-    forAll(eds, i)
-    {
-        inplaceRenumber(ptMap, eds[i]);
-    }
-
-    // Reinitialise the edgeMesh with sorted feature points and
-    // renumbered edges
-    edgeMesh::operator=(edgeMesh(pts, eds));
-
-    // Generate the featurePointNormals
-
-    labelListList featurePointNormals(nonFeatureStart_);
-
-    for (label i = 0; i < nonFeatureStart_; i++)
-    {
-        DynamicList<label> tmpFtPtNorms;
-
-        const labelList& ptEds = pointEdges()[i];
-
-        forAll(ptEds, j)
-        {
-            const labelList& ptEdNorms(edgeNormals[ptEds[j]]);
-
-            forAll(ptEdNorms, k)
-            {
-                if (findIndex(tmpFtPtNorms, ptEdNorms[k]) == -1)
-                {
-                    bool addNormal = true;
-
-                    // Check that the normal direction is unique at this feature
-                    forAll(tmpFtPtNorms, q)
-                    {
-                        if
-                        (
-                            (normals_[ptEdNorms[k]] & normals_[tmpFtPtNorms[q]])
-                          > cosNormalAngleTol_
-                        )
-                        {
-                            // Parallel to an existing normal, do not add
-                            addNormal = false;
-
-                            break;
-                        }
-                    }
-
-                    if (addNormal)
-                    {
-                        tmpFtPtNorms.append(ptEdNorms[k]);
-                    }
-                }
-            }
-        }
-
-        featurePointNormals[i] = tmpFtPtNorms;
-    }
-
-    featurePointNormals_ = featurePointNormals;
-}
-
-
-Foam::featureEdgeMesh::featureEdgeMesh
-(
-    const IOobject& io,
-    const pointField& pts,
-    const edgeList& eds,
-    label concaveStart,
-    label mixedStart,
-    label nonFeatureStart,
-    label internalStart,
-    label flatStart,
-    label openStart,
-    label multipleStart,
-    const vectorField& normals,
-    const vectorField& edgeDirections,
-    const labelListList& edgeNormals,
-    const labelListList& featurePointNormals,
-    const labelList& regionEdges
-)
-:
-    regIOobject(io),
-    edgeMesh(pts, eds),
-    concaveStart_(concaveStart),
-    mixedStart_(mixedStart),
-    nonFeatureStart_(nonFeatureStart),
-    internalStart_(internalStart),
-    flatStart_(flatStart),
-    openStart_(openStart),
-    multipleStart_(multipleStart),
-    normals_(normals),
-    edgeDirections_(edgeDirections),
-    edgeNormals_(edgeNormals),
-    featurePointNormals_(featurePointNormals),
-    regionEdges_(regionEdges),
-    edgeTree_(),
-    edgeTreesByType_()
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::featureEdgeMesh::~featureEdgeMesh()
-{}
-
-
-// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
-
-Foam::featureEdgeMesh::pointStatus Foam::featureEdgeMesh::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::featureEdgeMesh::edgeStatus Foam::featureEdgeMesh::classifyEdge
-(
-    const List<vector>& norms,
-    const labelList& edNorms,
-    const vector& fC0tofC1
-) const
-{
-    label nEdNorms = edNorms.size();
-
-    if (nEdNorms == 1)
-    {
-        return OPEN;
-    }
-    else if (nEdNorms == 2)
-    {
-        const vector n0(norms[edNorms[0]]);
-        const vector n1(norms[edNorms[1]]);
-
-        if ((n0 & n1) > cosNormalAngleTol_)
-        {
-            return FLAT;
-        }
-        else if ((fC0tofC1 & n0) > 0.0)
-        {
-            return INTERNAL;
-        }
-        else
-        {
-            return EXTERNAL;
-        }
-    }
-    else if (nEdNorms > 2)
-    {
-        return MULTIPLE;
-    }
-    else
-    {
-        // There is a problem - the edge has no normals
-        return NONE;
-    }
-}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-void Foam::featureEdgeMesh::nearestFeatureEdge
-(
-    const point& sample,
-    scalar searchDistSqr,
-    pointIndexHit& info
-) const
-{
-    info = edgeTree().findNearest
-    (
-        sample,
-        searchDistSqr
-    );
-}
-
-
-void Foam::featureEdgeMesh::nearestFeatureEdge
-(
-    const pointField& samples,
-    const scalarField& searchDistSqr,
-    List<pointIndexHit>& info
-) const
-{
-    info.setSize(samples.size());
-
-    forAll(samples, i)
-    {
-        nearestFeatureEdge
-        (
-            samples[i],
-            searchDistSqr[i],
-            info[i]
-        );
-    }
-}
-
-
-void Foam::featureEdgeMesh::nearestFeatureEdgeByType
-(
-    const point& sample,
-    const scalarField& searchDistSqr,
-    List<pointIndexHit>& info
-) const
-{
-    const PtrList<indexedOctree<treeDataEdge> >& edgeTrees = edgeTreesByType();
-
-    info.setSize(edgeTrees.size());
-
-    labelList sliceStarts(edgeTrees.size());
-
-    sliceStarts[0] = externalStart_;
-    sliceStarts[1] = internalStart_;
-    sliceStarts[2] = flatStart_;
-    sliceStarts[3] = openStart_;
-    sliceStarts[4] = multipleStart_;
-
-    forAll(edgeTrees, i)
-    {
-        info[i] = edgeTrees[i].findNearest
-        (
-            sample,
-            searchDistSqr[i]
-        );
-
-        // The index returned by the indexedOctree is local to the slice of
-        // edges it was supplied with, return the index to the value in the
-        // complete edge list
-
-        info[i].setIndex(info[i].index() + sliceStarts[i]);
-    }
-}
-
-
-const Foam::indexedOctree<Foam::treeDataEdge>&
-Foam::featureEdgeMesh::edgeTree() const
-{
-    if (edgeTree_.empty())
-    {
-        Random rndGen(17301893);
-
-        // Slightly extended bb. Slightly off-centred just so on symmetric
-        // geometry there are less face/edge aligned items.
-        treeBoundBox bb
-        (
-            treeBoundBox(points()).extend(rndGen, 1E-4)
-        );
-
-        bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-        bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-
-        labelList allEdges(identity(edges().size()));
-
-        edgeTree_.reset
-        (
-            new indexedOctree<treeDataEdge>
-            (
-                treeDataEdge
-                (
-                    false,          // cachebb
-                    edges(),        // edges
-                    points(),       // points
-                    allEdges        // selected edges
-                ),
-                bb,     // bb
-                8,      // maxLevel
-                10,     // leafsize
-                3.0     // duplicity
-            )
-        );
-    }
-
-    return edgeTree_();
-}
-
-
-const Foam::PtrList<Foam::indexedOctree<Foam::treeDataEdge> >&
-Foam::featureEdgeMesh::edgeTreesByType() const
-{
-    if (edgeTreesByType_.size() == 0)
-    {
-        edgeTreesByType_.setSize(nEdgeTypes);
-
-        Random rndGen(872141);
-
-        // Slightly extended bb. Slightly off-centred just so on symmetric
-        // geometry there are less face/edge aligned items.
-        treeBoundBox bb
-        (
-            treeBoundBox(points()).extend(rndGen, 1E-4)
-        );
-
-        bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-        bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-
-        labelListList sliceEdges(nEdgeTypes);
-
-        // External edges
-        sliceEdges[0] =
-            identity(internalStart_ - externalStart_) + externalStart_;
-
-        // Internal edges
-        sliceEdges[1] = identity(flatStart_ - internalStart_) + internalStart_;
-
-        // Flat edges
-        sliceEdges[2] = identity(openStart_ - flatStart_) + flatStart_;
-
-        // Open edges
-        sliceEdges[3] = identity(multipleStart_ - openStart_) + openStart_;
-
-        // Multiple edges
-        sliceEdges[4] =
-            identity(edges().size() - multipleStart_) + multipleStart_;
-
-        forAll(edgeTreesByType_, i)
-        {
-            edgeTreesByType_.set
-            (
-                i,
-                new indexedOctree<treeDataEdge>
-                (
-                    treeDataEdge
-                    (
-                        false,          // cachebb
-                        edges(),        // edges
-                        points(),       // points
-                        sliceEdges[i]   // selected edges
-                    ),
-                    bb,     // bb
-                    8,      // maxLevel
-                    10,     // leafsize
-                    3.0     // duplicity
-                )
-            );
-        }
-    }
-
-    return edgeTreesByType_;
-}
-
-
-void Foam::featureEdgeMesh::writeObj
-(
-    const fileName& prefix
-) const
-{
-    Pout<< nl << "Writing featureEdgeMesh components to " << prefix << endl;
-
-    label verti = 0;
-
-    edgeMesh::write(prefix + "_edgeMesh.obj");
-
-    OFstream convexFtPtStr(prefix + "_convexFeaturePts.obj");
-    Pout<< "Writing convex feature points to " << convexFtPtStr.name() << endl;
-
-    for(label i = 0; i < concaveStart_; i++)
-    {
-        meshTools::writeOBJ(convexFtPtStr, points()[i]);
-    }
-
-    OFstream concaveFtPtStr(prefix + "_concaveFeaturePts.obj");
-    Pout<< "Writing concave feature points to "
-        << concaveFtPtStr.name() << endl;
-
-    for(label i = concaveStart_; i < mixedStart_; i++)
-    {
-        meshTools::writeOBJ(concaveFtPtStr, points()[i]);
-    }
-
-    OFstream mixedFtPtStr(prefix + "_mixedFeaturePts.obj");
-    Pout<< "Writing mixed feature points to " << mixedFtPtStr.name() << endl;
-
-    for(label i = mixedStart_; i < nonFeatureStart_; i++)
-    {
-        meshTools::writeOBJ(mixedFtPtStr, points()[i]);
-    }
-
-    OFstream mixedFtPtStructureStr(prefix + "_mixedFeaturePtsStructure.obj");
-    Pout<< "Writing mixed feature point structure to "
-        << mixedFtPtStructureStr.name() << endl;
-
-    verti = 0;
-    for(label i = mixedStart_; i < nonFeatureStart_; i++)
-    {
-        const labelList& ptEds = pointEdges()[i];
-
-        forAll(ptEds, j)
-        {
-            const edge& e = edges()[ptEds[j]];
-
-            meshTools::writeOBJ(mixedFtPtStructureStr, points()[e[0]]); verti++;
-            meshTools::writeOBJ(mixedFtPtStructureStr, points()[e[1]]); verti++;
-            mixedFtPtStructureStr << "l " << verti-1 << ' ' << verti << endl;
-        }
-    }
-
-    OFstream externalStr(prefix + "_externalEdges.obj");
-    Pout<< "Writing external edges to " << externalStr.name() << endl;
-
-    verti = 0;
-    for (label i = externalStart_; i < internalStart_; i++)
-    {
-        const edge& e = edges()[i];
-
-        meshTools::writeOBJ(externalStr, points()[e[0]]); verti++;
-        meshTools::writeOBJ(externalStr, points()[e[1]]); verti++;
-        externalStr << "l " << verti-1 << ' ' << verti << endl;
-    }
-
-    OFstream internalStr(prefix + "_internalEdges.obj");
-    Pout<< "Writing internal edges to " << internalStr.name() << endl;
-
-    verti = 0;
-    for (label i = internalStart_; i < flatStart_; i++)
-    {
-        const edge& e = edges()[i];
-
-        meshTools::writeOBJ(internalStr, points()[e[0]]); verti++;
-        meshTools::writeOBJ(internalStr, points()[e[1]]); verti++;
-        internalStr << "l " << verti-1 << ' ' << verti << endl;
-    }
-
-    OFstream flatStr(prefix + "_flatEdges.obj");
-    Pout<< "Writing flat edges to " << flatStr.name() << endl;
-
-    verti = 0;
-    for (label i = flatStart_; i < openStart_; i++)
-    {
-        const edge& e = edges()[i];
-
-        meshTools::writeOBJ(flatStr, points()[e[0]]); verti++;
-        meshTools::writeOBJ(flatStr, points()[e[1]]); verti++;
-        flatStr << "l " << verti-1 << ' ' << verti << endl;
-    }
-
-    OFstream openStr(prefix + "_openEdges.obj");
-    Pout<< "Writing open edges to " << openStr.name() << endl;
-
-    verti = 0;
-    for (label i = openStart_; i < multipleStart_; i++)
-    {
-        const edge& e = edges()[i];
-
-        meshTools::writeOBJ(openStr, points()[e[0]]); verti++;
-        meshTools::writeOBJ(openStr, points()[e[1]]); verti++;
-        openStr << "l " << verti-1 << ' ' << verti << endl;
-    }
-
-    OFstream multipleStr(prefix + "_multipleEdges.obj");
-    Pout<< "Writing multiple edges to " << multipleStr.name() << endl;
-
-    verti = 0;
-    for (label i = multipleStart_; i < edges().size(); i++)
-    {
-        const edge& e = edges()[i];
-
-        meshTools::writeOBJ(multipleStr, points()[e[0]]); verti++;
-        meshTools::writeOBJ(multipleStr, points()[e[1]]); verti++;
-        multipleStr << "l " << verti-1 << ' ' << verti << endl;
-    }
-
-    OFstream regionStr(prefix + "_regionEdges.obj");
-    Pout<< "Writing region edges to " << regionStr.name() << endl;
-
-    verti = 0;
-    forAll(regionEdges_, i)
-    {
-        const edge& e = edges()[regionEdges_[i]];
-
-        meshTools::writeOBJ(regionStr, points()[e[0]]); verti++;
-        meshTools::writeOBJ(regionStr, points()[e[1]]); verti++;
-        regionStr << "l " << verti-1 << ' ' << verti << endl;
-    }
+    is >> *this;
+    return !is.bad();
 }
 
 
 bool Foam::featureEdgeMesh::writeData(Ostream& os) const
 {
-    os  << "// points, edges, concaveStart, mixedStart, nonFeatureStart, " << nl
-        << "// internalStart, flatStart, openStart, multipleStart, " << nl
-        << "// normals, edgeNormals, featurePointNormals, regionEdges" << nl
-        << endl;
-
-    os  << points() << nl
-        << edges() << nl
-        << concaveStart_ << token::SPACE
-        << mixedStart_ << token::SPACE
-        << nonFeatureStart_ << token::SPACE
-        << internalStart_ << token::SPACE
-        << flatStart_ << token::SPACE
-        << openStart_ << token::SPACE
-        << multipleStart_ << nl
-        << normals_ << nl
-        << edgeNormals_ << nl
-        << featurePointNormals_ << nl
-        << regionEdges_
-        << endl;
+    os << *this;
 
     return os.good();
 }
diff --git a/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.H b/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.H
index 4e1e360fa5..9e7982a7c6 100644
--- a/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.H
+++ b/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,27 +25,12 @@ Class
     Foam::featureEdgeMesh
 
 Description
+    edgeMesh + IO.
 
-    Feature points are a sorted subset at the start of the overall points list:
-        0 .. concaveStart_-1                : convex points (w.r.t normals)
-        concaveStart_-1 .. mixedStart_-1    : concave points
-        mixedStart_ .. nonFeatureStart_     : mixed internal/external points
-        nonFeatureStart_ .. size-1          : non-feature points
-
-    Feature edges are the edgeList of the edgeMesh and are sorted:
-        0 .. internalStart_-1           : external edges (convex w.r.t normals)
-        internalStart_ .. flatStart_-1  : internal edges (concave)
-        flatStart_ .. openStart_-1      : flat edges (neither concave or convex)
-                                          can arise from region interfaces on
-                                          flat surfaces
-        openStart_ .. multipleStart_-1  : open edges (e.g. from baffle surfaces)
-        multipleStart_ .. size-1        : multiply connected edges
-
-    The edge direction and feature edge and feature point adjacent normals
-    are stored.
+    See also extendedFeatureEdgeMesh type which stores additional classification
+    of features.
 
 SourceFiles
-    featureEdgeMeshI.H
     featureEdgeMesh.C
 
 \*---------------------------------------------------------------------------*/
@@ -54,12 +39,7 @@ SourceFiles
 #define featureEdgeMesh_H
 
 #include "edgeMesh.H"
-#include "surfaceFeatures.H"
-#include "objectRegistry.H"
-#include "IOdictionary.H"
-#include "indexedOctree.H"
-#include "treeDataEdge.H"
-#include "pointIndexHit.H"
+#include "regIOobject.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -67,7 +47,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                       Class featureEdgeMesh Declaration
+                           Class featureEdgeMesh Declaration
 \*---------------------------------------------------------------------------*/
 
 class featureEdgeMesh
@@ -78,281 +58,31 @@ class featureEdgeMesh
 
 public:
 
-    //- Runtime type information
     TypeName("featureEdgeMesh");
 
-    enum pointStatus
-    {
-        CONVEX,     // Fully convex point (w.r.t normals)
-        CONCAVE,    // Fully concave point
-        MIXED,      // A point surrounded by both convex and concave edges
-        NONFEATURE  // Not a feature point
-    };
-
-    enum edgeStatus
-    {
-        EXTERNAL,   // "Convex" edge
-        INTERNAL,   // "Concave" edge
-        FLAT,       // Neither concave or convex, on a flat surface
-        OPEN,       // i.e. only connected to one face
-        MULTIPLE,   // Multiply connected (connected to more than two faces)
-        NONE        // Not a classified feature edge (consistency with
-                    // surfaceFeatures)
-    };
-
-private:
-
-    // Static data
-
-        //- Angular closeness tolerance for treating normals as the same
-        static scalar cosNormalAngleTol_;
-
-        //- Index of the start of the convex feature points - static as 0
-        static label convexStart_;
-
-        //- Index of the start of the external feature edges - static as 0
-        static label externalStart_;
-
-
-    // Private data
-
-        //- Index of the start of the concave feature points
-        label concaveStart_;
-
-        //- Index of the start of the mixed type feature points
-        label mixedStart_;
-
-        //- Index of the start of the non-feature points
-        label nonFeatureStart_;
-
-        //- Index of the start of the internal feature edges
-        label internalStart_;
-
-        //- Index of the start of the flat feature edges
-        label flatStart_;
-
-        //- Index of the start of the open feature edges
-        label openStart_;
-
-        //- Index of the start of the multiply-connected feature edges
-        label multipleStart_;
-
-        //- Normals of the features, to be referred to by index by both feature
-        //  points and edges, unsorted
-        vectorField normals_;
-
-        //- Flat and open edges require the direction of the edge
-        vectorField edgeDirections_;
-
-        //- Indices of the normals that are adjacent to the feature edges
-        labelListList edgeNormals_;
-
-        //- Indices of the normals that are adjacent to the feature points
-        labelListList featurePointNormals_;
-
-        //- Feature edges which are on the boundary between regions
-        labelList regionEdges_;
-
-        //- Search tree for all edges
-        mutable autoPtr<indexedOctree<treeDataEdge> > edgeTree_;
-
-        //- Individual search trees for each type of edge
-        mutable PtrList<indexedOctree<treeDataEdge> > edgeTreesByType_;
-
-
-    // Private Member Functions
-
-        //- Classify the type of feature point.  Requires valid stored member
-        //  data for edges and normals.
-        pointStatus classifyFeaturePoint(label ptI) const;
-
-        //- Classify the type of feature edge.  Requires face centre 0 to face
-        //  centre 1 vector to distinguish internal from external
-        edgeStatus classifyEdge
-        (
-            const List<vector>& norms,
-            const labelList& edNorms,
-            const vector& fC0tofC1
-        ) const;
-
-
-public:
-
-    // Static data
-
-        //- Number of possible point types (i.e. number of slices)
-        static label nPointTypes;
-
-        //- Number of possible feature edge types (i.e. number of slices)
-        static label nEdgeTypes;
 
     // Constructors
 
         //- Construct (read) given an IOobject
         featureEdgeMesh(const IOobject&);
 
-        //- Construct as copy
-        explicit featureEdgeMesh(const IOobject&, const featureEdgeMesh&);
-
-        //- Construct by transferring components (points, edges)
+        //- Construct from featureEdgeMesh data
         featureEdgeMesh
         (
             const IOobject&,
-            const Xfer<pointField>&,
-            const Xfer<edgeList>&
+            const pointField&,
+            const edgeList&
         );
 
-        //- Construct (read) given surfaceFeatures, an objectRegistry and a
-        //  fileName to write to.  Extracts, classifies and reorders the data
-        //  from surfaceFeatures.
-        featureEdgeMesh
-        (
-            const surfaceFeatures& sFeat,
-            const objectRegistry& obr,
-            const fileName& sFeatFileName
-        );
-
-        //- Construct from all components
-        featureEdgeMesh
-        (
-            const IOobject& io,
-            const pointField& pts,
-            const edgeList& eds,
-            label concaveStart,
-            label mixedStart,
-            label nonFeatureStart,
-            label internalStart,
-            label flatStart,
-            label openStart,
-            label multipleStart,
-            const vectorField& normals,
-            const vectorField& edgeDirections,
-            const labelListList& edgeNormals,
-            const labelListList& featurePointNormals,
-            const labelList& regionEdges
-        );
+        //- Construct as copy
+        featureEdgeMesh(const IOobject&, const featureEdgeMesh&);
 
 
-    //- Destructor
-    ~featureEdgeMesh();
+        //- ReadData function required for regIOobject read operation
+        virtual bool readData(Istream&);
 
-
-    // Member Functions
-
-        // Find
-
-            //- Find nearest surface edge for the sample point.
-            void nearestFeatureEdge
-            (
-                const point& sample,
-                scalar searchDistSqr,
-                pointIndexHit& info
-            ) const;
-
-            //- Find nearest surface edge for each sample point.
-            void nearestFeatureEdge
-            (
-                const pointField& samples,
-                const scalarField& searchDistSqr,
-                List<pointIndexHit>& info
-            ) const;
-
-            //- Find the nearest point on each type of feature edge
-            void nearestFeatureEdgeByType
-            (
-                const point& sample,
-                const scalarField& searchDistSqr,
-                List<pointIndexHit>& info
-            ) const;
-
-        // Access
-
-            //- Return the index of the start of the convex feature points
-            inline label convexStart() const;
-
-            //- Return the index of the start of the concave feature points
-            inline label concaveStart() const;
-
-            //- Return the index of the start of the mixed type feature points
-            inline label mixedStart() const;
-
-            //- Return the index of the start of the non-feature points
-            inline label nonFeatureStart() const;
-
-            //- Return the index of the start of the external feature edges
-            inline label externalStart() const;
-
-            //- Return the index of the start of the internal feature edges
-            inline label internalStart() const;
-
-            //- Return the index of the start of the flat feature edges
-            inline label flatStart() const;
-
-            //- Return the index of the start of the open feature edges
-            inline label openStart() const;
-
-            //- Return the index of the start of the multiply-connected feature
-            //  edges
-            inline label multipleStart() const;
-
-            //- Return whether or not the point index is a feature point
-            inline bool featurePoint(label ptI) const;
-
-            //- Return the normals of the surfaces adjacent to the feature edges
-            //  and points
-            inline const vectorField& normals() const;
-
-            //- Return the edgeDirection vectors
-            inline const vectorField& edgeDirections() const;
-
-            //- Return the direction of edgeI, pointing away from ptI
-            inline vector edgeDirection(label edgeI, label ptI) const;
-
-            //- Return the indices of the normals that are adjacent to the
-            //  feature edges
-            inline const labelListList& edgeNormals() const;
-
-            //- Return the normal vectors for a given set of normal indices
-            inline vectorField edgeNormals(const labelList& edgeNormIs) const;
-
-            //- Return the normal vectors for a given edge
-            inline vectorField edgeNormals(label edgeI) const;
-
-            //- Return the indices of the normals that are adjacent to the
-            //  feature points
-            inline const labelListList& featurePointNormals() const;
-
-            //- Return the normal vectors for a given feature point
-            inline vectorField featurePointNormals(label ptI) const;
-
-            //- Return the feature edges which are on the boundary between
-            //  regions
-            inline const labelList& regionEdges() const;
-
-            //- Return the pointStatus of a specified point
-            inline pointStatus getPointStatus(label ptI) const;
-
-            //- Return the edgeStatus of a specified edge
-            inline edgeStatus getEdgeStatus(label edgeI) const;
-
-            //- Demand driven construction of octree for boundary edges
-            const indexedOctree<treeDataEdge>& edgeTree() const;
-
-            //- Demand driven construction of octree for boundary edges by type
-            const PtrList<indexedOctree<treeDataEdge> >&
-            edgeTreesByType() const;
-
-
-        // Write
-
-            //- Write all components of the featureEdgeMesh 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;
+        //- WriteData function required for regIOobject write operation
+        virtual bool writeData(Ostream&) const;
 };
 
 
@@ -362,10 +92,6 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#include "featureEdgeMeshI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #endif
 
 // ************************************************************************* //