diff --git a/applications/utilities/surface/surfaceFeatureExtract/Make/options b/applications/utilities/surface/surfaceFeatureExtract/Make/options
index 02e49d32c6..9c37e9e1be 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/Make/options
+++ b/applications/utilities/surface/surfaceFeatureExtract/Make/options
@@ -1,8 +1,10 @@
 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 8907ac46db..0522550729 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
@@ -30,10 +30,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"
@@ -108,9 +111,11 @@ int main(int argc, char *argv[])
     argList::validOptions.insert("minElem", "number of edges in feature");
     argList::validOptions.insert("subsetBox", "((x0 y0 z0)(x1 y1 z1))");
     argList::validOptions.insert("deleteBox", "((x0 y0 z0)(x1 y1 z1))");
-    argList args(argc, argv);
 
-    Pout<< "Feature line extraction is only valid on closed manifold surfaces." 
+#   include "setRootCase.H"
+#   include "createTime.H"
+
+    Pout<< "Feature line extraction is only valid on closed manifold surfaces."
         << endl;
 
 
@@ -159,7 +164,7 @@ int main(int argc, char *argv[])
 
         set = surfaceFeatures(surf, includedAngle);
 
-        Pout<< endl << "Writing initial features" << endl;    
+        Pout<< endl << "Writing initial features" << endl;
         set.write("initial.fSet");
         set.writeObj("initial");
     }
@@ -195,7 +200,7 @@ int main(int argc, char *argv[])
         minLen = readScalar(IStringStream(args.options()["minLen"])());
         Pout<< "Removing features of length < " << minLen << endl;
     }
-    
+
     label minElem = 0;
     if (args.options().found("minElem"))
     {
@@ -207,7 +212,7 @@ int main(int argc, char *argv[])
     if (minLen > 0 || minLen > 0)
     {
         set.trimFeatures(minLen, minElem);
-        Pout<< endl << "Removed small features" << endl;    
+        Pout<< endl << "Removed small features" << endl;
     }
 
 
@@ -254,11 +259,10 @@ int main(int argc, char *argv[])
 
     Pout<< endl << "Writing trimmed features to " << outFileName << endl;
     newSet.write(outFileName);
-    
+
     Pout<< endl << "Writing edge objs." << endl;
     newSet.writeObj("final");
 
-
     Pout<< nl
         << "Final feature set:" << nl
         << "    feature points : " << newSet.featurePoints().size() << nl
@@ -269,7 +273,22 @@ int main(int argc, char *argv[])
         << "        internal edges : " << newSet.nInternalEdges() << nl
         << endl;
 
-    Pout<< "End\n" << endl;
+    // Extracting and writing a featureEdgeMesh
+
+    Pout<< nl << "Writing featureEdgeMesh to constant/featureEdgeMesh."
+        << endl;
+
+    featureEdgeMesh feMesh
+    (
+        newSet,
+        runTime,
+        surfFileName.lessExt().name() + ".feMesh",
+        true
+    );
+
+    feMesh.write();
+
+    Pout<< nl << "End\n" << endl;
 
     return 0;
 }
diff --git a/src/edgeMesh/Make/files b/src/edgeMesh/Make/files
index 7ec66485e7..9dcdbca4ef 100644
--- a/src/edgeMesh/Make/files
+++ b/src/edgeMesh/Make/files
@@ -1,8 +1,6 @@
 primitiveEdgeMesh/primitiveEdgeMesh.C
 primitiveEdgeMesh/primitiveEdgeMeshIO.C
 
-edgeMesh/edgeMesh.C
-
 featureEdgeMesh/featureEdgeMesh.C
 
 LIB = $(FOAM_LIBBIN)/libedgeMesh
diff --git a/src/edgeMesh/Make/options b/src/edgeMesh/Make/options
index 280bd02ef2..d64c26ecad 100644
--- a/src/edgeMesh/Make/options
+++ b/src/edgeMesh/Make/options
@@ -1,7 +1,9 @@
 EXE_INC = \
+    -I$(LIB_SRC)/triSurface/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
 
 EXE_LIBS = \
+    -ltriSurface \
     -lmeshTools
 
 
