Merge branch 'defect-release-v1.1' into developmentPublicRepo

meshLibrary/utilities/containers/LongList/LongList.H

@@ -41,7 +41,6 @@ SourceFiles
 #define LongList_H
 
 #include "label.H"
-#include "longLong.H"
 #include "bool.H"
 #include "IOstreams.H"
 #include "error.H"

meshLibrary/utilities/meshes/polyMeshGenChecks/polyMeshGenChecksGeometry.C

@@ -367,7 +367,7 @@ bool checkFaceAreas
     const boolList* changedFacePtr
 )
 {
-    const scalarField magFaceAreas = mag(mesh.addressingData().faceAreas());
+    const vectorField& faceAreas = mesh.addressingData().faceAreas();
 
     const labelList& own = mesh.owner();
     const labelList& nei = mesh.neighbour();
@@ -384,12 +384,14 @@ bool checkFaceAreas
         # ifdef USE_OMP
         # pragma omp for schedule(guided)
         # endif
-        forAll(magFaceAreas, faceI)
+        forAll(faceAreas, faceI)
         {
             if( changedFacePtr && !changedFacePtr->operator[](faceI) )
                 continue;
 
-            if( magFaceAreas[faceI] < minFaceArea )
+            const scalar magFaceArea = mag(faceAreas[faceI]);
+
+            if( magFaceArea < minFaceArea )
             {
                 if( report )
                 {
@@ -399,14 +401,14 @@ bool checkFaceAreas
                             << "internal face " << faceI << " between cells "
                             << own[faceI] << " and " << nei[faceI]
                             << ".  Face area magnitude = "
-                            << magFaceAreas[faceI] << endl;
+                            << magFaceArea << endl;
                     }
                     else
                     {
                         Pout<< "Zero or negative face area detected for "
                             << "boundary face " << faceI << " next to cell "
                             << own[faceI] << ".  Face area magnitude = "
-                            << magFaceAreas[faceI] << endl;
+                            << magFaceArea << endl;
                     }
                 }
 
@@ -419,8 +421,8 @@ bool checkFaceAreas
                 }
             }
 
-            localMinArea = Foam::min(localMinArea, magFaceAreas[faceI]);
-            localMaxArea = Foam::max(localMaxArea, magFaceAreas[faceI]);
+            localMinArea = Foam::min(localMinArea, magFaceArea);
+            localMaxArea = Foam::max(localMaxArea, magFaceArea);
         }
 
         # ifdef USE_OMP

meshLibrary/utilities/meshes/triSurf/triSurf.C

@@ -32,6 +32,8 @@ Description
 #include "gzstream.h"
 #include "triSurface.H"
 
+#include "helperFunctions.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -163,6 +165,192 @@ void triSurf::writeToFMS(const fileName& fName) const
     fStream << subsets;
 }
 
