diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
index f309a2700d..836f6e4cca 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
@@ -684,6 +684,8 @@ private:
             labelList& patchStarts,
             labelList& procNeighbours,
             pointField& cellCentres,
+            labelList& cellToDelaunayVertex,
+            labelListList& patchToDelaunayVertex,
             bool filterFaces
         );
 
@@ -691,6 +693,7 @@ private:
         void calcTetMesh
         (
             pointField& points,
+            labelList& pointToDelaunayVertex,
             faceList& faces,
             labelList& owner,
             labelList& neighbour,
@@ -812,14 +815,15 @@ private:
             labelList& patchSizes,
             labelList& patchStarts,
             labelList& procNeighbours,
+            labelListList& patchPointPairSlaves,
             bool includeEmptyPatches = false
         ) const;
 
         //- Create the cell centres to use for the mesh
-        void createCellCentres
-        (
-            pointField& cellCentres
-        ) const;
+        void createCellCentres(pointField& cellCentres) const;
+
+        //- Extract all points in vertex-index order
+        tmp<pointField> allPoints() const;
 
         //- Sort the faces, owner and neighbour lists into
         //  upper-triangular order.  For internal faces only, use
@@ -837,6 +841,7 @@ private:
         (
             List<DynamicList<face> >& patchFaces,
             List<DynamicList<label> >& patchOwners,
+            List<DynamicList<label> >& patchPointPairSlaves,
             List<Pair<DynamicList<label> > >& patchSortingIndices
         ) const;
 
@@ -859,12 +864,12 @@ private:
             pointField& pts
         ) const;
 
-        //- Remove dual cells that are not used by any faces
-        void removeUnusedCells
+        //- Remove dual cells that are not used by any faces. Return compaction
+        //  map.
+        labelList removeUnusedCells
         (
             labelList& owner,
-            labelList& neighbour,
-            pointField& cellCentres
+            labelList& neighbour
         ) const;
 
         //- Disallow default bitwise copy construct
@@ -1008,11 +1013,11 @@ public:
                 faceList& faces,
                 labelList& owner,
                 labelList& neighbour,
-                wordList& patchTypes,
-                wordList& patchNames,
-                labelList& patchSizes,
-                labelList& patchStarts,
-                labelList& procNeighbours,
+                const wordList& patchTypes,
+                const wordList& patchNames,
+                const labelList& patchSizes,
+                const labelList& patchStarts,
+                const labelList& procNeighbours,
                 const pointField& cellCentres
             ) const;
 
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
index d74df4f061..5a5cbfe44d 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
@@ -41,6 +41,8 @@ void Foam::conformalVoronoiMesh::calcDualMesh
     labelList& patchStarts,
     labelList& procNeighbours,
     pointField& cellCentres,
+    labelList& cellToDelaunayVertex,
+    labelListList& patchToDelaunayVertex,
     bool filterFaces
 )
 {
@@ -274,14 +276,17 @@ void Foam::conformalVoronoiMesh::calcDualMesh
         patchSizes,
         patchStarts,
         procNeighbours,
+        patchToDelaunayVertex,  // from patch face to Delaunay vertex (slavePp)
         false
     );
 
     // deferredCollapseFaceSet(owner, neighbour, deferredCollapseFaces);
 
-    createCellCentres(cellCentres);
+    cellCentres = allPoints();
 
-    removeUnusedCells(owner, neighbour, cellCentres);
+    cellToDelaunayVertex = removeUnusedCells(owner, neighbour);
+
+    cellCentres = pointField(cellCentres, cellToDelaunayVertex);
 
     removeUnusedPoints(faces, points);
 
