Merge branch 'master' of /home/noisy3/OpenFOAM/OpenFOAM-dev

applications/utilities/mesh/manipulation/checkMesh/checkTopology.C

@@ -6,6 +6,7 @@
 #include "faceSet.H"
 #include "pointSet.H"
 #include "IOmanip.H"
+#include "emptyPolyPatch.H"
 
 Foam::label Foam::checkTopology
 (
@@ -21,6 +22,29 @@ Foam::label Foam::checkTopology
     // Check if the boundary definition is unique
     mesh.boundaryMesh().checkDefinition(true);
 
+    // Check that empty patches cover all sides of the mesh
+    {
+        label nEmpty = 0;
+        forAll(mesh.boundaryMesh(), patchI)
+        {
+            if (isA<emptyPolyPatch>(mesh.boundaryMesh()[patchI]))
+            {
+                nEmpty += mesh.boundaryMesh()[patchI].size();
+            }
+        }
+        reduce(nEmpty, sumOp<label>());
+        label nTotCells = returnReduce(mesh.cells().size(), sumOp<label>());
+
+        // These are actually warnings, not errors.
+        if (nEmpty % nTotCells)
+        {
+            Info<< " ***Total number of faces on empty patches"
+                << " is not divisible by the number of cells in the mesh."
+                << " Hence this mesh is not 1D or 2D."
+                << endl;
+        }
+    }
+
     // Check if the boundary processor patches are correct
     mesh.boundaryMesh().checkParallelSync(true);
 
@@ -41,6 +65,8 @@ Foam::label Foam::checkTopology
         noFailedChecks++;
     }
 
+
+
     {
         pointSet points(mesh, "unusedPoints", mesh.nPoints()/100);
         if (mesh.checkPoints(true, &points))
@@ -74,6 +100,22 @@ Foam::label Foam::checkTopology
         }
     }
 
+    {
+        faceSet faces(mesh, "outOfRangeFaces", mesh.nFaces()/100);
+        if (mesh.checkFaceVertices(true, &faces))
+        {
+            noFailedChecks++;
+
+            label nFaces = returnReduce(faces.size(), sumOp<label>());
+
+            Info<< "  <<Writing " << nFaces
+                << " faces with out-of-range or duplicate vertices to set "
+                << faces.name() << endl;
+            faces.instance() = mesh.pointsInstance();
+            faces.write();
+        }
+    }
+
     if (allTopology)
     {
         cellSet cells(mesh, "zipUpCells", mesh.nCells()/100);
@@ -91,22 +133,6 @@ Foam::label Foam::checkTopology
         }
     }
 
-    {
-        faceSet faces(mesh, "outOfRangeFaces", mesh.nFaces()/100);
-        if (mesh.checkFaceVertices(true, &faces))
-        {
-            noFailedChecks++;
-
-            label nFaces = returnReduce(faces.size(), sumOp<label>());
-
-            Info<< "  <<Writing " << nFaces
-                << " faces with out-of-range or duplicate vertices to set "
-                << faces.name() << endl;
-            faces.instance() = mesh.pointsInstance();
-            faces.write();
-        }
-    }
-
     if (allTopology)
     {
         faceSet faces(mesh, "edgeFaces", mesh.nFaces()/100);

src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C

@@ -81,30 +81,6 @@ void Foam::addPatchCellLayer::addVertex
             f[fp++] = pointI;
         }
     }
-
-    // Check for duplicates.
-    if (debug)
-    {
-        label n = 0;
-        for (label i = 0; i < fp; i++)
-        {
-            if (f[i] == pointI)
-            {
-                n++;
-
-                if (n == 2)
-                {
-                    f.setSize(fp);
-                    FatalErrorIn
-                    (
-                        "addPatchCellLayer::addVertex(const label, face&"
-                        ", label&)"
-                    )   << "Point " << pointI << " already present in face "
-                        << f << abort(FatalError);
-                }
-            }
-        }
-    }
 }
 
 