diff --git a/src/edgeMesh/edgeMesh/edgeMesh.C b/src/edgeMesh/edgeMesh/edgeMesh.C
deleted file mode 100644
index f22b3163e5..0000000000
--- a/src/edgeMesh/edgeMesh/edgeMesh.C
+++ /dev/null
@@ -1,109 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "edgeMesh.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-defineTypeNameAndDebug(edgeMesh, 0);
-
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::edgeMesh::edgeMesh(const IOobject& io)
-:
-    regIOobject(io),
-    primitiveEdgeMesh(pointField(0), edgeList(0))
-{
-    if
-    (
-        io.readOpt() == IOobject::MUST_READ
-     || (io.readOpt() == IOobject::READ_IF_PRESENT && headerOk())
-    )
-    {
-        readStream(typeName) >> *this;
-        close();
-    }
-
-    if (debug)
-    {
-        Pout<< "edgeMesh::edgeMesh :"
-            << " constructed from IOobject :"
-            << " points:" << points().size()
-            << " edges:" << edges().size()
-            << endl;
-    }
-}
-
-
-//- Construct from components
-Foam::edgeMesh::edgeMesh
-(
-    const IOobject& io,
-    const pointField& points,
-    const edgeList& edges
-)
-:
-    regIOobject(io),
-    primitiveEdgeMesh(points, edges)
-{}
-
-
-// Construct as copy
-Foam::edgeMesh::edgeMesh
-(
-    const IOobject& io,
-    const edgeMesh& em
-)
-:
-    regIOobject(io),
-    primitiveEdgeMesh(em)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-bool Foam::edgeMesh::readData(Istream& is)
-{
-    is >> *this;
-    return !is.bad();
-}
-
-
-bool Foam::edgeMesh::writeData(Ostream& os) const
-{
-    os << *this;
-
-    return os.good();
-}
-
-
-// ************************************************************************* //
diff --git a/src/edgeMesh/edgeMesh/edgeMesh.H b/src/edgeMesh/edgeMesh/edgeMesh.H
deleted file mode 100644
index ec3d7b7932..0000000000
--- a/src/edgeMesh/edgeMesh/edgeMesh.H
+++ /dev/null
@@ -1,95 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Class
-    Foam::edgeMesh
-
-Description
-    edges (lines), readable from file
-
-SourceFiles
-    edgeMesh.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef edgeMesh_H
-#define edgeMesh_H
-
-#include "primitiveEdgeMesh.H"
-#include "regIOobject.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                            Class edgeMesh Declaration
-\*---------------------------------------------------------------------------*/
-
-class edgeMesh
-:
-    public regIOobject,
-    public primitiveEdgeMesh
-{
-
-public:
-
-    TypeName("edgeMesh");
-
-
-    // Constructors
-
-        //- Construct (read) given an IOobject
-        edgeMesh(const IOobject&);
-
-        //- Construct from edgeMesh data
-        edgeMesh
-        (
-            const IOobject&,
-            const pointField&,
-            const edgeList&
-        );
-
-        //- Construct as copy
-        edgeMesh(const IOobject&, const edgeMesh&);
-
-
-        //- ReadData function required for regIOobject read operation
-        virtual bool readData(Istream&);
-
-        //- WriteData function required for regIOobject write operation
-        virtual bool writeData(Ostream&) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.C b/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.C
index 9045ae5cd9..eb893e0b04 100644
--- a/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.C
+++ b/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.C
@@ -25,22 +25,458 @@ 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"
 
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(featureEdgeMesh, 0);
+
+}
+
+
+Foam::scalar Foam::featureEdgeMesh::cosNormalAngleTol_ = Foam::cos
+(
+    0.1*mathematicalConstant::pi/180.0
+);
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::featureEdgeMesh::featureEdgeMesh(const IOobject& io)
 :
