From f07935e6c62c40c4eb0c07ade11d254417ce3843 Mon Sep 17 00:00:00 2001
From: laurence <laurence>
Date: Tue, 13 Dec 2011 17:50:22 +0000
Subject: [PATCH] ENH: cvMesh: 3D CV Mesh Generator

    - Remove cvMesh from constructor arguments for faceAreaWeightModel and cvControls
    - Add I.H files for indexedCell and indexedVertex classes
    - Separate feature point functions; put into conformalVoronoiMeshFeaturePoints.C
    - Add feature point specialisations for 2 external and 1 internal edge
    - Add a struct for feature point types
    - Add a writeCellCentres function
---
 .../cvMesh/conformalVoronoiMesh/Make/files    |   1 +
 .../conformalVoronoiMesh.C                    | 784 +----------------
 .../conformalVoronoiMesh.H                    |  31 +-
 .../conformalVoronoiMeshConformToSurface.C    |  11 +-
 ...alVoronoiMeshFeaturePointSpecialisations.C | 573 +++++++------
 .../conformalVoronoiMeshFeaturePoints.C       | 801 ++++++++++++++++++
 .../conformalVoronoiMeshIO.C                  |  30 +
 .../conformalVoronoiMesh/indexedCell.H        | 290 ++-----
 .../conformalVoronoiMesh/indexedCellI.H       | 288 +++++++
 .../conformalVoronoiMesh/indexedVertex.H      | 288 ++-----
 .../conformalVoronoiMesh/indexedVertexI.H     | 321 +++++++
 .../pointFeatureEdgesTypes.H                  |  87 ++
 .../cvControls/cvControls.C                   |   7 +-
 .../cvControls/cvControls.H                   |  10 +-
 .../cvControls/cvControlsI.H                  |   4 +
 .../faceAreaWeightModel/faceAreaWeightModel.C |   9 +-
 .../faceAreaWeightModel/faceAreaWeightModel.H |  15 +-
 .../piecewiseLinearRamp/piecewiseLinearRamp.C |   5 +-
 .../piecewiseLinearRamp/piecewiseLinearRamp.H |   3 +-
 .../utilities/mesh/generation/cvMesh/cvMesh.C |   4 +-
 20 files changed, 2028 insertions(+), 1534 deletions(-)
 create mode 100644 applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C
 create mode 100644 applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCellI.H
 create mode 100644 applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertexI.H
 create mode 100644 applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/pointFeatureEdgesTypes.H

diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/files b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/files
index 29305c3320..217ec9d12e 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/files
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/files
@@ -4,6 +4,7 @@ conformalVoronoiMesh/conformalVoronoiMesh.C
 conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
 conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
 conformalVoronoiMesh/conformalVoronoiMeshIO.C
+conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C
 conformalVoronoiMesh/conformalVoronoiMeshFeaturePointSpecialisations.C
 
 cvControls/cvControls.C
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index dd58fc5141..6ed54e3361 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@@ -507,734 +507,6 @@ void Foam::conformalVoronoiMesh::insertEdgePointGroups
 }
 
 
-void Foam::conformalVoronoiMesh::createEdgePointGroup
-(
-    const extendedFeatureEdgeMesh& feMesh,
-    const pointIndexHit& edHit,
-    DynamicList<Foam::point>& pts,
-    DynamicList<label>& indices,
-    DynamicList<label>& types
-)
-{
-    label edgeI = edHit.index();
-
-    extendedFeatureEdgeMesh::edgeStatus edStatus = feMesh.getEdgeStatus(edgeI);
-
-    switch (edStatus)
-    {
-        case extendedFeatureEdgeMesh::EXTERNAL:
-        {
-            createExternalEdgePointGroup(feMesh, edHit, pts, indices, types);
-            break;
-        }
-        case extendedFeatureEdgeMesh::INTERNAL:
-        {
-            createInternalEdgePointGroup(feMesh, edHit, pts, indices, types);
-            break;
-        }
-        case extendedFeatureEdgeMesh::FLAT:
-        {
-            createFlatEdgePointGroup(feMesh, edHit, pts, indices, types);
-            break;
-        }
-        case extendedFeatureEdgeMesh::OPEN:
-        {
-            createOpenEdgePointGroup(feMesh, edHit, pts, indices, types);
-            break;
-        }
-        case extendedFeatureEdgeMesh::MULTIPLE:
-        {
-            createMultipleEdgePointGroup(feMesh, edHit, pts, indices, types);
-            break;
-        }
-        case extendedFeatureEdgeMesh::NONE:
-        {
-            break;
-        }
-    }
-}
-
-
-void Foam::conformalVoronoiMesh::createExternalEdgePointGroup
-(
-    const extendedFeatureEdgeMesh& feMesh,
-    const pointIndexHit& edHit,
-    DynamicList<Foam::point>& pts,
-    DynamicList<label>& indices,
-    DynamicList<label>& types
-)
-{
-    const Foam::point& edgePt = edHit.hitPoint();
-
-    scalar ppDist = pointPairDistance(edgePt);
-
-    const vectorField& feNormals = feMesh.normals();
-    const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
-
-    // As this is an external edge, there are two normals by definition
-    const vector& nA = feNormals[edNormalIs[0]];
-    const vector& nB = feNormals[edNormalIs[1]];
-
-    if (mag(nA & nB) > 1 - SMALL)
-    {
-        // The normals are nearly parallel, so this is too sharp a feature to
-        // conform to.
-
-        return;
-    }
-
-    // Normalised distance of reference point from edge point
-    vector refVec((nA + nB)/(1 + (nA & nB)));
-
-    if (magSqr(refVec) > sqr(5.0))
-    {
-        // Limit the size of the conformation
-        ppDist *= 5.0/mag(refVec);
-
-        // Pout<< nl << "createExternalEdgePointGroup limit "
-        //     << "edgePt " << edgePt << nl
-        //     << "refVec " << refVec << nl
-        //     << "mag(refVec) " << mag(refVec) << nl
-        //     << "ppDist " << ppDist << nl
-        //     << "nA " << nA << nl
-        //     << "nB " << nB << nl
-        //     << "(nA & nB) " << (nA & nB) << nl
-        //     << endl;
-    }
-
-    // Convex. So refPt will be inside domain and hence a master point
-    Foam::point refPt = edgePt - ppDist*refVec;
-
-    // Result when the points are eventually inserted.
-    // Add number_of_vertices() at insertion of first vertex to all numbers:
-    // pt           index type
-    // refPt        0     1
-    // reflectedA   1     0
-    // reflectedB   2     0
-
-    // Insert the master point pairing the the first slave
-
-    pts.append(refPt);
-    indices.append(0);
-    types.append(1);
-
-    // Insert the slave points by reflecting refPt in both faces.
-    // with each slave refering to the master
-
-    Foam::point reflectedA = refPt + 2*ppDist*nA;
-    pts.append(reflectedA);
-    indices.append(0);
-    types.append(-1);
-
-    Foam::point reflectedB = refPt + 2*ppDist*nB;
-    pts.append(reflectedB);
-    indices.append(0);
-    types.append(-2);
-}
-
-
-void Foam::conformalVoronoiMesh::createInternalEdgePointGroup
-(
-    const extendedFeatureEdgeMesh& feMesh,
-    const pointIndexHit& edHit,
-    DynamicList<Foam::point>& pts,
-    DynamicList<label>& indices,
-    DynamicList<label>& types
-)
-{
-    const Foam::point& edgePt = edHit.hitPoint();
-
-    scalar ppDist = pointPairDistance(edgePt);
-
-    const vectorField& feNormals = feMesh.normals();
-    const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
-
-    // As this is an external edge, there are two normals by definition
-    const vector& nA = feNormals[edNormalIs[0]];
-    const vector& nB = feNormals[edNormalIs[1]];
-
-    if (mag(nA & nB) > 1 - SMALL)
-    {
-        // The normals are nearly parallel, so this is too sharp a feature to
-        // conform to.
-
-        return;
-    }
-
-    // Normalised distance of reference point from edge point
-    vector refVec((nA + nB)/(1 + (nA & nB)));
-
-    if (magSqr(refVec) > sqr(5.0))
-    {
-        // Limit the size of the conformation
-        ppDist *= 5.0/mag(refVec);
-
-        // Pout<< nl << "createInternalEdgePointGroup limit "
-        //     << "edgePt " << edgePt << nl
-        //     << "refVec " << refVec << nl
-        //     << "mag(refVec) " << mag(refVec) << nl
-        //     << "ppDist " << ppDist << nl
-        //     << "nA " << nA << nl
-        //     << "nB " << nB << nl
-        //     << "(nA & nB) " << (nA & nB) << nl
-        //     << endl;
-    }
-
-    // Concave. master and reflected points inside the domain.
-    Foam::point refPt = edgePt - ppDist*refVec;
-
-    // Generate reflected master to be outside.
-    Foam::point reflMasterPt = refPt + 2*(edgePt - refPt);
-
-    // Reflect reflMasterPt in both faces.
-    Foam::point reflectedA = reflMasterPt - 2*ppDist*nA;
-
-    Foam::point reflectedB = reflMasterPt - 2*ppDist*nB;
-
-    scalar totalAngle =
-        radToDeg(constant::mathematical::pi + acos(nA & nB));
-
-    // Number of quadrants the angle should be split into
-    int nQuads = int(totalAngle/cvMeshControls().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 = 2 + nAddPoints;
-
-    // Add number_of_vertices() at insertion of first vertex to all numbers:
-    // Result for nAddPoints 1 when the points are eventually inserted
-    // pt           index type
-    // reflectedA   0     2
-    // reflectedB   1     2
-    // reflMasterPt 2     0
-
-    // Result for nAddPoints 1 when the points are eventually inserted
-    // pt           index type
-    // reflectedA   0     3
-    // reflectedB   1     3
-    // refPt        2     3
-    // reflMasterPt 3     0
-
-    // Result for nAddPoints 2 when the points are eventually inserted
-    // pt           index type
-    // reflectedA   0     4
-    // reflectedB   1     4
-    // reflectedAa  2     4
-    // reflectedBb  3     4
-    // reflMasterPt 4     0
-
-    // Master A is inside.
-    pts.append(reflectedA);
-    indices.append(0);
-    types.append(reflectedMaster--);
-
-    // Master B is inside.
-    pts.append(reflectedB);
-    indices.append(0);
-    types.append(reflectedMaster--);
-
-    if (nAddPoints == 1)
-    {
-        // One additinal point is the reflection of the slave point,
-        // i.e. the original reference point
-        pts.append(refPt);
-        indices.append(0);
-        types.append(reflectedMaster--);
-    }
-    else if (nAddPoints == 2)
-    {
-        Foam::point reflectedAa = refPt + ppDist*nB;
-        pts.append(reflectedAa);
-        indices.append(0);
-        types.append(reflectedMaster--);
-
-        Foam::point reflectedBb = refPt + ppDist*nA;
-        pts.append(reflectedBb);
-        indices.append(0);
-        types.append(reflectedMaster--);
-    }
-
-    // Slave is outside.
-    pts.append(reflMasterPt);
-    indices.append(0);
-    types.append(-(nAddPoints + 2));
-}
-
-
-void Foam::conformalVoronoiMesh::createFlatEdgePointGroup
-(
-    const extendedFeatureEdgeMesh& feMesh,
-    const pointIndexHit& edHit,
-    DynamicList<Foam::point>& pts,
-    DynamicList<label>& indices,
-    DynamicList<label>& types
-)
-{
-    const Foam::point& edgePt = edHit.hitPoint();
-
-    scalar ppDist = pointPairDistance(edgePt);
-
-    const vectorField& feNormals = feMesh.normals();
-    const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
-
-    // As this is a flat edge, there are two normals by definition
-    const vector& nA = feNormals[edNormalIs[0]];
-    const vector& nB = feNormals[edNormalIs[1]];
-
-    // Average normal to remove any bias to one side, although as this
-    // is a flat edge, the normals should be essentially the same
-    vector n = 0.5*(nA + nB);
-
-    // Direction along the surface to the control point, sense of edge
-    // direction not important, as +s and -s can be used because this
-    // is a flat edge
-    vector s = 2.0*ppDist*(feMesh.edgeDirections()[edHit.index()] ^ n);
-
-    createPointPair(ppDist, edgePt + s, n, pts, indices, types);
-    createPointPair(ppDist, edgePt - s, n, pts, indices, types);
-}
-
-
-void Foam::conformalVoronoiMesh::createOpenEdgePointGroup
-(
-    const extendedFeatureEdgeMesh& feMesh,
-    const pointIndexHit& edHit,
-    DynamicList<Foam::point>& pts,
-    DynamicList<label>& indices,
-    DynamicList<label>& types
-)
-{
-    Info<< "NOT INSERTING OPEN EDGE POINT GROUP, NOT IMPLEMENTED" << endl;
-}
-
-
-void Foam::conformalVoronoiMesh::createMultipleEdgePointGroup
-(
-    const extendedFeatureEdgeMesh& feMesh,
-    const pointIndexHit& edHit,
-    DynamicList<Foam::point>& pts,
-    DynamicList<label>& indices,
-    DynamicList<label>& types
-)
-{
-    Info<< "NOT INSERTING MULTIPLE EDGE POINT GROUP, NOT IMPLEMENTED"
-        << endl;
-}
-
-
-void Foam::conformalVoronoiMesh::insertFeaturePoints()
-{
-    Info<< nl << "Conforming to feature points" << endl;
-
-    DynamicList<Foam::point> pts;
-    DynamicList<label> indices;
-    DynamicList<label> types;
-
-    label preFeaturePointSize = number_of_vertices();
-
-    createConvexFeaturePoints(pts, indices, types);
-
-    createConcaveFeaturePoints(pts, indices, types);
-
-    createMixedFeaturePoints(pts, indices, types);
-
-    // Insert the created points, distributing to the appropriate processor
-    insertPoints(pts, indices, types, true);
-
-    if (cvMeshControls().objOutput())
-    {
-        writePoints("featureVertices.obj", pts);
-    }
-
-    label nFeatureVertices = number_of_vertices() - preFeaturePointSize;
-
-    if (Pstream::parRun())
-    {
-        reduce(nFeatureVertices, sumOp<label>());
-    }
-
-    if (nFeatureVertices > 0)
-    {
-        Info<< "    Inserted " << nFeatureVertices
-            << " feature vertices" << endl;
-    }
-
-    featureVertices_.clear();
-
-    forAll(pts, pI)
-    {
-        featureVertices_.push_back
-        (
-            Vb(toPoint(pts[pI]), indices[pI], types[pI])
-        );
-    }
-
-    constructFeaturePointLocations();
-}
-
-
-void Foam::conformalVoronoiMesh::createConvexFeaturePoints
-(
-    DynamicList<Foam::point>& pts,
-    DynamicList<label>& indices,
-    DynamicList<label>& types
-)
-{
-    const PtrList<extendedFeatureEdgeMesh>& feMeshes
-    (
-        geometryToConformTo_.features()
-    );
-
-    forAll(feMeshes, i)
-    {
-        const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
-
-        for
-        (
-            label ptI = feMesh.convexStart();
-            ptI < feMesh.concaveStart();
-            ptI++
-        )
-        {
-            const Foam::point& featPt = feMesh.points()[ptI];
-
-            if (!positionOnThisProc(featPt))
-            {
-                continue;
-            }
-
-            vectorField featPtNormals = feMesh.featurePointNormals(ptI);
-
-            vector cornerNormal = sum(featPtNormals);
-            cornerNormal /= mag(cornerNormal);
-
-            Foam::point internalPt =
-                featPt - pointPairDistance(featPt)*cornerNormal;
-
-            // Result when the points are eventually inserted (example n = 4)
-            // Add number_of_vertices() at insertion of first vertex to all
-            // numbers:
-            // pt           index type
-            // internalPt   0     1
-            // externalPt0  1     0
-            // externalPt1  2     0
-            // externalPt2  3     0
-            // externalPt3  4     0
-
-            pts.append(internalPt);
-            indices.append(0);
-            types.append(1);
-
-            label internalPtIndex = -1;
-
-            forAll(featPtNormals, nI)
-            {
-                const vector& n = featPtNormals[nI];
-
-                plane planeN = plane(featPt, n);
-
-                Foam::point externalPt =
-                    internalPt + 2.0 * planeN.distance(internalPt) * n;
-
-                pts.append(externalPt);
-                indices.append(0);
-                types.append(internalPtIndex--);
-            }
-        }
-    }
-}
-
-
-void Foam::conformalVoronoiMesh::createConcaveFeaturePoints
-(
-    DynamicList<Foam::point>& pts,
-    DynamicList<label>& indices,
-    DynamicList<label>& types
-)
-{
-    const PtrList<extendedFeatureEdgeMesh>& feMeshes
-    (
-        geometryToConformTo_.features()
-    );
-
-    forAll(feMeshes, i)
-    {
-        const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
-
-        for
-        (
-            label ptI = feMesh.concaveStart();
-            ptI < feMesh.mixedStart();
-            ptI++
-        )
-        {
-            const Foam::point& featPt = feMesh.points()[ptI];
-
-            if (!positionOnThisProc(featPt))
-            {
-                continue;
-            }
-
-            vectorField featPtNormals = feMesh.featurePointNormals(ptI);
-
-            vector cornerNormal = sum(featPtNormals);
-            cornerNormal /= mag(cornerNormal);
-
-            Foam::point externalPt =
-                featPt + pointPairDistance(featPt)*cornerNormal;
-
-            label externalPtIndex = featPtNormals.size();
-
-            // Result when the points are eventually inserted (example n = 5)
-            // Add number_of_vertices() at insertion of first vertex to all
-            // numbers:
-            // pt           index type
-            // internalPt0  0     5
-            // internalPt1  1     5
-            // internalPt2  2     5
-            // internalPt3  3     5
-            // internalPt4  4     5
-            // externalPt   5     4
-
-            forAll(featPtNormals, nI)
-            {
-                const vector& n = featPtNormals[nI];
-
-                plane planeN = plane(featPt, n);
-
-                Foam::point internalPt =
-                    externalPt - 2.0 * planeN.distance(externalPt) * n;
-
-                pts.append(internalPt);
-                indices.append(0);
-                types.append(externalPtIndex--);
-            }
-
-            pts.append(externalPt);
-            indices.append(0);
-            types.append(-1);
-        }
-    }
-}
-
-
-void Foam::conformalVoronoiMesh::createMixedFeaturePoints
-(
-    DynamicList<Foam::point>& pts,
-    DynamicList<label>& indices,
-    DynamicList<label>& types
-)
-{
-    if (cvMeshControls().mixedFeaturePointPPDistanceCoeff() < 0)
-    {
-        Info<<nl << "Skipping specialised handling for mixed feature points"
-            << endl;
-        return;
-    }
-
-    const PtrList<extendedFeatureEdgeMesh>& feMeshes
-    (
-        geometryToConformTo_.features()
-    );
-
-    forAll(feMeshes, i)
-    {
-        const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
-
-        for
-        (
-            label ptI = feMesh.mixedStart();
-            ptI < feMesh.nonFeatureStart();
-            ptI++
-        )
-        {
-            if
-            (
-                !createSpecialisedFeaturePoint(feMesh, ptI, pts, indices, types)
-            )
-            {
-                // Specialisations available for some mixed feature points.  For
-                // non-specialised feature points, inserting mixed internal and
-                // external edge groups at feature point.
-
-                const labelList& pEds(feMesh.pointEdges()[ptI]);
-
-                // Skipping unsupported mixed feature point types
-
-                bool skipEdge = false;
-
-                forAll(pEds, e)
-                {
-                    label edgeI = pEds[e];
-
-                    extendedFeatureEdgeMesh::edgeStatus edStatus =
-                        feMesh.getEdgeStatus(edgeI);
-
-                    if
-                    (
-                        edStatus == extendedFeatureEdgeMesh::OPEN
-                     || edStatus == extendedFeatureEdgeMesh::MULTIPLE
-                    )
-                    {
-                        Info<< "Edge type " << edStatus
-                            << " found for mixed feature point " << ptI
-                            << ". Not supported."
-                            << endl;
-
-                        skipEdge = true;
-                    }
-                }
-
-                if (skipEdge)
-                {
-                    Info<< "Skipping point " << ptI << nl << endl;
-
-                    continue;
-                }
-
-                const Foam::point& pt(feMesh.points()[ptI]);
-
-                if (!positionOnThisProc(pt))
-                {
-                    continue;
-                }
-
-                scalar edgeGroupDistance = mixedFeaturePointDistance(pt);
-
-                forAll(pEds, e)
-                {
-                    label edgeI = pEds[e];
-
-                    Foam::point edgePt =
-                        pt + edgeGroupDistance*feMesh.edgeDirection(edgeI, ptI);
-
-                    pointIndexHit edgeHit(true, edgePt, edgeI);
-
-                    createEdgePointGroup(feMesh, edgeHit, pts, indices, types);
-                }
-            }
-        }
-    }
-}
-
-
-void Foam::conformalVoronoiMesh::constructFeaturePointLocations()
-{
-    DynamicList<Foam::point> ftPtLocs;
-
-    const PtrList<extendedFeatureEdgeMesh>& feMeshes
-    (
-        geometryToConformTo_.features()
-    );
-
-    forAll(feMeshes, i)
-    {
-        const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
-
-
-        if (cvMeshControls().mixedFeaturePointPPDistanceCoeff() < 0)
-        {
-            // Ignoring mixed feature points
-            for
-            (
-                label ptI = feMesh.convexStart();
-                ptI < feMesh.mixedStart();
-                ptI++
-            )
-            {
-                ftPtLocs.append(feMesh.points()[ptI]);
-            }
-        }
-        else
-        {
-            for
-            (
-                label ptI = feMesh.convexStart();
-                ptI < feMesh.nonFeatureStart();
-                ptI++
-            )
-            {
-                ftPtLocs.append(feMesh.points()[ptI]);
-            }
-        }
-    }
-
-    featurePointLocations_.transfer(ftPtLocs);
-}
-
-
-void Foam::conformalVoronoiMesh::reinsertFeaturePoints(bool distribute)
-{
-    Info<< nl << "Reinserting stored feature points" << endl;
-
-    label preReinsertionSize(number_of_vertices());
-
-    if (distribute)
-    {
-        DynamicList<Foam::point> pointsToInsert;
-        DynamicList<label> indices;
-        DynamicList<label> types;
-
-        for
-        (
-            std::list<Vb>::iterator vit=featureVertices_.begin();
-            vit != featureVertices_.end();
-            ++vit
-        )
-        {
-            pointsToInsert.append(topoint(vit->point()));
-            indices.append(vit->index());
-            types.append(vit->type());
-        }
-
-        // Insert distributed points
-        insertPoints(pointsToInsert, indices, types, true);
-
-        // Save points in new distribution
-        featureVertices_.clear();
-
-        forAll(pointsToInsert, pI)
-        {
-            featureVertices_.push_back
-            (
-                Vb(toPoint(pointsToInsert[pI]), indices[pI], types[pI])
-            );
-        }
-    }
-    else
-    {
-        for
-        (
-            std::list<Vb>::iterator vit=featureVertices_.begin();
-            vit != featureVertices_.end();
-            ++vit
-        )
-        {
-            // Assuming that all of the reinsertions are pair points, and that
-            // the index and type are relative, i.e. index 0 and type relative
-            // to it.
-            insertPoint
-            (
-                vit->point(),
-                vit->index() + number_of_vertices(),
-                vit->type() + number_of_vertices()
-            );
-        }
-    }
-
-    Info<< "    Reinserted "
-        << returnReduce
-        (
-            label(number_of_vertices()) - preReinsertionSize,
-            sumOp<label>()
-        )
-        << " vertices" << endl;
-}
-
-
 const Foam::indexedOctree<Foam::treeDataPoint>&
 Foam::conformalVoronoiMesh::featurePointTree() const
 {
@@ -1277,33 +549,17 @@ bool Foam::conformalVoronoiMesh::nearFeaturePt(const Foam::point& pt) const
 }
 
 
-void Foam::conformalVoronoiMesh::reset()
+void Foam::conformalVoronoiMesh::reset(const bool distribute)
 {
     this->clear();
 
-    insertBoundingPoints();
-}
+    reinsertBoundingPoints();
 
+    // Reinsert feature points, distributing them as necessary.
+    reinsertFeaturePoints(distribute);
+    //insertFeaturePoints();
 
-void Foam::conformalVoronoiMesh::insertBoundingPoints()
-{
-    pointField farPts = geometryToConformTo_.globalBounds().points();
-
-    // Shift corners of bounds relative to origin
-    farPts -= geometryToConformTo_.globalBounds().midpoint();
-
-    // Scale the box up
-    farPts *= 10.0;
-
-    // Shift corners of bounds back to be relative to midpoint
-    farPts += geometryToConformTo_.globalBounds().midpoint();
-
-    limitBounds_ = treeBoundBox(farPts);
-
-    forAll(farPts, fPI)
-    {
-        insertPoint(farPts[fPI], Vb::vtFar);
-    }
+    startOfInternalPoints_ = number_of_vertices();
 }
 
 
@@ -1466,13 +722,8 @@ bool Foam::conformalVoronoiMesh::distributeBackground()
             cellVertexTypes
         );
 
