src/meshTools/primitiveMeshGeometry/primitiveMeshGeometry.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your 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, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "primitiveMeshGeometry.H"
  27 #include "pyramidPointFaceRef.H"
  28 #include "unitConversion.H"
  29
  30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  31
  32 namespace Foam
  33 {
  34     defineTypeNameAndDebug(primitiveMeshGeometry, 0);
  35 }
  36
  37
  38 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  39
  40 void Foam::primitiveMeshGeometry::updateFaceCentresAndAreas
  41 (
  42     const pointField& p,
  43     const labelList& changedFaces
  44 )
  45 {
  46     const faceList& fs = mesh_.faces();
  47
  48     forAll(changedFaces, i)
  49     {
  50         label facei = changedFaces[i];
  51
  52         const labelList& f = fs[facei];
  53         label nPoints = f.size();
  54
  55         // If the face is a triangle, do a direct calculation for efficiency
  56         // and to avoid round-off error-related problems
  57         if (nPoints == 3)
  58         {
  59             faceCentres_[facei] = (1.0/3.0)*(p[f[0]] + p[f[1]] + p[f[2]]);
  60             faceAreas_[facei] = 0.5*((p[f[1]] - p[f[0]])^(p[f[2]] - p[f[0]]));
  61         }
  62         else
  63         {
  64             vector sumN = vector::zero;
  65             scalar sumA = 0.0;
  66             vector sumAc = vector::zero;
  67
  68             point fCentre = p[f[0]];
  69             for (label pi = 1; pi < nPoints; pi++)
  70             {
  71                 fCentre += p[f[pi]];
  72             }
  73
  74             fCentre /= nPoints;
  75
  76             for (label pi = 0; pi < nPoints; pi++)
  77             {
  78                 const point& nextPoint = p[f[(pi + 1) % nPoints]];
  79
  80                 vector c = p[f[pi]] + nextPoint + fCentre;
  81                 vector n = (nextPoint - p[f[pi]])^(fCentre - p[f[pi]]);
  82                 scalar a = mag(n);
  83
  84                 sumN += n;
  85                 sumA += a;
  86                 sumAc += a*c;
  87             }
  88
  89             faceCentres_[facei] = (1.0/3.0)*sumAc/(sumA + VSMALL);
  90             faceAreas_[facei] = 0.5*sumN;
  91         }
  92     }
  93 }
  94
  95
  96 void Foam::primitiveMeshGeometry::updateCellCentresAndVols
  97 (
  98     const labelList& changedCells,
  99     const labelList& changedFaces
 100 )
 101 {
 102     // Clear the fields for accumulation
 103     UIndirectList<vector>(cellCentres_, changedCells) = vector::zero;
 104     UIndirectList<scalar>(cellVolumes_, changedCells) = 0.0;
 105
 106     const labelList& own = mesh_.faceOwner();
 107     const labelList& nei = mesh_.faceNeighbour();
 108
 109     // first estimate the approximate cell centre as the average of face centres
 110
 111     vectorField cEst(mesh_.nCells());
 112     UIndirectList<vector>(cEst, changedCells) = vector::zero;
 113     scalarField nCellFaces(mesh_.nCells());
 114     UIndirectList<scalar>(nCellFaces, changedCells) = 0.0;
 115
 116     forAll(changedFaces, i)
 117     {
 118         label faceI = changedFaces[i];
 119         cEst[own[faceI]] += faceCentres_[faceI];
 120         nCellFaces[own[faceI]] += 1;
 121
 122         if (mesh_.isInternalFace(faceI))
 123         {
 124             cEst[nei[faceI]] += faceCentres_[faceI];
 125             nCellFaces[nei[faceI]] += 1;
 126         }
 127     }
 128
 129     forAll(changedCells, i)
 130     {
 131         label cellI = changedCells[i];
 132         cEst[cellI] /= nCellFaces[cellI];
 133     }
 134
 135     forAll(changedFaces, i)
 136     {
 137         label faceI = changedFaces[i];
 138
 139         // Calculate 3*face-pyramid volume
 140         scalar pyr3Vol = max
 141         (
 142             faceAreas_[faceI] & (faceCentres_[faceI] - cEst[own[faceI]]),
 143             VSMALL
 144         );
 145
 146         // Calculate face-pyramid centre
 147         vector pc = (3.0/4.0)*faceCentres_[faceI] + (1.0/4.0)*cEst[own[faceI]];
 148
 149         // Accumulate volume-weighted face-pyramid centre
 150         cellCentres_[own[faceI]] += pyr3Vol*pc;
 151
 152         // Accumulate face-pyramid volume
 153         cellVolumes_[own[faceI]] += pyr3Vol;
 154
 155         if (mesh_.isInternalFace(faceI))
 156         {
 157             // Calculate 3*face-pyramid volume
 158             scalar pyr3Vol = max
 159             (
 160                 faceAreas_[faceI] & (cEst[nei[faceI]] - faceCentres_[faceI]),
 161                 VSMALL
 162             );
 163
 164             // Calculate face-pyramid centre
 165             vector pc =
 166                 (3.0/4.0)*faceCentres_[faceI]
 167               + (1.0/4.0)*cEst[nei[faceI]];
 168
 169             // Accumulate volume-weighted face-pyramid centre
 170             cellCentres_[nei[faceI]] += pyr3Vol*pc;
 171
 172             // Accumulate face-pyramid volume
 173             cellVolumes_[nei[faceI]] += pyr3Vol;
 174         }
 175     }
 176
 177     forAll(changedCells, i)
 178     {
 179         label cellI = changedCells[i];
 180
 181         cellCentres_[cellI] /= cellVolumes_[cellI];
 182         cellVolumes_[cellI] *= (1.0/3.0);
 183     }
 184 }
 185
 186
 187 Foam::labelList Foam::primitiveMeshGeometry::affectedCells
 188 (
 189     const labelList& changedFaces
 190 ) const
 191 {
 192     const labelList& own = mesh_.faceOwner();
 193     const labelList& nei = mesh_.faceNeighbour();
 194
 195     labelHashSet affectedCells(2*changedFaces.size());
 196
 197     forAll(changedFaces, i)
 198     {
 199         label faceI = changedFaces[i];
 200
 201         affectedCells.insert(own[faceI]);
 202
 203         if (mesh_.isInternalFace(faceI))
 204         {
 205             affectedCells.insert(nei[faceI]);
 206         }
 207     }
 208     return affectedCells.toc();
 209 }
 210
 211
 212
 213 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 214
 215 // Construct from components
 216 Foam::primitiveMeshGeometry::primitiveMeshGeometry
 217 (
 218     const primitiveMesh& mesh
 219 )
 220 :
 221     mesh_(mesh)
 222 {
 223     correct();
 224 }
 225
 226
 227 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 228
 229
 230 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 231
 232 //- Take over properties from mesh
 233 void Foam::primitiveMeshGeometry::correct()
 234 {
 235     faceAreas_ = mesh_.faceAreas();
 236     faceCentres_ = mesh_.faceCentres();
 237     cellCentres_ = mesh_.cellCentres();
 238     cellVolumes_ = mesh_.cellVolumes();
 239 }
 240
 241
 242 //- Recalculate on selected faces
 243 void Foam::primitiveMeshGeometry::correct
 244 (
 245     const pointField& p,
 246     const labelList& changedFaces
 247 )
 248 {
 249     // Update face quantities
 250     updateFaceCentresAndAreas(p, changedFaces);
 251     // Update cell quantities from face quantities
 252     updateCellCentresAndVols(affectedCells(changedFaces), changedFaces);
 253 }
 254
 255
 256 bool Foam::primitiveMeshGeometry::checkFaceDotProduct
 257 (
 258     const bool report,
 259     const scalar orthWarn,
 260     const primitiveMesh& mesh,
 261     const vectorField& cellCentres,
 262     const vectorField& faceAreas,
 263     const labelList& checkFaces,
 264     labelHashSet* setPtr
 265 )
 266 {
 267     // for all internal faces check theat the d dot S product is positive
 268
 269     const labelList& own = mesh.faceOwner();
 270     const labelList& nei = mesh.faceNeighbour();
 271
 272     // Severe nonorthogonality threshold
 273     const scalar severeNonorthogonalityThreshold = ::cos(degToRad(orthWarn));
 274
 275     scalar minDDotS = GREAT;
 276
 277     scalar sumDDotS = 0;
 278
 279     label severeNonOrth = 0;
 280
 281     label errorNonOrth = 0;
 282
 283     forAll(checkFaces, i)
 284     {
 285         label faceI = checkFaces[i];
 286
 287         if (mesh.isInternalFace(faceI))
 288         {
 289             vector d = cellCentres[nei[faceI]] - cellCentres[own[faceI]];
 290             const vector& s = faceAreas[faceI];
 291
 292             scalar dDotS = (d & s)/(mag(d)*mag(s) + VSMALL);
 293
 294             if (dDotS < severeNonorthogonalityThreshold)
 295             {
 296                 if (dDotS > SMALL)
 297                 {
 298                     if (report)
 299                     {
 300                         // Severe non-orthogonality but mesh still OK
 301                         Pout<< "Severe non-orthogonality for face " << faceI
 302                             << " between cells " << own[faceI]
 303                             << " and " << nei[faceI]
 304                             << ": Angle = " << radToDeg(::acos(dDotS))
 305                             << " deg." << endl;
 306                     }
 307
 308                     if (setPtr)
 309                     {
 310                         setPtr->insert(faceI);
 311                     }
 312
 313                     severeNonOrth++;
 314                 }
 315                 else
 316                 {
 317                     // Non-orthogonality greater than 90 deg
 318                     if (report)
 319                     {
 320                         WarningIn
 321                         (
 322                             "primitiveMeshGeometry::checkFaceDotProduct"
 323                             "(const bool, const scalar, const labelList&"
 324                             ", labelHashSet*)"
 325                         )   << "Severe non-orthogonality detected for face "
 326                             << faceI
 327                             << " between cells " << own[faceI] << " and "
 328                             << nei[faceI]
 329                             << ": Angle = " << radToDeg(::acos(dDotS))
 330                             << " deg." << endl;
 331                     }
 332
 333                     errorNonOrth++;
 334
 335                     if (setPtr)
 336                     {
 337                         setPtr->insert(faceI);
 338                     }
 339                 }
 340             }
 341
 342             if (dDotS < minDDotS)
 343             {
 344                 minDDotS = dDotS;
 345             }
 346
 347             sumDDotS += dDotS;
 348         }
 349     }
 350
 351     reduce(minDDotS, minOp<scalar>());
 352     reduce(sumDDotS, sumOp<scalar>());
 353     reduce(severeNonOrth, sumOp<label>());
 354     reduce(errorNonOrth, sumOp<label>());
 355
 356     label neiSize = nei.size();
 357     reduce(neiSize, sumOp<label>());
 358
 359     // Only report if there are some internal faces
 360     if (neiSize > 0)
 361     {
 362         if (report && minDDotS < severeNonorthogonalityThreshold)
 363         {
 364             Info<< "Number of non-orthogonality errors: " << errorNonOrth
 365                 << ". Number of severely non-orthogonal faces: "
 366                 << severeNonOrth  << "." << endl;
 367         }
 368     }
 369
 370     if (report)
 371     {
 372         if (neiSize > 0)
 373         {
 374             Info<< "Mesh non-orthogonality Max: "
 375                 << radToDeg(::acos(minDDotS))
 376                 << " average: " << radToDeg(::acos(sumDDotS/neiSize))
 377                 << endl;
 378         }
 379     }
 380
 381     if (errorNonOrth > 0)
 382     {
 383         if (report)
 384         {
 385             SeriousErrorIn
 386             (
 387                 "primitiveMeshGeometry::checkFaceDotProduct"
 388                 "(const bool, const scalar, const labelList&, labelHashSet*)"
 389             )   << "Error in non-orthogonality detected" << endl;
 390         }
 391
 392         return true;
 393     }
 394     else
 395     {
 396         if (report)
 397         {
 398             Info<< "Non-orthogonality check OK.\n" << endl;
 399         }
 400
 401         return false;
 402     }
 403 }
 404
 405
 406 bool Foam::primitiveMeshGeometry::checkFacePyramids
 407 (
 408     const bool report,
 409     const scalar minPyrVol,
 410     const primitiveMesh& mesh,
 411     const vectorField& cellCentres,
 412     const pointField& p,
 413     const labelList& checkFaces,
 414     labelHashSet* setPtr
 415 )
 416 {
 417     // check whether face area vector points to the cell with higher label
 418     const labelList& own = mesh.faceOwner();
 419     const labelList& nei = mesh.faceNeighbour();
 420
 421     const faceList& f = mesh.faces();
 422
 423     label nErrorPyrs = 0;
 424
 425     forAll(checkFaces, i)
 426     {
 427         label faceI = checkFaces[i];
 428
 429         // Create the owner pyramid - it will have negative volume
 430         scalar pyrVol = pyramidPointFaceRef
 431         (
 432             f[faceI],
 433             cellCentres[own[faceI]]
 434         ).mag(p);
 435
 436         if (pyrVol > -minPyrVol)
 437         {
 438             if (report)
 439             {
 440                 Pout<< "bool primitiveMeshGeometry::checkFacePyramids("
 441                     << "const bool, const scalar, const pointField&"
 442                     << ", const labelList&, labelHashSet*): "
 443                     << "face " << faceI << " points the wrong way. " << endl
 444                     << "Pyramid volume: " << -pyrVol
 445                     << " Face " << f[faceI] << " area: " << f[faceI].mag(p)
 446                     << " Owner cell: " << own[faceI] << endl
 447                     << "Owner cell vertex labels: "
 448                     << mesh.cells()[own[faceI]].labels(f)
 449                     << endl;
 450             }
 451
 452
 453             if (setPtr)
 454             {
 455                 setPtr->insert(faceI);
 456             }
 457
 458             nErrorPyrs++;
 459         }
 460
 461         if (mesh.isInternalFace(faceI))
 462         {
 463             // Create the neighbour pyramid - it will have positive volume
 464             scalar pyrVol =
 465                 pyramidPointFaceRef(f[faceI], cellCentres[nei[faceI]]).mag(p);
 466
 467             if (pyrVol < minPyrVol)
 468             {
 469                 if (report)
 470                 {
 471                     Pout<< "bool primitiveMeshGeometry::checkFacePyramids("
 472                         << "const bool, const scalar, const pointField&"
 473                         << ", const labelList&, labelHashSet*): "
 474                         << "face " << faceI << " points the wrong way. " << endl
 475                         << "Pyramid volume: " << -pyrVol
 476                         << " Face " << f[faceI] << " area: " << f[faceI].mag(p)
 477                         << " Neighbour cell: " << nei[faceI] << endl
 478                         << "Neighbour cell vertex labels: "
 479                         << mesh.cells()[nei[faceI]].labels(f)
 480                         << endl;
 481                 }
 482
 483                 if (setPtr)
 484                 {
 485                     setPtr->insert(faceI);
 486                 }
 487
 488                 nErrorPyrs++;
 489             }
 490         }
 491     }
 492
 493     reduce(nErrorPyrs, sumOp<label>());
 494
 495     if (nErrorPyrs > 0)
 496     {
 497         if (report)
 498         {
 499             SeriousErrorIn
 500             (
 501                 "primitiveMeshGeometry::checkFacePyramids("
 502                 "const bool, const scalar, const pointField&"
 503                 ", const labelList&, labelHashSet*)"
 504             )   << "Error in face pyramids: faces pointing the wrong way!"
 505                 << endl;
 506         }
 507
 508         return true;
 509     }
 510     else
 511     {
 512         if (report)
 513         {
 514             Info<< "Face pyramids OK.\n" << endl;
 515         }
 516
 517         return false;
 518     }
 519 }
 520
 521
 522 bool Foam::primitiveMeshGeometry::checkFaceSkewness
 523 (
 524     const bool report,
 525     const scalar internalSkew,
 526     const scalar boundarySkew,
 527     const primitiveMesh& mesh,
 528     const vectorField& cellCentres,
 529     const vectorField& faceCentres,
 530     const vectorField& faceAreas,
 531     const labelList& checkFaces,
 532     labelHashSet* setPtr
 533 )
 534 {
 535     // Warn if the skew correction vector is more than skew times
 536     // larger than the face area vector
 537
 538     const labelList& own = mesh.faceOwner();
 539     const labelList& nei = mesh.faceNeighbour();
 540
 541     scalar maxSkew = 0;
 542
 543     label nWarnSkew = 0;
 544
 545     forAll(checkFaces, i)
 546     {
 547         label faceI = checkFaces[i];
 548
 549         if (mesh.isInternalFace(faceI))
 550         {
 551             scalar dOwn = mag(faceCentres[faceI] - cellCentres[own[faceI]]);
 552             scalar dNei = mag(faceCentres[faceI] - cellCentres[nei[faceI]]);
 553
 554             point faceIntersection =
 555                 cellCentres[own[faceI]]*dNei/(dOwn+dNei)
 556               + cellCentres[nei[faceI]]*dOwn/(dOwn+dNei);
 557
 558             scalar skewness =
 559                 mag(faceCentres[faceI] - faceIntersection)
 560               / (
 561                     mag(cellCentres[nei[faceI]]-cellCentres[own[faceI]])
 562                   + VSMALL
 563                 );
 564
 565             // Check if the skewness vector is greater than the PN vector.
 566             // This does not cause trouble but is a good indication of a poor
 567             // mesh.
 568             if (skewness > internalSkew)
 569             {
 570                 if (report)
 571                 {
 572                     Pout<< "Severe skewness for face " << faceI
 573                         << " skewness = " << skewness << endl;
 574                 }
 575
 576                 if (setPtr)
 577                 {
 578                     setPtr->insert(faceI);
 579                 }
 580
 581                 nWarnSkew++;
 582             }
 583
 584             if (skewness > maxSkew)
 585             {
 586                 maxSkew = skewness;
 587             }
 588         }
 589         else
 590         {
 591             // Boundary faces: consider them to have only skewness error.
 592             // (i.e. treat as if mirror cell on other side)
 593
 594             vector faceNormal = faceAreas[faceI];
 595             faceNormal /= mag(faceNormal) + VSMALL;
 596
 597             vector dOwn = faceCentres[faceI] - cellCentres[own[faceI]];
 598
 599             vector dWall = faceNormal*(faceNormal & dOwn);
 600
 601             point faceIntersection = cellCentres[own[faceI]] + dWall;
 602
 603             scalar skewness =
 604                 mag(faceCentres[faceI] - faceIntersection)
 605                 /(2*mag(dWall) + VSMALL);
 606
 607             // Check if the skewness vector is greater than the PN vector.
 608             // This does not cause trouble but is a good indication of a poor
 609             // mesh.
 610             if (skewness > boundarySkew)
 611             {
 612                 if (report)
 613                 {
 614                     Pout<< "Severe skewness for boundary face " << faceI
 615                         << " skewness = " << skewness << endl;
 616                 }
 617
 618                 if (setPtr)
 619                 {
 620                     setPtr->insert(faceI);
 621                 }
 622
 623                 nWarnSkew++;
 624             }
 625
 626             if (skewness > maxSkew)
 627             {
 628                 maxSkew = skewness;
 629             }
 630         }
 631     }
 632
 633     reduce(maxSkew, maxOp<scalar>());
 634     reduce(nWarnSkew, sumOp<label>());
 635
 636     if (nWarnSkew > 0)
 637     {
 638         if (report)
 639         {
 640             WarningIn
 641             (
 642                 "primitiveMeshGeometry::checkFaceSkewness"
 643                 "(const bool, const scalar, const labelList&, labelHashSet*)"
 644             )   << "Large face skewness detected.  Max skewness = "
 645                 << 100*maxSkew
 646                 << " percent.\nThis may impair the quality of the result." << nl
 647                 << nWarnSkew << " highly skew faces detected."
 648                 << endl;
 649         }
 650
 651         return true;
 652     }
 653     else
 654     {
 655         if (report)
 656         {
 657             Info<< "Max skewness = " << 100*maxSkew
 658                 << " percent.  Face skewness OK.\n" << endl;
 659         }
 660
 661         return false;
 662     }
 663 }
 664
 665
 666 bool Foam::primitiveMeshGeometry::checkFaceWeights
 667 (
 668     const bool report,
 669     const scalar warnWeight,
 670     const primitiveMesh& mesh,
 671     const vectorField& cellCentres,
 672     const vectorField& faceCentres,
 673     const vectorField& faceAreas,
 674     const labelList& checkFaces,
 675     labelHashSet* setPtr
 676 )
 677 {
 678     // Warn if the delta factor (0..1) is too large.
 679
 680     const labelList& own = mesh.faceOwner();
 681     const labelList& nei = mesh.faceNeighbour();
 682
 683     scalar minWeight = GREAT;
 684
 685     label nWarnWeight = 0;
 686
 687     forAll(checkFaces, i)
 688     {
 689         label faceI = checkFaces[i];
 690
 691         if (mesh.isInternalFace(faceI))
 692         {
 693             const point& fc = faceCentres[faceI];
 694
 695             scalar dOwn = mag(faceAreas[faceI] & (fc-cellCentres[own[faceI]]));
 696             scalar dNei = mag(faceAreas[faceI] & (cellCentres[nei[faceI]]-fc));
 697
 698             scalar weight = min(dNei,dOwn)/(dNei+dOwn);
 699
 700             if (weight < warnWeight)
 701             {
 702                 if (report)
 703                 {
 704                     Pout<< "Small weighting factor for face " << faceI
 705                         << " weight = " << weight << endl;
 706                 }
 707
 708                 if (setPtr)
 709                 {
 710                     setPtr->insert(faceI);
 711                 }
 712
 713                 nWarnWeight++;
 714             }
 715
 716             minWeight = min(minWeight, weight);
 717         }
 718     }
 719
 720     reduce(minWeight, minOp<scalar>());
 721     reduce(nWarnWeight, sumOp<label>());
 722
 723     if (minWeight < warnWeight)
 724     {
 725         if (report)
 726         {
 727             WarningIn
 728             (
 729                 "primitiveMeshGeometry::checkFaceWeights"
 730                 "(const bool, const scalar, const labelList&, labelHashSet*)"
 731             )   << "Small interpolation weight detected.  Min weight = "
 732                 << minWeight << '.' << nl
 733                 << nWarnWeight << " faces with small weights detected."
 734                 << endl;
 735         }
 736
 737         return true;
 738     }
 739     else
 740     {
 741         if (report)
 742         {
 743             Info<< "Min weight = " << minWeight
 744                 << " percent.  Weights OK.\n" << endl;
 745         }
 746
 747         return false;
 748     }
 749 }
 750
 751
 752 // Check convexity of angles in a face. Allow a slight non-convexity.
 753 // E.g. maxDeg = 10 allows for angles < 190 (or 10 degrees concavity)
 754 // (if truly concave and points not visible from face centre the face-pyramid
 755 //  check in checkMesh will fail)
 756 bool Foam::primitiveMeshGeometry::checkFaceAngles
 757 (
 758     const bool report,
 759     const scalar maxDeg,
 760     const primitiveMesh& mesh,
 761     const vectorField& faceAreas,
 762     const pointField& p,
 763     const labelList& checkFaces,
 764     labelHashSet* setPtr
 765 )
 766 {
 767     if (maxDeg < -SMALL || maxDeg > 180+SMALL)
 768     {
 769         FatalErrorIn
 770         (
 771             "primitiveMeshGeometry::checkFaceAngles"
 772             "(const bool, const scalar, const pointField&, const labelList&"
 773             ", labelHashSet*)"
 774         )   << "maxDeg should be [0..180] but is now " << maxDeg
 775             << abort(FatalError);
 776     }
 777
 778     const scalar maxSin = Foam::sin(degToRad(maxDeg));
 779
 780     const faceList& fcs = mesh.faces();
 781
 782     scalar maxEdgeSin = 0.0;
 783
 784     label nConcave = 0;
 785
 786     label errorFaceI = -1;
 787
 788     forAll(checkFaces, i)
 789     {
 790         label faceI = checkFaces[i];
 791
 792         const face& f = fcs[faceI];
 793
 794         vector faceNormal = faceAreas[faceI];
 795         faceNormal /= mag(faceNormal) + VSMALL;
 796
 797         // Get edge from f[0] to f[size-1];
 798         vector ePrev(p[f.first()] - p[f.last()]);
 799         scalar magEPrev = mag(ePrev);
 800         ePrev /= magEPrev + VSMALL;
 801
 802         forAll(f, fp0)
 803         {
 804             // Get vertex after fp
 805             label fp1 = f.fcIndex(fp0);
 806
 807             // Normalized vector between two consecutive points
 808             vector e10(p[f[fp1]] - p[f[fp0]]);
 809             scalar magE10 = mag(e10);
 810             e10 /= magE10 + VSMALL;
 811
 812             if (magEPrev > SMALL && magE10 > SMALL)
 813             {
 814                 vector edgeNormal = ePrev ^ e10;
 815                 scalar magEdgeNormal = mag(edgeNormal);
 816
 817                 if (magEdgeNormal < maxSin)
 818                 {
 819                     // Edges (almost) aligned -> face is ok.
 820                 }
 821                 else
 822                 {
 823                     // Check normal
 824                     edgeNormal /= magEdgeNormal;
 825
 826                     if ((edgeNormal & faceNormal) < SMALL)
 827                     {
 828                         if (faceI != errorFaceI)
 829                         {
 830                             // Count only one error per face.
 831                             errorFaceI = faceI;
 832                             nConcave++;
 833                         }
 834
 835                         if (setPtr)
 836                         {
 837                             setPtr->insert(faceI);
 838                         }
 839
 840                         maxEdgeSin = max(maxEdgeSin, magEdgeNormal);
 841                     }
 842                 }
 843             }
 844
 845             ePrev = e10;
 846             magEPrev = magE10;
 847         }
 848     }
 849
 850     reduce(nConcave, sumOp<label>());
 851     reduce(maxEdgeSin, maxOp<scalar>());
 852
 853     if (report)
 854     {
 855         if (maxEdgeSin > SMALL)
 856         {
 857             scalar maxConcaveDegr =
 858                 radToDeg(Foam::asin(Foam::min(1.0, maxEdgeSin)));
 859
 860             Info<< "There are " << nConcave
 861                 << " faces with concave angles between consecutive"
 862                 << " edges. Max concave angle = " << maxConcaveDegr
 863                 << " degrees.\n" << endl;
 864         }
 865         else
 866         {
 867             Info<< "All angles in faces are convex or less than "  << maxDeg
 868                 << " degrees concave.\n" << endl;
 869         }
 870     }
 871
 872     if (nConcave > 0)
 873     {
 874         if (report)
 875         {
 876             WarningIn
 877             (
 878                 "primitiveMeshGeometry::checkFaceAngles"
 879                 "(const bool, const scalar,  const pointField&"
 880                 ", const labelList&, labelHashSet*)"
 881             )   << nConcave  << " face points with severe concave angle (> "
 882                 << maxDeg << " deg) found.\n"
 883                 << endl;
 884         }
 885
 886         return true;
 887     }
 888     else
 889     {
 890         return false;
 891     }
 892 }
 893
 894
 895 //// Check warpage of faces. Is calculated as the difference between areas of
 896 //// individual triangles and the overall area of the face (which ifself is
 897 //// is the average of the areas of the individual triangles).
 898 //bool Foam::primitiveMeshGeometry::checkFaceFlatness
 899 //(
 900 //    const bool report,
 901 //    const scalar warnFlatness,
 902 //    const primitiveMesh& mesh,
 903 //    const vectorField& faceAreas,
 904 //    const vectorField& faceCentres,
 905 //    const pointField& p,
 906 //    const labelList& checkFaces,
 907 //    labelHashSet* setPtr
 908 //)
 909 //{
 910 //    if (warnFlatness < 0 || warnFlatness > 1)
 911 //    {
 912 //        FatalErrorIn
 913 //        (
 914 //            "primitiveMeshGeometry::checkFaceFlatness"
 915 //            "(const bool, const scalar, const pointField&"
 916 //            ", const labelList&, labelHashSet*)"
 917 //        )   << "warnFlatness should be [0..1] but is now " << warnFlatness
 918 //            << abort(FatalError);
 919 //    }
 920 //
 921 //
 922 //    const faceList& fcs = mesh.faces();
 923 //
 924 //    // Areas are calculated as the sum of areas. (see
 925 //    // primitiveMeshFaceCentresAndAreas.C)
 926 //
 927 //    label nWarped = 0;
 928 //
 929 //    scalar minFlatness = GREAT;
 930 //    scalar sumFlatness = 0;
 931 //    label nSummed = 0;
 932 //
 933 //    forAll(checkFaces, i)
 934 //    {
 935 //        label faceI = checkFaces[i];
 936 //
 937 //        const face& f = fcs[faceI];
 938 //
 939 //        scalar magArea = mag(faceAreas[faceI]);
 940 //
 941 //        if (f.size() > 3 && magArea > VSMALL)
 942 //        {
 943 //            const point& fc = faceCentres[faceI];
 944 //
 945 //            // Calculate the sum of magnitude of areas and compare to
 946 //            // magnitude of sum of areas.
 947 //
 948 //            scalar sumA = 0.0;
 949 //
 950 //            forAll(f, fp)
 951 //            {
 952 //                const point& thisPoint = p[f[fp]];
 953 //                const point& nextPoint = p[f.nextLabel(fp)];
 954 //
 955 //                // Triangle around fc.
 956 //                vector n = 0.5*((nextPoint - thisPoint)^(fc - thisPoint));
 957 //                sumA += mag(n);
 958 //            }
 959 //
 960 //            scalar flatness = magArea / (sumA+VSMALL);
 961 //
 962 //            sumFlatness += flatness;
 963 //            nSummed++;
 964 //
 965 //            minFlatness = min(minFlatness, flatness);
 966 //
 967 //            if (flatness < warnFlatness)
 968 //            {
 969 //                nWarped++;
 970 //
 971 //                if (setPtr)
 972 //                {
 973 //                    setPtr->insert(faceI);
 974 //                }
 975 //            }
 976 //        }
 977 //    }
 978 //
 979 //
 980 //    reduce(nWarped, sumOp<label>());
 981 //    reduce(minFlatness, minOp<scalar>());
 982 //
 983 //    reduce(nSummed, sumOp<label>());
 984 //    reduce(sumFlatness, sumOp<scalar>());
 985 //
 986 //    if (report)
 987 //    {
 988 //        if (nSummed > 0)
 989 //        {
 990 //            Info<< "Face flatness (1 = flat, 0 = butterfly) : average = "
 991 //                << sumFlatness / nSummed << "  min = " << minFlatness << endl;
 992 //        }
 993 //
 994 //        if (nWarped> 0)
 995 //        {
 996 //            Info<< "There are " << nWarped
 997 //                << " faces with ratio between projected and actual area < "
 998 //                << warnFlatness
 999 //                << ".\nMinimum ratio (minimum flatness, maximum warpage) = "
