From c0695e04faff3c54386780c0476a7e8c2630baf6 Mon Sep 17 00:00:00 2001
From: laurence <laurence>
Date: Fri, 20 Apr 2012 17:07:10 +0100
Subject: [PATCH] ENH: collapseEdges: Use PointEdgeWave to propagate collapse
 information

---
 .../mesh/advanced/collapseEdges/Make/files    |    1 +
 .../mesh/advanced/collapseEdges/Make/options  |    7 +-
 .../advanced/collapseEdges/collapseEdges.C    | 1293 ++++++++++++-----
 .../pointEdgeCollapse/pointEdgeCollapse.C     |   51 +
 .../pointEdgeCollapse/pointEdgeCollapse.H     |  227 +++
 .../pointEdgeCollapse/pointEdgeCollapseI.H    |  313 ++++
 6 files changed, 1549 insertions(+), 343 deletions(-)
 create mode 100644 applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapse.C
 create mode 100644 applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapse.H
 create mode 100644 applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapseI.H

diff --git a/applications/utilities/mesh/advanced/collapseEdges/Make/files b/applications/utilities/mesh/advanced/collapseEdges/Make/files
index d89ca6e737..a15838abe8 100644
--- a/applications/utilities/mesh/advanced/collapseEdges/Make/files
+++ b/applications/utilities/mesh/advanced/collapseEdges/Make/files
@@ -1,3 +1,4 @@
 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 73f34b0f85..d1efa61fd5 100644
--- a/applications/utilities/mesh/advanced/collapseEdges/Make/options
+++ b/applications/utilities/mesh/advanced/collapseEdges/Make/options
@@ -1,7 +1,10 @@
 EXE_INC = \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -IpointEdgeCollapse
 
 EXE_LIBS = \
     -ldynamicMesh \
-    -lmeshTools
+    -lmeshTools \
+    -lfiniteVolume
diff --git a/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C b/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C
index 400191275c..899159b2a2 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
@@ -41,10 +41,10 @@ Description
 
 #include "argList.H"
 #include "Time.H"
-#include "edgeCollapser.H"
 #include "polyTopoChange.H"
 #include "polyTopoChanger.H"
 #include "polyMesh.H"
+#include "fvMesh.H"
 #include "mapPolyMesh.H"
 #include "mathematicalConstants.H"
 #include "PackedBoolList.H"
@@ -52,6 +52,9 @@ Description
 #include "unitConversion.H"
 #include "globalMeshData.H"
 #include "globalIndex.H"
+#include "PointEdgeWave.H"
+#include "pointEdgeCollapse.H"
+#include "motionSmoother.H"
 
 #include "OFstream.H"
 #include "meshTools.H"
@@ -60,6 +63,413 @@ using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+label findIndex
+(
+    const labelList& elems,
+    const label start,
+    const label size,
+    const label val
+)
+{
+    for (label i = start; i < size; i++)
+    {
+        if (elems[i] == val)
+        {
+            return i;
+        }
+    }
+    return -1;
+}
+
+
+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(f, 0, 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(f, 0, 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<DynamicList<label> > collapseStrings;
+
+    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)
+                        {
+                            Info<< "Cell:" << cellI
+                                << " uses faces:" << cFaces
+                                << " of which too many are marked for removal:"
+                                << endl
+                                << "   ";
+                            forAll(cFaces, j)
+                            {
+                                if (newFaces[cFaces[j]].size() < 3)
+                                {
+                                    Info<< ' '<< cFaces[j];
+                                }
+                            }
+                            Info<< 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;
+                }
+            }
+        }
+    }
+
+    // Print regions:
+//    printRegions();
+
+    return meshChanged;
+}
+
+
 // Get faceEdges in order of face points, i.e. faceEdges[0] is between
 // f[0] and f[1]
 labelList getSortedEdges
@@ -98,84 +508,96 @@ labelList getSortedEdges
 
 
 // Merges edges which are in straight line. I.e. edge split by point.