@@ -1641,11 +1617,11 @@ void Foam::addPatchCellLayer::setRefinement
                         {
                             label offset =
                                 addedPoints_[vStart].size() - numEdgeSideFaces;
-                            for (label ioff = offset; ioff > 0; ioff--)
+                            for (label ioff = offset-1; ioff >= 0; ioff--)
                             {
                                 addVertex
                                 (
-                                    addedPoints_[vStart][ioff-1],
+                                    addedPoints_[vStart][ioff],
                                     newFace,
                                     newFp
                                 );
@@ -1660,6 +1636,51 @@ void Foam::addPatchCellLayer::setRefinement
 
                         newFace.setSize(newFp);
 
+                        if (debug)
+                        {
+                            labelHashSet verts(2*newFace.size());
+                            forAll(newFace, fp)
+                            {
+                                if (!verts.insert(newFace[fp]))
+                                {
+                                    FatalErrorIn
+                                    (
+                                        "addPatchCellLayer::setRefinement(..)"
+                                    )   << "Duplicate vertex in face"
+                                        << " to be added." << nl
+                                        << "newFace:" << newFace << nl
+                                        << "points:"
+                                        <<  UIndirectList<point>
+                                            (
+                                                meshMod.points(),
+                                                newFace
+                                            ) << nl
+                                        << "Layer:" << i
+                                        << " out of:" << numEdgeSideFaces << nl
+                                        << "ExtrudeEdge:" << meshEdgeI
+                                        << " at:"
+                                        <<  mesh_.edges()[meshEdgeI].line
+                                            (
+                                                mesh_.points()
+                                            ) << nl
+                                        << "string:" << stringedVerts
+                                        << "stringpoints:"
+                                        << UIndirectList<point>
+                                            (
+                                                pp.localPoints(),
+                                                stringedVerts
+                                            ) << nl
+                                        << "stringNLayers:"
+                                        <<  UIndirectList<label>
+                                            (
+                                                nPointLayers,
+                                                stringedVerts
+                                            ) << nl
+                                        << abort(FatalError);
+                                }
+                            }
+                        }
+
                         label nbrFaceI = nbrFace
                         (
                             pp.edgeFaces(),

src/finiteVolume/fields/fvPatchFields/constraint/empty/emptyFvPatchField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -138,18 +138,21 @@ emptyFvPatchField<Type>::emptyFvPatchField
 template<class Type>
 void emptyFvPatchField<Type>::updateCoeffs()
 {
-    if
+    //- Check moved to checkMesh. Test here breaks down if multiple empty
-    (
+    //  patches.
-        this->patch().patch().size()
+    //if
-      % this->dimensionedInternalField().mesh().nCells()
+    //(
-    )
+    //    this->patch().patch().size()
-    {
+    //  % this->dimensionedInternalField().mesh().nCells()
-        FatalErrorIn("emptyFvPatchField<Type>::updateCoeffs()")
+    //)
-            << "This mesh contains patches of type empty but is not 1D or 2D\n"
+    //{
-               "    by virtue of the fact that the number of faces of this\n"
+    //    FatalErrorIn("emptyFvPatchField<Type>::updateCoeffs()")
-               "    empty patch is not divisible by the number of cells."
+    //        << "This mesh contains patches of type empty but is not"
-            << exit(FatalError);
+    //        << "1D or 2D\n"
-    }
+    //           "    by virtue of the fact that the number of faces of this\n"
+    //           "    empty patch is not divisible by the number of cells."
+    //        << exit(FatalError);
+    //}
 }
 
 

src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C

@@ -129,7 +129,8 @@ Foam::refinementFeatures::refinementFeatures
         }
 
         Info<< "Detected " << featurePoints.size()
-            << " featurePoints out of " << points.size() << endl;
+            << " featurePoints out of " << points.size()
+            << " on feature " << eMesh.name() << endl;
 
         pointTrees_.set
         (
@@ -164,26 +165,30 @@ void Foam::refinementFeatures::findNearestEdge
     forAll(edgeTrees_, featI)
     {
         const indexedOctree<treeDataEdge>& tree = edgeTrees_[featI];
-        forAll(samples, sampleI)
+
+        if (tree.shapes().size() > 0)
         {
-            const point& sample = samples[sampleI];
+            forAll(samples, sampleI)
-
-            scalar distSqr;
-            if (nearInfo[sampleI].hit())
             {
-                distSqr = magSqr(nearInfo[sampleI].hitPoint()-sample);
+                const point& sample = samples[sampleI];
-            }
-            else
-            {
-                distSqr = nearestDistSqr[sampleI];
-            }
 
-            pointIndexHit info = tree.findNearest(sample, distSqr);
+                scalar distSqr;
+                if (nearInfo[sampleI].hit())
+                {
+                    distSqr = magSqr(nearInfo[sampleI].hitPoint()-sample);
+                }
+                else
+                {
+                    distSqr = nearestDistSqr[sampleI];
+                }
 
-            if (info.hit())
+                pointIndexHit info = tree.findNearest(sample, distSqr);
-            {
+
-                nearInfo[sampleI] = info;
+                if (info.hit())
-                nearFeature[sampleI] = featI;
+                {
+                    nearInfo[sampleI] = info;
+                    nearFeature[sampleI] = featI;
+                }
             }
         }
     }
@@ -206,33 +211,37 @@ void Foam::refinementFeatures::findNearestPoint
     forAll(pointTrees_, featI)
     {
         const indexedOctree<treeDataPoint>& tree = pointTrees_[featI];
-        forAll(samples, sampleI)
+
+        if (tree.shapes().pointLabels().size() > 0)
         {
-            const point& sample = samples[sampleI];
+            forAll(samples, sampleI)
-
-            scalar distSqr;
-            if (nearFeature[sampleI] != -1)
             {
-                label nearFeatI = nearFeature[sampleI];
+                const point& sample = samples[sampleI];
-                const indexedOctree<treeDataPoint>& nearTree =
-                    pointTrees_[nearFeatI];
-                label featPointI =
-                    nearTree.shapes().pointLabels()[nearIndex[sampleI]];
-                const point& featPt =
-                    operator[](nearFeatI).points()[featPointI];
-                distSqr = magSqr(featPt-sample);
-            }
-            else
-            {
-                distSqr = nearestDistSqr[sampleI];
-            }
 
-            pointIndexHit info = tree.findNearest(sample, distSqr);
+                scalar distSqr;
+                if (nearFeature[sampleI] != -1)
+                {
+                    label nearFeatI = nearFeature[sampleI];
+                    const indexedOctree<treeDataPoint>& nearTree =
+                        pointTrees_[nearFeatI];
+                    label featPointI =
+                        nearTree.shapes().pointLabels()[nearIndex[sampleI]];
+                    const point& featPt =
+                        operator[](nearFeatI).points()[featPointI];
+                    distSqr = magSqr(featPt-sample);
+                }
+                else
+                {
+                    distSqr = nearestDistSqr[sampleI];
+                }
 
-            if (info.hit())
+                pointIndexHit info = tree.findNearest(sample, distSqr);
-            {
+
-                nearFeature[sampleI] = featI;
+                if (info.hit())
-                nearIndex[sampleI] = info.index();
+                {
+                    nearFeature[sampleI] = featI;
+                    nearIndex[sampleI] = info.index();
+                }
             }
         }
     }

src/regionModels/regionModel/regionModel1D/regionModel1D.C

@@ -114,6 +114,7 @@ void Foam::regionModels::regionModel1D::initialise()
             boundaryFaceCells_[localPyrolysisFaceI].transfer(cellIDs);
 
             localPyrolysisFaceI++;
+            nLayers_ = nCells;
         }
     }
 
@@ -266,6 +267,7 @@ Foam::regionModels::regionModel1D::regionModel1D(const fvMesh& mesh)
     boundaryFaceFaces_(),
     boundaryFaceCells_(),
     boundaryFaceOppositeFace_(),
+    nLayers_(0),
     nMagSfPtr_(NULL),
     moveMesh_(false)
 {}
@@ -283,6 +285,7 @@ Foam::regionModels::regionModel1D::regionModel1D
     boundaryFaceFaces_(regionMesh().nCells()),
     boundaryFaceCells_(regionMesh().nCells()),
     boundaryFaceOppositeFace_(regionMesh().nCells()),
+    nLayers_(0),
     nMagSfPtr_(NULL),
     moveMesh_(true)
 {

src/regionModels/regionModel/regionModel1D/regionModel1D.H

@@ -88,6 +88,9 @@ protected:
             //- Global boundary face IDs oppossite coupled patch
             labelList boundaryFaceOppositeFace_;
 
+            //- Number of layers in the region
+            label nLayers_;
+
 
         // Geometry
 
@@ -157,6 +160,9 @@ public:
 
                 //- Return the face area magnitudes / [m2]
                 inline const surfaceScalarField& nMagSf() const;
+
+                //- Return the number of layers in the region
+                inline label nLayers() const;
 };
 
 

src/regionModels/regionModel/regionModel1D/regionModel1DI.H

@@ -65,4 +65,10 @@ Foam::regionModels::regionModel1D::nMagSf() const
 }
 
 
+inline Foam::label Foam::regionModels::regionModel1D::nLayers() const
+{
+    return nLayers_;
+}
+
+
 // ************************************************************************* //