diff --git a/applications/utilities/mesh/generation/Allwmake b/applications/utilities/mesh/generation/Allwmake
index ef28f9aaf5..48ca530c96 100755
--- a/applications/utilities/mesh/generation/Allwmake
+++ b/applications/utilities/mesh/generation/Allwmake
@@ -4,12 +4,15 @@ set -x
 
 wmake blockMesh
 wmake all extrude
-wmake extrude2DMesh
+
+extrude2DMesh/Allwmake
+
 wmake snappyHexMesh
 
 if [ -d "$CGAL_ARCH_PATH" ]
 then
-    cd cvMesh && ./Allwmake
+    cvMesh/Allwmake
+    cv2DMesh/Allwmake
 fi
 
 # ----------------------------------------------------------------- end-of-file
diff --git a/applications/utilities/mesh/generation/cv2DMesh/Allwclean b/applications/utilities/mesh/generation/cv2DMesh/Allwclean
new file mode 100755
index 0000000000..d0ae53e415
--- /dev/null
+++ b/applications/utilities/mesh/generation/cv2DMesh/Allwclean
@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+set -x
+
+wclean libso conformalVoronoi2DMesh
+wclean
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/applications/utilities/mesh/generation/cv2DMesh/Allwmake b/applications/utilities/mesh/generation/cv2DMesh/Allwmake
new file mode 100755
index 0000000000..5486849957
--- /dev/null
+++ b/applications/utilities/mesh/generation/cv2DMesh/Allwmake
@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+set -x
+
+wmake libso conformalVoronoi2DMesh
+wmake
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/applications/utilities/mesh/generation/cv2DMesh/CV2D.C b/applications/utilities/mesh/generation/cv2DMesh/CV2D.C
index f6ab73cc3f..c95a61573e 100644
--- a/applications/utilities/mesh/generation/cv2DMesh/CV2D.C
+++ b/applications/utilities/mesh/generation/cv2DMesh/CV2D.C
@@ -30,12 +30,19 @@ License
 #include "uint.H"
 #include "ulong.H"
 
+namespace Foam
+{
+
+defineTypeNameAndDebug(CV2D, 0);
+
+}
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void Foam::CV2D::insertBoundingBox()
 {
     Info<< "insertBoundingBox: creating bounding mesh" << endl;
-    scalar bigSpan = 10*tols_.span;
+    scalar bigSpan = 10*meshControls().span();
     insertPoint(point2D(-bigSpan, -bigSpan), Vb::FAR_POINT);
     insertPoint(point2D(-bigSpan, bigSpan), Vb::FAR_POINT);
     insertPoint(point2D(bigSpan, -bigSpan), Vb::FAR_POINT);
@@ -125,22 +132,26 @@ Foam::CV2D::CV2D
         allGeometry_,
         cvMeshDict.subDict("surfaceConformation")
     ),
-    controls_(cvMeshDict),
+    controls_(cvMeshDict, qSurf_.globalBounds()),
     cellSizeControl_
     (
         allGeometry_,
         cvMeshDict.subDict("motionControl")
     ),
-    tols_(cvMeshDict, controls_.minCellSize, qSurf_.globalBounds()),
     z_
     (
-        (1.0/3.0)
-       *(qSurf_.globalBounds().min().z() + qSurf_.globalBounds().max().z())
+        point
+        (
+            cvMeshDict.subDict("surfaceConformation").lookup("locationInMesh")
+        ).z()
     ),
     startOfInternalPoints_(0),
     startOfSurfacePointPairs_(0),
-    startOfBoundaryConformPointPairs_(0)
+    startOfBoundaryConformPointPairs_(0),
+    featurePoints_()
 {
+    Info<< meshControls() << endl;
+
     insertBoundingBox();
     insertFeaturePoints();
 }
@@ -183,7 +194,7 @@ void Foam::CV2D::insertPoints
 
     Info<< nVert << " vertices inserted" << endl;
 
