diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index fd1e5c4d42..e146e488a0 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@@ -1030,6 +1030,46 @@ Foam::face Foam::conformalVoronoiMesh::buildDualFace
 }
 
 
+Foam::scalar Foam::conformalVoronoiMesh::minFilterLimit
+(
+    const Triangulation::Finite_edges_iterator& eit
+) const
+{
+    Cell_circulator ccStart = incident_cells(*eit);
+    Cell_circulator cc = ccStart;
+
+    scalar minFilterLimit = GREAT;
+
+    do
+    {
+        if (cc->cellIndex() < 0)
+        {
+            Cell_handle c = eit->first;
+            Vertex_handle vA = c->vertex(eit->second);
+            Vertex_handle vB = c->vertex(eit->third);
+
+            FatalErrorIn("Foam::conformalVoronoiMesh::buildDualFace")
+                << "Dual face uses circumcenter defined by a "
+                << "Delaunay tetrahedron with no internal "
+                << "or boundary points.  Defining Delaunay edge ends: "
+                << topoint(vA->point()) << " "
+                << topoint(vB->point()) << nl
+                << exit(FatalError);
+        }
+
+        if (cc->filterLimit() < minFilterLimit)
+        {
+            minFilterLimit = cc->filterLimit();
+        }
+
+        cc++;
+
+    } while (cc != ccStart);
+
+    return minFilterLimit;
+}
+
+
 bool Foam::conformalVoronoiMesh::ownerAndNeighbour
 (
     Vertex_handle vA,
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
index 9e9c45fee7..095ddc5a18 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
@@ -357,6 +357,14 @@ private:
             const Triangulation::Finite_edges_iterator& eit
         ) const;
 
+        //- Finds the minimum filterLimit of the dual vertices
+        //  (Delaunay cells) that form the dual face produced by the
+        //  supplied edge
+        scalar minFilterLimit
+        (
+            const Triangulation::Finite_edges_iterator& eit
+        ) const;
+
         //- Determines the owner and neighbour labels for dual cells
         //  corresponding to the dual face formed by the supplied
         //  Delaunay vertices.  If the dual face is a boundary face
@@ -552,6 +560,10 @@ private:
         //  elements damage the mesh quality, allowing backtracking.
         label checkPolyMeshQuality(const pointField& pts) const;
 
+        //- Index all of the the Delaunay cells and calculate their
+        //- dual points
+        void indexDualVertices(pointField& pts);
+
         //- Re-index all of the the Delaunay cells
         void reindexDualVertices(const Map<label>& dualPtIndexMap);
 
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
index b0f9adbf78..489401a77d 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
@@ -76,53 +76,28 @@ void Foam::conformalVoronoiMesh::calcDualMesh
         }
     }
 
-    // Indexing Delaunay cells, which are the dual vertices
+    indexDualVertices(points);
 
-    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;
-        }
+        // No-risk face filtering to get rid of zero area faces and
+        // establish if the mesh can be produced at all to the
+        // specified criteria
+
+        Info<< nl << "    Merging close points" << endl;
+
+        // There is no guarantee that a merge of close points is
+        // no-risk, but it seems to work using 1e-4 as the mergeClosenessCoeff
+        mergeCloseDualVertices(points);
     }
 
-    points.setSize(dualVertI);
-
-    // No-risk face filtering to get rid of zero area faces and
-    // establish if the mesh can be produced at all to the
-    // specified criteria
-
-    Info<< nl << "    Merging close points" << endl;
-
-    mergeCloseDualVertices(points);
-
     label nInitialBadQualityFaces = checkPolyMeshQuality(points);
 
     Info<< "Initial check before face collapse, found "
         << nInitialBadQualityFaces << " bad quality faces"
         << endl;
 
+    HashSet<labelPair, labelPair::Hash<> > deferredCollapseFaces;
+
     if (nInitialBadQualityFaces > 0)
     {
         Info<< "A mesh could not be produced to satisfy the specified quality "
@@ -132,24 +107,41 @@ void Foam::conformalVoronoiMesh::calcDualMesh
             << "be applied in areas where problems are occurring."
             << endl;
     }
