From f94ce803739daa02656492422b99629f41c7dde4 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 15 May 2012 16:24:01 +0100
Subject: [PATCH 01/16] BUG: cylindricalInletVelocity: use of tmp after storage
 has been reused

---
 .../cylindricalInletVelocityFvPatchVectorField.C              | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/finiteVolume/fields/fvPatchFields/derived/cylindricalInletVelocity/cylindricalInletVelocityFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/cylindricalInletVelocity/cylindricalInletVelocityFvPatchVectorField.C
index 0348035314..d41c428e74 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/cylindricalInletVelocity/cylindricalInletVelocityFvPatchVectorField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/cylindricalInletVelocity/cylindricalInletVelocityFvPatchVectorField.C
@@ -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) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -132,7 +132,7 @@ void Foam::cylindricalInletVelocityFvPatchVectorField::updateCoeffs()
     vector hatAxis = axis_/mag(axis_);
 
     const vectorField r(patch().Cf() - centre_);
-    tmp<vectorField> d =  r - (hatAxis & r)*hatAxis;
+    const vectorField d =  r - (hatAxis & r)*hatAxis;
 
     tmp<vectorField> tangVel
     (

From 41308b96961bfd09f515225e61b97e430b741487 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Wed, 16 May 2012 10:16:18 +0100
Subject: [PATCH 02/16] ENH: collapseEdges: parallel operation

---
 .../advanced/collapseEdges/collapseEdges.C    | 1363 +++++++++--------
 .../pointEdgeCollapse/pointEdgeCollapseI.H    |   37 +-
 2 files changed, 693 insertions(+), 707 deletions(-)

diff --git a/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C b/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C
index 899159b2a2..f970c1d836 100644
--- a/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C
+++ b/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C
@@ -28,27 +28,22 @@ Description
     - merge two edges if they are in line. Maximum angle provided as argument.
     - remove unused points.
 
-    Cannot remove cells. Can remove faces and points but does not check
+    Optionally removes cells. Can remove faces and points but does not check
     for nonsense resulting topology.
 
     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. Note: it should in fact use features
-    where if one point is on a feature it collapses to that one. Alas we don't
-    have features on a polyMesh.
+    both points on boundary random again.
 
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
 #include "Time.H"
 #include "polyTopoChange.H"
-#include "polyTopoChanger.H"
-#include "polyMesh.H"
 #include "fvMesh.H"
 #include "mapPolyMesh.H"
 #include "mathematicalConstants.H"
 #include "PackedBoolList.H"
-#include "SortableList.H"
 #include "unitConversion.H"
 #include "globalMeshData.H"
 #include "globalIndex.H"
@@ -56,32 +51,10 @@ Description
 #include "pointEdgeCollapse.H"
 #include "motionSmoother.H"
 
-#include "OFstream.H"
-#include "meshTools.H"
-
 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,
@@ -121,7 +94,7 @@ void filterFace
         }
         else
         {
-            if (findIndex(f, 0, newFp, collapsePoint) == -1)
+            if (findIndex(SubList<label>(f, newFp), collapsePoint) == -1)
             {
                 f[newFp++] = collapsePoint;
             }
@@ -147,7 +120,7 @@ void filterFace
         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)
         {
@@ -206,19 +179,50 @@ bool setRefinement
 
     boolList removedPoints(mesh.nPoints(), false);
 
-    // Create strings of edges
+    // Create strings of edges.
+    // Map from collapseIndex(=global master point) to set of points
     Map<DynamicList<label> > collapseStrings;
-
-    forAll(allPointInfo, pointI)
     {
-        const label collapseIndex = allPointInfo[pointI].collapseIndex();
-
-        if (collapseIndex != -1 && collapseIndex != -2)
+        // 1. Count elements per collapseIndex
+        Map<label> nPerIndex(mesh.nPoints()/10);
+        forAll(allPointInfo, pointI)
         {
-            collapseStrings(collapseIndex).append(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)
@@ -293,7 +297,7 @@ bool setRefinement
 
                         if (nFaces < 4)
                         {
-                            Info<< "Cell:" << cellI
+                            Pout<< "Cell:" << cellI
                                 << " uses faces:" << cFaces
                                 << " of which too many are marked for removal:"
                                 << endl
@@ -302,10 +306,10 @@ bool setRefinement
                             {
                                 if (newFaces[cFaces[j]].size() < 3)
                                 {
-                                    Info<< ' '<< cFaces[j];
+                                    Pout<< ' '<< cFaces[j];
                                 }
                             }
-                            Info<< endl;
+                            Pout<< endl;
 
                             cellRemoved[cellI] = true;
 
@@ -463,9 +467,6 @@ bool setRefinement
         }
     }
 
-    // Print regions:
-//    printRegions();
-
     return meshChanged;
 }
 
@@ -507,99 +508,6 @@ labelList getSortedEdges
 }
 
 
-// Merges edges which are in straight line. I.e. edge split by point.
-//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)
 label edgeMaster
 (
@@ -673,112 +581,22 @@ label edgeMaster
 }
 
 
-label collapseSmallEdges
+// 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 scalarList& freezeEdges,
-    const labelList& boundaryPoint,
-    const scalar minLen,
+    const globalIndex& globalStrings,
+    const labelList& collapseEdge,
     List<pointEdgeCollapse>& allPointInfo
 )
 {
-    const pointField& points = mesh.points();
     const edgeList& edges = mesh.edges();
-
-    // Store collapse direction in collapseEdge
-    //   -1 -> Do not collapse
-    //    0 -> Collapse to start point
-    //    1 -> Collapse to end point
-    labelList collapseEdge(edges.size(), -1);
-
-    forAll(edges, edgeI)
-    {
-        const edge& e = edges[edgeI];
-
-        if (e.mag(points) < minLen*freezeEdges[edgeI])
-        {
-            collapseEdge[edgeI] = 0;
-        }
-    }
-
-    // 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 = 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;
-//            }
-//        }
-//    }
+    const pointField& points = mesh.points();
 
     label nCollapsed = 0;
 
@@ -786,46 +604,43 @@ label collapseSmallEdges
     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);
-    }
+    // 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)
         {
-            const label master =
-                edgeMaster
-                (
-                    boundaryPoint,
-                    collapseEdge[edgeI],
-                    e
-                );
+            label masterPointI = e[collapseEdge[edgeI]];
 
-//            if (master != -1)
-            {
-                pointEdgeCollapse pec
-                (
-                    points[master],
-                    globalStrings.toGlobal(master)
-                );
+            const pointEdgeCollapse pec
+            (
+                points[masterPointI],
+                globalStrings.toGlobal(masterPointI)
+            );
 
-                allEdgeInfo[edgeI] = pec;
+            // 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.start());
 
-                initPointInfo.append(pec);
-                initPoints.append(e.end());
+            initPointInfo.append(pec);
+            initPoints.append(e.end());
 
-                nCollapsed++;
-            }
+            nCollapsed++;
         }
     }
 
@@ -843,223 +658,336 @@ label collapseSmallEdges
 }
 
 