-    edgeMesh(io),
-    normals_(),
-    edgeNormals_(edges().size()),
-    allEdges_(identity(edges().size()))
+    regIOobject(io),
+    primitiveEdgeMesh(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)
+{
+    if
+    (
+        io.readOpt() == IOobject::MUST_READ
+     || (io.readOpt() == IOobject::READ_IF_PRESENT && headerOk())
+    )
+    {
+        dictionary edgeMeshDict(readStream(typeName));
 
-{}
+        primitiveEdgeMesh::operator=
+        (
+            primitiveEdgeMesh
+            (
+                edgeMeshDict.lookup("points"),
+                edgeMeshDict.lookup("edges")
+            )
+        );
+
+        concaveStart_ = readLabel(edgeMeshDict.lookup("concaveStart"));
+
+        mixedStart_ = readLabel(edgeMeshDict.lookup("mixedStart"));
+
+        nonFeatureStart_ = readLabel(edgeMeshDict.lookup("nonFeatureStart"));
+
+        internalStart_ = readLabel(edgeMeshDict.lookup("internalStart"));
+
+        flatStart_ = readLabel(edgeMeshDict.lookup("flatStart"));
+
+        openStart_ = readLabel(edgeMeshDict.lookup("openStart"));
+
+        multipleStart_ = readLabel(edgeMeshDict.lookup("multipleStart"));
+
+        normals_ = vectorField(edgeMeshDict.lookup("normals"));
+
+        edgeDirections_ = vectorField(edgeMeshDict.lookup("edgeDirections"));
+
+        edgeNormals_ = labelListList(edgeMeshDict.lookup("edgeNormals"));
+
+        featurePointNormals_ = labelListList
+        (
+            edgeMeshDict.lookup("featurePointNormals")
+        );
+
+        regionEdges_ = labelList(edgeMeshDict.lookup("regionEdges"));
+
+        close();
+    }
+}
+
+Foam::featureEdgeMesh::featureEdgeMesh
+(
+    const surfaceFeatures& sFeat,
+    const objectRegistry& obr,
+    const fileName& sFeatFileName,
+    bool write
+)
+:
+    regIOobject
+    (
+        IOobject
+        (
+            sFeatFileName,
+            obr.time().constant(),
+            "featureEdgeMesh",
+            obr,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        )
+    ),
+    primitiveEdgeMesh(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)
+{
+    // 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(5);
+
+    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,"
+                "    bool write"
+                ")"
+            )
+            << 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 primitiveEdgeMesh
+    primitiveEdgeMesh::operator=(primitiveEdgeMesh(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,"
+                "    bool write"
+                ")"
+            )
+                << 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 primitiveEdgeMesh with sorted feature points and
+    // renumbered edges
+    primitiveEdgeMesh::operator=(primitiveEdgeMesh(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;
+
+    // Optionally write the primitiveEdgeMesh to file
+    if (write)
+    {
+        writeObj(sFeatFileName.lessExt());
+    }
+}
 
 
 Foam::featureEdgeMesh::featureEdgeMesh
@@ -48,17 +484,39 @@ Foam::featureEdgeMesh::featureEdgeMesh
     const IOobject& io,
     const pointField& pts,
     const edgeList& eds,
-    const vectorField& normals
+    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
 )
 :
-    edgeMesh(io, pts, eds),
+    regIOobject(io),
+    primitiveEdgeMesh(pts, eds),
+    concaveStart_(concaveStart),
+    mixedStart_(mixedStart),
+    nonFeatureStart_(nonFeatureStart),
+    internalStart_(internalStart),
+    flatStart_(flatStart),
+    openStart_(openStart),
+    multipleStart_(multipleStart),
     normals_(normals),
-    edgeNormals_(edges().size()),
-    allEdges_(identity(edges().size()))
+    edgeDirections_(edgeDirections),
+    edgeNormals_(edgeNormals),
+    featurePointNormals_(featurePointNormals),
+    regionEdges_(regionEdges)
+
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * * Destruct or  * * * * * * * * * * * * * * * //
 
 Foam::featureEdgeMesh::~featureEdgeMesh()
 {}
@@ -66,6 +524,94 @@ 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  * * * * * * * * * * * * * * //
 
@@ -75,7 +621,7 @@ void Foam::featureEdgeMesh::nearestFeatureEdge
     scalarField searchDistSqr,
     labelList& edgeLabel,
     pointField& edgePoint,
-    labelListList& adjacentNormals
+    List<vectorField>& adjacentNormals
 ) const
 {
     edgeLabel.setSize(samples.size());
@@ -98,9 +644,9 @@ void Foam::featureEdgeMesh::nearestFeatureEdge
         }
         else
         {
-            edgeLabel[i] = allEdges_[pHit.index()];
+            edgeLabel[i] = pHit.index();
             edgePoint[i] = pHit.hitPoint();
-            adjacentNormals[i] = edgeNormals_[edgeLabel[i]];
+            adjacentNormals[i] = edgeNormals(edgeLabel[i]);
         }
     }
 }
@@ -123,6 +669,8 @@ Foam::featureEdgeMesh::edgeTree() const
         bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
         bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
 