+void triSurf::topologyCheck()
+{
+    const pointField& pts = this->points();
+    const LongList<labelledTri>& trias = this->facets();
+
+    //- check for inf and nan points
+    # ifdef USE_OMP
+    # pragma omp parallel for schedule(dynamic, 100)
+    # endif
+    forAll(pts, pointI)
+    {
+        const point& p = pts[pointI];
+
+        if( help::isnan(p) || help::isinf(p) )
+        {
+            # ifdef USE_OMP
+            # pragma omp critical
+            # endif
+            {
+                FatalErrorIn
+                (
+                    "void triSurf::topologyCheck()"
+                ) << "Point " << pointI << " has invalid coordinates "
+                  << p << exit(FatalError);
+            }
+        }
+    }
+
+    //- check whether the nodes are within the scope
+    //- report duplicate nodes in the same triangle
+    # ifdef USE_OMP
+    # pragma omp parallel for schedule(dynamic, 100)
+    # endif
+    forAll(trias, triI)
+    {
+        const labelledTri& ltri = trias[triI];
+
+        forAll(ltri, pI)
+        {
+            if( ltri[pI] < 0 || (ltri[pI] >= pts.size()) )
+            {
+                # ifdef USE_OMP
+                # pragma omp critical
+                # endif
+                FatalErrorIn
+                (
+                    "void triSurf::topologyCheck()"
+                ) << "Point " << ltri[pI] << " in triangle " << ltri
+                  << " is out of scope 0 " << pts.size() << exit(FatalError);
+            }
+
+            if( ltri[pI] == ltri[(pI+1)%3] || ltri[pI] == ltri[(pI+2)%3] )
+            {
+                # ifdef USE_OMP
+                # pragma omp critical
+                # endif
+                WarningIn
+                (
+                    "void triSurf::topologyCheck()"
+                ) << "Triangle " << ltri << " has duplicated points. "
+                  << "This may cause problems in the meshing process!" << endl;
+            }
+        }
+    }
+
+    //- check feature edges
+    const edgeLongList& featureEdges = this->featureEdges();
+
+    # ifdef USE_OMP
+    # pragma omp parallel for schedule(dynamic, 100)
+    # endif
+    forAll(featureEdges, eI)
+    {
+        const edge& fe = featureEdges[eI];
+
+        forAll(fe, pI)
+        {
+            if( fe[pI] < 0 || (fe[pI] >= pts.size()) )
+            {
+                # ifdef USE_OMP
+                # pragma omp critical
+                # endif
+                FatalErrorIn
+                (
+                    "void triSurf::topologyCheck()"
+                ) << "Feature edge " << fe << " point " << fe[pI]
+                  << " is out of scope 0 " << pts.size() << exit(FatalError);
+            }
+        }
+
+        if( fe.start() == fe.end() )
+        {
+            # ifdef USE_OMP
+            # pragma omp critical
+            # endif
+            WarningIn
+            (
+                "void triSurf::topologyCheck()"
+            ) << "Feature edge " << fe << " has duplicated points. "
+              << "This may cause problems in the meshing process!" << endl;
+        }
+    }
+
+    //- check point subsets
+    DynList<label> subsetIds;
+    this->pointSubsetIndices(subsetIds);
+    forAll(subsetIds, i)
+    {
+        labelLongList elmts;
+        this->pointsInSubset(subsetIds[i], elmts);
+
+        forAll(elmts, elmtI)
+        {
+            const label elI = elmts[elmtI];
+
+            if( elI < 0 || elI >= pts.size() )
+            {
+                # ifdef USE_OMP
+                # pragma omp critical
+                # endif
+                FatalErrorIn
+                (
+                    "void triSurf::topologyCheck()"
+                ) << "Point " << elI << " in point subset "
+                  << this->pointSubsetName(subsetIds[i])
+                  << " is out of scope 0 " << pts.size() << exit(FatalError);
+            }
+        }
+    }
+
+    //- check face subsets
+    subsetIds.clear();
+    this->facetSubsetIndices(subsetIds);
+    forAll(subsetIds, i)
+    {
+        labelLongList elmts;
+        this->facetsInSubset(subsetIds[i], elmts);
+
+        forAll(elmts, elmtI)
+        {
+            const label elI = elmts[elmtI];
+
+            if( elI < 0 || elI >= trias.size() )
+            {
+                # ifdef USE_OMP
+                # pragma omp critical
+                # endif
+                FatalErrorIn
+                (
+                    "void triSurf::topologyCheck()"
+                ) << "Triangle " << elI << " in facet subset "
+                  << this->facetSubsetName(subsetIds[i])
+                  << " is out of scope 0 " << trias.size() << exit(FatalError);
+            }
+        }
+    }
+
+    //- check feature edge subsets
+    subsetIds.clear();
+    this->edgeSubsetIndices(subsetIds);
+    forAll(subsetIds, i)
+    {
+        labelLongList elmts;
+        this->edgesInSubset(subsetIds[i], elmts);
+
+        forAll(elmts, elmtI)
+        {
+            const label elI = elmts[elmtI];
+
+            if( elI < 0 || elI >= featureEdges.size() )
+            {
+                # ifdef USE_OMP
+                # pragma omp critical
+                # endif
+                FatalErrorIn
+                (
+                    "void triSurf::topologyCheck()"
+                ) << "Feature edge " << elI << " in edge subset "
+                  << this->edgeSubsetName(subsetIds[i])
+                  << " is out of scope 0 " << featureEdges.size()
+                  << exit(FatalError);
+            }
+        }
+    }
+}
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 triSurf::triSurf()
@@ -186,7 +374,9 @@ triSurf::triSurf
     triSurfFacets(triangles, patches),
     triSurfFeatureEdges(featureEdges),
     triSurfAddressing(triSurfPoints::points_, triSurfFacets::triangles_)
