diff --git a/applications/utilities/surface/surfaceFeatureExtract/Make/options b/applications/utilities/surface/surfaceFeatureExtract/Make/options
index 02e49d32c6..f87dafaa5d 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/Make/options
+++ b/applications/utilities/surface/surfaceFeatureExtract/Make/options
@@ -1,8 +1,9 @@
 EXE_INC = \
-    -I$(LIB_SRC)/cfdTools/general/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/edgeMesh/lnInclude \
     -I$(LIB_SRC)/triSurface/lnInclude
 
 EXE_LIBS = \
     -lmeshTools \
+    -ledgeMesh \
     -ltriSurface
diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
index 83d84b92c8..7ba05f9d0e 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
@@ -29,10 +29,13 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
+
 #include "triangle.H"
 #include "triSurface.H"
 #include "argList.H"
+#include "Time.H"
 #include "surfaceFeatures.H"
+#include "featureEdgeMesh.H"
 #include "treeBoundBox.H"
 #include "meshTools.H"
 #include "OFstream.H"
@@ -135,7 +138,8 @@ int main(int argc, char *argv[])
         "remove edges within specified bounding box"
     );
 
-    argList args(argc, argv);
+#   include "setRootCase.H"
+#   include "createTime.H"
 
     Info<< "Feature line extraction is only valid on closed manifold surfaces."
         << endl;
@@ -282,7 +286,6 @@ int main(int argc, char *argv[])
     Info<< endl << "Writing edge objs." << endl;
     newSet.writeObj("final");
 
-
     Info<< nl
         << "Final feature set:" << nl
         << "    feature points : " << newSet.featurePoints().size() << nl
@@ -293,6 +296,22 @@ int main(int argc, char *argv[])
         << "        internal edges : " << newSet.nInternalEdges() << nl
         << endl;
 
+    // Extracting and writing a featureEdgeMesh
+
+    Pout<< nl << "Writing featureEdgeMesh to constant/featureEdgeMesh."
+        << endl;
+
+    featureEdgeMesh feMesh
+    (
+        newSet,
+        runTime,
+        surfFileName.lessExt().name() + ".featureEdgeMesh"
+    );
+
+    feMesh.writeObj(surfFileName.lessExt().name());
+
+    feMesh.write();
+
     Info<< "End\n" << endl;
 
     return 0;
diff --git a/src/edgeMesh/Make/options b/src/edgeMesh/Make/options
index 7e207d0dbe..4796c18e22 100644
--- a/src/edgeMesh/Make/options
+++ b/src/edgeMesh/Make/options
@@ -1,5 +1,9 @@
 EXE_INC = \
-    -I$(LIB_SRC)/fileFormats/lnInclude
+    -I$(LIB_SRC)/fileFormats/lnInclude \
+    -I$(LIB_SRC)/triSurface/lnInclude \ 
+    -I$(LIB_SRC)/meshTools/lnInclude    
 
 LIB_LIBS = \
+    -ltriSurface \
+    -lmeshTools \ 
     -lfileFormats
diff --git a/src/edgeMesh/edgeMesh.C b/src/edgeMesh/edgeMesh.C
index d5c8f12ef8..7dd7015144 100644
--- a/src/edgeMesh/edgeMesh.C
+++ b/src/edgeMesh/edgeMesh.C
@@ -387,14 +387,4 @@ void Foam::edgeMesh::mergePoints(const scalar mergeDist)
 }
 
 
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
-void Foam::edgeMesh::operator=(const edgeMesh& rhs)
-{
-    points_ = rhs.points_;
-    edges_ = rhs.edges_;
-    pointEdgesPtr_.clear();
-}
-
-
 // ************************************************************************* //
diff --git a/src/edgeMesh/edgeMeshI.H b/src/edgeMesh/edgeMeshI.H
index d5e38343ec..c15fbd14ce 100644
--- a/src/edgeMesh/edgeMeshI.H
+++ b/src/edgeMesh/edgeMeshI.H
@@ -69,4 +69,14 @@ inline Foam::edgeList& Foam::edgeMesh::storedEdges()
 }
 
 
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+void Foam::edgeMesh::operator=(const edgeMesh& rhs)
+{
+    points_ = rhs.points_;
+    edges_ = rhs.edges_;
+    pointEdgesPtr_.clear();
+}
+
+
 // ************************************************************************* //
diff --git a/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.C b/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.C
index 34bcd4388e..da89dcf2d5 100644
--- a/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.C
+++ b/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.C
@@ -24,18 +24,56 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "featureEdgeMesh.H"
+#include "triSurface.H"
+#include "Random.H"
+#include "Time.H"
+#include "meshTools.H"
+#include "linePointRef.H"
+#include "ListListOps.H"
+#include "OFstream.H"
+#include "IFstream.H"
+#include "unitConversion.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 defineTypeNameAndDebug(Foam::featureEdgeMesh, 0);
 
+Foam::scalar Foam::featureEdgeMesh::cosNormalAngleTol_ =
+    Foam::cos(degToRad(0.1));
+
+
+Foam::label Foam::featureEdgeMesh::convexStart_ = 0;
+
+
+Foam::label Foam::featureEdgeMesh::externalStart_ = 0;
+
+
+Foam::label Foam::featureEdgeMesh::nPointTypes = 4;
+
+
+Foam::label Foam::featureEdgeMesh::nEdgeTypes = 5;
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::featureEdgeMesh::featureEdgeMesh(const IOobject& io)
 :
     regIOobject(io),
-    edgeMesh(pointField(0), edgeList(0))
+    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
     (
@@ -52,9 +90,39 @@ Foam::featureEdgeMesh::featureEdgeMesh(const IOobject& io)
                 << 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)
@@ -71,11 +139,25 @@ Foam::featureEdgeMesh::featureEdgeMesh(const IOobject& io)
 Foam::featureEdgeMesh::featureEdgeMesh
 (
     const IOobject& io,
-    const featureEdgeMesh& em
+    const featureEdgeMesh& fem
 )
 :
     regIOobject(io),
-    edgeMesh(em)
+    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_()
 {}
 
 
@@ -87,42 +169,851 @@ Foam::featureEdgeMesh::featureEdgeMesh
 )
 :
     regIOobject(io),
-    edgeMesh(pointLst, edgeLst)
-{
-    if
+    edgeMesh(pointLst, edgeLst),
+    concaveStart_(0),
+    mixedStart_(0),
+    nonFeatureStart_(0),
+    internalStart_(0),
+    flatStart_(0),
+    openStart_(0),
+    multipleStart_(0),
+    normals_(0),
+    edgeDirections_(0),
+    edgeNormals_(0),
+    featurePointNormals_(0),
+    regionEdges_(0),
+    edgeTree_(),
+    edgeTreesByType_()
+{}
+
+
+Foam::featureEdgeMesh::featureEdgeMesh
+(
+    const surfaceFeatures& sFeat,
+    const objectRegistry& obr,
+    const fileName& sFeatFileName
+)
+:
+    regIOobject
     (
-        io.readOpt() == IOobject::MUST_READ
-     || io.readOpt() == IOobject::MUST_READ_IF_MODIFIED
-     || (io.readOpt() == IOobject::READ_IF_PRESENT && headerOk())
-    )
+        IOobject
+        (
+            sFeatFileName,
+            obr.time().constant(),
+            "featureEdgeMesh",
+            obr,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        )
+    ),
+    edgeMesh(pointField(0), edgeList(0)),
+    concaveStart_(-1),
+    mixedStart_(-1),
+    nonFeatureStart_(-1),
+    internalStart_(-1),
+    flatStart_(-1),
+    openStart_(-1),
+    multipleStart_(-1),
+    normals_(0),
+    edgeDirections_(0),
+    edgeNormals_(0),
+    featurePointNormals_(0),
+    regionEdges_(0),
+    edgeTree_(),
+    edgeTreesByType_()
+{
+    // Extract and reorder the data from surfaceFeatures
+
+    // References to the surfaceFeatures data
+    const triSurface& surf(sFeat.surface());
+    const pointField& sFeatLocalPts(surf.localPoints());
+    const edgeList& sFeatEds(surf.edges());
+
+    // Filling the featureEdgeMesh with the raw geometrical data.
+
+    label nFeatEds = sFeat.featureEdges().size();
+
+    DynamicList<point> tmpPts;
+    edgeList eds(nFeatEds);
+    DynamicList<vector> norms;
+    vectorField edgeDirections(nFeatEds);
+    labelListList edgeNormals(nFeatEds);
+    DynamicList<label> regionEdges;
+
+    // Mapping between old and new indices, there is entry in the map for each
+    // of sFeatLocalPts, -1 means that this point hasn't been used (yet), >= 0
+    // corresponds to the index
+    labelList pointMap(sFeatLocalPts.size(), -1);
+
+    // Noting when the normal of a face has been used so not to duplicate
+    labelList faceMap(surf.size(), -1);
+
+    // Collecting the status of edge for subsequent sorting
+    List<edgeStatus> edStatus(nFeatEds, NONE);
+
+    forAll(sFeat.featurePoints(), i)
     {
-        if (readOpt() == IOobject::MUST_READ_IF_MODIFIED)
+        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)
         {
-            WarningIn("featureEdgeMesh::featureEdgeMesh(const IOobject&)")
-                << "Specified IOobject::MUST_READ_IF_MODIFIED but class"
-                << " does not support automatic rereading."
-                << endl;
+             tmpPts.append(sFeatLocalPts[fE.start()]);
+
+             pointMap[fE.start()] = tmpPts.size() - 1;
         }
 
+        eds[i].start() = pointMap[fE.start()];
 
-        readStream(typeName) >> *this;
-        close();
+        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;
     }
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-bool Foam::featureEdgeMesh::readData(Istream& is)
+Foam::featureEdgeMesh::edgeStatus Foam::featureEdgeMesh::classifyEdge
+(
+    const List<vector>& norms,
+    const labelList& edNorms,
+    const vector& fC0tofC1
+) const
 {
-    is  >> *this;
-    return !is.bad();
+    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;
+    }
 }
 
 
 bool Foam::featureEdgeMesh::writeData(Ostream& os) const
 {
-    os  << *this;
+    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/featureEdgeMesh/featureEdgeMesh.H b/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.H
index 613c9bdd94..93f034d9aa 100644
--- a/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.H
+++ b/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.H
@@ -25,9 +25,27 @@ Class
     Foam::featureEdgeMesh
 
 Description
-    features (lines), readable from file
+
+    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
+    featureEdgeMeshI.H
     featureEdgeMesh.C
 
 \*---------------------------------------------------------------------------*/
@@ -36,7 +54,12 @@ SourceFiles
 #define featureEdgeMesh_H
 
 #include "edgeMesh.H"
-#include "regIOobject.H"
+#include "surfaceFeatures.H"
+#include "objectRegistry.H"
+#include "IOdictionary.H"
+#include "indexedOctree.H"
+#include "treeDataEdge.H"
+#include "pointIndexHit.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -58,11 +81,115 @@ 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
-        explicit featureEdgeMesh(const IOobject&);
+        featureEdgeMesh(const IOobject&);
 
         //- Construct as copy
         explicit featureEdgeMesh(const IOobject&, const featureEdgeMesh&);
@@ -75,14 +202,157 @@ public:
             const Xfer<edgeList>&
         );
 
