diff --git a/applications/utilities/mesh/generation/CV2DMesher/CV2DIO.C b/applications/utilities/mesh/generation/CV2DMesher/CV2DIO.C
index 70fb8aa083..4ded03e842 100644
--- a/applications/utilities/mesh/generation/CV2DMesher/CV2DIO.C
+++ b/applications/utilities/mesh/generation/CV2DMesher/CV2DIO.C
@@ -80,7 +80,7 @@ void Foam::CV2D::writeTriangles(const fileName& fName, bool internalOnly) const
     {
         if
         (
-            !internalOnly 
+            !internalOnly
          || (
                 fit->vertex(0)->internalOrBoundaryPoint()
              || fit->vertex(1)->internalOrBoundaryPoint()
diff --git a/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index e3445e68e9..2b3dd78c30 100644
--- a/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@@ -1782,7 +1782,7 @@ void Foam::conformalVoronoiMesh::calcDualMesh
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    label dualVerti = 0;
+    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
@@ -1824,11 +1824,11 @@ void Foam::conformalVoronoiMesh::calcDualMesh
                     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()
+                            internalOrBoundaryPoint()
+                     || cit->vertex(vertex_triple_index(f, 1))->
+                            internalOrBoundaryPoint()
+                     || cit->vertex(vertex_triple_index(f, 2))->
+                            internalOrBoundaryPoint()
                     )
                     {
                         point neighDualVertex;
@@ -1874,9 +1874,9 @@ void Foam::conformalVoronoiMesh::calcDualMesh
                      || cit->vertex(3)->internalOrBoundaryPoint()
                     )
                     {
-                        cit->cellIndex() = dualVerti;
-                        points[dualVerti] = dualVertex;
-                        dualVerti++;
+                        cit->cellIndex() = dualVertI;
+                        points[dualVertI] = dualVertex;
+                        dualVertI++;
                     }
                 }
                 else if
@@ -1958,22 +1958,22 @@ void Foam::conformalVoronoiMesh::calcDualMesh
                         {
                             if (edgeIsShort[f] && !neighbourAlreadyIndexed[f])
                             {
-                                cit->neighbor(f)->cellIndex() = dualVerti;
+                                cit->neighbor(f)->cellIndex() = dualVertI;
                             }
                         }
 
-                        cit->cellIndex() = dualVerti;
+                        cit->cellIndex() = dualVertI;
 
-                        points[dualVerti] = dualVertex;
+                        points[dualVertI] = dualVertex;
 
-                        dualVerti++;
+                        dualVertI++;
                     }
                 }
             }
         }
     }
 
-    points.setSize(dualVerti);
+    points.setSize(dualVertI);
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~ dual cell indexing ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // assigns an index to the Delaunay vertices which will be the dual cell
@@ -2028,7 +2028,7 @@ void Foam::conformalVoronoiMesh::calcDualMesh
 
     neighbour.setSize(number_of_edges());
 
-    label dualFacei = 0;
+    label dualFaceI = 0;
 
     for
     (
@@ -2174,13 +2174,13 @@ void Foam::conformalVoronoiMesh::calcDualMesh
                         reverse(newDualFace);
                     }
 
-                    faces[dualFacei] = newDualFace;
+                    faces[dualFaceI] = newDualFace;
 
-                    owner[dualFacei] = dcOwn;
+                    owner[dualFaceI] = dcOwn;
 
-                    neighbour[dualFacei] = dcNei;
+                    neighbour[dualFaceI] = dcNei;
 
-                    dualFacei++;
+                    dualFaceI++;
                 }
             }
             // else
@@ -2191,7 +2191,7 @@ void Foam::conformalVoronoiMesh::calcDualMesh
         }
     }
 
-    label nInternalFaces = dualFacei;
+    label nInternalFaces = dualFaceI;
 
     faces.setSize(nInternalFaces);
 
@@ -2345,11 +2345,355 @@ void Foam::conformalVoronoiMesh::calcDualMesh
     {
         forAll(patchFaces[p], f)
         {
-            faces[dualFacei] = patchFaces[p][f];
+            faces[dualFaceI] = patchFaces[p][f];
 
-            owner[dualFacei] = patchOwners[p][f];
+            owner[dualFaceI] = patchOwners[p][f];
 
-            dualFacei++;
+            dualFaceI++;
+        }
+    }
+}
+
+
+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();
+
+    patchSizes.setSize(nPatches);
+
+    patchStarts.setSize(nPatches);
+
+    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 c1Index = c1->cellIndex();
+        label c2Index = c2->cellIndex();
+
+        label ownerCell = -1;
+        label neighbourCell = -1;
+
+        if (c1Index == -1 && c2Index == -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 (c1Index == -1 || c2Index == -1)
+        {
+            // Boundary face...
+            if (c1Index == -1)
+            {
+                //... with c1 outside
+                ownerCell = c2Index;
+            }
+            else
+            {
+                // ... with c2 outside
+                ownerCell = c1Index;
+
+                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 (c1Index < c2Index)
+            {
+                // ...with c1 as the ownerCell
+                ownerCell = c1Index;
+                neighbourCell = c2Index;
+
+                reverse(newFace);
+            }
+            else
+            {
+                // ...with c2 as the ownerCell
+                ownerCell = c2Index;
+                neighbourCell = c1Index;
+            }
+
+            faces[faceI] = newFace;
+            owner[faceI] = ownerCell;
+            neighbour[faceI] = neighbourCell;
+            faceI++;
+        }
+    }
+
+    label nInternalFaces = faceI;
+
+    faces.setSize(nInternalFaces);
+
+    owner.setSize(nInternalFaces);
+
+    neighbour.setSize(nInternalFaces);
+
+    // ~~~~~~~~ sort owner, reordinging neighbour and faces to match ~~~~~~~~~~~
+    // two stage sort for upper triangular order:  sort by owner first, then for
+    // each block of owners sort by neighbour
+
+    labelList sortingIndices;
+
+    // Stage 1
+
+    {
+        SortableList<label> sortedOwner(owner);
+
+        sortingIndices = sortedOwner.indices();
+    }
+
+    {
+        labelList copyOwner(owner.size());
+
+        forAll(sortingIndices, sI)
+        {
+            copyOwner[sI] = owner[sortingIndices[sI]];
+        }
+
+        owner = copyOwner;
+    }
+
+    {
+        labelList copyNeighbour(neighbour.size());
+
+        forAll(sortingIndices, sI)
+        {
+            copyNeighbour[sI] = neighbour[sortingIndices[sI]];
+        }
+
+        neighbour = copyNeighbour;
+    }
+
+    {
+        faceList copyFaces(faces.size());
+
+        forAll(sortingIndices, sI)
+        {
+            copyFaces[sI] = faces[sortingIndices[sI]];
+        }
+
+        faces = copyFaces;
+    }
+
+    // Stage 2
+
+    sortingIndices = -1;
+
+    DynamicList<label> ownerCellJumps;
+
+    // Force first owner entry to be a jump
+    ownerCellJumps.append(0);
+
+    for (label o = 1; o < owner.size(); o++)
+    {
+        if (owner[o] > owner[o-1])
+        {
+            ownerCellJumps.append(o);
+        }
+    }
+
+    forAll(ownerCellJumps, oCJ)
+    {
+        label start = ownerCellJumps[oCJ];
+
+        label length;
+
+        if (oCJ == ownerCellJumps.size() - 1)
+        {
+            length = owner.size() - start;
+        }
+        else
+        {
+            length = ownerCellJumps[oCJ + 1] - start;
+        }
+
+        SubList<label> neighbourBlock(neighbour, length, start);
+
+        SortableList<label> sortedNeighbourBlock(neighbourBlock);
+
+        forAll(sortedNeighbourBlock, sNB)
+        {
+            sortingIndices[start + sNB] =
+            sortedNeighbourBlock.indices()[sNB] + start;
+        }
+    }
+
+    // Perform sort
+
+    {
+        labelList copyOwner(owner.size());
+
+        forAll(sortingIndices, sI)
+        {
+            copyOwner[sI] = owner[sortingIndices[sI]];
+        }
+
+        owner = copyOwner;
+    }
+
+    {
+        labelList copyNeighbour(neighbour.size());
+
+        forAll(sortingIndices, sI)
+        {
+            copyNeighbour[sI] = neighbour[sortingIndices[sI]];
+        }
+
+        neighbour = copyNeighbour;
+    }
+
+    {
+        faceList copyFaces(faces.size());
+
+        forAll(sortingIndices, sI)
+        {
+            copyFaces[sI] = faces[sortingIndices[sI]];
+        }
+
+        faces = copyFaces;
+    }
+
+    // ~~~~~~~~ add patch information ~~~~~~~~~~~
+
+    label nBoundaryFaces = 0;
+
+    forAll(patchFaces, p)
+    {
+        patchSizes[p] = patchFaces[p].size();
+
+        patchStarts[p] = nInternalFaces + nBoundaryFaces;
+
+        nBoundaryFaces += patchSizes[p];
+    }
+
+    faces.setSize(nInternalFaces + nBoundaryFaces);
+
+    owner.setSize(nInternalFaces + nBoundaryFaces);
+
+    forAll(patchFaces, p)
+    {
+        forAll(patchFaces[p], f)
+        {
+            faces[faceI] = patchFaces[p][f];
+
+            owner[faceI] = patchOwners[p][f];
+
+            faceI++;
         }
     }
 }