-{}
+{
+    topologyCheck();
+}
 
 //- Read from file
 triSurf::triSurf(const fileName& fName)
@@ -197,6 +387,8 @@ triSurf::triSurf(const fileName& fName)
     triSurfAddressing(triSurfPoints::points_, triSurfFacets::triangles_)
 {
     readSurface(fName);
+
+    topologyCheck();
 }
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

meshLibrary/utilities/meshes/triSurf/triSurf.H

@@ -74,6 +74,8 @@ class triSurf
         inline LongList<labelledTri>& accessToFacets();
         inline geometricSurfacePatchList& accessToPatches();
 
+        void topologyCheck();
+
         //- Disallow default bitwise assignment
         void operator=(const triSurf&);
 

meshLibrary/utilities/meshes/triSurf/triSurfAddressing.C

@@ -29,6 +29,8 @@ Description
 #include "VRWGraphSMPModifier.H"
 #include "demandDrivenData.H"
 
+#include <set>
+
 # ifdef USE_OMP
 #include <omp.h>
 # endif
@@ -70,49 +72,42 @@ void triSurfAddressing::calculateEdges() const
         # ifdef USE_OMP
         # pragma omp for schedule(static)
         # endif
-        forAll(facets_, fI)
+        forAll(pFacets, pI)
         {
-            const labelledTri& tri = facets_[fI];
+            std::set<std::pair<label, label> > edgesAtPoint;
 
-            forAll(tri, pI)
+            forAllRow(pFacets, pI, pfI)
             {
-                const edge fe(tri[pI], tri[(pI+1)%3]);
-                const label s = fe.start();
-                const label e = fe.end();
-
-                DynList<label> edgeFaces;
-
-                bool store(true);
+                const label triI = pFacets(pI, pfI);
+                const labelledTri& tri = facets_[triI];
 
-                //- find all faces attached to this edge
-                //- store the edge in case the face faceI is the face
-                //- with the smallest label
-                forAllRow(pFacets, s, pfI)
+                label pos(-1);
+                forAll(tri, i)
                 {
-                    const label ofI = pFacets(s, pfI);
-                    const labelledTri& of = facets_[ofI];
-
-                    label pos(-1);
-                    for(label i=0;i<3;++i)
-                        if( of[i] == e )
-                        {
-                            pos = i;
-                            break;
-                        }
-                    if( pos < 0 )
-                        continue;
-                    if( ofI < fI )
+                    if( tri[i] == pI )
                     {
-                        store = false;
-                        break;
+                        pos = i;
+
+
+                        if( tri[(pos+1)%3] >= pI )
+                            edgesAtPoint.insert
+                            (
+                                std::make_pair(pI, tri[(pos+1)%3])
+                            );
+                        if( tri[(pos+2)%3] >= pI )
+                            edgesAtPoint.insert
+                            (
+                                std::make_pair(pI, tri[(pos+2)%3])
+                            );
                     }
-
-                    edgeFaces.append(ofI);
                 }
 
-                if( store )
-                    edgesHelper.append(fe);
+
             }
+
+            std::set<std::pair<label, label> >::const_iterator it;
+            for(it=edgesAtPoint.begin();it!=edgesAtPoint.end();++it)
+                edgesHelper.append(edge(it->first, it->second));
         }
 
         //- this enables other threads to see the number of edges
@@ -158,56 +153,58 @@ void triSurfAddressing::calculateFacetEdges() const
     const edgeLongList& edges = this->edges();
     const VRWGraph& pointFaces = this->pointFacets();
 