+        labelList allEdges(identity(edges().size()));
+
         edgeTree_.reset
         (
             new indexedOctree<treeDataEdge>
@@ -132,7 +680,7 @@ Foam::featureEdgeMesh::edgeTree() const
                     false,          // cachebb
                     edges(),        // edges
                     points(),       // points
-                    allEdges_       // selected edges
+                    allEdges       // selected edges
                 ),
                 bb,     // bb
                 8,      // maxLevel
@@ -146,4 +694,166 @@ Foam::featureEdgeMesh::edgeTree() const
 }
 
 
+void Foam::featureEdgeMesh::writeObj
+(
+    const fileName& prefix
+) const
+{
+    Pout<< nl << "Writing featureEdgeMesh components to " << prefix << endl;
+
+    primitiveEdgeMesh::writeObj(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 externalStr(prefix + "_externalEdges.obj");
+    Pout<< "Writing external edges to " << externalStr.name() << endl;
+
+    label verti = 0;
+    for (label i = 0; 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.writeKeyword("concaveStart") << concaveStart_ << token::END_STATEMENT
+        << nl << nl;
+
+    os.writeKeyword("mixedStart") << mixedStart_ << token::END_STATEMENT
+        << nl << nl;
+
+    os.writeKeyword("nonFeatureStart") << nonFeatureStart_
+        << token::END_STATEMENT << nl << nl;
+
+    os.writeKeyword("points") << points() << token::END_STATEMENT
+        << nl << nl;
+
+    os.writeKeyword("internalStart") << internalStart_ << token::END_STATEMENT
+        << nl << nl;
+
+    os.writeKeyword("flatStart") << flatStart_ << token::END_STATEMENT
+        << nl << nl;
+
+    os.writeKeyword("openStart") << openStart_ << token::END_STATEMENT
+        << nl << nl;
+
+    os.writeKeyword("multipleStart") << multipleStart_ << token::END_STATEMENT
+        << nl << nl;
+
+    os.writeKeyword("edges") << edges() << token::END_STATEMENT
+        << nl << nl;
+
+    os.writeKeyword("normals") << normals_ << token::END_STATEMENT
+        << nl << nl;
+
+    os.writeKeyword("edgeDirections") << edgeDirections_ << token::END_STATEMENT
+        << nl << nl;
+
+    os.writeKeyword("edgeNormals") << edgeNormals_ << token::END_STATEMENT
+        << nl << nl;
+
+    os.writeKeyword("featurePointNormals") << featurePointNormals_
+        << token::END_STATEMENT << nl << nl;
+
+    os.writeKeyword("regionEdges") << regionEdges_ << token::END_STATEMENT
+        << nl << nl;
+
+    return os.good();
+}
+
+
 // ************************************************************************* //
diff --git a/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.H b/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.H
index 9addecb458..8b8c72792a 100644
--- a/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.H
+++ b/src/edgeMesh/featureEdgeMesh/featureEdgeMesh.H
@@ -34,7 +34,7 @@ Description
         nonFeatureStart_ .. size-1          : non-feature points
 
     Feature edges are the edgeList of the edgeMesh and are sorted:
-        0 .. internalStart_-1           : external edges (convex)
+        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
@@ -42,6 +42,8 @@ Description
         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
@@ -52,7 +54,10 @@ SourceFiles
 #ifndef featureEdgeMesh_H
 #define featureEdgeMesh_H
 
-#include "edgeMesh.H"
+#include "primitiveEdgeMesh.H"
+#include "surfaceFeatures.H"
+#include "objectRegistry.H"
+#include "IOdictionary.H"
 #include "indexedOctree.H"
 #include "treeDataEdge.H"
 #include "pointIndexHit.H"
@@ -68,31 +73,41 @@ namespace Foam
 
 class featureEdgeMesh
 :
-    public edgeMesh
+    public regIOobject,
+    public primitiveEdgeMesh
 {
 
 public:
 
-        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
-        };
+    TypeName("featureEdgeMesh");
 
-        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 feature edge (consistency with surfaceFeatures)
-        };
+    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_;
+
+
     // Private data
 
         //- Index of the start of the concave feature points
@@ -120,14 +135,17 @@ private:
         //  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_;
 
-        //- Indices of all edges, pre-calculated and stored
-        labelField allEdges_;
+        //- Feature edges which are on the boundary between regions
+        labelList regionEdges_;
 
         //- Search tree for edges.
         mutable autoPtr<indexedOctree<treeDataEdge> > edgeTree_;
@@ -135,6 +153,19 @@ private:
 
     // 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;
+
         //- Disallow default bitwise copy construct
         featureEdgeMesh(const featureEdgeMesh&);
 
@@ -149,13 +180,35 @@ public:
         //- Construct (read) given an IOobject
         featureEdgeMesh(const IOobject&);
 
-        //- Construct from edgeMesh data and
+         //- 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,
+            bool write = false
+        );
+
+        //- Construct from components
         featureEdgeMesh
         (
             const IOobject& io,
             const pointField& pts,
             const edgeList& eds,
-            const vectorField& normals
+            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
         );
 
 
