From e81813a9c2550427f41c938df51b6a98eafe4c70 Mon Sep 17 00:00:00 2001
From: laurence <laurence>
Date: Tue, 11 Dec 2012 17:12:27 +0000
Subject: [PATCH] ENH: Add edgeCollapser utility and libraries.

edgeCollapser collapses small edges and faces.
Works in parallel by using PointEdgeWave.
---
 .../mesh/advanced/collapseEdges/Make/files    |    1 -
 .../mesh/advanced/collapseEdges/Make/options  |    3 +-
 .../mesh/advanced/collapseEdges/collapseDict  |   85 +
 .../advanced/collapseEdges/collapseEdges.C    | 1420 +-----------
 .../mesh/advanced/modifyMesh/modifyMesh.C     |   30 +-
 .../extrude/extrudeMesh/extrudeMesh.C         |   24 +-
 .../extrude2DMesh/extrude2DMeshApp.C          |   23 +-
 src/dynamicMesh/Make/files                    |    2 +
 .../polyMeshFilter/polyMeshFilter.C           | 1082 +++++++++
 .../polyMeshFilter/polyMeshFilter.H           |  238 ++
 .../polyTopoChange/edgeCollapser.C            | 1962 ++++++++++++++---
 .../polyTopoChange/edgeCollapser.H            |  292 ++-
 .../pointEdgeCollapse/pointEdgeCollapse.C     |    8 +-
 .../pointEdgeCollapse/pointEdgeCollapse.H     |    8 +-
 .../pointEdgeCollapse/pointEdgeCollapseI.H    |   76 +-
 15 files changed, 3428 insertions(+), 1826 deletions(-)
 create mode 100644 applications/utilities/mesh/advanced/collapseEdges/collapseDict
 create mode 100644 src/dynamicMesh/polyMeshFilter/polyMeshFilter.C
 create mode 100644 src/dynamicMesh/polyMeshFilter/polyMeshFilter.H
 rename {applications/utilities/mesh/advanced/collapseEdges => src/dynamicMesh/polyTopoChange/polyTopoChange}/pointEdgeCollapse/pointEdgeCollapse.C (88%)
 rename {applications/utilities/mesh/advanced/collapseEdges => src/dynamicMesh/polyTopoChange/polyTopoChange}/pointEdgeCollapse/pointEdgeCollapse.H (96%)
 rename {applications/utilities/mesh/advanced/collapseEdges => src/dynamicMesh/polyTopoChange/polyTopoChange}/pointEdgeCollapse/pointEdgeCollapseI.H (85%)

diff --git a/applications/utilities/mesh/advanced/collapseEdges/Make/files b/applications/utilities/mesh/advanced/collapseEdges/Make/files
index a15838abe8..d89ca6e737 100644
--- a/applications/utilities/mesh/advanced/collapseEdges/Make/files
+++ b/applications/utilities/mesh/advanced/collapseEdges/Make/files
@@ -1,4 +1,3 @@
 collapseEdges.C
-pointEdgeCollapse/pointEdgeCollapse.C
 
 EXE = $(FOAM_APPBIN)/collapseEdges
diff --git a/applications/utilities/mesh/advanced/collapseEdges/Make/options b/applications/utilities/mesh/advanced/collapseEdges/Make/options
index d1efa61fd5..987eae5ed7 100644
--- a/applications/utilities/mesh/advanced/collapseEdges/Make/options
+++ b/applications/utilities/mesh/advanced/collapseEdges/Make/options
@@ -1,8 +1,7 @@
 EXE_INC = \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -IpointEdgeCollapse
+    -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
     -ldynamicMesh \
