diff --git a/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C b/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C
index 7a07b63d8f..0515205c78 100644
--- a/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C
+++ b/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C
@@ -95,24 +95,56 @@ int main(int argc, char *argv[])
     const bool collapseFaces = args.optionFound("collapseFaces");
     const bool collapseFaceZone = args.optionFound("collapseFaceZone");
 
+    if (collapseFaces && collapseFaceZone)
+    {
+        FatalErrorIn("main(int, char*[])")
+            << "Both face zone collapsing and face collapsing have been"
+            << "selected. Choose only one of:" << nl
+            << "    -collapseFaces" << nl
+            << "    -collapseFaceZone <faceZoneName>"
+            << abort(FatalError);
+    }
+
+    labelIOList pointPriority
+    (
+        IOobject
+        (
+            "pointPriority",
+            runTime.timeName(),
+            runTime,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        labelList(mesh.nPoints(), labelMin)
+    );
+
     forAll(timeDirs, timeI)
     {
         runTime.setTime(timeDirs[timeI], timeI);
 
         Info<< "Time = " << runTime.timeName() << endl;
 
-        polyMeshFilter meshFilter(mesh);
+        autoPtr<polyMeshFilter> meshFilterPtr;
 
-        // newMesh will be empty until it is filtered
-        const autoPtr<fvMesh>& newMesh = meshFilter.filteredMesh();
+        label nBadFaces = 0;
 
-        // Filter small edges only. This reduces the number of faces so that
-        // the face filtering is sped up.
-        label nBadFaces = meshFilter.filterEdges(0);
         {
-            polyTopoChange meshMod(newMesh);
+            meshFilterPtr.set(new polyMeshFilter(mesh, pointPriority));
+            polyMeshFilter& meshFilter = meshFilterPtr();
 
-            meshMod.changeMesh(mesh, false);
+            // newMesh will be empty until it is filtered
+            const autoPtr<fvMesh>& newMesh = meshFilter.filteredMesh();
+
+            // Filter small edges only. This reduces the number of faces so that
+            // the face filtering is sped up.
+            nBadFaces = meshFilter.filterEdges(0);
+            {
+                polyTopoChange meshMod(newMesh);
+
+                meshMod.changeMesh(mesh, false);
+            }
+
+            pointPriority = meshFilter.pointPriority();
         }
 
         if (collapseFaceZone)
@@ -121,18 +153,30 @@ int main(int argc, char *argv[])
 
             const faceZone& fZone = mesh.faceZones()[faceZoneName];
 
+            meshFilterPtr.reset(new polyMeshFilter(mesh, pointPriority));
+            polyMeshFilter& meshFilter = meshFilterPtr();
+
+            const autoPtr<fvMesh>& newMesh = meshFilter.filteredMesh();
+
             // Filter faces. Pass in the number of bad faces that are present
             // from the previous edge filtering to use as a stopping criterion.
-            meshFilter.filterFaceZone(fZone);
+            meshFilter.filter(fZone);
             {
                 polyTopoChange meshMod(newMesh);
 
                 meshMod.changeMesh(mesh, false);
             }
+
+            pointPriority = meshFilter.pointPriority();
         }
 
         if (collapseFaces)
         {
+            meshFilterPtr.reset(new polyMeshFilter(mesh, pointPriority));
+            polyMeshFilter& meshFilter = meshFilterPtr();
+
+            const autoPtr<fvMesh>& newMesh = meshFilter.filteredMesh();
+
             // Filter faces. Pass in the number of bad faces that are present
             // from the previous edge filtering to use as a stopping criterion.
             meshFilter.filter(nBadFaces);
@@ -141,6 +185,8 @@ int main(int argc, char *argv[])
 
                 meshMod.changeMesh(mesh, false);
             }
+
+            pointPriority = meshFilter.pointPriority();
         }
 
         // Write resulting mesh
@@ -157,6 +203,7 @@ int main(int argc, char *argv[])
             << runTime.timeName() << nl << endl;
 
         mesh.write();
+        pointPriority.write();
     }
 
     Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DistributedDelaunayMesh.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DistributedDelaunayMesh.C
index 152a82c7f2..40c76aec9d 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DistributedDelaunayMesh.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DistributedDelaunayMesh.C
@@ -966,6 +966,15 @@ Foam::DistributedDelaunayMesh<Triangulation>::rangeInsertReferredWithInfo
             )   << "Point is outside affine hull! pt = " << pointToInsert
                 << endl;
         }
+        else if (lt == Triangulation::OUTSIDE_CONVEX_HULL)
+        {
+            // @todo Can this be optimised?
+            //
+            // Only want to insert if a connection is formed between
+            // pointToInsert and an internal or internal boundary point.
+            hint = Triangulation::insert(pointToInsert, c);
+            inserted = true;
+        }
         else
         {
             // Get the cells that conflict with p in a vector V,
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index 8224492921..10ea65f772 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@@ -41,9 +41,26 @@ License
 
 namespace Foam
 {
-defineTypeNameAndDebug(conformalVoronoiMesh, 0);
+    defineTypeNameAndDebug(conformalVoronoiMesh, 0);
+
+    template<>
+    const char* NamedEnum
+    <
+        conformalVoronoiMesh::dualMeshPointType,
+        5
+    >::names[] =
+    {
+        "internal",
+        "surface",
+        "featureEdge",
+        "featurePoint",
+        "constrained"
+    };
 }
 
+const Foam::NamedEnum<Foam::conformalVoronoiMesh::dualMeshPointType, 5>
+    Foam::conformalVoronoiMesh::dualMeshPointTypeNames_;
+
 
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
@@ -1600,7 +1617,6 @@ void Foam::conformalVoronoiMesh::move()
         printVertexInfo(Info);
     }
 
-    // Write the intermediate mesh, do not filter the dual faces.
     if (time().outputTime())
     {
         writeMesh(time().timeName());
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
index eb75909d1c..334330f970 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
@@ -115,6 +115,19 @@ public:
     typedef List<pointIndexHitAndFeature>        pointIndexHitAndFeatureList;
     typedef DynamicList<pointIndexHitAndFeature> pointIndexHitAndFeatureDynList;
 
+    // Static data
+
+        enum dualMeshPointType
+        {
+            internal     = 0,
+            surface      = 1,
+            featureEdge  = 2,
+            featurePoint = 3,
+            constrained  = 4
+        };
+
+        static const NamedEnum<dualMeshPointType, 5> dualMeshPointTypeNames_;
+
 
 private:
 
@@ -682,14 +695,12 @@ private:
         //- Merge vertices that are identical
         void mergeIdenticalDualVertices
         (
-            const pointField& pts,
-            const labelList& boundaryPts
+            const pointField& pts
         );
 
         label mergeIdenticalDualVertices
         (
             const pointField& pts,
-            const labelList& boundaryPts,
             Map<label>& dualPtIndexMap
         ) const;
 
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
index 4e70aab964..cb17a1367d 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
@@ -207,209 +207,216 @@ void Foam::conformalVoronoiMesh::checkCells()
 }
 
 
-void Foam::conformalVoronoiMesh::checkDuals()
-{
-    List<List<Point> > pointFieldList(Pstream::nProcs());
-
-    List<Point> duals(number_of_finite_cells());
-
-//    PackedBoolList bPoints(number_of_finite_cells());
-
-//    indexDualVertices(duals, bPoints);
-
-    label count = 0;//duals.size();
-
-    duals.setSize(number_of_finite_cells());
-
-    globalIndex gIndex(number_of_vertices());
-
-    for
-    (
-        Delaunay::Finite_cells_iterator cit = finite_cells_begin();
-        cit != finite_cells_end();
-        ++cit
-    )
-    {
-        if (cit->hasFarPoint())
-        {
-            continue;
-        }
-
-        duals[count++] = cit->circumcenter();
-
-//        List<labelPair> cellVerticesPair(4);
-//        List<Point> cellVertices(4);
+//void Foam::conformalVoronoiMesh::checkDuals()
+//{
+//    List<List<Point> > pointFieldList(Pstream::nProcs());
 //
-//        for (label vI = 0; vI < 4; ++vI)
+//    List<Point> duals(number_of_finite_cells());
+//
+//    typedef CGAL::Exact_predicates_exact_constructions_kernel       EK2;
+//    typedef CGAL::Regular_triangulation_euclidean_traits_3<EK2>     EK;
+//    typedef CGAL::Cartesian_converter<baseK::Kernel, EK2>  To_exact;
+//    typedef CGAL::Cartesian_converter<EK2, baseK::Kernel>  Back_from_exact;
+//
+////    PackedBoolList bPoints(number_of_finite_cells());
+//
+////    indexDualVertices(duals, bPoints);
+//
+//    label count = 0;//duals.size();
+//
+//    duals.setSize(number_of_finite_cells());
+//
+//    globalIndex gIndex(number_of_vertices());
+//
+//    for
+//    (
+//        Delaunay::Finite_cells_iterator cit = finite_cells_begin();
+//        cit != finite_cells_end();
+//        ++cit
+//    )
+//    {
+//        if (cit->hasFarPoint())
 //        {
-//            cellVerticesPair[vI] = labelPair
-//            (
-//                cit->vertex(vI)->procIndex(),
-//                cit->vertex(vI)->index()
-//            );
-//            cellVertices[vI] = cit->vertex(vI)->point();
+//            continue;
 //        }
 //
-//        labelList oldToNew;
-//        sortedOrder(cellVerticesPair, oldToNew);
-//        oldToNew = invert(oldToNew.size(), oldToNew);
-//        inplaceReorder(oldToNew, cellVerticesPair);
-//        inplaceReorder(oldToNew, cellVertices);
+//        duals[count++] = cit->circumcenter();
 //
-//        duals[count++] = CGAL::circumcenter
-//        (
-//            cellVertices[0],
-//            cellVertices[1],
-//            cellVertices[2],
-//            cellVertices[3]
-//        );
-
-//        To_exact to_exact;
-//        Back_from_exact back_from_exact;
-//        EK::Construct_circumcenter_3 exact_circumcenter =
-//            EK().construct_circumcenter_3_object();
+////        List<labelPair> cellVerticesPair(4);
+////        List<Point> cellVertices(4);
+////
+////        for (label vI = 0; vI < 4; ++vI)
+////        {
+////            cellVerticesPair[vI] = labelPair
+////            (
+////                cit->vertex(vI)->procIndex(),
+////                cit->vertex(vI)->index()
+////            );
+////            cellVertices[vI] = cit->vertex(vI)->point();
+////        }
+////
+////        labelList oldToNew;
+////        sortedOrder(cellVerticesPair, oldToNew);
+////        oldToNew = invert(oldToNew.size(), oldToNew);
+////        inplaceReorder(oldToNew, cellVerticesPair);
+////        inplaceReorder(oldToNew, cellVertices);
+////
+////        duals[count++] = CGAL::circumcenter
+////        (
+////            cellVertices[0],
+////            cellVertices[1],
+////            cellVertices[2],
+////            cellVertices[3]
+////        );
 //
-//        duals[count++] = topoint
-//        (
-//            back_from_exact
-//            (
-//                exact_circumcenter
-//                (
-//                    to_exact(cit->vertex(0)->point()),
-//                    to_exact(cit->vertex(1)->point()),
-//                    to_exact(cit->vertex(2)->point()),
-//                    to_exact(cit->vertex(3)->point())
-//                )
-//            )
-//        );
-    }
-
-    Pout<< "Duals Calculated " << count << endl;
-
-    duals.setSize(count);
-
-    pointFieldList[Pstream::myProcNo()] = duals;
-
-    Pstream::gatherList(pointFieldList);
-
-    if (Pstream::master())
-    {
-        Info<< "Checking on master processor the dual locations of each " << nl
-            << "processor point list against the master dual list." << nl
-            << "There are " << pointFieldList.size() << " processors" << nl
-            << "The size of each processor's dual list is:" << endl;
-
-        forAll(pointFieldList, pfI)
-        {
-            Info<< "    Proc " << pfI << " has " << pointFieldList[pfI].size()
-                << " duals" << endl;
-        }
-
-        label nNonMatches = 0;
-        label nNearMatches = 0;
-        label nExactMatches = 0;
-
-        forAll(pointFieldList[0], pI)
-        {
-            const Point& masterPoint = pointFieldList[0][pI];
-
-            bool foundMatch = false;
-            bool foundNearMatch = false;
-
-            scalar minCloseness = GREAT;
-            Point closestPoint(0, 0, 0);
-
-            forAll(pointFieldList, pfI)
-            {
-                if (pfI == 0)
-                {
-                    continue;
-                }
-
-//                label pfI = 1;
-
-                forAll(pointFieldList[pfI], pISlave)
-                {
-                    const Point& slavePoint
-                        = pointFieldList[pfI][pISlave];
-
-                    if (masterPoint == slavePoint)
-                    {
-                        foundMatch = true;
-                        break;
-                    }
-
-                    const scalar closeness = mag
-                    (
-                        topoint(masterPoint) - topoint(slavePoint)
-                    );
-
-                    if (closeness < 1e-12)
-                    {
-                        foundNearMatch = true;
-                    }
-                    else
-                    {
-                        if (closeness < minCloseness)
-                        {
-                            minCloseness = closeness;
-                            closestPoint = slavePoint;
-                        }
-                    }
-                }
-
-                if (!foundMatch)
-                {
-                    if (foundNearMatch)
-                    {
-                        CGAL::Gmpq x(CGAL::to_double(masterPoint.x()));
-                        CGAL::Gmpq y(CGAL::to_double(masterPoint.y()));
-                        CGAL::Gmpq z(CGAL::to_double(masterPoint.z()));
-
-                        std::cout<< "master = " << x << " " << y << " " << z
-                            << std::endl;
-
-                        CGAL::Gmpq xs(CGAL::to_double(closestPoint.x()));
-                        CGAL::Gmpq ys(CGAL::to_double(closestPoint.y()));
-                        CGAL::Gmpq zs(CGAL::to_double(closestPoint.z()));
-                        std::cout<< "slave  = " << xs << " " << ys << " " << zs
-                            << std::endl;
-
-                        nNearMatches++;
-                    }
-                    else
-                    {
-                        nNonMatches++;
-                        Info<< "    Closest point to " << masterPoint << " is "
-                            << closestPoint << nl
-                            << "    Separation is " << minCloseness << endl;
-
-                        CGAL::Gmpq x(CGAL::to_double(masterPoint.x()));
-                        CGAL::Gmpq y(CGAL::to_double(masterPoint.y()));
-                        CGAL::Gmpq z(CGAL::to_double(masterPoint.z()));
-
-                        std::cout<< "master = " << x << " " << y << " " << z
-                                 << std::endl;
-
-                        CGAL::Gmpq xs(CGAL::to_double(closestPoint.x()));
-                        CGAL::Gmpq ys(CGAL::to_double(closestPoint.y()));
-                        CGAL::Gmpq zs(CGAL::to_double(closestPoint.z()));
-                        std::cout<< "slave  = " << xs << " " << ys << " " << zs
-                                 << std::endl;
-                    }
-                }
-                else
-                {
-                    nExactMatches++;
-                }
-            }
-        }
-
-        Info<< "Found " << nNonMatches << " non-matching duals" << nl
-            << " and " << nNearMatches << " near matches"
-            << " and " << nExactMatches << " exact matches" << endl;
-    }
-}
+////        To_exact to_exact;
+////        Back_from_exact back_from_exact;
+////        EK::Construct_circumcenter_3 exact_circumcenter =
+////            EK().construct_circumcenter_3_object();
+////
+////        duals[count++] = topoint
+////        (
+////            back_from_exact
+////            (
+////                exact_circumcenter
+////                (
+////                    to_exact(cit->vertex(0)->point()),
+////                    to_exact(cit->vertex(1)->point()),
+////                    to_exact(cit->vertex(2)->point()),
+////                    to_exact(cit->vertex(3)->point())
+////                )
+////            )
+////        );
+//    }
+//
+//    Pout<< "Duals Calculated " << count << endl;
+//
+//    duals.setSize(count);
+//
+//    pointFieldList[Pstream::myProcNo()] = duals;
+//
+//    Pstream::gatherList(pointFieldList);
+//
+//    if (Pstream::master())
+//    {
+//        Info<< "Checking on master processor the dual locations of each" << nl
+//            << " processor point list against the master dual list." << nl
+//            << "There are " << pointFieldList.size() << " processors" << nl
+//            << "The size of each processor's dual list is:" << endl;
+//
+//        forAll(pointFieldList, pfI)
+//        {
+//            Info<< "    Proc " << pfI << " has " << pointFieldList[pfI].size()
+//                << " duals" << endl;
+//        }
+//
+//        label nNonMatches = 0;
+//        label nNearMatches = 0;
+//        label nExactMatches = 0;
+//
+//        forAll(pointFieldList[0], pI)
+//        {
+//            const Point& masterPoint = pointFieldList[0][pI];
+//
+//            bool foundMatch = false;
+//            bool foundNearMatch = false;
+//
+//            scalar minCloseness = GREAT;
+//            Point closestPoint(0, 0, 0);
+//
+//            forAll(pointFieldList, pfI)
+//            {
+//                if (pfI == 0)
+//                {
+//                    continue;
+//                }
+//
+////                label pfI = 1;
+//
+//                forAll(pointFieldList[pfI], pISlave)
+//                {
+//                    const Point& slavePoint
+//                        = pointFieldList[pfI][pISlave];
+//
+//                    if (masterPoint == slavePoint)
+//                    {
+//                        foundMatch = true;
+//                        break;
+//                    }
+//
+//                    const scalar closeness = mag
+//                    (
+//                        topoint(masterPoint) - topoint(slavePoint)
+//                    );
+//
+//                    if (closeness < 1e-12)
+//                    {
+//                        foundNearMatch = true;
+//                    }
+//                    else
+//                    {
+//                        if (closeness < minCloseness)
+//                        {
+//                            minCloseness = closeness;
+//                            closestPoint = slavePoint;
+//                        }
+//                    }
+//                }
+//
+//                if (!foundMatch)
+//                {
+//                    if (foundNearMatch)
+//                    {
+//                        CGAL::Gmpq x(CGAL::to_double(masterPoint.x()));
+//                        CGAL::Gmpq y(CGAL::to_double(masterPoint.y()));
+//                        CGAL::Gmpq z(CGAL::to_double(masterPoint.z()));
+//
+//                        std::cout<< "master = " << x << " " << y << " " << z
+//                            << std::endl;
+//
+//                        CGAL::Gmpq xs(CGAL::to_double(closestPoint.x()));
+//                        CGAL::Gmpq ys(CGAL::to_double(closestPoint.y()));
+//                        CGAL::Gmpq zs(CGAL::to_double(closestPoint.z()));
+//                        std::cout<< "slave  = " << xs << " " << ys << " "
+//                                 << zs
+//                            << std::endl;
+//
+//                        nNearMatches++;
+//                    }
+//                    else
+//                    {
+//                        nNonMatches++;
+//                        Info<< "Closest point to " << masterPoint << " is "
+//                            << closestPoint << nl
+//                            << "    Separation is " << minCloseness << endl;
+//
+//                        CGAL::Gmpq x(CGAL::to_double(masterPoint.x()));
+//                        CGAL::Gmpq y(CGAL::to_double(masterPoint.y()));
+//                        CGAL::Gmpq z(CGAL::to_double(masterPoint.z()));
+//
+//                        std::cout<< "master = " << x << " " << y << " " << z
+//                                 << std::endl;
+//
+//                        CGAL::Gmpq xs(CGAL::to_double(closestPoint.x()));
+//                        CGAL::Gmpq ys(CGAL::to_double(closestPoint.y()));
+//                        CGAL::Gmpq zs(CGAL::to_double(closestPoint.z()));
+//                        std::cout<< "slave  = " << xs << " " << ys << " "
+//                                 << zs
+//                                 << std::endl;
+//                    }
+//                }
+//                else
+//                {
+//                    nExactMatches++;
+//                }
+//            }
+//        }
+//
+//        Info<< "Found " << nNonMatches << " non-matching duals" << nl
+//            << " and " << nNearMatches << " near matches"
+//            << " and " << nExactMatches << " exact matches" << endl;
+//    }
+//}
 
 
 void Foam::conformalVoronoiMesh::checkVertices()
