From 3715ed220747c4e47b11456f0ac7681e3e2b1396 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@Germany>
Date: Tue, 30 Mar 2010 14:07:51 +0200
Subject: [PATCH 01/20] STYLE: had OpenCFD instead of OpenFOAM in COPYING

---
 COPYING | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/COPYING b/COPYING
index 81ee282d31..7e02780f0b 100644
--- a/COPYING
+++ b/COPYING
@@ -3,7 +3,7 @@
  OpenFOAM(R) is Copyright (C) 1991-2010 OpenCFD Ltd.
  Contact: OpenCFD (enquiries@OpenCFD.co.uk)
 
- You may use, distribute and copy the OpenCFD CFD Toolbox under the terms
+ You may use, distribute and copy the OpenFOAM CFD Toolbox under the terms
  of GNU General Public License version 3, which is displayed below, or
  (at your option) any later version.
 

From 883dbc2e2dd8eda98da439cedfb0077c85506e6c Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 30 Mar 2010 21:48:57 +0100
Subject: [PATCH 02/20] ENH: Add concave cell checking

---
 .../manipulation/checkMesh/checkGeometry.C    |  16 +++
 .../meshes/primitiveMesh/primitiveMesh.H      |  21 ++-
 .../primitiveMeshCheck/primitiveMeshCheck.C   | 135 +++++++++++++++++-
 3 files changed, 160 insertions(+), 12 deletions(-)

diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
index e16d90bbbd..50ff6918b0 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
@@ -291,6 +291,22 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
         }
     }
 
+    if (allGeometry)
+    {
+        cellSet cells(mesh, "concaveCells", mesh.nCells()/100);
+        if (mesh.checkConcaveCells(true, &cells))
+        {
+            noFailedChecks++;
+
+            label nCells = returnReduce(cells.size(), sumOp<label>());
+
+            Info<< "  <<Writing " << nCells
+                << " concave cells to set " << cells.name() << endl;
+            cells.instance() = mesh.pointsInstance();
+            cells.write();
+        }
+    }
+
 
     return noFailedChecks;
 }
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
index 1e3e2293d9..dc4ddd33c5 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
@@ -8,10 +8,10 @@
 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 2 of the License, or (at your
-    option) any later version.
+    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
@@ -19,8 +19,7 @@ License
     for more details.
 
     You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
     Foam::primitiveMesh
@@ -292,6 +291,9 @@ class primitiveMesh
             //- Skewness warning threshold
             static scalar skewThreshold_;
 
+            //- Threshold where faces are considered coplanar
+            static scalar planarCosAngle_;
+
 
 protected:
 
@@ -647,6 +649,13 @@ public:
                     labelHashSet* setPtr = NULL
                 ) const;
 
+                //- Check for concave cells
+                bool checkConcaveCells
+                (
+                    const bool report = false,
+                    labelHashSet* setPtr = NULL
+                ) const;
+
 
             //- Check mesh topology for correctness.
             //  Returns false for no error.
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
index a7eab05ab9..30b051d1e9 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
@@ -8,10 +8,10 @@
 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 2 of the License, or (at your
-    option) any later version.
+    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
@@ -19,8 +19,7 @@ License
     for more details.
 
     You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 \*---------------------------------------------------------------------------*/
 
@@ -38,6 +37,7 @@ Foam::scalar Foam::primitiveMesh::closedThreshold_  = 1.0e-6;
 Foam::scalar Foam::primitiveMesh::aspectThreshold_  = 1000;
 Foam::scalar Foam::primitiveMesh::nonOrthThreshold_ = 70;    // deg
 Foam::scalar Foam::primitiveMesh::skewThreshold_    = 4;
+Foam::scalar Foam::primitiveMesh::planarCosAngle_   = 1.0e-6;
 
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
@@ -2094,6 +2094,129 @@ bool Foam::primitiveMesh::checkCellDeterminant
 }
 
 