-    facetEdgesPtr_ = new VRWGraph(facets_.size());
+    facetEdgesPtr_ = new VRWGraph(facets_.size(), 3, -1);
     VRWGraph& faceEdges = *facetEdgesPtr_;
 
-    labelList nfe(facets_.size());
-
     # ifdef USE_OMP
-    const label nThreads = 3 * omp_get_num_procs();
-    # pragma omp parallel num_threads(nThreads)
+    # pragma omp parallel for schedule(dynamic, 100)
     # endif
+    forAll(edges, edgeI)
     {
-        # ifdef USE_OMP
-        # pragma omp for schedule(static)
-        # endif
-        forAll(facets_, fI)
-            nfe[fI] = 3;
-
-        # ifdef USE_OMP
-        # pragma omp barrier
-
-        # pragma omp master
-        # endif
-        VRWGraphSMPModifier(faceEdges).setSizeAndRowSize(nfe);
-
-        # ifdef USE_OMP
-        # pragma omp barrier
+        const edge ee = edges[edgeI];
+        const label pI = ee.start();
 
-        # pragma omp for schedule(static)
-        # endif
-        forAll(edges, edgeI)
+        forAllRow(pointFaces, pI, pfI)
         {
-            const edge ee = edges[edgeI];
-            const label pI = ee.start();
+            const label fI = pointFaces(pI, pfI);
 
-            forAllRow(pointFaces, pI, pfI)
+            const labelledTri& tri = facets_[fI];
+            forAll(tri, eI)
             {
-                const label fI = pointFaces(pI, pfI);
+                const edge e(tri[eI], tri[(eI+1)%3]);
 
-                const labelledTri& tri = facets_[fI];
-                forAll(tri, eI)
+                if( e == ee )
                 {
-                    const edge e(tri[eI], tri[(eI+1)%3]);
-                    if( e == ee )
-                    {
-                        faceEdges(fI, eI) = edgeI;
-                        break;
-                    }
+                    faceEdges(fI, eI) = edgeI;
                 }
             }
         }
     }
+
+    # ifdef DEBUGTriSurfAddressing
+    forAll(faceEdges, triI)
+    {
+        forAllRow(faceEdges, triI, feI)
+        {
+            if( facets_[triI][feI] < 0 || facets_[triI][feI] >= points_.size() )
+                FatalErrorIn
+                (
+                    "void triSurfAddressing::calculateFacetEdges() const"
+                ) << "Invalid entry in triangle " << triI
+                  << " " << facets_[triI] << abort(FatalError);
+
+            const label edgeI = faceEdges(triI, feI);
+
+            if( edgeI < 0 || edgeI >= edges.size() )
+                FatalErrorIn
+                (
+                    "void triSurfAddressing::calculateFacetEdges() const"
+                ) << "Invalid entry in face " << triI << " "
+                     << facets_[triI] << " edges "
+                     << faceEdges[triI] << abort(FatalError);
+        }
+    }
+    # endif
 }
 
 void triSurfAddressing::calculateEdgeFacets() const

meshLibrary/utilities/octrees/meshOctree/meshOctreeAddressing/meshOctreeAddressingCreation.C

@@ -39,6 +39,10 @@ Description
 
 //#define DEBUGVrt
 