diff --git a/applications/utilities/mesh/advanced/collapseEdges/collapseDict b/applications/utilities/mesh/advanced/collapseEdges/collapseDict
new file mode 100644
index 0000000000..78905475e3
--- /dev/null
+++ b/applications/utilities/mesh/advanced/collapseEdges/collapseDict
@@ -0,0 +1,85 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      http://www.openfoam.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+FoamFile
+{
+    version         2.0;
+    format          ascii;
+
+    root            "";
+    case            "";
+    instance        "";
+    local           "";
+
+    class           dictionary;
+    object          collapseDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+collapseEdgesCoeffs
+{
+    // Edges shorter than this absolute value will be merged
+    minimumEdgeLength   1e-8;
+
+    // The maximum angle between two edges that share a point attached to
+    // no other edges
+    maximumMergeAngle   30;
+
+    // The amount that minimumEdgeLength will be reduced by for each
+    // edge if that edge's collapse generates a poor quality face
+    reductionFactor                         0.5;
+}
+
+
+collapseFacesCoeffs
+{
+    // The initial face length factor
+    initialFaceLengthFactor                 0.5;
+
+    // The amount that initialFaceLengthFactor will be reduced by for each
+    // face if its collapse generates a poor quality face
+    reductionFactor                         $initialFaceLengthFactor;
+
+    // If the face can't be collapsed to an edge, and it has a span less than
+    // the target face length multiplied by this coefficient, collapse it
+    // to a point.
+    maxCollapseFaceToPointSideLengthCoeff   0.3;
+
+    // Allow early collapse of edges to a point
+    allowEarlyCollapseToPoint               on;
+
+    // Fraction to premultiply maxCollapseFaceToPointSideLengthCoeff by if
+    // allowEarlyCollapseToPoint is enabled
+    allowEarlyCollapseCoeff                 0.2;
+
+    // Defining how close to the midpoint (M) of the projected
+    // vertices line a projected vertex (X) can be before making this
+    // an invalid edge collapse
+    //
+    // X---X-g----------------M----X-----------g----X--X
+    //
+    // Only allow a collapse if all projected vertices are outwith
+    // guardFraction (g) of the distance form the face centre to the
+    // furthest vertex in the considered direction
+    guardFraction                           0.1;
+}
+
+
+meshQualityCoeffs
+{
+    // Name of the dictionary that has the mesh quality coefficients used
+    // by motionSmoother::checkMesh
+    meshQualityCoeffDict    meshQualityDict;
+
+    // Maximum number of outer iterations is mesh quality checking is enabled
+    maximumIterations       30;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C b/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C
index c0523bf43b..132a602b1a 100644
--- a/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C
+++ b/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,1425 +27,96 @@ Description
     - collapse short edges. Length of edges to collapse provided as argument.
     - merge two edges if they are in line. Maximum angle provided as argument.
     - remove unused points.
+    - collapse faces:
+        - with small areas to a single point
+        - that have a high aspect ratio (i.e. sliver face) to a single edge
 
-    Optionally removes cells. Can remove faces and points but does not check
-    for nonsense resulting topology.
+    Optionally checks the resulting mesh for bad faces and reduces the desired
+    face length factor for those faces attached to the bad faces.
 
     When collapsing an edge with one point on the boundary it will leave
     the boundary point intact. When both points inside it chooses random. When
     both points on boundary random again.
 
 Usage
-    - collapseEdges <edge length> <merge angle> [OPTION]
-
-    \param -allowCellCollapse \n
-    Allow collapsing of cells to a single point
-
-    \param -meshQuality \n
-    Only collapse if not exceeding user defined (in \a system/meshQualityDict)
-    quality settings
+    - collapseEdges [OPTION]
 
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
 #include "Time.H"
+#include "timeSelector.H"
 #include "polyTopoChange.H"
 #include "fvMesh.H"
-#include "mapPolyMesh.H"
-#include "mathematicalConstants.H"
-#include "PackedBoolList.H"
-#include "unitConversion.H"
-#include "globalMeshData.H"
-#include "globalIndex.H"
-#include "PointEdgeWave.H"
-#include "pointEdgeCollapse.H"
-#include "motionSmoother.H"
-#include "removePoints.H"
+#include "polyMeshFilter.H"
 
 using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-void filterFace
-(
-    const label faceI,
-    face& f,
-    const List<pointEdgeCollapse>& allPointInfo,
-    const Map<DynamicList<label> >& collapseStrings
-)
-{
-    label newFp = 0;
-
-    face oldFace = f;
-
-    forAll(f, fp)
-    {
-        label pointI = f[fp];
-
-        label collapsePoint = allPointInfo[pointI].collapseIndex();
-
-        if (collapseStrings.found(collapsePoint))
-        {
-            collapsePoint = collapseStrings[collapsePoint][0];
-        }
-
-        if (collapsePoint == -1)
-        {
-            WarningIn
-                (
-                    "filterFace"
-                    "(const label, face&, const List<pointEdgeCollapse>&)"
-                )
-                << "Point " << pointI << " was not visited by PointEdgeWave"
-                << endl;
-        }
-        else if (collapsePoint == -2)
-        {
-            f[newFp++] = pointI;
-        }
-        else
-        {
-            if (findIndex(SubList<label>(f, newFp), collapsePoint) == -1)
-            {
-                f[newFp++] = collapsePoint;
-            }
-        }
-    }
-
-
-    // Check for pinched face. Tries to correct
-    // - consecutive duplicate vertex. Removes duplicate vertex.
-    // - duplicate vertex with one other vertex in between (spike).
-    // Both of these should not really occur! and should be checked before
-    // collapsing edges.
-
-    const label size = newFp;
-
-    newFp = 2;
-
-    for (label fp = 2; fp < size; fp++)
-    {
-        label fp1 = fp-1;
-        label fp2 = fp-2;
-
-        label pointI = f[fp];
-
-        // Search for previous occurrence.
-        label index = findIndex(SubList<label>(f, fp), pointI);
-
-        if (index == fp1)
-        {
-            WarningIn
-            (
-                "Foam::edgeCollapser::filterFace(const label faceI, "
-                "face& f) const"
-            )   << "Removing consecutive duplicate vertex in face "
-                << f << endl;
-            // Don't store current pointI
-        }
-        else if (index == fp2)
-        {
-            WarningIn
-            (
-                "Foam::edgeCollapser::filterFace(const label faceI, "
-                "face& f) const"
-            )   << "Removing non-consecutive duplicate vertex in face "
-                << f << endl;
-            // Don't store current pointI and remove previous
-            newFp--;
-        }
-        else if (index != -1)
-        {
-            WarningIn
-            (
-                "Foam::edgeCollapser::filterFace(const label faceI, "
-                "face& f) const"
-            )   << "Pinched face " << f << endl;
-            f[newFp++] = pointI;
-        }
-        else
-        {
-            f[newFp++] = pointI;
-        }
-    }
-
-    f.setSize(newFp);
-}
-
-
-bool setRefinement
-(
-    const polyMesh& mesh,
-    polyTopoChange& meshMod,
-    const List<pointEdgeCollapse>& allPointInfo
-)
-{
-    const cellList& cells = mesh.cells();
-    const labelList& faceOwner = mesh.faceOwner();
-    const labelList& faceNeighbour = mesh.faceNeighbour();
-    const labelListList& pointFaces = mesh.pointFaces();
-    const pointZoneMesh& pointZones = mesh.pointZones();
-
-    globalIndex globalStrings(mesh.nPoints());
-
-    boolList removedPoints(mesh.nPoints(), false);
-
-    // Create strings of edges.
-    // Map from collapseIndex(=global master point) to set of points
-    Map<DynamicList<label> > collapseStrings;
-    {
-        // 1. Count elements per collapseIndex
-        Map<label> nPerIndex(mesh.nPoints()/10);
-        forAll(allPointInfo, pointI)
-        {
-            const label collapseIndex = allPointInfo[pointI].collapseIndex();
-
-            if (collapseIndex != -1 && collapseIndex != -2)
-            {
-                Map<label>::iterator fnd = nPerIndex.find(collapseIndex);
-                if (fnd != nPerIndex.end())
-                {
-                    fnd()++;
-                }
-                else
-                {
-                    nPerIndex.insert(collapseIndex, 1);
-                }
-            }
-        }
-
-        // 2. Size
-        collapseStrings.resize(2*nPerIndex.size());
-        forAllConstIter(Map<label>, nPerIndex, iter)
-        {
-            collapseStrings.insert(iter.key(), DynamicList<label>(iter()));
-        }
-
-        // 3. Fill
-        forAll(allPointInfo, pointI)
-        {
-            const label collapseIndex = allPointInfo[pointI].collapseIndex();
-
-            if (collapseIndex != -1 && collapseIndex != -2)
-            {
-                collapseStrings[collapseIndex].append(pointI);
-            }
-        }
-    }
-
-
-    bool meshChanged = false;
-
-    // Current faces (is also collapseStatus: f.size() < 3)
-    faceList newFaces(mesh.faces());
-
-    // Current cellCollapse status
-    boolList cellRemoved(mesh.nCells(), false);
-
-    label nUnvisited = 0;
-    label nUncollapsed = 0;
-    label nCollapsed = 0;
-
-    forAll(allPointInfo, pI)
-    {
-        const pointEdgeCollapse& pec = allPointInfo[pI];
-
-        if (pec.collapseIndex() == -1)
-        {
-            nUnvisited++;
-        }
-        else if (pec.collapseIndex() == -2)
-        {
-            nUncollapsed++;
-        }
-        else if (pec.collapseIndex() >= 0)
-        {
-            nCollapsed++;
-        }
-    }
-
-    label nPoints = allPointInfo.size();
-
-    reduce(nPoints, sumOp<label>());
-    reduce(nUnvisited, sumOp<label>());
-    reduce(nUncollapsed, sumOp<label>());
-    reduce(nCollapsed, sumOp<label>());
-
-    Info<< incrIndent;
-    Info<< indent << "Number of points : " << nPoints << nl
-        << indent << "Not visited      : " << nUnvisited << nl
-        << indent << "Not collapsed    : " << nUncollapsed << nl
-        << indent << "Collapsed        : " << nCollapsed << nl
-        << endl;
-    Info<< decrIndent;
-
-    do
-    {
-        // Update face collapse from edge collapses
-        forAll(newFaces, faceI)
-        {
-            filterFace(faceI, newFaces[faceI], allPointInfo, collapseStrings);
-        }
-
-        // Check if faces to be collapsed cause cells to become collapsed.
-        label nCellCollapsed = 0;
-
-        forAll(cells, cellI)
-        {
-            if (!cellRemoved[cellI])
-            {
-                const cell& cFaces = cells[cellI];
-
-                label nFaces = cFaces.size();
-
-                forAll(cFaces, i)
-                {
-                    label faceI = cFaces[i];
-
-                    if (newFaces[faceI].size() < 3)
-                    {
-                        --nFaces;
-
-                        if (nFaces < 4)
-                        {
-                            Pout<< "Cell:" << cellI
-                                << " uses faces:" << cFaces
-                                << " of which too many are marked for removal:"
-                                << endl
-                                << "   ";
-                            forAll(cFaces, j)
-                            {
-                                if (newFaces[cFaces[j]].size() < 3)
-                                {
-                                    Pout<< ' '<< cFaces[j];
-                                }
-                            }
-                            Pout<< endl;
-
-                            cellRemoved[cellI] = true;
-
-                            // Collapse all edges of cell to nothing
-                            //collapseEdges(cellEdges[cellI]);
-
-                            nCellCollapsed++;
-
-                            break;
-                        }
-                    }
-                }
-            }
-        }
-
-        if (nCellCollapsed == 0)
-        {
-            break;
-        }
-    } while (true);
-
-
-    // Keep track of faces that have been done already.
-    boolList doneFace(mesh.nFaces(), false);
-
-    {
-        // Mark points used.
-        boolList usedPoint(mesh.nPoints(), false);
-
-        forAll(cellRemoved, cellI)
-        {
-            if (cellRemoved[cellI])
-            {
-                meshMod.removeCell(cellI, -1);
-            }
-        }
-
-        // Remove faces
-        forAll(newFaces, faceI)
-        {
-            const face& f = newFaces[faceI];
-
-            if (f.size() < 3)
-            {
-                meshMod.removeFace(faceI, -1);
-                meshChanged = true;
-
-                // Mark face as been done.
-                doneFace[faceI] = true;
-            }
-            else
-            {
-                // Kept face. Mark vertices
-                forAll(f, fp)
-                {
-                    usedPoint[f[fp]] = true;
-                }
-            }
-        }
-
-        // Remove unused vertices that have not been marked for removal already
-        forAll(usedPoint, pointI)
-        {
-            if (!usedPoint[pointI])
-            {
-                removedPoints[pointI] = true;
-                meshMod.removePoint(pointI, -1);
-                meshChanged = true;
-            }
-        }
-    }
-
-    // Modify the point location of the remaining points
-    forAll(allPointInfo, pointI)
-    {
-        const label collapseIndex = allPointInfo[pointI].collapseIndex();
-        const point& collapsePoint = allPointInfo[pointI].collapsePoint();
-
-        if
-        (
-            removedPoints[pointI] == false
-         && collapseIndex != -1
-         && collapseIndex != -2
-        )
-        {
-            meshMod.modifyPoint
-            (
-                pointI,
-                collapsePoint,
-                pointZones.whichZone(pointI),
-                false
-            );
-        }
-    }
-
-
-    const polyBoundaryMesh& boundaryMesh = mesh.boundaryMesh();
-    const faceZoneMesh& faceZones = mesh.faceZones();
-
-    // Renumber faces that use points
-    forAll(allPointInfo, pointI)
-    {
-        if (removedPoints[pointI] == true)
-        {
-            const labelList& changedFaces = pointFaces[pointI];
-
-            forAll(changedFaces, changedFaceI)
-            {
-                label faceI = changedFaces[changedFaceI];
-
-                if (!doneFace[faceI])
-                {
-                    doneFace[faceI] = true;
-
-                    // Get current zone info
-                    label zoneID = faceZones.whichZone(faceI);
-
-                    bool zoneFlip = false;
-
-                    if (zoneID >= 0)
-                    {
-                        const faceZone& fZone = faceZones[zoneID];
-
-                        zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
-                    }
-
-                    // Get current connectivity
-                    label own = faceOwner[faceI];
-                    label nei = -1;
-                    label patchID = -1;
-
-                    if (mesh.isInternalFace(faceI))
-                    {
-                        nei = faceNeighbour[faceI];
-                    }
-                    else
-                    {
-                        patchID = boundaryMesh.whichPatch(faceI);
-                    }
-
-                    meshMod.modifyFace
-                    (
-                        newFaces[faceI],            // face
-                        faceI,                      // faceI to change
-                        own,                        // owner
-                        nei,                        // neighbour
-                        false,                      // flipFaceFlux
-                        patchID,                    // patch
-                        zoneID,
-                        zoneFlip
-                    );
-                    meshChanged = true;
-                }
-            }
-        }
-    }
-
-    return meshChanged;
-}
-
-
-// Get faceEdges in order of face points, i.e. faceEdges[0] is between
-// f[0] and f[1]
-labelList getSortedEdges
-(
-    const edgeList& edges,
-    const labelList& f,
-    const labelList& edgeLabels
-)
-{
-    labelList faceEdges(edgeLabels.size(), -1);
-
-    // Find starting pos in f for every edgeLabels
-    forAll(edgeLabels, i)
-    {
-        label edgeI = edgeLabels[i];
-
-        const edge& e = edges[edgeI];
-
-        label fp = findIndex(f, e[0]);
-        label fp1 = f.fcIndex(fp);
-
-        if (f[fp1] == e[1])
-        {
-            // EdgeI between fp -> fp1
-            faceEdges[fp] = edgeI;
-        }
-        else
-        {
-            // EdgeI between fp-1 -> fp
-            faceEdges[f.rcIndex(fp)] = edgeI;
-        }
-    }
-
-    return faceEdges;
-}
-
-
-// Return master point edge needs to be collapsed to (or -1)
-label edgeMaster
-(
-    const labelList& boundaryPoint,
-    const bool flipEdge,
-    const edge& e
-)
-{
-    label masterPoint = -1;
-
-    label e0 = e[0];
-    label e1 = e[1];
-
-    if (flipEdge)
-    {
-        e0 = e[1];
-        e1 = e[0];
-    }
-
-    // Check if one of the points is on a processor
-//    if
-//    (
-//        boundaryPoint[e0] > 0
-//     && boundaryPoint[e1] > 0
-//    )
-//    {
-//        if (boundaryPoint[e0] != boundaryPoint[e1])
-//        {
-//            return -1;
-//        }
-//    }
-//
-//    if (boundaryPoint[e0] > 0)
-//    {
-//        return e0;
-//    }
-//    else if (boundaryPoint[e1] > 0)
-//    {
-//        return e1;
-//    }
-
-    // Collapse edge to boundary point.
-    if (boundaryPoint[e0] == 0)
-    {
-        if (boundaryPoint[e1] == 0)
-        {
-            // Both points on boundary. Choose one to collapse to.
-            // Note: should look at feature edges/points!
-            masterPoint = e0;
-        }
-        else
-        {
-            masterPoint = e0;
-        }
-    }
-    else
-    {
-        if (boundaryPoint[e1] == 0)
-        {
-            masterPoint = e1;
-        }
-        else
-        {
-            // None on boundary. Choose arbitrary.
-            // Note: should look at geometry?
-            masterPoint = e0;
-        }
-    }
-
-    return masterPoint;
-}
-
-
-// Create consistent set of collapses.
-//  collapseEdge : per edge:
-//      -1 : do not collapse
-//       0 : collapse to start
-//       1 : collapse to end
-//  Note: collapseEdge has to be parallel consistent (in orientation)
-label syncCollapse
-(
-    const polyMesh& mesh,
-    const globalIndex& globalStrings,
-    const labelList& collapseEdge,
-    List<pointEdgeCollapse>& allPointInfo
-)
-{
-    const edgeList& edges = mesh.edges();
-    const pointField& points = mesh.points();
-
-    label nCollapsed = 0;
-
-    DynamicList<label> initPoints(mesh.nPoints());
-    DynamicList<pointEdgeCollapse> initPointInfo(mesh.nPoints());
-    allPointInfo.resize(mesh.nPoints());
-
-    // Initialise edges to no collapse
-    List<pointEdgeCollapse> allEdgeInfo
-    (
-        mesh.nEdges(),
-        pointEdgeCollapse(vector::zero, -1)
-    );
-
-    // Mark selected edges for collapse
-    forAll(edges, edgeI)
-    {
-        const edge& e = edges[edgeI];
-
-        if (collapseEdge[edgeI] != -1)
-        {
-            label masterPointI = e[collapseEdge[edgeI]];
-
-            const pointEdgeCollapse pec
-            (
-                points[masterPointI],
-                globalStrings.toGlobal(masterPointI)
-            );
-
-            // Mark as collapsable but with nonsense master so it gets
-            // overwritten and starts an update wave
-            allEdgeInfo[edgeI] = pointEdgeCollapse
-            (
-                points[masterPointI],
-                labelMax
-            );
-
-            initPointInfo.append(pec);
-            initPoints.append(e.start());
-
-            initPointInfo.append(pec);
-            initPoints.append(e.end());
-
-            nCollapsed++;
-        }
-    }
-
-    PointEdgeWave<pointEdgeCollapse> collapsePropagator
-    (
-        mesh,
-        initPoints,
-        initPointInfo,
-        allPointInfo,
-        allEdgeInfo,
-        mesh.globalData().nTotalPoints()  // Maximum number of iterations
-    );
-
-    return nCollapsed;
-}
-
-
-void syncCollapseEdge(const polyMesh& mesh, labelList& collapseEdge)
-{
-    // Check whether edge point order is reversed from mesh to coupledPatch
-    const globalMeshData& globalData = mesh.globalData();
-    const mapDistribute& map = globalData.globalEdgeSlavesMap();
-    const labelList& coupledMeshEdges = globalData.coupledPatchMeshEdges();
-    const indirectPrimitivePatch& coupledPatch = globalData.coupledPatch();
-    const PackedBoolList& cppOrientation = globalData.globalEdgeOrientation();
-    PackedBoolList meshToPatchSameOrientation(coupledMeshEdges.size(), true);
-
-    forAll(coupledMeshEdges, eI)
-    {
-        const label meshEdgeIndex = coupledMeshEdges[eI];
-
-        if (collapseEdge[meshEdgeIndex] != -1)
-        {
-            const edge& meshEdge = mesh.edges()[meshEdgeIndex];
-            const edge& coupledPatchEdge = coupledPatch.edges()[eI];
-
-            if
-            (
-                meshEdge[0] == coupledPatch.meshPoints()[coupledPatchEdge[1]]
-             && meshEdge[1] == coupledPatch.meshPoints()[coupledPatchEdge[0]]
-            )
-            {
-                meshToPatchSameOrientation[eI] = false;
-            }
-        }
-    }
-
-
-    labelList cppEdgeData(map.constructSize());
-
-    forAll(coupledMeshEdges, eI)
-    {
-        const label meshEdgeIndex = coupledMeshEdges[eI];
-
-        cppEdgeData[eI] = collapseEdge[meshEdgeIndex];
-
-        if
-        (
-            (collapseEdge[meshEdgeIndex] != -1)
-         && (meshToPatchSameOrientation[eI] != cppOrientation[eI])
-        )
-        {
-            cppEdgeData[eI] = 1-cppEdgeData[eI];
-        }
-    }
-
-
-    // Synchronise cppEdgeData
-    // Use minEqOp reduction, so that edge will only be collapsed on processor
-    // boundary if both processors agree to collapse it
-    globalData.syncData
-    (
-        cppEdgeData,
-        globalData.globalEdgeSlaves(),
-        globalData.globalEdgeTransformedSlaves(),
-        map,
-        minEqOp<label>()
-    );
-
-
-    forAll(coupledMeshEdges, eI)
-    {
-        const label meshEdgeIndex = coupledMeshEdges[eI];
-
-        collapseEdge[meshEdgeIndex] = cppEdgeData[eI];
-
-        if
-        (
-            (cppEdgeData[eI] != -1)
-         && (meshToPatchSameOrientation[eI] != cppOrientation[eI])
-        )
-        {
-            collapseEdge[meshEdgeIndex] = 1-collapseEdge[meshEdgeIndex];
-        }
-    }
-}
-
-
-// Mark (in collapseEdge) any edges to collapse
-label collapseSmallEdges
-(
-    const polyMesh& mesh,
-    const labelList& boundaryPoint,
-    const scalarField& minEdgeLen,
-    labelList& collapseEdge
-)
-{
-    // Work out which edges to collapse
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    label nCollapsed = 0;
-
-    forAll(mesh.edges(), edgeI)
-    {
-        if (collapseEdge[edgeI] == -1)
-        {
-            const edge& e = mesh.edges()[edgeI];
-
-            if (e.mag(mesh.points()) < minEdgeLen[edgeI])
-            {
-                label masterPointI = edgeMaster(boundaryPoint, false, e);
-                if (masterPointI == e[0])
-                {
-                    collapseEdge[edgeI] = 0;
-                }
-                else
-                {
-                    collapseEdge[edgeI] = 1;
-                }
-                nCollapsed++;
-            }
-        }
-    }
-    return nCollapsed;
-}
-
-
-// Mark (in collapseEdge) any edges to merge
-label mergeEdges
-(
-    const polyMesh& mesh,
-    const scalar maxCos,
-    const labelList& boundaryPoint,
-    const scalarField& minEdgeLen,
-    labelList& collapseEdge
-)
-{
-    const edgeList& edges = mesh.edges();
-    const labelListList& pointEdges = mesh.pointEdges();
-
-    // Point removal engine
-    removePoints pointRemover(mesh, false);
-
-    // Find out points that can be deleted
-    boolList pointCanBeDeleted;
-    label nTotRemove = pointRemover.countPointUsage(maxCos, pointCanBeDeleted);
-
-
-    // Rework point-to-remove into edge-to-collapse.
-
-    label nCollapsed = 0;
-
-    if (nTotRemove > 0)
-    {
-        forAll(pointEdges, pointI)
-        {
-            if (pointCanBeDeleted[pointI])
-            {
-                const labelList& pEdges = pointEdges[pointI];
-
-                if (pEdges.size() == 2)
-                {
-                    // Always the case?
-
-                    label e0 = pEdges[0];
-                    label e1 = pEdges[1];
-
-                    if
-                    (
-                        collapseEdge[e0] == -1
-                     && minEdgeLen[e0] >= 0
-                     && collapseEdge[e1] == -1
-                     && minEdgeLen[e1] >= 0
-                    )
-                    {
-                        // Get the two vertices on both sides of the point
-                        label leftV = edges[e0].otherVertex(pointI);
-                        label rightV = edges[e1].otherVertex(pointI);
-
-                        // Can collapse pointI onto either leftV or rightV.
-                        // Preferentially choose an internal point to hopefully
-                        // give less distortion
-
-                        if (boundaryPoint[leftV] == -1)
-                        {
-                            collapseEdge[e0] = findIndex(edges[e0], leftV);
-                        }
-                        else
-                        {
-                            collapseEdge[e1] = findIndex(edges[e1], rightV);
-                        }
-                    }
-                }
-            }
-        }
-    }
-    return nCollapsed;
-}
-
-
-// Make consistent set of collapses that does not collapse any cells
-label consistentCollapse
-(
-    const bool allowCellCollapse,
-    const polyMesh& mesh,
-    const globalIndex& globalStrings,
-    labelList& collapseEdge,
-    List<pointEdgeCollapse>& allPointInfo
-)
-{
-    // Make sure we don't collapse cells
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    while (true)
-    {
-        // Sync collapseEdge
-        syncCollapseEdge(mesh, collapseEdge);
-
-
-        // Get collapsed faces
-
-        label nAdditionalCollapsed = 0;
-
-        PackedBoolList isCollapsedFace(mesh.nFaces());
-        forAll(mesh.faceEdges(), faceI)
-        {
-            const labelList& fEdges = mesh.faceEdges()[faceI];
-
-            // Count number of remaining edges
-            label nEdges = fEdges.size();
-            forAll(fEdges, fEdgeI)
-            {
-                label edgeI = fEdges[fEdgeI];
-                if (collapseEdge[edgeI] != -1)
-                {
-                    nEdges--;
-                }
-            }
-
-            if (nEdges < 3)
-            {
-                // Face is collapsed.
-                isCollapsedFace[faceI] = 1;
-
-                if (nEdges == 1)
-                {
-                    // Cannot collapse face down to single edge.
-
-                    //- Option 1: collapse remaining edge as well. However
-                    //  if this edge is on the coupled processor patch this
-                    //  logic clashes with that of syncCollapseEdge
-                    //  (do not collapse if any not collapse)
-                    //forAll(fEdges, fEdgeI)
-                    //{
-                    //    label edgeI = fEdges[fEdgeI];
-                    //    if (collapseEdge[edgeI] == -1)
-                    //    {
-                    //        collapseEdge[edgeI] = 0;
-                    //        nAdditionalCollapsed++;
-                    //    }
-                    //}
-
-                    //- Option 2: uncollapse this face.
-                    forAll(fEdges, fEdgeI)
-                    {
-                        label edgeI = fEdges[fEdgeI];
-                        if (collapseEdge[edgeI] != -1)
-                        {
-                            collapseEdge[edgeI] = -1;
-                            nAdditionalCollapsed++;
-                        }
-                    }
-                }
-            }
-        }
-
-        //Pout<< "nAdditionalCollapsed : " << nAdditionalCollapsed << endl;
-
-
-        label nUncollapsed = 0;
-
-        if (!allowCellCollapse)
-        {
-            // Check collapsed cells
-
-            forAll(mesh.cells(), cellI)
-            {
-                const cell& cFaces = mesh.cells()[cellI];
-                label nFaces = cFaces.size();
-                forAll(cFaces, i)
-                {
-                    label faceI = cFaces[i];
-                    if (isCollapsedFace[faceI])
-                    {
-                        nFaces--;
-                        if (nFaces < 4)
-                        {
-                            // Unmark this face for collapse
-                            const labelList& fEdges = mesh.faceEdges()[faceI];
-
-                            forAll(fEdges, fEdgeI)
-                            {
-                                label edgeI = fEdges[fEdgeI];
-                                if (collapseEdge[edgeI] != -1)
-                                {
-                                    collapseEdge[edgeI] = -1;
-                                    nUncollapsed++;
-                                }
-                            }
-
-                            // Uncollapsed this face.
-                            isCollapsedFace[faceI] = 0;
-                            nFaces++;
-                        }
-                    }
-                }
-            }
-            //Pout<< "** disallowing cells : nUncollapsed : "
-            //    << nUncollapsed << endl;
-        }
-
-
-        if
-        (
-            returnReduce(nUncollapsed+nAdditionalCollapsed, sumOp<label>())
-            == 0
-        )
-        {
-            break;
-        }
-    }
-
-
-    // Create consistent set of collapses
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    // Note: requires collapseEdge to be synchronised. (above loop makes sure
-    //       of that)
-
-    return syncCollapse(mesh, globalStrings, collapseEdge, allPointInfo);
-}
-
-
-// Check mesh and mark points on faces in error
-// Returns boolList with points in error set
-PackedBoolList checkMeshQuality(const polyMesh& mesh)
-{
-    //mesh.checkMesh(true);
-
-    IOdictionary meshQualityDict
-    (
-        IOobject
-        (
-            "meshQualityDict",
-            mesh.time().system(),
-            mesh.time(),
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    );
-
-    labelHashSet badFaces(mesh.nFaces()/100);
-    DynamicList<label> checkFaces(mesh.nFaces());
-
-    const vectorField& fAreas = mesh.faceAreas();
-
-    scalar faceAreaLimit = SMALL;
-
-    forAll(fAreas, fI)
-    {
-        if (mag(fAreas[fI]) > faceAreaLimit)
-        {
-            checkFaces.append(fI);
-        }
-    }
-
-    motionSmoother::checkMesh
-    (
-        false,
-        mesh,
-        meshQualityDict,
-        checkFaces,
-        badFaces
-    );
-
-    label nBadFaces = badFaces.size();
-    reduce(nBadFaces, sumOp<label>());
-
-    Info<< nl << "Number of bad faces          : " << nBadFaces << endl;
-
-    PackedBoolList isErrorPoint(mesh.nPoints());
-    forAllConstIter(labelHashSet, badFaces, iter)
-    {
-        const face& f = mesh.faces()[iter.key()];
-
-        forAll(f, pI)
-        {
-            isErrorPoint[f[pI]] = 1;
-        }
-    }
-
-    syncTools::syncPointList
-    (
-        mesh,
-        isErrorPoint,
-        orEqOp<unsigned int>(),
-        0
-    );
-
-    return isErrorPoint;
-}
-
-
-// Check mesh with collapses (newMesh), updates minEdgeLen, nFrozenEdges
-void checkMeshAndFreezeEdges
-(
-    const polyMesh& newMesh,
-    const labelList& oldToNewMesh,
-    const polyMesh& oldMesh,
-    scalarField& minEdgeLen,
-    label& nFrozenEdges
-)
-{
-    PackedBoolList isErrorPoint = checkMeshQuality(newMesh);
-
-    forAll(oldMesh.edges(), edgeI)
-    {
-        const edge& e = oldMesh.edges()[edgeI];
-        label newStart = oldToNewMesh[e[0]];
-        label newEnd = oldToNewMesh[e[1]];
-
-        if
-        (
-            (newStart >= 0 && isErrorPoint[newStart])
-         || (newEnd >= 0 && isErrorPoint[newEnd])
-        )
-        {
-            // Gradual decrease? For now just hard disable
-            if (minEdgeLen[edgeI] > -GREAT/2)
-            {
-                minEdgeLen[edgeI] = -GREAT;
-                nFrozenEdges++;
-            }
-        }
-    }
-}
-
-
-// Mark boundary points
-// boundaryPoint:
-// + -1 : point not on boundary
-// +  0 : point on a real boundary
-// + >0 : point on a processor patch with that ID
-labelList findBoundaryPoints(const polyMesh& mesh)
-{
-    const faceList& faces = mesh.faces();
-    const polyBoundaryMesh& bMesh = mesh.boundaryMesh();
-
-
-    labelList boundaryPoint(mesh.nPoints(), -1);
-
-    // Get all points on a boundary
-    label nIntFaces = mesh.nInternalFaces();
-    for (label faceI = nIntFaces; faceI < mesh.nFaces(); faceI++)
-    {
-        const face& f = faces[faceI];
-
-        forAll(f, fp)
-        {
-           boundaryPoint[f[fp]] = 0;
-        }
-    }
-
-    // Get all processor boundary points and the processor patch label
-    // that they are on.
-    forAll(bMesh, patchI)
-    {
-        const polyPatch& patch = bMesh[patchI];
-
-        if (isA<coupledPolyPatch>(patch))
-        {
-            if (patchI == 0)
-            {
-                // We mark 'normal' boundary points with 0 so make sure this
-                // coupled patch is not 0.
-                FatalErrorIn("findBoundaryPoints(const polyMesh&)")
-                    << "Your patches should have non-coupled ones before any"
-                    << " coupled ones. Current patches " << bMesh.names()
-                    << exit(FatalError);
-            }
-
-            forAll(patch, fI)
-            {
-                const face& f = patch[fI];
-
-                forAll(f, fp)
-                {
-                    boundaryPoint[f[fp]] = patchI;
-                }
-            }
-        }
-    }
-
-    return boundaryPoint;
-}
-
-
 // Main program:
 int main(int argc, char *argv[])
 {
+    timeSelector::addOptions(true, false);
     argList::addNote
     (
-        "Merges small and in-line edges.\n"
-        "Collapses faces and optionally cells to a point."
+        "Collapses small edges to a point.\n"
+        "Optionally collapse small faces to a point and thin faces to an edge."
+    );
+
+    argList::addBoolOption
+    (
+        "collapseFaces",
+        "Collapse small and sliver faces as well as small edges"
     );
 
 #   include "addOverwriteOption.H"
-
-    argList::validArgs.append("edge length [m]");
-    argList::validArgs.append("merge angle (degrees)");
-    argList::addBoolOption
-    (
-        "allowCellCollapse",
-        "Allow collapsing of cells to a single point"
-    );
-    argList::addBoolOption
-    (
-        "meshQuality",
-        "Only collapse if not exceeding given meshQualityDict limits"
-    );
-
-
-
 #   include "setRootCase.H"
 #   include "createTime.H"
+
     runTime.functionObjects().off();
-#   include "createPolyMesh.H"
+    instantList timeDirs = timeSelector::selectIfPresent(runTime, args);
+
+#   include "createMesh.H"
+
     const word oldInstance = mesh.pointsInstance();
 
-    scalar minLen = args.argRead<scalar>(1);
-    const scalar angle = args.argRead<scalar>(2);
-
-    const bool allowCellCollapse = args.optionFound("allowCellCollapse");
     const bool overwrite = args.optionFound("overwrite");
-    const bool checkQuality = args.optionFound("meshQuality");
 
-    scalar maxCos = Foam::cos(degToRad(angle));
+    const bool collapseFaces = args.optionFound("collapseFaces");
 
-    Info<< "Merging:" << nl
-        << "    edges with length less than " << minLen << " meters" << nl
-        << "    edges split by a point with edges in line to within " << angle
-        << " degrees" << nl
-        << endl;
-
-    if (allowCellCollapse)
+    forAll(timeDirs, timeI)
     {
-        Info<< "Allowing collapse of cells down to a point." << nl
-            << endl;
-    }
-    else
-    {
-        Info<< "Disallowing collapse of cells down to a point." << nl
-            << endl;
-    }
+        runTime.setTime(timeDirs[timeI], timeI);
 
+        Info<< "Time = " << runTime.timeName() << endl;
 
-    if (checkQuality)
-    {
-        Info<< "Selectively disabling wanted collapses until resulting quality"
-            << " satisfies constraints in system/meshQualityDict" << nl
-            << endl;
-    }
+        polyMeshFilter meshFilter(mesh);
 
+        // newMesh will be empty until it is filtered
+        const autoPtr<fvMesh>& newMesh = meshFilter.filteredMesh();
 
-
-    // To mark master of collapes
-    globalIndex globalStrings(mesh.nPoints());
-
-
-    // Local collapse length. Any edge below this length gets (attempted)
-    // collapsed. Use either aa gradually decreasing value
-    // (so from minLen to e.g. 0.5*minLen) or a hard limit (GREAT)
-    scalarField minEdgeLen(mesh.nEdges(), minLen);
-    label nFrozenEdges = 0;
-
-
-
-    // Initial mesh check
-    // ~~~~~~~~~~~~~~~~~~
-    // Do not allow collapses in regions of error.
-    // Updates minEdgeLen, nFrozenEdges
-    if (checkQuality)
-    {
-        checkMeshAndFreezeEdges
-        (
-            mesh,
-            identity(mesh.nPoints()),
-            mesh,
-            minEdgeLen,
-            nFrozenEdges
-        );
-        Info<< "Initial frozen edges "
-            << returnReduce(nFrozenEdges, sumOp<label>())
-            << " out of " << returnReduce(mesh.nEdges(), sumOp<label>())
-            << endl;
-    }
-
-
-    // Mark points on boundary
-    const labelList boundaryPoint = findBoundaryPoints(mesh);
-
-    // Copy of current set of topology changes. Used to generate final mesh.
-    polyTopoChange savedMeshMod(mesh.boundaryMesh().size());
-
-    // Keep track of whether mesh has changed at all
-    bool meshChanged = false;
-
-
-
-    // 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.
-
-    do
-    {
-        // Per edge collapse status:
-        // -1 : not collapsed
-        //  0 : collapse to start
-        //  1 : collapse to end
-        labelList collapseEdge(mesh.nEdges(), -1);
-
-
-        // Work out which edges to collapse
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-        // This is by looking at minEdgeLen (to avoid frozen edges)
-        // and marking in collapseEdge.
-
-        // Small edges
-        label nSmallCollapsed = collapseSmallEdges
-        (
-            mesh,
-            boundaryPoint,
-            minEdgeLen,
-            collapseEdge
-        );
-        reduce(nSmallCollapsed, sumOp<label>());
-
-        Info<< indent << "Collapsing " << nSmallCollapsed
-            << " small edges" << endl;
-
-
-        // Inline edges
-        label nMerged = mergeEdges
-        (
-            mesh,
-            maxCos,
-            boundaryPoint,
-            minEdgeLen,
-            collapseEdge
-        );
-
-        reduce(nMerged, sumOp<label>());
-        Info<< indent << "Collapsing " << nMerged << " in line edges"
-            << endl;
-
-
-        // Merge edge collapses into consistent collapse-network. Make sure
-        // no cells get collapsed.
-        List<pointEdgeCollapse> allPointInfo;
-        label nLocalCollapse = consistentCollapse
-        (
-            allowCellCollapse,
-            mesh,
-            globalStrings,
-            collapseEdge,
-            allPointInfo
-        );
-
-        reduce(nLocalCollapse, sumOp<label>());
-        Info<< "nLocalCollapse = " << nLocalCollapse << endl;
-
-        if (nLocalCollapse == 0)
+        // Filter small edges only. This reduces the number of faces so that
+        // the face filtering is sped up.
+        label nBadFaces = meshFilter.filterEdges(0);
         {
-            break;
+            polyTopoChange meshMod(newMesh);
+
+            meshMod.changeMesh(mesh, false);
         }
 
-
-        // There are collapses so mesh will get changed
-        meshChanged = true;
-
-
-        // Apply collapses to current mesh
-        polyTopoChange meshMod(mesh);
-
-        // Insert mesh refinement into polyTopoChange.
-        setRefinement(mesh, meshMod, allPointInfo);
-
-        // Do all changes
-        Info<< indent << "Applying changes to the mesh" << nl
-            //<< decrIndent
-            << endl;
-
-        savedMeshMod = meshMod;
-
-        autoPtr<fvMesh> newMeshPtr;
-        autoPtr<mapPolyMesh> mapPtr = meshMod.makeMesh
-        (
-            newMeshPtr,
-            IOobject
-            (
-                mesh.name(),
-                mesh.instance(),
-                mesh.time(),  // register with runTime
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            mesh,
-            true        // parallel sync
-        );
-
-        fvMesh& newMesh = newMeshPtr();
-        const mapPolyMesh& map = mapPtr();
-
-        // Update fields
-        newMesh.updateMesh(map);
-        if (map.hasMotionPoints())
+        if (collapseFaces)
         {
-            newMesh.movePoints(map.preMotionPoints());
+            // 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);
+            {
+                polyTopoChange meshMod(newMesh);
+
+                meshMod.changeMesh(mesh, false);
+            }
         }
 
-
-
-
-        // If no checks needed exit.
-        if (!checkQuality)
-        {
-            break;
-        }
-
-        // Mesh check
-        // ~~~~~~~~~~~~~~~~~~
-        // Do not allow collapses in regions of error.
-        // Updates minEdgeLen, nFrozenEdges
-        label nOldFrozenEdges = returnReduce(nFrozenEdges, sumOp<label>());
-        checkMeshAndFreezeEdges
-        (
-            newMesh,
-            map.reversePointMap(),
-            mesh,
-            minEdgeLen,
-            nFrozenEdges
-        );
-        label nNewFrozenEdges = returnReduce(nFrozenEdges, sumOp<label>());
-
-        Info<< "Frozen edges "
-            << returnReduce(nFrozenEdges, sumOp<label>())
-            << " out of " << returnReduce(mesh.nEdges(), sumOp<label>())
-            << endl;
-
-        if (nNewFrozenEdges == nOldFrozenEdges)
-        {
-            break;
-        }
-
-    } while (true);
-
-
-    if (meshChanged)
-    {
-        // Apply changes to current mesh
-        autoPtr<mapPolyMesh> mapPtr = savedMeshMod.changeMesh(mesh, false);
-        const mapPolyMesh& map = mapPtr();
-
-        // Update fields
-        mesh.updateMesh(map);
-        if (map.hasMotionPoints())
-        {
-            mesh.movePoints(map.preMotionPoints());
-        }
-
-
         // Write resulting mesh
         if (!overwrite)
         {
@@ -1462,7 +133,6 @@ int main(int argc, char *argv[])
         mesh.write();
     }
 
-
     Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
         << "  ClockTime = " << runTime.elapsedClockTime() << " s"
         << nl << endl;
diff --git a/applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C b/applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C
index adb8bdcd5b..a7db08ba71 100644
--- a/applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C
+++ b/applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C
@@ -55,6 +55,7 @@ Description
 #include "edgeCollapser.H"
 #include "meshTools.H"
 #include "Pair.H"
+#include "globalIndex.H"
 
 using namespace Foam;
 
@@ -569,26 +570,47 @@ int main(int argc, char *argv[])
         // Mesh change engine
         edgeCollapser cutter(mesh);
 
-        pointField newPoints(mesh.points());
+        const edgeList& edges = mesh.edges();
+        const pointField& points = mesh.points();
+
+        pointField newPoints(points);
+
+        PackedBoolList collapseEdge(mesh.nEdges());
+        Map<point> collapsePointToLocation(mesh.nPoints());
 
         // Get new positions and construct collapse network
         forAllConstIter(Map<point>, edgeToPos, iter)
         {
             label edgeI = iter.key();
-            const edge& e = mesh.edges()[edgeI];
+            const edge& e = edges[edgeI];
+
+            collapseEdge[edgeI] = true;
+            collapsePointToLocation.set(e[1], points[e[0]]);
 
-            cutter.collapseEdge(edgeI, e[0]);
             newPoints[e[0]] = iter();
         }
 
         // Move master point to destination.
         mesh.movePoints(newPoints);
 
+        List<pointEdgeCollapse> allPointInfo;
+        const globalIndex globalPoints(mesh.nPoints());
+        labelList pointPriority(mesh.nPoints(), 0);
+
+        cutter.consistentCollapse
+        (
+            globalPoints,
+            pointPriority,
+            collapsePointToLocation,
+            collapseEdge,
+            allPointInfo
+        );
+
         // Topo change container
         polyTopoChange meshMod(mesh);
 
         // Insert
-        cutter.setRefinement(meshMod);
+        cutter.setRefinement(allPointInfo, meshMod);
 
         // Do changes
         autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh, false);
diff --git a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C
index 1f72f1d221..bafe305504 100644
--- a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C
+++ b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C
@@ -473,6 +473,7 @@ int main(int argc, char *argv[])
             )
             {
                 label nbrProcI = patchToNbrProc[patchI];
+
                 word name =
                         "procBoundary"
                       + Foam::name(Pstream::myProcNo())
@@ -755,6 +756,9 @@ int main(int argc, char *argv[])
         const edgeList& edges = mesh.edges();
         const pointField& points = mesh.points();
 
+        PackedBoolList collapseEdge(mesh.nEdges());
+        Map<point> collapsePointToLocation(mesh.nPoints());
+
         forAll(edges, edgeI)
         {
             const edge& e = edges[edgeI];
@@ -766,15 +770,29 @@ int main(int argc, char *argv[])
                 Info<< "Merging edge " << e << " since length " << d
                     << " << " << mergeDim << nl;
 
-                // Collapse edge to e[0]
-                collapser.collapseEdge(edgeI, e[0]);
+                collapseEdge[edgeI] = true;
+                collapsePointToLocation.set(e[1], points[e[0]]);
             }
         }
 
+        List<pointEdgeCollapse> allPointInfo;
+        const globalIndex globalPoints(mesh.nPoints());
+        labelList pointPriority(mesh.nPoints(), 0);
+
+        collapser.consistentCollapse
+        (
+            globalPoints,
+            pointPriority,
+            collapsePointToLocation,
+            collapseEdge,
+            allPointInfo
+        );
+
         // Topo change container
         polyTopoChange meshMod(mesh);
+
         // Put all modifications into meshMod
-        bool anyChange = collapser.setRefinement(meshMod);
+        bool anyChange = collapser.setRefinement(allPointInfo, meshMod);
 
         if (anyChange)
         {
diff --git a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMeshApp.C b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMeshApp.C
index f161523ffd..aa033f26d3 100644
--- a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMeshApp.C
+++ b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMeshApp.C
@@ -43,6 +43,7 @@ Note
 #include "edgeCollapser.H"
 #include "addPatchCellLayer.H"
 #include "patchToPoly2DMesh.H"
+#include "globalIndex.H"
 
 using namespace Foam;
 
@@ -260,6 +261,9 @@ int main(int argc, char *argv[])
         const boundBox& bb = mesh().bounds();
         const scalar mergeDim = 1e-4 * bb.minDim();
 
+        PackedBoolList collapseEdge(mesh().nEdges());
+        Map<point> collapsePointToLocation(mesh().nPoints());
+
         forAll(edges, edgeI)
         {
             const edge& e = edges[edgeI];
@@ -271,14 +275,27 @@ int main(int argc, char *argv[])
                 Info<< "Merging edge " << e << " since length " << d
                     << " << " << mergeDim << nl;
 
-                // Collapse edge to e[0]
-                collapser.collapseEdge(edgeI, e[0]);
+                collapseEdge[edgeI] = true;
+                collapsePointToLocation.set(e[1], points[e[0]]);
             }
         }
 
+        List<pointEdgeCollapse> allPointInfo;
+        const globalIndex globalPoints(mesh().nPoints());
+        labelList pointPriority(mesh().nPoints(), 0);
+
+        collapser.consistentCollapse
+        (
+            globalPoints,
+            pointPriority,
+            collapsePointToLocation,
+            collapseEdge,
+            allPointInfo
+        );
+
         polyTopoChange meshModCollapse(mesh());
 
-        collapser.setRefinement(meshModCollapse);
+        collapser.setRefinement(allPointInfo, meshModCollapse);
 
         // Create a mesh from topo changes.
         autoPtr<mapPolyMesh> morphMap
diff --git a/src/dynamicMesh/Make/files b/src/dynamicMesh/Make/files
index bb0ec3ddef..a3e8de6892 100644
--- a/src/dynamicMesh/Make/files
+++ b/src/dynamicMesh/Make/files
@@ -25,6 +25,7 @@ polyTopoChange/polyTopoChange/topoAction/topoActions.C
 polyTopoChange/polyTopoChanger/polyTopoChanger.C
 polyTopoChange/polyTopoChange/polyTopoChange.C
 polyTopoChange/polyTopoChange/addPatchCellLayer.C
+polyTopoChange/polyTopoChange/pointEdgeCollapse/pointEdgeCollapse.C
 polyTopoChange/polyTopoChange/edgeCollapser.C
 polyTopoChange/polyTopoChange/faceCollapser.C
 polyTopoChange/polyTopoChange/hexRef8.C
@@ -97,5 +98,6 @@ createShellMesh/createShellMesh.C
 
 extrudePatchMesh/extrudePatchMesh.C
 
+polyMeshFilter/polyMeshFilter.C
 
 LIB = $(FOAM_LIBBIN)/libdynamicMesh
diff --git a/src/dynamicMesh/polyMeshFilter/polyMeshFilter.C b/src/dynamicMesh/polyMeshFilter/polyMeshFilter.C
new file mode 100644
index 0000000000..4512d9e58f
--- /dev/null
+++ b/src/dynamicMesh/polyMeshFilter/polyMeshFilter.C
@@ -0,0 +1,1082 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "polyMeshFilter.H"
+#include "polyMesh.H"
+#include "fvMesh.H"
+#include "unitConversion.H"
+#include "edgeCollapser.H"
+#include "syncTools.H"
+#include "polyTopoChange.H"
+#include "globalIndex.H"
+#include "PackedBoolList.H"
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::fvMesh> Foam::polyMeshFilter::copyMesh(const fvMesh& mesh)
+{
+    polyTopoChange originalMeshToNewMesh(mesh);
+
+    autoPtr<fvMesh> meshCopy;
+    autoPtr<mapPolyMesh> mapPtr = originalMeshToNewMesh.makeMesh
+    (
+        meshCopy,
+        IOobject
+        (
+            mesh.name(),
+            mesh.polyMesh::instance(),
+            mesh.time(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        mesh,
+        true // parallel sync
+    );
+
+    const mapPolyMesh& map = mapPtr();
+
+    // Update fields
+    meshCopy().updateMesh(map);
+    if (map.hasMotionPoints())
+    {
+        meshCopy().movePoints(map.preMotionPoints());
+    }
+
+    return meshCopy;
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::polyMeshFilter::updatePointErrorCount
+(
+    const PackedBoolList& isErrorPoint,
+    const labelList& oldToNewMesh,
+    labelList& pointErrorCount
+) const
+{
+    forAll(mesh_.points(), pI)
+    {
+        if (isErrorPoint[oldToNewMesh[pI]])
+        {
+            pointErrorCount[pI]++;
+        }
+    }
+}
+
+
+void Foam::polyMeshFilter::checkMeshEdgesAndRelaxEdges
+(
+    const polyMesh& newMesh,
+    const labelList& oldToNewMesh,
+    const PackedBoolList& isErrorPoint,
+    const labelList& pointErrorCount
+)
+{
+    const edgeList& edges = mesh_.edges();
+
+    forAll(edges, edgeI)
+    {
+        const edge& e = edges[edgeI];
+        label newStart = oldToNewMesh[e[0]];
+        label newEnd = oldToNewMesh[e[1]];
+
+        if
+        (
+            pointErrorCount[e[0]] >= maxPointErrorCount_
+         || pointErrorCount[e[1]] >= maxPointErrorCount_
+        )
+        {
+            minEdgeLen_[edgeI] = 0;
+        }
+
+        if
+        (
+            (newStart >= 0 && isErrorPoint[newStart])
+         || (newEnd >= 0 && isErrorPoint[newEnd])
+        )
+        {
+            minEdgeLen_[edgeI] *= edgeReductionFactor_;
+        }
+    }
+
+    syncTools::syncEdgeList(mesh_, minEdgeLen_, minEqOp<scalar>(), 0.0);
+
+    for (label smoothIter = 0; smoothIter < maxSmoothIters_; ++smoothIter)
+    {
+        // Smooth minEdgeLen
+        forAll(mesh_.edges(), edgeI)
+        {
+            const edge& e = mesh_.edges()[edgeI];
+
+            scalar sumMinEdgeLen = 0;
+            label nEdges = 0;
+
+            forAll(e, pointI)
+            {
+                const labelList& pEdges = mesh_.pointEdges()[e[pointI]];
+
+                forAll(pEdges, pEdgeI)
+                {
+                    const label pEdge = pEdges[pEdgeI];
+                    sumMinEdgeLen += minEdgeLen_[pEdge];
+                    nEdges++;
+                }
+            }
+
+            minEdgeLen_[edgeI] = min
+            (
+                minEdgeLen_[edgeI],
+                sumMinEdgeLen/nEdges
+            );
+        }
+
+        syncTools::syncEdgeList(mesh_, minEdgeLen_, minEqOp<scalar>(), 0.0);
+    }
+}
+
+
+void Foam::polyMeshFilter::checkMeshFacesAndRelaxEdges
+(
+    const polyMesh& newMesh,
+    const labelList& oldToNewMesh,
+    const PackedBoolList& isErrorPoint,
+    const labelList& pointErrorCount
+)
+{
+    const faceList& faces = mesh_.faces();
+
+    forAll(faces, faceI)
+    {
+        const face& f = faces[faceI];
+
+        forAll(f, fpI)
+        {
+            const label ptIndex = oldToNewMesh[f[fpI]];
+
+            if (pointErrorCount[f[fpI]] >= maxPointErrorCount_)
+            {
+                faceFilterFactor_[faceI] = 0;
+            }
+
+            if (isErrorPoint[ptIndex])
+            {
+                faceFilterFactor_[faceI] *= faceReductionFactor_;
+
+                break;
+            }
+        }
+    }
+
+    syncTools::syncFaceList(mesh_, faceFilterFactor_, minEqOp<scalar>());
+
+    for (label smoothIter = 0; smoothIter < maxSmoothIters_; ++smoothIter)
+    {
+        // Smooth faceFilterFactor
+        forAll(faces, faceI)
+        {
+            const labelList& fEdges = mesh_.faceEdges()[faceI];
+
+            scalar sumFaceFilterFactors = 0;
+            label nFaces = 0;
+
+            // This is important: Only smooth around faces that share an
+            // edge with a bad face
+            bool skipFace = true;
+
+            forAll(fEdges, fEdgeI)
+            {
+                const labelList& eFaces = mesh_.edgeFaces()[fEdges[fEdgeI]];
+
+                forAll(eFaces, eFaceI)
+                {
+                    const label eFace = eFaces[eFaceI];
+
+                    const face& f = faces[eFace];
+
+                    forAll(f, fpI)
+                    {
+                        const label ptIndex = oldToNewMesh[f[fpI]];
+
+                        if (isErrorPoint[ptIndex])
+                        {
+                            skipFace = false;
+                            break;
+                        }
+                    }
+
+                    if (eFace != faceI)
+                    {
+                        sumFaceFilterFactors += faceFilterFactor_[eFace];
+                        nFaces++;
+                    }
+                }
+            }
+
+            if (skipFace)
+            {
+                continue;
+            }
+
+            faceFilterFactor_[faceI] = min
+            (
+                faceFilterFactor_[faceI],
+                sumFaceFilterFactors/nFaces
+            );
+        }
+
+        // Face filter factor needs to be synchronised!
+        syncTools::syncFaceList(mesh_, faceFilterFactor_, minEqOp<scalar>());
+    }
+}
+
+
+Foam::labelList Foam::polyMeshFilter::findBoundaryPoints
+(
+    const polyMesh& mesh//,
+//    labelIOList& boundaryIOPts
+) const
+{
+    const polyBoundaryMesh& bMesh = mesh.boundaryMesh();
+
+    labelList boundaryPoint(mesh.nPoints(), -1);
+
+    // Get all processor boundary points and the processor patch label
+    // that they are on.
+    forAll(bMesh, patchI)
+    {
+        const polyPatch& patch = bMesh[patchI];
+
+        if (!isA<coupledPolyPatch>(patch))
+        {
+            forAll(patch, fI)
+            {
+                const face& f = patch[fI];
+
+                forAll(f, fp)
+                {
+                    boundaryPoint[f[fp]] = 0; //boundaryIOPts[f[fp]];
+                }
+            }
+        }
+    }
+
+    syncTools::syncPointList
+    (
+        mesh,
+        boundaryPoint,
+        maxEqOp<label>(),
+        labelMin
+    );
+
+    return boundaryPoint;
+}
+
+
+void Foam::polyMeshFilter::printScalarFieldStats
+(
+    const string desc,
+    const scalarField& fld
+) const
+{
+    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>())
+        << decrIndent << endl;
+}
+
+
+void Foam::polyMeshFilter::mapOldMeshEdgeFieldToNewMesh
+(
+    const polyMesh& newMesh,
+    const labelList& pointMap,
+    scalarField& newMeshMinEdgeLen
+) const
+{
+    scalarField tmp(newMesh.nEdges());
+
+    const edgeList& newEdges = newMesh.edges();
+
+    forAll(newEdges, newEdgeI)
+    {
+        const edge& newEdge = newEdges[newEdgeI];
+        const label pStart = newEdge.start();
+        const label pEnd = newEdge.end();
+
+        tmp[newEdgeI] = min
+        (
+            newMeshMinEdgeLen[pointMap[pStart]],
+            newMeshMinEdgeLen[pointMap[pEnd]]
+        );
+    }
+
+    newMeshMinEdgeLen.transfer(tmp);
+
+    syncTools::syncEdgeList
+    (
+        newMesh,
+        newMeshMinEdgeLen,
+        maxEqOp<scalar>(),
+        0.0
+    );
+}
+
+
+void Foam::polyMeshFilter::mapOldMeshFaceFieldToNewMesh
+(
+    const polyMesh& newMesh,
+    const labelList& faceMap,
+    scalarField& newMeshFaceFilterFactor
+) const
+{
+    scalarField tmp(newMesh.nFaces());
+
+    forAll(faceMap, newFaceI)
+    {
+        const label oldFaceI = faceMap[newFaceI];
+
+        tmp[newFaceI] = newMeshFaceFilterFactor[oldFaceI];
+    }
+
+    newMeshFaceFilterFactor.transfer(tmp);
+
+    syncTools::syncFaceList
+    (
+        newMesh,
+        newMeshFaceFilterFactor,
+        maxEqOp<scalar>()
+    );
+}
+
+
+void Foam::polyMeshFilter::updateOldToNewPointMap
+(
+    const labelList& currToNew,
+    labelList& origToCurrentPointMap
+) const
+{
+    forAll(origToCurrentPointMap, origPointI)
+    {
+        label oldPointI = origToCurrentPointMap[origPointI];
+
+        if (oldPointI < currToNew.size())
+        {
+            label newPointI = currToNew[oldPointI];
+
+            if (newPointI != -1)
+            {
+                origToCurrentPointMap[origPointI] = newPointI;
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::polyMeshFilter::polyMeshFilter(const fvMesh& mesh)
+:
+    mesh_(mesh),
+    newMeshPtr_(),
+    dict_
+    (
+        IOobject
+        (
+            "collapseDict",
+            mesh.time().system(),
+            mesh.time(),
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    ),
+    collapseEdgesCoeffDict_(dict_.subDict("collapseEdgesCoeffs")),
+    collapseFacesCoeffDict_(dict_.subDict("collapseFacesCoeffs")),
+    meshQualityCoeffDict_(dict_.subDict("meshQualityCoeffs")),
+    minLen_(readScalar(collapseEdgesCoeffDict_.lookup("minimumEdgeLength"))),
+    maxCos_
+    (
+        ::cos
+        (
+            degToRad
+            (
+                readScalar(collapseEdgesCoeffDict_.lookup("maximumMergeAngle"))
+            )
+        )
+    ),
+    edgeReductionFactor_
+    (
+        readScalar(collapseEdgesCoeffDict_.lookup("reductionFactor"))
+    ),
+    maxIterations_
+    (
+        readLabel(meshQualityCoeffDict_.lookup("maximumIterations"))
+    ),
+    maxSmoothIters_
+    (
+        readLabel(meshQualityCoeffDict_.lookup("maximumSmoothingIterations"))
+    ),
+    initialFaceLengthFactor_
+    (
+        readScalar(collapseFacesCoeffDict_.lookup("initialFaceLengthFactor"))
+    ),
+    faceReductionFactor_
+    (
+        readScalar(collapseFacesCoeffDict_.lookup("reductionFactor"))
+    ),
+    maxPointErrorCount_
+    (
+        readLabel(meshQualityCoeffDict_.lookup("maxPointErrorCount"))
+    ),
+    minEdgeLen_(),
+    faceFilterFactor_()
+{
+    Info<< "Merging:" << nl
+        << "    edges with length less than " << minLen_ << " meters" << 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;
+
+    Info<< "Collapse faces with reduction factor = " << faceReductionFactor_
+        << endl;
+
+    Info<< "Selectively disabling wanted collapses until resulting quality"
+        << " satisfies constraints in system/meshQualityDict" << nl
+        << endl;
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::polyMeshFilter::~polyMeshFilter()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::label Foam::polyMeshFilter::filter(const label nOriginalBadFaces)
+{
+    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//,
+    //                    boundaryIOPts
+                );
+
+                edgeCollapser collapser(newMesh, collapseFacesCoeffDict_);
+
+                // Per face collapse status:
+                //    -1 : not collapsed
+                //  >= 0 : index of point in face to collapse to
+                List<Map<point> > faceCollapseToPoints
+                (
+                    newMesh.nFaces(),
+                    Map<point>()
+                );
+
+                {
+                    // 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
+
+        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
+        );
+    }
+
+    return nBadFaces;
+}
+
+
+Foam::label Foam::polyMeshFilter::filterEdges
+(
+    const label nOriginalBadFaces
+)
+{
+    label nBadFaces = labelMax/2;
+    label nPreviousBadFaces = labelMax;
+    label nOuterIterations = 0;
+
+    minEdgeLen_.resize(mesh_.nEdges(), minLen_);
+    faceFilterFactor_.resize(0);
+
+    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
+     && nBadFaces < nPreviousBadFaces
+    )
+    {
+        Info<< nl << "Outer Iteration = " << nOuterIterations++ << nl
+            << endl;
+
+        printScalarFieldStats("Edge Filter Factor", minEdgeLen_);
+
+        nPreviousBadFaces = nBadFaces;
+
+        // Reset the new mesh to the old mesh
+        newMeshPtr_ = copyMesh(mesh_);
+        fvMesh& newMesh = newMeshPtr_();
+
+        scalarField newMeshMinEdgeLen = minEdgeLen_;
+
+        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//,
+//                    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;
+                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
+
+        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
+        );
+    }
+
+    return nBadFaces;
+}
+
+
+const Foam::autoPtr<Foam::fvMesh>& Foam::polyMeshFilter::filteredMesh() const
+{
+    return newMeshPtr_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/dynamicMesh/polyMeshFilter/polyMeshFilter.H b/src/dynamicMesh/polyMeshFilter/polyMeshFilter.H
new file mode 100644
index 0000000000..e00877a8c1
--- /dev/null
+++ b/src/dynamicMesh/polyMeshFilter/polyMeshFilter.H
@@ -0,0 +1,238 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::polyMeshFilter
+
+Description
+
+SourceFiles
+    polyMeshFilter.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef polyMeshFilter_H
+#define polyMeshFilter_H
+
+#include "IOdictionary.H"
+#include "Time.H"
+#include "List.H"
+#include "autoPtr.H"
+#include "scalarField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class polyMesh;
+class fvMesh;
+class PackedBoolList;
+
+/*---------------------------------------------------------------------------*\
+                       Class polyMeshFilter Declaration
+\*---------------------------------------------------------------------------*/
+
+class polyMeshFilter
+{
+    // Private data
+
+        //- Reference to the original mesh
+        const fvMesh& mesh_;
+
+        //- Copy of the original mesh to perform the filtering on
+        autoPtr<fvMesh> newMeshPtr_;
+
+        //- Dictionary containing the coefficient sub-dictionaries
+        const IOdictionary dict_;
+
+        //- Coefficients for collapsing edges
+        const dictionary& collapseEdgesCoeffDict_;
+
+        //- Coefficients for collapsing faces
+        const dictionary& collapseFacesCoeffDict_;
+
+        //- Coefficients for controlling the mesh quality
+        const 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_;
+
+        //-
+        const label maxPointErrorCount_;
+
+
+        //- The minimum edge length for each edge
+        scalarField minEdgeLen_;
+
+        //- The face filter factor for each face
+        scalarField faceFilterFactor_;
+
+
+    // Private Member Functions
+
+        //- Increment pointErrorCount for points attached to a bad face
+        void updatePointErrorCount
+        (
+            const PackedBoolList& isErrorPoint,
+            const labelList& oldToNewMesh,
+            labelList& pointErrorCount
+        ) const;
+
+
+        //- Given the new points that are part of bad faces, and a map from the
+        //  old mesh points to the new mesh points, relax minEdgeLen_
+        void checkMeshEdgesAndRelaxEdges
+        (
+            const polyMesh& newMesh,
+            const labelList& oldToNewMesh,
+            const PackedBoolList& isErrorPoint,
+            const labelList& pointErrorCount
+        );
+
+        //- Given the new points that are part of bad faces, and a map from the
+        //  old mesh points to the new mesh points, relax faceFilterFactor_
+        void checkMeshFacesAndRelaxEdges
+        (
+            const polyMesh& newMesh,
+            const labelList& oldToNewMesh,
+            const PackedBoolList& isErrorPoint,
+            const labelList& pointErrorCount
+        );
+
+        // Mark boundary points
+        // boundaryPoint:
+        // + -1 : point not on boundary
+        // +  0 : point on a real boundary
+        // + >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;
+
+        //- Print min/mean/max data for a field
+        void printScalarFieldStats
+        (
+            const string desc,
+            const scalarField& fld
+        ) const;
+
+        //- Update minEdgeLen_ for the new mesh based upon the movement of the
+        //  old points to the new points
+        void mapOldMeshEdgeFieldToNewMesh
+        (
+            const polyMesh& newMesh,
+            const labelList& pointMap,
+            scalarField& newMeshMinEdgeLen
+        ) const;
+
+        //- Update faceFilterFactor_ for the new mesh based upon the movement
+        //  of the old faces to the new faces
+        void mapOldMeshFaceFieldToNewMesh
+        (
+            const polyMesh& newMesh,
+            const labelList& faceMap,
+            scalarField& newMeshFaceFilterFactor
+        ) const;
+
+        //- Maintain a map of the original mesh points to the latest version of
+        //  the filtered mesh.
+        void updateOldToNewPointMap
+        (
+            const labelList& currToNew,
+            labelList& origToCurrentPointMap
+        ) const;
+
+        //- Disallow default bitwise copy construct
+        polyMeshFilter(const polyMeshFilter&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const polyMeshFilter&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from fvMesh
+        explicit polyMeshFilter(const fvMesh& mesh);
+
+
+    //- Destructor
+    ~polyMeshFilter();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return reference to the filtered mesh. Does not check if the
+            //  mesh has actually been filtered.
+            const autoPtr<fvMesh>& filteredMesh() const;
+
+
+        // Edit
+
+            //- Return a copy of an fvMesh
+            static autoPtr<fvMesh> copyMesh(const fvMesh& mesh);
+
+            //- Filter edges and faces
+            label filter(const label nOriginalBadFaces);
+
+            //- Filter edges only.
+            label filterEdges(const label nOriginalBadFaces);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C
index 9cc4349803..8290308a50 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C
@@ -26,102 +26,1121 @@ License
 #include "edgeCollapser.H"
 #include "polyMesh.H"
 #include "polyTopoChange.H"
-#include "ListOps.H"
 #include "globalMeshData.H"
-#include "OFstream.H"
-#include "meshTools.H"
 #include "syncTools.H"
+#include "PointEdgeWave.H"
+#include "globalIndex.H"
+#include "removePoints.H"
+#include "motionSmoother.H"
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
-Foam::label Foam::edgeCollapser::findIndex
+Foam::label Foam::edgeCollapser::longestEdge
 (
-    const labelList& elems,
-    const label start,
-    const label size,
-    const label val
+    const face& f,
+    const pointField& pts
 )
 {
-    for (label i = start; i < size; i++)
+    const edgeList& eds = f.edges();
+
+    label longestEdgeI = -1;
+    scalar longestEdgeLength = -SMALL;
+
+    forAll(eds, edI)
     {
-        if (elems[i] == val)
+        scalar edgeLength = eds[edI].mag(pts);
+
+        if (edgeLength > longestEdgeLength)
         {
-            return i;
+            longestEdgeI = edI;
+            longestEdgeLength = edgeLength;
         }
     }
-    return -1;
+
+    return longestEdgeI;
+}
+
+
+Foam::HashSet<Foam::label> Foam::edgeCollapser::checkBadFaces
+(
+    const polyMesh& mesh,
+    const dictionary& meshQualityDict
+)
+{
+    labelHashSet badFaces(mesh.nFaces()/100);
+    DynamicList<label> checkFaces(mesh.nFaces());
+
+    const vectorField& fAreas = mesh.faceAreas();
+
+    scalar faceAreaLimit = SMALL;
+
+    forAll(fAreas, fI)
+    {
+        if (mag(fAreas[fI]) > faceAreaLimit)
+        {
+            checkFaces.append(fI);
+        }
+    }
+
+    Info<< endl;
+
+    motionSmoother::checkMesh
+    (
+        false,
+        mesh,
+        meshQualityDict,
+        checkFaces,
+        badFaces
+    );
+
+    return badFaces;
+}
+
+
+Foam::label Foam::edgeCollapser::checkMeshQuality
+(
+    const polyMesh& mesh,
+    const dictionary& meshQualityDict,
+    PackedBoolList& isErrorPoint
+)
+{
+    labelHashSet badFaces = edgeCollapser::checkBadFaces
+    (
+        mesh,
+        meshQualityDict
+    );
+
+    label nBadFaces = returnReduce(badFaces.size(), sumOp<label>());
+
+    forAllConstIter(labelHashSet, badFaces, iter)
+    {
+        const face& f = mesh.faces()[iter.key()];
+
+        forAll(f, pI)
+        {
+            isErrorPoint[f[pI]] = true;
+        }
+    }
+
+    syncTools::syncPointList
+    (
+        mesh,
+        isErrorPoint,
+        orEqOp<unsigned int>(),
+        0
+    );
+
+    return nBadFaces;
 }
 
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-// Changes region of connected set of points
-Foam::label Foam::edgeCollapser::changePointRegion
+Foam::labelList Foam::edgeCollapser::edgesFromPoints
 (
-    const label pointI,
-    const label oldRegion,
-    const label newRegion
-)
+    const label& faceI,
+    const labelList& pointLabels
+) const
 {
-    label nChanged = 0;
+    labelList edgeLabels(pointLabels.size() - 1, -1);
 
-    if (pointRegion_[pointI] == oldRegion)
+    const labelList& faceEdges = mesh_.faceEdges()[faceI];
+    const edgeList& edges = mesh_.edges();
+
+    label count = 0;
+
+    forAll(faceEdges, eI)
     {
-        pointRegion_[pointI] = newRegion;
-        nChanged++;
+        const label edgeI = faceEdges[eI];
+        const edge& e = edges[edgeI];
 
-        // Step to neighbouring point across edges with same region number
+        label match = 0;
 
-        const labelList& pEdges = mesh_.pointEdges()[pointI];
-
-        forAll(pEdges, i)
+        forAll(pointLabels, pI)
         {
-            label otherPointI = mesh_.edges()[pEdges[i]].otherVertex(pointI);
+            if (e[0] == pointLabels[pI])
+            {
+                match++;
+            }
 
-            nChanged += changePointRegion(otherPointI, oldRegion, newRegion);
+            if (e[1] == pointLabels[pI])
+            {
+                match++;
+            }
+        }
+
+        if (match == 2)
+        {
+            // Edge found
+            edgeLabels[count++] = edgeI;
         }
     }
-    return nChanged;
+
+    if (count != edgeLabels.size())
+    {
+        edgeLabels.setSize(count);
+    }
+
+    return edgeLabels;
 }
 
 
-bool Foam::edgeCollapser::pointRemoved(const label pointI) const
+void Foam::edgeCollapser::collapseToEdge
+(
+    const label faceI,
+    const pointField& pts,
+    const labelList& pointPriority,
+    const vector& collapseAxis,
+    const point& fC,
+    const labelList& facePtsNeg,
+    const labelList& facePtsPos,
+    const scalarList& dNeg,
+    const scalarList& dPos,
+    const scalar dShift,
+    PackedBoolList& collapseEdge,
+    Map<point>& collapsePointToLocation
+) const
 {
-    label region = pointRegion_[pointI];
+    const face& f = mesh_.faces()[faceI];
 
-    if (region == -1 || pointRegionMaster_[region] == pointI)
+    // Negative half
+
+    Foam::point collapseToPtA =
+        collapseAxis*(sum(dNeg)/dNeg.size() - dShift) + fC;
+
+    DynamicList<label> faceBoundaryPts(f.size());
+    DynamicList<label> faceFeaturePts(f.size());
+
+    forAll(facePtsNeg, fPtI)
     {
-        return false;
+        if (pointPriority[facePtsNeg[fPtI]] == 1)
+        {
+            faceFeaturePts.append(facePtsNeg[fPtI]);
+        }
+        else if (pointPriority[facePtsNeg[fPtI]] == 0)
+        {
+            faceBoundaryPts.append(facePtsNeg[fPtI]);
+        }
+    }
+
+    if (!faceBoundaryPts.empty() || !faceFeaturePts.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);
+
+            collapseToPtA = bFace.centre(pts);
+        }
+    }
+
+    faceFeaturePts.clear();
+    faceBoundaryPts.clear();
+
+    labelList faceEdgesNeg = edgesFromPoints(faceI, facePtsNeg);
+
+    forAll(faceEdgesNeg, edgeI)
+    {
+        collapseEdge[faceEdgesNeg[edgeI]] = true;
+    }
+
+    forAll(facePtsNeg, pI)
+    {
+        collapsePointToLocation.set(facePtsNeg[pI], collapseToPtA);
+    }
+
+
+    // Positive half
+
+    Foam::point collapseToPtB
+        = collapseAxis*(sum(dPos)/dPos.size() - dShift) + fC;
+
+    forAll(facePtsPos, fPtI)
+    {
+        if (pointPriority[facePtsPos[fPtI]] == 1)
+        {
+            faceFeaturePts.append(facePtsPos[fPtI]);
+        }
+        else if (pointPriority[facePtsPos[fPtI]] == 0)
+        {
+            // 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]);
+        }
+    }
+
+    if (!faceBoundaryPts.empty() || !faceFeaturePts.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);
+
+            collapseToPtB = bFace.centre(pts);
+        }
+    }
+
+    labelList faceEdgesPos = edgesFromPoints(faceI, facePtsPos);
+
+    forAll(faceEdgesPos, edgeI)
+    {
+        collapseEdge[faceEdgesPos[edgeI]] = true;
+    }
+
+    forAll(facePtsPos, pI)
+    {
+        collapsePointToLocation.set(facePtsPos[pI], collapseToPtB);
+    }
+}
+
+
+void Foam::edgeCollapser::collapseToPoint
+(
+    const label& faceI,
+    const pointField& pts,
+    const labelList& pointPriority,
+    const point& fC,
+    const labelList& facePts,
+    PackedBoolList& collapseEdge,
+    Map<point>& collapsePointToLocation
+) const
+{
+    const face& f = mesh_.faces()[faceI];
+
+    Foam::point collapseToPt = fC;
+
+    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)
+    {
+        const label edgeI = faceEdges[eI];
+        collapseEdge[edgeI] = true;
+    }
+
+    forAll(f, pI)
+    {
+        collapsePointToLocation.set(f[pI], collapseToPt);
+    }
+}
+
+
+void Foam::edgeCollapser::faceCollapseAxisAndAspectRatio
+(
+    const face& f,
+    const point& fC,
+    vector& collapseAxis,
+    scalar& aspectRatio
+) const
+{
+    const pointField& pts = mesh_.points();
+
+    tensor J = f.inertia(pts, fC);
+
+    // Find the dominant collapse direction by finding the eigenvector
+    // that corresponds to the normal direction, discarding it.  The
+    // eigenvector corresponding to the smaller of the two remaining
+    // eigenvalues is the dominant axis in a high aspect ratio face.
+
+    scalar magJ = mag(J);
+
+    scalar detJ = SMALL;
+
+    if (magJ > VSMALL)
+    {
+        // Normalise inertia tensor to remove problems with small values
+
+        J /= mag(J);
+        // J /= cmptMax(J);
+        // J /= max(eigenValues(J).x(), SMALL);
+
+        // Calculating determinant, including stabilisation for zero or
+        // small negative values
+
+        detJ = max(det(J), SMALL);
+    }
+
+    if (detJ < 1e-5)
+    {
+        collapseAxis = f.edges()[longestEdge(f, pts)].vec(pts);
+
+        // It is possible that all the points of a face are the same
+        if (magSqr(collapseAxis) > VSMALL)
+        {
+            collapseAxis /= mag(collapseAxis);
+        }
+
+        // Empirical correlation for high aspect ratio faces
+
+        aspectRatio = Foam::sqrt(0.35/detJ);
     }
     else
     {
-        return true;
+        vector eVals = eigenValues(J);
+
+        if (mag(eVals.y() - eVals.x()) < 100*SMALL)
+        {
+            // First two eigenvalues are the same: i.e. a square face
+
+            // Cannot necessarily determine linearly independent
+            // eigenvectors, or any at all, use longest edge direction.
+
+            collapseAxis = f.edges()[longestEdge(f, pts)].vec(pts);
+
+            collapseAxis /= mag(collapseAxis);
+
+            aspectRatio = 1.0;
+        }
+        else
+        {
+            // The maximum eigenvalue (z()) must be the direction of the
+            // normal, as it has the greatest value.  The minimum eigenvalue
+            // is the dominant collapse axis for high aspect ratio faces.
+
+            collapseAxis = eigenVector(J, eVals.x());
+
+            // The inertia calculation describes the mass distribution as a
+            // function of distance squared to the axis, so the square root of
+            // the ratio of face-plane moments gives a good indication of the
+            // aspect ratio.
+
+            aspectRatio = Foam::sqrt(eVals.y()/max(eVals.x(), SMALL));
+        }
     }
 }
 
 
-void Foam::edgeCollapser::filterFace(const label faceI, face& f) const
+Foam::scalarField Foam::edgeCollapser::calcTargetFaceSizes() const
+{
+    scalarField targetFaceSizes(mesh_.nFaces(), -1);
+
+    const scalarField& cellVolumes = mesh_.cellVolumes();
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+    const labelList& cellOwner = mesh_.faceOwner();
+    const labelList& cellNeighbour = mesh_.faceNeighbour();
+
+    const label nBoundaryFaces = mesh_.nFaces() - mesh_.nInternalFaces();
+
+    // Calculate face size from cell volumes for internal faces
+    for (label intFaceI = 0; intFaceI < mesh_.nInternalFaces(); ++intFaceI)
+    {
+        const scalar cellOwnerVol = cellVolumes[cellOwner[intFaceI]];
+        const scalar cellNeighbourVol = cellVolumes[cellNeighbour[intFaceI]];
+
+        scalar targetFaceSizeA = Foam::pow(cellOwnerVol, 1.0/3.0);
+        scalar targetFaceSizeB = Foam::pow(cellNeighbourVol, 1.0/3.0);
+
+        targetFaceSizes[intFaceI] = 0.5*(targetFaceSizeA + targetFaceSizeB);
+    }
+
+    scalarField neiCellVolumes(nBoundaryFaces, -1);
+
+    // Now do boundary faces
+    forAll(patches, patchI)
+    {
+        const polyPatch& patch = patches[patchI];
+
+        label bFaceI = patch.start() - mesh_.nInternalFaces();
+
+        if (patch.coupled())
+        {
+            // Processor boundary face: Need to get the cell volume on the other
+            // processor
+            const labelUList& faceCells = patch.faceCells();
+
+            forAll(faceCells, facei)
+            {
+                neiCellVolumes[bFaceI++] = cellVolumes[faceCells[facei]];
+            }
+        }
+        else
+        {
+            // Normal boundary face: Just use owner cell volume to calculate
+            // the target face size
+            forAll(patch, patchFaceI)
+            {
+                const label extFaceI = patchFaceI + patch.start();
+                const scalar cellOwnerVol = cellVolumes[cellOwner[extFaceI]];
+
+                targetFaceSizes[extFaceI] = Foam::pow(cellOwnerVol, 1.0/3.0);
+            }
+        }
+    }
+
+    syncTools::swapBoundaryFaceList(mesh_, neiCellVolumes);
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& patch = patches[patchI];
+
+        label bFaceI = patch.start() - mesh_.nInternalFaces();
+
+        if (patch.coupled())
+        {
+            forAll(patch, patchFaceI)
+            {
+                const label localFaceI = patchFaceI + patch.start();
+                const scalar cellOwnerVol = cellVolumes[cellOwner[localFaceI]];
+                const scalar cellNeighbourVol = neiCellVolumes[bFaceI++];
+
+                scalar targetFaceSizeA = Foam::pow(cellOwnerVol, 1.0/3.0);
+                scalar targetFaceSizeB = Foam::pow(cellNeighbourVol, 1.0/3.0);
+
+                targetFaceSizes[localFaceI]
+                    = 0.5*(targetFaceSizeA + targetFaceSizeB);
+            }
+        }
+    }
+
+    // Returns a characteristic length, not an area
+    return targetFaceSizes;
+}
+
+
+Foam::edgeCollapser::collapseType Foam::edgeCollapser::collapseFace
+(
+    const labelList& pointPriority,
+    const face& f,
+    const label faceI,
+    const scalar targetFaceSize,
+    PackedBoolList& collapseEdge,
+    Map<point>& collapsePointToLocation,
+    const scalarField& faceFilterFactor
+) const
+{
+    const scalar collapseSizeLimitCoeff = faceFilterFactor[faceI];
+
+    const pointField& pts = mesh_.points();
+
+    labelList facePts(f);
+
+    const Foam::point fC = f.centre(pts);
+
+    const scalar fA = f.mag(pts);
+
+    vector collapseAxis = vector::zero;
+    scalar aspectRatio = 1.0;
+
+    faceCollapseAxisAndAspectRatio(f, fC, collapseAxis, aspectRatio);
+
+    // The signed distance along the collapse axis passing through the
+    // face centre that each vertex projects to.
+
+    scalarField d(f.size());
+
+    forAll(f, fPtI)
+    {
+        const Foam::point& pt = pts[f[fPtI]];
+
+        d[fPtI] = (collapseAxis & (pt - fC));
+    }
+
+    // Sort the projected distances and the corresponding vertex
+    // indices along the collapse axis
+
+    labelList oldToNew;
+
+    sortedOrder(d, oldToNew);
+
+    oldToNew = invert(oldToNew.size(), oldToNew);
+
+    inplaceReorder(oldToNew, d);
+
+    inplaceReorder(oldToNew, facePts);
+
+    // Shift the points so that they are relative to the centre of the
+    // collapse line.
+
+    scalar dShift = -0.5*(d.first() + d.last());
+
+    d += dShift;
+
+    // Form two lists, one for each half of the set of points
+    // projected along the collapse axis.
+
+    // Middle value, index of first entry in the second half
+    label middle = -1;
+
+    forAll(d, dI)
+    {
+        if (d[dI] > 0)
+        {
+            middle = dI;
+
+            break;
+        }
+    }
+
+    if (middle == -1)
+    {
+//        SeriousErrorIn("collapseFace")
+//            << "middle == -1, " << f << " " << d
+//            << endl;//abort(FatalError);
+
+        return noCollapse;
+    }
+
+    // Negative half
+    SubList<scalar> dNeg(d, middle, 0);
+    SubList<label> facePtsNeg(facePts, middle, 0);
+
+    // Positive half
+    SubList<scalar> dPos(d, d.size() - middle, middle);
+    SubList<label> facePtsPos(facePts, d.size() - middle, middle);
+
+    // Defining how close to the midpoint (M) of the projected
+    // vertices line a projected vertex (X) can be before making this
+    // an invalid edge collapse
+    //
+    // X---X-g----------------M----X-----------g----X--X
+    //
+    // Only allow a collapse if all projected vertices are outwith
+    // guardFraction (g) of the distance form the face centre to the
+    // furthest vertex in the considered direction
+
+    if (dNeg.size() == 0 || dPos.size() == 0)
+    {
+        WarningIn
+        (
+            "Foam::conformalVoronoiMesh::collapseFace"
+        )
+            << "All points on one side of face centre, not collapsing."
+            << endl;
+    }
+
+//    Info<< "Face : " << f << nl
+//        << "    Collapse Axis: " << collapseAxis << nl
+//        << "    Aspect Ratio : " << aspectRatio << endl;
+
+    collapseType typeOfCollapse = noCollapse;
+
+    if (magSqr(collapseAxis) < VSMALL)
+    {
+        typeOfCollapse = toPoint;
+    }
+    else if (fA < aspectRatio*sqr(targetFaceSize*collapseSizeLimitCoeff))
+    {
+        if
+        (
+            allowEarlyCollapseToPoint_
+         && (d.last() - d.first())
+          < targetFaceSize
+           *allowEarlyCollapseCoeff_*maxCollapseFaceToPointSideLengthCoeff_
+        )
+        {
+            typeOfCollapse = toPoint;
+        }
+        else if
+        (
+            (dNeg.last() < guardFraction_*dNeg.first())
+         && (dPos.first() > guardFraction_*dPos.last())
+        )
+        {
+            typeOfCollapse = toEdge;
+        }
+        else if
+        (
+            (d.last() - d.first())
+          < targetFaceSize
+           *maxCollapseFaceToPointSideLengthCoeff_
+        )
+        {
+            // If the face can't be collapsed to an edge, and it has a
+            // small enough span, collapse it to a point.
+            typeOfCollapse = toPoint;
+        }
+    }
+
+    if (typeOfCollapse == toPoint)
+    {
+        collapseToPoint
+        (
+            faceI,
+            pts,
+            pointPriority,
+            fC,
+            facePts,
+            collapseEdge,
+            collapsePointToLocation
+        );
+    }
+    else if (typeOfCollapse == toEdge)
+    {
+        collapseToEdge
+        (
+            faceI,
+            pts,
+            pointPriority,
+            collapseAxis,
+            fC,
+            facePtsNeg,
+            facePtsPos,
+            dNeg,
+            dPos,
+            dShift,
+            collapseEdge,
+            collapsePointToLocation
+        );
+    }
+
+    return typeOfCollapse;
+}
+
+
+Foam::label Foam::edgeCollapser::edgeMaster
+(
+    const labelList& pointPriority,
+    const edge& e
+) const
+{
+    label masterPoint = -1;
+
+    label e0 = e.start();
+    label e1 = e.end();
+
+    // 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;
+}
+
+
+void Foam::edgeCollapser::checkBoundaryPointMergeEdges
+(
+    const label pointI,
+    const label otherPointI,
+    const labelList& pointPriority,
+    Map<point>& collapsePointToLocation
+) const
+{
+   const pointField& points = mesh_.points();
+
+   if
+   (
+       pointPriority[pointI] == 0
+    && pointPriority[otherPointI] < 0
+   )
+   {
+       collapsePointToLocation.set
+       (
+           otherPointI,
+           points[pointI]
+       );
+   }
+   else
+   {
+       collapsePointToLocation.set
+       (
+           pointI,
+           points[otherPointI]
+       );
+   }
+}
+
+
+Foam::label Foam::edgeCollapser::breakStringsAtEdges
+(
+    const PackedBoolList& markedEdges,
+    PackedBoolList& collapseEdge,
+    List<pointEdgeCollapse>& allPointInfo
+) const
+{
+    const edgeList& edges = mesh_.edges();
+    const labelListList& pointEdges = mesh_.pointEdges();
+
+    label nUncollapsed = 0;
+
+    forAll(edges, eI)
+    {
+        if (markedEdges[eI])
+        {
+            const edge& e = edges[eI];
+
+            const label startCollapseIndex
+                = allPointInfo[e.start()].collapseIndex();
+
+            if (startCollapseIndex != -1 && startCollapseIndex != -2)
+            {
+                const label endCollapseIndex
+                    = allPointInfo[e.end()].collapseIndex();
+
+                if
+                (
+                    !collapseEdge[eI]
+                 && startCollapseIndex == endCollapseIndex
+                )
+                {
+                    const labelList& ptEdgesStart = pointEdges[e.start()];
+
+                    forAll(ptEdgesStart, ptEdgeI)
+                    {
+                        const label edgeI = ptEdgesStart[ptEdgeI];
+
+                        const label nbrPointI
+                            = edges[edgeI].otherVertex(e.start());
+                        const label nbrIndex
+                            = allPointInfo[nbrPointI].collapseIndex();
+
+                        if
+                        (
+                            collapseEdge[edgeI]
+                         && nbrIndex == startCollapseIndex
+                        )
+                        {
+                            collapseEdge[edgeI] = false;
+                            nUncollapsed++;
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    return nUncollapsed;
+}
+
+
+void Foam::edgeCollapser::determineDuplicatePointsOnFace
+(
+    const face& f,
+    PackedBoolList& markedPoints,
+    labelHashSet& uniqueCollapses,
+    labelHashSet& duplicateCollapses,
+    List<pointEdgeCollapse>& allPointInfo
+) const
+{
+    uniqueCollapses.clear();
+    duplicateCollapses.clear();
+
+    forAll(f, fpI)
+    {
+        label index = allPointInfo[f[fpI]].collapseIndex();
+
+        // Check for consecutive duplicate
+        if (index != allPointInfo[f.prevLabel(fpI)].collapseIndex())
+        {
+            if (!uniqueCollapses.insert(index))
+            {
+                // Failed inserting so duplicate
+                duplicateCollapses.insert(index);
+            }
+        }
+    }
+
+    // Now duplicateCollapses contains duplicate collapse indices.
+    // Convert to points.
+    forAll(f, fpI)
+    {
+        label index = allPointInfo[f[fpI]].collapseIndex();
+        if (duplicateCollapses.found(index))
+        {
+            markedPoints[f[fpI]] = true;
+        }
+    }
+}
+
+
+Foam::label Foam::edgeCollapser::countEdgesOnFace
+(
+    const face& f,
+    List<pointEdgeCollapse>& allPointInfo
+) const
+{
+    label nEdges = 0;
+
+    forAll(f, fpI)
+    {
+        const label pointI = f[fpI];
+        const label newPointI = allPointInfo[pointI].collapseIndex();
+
+        if (newPointI == -2)
+        {
+            nEdges++;
+        }
+        else
+        {
+            const label prevPointI = f[f.fcIndex(fpI)];
+            const label prevNewPointI
+                = allPointInfo[prevPointI].collapseIndex();
+
+            if (newPointI != prevNewPointI)
+            {
+                nEdges++;
+            }
+        }
+    }
+
+    return nEdges;
+}
+
+
+bool Foam::edgeCollapser::isFaceCollapsed
+(
+    const face& f,
+    List<pointEdgeCollapse>& allPointInfo
+) const
+{
+    label nEdges = countEdgesOnFace(f, allPointInfo);
+
+    // Polygons must have 3 or more edges to be valid
+    if (nEdges < 3)
+    {
+        return true;
+    }
+
+    return false;
+}
+
+
+// Create consistent set of collapses.
+//  collapseEdge : per edge:
+//      -1 : do not collapse
+//       0 : collapse to start
+//       1 : collapse to end
+//  Note: collapseEdge has to be parallel consistent (in orientation)
+Foam::label Foam::edgeCollapser::syncCollapse
+(
+    const globalIndex& globalPoints,
+    const labelList& pointPriority,
+    const PackedBoolList& collapseEdge,
+    const Map<point>& collapsePointToLocation,
+    List<pointEdgeCollapse>& allPointInfo
+) const
+{
+    const edgeList& edges = mesh_.edges();
+
+    label nCollapsed = 0;
+
+    DynamicList<label> initPoints(mesh_.nPoints());
+    DynamicList<pointEdgeCollapse> initPointInfo(mesh_.nPoints());
+
+    allPointInfo.clear();
+    allPointInfo.setSize(mesh_.nPoints());
+
+    // Initialise edges to no collapse
+    List<pointEdgeCollapse> allEdgeInfo
+    (
+        mesh_.nEdges(),
+        pointEdgeCollapse(vector::zero, -1, -1)
+    );
+
+    // Mark selected edges for collapse
+    forAll(edges, edgeI)
+    {
+        if (collapseEdge[edgeI])
+        {
+            const edge& e = edges[edgeI];
+
+            label masterPointI = e.start();
+
+            // Choose the point on the edge with the highest priority.
+            if (pointPriority[e.end()] > pointPriority[e.start()])
+            {
+                masterPointI = e.end();
+            }
+
+            label masterPointPriority = pointPriority[masterPointI];
+
+            label index = globalPoints.toGlobal(masterPointI);
+
+            if (!collapsePointToLocation.found(masterPointI))
+            {
+                const label otherVertex = e.otherVertex(masterPointI);
+
+                if (!collapsePointToLocation.found(otherVertex))
+                {
+                    FatalErrorIn
+                    (
+                        "syncCollapse\n"
+                        "(\n"
+                        "   const polyMesh&,\n"
+                        "   const globalIndex&,\n"
+                        "   const labelList&,\n"
+                        "   const PackedBoolList&,\n"
+                        "   Map<point>&,\n"
+                        "   List<pointEdgeCollapse>&\n"
+                        ")\n"
+                    )   << masterPointI << " on edge " << edgeI << " " << e
+                        << " is not marked for collapse."
+                        << abort(FatalError);
+                }
+                else
+                {
+                    masterPointI = otherVertex;
+                    masterPointPriority = pointPriority[masterPointI];
+                    index = globalPoints.toGlobal(masterPointI);
+                }
+            }
+
+            const point& collapsePoint = collapsePointToLocation[masterPointI];
+
+            const pointEdgeCollapse pec
+            (
+                collapsePoint,
+                index,
+                masterPointPriority
+            );
+
+            // Mark as collapsable but with nonsense master so it gets
+            // overwritten and starts an update wave
+            allEdgeInfo[edgeI] = pointEdgeCollapse
+            (
+                collapsePoint,
+                labelMax,
+                labelMin
+            );
+
+            initPointInfo.append(pec);
+            initPoints.append(e.start());
+
+            initPointInfo.append(pec);
+            initPoints.append(e.end());
+
+            nCollapsed++;
+        }
+    }
+
+    PointEdgeWave<pointEdgeCollapse> collapsePropagator
+    (
+        mesh_,
+        initPoints,
+        initPointInfo,
+        allPointInfo,
+        allEdgeInfo,
+        mesh_.globalData().nTotalPoints()  // Maximum number of iterations
+    );
+
+    return nCollapsed;
+}
+
+
+void Foam::edgeCollapser::filterFace
+(
+    const Map<DynamicList<label> >& collapseStrings,
+    const List<pointEdgeCollapse>& allPointInfo,
+    face& f
+) const
 {
     label newFp = 0;
 
+    face oldFace = f;
+
     forAll(f, fp)
     {
         label pointI = f[fp];
 
-        label region = pointRegion_[pointI];
+        label collapseIndex = allPointInfo[pointI].collapseIndex();
 
-        if (region == -1)
+        // Do we have a local point for this index?
+        if (collapseStrings.found(collapseIndex))
         {
-            f[newFp++] = pointI;
+            label localPointI = collapseStrings[collapseIndex][0];
+
+            if (findIndex(SubList<label>(f, newFp), localPointI) == -1)
+            {
+                f[newFp++] = localPointI;
+            }
+        }
+        else if (collapseIndex == -1)
+        {
+            WarningIn
+            (
+                "filterFace"
+                "(const label, const Map<DynamicList<label> >&, face&)"
+            )   << "Point " << pointI << " was not visited by PointEdgeWave"
+                << endl;
         }
         else
         {
-            label master = pointRegionMaster_[region];
-
-            if (findIndex(f, 0, newFp, master) == -1)
-            {
-                f[newFp++] = master;
-            }
+            f[newFp++] = pointI;
         }
     }
 
@@ -144,7 +1163,7 @@ void Foam::edgeCollapser::filterFace(const label faceI, face& f) const
         label pointI = f[fp];
 
         // Search for previous occurrence.
-        label index = findIndex(f, 0, fp, pointI);
+        label index = findIndex(SubList<label>(f, fp), pointI);
 
         if (index == fp1)
         {
@@ -186,281 +1205,99 @@ void Foam::edgeCollapser::filterFace(const label faceI, face& f) const
 }
 
 
-// Debugging.
-void Foam::edgeCollapser::printRegions() const
-{
-    forAll(pointRegionMaster_, regionI)
-    {
-        label master = pointRegionMaster_[regionI];
-
-        if (master != -1)
-        {
-            Pout<< "Region:" << regionI << nl
-                << "    master:" << master
-                << ' ' << pointRegionMasterLocation_[regionI] << nl;
-
-            forAll(pointRegion_, pointI)
-            {
-                if (pointRegion_[pointI] == regionI && pointI != master)
-                {
-                    Pout<< "    slave:" << pointI
-                        << ' ' <<  mesh_.points()[pointI] << nl;
-                }
-            }
-        }
-    }
-}
-
-
-
-// Collapse list of edges
-void Foam::edgeCollapser::collapseEdges(const labelList& edgeLabels)
-{
-    const edgeList& edges = mesh_.edges();
-
-    forAll(edgeLabels, i)
-    {
-        label edgeI = edgeLabels[i];
-        const edge& e = edges[edgeI];
-
-        label region0 = pointRegion_[e[0]];
-        label region1 = pointRegion_[e[1]];
-
-        if (region0 == -1)
-        {
-            if (region1 == -1)
-            {
-                // Both unaffected. Choose ad lib.
-                collapseEdge(edgeI, e[0]);
-            }
-            else
-            {
-                // Collapse to whatever e[1] collapses
-                collapseEdge(edgeI, e[1]);
-            }
-        }
-        else
-        {
-            if (region1 == -1)
-            {
-                // Collapse to whatever e[0] collapses
-                collapseEdge(edgeI, e[0]);
-            }
-            else
-            {
-                // Both collapsed.
-                if (pointRegionMaster_[region0] == e[0])
-                {
-                    // e[0] is a master
-                    collapseEdge(edgeI, e[0]);
-                }
-                else if (pointRegionMaster_[region1] == e[1])
-                {
-                    // e[1] is a master
-                    collapseEdge(edgeI, e[1]);
-                }
-                else
-                {
-                    // Dont know
-                    collapseEdge(edgeI, e[0]);
-                }
-            }
-        }
-    }
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from mesh
 Foam::edgeCollapser::edgeCollapser(const polyMesh& mesh)
 :
     mesh_(mesh),
-    pointRegion_(mesh.nPoints(), -1),
-    pointRegionMasterLocation_(mesh.nPoints() / 100),
-    pointRegionMaster_(mesh.nPoints() / 100),
-    freeRegions_()
+    guardFraction_(0),
+    maxCollapseFaceToPointSideLengthCoeff_(0),
+    allowEarlyCollapseToPoint_(false),
+    allowEarlyCollapseCoeff_(0)
+{}
+
+
+Foam::edgeCollapser::edgeCollapser
+(
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    mesh_(mesh),
+    guardFraction_(readScalar(dict.lookup("guardFraction"))),
+    maxCollapseFaceToPointSideLengthCoeff_
+    (
+        readScalar(dict.lookup("maxCollapseFaceToPointSideLengthCoeff"))
+    ),
+    allowEarlyCollapseToPoint_
+    (
+        dict.lookupOrDefault<Switch>("allowEarlyCollapseToPoint", true)
+    ),
+    allowEarlyCollapseCoeff_
+    (
+        readScalar(dict.lookup("allowEarlyCollapseCoeff"))
+    )
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-bool Foam::edgeCollapser::unaffectedEdge(const label edgeI) const
-{
-    const edge& e = mesh_.edges()[edgeI];
-
-    return (pointRegion_[e[0]] == -1) && (pointRegion_[e[1]] == -1);
-}
-
-
-bool Foam::edgeCollapser::collapseEdge(const label edgeI, const label master)
-{
-    const pointField& points = mesh_.points();
-
-    const edge& e = mesh_.edges()[edgeI];
-
-    label pointRegion0 = pointRegion_[e[0]];
-    label pointRegion1 = pointRegion_[e[1]];
-
-    if (pointRegion0 == -1)
-    {
-        if (pointRegion1 == -1)
-        {
-            // Both endpoints not collapsed. Create new region.
-
-            label freeRegion = -1;
-
-            if (freeRegions_.size())
-            {
-                freeRegion = freeRegions_.removeHead();
-
-                if (pointRegionMaster_[freeRegion] != -1)
-                {
-                    FatalErrorIn
-                    ("edgeCollapser::collapseEdge(const label, const label)")
-                        << "Problem : freed region :" << freeRegion
-                        << " has already master "
-                        << pointRegionMaster_[freeRegion]
-                        << abort(FatalError);
-                }
-            }
-            else
-            {
-                // If no region found create one. This is the only place where
-                // new regions are created.
-                freeRegion = pointRegionMaster_.size();
-            }
-
-            pointRegion_[e[0]] = freeRegion;
-            pointRegion_[e[1]] = freeRegion;
-
-            pointRegionMaster_(freeRegion) = master;
-            pointRegionMasterLocation_(freeRegion) = points[master];
-        }
-        else
-        {
-            // e[1] is part of collapse network, e[0] not. Add e0 to e1 region.
-            pointRegion_[e[0]] = pointRegion1;
-
-            if
-            (
-                pointRegionMaster_[pointRegion1] == e[0]
-             || pointRegionMaster_[pointRegion1] == e[1]
-            )
-            {
-                pointRegionMaster_[pointRegion1] = master;
-                pointRegionMasterLocation_[pointRegion1] = points[master];
-            }
-        }
-    }
-    else
-    {
-        if (pointRegion1 == -1)
-        {
-            // e[0] is part of collapse network. Add e1 to e0 region
-            pointRegion_[e[1]] = pointRegion0;
-
-            if
-            (
-                pointRegionMaster_[pointRegion0] == e[0]
-             || pointRegionMaster_[pointRegion0] == e[1]
-            )
-            {
-                pointRegionMaster_[pointRegion0] = master;
-                pointRegionMasterLocation_[pointRegion0] = points[master];
-            }
-        }
-        else if (pointRegion0 != pointRegion1)
-        {
-            // Both part of collapse network. Merge the two regions.
-
-            // Use the smaller region number for the whole network.
-            label minRegion = min(pointRegion0, pointRegion1);
-            label maxRegion = max(pointRegion0, pointRegion1);
-
-            // Use minRegion as region for combined net, free maxRegion.
-            pointRegionMaster_[minRegion] = master;
-            pointRegionMaster_[maxRegion] = -1;
-            pointRegionMasterLocation_[minRegion] = points[master];
-            pointRegionMasterLocation_[maxRegion] = point(0, 0, 0);
-
-            freeRegions_.insert(maxRegion);
-
-            if (minRegion != pointRegion0)
-            {
-                changePointRegion(e[0], pointRegion0, minRegion);
-            }
-            if (minRegion != pointRegion1)
-            {
-                changePointRegion(e[1], pointRegion1, minRegion);
-            }
-        }
-    }
-
-    return true;
-}
-
-
-bool Foam::edgeCollapser::setRefinement(polyTopoChange& meshMod)
+bool Foam::edgeCollapser::setRefinement
+(
+    const List<pointEdgeCollapse>& allPointInfo,
+    polyTopoChange& meshMod
+) const
 {
     const cellList& cells = mesh_.cells();
     const labelList& faceOwner = mesh_.faceOwner();
     const labelList& faceNeighbour = mesh_.faceNeighbour();
     const labelListList& pointFaces = mesh_.pointFaces();
-    const labelListList& cellEdges = mesh_.cellEdges();
     const pointZoneMesh& pointZones = mesh_.pointZones();
 
-
-
     bool meshChanged = false;
 
-    // Synchronise pointRegionMasterLocation_
-    const globalMeshData& globalData = mesh_.globalData();
-    const mapDistribute& map = globalData.globalPointSlavesMap();
-    const indirectPrimitivePatch& coupledPatch = globalData.coupledPatch();
-    const labelList& meshPoints = coupledPatch.meshPoints();
-    const Map<label>& meshPointMap = coupledPatch.meshPointMap();
+    PackedBoolList removedPoints(mesh_.nPoints());
 
-
-    List<point> newPoints = coupledPatch.localPoints();
-
-    for (label pI = 0; pI < coupledPatch.nPoints(); ++pI)
+    // Create strings of edges.
+    // Map from collapseIndex(=global master point) to set of points
+    Map<DynamicList<label> > collapseStrings;
     {
-        const label pointRegionMaster = pointRegion_[meshPoints[pI]];
-
-        if (pointRegionMaster != -1)
+        // 1. Count elements per collapseIndex
+        Map<label> nPerIndex(mesh_.nPoints()/10);
+        forAll(allPointInfo, pointI)
         {
-            newPoints[pI]
-                = pointRegionMasterLocation_[pointRegionMaster];
-        }
-    }
+            label collapseIndex = allPointInfo[pointI].collapseIndex();
 
-    globalData.syncData
-    (
-        newPoints,
-        globalData.globalPointSlaves(),
-        globalData.globalPointTransformedSlaves(),
-        map,
-        minMagSqrEqOp<point>()
-    );
-
-        OFstream str1("newPoints_" + name(Pstream::myProcNo()) + ".obj");
-        forAll(pointRegion_, pI)
-        {
-            if (meshPointMap.found(pI))
+            if (collapseIndex != -1 && collapseIndex != -2)
             {
-                meshTools::writeOBJ(str1, newPoints[meshPointMap[pI]]);
+                Map<label>::iterator fnd = nPerIndex.find(collapseIndex);
+                if (fnd != nPerIndex.end())
+                {
+                    fnd()++;
+                }
+                else
+                {
+                    nPerIndex.insert(collapseIndex, 1);
+                }
             }
         }
 
-    for (label pI = 0; pI < coupledPatch.nPoints(); ++pI)
-    {
-        const label pointRegionMaster = pointRegion_[meshPoints[pI]];
-
-        if (pointRegionMaster != -1)
+        // 2. Size
+        collapseStrings.resize(2*nPerIndex.size());
+        forAllConstIter(Map<label>, nPerIndex, iter)
         {
-            pointRegionMasterLocation_[pointRegionMaster]
-                = newPoints[pI];
+            collapseStrings.insert(iter.key(), DynamicList<label>(iter()));
+        }
+
+        // 3. Fill
+        forAll(allPointInfo, pointI)
+        {
+            const label collapseIndex = allPointInfo[pointI].collapseIndex();
+
+            if (collapseIndex != -1 && collapseIndex != -2)
+            {
+                collapseStrings[collapseIndex].append(pointI);
+            }
         }
     }
 
@@ -470,12 +1307,49 @@ bool Foam::edgeCollapser::setRefinement(polyTopoChange& meshMod)
     // Current cellCollapse status
     boolList cellRemoved(mesh_.nCells(), false);
 
+    label nUnvisited = 0;
+    label nUncollapsed = 0;
+    label nCollapsed = 0;
+
+    forAll(allPointInfo, pI)
+    {
+        const pointEdgeCollapse& pec = allPointInfo[pI];
+
+        if (pec.collapseIndex() == -1)
+        {
+            nUnvisited++;
+        }
+        else if (pec.collapseIndex() == -2)
+        {
+            nUncollapsed++;
+        }
+        else
+        {
+            nCollapsed++;
+        }
+    }
+
+    label nPoints = allPointInfo.size();
+
+    reduce(nPoints, sumOp<label>());
+    reduce(nUnvisited, sumOp<label>());
+    reduce(nUncollapsed, sumOp<label>());
+    reduce(nCollapsed, sumOp<label>());
+
+    Info<< incrIndent;
+    Info<< indent << "Number of points : " << nPoints << nl
+        << indent << "Not visited      : " << nUnvisited << nl
+        << indent << "Not collapsed    : " << nUncollapsed << nl
+        << indent << "Collapsed        : " << nCollapsed << nl
+        << endl;
+    Info<< decrIndent;
+
     do
     {
         // Update face collapse from edge collapses
         forAll(newFaces, faceI)
         {
-            filterFace(faceI, newFaces[faceI]);
+            filterFace(collapseStrings, allPointInfo, newFaces[faceI]);
         }
 
         // Check if faces to be collapsed cause cells to become collapsed.
@@ -504,6 +1378,8 @@ bool Foam::edgeCollapser::setRefinement(polyTopoChange& meshMod)
                                 << " of which too many are marked for removal:"
                                 << endl
                                 << "   ";
+
+
                             forAll(cFaces, j)
                             {
                                 if (newFaces[cFaces[j]].size() < 3)
@@ -516,7 +1392,7 @@ bool Foam::edgeCollapser::setRefinement(polyTopoChange& meshMod)
                             cellRemoved[cellI] = true;
 
                             // Collapse all edges of cell to nothing
-                            collapseEdges(cellEdges[cellI]);
+//                            collapseEdges(cellEdges[cellI]);
 
                             nCellCollapsed++;
 
@@ -527,10 +1403,14 @@ bool Foam::edgeCollapser::setRefinement(polyTopoChange& meshMod)
             }
         }
 
+        reduce(nCellCollapsed, sumOp<label>());
+        Info<< indent << "Collapsing " << nCellCollapsed << " cells" << endl;
+
         if (nCellCollapsed == 0)
         {
             break;
         }
+
     } while (true);
 
 
@@ -541,7 +1421,6 @@ bool Foam::edgeCollapser::setRefinement(polyTopoChange& meshMod)
         // Mark points used.
         boolList usedPoint(mesh_.nPoints(), false);
 
-
         forAll(cellRemoved, cellI)
         {
             if (cellRemoved[cellI])
@@ -550,7 +1429,6 @@ bool Foam::edgeCollapser::setRefinement(polyTopoChange& meshMod)
             }
         }
 
-
         // Remove faces
         forAll(newFaces, faceI)
         {
@@ -577,55 +1455,28 @@ bool Foam::edgeCollapser::setRefinement(polyTopoChange& meshMod)
         // Remove unused vertices that have not been marked for removal already
         forAll(usedPoint, pointI)
         {
-            if (!usedPoint[pointI] && !pointRemoved(pointI))
+            if (!usedPoint[pointI])
             {
+                removedPoints[pointI] = true;
                 meshMod.removePoint(pointI, -1);
                 meshChanged = true;
             }
         }
     }
 
-
-
-    // Remove points.
-    forAll(pointRegion_, pointI)
-    {
-        if (pointRemoved(pointI))
-        {
-            meshMod.removePoint(pointI, -1);
-            meshChanged = true;
-        }
-    }
-
     // Modify the point location of the remaining points
-    forAll(pointRegion_, pointI)
+    forAll(allPointInfo, pointI)
     {
-        const label pointRegion = pointRegion_[pointI];
+        const label collapseIndex = allPointInfo[pointI].collapseIndex();
+        const point& collapsePoint = allPointInfo[pointI].collapsePoint();
 
         if
         (
-            !pointRemoved(pointI)
-         && meshPointMap.found(pointI)
+            removedPoints[pointI] == false
+         && collapseIndex != -1
+         && collapseIndex != -2
         )
         {
-            meshMod.modifyPoint
-            (
-                pointI,
-                newPoints[meshPointMap[pointI]],
-                pointZones.whichZone(pointI),
-                false
-            );
-        }
-        else if
-        (
-            pointRegion != -1
-         && !pointRemoved(pointI)
-         && !meshPointMap.found(pointI)
-        )
-        {
-            const point& collapsePoint
-                = pointRegionMasterLocation_[pointRegion];
-
             meshMod.modifyPoint
             (
                 pointI,
@@ -641,9 +1492,9 @@ bool Foam::edgeCollapser::setRefinement(polyTopoChange& meshMod)
     const faceZoneMesh& faceZones = mesh_.faceZones();
 
     // Renumber faces that use points
-    forAll(pointRegion_, pointI)
+    forAll(allPointInfo, pointI)
     {
-        if (pointRemoved(pointI))
+        if (removedPoints[pointI] == true)
         {
             const labelList& changedFaces = pointFaces[pointI];
 
@@ -692,28 +1543,451 @@ bool Foam::edgeCollapser::setRefinement(polyTopoChange& meshMod)
                         zoneID,
                         zoneFlip
                     );
+
                     meshChanged = true;
                 }
             }
         }
     }
 
-    // Print regions:
-//    printRegions();
-
     return meshChanged;
 }
 
 
-void Foam::edgeCollapser::updateMesh(const mapPolyMesh& map)
+void Foam::edgeCollapser::consistentCollapse
+(
+    const globalIndex& globalPoints,
+    const labelList& pointPriority,
+    const Map<point>& collapsePointToLocation,
+    PackedBoolList& collapseEdge,
+    List<pointEdgeCollapse>& allPointInfo,
+    const bool allowCellCollapse
+) const
 {
-    pointRegion_.setSize(mesh_.nPoints());
-    pointRegion_ = -1;
+    // Make sure we don't collapse cells
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    const faceList faces = mesh_.faces();
+    const edgeList& edges = mesh_.edges();
+    const labelListList& faceEdges = mesh_.faceEdges();
+    const labelListList& pointEdges = mesh_.pointEdges();
+    const cellList& cells = mesh_.cells();
 
-    // Reset count, do not remove underlying storage
-    pointRegionMaster_.clear();
-    pointRegionMasterLocation_.clear();
-    freeRegions_.clear();
+    labelHashSet uniqueCollapses;
+    labelHashSet duplicateCollapses;
+
+    while (true)
+    {
+        label nUncollapsed = 0;
+
+        syncTools::syncEdgeList
+        (
+            mesh_,
+            collapseEdge,
+            minEqOp<unsigned int>(),
+            0
+        );
+
+        // Create consistent set of collapses
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        // Note: requires collapseEdge to be synchronised.
+        syncCollapse
+        (
+            globalPoints,
+            pointPriority,
+            collapseEdge,
+            collapsePointToLocation,
+            allPointInfo
+        );
+
+        // Get collapsed faces
+
+        PackedBoolList isCollapsedFace(mesh_.nFaces());
+        PackedBoolList markedPoints(mesh_.nPoints());
+
+        forAll(faces, faceI)
+        {
+            const face& f = faces[faceI];
+
+            isCollapsedFace[faceI] = isFaceCollapsed(f, allPointInfo);
+
+            if (isCollapsedFace[faceI] < 1)
+            {
+                determineDuplicatePointsOnFace
+                (
+                    f,
+                    markedPoints,
+                    uniqueCollapses,
+                    duplicateCollapses,
+                    allPointInfo
+                );
+            }
+        }
+
+        // Synchronise the marked points
+        syncTools::syncPointList
+        (
+            mesh_,
+            markedPoints,
+            orEqOp<unsigned int>(),
+            0
+        );
+
+        // Mark all edges attached to the point for collapse
+        forAll(markedPoints, pointI)
+        {
+            if (markedPoints[pointI])
+            {
+                const label index = allPointInfo[pointI].collapseIndex();
+
+                const labelList& ptEdges = pointEdges[pointI];
+
+                forAll(ptEdges, ptEdgeI)
+                {
+                    const label edgeI = ptEdges[ptEdgeI];
+                    const label nbrPointI = edges[edgeI].otherVertex(pointI);
+                    const label nbrIndex
+                        = allPointInfo[nbrPointI].collapseIndex();
+
+                    if (collapseEdge[edgeI] && nbrIndex == index)
+                    {
+                        collapseEdge[edgeI] = false;
+                        nUncollapsed++;
+                    }
+                }
+            }
+        }
+
+        PackedBoolList markedEdges(mesh_.nEdges());
+
+        if (!allowCellCollapse)
+        {
+            // Check collapsed cells
+            forAll(cells, cellI)
+            {
+                const cell& cFaces = cells[cellI];
+
+                label nFaces = cFaces.size();
+
+                forAll(cFaces, fI)
+                {
+                    label faceI = cFaces[fI];
+
+                    if (isCollapsedFace[faceI])
+                    {
+                        nFaces--;
+                    }
+                }
+
+                if (nFaces < 4)
+                {
+                    forAll(cFaces, fI)
+                    {
+                        label faceI = cFaces[fI];
+
+                        const labelList& fEdges = faceEdges[faceI];
+
+                        // Unmark this face for collapse
+                        forAll(fEdges, fEdgeI)
+                        {
+                            label edgeI = fEdges[fEdgeI];
+
+                            if (collapseEdge[edgeI])
+                            {
+                                collapseEdge[edgeI] = false;
+                                nUncollapsed++;
+                            }
+
+                            markedEdges[edgeI] = true;
+                        }
+
+                        // Uncollapsed this face.
+                        isCollapsedFace[faceI] = false;
+                        nFaces++;
+                    }
+                }
+
+                if (nFaces < 4)
+                {
+                    FatalErrorIn
+                    (
+                        "consistentCollapse\n"
+                        "(\n"
+                        "   labelList&,\n"
+                        "   Map<point>&,\n"
+                        "   bool&,\n"
+                        ")"
+                    )   << "Cell " << cellI << " " << cFaces << nl
+                        << "is " << nFaces << ", "
+                        << "but cell collapse has been disabled."
+                        << abort(FatalError);
+                }
+            }
+        }
+
+        syncTools::syncEdgeList
+        (
+            mesh_,
+            markedEdges,
+            orEqOp<unsigned int>(),
+            0
+        );
+
+        nUncollapsed += breakStringsAtEdges
+        (
+            markedEdges,
+            collapseEdge,
+            allPointInfo
+        );
+
+        reduce(nUncollapsed, sumOp<label>());
+
+        Info<< "            Uncollapsed edges = " << nUncollapsed << " / "
+            << returnReduce(mesh_.nEdges(), sumOp<label>()) << endl;
+
+        if (nUncollapsed == 0)
+        {
+            break;
+        }
+    }
+}
+
+
+Foam::label Foam::edgeCollapser::markSmallEdges
+(
+    const scalarField& minEdgeLen,
+    const labelList& pointPriority,
+    PackedBoolList& collapseEdge,
+    Map<point>& collapsePointToLocation
+) const
+{
+    // Work out which edges to collapse
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    const pointField& points = mesh_.points();
+    const edgeList& edges = mesh_.edges();
+
+    label nCollapsed = 0;
+
+    forAll(edges, edgeI)
+    {
+        const edge& e = edges[edgeI];
+
+        if (!collapseEdge[edgeI])
+        {
+            if (e.mag(points) < minEdgeLen[edgeI])
+            {
+                collapseEdge[edgeI] = true;
+
+                label masterPointI = edgeMaster(pointPriority, e);
+
+                if (masterPointI == -1)
+                {
+                    const point average
+                        = 0.5*(points[e.start()] + points[e.end()]);
+
+                    collapsePointToLocation.set(e.start(), average);
+                }
+                else
+                {
+                    const point& collapsePt = points[masterPointI];
+
+                    collapsePointToLocation.set(masterPointI, collapsePt);
+                }
+
+
+                nCollapsed++;
+            }
+        }
+    }
+
+    return nCollapsed;
+}
+
+
+Foam::label Foam::edgeCollapser::markMergeEdges
+(
+    const scalar maxCos,
+    const labelList& pointPriority,
+    PackedBoolList& collapseEdge,
+    Map<point>& collapsePointToLocation
+) const
+{
+    const edgeList& edges = mesh_.edges();
+    const pointField& points = mesh_.points();
+    const labelListList& pointEdges = mesh_.pointEdges();
+
+    // Point removal engine
+    removePoints pointRemover(mesh_, false);
+
+    // Find out points that can be deleted
+    boolList pointCanBeDeleted;
+    label nTotRemove = pointRemover.countPointUsage(maxCos, pointCanBeDeleted);
+
+    // Rework point-to-remove into edge-to-collapse.
+
+    label nCollapsed = 0;
+
+    if (nTotRemove > 0)
+    {
+        forAll(pointEdges, pointI)
+        {
+            if (pointCanBeDeleted[pointI])
+            {
+                const labelList& pEdges = pointEdges[pointI];
+
+                if (pEdges.size() == 2)
+                {
+                    // Always the case?
+
+                    label e0 = pEdges[0];
+                    label e1 = pEdges[1];
+
+                    if (!collapseEdge[e0] && !collapseEdge[e1])
+                    {
+                        // Get lengths of both edges and choose the smallest
+                        scalar e0length = mag
+                        (
+                            points[edges[e0][0]] - points[edges[e0][1]]
+                        );
+
+                        scalar e1length = mag
+                        (
+                            points[edges[e1][0]] - points[edges[e1][1]]
+                        );
+
+                        if (e0length <= e1length)
+                        {
+                            if (edges[e0].start() == pointI)
+                            {
+                                collapseEdge[e0] = true;
+
+                                checkBoundaryPointMergeEdges
+                                (
+                                    pointI,
+                                    edges[e0].end(),
+                                    pointPriority,
+                                    collapsePointToLocation
+                                );
+                            }
+                            else
+                            {
+                                collapseEdge[e0] = true;
+
+                                checkBoundaryPointMergeEdges
+                                (
+                                    pointI,
+                                    edges[e0].start(),
+                                    pointPriority,
+                                    collapsePointToLocation
+                                );
+                            }
+                        }
+                        else
+                        {
+                            if (edges[e1].start() == pointI)
+                            {
+                                collapseEdge[e1] = true;
+
+                                checkBoundaryPointMergeEdges
+                                (
+                                    pointI,
+                                    edges[e1].end(),
+                                    pointPriority,
+                                    collapsePointToLocation
+                                );
+                            }
+                            else
+                            {
+                                collapseEdge[e1] = true;
+
+                                checkBoundaryPointMergeEdges
+                                (
+                                    pointI,
+                                    edges[e1].start(),
+                                    pointPriority,
+                                    collapsePointToLocation
+                                );
+                            }
+                        }
+
+//                        if (boundaryPoint[leftV] == -1)
+//                        {
+//                            collapseEdge[e0] = findIndex(edges[e0], leftV);
+//                        }
+//                        else
+//                        {
+//                            collapseEdge[e1] = findIndex(edges[e1], rightV);
+//                        }
+
+                        nCollapsed++;
+                    }
+                }
+            }
+        }
+    }
+
+    return nCollapsed;
+}
+
+
+Foam::labelPair Foam::edgeCollapser::markSmallSliverFaces
+(
+    const scalarField& faceFilterFactor,
+    const labelList& pointPriority,
+    PackedBoolList& collapseEdge,
+    Map<point>& collapsePointToLocation
+) const
+{
+    const faceList& faces = mesh_.faces();
+
+    const scalarField targetFaceSizes = calcTargetFaceSizes();
+
+    // Calculate number of faces that will be collapsed to a point or an edge
+    label nCollapseToPoint = 0;
+    label nCollapseToEdge = 0;
+
+    forAll(faces, fI)
+    {
+        const face& f = faces[fI];
+
+        if (faceFilterFactor[fI] == 0)
+        {
+            continue;
+        }
+
+        collapseType flagCollapseFace = collapseFace
+        (
+            pointPriority,
+            f,
+            fI,
+            targetFaceSizes[fI],
+            collapseEdge,
+            collapsePointToLocation,
+            faceFilterFactor
+        );
+
+        if (flagCollapseFace == noCollapse)
+        {
+            continue;
+        }
+        else if (flagCollapseFace == toPoint)
+        {
+            nCollapseToPoint++;
+        }
+        else if (flagCollapseFace == toEdge)
+        {
+            nCollapseToEdge++;
+        }
+        else
+        {
+            FatalErrorIn("collapseFaces(const polyMesh&, List<labelPair>&)")
+                << "Face is marked to be collapsed to " << flagCollapseFace
+                << ". Currently can only collapse to point/edge."
+                << abort(FatalError);
+        }
+    }
+
+    return labelPair(nCollapseToPoint, nCollapseToEdge);
 }
 
 
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.H b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.H
index 3a10a93197..946ee2a72c 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.H
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,10 +37,15 @@ SourceFiles
 #ifndef edgeCollapser_H
 #define edgeCollapser_H
 
-#include "labelList.H"
+#include "pointEdgeCollapse.H"
 #include "DynamicList.H"
-#include "point.H"
+#include "Field.H"
+#include "pointFieldFwd.H"
+#include "Map.H"
+#include "labelPair.H"
+#include "HashSet.H"
 #include "typeInfo.H"
+#include "Switch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -49,9 +54,11 @@ namespace Foam
 
 // Forward declaration of classes
 class polyMesh;
+class PackedBoolList;
 class polyTopoChange;
+class globalIndex;
 class face;
-class mapPolyMesh;
+class edge;
 
 /*---------------------------------------------------------------------------*\
                            Class edgeCollapser Declaration
@@ -59,59 +66,169 @@ class mapPolyMesh;
 
 class edgeCollapser
 {
+public:
+
+    // The type of collapse of a face
+    enum collapseType
+    {
+        noCollapse = 0,
+        toPoint    = 1,
+        toEdge     = 2
+    };
+
+
+private:
+
     // Private data
 
         //- Reference to mesh
         const polyMesh& mesh_;
 
-        //- For every point -1 or region number
-        labelList pointRegion_;
+        const scalar guardFraction_;
 
-        //- Actual location of the point to collapse to for every region master
-        //  point. This will be forced to be consistent across processors
-        DynamicList<point> pointRegionMasterLocation_;
+        const scalar maxCollapseFaceToPointSideLengthCoeff_;
 
-        //- -1 or master vertex for region number
-        DynamicList<label> pointRegionMaster_;
+        const Switch allowEarlyCollapseToPoint_;
 
-        //- Stack of free region numbers. Corresponds to -1 in pointRegionMaster
-        SLList<label> freeRegions_;
-
-
-    // Static Functions
-
-        //- Find val in list. Search starts at start, continues to size-1.
-        static label findIndex
-        (
-            const labelList&,
-            const label start,
-            const label size,
-            const label val
-        );
+        const scalar allowEarlyCollapseCoeff_;
 
 
     // Private Member Functions
 
-        //- Determine points connected through edgesToRemove_.
-        //  Note: Only routine that uses edgesToRemove!
-        label changePointRegion
+        //- Create an edgeList of edges in faceI which have both their points
+        //  in pointLabels
+        labelList edgesFromPoints
+        (
+            const label& faceI,
+            const labelList& pointLabels
+        ) const;
+
+        //- Collapse a face to an edge, marking the collapsed edges and new
+        //  locations for points that will move as a result of the collapse
+        void collapseToEdge
+        (
+            const label faceI,
+            const pointField& pts,
+            const labelList& pointPriority,
+            const vector& collapseAxis,
+            const point& fC,
+            const labelList& facePtsNeg,
+            const labelList& facePtsPos,
+            const scalarList& dNeg,
+            const scalarList& dPos,
+            const scalar dShift,
+            PackedBoolList& collapseEdge,
+            Map<point>& collapsePointToLocation
+        ) const;
+
+        //- Collapse a face to a point, marking the collapsed edges and new
+        //  locations for points that will move as a result of the collapse
+        void collapseToPoint
+        (
+            const label& faceI,
+            const pointField& pts,
+            const labelList& pointPriority,
+            const point& fC,
+            const labelList& facePts,
+            PackedBoolList& collapseEdge,
+            Map<point>& collapsePointToLocation
+        ) const;
+
+        //- Do an eigenvector analysis of the face to get its collapse axis
+        //  and aspect ratio
+        void faceCollapseAxisAndAspectRatio
+        (
+            const face& f,
+            const point& fC,
+            vector& collapseAxis,
+            scalar& aspectRatio
+        ) const;
+
+        //- Return the target length scale for each face
+        scalarField calcTargetFaceSizes() const;
+
+        //- Decides whether a face should be collapsed (and if so it it is to a
+        //  point or an edge)
+        collapseType collapseFace
+        (
+            const labelList& pointPriority,
+            const face& f,
+            const label faceI,
+            const scalar targetFaceSize,
+            PackedBoolList& collapseEdge,
+            Map<point>& collapsePointToLocation,
+            const scalarField& faceFilterFactor
+        ) const;
+
+        //- Return label of point that has the highest priority. This will be
+        //  the point on the edge that will be collapsed to.
+        label edgeMaster(const labelList& pointPriority, const edge& e) const;
+
+        //- Decides which points in an edge to collapse, based on their priority
+        void checkBoundaryPointMergeEdges
         (
             const label pointI,
-            const label oldRegion,
-            const label newRegion
-        );
+            const label otherPointI,
+            const labelList& pointPriority,
+            Map<point>& collapsePointToLocation
+        ) const;
 
-        //- Whether point is master of region or has been removed
-        bool pointRemoved(const label) const;
+        //- Helper function that breaks strings of collapses if an edge is not
+        //  labelled to collapse, but its points both collapse to the same
+        //  location
+        label breakStringsAtEdges
+        (
+            const PackedBoolList& markedEdges,
+            PackedBoolList& collapseEdge,
+            List<pointEdgeCollapse>& allPointInfo
+        ) const;
 
-        //- Renumber f with new vertices. Removes duplicates.
-        void filterFace(const label faceI, face&) const;
+        //- Prevent face pinching by finding points in a face that will be
+        //  collapsed to the same location, but that are not ordered
+        //  consecutively in the face
+        void determineDuplicatePointsOnFace
+        (
+            const face& f,
+            PackedBoolList& markedPoints,
+            labelHashSet& uniqueCollapses,
+            labelHashSet& duplicateCollapses,
+            List<pointEdgeCollapse>& allPointInfo
+        ) const;
 
-        //- Some debugging printing
-        void printRegions() const;
+        //- Count the number of edges on the face that will exist as a result
+        //  of the collapse
+        label countEdgesOnFace
+        (
+            const face& f,
+            List<pointEdgeCollapse>& allPointInfo
+        ) const;
 
-        //- Collapse list of edges. Tries to find master to collapse to.
-        void collapseEdges(const labelList& edgeLabels);
+        //- Does the face have fewer than 3 edges as a result of the potential
+        //  collapse
+        bool isFaceCollapsed
+        (
+            const face& f,
+            List<pointEdgeCollapse>& allPointInfo
+        ) const;
+
+        //- Given the collapse information, propagates the information using
+        //  PointEdgeWave. Result is a list of new point locations and indices
+        label syncCollapse
+        (
+            const globalIndex& globalPoints,
+            const labelList& boundaryPoint,
+            const PackedBoolList& collapseEdge,
+            const Map<point>& collapsePointToLocation,
+            List<pointEdgeCollapse>& allPointInfo
+        ) const;
+
+        //- Renumber f with new vertices. Removes consecutive duplicates
+        void filterFace
+        (
+            const Map<DynamicList<label> >& collapseStrings,
+            const List<pointEdgeCollapse>& allPointInfo,
+            face& f
+        ) const;
 
         //- Disallow default bitwise copy construct
         edgeCollapser(const edgeCollapser&);
@@ -128,43 +245,92 @@ public:
 
     // Constructors
 
-        //- Construct from mesh.
+        //- Construct from mesh
         edgeCollapser(const polyMesh& mesh);
 
+        //- Construct from mesh and dict
+        edgeCollapser(const polyMesh& mesh, const dictionary& dict);
+
 
     // Member Functions
 
-        // Access
+        // Check
 
-            //- For every point the region it belongs to or -1.
-            const labelList& pointRegion() const
-            {
-                return pointRegion_;
-            }
+            //- Find the longest edge in a face
+            static label longestEdge(const face& f, const pointField& pts);
 
-            //- For every region the master (i.e. the point the region will
-            //  be replaced by)
-            const DynamicList<label>& pointRegionMaster() const
-            {
-                return pointRegionMaster_;
-            }
+            //- Calls motionSmoother::checkMesh and returns a set of bad faces
+            static HashSet<label> checkBadFaces
+            (
+                const polyMesh& mesh,
+                const dictionary& meshQualityDict
+            );
 
-            //- Check that edge is not marked for anything
-            bool unaffectedEdge(const label edgeI) const;
+            // Check mesh and mark points on faces in error
+            // Returns boolList with points in error set
+            static label checkMeshQuality
+            (
+                const polyMesh& mesh,
+                const dictionary& meshQualityDict,
+                PackedBoolList& isErrorPoint
+            );
+
+            //- Ensure that the collapse is parallel consistent and update
+            //  allPointInfo.
+            //  Returns a list of edge collapses that is consistent across
+            //  coupled boundaries and a list of pointEdgeCollapses.
+            void consistentCollapse
+            (
+                const globalIndex& globalPoints,
+                const labelList& pointPriority,
+                const Map<point>& collapsePointToLocation,
+                PackedBoolList& collapseEdge,
+                List<pointEdgeCollapse>& allPointInfo,
+                const bool allowCellCollapse = false
+            ) const;
 
 
-        // Edit
+        // Query
 
-            //- Set edge to collapse and point to collapse it to.
-            //  Return true if collapse is valid.
-            //  (always true at the moment)
-            bool collapseEdge(const label edgeI, const label master);
+            //- Play commands into polyTopoChange to create mesh.
+            //  Return true if anything changed.
+            bool setRefinement
+            (
+                const List<pointEdgeCollapse>& allPointInfo,
+                polyTopoChange& meshMod
+            ) const;
 
-            //- Play commands into polyTopoChange to create mesh. Return true
-            //  if anything changed.
-            bool setRefinement(polyTopoChange&);
+            // Mark (in collapseEdge) any edges to collapse
+            label markSmallEdges
+            (
+                const scalarField& minEdgeLen,
+                const labelList& pointPriority,
+                PackedBoolList& collapseEdge,
+                Map<point>& collapsePointToLocation
+            ) const;
 
-            void updateMesh(const mapPolyMesh&);
+            // Mark (in collapseEdge) any edges to merge
+            label markMergeEdges
+            (
+                const scalar maxCos,
+                const labelList& pointPriority,
+                PackedBoolList& collapseEdge,
+                Map<point>& collapsePointToLocation
+            ) const;
+
+            //- Find small faces and sliver faces in the mesh and mark the
+            //  edges that need to be collapsed in order to remove these faces.
+            //  Also returns a map of new locations for points that will move
+            //  as a result of the collapse.
+            //  Use in conjuctions with edgeCollapser to synchronise the
+            //  collapses and modify the mesh
+            labelPair markSmallSliverFaces
+            (
+                const scalarField& faceFilterFactor,
+                const labelList& pointPriority,
+                PackedBoolList& collapseEdge,
+                Map<point>& collapsePointToLocation
+            ) const;
 };
 
 
diff --git a/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapse.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/pointEdgeCollapse/pointEdgeCollapse.C
similarity index 88%
rename from applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapse.C
rename to src/dynamicMesh/polyTopoChange/polyTopoChange/pointEdgeCollapse/pointEdgeCollapse.C
index 6830e1c8f5..15747ded56 100644
--- a/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapse.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/pointEdgeCollapse/pointEdgeCollapse.C
@@ -34,7 +34,9 @@ Foam::Ostream& Foam::operator<<
 )
 {
     return os
-        << wDist.collapsePoint_ << wDist.collapseIndex_;
+        << wDist.collapsePoint_
+        << wDist.collapseIndex_
+        << wDist.collapsePriority_;
 }
 
 Foam::Istream& Foam::operator>>
@@ -44,7 +46,9 @@ Foam::Istream& Foam::operator>>
 )
 {
     return is
-        >> wDist.collapsePoint_ >> wDist.collapseIndex_;
+        >> wDist.collapsePoint_
+        >> wDist.collapseIndex_
+        >> wDist.collapsePriority_;
 }
 
 
diff --git a/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapse.H b/src/dynamicMesh/polyTopoChange/polyTopoChange/pointEdgeCollapse/pointEdgeCollapse.H
similarity index 96%
rename from applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapse.H
rename to src/dynamicMesh/polyTopoChange/polyTopoChange/pointEdgeCollapse/pointEdgeCollapse.H
index cd8383b8b3..83ac672b33 100644
--- a/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapse.H
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/pointEdgeCollapse/pointEdgeCollapse.H
@@ -62,6 +62,9 @@ class pointEdgeCollapse
         //- Collapse string index
         label collapseIndex_;
 
+        //- Priority of the collapse
+        label collapsePriority_;
+
 
     // Private Member Functions
 
@@ -89,7 +92,8 @@ public:
         inline pointEdgeCollapse
         (
             const point& collapsePoint,
-            const label collapseIndex
+            const label collapseIndex,
+            const label collapsePriority
         );
 
 
@@ -101,6 +105,8 @@ public:
 
             inline label collapseIndex() const;
 
+            inline label collapsePriority() const;
+
 
         // Needed by meshWave
 
diff --git a/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapseI.H b/src/dynamicMesh/polyTopoChange/polyTopoChange/pointEdgeCollapse/pointEdgeCollapseI.H
similarity index 85%
rename from applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapseI.H
rename to src/dynamicMesh/polyTopoChange/polyTopoChange/pointEdgeCollapse/pointEdgeCollapseI.H
index eeb4d6baf2..3b01268a92 100644
--- a/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapseI.H
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/pointEdgeCollapse/pointEdgeCollapseI.H
@@ -43,48 +43,58 @@ inline bool Foam::pointEdgeCollapse::update
             << "problem." << abort(FatalError);
     }
 
+    if (!valid(td))
+    {
+        operator=(w2);
+        return true;
+    }
+
     if (w2.collapseIndex_ == -1 || collapseIndex_ == -1)
     {
         // Not marked for collapse; only happens on edges.
         return false;
     }
 
-    if (!valid(td))
+    if (w2.collapsePriority_ < collapsePriority_)
+    {
+        return false;
+    }
+    else if (w2.collapsePriority_ > collapsePriority_)
     {
         operator=(w2);
         return true;
     }
-    else
-    {
-        // Take over w2 if it is 'better'
 
-        if (w2.collapseIndex_ < collapseIndex_)
+    // Get overwritten by w2 if it has a higher priority
+    if (w2.collapseIndex_ < collapseIndex_)
+    {
+        operator=(w2);
+        return true;
+    }
+    else if (w2.collapseIndex_ == collapseIndex_)
+    {
+        bool identicalPoint = samePoint(w2.collapsePoint_);
+
+        bool nearer = magSqr(w2.collapsePoint_)
+                    < magSqr(collapsePoint_);
+
+        if (nearer)
         {
-            // Take over string index and coordinate from w2
             operator=(w2);
-            return true;
         }
-        else if (w2.collapseIndex_ == collapseIndex_)
-        {
-            bool identicalPoint = samePoint(w2.collapsePoint_);
-            bool nearer = magSqr(w2.collapsePoint_) < magSqr(collapsePoint_);
-            if (nearer)
-            {
-                operator=(w2);
-            }
-            if (identicalPoint)
-            {
-                return false;
-            }
-            else
-            {
-                return nearer;
-            }
-        }
-        else
+
+        if (identicalPoint)
         {
             return false;
         }
+        else
+        {
+            return nearer;
+        }
+    }
+    else
+    {
+        return false;
     }
 }
 
@@ -95,7 +105,8 @@ inline bool Foam::pointEdgeCollapse::update
 inline Foam::pointEdgeCollapse::pointEdgeCollapse()
 :
     collapsePoint_(GREAT, GREAT, GREAT),
-    collapseIndex_(-2)
+    collapseIndex_(-2),
+    collapsePriority_(-2)
 {}
 
 
@@ -103,11 +114,13 @@ inline Foam::pointEdgeCollapse::pointEdgeCollapse()
 inline Foam::pointEdgeCollapse::pointEdgeCollapse
 (
     const point& collapsePoint,
-    const label collapseIndex
+    const label collapseIndex,
+    const label collapsePriority
 )
 :
     collapsePoint_(collapsePoint),
-    collapseIndex_(collapseIndex)
+    collapseIndex_(collapseIndex),
+    collapsePriority_(collapsePriority)
 {}
 
 
@@ -125,6 +138,12 @@ inline Foam::label Foam::pointEdgeCollapse::collapseIndex() const
 }
 
 
+inline Foam::label Foam::pointEdgeCollapse::collapsePriority() const
+{
+    return collapsePriority_;
+}
+
+
 inline bool Foam::pointEdgeCollapse::samePoint(const point& pt) const
 {
     bool isLegal1 = (cmptMin(collapsePoint_) < 0.5*GREAT);
@@ -268,6 +287,7 @@ inline bool Foam::pointEdgeCollapse::operator==
 {
     return
         collapseIndex_ == rhs.collapseIndex_
+     && collapsePriority_ == rhs.collapsePriority_
      && samePoint(rhs.collapsePoint_);
 }