src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C

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