+# ifdef DEBUGVrt
+#include "OFstream.H"
+# endif
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -46,6 +50,76 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+# ifdef DEBUGVrt
+void writeOctreeToVTK
+(
+    const fileName& fName,
+    const meshOctree& octree,
+    const List<direction>& boxTypes,
+    const direction bType
+)
+{
+    OFstream file(fName);
+
+    //- write the header
+    file << "# vtk DataFile Version 3.0\n";
+    file << "vtk output\n";
+    file << "ASCII\n";
+    file << "DATASET UNSTRUCTURED_GRID\n";
+
+    label nBoxes(0);
+    forAll(boxTypes, leafI)
+        if( boxTypes[leafI] & bType )
+            ++nBoxes;
+
+    //- write points
+    file << "POINTS " << 8*nBoxes << " float\n";
+    forAll(boxTypes, leafI)
+    {
+        if( boxTypes[leafI] & bType )
+        {
+            FixedList<point, 8> vertices;
+            octree.returnLeaf(leafI).vertices(octree.rootBox(), vertices);
+
+            forAll(vertices, vI)
+            {
+                const point& p = vertices[vI];
+
+                file << p.x() << ' ' << p.y() << ' ' << p.z() << nl;
+            }
+        }
+    }
+
+    //- write boxes
+    file << "\nCELLS " << nBoxes
+         << " " << 9 * nBoxes << nl;
+
+    nBoxes = 0;
+    forAll(boxTypes, leafI)
+    {
+        if( boxTypes[leafI] & bType )
+        {
+            const label start = 8 * nBoxes;
+            file << 8 << " " << start << " " << start+1
+                      << " " << start+3 << " " << start+2
+                      << " " << start+4 << " " << start+5
+                      << " " << start+7 << " " << start+6 << nl;
+
+            ++nBoxes;
+        }
+    }
+
+    file << nl;
+
+    //- write cell types
+    file << "CELL_TYPES " << nBoxes << nl;
+    for(label i=0;i<nBoxes;++i)
+        file << 12 << nl;
+
+    file << nl;
+}
+# endif
+
 void meshOctreeAddressing::createOctreePoints() const
 {
     const VRWGraph& nodeLabels = this->nodeLabels();
@@ -56,7 +130,7 @@ void meshOctreeAddressing::createOctreePoints() const
 
     const label nLeaves = nodeLabels.size();
     # ifdef USE_OMP
-    # pragma omp parallel for schedule(guided)
+    # pragma omp parallel for schedule(guided, 100)
     # endif
     for(label cubeI=0;cubeI<nLeaves;++cubeI)
     {
@@ -66,10 +140,12 @@ void meshOctreeAddressing::createOctreePoints() const
         FixedList<point, 8> vertices;
         const meshOctreeCubeBasic& oc = octree_.returnLeaf(cubeI);
         oc.vertices(rootBox, vertices);
-        forAllRow(nodeLabels, cubeI, pI)
+
+        forAllRow(nodeLabels, cubeI, nI)
         {
-            forAllRow(nodeLabels, cubeI, nI)
-                octreePoints[nodeLabels(cubeI, nI)] = vertices[nI];
+            const label nodeI = nodeLabels(cubeI, nI);
+
+            octreePoints[nodeI] = vertices[nI];
         }
     }
 }
@@ -97,7 +173,7 @@ void meshOctreeAddressing::createNodeLabels() const
     nNodes_ = 0;
     DynList<label> numLocalNodes;
     # ifdef USE_OMP
-    # pragma omp parallel //num_threads(Foam::max(nodeLabels.size() / 1000, 1))
+    # pragma omp parallel
     # endif
     {
         # ifdef USE_OMP
@@ -162,7 +238,7 @@ void meshOctreeAddressing::createNodeLabels() const
                     }
                 }
 
-                if( minLeaf == leafI && validLeaf[7-nI] )
+                if( (minLeaf == leafI) && validLeaf[7-nI] )
                 {
                     forAll(pLeaves, plI)
                         if( validLeaf[plI] )
@@ -251,6 +327,74 @@ void meshOctreeAddressing::createNodeLabels() const
             nNodes_ = Foam::max(nNodes_, startNode);
         }
     }
+
+    # ifdef DEBUGVrt
+    List<direction> badLeaves(nodeLabels.size(), direction(0));
+    forAll(nodeLabels, leafI)
+        forAllRow(nodeLabels, leafI, i)
+            if( nodeLabels(leafI, i) < 0 )
+                badLeaves[leafI] |= 1;
+    writeOctreeToVTK("badLeaves.vtk", octree_, badLeaves, 1);
+
+    writeOctreeToVTK("meshCells.vtk", octree_, boxType, MESHCELL);
+    writeOctreeToVTK("boundaryCells.vtk", octree_, boxType, BOUNDARY);
+