@@ -292,6 +297,7 @@ void Foam::conformalVoronoiMesh::calcDualMesh
 void Foam::conformalVoronoiMesh::calcTetMesh
 (
     pointField& points,
+    labelList& pointToDelaunayVertex,
     faceList& faces,
     labelList& owner,
     labelList& neighbour,
@@ -306,6 +312,7 @@ void Foam::conformalVoronoiMesh::calcTetMesh
     label vertI = 0;
 
     points.setSize(number_of_vertices());
+    pointToDelaunayVertex.setSize(number_of_vertices());
 
     for
     (
@@ -318,11 +325,13 @@ void Foam::conformalVoronoiMesh::calcTetMesh
         {
             vertexMap[vit->index()] = vertI;
             points[vertI] = topoint(vit->point());
+            pointToDelaunayVertex[vertI] = vit->index();
             vertI++;
         }
     }
 
     points.setSize(vertI);
+    pointToDelaunayVertex.setSize(vertI);
 
     label cellI = 0;
 
@@ -1508,6 +1517,7 @@ Foam::labelHashSet Foam::conformalVoronoiMesh::checkPolyMeshQuality
     labelList patchStarts;
     labelList procNeighbours;
     pointField cellCentres;
+    labelListList patchToDelaunayVertex;
 
     timeCheck("Start of checkPolyMeshQuality");
 
@@ -1523,12 +1533,15 @@ Foam::labelHashSet Foam::conformalVoronoiMesh::checkPolyMeshQuality
         patchSizes,
         patchStarts,
         procNeighbours,
+        patchToDelaunayVertex,
         false
     );
 
-    createCellCentres(cellCentres);
+    //createCellCentres(cellCentres);
+    cellCentres = allPoints();
 
-    removeUnusedCells(owner, neighbour, cellCentres);
+    labelList cellToDelaunayVertex(removeUnusedCells(owner, neighbour));
+    cellCentres = pointField(cellCentres, cellToDelaunayVertex);
 
     polyMesh pMesh
     (
@@ -1887,6 +1900,7 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
     labelList& patchSizes,
     labelList& patchStarts,
     labelList& procNeighbours,
+    labelListList& patchPointPairSlaves,
     bool includeEmptyPatches
 ) const
 {
@@ -1966,6 +1980,10 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
 
     List<DynamicList<label> > patchOwners(nPatches, DynamicList<label>(0));
 
+    // Per patch face the index of the slave node of the point pair
+    List<DynamicList<label> > patchPPSlaves(nPatches, DynamicList<label>(0));
+
+
     faces.setSize(number_of_edges());
 
     owner.setSize(number_of_edges());
@@ -2101,6 +2119,16 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
 
                     patchFaces[patchIndex].append(newDualFace);
                     patchOwners[patchIndex].append(own);
+
+                    // Store the non-internal or boundary point
+                    if (vA->internalOrBoundaryPoint())
+                    {
+                        patchPPSlaves[patchIndex].append(vB->index());
+                    }
+                    else
+                    {
+                        patchPPSlaves[patchIndex].append(vA->index());
+                    }
                 }
                 else
                 {
@@ -2139,6 +2167,7 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
     (
         patchFaces,
         patchOwners,
+        patchPPSlaves,
         procPatchSortingIndex
     );
 
@@ -2154,6 +2183,13 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
         patchFaces,
         patchOwners
     );
+
+    // Return     patchPointPairSlaves.setSize(nPatches);
+    patchPointPairSlaves.setSize(nPatches);
+    forAll(patchPPSlaves, patchI)
+    {
+        patchPointPairSlaves[patchI].transfer(patchPPSlaves[patchI]);
+    }
 }
 
 
@@ -2162,7 +2198,7 @@ void Foam::conformalVoronoiMesh::createCellCentres
     pointField& cellCentres
 ) const
 {
-    cellCentres.setSize(number_of_vertices());
+    cellCentres.setSize(number_of_vertices(), point::max);
 
     label vertI = 0;
 
@@ -2185,6 +2221,28 @@ void Foam::conformalVoronoiMesh::createCellCentres
 }
 
 