-        // Remove the entire tessellation
-        reset();
-
-        // Reinsert feature points, distributing them as necessary.
-        reinsertFeaturePoints(true);
-
-        startOfInternalPoints_ = number_of_vertices();
+        // Reset the entire tessellation
+        reset(true);
 
         timeCheck("Distribution performed");
 
@@ -1909,7 +1160,7 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
     Delaunay(),
     runTime_(runTime),
     rndGen_(64293*Pstream::myProcNo()),
-    cvMeshControls_(*this, cvMeshDict),
+    cvMeshControls_(cvMeshDict),
     allGeometry_
     (
         IOobject
@@ -1965,8 +1216,7 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
     (
         faceAreaWeightModel::New
         (
-            cvMeshDict.subDict("motionControl"),
-            *this
+            cvMeshDict.subDict("motionControl")
         )
     ),
     decomposition_()
@@ -2003,6 +1253,8 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
 
     buildSurfaceConformation(rmCoarse);
 
+    is_valid(true);
+
     // The introduction of the surface conformation may have distorted the
     // balance of vertices, distribute if necessary.
     if (distributeBackground())
@@ -2386,9 +1638,10 @@ void Foam::conformalVoronoiMesh::move()
     // Remove the entire tessellation
     reset();
 
-    reinsertFeaturePoints();
-
-    startOfInternalPoints_ = number_of_vertices();
+    if (cvMeshControls().objOutput() && runTime_.outputTime())
+    {
+        writePoints("featurePoints_" + runTime_.timeName() + ".obj", false);
+    }
 
     timeCheck("Displacement calculated");
 
@@ -2400,6 +1653,11 @@ void Foam::conformalVoronoiMesh::move()
 
     conformToSurface();
 
+    if (cvMeshControls().objOutput() && runTime_.outputTime())
+    {
+        writePoints("points_" + runTime_.timeName() + ".obj", false);
+    }
+
     timeCheck("After conformToSurface");
 
     updateSizesAndAlignments(pointsToInsert);
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
index 53b56da7b2..927827a675 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
@@ -35,6 +35,10 @@ SourceFiles
     conformalVoronoiMeshI.H
     conformalVoronoiMesh.C
     conformalVoronoiMeshIO.C
+    conformalVoronoiMeshConformToSurface.C
+    conformalVoronoiMeshFeaturePoints.C
+    conformalVoronoiMeshFeaturePointSpecialisations.C
+    conformalVoronoiMeshCalcDualMesh.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -73,6 +77,7 @@ SourceFiles
 #include "processorPolyPatch.H"
 #include "zeroGradientFvPatchFields.H"
 #include "globalIndex.H"
+#include "pointFeatureEdgesTypes.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -420,12 +425,24 @@ private:
             DynamicList<label>& types
         );
 