-
-    HashSet<labelPair, labelPair::Hash<> > deferredCollapseFaces;
-
+    else
     {
-        // Risky and undo-able face filtering to reduce the face count
-        // as much as possible staying within the specified criteria
+        label nBadQualityFaces = 0;
 
-        Info<< nl << "    Smoothing surface" << endl;
+        do
+        {
+            // Reindexing the Delaunay cells and regenerating the
+            // points resets the mesh to the starting condition.
 
-        smoothSurface(points);
+            indexDualVertices(points);
 
-        Info<< nl << "    Collapsing unnecessary faces" << endl;
+            {
+                Info<< nl << "    Merging close points" << endl;
 
-        collapseFaces(points, deferredCollapseFaces);
+                mergeCloseDualVertices(points);
+            }
 
-        label nBadQualityFaces = checkPolyMeshQuality(points);
+            {
+                // Risky and undo-able face filtering to reduce the face count
+                // as much as possible, staying within the specified criteria
 
-        Info<< "Found " << nBadQualityFaces << " bad quality faces" << endl;
+                Info<< nl << "    Smoothing surface" << endl;
+
+                smoothSurface(points);
+
+                Info<< nl << "    Collapsing unnecessary faces" << endl;
+
+                collapseFaces(points, deferredCollapseFaces);
+            }
+
+            nBadQualityFaces = checkPolyMeshQuality(points);
+
+            Info<< "Found " << nBadQualityFaces << " bad quality faces" << endl;
+
+        } while (nBadQualityFaces > 0);
     }
 
     // Final dual face and owner neighbour construction
@@ -485,6 +477,12 @@ void Foam::conformalVoronoiMesh::smoothSurface(pointField& pts)
     {
         label ptI = cit->cellIndex();
 
+        if (cit->filterLimit() < 1.0)
+        {
+            // This vertex has been limited, skip
+            continue;
+        }
+
         // Only cells with indices > -1 are valid
         if (ptI > -1)
         {
@@ -580,6 +578,14 @@ Foam::label Foam::conformalVoronoiMesh::smoothSurfaceDualFaces
                 continue;
             }
 
+            scalar minFL = minFilterLimit(eit);
+
+            if (minFL < 1.0)
+            {
+                // A vertex on this face has been limited, skip
+                continue;
+            }
+
             pointIndexHit surfHit;
             label hitSurface;
 
@@ -749,6 +755,14 @@ Foam::label Foam::conformalVoronoiMesh::collapseFaces
                 continue;
             }
 
+            scalar minFL = minFilterLimit(eit);
+
+            if (minFL < 1.0)
+            {
+                // A vertex on this face has been limited, skip
+                continue;
+            }
+
             scalar targetFaceSize = averageAnyCellSize(vA, vB);
 
             faceCollapseMode mode = collapseFace
@@ -1256,31 +1270,138 @@ Foam::label Foam::conformalVoronoiMesh::checkPolyMeshQuality
 
     motionSmoother::checkMesh
     (
-        false,
+        false,  // report
         pMesh,
         cvMeshControls().cvMeshDict().subDict("meshQualityControls"),
         wrongFaces
     );
 
+    label nInvalidPolyhedra = 0;
+
     const cellList& cells = pMesh.cells();
 
     forAll(cells, cI)
     {
-        if (cells[cI].size() < 4)
+        if (cells[cI].size() < 4 && cells[cI].size() > 0)
         {
-            Info<< "cell " << cI
-                << " has " << cells[cI].size() << " faces."
-                << endl;
+            // Info<< "cell " << cI << " " << cells[cI]
+            //     << " has " << cells[cI].size() << " faces."
+            //     << endl;
+
+            nInvalidPolyhedra++;
+
+            forAll(cells[cI], cFI)
+            {
+                wrongFaces.insert(cells[cI][cFI]);
+            }
         }
     }
 
+    Info<< "Cells with fewer than 4 faces                              : "
+        << nInvalidPolyhedra << endl;
+
+    PackedBoolList ptToBeLimited(pts.size(), false);
+
+    forAllConstIter(labelHashSet, wrongFaces, iter)
+    {
+        const face f = pMesh.faces()[iter.key()];
+
+        forAll(f, fPtI)
+        {
+            ptToBeLimited[f[fPtI]] = true;
+        }
+    }
+
+    // // Limit connected cells
+
+    // labelHashSet limitCells(pMesh.nCells()/100);
+
+    // const labelListList& ptCells = pMesh.pointCells();
+
     // forAllConstIter(labelHashSet, wrongFaces, iter)
     // {
-    //     label faceI = iter.key();
+    //     const face f = pMesh.faces()[iter.key()];
 
-    //     Info<< faceI << " " << pMesh.faces()[faceI] << endl;
+    //     forAll(f, fPtI)
+    //     {
+    //         label ptI = f[fPtI];
+
+    //         const labelList& pC = ptCells[ptI];
+
+    //         forAll(pC, pCI)
+    //         {
+    //             limitCells.insert(pC[pCI]);
+    //         }
+    //     }
     // }
 
+    // const labelListList& cellPts = pMesh.cellPoints();
+
+    // forAllConstIter(labelHashSet, limitCells, iter)
+    // {
+    //     label cellI = iter.key();
+
+    //     const labelList& cP = cellPts[cellI];
+
+    //     forAll(cP, cPI)
+    //     {
+    //         ptToBeLimited[cP[cPI]] = true;
+    //     }
+    // }
+
+
+    // Limit Delaunay cell filter values
+
+    for
+    (
+        Triangulation::Finite_cells_iterator cit = finite_cells_begin();
+        cit != finite_cells_end();
+        ++cit
+    )
+    {
+        label cI = cit->cellIndex();
+
+        if (cI >= 0)
+        {
+            if (ptToBeLimited[cI] == true)
+            {
+                cit->filterLimit() *= 0.75;
+            }
+        }
+    }
+
+    // Determine the minimum minFilterLimit
+
+    scalar minMinFilterLimit = GREAT;
+
+    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()
+        )
+        {
+            scalar minFL = minFilterLimit(eit);
+
+            if (minFL < minMinFilterLimit)
+            {
+                minMinFilterLimit = minFL;
+            }
+        }
+    }
+
+    Info<< "minMinFilterLimit " << minMinFilterLimit << endl;
+
     return wrongFaces.size();
 
     // For parallel running:
@@ -1292,6 +1413,46 @@ Foam::label Foam::conformalVoronoiMesh::checkPolyMeshQuality
 }
 
 