+Foam::tmp<Foam::pointField> Foam::conformalVoronoiMesh::allPoints() const
+{
+    tmp<pointField> tpts(new pointField(number_of_vertices(), point::max));
+    pointField& pts = tpts();
+
+    for
+    (
+        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
+        vit != finite_vertices_end();
+        ++vit
+    )
+    {
+        if (vit->internalOrBoundaryPoint())
+        {
+            pts[vit->index()] = topoint(vit->point());
+        }
+    }
+
+    return tpts;
+}
+
+
 void Foam::conformalVoronoiMesh::sortFaces
 (
     faceList& faces,
@@ -2287,6 +2345,7 @@ void Foam::conformalVoronoiMesh::sortProcPatches
 (
     List<DynamicList<face> >& patchFaces,
     List<DynamicList<label> >& patchOwners,
+    List<DynamicList<label> >& patchPointPairSlaves,
     List<Pair<DynamicList<label> > >& patchSortingIndices
 ) const
 {
@@ -2299,6 +2358,7 @@ void Foam::conformalVoronoiMesh::sortProcPatches
     {
         faceList& faces = patchFaces[patchI];
         labelList& owner = patchOwners[patchI];
+        DynamicList<label>& slaves = patchPointPairSlaves[patchI];
 
         Pair<DynamicList<label> >& sortingIndices = patchSortingIndices[patchI];
 
@@ -2311,6 +2371,7 @@ void Foam::conformalVoronoiMesh::sortProcPatches
             (
                 faces.size() != primary.size()
              || owner.size() != primary.size()
+             || slaves.size() != primary.size()
             )
             {
                 FatalErrorIn
@@ -2326,6 +2387,7 @@ void Foam::conformalVoronoiMesh::sortProcPatches
                     << " for patch " << patchI << nl
                     << " faces.size() " << faces.size() << nl
                     << " owner.size() " << owner.size() << nl
+                    << " slaves.size() " << slaves.size() << nl
                     << " sortingIndices.first().size() "
                     << sortingIndices.first().size()
                     << exit(FatalError) << endl;
@@ -2344,6 +2406,7 @@ void Foam::conformalVoronoiMesh::sortProcPatches
             inplaceReorder(oldToNew, secondary);
             inplaceReorder(oldToNew, faces);
             inplaceReorder(oldToNew, owner);
+            inplaceReorder(oldToNew, slaves);
 
             // 2) in each block of primary sort by secondary
 
@@ -2401,6 +2464,7 @@ void Foam::conformalVoronoiMesh::sortProcPatches
 
             inplaceReorder(oldToNew, faces);
             inplaceReorder(oldToNew, owner);
+            inplaceReorder(oldToNew, slaves);
         }
     }
 }
@@ -2505,11 +2569,10 @@ void Foam::conformalVoronoiMesh::removeUnusedPoints
 }
 
 
-void Foam::conformalVoronoiMesh::removeUnusedCells
+Foam::labelList Foam::conformalVoronoiMesh::removeUnusedCells
 (
     labelList& owner,
-    labelList& neighbour,
-    pointField& cellCentres
+    labelList& neighbour
 ) const
 {
     Info<< nl << "Removing unused cells" << endl;
@@ -2543,9 +2606,7 @@ void Foam::conformalVoronoiMesh::removeUnusedCells
         }
     }
 
-    inplaceReorder(oldToNew, cellCentres);
-
-    cellCentres.setSize(cellI);
+    labelList newToOld(invert(cellI, oldToNew));
 
     // Find all of the unused cells, create a list of them, then
     // subtract one from each owner and neighbour entry for each of
@@ -2580,6 +2641,8 @@ void Foam::conformalVoronoiMesh::removeUnusedCells
             n -= findLower(unusedCells, n) + 1;
         }
     }