-label mergeEdges
-(
-    const polyMesh& mesh,
-    const scalar maxCos,
-    edgeCollapser& collapser
-)
-{
-    const pointField& points = mesh.points();
-    const edgeList& edges = mesh.edges();
-    const labelListList& pointEdges = mesh.pointEdges();
-    const labelList& region = collapser.pointRegion();
-    const labelList& master = collapser.pointRegionMaster();
-
-    label nCollapsed = 0;
-
-    forAll(pointEdges, pointI)
-    {
-        const labelList& pEdges = pointEdges[pointI];
-
-        if (pEdges.size() == 2)
-        {
-            const edge& leftE = edges[pEdges[0]];
-            const edge& rightE = edges[pEdges[1]];
-
-            // Get the two vertices on both sides of the point
-            label leftV = leftE.otherVertex(pointI);
-            label rightV = rightE.otherVertex(pointI);
-
-            // Collapse only if none of the points part of merge network
-            // or all of networks with different masters.
-            label midMaster = -1;
-            if (region[pointI] != -1)
-            {
-                midMaster = master[region[pointI]];
-            }
-
-            label leftMaster = -2;
-            if (region[leftV] != -1)
-            {
-                leftMaster = master[region[leftV]];
-            }
-
-            label rightMaster = -3;
-            if (region[rightV] != -1)
-            {
-                rightMaster = master[region[rightV]];
-            }
-
-            if
-            (
-                midMaster != leftMaster
-             && midMaster != rightMaster
-             && leftMaster != rightMaster
-            )
-            {
-                // Check if the two edge are in line
-                vector leftVec = points[pointI] - points[leftV];
-                leftVec /= mag(leftVec) + VSMALL;
-
-                vector rightVec = points[rightV] - points[pointI];
-                rightVec /= mag(rightVec) + VSMALL;
-
-                if ((leftVec & rightVec) > maxCos)
-                {
-                    // Collapse one (left) side of the edge. Make left vertex
-                    // the master.
-                    //if (collapser.unaffectedEdge(pEdges[0]))
-                    {
-                        collapser.collapseEdge(pEdges[0], leftV);
-                        nCollapsed++;
-                    }
-                }
-            }
-        }
-    }
-
-    return nCollapsed;
-}
+//label mergeEdges
+//(
+//    const polyMesh& mesh,
+//    const scalar maxCos,
+//    List<pointEdgeCollapse>& allPointInfo
+//)
+//{
+//    const pointField& points = mesh.points();
+//    const edgeList& edges = mesh.edges();
+//    const labelListList& pointEdges = mesh.pointEdges();
+//    const labelList& region = collapser.pointRegion();
+//    const labelList& master = collapser.pointRegionMaster();
+//
+//    label nCollapsed = 0;
+//
+//    forAll(pointEdges, pointI)
+//    {
+//        const labelList& pEdges = pointEdges[pointI];
+//
+//        if (pEdges.size() == 2)
+//        {
+//            const edge& leftE = edges[pEdges[0]];
+//            const edge& rightE = edges[pEdges[1]];
+//
+//            // Get the two vertices on both sides of the point
+//            label leftV = leftE.otherVertex(pointI);
+//            label rightV = rightE.otherVertex(pointI);
+//
+//            // Collapse only if none of the points part of merge network
+//            // or all of networks with different masters.
+//            label midMaster = -1;
+//            if (region[pointI] != -1)
+//            {
+//                midMaster = master[region[pointI]];
+//            }
+//
+//            label leftMaster = -2;
+//            if (region[leftV] != -1)
+//            {
+//                leftMaster = master[region[leftV]];
+//            }
+//
+//            label rightMaster = -3;
+//            if (region[rightV] != -1)
+//            {
+//                rightMaster = master[region[rightV]];
+//            }
+//
+//            if
+//            (
+//                midMaster != leftMaster
+//             && midMaster != rightMaster
+//             && leftMaster != rightMaster
+//            )
+//            {
+//                // Check if the two edge are in line
+//                vector leftVec = points[pointI] - points[leftV];
+//                leftVec /= mag(leftVec) + VSMALL;
+//
+//                vector rightVec = points[rightV] - points[pointI];
+//                rightVec /= mag(rightVec) + VSMALL;
+//
+//                if ((leftVec & rightVec) > maxCos)
+//                {
+//                    // Collapse one (left) side of the edge. Make left vertex
+//                    // the master.
+//                    //if (collapser.unaffectedEdge(pEdges[0]))
+//                    const edge& e = mesh.edges()[pEdges[0]];
+//
+//                    if
+//                    (
+//                        allPointInfo[e[0]].collapseIndex() < 0
+//                     && allPointInfo[e[1]].collapseIndex() < 0
+//                    )
+//                    {
+//                        //pointEdgeCollapse pec(mesh.points()[leftV], leftV);
+//
+//                        allPointInfo[e[0]] = pec;
+//                        allPointInfo[e[1]] = pec;
+//
+//                        //collapser.collapseEdge(pEdges[0], leftV);
+//                        nCollapsed++;
+//                    }
+//                }
+//            }
+//        }
+//    }
+//
+//    return nCollapsed;
+//}
 
 
 // Return master point edge needs to be collapsed to (or -1)
@@ -198,26 +620,26 @@ label edgeMaster
     }
 
     // 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;