-// Faces which have edges just larger than collapse length but faces which
-// are very small. This one tries to collapse them if it can be done with
-// edge collapse. For faces where a face gets replace by two edges use
-// collapseFaces
-//label collapseHighAspectFaces
-//(
-//    const polyMesh& mesh,
-//    const PackedBoolList& boundaryPoint,
-//    const Map<label>& processorPoints,
-//    const scalar areaFac,
-//    const scalar edgeRatio,
-//    edgeCollapser& collapser
-//)
-//{
-//    const pointField& points = mesh.points();
-//    const edgeList& edges = mesh.edges();
-//    const faceList& faces = mesh.faces();
-//    const labelListList& faceEdges = mesh.faceEdges();
-//
-//    scalarField magArea(mag(mesh.faceAreas()));
-//
-//    label maxIndex = findMax(magArea);
-//
-//    scalar minArea = areaFac * magArea[maxIndex];
-//
-//    Info<< "Max face area:" << magArea[maxIndex] << endl
-//        << "Collapse area factor:" << areaFac << endl
-//        << "Collapse area:" << minArea << endl;
-//
-//    label nCollapsed = 0;
-//
-//    forAll(faces, faceI)
-//    {
-//        if (magArea[faceI] < minArea)
-//        {
-//            const face& f = faces[faceI];
-//
-//            // Get the edges in face point order
-//            labelList fEdges(getSortedEdges(edges, f, faceEdges[faceI]));
-//
-//            SortableList<scalar> lengths(fEdges.size());
-//            forAll(fEdges, i)
-//            {
-//                lengths[i] = edges[fEdges[i]].mag(points);
-//            }
-//            lengths.sort();
-//
-//
-//            label edgeI = -1;
-//
-//            if (f.size() == 4)
-//            {
-//                // Compare second largest to smallest
-//                if (lengths[2] > edgeRatio*lengths[0])
-//                {
-//                    // Collapse smallest only. Triangle should be cleared
-//                    // next time around.
-//                    edgeI = fEdges[lengths.indices()[0]];
-//                }
-//            }
-//            else if (f.size() == 3)
-//            {
-//                // Compare second largest to smallest
-//                if (lengths[1] > edgeRatio*lengths[0])
-//                {
-//                    edgeI = fEdges[lengths.indices()[0]];
-//                }
-//            }
-//
-//
-//            if (edgeI != -1)
-//            {
-//                label master = edgeMaster
-//                    (
-//                        boundaryPoint,
-//                        processorPoints,
-//                        false,
-//                        edges[edgeI]
-//                    );
-//
-//                if (master != -1)// && collapser.unaffectedEdge(edgeI))
-//                {
-//                    collapser.collapseEdge(edgeI, master);
-//                    nCollapsed++;
-//                }
-//            }
-//        }
-//    }
-//
-//    return nCollapsed;
-//}
-
-
-void set(const labelList& elems, const bool val, boolList& status)
+void syncCollapseEdge(const polyMesh& mesh, labelList& collapseEdge)
 {
-    forAll(elems, i)
+    // 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)
     {
-        status[elems[i]] = val;
+        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];
+        }
     }
 }
 
 
-// Tries to simplify polygons to face of minSize (4=quad, 3=triangle)
-//label simplifyFaces
-//(
-//    const polyMesh& mesh,
-//    const PackedBoolList& boundaryPoint,
-//    const Map<label>& processorPoints,
-//    const label minSize,
-//    const scalar lenGap,
-//    edgeCollapser& collapser
-//)
-//{
-//    const pointField& points = mesh.points();
-//    const edgeList& edges = mesh.edges();
-//    const faceList& faces = mesh.faces();
-//    const cellList& cells = mesh.cells();
-//    const labelListList& faceEdges = mesh.faceEdges();
-//    const labelList& faceOwner = mesh.faceOwner();
-//    const labelList& faceNeighbour = mesh.faceNeighbour();
-//    const labelListList& pointCells = mesh.pointCells();
-//    const labelListList& cellEdges = mesh.cellEdges();
-//
-//    label nCollapsed = 0;
-//
-//    boolList protectedEdge(mesh.nEdges(), false);
-//
-//    forAll(faces, faceI)
-//    {
-//        const face& f = faces[faceI];
-//
-//        if
-//        (
-//            f.size() > minSize
-//         && cells[faceOwner[faceI]].size() >= 6
-//         && (
-//                mesh.isInternalFace(faceI)
-//             && cells[faceNeighbour[faceI]].size() >= 6
-//            )
-//        )
-//        {
-//            // Get the edges in face point order
-//            labelList fEdges(getSortedEdges(edges, f, faceEdges[faceI]));
-//
-//            SortableList<scalar> lengths(fEdges.size());
-//            forAll(fEdges, i)
-//            {
-//                lengths[i] = edges[fEdges[i]].mag(points);
-//            }
-//            lengths.sort();
-//
-//
-//            // Now find a gap in length between consecutive elements greater
-//            // than lenGap.
-//
-//            label gapPos = -1;
-//
-//            for (label i = f.size()-1-minSize; i >= 0; --i)
-//            {
-//                if (lengths[i+1] > lenGap*lengths[i])
-//                {
-//                    gapPos = i;
-//
-//                    break;
-//                }
-//            }
-//
-//            if (gapPos != -1)
-//            {
-//                //for (label i = gapPos; i >= 0; --i)
-//                label i = 0;  // Hack: collapse smallest edge only.
-//                {
-//                    label edgeI = fEdges[lengths.indices()[i]];
-//
-//                    if (!protectedEdge[edgeI])
-//                    {
-//                        const edge& e = edges[edgeI];
-//
-//                        label master
-//                       = edgeMaster(boundaryPoint, processorPoints, false, e);
-//
-//                        if (master != -1)
-//                        {
-//                            collapser.collapseEdge(edgeI, master);
-//
-//                            // Protect all other edges on all cells using edge
-//                            // points.
-//
-//                            const labelList& pCells0 = pointCells[e[0]];
-//
-//                            forAll(pCells0, i)
-//                            {
-//                              set(cellEdges[pCells0[i]], true, protectedEdge);
-//                            }
-//                            const labelList& pCells1 = pointCells[e[1]];
-//
-//                            forAll(pCells1, i)
-//                            {
-//                              set(cellEdges[pCells1[i]], true, protectedEdge);
-//                            }
-//
-//                            nCollapsed++;
-//                        }
-//                    }
-//                }
-//            }
-//        }
-//    }
-//
-//    return nCollapsed;
-//}
+// 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;
+}
 
 
-labelHashSet checkMeshQuality(const polyMesh& mesh)
+// Mark (in collapseEdge) any edges to merge
+label mergeEdges
+(
+    const polyMesh& mesh,
+    const scalar maxCos,
+    const labelList& boundaryPoint,
+    const scalarField& minEdgeLen,
+    labelList& collapseEdge
+)
+{
+    const pointField& points = mesh.points();
+    const edgeList& edges = mesh.edges();
+    const labelListList& pointEdges = mesh.pointEdges();
+
+    label nCollapsed = 0;
+
+    forAll(pointEdges, pointI)
+    {
+        const labelList& pEdges = pointEdges[pointI];
+
+        if (pEdges.size() == 2 && boundaryPoint[pointI] <= 0)
+        {
+            // Collapse only if none of the points part of merge network and
+            // none protected (minEdgeLen < 0)
+            label e0 = pEdges[0];
+            label e1 = pEdges[1];
+
+            if
+            (
+                collapseEdge[e0] == -1
+             && minEdgeLen[e0] >= 0
+             && collapseEdge[e1] == -1
+             && minEdgeLen[e1] >= 0
+            )
+            {
+                const edge& leftE = edges[e0];
+                const edge& rightE = edges[e1];
+
+                // Get the two vertices on both sides of the point
+                label leftV = leftE.otherVertex(pointI);
+                label rightV = rightE.otherVertex(pointI);
+
+                // 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 pointI
+                    // the master.
+                    collapseEdge[e0] = findIndex(leftE, pointI);
+                    nCollapsed++;
+                }
+            }
+        }
+    }
+    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);
-    labelHashSet freezePoints;
 
     IOdictionary meshQualityDict
     (
@@ -1102,35 +1030,61 @@ labelHashSet checkMeshQuality(const polyMesh& mesh)
 
     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)
         {
-            freezePoints.insert(f[pI]);
+            isErrorPoint[f[pI]] = 1;
         }
     }
 
-//    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++;
-//        }
-//    }
+    syncTools::syncPointList
+    (
+        mesh,
+        isErrorPoint,
+        orEqOp<unsigned int>(),
+        0
+    );
 
-    return freezePoints;
+    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++;
+            }
+        }
+    }
 }
 
 
