Merge branch 'master' into particleInteractions

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;

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)")
+        fvc::div(phi/fvc::interpolate(rho), rho/psi, "div(U,p)")
-      - p*fvc::div(phi/fvc::interpolate(rho))
+      - (rho/psi)*fvc::div(phi/fvc::interpolate(rho))
     );
 
     hEqn.relax();

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;

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;
 }

src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H

@@ -291,6 +291,9 @@ class primitiveMesh
             //- Skewness warning threshold
             static scalar skewThreshold_;
 
+            //- Threshold where faces are considered coplanar
+            static scalar planarCosAngle_;
+
 
 protected:
 
@@ -646,6 +649,13 @@ public:
                     labelHashSet* setPtr = NULL
                 ) const;
 
+                //- Check for concave cells by the planes of faces
+                bool checkConcaveCells
+                (
+                    const bool report = false,
+                    labelHashSet* setPtr = NULL
+                ) const;
+
 
             //- Check mesh topology for correctness.
             //  Returns false for no error.

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C

@@ -37,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  * * * * * * * * * * * //
@@ -2093,6 +2094,121 @@ 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();
+
+    label nConcaveCells = 0;
+
+    forAll (c, cellI)
+    {
+        const cell& cFaces = c[cellI];
+
+        bool concave = false;
+
+        forAll(cFaces, i)
+        {
+            if (concave)
+            {
+                break;
+            }
+
+            label fI = cFaces[i];
+
+            const point& fC = fCentres[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 the centre of any other face of the cell on the
+            // wrong side of the plane of this face?
+
+            forAll(cFaces, j)
+            {
+                if (j != i)
+                {
+                    label fJ = cFaces[j];
+
+                    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.
+
+                    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

src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.C

@@ -202,12 +202,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
     (

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C

@@ -1375,7 +1375,7 @@ void Foam::autoLayerDriver::growNoExtrusion
     pointField& patchDisp,
     labelList& patchNLayers,
     List<extrudeMode>& extrudeStatus
-)
+) const
 {
     Info<< nl << "Growing non-extrusion points by one layer ..." << endl;
 
@@ -1406,7 +1406,7 @@ void Foam::autoLayerDriver::growNoExtrusion
             {
                 if
                 (
-                    extrudeStatus[f[fp]] == NOEXTRUDE
+                    extrudeStatus[f[fp]] == EXTRUDE
                  && grownExtrudeStatus[f[fp]] != NOEXTRUDE
                 )
                 {
@@ -1419,6 +1419,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)
@@ -2711,8 +2736,6 @@ void Foam::autoLayerDriver::addLayers
 
             const labelList& meshPoints = patches[patchI].meshPoints();
 
-            //scalar maxThickness = -VGREAT;
-            //scalar minThickness = VGREAT;
             scalar sumThickness = 0;
             scalar sumNearWallThickness = 0;
 
@@ -2720,8 +2743,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];
@@ -2749,8 +2770,6 @@ void Foam::autoLayerDriver::addLayers
 
             if (totNPoints > 0)
             {
-                //reduce(maxThickness, maxOp<scalar>());
-                //reduce(minThickness, minOp<scalar>());
                 avgThickness =
                     returnReduce(sumThickness, sumOp<scalar>())
                   / totNPoints;
@@ -2766,8 +2785,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;
@@ -3147,6 +3164,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;
     }
 
@@ -3293,7 +3323,6 @@ void Foam::autoLayerDriver::doLayers
 
 
         // Balance
-        if (Pstream::parRun())
         if (Pstream::parRun() && preBalance)
         {
             Info<< nl

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H

@@ -232,13 +232,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

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C

@@ -245,27 +245,8 @@ void Foam::meshRefinement::getBafflePatches
 
     const pointField& cellCentres = mesh_.cellCentres();
 
-    // Build list of surfaces that are not to be baffled.
+    // Surfaces that need to be baffled
-    const wordList& cellZoneNames = surfaces_.cellZoneNames();
+    const labelList surfacesToBaffle(surfaces_.getUnnamedSurfaces());
-
-    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);
 
     ownPatch.setSize(mesh_.nFaces());
     ownPatch = -1;
@@ -1253,7 +1234,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.
@@ -1305,7 +1286,7 @@ void Foam::meshRefinement::findCellZoneTopo
 
                     label surfI = namedSurfaceIndex[faceI];
 
-                    if (surfI != -1)
+                    if (surfI != -1 && surfaceToCellZone[surfI] != -1)
                     {
                         bool changedCell = calcRegionToZone
                         (
@@ -1439,6 +1420,15 @@ void Foam::meshRefinement::makeConsistentFaceIndex
                 }
             }
         }
+        else
+        {
+            // Unzonify boundary faces
+            forAll(pp, i)
+            {
+                label faceI = pp.start()+i;
+                namedSurfaceIndex[faceI] = -1;
+            }
+        }
     }
 }
 
@@ -2042,14 +2032,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;
@@ -2125,31 +2111,34 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
         {
             label surfI = namedSurfaces[i];
 
-            label zoneI = cellZones.findZoneID(cellZoneNames[surfI]);
+            if (cellZoneNames[surfI] != word::null)
-
-            if (zoneI == -1)
             {
-                zoneI = cellZones.size();
+                label zoneI = cellZones.findZoneID(cellZoneNames[surfI]);
-                cellZones.setSize(zoneI+1);
+
-                cellZones.set
+                if (zoneI == -1)
-                (
+                {
-                    zoneI,
+                    zoneI = cellZones.size();
-                    new cellZone
+                    cellZones.setSize(zoneI+1);
+                    cellZones.set
                     (
-                        cellZoneNames[surfI],   //name
+                        zoneI,
-                        labelList(0),           //addressing
+                        new cellZone
-                        zoneI,                  //index
+                        (
-                        cellZones               //cellZoneMesh
+                            cellZoneNames[surfI],   //name
-                    )
+                            labelList(0),           //addressing
-                );
+                            zoneI,                  //index
-            }
+                            cellZones               //cellZoneMesh
+                        )
+                    );
+                }
 
-            if (debug)
+                if (debug)
-            {
+                {
-                Pout<< "Cells inside " << surfaces_.names()[surfI]
+                    Pout<< "Cells inside " << surfaces_.names()[surfI]
-                    << " will go into cellZone " << zoneI << endl;
+                        << " will go into cellZone " << zoneI << endl;
+                }
+                surfaceToCellZone[surfI] = zoneI;
             }
-            surfaceToCellZone[surfI] = zoneI;
         }
 
         // Check they are synced
@@ -2321,6 +2310,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
@@ -2333,8 +2327,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
         (
@@ -2349,6 +2345,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);
     }
 

