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

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | cfMesh: A library for mesh generation
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by the origina 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 "pyramidPointFaceRef.H"
  29 #include "tetrahedron.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 checkClosedBoundary(const polyMeshGen& mesh, const bool report)
  46 {
  47     // Loop through all boundary faces and sum up the face area vectors.
  48     // For a closed boundary, this should be zero in all vector components
  49     vector sumClosed(vector::zero);
  50     scalar sumMagClosedBoundary = 0;
  51
  52     const vectorField& areas = mesh.addressingData().faceAreas();
  53
  54     for(label faceI = mesh.nInternalFaces();faceI<areas.size();++faceI)
  55     {
  56         sumClosed += areas[faceI];
  57         sumMagClosedBoundary += mag(areas[faceI]);
  58     }
  59
  60     // Check the openness in the maximum direction (this is APPROXIMATE!)
  61     scalar maxOpen = max
  62     (
  63         sumClosed.component(vector::X),
  64         max
  65         (
  66             sumClosed.component(vector::Y),
  67             sumClosed.component(vector::Z)
  68         )
  69     );
  70
  71     reduce(sumClosed, sumOp<vector>());
  72     reduce(maxOpen, maxOp<scalar>());
  73
  74     if( maxOpen > SMALL*max(1.0, sumMagClosedBoundary) )
  75     {
  76         SeriousErrorIn
  77         (
  78             "bool checkClosedBoundary(const polyMeshGen&, const bool report)"
  79         )   << "Possible hole in boundary description" << endl;
  80
  81         Info<< "Boundary openness in x-direction = "
  82             << sumClosed.component(vector::X) << endl;
  83
  84         Info<< "Boundary openness in y-direction = "
  85             << sumClosed.component(vector::Y) << endl;
  86
  87         Info<< "Boundary openness in z-direction = "
  88             << sumClosed.component(vector::Z) << endl;
  89
  90         return true;
  91     }
  92     else
  93     {
  94         if( report )
  95         {
  96             Info<< "Boundary openness in x-direction = "
  97                 << sumClosed.component(vector::X) << endl;
  98
  99             Info<< "Boundary openness in y-direction = "
 100                 << sumClosed.component(vector::Y) << endl;
 101
 102             Info<< "Boundary openness in z-direction = "
 103                 << sumClosed.component(vector::Z) << endl;
 104
 105             Info<< "Boundary closed (OK)." << endl;
 106         }
 107
 108         return false;
 109     }
 110 }
 111
 112 bool checkClosedCells
 113 (
 114     const polyMeshGen& mesh,
 115     const bool report,
 116     const scalar aspectWarn,
 117     labelHashSet* setPtr
 118 )
 119 {
 120     // Check that all cells labels are valid
 121     const cellListPMG& cells = mesh.cells();
 122     const label nFaces = mesh.faces().size();
 123
 124     label nErrorClosed = 0;
 125
 126     # ifdef USE_OMP
 127     # pragma omp parallel for schedule(guided) reduction(+ : nErrorClosed)
 128     # endif
 129     forAll(cells, cI)
 130     {
 131         const cell& curCell = cells[cI];
 132
 133         if( min(curCell) < 0 || max(curCell) > nFaces )
 134         {
 135             WarningIn
 136             (
 137                 "bool checkClosedCells("
 138                 "const polyMeshGen&, const bool, const scalar, labelHashSet*)"
 139             )   << "Cell " << cI << " contains face labels out of range: "
 140                 << curCell << " Max face index = " << nFaces << endl;
 141
 142             if( setPtr )
 143             {
 144                 # ifdef USE_OMP
 145                 # pragma omp critical
 146                 # endif
 147                 setPtr->insert(cI);
 148             }
 149
 150             ++nErrorClosed;
 151         }
 152     }
 153
 154     if( nErrorClosed > 0 )
 155     {
 156         SeriousErrorIn
 157         (
 158             "bool checkClosedCells("
 159             "const polyMeshGen&, const bool, const scalar, labelHashSet*)"
 160         )  << nErrorClosed << " cells with invalid face labels found"
 161             << endl;
 162
 163         return true;
 164     }
 165
 166     // Loop through cell faces and sum up the face area vectors for each cell.
 167     // This should be zero in all vector components
 168     vectorField sumClosed(cells.size(), vector::zero);
 169
 170     scalarField sumMagClosed(cells.size(), 0.0);
 171
 172     const labelList& own = mesh.owner();
 173     const labelList& nei = mesh.neighbour();
 174     const vectorField& areas = mesh.addressingData().faceAreas();
 175
 176     forAll(own, faceI)
 177     {
 178         // Add to owner
 179         sumClosed[own[faceI]] += areas[faceI];
 180         sumMagClosed[own[faceI]] += mag(areas[faceI]);
 181
 182         if( nei[faceI] == -1 )
 183             continue;
 184
 185         // Subtract from neighbour
 186         sumClosed[nei[faceI]] -= areas[faceI];
 187         sumMagClosed[nei[faceI]] += mag(areas[faceI]);
 188     }
 189
 190     label nOpen = 0;
 191     scalar maxOpenCell = 0;
 192
 193     label nAspect = 0;
 194     scalar maxAspectRatio = 0;
 195
 196     const scalarField& vols = mesh.addressingData().cellVolumes();
 197
 198     // Check the sums
 199     forAll(sumClosed, cellI)
 200     {
 201         scalar maxOpen = max
 202         (
 203             sumClosed[cellI].component(vector::X),
 204             max
 205             (
 206                 sumClosed[cellI].component(vector::Y),
 207                 sumClosed[cellI].component(vector::Z)
 208             )
 209         );
 210
 211         maxOpenCell = max(maxOpenCell, maxOpen);
 212
 213         if( maxOpen > SMALL*max(1.0, sumMagClosed[cellI]) )
 214         {
 215             if( report )
 216             {
 217                 Pout<< "Cell " << cellI << " is not closed. "
 218                     << "Face area vectors sum up to " << mag(sumClosed[cellI])
 219                     << " directionwise " << sumClosed[cellI] << " or "
 220                     << mag(sumClosed[cellI])
 221                        /(mag(sumMagClosed[cellI]) + VSMALL)
 222                     << " of the area of all faces of the cell. " << endl
 223                     << "    Area magnitudes sum up to "
 224                     << sumMagClosed[cellI] << endl;
 225             }
 226
 227             if( setPtr )
 228             {
 229                 setPtr->insert(cellI);
 230             }
 231
 232             ++nOpen;
 233         }
 234
 235         scalar aspectRatio =
 236             1.0/6.0*sumMagClosed[cellI]/pow(vols[cellI], 2.0/3.0);
 237
 238         maxAspectRatio = max(maxAspectRatio, aspectRatio);
 239
 240         if( aspectRatio > aspectWarn )
 241         {
 242             if( report )
 243             {
 244                 Pout<< "High aspect ratio for cell " << cellI
 245                     << ": " << aspectRatio << endl;
 246             }
 247
 248             ++nAspect;
 249         }
 250     }
 251
 252     reduce(nOpen, sumOp<label>());
 253     reduce(maxOpenCell, maxOp<scalar>());
 254
 255     reduce(nAspect, sumOp<label>());
 256     reduce(maxAspectRatio, maxOp<scalar>());
 257
 258     if( nOpen > 0 )
 259     {
 260         SeriousErrorIn
 261         (
 262             "bool checkClosedCells("
 263             "const polyMeshGen&, const bool, const scalar, labelHashSet*)"
 264         )   << nOpen << " open cells found. Max cell openness: "
 265             << maxOpenCell << endl;
 266
 267         return true;
 268     }
 269
 270     if( nAspect > 0 )
 271     {
 272         SeriousErrorIn
 273         (
 274             "bool checkClosedCells("
 275             "const polyMeshGen&, const bool, const scalar, labelHashSet*)"
 276         )   << nAspect << " high aspect ratio cells found.  "
 277             << "Max aspect ratio: " << maxAspectRatio
 278             << endl;
 279
 280         return true;
 281     }
 282
 283     if( report )
 284     {
 285         Info<< "Max cell openness = " << maxOpenCell
 286             << "  Max aspect ratio = " << maxAspectRatio
 287             << ".  All cells OK.\n" << endl;
 288     }
 289
 290     return false;
 291 }
 292
 293 bool checkCellVolumes
 294 (
 295     const polyMeshGen& mesh,
 296     const bool report,
 297     labelHashSet* setPtr
 298 )
 299 {
 300     const scalarField& vols = mesh.addressingData().cellVolumes();
 301
 302     scalar minVolume = GREAT;
 303     scalar maxVolume = -GREAT;
 304
 305     label nNegVolCells = 0;
 306
 307     forAll(vols, cellI)
 308     {
 309         if( vols[cellI] < VSMALL )
 310         {
 311             if( report )
 312                 SeriousErrorIn
 313                 (
 314                     "bool checkCellVolumes("
 315                     "const polyMeshGen&, const bool, labelHashSet*)"
 316                 )   << "Zero or negative cell volume detected for cell "
 317                     << cellI << ".  Volume = " << vols[cellI] << endl;
 318
 319             if( setPtr )
 320                 setPtr->insert(cellI);
 321
 322             ++nNegVolCells;
 323         }
 324
 325         minVolume = min(minVolume, vols[cellI]);
 326         maxVolume = max(maxVolume, vols[cellI]);
 327     }
 328
 329     reduce(minVolume, minOp<scalar>());
 330     reduce(maxVolume, maxOp<scalar>());
 331     reduce(nNegVolCells, sumOp<label>());
 332
 333     if( minVolume < VSMALL )
 334     {
 335         SeriousErrorIn
 336         (
 337             "bool checkCellVolumes("
 338             "const polyMeshGen&, const bool, labelHashSet*)"
 339         )   << "Zero or negative cell volume detected.  "
 340             << "Minimum negative volume: "
 341             << minVolume << ".\nNumber of negative volume cells: "
 342             << nNegVolCells << ".  This mesh is invalid"
 343             << endl;
 344
 345         return true;
 346     }
 347     else
 348     {
 349         if( report )
 350         {
 351             Info<< "Min volume = " << minVolume
 352                 << ". Max volume = " << maxVolume
 353                 << ".  Total volume = " << sum(vols)
 354                 << ".  Cell volumes OK.\n" << endl;
 355         }
 356
 357         return false;
 358     }
 359 }
 360
 361 bool checkFaceAreas
 362 (
 363     const polyMeshGen& mesh,
 364     const bool report,
 365     const scalar minFaceArea,
 366     labelHashSet* setPtr,
 367     const boolList* changedFacePtr
 368 )
 369 {
 370     const vectorField& faceAreas = mesh.addressingData().faceAreas();
 371
 372     const labelList& own = mesh.owner();
 373     const labelList& nei = mesh.neighbour();
 374
 375     scalar minArea = VGREAT;
 376     scalar maxArea = -VGREAT;
 377
 378     # ifdef USE_OMP
 379     # pragma omp parallel if( own.size() > 100 )
 380     # endif
 381     {
 382         scalar localMaxArea(-VGREAT), localMinArea(VGREAT);
 383
 384         # ifdef USE_OMP
 385         # pragma omp for schedule(guided)
 386         # endif
 387         forAll(faceAreas, faceI)
 388         {
 389             if( changedFacePtr && !changedFacePtr->operator[](faceI) )
 390                 continue;
 391
 392             const scalar magFaceArea = mag(faceAreas[faceI]);
 393
 394             if( magFaceArea < minFaceArea )
 395             {
 396                 if( report )
 397                 {
 398                     if( nei[faceI] != -1 )
 399                     {
 400                         Pout<< "Zero or negative face area detected for "
 401                             << "internal face " << faceI << " between cells "
 402                             << own[faceI] << " and " << nei[faceI]
 403                             << ".  Face area magnitude = "
 404                             << magFaceArea << endl;
 405                     }
 406                     else
 407                     {
 408                         Pout<< "Zero or negative face area detected for "
 409                             << "boundary face " << faceI << " next to cell "
 410                             << own[faceI] << ".  Face area magnitude = "
 411                             << magFaceArea << endl;
 412                     }
 413                 }
 414
 415                 if( setPtr )
 416                 {
 417                     # ifdef USE_OMP
 418                     # pragma omp critical
 419                     # endif
 420                     setPtr->insert(faceI);
 421                 }
 422             }
 423
 424             localMinArea = Foam::min(localMinArea, magFaceArea);
 425             localMaxArea = Foam::max(localMaxArea, magFaceArea);
 426         }
 427
 428         # ifdef USE_OMP
 429         # pragma omp critical
 430         # endif
 431         {
 432             minArea = Foam::min(minArea, localMinArea);
 433             maxArea = Foam::max(maxArea, localMaxArea);
 434         }
 435     }
 436
 437     reduce(minArea, minOp<scalar>());
 438     reduce(maxArea, maxOp<scalar>());
 439
 440     if( minArea < VSMALL )
 441     {
 442         SeriousErrorIn
 443         (
 444             "bool checkFaceAreas("
 445             "const polyMeshGen&, const bool, const scalar,"
 446             " labelHashSet*, const boolList*)"
 447         )   << "Zero or negative face area detected.  Minimum negative area: "
 448             << minArea << ". This mesh is invalid"
 449             << endl;
 450
 451         return true;
 452     }
 453     else
 454     {
 455         if( report )
 456         {
 457             Info<< "Minumum face area = " << minArea
 458                 << ". Maximum face area = " << maxArea
 459                 << ".  Face area magnitudes OK.\n" << endl;
 460         }
 461
 462         return false;
 463     }
 464 }
 465
 466 bool checkCellPartTetrahedra
 467 (
 468     const polyMeshGen& mesh,
 469     const bool report,
 470     const scalar minPartTet,
 471     labelHashSet* setPtr,
 472     const boolList* changedFacePtr
 473 )
 474 {
 475     const pointFieldPMG& points = mesh.points();
 476     const faceListPMG& faces = mesh.faces();
 477     const labelList& owner = mesh.owner();
 478     const labelList& neighbour = mesh.neighbour();
 479
 480     const vectorField& fCentres = mesh.addressingData().faceCentres();
 481     const vectorField& cCentres = mesh.addressingData().cellCentres();
 482
 483     label nNegVolCells = 0;
 484
 485     # ifdef USE_OMP
 486     # pragma omp parallel for if( owner.size() > 100 ) \
 487     schedule(guided) reduction(+ : nNegVolCells)
 488     # endif
 489     forAll(owner, faceI)
 490     {
 491         if( changedFacePtr && !(*changedFacePtr)[faceI] )
 492             continue;
 493
 494         const face& f = faces[faceI];
 495
 496         bool badFace(false);
 497
 498         forAll(f, eI)
 499         {
 500             const tetrahedron<point, point> tetOwn
 501             (
 502                 fCentres[faceI],
 503                 points[f.nextLabel(eI)],
 504                 points[f[eI]],
 505                 cCentres[owner[faceI]]
 506             );
 507
 508             if( tetOwn.mag() < minPartTet )
 509             {
 510                 if( report )
 511                 {
 512                     # ifdef USE_OMP
 513                     # pragma omp critical
 514                     # endif
 515                     Pout<< "Zero or negative cell volume detected for cell "
 516                         << owner[faceI] << "." << endl;
 517                 }
 518
 519                 badFace = true;
 520             }
 521
 522             if( neighbour[faceI] < 0 )
 523                 continue;
 524
 525             const tetrahedron<point, point> tetNei
 526             (
 527                 fCentres[faceI],
 528                 points[f[eI]],
 529                 points[f.nextLabel(eI)],
 530                 cCentres[neighbour[faceI]]
 531             );
 532
 533             if( tetNei.mag() < minPartTet )
 534             {
 535                 if( report )
 536                 {
 537                     # ifdef USE_OMP
 538                     # pragma omp critical
 539                     # endif
 540                     Pout<< "Zero or negative cell volume detected for cell "
 541                         << neighbour[faceI] << "." << endl;
 542                 }
 543
 544                 badFace = true;
 545             }
 546         }
 547
 548         if( badFace )
 549         {
 550             if( setPtr )
 551             {
 552                 # ifdef USE_OMP
 553                 # pragma omp critical
 554                 # endif
 555                 setPtr->insert(faceI);
 556             }
 557
 558             ++nNegVolCells;
 559         }
 560     }
 561
 562     if( setPtr )
 563     {
 564         //- ensure that faces are selected at both sides
 565         const PtrList<processorBoundaryPatch>& procBnd = mesh.procBoundaries();
 566         forAll(procBnd, patchI)
 567         {
 568             const label start = procBnd[patchI].patchStart();
 569             const label size = procBnd[patchI].patchSize();
 570
 571             labelLongList sendData;
 572             for(label faceI=0;faceI<size;++faceI)
 573             {
 574                 if( setPtr->found(faceI+start) )
 575                     sendData.append(faceI);
 576             }
 577
 578             OPstream toOtherProc
 579             (
 580                 Pstream::blocking,
 581                 procBnd[patchI].neiProcNo(),
 582                 sendData.byteSize()
 583             );
 584
 585             toOtherProc << sendData;
 586         }
 587
 588         forAll(procBnd, patchI)
 589         {
 590             labelList receivedData;
 591
 592             IPstream fromOtherProc
 593             (
 594                 Pstream::blocking,
 595                 procBnd[patchI].neiProcNo()
 596             );
 597
 598             fromOtherProc >> receivedData;
 599
 600             const label start = procBnd[patchI].patchStart();
 601             forAll(receivedData, i)
 602                 setPtr->insert(start+receivedData[i]);
 603         }
 604     }
 605
 606     reduce(nNegVolCells, sumOp<label>());
 607
 608     if( nNegVolCells != 0 )
 609     {
 610         WarningIn
 611         (
 612             "bool checkCellPartTetrahedra("
 613             "const polyMeshGen&, const bool, const scalar,"
 614             " labelHashSet*, const boolList*)"
 615         )   << nNegVolCells << " zero or negative part tetrahedra detected."
 616             << endl;
 617
 618         return true;
 619     }
 620     else
 621     {
 622         if( report )
 623             Info<< "Part cell tetrahedra OK.\n" << endl;
 624
 625         return false;
 626     }
 627 }
 628
 629 void checkFaceDotProduct
 630 (
 631     const polyMeshGen& mesh,
 632     scalarField& faceDotProduct,
 633     const boolList* changedFacePtr
 634 )
 635 {
 636     //- for all internal faces check theat the d dot S product is positive
 637     const vectorField& centres = mesh.addressingData().cellCentres();
 638     const vectorField& areas = mesh.addressingData().faceAreas();
 639
 640     const labelList& own = mesh.owner();
 641     const labelList& nei = mesh.neighbour();
 642     const label nInternalFaces = mesh.nInternalFaces();
 643
 644     faceDotProduct.setSize(own.size());
 645     faceDotProduct = 1.0;
 646
 647     # ifdef USE_OMP
 648     # pragma omp parallel if( nInternalFaces > 1000 )
 649     # endif
 650     {
 651         # ifdef USE_OMP
 652         # pragma omp for schedule(dynamic, 100)
 653         # endif
 654         for(label faceI=0;faceI<nInternalFaces;++faceI)
 655         {
 656             if( changedFacePtr && !(*changedFacePtr)[faceI] )
 657                 continue;
 658
 659             const vector d = centres[nei[faceI]] - centres[own[faceI]];
 660             const vector& s = areas[faceI];
 661
 662             faceDotProduct[faceI] = (d & s)/(mag(d)*mag(s) + VSMALL);
 663         }
 664     }
 665
 666     if( Pstream::parRun() )
 667     {
 668         const PtrList<processorBoundaryPatch>& procBoundaries =
 669             mesh.procBoundaries();
 670
 671         forAll(procBoundaries, patchI)
 672         {
 673             const label start = procBoundaries[patchI].patchStart();
 674
 675             vectorField cCentres(procBoundaries[patchI].patchSize());
 676             forAll(cCentres, faceI)
 677                 cCentres[faceI] = centres[own[start+faceI]];
 678
 679             OPstream toOtherProc
 680             (
 681                 Pstream::blocking,
 682                 procBoundaries[patchI].neiProcNo(),
 683                 cCentres.byteSize()
 684             );
 685
 686             toOtherProc << cCentres;
 687         }
 688
 689         forAll(procBoundaries, patchI)
 690         {
 691             vectorField otherCentres;
 692             IPstream fromOtherProc
 693             (
 694                 Pstream::blocking,
 695                 procBoundaries[patchI].neiProcNo()
 696             );
 697
 698             fromOtherProc >> otherCentres;
 699
 700             //- calculate skewness at processor faces
 701             const label start = procBoundaries[patchI].patchStart();
 702             # ifdef USE_OMP
 703             # pragma omp parallel
 704             # endif
 705             {
 706                 # ifdef USE_OMP
 707                 # pragma omp for schedule(dynamic, 100)
 708                 # endif
 709                 forAll(otherCentres, fI)
 710                 {
 711                     const label faceI = start + fI;
 712
 713                     if( changedFacePtr && !(*changedFacePtr)[faceI] )
 714                         continue;
 715
 716                     const point& cOwn = centres[own[faceI]];
 717                     const point& cNei = otherCentres[fI];
 718                     const vector d = cNei - cOwn;
 719                     const vector& s = areas[faceI];
 720
 721                     faceDotProduct[faceI] = (d & s)/(mag(d)*mag(s) + VSMALL);
 722                 }
 723             }
 724         }
 725     }
 726 }
 727
 728 bool checkFaceDotProduct
 729 (
 730     const polyMeshGen& mesh,
 731     const bool report,
 732     const scalar nonOrthWarn,
 733     labelHashSet* setPtr,
 734     const boolList* changedFacePtr
 735 )
 736 {
 737     //- for all internal faces check theat the d dot S product is positive
 738     scalarField faceDotProduct;
 739     checkFaceDotProduct(mesh, faceDotProduct, changedFacePtr);
 740
 741     const labelList& own = mesh.owner();
 742     const labelList& nei = mesh.neighbour();
 743     const label nInternalFaces = mesh.nInternalFaces();
 744
 745     // Severe nonorthogonality threshold
 746     const scalar severeNonorthogonalityThreshold =
 747         ::cos(nonOrthWarn/180.0*M_PI);
 748
 749     scalar minDDotS = VGREAT;
 750     scalar sumDDotS = 0.0;
 751
 752     label severeNonOrth = 0;
 753     label errorNonOrth = 0;
 754     label counter = 0;
 755
 756     # ifdef USE_OMP
 757     # pragma omp parallel if( nInternalFaces > 1000 ) \
 758     reduction(+ : severeNonOrth, errorNonOrth, sumDDotS)
 759     # endif
 760     {
 761         scalar localMinDDotS(VGREAT);
 762         # ifdef USE_OMP
 763         # pragma omp for schedule(guided)
 764         # endif
 765         for(label faceI=0;faceI<nInternalFaces;++faceI)
 766         {
 767             if( changedFacePtr && !(*changedFacePtr)[faceI] )
 768                 continue;
 769
 770             const scalar dDotS = faceDotProduct[faceI];
 771
 772             if( dDotS < severeNonorthogonalityThreshold )
 773             {
 774                 if( dDotS > SMALL )
 775                 {
 776                     if( report )
 777                     {
 778                         // Severe non-orthogonality but mesh still OK
 779                         # ifdef USE_OMP
 780                         # pragma omp critical
 781                         # endif
 782                         Pout<< "Severe non-orthogonality for face " << faceI
 783                             << " between cells " << own[faceI]
 784                             << " and " << nei[faceI]
 785                             << ": Angle = "
 786                             << ::acos(dDotS)/M_PI*180.0
 787                             << " deg." << endl;
 788                     }
 789
 790                     if( setPtr )
 791                     {
 792                         # ifdef USE_OMP
 793                         # pragma omp critical
 794                         # endif
 795                         setPtr->insert(faceI);
 796                     }
 797
 798                     ++severeNonOrth;
 799                 }
 800                 else
 801                 {
 802                     ++errorNonOrth;
 803
 804                     if( setPtr )
 805                     {
 806                         # ifdef USE_OMP
 807                         # pragma omp critical
 808                         # endif
 809                         setPtr->insert(faceI);
 810                     }
 811                 }
 812             }
 813
 814             localMinDDotS = Foam::min(dDotS, localMinDDotS);
 815             sumDDotS += dDotS;
 816         }
 817
 818         # ifdef USE_OMP
 819         # pragma omp critical
 820         # endif
 821         minDDotS = Foam::min(minDDotS, localMinDDotS);
 822     }
 823
 824     counter += nInternalFaces;
 825
 826     if( Pstream::parRun() )
 827     {
 828         const PtrList<processorBoundaryPatch>& procBoundaries =
 829             mesh.procBoundaries();
 830
 831         const label start = procBoundaries[0].patchStart();
 832
 833         # ifdef USE_OMP
 834         # pragma omp parallel \
 835         reduction(+ : severeNonOrth, errorNonOrth, sumDDotS, counter)
 836         # endif
 837         {
 838             scalar localMinDDotS(VGREAT);
 839
 840             # ifdef USE_OMP
 841             # pragma omp for schedule(dynamic, 100)
 842             # endif
 843             for(label faceI=start;faceI<own.size();++faceI)
 844             {
 845                 if( changedFacePtr && !(*changedFacePtr)[start+faceI] )
 846                             continue;
 847
 848                 const scalar dDotS = faceDotProduct[faceI];
 849
 850                 if( dDotS < severeNonorthogonalityThreshold )
 851                 {
 852                     if( dDotS > SMALL )
 853                     {
 854                         if( report )
 855                         {
 856                             // Severe non-orthogonality but mesh still OK
 857                             # ifdef USE_OMP
 858                             # pragma omp critical
 859                             # endif
 860                             {
 861                                 const scalar angle
 862                                 (
 863                                     Foam::acos(dDotS) /
 864                                     M_PI * 180.0
 865                                 );
 866                                 Pout<< "Severe non-orthogonality for face "
 867                                     << start+faceI
 868                                     << ": Angle = "
 869                                     << angle
 870                                     << " deg." << endl;
 871                             }
 872                         }
 873
 874                         if( setPtr )
 875                         {
 876                             # ifdef USE_OMP
 877                             # pragma omp critical
 878                             # endif
 879                             setPtr->insert(start+faceI);
 880                         }
 881
 882                         ++severeNonOrth;
 883                     }
 884                     else
 885                     {
 886                         ++errorNonOrth;
 887
 888                         if( setPtr )
 889                         {
 890                             # ifdef USE_OMP
 891                             # pragma omp critical
 892                             # endif
 893                             setPtr->insert(start+faceI);
 894                         }
 895                     }
 896                 }
 897
 898                 localMinDDotS = Foam::min(dDotS, localMinDDotS);
 899                 sumDDotS += 0.5 * dDotS;
 900
 901                 ++counter;
 902             }
 903
 904             # ifdef USE_OMP
 905             # pragma omp critical
 906             # endif
 907             minDDotS = Foam::min(minDDotS, localMinDDotS);
 908         }
 909     }
 910
 911     reduce(minDDotS, minOp<scalar>());
 912     reduce(sumDDotS, sumOp<scalar>());
 913     reduce(severeNonOrth, sumOp<label>());
 914     reduce(errorNonOrth, sumOp<label>());
 915
 916     reduce(counter, sumOp<label>());
 917
 918     // Only report if there are some internal faces
 919     if( counter > 0 )
 920     {
 921         if( minDDotS < severeNonorthogonalityThreshold )
 922         {
 923             Info<< "Number of non-orthogonality errors: " << errorNonOrth
 924                 << ". Number of severely non-orthogonal faces: "
 925                 << severeNonOrth  << "." << endl;
 926         }
 927     }
 928
 929     if( report )
 930     {
 931         if( counter > 0 )
 932         {
 933             Info<< "Mesh non-orthogonality Max: "
 934                 << ::acos(minDDotS)/M_PI*180.0
 935                 << " average: " <<
 936                    ::acos(sumDDotS/counter)/M_PI*180.0
 937                 << endl;
 938         }
 939     }
 940
 941     if( errorNonOrth > 0 )
 942     {
 943         WarningIn
 944         (
 945             "checkFaceDotProduct("
 946             "const polyMeshGen&, const bool, const scalar,"
 947             " labelHashSet*, const boolList*)"
 948         )   << "Error in non-orthogonality detected" << endl;
 949
 950         return true;
 951     }
 952     else
 953     {
 954         if( report )
 955             Info<< "Non-orthogonality check OK.\n" << endl;
 956
 957         return false;
 958     }
 959 }
 960
 961 bool checkFacePyramids
 962 (
 963     const polyMeshGen& mesh,
 964     const bool report,
 965     const scalar minPyrVol,
 966     labelHashSet* setPtr,
 967     const boolList* changedFacePtr
 968 )
 969 {
 970     // check whether face area vector points to the cell with higher label
 971     const vectorField& ctrs = mesh.addressingData().cellCentres();
 972
 973     const labelList& owner = mesh.owner();
 974     const labelList& neighbour = mesh.neighbour();
 975
 976     const faceListPMG& faces = mesh.faces();
 977     const pointFieldPMG& points = mesh.points();
 978
 979     label nErrorPyrs = 0;
 980
 981     # ifdef USE_OMP
 982     # pragma omp parallel for schedule(guided) reduction(+ : nErrorPyrs)
 983     # endif
 984     forAll(faces, faceI)
 985     {
 986         if( changedFacePtr && !(*changedFacePtr)[faceI] )
 987             continue;
 988
 989         // Create the owner pyramid - it will have negative volume
 990         const scalar pyrVol = pyramidPointFaceRef
 991         (
 992             faces[faceI],
 993             ctrs[owner[faceI]]
 994         ).mag(points);
 995
 996         bool badFace(false);
 997
 998         if( pyrVol > -minPyrVol )
 999         {
1000             if( report )
1001             {
1002                 # ifdef USE_OMP
1003                 # pragma omp critical
1004                 # endif
1005                 Pout<< "bool checkFacePyramids("
1006                     << "const bool, const scalar, labelHashSet*) : "
1007                     << "face " << faceI << " points the wrong way. " << endl
1008                     << "Pyramid volume: " << -pyrVol
1009                     << " Face " << faces[faceI] << " area: "
1010                     << faces[faceI].mag(points)
1011                     << " Owner cell: " << owner[faceI] << endl
1012                     << "Owner cell vertex labels: "
1013                     << mesh.cells()[owner[faceI]].labels(faces)
1014                     << endl;
1015             }
1016
1017             badFace = true;
1018         }
1019
1020         if( neighbour[faceI] != -1 )
1021         {
1022             // Create the neighbour pyramid - it will have positive volume
1023             const scalar pyrVol =
1024                 pyramidPointFaceRef
1025                 (
1026                     faces[faceI],
1027                     ctrs[neighbour[faceI]]
1028                 ).mag(points);
1029
1030             if( pyrVol < minPyrVol )
1031             {
1032                 if( report )
1033                 {
1034                     # ifdef USE_OMP
1035                     # pragma omp critical
1036                     # endif
1037                     Pout<< "bool checkFacePyramids("
1038                         << "const bool, const scalar, labelHashSet*) : "
1039                         << "face " << faceI << " points the wrong way. " << endl
1040                         << "Pyramid volume: " << -pyrVol
1041                         << " Face " << faces[faceI] << " area: "
1042                         << faces[faceI].mag(points)
1043                         << " Neighbour cell: " << neighbour[faceI] << endl
1044                         << "Neighbour cell vertex labels: "
1045                         << mesh.cells()[neighbour[faceI]].labels(faces)
1046                         << endl;
1047                 }
1048
1049                 badFace = true;
1050             }
1051         }
1052
1053         if( badFace )
1054         {
1055             if( setPtr )
1056             {
1057                 # ifdef USE_OMP
1058                 # pragma omp critical
1059                 # endif
1060                 setPtr->insert(faceI);
1061             }
1062
1063             ++nErrorPyrs;
1064         }
1065     }
1066
1067     reduce(nErrorPyrs, sumOp<label>());
1068
1069     if( setPtr )
1070     {
1071         //- make sure that processor faces are marked on both sides
1072         const PtrList<processorBoundaryPatch>& procBoundaries =
1073             mesh.procBoundaries();
1074
1075         //- send and receive data where needed
1076         forAll(procBoundaries, patchI)
1077         {
1078             const label start = procBoundaries[patchI].patchStart();
1079             const label size = procBoundaries[patchI].patchSize();
1080
1081             labelLongList markedFaces;
1082             for(label faceI=0;faceI<size;++faceI)
1083             {
1084                 if( setPtr->found(start+faceI) )
1085                     markedFaces.append(faceI);
1086             }
1087
1088             OPstream toOtherProc
1089             (
1090                 Pstream::blocking,
1091                 procBoundaries[patchI].neiProcNo(),
1092                 markedFaces.byteSize()
1093             );
1094
1095             toOtherProc << markedFaces;
1096         }
1097
1098         forAll(procBoundaries, patchI)
1099         {
1100             labelList receivedData;
1101             IPstream fromOtheProc
1102             (
1103                 Pstream::blocking,
1104                 procBoundaries[patchI].neiProcNo()
1105             );
1106             fromOtheProc >> receivedData;
1107
1108             const label start = procBoundaries[patchI].patchStart();
1109             forAll(receivedData, i)
1110                 setPtr->insert(start+receivedData[i]);
1111         }
1112     }
1113
1114     if( nErrorPyrs > 0 )
1115     {
1116         if( Pstream::master() )
1117             WarningIn
1118             (
1119                 "bool checkFacePyramids("
1120                 "const polyMeshGen&, const bool, const scalar,"
1121                 " labelHashSet*, const boolList*)"
1122             )   << "Error in face pyramids: " << nErrorPyrs
1123                 << " faces pointing the wrong way!"
1124                 << endl;
1125
1126         return true;
1127     }
1128     else
1129     {
1130         if( report )
1131             Info<< "Face pyramids OK.\n" << endl;
1132
1133         return false;
1134     }
1135 }
1136
1137 void checkFaceSkewness
1138 (
1139     const polyMeshGen& mesh,
1140     scalarField& faceSkewness,
1141     const boolList* changedFacePtr
1142 )
1143 {
1144     //- Warn if the skew correction vector is more than skewWarning times
1145     //- larger than the face area vector
1146     const labelList& own = mesh.owner();
1147     const labelList& nei = mesh.neighbour();
1148     const label nInternalFaces = mesh.nInternalFaces();
1149
1150     const vectorField& centres = mesh.addressingData().cellCentres();
1151     const vectorField& fCentres = mesh.addressingData().faceCentres();
1152
1153     faceSkewness.setSize(own.size());
1154     faceSkewness = 0.0;
1155
1156     //- check internal faces
1157     # ifdef USE_OMP
1158     # pragma omp parallel for schedule(dynamic, 100)
1159     # endif
1160     for(label faceI=0;faceI<nInternalFaces;++faceI)
1161     {
1162         if( changedFacePtr && !changedFacePtr->operator[](faceI) )
1163             continue;
1164
1165         const scalar dOwn = mag(fCentres[faceI] - centres[own[faceI]]);
1166         const scalar dNei = mag(fCentres[faceI] - centres[nei[faceI]]);
1167
1168         const point faceIntersection =
1169             centres[own[faceI]]*dNei/(dOwn+dNei)
1170           + centres[nei[faceI]]*dOwn/(dOwn+dNei);
1171
1172         faceSkewness[faceI] =
1173             mag(fCentres[faceI] - faceIntersection)
1174             /(mag(centres[nei[faceI]] - centres[own[faceI]]) + VSMALL);
1175     }
1176
1177     if( Pstream::parRun() )
1178     {
1179         //- check parallel boundaries
1180         const PtrList<processorBoundaryPatch>& procBoundaries =
1181             mesh.procBoundaries();
1182
1183         forAll(procBoundaries, patchI)
1184         {
1185             const label start = procBoundaries[patchI].patchStart();
1186
1187             vectorField cCentres(procBoundaries[patchI].patchSize());
1188             forAll(cCentres, faceI)
1189                 cCentres[faceI] = centres[own[start+faceI]];
1190
1191             OPstream toOtherProc
1192             (
1193                 Pstream::blocking,
1194                 procBoundaries[patchI].neiProcNo(),
1195                 cCentres.byteSize()
1196             );
1197
1198             toOtherProc << cCentres;
1199         }
1200
1201         forAll(procBoundaries, patchI)
1202         {
1203             const label start = procBoundaries[patchI].patchStart();
1204
1205             vectorField otherCentres;
1206             IPstream fromOtheProc
1207             (
1208                 Pstream::blocking,
1209                 procBoundaries[patchI].neiProcNo()
1210             );
1211
1212             fromOtheProc >> otherCentres;
1213
1214             //- calculate skewness at processor faces
1215             # ifdef USE_OMP
1216             # pragma omp parallel
1217             # endif
1218             {
1219                 # ifdef USE_OMP
1220                 # pragma omp for schedule(dynamic, 100)
1221                 # endif
1222                 forAll(otherCentres, fI)
1223                 {
1224                     const label faceI = start + fI;
1225
1226                     if(
1227                         changedFacePtr &&
1228                         !changedFacePtr->operator[](faceI)
1229                     )
1230                         continue;
1231
1232                     const point& cOwn = centres[own[faceI]];
1233                     const point& cNei = otherCentres[fI];
1234                     const scalar dOwn = mag(fCentres[faceI] - cOwn);
1235                     const scalar dNei = mag(fCentres[faceI] - cNei);
1236
1237                     const point faceIntersection =
1238                         cOwn*dNei/(dOwn+dNei)
1239                       + cNei*dOwn/(dOwn+dNei);
1240
1241                     faceSkewness[faceI] =
1242                         mag(fCentres[faceI] - faceIntersection)
1243                         /(mag(cOwn - cNei) + VSMALL);
1244                 }
1245             }
1246         }
1247     }
1248
1249     //- boundary faces
1250     const faceListPMG& faces = mesh.faces();
1251     const pointFieldPMG& points = mesh.points();
1252     const PtrList<boundaryPatch>& boundaries = mesh.boundaries();
1253     forAll(boundaries, patchI)
1254     {
1255         label faceI = boundaries[patchI].patchStart();
1256         const label end = faceI + boundaries[patchI].patchSize();
1257
1258         for(;faceI<end;++faceI)
1259         {
1260             const vector d = fCentres[faceI] - centres[own[faceI]];
1261
1262             vector n = faces[faceI].normal(points);
1263             const scalar magn = mag(n);
1264             if( magn > VSMALL )
1265             {
1266                 n /= magn;
1267             }
1268             else
1269             {
1270                 continue;
1271             }
1272
1273             const vector dn = (n & d) * n;
1274
1275             faceSkewness[faceI] = mag(d - dn) / (mag(d) + VSMALL);
1276         }
1277     }
1278 }
1279
1280 bool checkFaceSkewness
1281 (
1282     const polyMeshGen& mesh,
1283     const bool report,
1284     const scalar warnSkew,
1285     labelHashSet* setPtr,
1286     const boolList* changedFacePtr
1287 )
1288 {
1289     scalarField faceSkewness;
1290     checkFaceSkewness(mesh, faceSkewness, changedFacePtr);
1291
1292     //- Warn if the skew correction vector is more than skewWarning times
1293     //- larger than the face area vector
1294     scalar maxSkew = 0.0;
1295     scalar sumSkew = 0.0;
1296
1297     label nWarnSkew = 0;
1298
1299     //- check faces
1300     # ifdef USE_OMP
1301     # pragma omp parallel \
1302     reduction(+ : sumSkew, nWarnSkew)
1303     # endif
1304     {
1305         scalar localMaxSkew(0.0);
1306
1307         # ifdef USE_OMP
1308         # pragma omp for schedule(dynamic, 100)
1309         # endif
1310         forAll(faceSkewness, faceI)
1311         {
1312             if( changedFacePtr && !changedFacePtr->operator[](faceI) )
1313                 continue;
1314
1315             const scalar skewness = faceSkewness[faceI];
1316
1317             // Check if the skewness vector is greater than the PN vector.
1318             if( skewness > warnSkew )
1319             {
1320                 if( report )
1321                 {
1322                     # ifdef USE_OMP
1323                     # pragma omp critical
1324                     # endif
1325                     Pout<< " Severe skewness for face " << faceI
1326                         << " skewness = " << skewness << endl;
1327                 }
1328
1329                 if( setPtr )
1330                 {
1331                     # ifdef USE_OMP
1332                     # pragma omp critical
1333                     # endif
1334                     setPtr->insert(faceI);
1335                 }
1336
1337                 ++nWarnSkew;
1338             }
1339
1340             localMaxSkew = Foam::max(localMaxSkew, skewness);
1341             sumSkew += skewness;
1342         }
1343
1344         # ifdef USE_OMP
1345         # pragma omp critical
1346         # endif
1347         maxSkew = Foam::max(maxSkew, localMaxSkew);
1348     }
1349
1350     reduce(maxSkew, maxOp<scalar>());
1351     reduce(sumSkew, sumOp<scalar>());
1352     reduce(nWarnSkew, sumOp<label>());
1353
1354     if( nWarnSkew > 0 )
1355     {
1356         WarningIn
1357         (
1358             "checkFaceSkewness("
1359             "const polyMeshGen&, const bool, const scalar,"
1360             "labelHashSet*, const boolList*)"
1361         )   << "Large face skewness detected.  Max skewness = " << maxSkew
1362             << " Average skewness = " << sumSkew/faceSkewness.size()
1363             << ".\nThis may impair the quality of the result." << nl
1364             << nWarnSkew << " highly skew faces detected."
1365             << endl;
1366
1367         return true;
1368     }
1369     else
1370     {
1371         if( report )
1372             Info<< "Max skewness = " << maxSkew
1373                 << " Average skewness = " << sumSkew/faceSkewness.size()
1374                 << ".  Face skewness OK.\n" << endl;
1375
1376         return false;
1377     }
1378 }
1379
1380 void checkFaceUniformity
1381 (
1382     const polyMeshGen& mesh,
1383     scalarField& faceUniformity,
1384     const boolList* changedFacePtr
1385 )
1386 {
1387     //- for all internal faces check the uniformity of the mesh at faces
1388     const vectorField& centres = mesh.addressingData().cellCentres();
1389     const vectorField& fCentres = mesh.addressingData().faceCentres();
1390
1391     const labelList& owner = mesh.owner();
1392     const labelList& neighbour = mesh.neighbour();
1393     const label nInternalFaces = mesh.nInternalFaces();
1394
1395     faceUniformity.setSize(owner.size());
1396     faceUniformity = 0.5;
1397
1398     # ifdef USE_OMP
1399     # pragma omp parallel for schedule(dynamic, 100)
1400     # endif
1401     for(label faceI=0;faceI<nInternalFaces;++faceI)
1402     {
1403         if( changedFacePtr && !changedFacePtr->operator[](faceI) )
1404             continue;
1405
1406         const scalar dOwn
1407         (
1408             Foam::mag(centres[owner[faceI]] - fCentres[faceI])
1409         );
1410         const scalar dNei
1411         (
1412             Foam::mag(centres[neighbour[faceI]] - fCentres[faceI])
1413         );
1414
1415         faceUniformity[faceI] = Foam::min(dOwn, dNei) / (dOwn + dNei);
1416     }
1417
1418     if( Pstream::parRun() )
1419     {
1420         const PtrList<processorBoundaryPatch>& procBoundaries =
1421             mesh.procBoundaries();
1422
1423         forAll(procBoundaries, patchI)
1424         {
1425             scalarField dst(procBoundaries[patchI].patchSize());
1426             const label start = procBoundaries[patchI].patchStart();
1427
1428             forAll(dst, faceI)
1429             {
1430                 const label fI = start + faceI;
1431                 dst[faceI] = Foam::mag(centres[owner[fI]] - fCentres[fI]);
1432             }
1433
1434             OPstream toOtherProc
1435             (
1436                 Pstream::blocking,
1437                 procBoundaries[patchI].neiProcNo(),
1438                 dst.byteSize()
1439             );
1440
1441             toOtherProc << dst;
1442         }
1443
1444         forAll(procBoundaries, patchI)
1445         {
1446             const label start = procBoundaries[patchI].patchStart();
1447
1448             scalarField otherDst;
1449             IPstream fromOtheProc
1450             (
1451                 Pstream::blocking,
1452                 procBoundaries[patchI].neiProcNo()
1453             );
1454
1455             fromOtheProc >> otherDst;
1456
1457             forAll(otherDst, faceI)
1458             {
1459                 const label fI = start + faceI;
1460                 const scalar dOwn =
1461                     Foam::mag(centres[owner[fI]] - fCentres[fI]);
1462                 const scalar dNei = otherDst[faceI];
1463                 faceUniformity[fI] = Foam::min(dOwn, dNei) / (dOwn + dNei);
1464             }
1465         }
1466     }
1467 }
1468
1469 bool checkFaceUniformity
1470 (
1471     const polyMeshGen& mesh,
1472     const bool report,
1473     const scalar warnUniform,
1474     labelHashSet* setPtr,
1475     const boolList* changedFacePtr
1476 )
1477 {
1478     scalarField faceUniformity;
1479     checkFaceUniformity(mesh, faceUniformity, changedFacePtr);
1480
1481     //- for all internal faces check the uniformity of the mesh at faces
1482     const label nInternalFaces = mesh.nInternalFaces();
1483
1484     scalar maxUniformity = 0.0;
1485     scalar minUniformity = VGREAT;
1486
1487     scalar sumUniformity = 0.0;
1488
1489     label severeNonUniform = 0;
1490
1491     # ifdef USE_OMP
1492     # pragma omp parallel \
1493     reduction(+ : sumUniformity, severeNonUniform)
1494     # endif
1495     {
1496         scalar localMinUniformity(VGREAT);
1497         scalar localMaxUniformity(0.0);
1498
1499         # ifdef USE_OMP
1500         # pragma omp for schedule(guided)
1501         # endif
1502         for(label faceI=0;faceI<nInternalFaces;++faceI)
1503         {
1504             if( changedFacePtr && !changedFacePtr->operator[](faceI) )
1505                 continue;
1506
1507             const scalar uniformity = faceUniformity[faceI];
1508
1509             if( uniformity < warnUniform )
1510             {
1511                 if( setPtr )
1512                 {
1513                     # ifdef USE_OMP
1514                     # pragma omp critical
1515                     # endif
1516                     setPtr->insert(faceI);
1517                 }
1518
1519                 ++severeNonUniform;
1520             }
1521
1522             localMaxUniformity = Foam::max(localMaxUniformity, uniformity);
1523             localMinUniformity = Foam::min(localMinUniformity, uniformity);
1524             sumUniformity += uniformity;
1525         }
1526
1527         # ifdef USE_OMP
1528         # pragma omp critical
1529         # endif
1530         {
1531             maxUniformity = Foam::max(maxUniformity, localMaxUniformity);
1532             minUniformity = Foam::min(minUniformity, localMinUniformity);
1533         }
1534     }
1535
1536     label counter = nInternalFaces;
1537
1538     if( Pstream::parRun() )
1539     {
1540         const PtrList<processorBoundaryPatch>& procBoundaries =
1541             mesh.procBoundaries();
1542
1543         forAll(procBoundaries, patchI)
1544         {
1545             const label start = procBoundaries[patchI].patchStart();
1546             const label size = procBoundaries[patchI].patchSize();
1547
1548             for(label fI=0;fI<size;++fI)
1549             {
1550                 const label faceI = start + fI;
1551
1552                 const scalar uniformity = faceUniformity[faceI];
1553
1554                 if( uniformity < warnUniform )
1555                 {
1556                     if( setPtr)
1557                         setPtr->insert(faceI);
1558
1559                     ++severeNonUniform;
1560                 }
1561
1562                 maxUniformity = Foam::max(maxUniformity, uniformity);
1563                 minUniformity = Foam::min(minUniformity, uniformity);
1564                 sumUniformity += 0.5 * uniformity;
1565             }
1566
1567             if( procBoundaries[patchI].owner() )
1568                 counter += size;
1569         }
1570     }
1571
1572     reduce(maxUniformity, maxOp<scalar>());
1573     reduce(minUniformity, minOp<scalar>());
1574     reduce(sumUniformity, sumOp<scalar>());
1575     reduce(severeNonUniform, sumOp<label>());
1576     reduce(counter, sumOp<label>());
1577
1578     // Only report if there are some internal faces
1579     if( counter > 0 )
1580     {
1581         if( minUniformity < warnUniform )
1582             Info<< "Number of severely non-uniform faces: "
1583                 << severeNonUniform << "." << endl;
1584     }
1585
1586     if( report )
1587     {
1588         if( counter > 0 )
1589             Info<< "Mesh non-uniformity Max: "
1590                 << maxUniformity
1591                 << " Min: " << minUniformity
1592                 << " average: " << sumUniformity/counter
1593                 << endl;
1594     }
1595
1596     return false;
1597 }
1598
1599 void checkVolumeUniformity
1600 (
1601     const polyMeshGen& mesh,
1602     scalarField& faceUniformity,
1603     const boolList* changedFacePtr
1604 );
1605
1606 bool checkVolumeUniformity
1607 (
1608     const polyMeshGen&,
1609     const bool report,
1610     const scalar warnUniform,
1611     labelHashSet* setPtr,
1612     const boolList* changedFacePtr
1613 )
1614 {
1615
1616     return false;
1617 }
1618
1619 bool checkFaceAngles
1620 (
1621     const polyMeshGen& mesh,
1622     const bool report,
1623     const scalar maxDeg,
1624     labelHashSet* setPtr,
1625     const boolList* changedFacePtr
1626 )
1627 {
1628     if( maxDeg < -SMALL || maxDeg > 180+SMALL )
1629     {
1630         FatalErrorIn
1631         (
1632             "bool checkFaceAngles("
1633             "const polyMeshGen&, const bool, const scalar,"
1634             " labelHashSet*, const boolList*)"
1635         )   << "maxDeg should be [0..180] but is now " << maxDeg
1636             << abort(FatalError);
1637     }
1638
1639     const scalar maxSin = Foam::sin(maxDeg/180.0*M_PI);
1640
1641     const pointFieldPMG& points = mesh.points();
1642     const faceListPMG& faces = mesh.faces();
1643     vectorField faceNormals(mesh.addressingData().faceAreas());
1644     faceNormals /= mag(faceNormals) + VSMALL;
1645
1646     scalar maxEdgeSin = 0.0;
1647
1648     label nConcave = 0;
1649
1650     # ifdef USE_OMP
1651     # pragma omp parallel reduction(+ : nConcave)
1652     # endif
1653     {
1654         scalar localMaxEdgeSin(0.0);
1655         label errorFaceI(-1);
1656
1657         # ifdef USE_OMP
1658         # pragma omp for schedule(guided)
1659         # endif
1660         forAll(faces, faceI)
1661         {
1662             if( changedFacePtr && !changedFacePtr->operator[](faceI) )
1663                 continue;
1664
1665             const face& f = faces[faceI];
1666
1667             // Get edge from f[0] to f[size-1];
1668             vector ePrev(points[f[0]] - points[f[f.size()-1]]);
1669             scalar magEPrev = mag(ePrev);
1670             ePrev /= magEPrev + VSMALL;
1671
1672             forAll(f, fp0)
1673             {
1674                 // Get vertex after fp
1675                 label fp1 = f.fcIndex(fp0);
1676
1677                 // Normalized vector between two consecutive points
1678                 vector e10(points[f[fp1]] - points[f[fp0]]);
1679                 scalar magE10 = mag(e10);
1680                 e10 /= magE10 + VSMALL;
1681
1682                 if( magEPrev > SMALL && magE10 > SMALL )
1683                 {
1684                     vector edgeNormal = ePrev ^ e10;
1685                     scalar magEdgeNormal = mag(edgeNormal);
1686
1687                     if( magEdgeNormal < maxSin )
1688                     {
1689                         // Edges (almost) aligned -> face is ok.
1690                     }
1691                     else
1692                     {
1693                         // Check normal
1694                         edgeNormal /= magEdgeNormal;
1695
1696                         if( (edgeNormal & faceNormals[faceI]) < SMALL )
1697                         {
1698                             # ifdef USE_OMP
1699                             # pragma omp critical
1700                             # endif
1701                             {
1702                                 if( faceI != errorFaceI )
1703                                 {
1704                                     // Count only one error per face.
1705                                     errorFaceI = faceI;
1706                                     ++nConcave;
1707                                 }
1708
1709                                 if( setPtr )
1710                                     setPtr->insert(faceI);
1711
1712                                 localMaxEdgeSin =
1713                                     Foam::max(localMaxEdgeSin, magEdgeNormal);
1714                             }
1715                         }
1716                     }
1717                 }
1718
1719                 ePrev = e10;
1720                 magEPrev = magE10;
1721             }
1722         }
1723
1724         # ifdef USE_OMP
1725         # pragma omp critical
1726         # endif
1727         maxEdgeSin = Foam::max(maxEdgeSin, localMaxEdgeSin);
1728     }
1729
1730     reduce(nConcave, sumOp<label>());
1731     reduce(maxEdgeSin, maxOp<scalar>());
1732
1733     if( report )
1734     {
1735         if( maxEdgeSin > SMALL )
1736         {
1737             scalar maxConcaveDegr =
1738                 Foam::asin(Foam::min(1.0, maxEdgeSin))
1739              * 180.0/M_PI;
1740
1741             Warning<< "There are " << nConcave
1742                 << " faces with concave angles between consecutive"
1743                 << " edges. Max concave angle = "
1744                 << maxConcaveDegr
1745                 << " degrees.\n" << endl;
1746         }
1747         else
1748         {
1749             Info<< "All angles in faces are convex or less than "  << maxDeg
1750                 << " degrees concave.\n" << endl;
1751         }
1752     }
1753
1754     if( nConcave > 0 )
1755     {
1756         WarningIn
1757         (
1758             "bool checkFaceAngles("
1759             "const polyMeshGen&, const bool, const scalar,"
1760             " labelHashSet*, const boolList*)"
1761         )   << nConcave  << " face points with severe concave angle (> "
1762             << maxDeg << " deg) found.\n"
1763             << endl;
1764
1765         return true;
1766     }
1767     else
1768     {
1769         return false;
1770     }
1771 }
1772
1773 // Check warpage of faces. Is calculated as the difference between areas of
1774 // individual triangles and the overall area of the face (which ifself is
1775 // is the average of the areas of the individual triangles).
1776 bool checkFaceFlatness
1777 (
1778     const polyMeshGen& mesh,
1779     const bool report,
1780     const scalar warnFlatness,
1781     labelHashSet* setPtr,
1782     const boolList* changedFacePtr
1783 )
1784 {
1785     if( warnFlatness < 0 || warnFlatness > 1 )
1786     {
1787         FatalErrorIn
1788         (
1789             "bool checkFaceFlatness("
1790             "const polyMeshGen&, const bool, const scalar,"
1791             " labelHashSet*, const boolList*)"
1792         )   << "warnFlatness should be [0..1] but is now " << warnFlatness
1793             << abort(FatalError);
1794     }
1795
1796     const pointFieldPMG& points = mesh.points();
1797     const faceListPMG& faces = mesh.faces();
1798     const vectorField& fctrs = mesh.addressingData().faceCentres();
1799
1800     // Areas are calculated as the sum of areas. (see
1801     // polyMeshGenAddressingFaceCentresAndAreas.C)
1802     scalarField magAreas(mag(mesh.addressingData().faceAreas()));
1803
1804     label nWarped = 0;
1805
1806     scalar minFlatness = GREAT;
1807     scalar sumFlatness = 0;
1808     label nSummed = 0;
1809
1810     # ifdef USE_OMP
1811     # pragma omp parallel if( faces.size() > 1000 ) \
1812     reduction(+ : nSummed, nWarped) reduction(+ : sumFlatness)
1813     # endif
1814     {
1815         scalar minFlatnessProc = VGREAT;
1816
1817         # ifdef USE_OMP
1818         # pragma omp for schedule(guided)
1819         # endif
1820         forAll(faces, faceI)
1821         {
1822             if( changedFacePtr && !(*changedFacePtr)[faceI] )
1823                 continue;
1824
1825             const face& f = faces[faceI];
1826
1827             if( f.size() > 3 && magAreas[faceI] > VSMALL )
1828             {
1829                 const point& fc = fctrs[faceI];
1830
1831                 //- Calculate the sum of magnitude of areas and compare
1832                 //- the magnitude of sum of areas.
1833
1834                 scalar sumA = 0.0;
1835
1836                 forAll(f, fp)
1837                 {
1838                     const point& thisPoint = points[f[fp]];
1839                     const point& nextPoint = points[f.nextLabel(fp)];
1840
1841                     // Triangle around fc.
1842                     vector n = 0.5*((nextPoint - thisPoint)^(fc - thisPoint));
1843                     sumA += mag(n);
1844                 }
1845
1846                 scalar flatness = magAreas[faceI] / (sumA+VSMALL);
1847
1848                 sumFlatness += flatness;
1849                 ++nSummed;
1850
1851                 minFlatnessProc = Foam::min(minFlatnessProc, flatness);
1852
1853                 if( flatness < warnFlatness )
1854                 {
1855                     ++nWarped;
1856
1857                     if( setPtr )
1858                     {
1859                         # ifdef USE_OMP
1860                         # pragma omp critical
1861                         # endif
1862                         setPtr->insert(faceI);
1863                     }
1864                 }
1865             }
1866         }
1867
1868         # ifdef USE_OMP
1869         # pragma omp critical
1870         # endif
1871         {
1872             minFlatness = Foam::min(minFlatness, minFlatnessProc);
1873         }
1874     }
1875
1876     if( Pstream::parRun() && setPtr )
1877     {
1878         //- make sure that processor faces are marked on both sides
1879         const PtrList<processorBoundaryPatch>& procBoundaries =
1880             mesh.procBoundaries();
1881
1882         List<DynList<label> > markedFaces(procBoundaries.size());
1883         forAll(procBoundaries, patchI)
1884         {
1885             const label start = procBoundaries[patchI].patchStart();
1886             const label size = procBoundaries[patchI].patchSize();
1887
1888             for(label i=0;i<size;++i)
1889                 if( setPtr->found(start+i) )
1890                     markedFaces[patchI].append(i);
1891         }
1892
1893         //- exchange list sizes
1894         forAll(procBoundaries, patchI)
1895         {
1896             OPstream toOtherProc
1897             (
1898                 Pstream::blocking,
1899                 procBoundaries[patchI].neiProcNo(),
1900                 sizeof(label)
1901             );
1902
1903             toOtherProc << markedFaces[patchI].size();
1904         }
1905
1906         labelList nMarkedOnOtherProcs(procBoundaries.size());
1907         forAll(procBoundaries, patchI)
1908         {
1909             IPstream fromOtheProc
1910             (
1911                 Pstream::blocking,
1912                 procBoundaries[patchI].neiProcNo(),
1913                 sizeof(label)
1914             );
1915
1916             fromOtheProc >> nMarkedOnOtherProcs[patchI];
1917         }
1918
1919         //- exchange data
1920         forAll(procBoundaries, patchI)
1921         {
1922             if( markedFaces[patchI].size() == 0 )
1923                 continue;
1924
1925             OPstream toOtherProc
1926             (
1927                 Pstream::blocking,
1928                 procBoundaries[patchI].neiProcNo(),
1929                 markedFaces[patchI].byteSize()
1930             );
1931
1932             toOtherProc << markedFaces[patchI];
1933         }
1934
1935         forAll(procBoundaries, patchI)
1936         {
1937             if( nMarkedOnOtherProcs[patchI] == 0 )
1938                 continue;
1939
1940             labelList receivedData;
1941             IPstream fromOtheProc
1942             (
1943                 Pstream::blocking,
1944                 procBoundaries[patchI].neiProcNo(),
1945                 nMarkedOnOtherProcs[patchI]*sizeof(label)
1946             );
1947             fromOtheProc >> receivedData;
1948
1949             const label start = procBoundaries[patchI].patchStart();
1950             forAll(receivedData, i)
1951                 if( !setPtr->found(start+receivedData[i]) )
1952                     setPtr->insert(start+receivedData[i]);
1953         }
1954     }
1955
1956     reduce(nWarped, sumOp<label>());
1957     reduce(minFlatness, minOp<scalar>());
1958
1959     reduce(nSummed, sumOp<label>());
1960     reduce(sumFlatness, sumOp<scalar>());
1961
1962     if( report )
1963     {
1964         if( nSummed > 0 )
1965         {
1966             Info<< "Face flatness (1 = flat, 0 = butterfly) : average = "
1967                 << sumFlatness / nSummed << "  min = " << minFlatness << endl;
1968         }
1969
1970         if( nWarped> 0 )
1971         {
1972             Info<< "There are " << nWarped
1973                 << " faces with ratio between projected and actual area < "
1974                 << warnFlatness
1975                 << ".\nMinimum ratio (minimum flatness, maximum warpage) = "
1976                 << minFlatness << nl << endl;
1977         }
1978         else
1979         {
1980             Info<< "All faces are flat in that the ratio between projected"
1981                 << " and actual area is > " << warnFlatness << nl << endl;
1982         }
1983     }
1984
1985     if( nWarped > 0 )
1986     {
1987         WarningIn
1988         (
1989             "bool checkFaceFlatness("
1990             "const polyMeshGen&, const bool, const scalar,"
1991             " labelHashSet*, const boolList*)"
1992         )   << nWarped  << " faces with severe warpage (flatness < "
1993             << warnFlatness << ") found.\n"
1994             << endl;
1995
1996         return true;
1997     }
1998     else
1999     {
2000         return false;
2001     }
2002 }
2003
2004 label findBadFacesRelaxed
2005 (
2006     const polyMeshGen& mesh,
2007     labelHashSet& badFaces,
2008     const bool report,
2009     const boolList* activeFacePtr
2010 )
2011 {
2012     badFaces.clear();
2013
2014     polyMeshGenChecks::checkFacePyramids
2015     (
2016         mesh,
2017         report,
2018         VSMALL,
2019         &badFaces,
2020         activeFacePtr
2021     );
2022
2023     polyMeshGenChecks::checkFaceAreas
2024     (
2025         mesh,
2026         report,
2027         VSMALL,
2028         &badFaces,
2029         activeFacePtr
2030     );
2031
2032     const label nBadFaces = returnReduce(badFaces.size(), sumOp<label>());
2033
2034     return nBadFaces;
2035 }
2036
2037 label findBadFaces
2038 (
2039     const polyMeshGen& mesh,
2040     labelHashSet& badFaces,
2041     const bool report,
2042     const boolList* activeFacePtr
2043 )
2044 {
2045     badFaces.clear();
2046
2047     polyMeshGenChecks::checkFacePyramids
2048     (
2049         mesh,
2050         report,
2051         VSMALL,
2052         &badFaces,
2053         activeFacePtr
2054     );
2055
2056     polyMeshGenChecks::checkFaceFlatness
2057     (
2058         mesh,
2059         report,
2060         0.8,
2061         &badFaces,
2062         activeFacePtr
2063     );
2064
2065     polyMeshGenChecks::checkCellPartTetrahedra
2066     (
2067         mesh,
2068         report,
2069         VSMALL,
2070         &badFaces,
2071         activeFacePtr
2072     );
2073
2074     polyMeshGenChecks::checkFaceAreas
2075     (
2076         mesh,
2077         report,
2078         VSMALL,
2079         &badFaces,
2080         activeFacePtr
2081     );
2082
2083     const label nBadFaces = returnReduce(badFaces.size(), sumOp<label>());
2084
2085     return nBadFaces;
2086 }
2087
2088 label findLowQualityFaces
2089 (
2090     const polyMeshGen& mesh,
2091     labelHashSet& badFaces,
2092     const bool report,
2093     const boolList* activeFacePtr
2094 )
2095 {
2096     badFaces.clear();
2097
2098     polyMeshGenChecks::checkFaceDotProduct
2099     (
2100         mesh,
2101         report,
2102         65.0,
2103         &badFaces,
2104         activeFacePtr
2105     );
2106
2107     polyMeshGenChecks::checkFaceSkewness
2108     (
2109         mesh,
2110         report,
2111         2.0,
2112         &badFaces,
2113         activeFacePtr
2114     );
2115
2116     const label nBadFaces = returnReduce(badFaces.size(), sumOp<label>());
2117
2118     return nBadFaces;
2119 }
2120
2121 label findWorstQualityFaces
2122 (
2123     const polyMeshGen& mesh,
2124     labelHashSet& badFaces,
2125     const bool report,
2126     const boolList* activeFacePtr,
2127     const scalar relativeThreshold
2128 )
2129 {
2130     badFaces.clear();
2131
2132     scalarField checkValues;
2133     polyMeshGenChecks::checkFaceDotProduct
2134     (
2135         mesh,
2136         checkValues,
2137         activeFacePtr
2138     );
2139
2140     scalar minNonOrtho = returnReduce(min(checkValues), minOp<scalar>());
2141     const scalar warnNonOrtho =
2142         minNonOrtho + relativeThreshold * (1.0 - minNonOrtho);
2143
2144     Info << "Worst non-orthogonality " << Foam::acos(minNonOrtho) * 180.0 / M_PI
2145          << " selecting faces with non-orthogonality greater than "
2146          << (Foam::acos(warnNonOrtho) * 180.0 / M_PI) << endl;
2147
2148     forAll(checkValues, faceI)
2149     {
2150         if
2151         (
2152             activeFacePtr && activeFacePtr->operator[](faceI) &&
2153             checkValues[faceI] < warnNonOrtho
2154         )
2155             badFaces.insert(faceI);
2156     }
2157
2158     polyMeshGenChecks::checkFaceSkewness
2159     (
2160         mesh,
2161         checkValues,
2162         activeFacePtr
2163     );
2164
2165     const scalar maxSkew = returnReduce(max(checkValues), maxOp<scalar>());
2166     const scalar warnSkew = (1.0 - relativeThreshold) * maxSkew;
2167     forAll(checkValues, faceI)
2168         if
2169         (
2170             activeFacePtr && activeFacePtr->operator[](faceI) &&
2171             checkValues[faceI] > warnSkew
2172         )
2173             badFaces.insert(faceI);
2174
2175     Info << "Maximum skewness in the mesh is " << maxSkew
2176          << " selecting faces with skewness greater than " << warnSkew << endl;
2177
2178     const label nBadFaces = returnReduce(badFaces.size(), sumOp<label>());
2179
2180     Info << "Selected " << nBadFaces
2181          << " out of " << returnReduce(checkValues.size(), sumOp<label>())
2182          << " faces" << endl;
2183
2184     return nBadFaces;
2185 }
2186
2187 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
2188
2189 } // End namespace polyMeshGenChecks
2190
2191 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
2192
2193 } // End namespace Foam
2194
2195 // ************************************************************************* //