-    if (controls_.writeInitialTriangulation)
+    if (meshControls().objOutput())
     {
         // Checking validity of triangulation
         assert(is_valid());
@@ -200,7 +211,7 @@ void Foam::CV2D::insertPoints(const fileName& pointFileName)
 
     if (pointsFile.good())
     {
-        insertPoints(point2DField(pointsFile), 0.5*controls_.minCellSize2);
+        insertPoints(point2DField(pointsFile), 0.5*meshControls().minCellSize2());
     }
     else
     {
@@ -220,16 +231,16 @@ void Foam::CV2D::insertGrid()
 
     scalar x0 = qSurf_.globalBounds().min().x();
     scalar xR = qSurf_.globalBounds().max().x() - x0;
-    int ni = int(xR/controls_.minCellSize) + 1;
+    int ni = int(xR/meshControls().minCellSize()) + 1;
     scalar deltax = xR/ni;
 
     scalar y0 = qSurf_.globalBounds().min().y();
     scalar yR = qSurf_.globalBounds().max().y() - y0;
-    int nj = int(yR/controls_.minCellSize) + 1;
+    int nj = int(yR/meshControls().minCellSize()) + 1;
     scalar deltay = yR/nj;
 
     Random rndGen(1321);
-    scalar pert = controls_.randomPurturbation*min(deltax, deltay);
+    scalar pert = meshControls().randomPerturbation()*min(deltax, deltay);
 
     for (int i=0; i<ni; i++)
     {
@@ -237,13 +248,13 @@ void Foam::CV2D::insertGrid()
         {
             point p(x0 + i*deltax, y0 + j*deltay, 0);
 
-            if (controls_.randomiseInitialGrid)
+            if (meshControls().randomiseInitialGrid())
             {
                 p.x() += pert*(rndGen.scalar01() - 0.5);
                 p.y() += pert*(rndGen.scalar01() - 0.5);
             }
 
-            if (qSurf_.wellInside(p, 0.5*controls_.minCellSize2))
+            if (qSurf_.wellInside(p, 0.5*meshControls().minCellSize2()))
             {
                 insert(Point(p.x(), p.y()))->index() = nVert++;
             }
@@ -252,7 +263,7 @@ void Foam::CV2D::insertGrid()
 
     Info<< nVert << " vertices inserted" << endl;
 
-    if (controls_.writeInitialTriangulation)
+    if (meshControls().objOutput())
     {
         // Checking validity of triangulation
         assert(is_valid());
@@ -267,21 +278,16 @@ void Foam::CV2D::insertSurfacePointPairs()
 {
     startOfSurfacePointPairs_ = number_of_vertices();
 
-    if (controls_.insertSurfaceNearestPointPairs)
+    if (meshControls().insertSurfaceNearestPointPairs())
     {
         insertSurfaceNearestPointPairs();
     }
 
-    if (controls_.writeNearestTriangulation)
-    {
-        writeFaces("near_allFaces.obj", false);
-        writeFaces("near_faces.obj", true);
-        writeTriangles("near_triangles.obj", true);
-    }
+    write("nearest");
 
     // Insertion of point-pairs for near-points may cause protrusions
     // so insertBoundaryConformPointPairs must be executed last
-    if (controls_.insertSurfaceNearPointPairs)
+    if (meshControls().insertSurfaceNearPointPairs())
     {
         insertSurfaceNearPointPairs();
     }
@@ -292,7 +298,7 @@ void Foam::CV2D::insertSurfacePointPairs()
 
 void Foam::CV2D::boundaryConform()
 {
-    if (!controls_.insertSurfaceNearestPointPairs)
+    if (!meshControls().insertSurfaceNearestPointPairs())
     {
         markNearBoundaryPoints();
     }
@@ -308,7 +314,7 @@ void Foam::CV2D::boundaryConform()
         fit->faceIndex() = Fb::SAVE_CHANGED;
     }
 
-    for (label iter=1; iter<=controls_.maxBoundaryConformingIter; iter++)
+    for (label iter=1; iter<=meshControls().maxBoundaryConformingIter(); iter++)
     {
         label nIntersections = insertBoundaryConformPointPairs
         (
@@ -347,6 +353,8 @@ void Foam::CV2D::boundaryConform()
     }
 
     Info<< nl;
+
+    write("boundary");
 }
 
 
@@ -412,7 +420,7 @@ void Foam::CV2D::newPoints(const scalar relaxation)
     scalarField sizes
     (
         number_of_vertices(),
-        controls_.minCellSize
+        meshControls().minCellSize()
     );
 
     Field<vector2D> alignments
@@ -439,14 +447,13 @@ void Foam::CV2D::newPoints(const scalar relaxation)
             qSurf_.findSurfaceNearest
             (
                 toPoint3D(vert),
-                tols_.span2,
+                meshControls().span2(),
                 pHit,
                 hitSurface
             );
 
             if (pHit.hit())
             {
-
                 vectorField norm(1);
                 allGeometry_[hitSurface].getNormal
                 (
@@ -456,23 +463,8 @@ void Foam::CV2D::newPoints(const scalar relaxation)
 
                 alignments[vit->index()] = toPoint2D(norm[0]);
 
-                scalar surfDist = mag(toPoint3D(vert) - pHit.hitPoint());
-
-                 /*if (surfDist < 0.2)
-                 {
-                     sizes[vit->index()] *= 0.4;
-                 }*/
-
-                if (surfDist < 0.2)
-                {
-                    sizes[vit->index()] *= (1 - 0.1)*surfDist/0.2 + 0.1;
-                }
+                sizes[vit->index()] = cellSizeControl_.cellSize(toPoint3D(vit->point()));
             }
-
-            // if (vert.x() > 0)
-            // {
-            //     sizes[vit->index()] *= 0.5;
-            // }
         }
     }
 
@@ -486,7 +478,7 @@ void Foam::CV2D::newPoints(const scalar relaxation)
 
     PackedBoolList pointToBeRetained(startOfSurfacePointPairs_, true);
 
-    DynamicList<point2D> pointsToInsert;
+    std::list<Point> pointsToInsert;
 
     for
     (
@@ -621,7 +613,7 @@ void Foam::CV2D::newPoints(const scalar relaxation)
                 )
                 {
                     // Point insertion
-                    pointsToInsert.append(0.5*(dVA + dVB));
+                    pointsToInsert.push_back(toPoint(0.5*(dVA + dVB)));
                 }
                 else if
                 (
@@ -640,7 +632,7 @@ void Foam::CV2D::newPoints(const scalar relaxation)
                      && pointToBeRetained[vB->index()] == true
                     )
                     {
-                        pointsToInsert.append(0.5*(dVA + dVB));
+                        pointsToInsert.push_back(toPoint(0.5*(dVA + dVB)));
                     }
 
                     if (vA->internalPoint())
@@ -675,6 +667,8 @@ void Foam::CV2D::newPoints(const scalar relaxation)
     // Relax the calculated displacement
     displacementAccumulator *= relaxation;
 
+    label numberOfNewPoints = pointsToInsert.size();
+
     for
     (
         Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
@@ -684,62 +678,66 @@ void Foam::CV2D::newPoints(const scalar relaxation)
     {
         if (vit->internalPoint())
         {
-            if (!pointToBeRetained[vit->index()])
+            if (pointToBeRetained[vit->index()])
             {
-                remove(vit);
-            }
-            else
-            {
-                movePoint
+                pointsToInsert.push_front
                 (
-                    vit,
-                    vit->point()
-                  + K::Vector_2
+                    toPoint
                     (
-                        displacementAccumulator[vit->index()].x(),
-                        displacementAccumulator[vit->index()].y()
+                        toPoint2D(vit->point())
+                      + displacementAccumulator[vit->index()]
                     )
                 );
             }
         }
     }
 
-    removeSurfacePointPairs();
+    // Clear the triangulation and reinsert the bounding box and feature points.
+    // This is faster than removing and moving points.
+    this->clear();
 
-    // Re-index internal points
+    insertBoundingBox();
+
+    reinsertFeaturePoints();
+
+    startOfInternalPoints_ = number_of_vertices();
 
     label nVert = startOfInternalPoints_;
 
+    Info<< "Inserting " << numberOfNewPoints << " new points" << endl;
+
+    // Use the range insert as it is faster than individually inserting points.
+    insert(pointsToInsert.begin(), pointsToInsert.end());
+
     for
     (
-        Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
+        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
         vit != finite_vertices_end();
         ++vit
     )
     {
-        if (vit->internalPoint())
+        if
+        (
+            vit->type() == Vb::INTERNAL_POINT
+         && vit->index() == Vb::INTERNAL_POINT
+        )
         {
             vit->index() = nVert++;
         }
     }
 
-    // Insert new points
-
-    Info<< "Inserting " << pointsToInsert.size() << " new points" << endl;
-
-    forAll(pointsToInsert, i)
-    {
-        insertPoint(pointsToInsert[i], Vb::INTERNAL_POINT);
-    }
-
     Info<< "    Total displacement = " << totalDisp << nl
         << "    Total distance = " << totalDist << nl
         << "    Points added = " << pointsToInsert.size()
         << endl;
 
+    write("internal");
+
     insertSurfacePointPairs();
+
     boundaryConform();
 
+
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Old Method
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -816,19 +814,19 @@ void Foam::CV2D::newPoints(const scalar relaxation)
     //         // in in the x-y directions
     //         vector2D cd0(1, 0);
 
-    //         if (controls_.relaxOrientation)
+    //         if (meshControls().relaxOrientation())
     //         {
     //             // Get the longest edge from the array and use as the primary
     //             // direction of the coordinate system of the "square" cell
     //             cd0 = edges[edgecd0i];
     //         }
 
-    //         if (controls_.nearWallAlignedDist > 0)
+    //         if (meshControls().nearWallAlignedDist() > 0)
     //         {
     //             pointIndexHit pHit = qSurf_.tree().findNearest
     //             (
     //                 toPoint3D(defVert0),
-    //                 controls_.nearWallAlignedDist2
+    //                 meshControls().nearWallAlignedDist2()
     //             );
 
     //             if (pHit.hit())
@@ -871,7 +869,7 @@ void Foam::CV2D::newPoints(const scalar relaxation)
     //                 // Set the weight for this edge contribution
     //                 scalar w = 1;
 
-    //                 if (controls_.squares)
+    //                 if (meshControls().squares())
     //                 {
     //                     w = magSqr(deltai.x()*ei.y() - deltai.y()*ei.x());
     //                     // alternative weights
@@ -1014,7 +1012,7 @@ void Foam::CV2D::extractPatches
 
 void Foam::CV2D::write() const
 {
-    if (controls_.writeFinalTriangulation)
+    if (meshControls().objOutput())
     {
         writeFaces("allFaces.obj", false);
         writeFaces("faces.obj", true);
@@ -1025,4 +1023,32 @@ void Foam::CV2D::write() const
 }
 
 
+void Foam::CV2D::write(const word& stage) const
+{
+    if (meshControls().objOutput())
+    {
+        Foam::mkDir(stage + "Faces");
+        Foam::mkDir(stage + "Triangles");
+
+        writeFaces
+        (
+            stage
+          + "Faces/allFaces_"
+          + runTime_.timeName()
+          + ".obj",
+            false
+        );
+
+        writeTriangles
+        (
+            stage
+          + "Triangles/allTriangles_"
+          + runTime_.timeName()
+          + ".obj",
+            false
+        );
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cv2DMesh/CV2D.H b/applications/utilities/mesh/generation/cv2DMesh/CV2D.H
index e295300b12..c875ed3a87 100644
--- a/applications/utilities/mesh/generation/cv2DMesh/CV2D.H
+++ b/applications/utilities/mesh/generation/cv2DMesh/CV2D.H
@@ -125,7 +125,7 @@ SourceFiles
 #include "Switch.H"
 #include "PackedBoolList.H"
 #include "EdgeMap.H"
-#include "controls.H"
+#include "cv2DControls.H"
 #include "tolerances.H"
 #include "meshTools.H"
 #include "triSurface.H"
@@ -165,14 +165,11 @@ private:
         conformationSurfaces qSurf_;
 
         //- Meshing controls
-        controls controls_;
+        cv2DControls controls_;
 
         //- The cell size control object
         cellSizeControlSurfaces cellSizeControl_;
 
-        //- Meshing tolerances
-        tolerances tols_;
-
         //- z-level
         scalar z_;
 
@@ -188,6 +185,9 @@ private:
         //  removing and insertin the surface point-pairs
         label startOfBoundaryConformPointPairs_;
 
+        //- Store the feature points
+        std::list<Vb> featurePoints_;
+
         //- Temporary storage for a dual-cell
         static const label maxNvert = 20;
         mutable point2D vertices[maxNvert+1];
@@ -210,6 +210,21 @@ private:
             const label type
         );
 
+        //- Insert point and return it's index
+        inline label insertPoint
+        (
+            const point2D& pt,
+            const label index,
+            const label type
+        );
+
+        inline label insertPoint
+        (
+            const Point& p,
+            const label index,
+            const label type
+        );
+
         inline bool insertMirrorPoint
         (
             const point2D& nearSurfPt,
@@ -228,9 +243,15 @@ private:
         //- Create the initial mesh from the bounding-box
         void insertBoundingBox();
 
+        //- Check if a point is within a line.
+        bool on2DLine(const point2D& p, const linePointRef& line);
+
         //- Insert point groups at the feature points.
         void insertFeaturePoints();
 
+        //- Re-insert point groups at the feature points.
+        void reinsertFeaturePoints();
+
         //- Insert point-pairs at the given set of points using the surface
         //  normals corresponding to the given set of surface triangles
         //  and write the inserted point locations to the given file.
@@ -239,6 +260,7 @@ private:
             const DynamicList<point2D>& nearSurfacePoints,
             const DynamicList<point2D>& surfacePoints,
             const DynamicList<label>& surfaceTris,
+            const DynamicList<label>& surfaceHits,
             const fileName fName
         );
 
@@ -266,7 +288,8 @@ private:
         (
             Triangulation::Finite_vertices_iterator& vit,
             const point2D& p,
-            const label trii
+            const label trii,
+            const label hitSurface
         );
 
         //- Insert point-pair at the best intersection point between the lines
@@ -292,9 +315,16 @@ private:
         void external_flip(Face_handle& f, int i);
         bool internal_flip(Face_handle& f, int i);
 
+        //- Write all the faces and all the triangles at a particular stage.
+        void write(const word& stage) const;
+
 
 public:
 
+        //- Runtime type information
+        ClassName("CV2D");
+
+
     // Constructors
 
         //- Construct for given surface
@@ -313,15 +343,13 @@ public:
 
         // Access
 
-            const controls& meshingControls() const
-            {
-                return controls_;
-            }
+            inline const cv2DControls& meshControls() const;
 
 
         // Conversion functions between point2D, point and Point
 
             inline const point2D& toPoint2D(const point&) const;
+            inline const point2DField toPoint2D(const pointField&) const;
             inline point toPoint3D(const point2D&) const;
 
             #ifdef CGAL_INEXACT
diff --git a/applications/utilities/mesh/generation/cv2DMesh/CV2DI.H b/applications/utilities/mesh/generation/cv2DMesh/CV2DI.H
index a389c32891..3cef5a4dbc 100644
--- a/applications/utilities/mesh/generation/cv2DMesh/CV2DI.H
+++ b/applications/utilities/mesh/generation/cv2DMesh/CV2DI.H
@@ -33,16 +33,40 @@ inline Foam::label Foam::CV2D::insertPoint
 {
     uint nVert = number_of_vertices();
 
-    Vertex_handle vh = insert(toPoint(p));
+    return insertPoint(toPoint(p), nVert, type);
+}
+
+
+inline Foam::label Foam::CV2D::insertPoint
+(
+    const point2D& p,
+    const label index,
+    const label type
+)
+{
+    return insertPoint(toPoint(p), index, type);
+}
+
+
+inline Foam::label Foam::CV2D::insertPoint
+(
+    const Point& p,
+    const label index,
+    const label type
+)
+{
+    uint nVert = number_of_vertices();
+
+    Vertex_handle vh = insert(p);
 
     if (nVert == number_of_vertices())
     {
         WarningIn("Foam::CV2D::insertPoint")
-            << "Failed to insert point " << p << endl;
+            << "Failed to insert point " << toPoint2D(p) << endl;
     }
     else
     {
-        vh->index() = nVert;
+        vh->index() = index;
         vh->type() = type;
     }
 
@@ -91,11 +115,26 @@ inline void Foam::CV2D::insertPointPair
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+inline const Foam::cv2DControls& Foam::CV2D::meshControls() const
+{
+    return controls_;
+}
+
 inline const Foam::point2D& Foam::CV2D::toPoint2D(const point& p) const
 {
     return reinterpret_cast<const point2D&>(p);
 }
 
+inline const Foam::point2DField Foam::CV2D::toPoint2D(const pointField& p) const
+{
+    point2DField temp(p.size());
+    forAll(temp, pointI)
+    {
+        temp[pointI] = point2D(p[pointI].x(), p[pointI].y());
+    }
+    return temp;
+}
+
 inline Foam::point Foam::CV2D::toPoint3D(const point2D& p) const
 {
     return point(p.x(), p.y(), z_);
@@ -151,8 +190,8 @@ inline void Foam::CV2D::movePoint(const Vertex_handle& vh, const Point& P)
     // move(vh, P);
     // vh->index() = i;
 
-    // vh->set_point(P);
-    // fast_restore_Delaunay(vh);
+     //vh->set_point(P);
+     //fast_restore_Delaunay(vh);
 }
 
 
diff --git a/applications/utilities/mesh/generation/cv2DMesh/Make/files b/applications/utilities/mesh/generation/cv2DMesh/Make/files
index 75e58f25a1..ddc8c9df92 100755
--- a/applications/utilities/mesh/generation/cv2DMesh/Make/files
+++ b/applications/utilities/mesh/generation/cv2DMesh/Make/files
@@ -1,8 +1,6 @@
 #include CGAL_FILES
 
 CV2D.C
-controls.C
-tolerances.C
 insertFeaturePoints.C
 insertSurfaceNearestPointPairs.C
 insertSurfaceNearPointPairs.C
@@ -11,5 +9,5 @@ CV2DIO.C
 shortEdgeFilter2D.C
 cv2DMesh.C
 
-EXE = $(FOAM_USER_APPBIN)/cv2DMesh
+EXE = $(FOAM_APPBIN)/cv2DMesh
 
diff --git a/applications/utilities/mesh/generation/cv2DMesh/Make/options b/applications/utilities/mesh/generation/cv2DMesh/Make/options
index 2344ae6a69..82659091b6 100755
--- a/applications/utilities/mesh/generation/cv2DMesh/Make/options
+++ b/applications/utilities/mesh/generation/cv2DMesh/Make/options
@@ -10,26 +10,25 @@ EXE_INC = \
     ${EXE_NDEBUG} \
     ${CGAL_INC} \
     -I$(FOAM_APP)/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/lnInclude \
+    -I$(FOAM_APP)/utilities/mesh/generation/extrude/extrudeModel/lnInclude \
     -IconformalVoronoi2DMesh/lnInclude \
+    -I$(FOAM_APP)/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/edgeMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(FOAM_APP)/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/lnInclude \
-    -I$(FOAM_APP)/utilities/mesh/generation/extrude/extrudeModel/lnInclude \
-    -I$(LIB_SRC)/triSurface/lnInclude -DFULLDEBUG -g -O0
+    -I$(LIB_SRC)/triSurface/lnInclude
 
 EXE_LIBS = \
-    -L$(FOAM_USER_LIBBIN) \
     $(CGAL_LIBS) \
+    -lextrude2DMesh \
+    -lextrudeModel \
+    -lcv2DMesh \
+    -lconformalVoronoiMesh \
     -lmeshTools \
     -lsurfMesh \
     -ledgeMesh \
     -ltriSurface \
     -ldynamicMesh \
-    -lextrude2DMesh \
-    -lextrudeModel \
-    -lconformalVoronoiMesh \
-    -lcv2DMesh \
     -ldecompositionMethods
diff --git a/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/Make/files b/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/Make/files
index df6cc21377..b05aa47e45 100755
--- a/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/Make/files
+++ b/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/Make/files
@@ -1,4 +1,4 @@
-patchTo2DpolyMesh/patchTo2DpolyMesh.C
+cv2DControls/cv2DControls.C
 
 LIB = $(FOAM_LIBBIN)/libcv2DMesh
 
diff --git a/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/Make/options b/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/Make/options
index 298044fab8..79be6f3a7d 100755
--- a/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/Make/options
+++ b/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/Make/options
@@ -1,14 +1,3 @@
-EXE_INC = \
-    -I$(FOAM_APP)/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/surfMesh/lnInclude \
-    -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/triSurface/lnInclude -DFULLDEBUG -g -O0
+EXE_INC =
 
-LIB_LIBS = \
-    -lmeshTools \
-    -lsurfMesh \
-    -ltriSurface \
-    -ldynamicMesh \
-    -lextrude2DMesh
+LIB_LIBS =
diff --git a/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/cv2DControls/cv2DControls.C b/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/cv2DControls/cv2DControls.C
new file mode 100644
index 0000000000..0b8e6fe288
--- /dev/null
+++ b/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/cv2DControls/cv2DControls.C
@@ -0,0 +1,149 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2007-2010 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*----------------------------------------------------------------------------*/
+
+#include "cv2DControls.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+
+Foam::cv2DControls::cv2DControls
+(
+    const dictionary& controlDict,
+    const boundBox& bb
+)
+:
+    dict_(controlDict),
+
+    motionControl_(controlDict.subDict("motionControl")),
+    conformationControl_(controlDict.subDict("surfaceConformation")),
+
+    minCellSize_(readScalar(motionControl_.lookup("minCellSize"))),
+    minCellSize2_(Foam::sqr(minCellSize_)),
+
+    maxQuadAngle_(readScalar(conformationControl_.lookup("maxQuadAngle"))),
+
+    nearWallAlignedDist_
+    (
+        readScalar(motionControl_.lookup("nearWallAlignedDist"))*minCellSize_
+    ),
+    nearWallAlignedDist2_(Foam::sqr(nearWallAlignedDist_)),
+
+    insertSurfaceNearestPointPairs_
+    (
+            conformationControl_.lookup("insertSurfaceNearestPointPairs")
+    ),
+    mirrorPoints_(conformationControl_.lookup("mirrorPoints")),
+    insertSurfaceNearPointPairs_
+    (
+        conformationControl_.lookup("insertSurfaceNearPointPairs")
+    ),
+
+    objOutput_(motionControl_.lookupOrDefault<Switch>("objOutput", false)),
+
+    randomiseInitialGrid_(conformationControl_.lookup("randomiseInitialGrid")),
+    randomPerturbation_
+    (
+        readScalar(conformationControl_.lookup("randomPerturbation"))
+    ),
+
+    maxBoundaryConformingIter_
+    (
+        readLabel(conformationControl_.lookup("maxBoundaryConformingIter"))
+    ),
+
+    span_
+    (
+        max(mag(bb.max().x()), mag(bb.min().x()))
+      + max(mag(bb.max().y()), mag(bb.min().y()))
+    ),
+    span2_(Foam::sqr(span_)),
+
+    minEdgeLen_
+    (
+        readScalar(conformationControl_.lookup("minEdgeLenCoeff"))
+       *minCellSize_
+    ),
+    minEdgeLen2_(Foam::sqr(minEdgeLen_)),
+
+    maxNotchLen_
+    (
+        readScalar(conformationControl_.lookup("maxNotchLenCoeff"))
+       *minCellSize_
+    ),
+    maxNotchLen2_(Foam::sqr(maxNotchLen_)),
+
+    minNearPointDist_
+    (
+        readScalar(conformationControl_.lookup("minNearPointDistCoeff"))
+       *minCellSize_
+    ),
+    minNearPointDist2_(Foam::sqr(minNearPointDist_)),
+
+    ppDist_
+    (
+        readScalar(conformationControl_.lookup("pointPairDistanceCoeff"))
+       *minCellSize_
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::cv2DControls::~cv2DControls()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::cv2DControls::write(Ostream& os) const
+{
+    os.indentLevel() = 1;
+    os.precision(2);
+    os.flags(ios_base::scientific);
+
+    os << nl << "Outputting CV2D Mesher controls:" << nl
+       << token::BEGIN_BLOCK << nl
+       << indent << "minCellSize2_         : " << minCellSize2_ << nl
+       << indent << "span_ / span2_        : " << span_ << " / " << span2_ << nl
+       << indent << "maxNotchLen2_         : " << maxNotchLen2_ << nl
+       << indent << "minNearPointDist2_    : " << minNearPointDist2_ << nl
+       << indent << "nearWallAlignedDist2_ : " << nearWallAlignedDist2_ << nl
+       << indent << "ppDist_               : " << ppDist_ << nl
+       << indent << "minEdgeLen2_          : " << minEdgeLen2_ << nl
+       << token::END_BLOCK << endl;
+}
+
+
+// * * * * * * * * * * * * * * IOStream operators  * * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const cv2DControls& s)
+{
+    s.write(os);
+    return os;
+}
+
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/cv2DControls/cv2DControls.H b/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/cv2DControls/cv2DControls.H
new file mode 100644
index 0000000000..98a9587a0b
--- /dev/null
+++ b/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/cv2DControls/cv2DControls.H
@@ -0,0 +1,257 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::cv2DControls
+
+Description
+    Controls for the 2D CV mesh generator.
+
+SourceFiles
+    cv2DControls.C
+    cv2DControlsI.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef cv2DControls_H
+#define cv2DControls_H
+
+#include "Switch.H"
+#include "dictionary.H"
+#include "boundBox.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class cv2DControls Declaration
+\*---------------------------------------------------------------------------*/
+
+class cv2DControls
+{
+    // Private data
+
+        //- Description of data_
+        const dictionary& dict_;
+
+        const dictionary& motionControl_;
+
+        const dictionary& conformationControl_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        cv2DControls(const cv2DControls&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const cv2DControls&);
+
+
+public:
+
+        // Controls
+
+            //- Minimum cell size below which protusions through the surface are
+            //  not split
+            scalar minCellSize_;
+
+            //- Square of minCellSize
+            scalar minCellSize2_;
+
+            //- Maximum quadrant angle allowed at a concave corner before
+            //  additional "mitering" lines are added
+            scalar maxQuadAngle_;
+
+            //- Near-wall region where cells are aligned with the wall
+            scalar nearWallAlignedDist_;
+
+            //- Square of nearWallAlignedDist
+            scalar nearWallAlignedDist2_;
+
+            //- Insert near-boundary point mirror or point-pairs
+            Switch insertSurfaceNearestPointPairs_;
+
+            //- Mirror near-boundary points rather than insert point-pairs
+            Switch mirrorPoints_;
+
+            //- Insert point-pairs vor dual-cell vertices very near the surface
+            Switch insertSurfaceNearPointPairs_;
+
+            Switch objOutput_;
+
+            Switch randomiseInitialGrid_;
+
+            scalar randomPerturbation_;
+
+            label maxBoundaryConformingIter_;
+
+
+        // Tolerances
+
+            //- Maximum cartesian span of the geometry
+            scalar span_;
+
+            //- Square of span
+            scalar span2_;
+
+            //- Minumum edge-length of the cell size below which protusions
+            //  through the surface are not split
+            scalar minEdgeLen_;
+
+            //- Square of minEdgeLen
+            scalar minEdgeLen2_;
+
+            //- Maximum notch size below which protusions into the surface are
+            //  not filled
+            scalar maxNotchLen_;
+
+            //- Square of maxNotchLen
+            scalar maxNotchLen2_;
+
+            //- The minimum distance alowed between a dual-cell vertex
+            // and the surface before a point-pair is introduced
+            scalar minNearPointDist_;
+
+            //- Square of minNearPoint
+            scalar minNearPointDist2_;
+
+            //- Distance between boundary conforming point-pairs
+            scalar ppDist_;
+
+            //- Square of ppDist
+            scalar ppDist2_;
+
+
+    // Constructors
+
+        cv2DControls
+        (
+            const dictionary& controlDict,
+            const boundBox& bb
+        );
+
+
+    //- Destructor
+    ~cv2DControls();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return the minimum cell size
+            inline scalar minCellSize() const;
+
+            //- Return the square of the minimum cell size
+            inline scalar minCellSize2() const;
+
+            //- Return the maximum quadrant angle
+            inline scalar maxQuadAngle() const;
+
+            //- Return number of layers to align with the nearest wall
+            inline scalar nearWallAlignedDist() const;
+
+            //- Return square of nearWallAlignedDist
+            inline scalar nearWallAlignedDist2() const;
+
+            //- Return insertSurfaceNearestPointPairs Switch
+            inline Switch insertSurfaceNearestPointPairs() const;
+
+            //- Return mirrorPoints Switch
+            inline Switch mirrorPoints() const;
+
+            //- Return insertSurfaceNearPointPairs Switch
+            inline Switch insertSurfaceNearPointPairs() const;
+
+            //- Return the objOutput Switch
+            inline Switch objOutput() const;
+
+            //- Return the randomise initial point layout Switch
+            inline Switch randomiseInitialGrid() const;
+
+            //- Return the random perturbation factor
+            inline scalar randomPerturbation() const;
+
+            //- Return the maximum number of boundary conformation iterations
+            inline label maxBoundaryConformingIter() const;
+
+            //- Return the span
+            inline scalar span() const;
+
+            //- Return the span squared
+            inline scalar span2() const;
+
+            //- Return the minEdgeLen
+            inline scalar minEdgeLen() const;
+
+            //- Return the minEdgeLen squared
+            inline scalar minEdgeLen2() const;
+
+            //- Return the maxNotchLen
+            inline scalar maxNotchLen() const;
+
+            //- Return the maxNotchLen squared
+            inline scalar maxNotchLen2() const;
+
+            //- Return the minNearPointDist
+            inline scalar minNearPointDist() const;
+
+            //- Return the minNearPointDist squared
+            inline scalar minNearPointDist2() const;
+
+            //- Return the ppDist
+            inline scalar ppDist() const;
+
+
+        // Write
+
+            //- Write controls to output stream.
+            void write(Ostream& os) const;
+
+            //- Ostream Operator
+            friend Ostream& operator<<
+            (
+                Ostream& os,
+                const cv2DControls& s
+            );
+
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "cv2DControlsI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* // 
diff --git a/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/cv2DControls/cv2DControlsI.H b/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/cv2DControls/cv2DControlsI.H
new file mode 100644
index 0000000000..3ff16735a8
--- /dev/null
+++ b/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/cv2DControls/cv2DControlsI.H
@@ -0,0 +1,153 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+
+inline Foam::scalar Foam::cv2DControls::minCellSize() const
+{
+    return minCellSize_;
+}
+
+
+inline Foam::scalar Foam::cv2DControls::minCellSize2() const
+{
+    return minCellSize2_;
+}
+
+
+inline Foam::scalar Foam::cv2DControls::maxQuadAngle() const
+{
+    return maxQuadAngle_;
+}
+
+
+inline Foam::scalar Foam::cv2DControls::nearWallAlignedDist() const
+{
+    return nearWallAlignedDist_;
+}
+
+
+inline Foam::scalar Foam::cv2DControls::nearWallAlignedDist2() const
+{
+    return nearWallAlignedDist2_;
+}
+
+
+inline Foam::Switch Foam::cv2DControls::insertSurfaceNearestPointPairs() const
+{
+    return insertSurfaceNearestPointPairs_;
+}
+
+
+inline Foam::Switch Foam::cv2DControls::mirrorPoints() const
+{
+    return mirrorPoints_;
+}
+
+
+inline Foam::Switch Foam::cv2DControls::insertSurfaceNearPointPairs() const
+{
+    return insertSurfaceNearPointPairs_;
+}
+
+
+inline Foam::Switch Foam::cv2DControls::objOutput() const
+{
+    return objOutput_;
+}
+
+
+inline Foam::Switch Foam::cv2DControls::randomiseInitialGrid() const
+{
+    return randomiseInitialGrid_;
+}
+
+
+inline Foam::scalar Foam::cv2DControls::randomPerturbation() const
+{
+    return randomPerturbation_;
+}
+
+
+inline Foam::label Foam::cv2DControls::maxBoundaryConformingIter() const
+{
+    return maxBoundaryConformingIter_;
+}
+
+
+inline Foam::scalar Foam::cv2DControls::span() const
+{
+    return span_;
+}
+
+
+inline Foam::scalar Foam::cv2DControls::span2() const
+{
+    return span2_;
+}
+
+
+inline Foam::scalar Foam::cv2DControls::minEdgeLen() const
+{
+    return minEdgeLen_;
+}
+
+
+inline Foam::scalar Foam::cv2DControls::minEdgeLen2() const
+{
+    return minEdgeLen2_;
+}
+
+
+inline Foam::scalar Foam::cv2DControls::maxNotchLen() const
+{
+    return maxNotchLen_;
+}
+
+
+inline Foam::scalar Foam::cv2DControls::maxNotchLen2() const
+{
+    return maxNotchLen2_;
+}
+
+
+inline Foam::scalar Foam::cv2DControls::minNearPointDist() const
+{
+    return minNearPointDist_;
+}
+
+
+inline Foam::scalar Foam::cv2DControls::minNearPointDist2() const
+{
+    return minNearPointDist2_;
+}
+
+
+inline Foam::scalar Foam::cv2DControls::ppDist() const
+{
+    return ppDist_;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cv2DMesh/controls.C b/applications/utilities/mesh/generation/cv2DMesh/controls.C
deleted file mode 100644
index fd90a2d2dd..0000000000
--- a/applications/utilities/mesh/generation/cv2DMesh/controls.C
+++ /dev/null
@@ -1,82 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2007-2010 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 3 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, see <http://www.gnu.org/licenses/>.
-
-\*----------------------------------------------------------------------------*/
-
-#include "controls.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::controls::controls(const dictionary& controlDict)
-:
-    dict_(controlDict),
-    minCellSize(readScalar(controlDict.lookup("minCellSize"))),
-    minCellSize2(Foam::sqr(minCellSize)),
-
-    featAngle(readScalar(controlDict.lookup("featureAngle"))),
-    maxQuadAngle(readScalar(controlDict.lookup("maxQuadAngle"))),
-    squares(controlDict.lookup("squares")),
-
-    nearWallAlignedDist
-    (
-        readScalar(controlDict.lookup("nearWallAlignedDist"))*minCellSize
-    ),
-    nearWallAlignedDist2(Foam::sqr(nearWallAlignedDist)),
-
-    relaxOrientation(controlDict.lookup("relaxOrientation")),
-
-    insertSurfaceNearestPointPairs
-    (
-        controlDict.lookup("insertSurfaceNearestPointPairs")
-    ),
-    mirrorPoints(controlDict.lookup("mirrorPoints")),
-    insertSurfaceNearPointPairs
-    (
-        controlDict.lookup("insertSurfaceNearPointPairs")
-    ),
-    writeInitialTriangulation(controlDict.lookup("writeInitialTriangulation")),
-    writeFeatureTriangulation(controlDict.lookup("writeFeatureTriangulation")),
-    writeNearestTriangulation(controlDict.lookup("writeNearestTriangulation")),
-    writeInsertedPointPairs(controlDict.lookup("writeInsertedPointPairs")),
-    writeFinalTriangulation(controlDict.lookup("writeFinalTriangulation")),
-    randomiseInitialGrid(controlDict.lookup("randomiseInitialGrid")),
-    randomPurturbation(readScalar(controlDict.lookup("randomPurturbation"))),
-    maxBoundaryConformingIter
-    (
-        readLabel(controlDict.lookup("maxBoundaryConformingIter"))
-    ),
-    relaxationFactorStart
-    (
-        readScalar(controlDict.lookup("relaxationFactorStart"))
-    ),
-    relaxationFactorEnd(readScalar(controlDict.lookup("relaxationFactorEnd")))
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::controls::~controls()
-{}
-
-
-// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cv2DMesh/controls.H b/applications/utilities/mesh/generation/cv2DMesh/controls.H
deleted file mode 100644
index 1e6c7d79cd..0000000000
--- a/applications/utilities/mesh/generation/cv2DMesh/controls.H
+++ /dev/null
@@ -1,164 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2011 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 3 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, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::controls
-
-Description
-    Controls for the 2D cv mesh generator.
-
-SourceFiles
-    controls.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef controls_H
-#define controls_H
-
-#include "Switch.H"
-#include "dictionary.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                          Class controls Declaration
-\*---------------------------------------------------------------------------*/
-
-class controls
-{
-    // Private data
-
-        //- Description of data_
-        const dictionary& dict_;
-
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        controls(const controls&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const controls&);
-
-
-public:
-
-        //- Minimum cell size below which protusions through the surface are
-        //  not split
-        scalar minCellSize;
-
-        //- Square of minCellSize
-        scalar minCellSize2;
-
-        //- The feature angle used to select corners to be
-        //  explicitly represented in the mesh.
-        //  0 = all features, 180 = no features
-        scalar featAngle;
-
-        //- Maximum quadrant angle allowed at a concave corner before
-        //  additional "mitering" lines are added
-        scalar maxQuadAngle;
-
-        //- Should the mesh be square-dominated or of unbiased hexagons
-        Switch squares;
-
-        //- Near-wall region where cells are aligned with the wall
-        scalar nearWallAlignedDist;
-
-        //- Square of nearWallAlignedDist
-        scalar nearWallAlignedDist2;
-
-        //- Chose if the cell orientation should relax during the iterations
-        //  or remain fixed to the x-y directions
-        Switch relaxOrientation;
-
-        //- Insert near-boundary point mirror or point-pairs
-        Switch insertSurfaceNearestPointPairs;
-
-        //- Mirror near-boundary points rather than insert point-pairs
-        Switch mirrorPoints;
-
-        //- Insert point-pairs vor dual-cell vertices very near the surface
-        Switch insertSurfaceNearPointPairs;
-
-        Switch writeInitialTriangulation;
-        Switch writeFeatureTriangulation;
-        Switch writeNearestTriangulation;
-        Switch writeInsertedPointPairs;
-        Switch writeFinalTriangulation;
-
-        Switch randomiseInitialGrid;
-        scalar randomPurturbation;
-
-        label maxBoundaryConformingIter;
-
-        //- Relaxation factor at the start of the iteration
-        scalar relaxationFactorStart;
-
-        //- Relaxation factor at the end of the iteration
-        scalar relaxationFactorEnd;
-
-
-    // Constructors
-
-        controls(const dictionary& controlDict);
-
-
-    //- Destructor
-    ~controls();
-
-
-    // Member Functions
-
-        // Access
-
-        // Check
-
-        // Edit
-
-        // Write
-
-
-    // Friend Functions
-
-    // Friend Operators
-
-    // IOstream Operators
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-//#include "controlsI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* // 
diff --git a/applications/utilities/mesh/generation/cv2DMesh/cv2DMesh.C b/applications/utilities/mesh/generation/cv2DMesh/cv2DMesh.C
index 4ed77e7933..804e6610b0 100644
--- a/applications/utilities/mesh/generation/cv2DMesh/cv2DMesh.C
+++ b/applications/utilities/mesh/generation/cv2DMesh/cv2DMesh.C
@@ -33,16 +33,16 @@ Description
 
 #include "CV2D.H"
 #include "argList.H"
-#include "IFstream.H"
 
 #include "MeshedSurfaces.H"
 #include "shortEdgeFilter2D.H"
 #include "extrude2DMesh.H"
 #include "polyMesh.H"
-#include "PatchTools.H"
-#include "patchTo2DpolyMesh.H"
+#include "patchToPoly2DMesh.H"
 #include "extrudeModel.H"
 #include "polyTopoChange.H"
+#include "edgeCollapser.H"
+#include "relaxationModel.H"
 
 using namespace Foam;
 
@@ -53,9 +53,10 @@ int main(int argc, char *argv[])
 {
     argList::noParallel();
     argList::validArgs.clear();
-    argList::validArgs.append("surface");
     argList::validOptions.insert("pointsFile", "<filename>");
 
+    #include "addOverwriteOption.H"
+
     #include "setRootCase.H"
     #include "createTime.H"
 
@@ -66,19 +67,16 @@ int main(int argc, char *argv[])
     dictionary extrusionDict(controlDict.subDict("extrusion"));
 
     Switch extrude(extrusionDict.lookup("extrude"));
-    label nIterations(readLabel(controlDict.lookup("nIterations")));
-    label sefDebug(shortEdgeFilterDict.lookupOrDefault<label>("debug", 0));
+    const bool overwrite = args.optionFound("overwrite");
 
-    // Read the surface to conform to
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-//    querySurface surf(args.args()[1]);
-//    surf.writeTreeOBJ();
-
-//    Info<< nl
-//        << "Read surface with " << surf.size() << " triangles from file "
-//        << args.args()[1] << nl << endl;
-
-   // surf.write("surface.obj");
+    autoPtr<relaxationModel> relax
+    (
+        relaxationModel::New
+        (
+            controlDict.subDict("motionControl"),
+            runTime
+        )
+    );
 
     // Read and triangulation
     // ~~~~~~~~~~~~~~~~~~~~~~
@@ -97,24 +95,13 @@ int main(int argc, char *argv[])
         mesh.boundaryConform();
     }
 
-    for (int iter=1; iter<=nIterations; iter++)
+    while (runTime.loop())
     {
-        Info<< nl
-            << "Relaxation iteration " << iter << nl
-            << "~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
+        Info<< nl << "Time = " << runTime.timeName() << endl;
 
-        scalar relax =
-           mesh.meshingControls().relaxationFactorStart
-         +
-         (
-            mesh.meshingControls().relaxationFactorEnd
-          - mesh.meshingControls().relaxationFactorStart
-         )
-         *scalar(iter)/scalar(nIterations);
+        Info<< "Relaxation = " << relax->relaxation() << endl;
 
-        Info<< "Relaxation = " << relax << endl;
-
-        mesh.newPoints(relax);
+        mesh.newPoints(relax->relaxation());
     }
 
     mesh.write();
@@ -122,22 +109,23 @@ int main(int argc, char *argv[])
     Info<< "Finished Delaunay in = "
         << runTime.cpuTimeIncrement() << " s." << endl;
 
+    Info<< "Begin filtering short edges:" << endl;
     shortEdgeFilter2D sef(mesh, shortEdgeFilterDict);
 
-    shortEdgeFilter2D::debug = sefDebug;
-
     sef.filter();
 
     Info<< "Meshed surface after edge filtering :" << endl;
     sef.fMesh().writeStats(Info);
 
-    Info<< "Write .obj file : MeshedSurface.obj" << endl;
+    Info<< "Write .obj file of the 2D mesh: MeshedSurface.obj" << endl;
     sef.fMesh().write("MeshedSurface.obj");
 
     Info<< "Finished filtering in = "
         << runTime.cpuTimeIncrement() << " s." << endl;
 
-    patchTo2DpolyMesh poly2DMesh
+    Info<< "Begin constructing a polyMesh:" << endl;
+
+    patchToPoly2DMesh poly2DMesh
     (
         sef.fMesh(),
         sef.patchNames(),
@@ -183,21 +171,70 @@ int main(int argc, char *argv[])
 
     if (extrude)
     {
-        // Point generator
-        autoPtr<extrudeModel> model(extrudeModel::New(extrusionDict));
+        Info<< "Begin extruding the polyMesh:" << endl;
 
-        extrude2DMesh extruder(pMesh, extrusionDict, model());
+        {
+            // Point generator
+            autoPtr<extrudeModel> model(extrudeModel::New(extrusionDict));
 
-        polyTopoChange meshMod(pMesh.boundaryMesh().size());
+            extrude2DMesh extruder(pMesh, extrusionDict, model());
 
-        extruder.addFrontBackPatches();
-        extruder.setRefinement(meshMod);
+            extruder.addFrontBackPatches();
 
-        autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(pMesh, false);
+            polyTopoChange meshMod(pMesh.boundaryMesh().size());
 
-        pMesh.updateMesh(morphMap);
+            extruder.setRefinement(meshMod);
+
+            autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(pMesh, false);
+
+            pMesh.updateMesh(morphMap);
+        }
+
+        {
+            edgeCollapser collapser(pMesh);
+
+            const edgeList& edges = pMesh.edges();
+            const pointField& points = pMesh.points();
+
+            const boundBox& bb = pMesh.bounds();
+            const scalar mergeDim = 1E-4 * bb.minDim();
+
+            forAll(edges, edgeI)
+            {
+                const edge& e = edges[edgeI];
+
+                scalar d = e.mag(points);
+
+                if (d < mergeDim)
+                {
+                    Info<< "Merging edge " << e << " since length " << d
+                        << " << " << mergeDim << endl;
+
+                    // Collapse edge to e[0]
+                    collapser.collapseEdge(edgeI, e[0]);
+                }
+            }
+
+            polyTopoChange meshModCollapse(pMesh);
+
+            collapser.setRefinement(meshModCollapse);
+
+            // Create a mesh from topo changes.
+            autoPtr<mapPolyMesh> morphMap =
+                meshModCollapse.changeMesh(pMesh, false);
+
+            pMesh.updateMesh(morphMap);
+        }
+
+        if (!overwrite)
+        {
+            runTime++;
+        }
+        else
+        {
+            pMesh.setInstance("constant");
+        }
 
-        pMesh.setInstance("constant");
     }
 
     pMesh.write();
diff --git a/applications/utilities/mesh/generation/cv2DMesh/cv2DMeshDict b/applications/utilities/mesh/generation/cv2DMesh/cv2DMeshDict
new file mode 100644
index 0000000000..0ef6d2a0c4
--- /dev/null
+++ b/applications/utilities/mesh/generation/cv2DMesh/cv2DMeshDict
@@ -0,0 +1,208 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+FoamFile
+{
+    version         2.0;
+    format          ascii;
+
+    root            "";
+    case            "";
+    instance        "";
+    local           "";
+
+    class           dictionary;
+    object          cv2DMeshDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+geometry
+{
+    laurence_clean_preciser.stl
+    {
+        name laurence_clean_preciser;
+        type closedTriSurfaceMesh;
+        //type triSurfaceMesh;
+    }
+//    refinementBox
+//    {
+//        type searchableBox;
+//        min (-0.5 0.35 -1000);
+//        max (-0.5 0.35 1000);
+//    }
+//    refinementSphere
+//    {
+//        type searchableSphere;
+//        centre (0.85 0.4 0.0);
+//        radius 0.01;
+//    }
+}
+
+surfaceConformation
+{
+    locationInMesh (-2.8 0.7 0.5);
+
+    pointPairDistanceCoeff    0.005;
+
+    minEdgeLenCoeff           0.005;
+
+    maxNotchLenCoeff          0.003;
+
+    minNearPointDistCoeff     0.0025;
+
+    maxQuadAngle 125;
+
+    // Insert near-boundary point mirror or point-pairs
+    insertSurfaceNearestPointPairs yes;
+
+    // Mirror near-boundary points rather than insert point-pairs
+    mirrorPoints   no;
+
+    // Insert point-pairs vor dual-cell vertices very near the surface
+    insertSurfaceNearPointPairs yes;
+
+    // Maximum number of iterations used in boundaryConform.
+    maxBoundaryConformingIter 5;
+
+    geometryToConformTo
+    {
+        laurence_clean_preciser
+        {
+            featureMethod extendedFeatureEdgeMesh;
+            extendedFeatureEdgeMesh "laurence_clean_preciser.extendedFeatureEdgeMesh";
+        }
+    }
+
+    additionalFeatures
+    {
+    }
+
+    // Choose if to randomise the initial grid created by insertGrid.
+    randomiseInitialGrid yes;
+
+    // Perturbation fraction, 1 = cell-size.
+    randomPerturbation   0.1;
+
+}
+
+
+motionControl
+{
+    defaultCellSize     0.05;
+
+    // Assign a priority to all requests for cell sizes, the highest overrules.
+    defaultPriority     0;
+
+    cellSizeControlGeometry
+    {
+        laurence_clean_preciser
+        {
+            priority              1;
+            mode                  bothSides;
+            cellSizeFunction      linearDistance;
+            linearDistanceCoeffs
+            {
+                distanceCellSize  0.05;
+                surfaceCellSize   0.01;
+                distance          0.5;
+            }
+            uniformCoeffs
+            {
+                cellSize          0.01;
+            }
+        }
+//        refinementBox
+//        {
+//            priority            1;
+//            mode                outside;
+//            cellSizeFunction    linearDistance;
+//            linearDistanceCoeffs
+//            {
+//                distanceCellSize 0.04;
+//                surfaceCellSize 0.005;
+//                distance        0.2;
+//            }
+//        }
+//        refinementSphere
+//        {
+//            priority            1;
+//            mode                outside;
+//            cellSizeFunction    linearDistance;
+//            linearDistanceCoeffs
+//            {
+//                distanceCellSize 0.04;
+//                surfaceCellSize 0.005;
+//                distance        0.2;
+//            }
+//        }
+    }
+
+    relaxationModel     adaptiveLinear;
+
+    adaptiveLinearCoeffs
+    {
+        relaxationStart 0.5;
+        relaxationEnd   0.0;
+    }
+
+    objOutput                   no;
+
+    // Near-wall region where cells are aligned with the wall specified as a number
+    // of cell layers
+    nearWallAlignedDist 3;
+
+}
+
+shortEdgeFilter
+{
+    // Factor to multiply the average of a face's edge lengths by.
+    // If an edge of that face is smaller than that value then delete it.
+    shortEdgeFilterFactor     0.2;
+
+    // Weighting for the lengths of edges that are attached to the boundaries.
+    // Used when calculating the length of an edge. Default 2.0.
+    edgeAttachedToBoundaryFactor 2.0;
+}
+
+extrusion
+{
+    extrude         on;
+
+    extrudeModel        linearDirection;
+    //extrudeModel        wedge;
+
+    patchInfo
+    {
+        //type empty;
+        //startFace
+    }
+
+    patchType           empty;
+    //patchType           wedge;
+
+    nLayers             1;
+
+    expansionRatio      1.0;    //0.9;
+
+    linearDirectionCoeffs
+    {
+        direction       (0 0 1);
+        thickness       0.1;
+    }
+
+    wedgeCoeffs
+    {
+        axisPt      (0 0 0);
+        axis        (1 0 0);
+        angle       10;
+    }
+
+    thickness       0.1;
+}
diff --git a/applications/utilities/mesh/generation/cv2DMesh/featuresMissingFrom2D.org b/applications/utilities/mesh/generation/cv2DMesh/featuresMissingFrom2D.org
new file mode 100644
index 0000000000..857c41f31c
--- /dev/null
+++ b/applications/utilities/mesh/generation/cv2DMesh/featuresMissingFrom2D.org
@@ -0,0 +1,3 @@
+* Displacement limiting
+
+http://en.wikipedia.org/wiki/Geometric_algebra
diff --git a/applications/utilities/mesh/generation/cv2DMesh/indexedVertex.H b/applications/utilities/mesh/generation/cv2DMesh/indexedVertex.H
index edb8d4c954..1559714f82 100644
--- a/applications/utilities/mesh/generation/cv2DMesh/indexedVertex.H
+++ b/applications/utilities/mesh/generation/cv2DMesh/indexedVertex.H
@@ -99,6 +99,13 @@ public:
         type_(INTERNAL_POINT)
     {}
 
+    indexedVertex(const Point& p, const int index, const int& type)
+    :
+        Vb(p),
+        index_(index),
+        type_(type)
+    {}
+
     indexedVertex(const Point& p, Face_handle f)
     :
         Vb(f, p),
@@ -245,6 +252,7 @@ public:
             || (v1.farPoint() || v1.ppSlave())
             || (v2.farPoint() || v2.ppSlave());
     }
+
 };
 
 
diff --git a/applications/utilities/mesh/generation/cv2DMesh/insertBoundaryConformPointPairs.C b/applications/utilities/mesh/generation/cv2DMesh/insertBoundaryConformPointPairs.C
index dec7a5d38c..7a8def28fe 100644
--- a/applications/utilities/mesh/generation/cv2DMesh/insertBoundaryConformPointPairs.C
+++ b/applications/utilities/mesh/generation/cv2DMesh/insertBoundaryConformPointPairs.C
@@ -31,12 +31,13 @@ void Foam::CV2D::insertPointPair
 (
     Triangulation::Finite_vertices_iterator& vit,
     const point2D& p,
-    const label trii
+    const label trii,
+    const label hitSurface
 )
 {
     if
     (
-        !controls_.mirrorPoints
+        !meshControls().mirrorPoints()
      || !insertMirrorPoint(toPoint2D(vit->point()), p)
     )
     {
@@ -48,7 +49,7 @@ void Foam::CV2D::insertPointPair
         );
 
         vectorField norm(1);
-        qSurf_.geometry()[trii].getNormal
+        qSurf_.geometry()[hitSurface].getNormal
         (
             List<pointIndexHit>(1, pHit),
             norm
@@ -56,7 +57,7 @@ void Foam::CV2D::insertPointPair
 
         insertPointPair
         (
-            tols_.ppDist,
+            meshControls().ppDist(),
             p,
             toPoint2D(norm[0])
         );
@@ -84,6 +85,7 @@ bool Foam::CV2D::insertPointPairAtIntersection
     bool found = false;
     point2D interPoint;
     label interTri = -1;
+    label interHitSurface = -1;
     scalar interDist2 = 0;
 
     Face_circulator fcStart = incident_faces(vit);
@@ -97,7 +99,7 @@ bool Foam::CV2D::insertPointPairAtIntersection
             pointIndexHit pHit;
             label hitSurface = -1;
 
-            qSurf_.findSurfaceAnyIntersection
+            qSurf_.findSurfaceNearestIntersection
             (
                 toPoint3D(defVert),
                 toPoint3D(vertices[vi]),
@@ -120,15 +122,16 @@ bool Foam::CV2D::insertPointPairAtIntersection
                     // vertex
                     if (boundaryTriangle(fc))
                     {
-                        mps2 = tols_.maxNotchLen2;
+                        mps2 = meshControls().maxNotchLen2();
                     }
 
                     if (dist2 > mps2)
                     {
                         found = true;
                         interPoint = toPoint2D(pHit.hitPoint());
-                        interTri = hitSurface;
+                        interTri = pHit.index();
                         interDist2 = dist2;
+                        interHitSurface = hitSurface;
                     }
                 }
             }
@@ -139,7 +142,7 @@ bool Foam::CV2D::insertPointPairAtIntersection
 
     if (found)
     {
-        insertPointPair(vit, interPoint, interTri);
+        insertPointPair(vit, interPoint, interTri, interHitSurface);
         return true;
     }
     else
@@ -217,7 +220,7 @@ Foam::label Foam::CV2D::insertBoundaryConformPointPairs
         // Convert triangle vertex to OpenFOAM point
         point2DFromPoint defVert = toPoint2D(vit->point());
 
-        scalar maxProtSize2 = tols_.maxNotchLen2;
+        scalar maxProtSize2 = meshControls().maxNotchLen2();
 
         if (vit->internalOrBoundaryPoint())
         {
@@ -239,7 +242,7 @@ Foam::label Foam::CV2D::insertBoundaryConformPointPairs
             }
 
             // If the dual-cell is very small reject refinement
-            if (areaT2 < tols_.minEdgeLen2) continue;
+            if (areaT2 < meshControls().minEdgeLen2()) continue;
 
             // Estimate the cell width
             scalar cellWidth = areaT2/perimeter;
@@ -248,25 +251,25 @@ Foam::label Foam::CV2D::insertBoundaryConformPointPairs
             // Check dimensions of dual-cell
             /*
             // Quick rejection of dual-cell refinement based on it's perimeter
-            if (perimeter < 2*tols_.minCellSize) continue;
+            if (perimeter < 2*meshControls().minCellSize()) continue;
 
             // Also check the area of the cell and reject refinement
             // if it is less than that allowed
-            if (areaT2 < tols_.minCellSize2) continue;
+            if (areaT2 < meshControls().minCellSize2()) continue;
 
             // Estimate the cell width and reject refinement if it is less than
             // that allowed
-            if (cellWidth < 0.5*tols_.minEdgeLen) continue;
+            if (cellWidth < 0.5*meshControls().minEdgeLen()) continue;
             */
 
             if
             (
-                perimeter > 2*controls_.minCellSize
-             && areaT2 > controls_.minCellSize2
-             && cellWidth > 0.5*tols_.minEdgeLen
+                perimeter > 2*meshControls().minCellSize()
+             && areaT2 > meshControls().minCellSize2()
+             && cellWidth > 0.5*meshControls().minEdgeLen()
             )
             {
-                maxProtSize2 = 0.25*tols_.maxNotchLen2;
+                maxProtSize2 = 0.25*meshControls().maxNotchLen2();
             }
         }
 
@@ -282,6 +285,7 @@ Foam::label Foam::CV2D::insertBoundaryConformPointPairs
 
 void Foam::CV2D::markNearBoundaryPoints()
 {
+    label count = 0;
     for
     (
         Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
@@ -299,7 +303,7 @@ void Foam::CV2D::markNearBoundaryPoints()
             qSurf_.findSurfaceNearest
             (
                 vert,
-                4*controls_.minCellSize2,
+                4*meshControls().minCellSize2(),
                 pHit,
                 hitSurface
             );
@@ -307,9 +311,12 @@ void Foam::CV2D::markNearBoundaryPoints()
             if (pHit.hit())
             {
                 vit->setNearBoundary();
+                ++count;
             }
         }
     }
+
+    Info<< count << " points marked as being near a boundary" << endl;
 }
 
 
diff --git a/applications/utilities/mesh/generation/cv2DMesh/insertFeaturePoints.C b/applications/utilities/mesh/generation/cv2DMesh/insertFeaturePoints.C
index a1512f66f2..5cc736dbfa 100644
--- a/applications/utilities/mesh/generation/cv2DMesh/insertFeaturePoints.C
+++ b/applications/utilities/mesh/generation/cv2DMesh/insertFeaturePoints.C
@@ -28,170 +28,378 @@ License
 #include "triSurfaceTools.H"
 #include "unitConversion.H"
 
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+
+bool Foam::CV2D::on2DLine(const point2D& p, const linePointRef& line)
+{
+    const point2D& a = toPoint2D(line.start());
+    const point2D& b = toPoint2D(line.end());
+
+    if
+    (
+        p.x() < min(a.x(), b.x())
+     || p.x() > max(a.x(), b.x())
+     || p.y() < min(a.y(), b.y())
+     || p.y() > max(a.y(), b.y())
+    )
+    {
+        return false;
+    }
+
+    return true;
+}
+
+
 // Create feature points/edges by creating a triplet in the corner.
 // (this triplet will have as its circumcentre the feature)
 void Foam::CV2D::insertFeaturePoints()
 {
-    //labelList featEdges(qSurf_.extractFeatures2D(controls_.featAngle));
+    featurePoints_.clear();
+    label nVert = number_of_vertices();
 
-  /*  const PtrList<extendedFeatureEdgeMesh>& feMeshes
+    const PtrList<extendedFeatureEdgeMesh>& feMeshes
     (
         qSurf_.features()
     );
 
-   // const pointField& localPts = qSurf_.localPoints();
+    if (feMeshes.empty())
+    {
+        WarningIn("CV2D::insertFeaturePoints")
+            << "Extended Feature Edge Mesh is empty so no feature points will "
+            << "be found." << nl
+            << "    Use: featureMethod extendedFeatureEdgeMesh;" << nl
+            << endl;
+    }
 
     forAll(feMeshes, i)
     {
         const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
+        const edgeList& edges = feMesh.edges();
+        const pointField& points = feMesh.points();
 
-        // Loop over convex points
-        for
-        (
-            label ptI = feMesh.convexStart();
-            ptI < feMesh.concaveStart();
-            ptI++
-        )
+        if (debug)
         {
+            label nConvex = feMesh.concaveStart() - feMesh.convexStart();
+            label nConcave = feMesh.mixedStart() - feMesh.concaveStart();
+            label nMixed = feMesh.nonFeatureStart() - feMesh.mixedStart();
+            label nExternal = feMesh.internalStart() - feMesh.externalStart();
+            label nInternal = feMesh.flatStart() - feMesh.internalStart();
+            label nFlat = feMesh.openStart() - feMesh.flatStart();
+            label nOpen = feMesh.multipleStart() - feMesh.openStart();
+            label nMultiple = edges.size() - feMesh.multipleStart();
 
-            const Foam::point& featPt = feMesh.points()[ptI];
+            Info<< "Inserting Feature Points:" << nl
+                << "    Convex points: " << nConvex << nl
+                << "   Concave points: " << nConcave << nl
+                << "     Mixed points: " << nMixed << nl
+                << "   External edges: " << nExternal << nl
+                << "   Internal edges: " << nInternal << nl
+                << "       Flat edges: " << nFlat << nl
+                << "       Open edges: " << nOpen << nl
+                << "   Multiple edges: " << nMultiple << endl;
         }
 
-        // Loop over concave points
-        for
-        (
-            label ptI = feMesh.concaveStart();
-            ptI < feMesh.mixedStart();
-            ptI++
-        )
-        {
+        // Args: (base point, normal)
+        // @todo allow user to input this
+        plane zPlane(vector(0, 0, z_), vector(0, 0, 1));
 
+        if (debug)
+        {
+            Info<< "    plane: " << zPlane << " " << z_ << endl;
         }
 
-        // Loop over mixed points
-        for
-        (
-            label ptI = feMesh.mixedStart();
-            ptI < feMesh.nonFeatureStart();
-            ptI++
-        )
+        forAll(edges, edgeI)
         {
+            const edge& e = feMesh.edges()[edgeI];
 
-        }
+            const point& ep0 = points[e.start()];
+            const point& ep1 = points[e.end()];
 
-        //label edgeI = featEdges[i];
-        //const edge& featEdge = qSurf_.edges()[edgeI];
+            const linePointRef line(ep0, ep1);
 
-        // Get the feature point as the mid-point of the edge and convert to 2D
-        //point2D featPt = toPoint2D(featEdge.centre(qSurf_.localPoints()));
+            scalar intersect = zPlane.lineIntersect(line);
 
-        // Pick up the two faces adjacent to the feature edge
-        const labelList& eFaces = qSurf_.edgeFaces()[edgeI];
+            point2D featPoint = toPoint2D(intersect * (ep1 - ep0) + ep0);
 
-        label faceA = eFaces[0];
-        vector2D nA = toPoint2D(qSurf_.faceNormals()[faceA]);
-
-        label faceB = eFaces[1];
-        vector2D nB = toPoint2D(qSurf_.faceNormals()[faceB]);
-
-        // Intersect planes parallel to faceA and faceB offset by ppDist.
-        plane planeA(toPoint3D(featPt - tols_.ppDist*nA), toPoint3D(nA));
-        plane planeB(toPoint3D(featPt - tols_.ppDist*nB), toPoint3D(nB));
-        plane::ray interLine(planeA.planeIntersect(planeB));
-
-        // The reference point is where this line intersects the z_ plane
-        point2D refPt = toPoint2D
-        (
-            interLine.refPoint()
-          + ((z_ - interLine.refPoint().z())/interLine.dir().z())
-           *interLine.dir()
-        );
-
-        point2D faceAVert = toPoint2D
-        (
-            localPts[triSurfaceTools::oppositeVertex(qSurf_, faceA, edgeI)]
-        );
-
-        // Determine convex or concave angle
-        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
-
-            point2D reflectedA = refPt + 2*((featPt - refPt) & nA)*nA;
-            insertPoint(reflectedA, masterPtIndex);
-
-            point2D reflectedB = refPt + 2*((featPt - refPt) & nB)*nB;
-            insertPoint(reflectedB, masterPtIndex);
-        }
-        else
-        {
-            // Concave. master and reflected points inside the domain.
-            // Generate reflected master to be outside.
-            point2D reflMasterPt = refPt + 2*(featPt - refPt);
-
-            // Reflect refPt in both faces.
-            point2D reflectedA =
-                reflMasterPt + 2*((featPt - reflMasterPt) & nA)*nA;
-
-            point2D reflectedB =
-                reflMasterPt + 2*((featPt - reflMasterPt) & nB)*nB;
-
-            // Total angle around the concave feature
-            scalar totalAngle =
-                radToDeg(constant::mathematical::pi + acos(mag(nA & nB)));
-
-            // 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)
+            if (on2DLine(featPoint, line))
             {
-                // One additinal point is the reflection of the slave point,
-                // i.e. the original reference point
-                insertPoint(refPt, reflectedMaster);
+                vector2DField fpn = toPoint2D(feMesh.edgeNormals(edgeI));
+
+                vector2D cornerNormal = sum(fpn);
+                cornerNormal /= mag(cornerNormal);
+
+                if (debug)
+                {
+                    Info<< nl << "    line: " << line << nl
+                        << "        vec: " << line.vec() << nl
+                        << "    featurePoint: " << featPoint << nl
+                        << "    line length: " << line.mag() << nl
+                        << "    intersect: " << intersect << endl;
+                }
+
+                if
+                (
+                    feMesh.getEdgeStatus(edgeI)
+                 == extendedFeatureEdgeMesh::EXTERNAL
+                )
+                {
+                    // Convex Point
+                    Foam::point2D internalPt =
+                        featPoint - meshControls().ppDist()*cornerNormal;
+
+                    if (debug)
+                    {
+                        Info<< "PREC: " << internalPt << nl
+                            << "    : " << featPoint << nl
+                            << "    : " << meshControls().ppDist() << nl
+                            << "    : " << cornerNormal << endl;
+                    }
+
+                    featurePoints_.push_back
+                    (
+                        Vb
+                        (
+                            toPoint(internalPt),
+                            nVert,
+                            nVert + 1
+                        )
+                    );
+                    label masterPtIndex = nVert++;
+
+                    forAll(fpn, nI)
+                    {
+                        const vector n3D(fpn[nI][0], fpn[nI][1], 0.0);
+
+                        plane planeN = plane(toPoint3D(featPoint), n3D);
+
+                        Foam::point2D externalPt =
+                            internalPt
+                          + (
+                                2.0
+                              * planeN.distance(toPoint3D(internalPt))
+                              * fpn[nI]
+                            );
+
+                        if (debug)
+                        {
+                            Info<< "PREC1: " << externalPt << nl
+                                << "    : "<< n3D << nl
+                                << "    : "<< planeN.distance(toPoint3D(internalPt)) << nl
+                                << "    : "<< planeN.normal() << nl
+                                << "    : "<< planeN.refPoint() << endl;
+                        }
+
+                        featurePoints_.push_back
+                        (
+                            Vb
+                            (
+                                toPoint(externalPt),
+                                nVert++,
+                                masterPtIndex
+                            )
+                        );
+
+                        if (debug)
+                        {
+                            Info<< "  side point: " << externalPt << endl;
+                        }
+                    }
+
+                    if (debug)
+                    {
+                        Info<< "Convex Point: " << featPoint << nl
+                            << "  corner norm: " << cornerNormal << nl
+                            << "    reference: " << internalPt << endl;
+                    }
+                }
+                else if
+                (
+                    feMesh.getEdgeStatus(edgeI)
+                 == extendedFeatureEdgeMesh::INTERNAL
+                )
+                {
+                    // Concave Point
+                    Foam::point2D externalPt =
+                        featPoint + meshControls().ppDist()*cornerNormal;
+
+                    Foam::point2D refPt =
+                        featPoint - meshControls().ppDist()*cornerNormal;
+
+                    label slavePointIndex = 0;
+
+                    scalar totalAngle =
+                        radToDeg
+                        (
+                            constant::mathematical::pi
+                          + acos(mag(fpn[0] & fpn[1]))
+                        );
+
+                    // Number of quadrants the angle should be split into
+                    int nQuads = int(totalAngle/meshControls().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 = nVert + 2 + nAddPoints;
+
+                    if (debug)
+                    {
+                        Info<< "Concave Point: " <<  featPoint << nl
+                            << "  corner norm: " << cornerNormal << nl
+                            << "     external: " << externalPt << nl
+                            << "    reference: " << refPt << nl
+                            << "        angle: " << totalAngle << nl
+                            << "       nQuads: " << nQuads << nl
+                            << "   nAddPoints: " << nAddPoints << endl;
+                    }
+
+                    forAll(fpn, nI)
+                    {
+                        const vector n3D(fpn[nI][0], fpn[nI][1], 0.0);
+
+                        plane planeN = plane(toPoint3D(featPoint), n3D);
+
+                        Foam::point2D internalPt =
+                            externalPt
+                          - (
+                                2.0
+                              * planeN.distance(toPoint3D(externalPt))
+                              * fpn[nI]
+                            );
+
+                        featurePoints_.push_back
+                        (
+                            Vb
+                            (
+                                toPoint(internalPt),
+                                nVert,
+                                reflectedMaster
+                            )
+                        );
+                        slavePointIndex = nVert++;
+
+                        if (debug)
+                        {
+                            Info<< "Internal Point: " <<  internalPt << endl;
+                        }
+                    }
+
+                    if (nAddPoints == 1)
+                    {
+                        // One additional point is the reflection of the slave point,
+                        // i.e. the original reference point
+                        featurePoints_.push_back
+                        (
+                            Vb
+                            (
+                                toPoint(refPt),
+                                nVert++,
+                                reflectedMaster
+                            )
+                        );
+
+                        if (debug)
+                        {
+                           Info<< "ref Point: " <<  refPt << endl;
+                        }
+                    }
+                    else if (nAddPoints == 2)
+                    {
+                       point2D reflectedAa =
+                           refPt - ((featPoint - externalPt) & fpn[1])*fpn[1];
+
+                       featurePoints_.push_back
+                       (
+                           Vb
+                           (
+                               toPoint(reflectedAa),
+                               nVert++,
+                               reflectedMaster
+                           )
+                       );
+
+                       point2D reflectedBb =
+                           refPt - ((featPoint - externalPt) & fpn[0])*fpn[0];
+
+                       featurePoints_.push_back
+                       (
+                           Vb
+                           (
+                               toPoint(reflectedBb),
+                               nVert++,
+                               reflectedMaster
+                           )
+                       );
+
+                       if (debug)
+                       {
+                           Info<< "refA Point: " <<  reflectedAa << nl
+                               << "refb Point: " <<  reflectedBb << endl;
+                       }
+                    }
+
+                    featurePoints_.push_back
+                    (
+                        Vb
+                        (
+                            toPoint(externalPt),
+                            nVert++,
+                            slavePointIndex
+                        )
+                    );
+                }
+                else
+                {
+                    WarningIn("void Foam::CV2D::insertFeaturePoints()")
+                        << "Feature Edge " << edges[edgeI] << nl
+                        << "    points(" << points[edges[edgeI].start()]
+                        << ", " << points[edges[edgeI].end()] << ")" << nl
+                        << "    is not labelled as either concave or convex, it"
+                        << " is labelled as (#2 = flat): "
+                        << feMesh.getEdgeStatus(edgeI) << endl;
+                }
             }
-            else if (nAddPoints == 2)
+            else
             {
-                point2D reflectedAa = refPt - ((featPt - reflMasterPt) & nB)*nB;
-                insertPoint(reflectedAa, reflectedMaster);
-
-                point2D reflectedBb = refPt - ((featPt - reflMasterPt) & nA)*nA;
-                insertPoint(reflectedBb, reflectedMaster);
+                WarningIn("void Foam::CV2D::insertFeaturePoints()")
+                      << "Point " << featPoint << " is not on the line "
+                      << line << endl;
             }
-
-            // Slave is outside.
-            insertPoint(reflMasterPt, reflectedAI);
         }
     }
 
-    if (controls_.writeFeatureTriangulation)
+    // Insert the feature points.
+    reinsertFeaturePoints();
+
+    if (meshControls().objOutput())
     {
         writePoints("feat_allPoints.obj", false);
         writeFaces("feat_allFaces.obj", false);
         writeFaces("feat_faces.obj", true);
         writeTriangles("feat_triangles.obj", true);
-    }*/
+    }
 }
 
 
+void Foam::CV2D::reinsertFeaturePoints()
+{
+    for
+    (
+        std::list<Vb>::iterator vit=featurePoints_.begin();
+        vit != featurePoints_.end();
+        ++vit
+    )
+    {
+        insertPoint
+        (
+            toPoint2D(vit->point()),
+            vit->index(),
+            vit->type()
+        );
+    }
+}
 // ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cv2DMesh/insertSurfaceNearPointPairs.C b/applications/utilities/mesh/generation/cv2DMesh/insertSurfaceNearPointPairs.C
index 81349b8df5..dc5dc16cac 100644
--- a/applications/utilities/mesh/generation/cv2DMesh/insertSurfaceNearPointPairs.C
+++ b/applications/utilities/mesh/generation/cv2DMesh/insertSurfaceNearPointPairs.C
@@ -24,7 +24,6 @@ License
 \*----------------------------------------------------------------------------*/
 
 #include "CV2D.H"
-#include "treeDataTriSurface.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -51,7 +50,7 @@ void Foam::CV2D::insertSurfaceNearPointPairs()
 
             // Check that the two triangle vertices are further apart than the
             // minimum cell size
-            if (magSqr(v1 - v0) > controls_.minCellSize2)
+            if (magSqr(v1 - v0) > meshControls().minCellSize2())
             {
                 point2D e0(toPoint2D(circumcenter(eit->first)));
 
@@ -63,7 +62,7 @@ void Foam::CV2D::insertSurfaceNearPointPairs()
                 // Calculate the length^2 of the edge normal to the surface
                 scalar edgeLen2 = magSqr(e0 - e1);
 
-                if (edgeLen2 < tols_.minNearPointDist2)
+                if (edgeLen2 < meshControls().minNearPointDist2())
                 {
                     pointIndexHit pHit;
                     label hitSurface = -1;
@@ -71,7 +70,7 @@ void Foam::CV2D::insertSurfaceNearPointPairs()
                     qSurf_.findSurfaceNearest
                     (
                         toPoint3D(e0),
-                        tols_.minEdgeLen2,
+                        meshControls().minEdgeLen2(),
                         pHit,
                         hitSurface
                     );
@@ -87,7 +86,7 @@ void Foam::CV2D::insertSurfaceNearPointPairs()
 
                         insertPointPair
                         (
-                            tols_.ppDist,
+                            meshControls().ppDist(),
                             toPoint2D(pHit.hitPoint()),
                             toPoint2D(norm[0])
                         );
@@ -95,7 +94,7 @@ void Foam::CV2D::insertSurfaceNearPointPairs()
                         nNearPoints++;
 
                         // Correct the edge iterator for the change in the 
-                        // number od edges following the point-pair insertion
+                        // number of edges following the point-pair insertion
                         eit = Finite_edges_iterator
                         (
                             finite_edges_end().base(),
diff --git a/applications/utilities/mesh/generation/cv2DMesh/insertSurfaceNearestPointPairs.C b/applications/utilities/mesh/generation/cv2DMesh/insertSurfaceNearestPointPairs.C
index 63ba5061d9..5f9b278ca7 100644
--- a/applications/utilities/mesh/generation/cv2DMesh/insertSurfaceNearestPointPairs.C
+++ b/applications/utilities/mesh/generation/cv2DMesh/insertSurfaceNearestPointPairs.C
@@ -55,7 +55,7 @@ bool Foam::CV2D::dualCellSurfaceIntersection
                 return true;
             }
 
-            if (magSqr(e1 - e0) > tols_.minEdgeLen2)
+            if (magSqr(e1 - e0) > meshControls().minEdgeLen2())
             {
                 if (qSurf_.findSurfaceAnyIntersection(e0, e1))
                 {
@@ -75,10 +75,11 @@ void Foam::CV2D::insertPointPairs
     const DynamicList<point2D>& nearSurfacePoints,
     const DynamicList<point2D>& surfacePoints,
     const DynamicList<label>& surfaceTris,
+    const DynamicList<label>& surfaceHits,
     const fileName fName
 )
 {
-    if (controls_.mirrorPoints)
+    if (meshControls().mirrorPoints())
     {
         forAll(surfacePoints, ppi)
         {
@@ -101,7 +102,7 @@ void Foam::CV2D::insertPointPairs
             );
 
             vectorField norm(1);
-            qSurf_.geometry()[surfaceTris[ppi]].getNormal
+            qSurf_.geometry()[surfaceHits[ppi]].getNormal
             (
                 List<pointIndexHit>(1, pHit),
                 norm
@@ -109,7 +110,7 @@ void Foam::CV2D::insertPointPairs
 
             insertPointPair
             (
-                tols_.ppDist,
+                meshControls().ppDist(),
                 surfacePoints[ppi],
                 toPoint2D(norm[0])
             );
@@ -118,7 +119,7 @@ void Foam::CV2D::insertPointPairs
 
     Info<< surfacePoints.size() << " point-pairs inserted" << endl;
 
-    if (controls_.writeInsertedPointPairs)
+    if (meshControls().objOutput())
     {
         OFstream str(fName);
         label vertI = 0;
@@ -140,15 +141,17 @@ void Foam::CV2D::insertSurfaceNearestPointPairs()
 {
     Info<< "insertSurfaceNearestPointPairs: ";
 
-    label nSurfacePointsEst = min
-    (
-        number_of_vertices(),
-        size_t(10*sqrt(scalar(number_of_vertices())))
-    );
+    label nSurfacePointsEst =
+        min
+        (
+            number_of_vertices(),
+            size_t(10*sqrt(scalar(number_of_vertices())))
+        );
 
     DynamicList<point2D> nearSurfacePoints(nSurfacePointsEst);
     DynamicList<point2D> surfacePoints(nSurfacePointsEst);
     DynamicList<label> surfaceTris(nSurfacePointsEst);
+    DynamicList<label> surfaceHits(nSurfacePointsEst);
 
     // Local references to surface mesh addressing
 //    const pointField& localPoints = qSurf_.localPoints();
@@ -173,7 +176,7 @@ void Foam::CV2D::insertSurfaceNearestPointPairs()
             qSurf_.findSurfaceNearest
             (
                 toPoint3D(vert),
-                4*controls_.minCellSize2,
+                4*meshControls().minCellSize2(),
                 pHit,
                 hitSurface
             );
@@ -216,11 +219,12 @@ void Foam::CV2D::insertSurfaceNearestPointPairs()
 //                    }
 //                }
 
-//                if (!internalFeatureEdge && dualCellSurfaceIntersection(vit))
+                if (dualCellSurfaceIntersection(vit)) //&& !internalFeatureEdge)
                 {
                     nearSurfacePoints.append(vert);
                     surfacePoints.append(toPoint2D(pHit.hitPoint()));
-                    surfaceTris.append(hitSurface);
+                    surfaceTris.append(pHit.index());
+                    surfaceHits.append(hitSurface);
                 }
             }
         }
@@ -231,6 +235,7 @@ void Foam::CV2D::insertSurfaceNearestPointPairs()
         nearSurfacePoints,
         surfacePoints,
         surfaceTris,
+        surfaceHits,
         "surfaceNearestIntersections.obj"
     );
 }
diff --git a/applications/utilities/mesh/generation/cv2DMesh/querySurface.C b/applications/utilities/mesh/generation/cv2DMesh/querySurface.C
deleted file mode 100644
index 6a780c7a0d..0000000000
--- a/applications/utilities/mesh/generation/cv2DMesh/querySurface.C
+++ /dev/null
@@ -1,237 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2007-2010 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 3 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, see <http://www.gnu.org/licenses/>.
-
-\*----------------------------------------------------------------------------*/
-
-#include "querySurface.H"
-#include "unitConversion.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::querySurface::querySurface(const fileName& surfaceFileName)
-:
-    triSurface(surfaceFileName),
-    rndGen_(12345),
-    bb_(localPoints()),
-    tree_
-    (
-        treeDataTriSurface
-        (
-            *this,
-            indexedOctree<treeDataTriSurface>::perturbTol()
-        ),
-        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(degToRad(featAngle));
-
-    const labelListList& edgeFaces = this->edgeFaces();
-    const pointField& points = this->points();
-    const edgeList& edges = this->edges();
-    const vectorField& faceNormals = this->faceNormals();
-    const labelList& meshPoints = this->meshPoints();
-
-    DynamicList<label> featEdges(edges.size());
-
-    forAll(edgeFaces, edgeI)
-    {
-        const edge& e = edges[edgeI];
-
-        point p =
-            points[meshPoints[e.end()]]
-          - points[meshPoints[e.start()]];
-
-        if (magSqr(p & 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();
-}
-
-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/cv2DMesh/querySurface.H b/applications/utilities/mesh/generation/cv2DMesh/querySurface.H
deleted file mode 100644
index b5d9ad0185..0000000000
--- a/applications/utilities/mesh/generation/cv2DMesh/querySurface.H
+++ /dev/null
@@ -1,149 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2007-2010 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 3 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, see <http://www.gnu.org/licenses/>.
-
-Class
-    querySurface
-
-Description
-    Searchable triSurface using an octree to speed-up queries.
-
-SourceFiles
-    querySurface.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef querySurface_H
-#define querySurface_H
-
-#include "triSurface.H"
-#include "treeDataTriSurface.H"
-#include "indexedOctree.H"
-#include "Random.H"
-#include "meshTools.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                           Class querySurface Declaration
-\*---------------------------------------------------------------------------*/
-
-class querySurface
-:
-    public triSurface
-{
-    // Private data
-
-        Random rndGen_;
-
-        // Bounding box of surface. Used for relative tolerances.
-        treeBoundBox bb_;
-
-        // Search engine on surface
-        indexedOctree<treeDataTriSurface> tree_;
-
-
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        querySurface(const querySurface&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const querySurface&);
-
-
-public:
-
-    // Constructors
-
-        //- Construct given file name of the surface
-        querySurface(const fileName& surfaceFileName);
-
-
-    // Destructor
-
-        ~querySurface();
-
-
-    // Member Functions
-
-        // Access
-
-            const treeBoundBox& bb() const
-            {
-                return bb_;
-            }
-
-            const indexedOctree<treeDataTriSurface>& tree() const
-            {
-                return tree_;
-            }
-
-
-        // Query
-
-            //- Extract feature edges/points(2D)
-            //  using the given feature angle in deg
-            labelList extractFeatures2D(const scalar featAngle) const;
-
-            //- Returns inside, outside or mixed (= close to surface)
-            indexedOctree<Foam::treeDataTriSurface>::volumeType insideOutside
-            (
-                const scalar searchSpan2,
-                const point& pt
-            ) const;
-
-            //- Check if point is inside surface
-            bool inside(const point& pt) const;
-
-            //- Check if point is outside surface
-            bool outside(const point& pt) const;
-
-            //- Check if point is inside surface by at least dist2
-            bool wellInside(const point& pt, const scalar dist2) const;
-
-            //- Check if point is outside surface by at least dist2
-            bool wellOutside(const point& pt, const scalar dist2) const;
-
-
-        // Write
-
-            void writeTreeOBJ() const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-//#include "querySurfaceI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cv2DMesh/shortEdgeFilter2D.C b/applications/utilities/mesh/generation/cv2DMesh/shortEdgeFilter2D.C
index 1ef698f3af..6e50760bb1 100644
--- a/applications/utilities/mesh/generation/cv2DMesh/shortEdgeFilter2D.C
+++ b/applications/utilities/mesh/generation/cv2DMesh/shortEdgeFilter2D.C
@@ -78,7 +78,12 @@ Foam::shortEdgeFilter2D::shortEdgeFilter2D
     Info<< "Meshed surface stats before edge filtering :" << endl;
     ms_.writeStats(Info);
 
-    ms_.write("MeshedSurface_preFilter.obj");
+    if (debug)
+    {
+        writeInfo(Info);
+
+        ms_.write("MeshedSurface_preFilter.obj");
+    }
 }
 
 
@@ -93,11 +98,6 @@ Foam::shortEdgeFilter2D::~shortEdgeFilter2D()
 void
 Foam::shortEdgeFilter2D::filter()
 {
-    if (debug)
-    {
-        writeInfo(Info);
-    }
-
     // These are global indices.
     const pointField& points = ms_.points();
     const edgeList& edges = ms_.edges();
@@ -164,11 +164,12 @@ Foam::shortEdgeFilter2D::filter()
         const label startVertex = e.start();
         const label endVertex = e.end();
 
-        scalar edgeLength = mag
-        (
-            points[meshPoints[e.start()]]
-           -points[meshPoints[e.end()]]
-        );
+        scalar edgeLength =
+            mag
+            (
+                points[meshPoints[e.start()]]
+              - points[meshPoints[e.end()]]
+            );
 
         if (edgeAttachedToBoundary[edgeI])
         {
@@ -185,11 +186,12 @@ Foam::shortEdgeFilter2D::filter()
             const edge& psE = edges[psEdges[psEdgeI]];
             if (edgeI != psEdges[psEdgeI])
             {
-                shortEdgeFilterValue += mag
-                (
-                    points[meshPoints[psE.start()]]
-                   -points[meshPoints[psE.end()]]
-                );
+                shortEdgeFilterValue +=
+                    mag
+                    (
+                        points[meshPoints[psE.start()]]
+                       -points[meshPoints[psE.end()]]
+                    );
             }
         }
 
@@ -198,11 +200,12 @@ Foam::shortEdgeFilter2D::filter()
             const edge& peE = edges[peEdges[peEdgeI]];
             if (edgeI != peEdges[peEdgeI])
             {
-                shortEdgeFilterValue += mag
-                (
-                    points[meshPoints[peE.start()]]
-                   -points[meshPoints[peE.end()]]
-                );
+                shortEdgeFilterValue +=
+                    mag
+                    (
+                        points[meshPoints[peE.start()]]
+                       -points[meshPoints[peE.end()]]
+                    );
             }
         }
 
diff --git a/applications/utilities/mesh/generation/cv2DMesh/shortEdgeFilter2D.H b/applications/utilities/mesh/generation/cv2DMesh/shortEdgeFilter2D.H
index 0b63de73f7..adf36c29f0 100644
--- a/applications/utilities/mesh/generation/cv2DMesh/shortEdgeFilter2D.H
+++ b/applications/utilities/mesh/generation/cv2DMesh/shortEdgeFilter2D.H
@@ -43,7 +43,6 @@ SourceFiles
 namespace Foam
 {
 
-
 /*---------------------------------------------------------------------------*\
                      Class shortEdgeFilter2D Declaration
 \*---------------------------------------------------------------------------*/
diff --git a/applications/utilities/mesh/generation/cv2DMesh/tolerances.C b/applications/utilities/mesh/generation/cv2DMesh/tolerances.C
deleted file mode 100644
index ba033c907d..0000000000
--- a/applications/utilities/mesh/generation/cv2DMesh/tolerances.C
+++ /dev/null
@@ -1,68 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2007-2010 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 3 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, see <http://www.gnu.org/licenses/>.
-
-\*----------------------------------------------------------------------------*/
-
-#include "tolerances.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::tolerances::tolerances
-(
-    const dictionary& controlDict,
-    const scalar minCellSize,
-    const boundBox& bb
-)
-:
-    span
-    (
-        max(mag(bb.max().x()), mag(bb.min().x()))
-      + max(mag(bb.max().y()), mag(bb.min().y()))
-    ),
-    span2(Foam::sqr(span)),
-
-    minEdgeLen(readScalar(controlDict.lookup("minEdgeLenCoeff"))*minCellSize),
-    minEdgeLen2(Foam::sqr(minEdgeLen)),
-
-    maxNotchLen(readScalar(controlDict.lookup("maxNotchLenCoeff"))*minCellSize),
-    maxNotchLen2(Foam::sqr(maxNotchLen)),
-
-    minNearPointDist
-    (
-        readScalar(controlDict.lookup("minNearPointDistCoeff"))*minCellSize
-    ),
-    minNearPointDist2(Foam::sqr(minNearPointDist)),
-
-    ppDist(readScalar(controlDict.lookup("ppDistCoeff"))*minCellSize),
-    ppDist2(Foam::sqr(ppDist))
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::tolerances::~tolerances()
-{}
-
-
-
-// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cv2DMesh/tolerances.H b/applications/utilities/mesh/generation/cv2DMesh/tolerances.H
deleted file mode 100644
index 7d99c3a61b..0000000000
--- a/applications/utilities/mesh/generation/cv2DMesh/tolerances.H
+++ /dev/null
@@ -1,138 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2011 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 3 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, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::tolerances
-
-Description
-    Tolerances for the CV 2D mesher.
-
-SourceFiles
-    tolerances.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef tolerances_H
-#define tolerances_H
-
-#include "dictionary.H"
-#include "boundBox.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                         Class tolerances Declaration
-\*---------------------------------------------------------------------------*/
-
-class tolerances
-{
-    // Private data
-
-        //- Description of data_
-
-
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        tolerances(const tolerances&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const tolerances&);
-
-
-public:
-
-        //- Maximum cartesian span of the geometry
-        scalar span;
-
-        //- Square of span
-        scalar span2;
-
-        //- Minumum edge-length of the cell size below which protusions
-        //  through the surface are not split
-        scalar minEdgeLen;
-
-        //- Square of minEdgeLen
-        scalar minEdgeLen2;
-
-        //- Maximum notch size below which protusions into the surface are
-        //  not filled
-        scalar maxNotchLen;
-
-        //- Square of maxNotchLen
-        scalar maxNotchLen2;
-
-        //- The minimum distance alowed between a dual-cell vertex
-        // and the surface before a point-pair is introduced
-        scalar minNearPointDist;
-
-        //- Square of minNearPoint
-        scalar minNearPointDist2;
-
-        //- Distance between boundary conforming point-pairs
-        scalar ppDist;
-
-        //- Square of ppDist
-        scalar ppDist2;
-
-
-    // Constructors
-
-        //- Construct null
-        tolerances
-        (
-            const dictionary& controlDict,
-            scalar minCellSize,
-            const boundBox&
-        );
-
-
-    //- Destructor
-    ~tolerances();
-
-
-    // Member Functions
-
-        // Access
-
-        // Check
-
-        // Edit
-
-        // Write
-
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* // 
diff --git a/applications/utilities/mesh/generation/cvMesh/changesLaurence b/applications/utilities/mesh/generation/cvMesh/changesLaurence
new file mode 100644
index 0000000000..5f55eff462
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/changesLaurence
@@ -0,0 +1,5 @@
+wmake/rules/General/CGAL
+
+    -lboost_thread
+    -lboost_thread-mt
+
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index 6228f02e4c..dd58fc5141 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@@ -1958,7 +1958,7 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
         relaxationModel::New
         (
             cvMeshDict.subDict("motionControl"),
-            *this
+            runTime_
         )
     ),
     faceAreaWeightModel_
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/relaxationModel/adaptiveLinear/adaptiveLinear.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/relaxationModel/adaptiveLinear/adaptiveLinear.C
index e5589a92e1..a1967da49d 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/relaxationModel/adaptiveLinear/adaptiveLinear.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/relaxationModel/adaptiveLinear/adaptiveLinear.C
@@ -41,13 +41,13 @@ addToRunTimeSelectionTable(relaxationModel, adaptiveLinear, dictionary);
 adaptiveLinear::adaptiveLinear
 (
     const dictionary& relaxationDict,
-    const conformalVoronoiMesh& cvMesh
+    const Time& runTime
 )
 :
-    relaxationModel(typeName, relaxationDict, cvMesh),
+    relaxationModel(typeName, relaxationDict, runTime),
     relaxationStart_(readScalar(coeffDict().lookup("relaxationStart"))),
     relaxationEnd_(readScalar(coeffDict().lookup("relaxationEnd"))),
-    lastTimeValue_(cvMesh_.time().timeOutputValue()),
+    lastTimeValue_(runTime_.time().timeOutputValue()),
     relaxation_(relaxationStart_)
 {}
 
@@ -56,24 +56,24 @@ adaptiveLinear::adaptiveLinear
 
 scalar adaptiveLinear::relaxation()
 {
-    if (cvMesh_.time().timeOutputValue() > lastTimeValue_)
+    if (runTime_.time().timeOutputValue() > lastTimeValue_)
     {
-        scalar currentRelxation = relaxation_;
+        scalar currentRelaxation = relaxation_;
 
         relaxation_ -=
             (relaxation_ - relaxationEnd_)
            /(
                 (
-                    cvMesh_.time().endTime().value()
-                  - cvMesh_.time().timeOutputValue()
+                    runTime_.time().endTime().value()
+                  - runTime_.time().timeOutputValue()
                 )
-               /(cvMesh_.time().timeOutputValue() - lastTimeValue_)
+               /(runTime_.time().timeOutputValue() - lastTimeValue_)
               + 1
             );
 
-        lastTimeValue_ = cvMesh_.time().timeOutputValue();
+        lastTimeValue_ = runTime_.time().timeOutputValue();
 
-        return currentRelxation;
+        return currentRelaxation;
     }
 
     return relaxation_;
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/relaxationModel/adaptiveLinear/adaptiveLinear.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/relaxationModel/adaptiveLinear/adaptiveLinear.H
index 11a79c39d3..ecb4c32d4d 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/relaxationModel/adaptiveLinear/adaptiveLinear.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/relaxationModel/adaptiveLinear/adaptiveLinear.H
@@ -83,7 +83,7 @@ public:
         adaptiveLinear
         (
             const dictionary& relaxationDict,
-            const conformalVoronoiMesh& cvMesh
+            const Time& runTime
         );
 
 
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/relaxationModel/rampHoldFall/rampHoldFall.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/relaxationModel/rampHoldFall/rampHoldFall.C
index 17ffeab33f..92b42313de 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/relaxationModel/rampHoldFall/rampHoldFall.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/relaxationModel/rampHoldFall/rampHoldFall.C
@@ -41,10 +41,10 @@ addToRunTimeSelectionTable(relaxationModel, rampHoldFall, dictionary);
 rampHoldFall::rampHoldFall
 (
     const dictionary& relaxationDict,
-    const conformalVoronoiMesh& cvMesh
+    const Time& runTime
 )
 :
-    relaxationModel(typeName, relaxationDict, cvMesh),
+    relaxationModel(typeName, relaxationDict, runTime),
     rampStartRelaxation_(readScalar(coeffDict().lookup("rampStartRelaxation"))),
     holdRelaxation_(readScalar(coeffDict().lookup("holdRelaxation"))),
     fallEndRelaxation_(readScalar(coeffDict().lookup("fallEndRelaxation"))),
@@ -62,10 +62,10 @@ rampHoldFall::rampHoldFall
 
 scalar rampHoldFall::relaxation()
 {
-    scalar t = cvMesh_.time().timeOutputValue();
+    scalar t = runTime_.time().timeOutputValue();
 
-    scalar tStart = cvMesh_.time().startTime().value();
-    scalar tEnd = cvMesh_.time().endTime().value();
+    scalar tStart = runTime_.time().startTime().value();
+    scalar tEnd = runTime_.time().endTime().value();
     scalar tSpan = tEnd - tStart;
 
     if (tSpan < VSMALL)
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/relaxationModel/rampHoldFall/rampHoldFall.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/relaxationModel/rampHoldFall/rampHoldFall.H
index 6c65e15880..db814749f1 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/relaxationModel/rampHoldFall/rampHoldFall.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/relaxationModel/rampHoldFall/rampHoldFall.H
@@ -88,7 +88,7 @@ public:
         rampHoldFall
         (
             const dictionary& relaxationDict,
-            const conformalVoronoiMesh& cvMesh
+            const Time& runTime
         );
 
 
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/relaxationModel/relaxationModel/relaxationModel.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/relaxationModel/relaxationModel/relaxationModel.C
index dd7aee8980..91daf73f35 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/relaxationModel/relaxationModel/relaxationModel.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/relaxationModel/relaxationModel/relaxationModel.C
@@ -43,11 +43,11 @@ relaxationModel::relaxationModel
 (
     const word& type,
     const dictionary& relaxationDict,
-    const conformalVoronoiMesh& cvMesh
+    const Time& runTime
 )
 :
     dictionary(relaxationDict),
-    cvMesh_(cvMesh),
+    runTime_(runTime),
     coeffDict_(subDict(type + "Coeffs"))
 {}
 
@@ -57,7 +57,7 @@ relaxationModel::relaxationModel
 autoPtr<relaxationModel> relaxationModel::New
 (
     const dictionary& relaxationDict,
-    const conformalVoronoiMesh& cvMesh
+    const Time& runTime
 )
 {
     word relaxationModelTypeName
@@ -85,7 +85,7 @@ autoPtr<relaxationModel> relaxationModel::New
             << exit(FatalError);
     }
 
-    return autoPtr<relaxationModel>(cstrIter()(relaxationDict, cvMesh));
+    return autoPtr<relaxationModel>(cstrIter()(relaxationDict, runTime));
 }
 
 
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/relaxationModel/relaxationModel/relaxationModel.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/relaxationModel/relaxationModel/relaxationModel.H
index d60390ae9a..5517098f05 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/relaxationModel/relaxationModel/relaxationModel.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/relaxationModel/relaxationModel/relaxationModel.H
@@ -37,7 +37,7 @@ SourceFiles
 #define relaxationModel_H
 
 #include "point.H"
-#include "conformalVoronoiMesh.H"
+#include "Time.H"
 #include "dictionary.H"
 #include "autoPtr.H"
 #include "runTimeSelectionTables.H"
@@ -61,7 +61,7 @@ protected:
     // Protected data
 
         //- Reference to the conformalVoronoiMesh holding this cvControls object
-        const conformalVoronoiMesh& cvMesh_;
+        const Time& runTime_;
 
         //- Method coeffs dictionary
         dictionary coeffDict_;
@@ -92,9 +92,9 @@ public:
             dictionary,
             (
                 const dictionary& relaxationDict,
-                const conformalVoronoiMesh& cvMesh
+                const Time& runTime
             ),
-            (relaxationDict, cvMesh)
+            (relaxationDict, runTime)
         );
 
 
@@ -105,7 +105,7 @@ public:
         (
             const word& type,
             const dictionary& relaxationDict,
-            const conformalVoronoiMesh& cvMesh
+            const Time& runTime
         );
 
 
@@ -115,7 +115,7 @@ public:
         static autoPtr<relaxationModel> New
         (
             const dictionary& relaxationDict,
-            const conformalVoronoiMesh& cvMesh
+            const Time& runTime
         );
 
 
diff --git a/applications/utilities/mesh/generation/extrude/Allwmake b/applications/utilities/mesh/generation/extrude/Allwmake
index 091ec04743..5f09893321 100755
--- a/applications/utilities/mesh/generation/extrude/Allwmake
+++ b/applications/utilities/mesh/generation/extrude/Allwmake
@@ -4,6 +4,8 @@ set -x
 
 wmake libso extrudeModel
 wmake extrudeMesh
+
+wmake libso extrudeToRegionMesh/createShellMesh
 wmake extrudeToRegionMesh
 
 
diff --git a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/Make/files b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/Make/files
index 28ac7de544..22606204b5 100644
--- a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/Make/files
+++ b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/Make/files
@@ -1,4 +1,3 @@
-createShellMesh.C
 extrudeToRegionMesh.C
 
 EXE = $(FOAM_APPBIN)/extrudeToRegionMesh
diff --git a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/Make/options b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/Make/options
index 201089b74f..5732f9761d 100644
--- a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/Make/options
+++ b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/Make/options
@@ -1,11 +1,13 @@
 EXE_INC = \
     -I../extrudeModel/lnInclude \
+    -IcreateShellMesh/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude
 
 EXE_LIBS = \
     -lextrudeModel \
+    -lcreateShellMesh \
     -lfiniteVolume \
     -lmeshTools \
     -ldynamicMesh
diff --git a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh/Make/files b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh/Make/files
new file mode 100644
index 0000000000..4aa716f022
--- /dev/null
+++ b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh/Make/files
@@ -0,0 +1,3 @@
+createShellMesh.C
+
+LIB = $(FOAM_LIBBIN)/libcreateShellMesh
diff --git a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh/Make/options b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh/Make/options
new file mode 100644
index 0000000000..295f0b0bc9
--- /dev/null
+++ b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh/Make/options
@@ -0,0 +1,9 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude -DFULLDEBUG -g -O0
+
+LIB_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -ldynamicMesh
diff --git a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh.C b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh/createShellMesh.C
similarity index 99%
rename from applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh.C
rename to applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh/createShellMesh.C
index 6bfe6ca62e..13aefeac10 100644
--- a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh.C
+++ b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh/createShellMesh.C
@@ -509,7 +509,6 @@ void Foam::createShellMesh::setRefinement
         }
     }
 
-
     // Introduce original points
     // ~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -526,13 +525,12 @@ void Foam::createShellMesh::setRefinement
         );
         pointToPointMap.append(pointI);
 
-        //Pout<< "Added bottom point " << addedPointI
-        //    << " at " << patch_.localPoints()[pointI]
-        //    << "  from point " << pointI
-        //    << endl;
+//        Pout<< "Added bottom point " << pointToPointMap[pointI]
+//            << " at " << patch_.localPoints()[pointI]
+//            << "  from point " << pointI
+//            << endl;
     }
 
-
     // Introduce new points (one for every region)
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -543,6 +541,7 @@ void Foam::createShellMesh::setRefinement
 
         point pt = patch_.localPoints()[pointI];
         point disp = firstLayerDisp[regionI];
+
         for (label layerI = 0; layerI < nLayers; layerI++)
         {
             pt += disp;
@@ -676,7 +675,7 @@ void Foam::createShellMesh::setRefinement
             {
                 FatalErrorIn("createShellMesh::setRefinement(..)")
                     << "external/feature edge:" << edgeI
-                    << " has " << eFaces.size() << " connected extruded faces "
+                    << " has " << eFaces.size() << " connected extruded faces"
                     << " but only " << ePatches.size()
                     << " boundary faces defined." << exit(FatalError);
             }
diff --git a/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh.H b/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh/createShellMesh.H
similarity index 100%
rename from applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh.H
rename to applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh/createShellMesh.H
diff --git a/applications/utilities/mesh/generation/extrude2DMesh/Allwclean b/applications/utilities/mesh/generation/extrude2DMesh/Allwclean
new file mode 100755
index 0000000000..44537ed68b
--- /dev/null
+++ b/applications/utilities/mesh/generation/extrude2DMesh/Allwclean
@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+set -x
+
+wclean libso extrude2DMesh
+wclean
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/applications/utilities/mesh/generation/extrude2DMesh/Allwmake b/applications/utilities/mesh/generation/extrude2DMesh/Allwmake
new file mode 100755
index 0000000000..6f54612ffe
--- /dev/null
+++ b/applications/utilities/mesh/generation/extrude2DMesh/Allwmake
@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+set -x
+
+wmake libso extrude2DMesh
+wmake
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/applications/utilities/mesh/generation/extrude2DMesh/Make/files b/applications/utilities/mesh/generation/extrude2DMesh/Make/files
index 8cc4bd86cd..a30e03b5b7 100644
--- a/applications/utilities/mesh/generation/extrude2DMesh/Make/files
+++ b/applications/utilities/mesh/generation/extrude2DMesh/Make/files
@@ -1,4 +1,3 @@
-extrude2DMesh.C
 extrude2DMeshApp.C
 
 EXE = $(FOAM_APPBIN)/extrude2DMesh
diff --git a/applications/utilities/mesh/generation/extrude2DMesh/Make/options b/applications/utilities/mesh/generation/extrude2DMesh/Make/options
index ef62d4a7c3..8adfe67b36 100644
--- a/applications/utilities/mesh/generation/extrude2DMesh/Make/options
+++ b/applications/utilities/mesh/generation/extrude2DMesh/Make/options
@@ -1,10 +1,12 @@
 EXE_INC = \
-    /* -DFULLDEBUG -g -O0 */ \
-    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
+    -Iextrude2DMesh/lnInclude \
+    -I../extrude/extrudeModel/lnInclude
 
 EXE_LIBS = \
-    -lmeshTools \
+    -lsurfMesh \
     -ldynamicMesh \
-    -lfiniteVolume
+    -lextrude2DMesh \
+    -lextrudeModel
diff --git a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh.C b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh.C
deleted file mode 100644
index 69dc8d66ba..0000000000
--- a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh.C
+++ /dev/null
@@ -1,263 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "extrude2DMesh.H"
-#include "polyMesh.H"
-#include "polyTopoChange.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-defineTypeNameAndDebug(extrude2DMesh, 0);
-
-}
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-// Construct from mesh
-Foam::extrude2DMesh::extrude2DMesh(const polyMesh& mesh)
-:
-    mesh_(mesh)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::extrude2DMesh::setRefinement
-(
-    const direction extrudeDir,
-    const scalar thickness,
-    const label frontPatchI,
-    polyTopoChange& meshMod
-) const
-{
-    for (label cellI = 0; cellI < mesh_.nCells(); cellI++)
-    {
-        meshMod.addCell
-        (
-            -1,     //masterPointID,
-            -1,     //masterEdgeID,
-            -1,     //masterFaceID,
-            cellI,  //masterCellID,
-            mesh_.cellZones().whichZone(cellI)  //zoneID
-        );
-    }
-
-
-    // Generate points
-    // ~~~~~~~~~~~~~~~
-
-    forAll(mesh_.points(), pointI)
-    {
-        meshMod.addPoint
-        (
-            mesh_.points()[pointI],
-            pointI,
-            -1,             // zoneID
-            true            // inCell
-        );
-    }
-
-    //Info<< "Adding offsetted points." << nl << endl;
-    forAll(mesh_.points(), pointI)
-    {
-        point newPoint(mesh_.points()[pointI]);
-        newPoint[extrudeDir] += thickness;
-
-        meshMod.addPoint
-        (
-            newPoint,
-            pointI,
-            -1,             // zoneID
-            true            // inCell
-        );
-    }
-
-
-    // Generate faces
-    // ~~~~~~~~~~~~~~
-
-    const faceList& faces = mesh_.faces();
-    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-
-    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
-    {
-        label zoneID = mesh_.faceZones().whichZone(faceI);
-        bool zoneFlip = false;
-        if (zoneID != -1)
-        {
-            const faceZone& fZone = mesh_.faceZones()[zoneID];
-            zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
-        }
-
-        face newFace(4);
-        const face& f = faces[faceI];
-        newFace[0] = f[0];
-        newFace[1] = f[1];
-        newFace[2] = f[1]+mesh_.nPoints();
-        newFace[3] = f[0]+mesh_.nPoints();
-
-        meshMod.addFace
-        (
-            newFace,
-            mesh_.faceOwner()[faceI],       // own
-            mesh_.faceNeighbour()[faceI],   // nei
-            -1,                             // masterPointID
-            -1,                             // masterEdgeID
-            faceI,                          // masterFaceID
-            false,                          // flipFaceFlux
-            -1,                             // patchID
-            zoneID,                         // zoneID
-            zoneFlip                        // zoneFlip
-        );
-    }
-
-    forAll(patches, patchI)
-    {
-        label startFaceI = patches[patchI].start();
-        label endFaceI = startFaceI + patches[patchI].size();
-
-        for (label faceI = startFaceI; faceI < endFaceI; faceI++)
-        {
-            label zoneID = mesh_.faceZones().whichZone(faceI);
-            bool zoneFlip = false;
-            if (zoneID != -1)
-            {
-                const faceZone& fZone = mesh_.faceZones()[zoneID];
-                zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
-            }
-
-            face newFace(4);
-            const face& f = faces[faceI];
-            newFace[0] = f[0];
-            newFace[1] = f[1];
-            newFace[2] = f[1]+mesh_.nPoints();
-            newFace[3] = f[0]+mesh_.nPoints();
-
-            meshMod.addFace
-            (
-                newFace,
-                mesh_.faceOwner()[faceI],       // own
-                -1,                             // nei
-                -1,                             // masterPointID
-                -1,                             // masterEdgeID
-                faceI,                          // masterFaceID
-                false,                          // flipFaceFlux
-                patchI,                         // patchID
-                zoneID,                         // zoneID
-                zoneFlip                        // zoneFlip
-            );
-        }
-    }
-
-
-    // Generate front and back faces
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    forAll(mesh_.cells(), cellI)
-    {
-        const cell& cFaces = mesh_.cells()[cellI];
-
-        // Make a loop out of faces.
-        const face& f = faces[cFaces[0]];
-
-        face frontFace(cFaces.size());
-        frontFace[0] = f[0];
-
-        label nextPointI = f[1];
-        label nextFaceI = cFaces[0];
-
-        for (label i = 1; i < frontFace.size(); i++)
-        {
-            frontFace[i] = nextPointI;
-
-            // Find face containing pointI
-            forAll(cFaces, cFaceI)
-            {
-                label faceI = cFaces[cFaceI];
-                if (faceI != nextFaceI)
-                {
-                    const face& f = faces[faceI];
-
-                    if (f[0] == nextPointI)
-                    {
-                        nextPointI = f[1];
-                        nextFaceI = faceI;
-                        break;
-                    }
-                    else if (f[1] == nextPointI)
-                    {
-                        nextPointI = f[0];
-                        nextFaceI = faceI;
-                        break;
-                    }
-                }
-            }
-        }
-
-
-        // Add back face.
-        meshMod.addFace
-        (
-            frontFace.reverseFace(),
-            cellI,                          // own
-            -1,                             // nei
-            -1,                             // masterPointID
-            -1,                             // masterEdgeID
-            cFaces[0],                      // masterFaceID
-            false,                          // flipFaceFlux
-            frontPatchI,                    // patchID
-            -1,                             // zoneID
-            false                           // zoneFlip
-        );
-
-        // Offset to create front face.
-        forAll(frontFace, fp)
-        {
-            frontFace[fp] += mesh_.nPoints();
-        }
-        meshMod.addFace
-        (
-            frontFace,
-            cellI,                          // own
-            -1,                             // nei
-            -1,                             // masterPointID
-            -1,                             // masterEdgeID
-            cFaces[0],                      // masterFaceID
-            false,                          // flipFaceFlux
-            frontPatchI,                    // patchID
-            -1,                             // zoneID
-            false                           // zoneFlip
-        );
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/Make/files b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/Make/files
new file mode 100644
index 0000000000..95083a24fb
--- /dev/null
+++ b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/Make/files
@@ -0,0 +1,5 @@
+extrude2DMesh/extrude2DMesh.C
+
+patchToPoly2DMesh/patchToPoly2DMesh.C
+
+LIB = $(FOAM_LIBBIN)/libextrude2DMesh
diff --git a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/Make/options b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/Make/options
new file mode 100644
index 0000000000..8ace7c50da
--- /dev/null
+++ b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/Make/options
@@ -0,0 +1,13 @@
+EXE_INC = \
+    /* -DFULLDEBUG -g -O0 */ \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
+    -I$(FOAM_APP)/utilities/mesh/generation/extrude/extrudeModel/lnInclude
+
+LIB_LIBS = \
+    -lmeshTools \
+    -ldynamicMesh \
+    -lsurfMesh \
+    -lfiniteVolume
diff --git a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/extrude2DMesh/extrude2DMesh.C b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/extrude2DMesh/extrude2DMesh.C
new file mode 100644
index 0000000000..ff1a4aeacc
--- /dev/null
+++ b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/extrude2DMesh/extrude2DMesh.C
@@ -0,0 +1,588 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "extrude2DMesh.H"
+#include "polyMesh.H"
+#include "polyTopoChange.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(extrude2DMesh, 0);
+
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::extrude2DMesh::check2D() const
+{
+    const faceList& faces = mesh_.faces();
+    forAll(faces, faceI)
+    {
+        if (faces[faceI].size() != 2)
+        {
+            FatalErrorIn("extrude2DMesh.C")
+                << "Face " << faceI << " size " << faces[faceI].size()
+                << " is not of size 2 so mesh is not proper two-dimensional."
+                << exit(FatalError);
+        }
+    }
+}
+
+
+//void Foam::extrude2DMesh::findExtrudeDirection()
+//{
+//    scalar minRange = GREAT;
+
+//    for (direction dir = 0; dir < 3; dir++)
+//    {
+//        scalarField cmpts(mesh_.points().component(dir));
+
+//        scalar range = max(cmpts)-min(cmpts);
+
+//        Info<< "Direction:" << dir << " range:" << range << endl;
+
+//        if (range < minRange)
+//        {
+//            minRange = range;
+//            extrudeDir_ = dir;
+//        }
+//    }
+
+//    Info<< "Extruding in direction " << extrudeDir_
+//        << " with thickness " << thickness_ << nl
+//        << endl;
+//}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::extrude2DMesh::extrude2DMesh
+(
+    polyMesh& mesh,
+    const dictionary& dict,
+    const extrudeModel& model
+)
+:
+    mesh_(mesh),
+    dict_(dict),
+    //patchDict_(dict.subDict("patchInfo")),
+    model_(model),
+    modelType_(dict.lookup("extrudeModel")),
+    patchType_(dict.lookup("patchType")),
+    frontPatchI_(-1),
+    backPatchI_(-1)
+{
+    check2D();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::extrude2DMesh::~extrude2DMesh()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::extrude2DMesh::addFrontBackPatches()
+{
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+    frontPatchI_ = patches.findPatchID("front");
+    backPatchI_ = patches.findPatchID("back");
+
+    // Add patch.
+    List<polyPatch*> newPatches(patches.size() + 2);
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        newPatches[patchI] =
+            pp.clone
+            (
+                patches,
+                newPatches.size(),
+                pp.size(),
+                pp.start()
+            ).ptr();
+    }
+
+    if (frontPatchI_ == -1)
+    {
+        frontPatchI_ = patches.size();
+
+        newPatches[frontPatchI_] =
+            polyPatch::New
+            (
+                patchType_,
+                "front",
+                0,
+                mesh_.nFaces(),
+                frontPatchI_,
+                patches
+            ).ptr();
+
+//        newPatches[frontPatchI_] = polyPatch::New
+//        (
+//            "front",
+//            patchDict_,
+//            frontPatchI_,
+//            patches
+//        ).ptr();
+
+        Info<< "Adding patch " << newPatches[frontPatchI_]->name()
+            << " at index " << frontPatchI_
+            << " for front faces." << nl << endl;
+    }
+
+    if (backPatchI_ == -1)
+    {
+        backPatchI_ = patches.size() + 1;
+
+        newPatches[backPatchI_] =
+            polyPatch::New
+            (
+                patchType_,
+                "back",
+                0,
+                mesh_.nFaces(),
+                backPatchI_,
+                patches
+            ).ptr();
+
+//        newPatches[frontPatchI_] = polyPatch::New
+//        (
+//            "back",
+//            patchDict_,
+//            backPatchI_,
+//            patches
+//        ).ptr();
+
+        Info<< "Adding patch " << newPatches[backPatchI_]->name()
+            << " at index " << backPatchI_
+            << " for back faces." << nl << endl;
+    }
+
+    mesh_.removeBoundary();
+    mesh_.addPatches(newPatches);
+}
+
+
+void Foam::extrude2DMesh::setRefinement
+(
+    polyTopoChange& meshMod
+)
+{
+    const label nLayers = model_.nLayers();
+    const pointField& points = mesh_.points();
+    label nFaces = 0;
+
+    for (label layer = 0; layer < nLayers; ++layer)
+    {
+        label offset = layer * mesh_.nCells();
+
+        forAll(mesh_.cells(), cellI)
+        {
+            meshMod.addCell
+            (
+                -1,     //masterPointID,
+                -1,     //masterEdgeID,
+                -1,     //masterFaceID,
+                cellI + offset,  //masterCellID,
+                mesh_.cellZones().whichZone(cellI)  //zoneID
+            );
+        }
+    }
+
+
+    // Generate points
+    // ~~~~~~~~~~~~~~~
+
+    for (label layer = 0; layer <= nLayers; ++layer)
+    {
+        label offset = layer * points.size();
+
+        forAll(points, pointI)
+        {
+            // Don't need the surface normal for either linearDirection or
+            // wedge. Will need to add to be able to use others.
+            point newPoint = model_
+            (
+                points[pointI],
+                vector(),
+                layer
+            );
+
+            meshMod.addPoint
+            (
+                newPoint,
+                pointI + offset,
+                -1,             // zoneID
+                true            // inCell
+            );
+        }
+
+        Pout<< "Added " << points.size() << " points to layer "
+            << layer << endl;
+    }
+
+
+    // Generate faces
+    // ~~~~~~~~~~~~~~
+
+    const faceList& faces = mesh_.faces();
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+    for (label layer = 0; layer < nLayers; ++layer)
+    {
+        label currentLayerOffset = layer * mesh_.nPoints();
+        label nextLayerOffset = currentLayerOffset + mesh_.nPoints();
+
+        for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+        {
+            label zoneID = mesh_.faceZones().whichZone(faceI);
+            bool zoneFlip = false;
+            if (zoneID != -1)
+            {
+                const faceZone& fZone = mesh_.faceZones()[zoneID];
+                zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
+            }
+
+            face newFace(4);
+            const face& f = faces[faceI];
+            newFace[0] = f[0] + currentLayerOffset;
+            newFace[1] = f[1] + currentLayerOffset;
+            newFace[2] = f[1] + nextLayerOffset;
+            newFace[3] = f[0] + nextLayerOffset;
+
+//{
+//    vector n = newFace.normal(pointField(meshMod.points()));
+//    label own = mesh_.faceOwner()[faceI];
+//    const labelList& ownPoints = mesh_.cellPoints()[own];
+//    point ownCc = sum(pointField(mesh_.points(), ownPoints))/ownPoints.size();
+//    label nei = mesh_.faceNeighbour()[faceI];
+//    const labelList& neiPoints = mesh_.cellPoints()[nei];
+//    point neiCc = sum(pointField(mesh_.points(), neiPoints))/neiPoints.size();
+//    vector d = neiCc - ownCc;
+
+//    Pout<< "face:" << faceI << " at:" << f.centre(mesh_.points()) << endl
+//        << "    own:" << own << " at:" << ownCc << endl
+//        << "    nei:" << nei << " at:" << neiCc << endl
+//        << "    sign:" << (n & d) << endl
+//        << endl;
+//}
+
+            label offset = layer * mesh_.nCells();
+
+            meshMod.addFace
+            (
+                newFace,
+                mesh_.faceOwner()[faceI] + offset,       // own
+                mesh_.faceNeighbour()[faceI] + offset,   // nei
+                -1,                             // masterPointID
+                -1,                             // masterEdgeID
+                nFaces++,    // masterFaceID
+                false,                          // flipFaceFlux
+                -1,                             // patchID
+                zoneID,                         // zoneID
+                zoneFlip                        // zoneFlip
+            );
+
+            if (debug)
+            {
+                Info<< newFace << " "
+                    << mesh_.faceOwner()[faceI] + offset << " "
+                    << mesh_.faceNeighbour()[faceI] + offset << " "
+                    << nFaces - 1
+                    << endl;
+            }
+        }
+    }
+
+    forAll(patches, patchI)
+    {
+        for (label layer=0; layer < nLayers; layer++)
+        {
+            label currentLayerOffset = layer*mesh_.nPoints();
+            label nextLayerOffset = currentLayerOffset + mesh_.nPoints();
+
+            label startFaceI = patches[patchI].start();
+            label endFaceI = startFaceI + patches[patchI].size();
+
+            for (label faceI = startFaceI; faceI < endFaceI; faceI++)
+            {
+                label zoneID = mesh_.faceZones().whichZone(faceI);
+                bool zoneFlip = false;
+                if (zoneID != -1)
+                {
+                    const faceZone& fZone = mesh_.faceZones()[zoneID];
+                    zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
+                }
+
+                face newFace(4);
+                const face& f = faces[faceI];
+                newFace[0] = f[0] + currentLayerOffset;
+                newFace[1] = f[1] + currentLayerOffset;
+                newFace[2] = f[1] + nextLayerOffset;
+                newFace[3] = f[0] + nextLayerOffset;
+
+                label offset = layer * mesh_.nCells();
+
+                meshMod.addFace
+                (
+                    newFace,
+                    mesh_.faceOwner()[faceI] + offset,       // own
+                    -1,                                      // nei
+                    -1,                                      // masterPointID
+                    -1,                                      // masterEdgeID
+                    nFaces++,                                // masterFaceID
+                    false,                                   // flipFaceFlux
+                    patchI,                                  // patchID
+                    zoneID,                                  // zoneID
+                    zoneFlip                                 // zoneFlip
+                );
+
+                if (debug)
+                {
+                    Info<< newFace << " "
+                        << mesh_.faceOwner()[faceI] + offset << " "
+                        << nFaces - 1
+                        << endl;
+                }
+            }
+        }
+    }
+
+    // Add extra internal faces that need special treatment for owners and
+    // neighbours.
+    forAll(mesh_.cells(), cellI)
+    {
+        const cell& cFaces = mesh_.cells()[cellI];
+
+        face frontFace(cFaces.size());
+
+        // Make a loop out of faces.
+        label nextFaceI = cFaces[0];
+
+        const face& f = faces[nextFaceI];
+
+        label nextPointI;
+        if (mesh_.faceOwner()[nextFaceI] == cellI)
+        {
+            frontFace[0] = f[0];
+            nextPointI = f[1];
+        }
+        else
+        {
+            frontFace[0] = f[1];
+            nextPointI = f[0];
+        }
+
+
+        for (label i = 1; i < frontFace.size(); i++)
+        {
+            frontFace[i] = nextPointI;
+
+            // Find face containing pointI
+            forAll(cFaces, cFaceI)
+            {
+                label faceI = cFaces[cFaceI];
+                if (faceI != nextFaceI)
+                {
+                    const face& f = faces[faceI];
+
+                    if (f[0] == nextPointI)
+                    {
+                        nextPointI = f[1];
+                        nextFaceI = faceI;
+                        break;
+                    }
+                    else if (f[1] == nextPointI)
+                    {
+                        nextPointI = f[0];
+                        nextFaceI = faceI;
+                        break;
+                    }
+                }
+            }
+        }
+
+        for (label layer = 0; layer < nLayers - 1; ++layer)
+        {
+            // Offset to create front face.
+            forAll(frontFace, fp)
+            {
+                frontFace[fp] += mesh_.nPoints();
+            }
+
+            label offset = layer * mesh_.nCells();
+
+            label nei = -1;
+            if (layer != nLayers - 1)
+            {
+                nei = cellI + offset + mesh_.nCells();
+            }
+
+            meshMod.addFace
+            (
+                frontFace,
+                cellI + offset,                 // own
+                nei,                            // nei
+                -1,                             // masterPointID
+                -1,                             // masterEdgeID
+                nFaces++,                       // masterFaceID
+                false,                          // flipFaceFlux
+                -1,                             // patchID
+                -1,                             // zoneID
+                false                           // zoneFlip
+            );
+
+            if (debug)
+            {
+                Info<< frontFace << " "
+                    << cellI + offset << " "
+                    << nei << " "
+                    << nFaces - 1
+                    << endl;
+            }
+        }
+    }
+
+    // Generate front and back faces
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    forAll(mesh_.cells(), cellI)
+    {
+        const cell& cFaces = mesh_.cells()[cellI];
+
+        face frontFace(cFaces.size());
+
+        // Make a loop out of faces.
+        label nextFaceI = cFaces[0];
+
+        const face& f = faces[nextFaceI];
+
+        label nextPointI;
+        if (mesh_.faceOwner()[nextFaceI] == cellI)
+        {
+            frontFace[0] = f[0];
+            nextPointI = f[1];
+        }
+        else
+        {
+            frontFace[0] = f[1];
+            nextPointI = f[0];
+        }
+
+
+        for (label i = 1; i < frontFace.size(); i++)
+        {
+            frontFace[i] = nextPointI;
+
+            // Find face containing pointI
+            forAll(cFaces, cFaceI)
+            {
+                label faceI = cFaces[cFaceI];
+                if (faceI != nextFaceI)
+                {
+                    const face& f = faces[faceI];
+
+                    if (f[0] == nextPointI)
+                    {
+                        nextPointI = f[1];
+                        nextFaceI = faceI;
+                        break;
+                    }
+                    else if (f[1] == nextPointI)
+                    {
+                        nextPointI = f[0];
+                        nextFaceI = faceI;
+                        break;
+                    }
+                }
+            }
+        }
+
+        // Add back face.
+        meshMod.addFace
+        (
+            frontFace.reverseFace(),
+            cellI,                          // own
+            -1,                             // nei
+            -1,                             // masterPointID
+            -1,                             // masterEdgeID
+            nFaces++,                      // masterFaceID
+            false,                          // flipFaceFlux
+            backPatchI_,                    // patchID
+            -1,                             // zoneID
+            false                           // zoneFlip
+        );
+
+        if (debug)
+        {
+            Info<< nl<<frontFace.reverseFace() << " "
+                << cellI << " "
+                << nFaces - 1
+                << endl;
+        }
+
+        // Offset to create front face.
+        forAll(frontFace, fp)
+        {
+            frontFace[fp] += mesh_.nPoints()* (nLayers);
+        }
+
+        label offset = (nLayers - 1) * mesh_.nCells();
+
+        meshMod.addFace
+        (
+            frontFace,
+            cellI + offset,                 // own
+            -1,                             // nei
+            -1,                             // masterPointID
+            -1,                             // masterEdgeID
+            nFaces++, // masterFaceID
+            false,                          // flipFaceFlux
+            frontPatchI_,                   // patchID
+            -1,                             // zoneID
+            false                           // zoneFlip
+        );
+
+        if (debug)
+        {
+            Info<< frontFace << " "
+                << cellI + offset << " "
+                << nFaces - 1
+                << endl;
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh.H b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/extrude2DMesh/extrude2DMesh.H
similarity index 67%
rename from applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh.H
rename to applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/extrude2DMesh/extrude2DMesh.H
index cdeb787946..0733b5e236 100644
--- a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh.H
+++ b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/extrude2DMesh/extrude2DMesh.H
@@ -25,8 +25,8 @@ Class
     Foam::extrude2DMesh
 
 Description
-    Given 2D mesh insert all the topology changes to extrude. Does not work
-    in parallel
+    Given a 2D mesh insert all the topology changes to extrude. Does not work
+    in parallel.
 
 SourceFiles
     extrude2DMesh.C
@@ -36,9 +36,12 @@ SourceFiles
 #ifndef extrude2DMesh_H
 #define extrude2DMesh_H
 
+#include "typeInfo.H"
 #include "label.H"
 #include "scalar.H"
-#include "typeInfo.H"
+#include "labelList.H"
+#include "dictionary.H"
+#include "extrudeModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -50,6 +53,7 @@ class polyMesh;
 class polyTopoChange;
 class mapPolyMesh;
 class mapDistributePolyMesh;
+class polyBoundaryMesh;
 
 /*---------------------------------------------------------------------------*\
                            Class extrude2DMesh Declaration
@@ -60,8 +64,34 @@ class extrude2DMesh
     // Private data
 
         //- Reference to 2D mesh
-        const polyMesh& mesh_;
+        polyMesh& mesh_;
 
+        const dictionary dict_;
+
+        //const dictionary patchDict_;
+
+        const extrudeModel& model_;
+
+        const word modelType_;
+
+        const word patchType_;
+
+        label frontPatchI_;
+        label backPatchI_;
+
+    // Private Member Functions
+
+        //- Check the mesh is 2D
+        void check2D() const;
+
+        //- Find extrusion direction
+        //void findExtrudeDirection();
+
+        //- Disallow default bitwise copy construct
+        extrude2DMesh(const extrude2DMesh&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const extrude2DMesh&);
 
 public:
 
@@ -70,21 +100,25 @@ public:
 
 
     // Constructors
+        extrude2DMesh
+        (
+            polyMesh&,
+            const dictionary& dict,
+            const extrudeModel& model
+        );
 
-        //- Construct from mesh
-        extrude2DMesh(const polyMesh&);
+
+    //- Destructor
+    ~extrude2DMesh();
 
 
     // Member Functions
 
+        //- Add front and back patches
+        void addFrontBackPatches();
+
         //- Play commands into polyTopoChange to extrude mesh.
-        void setRefinement
-        (
-            const direction extrudeDir,
-            const scalar thickness,
-            const label frontPatchI,
-            polyTopoChange&
-        ) const;
+        void setRefinement(polyTopoChange&);
 
         //- Force recalculation of locally stored data on topological change
         void updateMesh(const mapPolyMesh&)
@@ -93,6 +127,17 @@ public:
         //- Force recalculation of locally stored data for mesh distribution
         void distribute(const mapDistributePolyMesh&)
         {}
+
+        label frontPatchI() const
+        {
+            return frontPatchI_;
+        }
+
+        label backPatchI() const
+        {
+            return backPatchI_;
+        }
+
 };
 
 
diff --git a/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/patchTo2DpolyMesh/patchTo2DpolyMesh.C b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.C
similarity index 92%
rename from applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/patchTo2DpolyMesh/patchTo2DpolyMesh.C
rename to applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.C
index c0121533bb..30c7a0334a 100644
--- a/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/patchTo2DpolyMesh/patchTo2DpolyMesh.C
+++ b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.C
@@ -23,9 +23,8 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "patchTo2DpolyMesh.H"
+#include "patchToPoly2DMesh.H"
 #include "PatchTools.H"
-#include "polyMesh.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -35,7 +34,7 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::patchTo2DpolyMesh::flipFaceOrder()
+void Foam::patchToPoly2DMesh::flipFaceOrder()
 {
     const edgeList& edges = patch_.edges();
     const faceList& localFaces = patch_.localFaces();
@@ -69,7 +68,7 @@ void Foam::patchTo2DpolyMesh::flipFaceOrder()
 }
 
 
-void Foam::patchTo2DpolyMesh::createNeighbours()
+void Foam::patchToPoly2DMesh::createNeighbours()
 {
     const edgeList& edges = patch_.edges();
     const labelListList& edgeFaces = patch_.edgeFaces();
@@ -102,7 +101,7 @@ void Foam::patchTo2DpolyMesh::createNeighbours()
 }
 
 
-Foam::labelList Foam::patchTo2DpolyMesh::internalFaceOrder()
+Foam::labelList Foam::patchToPoly2DMesh::internalFaceOrder()
 {
     const labelListList& cellFaces = patch_.faceEdges();
 
@@ -151,7 +150,7 @@ Foam::labelList Foam::patchTo2DpolyMesh::internalFaceOrder()
 }
 
 
-void Foam::patchTo2DpolyMesh::addPatchFacesToFaces()
+void Foam::patchToPoly2DMesh::addPatchFacesToFaces()
 {
     const labelList& meshPoints = patch_.meshPoints();
 
@@ -175,7 +174,7 @@ void Foam::patchTo2DpolyMesh::addPatchFacesToFaces()
 }
 
 
-void Foam::patchTo2DpolyMesh::addPatchFacesToOwner()
+void Foam::patchToPoly2DMesh::addPatchFacesToOwner()
 {
     const label nInternalEdges = patch_.nInternalEdges();
     const faceList& faces = patch_.faces();
@@ -231,8 +230,8 @@ void Foam::patchTo2DpolyMesh::addPatchFacesToOwner()
             }
             else if
             (
-                e[0] == patch_.edges()[bEdgeI][1]
-             && e[1] == patch_.edges()[bEdgeI][0]
+                e[0] == meshPoints[patch_.edges()[bEdgeI][1]]
+             && e[1] == meshPoints[patch_.edges()[bEdgeI][0]]
             )
             {
                 Info<< "Warning: Wrong orientation." << endl;
@@ -257,7 +256,7 @@ void Foam::patchTo2DpolyMesh::addPatchFacesToOwner()
 }
 
 
-void Foam::patchTo2DpolyMesh::createPolyMeshComponents()
+void Foam::patchToPoly2DMesh::createPolyMeshComponents()
 {
     flipFaceOrder();
 
@@ -282,7 +281,7 @@ void Foam::patchTo2DpolyMesh::createPolyMeshComponents()
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 //- Construct from a primitivePatch
-Foam::patchTo2DpolyMesh::patchTo2DpolyMesh
+Foam::patchToPoly2DMesh::patchToPoly2DMesh
 (
     const MeshedSurface<face>& patch,
     const wordList& patchNames,
@@ -304,13 +303,13 @@ Foam::patchTo2DpolyMesh::patchTo2DpolyMesh
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-Foam::patchTo2DpolyMesh::~patchTo2DpolyMesh()
+Foam::patchToPoly2DMesh::~patchToPoly2DMesh()
 {}
 
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-void Foam::patchTo2DpolyMesh::createMesh()
+void Foam::patchToPoly2DMesh::createMesh()
 {
     createPolyMeshComponents();
 
diff --git a/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/patchTo2DpolyMesh/patchTo2DpolyMesh.H b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.H
similarity index 89%
rename from applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/patchTo2DpolyMesh/patchTo2DpolyMesh.H
rename to applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.H
index 4b831b574e..a07f73f0d4 100644
--- a/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/patchTo2DpolyMesh/patchTo2DpolyMesh.H
+++ b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/patchToPoly2DMesh/patchToPoly2DMesh.H
@@ -22,25 +22,21 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::patchTo2DpolyMesh
+    Foam::patchToPoly2DMesh
 
 Description
     Convert a primitivePatch into a 2D polyMesh.
 
 SourceFiles
-    patchTo2DpolyMeshI.H
-    patchTo2DpolyMesh.C
-    patchTo2DpolyMeshIO.C
+    patchToPoly2DMesh.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef patchTo2DpolyMesh_H
-#define patchTo2DpolyMesh_H
+#ifndef patchToPoly2DMesh_H
+#define patchToPoly2DMesh_H
 
 #include "EdgeMap.H"
-#include "polyMesh.H"
 #include "MeshedSurface.H"
-#include "Time.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -48,10 +44,10 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                     Class patchTo2DpolyMesh Declaration
+                     Class patchToPoly2DMesh Declaration
 \*---------------------------------------------------------------------------*/
 
-class patchTo2DpolyMesh
+class patchToPoly2DMesh
 {
     // Private data
 
@@ -87,10 +83,10 @@ class patchTo2DpolyMesh
     // Private Member Functions
 
         //- Disallow default bitwise copy construct
-        patchTo2DpolyMesh(const patchTo2DpolyMesh&);
+        patchToPoly2DMesh(const patchToPoly2DMesh&);
 
         //- Disallow default bitwise assignment
-        void operator=(const patchTo2DpolyMesh&);
+        void operator=(const patchToPoly2DMesh&);
 
 
 public:
@@ -98,7 +94,7 @@ public:
     // Constructors
 
         //- Construct from a primitivePatch
-        patchTo2DpolyMesh
+        patchToPoly2DMesh
         (
             const MeshedSurface<face>& patch,
             const wordList& patchNames,
@@ -108,7 +104,7 @@ public:
 
 
     //- Destructor
-    ~patchTo2DpolyMesh();
+    ~patchToPoly2DMesh();
 
 
     // Member Functions
diff --git a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMeshApp.C b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMeshApp.C
index 364ee98b79..7cab3da3dd 100644
--- a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMeshApp.C
+++ b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMeshApp.C
@@ -28,13 +28,6 @@ Description
     Takes 2D mesh (all faces 2 points only, no front and back faces) and
     creates a 3D mesh by extruding with specified thickness.
 
-Usage
-
-    - extrude2DMesh thickness
-
-    \param thickness \n
-    Thickness (in metre) of slab.
-
 Note
     Not sure about the walking of the faces to create the front and back faces.
 
@@ -43,155 +36,266 @@ Note
 #include "argList.H"
 #include "Time.H"
 #include "polyMesh.H"
-#include "polyTopoChange.H"
 #include "extrude2DMesh.H"
-#include "emptyPolyPatch.H"
+#include "extrudeModel.H"
+#include "polyTopoChange.H"
+#include "MeshedSurface.H"
+#include "edgeCollapser.H"
+#include "addPatchCellLayer.H"
+#include "patchToPoly2DMesh.H"
 
 using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+enum ExtrudeMode
+{
+    POLYMESH2D,
+    MESHEDSURFACE
+};
+
+namespace Foam
+{
+    template<>
+    const char* NamedEnum<ExtrudeMode, 2>::names[] =
+    {
+        "polyMesh2D",
+        "MeshedSurface"
+    };
+}
+
+static const NamedEnum<ExtrudeMode, 2> ExtrudeModeNames;
+
+
+//pointField moveInitialPoints
+//(
+//    primitiveFacePatch& fMesh,
+//    const extrudeModel& model
+//)
+//{
+//    pointField layer0Points(fMesh.nPoints());
+//    pointField layer1Points(fMesh.nPoints());
+//    pointField displacement(fMesh.nPoints());
+
+//    forAll(layer0Points, pointI)
+//    {
+//        const labelList& meshPoints = fMesh.meshPoints();
+//        label meshPointI = meshPoints[pointI];
+
+//        layer0Points[meshPointI] = model
+//        (
+//            fMesh.points()[meshPointI],
+//            fMesh.pointNormals()[pointI],
+//            0
+//        );
+
+//        layer1Points[meshPointI] = model
+//        (
+//            fMesh.points()[meshPointI],
+//            fMesh.pointNormals()[pointI],
+//            1
+//        );
+
+//        displacement[pointI] =
+//            layer1Points[meshPointI]
+//          - layer0Points[meshPointI];
+//    }
+
+//    fMesh.movePoints(layer0Points);
+
+//    return displacement;
+//}
+
+
 // Main program:
 
 int main(int argc, char *argv[])
 {
-#   include "addOverwriteOption.H"
-    argList::validArgs.append("thickness");
+    #include "addOverwriteOption.H"
 
-#   include "setRootCase.H"
-#   include "createTime.H"
-    runTime.functionObjects().off();
-#   include "createPolyMesh.H"
-    const word oldInstance = mesh.pointsInstance();
+    argList::validArgs.append("surfaceFormat");
 
-    const scalar thickness = args.argRead<scalar>(1);
-    const bool overwrite   = args.optionFound("overwrite");
+    #include "setRootCase.H"
 
+    Info<< "Create time\n" << endl;
 
-    // Check that mesh is 2D
-    // ~~~~~~~~~~~~~~~~~~~~~
+    Time runTimeExtruded
+    (
+        Time::controlDictName,
+        args.rootPath(),
+        args.caseName()
+    );
 
-    const faceList& faces = mesh.faces();
-    forAll(faces, faceI)
+    runTimeExtruded.functionObjects().off();
+
+    const bool overwrite = args.optionFound("overwrite");
+
+    const ExtrudeMode surfaceFormat = ExtrudeModeNames[args[1]];
+
+    Info<< "Extruding from " << ExtrudeModeNames[surfaceFormat] << endl;
+
+    IOdictionary extrude2DMeshDict
+    (
+        IOobject
+        (
+            "extrude2DMeshDict",
+            runTimeExtruded.system(),
+            runTimeExtruded,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false
+        )
+    );
+
+    // Point generator
+    autoPtr<extrudeModel> model(extrudeModel::New(extrude2DMeshDict));
+
+    autoPtr<MeshedSurface<face> > fMesh;
+
+    autoPtr<polyMesh> mesh;
+
+    autoPtr<polyTopoChange> meshMod;
+
+    labelListList extrudeEdgePatches;
+
+    if (surfaceFormat == MESHEDSURFACE)
     {
-        if (faces[faceI].size() != 2)
+        fMesh.set(new MeshedSurface<face>("MeshedSurface.obj"));
+
+        EdgeMap<label> edgeRegionMap;
+        wordList patchNames(1, "default");
+        labelList patchSizes(1, fMesh->nEdges() - fMesh->nInternalEdges());
+
+        const edgeList& edges = fMesh->edges();
+        forAll(edges, edgeI)
         {
-            FatalErrorIn(args.executable())
-                << "Face " << faceI << " size " << faces[faceI].size()
-                << " is not of size 2 so mesh is not proper two-dimensional."
-                << exit(FatalError);
+            if (!fMesh->isInternalEdge(edgeI))
+            {
+                edgeRegionMap.insert(edges[edgeI], 0);
+            }
         }
-    }
 
+        patchToPoly2DMesh poly2DMesh
+        (
+            fMesh(),
+            patchNames,
+            patchSizes,
+            edgeRegionMap
+        );
 
-    // Find extrude direction
-    // ~~~~~~~~~~~~~~~~~~~~~~
+        poly2DMesh.createMesh();
 
-    scalar minRange = GREAT;
-    direction extrudeDir = 4;   //illegal value.
+        mesh.set
+        (
+            new polyMesh
+            (
+                IOobject
+                (
+                    polyMesh::defaultRegion,
+                    runTimeExtruded.constant(),
+                    runTimeExtruded,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                xferMove(poly2DMesh.points()),
+                xferMove(poly2DMesh.faces()),
+                xferMove(poly2DMesh.owner()),
+                xferMove(poly2DMesh.neighbour())
+            )
+        );
 
-    for (direction dir = 0; dir < 3; dir++)
-    {
-        scalarField cmpts(mesh.points().component(dir));
-
-        scalar range = max(cmpts)-min(cmpts);
-
-        Info<< "Direction:" << dir << " range:" << range << endl;
-
-        if (range < minRange)
-        {
-            minRange = range;
-            extrudeDir = dir;
-        }
-    }
-
-    Info<< "Extruding in direction " << extrudeDir
-        << " with thickness " << thickness << nl
-        << endl;
-
-
-
-    const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
-
-    // Add front and back patch
-    // ~~~~~~~~~~~~~~~~~~~~~~~~
-
-    label frontPatchI = patches.findPatchID("frontAndBack");
-
-    if (frontPatchI == -1)
-    {
-        // Add patch.
-        List<polyPatch*> newPatches(patches.size()+1);
+        Info<< "Constructing patches." << endl;
+        List<polyPatch*> patches(poly2DMesh.patchNames().size());
 
         forAll(patches, patchI)
         {
-            const polyPatch& pp = patches[patchI];
-
-            newPatches[patchI] = pp.clone
+            patches[patchI] = new polyPatch
             (
-                patches,
-                newPatches.size(),
-                pp.size(),
-                pp.start()
-            ).ptr();
+                poly2DMesh.patchNames()[patchI],
+                poly2DMesh.patchSizes()[patchI],
+                poly2DMesh.patchStarts()[patchI],
+                patchI,
+                mesh->boundaryMesh()
+            );
         }
 
-        frontPatchI = patches.size();
-
-        newPatches[frontPatchI] = new emptyPolyPatch
-        (
-            "frontAndBack",
-            0,
-            mesh.nFaces(),
-            frontPatchI,
-            patches
-        );
-
-        Info<< "Adding empty patch " << newPatches[frontPatchI]->name()
-            << " at index " << frontPatchI
-            << " for front and back faces." << nl << endl;
-
-        mesh.removeBoundary();
-        mesh.addPatches(newPatches);
+        mesh->addPatches(patches);
+    }
+    else if (surfaceFormat == POLYMESH2D)
+    {
+        mesh.set
+        (
+            new polyMesh
+            (
+                Foam::IOobject
+                (
+                    Foam::polyMesh::defaultRegion,
+                    runTimeExtruded.timeName(),
+                    runTimeExtruded,
+                    Foam::IOobject::MUST_READ
+                )
+            )
+        );
     }
-
-
-
-    // Topo changes container. Initialise with number of patches.
-    polyTopoChange meshMod(mesh.boundaryMesh().size());
 
     // Engine to extrude mesh
-    extrude2DMesh extruder(mesh);
+    extrude2DMesh extruder(mesh(), extrude2DMeshDict, model());
 
-    // Insert changes into meshMod
-    extruder.setRefinement
-    (
-        extrudeDir,
-        thickness,
-        frontPatchI,
-        meshMod
-    );
+    extruder.addFrontBackPatches();
+
+    meshMod.set(new polyTopoChange(mesh().boundaryMesh().size()));
+
+    extruder.setRefinement(meshMod());
 
     // Create a mesh from topo changes.
-    autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh, false);
+    autoPtr<mapPolyMesh> morphMap = meshMod->changeMesh(mesh(), false);
 
-    mesh.updateMesh(morphMap);
+    mesh->updateMesh(morphMap);
 
-    if (!overwrite)
     {
-        runTime++;
-    }
-    else
-    {
-        mesh.setInstance(oldInstance);
+        edgeCollapser collapser(mesh());
+
+        const edgeList& edges = mesh->edges();
+        const pointField& points = mesh->points();
+
+        const boundBox& bb = mesh->bounds();
+        const scalar mergeDim = 1E-4 * bb.minDim();
+
+        forAll(edges, edgeI)
+        {
+            const edge& e = edges[edgeI];
+
+            scalar d = e.mag(points);
+
+            if (d < mergeDim)
+            {
+                Info<< "Merging edge " << e << " since length " << d
+                    << " << " << mergeDim << nl;
+
+                // Collapse edge to e[0]
+                collapser.collapseEdge(edgeI, e[0]);
+            }
+        }
+
+        polyTopoChange meshModCollapse(mesh());
+
+        collapser.setRefinement(meshModCollapse);
+
+        // Create a mesh from topo changes.
+        autoPtr<mapPolyMesh> morphMap = meshModCollapse.changeMesh(mesh(), false);
+
+        mesh->updateMesh(morphMap);
     }
 
+
+    mesh->setInstance(runTimeExtruded.constant());
+
     // Take over refinement levels and write to new time directory.
-    Pout<< "Writing extruded mesh to time " << runTime.timeName() << nl
+    Pout<< "Writing extruded mesh to time = " << runTimeExtruded.timeName() << nl
         << endl;
 
-    mesh.write();
+    mesh().write();
 
     Pout<< "End\n" << endl;
 
diff --git a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMeshDict b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMeshDict
new file mode 100644
index 0000000000..5e53c29c88
--- /dev/null
+++ b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMeshDict
@@ -0,0 +1,46 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+FoamFile
+{
+    version         2.0;
+    format          ascii;
+
+    root            "";
+    case            "";
+    instance        "";
+    local           "";
+
+    class           dictionary;
+    object          extrude2DMeshDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+extrudeModel        linearDirection;
+//extrudeModel        wedge;
+
+patchType           empty;
+//patchType           wedge;
+
+nLayers             1;
+
+expansionRatio      1.0;
+
+linearDirectionCoeffs
+{
+    direction       (0 0 1);
+    thickness       0.1;
+}
+
+wedgeCoeffs
+{
+    axisPt      (0 0 0);
+    axis        (1 0 0);
+    angle       10;
+}
diff --git a/applications/utilities/surface/surfaceConvert/surfaceConvert.C b/applications/utilities/surface/surfaceConvert/surfaceConvert.C
index 504262012c..66af664ad7 100644
--- a/applications/utilities/surface/surfaceConvert/surfaceConvert.C
+++ b/applications/utilities/surface/surfaceConvert/surfaceConvert.C
@@ -49,6 +49,7 @@ Note
 #include "triSurface.H"
 #include "OFstream.H"
 #include "OSspecific.H"
+#include "Time.H"
 
 using namespace Foam;
 
@@ -83,9 +84,25 @@ int main(int argc, char *argv[])
         "factor",
         "geometry scaling factor - default is 1"
     );
+    argList::addOption
+    (
+        "writePrecision",
+        "label",
+        "write to output with the specified precision"
+    );
 
     argList args(argc, argv);
 
+    if (args.optionFound("writePrecision"))
+    {
+        label writePrecision = args.optionRead<label>("writePrecision");
+
+        IOstream::defaultPrecision(writePrecision);
+        Sout.precision(writePrecision);
+
+        Info<< "Output write precision set to " << writePrecision << endl;
+    }
+
     const fileName importName = args[1];
     const fileName exportName = args[2];
 
diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
index 6998c2b757..5f6ee05fb1 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
@@ -95,6 +95,25 @@ void deleteBox
 }
 
 
+bool onLine(const point& p, const linePointRef& line)
+{
+    const point& a = line.start();
+    const point& b = line.end();
+
+    if
+    (
+        ( p.x() < min(a.x(), b.x()) || p.x() > max(a.x(), b.x()) )
+     || ( p.y() < min(a.y(), b.y()) || p.y() > max(a.y(), b.y()) )
+     || ( p.z() < min(a.z(), b.z()) || p.z() > max(a.z(), b.z()) )
+    )
+    {
+        return false;
+    }
+
+    return true;
+}
+
+
 // Deletes all edges inside/outside bounding box from set.
 void deleteEdges
 (
@@ -116,7 +135,9 @@ void deleteEdges
         // If edge does not intersect the plane, delete.
         scalar intersect = cutPlane.lineIntersect(line);
 
-        if (mag(intersect) > line.mag())
+        point featPoint = intersect * (p1 - p0) + p0;
+
+        if (!onLine(featPoint, line))
         {
             edgeStat[edgeI] = surfaceFeatures::NONE;
         }
@@ -304,7 +325,7 @@ int main(int argc, char *argv[])
     (
         "plane",
         "(nx ny nz)(z0 y0 z0)",
-        "used to create feature edges for 2D meshing"
+        "use a plane to create feature edges for 2D meshing"
     );
 
 #   ifdef ENABLE_CURVATURE
diff --git a/wmake/rules/General/CGAL b/wmake/rules/General/CGAL
index 126640dc0c..3586a1de38 100644
--- a/wmake/rules/General/CGAL
+++ b/wmake/rules/General/CGAL
@@ -9,4 +9,5 @@ CGAL_LIBS = \
     -L$(MPFR_ARCH_PATH)/lib \
     -L$(BOOST_ARCH_PATH)/lib \
     -lmpfr \
-    -lboost_thread
+    -I/usr/lib64 \
+    -lboost_thread-mt