src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C

@@ -77,8 +77,18 @@ Foam::refinementSurfaces::refinementSurfaces
         if (dict.found("faceZone"))
         {
             dict.lookup("faceZone") >> faceZoneNames_[surfI];
-            dict.lookup("cellZone") >> cellZoneNames_[surfI];
+            if (dict.readIfPresent("cellZone", cellZoneNames_[surfI]))
-            dict.lookup("zoneInside") >> zoneInside_[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
@@ -314,8 +324,21 @@ Foam::refinementSurfaces::refinementSurfaces
             if (dict.found("faceZone"))
             {
                 dict.lookup("faceZone") >> faceZoneNames_[surfI];
-                dict.lookup("cellZone") >> cellZoneNames_[surfI];
+                if (dict.readIfPresent("cellZone", cellZoneNames_[surfI]))
-                dict.lookup("zoneInside") >> zoneInside_[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
@@ -475,18 +498,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(cellZoneNames_, surfI)
-    forAll(named, i)
     {
-        label surfI = named[i];
+        if (cellZoneNames_[surfI].size())
-
-        if (allGeometry_[surfaces_[surfI]].hasVolumeType())
         {
-            closed[closedI++] = surfI;
+            if (allGeometry_[surfaces_[surfI]].hasVolumeType())
+            {
+                closed[closedI++] = surfI;
+            }
         }
     }
     closed.setSize(closedI);

src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H

@@ -146,7 +146,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_;
@@ -158,7 +159,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

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