diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
index 496f980a11..4e0f1d160b 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
@@ -310,14 +310,13 @@ class autoSnapDriver
                 //- Find point on nearest feature edge (within searchDist).
                 //  Return point and feature
                 //  and store feature-edge to mesh-point and vice versa
-                pointIndexHit findNearFeatureEdge
+                Tuple2<label, pointIndexHit> findNearFeatureEdge
                 (
                     const indirectPrimitivePatch& pp,
                     const scalarField& snapDist,
                     const label pointI,
                     const point& estimatedPt,
 
-                    label& featI,
                     List<List<DynamicList<point> > >&,
                     List<List<DynamicList<pointConstraint> > >&,
                     vectorField&,
@@ -330,7 +329,7 @@ class autoSnapDriver
                 //  If another mesh point already referring to this feature
                 //  point and further away, reset that one to a near feature
                 //  edge (using findNearFeatureEdge above)
-                labelPair findNearFeaturePoint
+                Tuple2<label, pointIndexHit> findNearFeaturePoint
                 (
                     const indirectPrimitivePatch& pp,
                     const scalarField& snapDist,
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
index 67fee3ef96..03615a76bd 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
@@ -170,7 +170,7 @@ void Foam::autoSnapDriver::smoothAndConstrain
 
                     if (isMasterEdge[meshEdges[edgeI]])
                     {
-                        label nbrPointI = edges[pEdges[i]].otherVertex(pointI);
+                        label nbrPointI = edges[edgeI].otherVertex(pointI);
                         if (constraints[nbrPointI].first() >= nConstraints)
                         {
                             dispSum[pointI] += disp[nbrPointI];
@@ -213,80 +213,6 @@ void Foam::autoSnapDriver::smoothAndConstrain
         }
     }
 }
-//XXXXXX
-//TODO: make proper parallel so coupled edges don't have double influence
-//void Foam::autoSnapDriver::smoothAndConstrain2
-//(
-//    const bool applyConstraints,
-//    const indirectPrimitivePatch& pp,
-//    const List<pointConstraint>& constraints,
-//    vectorField& disp
-//) const
-//{
-//    const fvMesh& mesh = meshRefiner_.mesh();
-//
-//    for (label avgIter = 0; avgIter < 20; avgIter++)
-//    {
-//        vectorField dispSum(pp.nPoints(), vector::zero);
-//        labelList dispCount(pp.nPoints(), 0);
-//
-//        const labelListList& pointEdges = pp.pointEdges();
-//        const edgeList& edges = pp.edges();
-//
-//        forAll(pointEdges, pointI)
-//        {
-//            const labelList& pEdges = pointEdges[pointI];
-//
-//            forAll(pEdges, i)
-//            {
-//                label nbrPointI = edges[pEdges[i]].otherVertex(pointI);
-//                dispSum[pointI] += disp[nbrPointI];
-//                dispCount[pointI]++;
-//            }
-//        }
-//
-//        syncTools::syncPointList
-//        (
-//            mesh,
-//            pp.meshPoints(),
-//            dispSum,
-//            plusEqOp<point>(),
-//            vector::zero,
-//            mapDistribute::transform()
-//        );
-//        syncTools::syncPointList
-//        (
-//            mesh,
-//            pp.meshPoints(),
-//            dispCount,
-//            plusEqOp<label>(),
-//            0,
-//            mapDistribute::transform()
-//        );
-//
-//        // Constraints
-//        forAll(constraints, pointI)
-//        {
-//            if (dispCount[pointI] > 0)// && constraints[pointI].first() <= 1)
-//            {
-//                // Mix my displacement with neighbours' displacement
-//                disp[pointI] =
-//                    0.5
-//                   *(disp[pointI] + dispSum[pointI]/dispCount[pointI]);
-//
-//                if (applyConstraints)
-//                {
-//                    disp[pointI] = transform
-//                    (
-//                        constraints[pointI].constraintTransformation(),
-//                        disp[pointI]
-//                    );
-//                }
-//            }
-//        }
-//    }
-//}
-//XXXXXX
 
 
 void Foam::autoSnapDriver::calcNearestFace
@@ -716,7 +642,7 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
         pNormals = List<point>(pNormals, visitOrder);
         pDisp = List<point>(pDisp, visitOrder);
         pFc = List<point>(pFc, visitOrder);
-        pFid = UIndirectList<label>(pFid, visitOrder);
+        pFid = UIndirectList<label>(pFid, visitOrder)();
     }
 }
 