@@ -2368,7 +2712,7 @@ void Foam::conformalVoronoiMesh::move()
 
     pointField dualVertices(number_of_cells());
 
-    label dualVerti = 0;
+    label dualVertI = 0;
 
     // Find the dual point of each tetrahedron and assign it an index.
     for
@@ -2388,15 +2732,15 @@ void Foam::conformalVoronoiMesh::move()
          || cit->vertex(3)->internalOrBoundaryPoint()
         )
         {
-            cit->cellIndex() = dualVerti;
+            cit->cellIndex() = dualVertI;
 
-            dualVertices[dualVerti] = topoint(dual(cit));
+            dualVertices[dualVertI] = topoint(dual(cit));
 
-            dualVerti++;
+            dualVertI++;
         }
     }
 
-    dualVertices.setSize(dualVerti);
+    dualVertices.setSize(dualVertI);
 
     timeCheck();
 
diff --git a/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
index 0830fb9103..9b84aafabe 100644
--- a/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
+++ b/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
@@ -413,6 +413,18 @@ private:
             labelList& patchStarts
         );
 
+        //- Tet mesh calculation
+        void calcTetMesh
+        (
+            pointField& points,
+            faceList& faces,
+            labelList& owner,
+            labelList& neighbour,
+            wordList& patchNames,
+            labelList& patchSizes,
+            labelList& patchStarts
+        );
+
         //- Disallow default bitwise copy construct
         conformalVoronoiMesh(const conformalVoronoiMesh&);
 
@@ -493,11 +505,26 @@ public:
             //  pointField that may be used to restart the meshing process
             void writeInternalDelaunayVertices(bool writeToConstant) const;
 
-            //- Write polyMesh
+            //- Prepare data and call writeMesh for polyMesh and
+            //  tetDualMesh
             void writeMesh(bool writeToConstant = true);
 
