diff --git a/applications/utilities/mesh/generation/cvMesh/Make/options b/applications/utilities/mesh/generation/cvMesh/Make/options
index 529425600e..8a3f321499 100644
--- a/applications/utilities/mesh/generation/cvMesh/Make/options
+++ b/applications/utilities/mesh/generation/cvMesh/Make/options
@@ -26,5 +26,3 @@ EXE_LIBS = \
     -ldynamicMesh \
     -lmpfr \
     -lboost_thread
-
-
diff --git a/src/mesh/conformalVoronoiMesh/Make/files b/src/mesh/conformalVoronoiMesh/Make/files
index 8a6f6a2a3c..9f511dbff7 100644
--- a/src/mesh/conformalVoronoiMesh/Make/files
+++ b/src/mesh/conformalVoronoiMesh/Make/files
@@ -1,6 +1,8 @@
 #include CGAL_FILES
 
 conformalVoronoiMesh/conformalVoronoiMesh.C
+conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
+conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
 conformalVoronoiMesh/conformalVoronoiMeshIO.C
 
 cvControls/cvControls.C
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index 3c85a22075..801a1b07b4 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@@ -954,1405 +954,6 @@ void Foam::conformalVoronoiMesh::setVertexSizeAndAlignment()
 }
 
 
-void Foam::conformalVoronoiMesh::conformToSurface()
-{
-    reconformationMode reconfMode = reconformationControl();
-
-    if (reconfMode == rmNone)
-    {
-        // Reinsert stored surface conformation
-        reinsertSurfaceConformation();
-    }
-    else
-    {
-        // Rebuild, insert and store new surface conformation
-        buildSurfaceConformation(reconfMode);
-    }
-}
-
-Foam::conformalVoronoiMesh::reconformationMode
-Foam::conformalVoronoiMesh::reconformationControl() const
-{
-    if (!runTime_.run())
-    {
-        Info<< nl << "    Rebuilding surface conformation for final output"
-            << endl;
-
-        return rmFine;
-    }
-    else if
-    (
-        runTime_.timeIndex()
-      % cvMeshControls().surfaceConformationRebuildFrequency() == 0
-    )
-    {
-        Info<< nl << "    Rebuilding surface conformation for more iterations"
-            << endl;
-
-        return rmCoarse;
-    }
-
-    return rmNone;
-}
-
-void Foam::conformalVoronoiMesh::buildSurfaceConformation
-(
-    reconformationMode reconfMode
-)
-{
-    if (reconfMode == rmCoarse)
-    {
-        Info<< nl << "    Build coarse surface conformation" << endl;
-    }
-    else if (reconfMode == rmFine)
-    {
-        Info<< nl << "    Build fine surface conformation" << endl;
-    }
-    else if (reconfMode == rmNone)
-    {
-        WarningIn("buildSurfaceConformation(reconformationMode reconfMode)")
-            << "reconformationMode rmNone specified, not building conformation"
-            << endl;
-
-        return;
-    }
-    else
-    {
-        WarningIn("buildSurfaceConformation(reconformationMode reconfMode)")
-            << "Unknown reconformationMode " << reconfMode << endl;
-
-        return;
-    }
-
-    timeCheck();
-
-    startOfSurfacePointPairs_ = number_of_vertices();
-
-    // Initialise containers to store the edge conformation locations
-    DynamicList<point> newEdgeLocations;
-
-    pointField existingEdgeLocations(0);
-
-    autoPtr<indexedOctree<treeDataPoint> > edgeLocationTree;
-
-    // Initialise the edgeLocationTree
-    buildEdgeLocationTree(edgeLocationTree, existingEdgeLocations);
-
-    label initialTotalHits = 0;
-
-    // Initial surface protrusion conformation - nearest surface point
-    {
-        Info<< "    EDGE DISTANCE COEFFS HARD-CODED." << endl;
-        scalar edgeSearchDistCoeffSqr = sqr(1.1);
-        scalar surfacePtReplaceDistCoeffSqr = sqr(0.5);
-
-        DynamicList<pointIndexHit> surfaceHits;
-        DynamicList<label> hitSurfaces;
-
-        DynamicList<pointIndexHit> featureEdgeHits;
-        DynamicList<label> featureEdgeFeaturesHit;
-
-        for
-        (
-            Triangulation::Finite_vertices_iterator vit =
-            finite_vertices_begin();
-            vit != finite_vertices_end();
-            vit++
-        )
-        {
-            if (vit->internalPoint())
-            {
-                point vert(topoint(vit->point()));
-                scalar searchDistanceSqr = surfaceSearchDistanceSqr(vert);
-                pointIndexHit surfHit;
-                label hitSurface;
-
-                geometryToConformTo_.findSurfaceNearest
-                (
-                    vert,
-                    searchDistanceSqr,
-                    surfHit,
-                    hitSurface
-                );
-
-                if (surfHit.hit())
-                {
-                    vit->setNearBoundary();
-
-                    if (dualCellSurfaceAnyIntersection(vit))
-                    {
-                        addSurfaceAndEdgeHits
-                        (
-                            vit,
-                            vert,
-                            surfHit,
-                            hitSurface,
-                            surfacePtReplaceDistCoeffSqr,
-                            edgeSearchDistCoeffSqr,
-                            surfaceHits,
-                            hitSurfaces,
-                            featureEdgeHits,
-                            featureEdgeFeaturesHit,
-                            newEdgeLocations,
-                            existingEdgeLocations,
-                            edgeLocationTree
-                        );
-                    }
-                }
-            }
-        }
-
-        Info<< nl <<"    Initial conformation" << nl
-            << "        Number of vertices " << number_of_vertices() << nl
-            << "        Number of surface hits " << surfaceHits.size() << nl
-            << "        Number of edge hits " << featureEdgeHits.size()
-            << endl;
-
-        insertSurfacePointPairs
-        (
-            surfaceHits,
-            hitSurfaces,
-            "surfaceConformationLocations_initial.obj"
-        );
-
-        insertEdgePointGroups
-        (
-            featureEdgeHits,
-            featureEdgeFeaturesHit,
-            "edgeConformationLocations_initial.obj"
-        );
-
-        initialTotalHits = surfaceHits.size() + featureEdgeHits.size();
-    }
-
-    label iterationNo = 0;
-
-    label maxIterations = 10;
-
-    Info<< nl << "    MAX ITERATIONS HARD CODED TO "<< maxIterations << endl;
-
-    scalar iterationToIntialHitRatioLimit = 0.01;
-
-    label hitLimit = label(iterationToIntialHitRatioLimit*initialTotalHits);
-
-    Info<< "    STOPPING ITERATIONS WHEN TOTAL NUMBER OF HITS DROPS BELOW "
-        << iterationToIntialHitRatioLimit << " (HARD CODED) OF INITIAL HITS ("
-        << hitLimit << ")"
-        << endl;
-
-    // Set totalHits to a large enough positive value to enter the while loop on
-    // the first iteration
-    label totalHits = initialTotalHits;
-
-    while
-    (
-        totalHits > 0
-     && totalHits > hitLimit
-     && iterationNo < maxIterations
-    )
-    {
-        Info<< "    EDGE DISTANCE COEFFS HARD-CODED." << endl;
-        scalar edgeSearchDistCoeffSqr = sqr(1.25);
-        scalar surfacePtReplaceDistCoeffSqr = sqr(0.7);
-
-        DynamicList<pointIndexHit> surfaceHits;
-        DynamicList<label> hitSurfaces;
-
-        DynamicList<pointIndexHit> featureEdgeHits;
-        DynamicList<label> featureEdgeFeaturesHit;
-
-        for
-        (
-            Triangulation::Finite_vertices_iterator vit =
-            finite_vertices_begin();
-            vit != finite_vertices_end();
-            vit++
-        )
-        {
-            // The initial surface conformation has already identified the
-            // nearBoundary set of vertices.  Previously inserted boundary
-            // points can also generate protrusions and must be assessed too.
-
-            if (vit->nearBoundary() || vit->ppMaster())
-            {
-                point vert(topoint(vit->point()));
-                pointIndexHit surfHit;
-                label hitSurface;
-
-                dualCellLargestSurfaceProtrusion(vit, surfHit, hitSurface);
-
-                if (surfHit.hit())
-                {
-                    addSurfaceAndEdgeHits
-                    (
-                        vit,
-                        vert,
-                        surfHit,
-                        hitSurface,
-                        surfacePtReplaceDistCoeffSqr,
-                        edgeSearchDistCoeffSqr,
-                        surfaceHits,
-                        hitSurfaces,
-                        featureEdgeHits,
-                        featureEdgeFeaturesHit,
-                        newEdgeLocations,
-                        existingEdgeLocations,
-                        edgeLocationTree
-                    );
-                }
-            }
-        }
-
-        Info<< nl <<"    Conformation iteration " << iterationNo << nl
-            << "        Number of vertices " << number_of_vertices() << nl
-            << "        Number of surface hits " << surfaceHits.size() << nl
-            << "        Number of edge hits " << featureEdgeHits.size()
-            << endl;
-
-        totalHits = surfaceHits.size() + featureEdgeHits.size();
-
-        if (totalHits > 0)
-        {
-            insertSurfacePointPairs
-            (
-                surfaceHits,
-                hitSurfaces,
-                fileName
-                (
-                    "surfaceConformationLocations_" + name(iterationNo) + ".obj"
-                )
-            );
-
-            insertEdgePointGroups
-            (
-                featureEdgeHits,
-                featureEdgeFeaturesHit,
-                "edgeConformationLocations_" + name(iterationNo) + ".obj"
-            );
-        }
-
-        iterationNo++;
-
-        if (iterationNo == maxIterations)
-        {
-            WarningIn("conformalVoronoiMesh::conformToSurface()")
-                << "Maximum surface conformation iterations ("
-                << maxIterations <<  ") reached." << endl;
-        }
-
-        if (totalHits < hitLimit)
-        {
-            Info<< nl << "    Total hits (" << totalHits
-                << ") less than limit (" << hitLimit
-                << "), stopping iterations" << endl;
-        }
-    }
-
-    // Info<< nl << "    After iterations, check penetrations" << endl;
-
-    // for
-    // (
-    //     Triangulation::Finite_vertices_iterator vit =
-    //     finite_vertices_begin();
-    //     vit != finite_vertices_end();
-    //     vit++
-    // )
-    // {
-    //     if (vit->internalOrBoundaryPoint())
-    //     {
-    //         point vert(topoint(vit->point()));
-    //         pointIndexHit surfHit;
-    //         label hitSurface;
-
-    //         dualCellLargestSurfaceProtrusion(vit, surfHit, hitSurface);
-
-    //         if (surfHit.hit())
-    //         {
-    //             Info<< nl << "Residual penetration: " << nl
-    //                 << vit->index() << nl
-    //                 << vit->type() << nl
-    //                 << vit->ppMaster() << nl
-    //                 << "nearFeaturePt "
-    //                 << nearFeaturePt(surfHit.hitPoint()) << nl
-    //                 << vert << nl
-    //                 << surfHit.hitPoint()
-    //                 << endl;
-    //         }
-    //     }
-    // }
-
-    storeSurfaceConformation();
-}
-
-
-bool Foam::conformalVoronoiMesh::dualCellSurfaceAnyIntersection
-(
-    const Triangulation::Finite_vertices_iterator& vit
-) const
-{
-    std::list<Facet> facets;
-    incident_facets(vit, std::back_inserter(facets));
-
-    for
-    (
-        std::list<Facet>::iterator fit=facets.begin();
-        fit != facets.end();
-        ++fit
-    )
-    {
-        if
-        (
-            is_infinite(fit->first)
-         || is_infinite(fit->first->neighbor(fit->second))
-        )
-        {
-            return true;
-        }
-
-        point dE0 = topoint(dual(fit->first));
-
-        // If edge end is outside bounding box then edge cuts boundary
-        if (!geometryToConformTo_.bounds().contains(dE0))
-        {
-            return true;
-        }
-
-        point dE1 = topoint(dual(fit->first->neighbor(fit->second)));
-
-        // If other edge end is outside bounding box then edge cuts boundary
-        if (!geometryToConformTo_.bounds().contains(dE1))
-        {
-            return true;
-        }
-
-        // Check for the edge passing through a surface
-        if (geometryToConformTo_.findSurfaceAnyIntersection(dE0, dE1))
-        {
-            return true;
-        }
-    }
-
-    return false;
-}
-
-
-void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
-(
-    const Triangulation::Finite_vertices_iterator& vit,
-    pointIndexHit& surfHitLargest,
-    label& hitSurfaceLargest
-) const
-{
-    std::list<Facet> facets;
-    incident_facets(vit, std::back_inserter(facets));
-
-    point vert(topoint(vit->point()));
-
-    scalar maxProtrusionDistance = maxSurfaceProtrusion(vert);
-
-    for
-    (
-        std::list<Facet>::iterator fit=facets.begin();
-        fit != facets.end();
-        ++fit
-    )
-    {
-        if
-        (
-            !is_infinite(fit->first)
-         && !is_infinite(fit->first->neighbor(fit->second))
-        )
-        {
-            point edgeMid =
-                0.5
-               *(
-                    topoint(dual(fit->first))
-                  + topoint(dual(fit->first->neighbor(fit->second)))
-                );
-
-            pointIndexHit surfHit;
-            label hitSurface;
-
-            geometryToConformTo_.findSurfaceAnyIntersection
-            (
-                vert,
-                edgeMid,
-                surfHit,
-                hitSurface
-            );
-
-            if (surfHit.hit())
-            {
-                vectorField norm(1);
-
-                allGeometry_[hitSurface].getNormal
-                (
-                    List<pointIndexHit>(1, surfHit),
-                    norm
-                );
-
-                const vector& n = norm[0];
-
-                scalar normalProtrusionDistance =
-                    (edgeMid - surfHit.hitPoint()) & n;
-
-                if (normalProtrusionDistance > maxProtrusionDistance)
-                {
-                    surfHitLargest = surfHit;
-                    hitSurfaceLargest = hitSurface;
-
-                    maxProtrusionDistance = normalProtrusionDistance;
-                }
-            }
-        }
-    }
-}
-
-
-void Foam::conformalVoronoiMesh::limitDisplacement
-(
-    const Triangulation::Finite_vertices_iterator& vit,
-    vector& displacement
-) const
-{
-    point pt = topoint(vit->point());
-    point dispPt = pt + displacement;
-
-    bool limit = false;
-
-    pointIndexHit surfHit;
-    label hitSurface;
-
-    if (!geometryToConformTo_.bounds().contains(dispPt))
-    {
-        // If dispPt is outside bounding box then displacement cuts boundary
-        limit = true;
-
-        // Info<< "    bb limit" << endl;
-    }
-    else if (geometryToConformTo_.findSurfaceAnyIntersection(pt, dispPt))
-    {
-        // Full surface penetration test
-        limit = true;
-
-        // Info<< "    intersection limit" << endl;
-    }
-    else
-    {
-        // Testing if the displaced position is too close to the surface.
-        // Within twice the local surface point pair insertion distance is
-        // considered "too close"
-
-        scalar searchDistanceSqr = sqr
-        (
-            2*vit->targetCellSize()
-           *cvMeshControls().pointPairDistanceCoeff()
-        );
-
-        geometryToConformTo_.findSurfaceNearest
-        (
-            dispPt,
-            searchDistanceSqr,
-            surfHit,
-            hitSurface
-        );
-
-        if (surfHit.hit())
-        {
-            // Info<< "    proximity limit" << endl;
-
-            limit = true;
-
-            if (magSqr(pt - surfHit.hitPoint()) <= searchDistanceSqr)
-            {
-                // Info<< "    Cannot limit displacement, point " << pt
-                //     << " closer than tolerance" << endl;
-
-                return;
-            }
-        }
-    }
-
-    if (limit)
-    {
-        // Halve the displacement and call this function again.  Will continue
-        // recursively until the displacement is small enough.
-
-        displacement *= 0.5;
-
-        // Info<< "    Limiting displacement of point " << pt << endl;
-
-        limitDisplacement(vit, displacement);
-    }
-}
-
-
-bool Foam::conformalVoronoiMesh::nearFeatureEdgeLocation
-(
-    const point& pt,
-    DynamicList<point>& newEdgeLocations,
-    pointField& existingEdgeLocations,
-    autoPtr<indexedOctree<treeDataPoint> >& edgeLocationTree
-) const
-{
-    scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
-
-    // 0.01 and 1000 determined from speed tests, varying the indexedOctree
-    // rebuild frequency and balance of additions to queries.
-
-    if
-    (
-        newEdgeLocations.size()
-     >= max(0.01*existingEdgeLocations.size(), 1000)
-    )
-    {
-        existingEdgeLocations.append(newEdgeLocations);
-
-        buildEdgeLocationTree(edgeLocationTree, existingEdgeLocations);
-
-        newEdgeLocations.clear();
-    }
-    else
-    {
-        // Search for the nearest point in newEdgeLocations.
-        // Searching here first, because the intention is that the value of
-        // newEdgeLocationsSizeLimit should make this faster by design.
-
-        if (min(magSqr(newEdgeLocations - pt)) <= exclusionRangeSqr)
-        {
-            return true;
-        }
-    }
-
-    // Searching for the nearest point in existingEdgeLocations using the
-    // indexedOctree
-
-    pointIndexHit info = edgeLocationTree().findNearest(pt, exclusionRangeSqr);
-
-    return info.hit();
-}
-
-
-void Foam::conformalVoronoiMesh::buildEdgeLocationTree
-(
-    autoPtr<indexedOctree<treeDataPoint> >& edgeLocationTree,
-    const pointField& existingEdgeLocations
-) const
-{
-    treeBoundBox overallBb(geometryToConformTo_.bounds());
-
-    Random rndGen(72953);
-
-    overallBb.extend(rndGen, 1E-4);
-    overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-    overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-
-    edgeLocationTree.reset
-    (
-        new indexedOctree<treeDataPoint>
-        (
-            treeDataPoint(existingEdgeLocations),
-            overallBb,  // overall search domain
-            10,         // max levels
-            10.0,       // maximum ratio of cubes v.s. cells
-            100.0       // max. duplicity; n/a since no bounding boxes.
-        )
-    );
-}
-
-
-void Foam::conformalVoronoiMesh::buildSizeAndAlignmentTree() const
-{
-    treeBoundBox overallBb(geometryToConformTo_.bounds());
-
-    Random rndGen(627391);
-
-    overallBb.extend(rndGen, 1E-4);
-    overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-    overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-
-    sizeAndAlignmentTree_.reset
-    (
-        new indexedOctree<treeDataPoint>
-        (
-            treeDataPoint(sizeAndAlignmentLocations_),
-            overallBb,  // overall search domain
-            10,         // max levels
-            10.0,       // maximum ratio of cubes v.s. cells
-            100.0       // max. duplicity; n/a since no bounding boxes.
-        )
-    );
-}
-
-
-void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
-(
-    const Triangulation::Finite_vertices_iterator& vit,
-    const point& vert,
-    const pointIndexHit& surfHit,
-    label hitSurface,
-    scalar surfacePtReplaceDistCoeffSqr,
-    scalar edgeSearchDistCoeffSqr,
-    DynamicList<pointIndexHit>& surfaceHits,
-    DynamicList<label>& hitSurfaces,
-    DynamicList<pointIndexHit>& featureEdgeHits,
-    DynamicList<label>& featureEdgeFeaturesHit,
-    DynamicList<point>& newEdgeLocations,
-    pointField& existingEdgeLocations,
-    autoPtr<indexedOctree<treeDataPoint> >& edgeLocationTree
-) const
-{
-    bool keepSurfacePoint = true;
-
-    if (nearFeaturePt(surfHit.hitPoint()))
-    {
-        keepSurfacePoint = false;
-
-        if (vit->index() < startOfInternalPoints_)
-        {
-            surfaceHits.append(surfHit);
-
-            hitSurfaces.append(hitSurface);
-        }
-    }
-
-    List<pointIndexHit> edHits;
-
-    labelList featuresHit;
-
-    scalar targetCellSizeSqr = sqr(targetCellSize(vert));
-
-    geometryToConformTo_.findEdgeNearestByType
-    (
-        surfHit.hitPoint(),
-        edgeSearchDistCoeffSqr*targetCellSizeSqr,
-        edHits,
-        featuresHit
-    );
-
-    // Gather edge locations but do not add them to newEdgeLocations inside the
-    // loop as they will prevent nearby edge locations of different types being
-    // conformed to.
-
-    DynamicList<point> currentEdgeLocations;
-
-    forAll(edHits, i)
-    {
-        const pointIndexHit& edHit(edHits[i]);
-
-        label featureHit = featuresHit[i];
-
-        if (edHit.hit())
-        {
-            if (!nearFeaturePt(edHit.hitPoint()))
-            {
-                if
-                (
-                    magSqr(edHit.hitPoint() - surfHit.hitPoint())
-                  < surfacePtReplaceDistCoeffSqr*targetCellSizeSqr
-                )
-                {
-                    // If the point is within a given distance of a feature
-                    // edge, give control to edge control points instead, this
-                    // will prevent "pits" forming.
-
-                    keepSurfacePoint = false;
-                }
-
-                if
-                (
-                    !nearFeatureEdgeLocation
-                    (
-                        edHit.hitPoint(),
-                        newEdgeLocations,
-                        existingEdgeLocations,
-                        edgeLocationTree
-                    )
-                )
-                {
-                    // Do not place edge control points too close to a feature
-                    // point or existing edge control points
-
-                    featureEdgeHits.append(edHit);
-                    featureEdgeFeaturesHit.append(featureHit);
-
-                    currentEdgeLocations.append(edHit.hitPoint());
-                }
-            }
-        }
-    }
-
-    newEdgeLocations.append(currentEdgeLocations);
-
-    if (keepSurfacePoint)
-    {
-        surfaceHits.append(surfHit);
-
-        hitSurfaces.append(hitSurface);
-    }
-}
-
-
-void Foam::conformalVoronoiMesh::storeSurfaceConformation()
-{
-    Info<< nl << "    Storing surface conformation" << endl;
-
-    surfaceConformationVertices_.setSize
-    (
-        number_of_vertices() - startOfSurfacePointPairs_
-    );
-
-    label surfPtI = 0;
-
-    for
-    (
-        Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
-        vit != finite_vertices_end();
-        vit++
-    )
-    {
-        if (vit->index() >= startOfSurfacePointPairs_)
-        {
-            if (!vit->pairPoint())
-            {
-                FatalErrorIn("storeSurfaceConformation()")
-                    << "Trying to store a vertex that is not a surface point"
-                    << exit(FatalError);
-            }
-
-            surfaceConformationVertices_[surfPtI] = Vb(vit->point());
-
-            surfaceConformationVertices_[surfPtI].index() =
-            vit->index() - startOfSurfacePointPairs_;
-
-            surfaceConformationVertices_[surfPtI].type() =
-            vit->type() - startOfSurfacePointPairs_;
-
-            surfPtI++;
-        }
-    }
-
-    Info<< "    Stored " << surfaceConformationVertices_.size()
-        << " vertices" << endl;
-}
-
-
-void Foam::conformalVoronoiMesh::reinsertSurfaceConformation()
-{
-    Info<< nl << "    Reinserting stored surface conformation" << endl;
-
-    startOfSurfacePointPairs_ = number_of_vertices();
-
-    forAll(surfaceConformationVertices_, v)
-    {
-        insertVb(surfaceConformationVertices_[v], startOfSurfacePointPairs_);
-    }
-
-    Info<< "    Reinserted " << number_of_vertices() - startOfSurfacePointPairs_
-        << " vertices" << endl;
-}
-
-
-void Foam::conformalVoronoiMesh::calcDualMesh
-(
-    pointField& points,
-    faceList& faces,
-    labelList& owner,
-    labelList& neighbour,
-    wordList& patchNames,
-    labelList& patchSizes,
-    labelList& patchStarts
-)
-{
-    timeCheck();
-
-    setVertexSizeAndAlignment();
-
-    timeCheck();
-
-    // ------>  OLD, FOR REFERENCE
-        // // ~~~~~~~~~~~ removing short edges by indexing dual vertices ~~~~~~~~~~~~~~
-
-        // for
-        // (
-        //     Triangulation::Finite_cells_iterator cit = finite_cells_begin();
-        //     cit != finite_cells_end();
-        //     ++cit
-        // )
-        // {
-        //     cit->cellIndex() = -1;
-        // }
-
-        // points.setSize(number_of_cells());
-
-        // // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-        // label dualVertI = 0;
-
-        // // Scanning by number of short (dual) edges (nSE) attached to the
-        // // circumcentre of each Delaunay tet.  A Delaunay tet may only have four
-        // // dual edges emanating from its circumcentre, assigning positions and
-        // // indices to those with 4 short edges attached first, then >= 3, then >= 2
-        // // etc.
-        // for (label nSE = 4; nSE >= 0; nSE--)
-        // {
-        //     Info<< nl << "Scanning for dual vertices with >= "
-        //         << nSE
-        //         << " short edges attached." << endl;
-
-        //     for
-        //     (
-        //         Triangulation::Finite_cells_iterator cit = finite_cells_begin();
-        //         cit != finite_cells_end();
-        //         ++cit
-        //     )
-        //     {
-        //         // If the Delaunay tet has an index already then it has either
-        //         // evaluated itself and taken action or has had its index dictated
-        //         // by a neighbouring tet with more short edges attached.
-
-        //         if (cit->cellIndex() == -1)
-        //         {
-        //             point dualVertex = topoint(dual(cit));
-
-        //             label shortEdges = 0;
-
-        //             List<bool> edgeIsShort(4, false);
-
-        //             List<bool> neighbourAlreadyIndexed(4, false);
-
-        //             // Loop over the four facets of the Delaunay tet
-        //             for (label f = 0; f < 4; f++)
-        //             {
-        //                 // Check that at least one of the vertices of the facet is
-        //                 // an internal or boundary point
-        //                 if
-        //                 (
-        //                     cit->vertex(vertex_triple_index(f, 0))->
-        //                         internalOrBoundaryPoint()
-        //                  || cit->vertex(vertex_triple_index(f, 1))->
-        //                         internalOrBoundaryPoint()
-        //                  || cit->vertex(vertex_triple_index(f, 2))->
-        //                         internalOrBoundaryPoint()
-        //                 )
-        //                 {
-        //                     point neighDualVertex;
-
-        //                     label cNI = cit->neighbor(f)->cellIndex();
-
-        //                     if (cNI == -1)
-        //                     {
-        //                         neighDualVertex = topoint(dual(cit->neighbor(f)));
-        //                     }
-        //                     else
-        //                     {
-        //                         neighDualVertex = points[cNI];
-        //                     }
-
-        //                     if
-        //                     (
-        //                         magSqr(dualVertex - neighDualVertex)
-        //                       < sqr
-        //                         (
-        //                             minimumEdgeLength
-        //                             (
-        //                                 0.5*(dualVertex + neighDualVertex)
-        //                             )
-        //                         )
-        //                     )
-        //                     {
-        //                         edgeIsShort[f] = true;
-
-        //                         if (cNI > -1)
-        //                         {
-        //                             neighbourAlreadyIndexed[f] = true;
-        //                         }
-
-        //                         shortEdges++;
-        //                     }
-        //                 }
-        //             }
-
-        //             if (nSE == 0 && shortEdges == 0)
-        //             {
-        //                 // Final iteration and no short edges are found, index
-        //                 // remaining dual vertices.
-
-        //                 if
-        //                 (
-        //                     cit->vertex(0)->internalOrBoundaryPoint()
-        //                  || cit->vertex(1)->internalOrBoundaryPoint()
-        //                  || cit->vertex(2)->internalOrBoundaryPoint()
-        //                  || cit->vertex(3)->internalOrBoundaryPoint()
-        //                 )
-        //                 {
-        //                     cit->cellIndex() = dualVertI;
-        //                     points[dualVertI] = dualVertex;
-        //                     dualVertI++;
-        //                 }
-        //             }
-        //             else if
-        //             (
-        //                 shortEdges >= nSE
-        //             )
-        //             {
-        //                 // Info<< neighbourAlreadyIndexed << ' '
-        //                 //     << edgeIsShort << endl;
-
-        //                 label numUnindexedNeighbours = 1;
-
-        //                 for (label f = 0; f < 4; f++)
-        //                 {
-        //                     if (edgeIsShort[f] && !neighbourAlreadyIndexed[f])
-        //                     {
-        //                         dualVertex += topoint(dual(cit->neighbor(f)));
-
-        //                         numUnindexedNeighbours++;
-        //                     }
-        //                 }
-
-        //                 dualVertex /= numUnindexedNeighbours;
-
-        //                 label nearestExistingIndex = -1;
-
-        //                 point nearestIndexedNeighbourPos = vector::zero;
-
-        //                 scalar minDistSqrToNearestIndexedNeighbour = VGREAT;
-
-        //                 for (label f = 0; f < 4; f++)
-        //                 {
-        //                     if (edgeIsShort[f] && neighbourAlreadyIndexed[f])
-        //                     {
-        //                         label cNI = cit->neighbor(f)->cellIndex();
-
-        //                         point indexedNeighbourPos = points[cNI];
-
-        //                         if
-        //                         (
-        //                             magSqr(indexedNeighbourPos - dualVertex)
-        //                           < minDistSqrToNearestIndexedNeighbour
-        //                         )
-        //                         {
-        //                             nearestExistingIndex = cNI;
-
-        //                             nearestIndexedNeighbourPos =
-        //                             indexedNeighbourPos;
-
-        //                             minDistSqrToNearestIndexedNeighbour =
-        //                             magSqr(indexedNeighbourPos - dualVertex);
-        //                         }
-        //                     }
-        //                 }
-
-        //                 if
-        //                 (
-        //                     nearestExistingIndex > -1
-        //                  && minDistSqrToNearestIndexedNeighbour
-        //                   < sqr
-        //                     (
-        //                         minimumEdgeLength
-        //                         (
-        //                             0.5*(nearestIndexedNeighbourPos + dualVertex)
-        //                         )
-        //                     )
-        //                 )
-        //                 {
-        //                     points[nearestExistingIndex] =
-        //                     0.5*(dualVertex + nearestIndexedNeighbourPos);
-
-        //                     for (label f = 0; f < 4; f++)
-        //                     {
-        //                         if (edgeIsShort[f] && !neighbourAlreadyIndexed[f])
-        //                         {
-        //                             cit->neighbor(f)->cellIndex() =
-        //                             nearestExistingIndex;
-        //                         }
-        //                     }
-
-        //                     cit->cellIndex() = nearestExistingIndex;
-        //                 }
-        //                 else
-        //                 {
-        //                     for (label f = 0; f < 4; f++)
-        //                     {
-        //                         if (edgeIsShort[f] && !neighbourAlreadyIndexed[f])
-        //                         {
-        //                             cit->neighbor(f)->cellIndex() = dualVertI;
-        //                         }
-        //                     }
-
-        //                     cit->cellIndex() = dualVertI;
-
-        //                     points[dualVertI] = dualVertex;
-
-        //                     dualVertI++;
-        //                 }
-        //             }
-        //         }
-        //     }
-        // }
-
-        // points.setSize(dualVertI);
-
-    // <------ OLD, FOR REFERENCE
-
-    // Dual cell indexing
-
-    // Assign an index to the Delaunay vertices which will be the dual cell
-    // index used for owner neighbour assignment.
-
-    // The indices of the points are reset *which **destroys** the point-pair
-    // matching*, so the type of each vertex are reset to avoid any ambiguity.
-
-    label dualCelli = 0;
-
-    for
-    (
-        Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
-        vit != finite_vertices_end();
-        ++vit
-    )
-    {
-        if (vit->internalOrBoundaryPoint())
-        {
-            vit->type() = Vb::ptInternalPoint;
-            vit->index() = dualCelli;
-            dualCelli++;
-        }
-        else
-        {
-            vit->type() = Vb::ptFarPoint;
-            vit->index() = -1;
-        }
-    }
-
-    // Dual face filtering
-
-    // Indexing Delaunay cells, which are Dual vertices
-
-    label dualVertI = 0;
-
-    points.setSize(number_of_cells());
-
-    for
-    (
-        Triangulation::Finite_cells_iterator cit = finite_cells_begin();
-        cit != finite_cells_end();
-        ++cit
-    )
-    {
-        if
-        (
-            cit->vertex(0)->internalOrBoundaryPoint()
-         || cit->vertex(1)->internalOrBoundaryPoint()
-         || cit->vertex(2)->internalOrBoundaryPoint()
-         || cit->vertex(3)->internalOrBoundaryPoint()
-        )
-        {
-            cit->cellIndex() = dualVertI;
-            points[dualVertI] = topoint(dual(cit));
-            dualVertI++;
-        }
-        else
-        {
-            cit->cellIndex() = -1;
-        }
-    }
-
-    points.setSize(dualVertI);
-
-    // Assess close points to be merged
-
-    {
-        Info<< nl << "    Merging close points" << endl;
-
-        label nPtsMerged = 0;
-
-        do
-        {
-            Map<label> dualPtIndexMap;
-
-            nPtsMerged = mergeCloseDualVertices(points, dualPtIndexMap);
-
-            Info<< "        Merged " << nPtsMerged << " points" << endl;
-
-            reindexDualVertices(dualPtIndexMap);
-
-        } while (nPtsMerged > 0);
-    }
-
-    // Smooth the surface of the mesh
-
-    Info<< nl << "    Smoothing surface" << endl;
-
-    label nSmoothedVertices = 0;
-
-    do
-    {
-        Map<label> dualPtIndexMap;
-
-        nSmoothedVertices = smoothSurfaceDualFaces(points, dualPtIndexMap);
-
-        Info<< "        Smoothed " << nSmoothedVertices
-            << " points (0 HARD CODED)"
-            << endl;
-
-        reindexDualVertices(dualPtIndexMap);
-
-    } while (nSmoothedVertices > 0);
-
-    {
-        Info<< nl << "    Merging close points" << endl;
-
-        label nPtsMerged = 0;
-
-        do
-        {
-            Map<label> dualPtIndexMap;
-
-            nPtsMerged = mergeCloseDualVertices(points, dualPtIndexMap);
-
-            Info<< "        Merged " << nPtsMerged << " points" << endl;
-
-            reindexDualVertices(dualPtIndexMap);
-
-        } while (nPtsMerged > 0);
-    }
-
-    // Assess faces for collapse
-
-    Info<< nl << "    Collapsing unnecessary faces" << endl;
-
-    label nCollapsedFaces = 0;
-
-    do
-    {
-        Map<label> dualPtIndexMap;
-
-        nCollapsedFaces = collapseFaces(points, dualPtIndexMap);
-
-        reindexDualVertices(dualPtIndexMap);
-
-        Info<< "        Collapsed " << nCollapsedFaces << " faces" << endl;
-
-    } while (nCollapsedFaces > 0);
-
-    // Final dual face and owner neighbour construction
-
-    createFacesOwnerNeighbourAndPatches
-    (
-        points,
-        faces,
-        owner,
-        neighbour,
-        patchNames,
-        patchSizes,
-        patchStarts,
-        false
-    );
-}
-
-
-void Foam::conformalVoronoiMesh::calcTetMesh
-(
-    pointField& points,
-    faceList& faces,
-    labelList& owner,
-    labelList& neighbour,
-    wordList& patchNames,
-    labelList& patchSizes,
-    labelList& patchStarts
-)
-{
-    labelList vertexMap(number_of_vertices());
-
-    label vertI = 0;
-
-    points.setSize(number_of_vertices());
-
-    for
-    (
-        Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
-        vit != finite_vertices_end();
-        ++vit
-    )
-    {
-        if (vit->internalPoint() || vit->pairPoint())
-        {
-            vertexMap[vit->index()] = vertI;
-            points[vertI] =  topoint(vit->point());
-            vertI++;
-        }
-    }
-
-    points.setSize(vertI);
-
-    label cellI = 0;
-
-    for
-    (
-        Triangulation::Finite_cells_iterator cit = finite_cells_begin();
-        cit != finite_cells_end();
-        ++cit
-    )
-    {
-        if
-        (
-            cit->vertex(0)->internalOrBoundaryPoint()
-         || cit->vertex(1)->internalOrBoundaryPoint()
-         || cit->vertex(2)->internalOrBoundaryPoint()
-         || cit->vertex(3)->internalOrBoundaryPoint()
-        )
-        {
-             cit->cellIndex() = cellI++;
-        }
-        else
-        {
-            cit->cellIndex() = -1;
-        }
-    }
-
-    patchNames = geometryToConformTo_.patchNames();
-
-    patchNames.setSize(patchNames.size() + 1);
-
-    patchNames[patchNames.size() - 1] = "cvMesh_defaultPatch";
-
-    label nPatches = patchNames.size();
-
-    List<DynamicList<face> > patchFaces(nPatches, DynamicList<face>(0));
-
-    List<DynamicList<label> > patchOwners(nPatches, DynamicList<label>(0));
-
-    faces.setSize(number_of_facets());
-
-    owner.setSize(number_of_facets());
-
-    neighbour.setSize(number_of_facets());
-
-    label faceI = 0;
-
-    labelList verticesOnTriFace(3, -1);
-
-    face newFace(verticesOnTriFace);
-
-    for
-    (
-        Triangulation::Finite_facets_iterator fit = finite_facets_begin();
-        fit != finite_facets_end();
-        ++fit
-    )
-    {
-        const Cell_handle c1(fit->first);
-        const int oppositeVertex = fit->second;
-        const Cell_handle c2(c1->neighbor(oppositeVertex));
-
-        label c1I = c1->cellIndex();
-        label c2I = c2->cellIndex();
-
-        label ownerCell = -1;
-        label neighbourCell = -1;
-
-        if (c1I == -1 && c2I == -1)
-        {
-            // Both tets are outside, skip
-            continue;
-        }
-
-        for (label i = 0; i < 3; i++)
-        {
-            verticesOnTriFace[i] = vertexMap
-            [
-                c1->vertex(vertex_triple_index(oppositeVertex, i))->index()
-            ];
-        }
-
-        newFace = face(verticesOnTriFace);
-
-        if (c1I == -1 || c2I == -1)
-        {
-            // Boundary face...
-            if (c1I == -1)
-            {
-                //... with c1 outside
-                ownerCell = c2I;
-            }
-            else
-            {
-                // ... with c2 outside
-                ownerCell = c1I;
-
-                reverse(newFace);
-            }
-
-            label patchIndex = geometryToConformTo_.findPatch
-            (
-                newFace.centre(points)
-            );
-
-            if (patchIndex == -1)
-            {
-                patchIndex = patchNames.size() - 1;
-
-                WarningIn("Foam::conformalVoronoiMesh::calcTetMesh")
-                    << "Tet face centre at  " << nl
-                    << "    " << newFace.centre(points) << nl
-                    << "    did not find a surface patch. Adding to "
-                    << patchNames[patchIndex]
-                    << endl;
-            }
-
-            patchFaces[patchIndex].append(newFace);
-            patchOwners[patchIndex].append(ownerCell);
-        }
-        else
-        {
-            // Internal face...
-            if (c1I < c2I)
-            {
-                // ...with c1 as the ownerCell
-                ownerCell = c1I;
-                neighbourCell = c2I;
-
-                reverse(newFace);
-            }
-            else
-            {
-                // ...with c2 as the ownerCell
-                ownerCell = c2I;
-                neighbourCell = c1I;
-            }
-
-            faces[faceI] = newFace;
-            owner[faceI] = ownerCell;
-            neighbour[faceI] = neighbourCell;
-            faceI++;
-        }
-    }
-
-    label nInternalFaces = faceI;
-
-    faces.setSize(nInternalFaces);
-    owner.setSize(nInternalFaces);
-    neighbour.setSize(nInternalFaces);
-
-    sortFaces(faces, owner, neighbour);
-
-    addPatches
-    (
-        nInternalFaces,
-        faces,
-        owner,
-        patchNames,
-        patchSizes,
-        patchStarts,
-        patchFaces,
-        patchOwners,
-        false
-    );
-}
-
-
 Foam::face Foam::conformalVoronoiMesh::buildDualFace
 (
     const Triangulation::Finite_edges_iterator& eit
@@ -2404,955 +1005,6 @@ Foam::face Foam::conformalVoronoiMesh::buildDualFace
 }
 
 
-Foam::label Foam::conformalVoronoiMesh::mergeCloseDualVertices
-(
-    const pointField& pts,
-    Map<label>& dualPtIndexMap
-) const
-{
-    label nPtsMerged = 0;
-
-    scalar closenessTolerance = 1e-4;
-
-    for
-    (
-        Triangulation::Finite_facets_iterator fit = finite_facets_begin();
-        fit != finite_facets_end();
-        ++fit
-    )
-    {
-        const Cell_handle c1(fit->first);
-        const int oppositeVertex = fit->second;
-        const Cell_handle c2(c1->neighbor(oppositeVertex));
-
-        label& c1I = c1->cellIndex();
-        label& c2I = c2->cellIndex();
-
-        if (dualPtIndexMap.found(c1I) || dualPtIndexMap.found(c2I))
-        {
-            // One of the points of this edge has already been
-            // merged this sweep, leave for next sweep
-
-            continue;
-        }
-
-        if(c1I != -1 && c2I != -1 && (c1I != c2I))
-        {
-            if
-            (
-                magSqr(pts[c1I] - pts[c2I])
-              < sqr(averageAnyCellSize(fit)*closenessTolerance)
-            )
-            {
-                dualPtIndexMap.insert(c1I, c1I);
-                dualPtIndexMap.insert(c2I, c1I);
-
-                nPtsMerged++;
-            }
-        }
-    }
-
-    return nPtsMerged;
-}
-
-Foam::label Foam::conformalVoronoiMesh::smoothSurfaceDualFaces
-(
-    pointField& pts,
-    Map<label>& dualPtIndexMap
-) const
-{
-    label nSmoothedVertices = 0;
-
-    const scalar cosPerpendicularToleranceAngle = cos(degToRad(80));
-
-    for
-    (
-        Triangulation::Finite_edges_iterator eit = finite_edges_begin();
-        eit != finite_edges_end();
-        ++eit
-    )
-    {
-        Cell_circulator ccStart = incident_cells(*eit);
-        Cell_circulator cc = ccStart;
-
-        do
-        {
-            if (dualPtIndexMap.found(cc->cellIndex()))
-            {
-                // One of the points of this face has already been
-                // collapsed this sweep, leave for next sweep
-                continue;
-            }
-
-        } while (++cc != ccStart);
-
-        // Cell_handle c = eit->first;
-        // Vertex_handle vA = c->vertex(eit->second);
-        // Vertex_handle vB = c->vertex(eit->third);
-
-        if (isBoundaryDualFace(eit))
-        {
-            face dualFace = buildDualFace(eit);
-
-            if (dualFace.size() < 3)
-            {
-                // This face has been collapsed already
-                continue;
-            }
-
-            pointIndexHit surfHit;
-            label hitSurface;
-
-            geometryToConformTo_.findSurfaceNearest
-            (
-                dualFace.centre(pts),
-                cvMeshControls().spanSqr(),
-                surfHit,
-                hitSurface
-            );
-
-            vectorField norm(1);
-
-            allGeometry_[hitSurface].getNormal
-            (
-                List<pointIndexHit>(1, surfHit),
-                norm
-            );
-
-            const vector& surfaceNormal = norm[0];
-
-            vector faceNormal = dualFace.normal(pts);
-
-            faceNormal /= mag(faceNormal);
-
-            if ((faceNormal & surfaceNormal) < cosPerpendicularToleranceAngle)
-            {
-                collapseFaceToEdge(dualFace, pts, dualPtIndexMap);
-
-                // forAll(dualFace, fPtI)
-                // {
-                //     label ptI = dualFace[fPtI];
-
-                //     point& pt = pts[ptI];
-
-                //     pointIndexHit surfHit;
-                //     label hitSurface;
-
-                //     geometryToConformTo_.findSurfaceNearest
-                //     (
-                //         pt,
-                //         cvMeshControls().spanSqr(),
-                //         surfHit,
-                //         hitSurface
-                //     );
-
-                //     if (surfHit.hit())
-                //     {
-                //         pt = surfHit.hitPoint();
-
-                //         // dualPtIndexMap.insert(ptI, ptI);
-                //     }
-                // }
-            }
-        }
-    }
-
-    return nSmoothedVertices;
-}
-
-
-Foam::label Foam::conformalVoronoiMesh::collapseFaces
-(
-    pointField& pts,
-    Map<label>& dualPtIndexMap
-) const
-{
-    label nCollapsedFaces = 0;
-
-    scalar smallEdgeLengthCoeff = 1e-3;
-    scalar smallFaceAreaCoeff = sqr(smallEdgeLengthCoeff);
-    scalar collapseToEdgeCoeff = 0.02;
-    scalar longestEdgeLengthRatio = 0.35;
-
-    for
-    (
-        Triangulation::Finite_edges_iterator eit = finite_edges_begin();
-        eit != finite_edges_end();
-        ++eit
-    )
-    {
-        Cell_circulator ccStart = incident_cells(*eit);
-        Cell_circulator cc = ccStart;
-
-        do
-        {
-            if (dualPtIndexMap.found(cc->cellIndex()))
-            {
-                // One of the points of this face has already been
-                // collapsed this sweep, leave for next sweep
-                continue;
-            }
-
-        } while (++cc != ccStart);
-
-        Cell_handle c = eit->first;
-        Vertex_handle vA = c->vertex(eit->second);
-        Vertex_handle vB = c->vertex(eit->third);
-
-        if
-        (
-            vA->internalOrBoundaryPoint()
-         || vB->internalOrBoundaryPoint()
-        )
-        {
-            face dualFace = buildDualFace(eit);
-
-            if (dualFace.size() < 3)
-            {
-                // This face has been collapsed already
-                continue;
-            }
-
-            scalar area = dualFace.mag(pts);
-
-            scalar targetFaceSize = averageAnyCellSize(vA, vB);
-            scalar targetArea = sqr(targetFaceSize);
-
-            if (area < smallFaceAreaCoeff*targetArea)
-            {
-                // Collapse the dual face
-
-                // Determine if the face should be collapsed to a line or a
-                // point
-
-                const edgeList& eds = dualFace.edges();
-
-                label longestEdgeI = -1;
-
-                scalar longestEdgeLength = -SMALL;
-
-                scalar perimeter = 0.0;
-
-                forAll(eds, edI)
-                {
-                    scalar edgeLength = eds[edI].mag(pts);
-
-                    perimeter += edgeLength;
-
-                    if (edgeLength > longestEdgeLength)
-                    {
-                        longestEdgeI = edI;
-
-                        longestEdgeLength = edgeLength;
-                    }
-                }
-
-                if
-                (
-                    longestEdgeLength > collapseToEdgeCoeff*targetFaceSize
-                 && longestEdgeLength/perimeter > longestEdgeLengthRatio
-                )
-                {
-                    // Collapse to edge
-
-                    // Start at either end of the longest edge and consume the
-                    // rest of the points of the face
-
-                    const edge& longestEd = eds[longestEdgeI];
-
-                    label longestEdStartPtI = longestEd.start();
-                    label longestEdEndPtI = longestEd.end();
-
-                    label revEdI = longestEdgeI;
-                    label fwdEdI = longestEdgeI;
-
-                    point revPt = pts[longestEdStartPtI];
-                    point fwdPt = pts[longestEdEndPtI];
-
-                    dualPtIndexMap.insert(longestEdStartPtI, longestEdStartPtI);
-                    dualPtIndexMap.insert(longestEdEndPtI, longestEdEndPtI);
-
-                    for (label fcI = 1; fcI <= label(eds.size()/2); fcI++)
-                    {
-                        revEdI = eds.rcIndex(revEdI);
-                        fwdEdI = eds.fcIndex(fwdEdI);
-
-                        const edge& revEd = eds[revEdI];
-                        const edge& fwdEd = eds[fwdEdI];
-
-                        if (fcI < label(eds.size()/2))
-                        {
-                            revPt += pts[revEd.start()];
-                            fwdPt += pts[fwdEd.end()];
-
-                            dualPtIndexMap.insert
-                            (
-                                revEd.start(),
-                                longestEdStartPtI
-                            );
-
-                            dualPtIndexMap.insert
-                            (
-                                fwdEd.end(),
-                                longestEdEndPtI
-                            );
-                        }
-                        else
-                        {
-                            // Final circulation
-
-                            if
-                            (
-                                eds.size() % 2 == 1
-                                && revEd.start() == fwdEd.end()
-                            )
-                            {
-                                // Odd number of edges, give final point to
-                                // the edge direction that has the shorter
-                                // final edge
-
-                                if (fwdEd.mag(pts) < revEd.mag(pts))
-                                {
-                                    fwdPt += pts[fwdEd.end()];
-
-                                    dualPtIndexMap.insert
-                                    (
-                                        fwdEd.end(),
-                                        longestEdEndPtI
-                                    );
-
-                                    revPt /= fcI;
-                                    fwdPt /= (fcI + 1);
-                                }
-                                else
-                                {
-                                    revPt += pts[revEd.start()];
-
-                                    dualPtIndexMap.insert
-                                    (
-                                        revEd.start(),
-                                        longestEdStartPtI
-                                    );
-
-                                    revPt /= (fcI + 1);
-                                    fwdPt /= fcI;
-                                }
-                            }
-                            else if
-                            (
-                                eds.size() % 2 == 0
-                             && revEd.start() == fwdEd.start()
-                             && revEd.end() == fwdEd.end()
-                            )
-                            {
-                                // Even number of edges
-
-                                revPt /= fcI;
-                                fwdPt /= fcI;
-                            }
-                            else
-                            {
-                                FatalErrorIn
-                                (
-                                    "Foam::conformalVoronoiMesh::collapseFace"
-                                )
-                                    << "Face circulation failed for face "
-                                    << dualFace << nl
-                                    << exit(FatalError);
-                            }
-                        }
-                    }
-
-                    // Move the position of the accumulated points
-                    pts[longestEdStartPtI] = revPt;
-                    pts[longestEdEndPtI] = fwdPt;
-
-                    nCollapsedFaces++;
-                }
-                else if
-                (
-                    longestEdgeLength <= collapseToEdgeCoeff*targetFaceSize
-                )
-                {
-                    // Collapse to point
-
-                    point resultantPt = vector::zero;
-
-                    label collapseToPtI = dualFace[0];
-
-                    forAll(dualFace, fPtI)
-                    {
-                        label ptI = dualFace[fPtI];
-
-                        resultantPt += pts[ptI];
-
-                        dualPtIndexMap.insert(ptI, collapseToPtI);
-                    }
-
-                    resultantPt /= dualFace.size();
-
-                    pts[collapseToPtI] = resultantPt;
-
-                    nCollapsedFaces++;
-                }
-            }
-        }
-    }
-
-    return nCollapsedFaces;
-}
-
-
-void Foam::conformalVoronoiMesh::collapseFaceToEdge
-(
-    const face& f,
-    pointField& pts,
-    Map<label>& dualPtIndexMap
-) const
-{
-    const vector fC = f.centre(pts);
-
-    const edgeList& faceEds = f.edges();
-
-    const scalar fArea = f.mag(pts);
-
-    scalar averageQuality = 0.0;
-
-    vector faceCollapseDirection = vector::zero;
-
-    forAll(faceEds, edI)
-    {
-        const edge e = faceEds[edI];
-
-        triPointRef tri(pts[e.start()], pts[e.end()], fC);
-
-        scalar triArea = tri.mag();
-
-        if (triArea < SMALL)
-        {
-            // Sliver triangle, do not consider
-
-            continue;
-        }
-
-        // 0.7697997 converts from a quality being 1 for an area ratio
-        // of an equilateral triangle (4*pi/(3*sqrt(3))) to the
-        // quality being 1 for an area ratio of a quarter of a square
-        // (pi), i.e.  (4*pi/(3*sqrt(3)))/pi = 0.7697997
-
-        scalar triQuality = tri.quality()*0.7697997;
-
-        // Weight the quality by the tri area
-        averageQuality += triQuality*triArea;
-
-        // find the longest edge vector, and one of the other edges
-        // that isn't the longest
-        vector longestEdge = tri.b() - tri.a();
-        vector otherEdge = tri.c() - tri.a();
-        scalar longestEdgeMagSqr = magSqr(longestEdge);
-
-        if (magSqr(tri.c() - tri.b()) > longestEdgeMagSqr)
-        {
-            longestEdge = tri.c() - tri.b();
-            otherEdge = tri.a() - tri.b();
-            longestEdgeMagSqr = magSqr(longestEdgeMagSqr);
-        }
-
-        if (magSqr(tri.a() - tri.c()) > longestEdgeMagSqr)
-        {
-            longestEdge = tri.a() - tri.c();
-            otherEdge = tri.b() - tri.c();
-            longestEdgeMagSqr = magSqr(longestEdgeMagSqr);
-        }
-
-        vector triCollapseDirection =
-            (otherEdge & longestEdge)*longestEdge/longestEdgeMagSqr
-          - otherEdge;
-
-        triCollapseDirection /= mag(triCollapseDirection);
-
-        // Add the triCollapseDirection so that it is constructive
-
-        if
-        (
-            magSqr(faceCollapseDirection + triCollapseDirection)
-          > magSqr(faceCollapseDirection)
-        )
-        {
-            faceCollapseDirection += triCollapseDirection;
-            // faceCollapseDirection += triCollapseDirection/triQuality;
-        }
-        else
-        {
-            faceCollapseDirection -= triCollapseDirection;
-            // faceCollapseDirection -= triCollapseDirection/triQuality;
-        }
-    }
-
-    faceCollapseDirection /= mag(faceCollapseDirection);
-
-    averageQuality /= fArea;
-
-    vector collapseAxis = (faceCollapseDirection ^ (f.normal(pts)/fArea));
-
-    Info<< nl << "# averageQuality " << averageQuality << endl;
-
-    forAll(f, fPtI)
-    {
-        meshTools::writeOBJ(Info, pts[f[fPtI]]);
-    }
-
-    collapseAxis *= 3.0*mag(pts[f[0]] - fC);
-    faceCollapseDirection *= 3.0*mag(pts[f[0]] - fC);
-
-    meshTools::writeOBJ(Info, fC + faceCollapseDirection);
-    meshTools::writeOBJ(Info, fC - faceCollapseDirection);
-    meshTools::writeOBJ(Info, fC + collapseAxis);
-    meshTools::writeOBJ(Info, fC - collapseAxis);
-}
-
-
-void Foam::conformalVoronoiMesh::reindexDualVertices
-(
-    const Map<label>& dualPtIndexMap
-)
-{
-    for
-    (
-        Triangulation::Finite_cells_iterator cit = finite_cells_begin();
-        cit != finite_cells_end();
-        ++cit
-    )
-    {
-        if (dualPtIndexMap.found(cit->cellIndex()))
-        {
-            cit->cellIndex() = dualPtIndexMap[cit->cellIndex()];
-        }
-    }
-}
-
-
-void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
-(
-    pointField& points,
-    faceList& faces,
-    labelList& owner,
-    labelList& neighbour,
-    wordList& patchNames,
-    labelList& patchSizes,
-    labelList& patchStarts,
-    bool includeEmptyPatches
-) const
-{
-    patchNames = geometryToConformTo_.patchNames();
-
-    patchNames.setSize(patchNames.size() + 1);
-
-    patchNames[patchNames.size() - 1] = "cvMesh_defaultPatch";
-
-    label nPatches = patchNames.size();
-
-    List<DynamicList<face> > patchFaces(nPatches, DynamicList<face>(0));
-
-    List<DynamicList<label> > patchOwners(nPatches, DynamicList<label>(0));
-
-    faces.setSize(number_of_edges());
-
-    owner.setSize(number_of_edges());
-
-    neighbour.setSize(number_of_edges());
-
-    label dualFaceI = 0;
-
-    for
-    (
-        Triangulation::Finite_edges_iterator eit = finite_edges_begin();
-        eit != finite_edges_end();
-        ++eit
-    )
-    {
-        Cell_handle c = eit->first;
-        Vertex_handle vA = c->vertex(eit->second);
-        Vertex_handle vB = c->vertex(eit->third);
-
-        if
-        (
-            vA->internalOrBoundaryPoint()
-         || vB->internalOrBoundaryPoint()
-        )
-        {
-            face newDualFace = buildDualFace(eit);
-
-            if (newDualFace.size() >= 3)
-            {
-
-                label dcA = vA->index();
-
-                if (!vA->internalOrBoundaryPoint())
-                {
-                    dcA = -1;
-                }
-
-                label dcB = vB->index();
-
-                if (!vB->internalOrBoundaryPoint())
-                {
-                    dcB = -1;
-                }
-
-                label dcOwn = -1;
-                label dcNei = -1;
-
-                if (dcA == -1 && dcB == -1)
-                {
-                    FatalErrorIn("calcDualMesh")
-                    << "Attempting to create a face joining "
-                        << "two external dual cells "
-                        << exit(FatalError);
-                }
-                else if (dcA == -1 || dcB == -1)
-                {
-                    // boundary face, find which is the owner
-
-                    if (dcA == -1)
-                    {
-                        dcOwn = dcB;
-
-                        // reverse face order to correctly orientate normal
-                        reverse(newDualFace);
-                    }
-                    else
-                    {
-                        dcOwn = dcA;
-                    }
-
-                    // Find which patch this face is on by finding the
-                    // intersection with the surface of the Delaunay edge
-                    // generating the face and identify the region of the
-                    // intersection.
-
-                    point ptA = topoint(vA->point());
-
-                    point ptB = topoint(vB->point());
-
-                    label patchIndex = geometryToConformTo_.findPatch(ptA, ptB);
-
-                    if (patchIndex == -1)
-                    {
-                        patchIndex = patchNames.size() - 1;
-
-                        WarningIn("Foam::conformalVoronoiMesh::calcDualMesh")
-                            << "Dual face found between Dv pair " << nl
-                            << "    " << ptA << nl
-                            << "    " << ptB << nl
-                            << "    that is not on a surface patch. Adding to "
-                            << patchNames[patchIndex]
-                            << endl;
-                    }
-
-                    patchFaces[patchIndex].append(newDualFace);
-                    patchOwners[patchIndex].append(dcOwn);
-                }
-                else
-                {
-                    // internal face, find the lower cell to be the owner
-
-                    if (dcB > dcA)
-                    {
-                        dcOwn = dcA;
-                        dcNei = dcB;
-                    }
-                    else
-                    {
-                        dcOwn = dcB;
-                        dcNei = dcA;
-
-                        // reverse face order to correctly orientate normal
-                        reverse(newDualFace);
-                    }
-
-                    faces[dualFaceI] = newDualFace;
-
-                    owner[dualFaceI] = dcOwn;
-
-                    neighbour[dualFaceI] = dcNei;
-
-                    dualFaceI++;
-                }
-            }
-        }
-    }
-
-    label nInternalFaces = dualFaceI;
-
-    faces.setSize(nInternalFaces);
-    owner.setSize(nInternalFaces);
-    neighbour.setSize(nInternalFaces);
-
-    timeCheck();
-
-    sortFaces(faces, owner, neighbour);
-
-    timeCheck();
-
-    addPatches
-    (
-        nInternalFaces,
-        faces,
-        owner,
-        patchNames,
-        patchSizes,
-        patchStarts,
-        patchFaces,
-        patchOwners,
-        false
-    );
-
-    removeUnusedPoints(faces, points);
-
-    timeCheck();
-}
-
-
-void Foam::conformalVoronoiMesh::sortFaces
-(
-    faceList& faces,
-    labelList& owner,
-    labelList& neighbour
-) const
-{
-    // Upper triangular order:
-    // + owner is sorted in ascending cell order
-    // + within each block of equal value for owner, neighbour is sorted in
-    //   ascending cell order.
-    // + faces sorted to correspond
-    // e.g.
-    // owner | neighbour
-    // 0     | 2
-    // 0     | 23
-    // 0     | 71
-    // 1     | 23
-    // 1     | 24
-    // 1     | 91
-
-    // Two stage sort:
-    // 1) sort by owner
-
-    Info<< nl
-        << "    Sorting faces, owner and neighbour into upper triangular order"
-        << endl;
-
-    labelList oldToNew;
-
-    sortedOrder(owner, oldToNew);
-
-    oldToNew = invert(oldToNew.size(), oldToNew);
-
-    inplaceReorder(oldToNew, faces);
-    inplaceReorder(oldToNew, owner);
-    inplaceReorder(oldToNew, neighbour);
-
-    // 2) in each block of owners sort by neighbour
-
-    // Reset map.  Elements that are not sorted will retain their -1
-    // value, which will mean that they are ignored by inplaceReorder
-
-    oldToNew = -1;
-
-    label ownerBlockStart = 0;
-
-    for (label o = 1; o < owner.size(); o++)
-    {
-        label blockLength = -1;
-
-        if (owner[o] > owner[o-1])
-        {
-            blockLength = o - ownerBlockStart;
-        }
-        else if (o == owner.size() - 1)
-        {
-            // If the last element is not a jump in owner, then it
-            // needs to trigger a sort of the last block, but with a
-            // block length that is one element longer so that it
-            // sorts itself.
-
-            // If it is a jump in owner, then it will form a block of
-            // length one, and so will not need sorted.
-
-            blockLength = o - ownerBlockStart + 1;
-        }
-
-        if (blockLength >= 1)
-        {
-            labelList blockIndices =
-                identity(blockLength) + ownerBlockStart;
-
-            SubList<label> neighbourBlock
-            (
-                neighbour,
-                blockLength,
-                ownerBlockStart
-            );
-
-            sortedOrder(neighbourBlock, blockIndices);
-
-            blockIndices = invert(blockIndices.size(), blockIndices);
-
-            forAll(blockIndices, b)
-            {
-                oldToNew[ownerBlockStart + b] =
-                blockIndices[b] + ownerBlockStart;
-            }
-
-            ownerBlockStart = o;
-        }
-    }
-
-    // owner does not need re-sorted
-    inplaceReorder(oldToNew, faces);
-    inplaceReorder(oldToNew, neighbour);
-}
-
-
-void Foam::conformalVoronoiMesh::addPatches
-(
-    const label nInternalFaces,
-    faceList& faces,
-    labelList& owner,
-    wordList& patchNames,
-    labelList& patchSizes,
-    labelList& patchStarts,
-    List<DynamicList<face> >& patchFaces,
-    List<DynamicList<label> >& patchOwners,
-    bool includeEmptyPatches
-) const
-{
-    label nTotalPatches = patchNames.size();
-
-    label nValidPatches = 0;
-
-    PackedBoolList validPatch(nTotalPatches, false);
-
-    wordList allPatchNames = patchNames;
-
-    patchSizes.setSize(nTotalPatches);
-    patchStarts.setSize(nTotalPatches);
-
-    label nBoundaryFaces = 0;
-
-    forAll(patchFaces, p)
-    {
-        // Check if the patch has any faces.  Never create an empty
-        // default patch.
-
-        if
-        (
-            patchFaces[p].size()
-         || (includeEmptyPatches && (p != nTotalPatches - 1))
-        )
-        {
-            patchNames[nValidPatches] = allPatchNames[p];
-            patchSizes[nValidPatches] = patchFaces[p].size();
-            patchStarts[nValidPatches] = nInternalFaces + nBoundaryFaces;
-
-            nBoundaryFaces += patchSizes[p];
-
-            nValidPatches++;
-
-            validPatch[p] = 1;
-        }
-        else
-        {
-            // Warn if a patch is empty and includeEmptyPatches is
-            // false, unless it is the default patch.
-
-            if (p != nTotalPatches - 1)
-            {
-                WarningIn("void addPatches")
-                    << "Patch " << patchNames[p]
-                    << " has no faces, not creating." << endl;
-            }
-        }
-    }
-
-    patchNames.setSize(nValidPatches);
-    patchSizes.setSize(nValidPatches);
-    patchStarts.setSize(nValidPatches);
-
-    faces.setSize(nInternalFaces + nBoundaryFaces);
-    owner.setSize(nInternalFaces + nBoundaryFaces);
-
-    label faceI = nInternalFaces;
-
-    forAll(patchFaces, p)
-    {
-        if (validPatch[p])
-        {
-            forAll(patchFaces[p], f)
-            {
-                faces[faceI] = patchFaces[p][f];
-                owner[faceI] = patchOwners[p][f];
-
-                faceI++;
-            }
-        }
-    }
-}
-
-
-void Foam::conformalVoronoiMesh::removeUnusedPoints
-(
-    faceList& faces,
-    pointField& pts
-) const
-{
-    Info<< nl << "    Removing unused points after filtering" << endl;
-
-    PackedBoolList ptUsed(pts.size(), false);
-
-    // Scan all faces to find all of the points that are used
-
-    forAll(faces, fI)
-    {
-        const face& f = faces[fI];
-
-        forAll(f, fPtI)
-        {
-            ptUsed[f[fPtI]] = true;
-        }
-    }
-
-    label pointI = 0;
-    labelList oldToNew(pts.size(), -1);
-
-    // Move all of the used faces to the start of the pointField and
-    // truncate it
-
-    forAll(ptUsed, ptUI)
-    {
-        if (ptUsed[ptUI] == true)
-        {
-            oldToNew[ptUI] = pointI++;
-        }
-    }
-
-    inplaceReorder(oldToNew, pts);
-
-    Info<< "        Removing "
-        << pts.size() - pointI
-        << " unused points" << endl;
-
-    pts.setSize(pointI);
-
-    // Renumber the faces to use the new point numbers
-
-    forAll(faces, fI)
-    {
-        inplaceRenumber(oldToNew, faces[fI]);
-    }
-}
-
-
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 void Foam::conformalVoronoiMesh::move()
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
index 783089034e..e451844f36 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
@@ -343,6 +343,12 @@ private:
         //- Set the size and alignment data for each vertex
         void setVertexSizeAndAlignment();
 
+        //- Builds a dual face by circulating around the supplied edge.
+        face buildDualFace
+        (
+            const Triangulation::Finite_edges_iterator& eit
+        ) const;
+
         //- Insert the necessary point pairs to conform to the surface, either
         //  from stored results, or trigger a re-conformation
         void conformToSurface();
@@ -450,12 +456,6 @@ private:
             labelList& patchStarts
         );
 
-        //- Builds a dual face by circulating around the supplied edge
-        face buildDualFace
-        (
-            const Triangulation::Finite_edges_iterator& eit
-        ) const;
-
         //- Determines if the dual face constructed by the Delaunay
         //  edge is a boundary face
         inline bool isBoundaryDualFace
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
new file mode 100644
index 0000000000..2fbdcbe178
--- /dev/null
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
@@ -0,0 +1,1636 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "conformalVoronoiMesh.H"
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+void Foam::conformalVoronoiMesh::calcDualMesh
+(
+    pointField& points,
+    faceList& faces,
+    labelList& owner,
+    labelList& neighbour,
+    wordList& patchNames,
+    labelList& patchSizes,
+    labelList& patchStarts
+)
+{
+    timeCheck();
+
+    setVertexSizeAndAlignment();
+
+    timeCheck();
+
+    // ------>  OLD, FOR REFERENCE
+        // // ~~~~~~~~~~~ removing short edges by indexing dual vertices ~~~~~~~~~~~~~~
+
+        // for
+        // (
+        //     Triangulation::Finite_cells_iterator cit = finite_cells_begin();
+        //     cit != finite_cells_end();
+        //     ++cit
+        // )
+        // {
+        //     cit->cellIndex() = -1;
+        // }
+
+        // points.setSize(number_of_cells());
+
+        // // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        // label dualVertI = 0;
+
+        // // Scanning by number of short (dual) edges (nSE) attached to the
+        // // circumcentre of each Delaunay tet.  A Delaunay tet may only have four
+        // // dual edges emanating from its circumcentre, assigning positions and
+        // // indices to those with 4 short edges attached first, then >= 3, then >= 2
+        // // etc.
+        // for (label nSE = 4; nSE >= 0; nSE--)
+        // {
+        //     Info<< nl << "Scanning for dual vertices with >= "
+        //         << nSE
+        //         << " short edges attached." << endl;
+
+        //     for
+        //     (
+        //         Triangulation::Finite_cells_iterator cit = finite_cells_begin();
+        //         cit != finite_cells_end();
+        //         ++cit
+        //     )
+        //     {
+        //         // If the Delaunay tet has an index already then it has either
+        //         // evaluated itself and taken action or has had its index dictated
+        //         // by a neighbouring tet with more short edges attached.
+
+        //         if (cit->cellIndex() == -1)
+        //         {
+        //             point dualVertex = topoint(dual(cit));
+
+        //             label shortEdges = 0;
+
+        //             List<bool> edgeIsShort(4, false);
+
+        //             List<bool> neighbourAlreadyIndexed(4, false);
+
+        //             // Loop over the four facets of the Delaunay tet
+        //             for (label f = 0; f < 4; f++)
+        //             {
+        //                 // Check that at least one of the vertices of the facet is
+        //                 // an internal or boundary point
+        //                 if
+        //                 (
+        //                     cit->vertex(vertex_triple_index(f, 0))->
+        //                         internalOrBoundaryPoint()
+        //                  || cit->vertex(vertex_triple_index(f, 1))->
+        //                         internalOrBoundaryPoint()
+        //                  || cit->vertex(vertex_triple_index(f, 2))->
+        //                         internalOrBoundaryPoint()
+        //                 )
+        //                 {
+        //                     point neighDualVertex;
+
+        //                     label cNI = cit->neighbor(f)->cellIndex();
+
+        //                     if (cNI == -1)
+        //                     {
+        //                         neighDualVertex = topoint(dual(cit->neighbor(f)));
+        //                     }
+        //                     else
+        //                     {
+        //                         neighDualVertex = points[cNI];
+        //                     }
+
+        //                     if
+        //                     (
+        //                         magSqr(dualVertex - neighDualVertex)
+        //                       < sqr
+        //                         (
+        //                             minimumEdgeLength
+        //                             (
+        //                                 0.5*(dualVertex + neighDualVertex)
+        //                             )
+        //                         )
+        //                     )
+        //                     {
+        //                         edgeIsShort[f] = true;
+
+        //                         if (cNI > -1)
+        //                         {
+        //                             neighbourAlreadyIndexed[f] = true;
+        //                         }
+
+        //                         shortEdges++;
+        //                     }
+        //                 }
+        //             }
+
+        //             if (nSE == 0 && shortEdges == 0)
+        //             {
+        //                 // Final iteration and no short edges are found, index
+        //                 // remaining dual vertices.
+
+        //                 if
+        //                 (
+        //                     cit->vertex(0)->internalOrBoundaryPoint()
+        //                  || cit->vertex(1)->internalOrBoundaryPoint()
+        //                  || cit->vertex(2)->internalOrBoundaryPoint()
+        //                  || cit->vertex(3)->internalOrBoundaryPoint()
+        //                 )
+        //                 {
+        //                     cit->cellIndex() = dualVertI;
+        //                     points[dualVertI] = dualVertex;
+        //                     dualVertI++;
+        //                 }
+        //             }
+        //             else if
+        //             (
+        //                 shortEdges >= nSE
+        //             )
+        //             {
+        //                 // Info<< neighbourAlreadyIndexed << ' '
+        //                 //     << edgeIsShort << endl;
+
+        //                 label numUnindexedNeighbours = 1;
+
+        //                 for (label f = 0; f < 4; f++)
+        //                 {
+        //                     if (edgeIsShort[f] && !neighbourAlreadyIndexed[f])
+        //                     {
+        //                         dualVertex += topoint(dual(cit->neighbor(f)));
+
+        //                         numUnindexedNeighbours++;
+        //                     }
+        //                 }
+
+        //                 dualVertex /= numUnindexedNeighbours;
+
+        //                 label nearestExistingIndex = -1;
+
+        //                 point nearestIndexedNeighbourPos = vector::zero;
+
+        //                 scalar minDistSqrToNearestIndexedNeighbour = VGREAT;
+
+        //                 for (label f = 0; f < 4; f++)
+        //                 {
+        //                     if (edgeIsShort[f] && neighbourAlreadyIndexed[f])
+        //                     {
+        //                         label cNI = cit->neighbor(f)->cellIndex();
+
+        //                         point indexedNeighbourPos = points[cNI];
+
+        //                         if
+        //                         (
+        //                             magSqr(indexedNeighbourPos - dualVertex)
+        //                           < minDistSqrToNearestIndexedNeighbour
+        //                         )
+        //                         {
+        //                             nearestExistingIndex = cNI;
+
+        //                             nearestIndexedNeighbourPos =
+        //                             indexedNeighbourPos;
+
+        //                             minDistSqrToNearestIndexedNeighbour =
+        //                             magSqr(indexedNeighbourPos - dualVertex);
+        //                         }
+        //                     }
+        //                 }
+
+        //                 if
+        //                 (
+        //                     nearestExistingIndex > -1
+        //                  && minDistSqrToNearestIndexedNeighbour
+        //                   < sqr
+        //                     (
+        //                         minimumEdgeLength
+        //                         (
+        //                             0.5*(nearestIndexedNeighbourPos + dualVertex)
+        //                         )
+        //                     )
+        //                 )
+        //                 {
+        //                     points[nearestExistingIndex] =
+        //                     0.5*(dualVertex + nearestIndexedNeighbourPos);
+
+        //                     for (label f = 0; f < 4; f++)
+        //                     {
+        //                         if (edgeIsShort[f] && !neighbourAlreadyIndexed[f])
+        //                         {
+        //                             cit->neighbor(f)->cellIndex() =
+        //                             nearestExistingIndex;
+        //                         }
+        //                     }
+
+        //                     cit->cellIndex() = nearestExistingIndex;
+        //                 }
+        //                 else
+        //                 {
+        //                     for (label f = 0; f < 4; f++)
+        //                     {
+        //                         if (edgeIsShort[f] && !neighbourAlreadyIndexed[f])
+        //                         {
+        //                             cit->neighbor(f)->cellIndex() = dualVertI;
+        //                         }
+        //                     }
+
+        //                     cit->cellIndex() = dualVertI;
+
+        //                     points[dualVertI] = dualVertex;
+
+        //                     dualVertI++;
+        //                 }
+        //             }
+        //         }
+        //     }
+        // }
+
+        // points.setSize(dualVertI);
+
+    // <------ OLD, FOR REFERENCE
+
+    // Dual cell indexing
+
+    // Assign an index to the Delaunay vertices which will be the dual cell
+    // index used for owner neighbour assignment.
+
+    // The indices of the points are reset *which **destroys** the point-pair
+    // matching*, so the type of each vertex are reset to avoid any ambiguity.
+
+    label dualCelli = 0;
+
+    for
+    (
+        Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
+        vit != finite_vertices_end();
+        ++vit
+    )
+    {
+        if (vit->internalOrBoundaryPoint())
+        {
+            vit->type() = Vb::ptInternalPoint;
+            vit->index() = dualCelli;
+            dualCelli++;
+        }
+        else
+        {
+            vit->type() = Vb::ptFarPoint;
+            vit->index() = -1;
+        }
+    }
+
+    // Dual face filtering
+
+    // Indexing Delaunay cells, which are Dual vertices
+
+    label dualVertI = 0;
+
+    points.setSize(number_of_cells());
+
+    for
+    (
+        Triangulation::Finite_cells_iterator cit = finite_cells_begin();
+        cit != finite_cells_end();
+        ++cit
+    )
+    {
+        if
+        (
+            cit->vertex(0)->internalOrBoundaryPoint()
+         || cit->vertex(1)->internalOrBoundaryPoint()
+         || cit->vertex(2)->internalOrBoundaryPoint()
+         || cit->vertex(3)->internalOrBoundaryPoint()
+        )
+        {
+            cit->cellIndex() = dualVertI;
+            points[dualVertI] = topoint(dual(cit));
+            dualVertI++;
+        }
+        else
+        {
+            cit->cellIndex() = -1;
+        }
+    }
+
+    points.setSize(dualVertI);
+
+    // Assess close points to be merged
+
+    {
+        Info<< nl << "    Merging close points" << endl;
+
+        label nPtsMerged = 0;
+
+        do
+        {
+            Map<label> dualPtIndexMap;
+
+            nPtsMerged = mergeCloseDualVertices(points, dualPtIndexMap);
+
+            Info<< "        Merged " << nPtsMerged << " points" << endl;
+
+            reindexDualVertices(dualPtIndexMap);
+
+        } while (nPtsMerged > 0);
+    }
+
+    // Smooth the surface of the mesh
+
+    Info<< nl << "    Smoothing surface" << endl;
+
+    label nSmoothedVertices = 0;
+
+    do
+    {
+        Map<label> dualPtIndexMap;
+
+        nSmoothedVertices = smoothSurfaceDualFaces(points, dualPtIndexMap);
+
+        Info<< "        Smoothed " << nSmoothedVertices
+            << " points (0 HARD CODED)"
+            << endl;
+
+        reindexDualVertices(dualPtIndexMap);
+
+    } while (nSmoothedVertices > 0);
+
+    {
+        Info<< nl << "    Merging close points" << endl;
+
+        label nPtsMerged = 0;
+
+        do
+        {
+            Map<label> dualPtIndexMap;
+
+            nPtsMerged = mergeCloseDualVertices(points, dualPtIndexMap);
+
+            Info<< "        Merged " << nPtsMerged << " points" << endl;
+
+            reindexDualVertices(dualPtIndexMap);
+
+        } while (nPtsMerged > 0);
+    }
+
+    // Assess faces for collapse
+
+    Info<< nl << "    Collapsing unnecessary faces" << endl;
+
+    label nCollapsedFaces = 0;
+
+    do
+    {
+        Map<label> dualPtIndexMap;
+
+        nCollapsedFaces = collapseFaces(points, dualPtIndexMap);
+
+        reindexDualVertices(dualPtIndexMap);
+
+        Info<< "        Collapsed " << nCollapsedFaces << " faces" << endl;
+
+    } while (nCollapsedFaces > 0);
+
+    // Final dual face and owner neighbour construction
+
+    createFacesOwnerNeighbourAndPatches
+    (
+        points,
+        faces,
+        owner,
+        neighbour,
+        patchNames,
+        patchSizes,
+        patchStarts,
+        false
+    );
+}
+
+
+void Foam::conformalVoronoiMesh::calcTetMesh
+(
+    pointField& points,
+    faceList& faces,
+    labelList& owner,
+    labelList& neighbour,
+    wordList& patchNames,
+    labelList& patchSizes,
+    labelList& patchStarts
+)
+{
+    labelList vertexMap(number_of_vertices());
+
+    label vertI = 0;
+
+    points.setSize(number_of_vertices());
+
+    for
+    (
+        Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
+        vit != finite_vertices_end();
+        ++vit
+    )
+    {
+        if (vit->internalPoint() || vit->pairPoint())
+        {
+            vertexMap[vit->index()] = vertI;
+            points[vertI] =  topoint(vit->point());
+            vertI++;
+        }
+    }
+
+    points.setSize(vertI);
+
+    label cellI = 0;
+
+    for
+    (
+        Triangulation::Finite_cells_iterator cit = finite_cells_begin();
+        cit != finite_cells_end();
+        ++cit
+    )
+    {
+        if
+        (
+            cit->vertex(0)->internalOrBoundaryPoint()
+         || cit->vertex(1)->internalOrBoundaryPoint()
+         || cit->vertex(2)->internalOrBoundaryPoint()
+         || cit->vertex(3)->internalOrBoundaryPoint()
+        )
+        {
+             cit->cellIndex() = cellI++;
+        }
+        else
+        {
+            cit->cellIndex() = -1;
+        }
+    }
+
+    patchNames = geometryToConformTo_.patchNames();
+
+    patchNames.setSize(patchNames.size() + 1);
+
+    patchNames[patchNames.size() - 1] = "cvMesh_defaultPatch";
+
+    label nPatches = patchNames.size();
+
+    List<DynamicList<face> > patchFaces(nPatches, DynamicList<face>(0));
+
+    List<DynamicList<label> > patchOwners(nPatches, DynamicList<label>(0));
+
+    faces.setSize(number_of_facets());
+
+    owner.setSize(number_of_facets());
+
+    neighbour.setSize(number_of_facets());
+
+    label faceI = 0;
+
+    labelList verticesOnTriFace(3, -1);
+
+    face newFace(verticesOnTriFace);
+
+    for
+    (
+        Triangulation::Finite_facets_iterator fit = finite_facets_begin();
+        fit != finite_facets_end();
+        ++fit
+    )
+    {
+        const Cell_handle c1(fit->first);
+        const int oppositeVertex = fit->second;
+        const Cell_handle c2(c1->neighbor(oppositeVertex));
+
+        label c1I = c1->cellIndex();
+        label c2I = c2->cellIndex();
+
+        label ownerCell = -1;
+        label neighbourCell = -1;
+
+        if (c1I == -1 && c2I == -1)
+        {
+            // Both tets are outside, skip
+            continue;
+        }
+
+        for (label i = 0; i < 3; i++)
+        {
+            verticesOnTriFace[i] = vertexMap
+            [
+                c1->vertex(vertex_triple_index(oppositeVertex, i))->index()
+            ];
+        }
+
+        newFace = face(verticesOnTriFace);
+
+        if (c1I == -1 || c2I == -1)
+        {
+            // Boundary face...
+            if (c1I == -1)
+            {
+                //... with c1 outside
+                ownerCell = c2I;
+            }
+            else
+            {
+                // ... with c2 outside
+                ownerCell = c1I;
+
+                reverse(newFace);
+            }
+
+            label patchIndex = geometryToConformTo_.findPatch
+            (
+                newFace.centre(points)
+            );
+
+            if (patchIndex == -1)
+            {
+                patchIndex = patchNames.size() - 1;
+
+                WarningIn("Foam::conformalVoronoiMesh::calcTetMesh")
+                    << "Tet face centre at  " << nl
+                    << "    " << newFace.centre(points) << nl
+                    << "    did not find a surface patch. Adding to "
+                    << patchNames[patchIndex]
+                    << endl;
+            }
+
+            patchFaces[patchIndex].append(newFace);
+            patchOwners[patchIndex].append(ownerCell);
+        }
+        else
+        {
+            // Internal face...
+            if (c1I < c2I)
+            {
+                // ...with c1 as the ownerCell
+                ownerCell = c1I;
+                neighbourCell = c2I;
+
+                reverse(newFace);
+            }
+            else
+            {
+                // ...with c2 as the ownerCell
+                ownerCell = c2I;
+                neighbourCell = c1I;
+            }
+
+            faces[faceI] = newFace;
+            owner[faceI] = ownerCell;
+            neighbour[faceI] = neighbourCell;
+            faceI++;
+        }
+    }
+
+    label nInternalFaces = faceI;
+
+    faces.setSize(nInternalFaces);
+    owner.setSize(nInternalFaces);
+    neighbour.setSize(nInternalFaces);
+
+    sortFaces(faces, owner, neighbour);
+
+    addPatches
+    (
+        nInternalFaces,
+        faces,
+        owner,
+        patchNames,
+        patchSizes,
+        patchStarts,
+        patchFaces,
+        patchOwners,
+        false
+    );
+}
+
+
+Foam::label Foam::conformalVoronoiMesh::mergeCloseDualVertices
+(
+    const pointField& pts,
+    Map<label>& dualPtIndexMap
+) const
+{
+    label nPtsMerged = 0;
+
+    scalar closenessTolerance = 1e-4;
+
+    for
+    (
+        Triangulation::Finite_facets_iterator fit = finite_facets_begin();
+        fit != finite_facets_end();
+        ++fit
+    )
+    {
+        const Cell_handle c1(fit->first);
+        const int oppositeVertex = fit->second;
+        const Cell_handle c2(c1->neighbor(oppositeVertex));
+
+        label& c1I = c1->cellIndex();
+        label& c2I = c2->cellIndex();
+
+        if (dualPtIndexMap.found(c1I) || dualPtIndexMap.found(c2I))
+        {
+            // One of the points of this edge has already been
+            // merged this sweep, leave for next sweep
+
+            continue;
+        }
+
+        if(c1I != -1 && c2I != -1 && (c1I != c2I))
+        {
+            if
+            (
+                magSqr(pts[c1I] - pts[c2I])
+              < sqr(averageAnyCellSize(fit)*closenessTolerance)
+            )
+            {
+                dualPtIndexMap.insert(c1I, c1I);
+                dualPtIndexMap.insert(c2I, c1I);
+
+                nPtsMerged++;
+            }
+        }
+    }
+
+    return nPtsMerged;
+}
+
+Foam::label Foam::conformalVoronoiMesh::smoothSurfaceDualFaces
+(
+    pointField& pts,
+    Map<label>& dualPtIndexMap
+) const
+{
+    label nSmoothedVertices = 0;
+
+    const scalar cosPerpendicularToleranceAngle = cos(degToRad(80));
+
+    for
+    (
+        Triangulation::Finite_edges_iterator eit = finite_edges_begin();
+        eit != finite_edges_end();
+        ++eit
+    )
+    {
+        Cell_circulator ccStart = incident_cells(*eit);
+        Cell_circulator cc = ccStart;
+
+        do
+        {
+            if (dualPtIndexMap.found(cc->cellIndex()))
+            {
+                // One of the points of this face has already been
+                // collapsed this sweep, leave for next sweep
+                continue;
+            }
+
+        } while (++cc != ccStart);
+
+        if (isBoundaryDualFace(eit))
+        {
+            face dualFace = buildDualFace(eit);
+
+            if (dualFace.size() < 3)
+            {
+                // This face has been collapsed already
+                continue;
+            }
+
+            pointIndexHit surfHit;
+            label hitSurface;
+
+            geometryToConformTo_.findSurfaceNearest
+            (
+                dualFace.centre(pts),
+                cvMeshControls().spanSqr(),
+                surfHit,
+                hitSurface
+            );
+
+            vectorField norm(1);
+
+            allGeometry_[hitSurface].getNormal
+            (
+                List<pointIndexHit>(1, surfHit),
+                norm
+            );
+
+            const vector& surfaceNormal = norm[0];
+
+            // Orient the face correctly before calculating the normal
+            Cell_handle c = eit->first;
+            Vertex_handle vA = c->vertex(eit->second);
+            // Vertex_handle vB = c->vertex(eit->third);
+
+            if (!vA->internalOrBoundaryPoint())
+            {
+                reverse(dualFace);
+            }
+
+            vector faceNormal = dualFace.normal(pts);
+
+            faceNormal /= mag(faceNormal);
+
+            if ((faceNormal & surfaceNormal) < cosPerpendicularToleranceAngle)
+            {
+                collapseFaceToEdge(dualFace, pts, dualPtIndexMap);
+
+                // Force all points of the face to be on the surface
+
+                // forAll(dualFace, fPtI)
+                // {
+                //     label ptI = dualFace[fPtI];
+
+                //     point& pt = pts[ptI];
+
+                //     pointIndexHit surfHit;
+                //     label hitSurface;
+
+                //     geometryToConformTo_.findSurfaceNearest
+                //     (
+                //         pt,
+                //         cvMeshControls().spanSqr(),
+                //         surfHit,
+                //         hitSurface
+                //     );
+
+                //     if (surfHit.hit())
+                //     {
+                //         pt = surfHit.hitPoint();
+
+                //         // dualPtIndexMap.insert(ptI, ptI);
+                //     }
+                // }
+            }
+        }
+    }
+
+    return nSmoothedVertices;
+}
+
+
+Foam::label Foam::conformalVoronoiMesh::collapseFaces
+(
+    pointField& pts,
+    Map<label>& dualPtIndexMap
+) const
+{
+    label nCollapsedFaces = 0;
+
+    scalar smallEdgeLengthCoeff = 1e-3;
+    scalar smallFaceAreaCoeff = sqr(smallEdgeLengthCoeff);
+    scalar collapseToEdgeCoeff = 0.02;
+    scalar longestEdgeLengthRatio = 0.35;
+
+    for
+    (
+        Triangulation::Finite_edges_iterator eit = finite_edges_begin();
+        eit != finite_edges_end();
+        ++eit
+    )
+    {
+        Cell_circulator ccStart = incident_cells(*eit);
+        Cell_circulator cc = ccStart;
+
+        do
+        {
+            if (dualPtIndexMap.found(cc->cellIndex()))
+            {
+                // One of the points of this face has already been
+                // collapsed this sweep, leave for next sweep
+                continue;
+            }
+
+        } while (++cc != ccStart);
+
+        Cell_handle c = eit->first;
+        Vertex_handle vA = c->vertex(eit->second);
+        Vertex_handle vB = c->vertex(eit->third);
+
+        if
+        (
+            vA->internalOrBoundaryPoint()
+         || vB->internalOrBoundaryPoint()
+        )
+        {
+            face dualFace = buildDualFace(eit);
+
+            if (dualFace.size() < 3)
+            {
+                // This face has been collapsed already
+                continue;
+            }
+
+            scalar area = dualFace.mag(pts);
+
+            scalar targetFaceSize = averageAnyCellSize(vA, vB);
+            scalar targetArea = sqr(targetFaceSize);
+
+            if (area < smallFaceAreaCoeff*targetArea)
+            {
+                // Collapse the dual face
+
+                // Determine if the face should be collapsed to a line or a
+                // point
+
+                const edgeList& eds = dualFace.edges();
+
+                label longestEdgeI = -1;
+
+                scalar longestEdgeLength = -SMALL;
+
+                scalar perimeter = 0.0;
+
+                forAll(eds, edI)
+                {
+                    scalar edgeLength = eds[edI].mag(pts);
+
+                    perimeter += edgeLength;
+
+                    if (edgeLength > longestEdgeLength)
+                    {
+                        longestEdgeI = edI;
+
+                        longestEdgeLength = edgeLength;
+                    }
+                }
+
+                if
+                (
+                    longestEdgeLength > collapseToEdgeCoeff*targetFaceSize
+                 && longestEdgeLength/perimeter > longestEdgeLengthRatio
+                )
+                {
+                    // Collapse to edge
+
+                    // Start at either end of the longest edge and consume the
+                    // rest of the points of the face
+
+                    const edge& longestEd = eds[longestEdgeI];
+
+                    label longestEdStartPtI = longestEd.start();
+                    label longestEdEndPtI = longestEd.end();
+
+                    label revEdI = longestEdgeI;
+                    label fwdEdI = longestEdgeI;
+
+                    point revPt = pts[longestEdStartPtI];
+                    point fwdPt = pts[longestEdEndPtI];
+
+                    dualPtIndexMap.insert(longestEdStartPtI, longestEdStartPtI);
+                    dualPtIndexMap.insert(longestEdEndPtI, longestEdEndPtI);
+
+                    for (label fcI = 1; fcI <= label(eds.size()/2); fcI++)
+                    {
+                        revEdI = eds.rcIndex(revEdI);
+                        fwdEdI = eds.fcIndex(fwdEdI);
+
+                        const edge& revEd = eds[revEdI];
+                        const edge& fwdEd = eds[fwdEdI];
+
+                        if (fcI < label(eds.size()/2))
+                        {
+                            revPt += pts[revEd.start()];
+                            fwdPt += pts[fwdEd.end()];
+
+                            dualPtIndexMap.insert
+                            (
+                                revEd.start(),
+                                longestEdStartPtI
+                            );
+
+                            dualPtIndexMap.insert
+                            (
+                                fwdEd.end(),
+                                longestEdEndPtI
+                            );
+                        }
+                        else
+                        {
+                            // Final circulation
+
+                            if
+                            (
+                                eds.size() % 2 == 1
+                                && revEd.start() == fwdEd.end()
+                            )
+                            {
+                                // Odd number of edges, give final point to
+                                // the edge direction that has the shorter
+                                // final edge
+
+                                if (fwdEd.mag(pts) < revEd.mag(pts))
+                                {
+                                    fwdPt += pts[fwdEd.end()];
+
+                                    dualPtIndexMap.insert
+                                    (
+                                        fwdEd.end(),
+                                        longestEdEndPtI
+                                    );
+
+                                    revPt /= fcI;
+                                    fwdPt /= (fcI + 1);
+                                }
+                                else
+                                {
+                                    revPt += pts[revEd.start()];
+
+                                    dualPtIndexMap.insert
+                                    (
+                                        revEd.start(),
+                                        longestEdStartPtI
+                                    );
+
+                                    revPt /= (fcI + 1);
+                                    fwdPt /= fcI;
+                                }
+                            }
+                            else if
+                            (
+                                eds.size() % 2 == 0
+                             && revEd.start() == fwdEd.start()
+                             && revEd.end() == fwdEd.end()
+                            )
+                            {
+                                // Even number of edges
+
+                                revPt /= fcI;
+                                fwdPt /= fcI;
+                            }
+                            else
+                            {
+                                FatalErrorIn
+                                (
+                                    "Foam::conformalVoronoiMesh::collapseFace"
+                                )
+                                    << "Face circulation failed for face "
+                                    << dualFace << nl
+                                    << exit(FatalError);
+                            }
+                        }
+                    }
+
+                    // Move the position of the accumulated points
+                    pts[longestEdStartPtI] = revPt;
+                    pts[longestEdEndPtI] = fwdPt;
+
+                    nCollapsedFaces++;
+                }
+                else if
+                (
+                    longestEdgeLength <= collapseToEdgeCoeff*targetFaceSize
+                )
+                {
+                    // Collapse to point
+
+                    point resultantPt = vector::zero;
+
+                    label collapseToPtI = dualFace[0];
+
+                    forAll(dualFace, fPtI)
+                    {
+                        label ptI = dualFace[fPtI];
+
+                        resultantPt += pts[ptI];
+
+                        dualPtIndexMap.insert(ptI, collapseToPtI);
+                    }
+
+                    resultantPt /= dualFace.size();
+
+                    pts[collapseToPtI] = resultantPt;
+
+                    nCollapsedFaces++;
+                }
+            }
+        }
+    }
+
+    return nCollapsedFaces;
+}
+
+
+void Foam::conformalVoronoiMesh::collapseFaceToEdge
+(
+    const face& f,
+    pointField& pts,
+    Map<label>& dualPtIndexMap
+) const
+{
+    const vector fC = f.centre(pts);
+
+    vector fN = f.normal(pts);
+
+    const scalar fA = mag(fN);
+
+    // Face normal is now a unit vector
+    fN /= fA;
+
+    tensor J = f.inertia(pts, fC);
+
+    // Find the dominant collapse direction by finding the eigenvector
+    // that corresponds to the normal direction, discarding it.  The
+    // eigenvector corresponding to the smaller of the two remaining
+    // eigenvalues is the dominant axis in a high aspect ratio face.
+
+    // Normalise inertia tensor to remove problems with small values
+
+    J /= cmptMax(J);
+    // J /= max(eigenValues(J).x(), SMALL);
+
+    vector eVals = eigenValues(J);
+
+    // ***************************************************************
+    // The maximum eigenvalue (z()) must be the direction of the
+    // normal, as it has the greatest value.  The minimum eigenvalue
+    // is the dominant collapse axis for high aspect ratio faces.
+    // These assumption can be tested with the code below; they are
+    // only likely to be invalid for very warped faces.
+
+    // direction normalMatchEigenValueCmpt = -1;
+
+    // scalar normalMatchMaxMagDotProd = -SMALL;
+
+    // for (direction cmpt = 0; cmpt < vector::nComponents; cmpt++)
+    // {
+    //     vector eVec = eigenVector(J, eVals.component(cmpt));
+
+    //     scalar magDotProd = mag(eVec & fN);
+
+    //     if (magDotProd > normalMatchMaxMagDotProd)
+    //     {
+    //         normalMatchMaxMagDotProd = magDotProd;
+
+    //         normalMatchEigenValueCmpt = cmpt;
+    //     }
+    // }
+
+    // direction faceAxisACmpt = (normalMatchEigenValueCmpt + 1) % 3;
+    // vector faceAxisA = eigenVector(J, eVals.component(faceAxisACmpt));
+
+    // direction faceAxisBCmpt = (normalMatchEigenValueCmpt + 2) % 3;
+    // vector faceAxisB = eigenVector(J, eVals.component(faceAxisBCmpt));
+
+    // Info<< nl << "# Normal match " << normalMatchEigenValueCmpt << nl
+    //     << "# " << eVals.component(faceAxisACmpt) << nl
+    //     << "# " << eVals.component(faceAxisBCmpt) << nl
+    //     << "# Aspect ratio = " << sqrt
+    //     (
+    //         eVals.component(faceAxisBCmpt)
+    //        /eVals.component(faceAxisACmpt)
+    //     )
+    //     << endl;
+
+    // scalar scale = 2.0*mag(fC - pts[f[0]])/eVals.component(faceAxisACmpt);
+
+    // meshTools::writeOBJ(Info, fC);
+    // meshTools::writeOBJ
+    // (
+    //     Info,
+    //     fC + scale*eVals.component(faceAxisACmpt)*faceAxisA
+    // );
+    // meshTools::writeOBJ
+    // (
+    //     Info,
+    //     fC + scale*eVals.component(faceAxisBCmpt)*faceAxisB
+    // );
+
+    // Info<< "f 1 2" << nl << "f 1 3" << endl;
+
+    // forAll(f, fPtI)
+    // {
+    //     meshTools::writeOBJ(Info, pts[f[fPtI]]);
+    // }
+
+    // Info<< "f";
+
+    // forAll(f, fPtI)
+    // {
+    //     Info << " " << fPtI + 4;
+    // }
+
+    // Info << endl;
+
+    // ***************************************************************
+
+    // The inertia calculation describes the mass distribution as a
+    // function of distance squared to the axis, so the square root of
+    // the ratio of face-plane moments gives a good indication of the
+    // aspect ratio.
+    scalar aspectRatio = sqrt(eVals.y()/max(eVals.x(), SMALL));
+
+    vector collapseAxis = eigenVector(J, eVals.x());
+
+    // The signed distance along the collapse axis passing through the
+    // face centre that each vertex projects to.
+
+    List<scalar> d(f.size());
+    forAll(f, fPtI)
+    {
+        const point& pt = pts[f[fPtI]];
+
+        d[fPtI] = (collapseAxis & (pt - fC));
+    }
+
+    sort(d);
+
+    Info<< nl << "# Aspect ratio = " << aspectRatio << endl;
+
+    Info<< nl << "# eigenvalues = " << eVals << endl;
+
+    scalar scale = 2.0*mag(fC - pts[f[0]]);
+
+    meshTools::writeOBJ(Info, fC);
+    meshTools::writeOBJ(Info, fC + scale*collapseAxis);
+
+    Info<< "f 1 2" << endl;
+
+    forAll(f, fPtI)
+    {
+        meshTools::writeOBJ(Info, pts[f[fPtI]]);
+    }
+
+    Info<< "f";
+
+    forAll(f, fPtI)
+    {
+        Info << " " << fPtI + 3;
+    }
+
+    Info<< nl << "# " << d << endl;
+
+    Info<< "# " << d.first() << " " << d.last() << endl;
+
+    forAll(d, dI)
+    {
+        meshTools::writeOBJ(Info, fC + d[dI]*collapseAxis);
+    }
+}
+
+
+void Foam::conformalVoronoiMesh::reindexDualVertices
+(
+    const Map<label>& dualPtIndexMap
+)
+{
+    for
+    (
+        Triangulation::Finite_cells_iterator cit = finite_cells_begin();
+        cit != finite_cells_end();
+        ++cit
+    )
+    {
+        if (dualPtIndexMap.found(cit->cellIndex()))
+        {
+            cit->cellIndex() = dualPtIndexMap[cit->cellIndex()];
+        }
+    }
+}
+
+
+void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
+(
+    pointField& points,
+    faceList& faces,
+    labelList& owner,
+    labelList& neighbour,
+    wordList& patchNames,
+    labelList& patchSizes,
+    labelList& patchStarts,
+    bool includeEmptyPatches
+) const
+{
+    patchNames = geometryToConformTo_.patchNames();
+
+    patchNames.setSize(patchNames.size() + 1);
+
+    patchNames[patchNames.size() - 1] = "cvMesh_defaultPatch";
+
+    label nPatches = patchNames.size();
+
+    List<DynamicList<face> > patchFaces(nPatches, DynamicList<face>(0));
+
+    List<DynamicList<label> > patchOwners(nPatches, DynamicList<label>(0));
+
+    faces.setSize(number_of_edges());
+
+    owner.setSize(number_of_edges());
+
+    neighbour.setSize(number_of_edges());
+
+    label dualFaceI = 0;
+
+    for
+    (
+        Triangulation::Finite_edges_iterator eit = finite_edges_begin();
+        eit != finite_edges_end();
+        ++eit
+    )
+    {
+        Cell_handle c = eit->first;
+        Vertex_handle vA = c->vertex(eit->second);
+        Vertex_handle vB = c->vertex(eit->third);
+
+        if
+        (
+            vA->internalOrBoundaryPoint()
+         || vB->internalOrBoundaryPoint()
+        )
+        {
+            face newDualFace = buildDualFace(eit);
+
+            if (newDualFace.size() >= 3)
+            {
+                label dcA = vA->index();
+
+                if (!vA->internalOrBoundaryPoint())
+                {
+                    dcA = -1;
+                }
+
+                label dcB = vB->index();
+
+                if (!vB->internalOrBoundaryPoint())
+                {
+                    dcB = -1;
+                }
+
+                label dcOwn = -1;
+                label dcNei = -1;
+
+                if (dcA == -1 && dcB == -1)
+                {
+                    FatalErrorIn("calcDualMesh")
+                        << "Attempting to create a face joining "
+                        << "two unindexed dual cells "
+                        << exit(FatalError);
+                }
+                else if (dcA == -1 || dcB == -1)
+                {
+                    // boundary face, find which is the owner
+
+                    if (dcA == -1)
+                    {
+                        dcOwn = dcB;
+
+                        // reverse face order to correctly orientate normal
+                        reverse(newDualFace);
+                    }
+                    else
+                    {
+                        dcOwn = dcA;
+                    }
+
+                    // Find which patch this face is on by finding the
+                    // intersection with the surface of the Delaunay edge
+                    // generating the face and identify the region of the
+                    // intersection.
+
+                    point ptA = topoint(vA->point());
+
+                    point ptB = topoint(vB->point());
+
+                    label patchIndex = geometryToConformTo_.findPatch(ptA, ptB);
+
+                    if (patchIndex == -1)
+                    {
+                        patchIndex = patchNames.size() - 1;
+
+                        WarningIn("Foam::conformalVoronoiMesh::calcDualMesh")
+                            << "Dual face found between Dv pair " << nl
+                            << "    " << ptA << nl
+                            << "    " << ptB << nl
+                            << "    that is not on a surface patch. Adding to "
+                            << patchNames[patchIndex]
+                            << endl;
+                    }
+
+                    patchFaces[patchIndex].append(newDualFace);
+                    patchOwners[patchIndex].append(dcOwn);
+                }
+                else
+                {
+                    // internal face, find the lower cell to be the owner
+
+                    if (dcB > dcA)
+                    {
+                        dcOwn = dcA;
+                        dcNei = dcB;
+                    }
+                    else
+                    {
+                        dcOwn = dcB;
+                        dcNei = dcA;
+
+                        // reverse face order to correctly orientate normal
+                        reverse(newDualFace);
+                    }
+
+                    faces[dualFaceI] = newDualFace;
+
+                    owner[dualFaceI] = dcOwn;
+
+                    neighbour[dualFaceI] = dcNei;
+
+                    dualFaceI++;
+                }
+            }
+        }
+    }
+
+    label nInternalFaces = dualFaceI;
+
+    faces.setSize(nInternalFaces);
+    owner.setSize(nInternalFaces);
+    neighbour.setSize(nInternalFaces);
+
+    timeCheck();
+
+    sortFaces(faces, owner, neighbour);
+
+    timeCheck();
+
+    addPatches
+    (
+        nInternalFaces,
+        faces,
+        owner,
+        patchNames,
+        patchSizes,
+        patchStarts,
+        patchFaces,
+        patchOwners,
+        false
+    );
+
+    removeUnusedPoints(faces, points);
+
+    timeCheck();
+}
+
+
+void Foam::conformalVoronoiMesh::sortFaces
+(
+    faceList& faces,
+    labelList& owner,
+    labelList& neighbour
+) const
+{
+    // Upper triangular order:
+    // + owner is sorted in ascending cell order
+    // + within each block of equal value for owner, neighbour is sorted in
+    //   ascending cell order.
+    // + faces sorted to correspond
+    // e.g.
+    // owner | neighbour
+    // 0     | 2
+    // 0     | 23
+    // 0     | 71
+    // 1     | 23
+    // 1     | 24
+    // 1     | 91
+
+    // Two stage sort:
+    // 1) sort by owner
+
+    Info<< nl
+        << "    Sorting faces, owner and neighbour into upper triangular order"
+        << endl;
+
+    labelList oldToNew;
+
+    sortedOrder(owner, oldToNew);
+
+    oldToNew = invert(oldToNew.size(), oldToNew);
+
+    inplaceReorder(oldToNew, faces);
+    inplaceReorder(oldToNew, owner);
+    inplaceReorder(oldToNew, neighbour);
+
+    // 2) in each block of owners sort by neighbour
+
+    // Reset map.  Elements that are not sorted will retain their -1
+    // value, which will mean that they are ignored by inplaceReorder
+
+    oldToNew = -1;
+
+    label ownerBlockStart = 0;
+
+    for (label o = 1; o < owner.size(); o++)
+    {
+        label blockLength = -1;
+
+        if (owner[o] > owner[o-1])
+        {
+            blockLength = o - ownerBlockStart;
+        }
+        else if (o == owner.size() - 1)
+        {
+            // If the last element is not a jump in owner, then it
+            // needs to trigger a sort of the last block, but with a
+            // block length that is one element longer so that it
+            // sorts itself.
+
+            // If it is a jump in owner, then it will form a block of
+            // length one, and so will not need sorted.
+
+            blockLength = o - ownerBlockStart + 1;
+        }
+
+        if (blockLength >= 1)
+        {
+            labelList blockIndices =
+                identity(blockLength) + ownerBlockStart;
+
+            SubList<label> neighbourBlock
+            (
+                neighbour,
+                blockLength,
+                ownerBlockStart
+            );
+
+            sortedOrder(neighbourBlock, blockIndices);
+
+            blockIndices = invert(blockIndices.size(), blockIndices);
+
+            forAll(blockIndices, b)
+            {
+                oldToNew[ownerBlockStart + b] =
+                blockIndices[b] + ownerBlockStart;
+            }
+
+            ownerBlockStart = o;
+        }
+    }
+
+    // owner does not need re-sorted
+    inplaceReorder(oldToNew, faces);
+    inplaceReorder(oldToNew, neighbour);
+}
+
+
+void Foam::conformalVoronoiMesh::addPatches
+(
+    const label nInternalFaces,
+    faceList& faces,
+    labelList& owner,
+    wordList& patchNames,
+    labelList& patchSizes,
+    labelList& patchStarts,
+    List<DynamicList<face> >& patchFaces,
+    List<DynamicList<label> >& patchOwners,
+    bool includeEmptyPatches
+) const
+{
+    label nTotalPatches = patchNames.size();
+
+    label nValidPatches = 0;
+
+    PackedBoolList validPatch(nTotalPatches, false);
+
+    wordList allPatchNames = patchNames;
+
+    patchSizes.setSize(nTotalPatches);
+    patchStarts.setSize(nTotalPatches);
+
+    label nBoundaryFaces = 0;
+
+    forAll(patchFaces, p)
+    {
+        // Check if the patch has any faces.  Never create an empty
+        // default patch.
+
+        if
+        (
+            patchFaces[p].size()
+         || (includeEmptyPatches && (p != nTotalPatches - 1))
+        )
+        {
+            patchNames[nValidPatches] = allPatchNames[p];
+            patchSizes[nValidPatches] = patchFaces[p].size();
+            patchStarts[nValidPatches] = nInternalFaces + nBoundaryFaces;
+
+            nBoundaryFaces += patchSizes[p];
+
+            nValidPatches++;
+
+            validPatch[p] = 1;
+        }
+        else
+        {
+            // Warn if a patch is empty and includeEmptyPatches is
+            // false, unless it is the default patch.
+
+            if (p != nTotalPatches - 1)
+            {
+                WarningIn("void addPatches")
+                    << "Patch " << patchNames[p]
+                    << " has no faces, not creating." << endl;
+            }
+        }
+    }
+
+    patchNames.setSize(nValidPatches);
+    patchSizes.setSize(nValidPatches);
+    patchStarts.setSize(nValidPatches);
+
+    faces.setSize(nInternalFaces + nBoundaryFaces);
+    owner.setSize(nInternalFaces + nBoundaryFaces);
+
+    label faceI = nInternalFaces;
+
+    forAll(patchFaces, p)
+    {
+        if (validPatch[p])
+        {
+            forAll(patchFaces[p], f)
+            {
+                faces[faceI] = patchFaces[p][f];
+                owner[faceI] = patchOwners[p][f];
+
+                faceI++;
+            }
+        }
+    }
+}
+
+
+void Foam::conformalVoronoiMesh::removeUnusedPoints
+(
+    faceList& faces,
+    pointField& pts
+) const
+{
+    Info<< nl << "    Removing unused points after filtering" << endl;
+
+    PackedBoolList ptUsed(pts.size(), false);
+
+    // Scan all faces to find all of the points that are used
+
+    forAll(faces, fI)
+    {
+        const face& f = faces[fI];
+
+        forAll(f, fPtI)
+        {
+            ptUsed[f[fPtI]] = true;
+        }
+    }
+
+    label pointI = 0;
+    labelList oldToNew(pts.size(), -1);
+
+    // Move all of the used faces to the start of the pointField and
+    // truncate it
+
+    forAll(ptUsed, ptUI)
+    {
+        if (ptUsed[ptUI] == true)
+        {
+            oldToNew[ptUI] = pointI++;
+        }
+    }
+
+    inplaceReorder(oldToNew, pts);
+
+    Info<< "        Removing "
+        << pts.size() - pointI
+        << " unused points" << endl;
+
+    pts.setSize(pointI);
+
+    // Renumber the faces to use the new point numbers
+
+    forAll(faces, fI)
+    {
+        inplaceRenumber(oldToNew, faces[fI]);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
new file mode 100644
index 0000000000..0dbff5886b
--- /dev/null
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
@@ -0,0 +1,830 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "conformalVoronoiMesh.H"
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+void Foam::conformalVoronoiMesh::conformToSurface()
+{
+    reconformationMode reconfMode = reconformationControl();
+
+    if (reconfMode == rmNone)
+    {
+        // Reinsert stored surface conformation
+        reinsertSurfaceConformation();
+    }
+    else
+    {
+        // Rebuild, insert and store new surface conformation
+        buildSurfaceConformation(reconfMode);
+    }
+}
+
+Foam::conformalVoronoiMesh::reconformationMode
+Foam::conformalVoronoiMesh::reconformationControl() const
+{
+    if (!runTime_.run())
+    {
+        Info<< nl << "    Rebuilding surface conformation for final output"
+            << endl;
+
+        return rmFine;
+    }
+    else if
+    (
+        runTime_.timeIndex()
+      % cvMeshControls().surfaceConformationRebuildFrequency() == 0
+    )
+    {
+        Info<< nl << "    Rebuilding surface conformation for more iterations"
+            << endl;
+
+        return rmCoarse;
+    }
+
+    return rmNone;
+}
+
+void Foam::conformalVoronoiMesh::buildSurfaceConformation
+(
+    reconformationMode reconfMode
+)
+{
+    if (reconfMode == rmCoarse)
+    {
+        Info<< nl << "    Build coarse surface conformation" << endl;
+    }
+    else if (reconfMode == rmFine)
+    {
+        Info<< nl << "    Build fine surface conformation" << endl;
+    }
+    else if (reconfMode == rmNone)
+    {
+        WarningIn("buildSurfaceConformation(reconformationMode reconfMode)")
+            << "reconformationMode rmNone specified, not building conformation"
+            << endl;
+
+        return;
+    }
+    else
+    {
+        WarningIn("buildSurfaceConformation(reconformationMode reconfMode)")
+            << "Unknown reconformationMode " << reconfMode << endl;
+
+        return;
+    }
+
+    timeCheck();
+
+    startOfSurfacePointPairs_ = number_of_vertices();
+
+    // Initialise containers to store the edge conformation locations
+    DynamicList<point> newEdgeLocations;
+
+    pointField existingEdgeLocations(0);
+
+    autoPtr<indexedOctree<treeDataPoint> > edgeLocationTree;
+
+    // Initialise the edgeLocationTree
+    buildEdgeLocationTree(edgeLocationTree, existingEdgeLocations);
+
+    label initialTotalHits = 0;
+
+    // Initial surface protrusion conformation - nearest surface point
+    {
+        Info<< "    EDGE DISTANCE COEFFS HARD-CODED." << endl;
+        scalar edgeSearchDistCoeffSqr = sqr(1.1);
+        scalar surfacePtReplaceDistCoeffSqr = sqr(0.5);
+
+        DynamicList<pointIndexHit> surfaceHits;
+        DynamicList<label> hitSurfaces;
+
+        DynamicList<pointIndexHit> featureEdgeHits;
+        DynamicList<label> featureEdgeFeaturesHit;
+
+        for
+        (
+            Triangulation::Finite_vertices_iterator vit =
+            finite_vertices_begin();
+            vit != finite_vertices_end();
+            vit++
+        )
+        {
+            if (vit->internalPoint())
+            {
+                point vert(topoint(vit->point()));
+                scalar searchDistanceSqr = surfaceSearchDistanceSqr(vert);
+                pointIndexHit surfHit;
+                label hitSurface;
+
+                geometryToConformTo_.findSurfaceNearest
+                (
+                    vert,
+                    searchDistanceSqr,
+                    surfHit,
+                    hitSurface
+                );
+
+                if (surfHit.hit())
+                {
+                    vit->setNearBoundary();
+
+                    if (dualCellSurfaceAnyIntersection(vit))
+                    {
+                        addSurfaceAndEdgeHits
+                        (
+                            vit,
+                            vert,
+                            surfHit,
+                            hitSurface,
+                            surfacePtReplaceDistCoeffSqr,
+                            edgeSearchDistCoeffSqr,
+                            surfaceHits,
+                            hitSurfaces,
+                            featureEdgeHits,
+                            featureEdgeFeaturesHit,
+                            newEdgeLocations,
+                            existingEdgeLocations,
+                            edgeLocationTree
+                        );
+                    }
+                }
+            }
+        }
+
+        Info<< nl <<"    Initial conformation" << nl
+            << "        Number of vertices " << number_of_vertices() << nl
+            << "        Number of surface hits " << surfaceHits.size() << nl
+            << "        Number of edge hits " << featureEdgeHits.size()
+            << endl;
+
+        insertSurfacePointPairs
+        (
+            surfaceHits,
+            hitSurfaces,
+            "surfaceConformationLocations_initial.obj"
+        );
+
+        insertEdgePointGroups
+        (
+            featureEdgeHits,
+            featureEdgeFeaturesHit,
+            "edgeConformationLocations_initial.obj"
+        );
+
+        initialTotalHits = surfaceHits.size() + featureEdgeHits.size();
+    }
+
+    label iterationNo = 0;
+
+    label maxIterations = 10;
+
+    Info<< nl << "    MAX ITERATIONS HARD CODED TO "<< maxIterations << endl;
+
+    scalar iterationToIntialHitRatioLimit = 0.01;
+
+    label hitLimit = label(iterationToIntialHitRatioLimit*initialTotalHits);
+
+    Info<< "    STOPPING ITERATIONS WHEN TOTAL NUMBER OF HITS DROPS BELOW "
+        << iterationToIntialHitRatioLimit << " (HARD CODED) OF INITIAL HITS ("
+        << hitLimit << ")"
+        << endl;
+
+    // Set totalHits to a large enough positive value to enter the while loop on
+    // the first iteration
+    label totalHits = initialTotalHits;
+
+    while
+    (
+        totalHits > 0
+     && totalHits > hitLimit
+     && iterationNo < maxIterations
+    )
+    {
+        Info<< "    EDGE DISTANCE COEFFS HARD-CODED." << endl;
+        scalar edgeSearchDistCoeffSqr = sqr(1.25);
+        scalar surfacePtReplaceDistCoeffSqr = sqr(0.7);
+
+        DynamicList<pointIndexHit> surfaceHits;
+        DynamicList<label> hitSurfaces;
+
+        DynamicList<pointIndexHit> featureEdgeHits;
+        DynamicList<label> featureEdgeFeaturesHit;
+
+        for
+        (
+            Triangulation::Finite_vertices_iterator vit =
+            finite_vertices_begin();
+            vit != finite_vertices_end();
+            vit++
+        )
+        {
+            // The initial surface conformation has already identified the
+            // nearBoundary set of vertices.  Previously inserted boundary
+            // points can also generate protrusions and must be assessed too.
+
+            if (vit->nearBoundary() || vit->ppMaster())
+            {
+                point vert(topoint(vit->point()));
+                pointIndexHit surfHit;
+                label hitSurface;
+
+                dualCellLargestSurfaceProtrusion(vit, surfHit, hitSurface);
+
+                if (surfHit.hit())
+                {
+                    addSurfaceAndEdgeHits
+                    (
+                        vit,
+                        vert,
+                        surfHit,
+                        hitSurface,
+                        surfacePtReplaceDistCoeffSqr,
+                        edgeSearchDistCoeffSqr,
+                        surfaceHits,
+                        hitSurfaces,
+                        featureEdgeHits,
+                        featureEdgeFeaturesHit,
+                        newEdgeLocations,
+                        existingEdgeLocations,
+                        edgeLocationTree
+                    );
+                }
+            }
+        }
+
+        Info<< nl <<"    Conformation iteration " << iterationNo << nl
+            << "        Number of vertices " << number_of_vertices() << nl
+            << "        Number of surface hits " << surfaceHits.size() << nl
+            << "        Number of edge hits " << featureEdgeHits.size()
+            << endl;
+
+        totalHits = surfaceHits.size() + featureEdgeHits.size();
+
+        if (totalHits > 0)
+        {
+            insertSurfacePointPairs
+            (
+                surfaceHits,
+                hitSurfaces,
+                fileName
+                (
+                    "surfaceConformationLocations_" + name(iterationNo) + ".obj"
+                )
+            );
+
+            insertEdgePointGroups
+            (
+                featureEdgeHits,
+                featureEdgeFeaturesHit,
+                "edgeConformationLocations_" + name(iterationNo) + ".obj"
+            );
+        }
+
+        iterationNo++;
+
+        if (iterationNo == maxIterations)
+        {
+            WarningIn("conformalVoronoiMesh::conformToSurface()")
+                << "Maximum surface conformation iterations ("
+                << maxIterations <<  ") reached." << endl;
+        }
+
+        if (totalHits < hitLimit)
+        {
+            Info<< nl << "    Total hits (" << totalHits
+                << ") less than limit (" << hitLimit
+                << "), stopping iterations" << endl;
+        }
+    }
+
+    // Info<< nl << "    After iterations, check penetrations" << endl;
+
+    // for
+    // (
+    //     Triangulation::Finite_vertices_iterator vit =
+    //     finite_vertices_begin();
+    //     vit != finite_vertices_end();
+    //     vit++
+    // )
+    // {
+    //     if (vit->internalOrBoundaryPoint())
+    //     {
+    //         point vert(topoint(vit->point()));
+    //         pointIndexHit surfHit;
+    //         label hitSurface;
+
+    //         dualCellLargestSurfaceProtrusion(vit, surfHit, hitSurface);
+
+    //         if (surfHit.hit())
+    //         {
+    //             Info<< nl << "Residual penetration: " << nl
+    //                 << vit->index() << nl
+    //                 << vit->type() << nl
+    //                 << vit->ppMaster() << nl
+    //                 << "nearFeaturePt "
+    //                 << nearFeaturePt(surfHit.hitPoint()) << nl
+    //                 << vert << nl
+    //                 << surfHit.hitPoint()
+    //                 << endl;
+    //         }
+    //     }
+    // }
+
+    storeSurfaceConformation();
+}
+
+
+bool Foam::conformalVoronoiMesh::dualCellSurfaceAnyIntersection
+(
+    const Triangulation::Finite_vertices_iterator& vit
+) const
+{
+    std::list<Facet> facets;
+    incident_facets(vit, std::back_inserter(facets));
+
+    for
+    (
+        std::list<Facet>::iterator fit=facets.begin();
+        fit != facets.end();
+        ++fit
+    )
+    {
+        if
+        (
+            is_infinite(fit->first)
+         || is_infinite(fit->first->neighbor(fit->second))
+        )
+        {
+            return true;
+        }
+
+        point dE0 = topoint(dual(fit->first));
+
+        // If edge end is outside bounding box then edge cuts boundary
+        if (!geometryToConformTo_.bounds().contains(dE0))
+        {
+            return true;
+        }
+
+        point dE1 = topoint(dual(fit->first->neighbor(fit->second)));
+
+        // If other edge end is outside bounding box then edge cuts boundary
+        if (!geometryToConformTo_.bounds().contains(dE1))
+        {
+            return true;
+        }
+
+        // Check for the edge passing through a surface
+        if (geometryToConformTo_.findSurfaceAnyIntersection(dE0, dE1))
+        {
+            return true;
+        }
+    }
+
+    return false;
+}
+
+
+void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
+(
+    const Triangulation::Finite_vertices_iterator& vit,
+    pointIndexHit& surfHitLargest,
+    label& hitSurfaceLargest
+) const
+{
+    std::list<Facet> facets;
+    incident_facets(vit, std::back_inserter(facets));
+
+    point vert(topoint(vit->point()));
+
+    scalar maxProtrusionDistance = maxSurfaceProtrusion(vert);
+
+    for
+    (
+        std::list<Facet>::iterator fit=facets.begin();
+        fit != facets.end();
+        ++fit
+    )
+    {
+        if
+        (
+            !is_infinite(fit->first)
+         && !is_infinite(fit->first->neighbor(fit->second))
+        )
+        {
+            point edgeMid =
+                0.5
+               *(
+                    topoint(dual(fit->first))
+                  + topoint(dual(fit->first->neighbor(fit->second)))
+                );
+
+            pointIndexHit surfHit;
+            label hitSurface;
+
+            geometryToConformTo_.findSurfaceAnyIntersection
+            (
+                vert,
+                edgeMid,
+                surfHit,
+                hitSurface
+            );
+
+            if (surfHit.hit())
+            {
+                vectorField norm(1);
+
+                allGeometry_[hitSurface].getNormal
+                (
+                    List<pointIndexHit>(1, surfHit),
+                    norm
+                );
+
+                const vector& n = norm[0];
+
+                scalar normalProtrusionDistance =
+                    (edgeMid - surfHit.hitPoint()) & n;
+
+                if (normalProtrusionDistance > maxProtrusionDistance)
+                {
+                    surfHitLargest = surfHit;
+                    hitSurfaceLargest = hitSurface;
+
+                    maxProtrusionDistance = normalProtrusionDistance;
+                }
+            }
+        }
+    }
+}
+
+
+void Foam::conformalVoronoiMesh::limitDisplacement
+(
+    const Triangulation::Finite_vertices_iterator& vit,
+    vector& displacement
+) const
+{
+    point pt = topoint(vit->point());
+    point dispPt = pt + displacement;
+
+    bool limit = false;
+
+    pointIndexHit surfHit;
+    label hitSurface;
+
+    if (!geometryToConformTo_.bounds().contains(dispPt))
+    {
+        // If dispPt is outside bounding box then displacement cuts boundary
+        limit = true;
+
+        // Info<< "    bb limit" << endl;
+    }
+    else if (geometryToConformTo_.findSurfaceAnyIntersection(pt, dispPt))
+    {
+        // Full surface penetration test
+        limit = true;
+
+        // Info<< "    intersection limit" << endl;
+    }
+    else
+    {
+        // Testing if the displaced position is too close to the surface.
+        // Within twice the local surface point pair insertion distance is
+        // considered "too close"
+
+        scalar searchDistanceSqr = sqr
+        (
+            2*vit->targetCellSize()
+           *cvMeshControls().pointPairDistanceCoeff()
+        );
+
+        geometryToConformTo_.findSurfaceNearest
+        (
+            dispPt,
+            searchDistanceSqr,
+            surfHit,
+            hitSurface
+        );
+
+        if (surfHit.hit())
+        {
+            // Info<< "    proximity limit" << endl;
+
+            limit = true;
+
+            if (magSqr(pt - surfHit.hitPoint()) <= searchDistanceSqr)
+            {
+                // Info<< "    Cannot limit displacement, point " << pt
+                //     << " closer than tolerance" << endl;
+
+                return;
+            }
+        }
+    }
+
+    if (limit)
+    {
+        // Halve the displacement and call this function again.  Will continue
+        // recursively until the displacement is small enough.
+
+        displacement *= 0.5;
+
+        // Info<< "    Limiting displacement of point " << pt << endl;
+
+        limitDisplacement(vit, displacement);
+    }
+}
+
+
+bool Foam::conformalVoronoiMesh::nearFeatureEdgeLocation
+(
+    const point& pt,
+    DynamicList<point>& newEdgeLocations,
+    pointField& existingEdgeLocations,
+    autoPtr<indexedOctree<treeDataPoint> >& edgeLocationTree
+) const
+{
+    scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
+
+    // 0.01 and 1000 determined from speed tests, varying the indexedOctree
+    // rebuild frequency and balance of additions to queries.
+
+    if
+    (
+        newEdgeLocations.size()
+     >= max(0.01*existingEdgeLocations.size(), 1000)
+    )
+    {
+        existingEdgeLocations.append(newEdgeLocations);
+
+        buildEdgeLocationTree(edgeLocationTree, existingEdgeLocations);
+
+        newEdgeLocations.clear();
+    }
+    else
+    {
+        // Search for the nearest point in newEdgeLocations.
+        // Searching here first, because the intention is that the value of
+        // newEdgeLocationsSizeLimit should make this faster by design.
+
+        if (min(magSqr(newEdgeLocations - pt)) <= exclusionRangeSqr)
+        {
+            return true;
+        }
+    }
+
+    // Searching for the nearest point in existingEdgeLocations using the
+    // indexedOctree
+
+    pointIndexHit info = edgeLocationTree().findNearest(pt, exclusionRangeSqr);
+
+    return info.hit();
+}
+
+
+void Foam::conformalVoronoiMesh::buildEdgeLocationTree
+(
+    autoPtr<indexedOctree<treeDataPoint> >& edgeLocationTree,
+    const pointField& existingEdgeLocations
+) const
+{
+    treeBoundBox overallBb(geometryToConformTo_.bounds());
+
+    Random rndGen(72953);
+
+    overallBb.extend(rndGen, 1E-4);
+    overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+    overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+
+    edgeLocationTree.reset
+    (
+        new indexedOctree<treeDataPoint>
+        (
+            treeDataPoint(existingEdgeLocations),
+            overallBb,  // overall search domain
+            10,         // max levels
+            10.0,       // maximum ratio of cubes v.s. cells
+            100.0       // max. duplicity; n/a since no bounding boxes.
+        )
+    );
+}
+
+
+void Foam::conformalVoronoiMesh::buildSizeAndAlignmentTree() const
+{
+    treeBoundBox overallBb(geometryToConformTo_.bounds());
+
+    Random rndGen(627391);
+
+    overallBb.extend(rndGen, 1E-4);
+    overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+    overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+
+    sizeAndAlignmentTree_.reset
+    (
+        new indexedOctree<treeDataPoint>
+        (
+            treeDataPoint(sizeAndAlignmentLocations_),
+            overallBb,  // overall search domain
+            10,         // max levels
+            10.0,       // maximum ratio of cubes v.s. cells
+            100.0       // max. duplicity; n/a since no bounding boxes.
+        )
+    );
+}
+
+
+void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
+(
+    const Triangulation::Finite_vertices_iterator& vit,
+    const point& vert,
+    const pointIndexHit& surfHit,
+    label hitSurface,
+    scalar surfacePtReplaceDistCoeffSqr,
+    scalar edgeSearchDistCoeffSqr,
+    DynamicList<pointIndexHit>& surfaceHits,
+    DynamicList<label>& hitSurfaces,
+    DynamicList<pointIndexHit>& featureEdgeHits,
+    DynamicList<label>& featureEdgeFeaturesHit,
+    DynamicList<point>& newEdgeLocations,
+    pointField& existingEdgeLocations,
+    autoPtr<indexedOctree<treeDataPoint> >& edgeLocationTree
+) const
+{
+    bool keepSurfacePoint = true;
+
+    if (nearFeaturePt(surfHit.hitPoint()))
+    {
+        keepSurfacePoint = false;
+
+        if (vit->index() < startOfInternalPoints_)
+        {
+            surfaceHits.append(surfHit);
+
+            hitSurfaces.append(hitSurface);
+        }
+    }
+
+    List<pointIndexHit> edHits;
+
+    labelList featuresHit;
+
+    scalar targetCellSizeSqr = sqr(targetCellSize(vert));
+
+    geometryToConformTo_.findEdgeNearestByType
+    (
+        surfHit.hitPoint(),
+        edgeSearchDistCoeffSqr*targetCellSizeSqr,
+        edHits,
+        featuresHit
+    );
+
+    // Gather edge locations but do not add them to newEdgeLocations inside the
+    // loop as they will prevent nearby edge locations of different types being
+    // conformed to.
+
+    DynamicList<point> currentEdgeLocations;
+
+    forAll(edHits, i)
+    {
+        const pointIndexHit& edHit(edHits[i]);
+
+        label featureHit = featuresHit[i];
+
+        if (edHit.hit())
+        {
+            if (!nearFeaturePt(edHit.hitPoint()))
+            {
+                if
+                (
+                    magSqr(edHit.hitPoint() - surfHit.hitPoint())
+                  < surfacePtReplaceDistCoeffSqr*targetCellSizeSqr
+                )
+                {
+                    // If the point is within a given distance of a feature
+                    // edge, give control to edge control points instead, this
+                    // will prevent "pits" forming.
+
+                    keepSurfacePoint = false;
+                }
+
+                if
+                (
+                    !nearFeatureEdgeLocation
+                    (
+                        edHit.hitPoint(),
+                        newEdgeLocations,
+                        existingEdgeLocations,
+                        edgeLocationTree
+                    )
+                )
+                {
+                    // Do not place edge control points too close to a feature
+                    // point or existing edge control points
+
+                    featureEdgeHits.append(edHit);
+                    featureEdgeFeaturesHit.append(featureHit);
+
+                    currentEdgeLocations.append(edHit.hitPoint());
+                }
+            }
+        }
+    }
+
+    newEdgeLocations.append(currentEdgeLocations);
+
+    if (keepSurfacePoint)
+    {
+        surfaceHits.append(surfHit);
+
+        hitSurfaces.append(hitSurface);
+    }
+}
+
+
+void Foam::conformalVoronoiMesh::storeSurfaceConformation()
+{
+    Info<< nl << "    Storing surface conformation" << endl;
+
+    surfaceConformationVertices_.setSize
+    (
+        number_of_vertices() - startOfSurfacePointPairs_
+    );
+
+    label surfPtI = 0;
+
+    for
+    (
+        Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
+        vit != finite_vertices_end();
+        vit++
+    )
+    {
+        if (vit->index() >= startOfSurfacePointPairs_)
+        {
+            if (!vit->pairPoint())
+            {
+                FatalErrorIn("storeSurfaceConformation()")
+                    << "Trying to store a vertex that is not a surface point"
+                    << exit(FatalError);
+            }
+
+            surfaceConformationVertices_[surfPtI] = Vb(vit->point());
+
+            surfaceConformationVertices_[surfPtI].index() =
+            vit->index() - startOfSurfacePointPairs_;
+
+            surfaceConformationVertices_[surfPtI].type() =
+            vit->type() - startOfSurfacePointPairs_;
+
+            surfPtI++;
+        }
+    }
+
+    Info<< "    Stored " << surfaceConformationVertices_.size()
+        << " vertices" << endl;
+}
+
+
+void Foam::conformalVoronoiMesh::reinsertSurfaceConformation()
+{
+    Info<< nl << "    Reinserting stored surface conformation" << endl;
+
+    startOfSurfacePointPairs_ = number_of_vertices();
+
+    forAll(surfaceConformationVertices_, v)
+    {
+        insertVb(surfaceConformationVertices_[v], startOfSurfacePointPairs_);
+    }
+
+    Info<< "    Reinserted " << number_of_vertices() - startOfSurfacePointPairs_
+        << " vertices" << endl;
+}
+
+
+// ************************************************************************* //