@@ -1165,14 +1119,21 @@ labelList findBoundaryPoints(const polyMesh& mesh)
     {
         const polyPatch& patch = bMesh[patchI];
 
-        if (isA<processorPolyPatch>(patch))
+        if (isA<coupledPolyPatch>(patch))
         {
-            const processorPolyPatch& pPatch =
-                refCast<const processorPolyPatch>(patch);
-
-            forAll(pPatch, fI)
+            if (patchI == 0)
             {
-                const face& f = pPatch[fI];
+                // 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)
                 {
@@ -1189,11 +1150,28 @@ labelList findBoundaryPoints(const polyMesh& mesh)
 // Main program:
 int main(int argc, char *argv[])
 {
+    argList::addNote
+    (
+        "Merges small and in-line edges.\n"
+        "Collapses faces and optionally cells to a point."
+    );
+
 #   include "addOverwriteOption.H"
 
     argList::validArgs.append("edge length [m]");
     argList::validArgs.append("merge angle (degrees)");
-    argList::addOption("minLenFactor", "scalar", "edge length factor");
+    argList::addBoolOption
+    (
+        "allowCellCollapse",
+        "Allow collapsing of cells to a single point"
+    );
+    argList::addBoolOption
+    (
+        "checkMeshQuality",
+        "Only collapse if not exceeding given meshQualityDict limits"
+    );
+
+
 
 #   include "setRootCase.H"
 #   include "createTime.H"
@@ -1201,12 +1179,12 @@ int main(int argc, char *argv[])
 #   include "createPolyMesh.H"
     const word oldInstance = mesh.pointsInstance();
 
-    scalar minLen  = args.argRead<scalar>(1);
-    const scalar angle   = args.argRead<scalar>(2);
-    const scalar minLenFactor
-        = args.optionLookupOrDefault<scalar>("minLenFactor", 0.5);
+    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("checkMeshQuality");
 
     scalar maxCos = Foam::cos(degToRad(angle));
 
@@ -1216,190 +1194,216 @@ 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;
+    if (allowCellCollapse)
+    {
+        Info<< "Allowing collapse of cells down to a point." << nl
+            << endl;
+    }
+    else
+    {
+        Info<< "Disallowing collapse of cells down to a point." << nl
+            << endl;
+    }
 
+
+    if (checkQuality)
+    {
+        Info<< "Selectively disabling wanted collapses until resulting quality"
+            << " satisfies constraints in system/meshQualityDict" << nl
+            << endl;
+    }
+
+
+
+    // 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;
 
-    checkMeshQuality(mesh);
 
-    autoPtr<fvMesh> fvMeshPtr;
 
-    scalarList freezeEdges(mesh.nEdges(), 1.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.
 
     do
     {
-        label nIterations = 0;
-        label nFrozenEdges = 0;
-
-        fvMeshPtr.reset(new fvMesh(mesh));
-        fvMesh& fvMeshRef = fvMeshPtr();
-
-        // Contains new point label for original points
-        labelList pointMap(identity(mesh.nPoints()));
-
-        scalarList tmpFreezeEdges = freezeEdges;
-
-        autoPtr<mapPolyMesh> morphMap;
-
-        while (true)
-        {
-            Info<< "Iteration " << nIterations << incrIndent << endl;
-
-            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;
+        // 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
+        );
 
-            label nMerged = 0;
-
-            // Remove midpoints on straight edges.
-            if (nSmallCollapsed == 0)
-            {
-                //nMerged = mergeEdges(fvMeshRef, maxCos, allPointInfo);
-            }
-
-            reduce(nMerged, sumOp<label>());
-
-            Info<< indent << "Collapsing " << nMerged << " in line edges"
-                << 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,
-//                        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)
-            {
-                labelHashSet freezePoints = checkMeshQuality(fvMeshRef);
-
-                label nFreezePoints = freezePoints.size();
-                reduce(nFreezePoints, sumOp<label>());
-
-                nFrozenEdges = nFreezePoints;
-
-                Info<< "Number of frozen points      : " << nFreezePoints
-                    << endl;
-
-                break;
-            }
-
-            if (morphMap().hasMotionPoints())
-            {
-                fvMeshRef.movePoints(morphMap().preMotionPoints());
-            }
-
-            meshChanged = true;
-
-            nIterations++;
-        }
-
-        if (nFrozenEdges > 0)
-        {
-            minLen *= minLenFactor;
-        }
-
-        reduce(nFrozenEdges, sumOp<label>());
-
-        Info<< "Number of frozen edges       : " << nFrozenEdges << nl
+        reduce(nMerged, sumOp<label>());
+        Info<< indent << "Collapsing " << nMerged << " in line edges"
             << endl;
 
-        if (nFrozenEdges == 0)
+
+        // 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)
+        {
+            break;
+        }
+
+
+        // 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())
+        {
+            newMesh.movePoints(map.preMotionPoints());
+        }
+
+
+
+
+        // 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;
         }
@@ -1409,6 +1413,18 @@ int main(int argc, char *argv[])
 
     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)
         {
@@ -1416,16 +1432,15 @@ int main(int argc, char *argv[])
         }
         else
         {
-            fvMeshPtr().setInstance(oldInstance);
+            mesh.setInstance(oldInstance);
         }
 
         Info<< nl << "Writing collapsed mesh to time "
             << runTime.timeName() << nl << endl;
 
-        fvMeshPtr().write();
+        mesh.write();
     }
 
-    Info<< "Final minimum length : " << minLen << " m" << nl << endl;
 
     Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
         << "  ClockTime = " << runTime.elapsedClockTime() << " s"
diff --git a/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapseI.H b/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapseI.H
index a064bcf781..eeb4d6baf2 100644
--- a/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapseI.H
+++ b/applications/utilities/mesh/advanced/collapseEdges/pointEdgeCollapse/pointEdgeCollapseI.H
@@ -43,7 +43,7 @@ inline bool Foam::pointEdgeCollapse::update
             << "problem." << abort(FatalError);
     }
 
-    if (w2.collapseIndex_ == -1)
+    if (w2.collapseIndex_ == -1 || collapseIndex_ == -1)
     {
         // Not marked for collapse; only happens on edges.
         return false;
@@ -56,11 +56,11 @@ inline bool Foam::pointEdgeCollapse::update
     }
     else
     {
-        // Same coordinate. Same string?
+        // Take over w2 if it is 'better'
+
         if (w2.collapseIndex_ < collapseIndex_)
         {
-            // Take over string index from w2 (and also coordinate but this
-            // was same)
+            // Take over string index and coordinate from w2
             operator=(w2);
             return true;
         }
@@ -85,35 +85,6 @@ inline bool Foam::pointEdgeCollapse::update
         {
             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;
-//            }
-//        }
     }
 }
 

From 62a203c9db6fb0052a8aeb5606650a32e811a42e Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Wed, 16 May 2012 11:12:49 +0100
Subject: [PATCH 03/16] ENH: checkMesh: add meshQualityDict option

---
 .../mesh/manipulation/checkMesh/Make/files    |  1 +
 .../mesh/manipulation/checkMesh/Make/options  |  6 +-
 .../mesh/manipulation/checkMesh/checkMesh.C   | 87 +++++++++++++++++--
 .../manipulation/checkMesh/checkMeshQuality.C | 34 ++++++++
 .../manipulation/checkMesh/checkMeshQuality.H |  6 ++
 5 files changed, 125 insertions(+), 9 deletions(-)
 create mode 100644 applications/utilities/mesh/manipulation/checkMesh/checkMeshQuality.C
 create mode 100644 applications/utilities/mesh/manipulation/checkMesh/checkMeshQuality.H

diff --git a/applications/utilities/mesh/manipulation/checkMesh/Make/files b/applications/utilities/mesh/manipulation/checkMesh/Make/files
index f0b7c16694..1a3130a23c 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/Make/files
+++ b/applications/utilities/mesh/manipulation/checkMesh/Make/files
@@ -1,6 +1,7 @@
 printMeshStats.C
 checkTopology.C
 checkGeometry.C
+checkMeshQuality.C
 checkMesh.C
 
 EXE = $(FOAM_APPBIN)/checkMesh
diff --git a/applications/utilities/mesh/manipulation/checkMesh/Make/options b/applications/utilities/mesh/manipulation/checkMesh/Make/options
index 54c035b8f5..70c838b774 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/Make/options
+++ b/applications/utilities/mesh/manipulation/checkMesh/Make/options
@@ -1,5 +1,7 @@
 EXE_INC = \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude
 
 EXE_LIBS = \
-    -lmeshTools
+    -lmeshTools \
+    -ldynamicMesh
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkMesh.C b/applications/utilities/mesh/manipulation/checkMesh/checkMesh.C
index 3f0ae4589c..dc400f5f41 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkMesh.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkMesh.C
@@ -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) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,6 +27,21 @@ Application
 Description
     Checks validity of a mesh
 
+Usage
+    - checkMesh [OPTION]
+
+    \param -allGeometry \n
+    Checks all (including non finite-volume specific) geometry
+
+    \param -allTopology \n
+    Checks all (including non finite-volume specific) addressing
+
+    \param -meshQualityDict \n
+    Checks against user defined (in \a system/meshQualityDict) quality settings
+
+    \param -region \<name\> \n
+    Specify an alternative mesh region.
+
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
@@ -39,6 +54,7 @@ Description
 #include "printMeshStats.H"
 #include "checkTopology.H"
 #include "checkGeometry.H"
+#include "checkMeshQuality.H"
 
 using namespace Foam;
 
@@ -63,6 +79,11 @@ int main(int argc, char *argv[])
         "allTopology",
         "include extra topology checks"
     );
