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