-            //- Write dual points and faces as .obj file
-            void writeDual
+            //- Write mesh to disk
+            void writeMesh
+            (
+                const word& meshName,
+                bool writeToConstant,
+                pointField& points,
+                faceList& faces,
+                labelList& owner,
+                labelList& neighbour,
+                wordList& patchNames,
+                labelList& patchSizes,
+                labelList& patchStarts
+            ) const;
+
+            //- Write points and faces as .obj file
+            void writeObjMesh
             (
                 const pointField& points,
                 const faceList& faces,
diff --git a/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C b/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
index 86c1817c9f..cd86db685d 100644
--- a/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
+++ b/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
@@ -132,35 +132,98 @@ void Foam::conformalVoronoiMesh::writeInternalDelaunayVertices
 
 void Foam::conformalVoronoiMesh::writeMesh(bool writeToConstant)
 {
-    pointField points(0);
-    faceList faces(0);
-    labelList owner(0);
-    labelList neighbour(0);
-    wordList patchNames(0);
-    labelList patchSizes(0);
-    labelList patchStarts(0);
-
     writeInternalDelaunayVertices(writeToConstant);
 
-    calcDualMesh
-    (
-        points,
-        faces,
-        owner,
-        neighbour,
-        patchNames,
-        patchSizes,
-        patchStarts
-    );
+    // {
+    //     pointField points;
+    //     faceList faces;
+    //     labelList owner;
+    //     labelList neighbour;
+    //     wordList patchNames;
+    //     labelList patchSizes;
+    //     labelList patchStarts;
 
+    //     calcTetMesh
+    //     (
+    //         points,
+    //         faces,
+    //         owner,
+    //         neighbour,
+    //         patchNames,
+    //         patchSizes,
+    //         patchStarts
+    //     );
+
+    //     writeMesh
+    //     (
+    //         "tetDualMesh",
+    //         writeToConstant,
+    //         points,
+    //         faces,
+    //         owner,
+    //         neighbour,
+    //         patchNames,
+    //         patchSizes,
+    //         patchStarts
+    //     );
+    // }
+
+    {
+        pointField points;
+        faceList faces;
+        labelList owner;
+        labelList neighbour;
+        wordList patchNames;
+        labelList patchSizes;
+        labelList patchStarts;
+
+        calcDualMesh
+        (
+            points,
+            faces,
+            owner,
+            neighbour,
+            patchNames,
+            patchSizes,
+            patchStarts
+        );
+
+        writeMesh
+        (
+            Foam::polyMesh::defaultRegion,
+            writeToConstant,
+            points,
+            faces,
+            owner,
+            neighbour,
+            patchNames,
+            patchSizes,
+            patchStarts
+        );
+    }
+}
+
+void Foam::conformalVoronoiMesh::writeMesh
+(
+    const word& meshName,
+    bool writeToConstant,
+    pointField& points,
+    faceList& faces,
+    labelList& owner,
+    labelList& neighbour,
+    wordList& patchNames,
+    labelList& patchSizes,
+    labelList& patchStarts
+) const
+{
     if(cvMeshControls().objOutput())
     {
-        writeDual(points, faces, "dualMesh.obj");
+        writeObjMesh(points, faces, word(meshName + ".obj"));
     }
 
     IOobject io
     (
-        Foam::polyMesh::defaultRegion,
+        meshName,
         runTime_.timeName(),
         runTime_,
         IOobject::NO_READ,
@@ -169,11 +232,11 @@ void Foam::conformalVoronoiMesh::writeMesh(bool writeToConstant)
 
     if (writeToConstant)
     {
-        Info<< nl << "    Writing polyMesh to constant." << endl;
+        Info<< nl << "    Writing mesh to constant." << endl;
 
         io = IOobject
         (
-            Foam::polyMesh::defaultRegion,
+            meshName,
             runTime_.constant(),
             runTime_,
             IOobject::NO_READ,
@@ -182,7 +245,7 @@ void Foam::conformalVoronoiMesh::writeMesh(bool writeToConstant)
     }
     else
     {
-        Info<< nl << "    Writing polyMesh to time directory "
+        Info<< nl << "    Writing mesh to time directory "
             << runTime_.timeName() << endl;
     }
 
@@ -213,21 +276,21 @@ void Foam::conformalVoronoiMesh::writeMesh(bool writeToConstant)
 
     if (!pMesh.write())
     {
-        FatalErrorIn("Foam::conformalVoronoiMesh::writeMesh()")
-        << "Failed writing polyMesh."
+        FatalErrorIn("Foam::conformalVoronoiMesh::writeMesh")
+            << "Failed writing polyMesh."
             << exit(FatalError);
     }
 }
 
 
-void Foam::conformalVoronoiMesh::writeDual
+void Foam::conformalVoronoiMesh::writeObjMesh
 (
     const pointField& points,
     const faceList& faces,
     const fileName& fName
 ) const
 {
-    Info<< nl << "    Writing dual points and faces to " << fName << endl;
+    Info<< nl << "    Writing points and faces to " << fName << endl;
 
     OFstream str(fName);
 
diff --git a/src/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C b/src/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
index ffe497d1ef..e00e9cad57 100644
--- a/src/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
+++ b/src/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
@@ -609,10 +609,7 @@ void Foam::conformationSurfaces::findEdgeNearestByType
 }
 
 
-void Foam::conformationSurfaces::writeFeatureObj
-(
-    const fileName& prefix
-) const
+void Foam::conformationSurfaces::writeFeatureObj(const fileName& prefix) const
 {
     OFstream ftStr(prefix + "_allFeatures.obj");
     Pout<< nl << "Writing all features to " << ftStr.name() << endl;
@@ -665,7 +662,37 @@ Foam::label Foam::conformationSurfaces::findPatch
 
     return
         surfLocalRegion[0] + patchOffsets_[allGeometryToSurfaces_[hitSurface]];
+}
 
+
+Foam::label Foam::conformationSurfaces::findPatch(const point& pt) const
+{
+    pointIndexHit surfHit;
+    label hitSurface;
+
+    findSurfaceNearest
+    (
+        pt,
+        cvMesh_.cvMeshControls().spanSqr(),
+        surfHit,
+        hitSurface
+    );
+
+    if (!surfHit.hit())
+    {
+        return -1;
+    }
+
+    labelList surfLocalRegion;
+
+    allGeometry_[hitSurface].getRegion
+    (
+        List<pointIndexHit>(1, surfHit),
+        surfLocalRegion
+    );
+
+    return
+        surfLocalRegion[0] + patchOffsets_[allGeometryToSurfaces_[hitSurface]];
 }
 
 
diff --git a/src/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H b/src/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H
index c0418608ab..b717c98fd4 100644
--- a/src/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H
+++ b/src/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H
@@ -255,7 +255,11 @@ public:
 
             //- Find which patch is intersected by the line from one point to
             //  another
-            label findPatch(const point& dVA, const point& dVB) const;
+            label findPatch(const point& ptA, const point& ptB) const;
+
+            //- Find which patch is closest to the point
+            label findPatch(const point& pt) const;
+
 
         // Write