+    argList::addBoolOption
+    (
+        "meshQualityDict",
+        "read user-defined mesh quality criterions from system/meshQualityDict"
+    );
 
 #   include "setRootCase.H"
 #   include "createTime.H"
@@ -72,6 +93,46 @@ int main(int argc, char *argv[])
     const bool noTopology  = args.optionFound("noTopology");
     const bool allGeometry = args.optionFound("allGeometry");
     const bool allTopology = args.optionFound("allTopology");
+    const bool meshQualityDict = args.optionFound("meshQualityDict");
+
+    if (noTopology)
+    {
+        Info<< "Disabling all topology checks." << nl << endl;
+    }
+    if (allTopology)
+    {
+        Info<< "Enabling all (cell, face, edge, point) topology checks."
+            << nl << endl;
+    }
+    if (allGeometry)
+    {
+        Info<< "Enabling all geometry checks." << nl << endl;
+    }
+    if (meshQualityDict)
+    {
+        Info<< "Enabling user-defined geometry checks." << nl << endl;
+    }
+
+
+    autoPtr<IOdictionary> qualDict;
+    if (meshQualityDict)
+    {
+        qualDict.reset
+        (
+            new IOdictionary
+            (
+                IOobject
+                (
+                    "meshQualityDict",
+                    mesh.time().system(),
+                    mesh,
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE
+                )
+           )
+        );
+    }
+
 
     forAll(timeDirs, timeI)
     {
@@ -96,25 +157,31 @@ int main(int argc, char *argv[])
 
             printMeshStats(mesh, allTopology);
 
-            label noFailedChecks = 0;
+            label nFailedChecks = 0;
 
             if (!noTopology)
             {
-                noFailedChecks += checkTopology(mesh, allTopology, allGeometry);
+                nFailedChecks += checkTopology(mesh, allTopology, allGeometry);
             }
 
-            noFailedChecks += checkGeometry(mesh, allGeometry);
+            nFailedChecks += checkGeometry(mesh, allGeometry);
 
-            // Note: no reduction in noFailedChecks necessary since is
+            if (meshQualityDict)
+            {
+                nFailedChecks += checkMeshQuality(mesh, qualDict());
+            }
+
+
+            // Note: no reduction in nFailedChecks necessary since is
             //       counter of checks, not counter of failed cells,faces etc.
 
-            if (noFailedChecks == 0)
+            if (nFailedChecks == 0)
             {
                 Info<< "\nMesh OK.\n" << endl;
             }
             else
             {
-                Info<< "\nFailed " << noFailedChecks << " mesh checks.\n"
+                Info<< "\nFailed " << nFailedChecks << " mesh checks.\n"
                     << endl;
             }
         }
@@ -124,6 +191,12 @@ int main(int argc, char *argv[])
 
             label nFailedChecks = checkGeometry(mesh, allGeometry);
 