@@ -174,15 +227,48 @@ public:
                 scalarField searchDistSqr,
                 labelList& edgeLabel,
                 pointField& edgePoint,
-                labelListList& adjacentNormals
+                List<vectorField>& adjacentNormals
             ) const;
 
         // Access
 
+            //- 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 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 contruction of octree for boundary edges
             const indexedOctree<treeDataEdge>& edgeTree() const;
 
@@ -191,6 +277,15 @@ public:
         // Edit
 
         // Write
+
+            //- Write all components of the featureEdgeMesh as obj files
+            void writeObj(const fileName& fName) const;
+
+            // //- ReadData function required for regIOobject read operation
+            // virtual bool readData(Istream&);
+
+            //- WriteData function required for regIOobject write operation
+            virtual bool writeData(Ostream&) const;
 };
 
 
diff --git a/src/edgeMesh/featureEdgeMesh/featureEdgeMeshI.H b/src/edgeMesh/featureEdgeMesh/featureEdgeMeshI.H
index 42bc929311..06cd6e8040 100644
--- a/src/edgeMesh/featureEdgeMesh/featureEdgeMeshI.H
+++ b/src/edgeMesh/featureEdgeMesh/featureEdgeMeshI.H
@@ -26,34 +26,163 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+bool Foam::featureEdgeMesh::featurePoint(label ptI) const
+{
+    return ptI < nonFeatureStart_;
+}
 
 
-// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
+const Foam::vectorField& Foam::featureEdgeMesh::normals() const
+{
+    return normals_;
+}
+
+const Foam::vectorField& Foam::featureEdgeMesh::edgeDirections() const
+{
+    return edgeDirections_;
+}
 
 
-// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+const Foam::labelListList& Foam::featureEdgeMesh::edgeNormals() const
+{
+    return edgeNormals_;
+}
 
 
-// * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * * //
+Foam::vectorField Foam::featureEdgeMesh::edgeNormals
+(
+    const labelList& edgeNormIs) const
+{
+    vectorField norms(edgeNormIs.size());
+
+    forAll(edgeNormIs, i)
+    {
+        norms[i] = normals_[edgeNormIs[i]];
+    }
+
+    return norms;
+}
 
 
-// * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
+Foam::vectorField Foam::featureEdgeMesh::edgeNormals(label edgeI) const
+{
+    return edgeNormals(edgeNormals_[edgeI]);
+}
+
+
+const Foam::labelListList& Foam::featureEdgeMesh::featurePointNormals() const
+{
+    return featurePointNormals_;
+}
+
+
+Foam::vectorField Foam::featureEdgeMesh::featurePointNormals(label ptI) const
+{
+    if (!featurePoint(ptI))
+    {
+        WarningIn("vectorField featureEdgeMesh::featurePointNormals")
+            << "Requesting the normals of a non-feature point. "
+            << "Returned zero length vectorField."
+            << endl;
+
+        return vectorField(0);
+    }
+
+    labelList featPtNormIs(featurePointNormals_[ptI]);
+
+    vectorField norms(featPtNormIs.size());
+
+    forAll(featPtNormIs, i)
+    {
+        norms[i] = normals_[featPtNormIs[i]];
+    }
+
+    return norms;
+}
+
+
+const Foam::labelList& Foam::featureEdgeMesh::regionEdges() const
+{
+    return regionEdges_;
+}
+
+
+Foam::featureEdgeMesh::pointStatus Foam::featureEdgeMesh::getPointStatus
+(
+    label ptI
+) const
+{
+    if (ptI < concaveStart_)
+    {
+        return CONVEX;
+    }
+    else if (ptI < mixedStart_)
+    {
+        return CONCAVE;
+    }
+    else if (ptI < nonFeatureStart_)
+    {
+        return MIXED;
+    }
+    else
+    {
+        return NONFEATURE;
+    }
+}
+
+
+Foam::featureEdgeMesh::edgeStatus Foam::featureEdgeMesh::getEdgeStatus
+(
+    label edgeI
+) const
+{
+    if (edgeI < internalStart_)
+    {
+        return EXTERNAL;
+    }
+    else if (edgeI < flatStart_)
+    {
+        return INTERNAL;
+    }
+    else if (edgeI < openStart_)
+    {
+        return FLAT;
+    }
+    else if (edgeI < multipleStart_)
+    {
+        return OPEN;
+    }
+    else
+    {
+        return MULTIPLE;
+    }
+}
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 
 // ************************************************************************* //
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/src/edgeMesh/primitiveEdgeMesh/primitiveEdgeMesh.C b/src/edgeMesh/primitiveEdgeMesh/primitiveEdgeMesh.C
index 807ea3bea3..9e9bfd8749 100644
--- a/src/edgeMesh/primitiveEdgeMesh/primitiveEdgeMesh.C
+++ b/src/edgeMesh/primitiveEdgeMesh/primitiveEdgeMesh.C
@@ -27,6 +27,9 @@ License
 #include "primitiveEdgeMesh.H"
 #include "mergePoints.H"
 #include "StaticHashTable.H"