+        //- Fill the pointFeatureEdgesType struct with the types of feature
+        //  edges that are attached to the point.
+        List<extendedFeatureEdgeMesh::edgeStatus> calcPointFeatureEdgesTypes
+        (
+            const extendedFeatureEdgeMesh& feMesh,
+            const labelList& pEds,
+            pointFeatureEdgesTypes& pFEdgeTypes
+        );
+
         //- Create feature point groups if a specialisation exists for the
         //  structure
         bool createSpecialisedFeaturePoint
         (
             const extendedFeatureEdgeMesh& feMesh,
-            label ptI,
+            const labelList& pEds,
+            const pointFeatureEdgesTypes& pFEdgeTypes,
+            const List<extendedFeatureEdgeMesh::edgeStatus>& allEdStat,
+            const label ptI,
             DynamicList<Foam::point>& pts,
             DynamicList<label>& indices,
             DynamicList<label>& types
@@ -443,13 +460,18 @@ private:
         //- Check if a location is in exclusion range around a feature point
         bool nearFeaturePt(const Foam::point& pt) const;
 
-        //- Clear the entire tesselation and insert bounding points
-        void reset();
+        //- Clear the entire tesselation
+        //  Reinsert bounding points, feature points and recalculate
+        //  startOfInternalPoints_
+        void reset(const bool distribute = false);
 
         //- Insert far points in a large bounding box to avoid dual edges
         //  spanning huge distances
         void insertBoundingPoints();
 
+        //- Reinsert the bounding points
+        void reinsertBoundingPoints();
+
         //- Insert the initial points into the triangulation, based on the
         //  initialPointsMethod
         void insertInitialPoints();
@@ -1045,6 +1067,9 @@ public:
             //  actual cell size (cbrt(actual cell volume)).
             void writeCellSizes(const fvMesh& mesh) const;
 
+            //- Calculate and write the cell centres.
+            void writeCellCentres(const fvMesh& mesh) const;
+
             //- Find the cellSet of the boundary cells which have points that
             //  protrude out of the surface beyond a tolerance.
             void findRemainingProtrusionSet(const fvMesh& mesh) const;
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
index e65a5beefa..d488626825 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
@@ -203,10 +203,13 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
 
                 if (surfHit.hit())
                 {
-                    vit->setNearBoundary();
-
                     if (dualCellSurfaceAnyIntersection(vit))
                     {
+                        // This used to be just before this if statement.
+                        // Moved because a point is only near the boundary if
+                        // the dual cell intersects the surface.
+                        vit->setNearBoundary();
+
                         // meshTools::writeOBJ(Pout, vert);
                         // meshTools::writeOBJ(Pout, surfHit.hitPoint());
                         // Pout<< "l cr0 cr1" << endl;
@@ -497,8 +500,8 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAnyIntersection
 
         if (Pstream::parRun())
         {
-            Foam::point a = dE0;
-            Foam::point b = dE1;
+            Foam::point& a = dE0;
+            Foam::point& b = dE1;
 
             bool inProc = clipLineToProc(a, b);
 
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePointSpecialisations.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePointSpecialisations.C
index 95040bf865..64e46131d2 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePointSpecialisations.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePointSpecialisations.C
@@ -27,30 +27,14 @@ License
 
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
-bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
+Foam::List<Foam::extendedFeatureEdgeMesh::edgeStatus>
+Foam::conformalVoronoiMesh::calcPointFeatureEdgesTypes
 (
     const extendedFeatureEdgeMesh& feMesh,
-    label ptI,
-    DynamicList<Foam::point>& pts,
-    DynamicList<label>& indices,
-    DynamicList<label>& types
+    const labelList& pEds,
+    pointFeatureEdgesTypes& pFEdgeTypes
 )
 {
-    const labelList& pEds(feMesh.pointEdges()[ptI]);
-
-    if (pEds.size() != 3)
-    {
-        // Only three edge specialisations available
-
-        return false;
-    }
-
-    label nExternal = 0;
-    label nInternal = 0;
-    label nFlat = 0;
-    label nOpen = 0;
-    label nMultiple = 0;
-
     List<extendedFeatureEdgeMesh::edgeStatus> allEdStat(pEds.size());
 
     forAll(pEds, i)
@@ -65,293 +49,322 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
         {
             case extendedFeatureEdgeMesh::EXTERNAL:
             {
-                nExternal++;
+                pFEdgeTypes.nExternal++;
                 break;
             }
             case extendedFeatureEdgeMesh::INTERNAL:
             {
-                nInternal++;
+                pFEdgeTypes.nInternal++;
                 break;
             }
             case extendedFeatureEdgeMesh::FLAT:
             {
-                nFlat++;
+                pFEdgeTypes.nFlat++;
                 break;
             }
             case extendedFeatureEdgeMesh::OPEN:
             {
-                nOpen++;
+                pFEdgeTypes.nOpen++;
                 break;
             }
             case extendedFeatureEdgeMesh::MULTIPLE:
             {
-                nMultiple++;
+                pFEdgeTypes.nMultiple++;
                 break;
             }
             case extendedFeatureEdgeMesh::NONE:
             {
+                pFEdgeTypes.nNonFeature++;
                 break;
             }
         }
     }
 
-    if (nExternal == 2 && nInternal == 1)
+    if (debug)
     {
-        // if (!positionOnThisProc(pt))
-        // {
-        //     continue;
-        // }
-
-        // Info<< "nExternal == 2 && nInternal == 1" << endl;
-
-        // const Foam::point& featPt = feMesh.points()[ptI];
-
-        // scalar ppDist = pointPairDistance(featPt);
-
-        // const vectorField& normals = feMesh.normals();
-
-        // const labelListList& edgeNormals = feMesh.edgeNormals();
-
-        // label concaveEdgeI = pEds
-        // [
-        //     findIndex(allEdStat, extendedFeatureEdgeMesh::INTERNAL)
-        // ];
-
-        // // // Find which planes are joined to the concave edge
-
-        // // List<label> concaveEdgePlanes(2,label(-1));
-
-        // // label concaveEdgeI = concaveEdges[0];
-
-        // // // Pick up the two faces adjacent to the concave feature edge
-        // // const labelList& eFaces = qSurf_.edgeFaces()[concaveEdgeI];
-
-        // // label faceA = eFaces[0];
-
-        // // vector nA = qSurf_.faceNormals()[faceA];
-
-        // // scalar maxNormalDotProduct = -SMALL;
-
-        // // forAll(uniquePlaneNormals, uPN)
-        // // {
-        // //     scalar normalDotProduct = nA & uniquePlaneNormals[uPN];
-
-        // //     if (normalDotProduct > maxNormalDotProduct)
-        // //     {
-        // //         maxNormalDotProduct = normalDotProduct;
-
-        // //         concaveEdgePlanes[0] = uPN;
-        // //     }
-        // // }
-
-        // // label faceB = eFaces[1];
-        // // vector nB = qSurf_.faceNormals()[faceB];
-
-        // // maxNormalDotProduct = -SMALL;
-
-        // // forAll(uniquePlaneNormals, uPN)
-        // // {
-        // //     scalar normalDotProduct = nB & uniquePlaneNormals[uPN];
-
-        // //     if (normalDotProduct > maxNormalDotProduct)
-        // //     {
-        // //         maxNormalDotProduct = normalDotProduct;
-
-        // //         concaveEdgePlanes[1] = uPN;
-        // //     }
-        // // }
-
-        // // const vector& concaveEdgePlaneANormal =
-        // // uniquePlaneNormals[concaveEdgePlanes[0]];
-
-        // // const vector& concaveEdgePlaneBNormal =
-        // // uniquePlaneNormals[concaveEdgePlanes[1]];
-
-        // // label convexEdgesPlaneI;
-
-        // // if (findIndex(concaveEdgePlanes, 0) == -1)
-        // // {
-        // //     convexEdgesPlaneI = 0;
-        // // }
-        // // else if (findIndex(concaveEdgePlanes, 1) == -1)
-        // // {
-        // //     convexEdgesPlaneI = 1;
-        // // }
-        // // else
-        // // {
-        // //     convexEdgesPlaneI = 2;
-        // // }
-
-        // // const vector& convexEdgesPlaneNormal =
-        // // uniquePlaneNormals[convexEdgesPlaneI];
-
-        // // const edge& concaveEdge = edges[concaveEdgeI];
-
-        // // // Check direction of edge, if the feature point is at the end()
-        // // // the reverse direction.
-
-        // // scalar edgeDirection = 1.0;
-
-        // // if (ptI == concaveEdge.end())
-        // // {
-        // //     edgeDirection *= -1.0;
-        // // }
-
-
-        // const vector& concaveEdgePlaneANormal =
-        //     normals[edgeNormals[concaveEdgeI][0]];
-
-        // const vector& concaveEdgePlaneBNormal =
-        //     normals[edgeNormals[concaveEdgeI][1]];
-
-        // // Intersect planes parallel to the concave edge planes offset
-        // // by ppDist and the plane defined by featPt and the edge
-        // // vector.
-        // plane planeA
-        // (
-        //     featPt + ppDist*concaveEdgePlaneANormal,
-        //     concaveEdgePlaneANormal
-        // );
-
-        // plane planeB
-        // (
-        //     featPt + ppDist*concaveEdgePlaneBNormal,
-        //     concaveEdgePlaneBNormal
-        // );
-
-        // const vector& concaveEdgeDir = feMesh.edgeDirection
-        // (
-        //     concaveEdgeI,
-        //     ptI
-        // );
-
-        // Foam::point concaveEdgeLocalFeatPt = featPt + ppDist*concaveEdgeDir;
-
-        // // Finding the nearest point on the intersecting line to the
-        // // edge point.  Floating point errors often encountered using
-        // // planePlaneIntersect
-
-        // plane planeF(concaveEdgeLocalFeatPt, concaveEdgeDir);
-
-        // Foam::point concaveEdgeExternalPt = planeF.planePlaneIntersect
-        // (
-        //     planeA,
-        //     planeB
-        // );
-
-        // label concaveEdgeExternalPtI = number_of_vertices() + 4;
-
-        // // Redefine planes to be on the feature surfaces to project
-        // // through
-
-        // planeA = plane(featPt, concaveEdgePlaneANormal);
-
-        // planeB = plane(featPt, concaveEdgePlaneBNormal);
-
-        // Foam::point internalPtA =
-        //     concaveEdgeExternalPt
-        //   - 2*planeA.distance(concaveEdgeExternalPt)
-        //     *concaveEdgePlaneANormal;
-
-        // label internalPtAI = insertPoint
-        // (
-        //     internalPtA,
-        //     concaveEdgeExternalPtI
-        // );
-
-        // Foam::point internalPtB =
-        //     concaveEdgeExternalPt
-        //   - 2*planeB.distance(concaveEdgeExternalPt)
-        //    *concaveEdgePlaneBNormal;
-
-        // label internalPtBI = insertPoint
-        // (
-        //     internalPtB,
-        //     concaveEdgeExternalPtI
-        // );
-
-        // // TEMPORARY VARIABLE TO TEST
-        // vector convexEdgesPlaneNormal = -concaveEdgeDir;
-
-        // plane planeC(featPt, convexEdgesPlaneNormal);
-        // Foam::point externalPtD =
-        //     internalPtA
-        //   + 2*planeC.distance(internalPtA)*convexEdgesPlaneNormal;
-
-        // insertPoint(externalPtD, internalPtAI);
-
-        // Foam::point externalPtE =
-        //     internalPtB
-        //   + 2*planeC.distance(internalPtB)*convexEdgesPlaneNormal;
-
-        // insertPoint(externalPtE, internalPtBI);
-
-        // insertPoint(concaveEdgeExternalPt, internalPtAI);
-
-        // Info<< nl << "# featPt " << endl;
-        // meshTools::writeOBJ(Info, featPt);
-        // Info<< "# internalPtA" << endl;
-        // meshTools::writeOBJ(Info, internalPtA);
-        // Info<< "# internalPtB" << endl;
-        // meshTools::writeOBJ(Info, internalPtB);
-        // Info<< "# externalPtD" << endl;
-        // meshTools::writeOBJ(Info, externalPtD);
-        // Info<< "# externalPtE" << endl;
-        // meshTools::writeOBJ(Info, externalPtE);
-
-        // scalar totalAngle = radToDeg
-        // (
-        //     constant::mathematical::pi +
-        //     acos(mag(concaveEdgePlaneANormal & concaveEdgePlaneBNormal))
-        // );
-
-        // if (totalAngle > cvMeshControls().maxQuadAngle())
-        // {
-        //     // Add additional mitering points
-
-        //     scalar angleSign = 1.0;
-
-        //     if
-        //     (
-        //         geometryToConformTo_.outside
-        //         (
-        //             featPt - convexEdgesPlaneNormal*ppDist
-        //         )
-        //     )
-        //     {
-        //         angleSign = -1.0;
-        //     }
-
-        //     scalar phi =
-        //         angleSign*acos(concaveEdgeDir & -convexEdgesPlaneNormal);
-
-        //     scalar guard =
-        //     (
-        //         1 + sin(phi)*ppDist/mag
-        //         (
-        //             concaveEdgeLocalFeatPt - concaveEdgeExternalPt
-        //         )
-        //     )/cos(phi) - 1;
-
-        //     Foam::point internalPtF =
-        //         concaveEdgeExternalPt
-        //       + (2 + guard)*(concaveEdgeLocalFeatPt - concaveEdgeExternalPt);
-
-        //     label internalPtFI =
-        //         insertPoint(internalPtF, number_of_vertices() + 1);
-
-        //     Foam::point externalPtG =
-        //         internalPtF
-        //       + 2*planeC.distance(internalPtF) * convexEdgesPlaneNormal;
-
-        //     insertPoint(externalPtG, internalPtFI);
-        // }
-
-        // return true;
-
-        return false;
+        Info<< pFEdgeTypes << endl;
     }
-    else if (nExternal == 1 && nInternal == 2)
+
+    return allEdStat;
+}
+
+
+bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
+(
+    const extendedFeatureEdgeMesh& feMesh,
+    const labelList& pEds,
+    const pointFeatureEdgesTypes& pFEdgesTypes,
+    const List<extendedFeatureEdgeMesh::edgeStatus>& allEdStat,
+    const label ptI,
+    DynamicList<Foam::point>& pts,
+    DynamicList<label>& indices,
+    DynamicList<label>& types
+)
+{
+    if
+    (
+        pFEdgesTypes.nExternal == 2
+     && pFEdgesTypes.nInternal == 1
+     && pEds.size() == 3
+    )
+    {
+        Info<< "nExternal == 2 && nInternal == 1" << endl;
+
+        const Foam::point& featPt = feMesh.points()[ptI];
+
+        if (!positionOnThisProc(featPt))
+        {
+            return false;
+        }
+
+        const scalar ppDist = pointPairDistance(featPt);
+
+        const vectorField& normals = feMesh.normals();
+
+        const labelListList& edgeNormals = feMesh.edgeNormals();
+
+        label concaveEdgeI = -1;
+        labelList convexEdgesI(2, -1);
+        label nConvex = 0;
+
+        forAll(pEds, i)
+        {
+            const extendedFeatureEdgeMesh::edgeStatus& eS = allEdStat[i];
+
+            if (eS == extendedFeatureEdgeMesh::INTERNAL)
+            {
+                concaveEdgeI = pEds[i];
+            }
+            else if (eS == extendedFeatureEdgeMesh::EXTERNAL)
+            {
+                convexEdgesI[nConvex++] = pEds[i];
+            }
+            else if (eS == extendedFeatureEdgeMesh::FLAT)
+            {
+                WarningIn("Foam::conformalVoronoiMesh::"
+                    "createSpecialisedFeaturePoint")
+                    << "Edge " << eS << " is flat"
+                    << endl;
+            }
+            else
+            {
+                FatalErrorIn("Foam::conformalVoronoiMesh::"
+                    "createSpecialisedFeaturePoint")
+                    << "Edge " << eS << " not concave/convex"
+                    << exit(FatalError);
+            }
+        }
+
+        const vector& concaveEdgePlaneANormal =
+            normals[edgeNormals[concaveEdgeI][0]];
+
+        const vector& concaveEdgePlaneBNormal =
+            normals[edgeNormals[concaveEdgeI][1]];
+
+         // Intersect planes parallel to the concave edge planes offset
+         // by ppDist and the plane defined by featPt and the edge vector.
+         plane planeA
+         (
+             featPt + ppDist*concaveEdgePlaneANormal,
+             concaveEdgePlaneANormal
+         );
+
+         plane planeB
+         (
+             featPt + ppDist*concaveEdgePlaneBNormal,
+             concaveEdgePlaneBNormal
+         );
+
+        const vector& concaveEdgeDir = feMesh.edgeDirection
+        (
+            concaveEdgeI,
+            ptI
+        );
+
+        // Todo,needed later but want to get rid of this.
+        const Foam::point concaveEdgeLocalFeatPt = featPt + ppDist*concaveEdgeDir;
+
+        // Finding the nearest point on the intersecting line to the edge
+        // point. Floating point errors often occur using planePlaneIntersect
+
+        plane planeF(concaveEdgeLocalFeatPt, concaveEdgeDir);
+
+        const Foam::point concaveEdgeExternalPt = planeF.planePlaneIntersect
+        (
+            planeA,
+            planeB
+        );
+
+        // Redefine planes to be on the feature surfaces to project through
+
+        planeA = plane(featPt, concaveEdgePlaneANormal);
+
+        planeB = plane(featPt, concaveEdgePlaneBNormal);
+
+        const Foam::point internalPtA =
+            concaveEdgeExternalPt
+          - 2*planeA.distance(concaveEdgeExternalPt)
+            *concaveEdgePlaneANormal;
+
+        pts.append(internalPtA);
+        indices.append(0);
+        types.append(4);
+
+        const Foam::point internalPtB =
+            concaveEdgeExternalPt
+          - 2*planeB.distance(concaveEdgeExternalPt)
+            *concaveEdgePlaneBNormal;
+
+        pts.append(internalPtB);
+        indices.append(0);
+        types.append(3);
+
+        // Add the external points
+
+        Foam::point externalPtD;
+        Foam::point externalPtE;
+
+        vector convexEdgePlaneCNormal;
+        vector convexEdgePlaneDNormal;
+
+        const labelList& concaveEdgeNormals = edgeNormals[concaveEdgeI];
+        const labelList& convexEdgeANormals = edgeNormals[convexEdgesI[0]];
+        const labelList& convexEdgeBNormals = edgeNormals[convexEdgesI[1]];
+
+        forAll(concaveEdgeNormals, edgeNormalI)
+        {
+            forAll(convexEdgeANormals, edgeAnormalI)
+            {
+                const vector& concaveNormal
+                    = normals[concaveEdgeNormals[edgeNormalI]];
+                const vector& convexNormal
+                    = normals[convexEdgeANormals[edgeAnormalI]];
+
+                const scalar orientation = concaveNormal & convexNormal;
+
+                if (orientation <= 0.0)
+                {
+                    convexEdgePlaneCNormal = convexNormal;
+
+                    plane planeC(featPt, convexEdgePlaneCNormal);
+
+                    externalPtD =
+                        internalPtA
+                      + 2.0*planeC.distance(internalPtA)*convexEdgePlaneCNormal;
+
+                    pts.append(externalPtD);
+                    indices.append(0);
+                    types.append(-2);
+                }
+            }
+
+            forAll(convexEdgeBNormals, edgeBnormalI)
+            {
+                const vector& concaveNormal
+                    = normals[concaveEdgeNormals[edgeNormalI]];
+                const vector& convexNormal
+                    = normals[convexEdgeBNormals[edgeBnormalI]];
+
+                const scalar orientation = concaveNormal & convexNormal;
+
+                if (orientation <= 0.0)
+                {
+                    convexEdgePlaneDNormal = convexNormal;
+
+                    plane planeD(featPt, convexEdgePlaneDNormal);
+
+                    externalPtE =
+                        internalPtB
+                      + 2.0*planeD.distance(internalPtB)*convexEdgePlaneDNormal;
+
+                    pts.append(externalPtE);
+                    indices.append(0);
+                    types.append(-2);
+                }
+            }
+        }
+
+        pts.append(concaveEdgeExternalPt);
+        indices.append(0);
+        types.append(-4);
+
+        const scalar totalAngle = radToDeg
+        (
+            constant::mathematical::pi +
+            acos(mag(concaveEdgePlaneANormal & concaveEdgePlaneBNormal))
+        );
+
+        if (totalAngle > cvMeshControls().maxQuadAngle())
+        {
+            // Add additional mitering points
+            //scalar angleSign = 1.0;
+
+            vector convexEdgesPlaneNormal = convexEdgePlaneCNormal;
+
+            plane planeM(featPt, convexEdgesPlaneNormal);
+
+//            if
+//            (
+//                geometryToConformTo_.outside
+//                (
+//                    featPt - convexEdgesPlaneNormal*ppDist
+//                )
+//            )
+//            {
+//                angleSign = -1.0;
+//            }
+
+//            scalar phi =
+//                angleSign*acos(concaveEdgeDir & -convexEdgesPlaneNormal);
+//
+//            scalar guard =
+//            (
+//                1.0 + sin(phi)*ppDist/mag
+//                (
+//                    concaveEdgeLocalFeatPt - concaveEdgeExternalPt
+//                )
+//            )/cos(phi) - 1.0;
+
+            const Foam::point internalPtF =
+                concaveEdgeExternalPt
+            //+ (2.0 + guard)*(concaveEdgeLocalFeatPt - concaveEdgeExternalPt);
+              + 2.0*(concaveEdgeLocalFeatPt - concaveEdgeExternalPt);
+
+            pts.append(internalPtF);
+            indices.append(0);
+            types.append(1);
+
+            const Foam::point externalPtG =
+                internalPtF
+              //+ 2.0*planeM.distance(internalPtF) * convexEdgesPlaneNormal;
+              + 2.0*planeM.distance(internalPtF)*convexEdgesPlaneNormal;
+
+            pts.append(externalPtG);
+            indices.append(0);
+            types.append(-1);
+
+            if (debug)
+            {
+                Info<< nl << "# featPt " << endl;
+                meshTools::writeOBJ(Info, featPt);
+                Info<< "# internalPtA" << endl;
+                meshTools::writeOBJ(Info, internalPtA);
+                Info<< "# internalPtB" << endl;
+                meshTools::writeOBJ(Info, internalPtB);
+                Info<< "# externalPtD" << endl;
+                meshTools::writeOBJ(Info, externalPtD);
+                Info<< "# externalPtE" << endl;
+                meshTools::writeOBJ(Info, externalPtE);
+                Info<< "# internalPtF" << endl;
+                meshTools::writeOBJ(Info, internalPtF);
+                Info<< "# externalPtG" << endl;
+                meshTools::writeOBJ(Info, externalPtG);
+            }
+        }
+
+        return true;
+    }
+    else if (pFEdgesTypes.nExternal == 1 && pFEdgesTypes.nInternal == 2)
     {
         // Info<< "nExternal == 1 && nInternal == 2" << endl;
 
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C
new file mode 100644
index 0000000000..92041e715a
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C
@@ -0,0 +1,801 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "conformalVoronoiMesh.H"
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+void Foam::conformalVoronoiMesh::insertBoundingPoints()
+{
+    pointField farPts = geometryToConformTo_.globalBounds().points();
+
+    // Shift corners of bounds relative to origin
+    farPts -= geometryToConformTo_.globalBounds().midpoint();
+
+    // Scale the box up
+    farPts *= 10.0;
+
+    // Shift corners of bounds back to be relative to midpoint
+    farPts += geometryToConformTo_.globalBounds().midpoint();
+
+    limitBounds_ = treeBoundBox(farPts);
+
+    forAll(farPts, fPI)
+    {
+        insertPoint(farPts[fPI], Vb::vtFar);
+    }
+}
+
+
+void Foam::conformalVoronoiMesh::createEdgePointGroup
+(
+    const extendedFeatureEdgeMesh& feMesh,
+    const pointIndexHit& edHit,
+    DynamicList<Foam::point>& pts,
+    DynamicList<label>& indices,
+    DynamicList<label>& types
+)
+{
+    label edgeI = edHit.index();
+
+    extendedFeatureEdgeMesh::edgeStatus edStatus = feMesh.getEdgeStatus(edgeI);
+
+    switch (edStatus)
+    {
+        case extendedFeatureEdgeMesh::EXTERNAL:
+        {
+            createExternalEdgePointGroup(feMesh, edHit, pts, indices, types);
+            break;
+        }
+        case extendedFeatureEdgeMesh::INTERNAL:
+        {
+            createInternalEdgePointGroup(feMesh, edHit, pts, indices, types);
+            break;
+        }
+        case extendedFeatureEdgeMesh::FLAT:
+        {
+            createFlatEdgePointGroup(feMesh, edHit, pts, indices, types);
+            break;
+        }
+        case extendedFeatureEdgeMesh::OPEN:
+        {
+            createOpenEdgePointGroup(feMesh, edHit, pts, indices, types);
+            break;
+        }
+        case extendedFeatureEdgeMesh::MULTIPLE:
+        {
+            createMultipleEdgePointGroup(feMesh, edHit, pts, indices, types);
+            break;
+        }
+        case extendedFeatureEdgeMesh::NONE:
+        {
+            break;
+        }
+    }
+}
+
+
+void Foam::conformalVoronoiMesh::createExternalEdgePointGroup
+(
+    const extendedFeatureEdgeMesh& feMesh,
+    const pointIndexHit& edHit,
+    DynamicList<Foam::point>& pts,
+    DynamicList<label>& indices,
+    DynamicList<label>& types
+)
+{
+    const Foam::point& edgePt = edHit.hitPoint();
+
+    scalar ppDist = pointPairDistance(edgePt);
+
+    const vectorField& feNormals = feMesh.normals();
+    const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
+
+    // As this is an external edge, there are two normals by definition
+    const vector& nA = feNormals[edNormalIs[0]];
+    const vector& nB = feNormals[edNormalIs[1]];
+
+    if (mag(nA & nB) > 1 - SMALL)
+    {
+        // The normals are nearly parallel, so this is too sharp a feature to
+        // conform to.
+
+        return;
+    }
+
+    // Normalised distance of reference point from edge point
+    vector refVec((nA + nB)/(1 + (nA & nB)));
+
+    if (magSqr(refVec) > sqr(5.0))
+    {
+        // Limit the size of the conformation
+        ppDist *= 5.0/mag(refVec);
+
+        // Pout<< nl << "createExternalEdgePointGroup limit "
+        //     << "edgePt " << edgePt << nl
+        //     << "refVec " << refVec << nl
+        //     << "mag(refVec) " << mag(refVec) << nl
+        //     << "ppDist " << ppDist << nl
+        //     << "nA " << nA << nl
+        //     << "nB " << nB << nl
+        //     << "(nA & nB) " << (nA & nB) << nl
+        //     << endl;
+    }
+
+    // Convex. So refPt will be inside domain and hence a master point
+    Foam::point refPt = edgePt - ppDist*refVec;
+
+    // Result when the points are eventually inserted.
+    // Add number_of_vertices() at insertion of first vertex to all numbers:
+    // pt           index type
+    // refPt        0     1
+    // reflectedA   1     0
+    // reflectedB   2     0
+
+    // Insert the master point pairing the the first slave
+
+    pts.append(refPt);
+    indices.append(0);
+    types.append(1);
+
+    // Insert the slave points by reflecting refPt in both faces.
+    // with each slave refering to the master
+
+    Foam::point reflectedA = refPt + 2*ppDist*nA;
+    pts.append(reflectedA);
+    indices.append(0);
+    types.append(-1);
+
+    Foam::point reflectedB = refPt + 2*ppDist*nB;
+    pts.append(reflectedB);
+    indices.append(0);
+    types.append(-2);
+}
+
+
+void Foam::conformalVoronoiMesh::createInternalEdgePointGroup
+(
+    const extendedFeatureEdgeMesh& feMesh,
+    const pointIndexHit& edHit,
+    DynamicList<Foam::point>& pts,
+    DynamicList<label>& indices,
+    DynamicList<label>& types
+)
+{
+    const Foam::point& edgePt = edHit.hitPoint();
+
+    scalar ppDist = pointPairDistance(edgePt);
+
+    const vectorField& feNormals = feMesh.normals();
+    const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
+
+    // As this is an external edge, there are two normals by definition
+    const vector& nA = feNormals[edNormalIs[0]];
+    const vector& nB = feNormals[edNormalIs[1]];
+
+    if (mag(nA & nB) > 1 - SMALL)
+    {
+        // The normals are nearly parallel, so this is too sharp a feature to
+        // conform to.
+
+        return;
+    }
+
+    // Normalised distance of reference point from edge point
+    vector refVec((nA + nB)/(1 + (nA & nB)));
+
+    if (magSqr(refVec) > sqr(5.0))
+    {
+        // Limit the size of the conformation
+        ppDist *= 5.0/mag(refVec);
+
+        // Pout<< nl << "createInternalEdgePointGroup limit "
+        //     << "edgePt " << edgePt << nl
+        //     << "refVec " << refVec << nl
+        //     << "mag(refVec) " << mag(refVec) << nl
+        //     << "ppDist " << ppDist << nl
+        //     << "nA " << nA << nl
+        //     << "nB " << nB << nl
+        //     << "(nA & nB) " << (nA & nB) << nl
+        //     << endl;
+    }
+
+    // Concave. master and reflected points inside the domain.
+    Foam::point refPt = edgePt - ppDist*refVec;
+
+    // Generate reflected master to be outside.
+    Foam::point reflMasterPt = refPt + 2*(edgePt - refPt);
+
+    // Reflect reflMasterPt in both faces.
+    Foam::point reflectedA = reflMasterPt - 2*ppDist*nA;
+
+    Foam::point reflectedB = reflMasterPt - 2*ppDist*nB;
+
+    scalar totalAngle =
+        radToDeg(constant::mathematical::pi + acos(nA & nB));
+
+    // Number of quadrants the angle should be split into
+    int nQuads = int(totalAngle/cvMeshControls().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 = 2 + nAddPoints;
+
+    // Add number_of_vertices() at insertion of first vertex to all numbers:
+    // Result for nAddPoints 1 when the points are eventually inserted
+    // pt           index type
+    // reflectedA   0     2
+    // reflectedB   1     2
+    // reflMasterPt 2     0
+
+    // Result for nAddPoints 1 when the points are eventually inserted
+    // pt           index type
+    // reflectedA   0     3
+    // reflectedB   1     3
+    // refPt        2     3
+    // reflMasterPt 3     0
+
+    // Result for nAddPoints 2 when the points are eventually inserted
+    // pt           index type
+    // reflectedA   0     4
+    // reflectedB   1     4
+    // reflectedAa  2     4
+    // reflectedBb  3     4
+    // reflMasterPt 4     0
+
+    // Master A is inside.
+    pts.append(reflectedA);
+    indices.append(0);
+    types.append(reflectedMaster--);
+
+    // Master B is inside.
+    pts.append(reflectedB);
+    indices.append(0);
+    types.append(reflectedMaster--);
+
+    if (nAddPoints == 1)
+    {
+        // One additinal point is the reflection of the slave point,
+        // i.e. the original reference point
+        pts.append(refPt);
+        indices.append(0);
+        types.append(reflectedMaster--);
+    }
+    else if (nAddPoints == 2)
+    {
+        Foam::point reflectedAa = refPt + ppDist*nB;
+        pts.append(reflectedAa);
+        indices.append(0);
+        types.append(reflectedMaster--);
+
+        Foam::point reflectedBb = refPt + ppDist*nA;
+        pts.append(reflectedBb);
+        indices.append(0);
+        types.append(reflectedMaster--);
+    }
+
+    // Slave is outside.
+    pts.append(reflMasterPt);
+    indices.append(0);
+    types.append(-(nAddPoints + 2));
+}
+
+
+void Foam::conformalVoronoiMesh::createFlatEdgePointGroup
+(
+    const extendedFeatureEdgeMesh& feMesh,
+    const pointIndexHit& edHit,
+    DynamicList<Foam::point>& pts,
+    DynamicList<label>& indices,
+    DynamicList<label>& types
+)
+{
+    const Foam::point& edgePt = edHit.hitPoint();
+
+    const scalar ppDist = pointPairDistance(edgePt);
+
+    const vectorField& feNormals = feMesh.normals();
+    const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
+
+    // As this is a flat edge, there are two normals by definition
+    const vector& nA = feNormals[edNormalIs[0]];
+    const vector& nB = feNormals[edNormalIs[1]];
+
+    // Average normal to remove any bias to one side, although as this
+    // is a flat edge, the normals should be essentially the same
+    vector n = 0.5*(nA + nB);
+
+    // Direction along the surface to the control point, sense of edge
+    // direction not important, as +s and -s can be used because this
+    // is a flat edge
+    vector s = ppDist*(feMesh.edgeDirections()[edHit.index()] ^ n);
+
+    createPointPair(ppDist, edgePt + s, n, pts, indices, types);
+    createPointPair(ppDist, edgePt - s, n, pts, indices, types);
+}
+
+
+void Foam::conformalVoronoiMesh::createOpenEdgePointGroup
+(
+    const extendedFeatureEdgeMesh& feMesh,
+    const pointIndexHit& edHit,
+    DynamicList<Foam::point>& pts,
+    DynamicList<label>& indices,
+    DynamicList<label>& types
+)
+{
+    Info<< "NOT INSERTING OPEN EDGE POINT GROUP, NOT IMPLEMENTED" << endl;
+}
+
+
+void Foam::conformalVoronoiMesh::createMultipleEdgePointGroup
+(
+    const extendedFeatureEdgeMesh& feMesh,
+    const pointIndexHit& edHit,
+    DynamicList<Foam::point>& pts,
+    DynamicList<label>& indices,
+    DynamicList<label>& types
+)
+{
+    Info<< "NOT INSERTING MULTIPLE EDGE POINT GROUP, NOT IMPLEMENTED" << endl;
+}
+
+
+void Foam::conformalVoronoiMesh::reinsertBoundingPoints()
+{
+    tmp<pointField> farPts = limitBounds_.points();
+
+    forAll(farPts(), fPI)
+    {
+        insertPoint(farPts()[fPI], Vb::vtFar);
+    }
+}
+
+
+void Foam::conformalVoronoiMesh::reinsertFeaturePoints(bool distribute)
+{
+    Info<< nl << "Reinserting stored feature points" << endl;
+
+    label preReinsertionSize(number_of_vertices());
+
+    if (distribute)
+    {
+        DynamicList<Foam::point> pointsToInsert;
+        DynamicList<label> indices;
+        DynamicList<label> types;
+
+        for
+        (
+            std::list<Vb>::iterator vit=featureVertices_.begin();
+            vit != featureVertices_.end();
+            ++vit
+        )
+        {
+            pointsToInsert.append(topoint(vit->point()));
+            indices.append(vit->index());
+            types.append(vit->type());
+        }
+
+        // Insert distributed points
+        insertPoints(pointsToInsert, indices, types, true);
+
+        // Save points in new distribution
+        featureVertices_.clear();
+
+        forAll(pointsToInsert, pI)
+        {
+            featureVertices_.push_back
+            (
+                Vb(toPoint(pointsToInsert[pI]), indices[pI], types[pI])
+            );
+        }
+    }
+    else
+    {
+        for
+        (
+            std::list<Vb>::iterator vit=featureVertices_.begin();
+            vit != featureVertices_.end();
+            ++vit
+        )
+        {
+            // Assuming that all of the reinsertions are pair points, and that
+            // the index and type are relative, i.e. index 0 and type relative
+            // to it.
+            insertPoint
+            (
+                vit->point(),
+                vit->index() + number_of_vertices(),
+                vit->type() + number_of_vertices()
+            );
+        }
+    }
+
+    Info<< "    Reinserted "
+        << returnReduce
+        (
+            label(number_of_vertices()) - preReinsertionSize,
+            sumOp<label>()
+        )
+        << " vertices" << endl;
+}
+
+
+void Foam::conformalVoronoiMesh::createConvexFeaturePoints
+(
+    DynamicList<Foam::point>& pts,
+    DynamicList<label>& indices,
+    DynamicList<label>& types
+)
+{
+    const PtrList<extendedFeatureEdgeMesh>& feMeshes
+    (
+        geometryToConformTo_.features()
+    );
+
+    forAll(feMeshes, i)
+    {
+        const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
+
+        for
+        (
+            label ptI = feMesh.convexStart();
+            ptI < feMesh.concaveStart();
+            ptI++
+        )
+        {
+            const Foam::point& featPt = feMesh.points()[ptI];
+
+            if (!positionOnThisProc(featPt))
+            {
+                continue;
+            }
+
+            vectorField featPtNormals = feMesh.featurePointNormals(ptI);
+
+            vector cornerNormal = sum(featPtNormals);
+            cornerNormal /= mag(cornerNormal);
+
+            Foam::point internalPt =
+                featPt - pointPairDistance(featPt)*cornerNormal;
+
+            // Result when the points are eventually inserted (example n = 4)
+            // Add number_of_vertices() at insertion of first vertex to all
+            // numbers:
+            // pt           index type
+            // internalPt   0     1
+            // externalPt0  1     0
+            // externalPt1  2     0
+            // externalPt2  3     0
+            // externalPt3  4     0
+
+            pts.append(internalPt);
+            indices.append(0);
+            types.append(1);
+
+            label internalPtIndex = -1;
+
+            forAll(featPtNormals, nI)
+            {
+                const vector& n = featPtNormals[nI];
+
+                plane planeN = plane(featPt, n);
+
+                Foam::point externalPt =
+                    internalPt + 2.0 * planeN.distance(internalPt) * n;
+
+                pts.append(externalPt);
+                indices.append(0);
+                types.append(internalPtIndex--);
+            }
+        }
+    }
+}
+
+
+void Foam::conformalVoronoiMesh::createConcaveFeaturePoints
+(
+    DynamicList<Foam::point>& pts,
+    DynamicList<label>& indices,
+    DynamicList<label>& types
+)
+{
+    const PtrList<extendedFeatureEdgeMesh>& feMeshes
+    (
+        geometryToConformTo_.features()
+    );
+
+    forAll(feMeshes, i)
+    {
+        const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
+
+        for
+        (
+            label ptI = feMesh.concaveStart();
+            ptI < feMesh.mixedStart();
+            ptI++
+        )
+        {
+            const Foam::point& featPt = feMesh.points()[ptI];
+
+            if (!positionOnThisProc(featPt))
+            {
+                continue;
+            }
+
+            vectorField featPtNormals = feMesh.featurePointNormals(ptI);
+
+            vector cornerNormal = sum(featPtNormals);
+            cornerNormal /= mag(cornerNormal);
+
+            Foam::point externalPt =
+                featPt + pointPairDistance(featPt)*cornerNormal;
+
+            label externalPtIndex = featPtNormals.size();
+
+            // Result when the points are eventually inserted (example n = 5)
+            // Add number_of_vertices() at insertion of first vertex to all
+            // numbers:
+            // pt           index type
+            // internalPt0  0     5
+            // internalPt1  1     5
+            // internalPt2  2     5
+            // internalPt3  3     5
+            // internalPt4  4     5
+            // externalPt   5     4
+
+            forAll(featPtNormals, nI)
+            {
+                const vector& n = featPtNormals[nI];
+
+                plane planeN = plane(featPt, n);
+
+                Foam::point internalPt =
+                    externalPt - 2.0 * planeN.distance(externalPt) * n;
+
+                pts.append(internalPt);
+                indices.append(0);
+                types.append(externalPtIndex--);
+            }
+
+            pts.append(externalPt);
+            indices.append(0);
+            types.append(-1);
+        }
+    }
+}
+
+
+void Foam::conformalVoronoiMesh::createMixedFeaturePoints
+(
+    DynamicList<Foam::point>& pts,
+    DynamicList<label>& indices,
+    DynamicList<label>& types
+)
+{
+    if (cvMeshControls().mixedFeaturePointPPDistanceCoeff() < 0)
+    {
+        Info<<nl << "Skipping specialised handling for mixed feature points"
+            << endl;
+        return;
+    }
+
+    const PtrList<extendedFeatureEdgeMesh>& feMeshes
+    (
+        geometryToConformTo_.features()
+    );
+
+    forAll(feMeshes, i)
+    {
+        const extendedFeatureEdgeMesh& feMesh = feMeshes[i];
+        const labelListList& pointsEdges = feMesh.pointEdges();
+        const pointField& points = feMesh.points();
+
+        for
+        (
+            label ptI = feMesh.mixedStart();
+            ptI < feMesh.nonFeatureStart();
+            ptI++
+        )
+        {
+            const labelList& pEds = pointsEdges[ptI];
+
+            pointFeatureEdgesTypes pFEdgeTypes(ptI);
+
+            const List<extendedFeatureEdgeMesh::edgeStatus> allEdStat
+                = calcPointFeatureEdgesTypes(feMesh, pEds, pFEdgeTypes);
+
+            bool specialisedSuccess = false;
+
+            if (cvMeshControls().specialiseFeaturePoints())
+            {
+                specialisedSuccess = createSpecialisedFeaturePoint
+                (
+                    feMesh, pEds, pFEdgeTypes, allEdStat, ptI,
+                    pts, indices, types
+                );
+            }
+
+            if (!specialisedSuccess)
+            {
+                // Specialisations available for some mixed feature points.  For
+                // non-specialised feature points, inserting mixed internal and
+                // external edge groups at feature point.
+
+                // Skipping unsupported mixed feature point types
+                bool skipEdge = false;
+
+                forAll(pEds, e)
+                {
+                    const label edgeI = pEds[e];
+
+                    const extendedFeatureEdgeMesh::edgeStatus edStatus =
+                        feMesh.getEdgeStatus(edgeI);
+
+                    if
+                    (
+                        edStatus == extendedFeatureEdgeMesh::OPEN
+                     || edStatus == extendedFeatureEdgeMesh::MULTIPLE
+                    )
+                    {
+                        Info<< "Edge type " << edStatus
+                            << " found for mixed feature point " << ptI
+                            << ". Not supported."
+                            << endl;
+
+                        skipEdge = true;
+                    }
+                }
+
+                if (skipEdge)
+                {
+                    Info<< "Skipping point " << ptI << nl << endl;
+
+                    continue;
+                }
+
+                const Foam::point& pt = points[ptI];
+
+                if (!positionOnThisProc(pt))
+                {
+                    continue;
+                }
+
+                const scalar edgeGroupDistance = mixedFeaturePointDistance(pt);
+
+                forAll(pEds, e)
+                {
+                    const label edgeI = pEds[e];
+
+                    const Foam::point edgePt =
+                        pt + edgeGroupDistance*feMesh.edgeDirection(edgeI, ptI);
+
+                    const pointIndexHit edgeHit(true, edgePt, edgeI);
+
+                    createEdgePointGroup(feMesh, edgeHit, pts, indices, types);
+                }
+            }
+        }
+    }
+}
+
+
+void Foam::conformalVoronoiMesh::insertFeaturePoints()
+{
+    Info<< nl << "Conforming to feature points" << endl;
+
+    DynamicList<Foam::point> pts;
+    DynamicList<label> indices;
+    DynamicList<label> types;
+
+    const label preFeaturePointSize = number_of_vertices();
+
+    createConvexFeaturePoints(pts, indices, types);
+
+    createConcaveFeaturePoints(pts, indices, types);
+
+    createMixedFeaturePoints(pts, indices, types);
+
+    // Insert the created points, distributing to the appropriate processor
+    insertPoints(pts, indices, types, true);
+
+    if (cvMeshControls().objOutput())
+    {
+        writePoints("featureVertices.obj", pts);
+    }
+
+    label nFeatureVertices = number_of_vertices() - preFeaturePointSize;
+
+    if (Pstream::parRun())
+    {
+        reduce(nFeatureVertices, sumOp<label>());
+    }
+
+    if (nFeatureVertices > 0)
+    {
+        Info<< "    Inserted " << nFeatureVertices
+            << " feature vertices" << endl;
+    }
+
+    featureVertices_.clear();
+
+    forAll(pts, pI)
+    {
+        featureVertices_.push_back
+        (
+            Vb(toPoint(pts[pI]), indices[pI], types[pI])
+        );
+    }
+
+    constructFeaturePointLocations();
+}
+
+
+void Foam::conformalVoronoiMesh::constructFeaturePointLocations()
+{
+    DynamicList<Foam::point> ftPtLocs;
+
+    const PtrList<extendedFeatureEdgeMesh>& feMeshes
+    (
+        geometryToConformTo_.features()
+    );
+
+    forAll(feMeshes, i)
+    {
+        const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
+
+
+        if (cvMeshControls().mixedFeaturePointPPDistanceCoeff() < 0)
+        {
+            // Ignoring mixed feature points
+            for
+            (
+                label ptI = feMesh.convexStart();
+                ptI < feMesh.mixedStart();
+                ptI++
+            )
+            {
+                ftPtLocs.append(feMesh.points()[ptI]);
+            }
+        }
+        else
+        {
+            for
+            (
+                label ptI = feMesh.convexStart();
+                ptI < feMesh.nonFeatureStart();
+                ptI++
+            )
+            {
+                ftPtLocs.append(feMesh.points()[ptI]);
+            }
+        }
+    }
+
+    featurePointLocations_.transfer(ftPtLocs);
+}
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
index 6430eab8ce..421bf98248 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
@@ -620,6 +620,8 @@ void Foam::conformalVoronoiMesh::writeMesh
 
     writeCellSizes(mesh);
 
+    writeCellCentres(mesh);
+
     findRemainingProtrusionSet(mesh);
 }
 
@@ -805,6 +807,34 @@ void Foam::conformalVoronoiMesh::writeCellSizes
 }
 
 
+void Foam::conformalVoronoiMesh::writeCellCentres
+(
+    const fvMesh& mesh
+) const
+{
+    Info<< "Writing components of cellCentre positions to volScalarFields"
+        << " ccx, ccy, ccz in " <<  runTime_.timeName() << endl;
+
+    for (direction i=0; i<vector::nComponents; i++)
+    {
+        volScalarField cci
+        (
+            IOobject
+            (
+                "cc" + word(vector::componentNames[i]),
+                runTime_.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            mesh.C().component(i)
+        );
+
+        cci.write();
+    }
+}
+
+
 void Foam::conformalVoronoiMesh::findRemainingProtrusionSet
 (
     const fvMesh& mesh
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H
index 7e5ea9112e..7ab4f10853 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H
@@ -33,6 +33,10 @@ Description
     An indexed form of CGAL::Triangulation_cell_base_3<K> used to keep
     track of the Delaunay cells (tets) in the tessellation.
 
+SourceFiles
+    indexedCellI.H
+    indexedCell.C
+
 \*---------------------------------------------------------------------------*/
 
 #ifndef indexedCell_H
@@ -119,248 +123,72 @@ public:
     };
 
 
-    indexedCell()
-    :
-        Cb(),
-        index_(ctFar),
-        filterCount_(0)
-    {}
+    // Constructors
 
+        inline indexedCell();
 
-    indexedCell
-    (
-        Vertex_handle v0, Vertex_handle v1, Vertex_handle v2, Vertex_handle v3
-    )
-    :
-        Cb(v0, v1, v2, v3),
-        index_(ctFar),
-        filterCount_(0)
-    {}
-
-
-    indexedCell
-    (
-        Vertex_handle v0,
-        Vertex_handle v1,
-        Vertex_handle v2,
-        Vertex_handle v3,
-        Cell_handle n0,
-        Cell_handle n1,
-        Cell_handle n2,
-        Cell_handle n3
-    )
-    :
-        Cb(v0, v1, v2, v3, n0, n1, n2, n3),
-        index_(ctFar),
-        filterCount_(0)
-    {}
-
-
-    int& cellIndex()
-    {
-        return index_;
-    }
-
-
-    int cellIndex() const
-    {
-        return index_;
-    }
-
-    inline bool farCell() const
-    {
-        return index_ == ctFar;
-    }
-
-    inline int& filterCount()
-    {
-        return filterCount_;
-    }
-
-
-    inline int filterCount() const
-    {
-        return filterCount_;
-    }
-
-
-    //- Is the Delaunay cell real, i.e. any real vertex
-    inline bool real() const
-    {
-        return
+        inline indexedCell
         (
-            this->vertex(0)->real()
-         || this->vertex(1)->real()
-         || this->vertex(2)->real()
-         || this->vertex(3)->real()
-        );
-    }
-
-
-    //- Does the Dual vertex form part of a processor patch
-    inline bool parallelDualVertex() const
-    {
-        return
-        (
-            this->vertex(0)->referred()
-         || this->vertex(1)->referred()
-         || this->vertex(2)->referred()
-         || this->vertex(3)->referred()
-        );
-    }
-
-    //- Does the Dual vertex form part of a processor patch
-    inline Foam::label dualVertexMasterProc() const
-    {
-        if (!parallelDualVertex())
-        {
-            return -1;
-        }
-
-        // The master processor is the lowest numbered of the four on this tet.
-
-        int masterProc = Foam::Pstream::nProcs() + 1;
-
-        for (int i = 0; i < 4; i++)
-        {
-            if (this->vertex(i)->referred())
-            {
-                masterProc = min(masterProc, this->vertex(i)->procIndex());
-            }
-            else
-            {
-                masterProc = min(masterProc, Foam::Pstream::myProcNo());
-            }
-        }
-
-        return masterProc;
-    }
-
-
-    inline Foam::FixedList<Foam::label, 4> processorsAttached() const
-    {
-        if (!parallelDualVertex())
-        {
-            return Foam::FixedList<Foam::label, 4>(Foam::Pstream::myProcNo());
-        }
-
-        Foam::FixedList<Foam::label, 4> procsAttached
-        (
-            Foam::Pstream::myProcNo()
+            Vertex_handle v0,
+            Vertex_handle v1,
+            Vertex_handle v2,
+            Vertex_handle v3
         );
 
-        for (int i = 0; i < 4; i++)
-        {
-            if (this->vertex(i)->referred())
-            {
-                procsAttached[i] = this->vertex(i)->procIndex();
-            }
-        }
-
-        return procsAttached;
-    }
-
-
-    // Using the globalIndex object, return a list of four (sorted) global
-    // Delaunay vertex indices that uniquely identify this tet in parallel
-    inline Foam::tetCell vertexGlobalIndices
-    (
-        const Foam::globalIndex& globalDelaunayVertexIndices
-    ) const
-    {
-        // tetVertexGlobalIndices
-        Foam::tetCell tVGI;
-
-        for (int i = 0; i < 4; i++)
-        {
-            Vertex_handle v = this->vertex(i);
-
-            // Finding the global index of each Delaunay vertex
-
-            if (v->referred())
-            {
-                // Referred vertices may have negative indices
-
-                tVGI[i] = globalDelaunayVertexIndices.toGlobal
-                (
-                    v->procIndex(),
-                    Foam::mag(v->index())
-                );
-            }
-            else
-            {
-                tVGI[i] = globalDelaunayVertexIndices.toGlobal
-                (
-                    Foam::Pstream::myProcNo(),
-                    v->index()
-                );
-            }
-        }
-
-        // bubble sort
-        for (int i = 0; i < tVGI.size(); i++)
-        {
-            for (int j = tVGI.size() - 1 ; j > i; j--)
-            {
-                if (tVGI[j - 1] > tVGI[j])
-                {
-                    Foam::Swap(tVGI[j - 1], tVGI[j]);
-                }
-            }
-        }
-
-        return tVGI;
-    }
-
-
-    // Is the Delaunay cell part of the final dual mesh, i.e. any vertex form
-    // part of the internal or boundary definition
-    inline bool internalOrBoundaryDualVertex() const
-    {
-        return
+        inline indexedCell
         (
-            this->vertex(0)->internalOrBoundaryPoint()
-         || this->vertex(1)->internalOrBoundaryPoint()
-         || this->vertex(2)->internalOrBoundaryPoint()
-         || this->vertex(3)->internalOrBoundaryPoint()
+            Vertex_handle v0,
+            Vertex_handle v1,
+            Vertex_handle v2,
+            Vertex_handle v3,
+            Cell_handle n0,
+            Cell_handle n1,
+            Cell_handle n2,
+            Cell_handle n3
         );
-    }
 
 
-    // Is the Delaunay cell real or referred (or mixed), i.e. all vertices form
-    // part of the real or referred internal or boundary definition
-    inline bool anyInternalOrBoundaryDualVertex() const
-    {
-        return
+    // Member Functions
+
+        inline int& cellIndex();
+
+        inline int cellIndex() const;
+
+        inline bool farCell() const;
+
+        inline int& filterCount();
+
+        inline int filterCount() const;
+
+        //- Is the Delaunay cell real, i.e. any real vertex
+        inline bool real() const;
+
+        //- Does the Dual vertex form part of a processor patch
+        inline bool parallelDualVertex() const;
+
+        //- Does the Dual vertex form part of a processor patch
+        inline Foam::label dualVertexMasterProc() const;
+
+        inline Foam::FixedList<Foam::label, 4> processorsAttached() const;
+
+        //- Using the globalIndex object, return a list of four (sorted) global
+        //  Delaunay vertex indices that uniquely identify this tet in parallel
+        inline Foam::tetCell vertexGlobalIndices
         (
-            this->vertex(0)->anyInternalOrBoundaryPoint()
-         || this->vertex(1)->anyInternalOrBoundaryPoint()
-         || this->vertex(2)->anyInternalOrBoundaryPoint()
-         || this->vertex(3)->anyInternalOrBoundaryPoint()
-        );
-    }
+            const Foam::globalIndex& globalDelaunayVertexIndices
+        ) const;
 
+        //- Is the Delaunay cell part of the final dual mesh, i.e. any vertex
+        //  form part of the internal or boundary definition
+        inline bool internalOrBoundaryDualVertex() const;
 
-    // A dual vertex on the boundary will result from a Delaunay cell with
-    // least one Delaunay vertex outside and at least one inside
-    inline bool boundaryDualVertex() const
-    {
-        return
-        (
-            (
-               this->vertex(0)->internalOrBoundaryPoint()
-            || this->vertex(1)->internalOrBoundaryPoint()
-            || this->vertex(2)->internalOrBoundaryPoint()
-            || this->vertex(3)->internalOrBoundaryPoint()
-            )
-         && (
-               !this->vertex(0)->internalOrBoundaryPoint()
-            || !this->vertex(1)->internalOrBoundaryPoint()
-            || !this->vertex(2)->internalOrBoundaryPoint()
-            || !this->vertex(3)->internalOrBoundaryPoint()
-            )
-        );
-    }
+        //- Is the Delaunay cell real or referred (or mixed), i.e. all vertices
+        //  form part of the real or referred internal or boundary definition
+        inline bool anyInternalOrBoundaryDualVertex() const;
+
+        //- A dual vertex on the boundary will result from a Delaunay cell with
+        //  least one Delaunay vertex outside and at least one inside
+        inline bool boundaryDualVertex() const;
 
 
     // Info
@@ -387,6 +215,8 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#include "indexedCellI.H"
+
 #ifdef NoRepository
 #   include "indexedCell.C"
 #endif
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCellI.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCellI.H
new file mode 100644
index 0000000000..5fbd378f9d
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCellI.H
@@ -0,0 +1,288 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+    As a special exception, you have permission to link this program with the
+    CGAL library and distribute executables, as long as you follow the
+    requirements of the GNU GPL in regard to all of the software in the
+    executable aside from CGAL.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Gt, class Cb>
+CGAL::indexedCell<Gt, Cb>::indexedCell()
+:
+    Cb(),
+    index_(ctFar),
+    filterCount_(0)
+{}
+
+
+template<class Gt, class Cb>
+CGAL::indexedCell<Gt, Cb>::indexedCell
+(
+    Vertex_handle v0, Vertex_handle v1, Vertex_handle v2, Vertex_handle v3
+)
+:
+    Cb(v0, v1, v2, v3),
+    index_(ctFar),
+    filterCount_(0)
+{}
+
+
+template<class Gt, class Cb>
+CGAL::indexedCell<Gt, Cb>::indexedCell
+(
+    Vertex_handle v0,
+    Vertex_handle v1,
+    Vertex_handle v2,
+    Vertex_handle v3,
+    Cell_handle n0,
+    Cell_handle n1,
+    Cell_handle n2,
+    Cell_handle n3
+)
+:
+    Cb(v0, v1, v2, v3, n0, n1, n2, n3),
+    index_(ctFar),
+    filterCount_(0)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Gt, class Cb>
+int& CGAL::indexedCell<Gt, Cb>::cellIndex()
+{
+    return index_;
+}
+
+
+template<class Gt, class Cb>
+int CGAL::indexedCell<Gt, Cb>::cellIndex() const
+{
+    return index_;
+}
+
+
+template<class Gt, class Cb>
+inline bool CGAL::indexedCell<Gt, Cb>::farCell() const
+{
+    return index_ == ctFar;
+}
+
+
+template<class Gt, class Cb>
+inline int& CGAL::indexedCell<Gt, Cb>::filterCount()
+{
+    return filterCount_;
+}
+
+
+template<class Gt, class Cb>
+inline int CGAL::indexedCell<Gt, Cb>::filterCount() const
+{
+    return filterCount_;
+}
+
+
+template<class Gt, class Cb>
+inline bool CGAL::indexedCell<Gt, Cb>::real() const
+{
+    return
+    (
+        this->vertex(0)->real()
+     || this->vertex(1)->real()
+     || this->vertex(2)->real()
+     || this->vertex(3)->real()
+    );
+}
+
+
+template<class Gt, class Cb>
+inline bool CGAL::indexedCell<Gt, Cb>::parallelDualVertex() const
+{
+    return
+    (
+        this->vertex(0)->referred()
+     || this->vertex(1)->referred()
+     || this->vertex(2)->referred()
+     || this->vertex(3)->referred()
+    );
+}
+
+
+template<class Gt, class Cb>
+inline Foam::label CGAL::indexedCell<Gt, Cb>::dualVertexMasterProc() const
+{
+    if (!parallelDualVertex())
+    {
+        return -1;
+    }
+
+    // The master processor is the lowest numbered of the four on this tet.
+
+    int masterProc = Foam::Pstream::nProcs() + 1;
+
+    for (int i = 0; i < 4; i++)
+    {
+        if (this->vertex(i)->referred())
+        {
+            masterProc = min(masterProc, this->vertex(i)->procIndex());
+        }
+        else
+        {
+            masterProc = min(masterProc, Foam::Pstream::myProcNo());
+        }
+    }
+
+    return masterProc;
+}
+
+
+template<class Gt, class Cb>
+inline Foam::FixedList<Foam::label, 4>
+CGAL::indexedCell<Gt, Cb>::processorsAttached() const
+{
+    if (!parallelDualVertex())
+    {
+        return Foam::FixedList<Foam::label, 4>(Foam::Pstream::myProcNo());
+    }
+
+    Foam::FixedList<Foam::label, 4> procsAttached
+    (
+        Foam::Pstream::myProcNo()
+    );
+
+    for (int i = 0; i < 4; i++)
+    {
+        if (this->vertex(i)->referred())
+        {
+            procsAttached[i] = this->vertex(i)->procIndex();
+        }
+    }
+
+    return procsAttached;
+}
+
+
+template<class Gt, class Cb>
+inline Foam::tetCell CGAL::indexedCell<Gt, Cb>::vertexGlobalIndices
+(
+    const Foam::globalIndex& globalDelaunayVertexIndices
+) const
+{
+    // tetVertexGlobalIndices
+    Foam::tetCell tVGI;
+
+    for (int i = 0; i < 4; i++)
+    {
+        Vertex_handle v = this->vertex(i);
+
+        // Finding the global index of each Delaunay vertex
+
+        if (v->referred())
+        {
+            // Referred vertices may have negative indices
+
+            tVGI[i] = globalDelaunayVertexIndices.toGlobal
+            (
+                v->procIndex(),
+                Foam::mag(v->index())
+            );
+        }
+        else
+        {
+            tVGI[i] = globalDelaunayVertexIndices.toGlobal
+            (
+                Foam::Pstream::myProcNo(),
+                v->index()
+            );
+        }
+    }
+
+    // bubble sort
+    for (int i = 0; i < tVGI.size(); i++)
+    {
+        for (int j = tVGI.size() - 1 ; j > i; j--)
+        {
+            if (tVGI[j - 1] > tVGI[j])
+            {
+                Foam::Swap(tVGI[j - 1], tVGI[j]);
+            }
+        }
+    }
+
+    return tVGI;
+}
+
+
+template<class Gt, class Cb>
+inline bool CGAL::indexedCell<Gt, Cb>::internalOrBoundaryDualVertex() const
+{
+    return
+    (
+        this->vertex(0)->internalOrBoundaryPoint()
+     || this->vertex(1)->internalOrBoundaryPoint()
+     || this->vertex(2)->internalOrBoundaryPoint()
+     || this->vertex(3)->internalOrBoundaryPoint()
+    );
+}
+
+
+template<class Gt, class Cb>
+inline bool CGAL::indexedCell<Gt, Cb>::anyInternalOrBoundaryDualVertex() const
+{
+    return
+    (
+        this->vertex(0)->anyInternalOrBoundaryPoint()
+     || this->vertex(1)->anyInternalOrBoundaryPoint()
+     || this->vertex(2)->anyInternalOrBoundaryPoint()
+     || this->vertex(3)->anyInternalOrBoundaryPoint()
+    );
+}
+
+
+template<class Gt, class Cb>
+inline bool CGAL::indexedCell<Gt, Cb>::boundaryDualVertex() const
+{
+    return
+    (
+        (
+           this->vertex(0)->internalOrBoundaryPoint()
+        || this->vertex(1)->internalOrBoundaryPoint()
+        || this->vertex(2)->internalOrBoundaryPoint()
+        || this->vertex(3)->internalOrBoundaryPoint()
+        )
+     && (
+           !this->vertex(0)->internalOrBoundaryPoint()
+        || !this->vertex(1)->internalOrBoundaryPoint()
+        || !this->vertex(2)->internalOrBoundaryPoint()
+        || !this->vertex(3)->internalOrBoundaryPoint()
+        )
+    );
+}
+
+
+// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H
index 0dae13edb3..a9b997f6cb 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H
@@ -33,6 +33,10 @@ Description
     An indexed form of CGAL::Triangulation_vertex_base_3<K> used to keep
     track of the Delaunay vertices in the tessellation.
 
+SourceFiles
+    indexedVertexI.H
+    indexedVertex.C
+
 \*---------------------------------------------------------------------------*/
 
 #ifndef indexedVertex_H
@@ -116,267 +120,99 @@ public:
     };
 
 
-    indexedVertex()
-    :
-        Vb(),
-        index_(vtInternal),
-        type_(vtInternal),
-        alignment_(),
-        targetCellSize_(0.0)
-    {}
+    // Constructors
 
+        inline indexedVertex();
 
-    indexedVertex(const Point& p)
-    :
-        Vb(p),
-        index_(vtInternal),
-        type_(vtInternal),
-        alignment_(),
-        targetCellSize_(0.0)
-    {}
+        inline indexedVertex(const Point& p);
 
+        inline indexedVertex(const Point& p, int index, int type);
 
-    indexedVertex(const Point& p, int index, int type)
-    :
-        Vb(p),
-        index_(index),
-        type_(type),
-        alignment_(),
-        targetCellSize_(0.0)
-    {}
+        inline indexedVertex(const Point& p, Cell_handle f);
 
+        inline indexedVertex(Cell_handle f);
 
-    indexedVertex(const Point& p, Cell_handle f)
-    :
-        Vb(f, p),
-        index_(vtInternal),
-        type_(vtInternal),
-        alignment_(),
-        targetCellSize_(0.0)
-    {}
 
+    // Member Functions
 
-    indexedVertex(Cell_handle f)
-    :
-        Vb(f),
-        index_(vtInternal),
-        type_(vtInternal),
-        alignment_(),
-        targetCellSize_(0.0)
-    {}
+        inline int& index();
 
+        inline int index() const;
 
-    int& index()
-    {
-        return index_;
-    }
+        inline int& type();
 
+        inline int type() const;
 
-    int index() const
-    {
-        return index_;
-    }
+        inline Foam::tensor& alignment();
 
+        inline const Foam::tensor& alignment() const;
 
-    int& type()
-    {
-        return type_;
-    }
+        inline Foam::scalar& targetCellSize();
 
+        inline Foam::scalar targetCellSize() const;
 
-    int type() const
-    {
-        return type_;
-    }
+        inline bool uninitialised() const;
 
+        //- Is point a far-point
+        inline bool farPoint() const;
 
-    inline Foam::tensor& alignment()
-    {
-        return alignment_;
-    }
+        //- Is point internal, i.e. not on boundary
+        inline bool internalPoint() const;
 
+        //- Is point internal, i.e. not on boundary, external query.
+        inline static bool internalPoint(int type);
 
-    inline const Foam::tensor& alignment() const
-    {
-        return alignment_;
-    }
+        // is this a referred vertex
+        inline bool referred() const;
 
+        // is this a referred internal or boundary vertex
+        inline bool referredInternalOrBoundaryPoint() const;
 
-    inline Foam::scalar& targetCellSize()
-    {
-        return targetCellSize_;
-    }
+        // is this a referred external (pair slave) vertex
+        inline bool referredExternal() const;
 
+        // Is this a "real" point on this processor, i.e. is it internal or part
+        // of the boundary description, and not a "far" or "referred" point
+        inline bool real() const;
 
-    inline Foam::scalar targetCellSize() const
-    {
-        return targetCellSize_;
-    }
+        // For referred vertices, what is the original processor index
+        inline int procIndex() const;
 
+        inline static int encodeProcIndex(int procI);
 
-    inline bool uninitialised() const
-    {
-        return type_ == vtInternal && index_ == vtInternal;
-    }
+        //- Set the point to be internal
+        inline void setInternal();
 
+        //- Is point internal and near the boundary
+        inline bool nearBoundary() const;
 
-    //- Is point a far-point
-    inline bool farPoint() const
-    {
-        return type_ == vtFar;
-    }
+        //- Set the point to be near the boundary
+        inline void setNearBoundary();
 
+        //- Either master or slave of pointPair.
+        inline bool pairPoint() const;
 
-    //- Is point internal, i.e. not on boundary
-    inline bool internalPoint() const
-    {
-        return internalPoint(type_);
-    }
+        //- Master of a pointPair is the lowest numbered one.
+        inline bool ppMaster() const;
 
+        //- Master of a pointPair is the lowest numbered one, external query.
+        inline static bool ppMaster(int index, int type);
 
-    //- Is point internal, i.e. not on boundary, external query.
-    inline static bool internalPoint(int type)
-    {
-        return type <= vtInternal;
-    }
+        //- Slave of a pointPair is the highest numbered one.
+        inline bool ppSlave() const;
 
+        //- Either original internal point or master of pointPair.
+        inline bool internalOrBoundaryPoint() const;
 
-    // is this a referred vertex
-    inline bool referred() const
-    {
-        return (type_ < 0 && type_ > vtFar);
-    }
+        //- Either original internal point or master of pointPair.
+        //  External query.
+        inline static bool internalOrBoundaryPoint(int index, int type);
 
+        //- Is point near the boundary or part of the boundary definition
+        inline bool nearOrOnBoundary() const;
 
-    // is this a referred internal or boundary vertex
-    inline bool referredInternalOrBoundaryPoint() const
-    {
-        return referred() && index_ >= 0;
-    }
-
-
-    // is this a referred external (pair slave) vertex
-    inline bool referredExternal() const
-    {
-        return referred() && index_ < 0;
-    }
-
-
-    // is this a "real" point on this processor, i.e. is it internal or part of
-    // the boundary description, and not a "far" or "referred" point
-    inline bool real() const
-    {
-        return internalPoint() || pairPoint();
-    }
-
-
-    // For referred vertices, what is the original processor index
-    inline int procIndex() const
-    {
-        if (referred())
-        {
-            return -(type_ + 1);
-        }
-        else
-        {
-            return -1;
-        }
-    }
-
-
-    inline static int encodeProcIndex(int procI)
-    {
-        return -(procI + 1);
-    }
-
-
-    //- Set the point to be internal
-    inline void setInternal()
-    {
-        type_ = vtInternal;
-    }
-
-
-    //- Is point internal and near the boundary
-    inline bool nearBoundary() const
-    {
-        return type_ == vtNearBoundary;
-    }
-
-
-    //- Set the point to be near the boundary
-    inline void setNearBoundary()
-    {
-        type_ = vtNearBoundary;
-    }
-
-
-    //- Either master or slave of pointPair.
-    inline bool pairPoint() const
-    {
-        return type_ >= 0;
-    }
-
-
-    //- Master of a pointPair is the lowest numbered one.
-    inline bool ppMaster() const
-    {
-        return ppMaster(index_, type_);
-    }
-
-
-    //- Master of a pointPair is the lowest numbered one, external query.
-    inline static bool ppMaster(int index, int type)
-    {
-        if (index >= 0 && type > index)
-        {
-            return true;
-        }
-
-        return false;
-    }
-
-
-    //- Slave of a pointPair is the highest numbered one.
-    inline bool ppSlave() const
-    {
-        if (type_ >= 0 && type_ < index_)
-        {
-            return true;
-        }
-        else
-        {
-            return false;
-        }
-    }
-
-
-    //- Either original internal point or master of pointPair.
-    inline bool internalOrBoundaryPoint() const
-    {
-        return internalOrBoundaryPoint(index_, type_);
-    }
-
-
-    //- Either original internal point or master of pointPair, external query.
-    inline static bool internalOrBoundaryPoint(int index, int type)
-    {
-        return internalPoint(type) || ppMaster(index, type);
-    }
-
-
-    //- Is point near the boundary or part of the boundary definition
-    inline bool nearOrOnBoundary() const
-    {
-        return pairPoint() || nearBoundary();
-    }
-
-
-    //- Either a real or referred internal or boundary point
-    inline bool anyInternalOrBoundaryPoint() const
-    {
-        return internalOrBoundaryPoint() || referredInternalOrBoundaryPoint();
-    }
+        //- Either a real or referred internal or boundary point
+        inline bool anyInternalOrBoundaryPoint() const;
 
 
     // inline void operator=(const Delaunay::Finite_vertices_iterator vit)
@@ -414,6 +250,8 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#include "indexedVertexI.H"
+
 #ifdef NoRepository
 #   include "indexedVertex.C"
 #endif
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertexI.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertexI.H
new file mode 100644
index 0000000000..86479321a9
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertexI.H
@@ -0,0 +1,321 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+    As a special exception, you have permission to link this program with the
+    CGAL library and distribute executables, as long as you follow the
+    requirements of the GNU GPL in regard to all of the software in the
+    executable aside from CGAL.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Gt, class Vb>
+inline CGAL::indexedVertex<Gt, Vb>::indexedVertex()
+:
+    Vb(),
+    index_(vtInternal),
+    type_(vtInternal),
+    alignment_(),
+    targetCellSize_(0.0)
+{}
+
+
+template<class Gt, class Vb>
+inline CGAL::indexedVertex<Gt, Vb>::indexedVertex(const Point& p)
+:
+    Vb(p),
+    index_(vtInternal),
+    type_(vtInternal),
+    alignment_(),
+    targetCellSize_(0.0)
+{}
+
+
+template<class Gt, class Vb>
+inline CGAL::indexedVertex<Gt, Vb>::indexedVertex
+(
+    const Point& p,
+    int index,
+    int type
+)
+:
+    Vb(p),
+    index_(index),
+    type_(type),
+    alignment_(),
+    targetCellSize_(0.0)
+{}
+
+
+template<class Gt, class Vb>
+inline CGAL::indexedVertex<Gt, Vb>::indexedVertex(const Point& p, Cell_handle f)
+:
+    Vb(f, p),
+    index_(vtInternal),
+    type_(vtInternal),
+    alignment_(),
+    targetCellSize_(0.0)
+{}
+
+
+template<class Gt, class Vb>
+inline CGAL::indexedVertex<Gt, Vb>::indexedVertex(Cell_handle f)
+:
+    Vb(f),
+    index_(vtInternal),
+    type_(vtInternal),
+    alignment_(),
+    targetCellSize_(0.0)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Gt, class Vb>
+inline int& CGAL::indexedVertex<Gt, Vb>::index()
+{
+    return index_;
+}
+
+
+template<class Gt, class Vb>
+inline int CGAL::indexedVertex<Gt, Vb>::index() const
+{
+    return index_;
+}
+
+
+template<class Gt, class Vb>
+inline int& CGAL::indexedVertex<Gt, Vb>::type()
+{
+    return type_;
+}
+
+
+template<class Gt, class Vb>
+inline int CGAL::indexedVertex<Gt, Vb>::type() const
+{
+    return type_;
+}
+
+
+template<class Gt, class Vb>
+inline Foam::tensor& CGAL::indexedVertex<Gt, Vb>::alignment()
+{
+    return alignment_;
+}
+
+
+template<class Gt, class Vb>
+inline const Foam::tensor& CGAL::indexedVertex<Gt, Vb>::alignment() const
+{
+    return alignment_;
+}
+
+
+template<class Gt, class Vb>
+inline Foam::scalar& CGAL::indexedVertex<Gt, Vb>::targetCellSize()
+{
+    return targetCellSize_;
+}
+
+
+template<class Gt, class Vb>
+inline Foam::scalar CGAL::indexedVertex<Gt, Vb>::targetCellSize() const
+{
+    return targetCellSize_;
+}
+
+
+template<class Gt, class Vb>
+inline bool CGAL::indexedVertex<Gt, Vb>::uninitialised() const
+{
+    return type_ == vtInternal && index_ == vtInternal;
+}
+
+
+template<class Gt, class Vb>
+inline bool CGAL::indexedVertex<Gt, Vb>::farPoint() const
+{
+    return type_ == vtFar;
+}
+
+
+template<class Gt, class Vb>
+inline bool CGAL::indexedVertex<Gt, Vb>::internalPoint() const
+{
+    return internalPoint(type_);
+}
+
+
+template<class Gt, class Vb>
+inline bool CGAL::indexedVertex<Gt, Vb>::internalPoint(int type)
+{
+    return type <= vtInternal;
+}
+
+
+template<class Gt, class Vb>
+inline bool CGAL::indexedVertex<Gt, Vb>::referred() const
+{
+    return (type_ < 0 && type_ > vtFar);
+}
+
+
+template<class Gt, class Vb>
+inline bool CGAL::indexedVertex<Gt, Vb>::referredInternalOrBoundaryPoint() const
+{
+    return referred() && index_ >= 0;
+}
+
+
+template<class Gt, class Vb>
+inline bool CGAL::indexedVertex<Gt, Vb>::referredExternal() const
+{
+    return referred() && index_ < 0;
+}
+
+
+template<class Gt, class Vb>
+inline bool CGAL::indexedVertex<Gt, Vb>::real() const
+{
+    return internalPoint() || pairPoint();
+}
+
+
+template<class Gt, class Vb>
+inline int CGAL::indexedVertex<Gt, Vb>::procIndex() const
+{
+    if (referred())
+    {
+        return -(type_ + 1);
+    }
+    else
+    {
+        return -1;
+    }
+}
+
+
+template<class Gt, class Vb>
+inline int CGAL::indexedVertex<Gt, Vb>::encodeProcIndex(int procI)
+{
+    return -(procI + 1);
+}
+
+
+template<class Gt, class Vb>
+inline void CGAL::indexedVertex<Gt, Vb>::setInternal()
+{
+    type_ = vtInternal;
+}
+
+
+template<class Gt, class Vb>
+inline bool CGAL::indexedVertex<Gt, Vb>::nearBoundary() const
+{
+    return type_ == vtNearBoundary;
+}
+
+
+template<class Gt, class Vb>
+inline void CGAL::indexedVertex<Gt, Vb>::setNearBoundary()
+{
+    type_ = vtNearBoundary;
+}
+
+
+template<class Gt, class Vb>
+inline bool CGAL::indexedVertex<Gt, Vb>::pairPoint() const
+{
+    return type_ >= 0;
+}
+
+
+template<class Gt, class Vb>
+inline bool CGAL::indexedVertex<Gt, Vb>::ppMaster() const
+{
+    return ppMaster(index_, type_);
+}
+
+
+template<class Gt, class Vb>
+inline bool CGAL::indexedVertex<Gt, Vb>::ppMaster(int index, int type)
+{
+    if (index >= 0 && type > index)
+    {
+        return true;
+    }
+
+    return false;
+}
+
+
+template<class Gt, class Vb>
+inline bool CGAL::indexedVertex<Gt, Vb>::ppSlave() const
+{
+    if (type_ >= 0 && type_ < index_)
+    {
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template<class Gt, class Vb>
+inline bool CGAL::indexedVertex<Gt, Vb>::internalOrBoundaryPoint() const
+{
+    return internalOrBoundaryPoint(index_, type_);
+}
+
+
+template<class Gt, class Vb>
+inline bool CGAL::indexedVertex<Gt, Vb>::internalOrBoundaryPoint
+(
+    int index,
+    int type
+)
+{
+    return internalPoint(type) || ppMaster(index, type);
+}
+
+
+template<class Gt, class Vb>
+inline bool CGAL::indexedVertex<Gt, Vb>::nearOrOnBoundary() const
+{
+    return pairPoint() || nearBoundary();
+}
+
+
+template<class Gt, class Vb>
+inline bool CGAL::indexedVertex<Gt, Vb>::anyInternalOrBoundaryPoint() const
+{
+    return internalOrBoundaryPoint() || referredInternalOrBoundaryPoint();
+}
+
+
+// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/pointFeatureEdgesTypes.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/pointFeatureEdgesTypes.H
new file mode 100644
index 0000000000..3009f385b2
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/pointFeatureEdgesTypes.H
@@ -0,0 +1,87 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+    As a special exception, you have permission to link this program with the
+    CGAL library and distribute executables, as long as you follow the
+    requirements of the GNU GPL in regard to all of the software in the
+    executable aside from CGAL.
+
+Class
+    pointFeatureEdgesTypes
+
+Description
+    struct for holding information on the types of feature edges attached to
+    feature points
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef pointFeatureEdgesTypes_H
+#define pointFeatureEdgesTypes_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class pointFeatureEdgesTypes Declaration
+\*---------------------------------------------------------------------------*/
+
+//- Hold the types of feature edges attached to the point.
+struct pointFeatureEdgesTypes
+{
+    label ptI;
+    label nExternal, nInternal;
+    label nFlat, nOpen, nMultiple, nNonFeature;
+
+    pointFeatureEdgesTypes(const label pointLabel)
+    {
+        ptI = pointLabel;
+        nExternal = nInternal = 0;
+        nFlat = nOpen = nMultiple = nNonFeature = 0;
+    }
+
+    friend Ostream& operator<<(Ostream& os, const pointFeatureEdgesTypes& p)
+    {
+        os  << "E " << "I " << "F " << "O " << "M " << "N "
+            << "Pt" << nl
+            << p.nExternal << " " << p.nInternal << " "
+            << p.nFlat << " " << p.nOpen << " "
+            << p.nMultiple << " " << p.nNonFeature << " "
+            << p.ptI
+            << endl;
+
+        return os;
+    }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.C
index ce1994e7bd..1058b6f133 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.C
@@ -30,11 +30,9 @@ License
 
 Foam::cvControls::cvControls
 (
-    const conformalVoronoiMesh& cvMesh,
     const dictionary& cvMeshDict
 )
 :
-    cvMesh_(cvMesh),
     cvMeshDict_(cvMeshDict)
 {
     // Surface conformation controls
@@ -61,6 +59,11 @@ Foam::cvControls::cvControls
         surfDict.lookup("featureEdgeExclusionDistanceCoeff")
     );
 
+    specialiseFeaturePoints_ = Switch
+    (
+        surfDict.lookupOrDefault<Switch>("specialiseFeaturePoints", false)
+    );
+
     surfaceSearchDistanceCoeff_ = readScalar
     (
         surfDict.lookup("surfaceSearchDistanceCoeff")
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.H
index 8b00b9a163..28e92c05bd 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.H
@@ -54,9 +54,6 @@ class cvControls
 {
     // Private data
 
-        //- Reference to the conformalVoronoiMesh holding this cvControls object
-        const conformalVoronoiMesh& cvMesh_;
-
         //- Reference to the cvMeshDict
         const dictionary& cvMeshDict_;
 
@@ -80,6 +77,9 @@ class cvControls
             //  fraction of the local target cell size
             scalar featureEdgeExclusionDistanceCoeff_;
 
+            //- Switch for using specialised feature points
+            Switch specialiseFeaturePoints_;
+
             //- Surface search distance coefficient - fraction of the local
             //  target cell size
             scalar surfaceSearchDistanceCoeff_;
@@ -291,7 +291,6 @@ public:
         //- Construct from references to conformalVoronoiMesh and dictionary
         cvControls
         (
-            const conformalVoronoiMesh& cvMesh,
             const dictionary& cvMeshDict
         );
 
@@ -319,6 +318,9 @@ public:
             //- Return the featureEdgeExclusionDistanceCoeff
             inline scalar featureEdgeExclusionDistanceCoeff() const;
 
+            //- Return whether to use specialised feature points
+            inline Switch specialiseFeaturePoints() const;
+
             //- Return the surfaceSearchDistanceCoeff
             inline scalar surfaceSearchDistanceCoeff() const;
 
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControlsI.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControlsI.H
index 5fd3364019..53e79ca1cd 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControlsI.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControlsI.H
@@ -54,6 +54,10 @@ inline Foam::scalar Foam::cvControls::featureEdgeExclusionDistanceCoeff() const
     return featureEdgeExclusionDistanceCoeff_;
 }
 
+inline Foam::Switch Foam::cvControls::specialiseFeaturePoints() const
+{
+    return specialiseFeaturePoints_;
+}
 
 inline Foam::scalar Foam::cvControls::surfaceSearchDistanceCoeff() const
 {
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/faceAreaWeightModel/faceAreaWeightModel.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/faceAreaWeightModel/faceAreaWeightModel.C
index ec7a8ed29d..76d68aa073 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/faceAreaWeightModel/faceAreaWeightModel.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/faceAreaWeightModel/faceAreaWeightModel.C
@@ -42,12 +42,10 @@ defineRunTimeSelectionTable(faceAreaWeightModel, dictionary);
 faceAreaWeightModel::faceAreaWeightModel
 (
     const word& type,
-    const dictionary& relaxationDict,
-    const conformalVoronoiMesh& cvMesh
+    const dictionary& relaxationDict
 )
 :
     dictionary(relaxationDict),
-    cvMesh_(cvMesh),
     coeffDict_(subDict(type + "Coeffs"))
 {}
 
@@ -56,8 +54,7 @@ faceAreaWeightModel::faceAreaWeightModel
 
 autoPtr<faceAreaWeightModel> faceAreaWeightModel::New
 (
-    const dictionary& relaxationDict,
-    const conformalVoronoiMesh& cvMesh
+    const dictionary& relaxationDict
 )
 {
     word faceAreaWeightModelTypeName
@@ -85,7 +82,7 @@ autoPtr<faceAreaWeightModel> faceAreaWeightModel::New
             << exit(FatalError);
     }
 
-    return autoPtr<faceAreaWeightModel>(cstrIter()(relaxationDict, cvMesh));
+    return autoPtr<faceAreaWeightModel>(cstrIter()(relaxationDict));
 }
 
 
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/faceAreaWeightModel/faceAreaWeightModel.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/faceAreaWeightModel/faceAreaWeightModel.H
index 0eaf7b14f6..3e88f5bb74 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/faceAreaWeightModel/faceAreaWeightModel.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/faceAreaWeightModel/faceAreaWeightModel.H
@@ -39,7 +39,6 @@ SourceFiles
 #define faceAreaWeightModel_H
 
 #include "point.H"
-#include "conformalVoronoiMesh.H"
 #include "dictionary.H"
 #include "autoPtr.H"
 #include "runTimeSelectionTables.H"
@@ -62,9 +61,6 @@ protected:
 
     // Protected data
 
-        //- Reference to the conformalVoronoiMesh holding this cvControls object
-        const conformalVoronoiMesh& cvMesh_;
-
         //- Method coeffs dictionary
         dictionary coeffDict_;
 
@@ -93,10 +89,9 @@ public:
             faceAreaWeightModel,
             dictionary,
             (
-                const dictionary& faceAreaWeightDict,
-                const conformalVoronoiMesh& cvMesh
+                const dictionary& faceAreaWeightDict
             ),
-            (faceAreaWeightDict, cvMesh)
+            (faceAreaWeightDict)
         );
 
 
@@ -106,8 +101,7 @@ public:
         faceAreaWeightModel
         (
             const word& type,
-            const dictionary& faceAreaWeightDict,
-            const conformalVoronoiMesh& cvMesh
+            const dictionary& faceAreaWeightDict
         );
 
 
@@ -116,8 +110,7 @@ public:
         //- Return a reference to the selected faceAreaWeightModel
         static autoPtr<faceAreaWeightModel> New
         (
-            const dictionary& faceAreaWeightDict,
-            const conformalVoronoiMesh& cvMesh
+            const dictionary& faceAreaWeightDict
         );
 
 
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/piecewiseLinearRamp/piecewiseLinearRamp.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/piecewiseLinearRamp/piecewiseLinearRamp.C
index 2f7df3d353..0fa6025c99 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/piecewiseLinearRamp/piecewiseLinearRamp.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/piecewiseLinearRamp/piecewiseLinearRamp.C
@@ -45,11 +45,10 @@ addToRunTimeSelectionTable
 
 piecewiseLinearRamp::piecewiseLinearRamp
 (
-    const dictionary& faceAreaWeightDict,
-    const conformalVoronoiMesh& cvMesh
+    const dictionary& faceAreaWeightDict
 )
 :
-    faceAreaWeightModel(typeName, faceAreaWeightDict, cvMesh),
+    faceAreaWeightModel(typeName, faceAreaWeightDict),
     lAF_(readScalar(coeffDict().lookup("lowerAreaFraction"))),
     uAF_(readScalar(coeffDict().lookup("upperAreaFraction")))
 {}
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/piecewiseLinearRamp/piecewiseLinearRamp.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/piecewiseLinearRamp/piecewiseLinearRamp.H
index 0c0ba587f3..14f6a929bb 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/piecewiseLinearRamp/piecewiseLinearRamp.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/piecewiseLinearRamp/piecewiseLinearRamp.H
@@ -71,8 +71,7 @@ public:
         //- Construct from components
         piecewiseLinearRamp
         (
-            const dictionary& faceAreaWeightDict,
-            const conformalVoronoiMesh& cvMesh
+            const dictionary& faceAreaWeightDict
         );
 
 
diff --git a/applications/utilities/mesh/generation/cvMesh/cvMesh.C b/applications/utilities/mesh/generation/cvMesh/cvMesh.C
index fa95c852c1..daf03aa376 100644
--- a/applications/utilities/mesh/generation/cvMesh/cvMesh.C
+++ b/applications/utilities/mesh/generation/cvMesh/cvMesh.C
@@ -59,6 +59,8 @@ int main(int argc, char *argv[])
         )
     );
 
+    conformalVoronoiMesh::debug = true;
+
     conformalVoronoiMesh mesh(runTime, cvMeshDict);
 
     while (runTime.loop())
@@ -72,7 +74,7 @@ int main(int argc, char *argv[])
             << nl << endl;
     }
 
-    mesh.writeMesh(runTime.constant());
+    mesh.writeMesh(runTime.constant(), true);
 
     Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
         << "  ClockTime = " << runTime.elapsedClockTime() << " s"