@@ -1412,14 +1338,14 @@ Foam::labelPair Foam::autoSnapDriver::findDiagonalAttraction
 }
 
 
-Foam::pointIndexHit Foam::autoSnapDriver::findNearFeatureEdge
+Foam::Tuple2<Foam::label, Foam::pointIndexHit>
+Foam::autoSnapDriver::findNearFeatureEdge
 (
     const indirectPrimitivePatch& pp,
     const scalarField& snapDist,
     const label pointI,
     const point& estimatedPt,
 
-    label& featI,
     List<List<DynamicList<point> > >& edgeAttractors,
     List<List<DynamicList<pointConstraint> > >& edgeConstraints,
     vectorField& patchAttraction,
@@ -1430,16 +1356,18 @@ Foam::pointIndexHit Foam::autoSnapDriver::findNearFeatureEdge
 
     labelList nearEdgeFeat;
     List<pointIndexHit> nearEdgeInfo;
+    vectorField nearNormal;
     features.findNearestEdge
     (
         pointField(1, estimatedPt),
         scalarField(1, sqr(snapDist[pointI])),
         nearEdgeFeat,
-        nearEdgeInfo
+        nearEdgeInfo,
+        nearNormal
     );
 
     const pointIndexHit& nearInfo = nearEdgeInfo[0];
-    featI = nearEdgeFeat[0];
+    label featI = nearEdgeFeat[0];
 
     if (nearInfo.hit())
     {
@@ -1449,12 +1377,7 @@ Foam::pointIndexHit Foam::autoSnapDriver::findNearFeatureEdge
         (
             nearInfo.hitPoint()
         );
-        pointConstraint c;
-        const edge e = features[featI].edges()[nearInfo.index()];
-        vector eVec = e.vec(features[featI].points());
-        eVec /= mag(eVec)+VSMALL;
-        c.first() = 2;
-        c.second() = eVec;
+        pointConstraint c(Tuple2<label, vector>(2, nearNormal[0]));
         edgeConstraints[featI][nearInfo.index()].append(c);
 
         // Store for later use
@@ -1462,11 +1385,12 @@ Foam::pointIndexHit Foam::autoSnapDriver::findNearFeatureEdge
             nearInfo.hitPoint()-pp.localPoints()[pointI];
         patchConstraints[pointI] = c;
     }
-    return nearInfo;
+    return Tuple2<label, pointIndexHit>(featI, nearInfo);
 }
 
 
-Foam::labelPair Foam::autoSnapDriver::findNearFeaturePoint
+Foam::Tuple2<Foam::label, Foam::pointIndexHit>
+Foam::autoSnapDriver::findNearFeaturePoint
 (
     const indirectPrimitivePatch& pp,
     const scalarField& snapDist,
@@ -1487,26 +1411,23 @@ Foam::labelPair Foam::autoSnapDriver::findNearFeaturePoint
     const refinementFeatures& features = meshRefiner_.features();
 
     labelList nearFeat;
-    labelList nearIndex;
+    List<pointIndexHit> nearInfo;
     features.findNearestPoint
     (
         pointField(1, estimatedPt),
         scalarField(1, sqr(snapDist[pointI])),
         nearFeat,
-        nearIndex
+        nearInfo
     );
 
     label featI = nearFeat[0];
-    label featPointI = -1;
 
     if (featI != -1)
     {
         const point& pt = pp.localPoints()[pointI];
 
-        const treeDataPoint& shapes =
-            features.pointTrees()[featI].shapes();
-        featPointI = shapes.pointLabels()[nearIndex[0]];
-        const point& featPt = shapes.points()[featPointI];
+        label featPointI = nearInfo[0].index();
+        const point& featPt = nearInfo[0].hitPoint();
         scalar distSqr = magSqr(featPt-pt);
 
         // Check if already attracted
@@ -1537,7 +1458,6 @@ Foam::labelPair Foam::autoSnapDriver::findNearFeaturePoint
                 patchAttraction[oldPointI] = vector::zero;
                 patchConstraints[oldPointI] = pointConstraint();
 
-                label edgeFeatI;
                 findNearFeatureEdge
                 (
                     pp,
@@ -1545,7 +1465,6 @@ Foam::labelPair Foam::autoSnapDriver::findNearFeaturePoint
                     oldPointI,
                     pp.localPoints()[oldPointI],
 
-                    edgeFeatI,
                     edgeAttractors,
                     edgeConstraints,
                     patchAttraction,
@@ -1566,7 +1485,7 @@ Foam::labelPair Foam::autoSnapDriver::findNearFeaturePoint
         }
     }
 
-    return labelPair(featI, featPointI);
+    return Tuple2<label, pointIndexHit>(featI, nearInfo[0]);
 }
 
 
@@ -1637,7 +1556,6 @@ void Foam::autoSnapDriver::determineFeatures
             << featurePointStr().name() << endl;
     }
 
-    const refinementFeatures& features = meshRefiner_.features();
 
     forAll(pp.localPoints(), pointI)
     {
@@ -1664,6 +1582,7 @@ void Foam::autoSnapDriver::determineFeatures
             constraint
         );
 
+
         if
         (
             (constraint.first() > patchConstraints[pointI].first())
@@ -1696,15 +1615,14 @@ void Foam::autoSnapDriver::determineFeatures
                         // Behave like when having two surface normals so
                         // attract to nearest feature edge (with a guess for
                         // the multipatch point as starting point)
-                        label featI = -1;
-                        pointIndexHit nearInfo = findNearFeatureEdge
+                        Tuple2<label, pointIndexHit> nearInfo =
+                        findNearFeatureEdge
                         (
                             pp,
                             snapDist,
                             pointI,
                             multiPatchPt.hitPoint(),        //estimatedPt
 
-                            featI,
                             edgeAttractors,
                             edgeConstraints,
 
@@ -1712,14 +1630,15 @@ void Foam::autoSnapDriver::determineFeatures
                             patchConstraints
                         );
 
-                        if (nearInfo.hit())
+                        const pointIndexHit& info = nearInfo.second();
+                        if (info.hit())
                         {
                             // Dump
                             if (featureEdgeStr.valid())
                             {
                                 featureEdgeStr().write
                                 (
-                                    linePointRef(pt, nearInfo.hitPoint())
+                                    linePointRef(pt, info.hitPoint())
                                 );
                             }
                         }
@@ -1729,7 +1648,7 @@ void Foam::autoSnapDriver::determineFeatures
                             {
                                 missedEdgeStr().write
                                 (
-                                    linePointRef(pt, nearInfo.missPoint())
+                                    linePointRef(pt, info.missPoint())
                                 );
                             }
                         }
@@ -1746,15 +1665,13 @@ void Foam::autoSnapDriver::determineFeatures
                 // Determine nearest point on feature edge. Store constraint
                 // (calculated from feature edge, alternative would be to
                 //  use constraint calculated from both surfaceNormals)
-                label featI = -1;
-                pointIndexHit nearInfo = findNearFeatureEdge
+                Tuple2<label, pointIndexHit> nearInfo = findNearFeatureEdge
                 (
                     pp,
                     snapDist,
                     pointI,
                     estimatedPt,
 
-                    featI,
                     edgeAttractors,
                     edgeConstraints,
 
@@ -1762,14 +1679,15 @@ void Foam::autoSnapDriver::determineFeatures
                     patchConstraints
                 );
 
-                if (nearInfo.hit())
+                // Dump to obj
+                const pointIndexHit& info = nearInfo.second();
+                if (info.hit())
                 {
-                    // Dump
                     if (featureEdgeStr.valid())
                     {
                         featureEdgeStr().write
                         (
-                            linePointRef(pt, nearInfo.hitPoint())
+                            linePointRef(pt, info.hitPoint())
                         );
                     }
                 }
@@ -1779,7 +1697,7 @@ void Foam::autoSnapDriver::determineFeatures
                     {
                         missedEdgeStr().write
                         (
-                            linePointRef(pt, nearInfo.missPoint())
+                            linePointRef(pt, info.missPoint())
                         );
                     }
                 }
@@ -1789,7 +1707,7 @@ void Foam::autoSnapDriver::determineFeatures
                 // Mark point on the nearest feature point.
                 const point estimatedPt(pt + patchAttraction[pointI]);
 
-                labelPair nearInfo = findNearFeaturePoint
+                Tuple2<label, pointIndexHit> nearInfo = findNearFeaturePoint
                 (
                     pp,
                     snapDist,
@@ -1807,18 +1725,10 @@ void Foam::autoSnapDriver::determineFeatures
                     patchConstraints
                 );
 
-                if (nearInfo.first() != -1)
+                const pointIndexHit& info = nearInfo.second();
+                if (info.hit() && featurePointStr.valid())
                 {
-                    // Dump
-                    if (featurePointStr.valid())
-                    {
-                        const treeDataPoint& shapes =
-                            features.pointTrees()[nearInfo.first()].shapes();
-                        const point& featPt =
-                            shapes.points()[nearInfo.second()];
-
-                        featurePointStr().write(linePointRef(pt, featPt));
-                    }
+                    featurePointStr().write(linePointRef(pt, info.hitPoint()));
                 }
             }
         }
@@ -2029,22 +1939,20 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
 
             if (pointStatus[pointI] == 0)   // baffle edge
             {
-                label featI;
-                const pointIndexHit nearInfo = findNearFeatureEdge
+                Tuple2<label, pointIndexHit> nearInfo = findNearFeatureEdge
                 (
                     pp,
                     snapDist,
                     pointI,
                     pt,
 
-                    featI,
                     edgeAttractors,
                     edgeConstraints,
                     rawPatchAttraction,
                     rawPatchConstraints
                 );
 
-                if (!nearInfo.hit())
+                if (!nearInfo.second().hit())
                 {
                     //Pout<< "*** Failed to find close edge to point " << pt
                     //    << endl;
@@ -2053,23 +1961,21 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
             else if (pointStatus[pointI] == 1)   // baffle point
             {
                 labelList nearFeat;
-                labelList nearIndex;
+                List<pointIndexHit> nearInfo;
                 features.findNearestPoint
                 (
                     pointField(1, pt),
                     scalarField(1, sqr(snapDist[pointI])),
                     nearFeat,
-                    nearIndex
+                    nearInfo
                 );
 
                 label featI = nearFeat[0];
 
                 if (featI != -1)
                 {
-                    const treeDataPoint& shapes =
-                        features.pointTrees()[featI].shapes();
-                    label featPointI = shapes.pointLabels()[nearIndex[0]];
-                    const point& featPt = shapes.points()[featPointI];
+                    label featPointI = nearInfo[0].index();
+                    const point& featPt = nearInfo[0].hitPoint();
                     scalar distSqr = magSqr(featPt-pt);
 
                     // Check if already attracted
@@ -2099,7 +2005,6 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                             // The current point is closer so wins. Reset
                             // the old point to attract to nearest edge
                             // instead.
-                            label edgeFeatI;
                             findNearFeatureEdge
                             (
                                 pp,
@@ -2107,7 +2012,6 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                                 oldPointI,
                                 pp.localPoints()[oldPointI],
 
-                                edgeFeatI,
                                 edgeAttractors,
                                 edgeConstraints,
                                 rawPatchAttraction,
@@ -2130,7 +2034,6 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                     //    << " for baffle-feature-point " << pt
                     //    << endl;
 
-                    label featI;
                     findNearFeatureEdge
                     (
                         pp,
@@ -2138,7 +2041,6 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                         pointI,
                         pt,                     // starting point
 
-                        featI,
                         edgeAttractors,
                         edgeConstraints,
                         rawPatchAttraction,
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C b/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C
index a5109bc3c8..f2b730bc9d 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C
+++ b/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C
@@ -26,6 +26,7 @@ License
 #include "refinementFeatures.H"
 #include "Time.H"
 #include "Tuple2.H"
+#include "DynamicField.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -41,7 +42,40 @@ void Foam::refinementFeatures::read
 
         fileName featFileName(dict.lookup("file"));
 
+
+        // Try reading extendedEdgeMesh first
+
+        IOobject extFeatObj
+        (
+            featFileName,                       // name
+            io.time().constant(),               // instance
+            "extendedFeatureEdgeMesh",          // local
+            io.time(),                          // registry
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false
+        );
+
+        const fileName fName(extFeatObj.filePath());
+
+        if (!fName.empty() && extendedEdgeMesh::canRead(fName))
         {
+            autoPtr<extendedEdgeMesh> eMeshPtr = extendedEdgeMesh::New
+            (
+                fName
+            );
+
+            Info<< "Read extendedFeatureEdgeMesh " << extFeatObj.name()
+                << nl << incrIndent;
+            eMeshPtr().writeStats(Info);
+            Info<< decrIndent << endl;
+
+            set(featI, new extendedFeatureEdgeMesh(extFeatObj, eMeshPtr()));
+        }
+        else
+        {
+            // Try reading edgeMesh
+
             IOobject featObj
             (
                 featFileName,                       // name
@@ -53,21 +87,106 @@ void Foam::refinementFeatures::read
                 false
             );
 
-            autoPtr<edgeMesh> eMeshPtr = edgeMesh::New(featObj.filePath());
+            const fileName fName(featObj.filePath());
 
-            set
-            (
-                featI,
-                new featureEdgeMesh
+            if (fName.empty())
+            {
+                FatalIOErrorIn
                 (
-                    featObj,
-                    eMeshPtr->points(),
-                    eMeshPtr->edges()
-                )
+                    "refinementFeatures::read"
+                    "(const objectRegistry&"
+                    ", const PtrList<dictionary>&)",
+                    dict
+                )   << "Could not open " << featObj.objectPath()
+                    << exit(FatalIOError);
+            }
+
+
+            // Read as edgeMesh
+            autoPtr<edgeMesh> eMeshPtr = edgeMesh::New(fName);
+            const edgeMesh& eMesh = eMeshPtr();
+
+            Info<< "Read edgeMesh " << featObj.name() << nl
+                << incrIndent;
+            eMesh.writeStats(Info);
+            Info<< decrIndent << endl;
+
+
+            // Analyse for feature points. These are all classified as mixed
+            // points for lack of anything better
+            const labelListList& pointEdges = eMesh.pointEdges();
+
+            labelList oldToNew(eMesh.points().size(), -1);
+            DynamicField<point> newPoints(eMesh.points().size());
+            forAll(pointEdges, pointI)
+            {
+                if (pointEdges[pointI].size() > 2)
+                {
+                    oldToNew[pointI] = newPoints.size();
+                    newPoints.append(eMesh.points()[pointI]);
+                }
+                //else if (pointEdges[pointI].size() == 2)
+                //MEJ: do something based on a feature angle?
+            }
+            label nFeatures = newPoints.size();
+            forAll(oldToNew, pointI)
+            {
+                if (oldToNew[pointI] == -1)
+                {
+                    oldToNew[pointI] = newPoints.size();
+                    newPoints.append(eMesh.points()[pointI]);
+                }
+            }
+
+
+            const edgeList& edges = eMesh.edges();
+            edgeList newEdges(edges.size());
+            forAll(edges, edgeI)
+            {
+                const edge& e = edges[edgeI];
+                newEdges[edgeI] = edge
+                (
+                    oldToNew[e[0]],
+                    oldToNew[e[1]]
+                );
+            }
+
+            // Construct an extendedEdgeMesh with
+            // - all points on more than 2 edges : mixed feature points
+            // - all edges : external edges
+
+            extendedEdgeMesh eeMesh
+            (
+                newPoints,          // pts
+                newEdges,           // eds
+                0,                  // (point) concaveStart
+                0,                  // (point) mixedStart
+                nFeatures,          // (point) nonFeatureStart
+                edges.size(),       // (edge) internalStart
+                edges.size(),       // (edge) flatStart
+                edges.size(),       // (edge) openStart
+                edges.size(),       // (edge) multipleStart
+                vectorField(0),     // normals
+                List<extendedEdgeMesh::sideVolumeType>(0),// normalVolumeTypes
+                vectorField(0),     // edgeDirections
+                labelListList(0),   // normalDirections
+                labelListList(0),   // edgeNormals
+                labelListList(0),   // featurePointNormals
+                labelListList(0),   // featurePointEdges
+                labelList(0)        // regionEdges
             );
+
+
+            Info<< "Constructed extendedFeatureEdgeMesh " << featObj.name()
+                << nl << incrIndent;
+            eeMesh.writeStats(Info);
+            Info<< decrIndent << endl;
+
+
+            set(featI, new extendedFeatureEdgeMesh(featObj, eeMesh));
         }
 
-        const featureEdgeMesh& eMesh = operator[](featI);
+        const edgeMesh& eMesh = operator[](featI);
 
         //eMesh.mergePoints(meshRefiner_.mergeDistance());
 
@@ -145,7 +264,7 @@ void Foam::refinementFeatures::buildTrees
     const labelList& featurePoints
 )
 {
-    const featureEdgeMesh& eMesh = operator[](featI);
+    const edgeMesh& eMesh = operator[](featI);
     const pointField& points = eMesh.points();
     const edgeList& edges = eMesh.edges();
 
@@ -277,7 +396,7 @@ Foam::refinementFeatures::refinementFeatures
     const PtrList<dictionary>& featDicts
 )
 :
-    PtrList<featureEdgeMesh>(featDicts.size()),
+    PtrList<extendedFeatureEdgeMesh>(featDicts.size()),
     distances_(featDicts.size()),
     levels_(featDicts.size()),
     edgeTrees_(featDicts.size()),
@@ -289,7 +408,7 @@ Foam::refinementFeatures::refinementFeatures
     // Search engines
     forAll(*this, i)
     {
-        const featureEdgeMesh& eMesh = operator[](i);
+        const extendedEdgeMesh& eMesh = operator[](i);
         const labelListList& pointEdges = eMesh.pointEdges();
 
         DynamicList<label> featurePoints;
@@ -303,81 +422,82 @@ Foam::refinementFeatures::refinementFeatures
 
         Info<< "Detected " << featurePoints.size()
             << " featurePoints out of " << pointEdges.size()
-            << " points on feature " << eMesh.name() << endl;
+            << " points on feature " << operator[](i).name()
+            << endl;
 
         buildTrees(i, featurePoints);
     }
 }
 
 
-Foam::refinementFeatures::refinementFeatures
-(
-    const objectRegistry& io,
-    const PtrList<dictionary>& featDicts,
-    const scalar minCos
-)
-:
-    PtrList<featureEdgeMesh>(featDicts.size()),
-    distances_(featDicts.size()),
-    levels_(featDicts.size()),
-    edgeTrees_(featDicts.size()),
-    pointTrees_(featDicts.size())
-{
-    // Read features
-    read(io, featDicts);
-
-    // Search engines
-    forAll(*this, i)
-    {
-        const featureEdgeMesh& eMesh = operator[](i);
-        const pointField& points = eMesh.points();
-        const edgeList& edges = eMesh.edges();
-        const labelListList& pointEdges = eMesh.pointEdges();
-
-        DynamicList<label> featurePoints;
-        forAll(pointEdges, pointI)
-        {
-            const labelList& pEdges = pointEdges[pointI];
-            if (pEdges.size() > 2)
-            {
-                featurePoints.append(pointI);
-            }
-            else if (pEdges.size() == 2)
-            {
-                // Check the angle
-                const edge& e0 = edges[pEdges[0]];
-                const edge& e1 = edges[pEdges[1]];
-
-                const point& p = points[pointI];
-                const point& p0 = points[e0.otherVertex(pointI)];
-                const point& p1 = points[e1.otherVertex(pointI)];
-
-                vector v0 = p-p0;
-                scalar v0Mag = mag(v0);
-
-                vector v1 = p1-p;
-                scalar v1Mag = mag(v1);
-
-                if
-                (
-                    v0Mag > SMALL
-                 && v1Mag > SMALL
-                 && ((v0/v0Mag & v1/v1Mag) < minCos)
-                )
-                {
-                    featurePoints.append(pointI);
-                }
-            }
-        }
-
-        Info<< "Detected " << featurePoints.size()
-            << " featurePoints out of " << points.size()
-            << " points on feature " << eMesh.name()
-            << " when using feature cos " << minCos << endl;
-
-        buildTrees(i, featurePoints);
-    }
-}
+//Foam::refinementFeatures::refinementFeatures
+//(
+//    const objectRegistry& io,
+//    const PtrList<dictionary>& featDicts,
+//    const scalar minCos
+//)
+//:
+//    PtrList<extendedFeatureEdgeMesh>(featDicts.size()),
+//    distances_(featDicts.size()),
+//    levels_(featDicts.size()),
+//    edgeTrees_(featDicts.size()),
+//    pointTrees_(featDicts.size())
+//{
+//    // Read features
+//    read(io, featDicts);
+//
+//    // Search engines
+//    forAll(*this, i)
+//    {
+//        const edgeMesh& eMesh = operator[](i);
+//        const pointField& points = eMesh.points();
+//        const edgeList& edges = eMesh.edges();
+//        const labelListList& pointEdges = eMesh.pointEdges();
+//
+//        DynamicList<label> featurePoints;
+//        forAll(pointEdges, pointI)
+//        {
+//            const labelList& pEdges = pointEdges[pointI];
+//            if (pEdges.size() > 2)
+//            {
+//                featurePoints.append(pointI);
+//            }
+//            else if (pEdges.size() == 2)
+//            {
+//                // Check the angle
+//                const edge& e0 = edges[pEdges[0]];
+//                const edge& e1 = edges[pEdges[1]];
+//
+//                const point& p = points[pointI];
+//                const point& p0 = points[e0.otherVertex(pointI)];
+//                const point& p1 = points[e1.otherVertex(pointI)];
+//
+//                vector v0 = p-p0;
+//                scalar v0Mag = mag(v0);
+//
+//                vector v1 = p1-p;
+//                scalar v1Mag = mag(v1);
+//
+//                if
+//                (
+//                    v0Mag > SMALL
+//                 && v1Mag > SMALL
+//                 && ((v0/v0Mag & v1/v1Mag) < minCos)
+//                )
+//                {
+//                    featurePoints.append(pointI);
+//                }
+//            }
+//        }
+//
+//        Info<< "Detected " << featurePoints.size()
+//            << " featurePoints out of " << points.size()
+//            << " points on feature " << i   //eMesh.name()
+//            << " when using feature cos " << minCos << endl;
+//
+//        buildTrees(i, featurePoints);
+//    }
+//}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -387,12 +507,16 @@ void Foam::refinementFeatures::findNearestEdge
     const pointField& samples,
     const scalarField& nearestDistSqr,
     labelList& nearFeature,
-    List<pointIndexHit>& nearInfo
+    List<pointIndexHit>& nearInfo,
+    vectorField& nearNormal
 ) const
 {
     nearFeature.setSize(samples.size());
     nearFeature = -1;
     nearInfo.setSize(samples.size());
+    nearInfo = pointIndexHit();
+    nearNormal.setSize(samples.size());
+    nearNormal = vector::zero;
 
     forAll(edgeTrees_, featI)
     {
@@ -418,8 +542,18 @@ void Foam::refinementFeatures::findNearestEdge
 
                 if (info.hit())
                 {
-                    nearInfo[sampleI] = info;
                     nearFeature[sampleI] = featI;
+                    nearInfo[sampleI] = pointIndexHit
+                    (
+                        info.hit(),
+                        info.hitPoint(),
+                        tree.shapes().edgeLabels()[info.index()]
+                    );
+
+                    const treeDataEdge& td = tree.shapes();
+                    const edge& e = td.edges()[nearInfo[sampleI].index()];
+                    nearNormal[sampleI] =  e.vec(td.points());
+                    nearNormal[sampleI] /= mag(nearNormal[sampleI])+VSMALL;
                 }
             }
         }
@@ -427,18 +561,71 @@ void Foam::refinementFeatures::findNearestEdge
 }
 
 
+//void Foam::refinementFeatures::findNearestPoint
+//(
+//    const pointField& samples,
+//    const scalarField& nearestDistSqr,
+//    labelList& nearFeature,
+//    labelList& nearIndex
+//) const
+//{
+//    nearFeature.setSize(samples.size());
+//    nearFeature = -1;
+//    nearIndex.setSize(samples.size());
+//    nearIndex = -1;
+//
+//    forAll(pointTrees_, featI)
+//    {
+//        const indexedOctree<treeDataPoint>& tree = pointTrees_[featI];
+//
+//        if (tree.shapes().pointLabels().size() > 0)
+//        {
+//            forAll(samples, sampleI)
+//            {
+//                const point& sample = samples[sampleI];
+//
+//                scalar distSqr;
+//                if (nearFeature[sampleI] != -1)
+//                {
+//                    label nearFeatI = nearFeature[sampleI];
+//                    const indexedOctree<treeDataPoint>& nearTree =
+//                        pointTrees_[nearFeatI];
+//                    label featPointI =
+//                        nearTree.shapes().pointLabels()[nearIndex[sampleI]];
+//                    const point& featPt =
+//                        operator[](nearFeatI).points()[featPointI];
+//                    distSqr = magSqr(featPt-sample);
+//                }
+//                else
+//                {
+//                    distSqr = nearestDistSqr[sampleI];
+//                }
+//
+//                pointIndexHit info = tree.findNearest(sample, distSqr);
+//
+//                if (info.hit())
+//                {
+//                    nearFeature[sampleI] = featI;
+//                    nearIndex[sampleI] = info.index();
+//                }
+//            }
+//        }
+//    }
+//}
+
+
 void Foam::refinementFeatures::findNearestPoint
 (
     const pointField& samples,
     const scalarField& nearestDistSqr,
     labelList& nearFeature,
-    labelList& nearIndex
+    List<pointIndexHit>& nearInfo
 ) const
 {
     nearFeature.setSize(samples.size());
     nearFeature = -1;
-    nearIndex.setSize(samples.size());
-    nearIndex = -1;
+    nearInfo.setSize(samples.size());
+    nearInfo = pointIndexHit();
 
     forAll(pointTrees_, featI)
     {
@@ -453,14 +640,7 @@ void Foam::refinementFeatures::findNearestPoint
                 scalar distSqr;
                 if (nearFeature[sampleI] != -1)
                 {
-                    label nearFeatI = nearFeature[sampleI];
-                    const indexedOctree<treeDataPoint>& nearTree =
-                        pointTrees_[nearFeatI];
-                    label featPointI =
-                        nearTree.shapes().pointLabels()[nearIndex[sampleI]];
-                    const point& featPt =
-                        operator[](nearFeatI).points()[featPointI];
-                    distSqr = magSqr(featPt-sample);
+                    distSqr = magSqr(nearInfo[sampleI].hitPoint()-sample);
                 }
                 else
                 {
@@ -472,7 +652,12 @@ void Foam::refinementFeatures::findNearestPoint
                 if (info.hit())
                 {
                     nearFeature[sampleI] = featI;
-                    nearIndex[sampleI] = info.index();
+                    nearInfo[sampleI] = pointIndexHit
+                    (
+                        info.hit(),
+                        info.hitPoint(),
+                        tree.shapes().pointLabels()[info.index()]
+                    );
                 }
             }
         }
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.H b/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.H
index aa2807bac1..4cfaf20ffa 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.H
+++ b/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,7 +35,7 @@ SourceFiles
 #ifndef refinementFeatures_H
 #define refinementFeatures_H
 
-#include "featureEdgeMesh.H"
+#include "extendedFeatureEdgeMesh.H"
 #include "indexedOctree.H"
 #include "treeDataEdge.H"
 #include "treeDataPoint.H"
@@ -51,7 +51,7 @@ namespace Foam
 
 class refinementFeatures
 :
-    public PtrList<featureEdgeMesh>
+    public PtrList<extendedFeatureEdgeMesh>
 {
 private:
 
@@ -86,6 +86,20 @@ private:
             labelList& maxLevel
         ) const;
 
+
+protected:
+
+        const PtrList<indexedOctree<treeDataEdge> >& edgeTrees() const
+        {
+            return edgeTrees_;
+        }
+
+        const PtrList<indexedOctree<treeDataPoint> >& pointTrees() const
+        {
+            return pointTrees_;
+        }
+
+
 public:
 
     // Constructors
@@ -97,15 +111,6 @@ public:
             const PtrList<dictionary>& featDicts
         );
 
-        //- Construct from description and do geometric analysis to determine
-        //  feature points
-        refinementFeatures
-        (
-            const objectRegistry& io,
-            const PtrList<dictionary>& featDicts,
-            const scalar minCos
-        );
-
 
     // Member Functions
 
@@ -123,42 +128,38 @@ public:
                 return distances_;
             }
 
-            const PtrList<indexedOctree<treeDataEdge> >& edgeTrees() const
-            {
-                return edgeTrees_;
-            }
-
-            const PtrList<indexedOctree<treeDataPoint> >& pointTrees() const
-            {
-                return pointTrees_;
-            }
-
 
         // Query
 
             //- Highest distance of all features
             scalar maxDistance() const;
 
-            //- Find nearest point on nearest feature edge
+            //- Find nearest point on nearest feature edge. Sets
+            //  - nearFeature: index of feature mesh
+            //  - nearInfo   : location on feature edge and edge index
+            //                 (note: not feature edge index but index into
+            //                  edges() directly)
+            //  - nearNormal : local feature edge normal
             void findNearestEdge
             (
                 const pointField& samples,
                 const scalarField& nearestDistSqr,
                 labelList& nearFeature,
-                List<pointIndexHit>& nearInfo
+                List<pointIndexHit>& nearInfo,
+                vectorField& nearNormal
             ) const;
 
-            //- Find nearest feature point. Is an index into feature points
-            //  which itself is an index into the edgeMesh points.
-            //  So the point index is
-            //      pointTrees()[nearFeature].shapes().pointLabels()[nearIndex]
-            //  Wip.
+            //- Find nearest feature point. Sets
+            //  - nearFeature: index of feature mesh
+            //  - nearInfo   : location on feature point and point index.
+            //                 (note: not index into shapes().pointLabels() but
+            //                  index into points() directly)
             void findNearestPoint
             (
                 const pointField& samples,
                 const scalarField& nearestDistSqr,
                 labelList& nearFeature,
-                labelList& nearIndex
+                List<pointIndexHit>& nearInfo
             ) const;
 
             //- Find shell level higher than ptLevel