diff --git a/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C b/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C
index 56533c0b00..0cdd440d1a 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C
+++ b/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C
@@ -26,21 +26,14 @@ License
 #include "refinementFeatures.H"
 #include "Time.H"
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-Foam::refinementFeatures::refinementFeatures
+void Foam::refinementFeatures::read
 (
     const objectRegistry& io,
     const PtrList<dictionary>& featDicts
 )
-:
-    PtrList<featureEdgeMesh>(featDicts.size()),
-    levels_(featDicts.size()),
-    edgeTrees_(featDicts.size()),
-    pointTrees_(featDicts.size())
 {
-    // Read features
-
     forAll(featDicts, i)
     {
         const dictionary& dict = featDicts[i];
@@ -75,50 +68,87 @@ Foam::refinementFeatures::refinementFeatures
             << " (" << eMesh.points().size() << " points, "
             << eMesh.edges().size() << " edges)." << endl;
     }
+}
 
 
+void Foam::refinementFeatures::buildTrees
+(
+    const label featI,
+    const labelList& featurePoints
+)
+{
+    const featureEdgeMesh& eMesh = operator[](featI);
+    const pointField& points = eMesh.points();
+    const edgeList& edges = eMesh.edges();
+
+    // Calculate bb of all points
+    treeBoundBox bb(points);
+
+    // Random number generator. Bit dodgy since not exactly random ;-)
+    Random rndGen(65431);
+
+    // Slightly extended bb. Slightly off-centred just so on symmetric
+    // geometry there are less face/edge aligned items.
+    bb = bb.extend(rndGen, 1e-4);
+    bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+    bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+
+    edgeTrees_.set
+    (
+        featI,
+        new indexedOctree<treeDataEdge>
+        (
+            treeDataEdge
+            (
+                false,                  // do not cache bb
+                edges,
+                points,
+                identity(edges.size())
+            ),
+            bb,     // overall search domain
+            8,      // maxLevel
+            10,     // leafsize
+            3.0     // duplicity
+        )
+    );
+
+    pointTrees_.set
+    (
+        featI,
+        new indexedOctree<treeDataPoint>
+        (
+            treeDataPoint(points, featurePoints),
+            bb,     // overall search domain
+            8,      // maxLevel
+            10,     // leafsize
+            3.0     // duplicity
+        )
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::refinementFeatures::refinementFeatures
+(
+    const objectRegistry& io,
+    const PtrList<dictionary>& featDicts
+)
+:
+    PtrList<featureEdgeMesh>(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();
-
-        // Calculate bb of all points
-        treeBoundBox bb(points);
-
-        // Random number generator. Bit dodgy since not exactly random ;-)
-        Random rndGen(65431);
-
-        // Slightly extended bb. Slightly off-centred just so on symmetric
-        // geometry there are less face/edge aligned items.
-        bb = bb.extend(rndGen, 1e-4);
-        bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-        bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-
-        edgeTrees_.set
-        (
-            i,
-            new indexedOctree<treeDataEdge>
-            (
-                treeDataEdge
-                (
-                    false,                  // do not cache bb
-                    edges,
-                    points,
-                    identity(edges.size())
-                ),
-                bb,     // overall search domain
-                8,      // maxLevel
-                10,     // leafsize
-                3.0     // duplicity
-            )
-        );
-
-
-        // Detect feature points from edges.
         const labelListList& pointEdges = eMesh.pointEdges();
+
         DynamicList<label> featurePoints;
         forAll(pointEdges, pointI)
         {
@@ -129,21 +159,79 @@ Foam::refinementFeatures::refinementFeatures
         }
 
         Info<< "Detected " << featurePoints.size()
-            << " featurePoints out of " << points.size()
+            << " featurePoints out of " << pointEdges.size()
             << " on feature " << eMesh.name() << endl;
 
-        pointTrees_.set
-        (
-            i,
-            new indexedOctree<treeDataPoint>
-            (
-                treeDataPoint(points, featurePoints),
-                bb,     // overall search domain
-                8,      // maxLevel
-                10,     // leafsize
-                3.0     // duplicity
-            )
-        );
+        buildTrees(i, featurePoints);
+    }
+}
+
+
+Foam::refinementFeatures::refinementFeatures
+(
+    const objectRegistry& io,
+    const PtrList<dictionary>& featDicts,
+    const scalar minCos
+)
+:
+    PtrList<featureEdgeMesh>(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()
+            << " on feature " << eMesh.name()
+            << " when using feature cos " << minCos << endl;
+
+        buildTrees(i, featurePoints);
     }
 }
 
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.H b/src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.H
index 3e56be043c..5442e140dd 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 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -69,17 +69,33 @@ private:
 
     // Private Member Functions
 
+        //- Read set of feature edge meshes
+        void read(const objectRegistry&, const PtrList<dictionary>&);
+
+        //- Build edge tree and feature point tree
+        void buildTrees(const label, const labelList&);
+
+
 public:
 
     // Constructors
 
-        //- Construct from components
+        //- Construct from description
         refinementFeatures
         (
             const objectRegistry& io,
             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
 
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
index 52b2cd1b5b..38afad6fa2 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
+++ b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
@@ -917,6 +917,127 @@ void Foam::refinementSurfaces::findNearestIntersection
 }
 
 
+void Foam::refinementSurfaces::findNearestIntersection
+(
+    const labelList& surfacesToTest,
+    const pointField& start,
+    const pointField& end,
+
+    labelList& surface1,
+    List<pointIndexHit>& hit1,
+    labelList& region1,
+    vectorField& normal1,
+
+    labelList& surface2,
+    List<pointIndexHit>& hit2,
+    labelList& region2,
+    vectorField& normal2
+) const
+{
+    // 1. intersection from start to end
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Initialize arguments
+    surface1.setSize(start.size());
+    surface1 = -1;
+    hit1.setSize(start.size());
+    region1.setSize(start.size());
+    normal1.setSize(start.size());
+
+    // Current end of segment to test.
+    pointField nearest(end);
+    // Work array
+    List<pointIndexHit> nearestInfo(start.size());
+    labelList region;
+    vectorField normal;
+
+    forAll(surfacesToTest, testI)
+    {
+        label surfI = surfacesToTest[testI];
+        const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
+
+        // See if any intersection between start and current nearest
+        geom.findLine(start, nearest, nearestInfo);
+        geom.getRegion(nearestInfo, region);
+        geom.getNormal(nearestInfo, normal);
+
+        forAll(nearestInfo, pointI)
+        {
+            if (nearestInfo[pointI].hit())
+            {
+                hit1[pointI] = nearestInfo[pointI];
+                surface1[pointI] = surfI;
+                region1[pointI] = region[pointI];
+                normal1[pointI] = normal[pointI];
+                nearest[pointI] = hit1[pointI].hitPoint();
+            }
+        }
+    }
+
+
+    // 2. intersection from end to last intersection
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Find the nearest intersection from end to start. Note that we initialize
+    // to the first intersection (if any).
+    surface2 = surface1;
+    hit2 = hit1;
+    region2 = region1;
+    normal2 = normal1;
+
+    // Set current end of segment to test.
+    forAll(nearest, pointI)
+    {
+        if (hit1[pointI].hit())
+        {
+            nearest[pointI] = hit1[pointI].hitPoint();
+        }
+        else
+        {
+            // Disable testing by setting to end.
+            nearest[pointI] = end[pointI];
+        }
+    }
+
+    forAll(surfacesToTest, testI)
+    {
+        label surfI = surfacesToTest[testI];
+        const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
+
+        // See if any intersection between end and current nearest
+        geom.findLine(end, nearest, nearestInfo);
+        geom.getRegion(nearestInfo, region);
+        geom.getNormal(nearestInfo, normal);
+
+        forAll(nearestInfo, pointI)
+        {
+            if (nearestInfo[pointI].hit())
+            {
+                hit2[pointI] = nearestInfo[pointI];
+                surface2[pointI] = surfI;
+                region2[pointI] = region[pointI];
+                normal2[pointI] = normal[pointI];
+                nearest[pointI] = hit2[pointI].hitPoint();
+            }
+        }
+    }
+
+
+    // Make sure that if hit1 has hit something, hit2 will have at least the
+    // same point (due to tolerances it might miss its end point)
+    forAll(hit1, pointI)
+    {
+        if (hit1[pointI].hit() && !hit2[pointI].hit())
+        {
+            hit2[pointI] = hit1[pointI];
+            surface2[pointI] = surface1[pointI];
+            region2[pointI] = region1[pointI];
+            normal2[pointI] = normal1[pointI];
+        }
+    }
+}
+
+
 void Foam::refinementSurfaces::findAnyIntersection
 (
     const pointField& start,
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
index b0edaec614..2f685ed4ff 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
+++ b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
@@ -317,6 +317,24 @@ public:
                 labelList& region2
             ) const;
 
+            //- findNearestIntersection but also get normals
+            void findNearestIntersection
+            (
+                const labelList& surfacesToTest,
+                const pointField& start,
+                const pointField& end,
+
+                labelList& surface1,
+                List<pointIndexHit>& hit1,
+                labelList& region1,
+                vectorField& normal1,
+
+                labelList& surface2,
+                List<pointIndexHit>& hit2,
+                labelList& region2,
+                vectorField& normal2
+            ) const;
+
             //- Used for debugging only: find intersection of edge.
             void findAnyIntersection
             (