+bool Foam::primitiveMesh::checkConcaveCells
+(
+    const bool report,
+    labelHashSet* setPtr
+) const
+{
+    if (debug)
+    {
+        Info<< "bool primitiveMesh::checkConcaveCells(const bool"
+            << ", labelHashSet*) const: "
+            << "checking for concave cells" << endl;
+    }
+
+    const cellList& c = cells();
+    const labelList& fOwner = faceOwner();
+    const vectorField& fAreas = faceAreas();
+    const pointField& fCentres = faceCentres();
+    const pointField& pts = points();
+
+    label nConcaveCells = 0;
+
+    forAll (c, cellI)
+    {
+        const cell& cFaces = c[cellI];
+
+        const labelList& cPoints = cellPoints()[cellI];
+
+        bool concave = false;
+
+        forAll(cFaces, i)
+        {
+            if (concave)
+            {
+                break;
+            }
+
+            label fI = cFaces[i];
+
+            const point& fC = fCentres[fI];
+
+            const face& f = faces()[fI];
+
+            vector fN = fAreas[fI];
+
+            fN /= max(mag(fN), VSMALL);
+
+            // Flip normal if required so that it is always pointing out of
+            // the cell
+            if (fOwner[fI] != cellI)
+            {
+                fN *= -1;
+            }
+
+            // Is any vertex of the cell on the wrong side of the
+            // plane of this face?
+
+            forAll(cPoints, cPtI)
+            {
+                label ptI = cPoints[cPtI];
+
+                // Skip points that are on this face
+                if (findIndex(f, ptI) > -1)
+                {
+                    continue;
+                }
+
+                const point& pt = pts[ptI];
+
+                // If the cell is concave, the point will be on the
+                // positive normal side of the plane of f, defined by
+                // its centre and normal, and the angle between (pt -
+                // fC) and fN will be less than 90 degrees, so the dot
+                // product will be positive.
+
+                vector pC = (pt - fC);
+
+                pC /= max(mag(pC), VSMALL);
+
+                if ((pC & fN) > -planarCosAngle_)
+                {
+                    // Concave or planar face
+
+                    concave = true;
+
+                    if (setPtr)
+                    {
+                        setPtr->insert(cellI);
+                    }
+
+                    nConcaveCells++;
+
+                    break;
+                }
+            }
+        }
+    }
+
+    reduce(nConcaveCells, sumOp<label>());
+
+    if (nConcaveCells > 0)
+    {
+        if (debug || report)
+        {
+            Info<< " ***Concave cells found, number of cells: "
+                << nConcaveCells << endl;
+        }
+
+        return true;
+    }
+    else
+    {
+        if (debug || report)
+        {
+            Info<< "    Concave cell check OK." << endl;
+        }
+
+        return false;
+    }
+
+    return false;
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 bool Foam::primitiveMesh::checkTopology(const bool report) const

From d8c1aa3afd45990ee7ca21fd46c44aad9342f121 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 30 Mar 2010 21:49:37 +0100
Subject: [PATCH 03/20] BUG: nGrow of layer disabling not handled correctly.

Now also moved to within layer iteration loop.
---
 .../autoHexMeshDriver/autoLayerDriver.C       | 51 +++++++++++++++----
 .../autoHexMeshDriver/autoLayerDriver.H       |  4 +-
 2 files changed, 42 insertions(+), 13 deletions(-)

diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
index 021568e461..4841b49b96 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
@@ -1376,7 +1376,7 @@ void Foam::autoLayerDriver::growNoExtrusion
     pointField& patchDisp,
     labelList& patchNLayers,
     List<extrudeMode>& extrudeStatus
-)
+) const
 {
     Info<< nl << "Growing non-extrusion points by one layer ..." << endl;
 
@@ -1407,7 +1407,7 @@ void Foam::autoLayerDriver::growNoExtrusion
             {
                 if
                 (
-                    extrudeStatus[f[fp]] == NOEXTRUDE
+                    extrudeStatus[f[fp]] == EXTRUDE
                  && grownExtrudeStatus[f[fp]] != NOEXTRUDE
                 )
                 {
@@ -1420,6 +1420,31 @@ void Foam::autoLayerDriver::growNoExtrusion
 
     extrudeStatus.transfer(grownExtrudeStatus);
 
+
+    // Synchronise since might get called multiple times.
+    // Use the fact that NOEXTRUDE is the minimum value.
+    {
+        labelList status(extrudeStatus.size());
+        forAll(status, i)
+        {
+            status[i] = extrudeStatus[i];
+        }
+        syncTools::syncPointList
+        (
+            meshRefiner_.mesh(),
+            pp.meshPoints(),
+            status,
+            minEqOp<label>(),
+            labelMax,           // null value
+            false               // no separation
+        );
+        forAll(status, i)
+        {
+            extrudeStatus[i] = extrudeMode(status[i]);
+        }
+    }
+
+
     forAll(extrudeStatus, patchPointI)
     {
         if (extrudeStatus[patchPointI] == NOEXTRUDE)
@@ -2712,8 +2737,6 @@ void Foam::autoLayerDriver::addLayers
 
             const labelList& meshPoints = patches[patchI].meshPoints();
 
-            //scalar maxThickness = -VGREAT;
-            //scalar minThickness = VGREAT;
             scalar sumThickness = 0;
             scalar sumNearWallThickness = 0;
 
@@ -2721,8 +2744,6 @@ void Foam::autoLayerDriver::addLayers
             {
                 label ppPointI = pp().meshPointMap()[meshPoints[patchPointI]];
 
-                //maxThickness = max(maxThickness, thickness[ppPointI]);
-                //minThickness = min(minThickness, thickness[ppPointI]);
                 sumThickness += thickness[ppPointI];
 
                 label nLay = patchNLayers[ppPointI];
@@ -2750,8 +2771,6 @@ void Foam::autoLayerDriver::addLayers
 
             if (totNPoints > 0)
             {
-                //reduce(maxThickness, maxOp<scalar>());
-                //reduce(minThickness, minOp<scalar>());
                 avgThickness =
                     returnReduce(sumThickness, sumOp<scalar>())
                   / totNPoints;
@@ -2767,8 +2786,6 @@ void Foam::autoLayerDriver::addLayers
                 << " " << setw(6) << layerParams.numLayers()[patchI]
                 << " " << setw(8) << avgNearWallThickness
                 << "  " << setw(8) << avgThickness
-                //<< " " << setw(8) << minThickness
-                //<< " " << setw(8) << maxThickness
                 << endl;
         }
         Info<< endl;
@@ -3148,6 +3165,19 @@ void Foam::autoLayerDriver::addLayers
         meshMover().movePoints(oldPoints);
         meshMover().correct();
 
+
+        // Grow out region of non-extrusion
+        for (label i = 0; i < layerParams.nGrow(); i++)
+        {
+            growNoExtrusion
+            (
+                pp,
+                patchDisp,
+                patchNLayers,
+                extrudeStatus
+            );
+        }
+
         Info<< endl;
     }
 
@@ -3294,7 +3324,6 @@ void Foam::autoLayerDriver::doLayers
 
 
         // Balance
-        if (Pstream::parRun())
         if (Pstream::parRun() && preBalance)
         {
             Info<< nl
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
index 8f1e27334b..5061a82d8e 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
@@ -233,13 +233,13 @@ class autoLayerDriver
                 ) const;
 
                 //- Grow no-extrusion layer.
-                static void growNoExtrusion
+                void growNoExtrusion
                 (
                     const indirectPrimitivePatch& pp,
                     pointField& patchDisp,
                     labelList& patchNLayers,
                     List<extrudeMode>& extrudeStatus
-                );
+                ) const;
 
                 //- Calculate pointwise wanted and minimum thickness.
                 //  thickness: wanted thickness

From 0b1fc8e0cd3dce01eac990b8dfa9695f559e6c25 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 30 Mar 2010 21:50:40 +0100
Subject: [PATCH 04/20] ENH: Changed behaviour of surfaces with zones.

- if allowFreeStanding = false it also disables faceZones on boundaries.
- can have surfaces with faceZone but no cellZone.
---
 .../meshRefinement/meshRefinementBaffles.C    | 97 +++++++++----------
 .../refinementSurfaces/refinementSurfaces.C   | 46 ++++++---
 .../refinementSurfaces/refinementSurfaces.H   |  5 +-
 3 files changed, 85 insertions(+), 63 deletions(-)

diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
index 8d6627882c..1b2a6391d0 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -246,27 +246,8 @@ void Foam::meshRefinement::getBafflePatches
 
     const pointField& cellCentres = mesh_.cellCentres();
 
-    // Build list of surfaces that are not to be baffled.
-    const wordList& cellZoneNames = surfaces_.cellZoneNames();
-
-    labelList surfacesToBaffle(cellZoneNames.size());
-    label baffleI = 0;
-    forAll(cellZoneNames, surfI)
-    {
-        if (cellZoneNames[surfI].size())
-        {
-            if (debug)
-            {
-                Pout<< "getBafflePatches : Not baffling surface "
-                    << surfaces_.names()[surfI] << endl;
-            }
-        }
-        else
-        {
-            surfacesToBaffle[baffleI++] = surfI;
-        }
-    }
-    surfacesToBaffle.setSize(baffleI);
+    // Surfaces that need to be baffled
+    const labelList surfacesToBaffle(surfaces_.getUnnamedSurfaces());
 
     ownPatch.setSize(mesh_.nFaces());
     ownPatch = -1;
@@ -1254,7 +1235,7 @@ void Foam::meshRefinement::findCellZoneTopo
         {
             label surfI = namedSurfaceIndex[faceI];
 
-            if (surfI != -1)
+            if (surfI != -1 && surfaceToCellZone[surfI] != -1)
             {
                 // Calculate region to zone from cellRegions on either side
                 // of internal face.
@@ -1306,7 +1287,7 @@ void Foam::meshRefinement::findCellZoneTopo
 
                     label surfI = namedSurfaceIndex[faceI];
 
-                    if (surfI != -1)
+                    if (surfI != -1 && surfaceToCellZone[surfI] != -1)
                     {
                         bool changedCell = calcRegionToZone
                         (
@@ -1440,6 +1421,15 @@ void Foam::meshRefinement::makeConsistentFaceIndex
                 }
             }
         }
+        else
+        {
+            // Unzonify boundary faces
+            forAll(pp, i)
+            {
+                label faceI = pp.start()+i;
+                namedSurfaceIndex[faceI] = -1;
+            }
+        }
     }
 }
 
@@ -2043,14 +2033,10 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
 
     labelList namedSurfaces(surfaces_.getNamedSurfaces());
 
-    boolList isNamedSurface(cellZoneNames.size(), false);
-
     forAll(namedSurfaces, i)
     {
         label surfI = namedSurfaces[i];
 
-        isNamedSurface[surfI] = true;
-
         Info<< "Surface : " << surfaces_.names()[surfI] << nl
             << "    faceZone : " << faceZoneNames[surfI] << nl
             << "    cellZone : " << cellZoneNames[surfI] << endl;
@@ -2126,31 +2112,34 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
         {
             label surfI = namedSurfaces[i];
 
-            label zoneI = cellZones.findZoneID(cellZoneNames[surfI]);
-
-            if (zoneI == -1)
+            if (cellZoneNames[surfI] != word::null)
             {
-                zoneI = cellZones.size();
-                cellZones.setSize(zoneI+1);
-                cellZones.set
-                (
-                    zoneI,
-                    new cellZone
+                label zoneI = cellZones.findZoneID(cellZoneNames[surfI]);
+
+                if (zoneI == -1)
+                {
+                    zoneI = cellZones.size();
+                    cellZones.setSize(zoneI+1);
+                    cellZones.set
                     (
-                        cellZoneNames[surfI],   //name
-                        labelList(0),           //addressing
-                        zoneI,                  //index
-                        cellZones               //cellZoneMesh
-                    )
-                );
-            }
+                        zoneI,
+                        new cellZone
+                        (
+                            cellZoneNames[surfI],   //name
+                            labelList(0),           //addressing
+                            zoneI,                  //index
+                            cellZones               //cellZoneMesh
+                        )
+                    );
+                }
 
-            if (debug)
-            {
-                Pout<< "Cells inside " << surfaces_.names()[surfI]
-                    << " will go into cellZone " << zoneI << endl;
+                if (debug)
+                {
+                    Pout<< "Cells inside " << surfaces_.names()[surfI]
+                        << " will go into cellZone " << zoneI << endl;
+                }
+                surfaceToCellZone[surfI] = zoneI;
             }
-            surfaceToCellZone[surfI] = zoneI;
         }
 
         // Check they are synced
@@ -2322,6 +2311,11 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
 
     if (closedNamedSurfaces.size())
     {
+        Info<< "Found " << closedNamedSurfaces.size()
+            << " closed, named surfaces. Assigning cells in/outside"
+            << " these surfaces to the corresponding cellZone."
+            << nl << endl;
+
         findCellZoneGeometric
         (
             closedNamedSurfaces,   // indices of closed surfaces
@@ -2334,8 +2328,10 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
     // Set using walking
     // ~~~~~~~~~~~~~~~~~
 
-    //if (returnReduce(nSet, sumOp<label>()) < mesh_.globalData().nTotalCells())
+    if (!allowFreeStandingZoneFaces)
     {
+        Info<< "Walking to assign cellZones." << nl << endl;
+
         // Topological walk
         findCellZoneTopo
         (
@@ -2350,6 +2346,9 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
     // Make sure namedSurfaceIndex is unset inbetween same cell cell zones.
     if (!allowFreeStandingZoneFaces)
     {
+        Info<< "Only keeping zone faces inbetween different cellZones."
+            << nl << endl;
+
         makeConsistentFaceIndex(cellToZone, namedSurfaceIndex);
     }
 
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
index a23c1ef56c..2b1d12c6da 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
+++ b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
@@ -78,8 +78,18 @@ Foam::refinementSurfaces::refinementSurfaces
         if (dict.found("faceZone"))
         {
             dict.lookup("faceZone") >> faceZoneNames_[surfI];
-            dict.lookup("cellZone") >> cellZoneNames_[surfI];
-            dict.lookup("zoneInside") >> zoneInside_[surfI];
+            if (dict.readIfPresent("cellZone", cellZoneNames_[surfI]))
+            {
+                dict.lookup("zoneInside") >> zoneInside_[surfI];
+            }
+            else if (dict.found("zoneInside"))
+            {
+                IOWarningIn("refinementSurfaces::refinementSurfaces(..)", dict)
+                    << "Unused entry zoneInside for faceZone "
+                    << faceZoneNames_[surfI]
+                    << " since no cellZone specified."
+                    << endl;
+            }
         }
 
         // Global perpendicular angle
@@ -315,8 +325,21 @@ Foam::refinementSurfaces::refinementSurfaces
             if (dict.found("faceZone"))
             {
                 dict.lookup("faceZone") >> faceZoneNames_[surfI];
-                dict.lookup("cellZone") >> cellZoneNames_[surfI];
-                dict.lookup("zoneInside") >> zoneInside_[surfI];
+                if (dict.readIfPresent("cellZone", cellZoneNames_[surfI]))
+                {
+                    dict.lookup("zoneInside") >> zoneInside_[surfI];
+                }
+                else if (dict.found("zoneInside"))
+                {
+                    IOWarningIn
+                    (
+                        "refinementSurfaces::refinementSurfaces(..)",
+                        dict
+                    )   << "Unused entry zoneInside for faceZone "
+                        << faceZoneNames_[surfI]
+                        << " since no cellZone specified."
+                        << endl;
+                }
             }
 
             // Global perpendicular angle
@@ -476,18 +499,17 @@ Foam::labelList Foam::refinementSurfaces::getNamedSurfaces() const
 // Get indices of closed named surfaces
 Foam::labelList Foam::refinementSurfaces::getClosedNamedSurfaces() const
 {
-    labelList named(getNamedSurfaces());
+    labelList closed(cellZoneNames_.size());
 
-    labelList closed(named.size());
     label closedI = 0;
-
-    forAll(named, i)
+    forAll(cellZoneNames_, surfI)
     {
-        label surfI = named[i];
-
-        if (allGeometry_[surfaces_[surfI]].hasVolumeType())
+        if (cellZoneNames_[surfI].size())
         {
-            closed[closedI++] = surfI;
+            if (allGeometry_[surfaces_[surfI]].hasVolumeType())
+            {
+                closed[closedI++] = surfI;
+            }
         }
     }
     closed.setSize(closedI);
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
index b1188f1dff..b70d6252f9 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
+++ b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
@@ -147,7 +147,8 @@ public:
                 return faceZoneNames_;
             }
 
-            //- Per 'interface' surface : name of cellZone to put cells into
+            //- Per 'interface' surface : empty or name of cellZone to put
+            //  cells into
             const wordList& cellZoneNames() const
             {
                 return cellZoneNames_;
@@ -159,7 +160,7 @@ public:
             //- Get indices of named surfaces (surfaces with faceZoneName)
             labelList getNamedSurfaces() const;
 
-            //- Get indices of closed named surfaces
+            //- Get indices of closed surfaces with a cellZone
             labelList getClosedNamedSurfaces() const;
 
             //- From local region number to global region number

From 1b46e0515851e40647d3134779e7d8d4151ad2d3 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 30 Mar 2010 21:52:48 +0100
Subject: [PATCH 05/20] BUG: allowed sliver triangle to be collapsed down to
 illegal 2 vertex face.

Added simple check to disable collapsing faces to <3 vertices.
---
 .../polyTopoChange/removePoints.C             | 35 ++++++++++++++++++-
 1 file changed, 34 insertions(+), 1 deletion(-)

diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.C
index 78108a6523..382391e3a5 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.C
@@ -203,12 +203,45 @@ Foam::label Foam::removePoints::countPointUsage
             pointCanBeDeleted[pointI] = true;
             nDeleted++;
         }
-
     }
     edge0.clear();
     edge1.clear();
 
 
+    // Protect any points on faces that would collapse down to nothing
+    // No particular intelligence so might protect too many points
+    forAll(mesh_.faces(), faceI)
+    {
+        const face& f = mesh_.faces()[faceI];
+
+        label nCollapse = 0;
+        forAll(f, fp)
+        {
+            if (pointCanBeDeleted[f[fp]])
+            {
+                nCollapse++;
+            }
+        }
+
+        if ((f.size() - nCollapse) < 3)
+        {
+            // Just unmark enough points
+            forAll(f, fp)
+            {
+                if (pointCanBeDeleted[f[fp]])
+                {
+                    pointCanBeDeleted[f[fp]] = false;
+                    --nCollapse;
+                    if (nCollapse == 0)
+                    {
+                        break;
+                    }
+                }
+            }
+        }
+    }
+
+
     // Point can be deleted only if all processors want to delete it
     syncTools::syncPointList
     (

From a5d03f9ba39664fe849c6b52496a8083bba4bc03 Mon Sep 17 00:00:00 2001
From: henry <Henry Weller h.weller@opencfd.co.uk>
Date: Wed, 31 Mar 2010 14:35:51 +0100
Subject: [PATCH 06/20] rhoSimpleFoam: Improved stability by limiting the
 changes in rho

and formulating the enthalpy sources consistently.
---
 .../compressible/rhoSimpleFoam/createFields.H     |  9 +++++++--
 .../solvers/compressible/rhoSimpleFoam/hEqn.H     |  4 ++--
 .../solvers/compressible/rhoSimpleFoam/pEqn.H     | 15 +++++++++------
 3 files changed, 18 insertions(+), 10 deletions(-)

diff --git a/applications/solvers/compressible/rhoSimpleFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/createFields.H
index 690b53760d..122280cfac 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/createFields.H
@@ -44,9 +44,14 @@
     scalar pRefValue = 0.0;
     setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue);
 
-    dimensionedScalar pMin
+    dimensionedScalar rhoMax
     (
-        mesh.solutionDict().subDict("SIMPLE").lookup("pMin")
+        mesh.solutionDict().subDict("SIMPLE").lookup("rhoMax")
+    );
+
+    dimensionedScalar rhoMin
+    (
+        mesh.solutionDict().subDict("SIMPLE").lookup("rhoMin")
     );
 
     Info<< "Creating turbulence model\n" << endl;
diff --git a/applications/solvers/compressible/rhoSimpleFoam/hEqn.H b/applications/solvers/compressible/rhoSimpleFoam/hEqn.H
index e299d99f83..57395e977e 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/hEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/hEqn.H
@@ -5,8 +5,8 @@
       - fvm::Sp(fvc::div(phi), h)
       - fvm::laplacian(turbulence->alphaEff(), h)
      ==
-        fvc::div(phi/fvc::interpolate(rho), p, "div(U,p)")
-      - p*fvc::div(phi/fvc::interpolate(rho))
+        fvc::div(phi/fvc::interpolate(rho), rho/psi, "div(U,p)")
+      - (rho/psi)*fvc::div(phi/fvc::interpolate(rho))
     );
 
     hEqn.relax();
diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
index 38922c99c9..329ebbd27f 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
@@ -1,4 +1,7 @@
 rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+rho.relax();
 
 volScalarField rUA = 1.0/UEqn().A();
 U = rUA*UEqn().H();
@@ -82,15 +85,9 @@ else
 // Explicitly relax pressure for momentum corrector
 p.relax();
 
-rho = thermo.rho();
-rho.relax();
-Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
-
 U -= rUA*fvc::grad(p);
 U.correctBoundaryConditions();
 
-bound(p, pMin);
-
 // For closed-volume cases adjust the pressure and density levels
 // to obey overall mass continuity
 if (closedVolume)
@@ -98,3 +95,9 @@ if (closedVolume)
     p += (initialMass - fvc::domainIntegrate(psi*p))
         /fvc::domainIntegrate(psi);
 }
+
+rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+rho.relax();
+Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;

From 0dd8a2b8a4611ba7ef5519326c6d02597b1458d6 Mon Sep 17 00:00:00 2001
From: henry <Henry Weller h.weller@opencfd.co.uk>
Date: Wed, 31 Mar 2010 14:36:09 +0100
Subject: [PATCH 07/20] kOmega model: Corrected date of publication containing
 the description of the model.

---
 src/turbulenceModels/incompressible/RAS/kOmega/kOmega.H | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.H b/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.H
index cf8f8a07ce..d2759d2927 100644
--- a/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.H
+++ b/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.H
@@ -33,7 +33,7 @@ Description
         "Turbulence Modeling for CFD"
         D. C. Wilcox,
         DCW Industries, Inc., La Canada,
-        California, 1998.
+        California, 1988.
 
         See also:
         http://www.cfd-online.com/Wiki/Wilcox's_k-omega_model

From cf5ae6d629843f86b18d900b4a53fc7c62ad5ccd Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Wed, 31 Mar 2010 15:53:07 +0100
Subject: [PATCH 08/20] ENH: Use concavity check based on face planes.

A cell is concave if there is a face whose face-centre is outside the plane of any other face on the cell.
Another option is to check for any vertices being         outside      ,,    ,,
but this is not consistent with current tracking.
---
 .../primitiveMeshCheck/primitiveMeshCheck.C   | 64 ++++++++-----------
 1 file changed, 28 insertions(+), 36 deletions(-)

diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
index 30b051d1e9..f28e9affd0 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
@@ -2111,7 +2111,6 @@ bool Foam::primitiveMesh::checkConcaveCells
     const labelList& fOwner = faceOwner();
     const vectorField& fAreas = faceAreas();
     const pointField& fCentres = faceCentres();
-    const pointField& pts = points();
 
     label nConcaveCells = 0;
 
@@ -2119,8 +2118,6 @@ bool Foam::primitiveMesh::checkConcaveCells
     {
         const cell& cFaces = c[cellI];
 
-        const labelList& cPoints = cellPoints()[cellI];
-
         bool concave = false;
 
         forAll(cFaces, i)
@@ -2134,8 +2131,6 @@ bool Foam::primitiveMesh::checkConcaveCells
 
             const point& fC = fCentres[fI];
 
-            const face& f = faces()[fI];
-
             vector fN = fAreas[fI];
 
             fN /= max(mag(fN), VSMALL);
@@ -2147,45 +2142,42 @@ bool Foam::primitiveMesh::checkConcaveCells
                 fN *= -1;
             }
 
-            // Is any vertex of the cell on the wrong side of the
-            // plane of this face?
+            // Is the centre of any other face of the cell on the
+            // wrong side of the plane of this face?
 
-            forAll(cPoints, cPtI)
+            forAll(cFaces, j)
             {
-                label ptI = cPoints[cPtI];
-
-                // Skip points that are on this face
-                if (findIndex(f, ptI) > -1)
+                if (j != i)
                 {
-                    continue;
-                }
+                    label fJ = cFaces[j];
 
-                const point& pt = pts[ptI];
+                    const point& pt = fCentres[fJ];
 
-                // If the cell is concave, the point will be on the
-                // positive normal side of the plane of f, defined by
-                // its centre and normal, and the angle between (pt -
-                // fC) and fN will be less than 90 degrees, so the dot
-                // product will be positive.
+                    // If the cell is concave, the point will be on the
+                    // positive normal side of the plane of f, defined by
+                    // its centre and normal, and the angle between (pt -
+                    // fC) and fN will be less than 90 degrees, so the dot
+                    // product will be positive.
 
-                vector pC = (pt - fC);
+                    vector pC = (pt - fC);
 
-                pC /= max(mag(pC), VSMALL);
+                    pC /= max(mag(pC), VSMALL);
 
-                if ((pC & fN) > -planarCosAngle_)
-                {
-                    // Concave or planar face
-
-                    concave = true;
-
-                    if (setPtr)
+                    if ((pC & fN) > -planarCosAngle_)
                     {
-                        setPtr->insert(cellI);
+                        // Concave or planar face
+
+                        concave = true;
+
+                        if (setPtr)
+                        {
+                            setPtr->insert(cellI);
+                        }
+
+                        nConcaveCells++;
+
+                        break;
                     }
-
-                    nConcaveCells++;
-
-                    break;
                 }
             }
         }
@@ -2197,8 +2189,8 @@ bool Foam::primitiveMesh::checkConcaveCells
     {
         if (debug || report)
         {
-            Info<< " ***Concave cells found, number of cells: "
-                << nConcaveCells << endl;
+            Info<< " ***Concave cells (using face planes) found,"
+                << " number of cells: " << nConcaveCells << endl;
         }
 
         return true;

From bb1a2476fddd7670d81c7bdb0e48170cd3f23fea Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Wed, 31 Mar 2010 18:10:05 +0100
Subject: [PATCH 09/20] ENH: Change concave cell checking to be based on face
 centres and face planes, rather than a vertex-by-vertex check.  Assesses
 concavity as relevant to the tracking.

---
 .../meshes/primitiveMesh/primitiveMesh.H      |  2 +-
 .../primitiveMeshCheck/primitiveMeshCheck.C   | 60 ++++++++-----------
 2 files changed, 27 insertions(+), 35 deletions(-)

diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
index b1b54c9aa5..718df480e0 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
@@ -649,7 +649,7 @@ public:
                     labelHashSet* setPtr = NULL
                 ) const;
 
-                //- Check for concave cells
+                //- Check for concave cells by the planes of faces
                 bool checkConcaveCells
                 (
                     const bool report = false,
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
index 30b051d1e9..7dcfceaf10 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
@@ -2111,7 +2111,6 @@ bool Foam::primitiveMesh::checkConcaveCells
     const labelList& fOwner = faceOwner();
     const vectorField& fAreas = faceAreas();
     const pointField& fCentres = faceCentres();
-    const pointField& pts = points();
 
     label nConcaveCells = 0;
 
@@ -2119,8 +2118,6 @@ bool Foam::primitiveMesh::checkConcaveCells
     {
         const cell& cFaces = c[cellI];
 
-        const labelList& cPoints = cellPoints()[cellI];
-
         bool concave = false;
 
         forAll(cFaces, i)
@@ -2134,8 +2131,6 @@ bool Foam::primitiveMesh::checkConcaveCells
 
             const point& fC = fCentres[fI];
 
-            const face& f = faces()[fI];
-
             vector fN = fAreas[fI];
 
             fN /= max(mag(fN), VSMALL);
@@ -2147,45 +2142,42 @@ bool Foam::primitiveMesh::checkConcaveCells
                 fN *= -1;
             }
 
-            // Is any vertex of the cell on the wrong side of the
-            // plane of this face?
+            // Is the centre of any other face of the cell on the
+            // wrong side of the plane of this face?
 
-            forAll(cPoints, cPtI)
+            forAll(cFaces, j)
             {
-                label ptI = cPoints[cPtI];
-
-                // Skip points that are on this face
-                if (findIndex(f, ptI) > -1)
+                if (j != i)
                 {
-                    continue;
-                }
+                    label fJ = cFaces[j];
 
-                const point& pt = pts[ptI];
+                    const point& pt = fCentres[fJ];
 
-                // If the cell is concave, the point will be on the
-                // positive normal side of the plane of f, defined by
-                // its centre and normal, and the angle between (pt -
-                // fC) and fN will be less than 90 degrees, so the dot
-                // product will be positive.
+                    // If the cell is concave, the point will be on the
+                    // positive normal side of the plane of f, defined by
+                    // its centre and normal, and the angle between (pt -
+                    // fC) and fN will be less than 90 degrees, so the dot
+                    // product will be positive.
 
-                vector pC = (pt - fC);
+                    vector pC = (pt - fC);
 
-                pC /= max(mag(pC), VSMALL);
+                    pC /= max(mag(pC), VSMALL);
 
-                if ((pC & fN) > -planarCosAngle_)
-                {
-                    // Concave or planar face
-
-                    concave = true;
-
-                    if (setPtr)
+                    if ((pC & fN) > -planarCosAngle_)
                     {
-                        setPtr->insert(cellI);
+                        // Concave or planar face
+
+                        concave = true;
+
+                        if (setPtr)
+                        {
+                            setPtr->insert(cellI);
+                        }
+
+                        nConcaveCells++;
+
+                        break;
                     }
-
-                    nConcaveCells++;
-
-                    break;
                 }
             }
         }

From fd7e9c9ae49e31ff403f5200d852cfb46e0ba965 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 1 Apr 2010 14:20:27 +0100
Subject: [PATCH 10/20] ENH: Added option for explicitly providing regions.
 Restructured code.

---
 .../splitMeshRegions/splitMeshRegions.C       | 668 ++++++++++--------
 1 file changed, 376 insertions(+), 292 deletions(-)

diff --git a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
index 88ca441e98..7da3c22269 100644
--- a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
+++ b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
@@ -39,17 +39,29 @@ Description
     - mesh with cells put into cellZones (-makeCellZones)
 
     Note:
+    - cellZonesOnly does not do a walk and uses the cellZones only. Use
+    this if you don't mind having disconnected domains in a single region.
+    This option requires all cells to be in one (and one only) cellZone.
+
+    - cellZonesFileOnly behaves like -cellZonesOnly but reads the cellZones
+    from the specified file. This allows one to explicitly specify the region
+    distribution and still have multiple cellZones per region.
+
+    - useCellZonesOnly does not do a walk and uses the cellZones only. Use
+    this if you don't mind having disconnected domains in a single region.
+    This option requires all cells to be in one (and one only) cellZone.
+
+
     - Should work in parallel.
     cellZones can differ on either side of processor boundaries in which case
     the faces get moved from processor patch to directMapped patch. Not
-    ery well tested.
+    very well tested.
+
     - If a cell zone gets split into more than one region it can detect
     the largest matching region (-sloppyCellZones). This will accept any
     region that covers more than 50% of the zone. It has to be a subset
     so cannot have any cells in any other zone.
-    - useCellZonesOnly does not do a walk and uses the cellZones only. Use
-    this if you don't mind having disconnected domains in a single region.
-    This option requires all cells to be in one (and one only) cellZone.
+
     - writes maps like decomposePar back to original mesh:
         - pointRegionAddressing : for every point in this region the point in
         the original mesh
@@ -1137,63 +1149,6 @@ EdgeMap<label> addRegionPatches
 }
 
 
-//// Checks if regionI in cellRegion is subset of existing cellZone. Returns -1
-//// if no zone found, zone otherwise
-//label findCorrespondingSubZone
-//(
-//    const cellZoneMesh& cellZones,
-//    const labelList& existingZoneID,
-//    const labelList& cellRegion,
-//    const label regionI
-//)
-//{
-//    // Zone corresponding to region. No corresponding zone.
-//    label zoneI = labelMax;
-//
-//    labelList regionCells = findIndices(cellRegion, regionI);
-//
-//    if (regionCells.empty())
-//    {
-//        // My local portion is empty. Maps to any empty cellZone. Mark with
-//        // special value which can get overwritten by other processors.
-//        zoneI = -1;
-//    }
-//    else
-//    {
-//        // Get zone for first element.
-//        zoneI = existingZoneID[regionCells[0]];
-//
-//        if (zoneI == -1)
-//        {
-//            zoneI = labelMax;
-//        }
-//        else
-//        {
-//            // 1. All regionCells in zoneI?
-//            forAll(regionCells, i)
-//            {
-//                if (existingZoneID[regionCells[i]] != zoneI)
-//                {
-//                    zoneI = labelMax;
-//                    break;
-//                }
-//            }
-//        }
-//    }
-//
-//    // Determine same zone over all processors.
-//    reduce(zoneI, maxOp<label>());
-//
-//    if (zoneI == labelMax)
-//    {
-//        // Cells in region that are not in zoneI
-//        zoneI = -1;
-//    }
-//
-//    return zoneI;
-//}
-
-
 // Find region that covers most of cell zone
 label findCorrespondingRegion
 (
@@ -1314,62 +1269,20 @@ label findCorrespondingRegion
 //}
 
 
-// Main program:
-
-int main(int argc, char *argv[])
+// Get zone per cell
+// - non-unique zoning
+// - coupled zones
+void getZoneID
+(
+    const polyMesh& mesh,
+    const cellZoneMesh& cellZones,
+    labelList& zoneID,
+    labelList& neiZoneID
+)
 {
-#   include "addOverwriteOption.H"
-    argList::addBoolOption("cellZones");
-    argList::addBoolOption("cellZonesOnly");
-    argList::addOption("blockedFaces", "faceSet");
-    argList::addBoolOption("makeCellZones");
-    argList::addBoolOption("largestOnly");
-    argList::addOption("insidePoint", "point");
-    argList::addBoolOption("detectOnly");
-    argList::addBoolOption("sloppyCellZones");
-
-#   include "setRootCase.H"
-#   include "createTime.H"
-    runTime.functionObjects().off();
-#   include "createMesh.H"
-    const word oldInstance = mesh.pointsInstance();
-
-    word blockedFacesName;
-    if (args.optionReadIfPresent("blockedFaces", blockedFacesName))
-    {
-        Info<< "Reading blocked internal faces from faceSet "
-            << blockedFacesName << nl << endl;
-    }
-
-    const bool makeCellZones    = args.optionFound("makeCellZones");
-    const bool largestOnly      = args.optionFound("largestOnly");
-    const bool insidePoint      = args.optionFound("insidePoint");
-    const bool useCellZones     = args.optionFound("cellZones");
-    const bool useCellZonesOnly = args.optionFound("cellZonesOnly");
-    const bool overwrite        = args.optionFound("overwrite");
-    const bool detectOnly       = args.optionFound("detectOnly");
-    const bool sloppyCellZones  = args.optionFound("sloppyCellZones");
-
-    if (insidePoint && largestOnly)
-    {
-        FatalErrorIn(args.executable())
-            << "You cannot specify both -largestOnly"
-            << " (keep region with most cells)"
-            << " and -insidePoint (keep region containing point)"
-            << exit(FatalError);
-    }
-
-
-    const cellZoneMesh& cellZones = mesh.cellZones();
-
-
-    // Collect zone per cell
-    // ~~~~~~~~~~~~~~~~~~~~~
-    // - non-unique zoning
-    // - coupled zones
-
     // Existing zoneID
-    labelList zoneID(mesh.nCells(), -1);
+    zoneID.setSize(mesh.nCells());
+    zoneID = -1;
 
     forAll(cellZones, zoneI)
     {
@@ -1384,7 +1297,7 @@ int main(int argc, char *argv[])
             }
             else
             {
-                FatalErrorIn(args.executable())
+                FatalErrorIn("getZoneID(..)")
                     << "Cell " << cellI << " with cell centre "
                     << mesh.cellCentres()[cellI]
                     << " is multiple zones. This is not allowed." << endl
@@ -1396,175 +1309,42 @@ int main(int argc, char *argv[])
     }
 
     // Neighbour zoneID.
-    labelList neiZoneID(mesh.nFaces()-mesh.nInternalFaces());
+    neiZoneID.setSize(mesh.nFaces()-mesh.nInternalFaces());
 
     forAll(neiZoneID, i)
     {
         neiZoneID[i] = zoneID[mesh.faceOwner()[i+mesh.nInternalFaces()]];
     }
     syncTools::swapBoundaryFaceList(mesh, neiZoneID, false);
+}
 
 
-    // Determine connected regions
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+void matchRegions
+(
+    const bool sloppyCellZones,
+    const polyMesh& mesh,
 
-    // Mark additional faces that are blocked
-    boolList blockedFace;
+    const label nCellRegions,
+    const labelList& cellRegion,
 
-    // Read from faceSet
-    if (blockedFacesName.size())
-    {
-        faceSet blockedFaceSet(mesh, blockedFacesName);
-        Info<< "Read " << returnReduce(blockedFaceSet.size(), sumOp<label>())
-            << " blocked faces from set " << blockedFacesName << nl << endl;
+    labelList& regionToZone,
+    wordList& regionNames,
+    labelList& zoneToRegion
+)
+{
+    const cellZoneMesh& cellZones = mesh.cellZones();
 
-        blockedFace.setSize(mesh.nFaces(), false);
-
-        forAllConstIter(faceSet, blockedFaceSet, iter)
-        {
-            blockedFace[iter.key()] = true;
-        }
-    }
-
-    // Imply from differing cellZones
-    if (useCellZones)
-    {
-        blockedFace.setSize(mesh.nFaces(), false);
-
-        for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
-        {
-            label own = mesh.faceOwner()[faceI];
-            label nei = mesh.faceNeighbour()[faceI];
-
-            if (zoneID[own] != zoneID[nei])
-            {
-                blockedFace[faceI] = true;
-            }
-        }
-
-        // Different cellZones on either side of processor patch.
-        forAll(neiZoneID, i)
-        {
-            label faceI = i+mesh.nInternalFaces();
-
-            if (zoneID[mesh.faceOwner()[faceI]] != neiZoneID[i])
-            {
-                blockedFace[faceI] = true;
-            }
-        }
-    }
-
-
-    // Determine per cell the region it belongs to
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    // cellRegion is the labelList with the region per cell.
-    labelList cellRegion;
-    label nCellRegions = 0;
-    if (useCellZonesOnly)
-    {
-        label unzonedCellI = findIndex(zoneID, -1);
-        if (unzonedCellI != -1)
-        {
-            FatalErrorIn(args.executable())
-                << "For the cellZonesOnly option all cells "
-                << "have to be in a cellZone." << endl
-                << "Cell " << unzonedCellI
-                << " at" << mesh.cellCentres()[unzonedCellI]
-                << " is not in a cellZone. There might be more unzoned cells."
-                << exit(FatalError);
-        }
-        cellRegion = zoneID;
-        nCellRegions = gMax(cellRegion)+1;
-    }
-    else
-    {
-        // Do a topological walk to determine regions
-        regionSplit regions(mesh, blockedFace);
-        nCellRegions = regions.nRegions();
-        cellRegion.transfer(regions);
-    }
-
-    Info<< endl << "Number of regions:" << nCellRegions << nl << endl;
-
-
-    // Write to manual decomposition option
-    {
-        labelIOList cellToRegion
-        (
-            IOobject
-            (
-                "cellToRegion",
-                mesh.facesInstance(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            cellRegion
-        );
-        cellToRegion.write();
-
-        Info<< "Writing region per cell file (for manual decomposition) to "
-            << cellToRegion.objectPath() << nl << endl;
-    }
-    // Write for postprocessing
-    {
-        volScalarField cellToRegion
-        (
-            IOobject
-            (
-                "cellToRegion",
-                runTime.timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            mesh,
-            dimensionedScalar("zero", dimless, 0),
-            zeroGradientFvPatchScalarField::typeName
-        );
-        forAll(cellRegion, cellI)
-        {
-            cellToRegion[cellI] = cellRegion[cellI];
-        }
-        cellToRegion.write();
-
-        Info<< "Writing region per cell as volScalarField to "
-            << cellToRegion.objectPath() << nl << endl;
-    }
-
-
-    // Sizes per region
-    // ~~~~~~~~~~~~~~~~
-
-    labelList regionSizes(nCellRegions, 0);
-
-    forAll(cellRegion, cellI)
-    {
-        regionSizes[cellRegion[cellI]]++;
-    }
-    forAll(regionSizes, regionI)
-    {
-        reduce(regionSizes[regionI], sumOp<label>());
-    }
-
-    Info<< "Region\tCells" << nl
-        << "------\t-----" << endl;
-
-    forAll(regionSizes, regionI)
-    {
-        Info<< regionI << '\t' << regionSizes[regionI] << nl;
-    }
-    Info<< endl;
+    regionToZone.setSize(nCellRegions, -1);
+    regionNames.setSize(nCellRegions);
+    zoneToRegion.setSize(cellZones.size(), -1);
 
+    // Get current per cell zoneID
+    labelList zoneID(mesh.nCells(), -1);
+    labelList neiZoneID(mesh.nFaces()-mesh.nInternalFaces());
+    getZoneID(mesh, cellZones, zoneID, neiZoneID);
 
     // Sizes per cellzone
-    // ~~~~~~~~~~~~~~~~~~
-
     labelList zoneSizes(cellZones.size(), 0);
-    if (useCellZones || makeCellZones || sloppyCellZones)
     {
         List<wordList> zoneNames(Pstream::nProcs());
         zoneNames[Pstream::myProcNo()] = cellZones.names();
@@ -1575,7 +1355,7 @@ int main(int argc, char *argv[])
         {
             if (zoneNames[procI] != zoneNames[0])
             {
-                FatalErrorIn(args.executable())
+                FatalErrorIn("matchRegions(..)")
                     << "cellZones not synchronised across processors." << endl
                     << "Master has cellZones " << zoneNames[0] << endl
                     << "Processor " << procI
@@ -1595,16 +1375,6 @@ int main(int argc, char *argv[])
     }
 
 
-    // Whether region corresponds to a cellzone
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Region per zone
-    labelList regionToZone(nCellRegions, -1);
-    // Name of region
-    wordList regionNames(nCellRegions);
-    // Zone to region
-    labelList zoneToRegion(cellZones.size(), -1);
-
     if (sloppyCellZones)
     {
         Info<< "Trying to match regions to existing cell zones;"
@@ -1624,7 +1394,7 @@ int main(int argc, char *argv[])
             if (regionI != -1)
             {
                 Info<< "Sloppily matched region " << regionI
-                    << " size " << regionSizes[regionI]
+                    //<< " size " << regionSizes[regionI]
                     << " to zone " << zoneI << " size " << zoneSizes[zoneI]
                     << endl;
                 zoneToRegion[zoneI] = regionI;
@@ -1664,6 +1434,330 @@ int main(int argc, char *argv[])
             regionNames[regionI] = "domain" + Foam::name(regionI);
         }
     }
+}
+
+
+void writeCellToRegion(const fvMesh& mesh, const labelList& cellRegion)
+{
+    // Write to manual decomposition option
+    {
+        labelIOList cellToRegion
+        (
+            IOobject
+            (
+                "cellToRegion",
+                mesh.facesInstance(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            cellRegion
+        );
+        cellToRegion.write();
+
+        Info<< "Writing region per cell file (for manual decomposition) to "
+            << cellToRegion.objectPath() << nl << endl;
+    }
+    // Write for postprocessing
+    {
+        volScalarField cellToRegion
+        (
+            IOobject
+            (
+                "cellToRegion",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh,
+            dimensionedScalar("zero", dimless, 0),
+            zeroGradientFvPatchScalarField::typeName
+        );
+        forAll(cellRegion, cellI)
+        {
+            cellToRegion[cellI] = cellRegion[cellI];
+        }
+        cellToRegion.write();
+
+        Info<< "Writing region per cell as volScalarField to "
+            << cellToRegion.objectPath() << nl << endl;
+    }
+}
+
+
+
+// Main program:
+
+int main(int argc, char *argv[])
+{
+#   include "addOverwriteOption.H"
+    argList::addBoolOption("cellZones");
+    argList::addBoolOption("cellZonesOnly");
+    argList::addOption("cellZonesFileOnly", "cellZonesName");
+    argList::addOption("blockedFaces", "faceSet");
+    argList::addBoolOption("makeCellZones");
+    argList::addBoolOption("largestOnly");
+    argList::addOption("insidePoint", "point");
+    argList::addBoolOption("detectOnly");
+    argList::addBoolOption("sloppyCellZones");
+
+#   include "setRootCase.H"
+#   include "createTime.H"
+    runTime.functionObjects().off();
+#   include "createMesh.H"
+    const word oldInstance = mesh.pointsInstance();
+
+    word blockedFacesName;
+    if (args.optionReadIfPresent("blockedFaces", blockedFacesName))
+    {
+        Info<< "Reading blocked internal faces from faceSet "
+            << blockedFacesName << nl << endl;
+    }
+
+    const bool makeCellZones    = args.optionFound("makeCellZones");
+    const bool largestOnly      = args.optionFound("largestOnly");
+    const bool insidePoint      = args.optionFound("insidePoint");
+    const bool useCellZones     = args.optionFound("cellZones");
+    const bool useCellZonesOnly = args.optionFound("cellZonesOnly");
+    const bool useCellZonesFile = args.optionFound("cellZonesFileOnly");
+    const bool overwrite        = args.optionFound("overwrite");
+    const bool detectOnly       = args.optionFound("detectOnly");
+    const bool sloppyCellZones  = args.optionFound("sloppyCellZones");
+    if
+    (
+        (useCellZonesOnly || useCellZonesFile)
+     && (
+            blockedFacesName != word::null
+         || useCellZones
+        )
+    )
+    {
+        FatalErrorIn(args.executable())
+            << "You cannot specify both -cellZonesOnly or -cellZonesFileOnly"
+            << " (which specify complete split)"
+            << " in combination with -blockedFaces or -cellZones"
+            << " (which imply a split based on topology)"
+            << exit(FatalError);
+    }
+
+
+
+    if (insidePoint && largestOnly)
+    {
+        FatalErrorIn(args.executable())
+            << "You cannot specify both -largestOnly"
+            << " (keep region with most cells)"
+            << " and -insidePoint (keep region containing point)"
+            << exit(FatalError);
+    }
+
+
+    const cellZoneMesh& cellZones = mesh.cellZones();
+
+    // Existing zoneID
+    labelList zoneID(mesh.nCells(), -1);
+    // Neighbour zoneID.
+    labelList neiZoneID(mesh.nFaces()-mesh.nInternalFaces());
+    getZoneID(mesh, cellZones, zoneID, neiZoneID);
+
+
+
+    // Determine per cell the region it belongs to
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // cellRegion is the labelList with the region per cell.
+    labelList cellRegion;
+    // Region per zone
+    labelList regionToZone;
+    // Name of region
+    wordList regionNames;
+    // Zone to region
+    labelList zoneToRegion;
+
+    label nCellRegions = 0;
+    if (useCellZonesOnly)
+    {
+        Info<< "Using current cellZones to split mesh into regions."
+            << " This requires all"
+            << " cells to be in one and only one cellZone." << nl << endl;
+
+        label unzonedCellI = findIndex(zoneID, -1);
+        if (unzonedCellI != -1)
+        {
+            FatalErrorIn(args.executable())
+                << "For the cellZonesOnly option all cells "
+                << "have to be in a cellZone." << endl
+                << "Cell " << unzonedCellI
+                << " at" << mesh.cellCentres()[unzonedCellI]
+                << " is not in a cellZone. There might be more unzoned cells."
+                << exit(FatalError);
+        }
+        cellRegion = zoneID;
+        nCellRegions = gMax(cellRegion)+1;
+        regionToZone.setSize(nCellRegions);
+        regionNames.setSize(nCellRegions);
+        zoneToRegion.setSize(cellZones.size(), -1);
+        for (label regionI = 0; regionI < nCellRegions; regionI++)
+        {
+            regionToZone[regionI] = regionI;
+            zoneToRegion[regionI] = regionI;
+            regionNames[regionI] = cellZones[regionI].name();
+        }
+    }
+    else if (useCellZonesFile)
+    {
+        const word zoneFile = args.option("cellZonesFileOnly");
+        Info<< "Reading split from cellZones file " << zoneFile << endl
+            << "This requires all"
+            << " cells to be in one and only one cellZone." << nl << endl;
+
+        cellZoneMesh newCellZones
+        (
+            IOobject
+            (
+                zoneFile,
+                mesh.facesInstance(),
+                polyMesh::meshSubDir,
+                mesh,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh
+        );
+
+        labelList newZoneID(mesh.nCells(), -1);
+        labelList newNeiZoneID(mesh.nFaces()-mesh.nInternalFaces());
+        getZoneID(mesh, newCellZones, newZoneID, newNeiZoneID);
+
+        label unzonedCellI = findIndex(newZoneID, -1);
+        if (unzonedCellI != -1)
+        {
+            FatalErrorIn(args.executable())
+                << "For the cellZonesFileOnly option all cells "
+                << "have to be in a cellZone." << endl
+                << "Cell " << unzonedCellI
+                << " at" << mesh.cellCentres()[unzonedCellI]
+                << " is not in a cellZone. There might be more unzoned cells."
+                << exit(FatalError);
+        }
+        cellRegion = newZoneID;
+        nCellRegions = gMax(cellRegion)+1;
+        zoneToRegion.setSize(newCellZones.size(), -1);
+        regionToZone.setSize(nCellRegions);
+        regionNames.setSize(nCellRegions);
+        for (label regionI = 0; regionI < nCellRegions; regionI++)
+        {
+            regionToZone[regionI] = regionI;
+            zoneToRegion[regionI] = regionI;
+            regionNames[regionI] = newCellZones[regionI].name();
+        }
+    }
+    else
+    {
+        // Determine connected regions
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        // Mark additional faces that are blocked
+        boolList blockedFace;
+
+        // Read from faceSet
+        if (blockedFacesName.size())
+        {
+            faceSet blockedFaceSet(mesh, blockedFacesName);
+            Info<< "Read "
+                << returnReduce(blockedFaceSet.size(), sumOp<label>())
+                << " blocked faces from set " << blockedFacesName << nl << endl;
+
+            blockedFace.setSize(mesh.nFaces(), false);
+
+            forAllConstIter(faceSet, blockedFaceSet, iter)
+            {
+                blockedFace[iter.key()] = true;
+            }
+        }
+
+        // Imply from differing cellZones
+        if (useCellZones)
+        {
+            blockedFace.setSize(mesh.nFaces(), false);
+
+            for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
+            {
+                label own = mesh.faceOwner()[faceI];
+                label nei = mesh.faceNeighbour()[faceI];
+
+                if (zoneID[own] != zoneID[nei])
+                {
+                    blockedFace[faceI] = true;
+                }
+            }
+
+            // Different cellZones on either side of processor patch.
+            forAll(neiZoneID, i)
+            {
+                label faceI = i+mesh.nInternalFaces();
+
+                if (zoneID[mesh.faceOwner()[faceI]] != neiZoneID[i])
+                {
+                    blockedFace[faceI] = true;
+                }
+            }
+        }
+
+        // Do a topological walk to determine regions
+        regionSplit regions(mesh, blockedFace);
+        nCellRegions = regions.nRegions();
+        cellRegion.transfer(regions);
+
+        // Make up region names. If possible match them to existing zones.
+        matchRegions
+        (
+            sloppyCellZones,
+            mesh,
+            nCellRegions,
+            cellRegion,
+
+            regionToZone,
+            regionNames,
+            zoneToRegion
+        );
+    }
+
+    Info<< endl << "Number of regions:" << nCellRegions << nl << endl;
+
+
+    // Write decomposition to file
+    writeCellToRegion(mesh, cellRegion);
+
+
+
+    // Sizes per region
+    // ~~~~~~~~~~~~~~~~
+
+    labelList regionSizes(nCellRegions, 0);
+
+    forAll(cellRegion, cellI)
+    {
+        regionSizes[cellRegion[cellI]]++;
+    }
+    forAll(regionSizes, regionI)
+    {
+        reduce(regionSizes[regionI], sumOp<label>());
+    }
+
+    Info<< "Region\tCells" << nl
+        << "------\t-----" << endl;
+
+    forAll(regionSizes, regionI)
+    {
+        Info<< regionI << '\t' << regionSizes[regionI] << nl;
+    }
+    Info<< endl;
+
 
 
     // Print region to zone
@@ -1676,16 +1770,6 @@ int main(int argc, char *argv[])
     }
     Info<< endl;
 
-    //// Print zone to region
-    //Info<< "Zone\tName\tRegion" << nl
-    //    << "----\t----\t------" << endl;
-    //forAll(zoneToRegion, zoneI)
-    //{
-    //    Info<< zoneI << '\t' << cellZones[zoneI].name() << '\t'
-    //        << zoneToRegion[zoneI] << nl;
-    //}
-    //Info<< endl;
-
 
 
     // Since we're going to mess with patches make sure all non-processor ones

From fd952687fd2490214522ef5ce46cc78ad3414ff2 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 1 Apr 2010 14:21:25 +0100
Subject: [PATCH 11/20] ENH: Removed extraneous if

---
 .../fluid/compressibleMultiRegionCourantNo.H  | 26 +++++++++----------
 1 file changed, 12 insertions(+), 14 deletions(-)

diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H
index 8f1f654d2b..b5db5078f1 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H
@@ -1,18 +1,16 @@
     scalar CoNum = -GREAT;
-    if (fluidRegions.size())
+
+    forAll(fluidRegions, regionI)
     {
-        forAll(fluidRegions, regionI)
-        {
-            CoNum = max
+        CoNum = max
+        (
+            compressibleCourantNo
             (
-                compressibleCourantNo
-                (
-                    fluidRegions[regionI],
-                    runTime,
-                    rhoFluid[regionI],
-                    phiFluid[regionI]
-                ),
-                CoNum
-            );
-        }
+                fluidRegions[regionI],
+                runTime,
+                rhoFluid[regionI],
+                phiFluid[regionI]
+            ),
+            CoNum
+        );
     }

From 16cf8eaf9e05b7334ea40b6cc8027ba2540ab9ce Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 1 Apr 2010 14:22:06 +0100
Subject: [PATCH 12/20] ENH: Changed flux condition into gradient bc

---
 ...allHeatFluxTemperatureFvPatchScalarField.C | 29 ++++++++++---------
 ...allHeatFluxTemperatureFvPatchScalarField.H |  4 +--
 2 files changed, 17 insertions(+), 16 deletions(-)

diff --git a/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C b/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C
index bf48d47292..e2b1d536d2 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C
@@ -37,7 +37,7 @@ solidWallHeatFluxTemperatureFvPatchScalarField
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    fixedValueFvPatchScalarField(p, iF),
+    fixedGradientFvPatchScalarField(p, iF),
     q_(p.size(), 0.0),
     KName_("undefined-K")
 {}
@@ -52,7 +52,7 @@ solidWallHeatFluxTemperatureFvPatchScalarField
     const fvPatchFieldMapper& mapper
 )
 :
-    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
+    fixedGradientFvPatchScalarField(ptf, p, iF, mapper),
     q_(ptf.q_, mapper),
     KName_(ptf.KName_)
 {}
@@ -66,7 +66,7 @@ solidWallHeatFluxTemperatureFvPatchScalarField
     const dictionary& dict
 )
 :
-    fixedValueFvPatchScalarField(p, iF, dict),
+    fixedGradientFvPatchScalarField(p, iF, dict),
     q_("q", dict, p.size()),
     KName_(dict.lookup("K"))
 {}
@@ -78,7 +78,7 @@ solidWallHeatFluxTemperatureFvPatchScalarField
     const solidWallHeatFluxTemperatureFvPatchScalarField& tppsf
 )
 :
-    fixedValueFvPatchScalarField(tppsf),
+    fixedGradientFvPatchScalarField(tppsf),
     q_(tppsf.q_),
     KName_(tppsf.KName_)
 {}
@@ -91,7 +91,7 @@ solidWallHeatFluxTemperatureFvPatchScalarField
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    fixedValueFvPatchScalarField(tppsf, iF),
+    fixedGradientFvPatchScalarField(tppsf, iF),
     q_(tppsf.q_),
     KName_(tppsf.KName_)
 {}
@@ -104,7 +104,7 @@ void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::autoMap
     const fvPatchFieldMapper& m
 )
 {
-    fixedValueFvPatchScalarField::autoMap(m);
+    fixedGradientFvPatchScalarField::autoMap(m);
     q_.autoMap(m);
 }
 
@@ -115,7 +115,7 @@ void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::rmap
     const labelList& addr
 )
 {
-    fixedValueFvPatchScalarField::rmap(ptf, addr);
+    fixedGradientFvPatchScalarField::rmap(ptf, addr);
 
     const solidWallHeatFluxTemperatureFvPatchScalarField& hfptf =
         refCast<const solidWallHeatFluxTemperatureFvPatchScalarField>(ptf);
@@ -131,14 +131,14 @@ void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::updateCoeffs()
         return;
     }
 
-    const scalarField& Kw =
-        patch().lookupPatchField<volScalarField, scalar>(KName_);
+    const scalarField& Kw = patch().lookupPatchField<volScalarField, scalar>
+    (
+        KName_
+    );
 
-    const fvPatchScalarField& Tw = *this;
+    gradient() = q_/Kw;
 
-    operator==(q_/(patch().deltaCoeffs()*Kw) + Tw.patchInternalField());
-
-    fixedValueFvPatchScalarField::updateCoeffs();
+    fixedGradientFvPatchScalarField::updateCoeffs();
 }
 
 
@@ -147,9 +147,10 @@ void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::write
     Ostream& os
 ) const
 {
-    fixedValueFvPatchScalarField::write(os);
+    fixedGradientFvPatchScalarField::write(os);
     q_.writeEntry("q", os);
     os.writeKeyword("K") << KName_ << token::END_STATEMENT << nl;
+    this->writeEntry("value", os);
 }
 
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.H b/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.H
index 7558c4df7e..ffbcf52ed5 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.H
@@ -45,7 +45,7 @@ SourceFiles
 #ifndef solidWallHeatFluxTemperatureFvPatchScalarField_H
 #define solidWallHeatFluxTemperatureFvPatchScalarField_H
 
-#include "fixedValueFvPatchFields.H"
+#include "fixedGradientFvPatchFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -58,7 +58,7 @@ namespace Foam
 
 class solidWallHeatFluxTemperatureFvPatchScalarField
 :
-    public fixedValueFvPatchScalarField
+    public fixedGradientFvPatchScalarField
 {
     // Private data
 

From 0b54b4b12d7f8d6debc075f173b990e001f51915 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 1 Apr 2010 14:22:49 +0100
Subject: [PATCH 13/20] STYLE: mixed bc already present in
 compressible::turbulentTemperatureBaffleMixed

---
 .../chtMultiRegionFoam/Make/files             |   1 -
 ...ixedTemperatureCoupledFvPatchScalarField.C | 381 ------------------
 ...ixedTemperatureCoupledFvPatchScalarField.H | 192 ---------
 3 files changed, 574 deletions(-)
 delete mode 100644 applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C
 delete mode 100644 applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.H

diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files
index dcec73b268..68d70ae718 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files
@@ -1,7 +1,6 @@
 regionProperties/regionProperties.C
 
 derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C
-derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C
 
 fluid/compressibleCourantNo.C
 solid/solidRegionDiffNo.C
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C
deleted file mode 100644
index 42a73343cc..0000000000
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C
+++ /dev/null
@@ -1,381 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "solidWallMixedTemperatureCoupledFvPatchScalarField.H"
-#include "addToRunTimeSelectionTable.H"
-#include "fvPatchFieldMapper.H"
-#include "volFields.H"
-#include "directMappedPatchBase.H"
-#include "regionProperties.H"
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-bool Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::interfaceOwner
-(
-    const polyMesh& nbrRegion
-) const
-{
-    const fvMesh& myRegion = patch().boundaryMesh().mesh();
-
-    const regionProperties& props =
-        myRegion.objectRegistry::parent().lookupObject<regionProperties>
-        (
-            "regionProperties"
-        );
-
-    label myIndex = findIndex(props.fluidRegionNames(), myRegion.name());
-    if (myIndex == -1)
-    {
-        label i = findIndex(props.solidRegionNames(), myRegion.name());
-
-        if (i == -1)
-        {
-            FatalErrorIn
-            (
-                "solidWallMixedTemperatureCoupledFvPatchScalarField"
-                "::interfaceOwner(const polyMesh&) const"
-            )   << "Cannot find region " << myRegion.name()
-                << " neither in fluids " << props.fluidRegionNames()
-                << " nor in solids " << props.solidRegionNames()
-                << exit(FatalError);
-        }
-        myIndex = props.fluidRegionNames().size() + i;
-    }
-    label nbrIndex = findIndex(props.fluidRegionNames(), nbrRegion.name());
-    if (nbrIndex == -1)
-    {
-        label i = findIndex(props.solidRegionNames(), nbrRegion.name());
-
-        if (i == -1)
-        {
-            FatalErrorIn("coupleManager::interfaceOwner(const polyMesh&) const")
-                << "Cannot find region " << nbrRegion.name()
-                << " neither in fluids " << props.fluidRegionNames()
-                << " nor in solids " << props.solidRegionNames()
-                << exit(FatalError);
-        }
-        nbrIndex = props.fluidRegionNames().size() + i;
-    }
-
-    return myIndex < nbrIndex;
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::
-solidWallMixedTemperatureCoupledFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    mixedFvPatchScalarField(p, iF),
-    neighbourFieldName_("undefined-neighbourFieldName"),
-    KName_("undefined-K")
-{
-    this->refValue() = 0.0;
-    this->refGrad() = 0.0;
-    this->valueFraction() = 1.0;
-    this->fixesValue_ = true;
-}
-
-
-Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::
-solidWallMixedTemperatureCoupledFvPatchScalarField
-(
-    const solidWallMixedTemperatureCoupledFvPatchScalarField& ptf,
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const fvPatchFieldMapper& mapper
-)
-:
-    mixedFvPatchScalarField(ptf, p, iF, mapper),
-    neighbourFieldName_(ptf.neighbourFieldName_),
-    KName_(ptf.KName_),
-    fixesValue_(ptf.fixesValue_)
-{}
-
-
-Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::
-solidWallMixedTemperatureCoupledFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const dictionary& dict
-)
-:
-    mixedFvPatchScalarField(p, iF),
-    neighbourFieldName_(dict.lookup("neighbourFieldName")),
-    KName_(dict.lookup("K"))
-{
-    if (!isA<directMappedPatchBase>(this->patch().patch()))
-    {
-        FatalErrorIn
-        (
-            "solidWallMixedTemperatureCoupledFvPatchScalarField::"
-            "solidWallMixedTemperatureCoupledFvPatchScalarField\n"
-            "(\n"
-            "    const fvPatch& p,\n"
-            "    const DimensionedField<scalar, volMesh>& iF,\n"
-            "    const dictionary& dict\n"
-            ")\n"
-        )   << "\n    patch type '" << p.type()
-            << "' not type '" << directMappedPatchBase::typeName << "'"
-            << "\n    for patch " << p.name()
-            << " of field " << dimensionedInternalField().name()
-            << " in file " << dimensionedInternalField().objectPath()
-            << exit(FatalError);
-    }
-
-    fvPatchScalarField::operator=(scalarField("value", dict, p.size()));
-
-    if (dict.found("refValue"))
-    {
-        // Full restart
-        refValue() = scalarField("refValue", dict, p.size());
-        refGrad() = scalarField("refGradient", dict, p.size());
-        valueFraction() = scalarField("valueFraction", dict, p.size());
-        fixesValue_ = readBool(dict.lookup("fixesValue"));
-    }
-    else
-    {
-        // Start from user entered data. Assume fixedValue.
-        refValue() = *this;
-        refGrad() = 0.0;
-        valueFraction() = 1.0;
-        fixesValue_ = true;
-    }
-}
-
-
-Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::
-solidWallMixedTemperatureCoupledFvPatchScalarField
-(
-    const solidWallMixedTemperatureCoupledFvPatchScalarField& wtcsf,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    mixedFvPatchScalarField(wtcsf, iF),
-    neighbourFieldName_(wtcsf.neighbourFieldName_),
-    KName_(wtcsf.KName_),
-    fixesValue_(wtcsf.fixesValue_)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-const Foam::fvPatchScalarField&
-Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::K() const
-{
-    return this->patch().lookupPatchField<volScalarField, scalar>(KName_);
-}
-
-
-void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::updateCoeffs()
-{
-    if (updated())
-    {
-        return;
-    }
-
-    // Get the coupling information from the directMappedPatchBase
-    const directMappedPatchBase& mpp = refCast<const directMappedPatchBase>
-    (
-        patch().patch()
-    );
-    const polyMesh& nbrMesh = mpp.sampleMesh();
-    // Force recalculation of mapping and schedule
-    const mapDistribute& distMap = mpp.map();
-    (void)distMap.schedule();
-
-    tmp<scalarField> intFld = patchInternalField();
-
-    if (interfaceOwner(nbrMesh))
-    {
-        // Note: other side information could be cached - it only needs
-        // to be updated the first time round the iteration (i.e. when
-        // switching regions) but unfortunately we don't have this information.
-
-        const fvPatch& nbrPatch = refCast<const fvMesh>
-        (
-            nbrMesh
-        ).boundary()[mpp.samplePolyPatch().index()];
-
-
-        // Calculate the temperature by harmonic averaging
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-        const solidWallMixedTemperatureCoupledFvPatchScalarField& nbrField =
-        refCast<const solidWallMixedTemperatureCoupledFvPatchScalarField>
-        (
-            nbrPatch.lookupPatchField<volScalarField, scalar>
-            (
-                neighbourFieldName_
-            )
-        );
-
-        // Swap to obtain full local values of neighbour internal field
-        scalarField nbrIntFld = nbrField.patchInternalField();
-        mapDistribute::distribute
-        (
-            Pstream::defaultCommsType,
-            distMap.schedule(),
-            distMap.constructSize(),
-            distMap.subMap(),           // what to send
-            distMap.constructMap(),     // what to receive
-            nbrIntFld
-        );
-
-        // Swap to obtain full local values of neighbour K*delta
-        scalarField nbrKDelta = nbrField.K()*nbrPatch.deltaCoeffs();
-        mapDistribute::distribute
-        (
-            Pstream::defaultCommsType,
-            distMap.schedule(),
-            distMap.constructSize(),
-            distMap.subMap(),           // what to send
-            distMap.constructMap(),     // what to receive
-            nbrKDelta
-        );
-
-
-        tmp<scalarField> myKDelta = K()*patch().deltaCoeffs();
-
-        // Calculate common wall temperature. Reuse *this to store common value.
-        scalarField Twall
-        (
-            (myKDelta()*intFld() + nbrKDelta*nbrIntFld)
-          / (myKDelta() + nbrKDelta)
-        );
-        // Assign to me
-        fvPatchScalarField::operator=(Twall);
-        // Distribute back and assign to neighbour
-        mapDistribute::distribute
-        (
-            Pstream::defaultCommsType,
-            distMap.schedule(),
-            nbrField.size(),
-            distMap.constructMap(),     // reverse : what to send
-            distMap.subMap(),
-            Twall
-        );
-        const_cast<solidWallMixedTemperatureCoupledFvPatchScalarField&>
-        (
-            nbrField
-        ).fvPatchScalarField::operator=(Twall);
-    }
-
-
-    // Switch between fixed value (of harmonic avg) or gradient
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    label nFixed = 0;
-
-    // Like snGrad but bypass switching on refValue/refGrad.
-    tmp<scalarField> normalGradient = (*this-intFld())*patch().deltaCoeffs();
-
-    if (debug)
-    {
-        scalar Q = gSum(K()*patch().magSf()*normalGradient());
-
-        Info<< "solidWallMixedTemperatureCoupledFvPatchScalarField::"
-            << "updateCoeffs() :"
-            << " patch:" << patch().name()
-            << " heatFlux:" << Q
-            << " walltemperature "
-            << " min:" << gMin(*this)
-            << " max:" << gMax(*this)
-            << " avg:" << gAverage(*this)
-            << endl;
-    }
-
-    forAll(*this, i)
-    {
-        // if outgoing flux use fixed value.
-        if (normalGradient()[i] < 0.0)
-        {
-            this->refValue()[i] = operator[](i);
-            this->refGrad()[i] = 0.0;   // not used by me
-            this->valueFraction()[i] = 1.0;
-            nFixed++;
-        }
-        else
-        {
-            // Fixed gradient. Make sure to have valid refValue (even though
-            // I am not using it - other boundary conditions might)
-            this->refValue()[i] = operator[](i);
-            this->refGrad()[i] = normalGradient()[i];
-            this->valueFraction()[i] = 0.0;
-        }
-    }
-
-    reduce(nFixed, sumOp<label>());
-
-    fixesValue_ = (nFixed > 0);
-
-    if (debug)
-    {
-        label nTotSize = returnReduce(this->size(), sumOp<label>());
-
-        Info<< "solidWallMixedTemperatureCoupledFvPatchScalarField::"
-            << "updateCoeffs() :"
-            << " patch:" << patch().name()
-            << " out of:" << nTotSize
-            << " fixedBC:" << nFixed
-            << " gradient:" << nTotSize-nFixed << endl;
-    }
-
-    mixedFvPatchScalarField::updateCoeffs();
-}
-
-
-void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::write
-(
-    Ostream& os
-) const
-{
-    mixedFvPatchScalarField::write(os);
-    os.writeKeyword("neighbourFieldName")<< neighbourFieldName_
-        << token::END_STATEMENT << nl;
-    os.writeKeyword("K") << KName_ << token::END_STATEMENT << nl;
-    os.writeKeyword("fixesValue") << fixesValue_ << token::END_STATEMENT << nl;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-makePatchTypeField
-(
-    fvPatchScalarField,
-    solidWallMixedTemperatureCoupledFvPatchScalarField
-);
-
-} // End namespace Foam
-
-// ************************************************************************* //
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.H
deleted file mode 100644
index f4766e5dc3..0000000000
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.H
+++ /dev/null
@@ -1,192 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    solidWallMixedTemperatureCoupledFvPatchScalarField
-
-Description
-    Mixed boundary condition for temperature, to be used by the
-    conjugate heat transfer solver.
-    If my temperature is T1, neighbour is T2:
-
-    T1 > T2: my side becomes fixedValue T2 bc, other side becomes fixedGradient.
-    
-
-    Example usage:
-        myInterfacePatchName
-        {
-            type                solidWallMixedTemperatureCoupled;
-            neighbourFieldName  T;
-            K                   K;
-            value               uniform 300;
-        }
-
-    Needs to be on underlying directMapped(Wall)FvPatch.
-
-    Note: runs in parallel with arbitrary decomposition. Uses directMapped
-    functionality to calculate exchange.
-
-    Note: lags interface data so both sides use same data.
-    - problem: schedule to calculate average would interfere
-    with standard processor swaps.
-    - so: updateCoeffs sets both to same Twall. Only need to do
-    this for last outer iteration but don't have access to this.
-
-SourceFiles
-    solidWallMixedTemperatureCoupledFvPatchScalarField.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef solidWallMixedTemperatureCoupledFvPatchScalarField_H
-#define solidWallMixedTemperatureCoupledFvPatchScalarField_H
-
-#include "fvPatchFields.H"
-#include "mixedFvPatchFields.H"
-#include "fvPatch.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-        Class solidWallMixedTemperatureCoupledFvPatchScalarField Declaration
-\*---------------------------------------------------------------------------*/
-
-class solidWallMixedTemperatureCoupledFvPatchScalarField
-:
-    public mixedFvPatchScalarField
-{
-    // Private data
-
-        //- Name of field on the neighbour region
-        const word neighbourFieldName_;
-        
-        //- Name of thermal conductivity field
-        const word KName_;
-
-        bool fixesValue_;
-
-
-    // Private Member Functions
-
-        //- Am I or neighbour owner of interface
-        bool interfaceOwner(const polyMesh& nbrRegion) const;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("solidWallMixedTemperatureCoupled");
-
-
-    // Constructors
-
-        //- Construct from patch and internal field
-        solidWallMixedTemperatureCoupledFvPatchScalarField
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct from patch, internal field and dictionary
-        solidWallMixedTemperatureCoupledFvPatchScalarField
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const dictionary&
-        );
-
-        //- Construct by mapping given
-        //  solidWallMixedTemperatureCoupledFvPatchScalarField onto a new patch
-        solidWallMixedTemperatureCoupledFvPatchScalarField
-        (
-            const solidWallMixedTemperatureCoupledFvPatchScalarField&,
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const fvPatchFieldMapper&
-        );
-
-        //- Construct and return a clone
-        virtual tmp<fvPatchScalarField> clone() const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new solidWallMixedTemperatureCoupledFvPatchScalarField(*this)
-            );
-        }
-
-        //- Construct as copy setting internal field reference
-        solidWallMixedTemperatureCoupledFvPatchScalarField
-        (
-            const solidWallMixedTemperatureCoupledFvPatchScalarField&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct and return a clone setting internal field reference
-        virtual tmp<fvPatchScalarField> clone
-        (
-            const DimensionedField<scalar, volMesh>& iF
-        ) const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new solidWallMixedTemperatureCoupledFvPatchScalarField
-                (
-                    *this,
-                    iF
-                )
-            );
-        }
-
-
-    // Member functions
-
-        //- Get corresponding K field
-        const fvPatchScalarField& K() const;
-
-        //- Return true if this patch field fixes a value.
-        //  Needed to check if a level has to be specified while solving
-        //  Poissons equations.
-        virtual bool fixesValue() const
-        {
-            return fixesValue_;
-        }
-
-        //- Update the coefficients associated with the patch field
-        virtual void updateCoeffs();
-
-        //- Write
-        virtual void write(Ostream&) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

From 3f2b57f0f38262ca6bb5f29551a131b375f28d19 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 1 Apr 2010 14:23:58 +0100
Subject: [PATCH 14/20] BUG: Always do topological walking to correct
 cellZones.

---
 .../meshRefinement/meshRefinementBaffles.C    |  6 ++-
 .../refinementSurfaces/refinementSurfaces.C   | 38 +++++++++++++++++--
 2 files changed, 38 insertions(+), 6 deletions(-)

diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
index e2a2fee10a..9798046ca3 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -2327,9 +2327,11 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
     // Set using walking
     // ~~~~~~~~~~~~~~~~~
 
-    if (!allowFreeStandingZoneFaces)
+    //if (!allowFreeStandingZoneFaces)
     {
-        Info<< "Walking to assign cellZones." << nl << endl;
+        Info<< "Walking from location-in-mesh " << keepPoint
+            << " to assign cellZones "
+            << "- crossing a faceZone face changes cellZone" << nl << endl;
 
         // Topological walk
         findCellZoneTopo
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
index c6aae8ccae..ac01e052cd 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
+++ b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
@@ -77,11 +77,22 @@ Foam::refinementSurfaces::refinementSurfaces
         if (dict.found("faceZone"))
         {
             dict.lookup("faceZone") >> faceZoneNames_[surfI];
+            bool hasSide = dict.readIfPresent("zoneInside", zoneInside_[surfI]);
             if (dict.readIfPresent("cellZone", cellZoneNames_[surfI]))
             {
-                dict.lookup("zoneInside") >> zoneInside_[surfI];
+                if (hasSide && !allGeometry_[surfaces_[surfI]].hasVolumeType())
+                {
+                    IOWarningIn
+                    (
+                        "refinementSurfaces::refinementSurfaces(..)",
+                        dict
+                    )   << "Unused entry zoneInside for faceZone "
+                        << faceZoneNames_[surfI]
+                        << " since surface " << names_[surfI]
+                        << " is not closed." << endl;
+                }
             }
-            else if (dict.found("zoneInside"))
+            else if (hasSide)
             {
                 IOWarningIn("refinementSurfaces::refinementSurfaces(..)", dict)
                     << "Unused entry zoneInside for faceZone "
@@ -324,11 +335,30 @@ Foam::refinementSurfaces::refinementSurfaces
             if (dict.found("faceZone"))
             {
                 dict.lookup("faceZone") >> faceZoneNames_[surfI];
+                bool hasSide = dict.readIfPresent
+                (
+                    "zoneInside",
+                    zoneInside_[surfI]
+                );
                 if (dict.readIfPresent("cellZone", cellZoneNames_[surfI]))
                 {
-                    dict.lookup("zoneInside") >> zoneInside_[surfI];
+                    if
+                    (
+                        hasSide
+                    && !allGeometry_[surfaces_[surfI]].hasVolumeType()
+                    )
+                    {
+                        IOWarningIn
+                        (
+                            "refinementSurfaces::refinementSurfaces(..)",
+                            dict
+                        )   << "Unused entry zoneInside for faceZone "
+                            << faceZoneNames_[surfI]
+                            << " since surface " << names_[surfI]
+                            << " is not closed." << endl;
+                    }
                 }
-                else if (dict.found("zoneInside"))
+                else if (hasSide)
                 {
                     IOWarningIn
                     (

From 25e7d800bb085bb02cead0ca11a396cc61134033 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@Germany>
Date: Tue, 6 Apr 2010 13:46:37 +0200
Subject: [PATCH 15/20] BUG: fix inconsistent ListOps behaviour with
 DynamicList

---
 .../containers/Lists/ListOps/ListOpsTemplates.C | 17 ++++++++++++++++-
 .../primitives/Lists/stringListOpsTemplates.C   |  6 +++++-
 2 files changed, 21 insertions(+), 2 deletions(-)

diff --git a/src/OpenFOAM/containers/Lists/ListOps/ListOpsTemplates.C b/src/OpenFOAM/containers/Lists/ListOps/ListOpsTemplates.C
index d50fbcf0f8..65bd73563d 100644
--- a/src/OpenFOAM/containers/Lists/ListOps/ListOpsTemplates.C
+++ b/src/OpenFOAM/containers/Lists/ListOps/ListOpsTemplates.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,6 +37,9 @@ ListType Foam::renumber
     // Create copy
     ListType newLst(lst.size());
 
+    // ensure consistent addressable size (eg, DynamicList)
+    newLst.setSize(lst.size());
+
     forAll(lst, elemI)
     {
         if (lst[elemI] >= 0)
@@ -76,6 +79,9 @@ ListType Foam::reorder
     // Create copy
     ListType newLst(lst.size());
 
+    // ensure consistent addressable size (eg, DynamicList)
+    newLst.setSize(lst.size());
+
     forAll(lst, elemI)
     {
         if (oldToNew[elemI] >= 0)
@@ -101,6 +107,9 @@ void Foam::inplaceReorder
     // Create copy
     ListType newLst(lst.size());
 
+    // ensure consistent addressable size (eg, DynamicList)
+    newLst.setSize(lst.size());
+
     forAll(lst, elemI)
     {
         if (oldToNew[elemI] >= 0)
@@ -258,6 +267,9 @@ ListType Foam::subset
 
     ListType newLst(lst.size());
 
+    // ensure consistent addressable size (eg, DynamicList)
+    newLst.setSize(lst.size());
+
     label nElem = 0;
     forAll(lst, elemI)
     {
@@ -318,6 +330,9 @@ ListType Foam::subset
 
     ListType newLst(lst.size());
 
+    // ensure consistent addressable size (eg, DynamicList)
+    newLst.setSize(lst.size());
+
     label nElem = 0;
     forAll(lst, elemI)
     {
diff --git a/src/OpenFOAM/primitives/Lists/stringListOpsTemplates.C b/src/OpenFOAM/primitives/Lists/stringListOpsTemplates.C
index 9fdfee1b0a..695d211d61 100644
--- a/src/OpenFOAM/primitives/Lists/stringListOpsTemplates.C
+++ b/src/OpenFOAM/primitives/Lists/stringListOpsTemplates.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -57,8 +57,12 @@ StringListType Foam::subsetMatchingStrings
     const bool invert
 )
 {
+    // Create copy
     StringListType newLst(lst.size());
 
+    // ensure consistent addressable size (eg, DynamicList)
+    newLst.setSize(lst.size());
+
     label nElem = 0;
     forAll(lst, elemI)
     {

From 4c1ff609f928493b4b82a6bd35490a7a5a7a7639 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 6 Apr 2010 16:28:29 +0100
Subject: [PATCH 16/20] STYLE: unnecessary check

---
 .../solid/solidRegionDiffusionNo.H            | 26 +++++++++----------
 1 file changed, 12 insertions(+), 14 deletions(-)

diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H
index 20ca225676..c2e5668a68 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H
@@ -1,18 +1,16 @@
     scalar DiNum = -GREAT;
-    if (solidRegions.size())
+
+    forAll(solidRegions, regionI)
     {
-        forAll(solidRegions, regionI)
-        {
-            DiNum = max
+        DiNum = max
+        (
+            solidRegionDiffNo
             (
-                solidRegionDiffNo
-                (
-                    solidRegions[regionI],
-                    runTime,
-                    rhosCps[regionI],
-                    Ks[regionI]
-                ),
-                DiNum
-            );
-        }
+                solidRegions[regionI],
+                runTime,
+                rhosCps[regionI],
+                Ks[regionI]
+            ),
+            DiNum
+        );
     }

From ed8fb0aa99a82c7155edcb8ae0937730f0ac2313 Mon Sep 17 00:00:00 2001
From: henry <Henry Weller h.weller@opencfd.co.uk>
Date: Tue, 6 Apr 2010 16:28:59 +0100
Subject: [PATCH 17/20] UEqn: Removed "Sp" term from the momentum equation of
 the steady-state solvers.

It is not clear that this "correction" term has provided any benefit and there
are cases for which it has proved detrimental to stability and convergence.
---
 applications/solvers/compressible/rhoPorousSimpleFoam/UEqn.H     | 1 -
 applications/solvers/compressible/rhoSimpleFoam/UEqn.H           | 1 -
 applications/solvers/compressible/rhoSimplecFoam/UEqn.H          | 1 -
 .../solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H      | 1 -
 applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H       | 1 -
 .../solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/UEqn.H   | 1 -
 applications/solvers/incompressible/porousSimpleFoam/UEqn.H      | 1 -
 .../incompressible/MRFSimpleFoam/MRFSimpleFoam/MRFSimpleFoam.C   | 1 -
 8 files changed, 8 deletions(-)

diff --git a/applications/solvers/compressible/rhoPorousSimpleFoam/UEqn.H b/applications/solvers/compressible/rhoPorousSimpleFoam/UEqn.H
index 4150cd7503..5460944288 100644
--- a/applications/solvers/compressible/rhoPorousSimpleFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoPorousSimpleFoam/UEqn.H
@@ -3,7 +3,6 @@
     tmp<fvVectorMatrix> UEqn
     (
         fvm::div(phi, U)
-      - fvm::Sp(fvc::div(phi), U)
       + turbulence->divDevRhoReff(U)
     );
 
diff --git a/applications/solvers/compressible/rhoSimpleFoam/UEqn.H b/applications/solvers/compressible/rhoSimpleFoam/UEqn.H
index c41bc9b6c7..401e1203d2 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/UEqn.H
@@ -3,7 +3,6 @@
     tmp<fvVectorMatrix> UEqn
     (
         fvm::div(phi, U)
-      - fvm::Sp(fvc::div(phi), U)
       + turbulence->divDevRhoReff(U)
     );
 
diff --git a/applications/solvers/compressible/rhoSimplecFoam/UEqn.H b/applications/solvers/compressible/rhoSimplecFoam/UEqn.H
index c41bc9b6c7..401e1203d2 100644
--- a/applications/solvers/compressible/rhoSimplecFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoSimplecFoam/UEqn.H
@@ -3,7 +3,6 @@
     tmp<fvVectorMatrix> UEqn
     (
         fvm::div(phi, U)
-      - fvm::Sp(fvc::div(phi), U)
       + turbulence->divDevRhoReff(U)
     );
 
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H
index a813a2d8b7..bcf99e1668 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H
@@ -3,7 +3,6 @@
     tmp<fvVectorMatrix> UEqn
     (
         fvm::div(phi, U)
-      - fvm::Sp(fvc::div(phi), U)
       + turbulence->divDevReff(U)
     );
 
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H
index aa476b4dc8..0170d6e384 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H
@@ -3,7 +3,6 @@
     tmp<fvVectorMatrix> UEqn
     (
         fvm::div(phi, U)
-      - fvm::Sp(fvc::div(phi), U)
       + turbulence->divDevRhoReff(U)
     );
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/UEqn.H b/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/UEqn.H
index 1cbafe368a..7f317cc2e6 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/UEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/UEqn.H
@@ -2,7 +2,6 @@
     tmp<fvVectorMatrix> UEqn
     (
         fvm::div(phi, U)
-      - fvm::Sp(fvc::div(phi), U)
       + turb.divDevRhoReff(U)
     );
 
diff --git a/applications/solvers/incompressible/porousSimpleFoam/UEqn.H b/applications/solvers/incompressible/porousSimpleFoam/UEqn.H
index c1925c7708..688120809a 100644
--- a/applications/solvers/incompressible/porousSimpleFoam/UEqn.H
+++ b/applications/solvers/incompressible/porousSimpleFoam/UEqn.H
@@ -3,7 +3,6 @@
     tmp<fvVectorMatrix> UEqn
     (
         fvm::div(phi, U)
-      - fvm::Sp(fvc::div(phi), U)
       + turbulence->divDevReff(U)
     );
 
diff --git a/tutorials/incompressible/MRFSimpleFoam/MRFSimpleFoam/MRFSimpleFoam.C b/tutorials/incompressible/MRFSimpleFoam/MRFSimpleFoam/MRFSimpleFoam.C
index c727a74d65..0941f03528 100644
--- a/tutorials/incompressible/MRFSimpleFoam/MRFSimpleFoam/MRFSimpleFoam.C
+++ b/tutorials/incompressible/MRFSimpleFoam/MRFSimpleFoam/MRFSimpleFoam.C
@@ -64,7 +64,6 @@ int main(int argc, char *argv[])
             tmp<fvVectorMatrix> UEqn
             (
                 fvm::div(phi, U)
-              - fvm::Sp(fvc::div(phi), U)
               + turbulence->divDevReff(U)
             );
             mrfZones.addCoriolis(UEqn());

From 7783e806d85e9feff27b4f8a1bb51b0291a00c9f Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 6 Apr 2010 16:30:08 +0100
Subject: [PATCH 18/20] BUG: Was not crossing into unnamed cellZone regions

---
 .../meshRefinement/meshRefinementBaffles.C    | 30 +++++++++++++++++--
 1 file changed, 28 insertions(+), 2 deletions(-)

diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
index 9798046ca3..0d9cecd039 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -1234,7 +1234,7 @@ void Foam::meshRefinement::findCellZoneTopo
         {
             label surfI = namedSurfaceIndex[faceI];
 
-            if (surfI != -1 && surfaceToCellZone[surfI] != -1)
+            if (surfI != -1)
             {
                 // Calculate region to zone from cellRegions on either side
                 // of internal face.
@@ -1286,7 +1286,7 @@ void Foam::meshRefinement::findCellZoneTopo
 
                     label surfI = namedSurfaceIndex[faceI];
 
-                    if (surfI != -1 && surfaceToCellZone[surfI] != -1)
+                    if (surfI != -1)
                     {
                         bool changedCell = calcRegionToZone
                         (
@@ -2522,6 +2522,32 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
         mesh_.setInstance(oldInstance());
     }
 
+    // Print some stats (note: zones are synchronised)
+    if (mesh_.cellZones().size() > 0)
+    {
+        Info<< "CellZones:" << endl;
+        forAll(mesh_.cellZones(), zoneI)
+        {
+            const cellZone& cz = mesh_.cellZones()[zoneI];
+            Info<< "    " << cz.name()
+                << "\tsize:" << returnReduce(cz.size(), sumOp<label>())
+                << endl;
+        }
+        Info<< endl;
+    }
+    if (mesh_.faceZones().size() > 0)
+    {
+        Info<< "FaceZones:" << endl;
+        forAll(mesh_.faceZones(), zoneI)
+        {
+            const faceZone& fz = mesh_.faceZones()[zoneI];
+            Info<< "    " << fz.name()
+                << "\tsize:" << returnReduce(fz.size(), sumOp<label>())
+                << endl;
+        }
+        Info<< endl;
+    }
+
     // None of the faces has changed, only the zones. Still...
     updateMesh(map, labelList());
 

From 645ba8e612e7f5f7d9b51c8497f844262dacaa4d Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 6 Apr 2010 16:30:35 +0100
Subject: [PATCH 19/20] ENH: Split exchange into two step process
 (send+receive) so we can intersperse calculations.

---
 .../mapPolyMesh/mapDistribute/mapDistribute.H |  8 ++
 .../mapDistribute/mapDistributeTemplates.C    | 78 +++++++++++++++++--
 2 files changed, 78 insertions(+), 8 deletions(-)

diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H
index ab7193a682..0d47c7db9d 100644
--- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H
+++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H
@@ -70,6 +70,7 @@ namespace Foam
 
 class mapPolyMesh;
 class globalIndex;
+class PstreamBuffers;
 
 /*---------------------------------------------------------------------------*\
                            Class mapDistribute Declaration
@@ -318,6 +319,13 @@ public:
                 }
             }
 
+            //- Do all sends using PstreamBuffers
+            template<class T>
+            void send(PstreamBuffers&, const List<T>&) const;
+            //- Do all receives using PstreamBuffers
+            template<class T>
+            void receive(PstreamBuffers&, List<T>&) const;
+
             //- Correct for topo change.
             void updateMesh(const mapPolyMesh&)
             {
diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C
index e219bc95ab..8797519d4d 100644
--- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C
+++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C
@@ -185,7 +185,7 @@ void Foam::mapDistribute::distribute
     {
         if (!contiguous<T>())
         {
-            PstreamBuffers pBuffs(Pstream::nonBlocking);
+            PstreamBuffers pBufs(Pstream::nonBlocking);
 
             // Stream data into buffer
             for (label domain = 0; domain < Pstream::nProcs(); domain++)
@@ -195,13 +195,13 @@ void Foam::mapDistribute::distribute
                 if (domain != Pstream::myProcNo() && map.size())
                 {
                     // Put data into send buffer
-                    UOPstream toDomain(domain, pBuffs);
+                    UOPstream toDomain(domain, pBufs);
                     toDomain << UIndirectList<T>(field, map);
                 }
             }
 
             // Start receiving
-            pBuffs.finishedSends();
+            pBufs.finishedSends();
 
             {
                 // Set up 'send' to myself
@@ -231,7 +231,7 @@ void Foam::mapDistribute::distribute
 
                 if (domain != Pstream::myProcNo() && map.size())
                 {
-                    UIPstream str(domain, pBuffs);
+                    UIPstream str(domain, pBufs);
                     List<T> recvField(str);
 
                     if (recvField.size() != map.size())
@@ -551,7 +551,7 @@ void Foam::mapDistribute::distribute
     {
         if (!contiguous<T>())
         {
-            PstreamBuffers pBuffs(Pstream::nonBlocking);
+            PstreamBuffers pBufs(Pstream::nonBlocking);
 
             // Stream data into buffer
             for (label domain = 0; domain < Pstream::nProcs(); domain++)
@@ -561,13 +561,13 @@ void Foam::mapDistribute::distribute
                 if (domain != Pstream::myProcNo() && map.size())
                 {
                     // Put data into send buffer
-                    UOPstream toDomain(domain, pBuffs);
+                    UOPstream toDomain(domain, pBufs);
                     toDomain << UIndirectList<T>(field, map);
                 }
             }
 
             // Start receiving
-            pBuffs.finishedSends();
+            pBufs.finishedSends();
 
             {
                 // Set up 'send' to myself
@@ -597,7 +597,7 @@ void Foam::mapDistribute::distribute
 
                 if (domain != Pstream::myProcNo() && map.size())
                 {
-                    UIPstream str(domain, pBuffs);
+                    UIPstream str(domain, pBufs);
                     List<T> recvField(str);
 
                     if (recvField.size() != map.size())
@@ -757,4 +757,66 @@ void Foam::mapDistribute::distribute
 }
 
 
+template<class T>
+void Foam::mapDistribute::send(PstreamBuffers& pBufs, const List<T>& field)
+const
+{
+    // Stream data into buffer
+    for (label domain = 0; domain < Pstream::nProcs(); domain++)
+    {
+        const labelList& map = subMap_[domain];
+
+        if (map.size())
+        {
+            // Put data into send buffer
+            UOPstream toDomain(domain, pBufs);
+            toDomain << UIndirectList<T>(field, map);
+        }
+    }
+
+    // Start sending and receiving but do not block.
+    pBufs.finishedSends(false);
+}
+
+
+template<class T>
+void Foam::mapDistribute::receive(PstreamBuffers& pBufs, List<T>& field) const
+{
+    // Consume
+    field.setSize(constructSize_);
+
+    for (label domain = 0; domain < Pstream::nProcs(); domain++)
+    {
+        const labelList& map = constructMap_[domain];
+
+        if (map.size())
+        {
+            UIPstream str(domain, pBufs);
+            List<T> recvField(str);
+
+            if (recvField.size() != map.size())
+            {
+                FatalErrorIn
+                (
+                    "template<class T>\n"
+                    "void mapDistribute::receive\n"
+                    "(\n"
+                    "    PstreamBuffers&,\n"
+                    "    List<T>&\n"
+                    ")\n"
+                )   << "Expected from processor " << domain
+                    << " " << map.size() << " but received "
+                    << recvField.size() << " elements."
+                    << abort(FatalError);
+            }
+
+            forAll(map, i)
+            {
+                field[map[i]] = recvField[i];
+            }
+        }
+    }
+}
+
+
 // ************************************************************************* //

From a85b221bf72230ed5b62fbfcf4586e98e9875f11 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 6 Apr 2010 18:45:29 +0100
Subject: [PATCH 20/20] BUG: offsets not initialised consistently with offset.

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

diff --git a/src/meshTools/directMapped/directMappedPolyPatch/directMappedPatchBase.C b/src/meshTools/directMapped/directMappedPolyPatch/directMappedPatchBase.C
index 166d74ae2a..91c24ecf28 100644
--- a/src/meshTools/directMapped/directMappedPolyPatch/directMappedPatchBase.C
+++ b/src/meshTools/directMapped/directMappedPolyPatch/directMappedPatchBase.C
@@ -637,7 +637,7 @@ Foam::directMappedPatchBase::directMappedPatchBase
     samplePatch_(samplePatch),
     uniformOffset_(true),
     offset_(offset),
-    offsets_(0),
+    offsets_(pp.size(), offset_),
     sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
     mapPtr_(NULL)
 {}
@@ -670,7 +670,7 @@ Foam::directMappedPatchBase::directMappedPatchBase
     offsets_
     (
         uniformOffset_
-      ? pointField(patch_.size(), offset_)
+      ? pointField(pp.size(), offset_)
       : dict.lookup("offsets")
     ),
     sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),