@@ -578,7 +585,7 @@ void Foam::conformalVoronoiMesh::calcDualMesh
         Info<< nl << "Merging identical points" << endl;
 
         // There is no guarantee that a merge of close points is no-risk
-        mergeIdenticalDualVertices(points, boundaryPts);
+        mergeIdenticalDualVertices(points);
     }
 
     // Final dual face and owner neighbour construction
@@ -813,8 +820,7 @@ void Foam::conformalVoronoiMesh::calcTetMesh
 
 void Foam::conformalVoronoiMesh::mergeIdenticalDualVertices
 (
-    const pointField& pts,
-    const labelList& boundaryPts
+    const pointField& pts
 )
 {
     // Assess close points to be merged
@@ -829,7 +835,6 @@ void Foam::conformalVoronoiMesh::mergeIdenticalDualVertices
         nPtsMerged = mergeIdenticalDualVertices
         (
             pts,
-            boundaryPts,
             dualPtIndexMap
         );
 
@@ -851,7 +856,6 @@ void Foam::conformalVoronoiMesh::mergeIdenticalDualVertices
 Foam::label Foam::conformalVoronoiMesh::mergeIdenticalDualVertices
 (
     const pointField& pts,
-    const labelList& boundaryPts,
     Map<label>& dualPtIndexMap
 ) const
 {
@@ -883,6 +887,19 @@ Foam::label Foam::conformalVoronoiMesh::mergeIdenticalDualVertices
 
             if (p1 == p2)
             {
+//                if (c1->parallelDualVertex() || c2->parallelDualVertex())
+//                {
+//                    if (c1->vertexLowestProc() < c2->vertexLowestProc())
+//                    {
+//                        dualPtIndexMap.insert(c1I, c1I);
+//                        dualPtIndexMap.insert(c2I, c1I);
+//                    }
+//                    else
+//                    {
+//                        dualPtIndexMap.insert(c1I, c2I);
+//                        dualPtIndexMap.insert(c2I, c2I);
+//                    }
+//                }
                 if (c1I < c2I)
                 {
                     dualPtIndexMap.insert(c1I, c1I);
@@ -1338,13 +1355,13 @@ void Foam::conformalVoronoiMesh::checkCellSizing()
 
     timeCheck("Start of Cell Sizing");
 
-    labelList boundaryPts(number_of_finite_cells(), -1);
+    labelList boundaryPts(number_of_finite_cells(), internal);
     pointField ptsField;
 
     indexDualVertices(ptsField, boundaryPts);
 
     // Merge close dual vertices.
-    mergeIdenticalDualVertices(ptsField, boundaryPts);
+    mergeIdenticalDualVertices(ptsField);
 
     autoPtr<polyMesh> meshPtr = createPolyMeshFromPoints(ptsField);
     const polyMesh& pMesh = meshPtr();
@@ -1755,7 +1772,7 @@ void Foam::conformalVoronoiMesh::indexDualVertices
     boundaryPts.setSize
     (
         number_of_finite_cells() + nConstrainedVertices,
-        -1
+        internal
     );
 
     if (foamyHexMeshControls().guardFeaturePoints())
@@ -1774,7 +1791,7 @@ void Foam::conformalVoronoiMesh::indexDualVertices
                     topoint(vit->point());
 
                 boundaryPts[number_of_finite_cells() + nConstrainedVertices] =
-                    1;
+                    constrained;
 
                 nConstrainedVertices++;
             }
@@ -1974,15 +1991,40 @@ void Foam::conformalVoronoiMesh::indexDualVertices
 
             if (cit->boundaryDualVertex())
             {
-                if (cit->featureEdgeDualVertex())
+                if (cit->featurePointDualVertex())
                 {
-                    boundaryPts[cit->cellIndex()] = 1;
+                    boundaryPts[cit->cellIndex()] = featurePoint;
+                }
+                else if (cit->featureEdgeDualVertex())
+                {
+                    boundaryPts[cit->cellIndex()] = featureEdge;
                 }
                 else
                 {
-                    boundaryPts[cit->cellIndex()] = 0;
+                    boundaryPts[cit->cellIndex()] = surface;
                 }
             }
+            else if
+            (
+                cit->baffleBoundaryDualVertex()
+            )
+            {
+                boundaryPts[cit->cellIndex()] = surface;
+            }
+            else if
+            (
+                cit->vertex(0)->featureEdgePoint()
+             && cit->vertex(1)->featureEdgePoint()
+             && cit->vertex(2)->featureEdgePoint()
+             && cit->vertex(3)->featureEdgePoint()
+            )
+            {
+                boundaryPts[cit->cellIndex()] = featureEdge;
+            }
+            else
+            {
+                boundaryPts[cit->cellIndex()] = internal;
+            }
         }
         else
         {
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
index 4f02fb67f6..d733b33e35 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
@@ -958,7 +958,7 @@ void Foam::conformalVoronoiMesh::writeMesh
     {
         Info<< nl << "Filtering edges on polyMesh" << nl << endl;
 
-        meshFilter.reset(new polyMeshFilter(mesh));
+        meshFilter.reset(new polyMeshFilter(mesh, boundaryPts));
 
         // Filter small edges only. This reduces the number of faces so that
         // the face filtering is sped up.
@@ -974,9 +974,28 @@ void Foam::conformalVoronoiMesh::writeMesh
 
     if (foamyHexMeshControls().filterFaces())
     {
+        labelIOList boundaryPtsIO
+        (
+            IOobject
+            (
+                "pointPriority",
+                instance,
+                time(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            labelList(mesh.nPoints(), labelMin)
+        );
+
+        forAll(mesh.points(), ptI)
+        {
+            boundaryPtsIO[ptI] = mesh.pointZones().whichZone(ptI);
+        }
+
+
         Info<< nl << "Filtering faces on polyMesh" << nl << endl;
 
-        meshFilter.reset(new polyMeshFilter(mesh));
+        meshFilter.reset(new polyMeshFilter(mesh, boundaryPtsIO));
 
         meshFilter().filter(nInitialBadFaces);
         {
@@ -1005,159 +1024,44 @@ void Foam::conformalVoronoiMesh::writeMesh
             << endl;
     }
 
-
-//    volTensorField alignments
-//    (
-//        IOobject
-//        (
-//            "alignmentsField",
-//            runTime_.timeName(),
-//            runTime_,
-//            IOobject::NO_READ,
-//            IOobject::AUTO_WRITE
-//        ),
-//        mesh,
-//        tensor::zero
-//    );
-//
-//    forAll(mesh.cellCentres(), pI)
-//    {
-//        Vertex_handle nearV =
-//            nearest_vertex
-//            (
-//                toPoint<Point>(mesh.cellCentres()[pI])
-//            );
-//        alignments[pI] = nearV->alignment();
-//    }
-//    alignments.write();
-//
-//    {
-//        volVectorField alignmentx
-//        (
-//            IOobject
-//            (
-//                "alignmentsx",
-//                runTime_.timeName(),
-//                runTime_,
-//                IOobject::NO_READ,
-//                IOobject::AUTO_WRITE
-//            ),
-//            mesh,
-//            vector::zero
-//        );
-//        forAll(alignmentx, aI)
-//        {
-//            alignmentx[aI] = alignments[aI].x();
-//        }
-//        alignmentx.write();
-//    }
-//    {
-//        volVectorField alignmenty
-//        (
-//            IOobject
-//            (
-//                "alignmentsy",
-//                runTime_.timeName(),
-//                runTime_,
-//                IOobject::NO_READ,
-//                IOobject::AUTO_WRITE
-//            ),
-//            mesh,
-//            vector::zero
-//        );
-//        forAll(alignmenty, aI)
-//        {
-//            alignmenty[aI] = alignments[aI].y();
-//        }
-//        alignmenty.write();
-//    }
-//    {
-//        volVectorField alignmentz
-//        (
-//            IOobject
-//            (
-//                "alignmentsz",
-//                runTime_.timeName(),
-//                runTime_,
-//                IOobject::NO_READ,
-//                IOobject::AUTO_WRITE
-//            ),
-//            mesh,
-//            vector::zero
-//        );
-//        forAll(alignmentz, aI)
-//        {
-//            alignmentz[aI] = alignments[aI].z();
-//        }
-//        alignmentz.write();
-//    }
-
-
-    labelIOList boundaryIOPts
-    (
-        IOobject
-        (
-            "boundaryPoints",
-            instance,
-            runTime_,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        boundaryPts
-    );
-
-    // Dump list of boundary points
-    forAll(mesh.boundaryMesh(), patchI)
     {
-        const polyPatch& pp = mesh.boundaryMesh()[patchI];
+        pointScalarField boundaryPtsScalarField
+        (
+            IOobject
+            (
+                "boundaryPoints_collapsed",
+                instance,
+                time(),
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            pointMesh::New(mesh),
+            scalar(labelMin)
+        );
 
-        if (!isA<coupledPolyPatch>(pp))
+        labelIOList boundaryPtsIO
+        (
+            IOobject
+            (
+                "pointPriority",
+                instance,
+                time(),
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            labelList(mesh.nPoints(), labelMin)
+        );
+
+        forAll(mesh.points(), ptI)
         {
-            forAll(pp, fI)
-            {
-                const face& boundaryFace = pp[fI];
-
-                forAll(boundaryFace, pI)
-                {
-                    const label boundaryPointI = boundaryFace[pI];
-
-                    boundaryIOPts[boundaryPointI] = boundaryPts[boundaryPointI];
-                }
-            }
+            boundaryPtsScalarField[ptI] = mesh.pointZones().whichZone(ptI);
+            boundaryPtsIO[ptI] = mesh.pointZones().whichZone(ptI);
         }
+
+        boundaryPtsScalarField.write();
+        boundaryPtsIO.write();
     }
 
-    boundaryIOPts.write();
-
-//    forAllConstIter(labelHashSet, pointsInPatch, pI)
-//    {
-//        const Foam::point& ptMaster = mesh.points()[pI.key()];
-//
-//        forAllConstIter(labelHashSet, pointsInPatch, ptI)
-//        {
-//            if (ptI.key() != pI.key())
-//            {
-//                const Foam::point& ptSlave = mesh.points()[ptI.key()];
-//
-//                const scalar dist = mag(ptMaster - ptSlave);
-//                if (ptMaster == ptSlave)
-//                {
-//                    Pout<< "Point(" << pI.key() << ") " << ptMaster
-//                        << " == "
-//                        << "(" << ptI.key() << ") " << ptSlave
-//                        << endl;
-//                }
-//                else if (dist == 0)
-//                {
-//                    Pout<< "Point(" << pI.key() << ") " << ptMaster
-//                        << " ~= "
-//                        << "(" << ptI.key() << ") " << ptSlave
-//                        << endl;
-//                }
-//            }
-//        }
-//    }
-
 //    writeCellSizes(mesh);
 
 //    writeCellAlignments(mesh);
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCell.H b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCell.H
index 82bbe43b95..15b155d85b 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCell.H
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCell.H
@@ -176,6 +176,9 @@ public:
         //- Does the Delaunay cell have a far point
         inline bool hasFarPoint() const;
 
+        //- Does the Delaunay cell have a referred point
+        inline bool hasReferredPoint() const;
+
         //- Does the Delaunay cell have a feature point
         inline bool hasFeaturePoint() const;
 
@@ -216,6 +219,8 @@ public:
         //  least one Delaunay vertex outside and at least one inside
         inline bool boundaryDualVertex() const;
 
+        inline bool baffleBoundaryDualVertex() const;
+
         //- A dual vertex on a feature edge will result from this Delaunay cell
         inline bool featureEdgeDualVertex() const;
 
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCellI.H b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCellI.H
index 2f5d076138..6c55a8eb30 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCellI.H
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCellI.H
@@ -189,6 +189,19 @@ inline bool CGAL::indexedCell<Gt, Cb>::hasFarPoint() const
 }
 
 
+template<class Gt, class Cb>
+inline bool CGAL::indexedCell<Gt, Cb>::hasReferredPoint() const
+{
+    return
+    (
+        this->vertex(0)->referred()
+     || this->vertex(1)->referred()
+     || this->vertex(2)->referred()
+     || this->vertex(3)->referred()
+    );
+}
+
+
 template<class Gt, class Cb>
 inline bool CGAL::indexedCell<Gt, Cb>::hasFeaturePoint() const
 {
@@ -372,6 +385,14 @@ inline bool CGAL::indexedCell<Gt, Cb>::anyInternalOrBoundaryDualVertex() const
 template<class Gt, class Cb>
 inline bool CGAL::indexedCell<Gt, Cb>::boundaryDualVertex() const
 {
+//    return
+//    (
+//        this->vertex(0)->boundaryPoint()
+//     && this->vertex(1)->boundaryPoint()
+//     && this->vertex(2)->boundaryPoint()
+//     && this->vertex(3)->boundaryPoint()
+//    );
+
     return
     (
         (
@@ -387,6 +408,41 @@ inline bool CGAL::indexedCell<Gt, Cb>::boundaryDualVertex() const
         || this->vertex(3)->externalBoundaryPoint()
         )
     );
+
+//    Foam::label nBoundaryPoints = 0;
+//
+//    for (Foam::label i = 0; i < 4; ++i)
+//    {
+//        Vertex_handle v = this->vertex(i);
+//
+//        if (v->boundaryPoint())
+//        {
+//            nBoundaryPoints++;
+//        }
+//    }
+//
+//    return (nBoundaryPoints > 1);
+}
+
+
+template<class Gt, class Cb>
+inline bool CGAL::indexedCell<Gt, Cb>::baffleBoundaryDualVertex() const
+{
+    return
+    (
+        (
+           this->vertex(0)->internalBafflePoint()
+        || this->vertex(1)->internalBafflePoint()
+        || this->vertex(2)->internalBafflePoint()
+        || this->vertex(3)->internalBafflePoint()
+        )
+     && (
+           this->vertex(0)->externalBafflePoint()
+        || this->vertex(1)->externalBafflePoint()
+        || this->vertex(2)->externalBafflePoint()
+        || this->vertex(3)->externalBafflePoint()
+        )
+    );
 }
 
 
@@ -400,6 +456,31 @@ inline bool CGAL::indexedCell<Gt, Cb>::featureEdgeDualVertex() const
      && this->vertex(2)->featureEdgePoint()
      && this->vertex(3)->featureEdgePoint()
     );
+//        (
+//            this->vertex(0)->featureEdgePoint()
+//         || this->vertex(1)->featureEdgePoint()
+//         || this->vertex(2)->featureEdgePoint()
+//         || this->vertex(3)->featureEdgePoint()
+//        )
+//     && (
+//            (
+//                this->vertex(0)->featureEdgePoint()
+//             || this->vertex(0)->featurePoint()
+//            )
+//         && (
+//                this->vertex(1)->featureEdgePoint()
+//             || this->vertex(1)->featurePoint()
+//            )
+//         && (
+//                this->vertex(2)->featureEdgePoint()
+//             || this->vertex(2)->featurePoint()
+//            )
+//         && (
+//                this->vertex(3)->featureEdgePoint()
+//             || this->vertex(3)->featurePoint()
+//            )
+//        )
+//    );
 }
 
 
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex/indexedVertex.H b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex/indexedVertex.H
index b556be66cc..362e457d74 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex/indexedVertex.H
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex/indexedVertex.H
@@ -243,8 +243,10 @@ public:
         inline bool surfacePoint() const;
 
         inline bool internalBoundaryPoint() const;
+        inline bool internalBafflePoint() const;
 
         inline bool externalBoundaryPoint() const;
+        inline bool externalBafflePoint() const;
 
         inline bool constrained() const;
 
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex/indexedVertexEnum.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex/indexedVertexEnum.C
index d0ac971782..34c933e162 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex/indexedVertexEnum.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex/indexedVertexEnum.C
@@ -30,13 +30,17 @@ License
 
 template<>
 const char*
-Foam::NamedEnum<Foam::indexedVertexEnum::vertexType, 11>::names[] =
+Foam::NamedEnum<Foam::indexedVertexEnum::vertexType, 15>::names[] =
 {
     "Unassigned",
     "Internal",
     "InternalNearBoundary",
     "InternalSurface",
+    "InternalSurfaceBaffle",
+    "ExternalSurfaceBaffle",
     "InternalFeatureEdge",
+    "InternalFeatureEdgeBaffle",
+    "ExternalFeatureEdgeBaffle",
     "InternalFeaturePoint",
     "ExternalSurface",
     "ExternalFeatureEdge",
@@ -45,7 +49,7 @@ Foam::NamedEnum<Foam::indexedVertexEnum::vertexType, 11>::names[] =
     "Constrained"
 };
 
-const Foam::NamedEnum<Foam::indexedVertexEnum::vertexType, 11>
+const Foam::NamedEnum<Foam::indexedVertexEnum::vertexType, 15>
 Foam::indexedVertexEnum::vertexTypeNames_;
 
 
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex/indexedVertexEnum.H b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex/indexedVertexEnum.H
index 072470d673..c606f815d2 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex/indexedVertexEnum.H
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex/indexedVertexEnum.H
@@ -49,17 +49,21 @@ public:
 
     enum vertexType
     {
-        vtUnassigned            = 0,
-        vtInternal              = 1,
-        vtInternalNearBoundary  = 2,
-        vtInternalSurface       = 3,
-        vtInternalFeatureEdge   = 4,
-        vtInternalFeaturePoint  = 5,
-        vtExternalSurface       = 6,
-        vtExternalFeatureEdge   = 7,
-        vtExternalFeaturePoint  = 8,
-        vtFar                   = 9,
-        vtConstrained           = 10
+        vtUnassigned                = 0,
+        vtInternal                  = 1,
+        vtInternalNearBoundary      = 2,
+        vtInternalSurface           = 3,
+        vtInternalSurfaceBaffle     = 4,
+        vtExternalSurfaceBaffle     = 5,
+        vtInternalFeatureEdge       = 6,
+        vtInternalFeatureEdgeBaffle = 7,
+        vtExternalFeatureEdgeBaffle = 8,
+        vtInternalFeaturePoint      = 9,
+        vtExternalSurface           = 10,
+        vtExternalFeatureEdge       = 11,
+        vtExternalFeaturePoint      = 12,
+        vtFar                       = 13,
+        vtConstrained               = 14
     };
 
     enum vertexMotion
@@ -68,7 +72,7 @@ public:
         movable = 1
     };
 
-    static const Foam::NamedEnum<vertexType, 11> vertexTypeNames_;
+    static const Foam::NamedEnum<vertexType, 15> vertexTypeNames_;
 
     static const Foam::NamedEnum<vertexMotion, 2> vertexMotionNames_;
 
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex/indexedVertexI.H b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex/indexedVertexI.H
index 76a77da47d..c5930395c3 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex/indexedVertexI.H
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex/indexedVertexI.H
@@ -307,6 +307,16 @@ inline bool CGAL::indexedVertex<Gt, Vb>::internalBoundaryPoint() const
     return type_ >= vtInternalSurface && type_ <= vtInternalFeaturePoint;
 }
 
+template<class Gt, class Vb>
+inline bool CGAL::indexedVertex<Gt, Vb>::internalBafflePoint() const
+{
+    return
+    (
+        type_ == vtInternalSurfaceBaffle
+     || type_ == vtInternalFeatureEdgeBaffle
+    );
+}
+
 
 template<class Gt, class Vb>
 inline bool CGAL::indexedVertex<Gt, Vb>::externalBoundaryPoint() const
@@ -314,6 +324,16 @@ inline bool CGAL::indexedVertex<Gt, Vb>::externalBoundaryPoint() const
     return type_ >= vtExternalSurface && type_ <= vtExternalFeaturePoint;
 }
 
+template<class Gt, class Vb>
+inline bool CGAL::indexedVertex<Gt, Vb>::externalBafflePoint() const
+{
+    return
+    (
+        type_ == vtExternalSurfaceBaffle
+     || type_ == vtExternalFeatureEdgeBaffle
+    );
+}
+
 
 template<class Gt, class Vb>
 inline bool CGAL::indexedVertex<Gt, Vb>::constrained() const
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C
index 66dc10737d..6568ef1d02 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C
@@ -49,7 +49,6 @@ void rayShooting::splitLine
 {
     Foam::point midPoint(l.centre());
     const scalar localCellSize(cellShapeControls().cellSize(midPoint));
-    const scalar lineLength(l.mag());
 
     const scalar minDistFromSurfaceSqr
     (
@@ -64,6 +63,8 @@ void rayShooting::splitLine
     )
     {
         // Add extra points if line length is much bigger than local cell size
+//        const scalar lineLength(l.mag());
+//
 //        if (lineLength > 4.0*localCellSize)
 //        {
 //            splitLine
diff --git a/src/OpenFOAM/db/error/IOerror.C b/src/OpenFOAM/db/error/IOerror.C
index f68fd651f7..08f490c0bd 100644
--- a/src/OpenFOAM/db/error/IOerror.C
+++ b/src/OpenFOAM/db/error/IOerror.C
@@ -255,28 +255,31 @@ void Foam::IOerror::abort()
 
 Foam::Ostream& Foam::operator<<(Ostream& os, const IOerror& ioErr)
 {
-    os  << endl
-        << ioErr.title().c_str() << endl
-        << ioErr.message().c_str() << endl << endl;
-
-    os  << "file: " << ioErr.ioFileName().c_str();
-
-    if (ioErr.ioStartLineNumber() >= 0 && ioErr.ioEndLineNumber() >= 0)
+    if (!os.bad())
     {
-        os  << " from line " << ioErr.ioStartLineNumber()
-            << " to line " << ioErr.ioEndLineNumber() << '.';
-    }
-    else if (ioErr.ioStartLineNumber() >= 0)
-    {
-        os  << " at line " << ioErr.ioStartLineNumber() << '.';
-    }
+        os  << endl
+            << ioErr.title().c_str() << endl
+            << ioErr.message().c_str() << endl << endl;
 
-    if (IOerror::level >= 2 && ioErr.sourceFileLineNumber())
-    {
-        os  << endl << endl
-            << "    From function " << ioErr.functionName().c_str() << endl
-            << "    in file " << ioErr.sourceFileName().c_str()
-            << " at line " << ioErr.sourceFileLineNumber() << '.';
+        os  << "file: " << ioErr.ioFileName().c_str();
+
+        if (ioErr.ioStartLineNumber() >= 0 && ioErr.ioEndLineNumber() >= 0)
+        {
+            os  << " from line " << ioErr.ioStartLineNumber()
+                << " to line " << ioErr.ioEndLineNumber() << '.';
+        }
+        else if (ioErr.ioStartLineNumber() >= 0)
+        {
+            os  << " at line " << ioErr.ioStartLineNumber() << '.';
+        }
+
+        if (IOerror::level >= 2 && ioErr.sourceFileLineNumber())
+        {
+            os  << endl << endl
+                << "    From function " << ioErr.functionName().c_str() << endl
+                << "    in file " << ioErr.sourceFileName().c_str()
+                << " at line " << ioErr.sourceFileLineNumber() << '.';
+        }
     }
 
     return os;
diff --git a/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.C b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.C
index 0968c71ca3..1cec955dbd 100644
--- a/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.C
+++ b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.C
@@ -269,7 +269,8 @@ bool Foam::treeBoundBox::intersects
     const direction endBits = posBits(end);
     pt = start;
 
-    while (true)
+    // Allow maximum of 3 clips.
+    for (label i = 0; i < 4; ++i)
     {
         direction ptBits = posBits(pt);
 
@@ -380,6 +381,9 @@ bool Foam::treeBoundBox::intersects
             }
         }
     }
+
+    // Can end up here if the end point is on the edge of the boundBox
+    return true;
 }
 
 
diff --git a/src/dynamicMesh/Make/files b/src/dynamicMesh/Make/files
index 6429d3b4a5..f769607f94 100644
--- a/src/dynamicMesh/Make/files
+++ b/src/dynamicMesh/Make/files
@@ -99,6 +99,7 @@ createShellMesh/createShellMesh.C
 
 extrudePatchMesh/extrudePatchMesh.C
 
+polyMeshFilter/polyMeshFilterSettings.C
 polyMeshFilter/polyMeshFilter.C
 
 LIB = $(FOAM_LIBBIN)/libdynamicMesh
diff --git a/src/dynamicMesh/polyMeshFilter/polyMeshFilter.C b/src/dynamicMesh/polyMeshFilter/polyMeshFilter.C
index 0c18d84e4f..2958356dc9 100644
--- a/src/dynamicMesh/polyMeshFilter/polyMeshFilter.C
+++ b/src/dynamicMesh/polyMeshFilter/polyMeshFilter.C
@@ -78,6 +78,424 @@ Foam::autoPtr<Foam::fvMesh> Foam::polyMeshFilter::copyMesh(const fvMesh& mesh)
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+Foam::label Foam::polyMeshFilter::filterFacesLoop(const label nOriginalBadFaces)
+{
+    label nBadFaces = labelMax;
+    label nOuterIterations = 0;
+
+    // Maintain the number of times a point has been part of a bad face
+    labelList pointErrorCount(mesh_.nPoints(), 0);
+
+    PackedBoolList newErrorPoint(mesh_.nPoints());
+    edgeCollapser::checkMeshQuality
+    (
+        mesh_,
+        meshQualityCoeffDict(),
+        newErrorPoint
+    );
+
+    bool newBadFaces = true;
+
+    // Main loop
+    // ~~~~~~~~~
+    // It tries and do some collapses, checks the resulting mesh and
+    // 'freezes' some edges (by marking in minEdgeLen) and tries again.
+    // This will iterate ultimately to the situation where every edge is
+    // frozen and nothing gets collapsed.
+    while
+    (
+        nOuterIterations < maxIterations()
+     //&& nBadFaces > nOriginalBadFaces
+     && newBadFaces
+    )
+    {
+        Info<< nl << "Outer Iteration = " << nOuterIterations++ << nl
+            << endl;
+
+        printScalarFieldStats("Edge Filter Factor", minEdgeLen_);
+        printScalarFieldStats("Face Filter Factor", faceFilterFactor_);
+
+        // Reset the new mesh to the old mesh
+        newMeshPtr_ = copyMesh(mesh_);
+        fvMesh& newMesh = newMeshPtr_();
+
+        scalarField newMeshFaceFilterFactor = faceFilterFactor_;
+        pointPriority_.reset(new labelList(originalPointPriority_));
+
+        labelList origToCurrentPointMap(identity(newMesh.nPoints()));
+
+        {
+            label nInnerIterations = 0;
+            label nPrevLocalCollapse = labelMax;
+
+            Info<< incrIndent;
+
+            while (true)
+            {
+                Info<< nl << indent << "Inner iteration = "
+                    << nInnerIterations++ << nl << incrIndent << endl;
+
+                label nLocalCollapse = filterFaces
+                (
+                    newMesh,
+                    newMeshFaceFilterFactor,
+                    origToCurrentPointMap
+                );
+                Info<< decrIndent;
+
+                if
+                (
+                    nLocalCollapse >= nPrevLocalCollapse
+                 || nLocalCollapse == 0
+                )
+                {
+                    Info<< decrIndent;
+                    break;
+                }
+                else
+                {
+                    nPrevLocalCollapse = nLocalCollapse;
+                }
+            }
+        }
+
+        scalarField newMeshMinEdgeLen = minEdgeLen_;
+
+        {
+            label nInnerIterations = 0;
+            label nPrevLocalCollapse = labelMax;
+
+            Info<< incrIndent;
+
+            while (true)
+            {
+                Info<< nl << indent << "Inner iteration = "
+                    << nInnerIterations++ << nl << incrIndent << endl;
+
+                label nLocalCollapse = filterEdges
+                (
+                    newMesh,
+                    newMeshMinEdgeLen,
+                    origToCurrentPointMap
+                );
+                Info<< decrIndent;
+
+                if
+                (
+                    nLocalCollapse >= nPrevLocalCollapse
+                 || nLocalCollapse == 0
+                )
+                {
+                    Info<< decrIndent;
+                    break;
+                }
+                else
+                {
+                    nPrevLocalCollapse = nLocalCollapse;
+                }
+            }
+        }
+
+
+        // Mesh check
+        // ~~~~~~~~~~~~~~~~~~
+        // Do not allow collapses in regions of error.
+        // Updates minEdgeLen, nRelaxedEdges
+
+        if (controlMeshQuality())
+        {
+            PackedBoolList isErrorPoint(newMesh.nPoints());
+            nBadFaces = edgeCollapser::checkMeshQuality
+            (
+                newMesh,
+                meshQualityCoeffDict(),
+                isErrorPoint
+            );
+
+            Info<< nl << "    Number of bad faces     : " << nBadFaces << nl
+                << "    Number of marked points : "
+                << returnReduce(isErrorPoint.count(), sumOp<unsigned int>())
+                << endl;
+
+            updatePointErrorCount
+            (
+                isErrorPoint,
+                origToCurrentPointMap,
+                pointErrorCount
+            );
+
+            checkMeshEdgesAndRelaxEdges
+            (
+                newMesh,
+                origToCurrentPointMap,
+                isErrorPoint,
+                pointErrorCount
+            );
+
+            checkMeshFacesAndRelaxEdges
+            (
+                newMesh,
+                origToCurrentPointMap,
+                isErrorPoint,
+                pointErrorCount
+            );
+
+            newBadFaces = false;
+            forAll(mesh_.points(), pI)
+            {
+                if (isErrorPoint[origToCurrentPointMap[pI]])
+                {
+                    if (!newErrorPoint[pI])
+                    {
+                        newBadFaces = true;
+                        break;
+                    }
+                }
+            }
+
+            reduce(newBadFaces, orOp<bool>());
+        }
+        else
+        {
+            return -1;
+        }
+    }
+
+    return nBadFaces;
+}
+
+
+Foam::label Foam::polyMeshFilter::filterFaces
+(
+    polyMesh& newMesh,
+    scalarField& newMeshFaceFilterFactor,
+    labelList& origToCurrentPointMap
+)
+{
+    // Per edge collapse status
+    PackedBoolList collapseEdge(newMesh.nEdges());
+
+    Map<point> collapsePointToLocation(newMesh.nPoints());
+
+    edgeCollapser collapser(newMesh, collapseFacesCoeffDict());
+
+    {
+        // Collapse faces
+        labelPair nCollapsedPtEdge = collapser.markSmallSliverFaces
+        (
+            newMeshFaceFilterFactor,
+            pointPriority_(),
+            collapseEdge,
+            collapsePointToLocation
+        );
+
+        label nCollapsed = 0;
+        forAll(nCollapsedPtEdge, collapseTypeI)
+        {
+            nCollapsed += nCollapsedPtEdge[collapseTypeI];
+        }
+
+        reduce(nCollapsed, sumOp<label>());
+
+        label nToPoint = returnReduce(nCollapsedPtEdge.first(), sumOp<label>());
+        label nToEdge = returnReduce(nCollapsedPtEdge.second(), sumOp<label>());
+
+        Info<< indent
+            << "Collapsing " << nCollapsed << " faces "
+            << "(to point = " << nToPoint << ", to edge = " << nToEdge << ")"
+            << endl;
+
+        if (nCollapsed == 0)
+        {
+            return 0;
+        }
+    }
+
+    // Merge edge collapses into consistent collapse-network.
+    // Make sure no cells get collapsed.
+    List<pointEdgeCollapse> allPointInfo;
+    const globalIndex globalPoints(newMesh.nPoints());
+
+    collapser.consistentCollapse
+    (
+        globalPoints,
+        pointPriority_(),
+        collapsePointToLocation,
+        collapseEdge,
+        allPointInfo
+    );
+
+    label nLocalCollapse = collapseEdge.count();
+
+    reduce(nLocalCollapse, sumOp<label>());
+    Info<< nl << indent << "Collapsing " << nLocalCollapse
+        << " edges after synchronisation and PointEdgeWave" << endl;
+
+    if (nLocalCollapse == 0)
+    {
+        return 0;
+    }
+
+    {
+        // Apply collapses to current mesh
+        polyTopoChange newMeshMod(newMesh);
+
+        // Insert mesh refinement into polyTopoChange.
+        collapser.setRefinement(allPointInfo, newMeshMod);
+
+        Info<< indent << "Apply changes to the current mesh" << endl;
+
+        // Apply changes to current mesh
+        autoPtr<mapPolyMesh> newMapPtr = newMeshMod.changeMesh
+        (
+            newMesh,
+            false
+        );
+        const mapPolyMesh& newMap = newMapPtr();
+
+        // Update fields
+        newMesh.updateMesh(newMap);
+        if (newMap.hasMotionPoints())
+        {
+            newMesh.movePoints(newMap.preMotionPoints());
+        }
+
+        updatePointPriorities(newMesh, newMap.pointMap());
+
+        mapOldMeshFaceFieldToNewMesh
+        (
+            newMesh,
+            newMap.faceMap(),
+            newMeshFaceFilterFactor
+        );
+
+        updateOldToNewPointMap
+        (
+            newMap.reversePointMap(),
+            origToCurrentPointMap
+        );
+    }
+
+    return nLocalCollapse;
+}
+
+
+Foam::label Foam::polyMeshFilter::filterEdges
+(
+    polyMesh& newMesh,
+    scalarField& newMeshMinEdgeLen,
+    labelList& origToCurrentPointMap
+)
+{
+    // Per edge collapse status
+    PackedBoolList collapseEdge(newMesh.nEdges());
+
+    Map<point> collapsePointToLocation(newMesh.nPoints());
+
+    edgeCollapser collapser(newMesh, collapseFacesCoeffDict());
+
+    // Work out which edges to collapse
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // This is by looking at minEdgeLen (to avoid frozen edges)
+    // and marking in collapseEdge.
+    label nSmallCollapsed = collapser.markSmallEdges
+    (
+        newMeshMinEdgeLen,
+        pointPriority_(),
+        collapseEdge,
+        collapsePointToLocation
+    );
+
+    reduce(nSmallCollapsed, sumOp<label>());
+    Info<< indent << "Collapsing " << nSmallCollapsed
+        << " small edges" << endl;
+
+    // Merge inline edges
+    label nMerged = collapser.markMergeEdges
+    (
+        maxCos(),
+        pointPriority_(),
+        collapseEdge,
+        collapsePointToLocation
+    );
+
+    reduce(nMerged, sumOp<label>());
+    Info<< indent << "Collapsing " << nMerged << " in line edges"
+        << endl;
+
+    if (nMerged + nSmallCollapsed == 0)
+    {
+        return 0;
+    }
+
+    // Merge edge collapses into consistent collapse-network.
+    // Make sure no cells get collapsed.
+    List<pointEdgeCollapse> allPointInfo;
+    const globalIndex globalPoints(newMesh.nPoints());
+
+    collapser.consistentCollapse
+    (
+        globalPoints,
+        pointPriority_(),
+        collapsePointToLocation,
+        collapseEdge,
+        allPointInfo
+    );
+
+    label nLocalCollapse = collapseEdge.count();
+
+    reduce(nLocalCollapse, sumOp<label>());
+    Info<< nl << indent << "Collapsing " << nLocalCollapse
+        << " edges after synchronisation and PointEdgeWave" << endl;
+
+    if (nLocalCollapse == 0)
+    {
+        return 0;
+    }
+
+    // Apply collapses to current mesh
+    polyTopoChange newMeshMod(newMesh);
+
+    // Insert mesh refinement into polyTopoChange.
+    collapser.setRefinement(allPointInfo, newMeshMod);
+
+    Info<< indent << "Apply changes to the current mesh" << endl;
+
+    // Apply changes to current mesh
+    autoPtr<mapPolyMesh> newMapPtr = newMeshMod.changeMesh
+    (
+        newMesh,
+        false
+    );
+    const mapPolyMesh& newMap = newMapPtr();
+
+    // Update fields
+    newMesh.updateMesh(newMap);
+    if (newMap.hasMotionPoints())
+    {
+        newMesh.movePoints(newMap.preMotionPoints());
+    }
+
+    // Synchronise the factors
+    mapOldMeshEdgeFieldToNewMesh
+    (
+        newMesh,
+        newMap.pointMap(),
+        newMeshMinEdgeLen
+    );
+
+    updateOldToNewPointMap
+    (
+        newMap.reversePointMap(),
+        origToCurrentPointMap
+    );
+
+    updatePointPriorities(newMesh, newMap.pointMap());
+
+    return nLocalCollapse;
+}
+
+
 void Foam::polyMeshFilter::updatePointErrorCount
 (
     const PackedBoolList& isErrorPoint,
@@ -113,11 +531,11 @@ void Foam::polyMeshFilter::checkMeshEdgesAndRelaxEdges
 
         if
         (
-            pointErrorCount[e[0]] >= maxPointErrorCount_
-         || pointErrorCount[e[1]] >= maxPointErrorCount_
+            pointErrorCount[e[0]] >= maxPointErrorCount()
+         || pointErrorCount[e[1]] >= maxPointErrorCount()
         )
         {
-            minEdgeLen_[edgeI] = 0;
+            minEdgeLen_[edgeI] = -1;
         }
 
         if
@@ -126,13 +544,13 @@ void Foam::polyMeshFilter::checkMeshEdgesAndRelaxEdges
          || (newEnd >= 0 && isErrorPoint[newEnd])
         )
         {
-            minEdgeLen_[edgeI] *= edgeReductionFactor_;
+            minEdgeLen_[edgeI] *= edgeReductionFactor();
         }
     }
 
     syncTools::syncEdgeList(mesh_, minEdgeLen_, minEqOp<scalar>(), scalar(0.0));
 
-    for (label smoothIter = 0; smoothIter < maxSmoothIters_; ++smoothIter)
+    for (label smoothIter = 0; smoothIter < maxSmoothIters(); ++smoothIter)
     {
         // Smooth minEdgeLen
         forAll(mesh_.edges(), edgeI)
@@ -190,14 +608,14 @@ void Foam::polyMeshFilter::checkMeshFacesAndRelaxEdges
         {
             const label ptIndex = oldToNewMesh[f[fpI]];
 
-            if (pointErrorCount[f[fpI]] >= maxPointErrorCount_)
+            if (pointErrorCount[f[fpI]] >= maxPointErrorCount())
             {
-                faceFilterFactor_[faceI] = 0;
+                faceFilterFactor_[faceI] = -1;
             }
 
             if (isErrorPoint[ptIndex])
             {
-                faceFilterFactor_[faceI] *= faceReductionFactor_;
+                faceFilterFactor_[faceI] *= faceReductionFactor();
 
                 break;
             }
@@ -206,7 +624,7 @@ void Foam::polyMeshFilter::checkMeshFacesAndRelaxEdges
 
     syncTools::syncFaceList(mesh_, faceFilterFactor_, minEqOp<scalar>());
 
-    for (label smoothIter = 0; smoothIter < maxSmoothIters_; ++smoothIter)
+    for (label smoothIter = 0; smoothIter < maxSmoothIters(); ++smoothIter)
     {
         // Smooth faceFilterFactor
         forAll(faces, faceI)
@@ -267,45 +685,32 @@ void Foam::polyMeshFilter::checkMeshFacesAndRelaxEdges
 }
 
 
-Foam::labelList Foam::polyMeshFilter::findBoundaryPoints
+void Foam::polyMeshFilter::updatePointPriorities
 (
-    const polyMesh& mesh//,
-//    labelIOList& boundaryIOPts
-) const
+    const polyMesh& newMesh,
+    const labelList& pointMap
+)
 {
-    const polyBoundaryMesh& bMesh = mesh.boundaryMesh();
+    labelList newPointPriority(newMesh.nPoints(), labelMin);
+    const labelList& currPointPriority = pointPriority_();
 
-    labelList boundaryPoint(mesh.nPoints(), -1);
-
-    // Get all processor boundary points and the processor patch label
-    // that they are on.
-    forAll(bMesh, patchI)
+    forAll(newPointPriority, ptI)
     {
-        const polyPatch& patch = bMesh[patchI];
+        const label newPointToOldPoint = pointMap[ptI];
+        const label origPointPriority = currPointPriority[newPointToOldPoint];
 
-        if (!isA<coupledPolyPatch>(patch))
-        {
-            forAll(patch, fI)
-            {
-                const face& f = patch[fI];
-
-                forAll(f, fp)
-                {
-                    boundaryPoint[f[fp]] = 0; //boundaryIOPts[f[fp]];
-                }
-            }
-        }
+        newPointPriority[ptI] = max(origPointPriority, newPointPriority[ptI]);
     }
 
     syncTools::syncPointList
     (
-        mesh,
-        boundaryPoint,
+        newMesh,
+        newPointPriority,
         maxEqOp<label>(),
         labelMin
     );
 
-    return boundaryPoint;
+    pointPriority_.reset(new labelList(newPointPriority));
 }
 
 
@@ -315,12 +720,45 @@ void Foam::polyMeshFilter::printScalarFieldStats
     const scalarField& fld
 ) const
 {
+    scalar sum = 0;
+    scalar validElements = 0;
+    scalar min = GREAT;
+    scalar max = -GREAT;
+
+    forAll(fld, i)
+    {
+        const scalar fldElement = fld[i];
+
+        if (fldElement >= 0)
+        {
+            sum += fldElement;
+
+            if (fldElement < min)
+            {
+                min = fldElement;
+            }
+
+            if (fldElement > max)
+            {
+                max = fldElement;
+            }
+
+            validElements++;
+        }
+    }
+
+    reduce(sum, sumOp<scalar>());
+    reduce(min, minOp<scalar>());
+    reduce(max, maxOp<scalar>());
+    reduce(validElements, sumOp<label>());
+    const label totFieldSize = returnReduce(fld.size(), sumOp<label>());
+
     Info<< incrIndent << indent << desc
-        << ": min = " << returnReduce(min(fld), minOp<scalar>())
-        << " av = "
-        << returnReduce(sum(fld), sumOp<scalar>())
-          /returnReduce(fld.size(), sumOp<label>())
-        << " max = " << returnReduce(max(fld), maxOp<scalar>())
+        << ": min = " << min
+        << " av = " << sum/(validElements + SMALL)
+        << " max = " << max << nl
+        << indent
+        << "    " << validElements << " / " << totFieldSize << " elements used"
         << decrIndent << endl;
 }
 
@@ -415,119 +853,59 @@ void Foam::polyMeshFilter::updateOldToNewPointMap
 
 Foam::polyMeshFilter::polyMeshFilter(const fvMesh& mesh)
 :
-    mesh_(mesh),
-    newMeshPtr_(),
-    dict_
+    polyMeshFilterSettings
     (
-        IOobject
+        IOdictionary
         (
-            "collapseDict",
-            mesh.time().system(),
-            mesh.time(),
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    ),
-    controlMeshQuality_
-    (
-        dict_.lookupOrDefault<Switch>("controlMeshQuality", false)
-    ),
-    collapseEdgesCoeffDict_(dict_.subDict("collapseEdgesCoeffs")),
-    collapseFacesCoeffDict_(dict_.subOrEmptyDict("collapseFacesCoeffs")),
-    meshQualityCoeffDict_(dict_.subOrEmptyDict("controlMeshQualityCoeffs")),
-    minLen_(readScalar(collapseEdgesCoeffDict_.lookup("minimumEdgeLength"))),
-    maxCos_
-    (
-        ::cos
-        (
-            degToRad
+            IOobject
             (
-                readScalar(collapseEdgesCoeffDict_.lookup("maximumMergeAngle"))
+                "collapseDict",
+                mesh.time().system(),
+                mesh.time(),
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE
             )
         )
     ),
-    edgeReductionFactor_
-    (
-        meshQualityCoeffDict_.lookupOrDefault<scalar>("edgeReductionFactor", -1)
-    ),
-    maxIterations_
-    (
-        meshQualityCoeffDict_.lookupOrAddDefault<label>("maximumIterations", 1)
-    ),
-    maxSmoothIters_
-    (
-        meshQualityCoeffDict_.lookupOrAddDefault<label>
-        (
-            "maximumSmoothingIterations",
-            0
-        )
-    ),
-    initialFaceLengthFactor_
-    (
-        collapseFacesCoeffDict_.lookupOrAddDefault<scalar>
-        (
-            "initialFaceLengthFactor",
-            -1
-        )
-    ),
-    faceReductionFactor_
-    (
-        meshQualityCoeffDict_.lookupOrAddDefault<scalar>
-        (
-            "faceReductionFactor",
-            -1
-        )
-    ),
-    maxPointErrorCount_
-    (
-        meshQualityCoeffDict_.lookupOrAddDefault<label>("maxPointErrorCount", 0)
-    ),
+    mesh_(mesh),
+    newMeshPtr_(),
+    originalPointPriority_(mesh.nPoints(), labelMin),
+    pointPriority_(),
     minEdgeLen_(),
     faceFilterFactor_()
 {
-    Info<< "Merging:" << nl
-        << "    edges with length less than " << minLen_ << " metres" << nl
-        << "    edges split by a point with edges in line to within "
-        << radToDeg(::acos(maxCos_)) << " degrees" << nl
-        << "    Minimum edge length reduction factor = "
-        << edgeReductionFactor_ << nl
-        << endl;
+    writeSettings(Info);
+}
 
-    if (collapseFacesCoeffDict_.empty())
-    {
-        Info<< "Face collapsing is off" << endl;
-    }
-    else
-    {
-        Info<< "Face collapsing is on" << endl;
-        Info<< "    Initial face length factor = "<< initialFaceLengthFactor_
-            << endl;
-    }
 
-    Info<< "Control mesh quality = " << controlMeshQuality_.asText() << endl;
-
-    if (controlMeshQuality_)
-    {
-        Info<< "    Minimum edge length reduction factor = "
-            << edgeReductionFactor_ << nl
-            << "    Minimum face area reduction factor = "
-            << faceReductionFactor_ << endl;
-
-        Info<< "    Maximum number of collapse iterations = " << maxIterations_
-            << endl;
-
-        Info<< "    Maximum number of edge/face reduction factor smoothing "
-            << "iterations = " << maxSmoothIters_ << endl;
-
-        Info<< "    Maximum number of times a point can contribute to bad "
-            << "faces across " << nl
-            << "    collapse iterations = " << maxPointErrorCount_
-            << endl;
-    }
-
-    Info<< "Selectively disabling wanted collapses until resulting quality"
-        << " satisfies constraints in system/meshQualityDict" << nl
-        << endl;
+Foam::polyMeshFilter::polyMeshFilter
+(
+    const fvMesh& mesh,
+    const labelList& pointPriority
+)
+:
+    polyMeshFilterSettings
+    (
+        IOdictionary
+        (
+            IOobject
+            (
+                "collapseDict",
+                mesh.time().system(),
+                mesh.time(),
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE
+            )
+        )
+    ),
+    mesh_(mesh),
+    newMeshPtr_(),
+    originalPointPriority_(pointPriority),
+    pointPriority_(),
+    minEdgeLen_(),
+    faceFilterFactor_()
+{
+    writeSettings(Info);
 }
 
 
@@ -541,401 +919,27 @@ Foam::polyMeshFilter::~polyMeshFilter()
 
 Foam::label Foam::polyMeshFilter::filter(const label nOriginalBadFaces)
 {
-    label nBadFaces = labelMax;
-    label nOuterIterations = 0;
+    minEdgeLen_.resize(mesh_.nEdges(), minLen());
+    faceFilterFactor_.resize(mesh_.nFaces(), initialFaceLengthFactor());
 
-    minEdgeLen_.resize(mesh_.nEdges(), minLen_);
-    faceFilterFactor_.resize(mesh_.nFaces(), initialFaceLengthFactor_);
+    return filterFacesLoop(nOriginalBadFaces);
+}
 
-    // Maintain the number of times a point has been part of a bad face
-    labelList pointErrorCount(mesh_.nPoints(), 0);
 
-    PackedBoolList newErrorPoint(mesh_.nPoints());
-    edgeCollapser::checkMeshQuality
-    (
-        mesh_,
-        meshQualityCoeffDict_,
-        newErrorPoint
-    );
+Foam::label Foam::polyMeshFilter::filter(const faceZone& fZone)
+{
+    minEdgeLen_.resize(mesh_.nEdges(), minLen());
+    faceFilterFactor_.resize(mesh_.nFaces(), initialFaceLengthFactor());
 
-    bool newBadFaces = true;
-
-    // Main loop
-    // ~~~~~~~~~
-    // It tries and do some collapses, checks the resulting mesh and
-    // 'freezes' some edges (by marking in minEdgeLen) and tries again.
-    // This will iterate ultimately to the situation where every edge is
-    // frozen and nothing gets collapsed.
-    while
-    (
-        nOuterIterations < maxIterations_
-     //&& nBadFaces > nOriginalBadFaces
-     && newBadFaces
-    )
+    forAll(faceFilterFactor_, fI)
     {
-        Info<< nl << "Outer Iteration = " << nOuterIterations++ << nl
-            << endl;
-
-        printScalarFieldStats("Edge Filter Factor", minEdgeLen_);
-        printScalarFieldStats("Face Filter Factor", faceFilterFactor_);
-
-        // Reset the new mesh to the old mesh
-        newMeshPtr_ = copyMesh(mesh_);
-        fvMesh& newMesh = newMeshPtr_();
-
-        scalarField newMeshFaceFilterFactor = faceFilterFactor_;
-
-        labelList origToCurrentPointMap(identity(newMesh.nPoints()));
+        if (fZone.whichFace(fI) == -1)
         {
-
-            label nInnerIterations = 0;
-            label nPrevLocalCollapse = labelMax;
-
-            Info<< incrIndent;
-
-            while (true)
-            {
-                Info<< nl << indent << "Inner iteration = "
-                    << nInnerIterations++ << nl << incrIndent << endl;
-
-                // Per edge collapse status
-                PackedBoolList collapseEdge(newMesh.nEdges());
-
-                Map<point> collapsePointToLocation(newMesh.nPoints());
-
-                // Mark points on boundary
-                const labelList boundaryPoint = findBoundaryPoints(newMesh);
-
-                edgeCollapser collapser(newMesh, collapseFacesCoeffDict_);
-
-                {
-                    // Collapse faces
-                    labelPair nCollapsedPtEdge = collapser.markSmallSliverFaces
-                    (
-                        newMeshFaceFilterFactor,
-                        boundaryPoint,
-                        collapseEdge,
-                        collapsePointToLocation
-                    );
-
-                    label nCollapsed = 0;
-                    forAll(nCollapsedPtEdge, collapseTypeI)
-                    {
-                        nCollapsed += nCollapsedPtEdge[collapseTypeI];
-                    }
-
-                    reduce(nCollapsed, sumOp<label>());
-
-                    Info<< indent
-                        << "Collapsing " << nCollapsed << " faces"
-                        << " (to point = "
-                        << returnReduce
-                           (
-                               nCollapsedPtEdge.first(),
-                               sumOp<label>()
-                           )
-                        << ", to edge = "
-                        << returnReduce
-                           (
-                               nCollapsedPtEdge.second(),
-                               sumOp<label>()
-                           )
-                        << ")" << endl;
-
-                    if (nCollapsed == 0)
-                    {
-                        Info<< decrIndent;
-                        Info<< decrIndent;
-                        break;
-                    }
-                }
-
-                // Merge edge collapses into consistent collapse-network.
-                // Make sure no cells get collapsed.
-                List<pointEdgeCollapse> allPointInfo;
-                const globalIndex globalPoints(newMesh.nPoints());
-
-                collapser.consistentCollapse
-                (
-                    globalPoints,
-                    boundaryPoint,
-                    collapsePointToLocation,
-                    collapseEdge,
-                    allPointInfo
-                );
-
-                label nLocalCollapse = collapseEdge.count();
-
-                reduce(nLocalCollapse, sumOp<label>());
-                Info<< nl << indent << "Collapsing " << nLocalCollapse
-                    << " edges after synchronisation and PointEdgeWave" << endl;
-
-                if (nLocalCollapse >= nPrevLocalCollapse)
-                {
-                    Info<< decrIndent;
-                    Info<< decrIndent;
-                    break;
-                }
-                else
-                {
-                    nPrevLocalCollapse = nLocalCollapse;
-                }
-
-                {
-                    // Apply collapses to current mesh
-                    polyTopoChange newMeshMod(newMesh);
-
-                    // Insert mesh refinement into polyTopoChange.
-                    collapser.setRefinement(allPointInfo, newMeshMod);
-
-                    Info<< indent << "Apply changes to the current mesh"
-                        << decrIndent << endl;
-
-                    // Apply changes to current mesh
-                    autoPtr<mapPolyMesh> newMapPtr = newMeshMod.changeMesh
-                    (
-                        newMesh,
-                        false
-                    );
-                    const mapPolyMesh& newMap = newMapPtr();
-
-                    // Update fields
-                    newMesh.updateMesh(newMap);
-                    if (newMap.hasMotionPoints())
-                    {
-                        newMesh.movePoints(newMap.preMotionPoints());
-                    }
-
-                    // Relabel the boundary points
-        //            labelList newBoundaryPoints(newMesh.nPoints(), -1);
-        //            forAll(newBoundaryPoints, pI)
-        //            {
-        //                const label newToOldptI= map.pointMap()[pI];
-        //                newBoundaryPoints[pI] = boundaryIOPts[newToOldptI];
-        //            }
-        //            boundaryIOPts = newBoundaryPoints;
-
-
-                    mapOldMeshFaceFieldToNewMesh
-                    (
-                        newMesh,
-                        newMap.faceMap(),
-                        newMeshFaceFilterFactor
-                    );
-
-                    updateOldToNewPointMap
-                    (
-                        newMap.reversePointMap(),
-                        origToCurrentPointMap
-                    );
-                }
-            }
-        }
-
-
-        scalarField newMeshMinEdgeLen = minEdgeLen_;
-
-        label nInnerIterations = 0;
-        label nPrevLocalCollapse = labelMax;
-
-        while (true)
-        {
-            Info<< nl << indent << "Inner iteration = "
-                << nInnerIterations++ << nl << incrIndent << endl;
-
-            // Per edge collapse status
-            PackedBoolList collapseEdge(newMesh.nEdges());
-
-            Map<point> collapsePointToLocation(newMesh.nPoints());
-
-            // Mark points on boundary
-            const labelList boundaryPoint = findBoundaryPoints
-            (
-                newMesh//,
-//                    boundaryIOPts
-            );
-
-            edgeCollapser collapser(newMesh, collapseFacesCoeffDict_);
-
-            // Work out which edges to collapse
-            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-            // This is by looking at minEdgeLen (to avoid frozen edges)
-            // and marking in collapseEdge.
-            label nSmallCollapsed = collapser.markSmallEdges
-            (
-                newMeshMinEdgeLen,
-                boundaryPoint,
-                collapseEdge,
-                collapsePointToLocation
-            );
-
-            reduce(nSmallCollapsed, sumOp<label>());
-            Info<< indent << "Collapsing " << nSmallCollapsed
-                << " small edges" << endl;
-
-            // Merge inline edges
-            label nMerged = collapser.markMergeEdges
-            (
-                maxCos_,
-                boundaryPoint,
-                collapseEdge,
-                collapsePointToLocation
-            );
-
-            reduce(nMerged, sumOp<label>());
-            Info<< indent << "Collapsing " << nMerged << " in line edges"
-                << endl;
-
-            if (nMerged + nSmallCollapsed == 0)
-            {
-                Info<< decrIndent;
-                break;
-            }
-
-            // Merge edge collapses into consistent collapse-network.
-            // Make sure no cells get collapsed.
-            List<pointEdgeCollapse> allPointInfo;
-            const globalIndex globalPoints(newMesh.nPoints());
-
-            collapser.consistentCollapse
-            (
-                globalPoints,
-                boundaryPoint,
-                collapsePointToLocation,
-                collapseEdge,
-                allPointInfo
-            );
-
-            label nLocalCollapse = collapseEdge.count();
-
-            reduce(nLocalCollapse, sumOp<label>());
-            Info<< nl << indent << "Collapsing " << nLocalCollapse
-                << " edges after synchronisation and PointEdgeWave" << endl;
-
-            if (nLocalCollapse >= nPrevLocalCollapse)
-            {
-                Info<< decrIndent;
-                break;
-            }
-            else
-            {
-                nPrevLocalCollapse = nLocalCollapse;
-            }
-
-            // Apply collapses to current mesh
-            polyTopoChange newMeshMod(newMesh);
-
-            // Insert mesh refinement into polyTopoChange.
-            collapser.setRefinement(allPointInfo, newMeshMod);
-
-            Info<< indent << "Apply changes to the current mesh"
-                << decrIndent << endl;
-
-            // Apply changes to current mesh
-            autoPtr<mapPolyMesh> newMapPtr = newMeshMod.changeMesh
-            (
-                newMesh,
-                false
-            );
-            const mapPolyMesh& newMap = newMapPtr();
-
-            // Update fields
-            newMesh.updateMesh(newMap);
-            if (newMap.hasMotionPoints())
-            {
-                newMesh.movePoints(newMap.preMotionPoints());
-            }
-
-            // Relabel the boundary points
-//            labelList newBoundaryPoints(newMesh.nPoints(), -1);
-//            forAll(newBoundaryPoints, pI)
-//            {
-//                const label newToOldptI= map.pointMap()[pI];
-//                newBoundaryPoints[pI] = boundaryIOPts[newToOldptI];
-//            }
-//            boundaryIOPts = newBoundaryPoints;
-
-            // Synchronise the factors
-            mapOldMeshEdgeFieldToNewMesh
-            (
-                newMesh,
-                newMap.pointMap(),
-                newMeshMinEdgeLen
-            );
-
-            updateOldToNewPointMap
-            (
-                newMap.reversePointMap(),
-                origToCurrentPointMap
-            );
-        }
-
-
-        // Mesh check
-        // ~~~~~~~~~~~~~~~~~~
-        // Do not allow collapses in regions of error.
-        // Updates minEdgeLen, nRelaxedEdges
-
-        if (controlMeshQuality_)
-        {
-            PackedBoolList isErrorPoint(newMesh.nPoints());
-            nBadFaces = edgeCollapser::checkMeshQuality
-            (
-                newMesh,
-                meshQualityCoeffDict_,
-                isErrorPoint
-            );
-
-            Info<< nl << "    Number of bad faces     : " << nBadFaces << nl
-                << "    Number of marked points : "
-                << returnReduce(isErrorPoint.count(), sumOp<unsigned int>())
-                << endl;
-
-            updatePointErrorCount
-            (
-                isErrorPoint,
-                origToCurrentPointMap,
-                pointErrorCount
-            );
-
-            checkMeshEdgesAndRelaxEdges
-            (
-                newMesh,
-                origToCurrentPointMap,
-                isErrorPoint,
-                pointErrorCount
-            );
-
-            checkMeshFacesAndRelaxEdges
-            (
-                newMesh,
-                origToCurrentPointMap,
-                isErrorPoint,
-                pointErrorCount
-            );
-
-            newBadFaces = false;
-            forAll(mesh_.points(), pI)
-            {
-                if (isErrorPoint[origToCurrentPointMap[pI]])
-                {
-                    if (!newErrorPoint[pI])
-                    {
-                        newBadFaces = true;
-                        break;
-                    }
-                }
-            }
-
-            reduce(newBadFaces, orOp<bool>());
-        }
-        else
-        {
-            return -1;
+            faceFilterFactor_[fI] = -1;
         }
     }
 
-    return nBadFaces;
+    return filterFacesLoop(0);
 }
 
 
@@ -948,7 +952,7 @@ Foam::label Foam::polyMeshFilter::filterEdges
     label nPreviousBadFaces = labelMax;
     label nOuterIterations = 0;
 
-    minEdgeLen_.resize(mesh_.nEdges(), minLen_);
+    minEdgeLen_.resize(mesh_.nEdges(), minLen());
     faceFilterFactor_.resize(0);
 
     labelList pointErrorCount(mesh_.nPoints(), 0);
@@ -961,7 +965,7 @@ Foam::label Foam::polyMeshFilter::filterEdges
     // frozen and nothing gets collapsed.
     while
     (
-        nOuterIterations < maxIterations_
+        nOuterIterations < maxIterations()
      && nBadFaces > nOriginalBadFaces
      && nBadFaces < nPreviousBadFaces
     )
@@ -978,455 +982,35 @@ Foam::label Foam::polyMeshFilter::filterEdges
         fvMesh& newMesh = newMeshPtr_();
 
         scalarField newMeshMinEdgeLen = minEdgeLen_;
+        pointPriority_.reset(new labelList(originalPointPriority_));
 
         labelList origToCurrentPointMap(identity(newMesh.nPoints()));
 
-        label nInnerIterations = 0;
-        label nPrevLocalCollapse = labelMax;
-
         Info<< incrIndent;
 
-        while (true)
-        {
-            Info<< nl << indent << "Inner iteration = "
-                << nInnerIterations++ << nl << incrIndent << endl;
-
-            // Per edge collapse status
-            PackedBoolList collapseEdge(newMesh.nEdges());
-
-            Map<point> collapsePointToLocation(newMesh.nPoints());
-
-            // Mark points on boundary
-            const labelList boundaryPoint = findBoundaryPoints(newMesh);
-
-            edgeCollapser collapser(newMesh, collapseFacesCoeffDict_);
-
-            // Work out which edges to collapse
-            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-            // This is by looking at minEdgeLen (to avoid frozen edges)
-            // and marking in collapseEdge.
-            label nSmallCollapsed = collapser.markSmallEdges
-            (
-                newMeshMinEdgeLen,
-                boundaryPoint,
-                collapseEdge,
-                collapsePointToLocation
-            );
-
-            reduce(nSmallCollapsed, sumOp<label>());
-            Info<< indent << "Collapsing " << nSmallCollapsed
-                << " small edges" << endl;
-
-            // Merge inline edges
-            label nMerged = collapser.markMergeEdges
-            (
-                maxCos_,
-                boundaryPoint,
-                collapseEdge,
-                collapsePointToLocation
-            );
-
-            reduce(nMerged, sumOp<label>());
-            Info<< indent << "Collapsing " << nMerged << " in line edges"
-                << endl;
-
-            if (nMerged + nSmallCollapsed == 0)
-            {
-                Info<< decrIndent;
-                Info<< decrIndent;
-                break;
-            }
-
-            // Merge edge collapses into consistent collapse-network.
-            // Make sure no cells get collapsed.
-            List<pointEdgeCollapse> allPointInfo;
-            const globalIndex globalPoints(newMesh.nPoints());
-
-            collapser.consistentCollapse
-            (
-                globalPoints,
-                boundaryPoint,
-                collapsePointToLocation,
-                collapseEdge,
-                allPointInfo
-            );
-
-            label nLocalCollapse = collapseEdge.count();
-
-            reduce(nLocalCollapse, sumOp<label>());
-            Info<< nl << indent << "Collapsing " << nLocalCollapse
-                << " edges after synchronisation and PointEdgeWave" << endl;
-
-            if (nLocalCollapse >= nPrevLocalCollapse)
-            {
-                Info<< decrIndent;
-                Info<< decrIndent;
-                break;
-            }
-            else
-            {
-                nPrevLocalCollapse = nLocalCollapse;
-            }
-
-            // Apply collapses to current mesh
-            polyTopoChange newMeshMod(newMesh);
-
-            // Insert mesh refinement into polyTopoChange.
-            collapser.setRefinement(allPointInfo, newMeshMod);
-
-            Info<< indent << "Apply changes to the current mesh"
-                << decrIndent << endl;
-
-            // Apply changes to current mesh
-            autoPtr<mapPolyMesh> newMapPtr = newMeshMod.changeMesh
-            (
-                newMesh,
-                false
-            );
-            const mapPolyMesh& newMap = newMapPtr();
-
-            // Update fields
-            newMesh.updateMesh(newMap);
-            if (newMap.hasMotionPoints())
-            {
-                newMesh.movePoints(newMap.preMotionPoints());
-            }
-
-            // Relabel the boundary points
-//            labelList newBoundaryPoints(newMesh.nPoints(), -1);
-//            forAll(newBoundaryPoints, pI)
-//            {
-//                const label newToOldptI= map.pointMap()[pI];
-//                newBoundaryPoints[pI] = boundaryIOPts[newToOldptI];
-//            }
-//            boundaryIOPts = newBoundaryPoints;
-
-            // Synchronise the factors
-            mapOldMeshEdgeFieldToNewMesh
-            (
-                newMesh,
-                newMap.pointMap(),
-                newMeshMinEdgeLen
-            );
-
-            updateOldToNewPointMap
-            (
-                newMap.reversePointMap(),
-                origToCurrentPointMap
-            );
-        }
-
-        // Mesh check
-        // ~~~~~~~~~~~~~~~~~~
-        // Do not allow collapses in regions of error.
-        // Updates minEdgeLen, nRelaxedEdges
-
-        if (controlMeshQuality_)
-        {
-            PackedBoolList isErrorPoint(newMesh.nPoints());
-            nBadFaces = edgeCollapser::checkMeshQuality
-            (
-                newMesh,
-                meshQualityCoeffDict_,
-                isErrorPoint
-            );
-
-            Info<< nl << "    Number of bad faces     : " << nBadFaces << nl
-                << "    Number of marked points : "
-                << returnReduce(isErrorPoint.count(), sumOp<unsigned int>())
-                << endl;
-
-            updatePointErrorCount
-            (
-                isErrorPoint,
-                origToCurrentPointMap,
-                pointErrorCount
-            );
-
-            checkMeshEdgesAndRelaxEdges
-            (
-                newMesh,
-                origToCurrentPointMap,
-                isErrorPoint,
-                pointErrorCount
-            );
-        }
-        else
-        {
-            return -1;
-        }
-    }
-
-    return nBadFaces;
-}
-
-
-Foam::label Foam::polyMeshFilter::filterFaceZone(const faceZone& fZone)
-{
-    const label nOriginalBadFaces = 0;
-
-    label nBadFaces = labelMax;
-    label nOuterIterations = 0;
-
-    minEdgeLen_.resize(mesh_.nEdges(), minLen_);
-    faceFilterFactor_.resize(mesh_.nFaces(), initialFaceLengthFactor_);
-
-    // Maintain the number of times a point has been part of a bad face
-    labelList pointErrorCount(mesh_.nPoints(), 0);
-
-    // Main loop
-    // ~~~~~~~~~
-    // It tries and do some collapses, checks the resulting mesh and
-    // 'freezes' some edges (by marking in minEdgeLen) and tries again.
-    // This will iterate ultimately to the situation where every edge is
-    // frozen and nothing gets collapsed.
-    while
-    (
-        nOuterIterations < maxIterations_
-     && nBadFaces > nOriginalBadFaces
-    )
-    {
-        Info<< nl << "Outer Iteration = " << nOuterIterations++ << nl
-            << endl;
-
-        printScalarFieldStats("Edge Filter Factor", minEdgeLen_);
-        printScalarFieldStats("Face Filter Factor", faceFilterFactor_);
-
-        // Reset the new mesh to the old mesh
-        newMeshPtr_ = copyMesh(mesh_);
-        fvMesh& newMesh = newMeshPtr_();
-
-        scalarField newMeshFaceFilterFactor = faceFilterFactor_;
-
-        labelList origToCurrentPointMap(identity(newMesh.nPoints()));
-        {
-
-            label nInnerIterations = 0;
-            label nPrevLocalCollapse = labelMax;
-
-            Info<< incrIndent;
-
-            while (true)
-            {
-                Info<< nl << indent << "Inner iteration = "
-                    << nInnerIterations++ << nl << incrIndent << endl;
-
-                // Per edge collapse status
-                PackedBoolList collapseEdge(newMesh.nEdges());
-
-                Map<point> collapsePointToLocation(newMesh.nPoints());
-
-                // Mark points on boundary
-                const labelList boundaryPoint = findBoundaryPoints(newMesh);
-
-                edgeCollapser collapser(newMesh, collapseFacesCoeffDict_);
-
-                {
-                    // Collapse faces
-                    labelPair nCollapsedPtEdge = collapser.markFaceZoneEdges
-                    (
-                        fZone,
-                        newMeshFaceFilterFactor,
-                        boundaryPoint,
-                        collapseEdge,
-                        collapsePointToLocation
-                    );
-
-                    label nCollapsed = 0;
-                    forAll(nCollapsedPtEdge, collapseTypeI)
-                    {
-                        nCollapsed += nCollapsedPtEdge[collapseTypeI];
-                    }
-
-                    reduce(nCollapsed, sumOp<label>());
-
-                    Info<< indent
-                        << "Collapsing " << nCollapsed << " faces"
-                        << " (to point = "
-                        << returnReduce
-                           (
-                               nCollapsedPtEdge.first(),
-                               sumOp<label>()
-                           )
-                        << ", to edge = "
-                        << returnReduce
-                           (
-                               nCollapsedPtEdge.second(),
-                               sumOp<label>()
-                           )
-                        << ")" << endl;
-
-                    if (nCollapsed == 0)
-                    {
-                        Info<< decrIndent;
-                        Info<< decrIndent;
-                        break;
-                    }
-                }
-
-                // Merge edge collapses into consistent collapse-network.
-                // Make sure no cells get collapsed.
-                List<pointEdgeCollapse> allPointInfo;
-                const globalIndex globalPoints(newMesh.nPoints());
-
-                collapser.consistentCollapse
-                (
-                    globalPoints,
-                    boundaryPoint,
-                    collapsePointToLocation,
-                    collapseEdge,
-                    allPointInfo
-                );
-
-                label nLocalCollapse = collapseEdge.count();
-
-                reduce(nLocalCollapse, sumOp<label>());
-                Info<< nl << indent << "Collapsing " << nLocalCollapse
-                    << " edges after synchronisation and PointEdgeWave" << endl;
-
-                if (nLocalCollapse >= nPrevLocalCollapse)
-                {
-                    Info<< decrIndent;
-                    Info<< decrIndent;
-                    break;
-                }
-                else
-                {
-                    nPrevLocalCollapse = nLocalCollapse;
-                }
-
-                {
-                    // Apply collapses to current mesh
-                    polyTopoChange newMeshMod(newMesh);
-
-                    // Insert mesh refinement into polyTopoChange.
-                    collapser.setRefinement(allPointInfo, newMeshMod);
-
-                    Info<< indent << "Apply changes to the current mesh"
-                        << decrIndent << endl;
-
-                    // Apply changes to current mesh
-                    autoPtr<mapPolyMesh> newMapPtr = newMeshMod.changeMesh
-                    (
-                        newMesh,
-                        false
-                    );
-                    const mapPolyMesh& newMap = newMapPtr();
-
-                    // Update fields
-                    newMesh.updateMesh(newMap);
-                    if (newMap.hasMotionPoints())
-                    {
-                        newMesh.movePoints(newMap.preMotionPoints());
-                    }
-
-                    // Relabel the boundary points
-        //            labelList newBoundaryPoints(newMesh.nPoints(), -1);
-        //            forAll(newBoundaryPoints, pI)
-        //            {
-        //                const label newToOldptI= map.pointMap()[pI];
-        //                newBoundaryPoints[pI] = boundaryIOPts[newToOldptI];
-        //            }
-        //            boundaryIOPts = newBoundaryPoints;
-
-
-                    mapOldMeshFaceFieldToNewMesh
-                    (
-                        newMesh,
-                        newMap.faceMap(),
-                        newMeshFaceFilterFactor
-                    );
-
-                    updateOldToNewPointMap
-                    (
-                        newMap.reversePointMap(),
-                        origToCurrentPointMap
-                    );
-                }
-            }
-        }
-
-
-        scalarField newMeshMinEdgeLen = minEdgeLen_;
-
         label nInnerIterations = 0;
         label nPrevLocalCollapse = labelMax;
 
         while (true)
         {
-            Info<< nl << indent << "Inner iteration = "
-                << nInnerIterations++ << nl << incrIndent << endl;
+            Info<< nl
+                << indent << "Inner iteration = " << nInnerIterations++ << nl
+                << incrIndent << endl;
 
-            // Per edge collapse status
-            PackedBoolList collapseEdge(newMesh.nEdges());
-
-            Map<point> collapsePointToLocation(newMesh.nPoints());
-
-            // Mark points on boundary
-            const labelList boundaryPoint = findBoundaryPoints
-            (
-                newMesh//,
-//                    boundaryIOPts
-            );
-
-            edgeCollapser collapser(newMesh, collapseFacesCoeffDict_);
-
-            // Work out which edges to collapse
-            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-            // This is by looking at minEdgeLen (to avoid frozen edges)
-            // and marking in collapseEdge.
-            label nSmallCollapsed = collapser.markSmallEdges
+            label nLocalCollapse = filterEdges
             (
+                newMesh,
                 newMeshMinEdgeLen,
-                boundaryPoint,
-                collapseEdge,
-                collapsePointToLocation
+                origToCurrentPointMap
             );
 
-            reduce(nSmallCollapsed, sumOp<label>());
-            Info<< indent << "Collapsing " << nSmallCollapsed
-                << " small edges" << endl;
+            Info<< decrIndent;
 
-            // Merge inline edges
-            label nMerged = collapser.markMergeEdges
+            if
             (
-                maxCos_,
-                boundaryPoint,
-                collapseEdge,
-                collapsePointToLocation
-            );
-
-            reduce(nMerged, sumOp<label>());
-            Info<< indent << "Collapsing " << nMerged << " in line edges"
-                << endl;
-
-            if (nMerged + nSmallCollapsed == 0)
-            {
-                Info<< decrIndent;
-                break;
-            }
-
-            // Merge edge collapses into consistent collapse-network.
-            // Make sure no cells get collapsed.
-            List<pointEdgeCollapse> allPointInfo;
-            const globalIndex globalPoints(newMesh.nPoints());
-
-            collapser.consistentCollapse
-            (
-                globalPoints,
-                boundaryPoint,
-                collapsePointToLocation,
-                collapseEdge,
-                allPointInfo
-            );
-
-            label nLocalCollapse = collapseEdge.count();
-
-            reduce(nLocalCollapse, sumOp<label>());
-            Info<< nl << indent << "Collapsing " << nLocalCollapse
-                << " edges after synchronisation and PointEdgeWave" << endl;
-
-            if (nLocalCollapse >= nPrevLocalCollapse)
+                nLocalCollapse >= nPrevLocalCollapse
+             || nLocalCollapse == 0
+            )
             {
                 Info<< decrIndent;
                 break;
@@ -1435,68 +1019,20 @@ Foam::label Foam::polyMeshFilter::filterFaceZone(const faceZone& fZone)
             {
                 nPrevLocalCollapse = nLocalCollapse;
             }
-
-            // Apply collapses to current mesh
-            polyTopoChange newMeshMod(newMesh);
-
-            // Insert mesh refinement into polyTopoChange.
-            collapser.setRefinement(allPointInfo, newMeshMod);
-
-            Info<< indent << "Apply changes to the current mesh"
-                << decrIndent << endl;
-
-            // Apply changes to current mesh
-            autoPtr<mapPolyMesh> newMapPtr = newMeshMod.changeMesh
-            (
-                newMesh,
-                false
-            );
-            const mapPolyMesh& newMap = newMapPtr();
-
-            // Update fields
-            newMesh.updateMesh(newMap);
-            if (newMap.hasMotionPoints())
-            {
-                newMesh.movePoints(newMap.preMotionPoints());
-            }
-
-            // Relabel the boundary points
-//            labelList newBoundaryPoints(newMesh.nPoints(), -1);
-//            forAll(newBoundaryPoints, pI)
-//            {
-//                const label newToOldptI= map.pointMap()[pI];
-//                newBoundaryPoints[pI] = boundaryIOPts[newToOldptI];
-//            }
-//            boundaryIOPts = newBoundaryPoints;
-
-            // Synchronise the factors
-            mapOldMeshEdgeFieldToNewMesh
-            (
-                newMesh,
-                newMap.pointMap(),
-                newMeshMinEdgeLen
-            );
-
-            updateOldToNewPointMap
-            (
-                newMap.reversePointMap(),
-                origToCurrentPointMap
-            );
         }
 
-
         // Mesh check
         // ~~~~~~~~~~~~~~~~~~
         // Do not allow collapses in regions of error.
         // Updates minEdgeLen, nRelaxedEdges
 
-        if (controlMeshQuality_)
+        if (controlMeshQuality())
         {
             PackedBoolList isErrorPoint(newMesh.nPoints());
             nBadFaces = edgeCollapser::checkMeshQuality
             (
                 newMesh,
-                meshQualityCoeffDict_,
+                meshQualityCoeffDict(),
                 isErrorPoint
             );
 
@@ -1519,14 +1055,6 @@ Foam::label Foam::polyMeshFilter::filterFaceZone(const faceZone& fZone)
                 isErrorPoint,
                 pointErrorCount
             );
-
-            checkMeshFacesAndRelaxEdges
-            (
-                newMesh,
-                origToCurrentPointMap,
-                isErrorPoint,
-                pointErrorCount
-            );
         }
         else
         {
@@ -1544,4 +1072,11 @@ const Foam::autoPtr<Foam::fvMesh>& Foam::polyMeshFilter::filteredMesh() const
 }
 
 
+const Foam::autoPtr<Foam::labelList>&
+Foam::polyMeshFilter::pointPriority() const
+{
+    return pointPriority_;
+}
+
+
 // ************************************************************************* //
diff --git a/src/dynamicMesh/polyMeshFilter/polyMeshFilter.H b/src/dynamicMesh/polyMeshFilter/polyMeshFilter.H
index 0f26687608..95cddd0058 100644
--- a/src/dynamicMesh/polyMeshFilter/polyMeshFilter.H
+++ b/src/dynamicMesh/polyMeshFilter/polyMeshFilter.H
@@ -43,6 +43,7 @@ SourceFiles
 #include "List.H"
 #include "autoPtr.H"
 #include "scalarField.H"
+#include "polyMeshFilterSettings.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -59,6 +60,8 @@ class faceZone;
 \*---------------------------------------------------------------------------*/
 
 class polyMeshFilter
+:
+    private polyMeshFilterSettings
 {
     // Private data
 
@@ -68,52 +71,9 @@ class polyMeshFilter
         //- Copy of the original mesh to perform the filtering on
         autoPtr<fvMesh> newMeshPtr_;
 
-        //- Dictionary containing the coefficient sub-dictionaries
-        const IOdictionary dict_;
-
-        //- After collapsing, check the mesh quality and redo the collapsing
-        //  iteration if there are too many bad faces in the mesh
-        Switch controlMeshQuality_;
-
-        //- Coefficients for collapsing edges
-        const dictionary& collapseEdgesCoeffDict_;
-
-        //- Coefficients for collapsing faces
-        dictionary collapseFacesCoeffDict_;
-
-        //- Coefficients for controlling the mesh quality
-        dictionary meshQualityCoeffDict_;
-
-        //- Remove edges shorter than this length
-        const scalar minLen_;
-
-        //- Merge points that are only attached to two edges and have an angle
-        //  between the edge greater than this value
-        const scalar maxCos_;
-
-        //- The amount that the local minimum edge length will be reduced by if
-        //  the edge is part of a collapse string that generates poor quality
-        //  faces
-        const scalar edgeReductionFactor_;
-
-        //- Maximum number of outer iterations
-        const label maxIterations_;
-
-        //- Maximum number of smoothing iterations of minEdgeLen_ and
-        //  faceFilterFactor_
-        const label maxSmoothIters_;
-
-        //- Initialisation value of faceFilterFactor_
-        const scalar initialFaceLengthFactor_;
-
-        //- The amount that the local face size factor will be reduced by if
-        //  the face is part of a collapse string that generates poor quality
-        //  faces
-        const scalar faceReductionFactor_;
-
-        //- Maximum number of times a deleted point can be associated with the
-        //  creation of a bad face it is forced to be kept.
-        const label maxPointErrorCount_;
+        //-
+        labelList originalPointPriority_;
+        autoPtr<labelList> pointPriority_;
 
         //- The minimum edge length for each edge
         scalarField minEdgeLen_;
@@ -124,6 +84,22 @@ class polyMeshFilter
 
     // Private Member Functions
 
+        label filterFacesLoop(const label nOriginalBadFaces);
+
+        label filterFaces
+        (
+            polyMesh& newMesh,
+            scalarField& newMeshFaceFilterFactor,
+            labelList& origToCurrentPointMap
+        );
+
+        label filterEdges
+        (
+            polyMesh& newMesh,
+            scalarField& newMeshMinEdgeLen,
+            labelList& origToCurrentPointMap
+        );
+
         //- Increment pointErrorCount for points attached to a bad face
         void updatePointErrorCount
         (
@@ -160,7 +136,11 @@ class polyMeshFilter
         // + >0 : point on a processor patch with that ID
         // @todo Need to mark boundaryEdges as well, as an edge may have two
         //       boundary points but not itself lie on a boundary
-        labelList findBoundaryPoints(const polyMesh& mesh) const;
+        void updatePointPriorities
+        (
+            const polyMesh& newMesh,
+            const labelList& pointMap
+        );
 
         //- Print min/mean/max data for a field
         void printScalarFieldStats
@@ -213,6 +193,9 @@ public:
         //- Construct from fvMesh
         explicit polyMeshFilter(const fvMesh& mesh);
 
+        //- Construct from fvMesh and a label list of point priorities
+        polyMeshFilter(const fvMesh& mesh, const labelList& pointPriority);
+
 
     //- Destructor
     ~polyMeshFilter();
@@ -226,6 +209,8 @@ public:
             //  mesh has actually been filtered.
             const autoPtr<fvMesh>& filteredMesh() const;
 
+            const autoPtr<labelList>& pointPriority() const;
+
 
         // Edit
 
@@ -235,11 +220,11 @@ public:
             //- Filter edges and faces
             label filter(const label nOriginalBadFaces);
 
+            //- Filter all faces that are in the face zone indirectPatchFaces
+            label filter(const faceZone& fZone);
+
             //- Filter edges only.
             label filterEdges(const label nOriginalBadFaces);
-
-            //- Filter all faces that are in the face zone indirectPatchFaces
-            label filterFaceZone(const faceZone& fZone);
 };
 
 
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C
index 6ea1e748b0..cb6334306e 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C
@@ -182,53 +182,46 @@ void Foam::edgeCollapser::collapseToEdge
     Map<point>& collapsePointToLocation
 ) const
 {
-    const face& f = mesh_.faces()[faceI];
-
     // Negative half
 
-    Foam::point collapseToPtA =
-        collapseAxis*(sum(dNeg)/dNeg.size() - dShift) + fC;
+    Foam::point collapseToPtA(GREAT, GREAT, GREAT);
+        //collapseAxis*(sum(dNeg)/dNeg.size() - dShift) + fC;
 
-    DynamicList<label> faceBoundaryPts(f.size());
-    DynamicList<label> faceFeaturePts(f.size());
+    label maxPriority = labelMin;
+    DynamicList<label> maxPriorityPts(max(dNeg.size(), dPos.size()));
 
     forAll(facePtsNeg, fPtI)
     {
-        if (pointPriority[facePtsNeg[fPtI]] == 1)
+        const label facePointI = facePtsNeg[fPtI];
+        const label facePtPriority = pointPriority[facePointI];
+
+        if (facePtPriority > maxPriority)
         {
-            faceFeaturePts.append(facePtsNeg[fPtI]);
+            maxPriority = facePtPriority;
+            maxPriorityPts.clear();
+            maxPriorityPts.append(facePointI);
         }
-        else if (pointPriority[facePtsNeg[fPtI]] == 0)
+        else if (facePtPriority == maxPriority)
         {
-            faceBoundaryPts.append(facePtsNeg[fPtI]);
+            maxPriorityPts.append(facePointI);
         }
     }
 
-    if (!faceBoundaryPts.empty() || !faceFeaturePts.empty())
+    if (!maxPriorityPts.empty())
     {
-        if (!faceFeaturePts.empty())
-        {
-            collapseToPtA = pts[faceFeaturePts.first()];
-        }
-        else if (faceBoundaryPts.size() == 2)
-        {
-            collapseToPtA =
-                0.5
-               *(
-                    pts[faceBoundaryPts[0]]
-                  + pts[faceBoundaryPts[1]]
-                );
-        }
-        else if (faceBoundaryPts.size() <= f.size())
-        {
-            face bFace(faceBoundaryPts);
+        Foam::point averagePt(vector::zero);
 
-            collapseToPtA = bFace.centre(pts);
+        forAll(maxPriorityPts, ptI)
+        {
+            averagePt += pts[maxPriorityPts[ptI]];
         }
+
+        collapseToPtA = averagePt/maxPriorityPts.size();
+//        collapseToPtA = pts[maxPriorityPts.first()];
     }
 
-    faceFeaturePts.clear();
-    faceBoundaryPts.clear();
+    maxPriority = labelMin;
+    maxPriorityPts.clear();
 
     labelList faceEdgesNeg = edgesFromPoints(faceI, facePtsNeg);
 
@@ -244,47 +237,37 @@ void Foam::edgeCollapser::collapseToEdge
 
 
     // Positive half
-
-    Foam::point collapseToPtB
-        = collapseAxis*(sum(dPos)/dPos.size() - dShift) + fC;
+    Foam::point collapseToPtB(GREAT, GREAT, GREAT);
+//        = collapseAxis*(sum(dPos)/dPos.size() - dShift) + fC;
 
     forAll(facePtsPos, fPtI)
     {
-        if (pointPriority[facePtsPos[fPtI]] == 1)
+        const label facePointI = facePtsPos[fPtI];
+        const label facePtPriority = pointPriority[facePointI];
+
+        if (facePtPriority > maxPriority)
         {
-            faceFeaturePts.append(facePtsPos[fPtI]);
+            maxPriority = facePtPriority;
+            maxPriorityPts.clear();
+            maxPriorityPts.append(facePointI);
         }
-        else if (pointPriority[facePtsPos[fPtI]] == 0)
+        else if (facePtPriority == maxPriority)
         {
-            // If there is a point which is on the boundary,
-            // use it as the point to collapse others to, will
-            // use the first boundary point encountered if
-            // there are multiple boundary points.
-            faceBoundaryPts.append(facePtsPos[fPtI]);
+            maxPriorityPts.append(facePointI);
         }
     }
 
-    if (!faceBoundaryPts.empty() || !faceFeaturePts.empty())
+    if (!maxPriorityPts.empty())
     {
-        if (!faceFeaturePts.empty())
-        {
-            collapseToPtB = pts[faceFeaturePts.first()];
-        }
-        else if (faceBoundaryPts.size() == 2)
-        {
-            collapseToPtB =
-                0.5
-               *(
-                    pts[faceBoundaryPts[0]]
-                  + pts[faceBoundaryPts[1]]
-                );
-        }
-        else if (faceBoundaryPts.size() <= f.size())
-        {
-            face bFace(faceBoundaryPts);
+        Foam::point averagePt(vector::zero);
 
-            collapseToPtB = bFace.centre(pts);
+        forAll(maxPriorityPts, ptI)
+        {
+            averagePt += pts[maxPriorityPts[ptI]];
         }
+
+        collapseToPtB = averagePt/maxPriorityPts.size();
+//        collapseToPtB = pts[maxPriorityPts.first()];
     }
 
     labelList faceEdgesPos = edgesFromPoints(faceI, facePtsPos);
@@ -316,45 +299,79 @@ void Foam::edgeCollapser::collapseToPoint
 
     Foam::point collapseToPt = fC;
 
-    DynamicList<label> faceBoundaryPts(f.size());
-    DynamicList<label> faceFeaturePts(f.size());
+    label maxPriority = labelMin;
+    DynamicList<label> maxPriorityPts(f.size());
 
     forAll(facePts, fPtI)
     {
-        if (pointPriority[facePts[fPtI]] == 1)
+        const label facePointI = facePts[fPtI];
+        const label facePtPriority = pointPriority[facePointI];
+
+        if (facePtPriority > maxPriority)
         {
-            faceFeaturePts.append(facePts[fPtI]);
+            maxPriority = facePtPriority;
+            maxPriorityPts.clear();
+            maxPriorityPts.append(facePointI);
         }
-        else if (pointPriority[facePts[fPtI]] == 0)
+        else if (facePtPriority == maxPriority)
         {
-            faceBoundaryPts.append(facePts[fPtI]);
+            maxPriorityPts.append(facePointI);
         }
     }
 
-    if (!faceBoundaryPts.empty() || !faceFeaturePts.empty())
+    if (!maxPriorityPts.empty())
     {
-        if (!faceFeaturePts.empty())
-        {
-            collapseToPt = pts[faceFeaturePts.first()];
-        }
-        else if (faceBoundaryPts.size() == 2)
-        {
-            collapseToPt =
-                0.5
-               *(
-                    pts[faceBoundaryPts[0]]
-                  + pts[faceBoundaryPts[1]]
-                );
-        }
-        else if (faceBoundaryPts.size() <= f.size())
-        {
-            face bFace(faceBoundaryPts);
+        Foam::point averagePt(vector::zero);
 
-            collapseToPt = bFace.centre(pts);
+        forAll(maxPriorityPts, ptI)
+        {
+            averagePt += pts[maxPriorityPts[ptI]];
         }
+
+        collapseToPt = averagePt/maxPriorityPts.size();
+
+//        collapseToPt = pts[maxPriorityPts.first()];
     }
 
-    const labelList faceEdges = mesh_.faceEdges()[faceI];
+//    DynamicList<label> faceBoundaryPts(f.size());
+//    DynamicList<label> faceFeaturePts(f.size());
+//
+//    forAll(facePts, fPtI)
+//    {
+//        if (pointPriority[facePts[fPtI]] == 1)
+//        {
+//            faceFeaturePts.append(facePts[fPtI]);
+//        }
+//        else if (pointPriority[facePts[fPtI]] == 0)
+//        {
+//            faceBoundaryPts.append(facePts[fPtI]);
+//        }
+//    }
+//
+//    if (!faceBoundaryPts.empty() || !faceFeaturePts.empty())
+//    {
+//        if (!faceFeaturePts.empty())
+//        {
+//            collapseToPt = pts[faceFeaturePts.first()];
+//        }
+//        else if (faceBoundaryPts.size() == 2)
+//        {
+//            collapseToPt =
+//                0.5
+//               *(
+//                    pts[faceBoundaryPts[0]]
+//                  + pts[faceBoundaryPts[1]]
+//                );
+//        }
+//        else if (faceBoundaryPts.size() <= f.size())
+//        {
+//            face bFace(faceBoundaryPts);
+//
+//            collapseToPt = bFace.centre(pts);
+//        }
+//    }
+
+    const labelList& faceEdges = mesh_.faceEdges()[faceI];
 
     forAll(faceEdges, eI)
     {
@@ -740,35 +757,51 @@ Foam::label Foam::edgeCollapser::edgeMaster
 {
     label masterPoint = -1;
 
-    label e0 = e.start();
-    label e1 = e.end();
+    const label e0 = e.start();
+    const label e1 = e.end();
 
-    // Collapse edge to point with higher priority.
-    if (pointPriority[e0] >= 0)
+    const label e0Priority = pointPriority[e0];
+    const label e1Priority = pointPriority[e1];
+
+    if (e0Priority > e1Priority)
     {
-        if (pointPriority[e1] >= 0)
-        {
-            // Both points have high priority. Choose one to collapse to.
-            // Note: should look at feature edges/points!
-            masterPoint = e0;
-        }
-        else
-        {
-            masterPoint = e0;
-        }
+        masterPoint = e0;
     }
-    else
+    else if (e0Priority < e1Priority)
     {
-        if (pointPriority[e1] >= 0)
-        {
-            masterPoint = e1;
-        }
-        else
-        {
-            // None on boundary. Neither is a master.
-            return -1;
-        }
+        masterPoint = e1;
     }
+    else if (e0Priority == e1Priority)
+    {
+        masterPoint = e0;
+    }
+
+//    // Collapse edge to point with higher priority.
+//    if (pointPriority[e0] >= 0)
+//    {
+//        if (pointPriority[e1] >= 0)
+//        {
+//            // Both points have high priority. Choose one to collapse to.
+//            // Note: should look at feature edges/points!
+//            masterPoint = e0;
+//        }
+//        else
+//        {
+//            masterPoint = e0;
+//        }
+//    }
+//    else
+//    {
+//        if (pointPriority[e1] >= 0)
+//        {
+//            masterPoint = e1;
+//        }
+//        else
+//        {
+//            // None on boundary. Neither is a master.
+//            return -1;
+//        }
+//    }
 
     return masterPoint;
 }
@@ -784,7 +817,10 @@ void Foam::edgeCollapser::checkBoundaryPointMergeEdges
 {
    const pointField& points = mesh_.points();
 
-   if (pointPriority[pointI] >= 0 && pointPriority[otherPointI] < 0)
+   const label e0Priority = pointPriority[pointI];
+   const label e1Priority = pointPriority[otherPointI];
+
+   if (e0Priority > e1Priority)
    {
        collapsePointToLocation.set
        (
@@ -792,7 +828,7 @@ void Foam::edgeCollapser::checkBoundaryPointMergeEdges
            points[pointI]
        );
    }
-   else
+   else if (e0Priority < e1Priority)
    {
        collapsePointToLocation.set
        (
@@ -800,6 +836,22 @@ void Foam::edgeCollapser::checkBoundaryPointMergeEdges
            points[otherPointI]
        );
    }
+   else // e0Priority == e1Priority
+   {
+       collapsePointToLocation.set
+       (
+           pointI,
+           points[otherPointI]
+       );
+
+//       Foam::point averagePt
+//       (
+//           0.5*(points[otherPointI] + points[pointI])
+//       );
+//
+//       collapsePointToLocation.set(pointI, averagePt);
+//       collapsePointToLocation.set(otherPointI, averagePt);
+   }
 }
 
 
@@ -1963,7 +2015,7 @@ Foam::labelPair Foam::edgeCollapser::markSmallSliverFaces
     {
         const face& f = faces[fI];
 
-        if (faceFilterFactor[fI] == 0)
+        if (faceFilterFactor[fI] <= 0)
         {
             continue;
         }
@@ -2030,7 +2082,7 @@ Foam::labelPair Foam::edgeCollapser::markFaceZoneEdges
 
         const face& f = faces[fI];
 
-        if (faceFilterFactor[fI] == 0)
+        if (faceFilterFactor[fI] <= 0)
         {
             continue;
         }