+
+    return newToOld;
 }
 
 
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
index 4600237231..77569aa016 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
@@ -250,49 +250,12 @@ void Foam::conformalVoronoiMesh::writeMesh
 {
     writeInternalDelaunayVertices(instance);
 
-    if (cvMeshControls().writeTetDualMesh())
-    {
-        pointField points;
-        faceList faces;
-        labelList owner;
-        labelList neighbour;
-        wordList patchTypes;
-        wordList patchNames;
-        labelList patchSizes;
-        labelList patchStarts;
-        pointField cellCentres;
-
-        calcTetMesh
-        (
-            points,
-            faces,
-            owner,
-            neighbour,
-            patchTypes,
-            patchNames,
-            patchSizes,
-            patchStarts
-        );
-        labelList procNeighbours(patchNames.size(), -1);
-
-        Info<< nl << "Writing tetDualMesh to " << instance << endl;
-
-        writeMesh
-        (
-            "tetDualMesh",
-            instance,
-            points,
-            faces,
-            owner,
-            neighbour,
-            patchTypes,
-            patchNames,
-            patchSizes,
-            patchStarts,
-            procNeighbours,
-            cellCentres
-        );
-    }
+    // Per cell the Delaunay vertex
+    labelList cellToDelaunayVertex;
+    // Per patch, per face the Delaunay vertex
+    labelListList patchToDelaunayVertex;
+    // Per patch the start of the dual faces
+    labelList dualPatchStarts;
 
     {
         pointField points;
@@ -302,7 +265,6 @@ void Foam::conformalVoronoiMesh::writeMesh
         wordList patchTypes;
         wordList patchNames;
         labelList patchSizes;
-        labelList patchStarts;
         labelList procNeighbours;
         pointField cellCentres;
 
@@ -315,9 +277,11 @@ void Foam::conformalVoronoiMesh::writeMesh
             patchTypes,
             patchNames,
             patchSizes,
-            patchStarts,
+            dualPatchStarts,
             procNeighbours,
             cellCentres,
+            cellToDelaunayVertex,
+            patchToDelaunayVertex,
             filterFaces
         );
 
@@ -334,6 +298,207 @@ void Foam::conformalVoronoiMesh::writeMesh
             patchTypes,
             patchNames,
             patchSizes,
+            dualPatchStarts,
+            procNeighbours,
+            cellCentres
+        );
+    }
+
+
+
+    if (cvMeshControls().writeTetDualMesh())
+    {
+        // Determine map from Delaunay vertex to Dual mesh
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        
+        // From all Delaunay vertices to cell (positive index)
+        // or patch face (negative index)
+        labelList vertexToDualAddressing(number_of_vertices(), 0);
+        
+        forAll(cellToDelaunayVertex, cellI)
+        {
+            label vertI = cellToDelaunayVertex[cellI];
+        
+            if (vertexToDualAddressing[vertI] != 0)
+            {
+                FatalErrorIn("conformalVoronoiMesh::writeMesh(..)")
+                    << "Delaunay vertex " << vertI
+                    << " from cell " << cellI
+                    << " is already mapped to "
+                    << vertexToDualAddressing[vertI]
+                    << exit(FatalError);
+            }
+            vertexToDualAddressing[vertI] = cellI+1;
+        }
+        
+        forAll(patchToDelaunayVertex, patchI)
+        {
+            const labelList& patchVertices = patchToDelaunayVertex[patchI];
+            forAll(patchVertices, i)
+            {
+                label vertI = patchVertices[i];
+        
+                if (vertexToDualAddressing[vertI] > 0)
+                {
+                    FatalErrorIn("conformalVoronoiMesh::writeMesh(..)")
+                        << "Delaunay vertex " << vertI
+                        << " from patch " << patchI
+                        << " local index " << i
+                        << " is already mapped to cell "
+                        << vertexToDualAddressing[vertI]-1
+                        << exit(FatalError);
+                }
+        
+                // Vertex might be used by multiple faces. Which one to
+                // use? For now last one wins.
+                label dualFaceI = dualPatchStarts[patchI]+i;
+                vertexToDualAddressing[vertI] = -dualFaceI-1;
+            }
+        }
+
+
+        // Calculate tet mesh addressing
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        pointField points;
+        // From tet point back to Delaunay vertex index
+        labelList pointToDelaunayVertex;
+        faceList faces;
+        labelList owner;
+        labelList neighbour;
+        wordList patchTypes;
+        wordList patchNames;
+        labelList patchSizes;
+        labelList patchStarts;
+        pointField cellCentres;
+
+        calcTetMesh
+        (
+            points,
+            pointToDelaunayVertex,
+            faces,
+            owner,
+            neighbour,
+            patchTypes,
+            patchNames,
+            patchSizes,
+            patchStarts
+        );
+
+
+
+        // Calculate map from tet points to dual mesh cells/patch faces
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        
+        labelIOList pointDualAddressing
+        (
+            IOobject
+            (
+                "pointDualAddressing",
+                instance,
+                "tetDualMesh"/polyMesh::meshSubDir,
+                runTime_,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            UIndirectList<label>
+            (
+                vertexToDualAddressing,
+                pointToDelaunayVertex
+            )()
+        );
+        
+        label pointI = findIndex(pointDualAddressing, -1);
+        if (pointI != -1)
+        {
+            WarningIn
+            (
+                "conformalVoronoiMesh::writeMesh\n"
+                "(\n"
+                "    const fileName& instance,\n"
+                "    bool filterFaces\n"
+                ")\n"
+            )   << "Delaunay vertex " << pointI
+                << " does not have a corresponding dual cell." << endl;
+        }
+        
+        Info<< "Writing map from tetDualMesh points to Voronoi mesh to "
+            << pointDualAddressing.objectPath() << endl;
+        pointDualAddressing.write();
+        
+        
+        
+        // Write tet points corresponding to the Voronoi cell/face centre
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        {
+            // Read Voronoi mesh
+            fvMesh mesh
+            (
+                IOobject
+                (
+                    Foam::polyMesh::defaultRegion,
+                    instance,
+                    runTime_,
+                    IOobject::MUST_READ
+                )
+            );
+            pointIOField dualPoints
+            (
+                IOobject
+                (
+                    "dualPoints",
+                    instance,
+                    "tetDualMesh"/polyMesh::meshSubDir,
+                    runTime_,
+                    IOobject::NO_READ,
+                    IOobject::AUTO_WRITE,
+                    false
+                ),
+                points
+            );
+        
+            forAll(pointDualAddressing, pointI)
+            {
+                label index = pointDualAddressing[pointI];
+        
+                if (index > 0)
+                {
+                    label cellI = index-1;
+                    dualPoints[pointI] = mesh.cellCentres()[cellI];
+                }
+                else if (index < 0)
+                {
+                    label faceI = -index-1;
+                    if (faceI >= mesh.nInternalFaces())
+                    {
+                        dualPoints[pointI] = mesh.faceCentres()[faceI];
+                    }
+                }
+            }
+        
+            Info<< "Writing new tetDualMesh points mapped onto Voronoi mesh to "
+                << dualPoints.objectPath() << endl
+                << "Replace the polyMesh/points with these." << endl;
+            dualPoints.write();
+        }
+
+
+        labelList procNeighbours(patchNames.size(), -1);
+
+        Info<< nl << "Writing tetDualMesh to " << instance << endl;
+
+        writeMesh
+        (
+            "tetDualMesh",
+            instance,
+            points,
+            faces,
+            owner,
+            neighbour,
+            patchTypes,
+            patchNames,
+            patchSizes,
             patchStarts,
             procNeighbours,
             cellCentres
@@ -350,11 +515,11 @@ void Foam::conformalVoronoiMesh::writeMesh
     faceList& faces,
     labelList& owner,
     labelList& neighbour,
-    wordList& patchTypes,
-    wordList& patchNames,
-    labelList& patchSizes,
-    labelList& patchStarts,
-    labelList& procNeighbours,
+    const wordList& patchTypes,
+    const wordList& patchNames,
+    const labelList& patchSizes,
+    const labelList& patchStarts,
+    const labelList& procNeighbours,
     const pointField& cellCentres
 ) const
 {
@@ -425,12 +590,9 @@ void Foam::conformalVoronoiMesh::writeMesh
 
     mesh.addFvPatches(patches);
 
-    // Info<< "ADDPATCHES NOT IN PARALLEL" << endl;
-    // mesh.addFvPatches(patches, false);
-
     if (!mesh.write())
     {
-        FatalErrorIn("Foam::conformalVoronoiMesh::writeMesh")
+        FatalErrorIn("Foam::conformalVoronoiMesh::writeMesh(..)")
             << "Failed writing polyMesh."
             << exit(FatalError);
     }