+void Foam::conformalVoronoiMesh::indexDualVertices
+(
+    pointField& pts
+)
+{
+    // Indexing Delaunay cells, which are the dual vertices
+
+    label dualVertI = 0;
+
+    pts.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;
+            pts[dualVertI] = topoint(dual(cit));
+            dualVertI++;
+        }
+        else
+        {
+            cit->cellIndex() = -1;
+        }
+    }
+
+    pts.setSize(dualVertI);
+}
+
+
 void Foam::conformalVoronoiMesh::reindexDualVertices
 (
     const Map<label>& dualPtIndexMap
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H
index 54775ed4a7..91ac79c57a 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H
@@ -36,6 +36,7 @@ Description
 
 #include <CGAL/Triangulation_3.h>
 #include "indexedVertex.H"
+#include "scalar.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -43,7 +44,7 @@ namespace CGAL
 {
 
 /*---------------------------------------------------------------------------*\
-                       Class indexedCell Declaration
+                          Class indexedCell Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class Gt, class Cb=CGAL::Triangulation_cell_base_3<Gt> >
@@ -54,20 +55,15 @@ class indexedCell
     // Private data
 
         //- The index for this tetrahedral cell
-        //  -1: tetrahedral and changed and associated data needs updating
-        // >=0: index of cells, set by external update algorithm
         int index_;
 
+        //- The applicable fraction of the face filtering mechanism
+        //  controls
+        Foam::scalar filterLimit_;
+
 
 public:
 
-    enum cellTypes
-    {
-        ctUnchanged           = 0,
-        ctChanged             = -1,
-        ctSaveChanged         = -2
-    };
-
     typedef typename Cb::Vertex_handle                   Vertex_handle;
     typedef typename Cb::Cell_handle                     Cell_handle;
 
@@ -82,18 +78,22 @@ public:
     indexedCell()
     :
         Cb(),
-        index_(ctChanged)
+        index_(-1),
+        filterLimit_(1.0)
     {}
 
+
     indexedCell
     (
         Vertex_handle v0, Vertex_handle v1, Vertex_handle v2, Vertex_handle v3
     )
     :
         Cb(v0, v1, v2, v3),
-        index_(ctChanged)
+        index_(-1),
+        filterLimit_(1.0)
     {}
 
+
     indexedCell
     (
         Vertex_handle v0,
@@ -107,39 +107,33 @@ public:
     )
     :
         Cb(v0, v1, v2, v3, n0, n1, n2, n3),
-        index_(ctChanged)
+        index_(-1),
+        filterLimit_(1.0)
     {}
 
-    void set_vertex(int i, Vertex_handle v)
-    {
-        index_ = ctChanged;
-        Cb::set_vertex(i, v);
-    }
-
-    void set_vertices()
-    {
-        index_ = ctChanged;
-        Cb::set_vertices();
-    }
-
-    void set_vertices
-    (
-        Vertex_handle v0, Vertex_handle v1, Vertex_handle v2, Vertex_handle v3
-    )
-    {
-        index_ = ctChanged;
-        Cb::set_vertices(v0, v1, v2, v3);
-    }
 
     int& cellIndex()
     {
         return index_;
     }
 
+
     int cellIndex() const
     {
         return index_;
     }
+
+
+    inline Foam::scalar& filterLimit()
+    {
+        return filterLimit_;
+    }
+
+
+    inline Foam::scalar filterLimit() const
+    {
+        return filterLimit_;
+    }
 };
 
 
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H
index 93fcb12346..a9116054cd 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H
@@ -35,7 +35,6 @@ Description
 #define indexedVertex_H
 
 #include <CGAL/Triangulation_3.h>
-#include "List.H"
 #include "tensor.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -44,7 +43,7 @@ namespace CGAL
 {
 
 /*---------------------------------------------------------------------------*\
-                       Class indexedVertex Declaration
+                         Class indexedVertex Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class Gt, class Vb=CGAL::Triangulation_vertex_base_3<Gt> >
@@ -91,6 +90,7 @@ public:
         typedef indexedVertex<Gt,Vb2>                           Other;
     };
 
+
     indexedVertex()
     :
         Vb(),
@@ -100,6 +100,7 @@ public:
         targetCellSize_(0.0)
     {}
 
+
     indexedVertex(const Point& p)
     :
         Vb(p),
@@ -109,6 +110,7 @@ public:
         targetCellSize_(0.0)
     {}
 
+
     indexedVertex(const Point& p, Cell_handle f)
     :
         Vb(f, p),
@@ -118,6 +120,7 @@ public:
         targetCellSize_(0.0)
     {}
 
+
     indexedVertex(Cell_handle f)
     :
         Vb(f),
@@ -133,70 +136,83 @@ public:
         return index_;
     }
 
+
     int index() const
     {
         return index_;
     }
 
+
     int& type()
     {
         return type_;
     }
 
+
     int type() const
     {
         return type_;
     }
 
+
     inline Foam::tensor& alignment()
     {
         return alignment_;
     }
 
+
     inline const Foam::tensor& alignment() const
     {
         return alignment_;
     }
 
+
     inline Foam::scalar& targetCellSize()
     {
         return targetCellSize_;
     }
 
+
     inline Foam::scalar targetCellSize() const
     {
         return targetCellSize_;
     }
 
+
     inline bool uninitialised() const
     {
         return type_ == ptInternalPoint && index_ == ptInternalPoint;
     }
 
+
     //- Is point a far-point
     inline bool farPoint() const
     {
         return type_ == ptFarPoint;
     }
 
+
     //- Is point internal, i.e. not on boundary
     inline bool internalPoint() const
     {
         return type_ <= ptInternalPoint;
     }
 
+
     //- Set the point to be internal
     inline void setInternal()
     {
         type_ = ptInternalPoint;
     }
 
+
     //- Is point internal and near the boundary
     inline bool nearBoundary() const
     {
         return type_ == ptNearBoundaryPoint;
     }
 
+
     //- Set the point to be near the boundary
     inline void setNearBoundary()
     {
@@ -209,12 +225,14 @@ public:
         return type_ == ptMirrorPoint;
     }
 
+
     //- Either master or slave of pointPair.
     inline bool pairPoint() const
     {
         return type_ >= 0;
     }
 
+
     //- Master of a pointPair is the lowest numbered one.
     inline bool ppMaster() const
     {
@@ -228,6 +246,7 @@ public:
         }
     }
 
+
     //- Slave of a pointPair is the highest numbered one.
     inline bool ppSlave() const
     {
@@ -241,18 +260,21 @@ public:
         }
     }
 
+
     //- Either original internal point or master of pointPair.
     inline bool internalOrBoundaryPoint() const
     {
         return internalPoint() || ppMaster();
     }
 
+
     //- Is point near the boundary or part of the boundary definition
     inline bool nearOrOnBoundary() const
     {
         return pairPoint() || mirrorPoint() || nearBoundary();
     }
 
+
     //- Do the two given vertices consitute a boundary point-pair
     inline friend bool pointPair
     (
@@ -263,6 +285,7 @@ public:
         return v0.index_ == v1.type_ || v1.index_ == v0.type_;
     }
 
+
     //- Do the three given vertices consitute a boundary triangle
     inline friend bool boundaryTriangle
     (
@@ -276,6 +299,7 @@ public:
             || (v2.pairPoint() && pointPair(v0, v1));
     }
 
+
     //- Do the three given vertices consitute an outside triangle
     inline friend bool outsideTriangle
     (
diff --git a/src/mesh/conformalVoronoiMesh/log b/src/mesh/conformalVoronoiMesh/log
deleted file mode 100644
index 0e31b724c8..0000000000
--- a/src/mesh/conformalVoronoiMesh/log
+++ /dev/null
@@ -1,20 +0,0 @@
-/*---------------------------------------------------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-Build  : dev-79ee164b7a49
-Exec   : cvMesh
-Date   : Jan 12 2010
-Time   : 19:49:15
-Host   : rockall
-PID    : 14798
-Case   : /home/graham/OpenFOAM/OpenFOAM-dev/src/mesh/conformalVoronoiMesh
-nProcs : 1
-SigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-Create time
-