-    }
+//    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)
@@ -254,23 +676,15 @@ label edgeMaster
 label collapseSmallEdges
 (
     const polyMesh& mesh,
+    const scalarList& freezeEdges,
     const labelList& boundaryPoint,
     const scalar minLen,
-    edgeCollapser& collapser
+    List<pointEdgeCollapse>& allPointInfo
 )
 {
     const pointField& points = mesh.points();
     const edgeList& edges = mesh.edges();
 
-    // Collapse all edges that are too small. Choose intelligently which
-    // point to collapse edge to.
-
-    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();
-
     // Store collapse direction in collapseEdge
     //   -1 -> Do not collapse
     //    0 -> Collapse to start point
@@ -281,97 +695,111 @@ label collapseSmallEdges
     {
         const edge& e = edges[edgeI];
 
-        if (e.mag(points) < minLen)
+        if (e.mag(points) < minLen*freezeEdges[edgeI])
         {
             collapseEdge[edgeI] = 0;
         }
     }
 
     // Check whether edge point order is reversed from mesh to coupledPatch
-    PackedBoolList meshToPatchSameOrientation(coupledMeshEdges.size(), true);
-
-    forAll(coupledMeshEdges, eI)
-    {
-        const label meshEdgeIndex = coupledMeshEdges[eI];
-
-        if (collapseEdge[meshEdgeIndex] != -1)
-        {
-            const edge& meshEdge = 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(coupledMeshEdges.size(), -1);
-
-    forAll(coupledMeshEdges, eI)
-    {
-        const label meshEdgeIndex = coupledMeshEdges[eI];
-
-        if (collapseEdge[meshEdgeIndex] != -1)
-        {
-            if (meshToPatchSameOrientation[eI] == cppOrientation[eI])
-            {
-                cppEdgeData[eI] = 0;
-            }
-            else
-            {
-                cppEdgeData[eI] = 1;
-            }
-        }
-    }
-
-
-    // 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];
-
-        if (collapseEdge[meshEdgeIndex] != -1)
-        {
-            if (meshToPatchSameOrientation[eI] == cppOrientation[eI])
-            {
-                collapseEdge[meshEdgeIndex] = 0;
-            }
-            else
-            {
-                collapseEdge[meshEdgeIndex] = 1;
-            }
-        }
-    }
-
-
-    OFstream str1("collapsedPoints_" + name(Pstream::myProcNo()) + ".obj");
+//    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 = 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(coupledMeshEdges.size(), -1);
+//
+//    forAll(coupledMeshEdges, eI)
+//    {
+//        const label meshEdgeIndex = coupledMeshEdges[eI];
+//
+//        if (collapseEdge[meshEdgeIndex] != -1)
+//        {
+//            if (meshToPatchSameOrientation[eI] == cppOrientation[eI])
+//            {
+//                cppEdgeData[eI] = 0;
+//            }
+//            else
+//            {
+//                cppEdgeData[eI] = 1;
+//            }
+//        }
+//    }
+//
+//
+//    // 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];
+//
+//        if (collapseEdge[meshEdgeIndex] != -1)
+//        {
+//            if (meshToPatchSameOrientation[eI] == cppOrientation[eI])
+//            {
+//                collapseEdge[meshEdgeIndex] = 0;
+//            }
+//            else
+//            {
+//                collapseEdge[meshEdgeIndex] = 1;
+//            }
+//        }
+//    }
 
     label nCollapsed = 0;
 
+    DynamicList<label> initPoints(mesh.nPoints());
+    DynamicList<pointEdgeCollapse> initPointInfo(mesh.nPoints());
+    allPointInfo.resize(mesh.nPoints());
+
+    globalIndex globalStrings(mesh.nPoints());
+
+    List<pointEdgeCollapse> allEdgeInfo(mesh.nEdges());
+    forAll(allEdgeInfo, edgeI)
+    {
+        allEdgeInfo[edgeI] = pointEdgeCollapse(vector::zero, -1);
+    }
+
     forAll(edges, edgeI)
     {
+        const edge& e = edges[edgeI];
+
         if (collapseEdge[edgeI] != -1)
         {
-            const edge& e = edges[edgeI];
-
             const label master =
                 edgeMaster
                 (
@@ -380,24 +808,37 @@ label collapseSmallEdges
                     e
                 );
 
-            if (e[0] == master)
+//            if (master != -1)
             {
-                meshTools::writeOBJ(str1, points[e[1]], points[e[0]]);
-            }
-            else if (e[1] == master)
-            {
-                meshTools::writeOBJ(str1, points[e[0]], points[e[1]]);
-            }
+                pointEdgeCollapse pec
+                (
+                    points[master],
+                    globalStrings.toGlobal(master)
+                );
 
-            if (master != -1)
-            {
-                collapser.collapseEdge(edgeI, master);
+                allEdgeInfo[edgeI] = pec;
+
+                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;
 }
 
@@ -474,8 +915,13 @@ label collapseSmallEdges
 //
 //            if (edgeI != -1)
 //            {
-//                label master =
-//                    edgeMaster(boundaryPoint, processorPoints, false, edges[edgeI]);
+//                label master = edgeMaster
+//                    (
+//                        boundaryPoint,
+//                        processorPoints,
+//                        false,
+//                        edges[edgeI]
+//                    );
 //
 //                if (master != -1)// && collapser.unaffectedEdge(edgeI))
 //                {
@@ -575,7 +1021,8 @@ void set(const labelList& elems, const bool val, boolList& status)
 //                    {
 //                        const edge& e = edges[edgeI];
 //
-//                        label master = edgeMaster(boundaryPoint, processorPoints, false, e);
+//                        label master
+//                       = edgeMaster(boundaryPoint, processorPoints, false, e);
 //
 //                        if (master != -1)
 //                        {
@@ -588,13 +1035,13 @@ void set(const labelList& elems, const bool val, boolList& status)
 //
 //                            forAll(pCells0, i)
 //                            {
-//                                set(cellEdges[pCells0[i]], true, protectedEdge);
+//                              set(cellEdges[pCells0[i]], true, protectedEdge);
 //                            }
 //                            const labelList& pCells1 = pointCells[e[1]];
 //
 //                            forAll(pCells1, i)
 //                            {
-//                                set(cellEdges[pCells1[i]], true, protectedEdge);
+//                              set(cellEdges[pCells1[i]], true, protectedEdge);
 //                            }
 //
 //                            nCollapsed++;
@@ -609,14 +1056,144 @@ void set(const labelList& elems, const bool val, boolList& status)
 //}
 
 
-// Main program:
+labelHashSet checkMeshQuality(const polyMesh& mesh)
+{
+    //mesh.checkMesh(true);
+    labelHashSet freezePoints;
 
+    IOdictionary meshQualityDict
+    (
+        IOobject
+        (
+            "meshQualityControls",
+            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;
+
+    forAllConstIter(labelHashSet, badFaces, iter)
+    {
+        const face& f = mesh.faces()[iter.key()];
+
+        forAll(f, pI)
+        {
+            freezePoints.insert(f[pI]);
+        }
+    }
+
+//    const edgeList& edges = mesh.edges();
+//
+//    label nFrozenEdges = 0;
+//
+//    OFstream str("frozenEdges.obj");
+//
+//    label count = 0;
+//    forAll(edges, eI)
+//    {
+//        const edge& e = edges[eI];
+//
+//        if (freezePoints.found(e[0]) && freezePoints.found(e[1]))
+//        {
+//            freezeEdges[eI] = true;
+//            nFrozenEdges++;
+//        }
+//    }
+
+    return freezePoints;
+}
+
+
+// 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<processorPolyPatch>(patch))
+        {
+            const processorPolyPatch& pPatch =
+                refCast<const processorPolyPatch>(patch);
+
+            forAll(pPatch, fI)
+            {
+                const face& f = pPatch[fI];
+
+                forAll(f, fp)
+                {
+                    boundaryPoint[f[fp]] = patchI;
+                }
+            }
+        }
+    }
+
+    return boundaryPoint;
+}
+
+
+// Main program:
 int main(int argc, char *argv[])
 {
 #   include "addOverwriteOption.H"
 
     argList::validArgs.append("edge length [m]");
     argList::validArgs.append("merge angle (degrees)");
+    argList::addOption("minLenFactor", "scalar", "edge length factor");
 
 #   include "setRootCase.H"
 #   include "createTime.H"
@@ -624,8 +1201,11 @@ int main(int argc, char *argv[])
 #   include "createPolyMesh.H"
     const word oldInstance = mesh.pointsInstance();
 
-    const scalar minLen  = args.argRead<scalar>(1);
+    scalar minLen  = args.argRead<scalar>(1);
     const scalar angle   = args.argRead<scalar>(2);
+    const scalar minLenFactor
+        = args.optionLookupOrDefault<scalar>("minLenFactor", 0.5);
+
     const bool overwrite = args.optionFound("overwrite");
 
     scalar maxCos = Foam::cos(degToRad(angle));
@@ -636,171 +1216,196 @@ int main(int argc, char *argv[])
         << " degrees" << nl
         << endl;
 
+    Info<< "If an invalid mesh is generated then the edge length will be " << nl
+        << "multiplied by a factor of " << minLenFactor << " and collapsing "
+        << "will be reattempted" << nl << endl;
+
     bool meshChanged = false;
 
-    // Edge collapsing engine
-    edgeCollapser collapser(mesh);
+    checkMeshQuality(mesh);
 
-    label nIterations = 0;
+    autoPtr<fvMesh> fvMeshPtr;
 
-    while (true)
+    scalarList freezeEdges(mesh.nEdges(), 1.0);
+
+    do
     {
-        Info<< "Iteration " << nIterations << incrIndent << endl;
+        label nIterations = 0;
+        label nFrozenEdges = 0;
 
-        const faceList& faces = mesh.faces();
+        fvMeshPtr.reset(new fvMesh(mesh));
+        fvMesh& fvMeshRef = fvMeshPtr();
 
-        // boundaryPoint:
-        // + -1 : point not on boundary
-        // +  0 : point on a real boundary
-        // + >0 : point on a processor patch with that ID
-        labelList boundaryPoint(mesh.nPoints(), -1);
+        // Contains new point label for original points
+        labelList pointMap(identity(mesh.nPoints()));
 
-        // Get all points on a boundary
-        label nIntFaces = mesh.nInternalFaces();
-        for (label faceI = nIntFaces; faceI < mesh.nFaces(); faceI++)
+        scalarList tmpFreezeEdges = freezeEdges;
+
+        autoPtr<mapPolyMesh> morphMap;
+
+        while (true)
         {
-            const face& f = faces[faceI];
+            Info<< "Iteration " << nIterations << incrIndent << endl;
 
-            forAll(f, fp)
+            labelList boundaryPoint = findBoundaryPoints(fvMeshRef);
+
+            List<pointEdgeCollapse> allPointInfo;
+
+            // Collapse all edges that are too small.
+            label nSmallCollapsed =
+                collapseSmallEdges
+                (
+                    fvMeshRef,
+                    tmpFreezeEdges,
+                    boundaryPoint,
+                    minLen,
+                    allPointInfo
+                );
+
+            reduce(nSmallCollapsed, sumOp<label>());
+
+            Info<< indent << "Collapsing " << nSmallCollapsed
+                << " small edges" << endl;
+
+
+
+
+
+            label nMerged = 0;
+
+            // Remove midpoints on straight edges.
+            if (nSmallCollapsed == 0)
             {
-                boundaryPoint[f[fp]] = 0;
+                //nMerged = mergeEdges(fvMeshRef, maxCos, allPointInfo);
             }
-        }
 
-        // Get all processor boundary points and the processor patch label that
-        // they are on.
-        const polyBoundaryMesh& bMesh = mesh.boundaryMesh();
+            reduce(nMerged, sumOp<label>());
 
-        forAll(bMesh, patchI)
-        {
-            const polyPatch& patch = bMesh[patchI];
+            Info<< indent << "Collapsing " << nMerged << " in line edges"
+                << endl;
 
-            if (isA<processorPolyPatch>(patch))
+
+
+
+
+
+            label nSliversCollapsed = 0;
+
+            // Remove small sliver faces that can be collapsed to single edge
+//            if (nSmallCollapsed == 0 && nMerged == 0)
+//            {
+//                nSliversCollapsed =
+//                    collapseHighAspectFaces
+//                    (
+//                        mesh,
+//                        boundaryPoint,
+//                        processorPoints,
+//                        1E-9,
+//                        5,
+//                        collapser
+//                    );
+//            }
+
+            reduce(nSliversCollapsed, sumOp<label>());
+
+            Info<< indent << "Collapsing " << nSliversCollapsed
+                << " small high aspect ratio faces" << endl;
+
+
+            // Simplify faces to quads wherever possible
+            //if (nCollapsed == 0)
+            //{
+            //    nCollapsed =
+            //        simplifyFaces
+            //        (
+            //            mesh,
+            //            boundaryPoint,
+            //            4,              // minimum size of face
+            //            0.2,            // gap in edge lengths on face
+            //            collapser
+            //        );
+            //    Info<< "Collapsing " << nCollapsed << " polygonal faces"
+            //        << endl;
+            //}
+
+
+            label totalCollapsed =
+                nSmallCollapsed
+              + nMerged
+              + nSliversCollapsed;
+
+            polyTopoChange meshMod(fvMeshRef);
+
+            // Insert mesh refinement into polyTopoChange.
+            setRefinement(fvMeshRef, meshMod, allPointInfo);
+
+            // Do all changes
+            Info<< indent << "Applying changes to the mesh" << nl
+                << decrIndent << endl;
+
+            morphMap = meshMod.changeMesh(fvMeshRef, false);
+
+//            // Contains new point label for old points
+//            const labelList& reversePointMap = morphMap().reversePointMap();
+//
+//            forAll(pointMap, pI)
+//            {
+//                const label originalPoint = pI;
+//                const label currentPoint = pointMap[pI];
+//
+//                if (currentPoint < reversePointMap.size())
+//                {
+//                    const label newPoint = reversePointMap[currentPoint];
+//
+//                    if (newPoint != -1)
+//                    {
+//                        pointMap[originalPoint] = newPoint;
+//                    }
+//                }
+//            }
+
+            if (totalCollapsed == 0)
             {
-                const processorPolyPatch& pPatch =
-                    refCast<const processorPolyPatch>(patch);
+                labelHashSet freezePoints = checkMeshQuality(fvMeshRef);
 
-                forAll(pPatch, fI)
-                {
-                    const face& f = pPatch[fI];
+                label nFreezePoints = freezePoints.size();
+                reduce(nFreezePoints, sumOp<label>());
 
-                    forAll(f, fp)
-                    {
-                        boundaryPoint[f[fp]] = patchI;
-                    }
-                }
+                nFrozenEdges = nFreezePoints;
+
+                Info<< "Number of frozen points      : " << nFreezePoints
+                    << endl;
+
+                break;
             }
+
+            if (morphMap().hasMotionPoints())
+            {
+                fvMeshRef.movePoints(morphMap().preMotionPoints());
+            }
+
+            meshChanged = true;
+
+            nIterations++;
         }
 
-        // Collapse all edges that are too small.
-        label nSmallCollapsed =
-            collapseSmallEdges
-            (
-                mesh,
-                boundaryPoint,
-                minLen,
-                collapser
-            );
-
-        reduce(nSmallCollapsed, sumOp<label>());
-
-        Info<< indent << "Collapsing " << nSmallCollapsed
-            << " small edges" << endl;
-
-
-
-
-
-        label nMerged = 0;
-
-        // Remove midpoints on straight edges.
-        if (nSmallCollapsed == 0)
+        if (nFrozenEdges > 0)
         {
-            //nMerged = mergeEdges(mesh, maxCos, collapser);
+            minLen *= minLenFactor;
         }
 
-        reduce(nMerged, sumOp<label>());
+        reduce(nFrozenEdges, sumOp<label>());
 
-        Info<< indent << "Collapsing " << nMerged << " in line edges" << endl;
+        Info<< "Number of frozen edges       : " << nFrozenEdges << nl
+            << endl;
 
-
-
-
-
-
-        label nSliversCollapsed = 0;
-
-        // Remove small sliver faces that can be collapsed to single edge
-        if (nSmallCollapsed == 0 && nMerged == 0)
-        {
-//            nSliversCollapsed =
-//                collapseHighAspectFaces
-//                (
-//                    mesh,
-//                    boundaryPoint,
-//                    processorPoints,
-//                    1E-9,       // factor of largest face area
-//                    5,          // factor between smallest and largest edge on
-//                                // face
-//                    collapser
-//                );
-        }
-
-        reduce(nSliversCollapsed, sumOp<label>());
-
-        Info<< indent << "Collapsing " << nSliversCollapsed
-            << " small high aspect ratio faces" << endl;
-
-
-        // Simplify faces to quads wherever possible
-        //if (nCollapsed == 0)
-        //{
-        //    nCollapsed =
-        //        simplifyFaces
-        //        (
-        //            mesh,
-        //            boundaryPoint,
-        //            4,              // minimum size of face
-        //            0.2,            // gap in edge lengths on face
-        //            collapser
-        //        );
-        //    Info<< "Collapsing " << nCollapsed << " polygonal faces" << endl;
-        //}
-
-
-        label totalCollapsed =
-            nSmallCollapsed
-          + nMerged
-          + nSliversCollapsed;
-
-        if (totalCollapsed == 0)
+        if (nFrozenEdges == 0)
         {
             break;
         }
 
-        polyTopoChange meshMod(mesh);
+    } while (true);
 
-        // Insert mesh refinement into polyTopoChange.
-        collapser.setRefinement(meshMod);
-
-        // Do all changes
-        Info<< indent << "Applying changes to the mesh" << nl << endl;
-
-        autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh, false);
-
-        collapser.updateMesh(morphMap());
-
-        if (morphMap().hasMotionPoints())
-        {
-            mesh.movePoints(morphMap().preMotionPoints());
-        }
-
-        meshChanged = true;
-
-        Info<< decrIndent;
-
-        nIterations++;
-    }
 
     if (meshChanged)
     {
@@ -811,15 +1416,21 @@ int main(int argc, char *argv[])
         }
         else
         {
-            mesh.setInstance(oldInstance);
+            fvMeshPtr().setInstance(oldInstance);
         }
 
         Info<< nl << "Writing collapsed mesh to time "
             << runTime.timeName() << nl << endl;
 
-        mesh.write();
+        fvMeshPtr().write();
     }
 
+    Info<< "Final minimum length : " << minLen << " m" << nl << endl;
+
+    Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+        << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+        << nl << endl;
+
     Info<< "End\n" << endl;
 
     return 0;
diff --git a/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapse.C b/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapse.C
new file mode 100644
index 0000000000..6830e1c8f5
--- /dev/null
+++ b/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapse.C
@@ -0,0 +1,51 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "pointEdgeCollapse.H"
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<
+(
+    Foam::Ostream& os,
+    const Foam::pointEdgeCollapse& wDist
+)
+{
+    return os
+        << wDist.collapsePoint_ << wDist.collapseIndex_;
+}
+
+Foam::Istream& Foam::operator>>
+(
+    Foam::Istream& is,
+    Foam::pointEdgeCollapse& wDist
+)
+{
+    return is
+        >> wDist.collapsePoint_ >> wDist.collapseIndex_;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapse.H b/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapse.H
new file mode 100644
index 0000000000..cd8383b8b3
--- /dev/null
+++ b/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapse.H
@@ -0,0 +1,227 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::pointEdgeCollapse
+
+Description
+    Determines length of string of edges walked to point.
+
+SourceFiles
+    pointEdgeCollapseI.H
+    pointEdgeCollapse.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef pointEdgeCollapse_H
+#define pointEdgeCollapse_H
+
+#include "point.H"
+#include "tensor.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class polyPatch;
+class polyMesh;
+
+/*---------------------------------------------------------------------------*\
+                           Class pointEdgeCollapse Declaration
+\*---------------------------------------------------------------------------*/
+
+class pointEdgeCollapse
+{
+    // Private data
+
+        //- Collapse location
+        point collapsePoint_;
+
+        //- Collapse string index
+        label collapseIndex_;
+
+
+    // Private Member Functions
+
+        //- Evaluate distance to point.
+        template<class TrackingData>
+        inline bool update
+        (
+            const pointEdgeCollapse& w2,
+            const scalar tol,
+            TrackingData& td
+        );
+
+
+        //- Check for same coordinate
+        inline bool samePoint(const point& pt) const;
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        inline pointEdgeCollapse();
+
+        //- Construct from components
+        inline pointEdgeCollapse
+        (
+            const point& collapsePoint,
+            const label collapseIndex
+        );
+
+
+    // Member Functions
+
+        // Access
+
+            inline const point& collapsePoint() const;
+
+            inline label collapseIndex() const;
+
+
+        // Needed by meshWave
+
+            //- Check whether origin has been changed at all or
+            //  still contains original (invalid) value.
+            template<class TrackingData>
+            inline bool valid(TrackingData& td) const;
+
+            //- Convert origin to relative vector to leaving point
+            //  (= point coordinate)
+            template<class TrackingData>
+            inline void leaveDomain
+            (
+                const polyPatch& patch,
+                const label patchPointI,
+                const point& pos,
+                TrackingData& td
+            );
+
+            //- Convert relative origin to absolute by adding entering point
+            template<class TrackingData>
+            inline void enterDomain
+            (
+                const polyPatch& patch,
+                const label patchPointI,
+                const point& pos,
+                TrackingData& td
+            );
+
+            //- Apply rotation matrix to origin
+            template<class TrackingData>
+            inline void transform
+            (
+                const tensor& rotTensor,
+                TrackingData& td
+            );
+
+            //- Influence of edge on point
+            template<class TrackingData>
+            inline bool updatePoint
+            (
+                const polyMesh& mesh,
+                const label pointI,
+                const label edgeI,
+                const pointEdgeCollapse& edgeInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Influence of different value on same point.
+            //  Merge new and old info.
+            template<class TrackingData>
+            inline bool updatePoint
+            (
+                const polyMesh& mesh,
+                const label pointI,
+                const pointEdgeCollapse& newPointInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Influence of different value on same point.
+            //  No information about current position whatsoever.
+            template<class TrackingData>
+            inline bool updatePoint
+            (
+                const pointEdgeCollapse& newPointInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Influence of point on edge.
+            template<class TrackingData>
+            inline bool updateEdge
+            (
+                const polyMesh& mesh,
+                const label edgeI,
+                const label pointI,
+                const pointEdgeCollapse& pointInfo,
+                const scalar tol,
+                TrackingData& td
+            );
+
+            //- Same (like operator==)
+            template<class TrackingData>
+            inline bool equal(const pointEdgeCollapse&, TrackingData&)
+            const;
+
+
+    // Member Operators
+
+        //Note: Used to determine whether to call update.
+        inline bool operator==(const pointEdgeCollapse&) const;
+        inline bool operator!=(const pointEdgeCollapse&) const;
+
+
+    // IOstream Operators
+
+        friend Ostream& operator<<(Ostream&, const pointEdgeCollapse&);
+        friend Istream& operator>>(Istream&, pointEdgeCollapse&);
+};
+
+
+//- Data associated with pointEdgeCollapse type are contiguous
+template<>
+inline bool contiguous<pointEdgeCollapse>()
+{
+    return true;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "pointEdgeCollapseI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapseI.H b/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapseI.H
new file mode 100644
index 0000000000..3fb1468a88
--- /dev/null
+++ b/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapseI.H
@@ -0,0 +1,313 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "polyMesh.H"
+#include "transform.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// Update this with w2.
+template<class TrackingData>
+inline bool Foam::pointEdgeCollapse::update
+(
+    const pointEdgeCollapse& w2,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    if (!w2.valid(td))
+    {
+        FatalErrorIn("pointEdgeCollapse::update(..)")
+            << "problem." << abort(FatalError);
+    }
+
+    if (w2.collapseIndex_ == -1)
+    {
+        // Not marked for collapse; only happens on edges.
+        return false;
+    }
+
+    if (!valid(td))
+    {
+        operator=(w2);
+        return true;
+    }
+    else
+    {
+        // Same coordinate. Same string?
+        if (w2.collapseIndex_ < collapseIndex_)
+        {
+            // Take over string index from w2 (and also coordinate but this
+            // was same)
+            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
+        {
+            return false;
+        }
+
+//        if (samePoint(w2.collapsePoint_))
+//        {
+//            // Same coordinate. Same string?
+//            if (w2.collapseIndex_ < collapseIndex_)
+//            {
+//                // Take over string index from w2 (and also coordinate but this
+//                // was same)
+//                operator=(w2);
+//                return true;
+//            }
+//            else
+//            {
+//                return false;
+//            }
+//        }
+//        else
+//        {
+//            // Find nearest coordinate
+//            if (magSqr(w2.collapsePoint_) < magSqr(collapsePoint_))
+//            {
+//                operator=(w2);
+//                return true;
+//            }
+//            else
+//            {
+//                return false;
+//            }
+//        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Null constructor
+inline Foam::pointEdgeCollapse::pointEdgeCollapse()
+:
+    collapsePoint_(GREAT, GREAT, GREAT),
+    collapseIndex_(-2)
+{}
+
+
+// Construct from origin, distance
+inline Foam::pointEdgeCollapse::pointEdgeCollapse
+(
+    const point& collapsePoint,
+    const label collapseIndex
+)
+:
+    collapsePoint_(collapsePoint),
+    collapseIndex_(collapseIndex)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const Foam::point& Foam::pointEdgeCollapse::collapsePoint() const
+{
+    return collapsePoint_;
+}
+
+
+inline Foam::label Foam::pointEdgeCollapse::collapseIndex() const
+{
+    return collapseIndex_;
+}
+
+
+inline bool Foam::pointEdgeCollapse::samePoint(const point& pt) const
+{
+    bool isLegal1 = (cmptMin(collapsePoint_) < 0.5*GREAT);
+    bool isLegal2 = (cmptMin(pt) < 0.5*GREAT);
+
+    if (isLegal1 && isLegal2)
+    {
+        return mag(collapsePoint_ - pt) < 1e-9;//SMALL;
+    }
+    else
+    {
+        return isLegal1 == isLegal2;
+    }
+}
+
+
+template<class TrackingData>
+inline bool Foam::pointEdgeCollapse::valid(TrackingData& td) const
+{
+    return collapseIndex_ != -2;
+}
+
+
+template<class TrackingData>
+inline void Foam::pointEdgeCollapse::leaveDomain
+(
+    const polyPatch& patch,
+    const label patchPointI,
+    const point& coord,
+    TrackingData& td
+)
+{
+    collapsePoint_ -= coord;
+}
+
+
+template<class TrackingData>
+inline void Foam::pointEdgeCollapse::transform
+(
+    const tensor& rotTensor,
+    TrackingData& td
+)
+{
+    collapsePoint_ = Foam::transform(rotTensor, collapsePoint_);
+}
+
+
+// Update absolute geometric quantities. Note that distance (dist_)
+// is not affected by leaving/entering domain.
+template<class TrackingData>
+inline void Foam::pointEdgeCollapse::enterDomain
+(
+    const polyPatch& patch,
+    const label patchPointI,
+    const point& coord,
+    TrackingData& td
+)
+{
+    // back to absolute form
+    collapsePoint_ += coord;
+}
+
+
+// Update this with information from connected edge
+template<class TrackingData>
+inline bool Foam::pointEdgeCollapse::updatePoint
+(
+    const polyMesh& mesh,
+    const label pointI,
+    const label edgeI,
+    const pointEdgeCollapse& edgeInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    return update(edgeInfo, tol, td);
+}
+
+
+// Update this with new information on same point
+template<class TrackingData>
+inline bool Foam::pointEdgeCollapse::updatePoint
+(
+    const polyMesh& mesh,
+    const label pointI,
+    const pointEdgeCollapse& newPointInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    return update(newPointInfo, tol, td);
+}
+
+
+// Update this with new information on same point. No extra information.
+template<class TrackingData>
+inline bool Foam::pointEdgeCollapse::updatePoint
+(
+    const pointEdgeCollapse& newPointInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    return update(newPointInfo, tol, td);
+}
+
+
+// Update this with information from connected point
+template<class TrackingData>
+inline bool Foam::pointEdgeCollapse::updateEdge
+(
+    const polyMesh& mesh,
+    const label edgeI,
+    const label pointI,
+    const pointEdgeCollapse& pointInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    return update(pointInfo, tol, td);
+}
+
+
+template <class TrackingData>
+inline bool Foam::pointEdgeCollapse::equal
+(
+    const pointEdgeCollapse& rhs,
+    TrackingData& td
+) const
+{
+    return operator==(rhs);
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+inline bool Foam::pointEdgeCollapse::operator==
+(
+    const Foam::pointEdgeCollapse& rhs
+) const
+{
+    return
+        collapseIndex_ == rhs.collapseIndex_
+     && samePoint(rhs.collapsePoint_);
+}
+
+
+inline bool Foam::pointEdgeCollapse::operator!=
+(
+    const Foam::pointEdgeCollapse& rhs
+) const
+{
+    return !(*this == rhs);
+}
+
+
+// ************************************************************************* //