+            if (meshQualityDict)
+            {
+                nFailedChecks += checkMeshQuality(mesh, qualDict());
+            }
+
+
             if (nFailedChecks)
             {
                 Info<< "\nFailed " << nFailedChecks << " mesh checks.\n"
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkMeshQuality.C b/applications/utilities/mesh/manipulation/checkMesh/checkMeshQuality.C
new file mode 100644
index 0000000000..b00f02f8ef
--- /dev/null
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkMeshQuality.C
@@ -0,0 +1,34 @@
+#include "checkMeshQuality.H"
+#include "polyMesh.H"
+#include "cellSet.H"
+#include "faceSet.H"
+#include "motionSmoother.H"
+
+
+Foam::label Foam::checkMeshQuality
+(
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+{
+    label noFailedChecks = 0;
+
+    {
+        faceSet faces(mesh, "meshQualityFaces", mesh.nFaces()/100+1);
+        motionSmoother::checkMesh(false, mesh, dict, faces);
+
+        label nFaces = returnReduce(faces.size(), sumOp<label>());
+
+        if (nFaces > 0)
+        {
+            noFailedChecks++;
+
+            Info<< "  <<Writing " << nFaces
+                << " faces in error to set " << faces.name() << endl;
+            faces.instance() = mesh.pointsInstance();
+            faces.write();
+        }
+    }
+
+    return noFailedChecks;
+}
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkMeshQuality.H b/applications/utilities/mesh/manipulation/checkMesh/checkMeshQuality.H
new file mode 100644
index 0000000000..1e5b3489c9
--- /dev/null
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkMeshQuality.H
@@ -0,0 +1,6 @@
+#include "polyMesh.H"
+
+namespace Foam
+{
+    label checkMeshQuality(const polyMesh& mesh, const dictionary&);
+}

From e14d8d14d756120660063a8bef9ff93ce0b0dffb Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Wed, 16 May 2012 11:13:58 +0100
Subject: [PATCH 04/16] GIT: checkCvMesh: surpassed by checkMesh
 -meshQualityDict option

---
 .../generation/cvMesh/checkCvMesh/Make/files  |   3 -
 .../cvMesh/checkCvMesh/Make/options           |  14 --
 .../cvMesh/checkCvMesh/checkCvMesh.C          | 122 ------------------
 3 files changed, 139 deletions(-)
 delete mode 100644 applications/utilities/mesh/generation/cvMesh/checkCvMesh/Make/files
 delete mode 100644 applications/utilities/mesh/generation/cvMesh/checkCvMesh/Make/options
 delete mode 100644 applications/utilities/mesh/generation/cvMesh/checkCvMesh/checkCvMesh.C

diff --git a/applications/utilities/mesh/generation/cvMesh/checkCvMesh/Make/files b/applications/utilities/mesh/generation/cvMesh/checkCvMesh/Make/files
deleted file mode 100644
index 73871a7532..0000000000
--- a/applications/utilities/mesh/generation/cvMesh/checkCvMesh/Make/files
+++ /dev/null
@@ -1,3 +0,0 @@
-checkCvMesh.C
-
-EXE = $(FOAM_APPBIN)/checkCvMesh
diff --git a/applications/utilities/mesh/generation/cvMesh/checkCvMesh/Make/options b/applications/utilities/mesh/generation/cvMesh/checkCvMesh/Make/options
deleted file mode 100644
index ba68fd3819..0000000000
--- a/applications/utilities/mesh/generation/cvMesh/checkCvMesh/Make/options
+++ /dev/null
@@ -1,14 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/triSurface/lnInclude \
-    -I$(LIB_SRC)/mesh/autoMesh/lnInclude \
-    -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
-
-EXE_LIBS = \
-    -lfiniteVolume \
-    -ldynamicMesh \
-    -ltriSurface \
-    -lautoMesh \
-    -lmeshTools
-
diff --git a/applications/utilities/mesh/generation/cvMesh/checkCvMesh/checkCvMesh.C b/applications/utilities/mesh/generation/cvMesh/checkCvMesh/checkCvMesh.C
deleted file mode 100644
index 88695dbd39..0000000000
--- a/applications/utilities/mesh/generation/cvMesh/checkCvMesh/checkCvMesh.C
+++ /dev/null
@@ -1,122 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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/>.
-
-Application
-    checkCvMesh
-
-Description
-
-\*---------------------------------------------------------------------------*/
-
-#include "argList.H"
-#include "Time.H"
-#include "fvMesh.H"
-#include "autoSnapDriver.H"
-#include "faceSet.H"
-#include "motionSmoother.H"
-#include "timeSelector.H"
-
-
-using namespace Foam;
-
-
-int main(int argc, char *argv[])
-{
-    timeSelector::addOptions();
-
-#   include "addOverwriteOption.H"
-
-#   include "setRootCase.H"
-#   include "createTime.H"
-
-    instantList timeDirs = timeSelector::select0(runTime, args);
-
-#   include "createNamedPolyMesh.H"
-
-    runTime.functionObjects().off();
-
-    forAll(timeDirs, timeI)
-    {
-        runTime.setTime(timeDirs[timeI], timeI);
-
-        Info<< "Time = " << runTime.timeName()
-            << nl << endl;
-
-        mesh.readUpdate();
-
-        // Check patches and faceZones are synchronised
-        mesh.boundaryMesh().checkParallelSync(true);
-        meshRefinement::checkCoupledFaceZones(mesh);
-
-        // Read meshing dictionary
-        IOdictionary cvMeshDict
-        (
-           IOobject
-           (
-                "cvMeshDict",
-                runTime.system(),
-                mesh,
-                IOobject::MUST_READ_IF_MODIFIED,
-                IOobject::NO_WRITE
-           )
-        );
-
-        // mesh motion and mesh quality parameters
-        const dictionary& meshQualityDict
-            = cvMeshDict.subDict("meshQualityControls");
-
-
-        Info<< "Checking mesh ..." << endl;
-
-        faceSet wrongFaces(mesh, "wrongFaces", mesh.nFaces()/100);
-
-        motionSmoother::checkMesh(false, mesh, meshQualityDict, wrongFaces);
-
-        const label nInitErrors = returnReduce
-        (
-            wrongFaces.size(),
-            sumOp<label>()
-        );
-
-        Info<< "Detected " << nInitErrors << " illegal faces"
-            << " (concave, zero area or negative cell pyramid volume)"
-            << endl;
-
-        if (nInitErrors > 0)
-        {
-            Info<< "Writing " << nInitErrors
-                << " faces in error to set "
-                << wrongFaces.name() << endl;
-
-            wrongFaces.instance() = mesh.pointsInstance();
-            wrongFaces.write();
-        }
-
-        Info<< nl << "End of time " << runTime.timeName() << nl << endl;
-    }
-
-    Info<< "End\n" << endl;
-
-    return 0;
-}
-

From b1d2767282132ad6b3a469368bfe58ddd269bb6c Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Wed, 16 May 2012 11:16:46 +0100
Subject: [PATCH 05/16] ENH: collapseEdges: read meshQualityDict. parallel
 operation.

---
 .../advanced/collapseEdges/collapseEdges.C    |  2 +-
 .../advanced/collapseEdges/meshQualityDict    | 75 +++++++++++++++++++
 2 files changed, 76 insertions(+), 1 deletion(-)
 create mode 100644 applications/utilities/mesh/advanced/collapseEdges/meshQualityDict

diff --git a/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C b/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C
index f970c1d836..f1d61f14a7 100644
--- a/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C
+++ b/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C
@@ -993,7 +993,7 @@ PackedBoolList checkMeshQuality(const polyMesh& mesh)
     (
         IOobject
         (
-            "meshQualityControls",
+            "meshQualityDict",
             mesh.time().system(),
             mesh.time(),
             IOobject::MUST_READ,
diff --git a/applications/utilities/mesh/advanced/collapseEdges/meshQualityDict b/applications/utilities/mesh/advanced/collapseEdges/meshQualityDict
new file mode 100644
index 0000000000..49810f468b
--- /dev/null
+++ b/applications/utilities/mesh/advanced/collapseEdges/meshQualityDict
@@ -0,0 +1,75 @@
+/*--------------------------------*- 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          meshQualityDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//- Maximum non-orthogonality allowed. Set to 180 to disable.
+maxNonOrtho         65;
+
+//- Max skewness allowed. Set to <0 to disable.
+maxBoundarySkewness 50;
+
+//- Max skewness allowed. Set to <0 to disable.
+maxInternalSkewness 10;
+
+//- Max concaveness allowed. Is angle (in degrees) below which concavity
+//  is allowed. 0 is straight face, <0 would be convex face.
+//  Set to 180 to disable.
+maxConcave          80;
+
+//- Minimum pyramid volume. Is absolute volume of cell pyramid.
+//  Set to a sensible fraction of the smallest cell volume expected.
+//  Set to very negative number (e.g. -1E30) to disable.
+minVol              1e-20;
+
+//- Minimum quality of the tet formed by the face-centre
+//  and variable base point minimum decomposition triangles and
+//  the cell centre. This has to be a positive number for tracking
+//  to work. Set to very negative number (e.g. -1E30) to
+//  disable.
+//     <0 = inside out tet,
+//      0 = flat tet
+//      1 = regular tet
+minTetQuality       1e-30;
+
+//- Minimum face area. Set to <0 to disable.
+minArea             -1;
+
+//- Minimum face twist. Set to <-1 to disable. dot product of face normal
+//- and face centre triangles normal
+minTwist            0.0;
+
+//- minimum normalised cell determinant
+//- 1 = hex, <= 0 = folded or flattened illegal cell
+minDeterminant      0.001;
+
+//- minFaceWeight (0 -> 0.5)
+minFaceWeight       0.02;
+
+//- minVolRatio (0 -> 1)
+minVolRatio         0.01;
+
+//must be >0 for Fluent compatibility
+minTriangleTwist    -1;
+
+
+// ************************************************************************* //

From 30f6574ceeea4270e33fe01671de77e29a5c96ba Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Wed, 16 May 2012 14:15:10 +0100
Subject: [PATCH 06/16] ENH: vtkTopo.C: paraview3.12 handle wedge correctly

---
 .../dataConversion/foamToVTK/vtkTopo.C        | 35 ++++++++++---------
 1 file changed, 18 insertions(+), 17 deletions(-)

diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/vtkTopo.C b/applications/utilities/postProcessing/dataConversion/foamToVTK/vtkTopo.C
index 6edcb09f98..f0ef8cd88a 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/vtkTopo.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/vtkTopo.C
@@ -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) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -52,6 +52,7 @@ Foam::vtkTopo::vtkTopo(const polyMesh& mesh)
     const cellModel& tet = *(cellModeller::lookup("tet"));
     const cellModel& pyr = *(cellModeller::lookup("pyr"));
     const cellModel& prism = *(cellModeller::lookup("prism"));
+    const cellModel& wedge = *(cellModeller::lookup("wedge"));
     const cellModel& tetWedge = *(cellModeller::lookup("tetWedge"));
     const cellModel& hex = *(cellModeller::lookup("hex"));
 
@@ -77,7 +78,7 @@ Foam::vtkTopo::vtkTopo(const polyMesh& mesh)
             if
             (
                 model != hex
-                // && model != wedge    // See above.
+             && model != wedge    // See above.
              && model != prism
              && model != pyr
              && model != tet
@@ -164,21 +165,21 @@ Foam::vtkTopo::vtkTopo(const polyMesh& mesh)
 
             cellTypes_[cellI] = VTK_WEDGE;
         }
-//        else if (cellModel == wedge)
-//        {
-//            // Treat as squeezed hex
-//            vtkVerts.setSize(8);
-//            vtkVerts[0] = cellShape[0];
-//            vtkVerts[1] = cellShape[1];
-//            vtkVerts[2] = cellShape[2];
-//            vtkVerts[3] = cellShape[0];
-//            vtkVerts[4] = cellShape[3];
-//            vtkVerts[5] = cellShape[4];
-//            vtkVerts[6] = cellShape[5];
-//            vtkVerts[7] = cellShape[6];
-//
-//            cellTypes_[cellI] = VTK_HEXAHEDRON;
-//        }
+        else if (cellModel == wedge)
+        {
+            // Treat as squeezed hex
+            vtkVerts.setSize(8);
+            vtkVerts[0] = cellShape[0];
+            vtkVerts[1] = cellShape[1];
+            vtkVerts[2] = cellShape[2];
+            vtkVerts[3] = cellShape[2];
+            vtkVerts[4] = cellShape[3];
+            vtkVerts[5] = cellShape[4];
+            vtkVerts[6] = cellShape[5];
+            vtkVerts[7] = cellShape[6];
+
+            cellTypes_[cellI] = VTK_HEXAHEDRON;
+        }
         else if (cellModel == hex)
         {
             vtkVerts = cellShape;

From 9b30d01b963593bd7cc7b7386c548fae84b44895 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Wed, 16 May 2012 14:43:26 +0100
Subject: [PATCH 07/16] ENH: snappyHexMesh: removed checking code since now in
 checkMesh

---
 .../generation/snappyHexMesh/snappyHexMesh.C  | 40 -------------------
 1 file changed, 40 deletions(-)

diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
index 40adb83035..279808bff1 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
@@ -123,14 +123,6 @@ void writeMesh
 int main(int argc, char *argv[])
 {
 #   include "addOverwriteOption.H"
-
-    Foam::argList::addBoolOption
-    (
-        "checkOnly",
-        "check existing mesh against snappyHexMeshDict settings"
-    );
-
-
 #   include "setRootCase.H"
 #   include "createTime.H"
     runTime.functionObjects().off();
@@ -141,8 +133,6 @@ int main(int argc, char *argv[])
 
     const bool overwrite = args.optionFound("overwrite");
 
-    const bool checkOnly = args.optionFound("checkOnly");
-
     // Check patches and faceZones are synchronised
     mesh.boundaryMesh().checkParallelSync(true);
     meshRefinement::checkCoupledFaceZones(mesh);
@@ -184,36 +174,6 @@ int main(int argc, char *argv[])
     );
 
 
-    if (checkOnly)
-    {
-        Info<< "Checking initial mesh ..." << endl;
-        faceSet wrongFaces(mesh, "wrongFaces", mesh.nFaces()/100);
-        motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces);
-
-        const label nInitErrors = returnReduce
-        (
-            wrongFaces.size(),
-            sumOp<label>()
-        );
-
-        Info<< "Detected " << nInitErrors << " illegal faces"
-            << " (concave, zero area or negative cell pyramid volume)"
-            << endl;
-
-        if (nInitErrors > 0)
-        {
-            Info<< "Writing " << nInitErrors
-                << " faces in error to set "
-                << wrongFaces.name() << endl;
-            wrongFaces.instance() = mesh.pointsInstance();
-            wrongFaces.write();
-        }
-
-        Info<< "End\n" << endl;
-
-        return 0;
-    }
-
 
     // Read decomposePar dictionary
     IOdictionary decomposeDict

From cf89f77c83cdeec12aa107c976e36d1383e94bd0 Mon Sep 17 00:00:00 2001
From: laurence <laurence>
Date: Wed, 16 May 2012 15:35:25 +0100
Subject: [PATCH 08/16] BUG: correct path in processorPolyPatch dump when faces
 do not match

---
 .../polyPatches/constraint/processor/processorPolyPatch.C       | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C
index e9faf44971..be7f1cb3fa 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C
@@ -634,7 +634,7 @@ bool Foam::processorPolyPatch::order
             fileName str
             (
                 boundaryMesh().mesh().time().path()
-               /name()/name()+"_faces.obj"
+               /name() + "_faces.obj"
             );
             Pout<< "processorPolyPatch::order :"
                 << " Writing faces to OBJ file " << str.name() << endl;

From abbbe31f41b2fcbbe69e252e34dc7c54d9fe38d6 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 17 May 2012 11:58:35 +0100
Subject: [PATCH 09/16] ENH: collapseEdges: parallel operation of mergeEdges

---
 .../advanced/collapseEdges/collapseEdges.C    | 83 +++++++++++--------
 1 file changed, 47 insertions(+), 36 deletions(-)

diff --git a/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C b/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C
index f1d61f14a7..5fa068f730 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) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -50,6 +50,7 @@ Description
 #include "PointEdgeWave.H"
 #include "pointEdgeCollapse.H"
 #include "motionSmoother.H"
+#include "removePoints.H"
 
 using namespace Foam;
 
@@ -788,51 +789,61 @@ label mergeEdges
     labelList& collapseEdge
 )
 {
-    const pointField& points = mesh.points();
     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;
 
-    forAll(pointEdges, pointI)
+    if (nTotRemove > 0)
     {
-        const labelList& pEdges = pointEdges[pointI];
-
-        if (pEdges.size() == 2 && boundaryPoint[pointI] <= 0)
+        forAll(pointEdges, pointI)
         {
-            // Collapse only if none of the points part of merge network and
-            // none protected (minEdgeLen < 0)
-            label e0 = pEdges[0];
-            label e1 = pEdges[1];
-
-            if
-            (
-                collapseEdge[e0] == -1
-             && minEdgeLen[e0] >= 0
-             && collapseEdge[e1] == -1
-             && minEdgeLen[e1] >= 0
-            )
+            if (pointCanBeDeleted[pointI])
             {
-                const edge& leftE = edges[e0];
-                const edge& rightE = edges[e1];
+                const labelList& pEdges = pointEdges[pointI];
 
-                // Get the two vertices on both sides of the point
-                label leftV = leftE.otherVertex(pointI);
-                label rightV = rightE.otherVertex(pointI);
-
-                // 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)
+                if (pEdges.size() == 2)
                 {
-                    // Collapse one (left) side of the edge. Make pointI
-                    // the master.
-                    collapseEdge[e0] = findIndex(leftE, pointI);
-                    nCollapsed++;
+                    // 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);
+                        }
+                    }
                 }
             }
         }

From 142c28d4f90191c6bd88a42aa6925071e9f751b8 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 17 May 2012 12:30:48 +0100
Subject: [PATCH 10/16] GIT: meshQualityDict: put under git

---
 .../mesh/cvMesh/blob/system/meshQualityDict   | 73 +++++++++++++++++++
 1 file changed, 73 insertions(+)
 create mode 100644 tutorials/mesh/cvMesh/blob/system/meshQualityDict

diff --git a/tutorials/mesh/cvMesh/blob/system/meshQualityDict b/tutorials/mesh/cvMesh/blob/system/meshQualityDict
new file mode 100644
index 0000000000..fa5319e087
--- /dev/null
+++ b/tutorials/mesh/cvMesh/blob/system/meshQualityDict
@@ -0,0 +1,73 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+FoamFile
+{
+    version         2.0;
+    format          ascii;
+
+    root            "";
+    case            "";
+    instance        "";
+    local           "";
+
+    class           dictionary;
+    object          meshQualityDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//- Maximum non-orthogonality allowed. Set to 180 to disable.
+maxNonOrtho         65;
+
+//- Max skewness allowed. Set to <0 to disable.
+maxBoundarySkewness 50;
+maxInternalSkewness 10;
+
+//- Max concaveness allowed. Is angle (in degrees) below which concavity
+//  is allowed. 0 is straight face, <0 would be convex face.
+//  Set to 180 to disable.
+maxConcave          80;
+
+//- Minimum quality of the tet formed by the face-centre
+//  and variable base point minimum decomposition triangles and
+//  the cell centre. This has to be a positive number for tracking
+//  to work. Set to very negative number (e.g. -1E30) to
+//  disable.
+//     <0 = inside out tet,
+//      0 = flat tet
+//      1 = regular tet
+minTetQuality       1e-30;
+
+//- Minimum pyramid volume. Is absolute volume of cell pyramid.
+//  Set to a sensible fraction of the smallest cell volume expected.
+//  Set to very negative number (e.g. -1E30) to disable.
+minVol              1e-20;
+
+//- Minimum face area. Set to <0 to disable.
+minArea             -1;
+
+//- Minimum face twist. Set to <-1 to disable. dot product of face normal
+//- and face centre triangles normal
+minTwist            0.001;
+
+//- minimum normalised cell determinant
+//- 1 = hex, <= 0 = folded or flattened illegal cell
+minDeterminant      0.001;
+
+//- minFaceWeight (0 -> 0.5)
+minFaceWeight       0.02;
+
+//- minVolRatio (0 -> 1)
+minVolRatio         0.01;
+
+//must be >0 for Fluent compatibility
+minTriangleTwist    -1;
+
+
+// ************************************************************************* //

From ba7aff0333954d313e436ca56b515b28d15c0007 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 17 May 2012 12:32:05 +0100
Subject: [PATCH 11/16] STYLE: modifyMeshDict: improved comment

---
 .../mesh/advanced/modifyMesh/modifyMeshDict     | 17 +++++++++--------
 1 file changed, 9 insertions(+), 8 deletions(-)

diff --git a/applications/utilities/mesh/advanced/modifyMesh/modifyMeshDict b/applications/utilities/mesh/advanced/modifyMesh/modifyMeshDict
index cd2de0b8eb..b388da5dac 100644
--- a/applications/utilities/mesh/advanced/modifyMesh/modifyMeshDict
+++ b/applications/utilities/mesh/advanced/modifyMesh/modifyMeshDict
@@ -14,7 +14,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// Move points:
+// Move boundary points:
 // Every entry is two coordinates. First one is location of the point to move,
 // the second is the position to move to.
 pointsToMove
@@ -22,15 +22,15 @@ pointsToMove
     (( -0.17861 -0.45073 0.75276)( -0.18 -0.45073 0.75276))
 );
 
-// Split edge in two:
-// First coord is a point on the edge to cut, second is the position of the
-// newly introduced point
+// Split boundary edge in two:
+// First coord is a point on the (boundary) edge to cut, second is the
+// position of the newly introduced point
 edgesToSplit
 (
     (( -0.17692 -0.45312 0.74516)( -0.18 -0.45 0.742))
 );
 
-// Triangulate a face:
+// Triangulate a boundary face:
 // First coord is a point on the face to triangulate. It will introduce a
 // point on the face, triangulate and move the point to the second coordinate.
 facesToTriangulate
@@ -38,8 +38,8 @@ facesToTriangulate
     (( -0.039123 -0.45045 0.74083) (-0.03844 -0.45049 0.73572))
 );
 
-// Edges to collapse. First coord is point on the edge, second is coordinate
-// to collapse to.
+// Boundary edges to collapse. First coord is point on the edge, second
+// is coordinate to collapse to.
 edgesToCollapse
 (
     ((0.054975 0.099987 0.0044074)(0.054975 0.099987 0.0044074))
@@ -55,7 +55,8 @@ cellsToSplit
 );
 
 // Change patch:
-// Changes patchID of faces. Coord selects the face, label is the patch index.
+// Changes patchID of boundary faces. Coord selects the face, label is the
+// patch index.
 facesToRepatch
 (
     (( -0.039123 -0.45045 0.74083) 1)

From 61f9b0ec48d8777fbcd732eea445165bdb39b66f Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 17 May 2012 12:32:40 +0100
Subject: [PATCH 12/16] ENH: removePoints: updated comment

---
 .../polyTopoChange/polyTopoChange/removePoints.H            | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.H b/src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.H
index 10e31d034e..4a2bfcfc07 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.H
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.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) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -127,9 +127,7 @@ public:
             //  A point can be deleted if
             //  - it is not used by any edge.
             //  or
-            //  - is not used by an internal edge
-            //  - is used by only two boundary edges. (note that these two
-            //    edges will always be boundary ones!)
+            //  - is used by only two edges
             //  - these two edges are sufficiently in line (cos > minCos)
             //  - all processors agree that point can be deleted.
             label countPointUsage

From 255b1c4ab439735f7ef79e22976b7c6fcb4fd34d Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 17 May 2012 12:58:47 +0100
Subject: [PATCH 13/16] ENH: blob: updated tutorial

---
 tutorials/mesh/cvMesh/blob/Allrun             |  3 +-
 tutorials/mesh/cvMesh/blob/system/cvMeshDict  | 16 ++-----
 .../mesh/cvMesh/blob/system/snappyHexMeshDict | 48 +------------------
 3 files changed, 7 insertions(+), 60 deletions(-)

diff --git a/tutorials/mesh/cvMesh/blob/Allrun b/tutorials/mesh/cvMesh/blob/Allrun
index aafdcec484..860857ec95 100755
--- a/tutorials/mesh/cvMesh/blob/Allrun
+++ b/tutorials/mesh/cvMesh/blob/Allrun
@@ -6,7 +6,8 @@ cd ${0%/*} || exit 1    # run from this directory
 
 runApplication blockMesh
 runApplication cvMesh
-runApplication snappyHexMesh
+#runApplication collapseEdges 1e-10 180 -overwrite
+runApplication snappyHexMesh -overwrite
 runApplication checkMesh -constant -allGeometry -allTopology
 
 # ----------------------------------------------------------------- end-of-file
diff --git a/tutorials/mesh/cvMesh/blob/system/cvMeshDict b/tutorials/mesh/cvMesh/blob/system/cvMeshDict
index ca2fd3ccd6..27fff06883 100644
--- a/tutorials/mesh/cvMesh/blob/system/cvMeshDict
+++ b/tutorials/mesh/cvMesh/blob/system/cvMeshDict
@@ -92,6 +92,8 @@ surfaceConformation
 
     surfaceConformationRebuildFrequency 10;
 
+    specialiseFeaturePoints off;
+
     conformationControls
     {
         initial
@@ -230,18 +232,8 @@ polyMeshFiltering
 
 meshQualityControls
 {
-    maxNonOrtho         65;
-    maxBoundarySkewness 50;
-    maxInternalSkewness 10;
-    maxConcave          80;
-    minVol              1e-20;
-    minTetQuality       1e-30;
-    minArea             -1;
-    minTwist            0.0;
-    minDeterminant      0.001;
-    minFaceWeight       0.02;
-    minVolRatio         0.01;
-    minTriangleTwist    -1;
+    #include "meshQualityDict"
 }
 
+
 // ************************************************************************* //
diff --git a/tutorials/mesh/cvMesh/blob/system/snappyHexMeshDict b/tutorials/mesh/cvMesh/blob/system/snappyHexMeshDict
index 3ac70912b4..da3f1895b4 100644
--- a/tutorials/mesh/cvMesh/blob/system/snappyHexMeshDict
+++ b/tutorials/mesh/cvMesh/blob/system/snappyHexMeshDict
@@ -282,53 +282,7 @@ addLayersControls
 // where to undo.
 meshQualityControls
 {
-    //- Maximum non-orthogonality allowed. Set to 180 to disable.
-    maxNonOrtho 65;
-
-    //- Max skewness allowed. Set to <0 to disable.
-    maxBoundarySkewness 20;
-    maxInternalSkewness 4;
-
-    //- Max concaveness allowed. Is angle (in degrees) below which concavity
-    //  is allowed. 0 is straight face, <0 would be convex face.
-    //  Set to 180 to disable.
-    maxConcave 80;
-
-    //- Minimum pyramid volume. Is absolute volume of cell pyramid.
-    //  Set to a sensible fraction of the smallest cell volume expected.
-    //  Set to very negative number (e.g. -1E30) to disable.
-    minVol 1e-13;
-
-    //- Minimum quality of the tet formed by the face-centre
-    //  and variable base point minimum decomposition triangles and
-    //  the cell centre. This has to be a positive number for tracking
-    //  to work. Set to very negative number (e.g. -1E30) to
-    //  disable.
-    //     <0 = inside out tet,
-    //      0 = flat tet
-    //      1 = regular tet
-    minTetQuality -1;   //1e-30;
-
-    //- Minimum face area. Set to <0 to disable.
-    minArea -1;
-
-    //- Minimum face twist. Set to <-1 to disable. dot product of face normal
-    //- and face centre triangles normal
-    minTwist 0.02;
-
-    //- minimum normalised cell determinant
-    //- 1 = hex, <= 0 = folded or flattened illegal cell
-    minDeterminant 0.001;
-
-    //- minFaceWeight (0 -> 0.5)
-    minFaceWeight 0.02;
-
-    //- minVolRatio (0 -> 1)
-    minVolRatio 0.01;
-
-    //must be >0 for Fluent compatibility
-    minTriangleTwist -1;
-
+    #include "meshQualityDict"
 
     // Advanced
 

From 3f8aefc4707c28c7b9ccfd57e8c87821d1813eec Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Fri, 18 May 2012 12:16:01 +0100
Subject: [PATCH 14/16] ENH: mapDistribute: added wrapper for DynamicList

---
 .../mapPolyMesh/mapDistribute/mapDistribute.H | 11 +++++++++-
 .../mapDistribute/mapDistributeTemplates.C    | 21 ++++++++++++++++++-
 2 files changed, 30 insertions(+), 2 deletions(-)

diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H
index 633a07d840..5c17f0547f 100644
--- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H
+++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.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) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -554,6 +554,15 @@ public:
                 const int tag = UPstream::msgType()
             );
 
+            //- Distribute data using default commsType.
+            template<class T>
+            void distribute
+            (
+                DynamicList<T>& fld,
+                const bool dummyTransform = true,
+                const int tag = UPstream::msgType()
+            ) const;
+
             //- Distribute data using default commsType.
             template<class T>
             void distribute
diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C
index 037075c96e..577af21ed1 100644
--- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C
+++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C
@@ -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) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -887,6 +887,25 @@ void Foam::mapDistribute::applyInverseTransforms
 }
 
 
+//- Distribute data using default commsType.
+template<class T>
+void Foam::mapDistribute::distribute
+(
+    DynamicList<T>& fld,
+    const bool dummyTransform,
+    const int tag
+) const
+{
+    fld.shrink();
+
+    List<T>& fldList = static_cast<List<T>& >(fld);
+
+    distribute(fldList, dummyTransform, tag);
+
+    fld.setCapacity(fldList.size());
+}
+
+
 //- Distribute data using default commsType.
 template<class T>
 void Foam::mapDistribute::distribute

From 3fc8c16b0d46addb865b15968a2cf2af4e5e0941 Mon Sep 17 00:00:00 2001
From: sergio <sergio>
Date: Tue, 22 May 2012 16:29:41 +0100
Subject: [PATCH 15/16] BUG: Fix to source in energy Eq for solid

---
 .../porousSolid/setPorousRegionSolidFields.H                   | 3 +++
 .../chtMultiRegionSimpleFoam/porousSolid/solvePorousSolid.H    | 2 +-
 .../porousSolid/setPorousRegionSolidFields.H                   | 3 +++
 .../chtMultiRegionFoam/porousSolid/solvePorousSolid.H          | 2 +-
 4 files changed, 8 insertions(+), 2 deletions(-)

diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/setPorousRegionSolidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/setPorousRegionSolidFields.H
index bbbd079ef2..3ce6b03b9d 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/setPorousRegionSolidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/setPorousRegionSolidFields.H
@@ -15,6 +15,9 @@
     const volScalarField& kappa = tkappa();
     //const volSymmTensorField& K = tK();
 
+    tmp<volScalarField> trhoCp = cp*rho;
+    const volScalarField& rhoCp = trhoCp();
+
     volScalarField& T = thermo.T();
 
     IObasicSourceList& sources = solidHeatSources[i];
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/solvePorousSolid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/solvePorousSolid.H
index 2dd09b974f..5d4bb47b22 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/solvePorousSolid.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/solvePorousSolid.H
@@ -4,7 +4,7 @@
         tmp<fvScalarMatrix> TEqn
         (
           - fvm::laplacian(betav*kappa, T, "laplacian(K,T)")
-          + sources(rho, T)
+          + sources(rhoCp, T)
         );
 
         TEqn().relax();
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/setPorousRegionSolidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/setPorousRegionSolidFields.H
index 6f4834000d..2ade49d0ac 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/setPorousRegionSolidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/setPorousRegionSolidFields.H
@@ -15,6 +15,9 @@
     const volScalarField& kappa = tkappa();
     //const volSymmTensorField& K = tK();
 
+    tmp<volScalarField> trhoCp = cp*rho;
+    const volScalarField& rhoCp = trhoCp();
+
     volScalarField& T = thermo.T();
 
     IObasicSourceList& sources = solidHeatSources[i];
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/solvePorousSolid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/solvePorousSolid.H
index e8941cf2ad..dfe0bff675 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/solvePorousSolid.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/solvePorousSolid.H
@@ -10,7 +10,7 @@ if (finalIter)
         (
             fvm::ddt(betav*rho*cp, T)
           - fvm::laplacian(betav*kappa, T, "laplacian(K,T)")
-          + sources(rho, T)
+          + sources(rhoCp, T)
         );
 
         TEqn().relax();

From f22ba88d8ae2845a27edecce2f9ec1952627e519 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 22 May 2012 17:10:26 +0100
Subject: [PATCH 16/16] ENH: extrudeMesh: allow flipNormals for
 extrude-from-patch

---
 .../extrude/extrudeMesh/extrudeMesh.C         | 63 ++++++++++++++++---
 .../extrude/extrudeMesh/extrudeMeshDict       |  3 +-
 2 files changed, 58 insertions(+), 8 deletions(-)

diff --git a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C
index 1c769a4346..6cd54bb461 100644
--- a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C
+++ b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C
@@ -251,10 +251,10 @@ int main(int argc, char *argv[])
 
     if (mode == PATCH || mode == MESH)
     {
-        if (flipNormals)
+        if (flipNormals && mode == MESH)
         {
             FatalErrorIn(args.executable())
-                << "Flipping normals not supported for extrusions from patch."
+                << "Flipping normals not supported for extrusions from mesh."
                 << exit(FatalError);
         }
 
@@ -297,6 +297,60 @@ int main(int argc, char *argv[])
         addPatchCellLayer layerExtrude(mesh, (mode == MESH));
 
         const labelList meshFaces(patchFaces(patches, sourcePatches));
+
+        if (mode == PATCH && flipNormals)
+        {
+            // Cheat. Flip patch faces in mesh. This invalidates the
+            // mesh (open cells) but does produce the correct extrusion.
+            polyTopoChange meshMod(mesh);
+            forAll(meshFaces, i)
+            {
+                label meshFaceI = meshFaces[i];
+
+                label patchI = patches.whichPatch(meshFaceI);
+                label own = mesh.faceOwner()[meshFaceI];
+                label nei = -1;
+                if (patchI == -1)
+                {
+                    nei = mesh.faceNeighbour()[meshFaceI];
+                }
+
+                label zoneI = mesh.faceZones().whichZone(meshFaceI);
+                bool zoneFlip = false;
+                if (zoneI != -1)
+                {
+                    label index = mesh.faceZones()[zoneI].whichFace(meshFaceI);
+                    zoneFlip = mesh.faceZones()[zoneI].flipMap()[index];
+                }
+
+                meshMod.modifyFace
+                (
+                    mesh.faces()[meshFaceI].reverseFace(),  // modified face
+                    meshFaceI,                      // label of face
+                    own,                            // owner
+                    nei,                            // neighbour
+                    true,                           // face flip
+                    patchI,                         // patch for face
+                    zoneI,                          // zone for face
+                    zoneFlip                        // face flip in zone
+                );
+            }
+
+            // Change the mesh. No inflation.
+            autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);
+
+            // Update fields
+            mesh.updateMesh(map);
+
+            // Move mesh (since morphing does not do this)
+            if (map().hasMotionPoints())
+            {
+                mesh.movePoints(map().preMotionPoints());
+            }
+        }
+
+
+
         indirectPrimitivePatch extrudePatch
         (
             IndirectList<face>
@@ -471,11 +525,6 @@ int main(int argc, char *argv[])
             displacement[pointI] = extrudePt - layer0Points[pointI];
         }
 
-        if (flipNormals)
-        {
-            Info<< "Flipping faces." << nl << endl;
-            displacement = -displacement;
-        }
 
         // Check if wedge (has layer0 different from original patch points)
         // If so move the mesh to starting position.
diff --git a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMeshDict b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMeshDict
index f1985fe7f6..e075acbc4e 100644
--- a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMeshDict
+++ b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMeshDict
@@ -31,7 +31,8 @@ exposedPatchName movingWall;
 // If construct from surface:
 surface "movingWall.stl";
 
-// Flip surface normals before usage.
+// Flip surface normals before usage. Valid only for extrude from surface or
+// patch.
 flipNormals false;
 
 //- Linear extrusion in point-normal direction