+#include "demandDrivenData.H"
+#include "meshTools.H"
+#include "OFstream.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -98,6 +101,12 @@ Foam::primitiveEdgeMesh::primitiveEdgeMesh(const primitiveEdgeMesh& em)
 {}
 
 
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::primitiveEdgeMesh::~primitiveEdgeMesh()
+{}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::label Foam::primitiveEdgeMesh::regions(labelList& edgeRegion) const
@@ -244,13 +253,28 @@ void Foam::primitiveEdgeMesh::mergePoints(const scalar mergeDist)
 }
 
 
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
-void Foam::primitiveEdgeMesh::operator=(const primitiveEdgeMesh& rhs)
+void Foam::primitiveEdgeMesh::writeObj
+(
+    const fileName& fName
+) const
 {
-    points_ = rhs.points_;
-    edges_ = rhs.edges_;
-    pointEdgesPtr_.reset(NULL);
+    Pout<< nl << "Writing points and edges to " << fName << endl;
+
+    OFstream str(fName);
+
+    forAll(points_, p)
+    {
+        meshTools::writeOBJ(str, points_[p]);
+    }
+
+    forAll (edges_, e)
+    {
+        const edge& ed = edges_[e];
+
+        str<< "l " << ed.start() + 1 << ' ' << ed.end() + 1;
+
+        str<< nl;
+    }
 }
 
 
diff --git a/src/edgeMesh/primitiveEdgeMesh/primitiveEdgeMesh.H b/src/edgeMesh/primitiveEdgeMesh/primitiveEdgeMesh.H
index 9c33b509f4..1bafdbb7fe 100644
--- a/src/edgeMesh/primitiveEdgeMesh/primitiveEdgeMesh.H
+++ b/src/edgeMesh/primitiveEdgeMesh/primitiveEdgeMesh.H
@@ -85,6 +85,10 @@ public:
         primitiveEdgeMesh(const primitiveEdgeMesh&);
 
 
+    // Destructor
+
+        ~primitiveEdgeMesh();
+
 
     // Member Functions
 
@@ -101,6 +105,9 @@ public:
         //- Merge common points (points within mergeDist)
         void mergePoints(const scalar mergeDist);
 
+        //- Write the primitiveEdgeMesh as an obj file
+        void writeObj(const fileName& fName) const;
+
     // Member Operators
 
         inline void operator=(const primitiveEdgeMesh&);
diff --git a/src/edgeMesh/primitiveEdgeMesh/primitiveEdgeMeshI.H b/src/edgeMesh/primitiveEdgeMesh/primitiveEdgeMeshI.H
index 5d944bcbdc..be2cdc9a5b 100644
--- a/src/edgeMesh/primitiveEdgeMesh/primitiveEdgeMeshI.H
+++ b/src/edgeMesh/primitiveEdgeMesh/primitiveEdgeMeshI.H
@@ -48,4 +48,14 @@ inline const Foam::labelListList& Foam::primitiveEdgeMesh::pointEdges() const
 }
 
 
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+void Foam::primitiveEdgeMesh::operator=(const primitiveEdgeMesh& rhs)
+{
+    points_ = rhs.points_;
+    edges_ = rhs.edges_;
+    pointEdgesPtr_.reset(NULL);
+}
+
+
 // ************************************************************************* //