1000 //                << minFlatness << nl << endl;
1001 //        }
1002 //        else
1003 //        {
1004 //            Info<< "All faces are flat in that the ratio between projected"
1005 //                << " and actual area is > " << warnFlatness << nl << endl;
1006 //        }
1007 //    }
1008 //
1009 //    if (nWarped > 0)
1010 //    {
1011 //        if (report)
1012 //        {
1013 //            WarningIn
1014 //            (
1015 //                "primitiveMeshGeometry::checkFaceFlatness"
1016 //                "(const bool, const scalar, const pointField&"
1017 //                ", const labelList&, labelHashSet*)"
1018 //            )   << nWarped  << " faces with severe warpage (flatness < "
1019 //                << warnFlatness << ") found.\n"
1020 //                << endl;
1021 //        }
1022 //
1023 //        return true;
1024 //    }
1025 //    else
1026 //    {
1027 //        return false;
1028 //    }
1029 //}
1030
1031
1032 // Check twist of faces. Is calculated as the difference between areas of
1033 // individual triangles and the overall area of the face (which ifself is
1034 // is the average of the areas of the individual triangles).
1035 bool Foam::primitiveMeshGeometry::checkFaceTwist
1036 (
1037     const bool report,
1038     const scalar minTwist,
1039     const primitiveMesh& mesh,
1040     const vectorField& faceAreas,
1041     const vectorField& faceCentres,
1042     const pointField& p,
1043     const labelList& checkFaces,
1044     labelHashSet* setPtr
1045 )
1046 {
1047     if (minTwist < -1-SMALL || minTwist > 1+SMALL)
1048     {
1049         FatalErrorIn
1050         (
1051             "primitiveMeshGeometry::checkFaceTwist"
1052             "(const bool, const scalar, const primitiveMesh&, const pointField&"
1053             ", const labelList&, labelHashSet*)"
1054         )   << "minTwist should be [-1..1] but is now " << minTwist
1055             << abort(FatalError);
1056     }
1057
1058
1059     const faceList& fcs = mesh.faces();
1060
1061     // Areas are calculated as the sum of areas. (see
1062     // primitiveMeshFaceCentresAndAreas.C)
1063
1064     label nWarped = 0;
1065
1066     forAll(checkFaces, i)
1067     {
1068         label faceI = checkFaces[i];
1069
1070         const face& f = fcs[faceI];
1071
1072         scalar magArea = mag(faceAreas[faceI]);
1073
1074         if (f.size() > 3 && magArea > VSMALL)
1075         {
1076             const vector nf = faceAreas[faceI] / magArea;
1077
1078             const point& fc = faceCentres[faceI];
1079
1080             forAll(f, fpI)
1081             {
1082                 vector triArea
1083                 (
1084                     triPointRef
1085                     (
1086                         p[f[fpI]],
1087                         p[f.nextLabel(fpI)],
1088                         fc
1089                     ).normal()
1090                 );
1091
1092                 scalar magTri = mag(triArea);
1093
1094                 if (magTri > VSMALL && ((nf & triArea/magTri) < minTwist))
1095                 {
1096                     nWarped++;
1097
1098                     if (setPtr)
1099                     {
1100                         setPtr->insert(faceI);
1101                     }
1102                 }
1103             }
1104         }
1105     }
1106
1107
1108     reduce(nWarped, sumOp<label>());
1109
1110     if (report)
1111     {
1112         if (nWarped> 0)
1113         {
1114             Info<< "There are " << nWarped
1115                 << " faces with cosine of the angle"
1116                 << " between triangle normal and face normal less than "
1117                 << minTwist << nl << endl;
1118         }
1119         else
1120         {
1121             Info<< "All faces are flat in that the cosine of the angle"
1122                 << " between triangle normal and face normal less than "
1123                 << minTwist << nl << endl;
1124         }
1125     }
1126
1127     if (nWarped > 0)
1128     {
1129         if (report)
1130         {
1131             WarningIn
1132             (
1133                 "primitiveMeshGeometry::checkFaceTwist"
1134                 "(const bool, const scalar, const primitiveMesh&"
1135                 ", const pointField&, const labelList&, labelHashSet*)"
1136             )   << nWarped  << " faces with severe warpage "
1137                 << "(cosine of the angle between triangle normal and "
1138                 << "face normal < " << minTwist << ") found.\n"
1139                 << endl;
1140         }
1141
1142         return true;
1143     }
1144     else
1145     {
1146         return false;
1147     }
1148 }
1149
1150
1151 bool Foam::primitiveMeshGeometry::checkFaceArea
1152 (
1153     const bool report,
1154     const scalar minArea,
1155     const primitiveMesh& mesh,
1156     const vectorField& faceAreas,
1157     const labelList& checkFaces,
1158     labelHashSet* setPtr
1159 )
1160 {
1161     label nZeroArea = 0;
1162
1163     forAll(checkFaces, i)
1164     {
1165         label faceI = checkFaces[i];
1166
1167         if (mag(faceAreas[faceI]) < minArea)
1168         {
1169             if (setPtr)
1170             {
1171                 setPtr->insert(faceI);
1172             }
1173             nZeroArea++;
1174         }
1175     }
1176
1177
1178     reduce(nZeroArea, sumOp<label>());
1179
1180     if (report)
1181     {
1182         if (nZeroArea > 0)
1183         {
1184             Info<< "There are " << nZeroArea
1185                 << " faces with area < " << minArea << '.' << nl << endl;
1186         }
1187         else
1188         {
1189             Info<< "All faces have area > " << minArea << '.' << nl << endl;
1190         }
1191     }
1192
1193     if (nZeroArea > 0)
1194     {
1195         if (report)
1196         {
1197             WarningIn
1198             (
1199                 "primitiveMeshGeometry::checkFaceArea"
1200                 "(const bool, const scalar, const primitiveMesh&"
1201                 ", const pointField&, const labelList&, labelHashSet*)"
1202             )   << nZeroArea  << " faces with area < " << minArea
1203                 << " found.\n"
1204                 << endl;
1205         }
1206
1207         return true;
1208     }
1209     else
1210     {
1211         return false;
1212     }
1213 }
1214
1215
1216 bool Foam::primitiveMeshGeometry::checkCellDeterminant
1217 (
1218     const bool report,
1219     const scalar warnDet,
1220     const primitiveMesh& mesh,
1221     const vectorField& faceAreas,
1222     const labelList& checkFaces,
1223     const labelList& affectedCells,
1224     labelHashSet* setPtr
1225 )
1226 {
1227     const cellList& cells = mesh.cells();
1228
1229     scalar minDet = GREAT;
1230     scalar sumDet = 0.0;
1231     label nSumDet = 0;
1232     label nWarnDet = 0;
1233
1234     forAll(affectedCells, i)
1235     {
1236         const cell& cFaces = cells[affectedCells[i]];
1237
1238         tensor areaSum(tensor::zero);
1239         scalar magAreaSum = 0;
1240
1241         forAll(cFaces, cFaceI)
1242         {
1243             label faceI = cFaces[cFaceI];
1244
1245             scalar magArea = mag(faceAreas[faceI]);
1246
1247             magAreaSum += magArea;
1248             areaSum += faceAreas[faceI]*(faceAreas[faceI]/magArea);
1249         }
1250
1251         scalar scaledDet = det(areaSum/magAreaSum)/0.037037037037037;
1252
1253         minDet = min(minDet, scaledDet);
1254         sumDet += scaledDet;
1255         nSumDet++;
1256
1257         if (scaledDet < warnDet)
1258         {
1259             if (setPtr)
1260             {
1261                 // Insert all faces of the cell.
1262                 forAll(cFaces, cFaceI)
1263                 {
1264                     label faceI = cFaces[cFaceI];
1265                     setPtr->insert(faceI);
1266                 }
1267             }
1268             nWarnDet++;
1269         }
1270     }
1271
1272     reduce(minDet, minOp<scalar>());
1273     reduce(sumDet, sumOp<scalar>());
1274     reduce(nSumDet, sumOp<label>());
1275     reduce(nWarnDet, sumOp<label>());
1276
1277     if (report)
1278     {
1279         if (nSumDet > 0)
1280         {
1281             Info<< "Cell determinant (1 = uniform cube) : average = "
1282                 << sumDet / nSumDet << "  min = " << minDet << endl;
1283         }
1284
1285         if (nWarnDet > 0)
1286         {
1287             Info<< "There are " << nWarnDet
1288                 << " cells with determinant < " << warnDet << '.' << nl
1289                 << endl;
1290         }
1291         else
1292         {
1293             Info<< "All faces have determinant > " << warnDet << '.' << nl
1294                 << endl;
1295         }
1296     }
1297
1298     if (nWarnDet > 0)
1299     {
1300         if (report)
1301         {
1302             WarningIn
1303             (
1304                 "primitiveMeshGeometry::checkCellDeterminant"
1305                 "(const bool, const scalar, const primitiveMesh&"
1306                 ", const pointField&, const labelList&, const labelList&"
1307                 ", labelHashSet*)"
1308             )   << nWarnDet << " cells with determinant < " << warnDet
1309                 << " found.\n"
1310                 << endl;
1311         }
1312
1313         return true;
1314     }
1315     else
1316     {
1317         return false;
1318     }
1319 }
1320
1321
1322 bool Foam::primitiveMeshGeometry::checkFaceDotProduct
1323 (
1324     const bool report,
1325     const scalar orthWarn,
1326     const labelList& checkFaces,
1327     labelHashSet* setPtr
1328 ) const
1329 {
1330     return checkFaceDotProduct
1331     (
1332         report,
1333         orthWarn,
1334         mesh_,
1335         cellCentres_,
1336         faceAreas_,
1337         checkFaces,
1338         setPtr
1339     );
1340 }
1341
1342
1343 bool Foam::primitiveMeshGeometry::checkFacePyramids
1344 (
1345     const bool report,
1346     const scalar minPyrVol,
1347     const pointField& p,
1348     const labelList& checkFaces,
1349     labelHashSet* setPtr
1350 ) const
1351 {
1352     return checkFacePyramids
1353     (
1354         report,
1355         minPyrVol,
1356         mesh_,
1357         cellCentres_,
1358         p,
1359         checkFaces,
1360         setPtr
1361     );
1362 }
1363
1364
1365 bool Foam::primitiveMeshGeometry::checkFaceSkewness
1366 (
1367     const bool report,
1368     const scalar internalSkew,
1369     const scalar boundarySkew,
1370     const labelList& checkFaces,
1371     labelHashSet* setPtr
1372 ) const
1373 {
1374     return checkFaceSkewness
1375     (
1376         report,
1377         internalSkew,
1378         boundarySkew,
1379         mesh_,
1380         cellCentres_,
1381         faceCentres_,
1382         faceAreas_,
1383         checkFaces,
1384         setPtr
1385     );
1386 }
1387
1388
1389 bool Foam::primitiveMeshGeometry::checkFaceWeights
1390 (
1391     const bool report,
1392     const scalar warnWeight,
1393     const labelList& checkFaces,
1394     labelHashSet* setPtr
1395 ) const
1396 {
1397     return checkFaceWeights
1398     (
1399         report,
1400         warnWeight,
1401         mesh_,
1402         cellCentres_,
1403         faceCentres_,
1404         faceAreas_,
1405         checkFaces,
1406         setPtr
1407     );
1408 }
1409
1410
1411 bool Foam::primitiveMeshGeometry::checkFaceAngles
1412 (
1413     const bool report,
1414     const scalar maxDeg,
1415     const pointField& p,
1416     const labelList& checkFaces,
1417     labelHashSet* setPtr
1418 ) const
1419 {
1420     return checkFaceAngles
1421     (
1422         report,
1423         maxDeg,
1424         mesh_,
1425         faceAreas_,
1426         p,
1427         checkFaces,
1428         setPtr
1429     );
1430 }
1431
1432
1433 //bool Foam::primitiveMeshGeometry::checkFaceFlatness
1434 //(
1435 //    const bool report,
1436 //    const scalar warnFlatness,
1437 //    const pointField& p,
1438 //    const labelList& checkFaces,
1439 //    labelHashSet* setPtr
1440 //) const
1441 //{
1442 //    return checkFaceFlatness
1443 //    (
1444 //        report,
1445 //        warnFlatness,
1446 //        mesh_,
1447 //        faceAreas_,
1448 //        faceCentres_,
1449 //        p,
1450 //        checkFaces,
1451 //        setPtr
1452 //    );
1453 //}
1454
1455
1456 bool Foam::primitiveMeshGeometry::checkFaceTwist
1457 (
1458     const bool report,
1459     const scalar minTwist,
1460     const pointField& p,
1461     const labelList& checkFaces,
1462     labelHashSet* setPtr
1463 ) const
1464 {
1465     return checkFaceTwist
1466     (
1467         report,
1468         minTwist,
1469         mesh_,
1470         faceAreas_,
1471         faceCentres_,
1472         p,
1473         checkFaces,
1474         setPtr
1475     );
1476 }
1477
1478
1479 bool Foam::primitiveMeshGeometry::checkFaceArea
1480 (
1481     const bool report,
1482     const scalar minArea,
1483     const labelList& checkFaces,
1484     labelHashSet* setPtr
1485 ) const
1486 {
1487     return checkFaceArea
1488     (
1489         report,
1490         minArea,
1491         mesh_,
1492         faceAreas_,
1493         checkFaces,
1494         setPtr
1495     );
1496 }
1497
1498
1499 bool Foam::primitiveMeshGeometry::checkCellDeterminant
1500 (
1501     const bool report,
1502     const scalar warnDet,
1503     const labelList& checkFaces,
1504     const labelList& affectedCells,
1505     labelHashSet* setPtr
1506 ) const
1507 {
1508     return checkCellDeterminant
1509     (
1510         report,
1511         warnDet,
1512         mesh_,
1513         faceAreas_,
1514         checkFaces,
1515         affectedCells,
1516         setPtr
1517     );
1518 }
1519
1520
1521 // ************************************************************************* //