-        //- Give precedence to the regIOobject write
-        using regIOobject::write;
+        //- 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
+        );
 
-        //- ReadData function required for regIOobject read operation
-        virtual bool readData(Istream&);
+        //- 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
+        );
 
-        //- WriteData function required for regIOobject write operation
-        virtual bool writeData(Ostream&) const;
+
+    //- Destructor
+    ~featureEdgeMesh();
+
+
+    // Member Functions
+
+        // Find
+
+            //- Find nearest surface edge for the sample point.
+            void nearestFeatureEdge
+            (
+                const point& sample,
+                scalar searchDistSqr,
+                pointIndexHit& info
+            ) const;
+
+            //- Find nearest surface edge for each sample point.
+            void nearestFeatureEdge
+            (
+                const pointField& samples,
+                const scalarField& searchDistSqr,
+                List<pointIndexHit>& info
+            ) const;
+
+            //- Find the nearest point on each type of feature edge
+            void nearestFeatureEdgeByType
+            (
+                const point& sample,
+                const scalarField& searchDistSqr,
+                List<pointIndexHit>& info
+            ) const;
+
+        // Access
+
+            //- Return the index of the start of the convex feature points
+            inline label convexStart() const;
+
+            //- Return the index of the start of the concave feature points
+            inline label concaveStart() const;
+
+            //- Return the index of the start of the mixed type feature points
+            inline label mixedStart() const;
+
+            //- Return the index of the start of the non-feature points
+            inline label nonFeatureStart() const;
+
+            //- Return the index of the start of the external feature edges
+            inline label externalStart() const;
+
+            //- Return the index of the start of the internal feature edges
+            inline label internalStart() const;
+
+            //- Return the index of the start of the flat feature edges
+            inline label flatStart() const;
+
+            //- Return the index of the start of the open feature edges
+            inline label openStart() const;
+
+            //- Return the index of the start of the multiply-connected feature
+            //  edges
+            inline label multipleStart() const;
+
+            //- Return whether or not the point index is a feature point
+            inline bool featurePoint(label ptI) const;
+
+            //- Return the normals of the surfaces adjacent to the feature edges
+            //  and points
+            inline const vectorField& normals() const;
+
+            //- Return the edgeDirection vectors
+            inline const vectorField& edgeDirections() const;
+
+            //- Return the direction of edgeI, pointing away from ptI
+            inline vector edgeDirection(label edgeI, label ptI) const;
+
+            //- Return the indices of the normals that are adjacent to the
+            //  feature edges
+            inline const labelListList& edgeNormals() const;
+
+            //- Return the normal vectors for a given set of normal indices
+            inline vectorField edgeNormals(const labelList& edgeNormIs) const;
+
+            //- Return the normal vectors for a given edge
+            inline vectorField edgeNormals(label edgeI) const;
+
+            //- Return the indices of the normals that are adjacent to the
+            //  feature points
+            inline const labelListList& featurePointNormals() const;
+
+            //- Return the normal vectors for a given feature point
+            inline vectorField featurePointNormals(label ptI) const;
+
+            //- Return the feature edges which are on the boundary between
+            //  regions
+            inline const labelList& regionEdges() const;
+
+            //- Return the pointStatus of a specified point
+            inline pointStatus getPointStatus(label ptI) const;
+
+            //- Return the edgeStatus of a specified edge
+            inline edgeStatus getEdgeStatus(label edgeI) const;
+
+            //- Demand driven construction of octree for boundary edges
+            const indexedOctree<treeDataEdge>& edgeTree() const;
+
+            //- Demand driven construction of octree for boundary edges by type
+            const PtrList<indexedOctree<treeDataEdge> >&
+            edgeTreesByType() const;
+
+
+        // Write
+
+            //- 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;
 };
 
 
@@ -92,6 +362,10 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#include "featureEdgeMeshI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //