diff --git a/applications/utilities/mesh/generation/CV3DMesher/backupPreviousIdeas/featureConformationBackups/23_7_08_edgeAimingReconstuction/insertFeaturePoints.C b/applications/utilities/mesh/generation/CV3DMesher/backupPreviousIdeas/featureConformationBackups/23_7_08_edgeAimingReconstuction/insertFeaturePoints.C
new file mode 100644
index 0000000000..18f3db880a
--- /dev/null
+++ b/applications/utilities/mesh/generation/CV3DMesher/backupPreviousIdeas/featureConformationBackups/23_7_08_edgeAimingReconstuction/insertFeaturePoints.C
@@ -0,0 +1,178 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 "CV3D.H"
+#include "plane.H"
+#include "triSurfaceTools.H"
+#include "mathematicalConstants.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::CV3D::insertFeaturePoints()
+{
+    const edgeList& edges = qSurf_.edges();
+    const pointField& localPts = qSurf_.localPoints();
+
+    labelList featPoints(0);
+    labelListList featPointFeatEdges(0);
+
+    qSurf_.extractFeatures3D
+    (
+        controls_.featAngle,
+        featPoints,
+        featPointFeatEdges
+    );
+
+    forAll(featPoints, i)
+    {
+        label ptI = featPoints[i];
+        const point& featPt = localPts[ptI];
+
+        const labelList& featEdges = featPointFeatEdges[i];
+
+        forAll(featEdges, fE)
+        {
+            label edgeI = featEdges[fE];
+            const edge& featEdge = edges[edgeI];
+
+            // Check direction of edge, if the feature point is at the end()
+            // the reverse direction.
+
+            scalar edgeDirection = 1.0;
+
+            if (ptI == featEdge.end())
+            {
+                edgeDirection = -1.0;
+            }
+
+            point edgeLocalFeatPt = featPt
+              + 2.0*tols_.ppDist*edgeDirection
+              * featEdge.vec(localPts)/featEdge.mag(localPts);
+
+            // Pick up the two faces adjacent to the feature edge
+            const labelList& eFaces = qSurf_.edgeFaces()[edgeI];
+
+            label faceA = eFaces[0];
+            vector nA = qSurf_.faceNormals()[faceA];
+
+            label faceB = eFaces[1];
+            vector nB = qSurf_.faceNormals()[faceB];
+
+            // Intersect planes parallel to faceA and faceB offset by ppDist
+            // and the plane defined by edgeLocalFeatPt and the edge vector.
+            plane planeA(edgeLocalFeatPt - tols_.ppDist*nA, nA);
+            plane planeB(edgeLocalFeatPt - tols_.ppDist*nB, nB);
+
+            plane planeF(edgeLocalFeatPt, featEdge.vec(localPts));
+
+            point refPt = planeF.planePlaneIntersect(planeA,planeB);
+
+            point faceAVert =
+                localPts[triSurfaceTools::oppositeVertex(qSurf_, faceA, edgeI)];
+
+            // Determine convex or concave angle
+            if (((faceAVert - edgeLocalFeatPt) & nB) < 0)
+            {
+                // Convex. So refPt will be inside domain and hence a master point
+
+                // Insert the master point refering the the first slave
+                label masterPtIndex = insertPoint(refPt, number_of_vertices() + 1);
+
+                // Insert the slave points by reflecting refPt in both faces.
+                // with each slave refering to the master
+
+                point reflectedA = refPt + 2*((edgeLocalFeatPt - refPt) & nA)*nA;
+                insertPoint(reflectedA, masterPtIndex);
+
+                point reflectedB = refPt + 2*((edgeLocalFeatPt - refPt) & nB)*nB;
+                insertPoint(reflectedB, masterPtIndex);
+            }
+            else
+            {
+                // Concave. master and reflected points inside the domain.
+                // Generate reflected master to be outside.
+                point reflMasterPt = refPt + 2*(edgeLocalFeatPt - refPt);
+
+                // Reflect refPt in both faces.
+                point reflectedA =
+                    reflMasterPt + 2*((edgeLocalFeatPt - reflMasterPt) & nA)*nA;
+
+                point reflectedB =
+                    reflMasterPt + 2*((edgeLocalFeatPt - reflMasterPt) & nB)*nB;
+
+                scalar totalAngle =
+                180*(mathematicalConstant::pi + acos(mag(nA & nB)))
+                /mathematicalConstant::pi;
+
+                // Number of quadrants the angle should be split into
+                int nQuads = int(totalAngle/controls_.maxQuadAngle) + 1;
+
+                // The number of additional master points needed to obtain the
+                // required number of quadrants.
+                int nAddPoints = min(max(nQuads - 2, 0), 2);
+
+                // index of reflMaster
+                label reflectedMaster = number_of_vertices() + 2 + nAddPoints;
+
+                // Master A is inside.
+                label reflectedAI = insertPoint(reflectedA, reflectedMaster);
+
+                // Master B is inside.
+                insertPoint(reflectedB, reflectedMaster);
+
+                if (nAddPoints == 1)
+                {
+                    // One additinal point is the reflection of the slave point,
+                    // i.e. the original reference point
+                    insertPoint(refPt, reflectedMaster);
+                }
+                else if (nAddPoints == 2)
+                {
+                    point reflectedAa = refPt
+                      - ((edgeLocalFeatPt - reflMasterPt) & nB)*nB;
+                    insertPoint(reflectedAa, reflectedMaster);
+
+                    point reflectedBb = refPt
+                      - ((edgeLocalFeatPt - reflMasterPt) & nA)*nA;
+                    insertPoint(reflectedBb, reflectedMaster);
+                }
+
+                // Slave is outside.
+                insertPoint(reflMasterPt, reflectedAI);
+            }
+        }
+    }
+
+    if (controls_.writeFeatureTriangulation)
+    {
+        writePoints("feat_allPoints.obj", false);
+        writePoints("feat_points.obj", true);
+        writeTriangles("feat_triangles.obj", true);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/CV3DMesher/backupPreviousIdeas/featureConformationBackups/23_7_08_edgeAimingReconstuction/querySurface.C b/applications/utilities/mesh/generation/CV3DMesher/backupPreviousIdeas/featureConformationBackups/23_7_08_edgeAimingReconstuction/querySurface.C
new file mode 100644
index 0000000000..311c6350c1
--- /dev/null
+++ b/applications/utilities/mesh/generation/CV3DMesher/backupPreviousIdeas/featureConformationBackups/23_7_08_edgeAimingReconstuction/querySurface.C
@@ -0,0 +1,296 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 "querySurface.H"
+#include "mathematicalConstants.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::querySurface::querySurface(const fileName& surfaceFileName)
+:
+    triSurface(surfaceFileName),
+    rndGen_(12345),
+    bb_(localPoints()),
+    tree_
+    (
+        treeDataTriSurface(*this),
+        bb_.extend(rndGen_, 1e-3), // slightly randomize bb
+        8,                         // maxLevel
+        4, //10,                   // leafsize
+        10.0 //3.0                 // duplicity
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::querySurface::~querySurface()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::labelList Foam::querySurface::extractFeatures2D
+(
+    const scalar featAngle
+) const
+{
+    scalar featCos = cos(mathematicalConstant::pi*featAngle/180.0);
+
+    const labelListList& edgeFaces = this->edgeFaces();
+    const pointField& localPoints = this->localPoints();
+    const edgeList& edges = this->edges();
+    const vectorField& faceNormals = this->faceNormals();
+
+    DynamicList<label> featEdges(edges.size());
+
+    forAll(edgeFaces, edgeI)
+    {
+        const edge& e = edges[edgeI];
+
+        if (magSqr(e.vec(localPoints) & vector(1,1,0)) < SMALL)
+        {
+            const labelList& eFaces = edgeFaces[edgeI];
+
+            if
+            (
+                eFaces.size() == 2
+             && (faceNormals[eFaces[0]] & faceNormals[eFaces[1]]) < featCos
+            )
+            {
+                featEdges.append(edgeI);
+            }
+        }
+    }
+
+    return featEdges.shrink();
+}
+
+void Foam::querySurface::extractFeatures3D
+(
+    const scalar featAngle,
+    labelList& featPoints,
+    labelListList& featPointFeatEdges
+) const
+{
+    scalar featCos = cos(mathematicalConstant::pi*featAngle/180.0);
+
+    const labelListList& edgeFaces = this->edgeFaces();
+    const labelListList& pointEdges = this->pointEdges();
+    const pointField& localPoints = this->localPoints();
+    // const edgeList& edges = this->edges();
+    const vectorField& faceNormals = this->faceNormals();
+
+    DynamicList<label> tempFeatPoints(localPoints.size());
+
+    DynamicList<DynamicList<label> > tempFeatPointFeatEdges(localPoints.size());
+
+    forAll(pointEdges, ptI)
+    {
+        // const point& p = localPoints[ptI];
+
+        // const vector& pNormal = pointNormals[ptI];
+
+        const labelList& pEdges = pointEdges[ptI];
+
+        DynamicList<label> tempFeatEdges(0);
+
+        forAll(pEdges, pE)
+        {
+            label edgeI = pEdges[pE];
+
+            const labelList& eFaces = edgeFaces[edgeI];
+
+            if
+            (
+                eFaces.size() == 2
+                && (faceNormals[eFaces[0]] & faceNormals[eFaces[1]]) < featCos
+            )
+            {
+                tempFeatEdges.append(edgeI);
+            }
+        }
+
+        tempFeatEdges.shrink();
+
+        if(tempFeatEdges.size())
+        {
+            tempFeatPoints.append(ptI);
+
+            tempFeatPointFeatEdges.append(tempFeatEdges);
+        }
+    }
+
+    featPoints.transfer(tempFeatPoints.shrink());
+
+    tempFeatPointFeatEdges.shrink();
+
+    featPointFeatEdges.setSize(tempFeatPointFeatEdges.size());
+
+    forAll(featPointFeatEdges, fPFE)
+    {
+        featPointFeatEdges[fPFE].transfer(tempFeatPointFeatEdges[fPFE].shrink());
+    }
+}
+
+Foam::indexedOctree<Foam::treeDataTriSurface>::volumeType
+Foam::querySurface::insideOutside
+(
+    const scalar searchSpan2,
+    const point& pt
+) const
+{
+    if (!bb_.contains(pt))
+    {
+        return indexedOctree<treeDataTriSurface>::OUTSIDE;
+    }
+
+    pointIndexHit pHit = tree_.findNearest(pt, searchSpan2);
+
+    if (!pHit.hit())
+    {
+        return tree_.getVolumeType(pt);
+    }
+    else
+    {
+        return indexedOctree<treeDataTriSurface>::MIXED;
+    }
+}
+
+
+// Check if point is inside surface
+bool Foam::querySurface::inside(const point& pt) const
+{
+    if (!bb_.contains(pt))
+    {
+        return false;
+    }
+
+    return
+    (
+        tree_.getVolumeType(pt) == indexedOctree<treeDataTriSurface>::INSIDE
+    );
+}
+
+
+// Check if point is outside surface
+bool Foam::querySurface::outside(const point& pt) const
+{
+    if (!bb_.contains(pt))
+    {
+        return true;
+    }
+
+    return
+    (
+        tree_.getVolumeType(pt) == indexedOctree<treeDataTriSurface>::OUTSIDE
+    );
+}
+
+
+// Check if point is inside surface by at least dist2
+bool Foam::querySurface::wellInside(const point& pt, const scalar dist2) const
+{
+    if (!bb_.contains(pt))
+    {
+        return false;
+    }
+
+    pointIndexHit pHit = tree_.findNearest(pt, dist2);
+
+    if (pHit.hit())
+    {
+        return false;
+    }
+    else
+    {
+        return
+            tree_.getVolumeType(pt)
+         == indexedOctree<treeDataTriSurface>::INSIDE;
+    }
+}
+
+
+// Check if point is outside surface by at least dist2
+bool Foam::querySurface::wellOutside(const point& pt, const scalar dist2) const
+{
+    if (!bb_.contains(pt))
+    {
+        return true;
+    }
+
+    pointIndexHit pHit = tree_.findNearest(pt, dist2);
+
+    if (pHit.hit())
+    {
+        return false;
+    }
+    else
+    {
+        return
+            tree_.getVolumeType(pt)
+         == indexedOctree<treeDataTriSurface>::OUTSIDE;
+    }
+}
+
+
+void Foam::querySurface::writeTreeOBJ() const
+{
+    OFstream str("tree.obj");
+    label vertI = 0;
+
+    const List<indexedOctree<treeDataTriSurface>::node>& nodes = tree_.nodes();
+
+    forAll(nodes, nodeI)
+    {
+        const indexedOctree<treeDataTriSurface>::node& nod = nodes[nodeI];
+
+        const treeBoundBox& bb = nod.bb_;
+
+        const pointField points(bb.points());
+
+        label startVertI = vertI;
+
+        forAll(points, i)
+        {
+            meshTools::writeOBJ(str, points[i]);
+            vertI++;
+        }
+
+        const edgeList edges(treeBoundBox::edges);
+
+        forAll(edges, i)
+        {
+            const edge& e = edges[i];
+
+            str << "l " << e[0]+startVertI+1 << ' ' << e[1]+startVertI+1
+                << nl;
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/CV3DMesher/backupPreviousIdeas/hardCodedSimpleCubeForPolyTopoChange.H b/applications/utilities/mesh/generation/CV3DMesher/backupPreviousIdeas/hardCodedSimpleCubeForPolyTopoChange.H
new file mode 100644
index 0000000000..8fd0e855ff
--- /dev/null
+++ b/applications/utilities/mesh/generation/CV3DMesher/backupPreviousIdeas/hardCodedSimpleCubeForPolyTopoChange.H
@@ -0,0 +1,96 @@
+    IOobject io
+    (
+        Foam::polyMesh::defaultRegion,
+        runTime.constant(),
+        runTime,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    );
+
+    polyMesh pMesh(io);
+
+    IOobject io2
+    (
+        Foam::polyMesh::defaultRegion,
+        runTime.timeName(),
+        runTime,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    );
+
+    polyMesh pMesh2(io2, pointField(0), faceList(0), cellList(0));
+
+    polyTopoChange meshCreator(pMesh.boundaryMesh().size());
+
+    meshCreator.addPoint(point(-1,-1,-1), -1, -1, true);
+    meshCreator.addPoint(point(1,-1,-1), -1, -1, true);
+    meshCreator.addPoint(point(1,1,-1), -1, -1, true);
+    meshCreator.addPoint(point(-1,1,-1), -1, -1, true);
+    meshCreator.addPoint(point(-1,-1,1), -1, -1, true);
+    meshCreator.addPoint(point(1,-1,1), -1, -1, true);
+    meshCreator.addPoint(point(1,1,1), -1, -1, true);
+    meshCreator.addPoint(point(-1,1,1), -1, -1, true);
+    
+    meshCreator.addCell(-1, -1, -1, -1, -1);
+
+    labelList faceConstruct(4);
+
+    faceConstruct[0] = 1;
+    faceConstruct[1] = 2;
+    faceConstruct[2] = 6;
+    faceConstruct[3] = 5;
+
+    meshCreator.addFace(face(faceConstruct), 0, -1, -1, -1, -1, false, 0, -1, false);
+
+    faceConstruct[0] = 0;
+    faceConstruct[1] = 4;
+    faceConstruct[2] = 7;
+    faceConstruct[3] = 3;
+
+    meshCreator.addFace(face(faceConstruct), 0, -1, -1, -1, -1, false, 0, -1, false);
+
+    faceConstruct[0] = 2;
+    faceConstruct[1] = 3;
+    faceConstruct[2] = 7;
+    faceConstruct[3] = 6;
+
+    meshCreator.addFace(face(faceConstruct), 0, -1, -1, -1, -1, false, 0, -1, false);
+
+    faceConstruct[0] = 0;
+    faceConstruct[1] = 1;
+    faceConstruct[2] = 5;
+    faceConstruct[3] = 4;
+
+    meshCreator.addFace(face(faceConstruct), 0, -1, -1, -1, -1, false, 0, -1, false);
+
+    faceConstruct[0] = 0;
+    faceConstruct[1] = 3;
+    faceConstruct[2] = 2;
+    faceConstruct[3] = 1;
+
+    meshCreator.addFace(face(faceConstruct), 0, -1, -1, -1, -1, false, 0, -1, false);
+
+    faceConstruct[0] = 4;
+    faceConstruct[1] = 5;
+    faceConstruct[2] = 6;
+    faceConstruct[3] = 7;
+
+    meshCreator.addFace(face(faceConstruct), 0, -1, -1, -1, -1, false, 0, -1, false);
+
+    Info<< meshCreator.points() << endl;
+    Info<< meshCreator.faces() << endl;
+    Info<< meshCreator.faceOwner() << endl;
+    Info<< meshCreator.faceNeighbour() << endl;
+
+//     calcDualMesh(meshCreator);
+
+    autoPtr<mapPolyMesh> map = meshCreator.changeMesh(pMesh, false);
+
+    pMesh.updateMesh(map);
+
+    if (!pMesh.write())
+    {
+        FatalErrorIn("CV3D::writeMesh(const Time& runTime)")
+            << "Failed writing polyMesh."
+            << exit(FatalError);
+    }
\ No newline at end of file
diff --git a/applications/utilities/mesh/generation/CV3DMesher/calcDualMesh.C b/applications/utilities/mesh/generation/CV3DMesher/calcDualMesh.C
index 74c66438f0..018abb495d 100644
--- a/applications/utilities/mesh/generation/CV3DMesher/calcDualMesh.C
+++ b/applications/utilities/mesh/generation/CV3DMesher/calcDualMesh.C
@@ -391,8 +391,6 @@ void Foam::CV3D::calcDualMesh
 
         patchOwners[p].shrink();
 
-        // TODO SORT BOUNDARY FACES AND NEIGHBOURS?
-
         patchSizes[p] = patchFaces[p].size();
 
         patchStarts[p] = nInternalFaces + nBoundaryFaces;
diff --git a/applications/utilities/mesh/generation/CV3DMesher/insertFeaturePoints.C b/applications/utilities/mesh/generation/CV3DMesher/insertFeaturePoints.C
index 18f3db880a..af97d348cf 100644
--- a/applications/utilities/mesh/generation/CV3DMesher/insertFeaturePoints.C
+++ b/applications/utilities/mesh/generation/CV3DMesher/insertFeaturePoints.C
@@ -39,13 +39,17 @@ void Foam::CV3D::insertFeaturePoints()
     labelList featPoints(0);
     labelListList featPointFeatEdges(0);
 
-    qSurf_.extractFeatures3D
+    qSurf_.extractFeatures
     (
         controls_.featAngle,
         featPoints,
         featPointFeatEdges
     );
 
+    scalar planeErrorAngle = 0.1*controls_.featAngle;
+
+    scalar planeErrorAngleCos = cos(mathematicalConstant::pi*planeErrorAngle/180.0);
+
     forAll(featPoints, i)
     {
         label ptI = featPoints[i];
@@ -53,25 +57,17 @@ void Foam::CV3D::insertFeaturePoints()
 
         const labelList& featEdges = featPointFeatEdges[i];
 
+        DynamicList<label> convexEdges(0);
+        DynamicList<label> concaveEdges(0);
+
+        List<Pair<vector> > planeNormalPairs(featEdges.size());
+
+        // Classify the edges around the feature point.
         forAll(featEdges, fE)
         {
             label edgeI = featEdges[fE];
             const edge& featEdge = edges[edgeI];
 
-            // Check direction of edge, if the feature point is at the end()
-            // the reverse direction.
-
-            scalar edgeDirection = 1.0;
-
-            if (ptI == featEdge.end())
-            {
-                edgeDirection = -1.0;
-            }
-
-            point edgeLocalFeatPt = featPt
-              + 2.0*tols_.ppDist*edgeDirection
-              * featEdge.vec(localPts)/featEdge.mag(localPts);
-
             // Pick up the two faces adjacent to the feature edge
             const labelList& eFaces = qSurf_.edgeFaces()[edgeI];
 
@@ -81,91 +77,108 @@ void Foam::CV3D::insertFeaturePoints()
             label faceB = eFaces[1];
             vector nB = qSurf_.faceNormals()[faceB];
 
-            // Intersect planes parallel to faceA and faceB offset by ppDist
-            // and the plane defined by edgeLocalFeatPt and the edge vector.
-            plane planeA(edgeLocalFeatPt - tols_.ppDist*nA, nA);
-            plane planeB(edgeLocalFeatPt - tols_.ppDist*nB, nB);
-
-            plane planeF(edgeLocalFeatPt, featEdge.vec(localPts));
-
-            point refPt = planeF.planePlaneIntersect(planeA,planeB);
-
             point faceAVert =
                 localPts[triSurfaceTools::oppositeVertex(qSurf_, faceA, edgeI)];
 
             // Determine convex or concave angle
-            if (((faceAVert - edgeLocalFeatPt) & nB) < 0)
+            if (((faceAVert - featPt) & nB) < 0)
             {
-                // Convex. So refPt will be inside domain and hence a master point
-
-                // Insert the master point refering the the first slave
-                label masterPtIndex = insertPoint(refPt, number_of_vertices() + 1);
-
-                // Insert the slave points by reflecting refPt in both faces.
-                // with each slave refering to the master
-
-                point reflectedA = refPt + 2*((edgeLocalFeatPt - refPt) & nA)*nA;
-                insertPoint(reflectedA, masterPtIndex);
-
-                point reflectedB = refPt + 2*((edgeLocalFeatPt - refPt) & nB)*nB;
-                insertPoint(reflectedB, masterPtIndex);
+                // Convex feature edge
+                convexEdges.append(edgeI);
             }
             else
             {
-                // Concave. master and reflected points inside the domain.
-                // Generate reflected master to be outside.
-                point reflMasterPt = refPt + 2*(edgeLocalFeatPt - refPt);
-
-                // Reflect refPt in both faces.
-                point reflectedA =
-                    reflMasterPt + 2*((edgeLocalFeatPt - reflMasterPt) & nA)*nA;
-
-                point reflectedB =
-                    reflMasterPt + 2*((edgeLocalFeatPt - reflMasterPt) & nB)*nB;
-
-                scalar totalAngle =
-                180*(mathematicalConstant::pi + acos(mag(nA & nB)))
-                /mathematicalConstant::pi;
-
-                // Number of quadrants the angle should be split into
-                int nQuads = int(totalAngle/controls_.maxQuadAngle) + 1;
-
-                // The number of additional master points needed to obtain the
-                // required number of quadrants.
-                int nAddPoints = min(max(nQuads - 2, 0), 2);
-
-                // index of reflMaster
-                label reflectedMaster = number_of_vertices() + 2 + nAddPoints;
-
-                // Master A is inside.
-                label reflectedAI = insertPoint(reflectedA, reflectedMaster);
-
-                // Master B is inside.
-                insertPoint(reflectedB, reflectedMaster);
-
-                if (nAddPoints == 1)
-                {
-                    // One additinal point is the reflection of the slave point,
-                    // i.e. the original reference point
-                    insertPoint(refPt, reflectedMaster);
-                }
-                else if (nAddPoints == 2)
-                {
-                    point reflectedAa = refPt
-                      - ((edgeLocalFeatPt - reflMasterPt) & nB)*nB;
-                    insertPoint(reflectedAa, reflectedMaster);
-
-                    point reflectedBb = refPt
-                      - ((edgeLocalFeatPt - reflMasterPt) & nA)*nA;
-                    insertPoint(reflectedBb, reflectedMaster);
-                }
-
-                // Slave is outside.
-                insertPoint(reflMasterPt, reflectedAI);
+                // Concave feature edge
+                concaveEdges.append(edgeI);
             }
+
+            // Identify the normals of the faces attached to feature edges to
+            // identify the unique planes to be reconstructed.  The triangulated
+            // surface will usually not mean that two feature edges that should
+            // bound a plane are attached to the same face.
+
+            planeNormalPairs[fE].first() = nA;
+            planeNormalPairs[fE].second() = nB;
         }
+
+        convexEdges.shrink();
+        concaveEdges.shrink();
+
+        Info<< nl << convexEdges
+            << nl << concaveEdges
+            << nl << planeNormalPairs
+            << endl;
+
+        List<vector> uniquePlaneNormals(featEdges.size());
+        label uniquePlaneNormalI = 0;
+
+        List<Pair<bool> > planeMatchedStatus(featEdges.size(), Pair<bool>(false,false));
+
+        Info<< planeMatchedStatus << endl;
+
+        // Examine the plane normals to identify unique planes.
+        forAll(planeNormalPairs, nA)
+        {
+            const Pair<vector>& normalPairA = planeNormalPairs[nA];
+
+            if (!planeMatchedStatus[nA].first())
+            {
+                const vector& nAf = normalPairA.first();
+
+                scalar minNormalDotProduct = 1 + SMALL;
+
+                forAll(planeNormalPairs, nB)
+                {
+                    if (nA == nB)
+                    {
+                        continue;
+                    }
+
+                    const Pair<vector>& normalPairB = planeNormalPairs[nB];
+
+                    if (!planeMatchedStatus[nB].first())
+                    {
+                        const vector& nBf = normalPairB.first();
+
+                        scalar normalDotProduct = nAf & nBf;
+
+                        if (normalDotProduct < minNormalDotProduct)
+                        {
+                            minNormalDotProduct = normalDotProduct;
+                        }
+                    }
+
+                    if (!planeMatchedStatus[nB].second())
+                    {
+                        const vector& nBs = normalPairA.second();
+                    }
+                }
+            }
+
+            if (!planeMatchedStatus[nA].second())
+            {
+                const vector& nAs = normalPairA.second();
+            }
+
+            // if (minNormalDotProduct > planeErrorAngleCos)
+            // {
+            //     FatalErrorIn("insertFeaturePoints")
+            //         << "Could not find unique planes matching feature edges "
+            //         << "at point located at "
+            //         << featPt << nl
+            //         << exit(FatalError);
+            // }
+
+            // const vector& matchingPNB = planeNormalsB[matchingPB];
+
+            // uniquePlaneNormals[pA] =
+            // (pNA + matchingPNB)/(mag(pNA + matchingPNB) + VSMALL);
+        }
+
+        Info<< uniquePlaneNormals << endl;
     }
 
+
     if (controls_.writeFeatureTriangulation)
     {
         writePoints("feat_allPoints.obj", false);
diff --git a/applications/utilities/mesh/generation/CV3DMesher/querySurface.C b/applications/utilities/mesh/generation/CV3DMesher/querySurface.C
index 311c6350c1..4575abc1cb 100644
--- a/applications/utilities/mesh/generation/CV3DMesher/querySurface.C
+++ b/applications/utilities/mesh/generation/CV3DMesher/querySurface.C
@@ -53,43 +53,7 @@ Foam::querySurface::~querySurface()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::labelList Foam::querySurface::extractFeatures2D
-(
-    const scalar featAngle
-) const
-{
-    scalar featCos = cos(mathematicalConstant::pi*featAngle/180.0);
-
-    const labelListList& edgeFaces = this->edgeFaces();
-    const pointField& localPoints = this->localPoints();
-    const edgeList& edges = this->edges();
-    const vectorField& faceNormals = this->faceNormals();
-
-    DynamicList<label> featEdges(edges.size());
-
-    forAll(edgeFaces, edgeI)
-    {
-        const edge& e = edges[edgeI];
-
-        if (magSqr(e.vec(localPoints) & vector(1,1,0)) < SMALL)
-        {
-            const labelList& eFaces = edgeFaces[edgeI];
-
-            if
-            (
-                eFaces.size() == 2
-             && (faceNormals[eFaces[0]] & faceNormals[eFaces[1]]) < featCos
-            )
-            {
-                featEdges.append(edgeI);
-            }
-        }
-    }
-
-    return featEdges.shrink();
-}
-
-void Foam::querySurface::extractFeatures3D
+void Foam::querySurface::extractFeatures
 (
     const scalar featAngle,
     labelList& featPoints,
@@ -136,7 +100,7 @@ void Foam::querySurface::extractFeatures3D
 
         tempFeatEdges.shrink();
 
-        if(tempFeatEdges.size())
+        if(tempFeatEdges.size() > 2)
         {
             tempFeatPoints.append(ptI);
 
diff --git a/applications/utilities/mesh/generation/CV3DMesher/querySurface.H b/applications/utilities/mesh/generation/CV3DMesher/querySurface.H
index 8e4a945c40..44e89380bf 100644
--- a/applications/utilities/mesh/generation/CV3DMesher/querySurface.H
+++ b/applications/utilities/mesh/generation/CV3DMesher/querySurface.H
@@ -104,13 +104,9 @@ public:
 
         // Query
 
-            //- Extract feature edges/points(2D)
+            //- Extract feature edges/points
             //  using the given feature angle in deg
-            labelList extractFeatures2D(const scalar featAngle) const;
-
-            //- Extract feature edges/points(3D)
-            //  using the given feature angle in deg
-            void extractFeatures3D
+            void extractFeatures
             (
                 const scalar featAngle,
                 labelList& featPoints,
diff --git a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/G b/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/G
deleted file mode 100644
index cb9783649f..0000000000
--- a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/G
+++ /dev/null
@@ -1,56 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.5                                   |
-|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    object      G;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [1 0 -3 0 0 0 0];
-
-internalField   uniform 0;
-
-boundaryField
-{
-    floor
-    {
-        type            MarshakRadiation;
-        T               T;
-        emissivity      1;
-        value           uniform 0;
-    }
-
-    fixedWalls
-    {
-        type            MarshakRadiation;
-        T               T;
-        emissivity      1;
-        value           uniform 0;
-    }
-
-    ceiling
-    {
-        type            MarshakRadiation;
-        T               T;
-        emissivity      1;
-        value           uniform 0;
-    }
-
-    box
-    {
-        type            MarshakRadiation;
-        T               T;
-        emissivity      1;
-        value           uniform 0;
-    }
-}
-
-// ************************************************************************* //
diff --git a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/T b/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/T
deleted file mode 100644
index 19ee7d9d08..0000000000
--- a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/T
+++ /dev/null
@@ -1,47 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.5                                   |
-|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    object      T;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 0 0 1 0 0 0];
-
-internalField   uniform 300;
-
-boundaryField
-{
-    floor
-    {
-        type            fixedValue;
-        value           uniform 300.0;
-    }
-
-    ceiling
-    {
-        type            fixedValue;
-        value           uniform 300.0;
-    }
-
-    fixedWalls
-    {
-        type            zeroGradient;
-    }
-
-    box
-    {
-        type            fixedValue;
-        value           uniform 500.0;
-    }
-}
-
-// ************************************************************************* //
diff --git a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/U b/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/U
deleted file mode 100644
index 9a0eaf66b5..0000000000
--- a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/U
+++ /dev/null
@@ -1,48 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.5                                   |
-|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volVectorField;
-    object      U;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 1 -1 0 0 0 0];
-
-internalField   uniform (0 0 0);
-
-boundaryField
-{
-    floor           
-    {
-        type            fixedValue;
-        value           uniform (0 0 0);
-    }
-
-    ceiling         
-    {
-        type            fixedValue;
-        value           uniform (0 0 0);
-    }
-
-    fixedWalls      
-    {
-        type            fixedValue;
-        value           uniform (0 0 0);
-    }
-
-    box
-    {
-        type            fixedValue;
-        value           uniform (0 0 0);
-    }
-}
-
-// ************************************************************************* //
diff --git a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/epsilon b/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/epsilon
deleted file mode 100644
index 5f8de5792e..0000000000
--- a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/epsilon
+++ /dev/null
@@ -1,44 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.5                                   |
-|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    object      epsilon;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 2 -3 0 0 0 0];
-
-internalField   uniform 0.01;
-
-boundaryField
-{
-    floor           
-    {
-        type            zeroGradient;
-    }
-
-    ceiling         
-    {
-        type            zeroGradient;
-    }
-
-    fixedWalls      
-    {
-        type            zeroGradient;
-    }
-
-    box
-    {
-        type            zeroGradient;
-    }
-}
-
-// ************************************************************************* //
diff --git a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/k b/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/k
deleted file mode 100644
index a7d2d3f4f0..0000000000
--- a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/k
+++ /dev/null
@@ -1,44 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.5                                   |
-|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    object      k;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 2 -2 0 0 0 0];
-
-internalField   uniform 0.1;
-
-boundaryField
-{
-    floor           
-    {
-        type            zeroGradient;
-    }
-
-    ceiling         
-    {
-        type            zeroGradient;
-    }
-
-    fixedWalls      
-    {
-        type            zeroGradient;
-    }
-
-    box
-    {
-        type            zeroGradient;
-    }
-}
-
-// ************************************************************************* //
diff --git a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/p b/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/p
deleted file mode 100644
index f200e3eeb0..0000000000
--- a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/p
+++ /dev/null
@@ -1,48 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.5                                   |
-|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    object      p;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [1 -1 -2 0 0 0 0];
-
-internalField   uniform 100000;
-
-boundaryField
-{
-    floor
-    {
-        type            calculated;
-        value           uniform 100000;
-    }
-
-    ceiling
-    {
-        type            calculated;
-        value           uniform 100000;
-    }
-
-    fixedWalls
-    {
-        type            calculated;
-        value           uniform 100000;
-    }
-
-    box
-    {
-        type            calculated;
-        value           uniform 100000;
-    }
-}
-
-// ************************************************************************* //
diff --git a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/pd b/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/pd
deleted file mode 100644
index 1841d7882f..0000000000
--- a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/0/pd
+++ /dev/null
@@ -1,48 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.5                                   |
-|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    object      pd;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [1 -1 -2 0 0 0 0];
-
-internalField   uniform 0;
-
-boundaryField
-{
-    floor
-    {
-        type            fixedFluxBuoyantPressure;
-        value           uniform 0;
-    }
-
-    ceiling
-    {
-        type            fixedFluxBuoyantPressure;
-        value           uniform 0;
-    }
-
-    fixedWalls
-    {
-        type            fixedFluxBuoyantPressure;
-        value           uniform 0;
-    }
-
-    box
-    {
-        type            fixedFluxBuoyantPressure;
-        value           uniform 0;
-    }
-}
-
-// ************************************************************************* //
diff --git a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/constant/polyMesh/boundary b/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/constant/polyMesh/boundary
index add6399d3f..aebfe3f9a7 100644
--- a/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/constant/polyMesh/boundary
+++ b/tutorials/buoyantSimpleRadiationFoam/hotRadiationRoom/constant/polyMesh/boundary
@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.5                                   |
+|  \\    /   O peration     | Version:  dev                                   |
 |   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
@@ -10,6 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       polyBoundaryMesh;
+    location    "constant/polyMesh";
     object      boundary;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //