src/mesh/cfMesh/utilities/meshes/polyMeshGenChecks/polyMeshGenChecksTopology.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | cfMesh: A library for mesh generation
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by the original author
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of cfMesh.
  10
  11     cfMesh is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     cfMesh is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with cfMesh.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "polyMeshGenChecks.H"
  27 #include "polyMeshGenAddressing.H"
  28 #include "cell.H"
  29 #include "Map.H"
  30
  31 # ifdef USE_OMP
  32 #include <omp.h>
  33 # endif
  34
  35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  36
  37 namespace Foam
  38 {
  39
  40 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  41
  42 namespace polyMeshGenChecks
  43 {
  44
  45 bool checkPoints
  46 (
  47     const polyMeshGen& mesh,
  48     const bool report,
  49     labelHashSet* setPtr
  50 )
  51 {
  52     label nFaceErrors = 0;
  53     label nCellErrors = 0;
  54
  55     const VRWGraph& pf = mesh.addressingData().pointFaces();
  56
  57     forAll(pf, pointI)
  58     {
  59         if( pf.sizeOfRow(pointI) == 0 )
  60         {
  61             WarningIn
  62             (
  63                 "bool checkPoints"
  64                 "(const polyMeshGen&, const bool, labelHashSet*)"
  65             )   << "Point " << pointI << " not used by any faces." << endl;
  66
  67             if( setPtr )
  68                 setPtr->insert(pointI);
  69
  70             ++nFaceErrors;
  71         }
  72     }
  73
  74     const VRWGraph& pc = mesh.addressingData().pointCells();
  75
  76     forAll(pc, pointI)
  77     {
  78         if( pc.sizeOfRow(pointI) == 0 )
  79         {
  80             WarningIn
  81             (
  82                 "bool checkPoints"
  83                 "(const polyMeshGen&, const bool, labelHashSet*)"
  84             )   << "Point " << pointI << " not used by any cells." << endl;
  85
  86             if( setPtr )
  87                 setPtr->insert(pointI);
  88
  89             ++nCellErrors;
  90         }
  91     }
  92
  93     reduce(nFaceErrors, sumOp<label>());
  94     reduce(nCellErrors, sumOp<label>());
  95
  96     if( nFaceErrors > 0 || nCellErrors > 0 )
  97     {
  98         WarningIn
  99         (
 100             "bool checkPoints"
 101             "(const polyMeshGen&, const bool, labelHashSet*)"
 102         )   << "Error in point usage detected: " << nFaceErrors
 103             << " unused points found in the mesh.  This mesh is invalid."
 104             << endl;
 105
 106         return true;
 107     }
 108     else
 109     {
 110         if( report )
 111             Info << "Point usage check OK.\n" << endl;
 112
 113         return false;
 114     }
 115 }
 116
 117 bool checkUpperTriangular
 118 (
 119     const polyMeshGen& mesh,
 120     const bool report,
 121     labelHashSet* setPtr
 122 )
 123 {
 124     // Check whether internal faces are ordered in the upper triangular order
 125     const labelList& own = mesh.owner();
 126     const labelList& nei = mesh.neighbour();
 127
 128     const cellListPMG& cells = mesh.cells();
 129     const VRWGraph& cc = mesh.addressingData().cellCells();
 130
 131     const label internal = mesh.nInternalFaces();
 132
 133     labelList checkInternalFaces(internal, -1);
 134
 135     label nChecks = 0;
 136
 137     bool error = false;
 138
 139     // Loop through faceCells once more and make sure that for internal cell
 140     // the first label is smaller
 141     for(label faceI=0;faceI<internal;++faceI)
 142     {
 143         if( own[faceI] >= nei[faceI] )
 144         {
 145             if( report )
 146             {
 147                 Pout<< "bool checkUpperTriangular(const polyMeshGen&, "
 148                     << "const bool, labelHashSet*) : " << endl
 149                     << "face " << faceI
 150                     << " has the owner label greater than neighbour:" << endl
 151                     << own[faceI] << tab << nei[faceI] << endl;
 152             }
 153
 154             if( setPtr )
 155                 setPtr->insert(faceI);
 156
 157             error  = true;
 158         }
 159     }
 160
 161     // Loop through all cells. For each cell, find the face that is internal and
 162     // add it to the check list (upper triangular order).
 163     // Once the list is completed, check it against the faceCell list
 164     forAll(cells, cellI)
 165     {
 166         const labelList& curFaces = cells[cellI];
 167
 168         // Using the fact that cell neighbour always appear
 169         // in the increasing order
 170         boolList usedNbr(cc.sizeOfRow(cellI), false);
 171
 172         for(label nSweeps=0;nSweeps<usedNbr.size();++nSweeps)
 173         {
 174             // Find the lowest neighbour which is still valid
 175             label nextNei = -1;
 176             label minNei = cells.size();
 177
 178             forAllRow(cc, cellI, nbrI)
 179             {
 180                 const label neiI = cc(cellI, nbrI);
 181                 if( (neiI > cellI) && !usedNbr[nbrI] && (neiI < minNei) )
 182                 {
 183                     nextNei = nbrI;
 184                     minNei = neiI;
 185                 }
 186             }
 187
 188             if( nextNei > -1 )
 189             {
 190                 // Mark this neighbour as used
 191                 usedNbr[nextNei] = true;
 192
 193                 forAll(curFaces, faceI)
 194                 {
 195                     if( curFaces[faceI] < internal )
 196                     {
 197                         if( nei[curFaces[faceI]] == cc(cellI, nextNei) )
 198                         {
 199                             checkInternalFaces[nChecks] = curFaces[faceI];
 200                             ++nChecks;
 201
 202                             break;
 203                         }
 204                     }
 205                 }
 206             }
 207         }
 208     }
 209
 210     // Check list created. If everything is OK, the face label is equal to index
 211     forAll(checkInternalFaces, faceI)
 212     {
 213         if( checkInternalFaces[faceI] != faceI )
 214         {
 215             error = true;
 216
 217             Pout<< "bool checkUpperTriangular(const polyMeshGen&, const bool"
 218                 << ", labelHashSet*) : " << endl
 219                 << "face " << faceI << " out of position. Markup label: "
 220                 << checkInternalFaces[faceI] << ". All subsequent faces will "
 221                 << "also be out of position. Please check the mesh manually."
 222                 << endl;
 223
 224             if( setPtr )
 225                 setPtr->insert(faceI);
 226
 227             break;
 228         }
 229     }
 230
 231     reduce(error, orOp<bool>());
 232
 233     if( error )
 234     {
 235         WarningIn
 236         (
 237             "bool checkUpperTriangular(const polyMeshGen&, const bool"
 238             ", labelHashSet*)"
 239         )   << "Error in face ordering: faces not in upper triangular order!"
 240             << endl;
 241
 242         return true;
 243     }
 244     else
 245     {
 246         if( report )
 247             Info<< "Upper triangular ordering OK.\n" << endl;
 248
 249         return false;
 250     }
 251 }
 252
 253 bool checkCellsZipUp
 254 (
 255     const polyMeshGen& mesh,
 256     const bool report,
 257     labelHashSet* setPtr
 258 )
 259 {
 260     label nOpenCells = 0;
 261
 262     const faceListPMG& faces = mesh.faces();
 263     const cellListPMG& cells = mesh.cells();
 264
 265     # ifdef USE_OMP
 266     # pragma omp parallel for schedule(guided) reduction(+ : nOpenCells)
 267     # endif
 268     forAll(cells, cellI)
 269     {
 270         const labelList& c = cells[cellI];
 271
 272         DynList<edge> cellEdges;
 273         DynList<label> edgeUsage;
 274
 275         forAll(c, faceI)
 276         {
 277             const face& f = faces[c[faceI]];
 278
 279             forAll(f, eI)
 280             {
 281                 const edge e = f.faceEdge(eI);
 282
 283                 const label pos = cellEdges.containsAtPosition(e);
 284
 285                 if( pos < 0 )
 286                 {
 287                     cellEdges.append(e);
 288                     edgeUsage.append(1);
 289                 }
 290                 else
 291                 {
 292                     ++edgeUsage[pos];
 293                 }
 294             }
 295         }
 296
 297         DynList<edge> singleEdges;
 298
 299         forAll(edgeUsage, edgeI)
 300         {
 301             if( edgeUsage[edgeI] == 1 )
 302             {
 303                 singleEdges.append(cellEdges[edgeI]);
 304             }
 305             else if( edgeUsage[edgeI] != 2 )
 306             {
 307                 WarningIn
 308                 (
 309                     "bool checkCellsZipUp(const polyMeshGen&,"
 310                     "const bool, labelHashSet*)"
 311                 )   << "edge " << cellEdges[edgeI] << " in cell " << cellI
 312                     << " used " << edgeUsage[edgeI] << " times. " << endl
 313                     << "Should be 1 or 2 - serious error in mesh structure"
 314                     << endl;
 315
 316                 if( setPtr )
 317                 {
 318                     # ifdef USE_OMP
 319                     # pragma omp critical
 320                     # endif
 321                     setPtr->insert(cellI);
 322                 }
 323             }
 324         }
 325
 326         if( singleEdges.size() > 0 )
 327         {
 328             if( report )
 329             {
 330                 Pout<< "bool checkCellsZipUp(const polyMeshGen&, const bool"
 331                     << ", labelHashSet*) : " << endl
 332                     << "Cell " << cellI << " has got " << singleEdges.size()
 333                     << " unmatched edges: " << singleEdges << endl;
 334             }
 335
 336             if( setPtr )
 337             {
 338                 # ifdef USE_OMP
 339                 # pragma omp critical
 340                 # endif
 341                 setPtr->insert(cellI);
 342             }
 343
 344             ++nOpenCells;
 345         }
 346     }
 347
 348     reduce(nOpenCells, sumOp<label>());
 349
 350     if( nOpenCells > 0 )
 351     {
 352         WarningIn
 353         (
 354             "bool checkCellsZipUp(const polyMeshGen&,"
 355             " const bool, labelHashSet*)"
 356         )   << nOpenCells
 357             << " open cells found.  Please use the mesh zip-up tool. "
 358             << endl;
 359
 360         return true;
 361     }
 362     else
 363     {
 364         if( report )
 365             Info<< "Topological cell zip-up check OK.\n" << endl;
 366
 367         return false;
 368     }
 369 }
 370
 371 // Vertices of face within point range and unique.
 372 bool checkFaceVertices
 373 (
 374     const polyMeshGen& mesh,
 375     const bool report,
 376     labelHashSet* setPtr
 377 )
 378 {
 379     // Check that all vertex labels are valid
 380     const faceListPMG& faces = mesh.faces();
 381
 382     label nErrorFaces = 0;
 383     const label nPoints = mesh.points().size();
 384
 385     forAll(faces, fI)
 386     {
 387         const face& curFace = faces[fI];
 388
 389         if( min(curFace) < 0 || max(curFace) > nPoints )
 390         {
 391             WarningIn
 392             (
 393                 "bool checkFaceVertices("
 394                 "const polyMesgGen&, const bool, labelHashSet*)"
 395             )   << "Face " << fI << " contains vertex labels out of range: "
 396                 << curFace << " Max point index = " << nPoints-1 << endl;
 397
 398             if( setPtr )
 399                 setPtr->insert(fI);
 400
 401             ++nErrorFaces;
 402         }
 403
 404         // Uniqueness of vertices
 405         labelHashSet facePoints(2*curFace.size());
 406
 407         forAll(curFace, fp)
 408         {
 409             bool inserted = facePoints.insert(curFace[fp]);
 410
 411             if( !inserted )
 412             {
 413                 WarningIn
 414                 (
 415                     "bool checkFaceVertices("
 416                     "const polyMeshGen&, const bool, labelHashSet*)"
 417                 )   << "Face " << fI << " contains duplicate vertex labels: "
 418                     << curFace << endl;
 419
 420                 if( setPtr )
 421                     setPtr->insert(fI);
 422
 423                 ++nErrorFaces;
 424             }
 425         }
 426     }
 427
 428     reduce(nErrorFaces, sumOp<label>());
 429
 430     if( nErrorFaces > 0 )
 431     {
 432         SeriousErrorIn
 433         (
 434             "bool checkFaceVertices("
 435             "const polyMeshGen&, const bool, labelHashSet*)"
 436         )   << "const bool, labelHashSet*) const: "
 437             << nErrorFaces << " faces with invalid vertex labels found"
 438             << endl;
 439
 440         return true;
 441     }
 442     else
 443     {
 444         if( report )
 445             Info<< "Face vertices OK.\n" << endl;
 446
 447         return false;
 448     }
 449 }
 450
 451 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 452
 453 } // End namespace polyMeshGenChecks
 454
 455 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 456
 457 } // End namespace Foam
 458
 459 // ************************************************************************* //