src/dynamicMesh/motionSmoother/motionSmootherCheck.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2010 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
  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 "motionSmoother.H"
  27 #include <dynamicMesh/polyMeshGeometry.H>
  28 #include <OpenFOAM/IOmanip.H>
  29
  30 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  31
  32 bool Foam::motionSmoother::checkMesh
  33 (
  34     const bool report,
  35     const polyMesh& mesh,
  36     const dictionary& dict,
  37     const labelList& checkFaces,
  38     labelHashSet& wrongFaces
  39 )
  40 {
  41     List<labelPair> emptyBaffles;
  42     return checkMesh
  43     (
  44         report,
  45         mesh,
  46         dict,
  47         checkFaces,
  48         emptyBaffles,
  49         wrongFaces
  50     );
  51 }
  52
  53 bool Foam::motionSmoother::checkMesh
  54 (
  55     const bool report,
  56     const polyMesh& mesh,
  57     const dictionary& dict,
  58     const labelList& checkFaces,
  59     const List<labelPair>& baffles,
  60     labelHashSet& wrongFaces
  61 )
  62 {
  63     const scalar maxNonOrtho
  64     (
  65         readScalar(dict.lookup("maxNonOrtho", true))
  66     );
  67     const scalar minVol
  68     (
  69         readScalar(dict.lookup("minVol", true))
  70     );
  71     const scalar maxConcave
  72     (
  73         readScalar(dict.lookup("maxConcave", true))
  74     );
  75     const scalar minArea
  76     (
  77         readScalar(dict.lookup("minArea", true))
  78     );
  79     const scalar maxIntSkew
  80     (
  81         readScalar(dict.lookup("maxInternalSkewness", true))
  82     );
  83     const scalar maxBounSkew
  84     (
  85         readScalar(dict.lookup("maxBoundarySkewness", true))
  86     );
  87     const scalar minWeight
  88     (
  89         readScalar(dict.lookup("minFaceWeight", true))
  90     );
  91     const scalar minVolRatio
  92     (
  93         readScalar(dict.lookup("minVolRatio", true))
  94     );
  95     const scalar minTwist
  96     (
  97         readScalar(dict.lookup("minTwist", true))
  98     );
  99     const scalar minTriangleTwist
 100     (
 101         readScalar(dict.lookup("minTriangleTwist", true))
 102     );
 103     const scalar minDet
 104     (
 105         readScalar(dict.lookup("minDeterminant", true))
 106     );
 107
 108     label nWrongFaces = 0;
 109
 110     Info<< "Checking faces in error :" << endl;
 111     //Pout.setf(ios_base::left);
 112
 113     if (maxNonOrtho < 180.0-SMALL)
 114     {
 115         polyMeshGeometry::checkFaceDotProduct
 116         (
 117             report,
 118             maxNonOrtho,
 119             mesh,
 120             mesh.cellCentres(),
 121             mesh.faceAreas(),
 122             checkFaces,
 123             baffles,
 124             &wrongFaces
 125         );
 126
 127         label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
 128
 129         Info<< "    non-orthogonality > "
 130             << setw(3) << maxNonOrtho
 131             << " degrees                        : "
 132             << nNewWrongFaces-nWrongFaces << endl;
 133
 134         nWrongFaces = nNewWrongFaces;
 135     }
 136
 137     if (minVol > -GREAT)
 138     {
 139         polyMeshGeometry::checkFacePyramids
 140         (
 141             report,
 142             minVol,
 143             mesh,
 144             mesh.cellCentres(),
 145             mesh.points(),
 146             checkFaces,
 147             baffles,
 148             &wrongFaces
 149         );
 150
 151         label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
 152
 153         Info<< "    faces with face pyramid volume < "
 154             << setw(5) << minVol << "                 : "
 155             << nNewWrongFaces-nWrongFaces << endl;
 156
 157         nWrongFaces = nNewWrongFaces;
 158     }
 159
 160     if (maxConcave < 180.0-SMALL)
 161     {
 162         polyMeshGeometry::checkFaceAngles
 163         (
 164             report,
 165             maxConcave,
 166             mesh,
 167             mesh.faceAreas(),
 168             mesh.points(),
 169             checkFaces,
 170             &wrongFaces
 171         );
 172
 173         label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
 174
 175         Info<< "    faces with concavity > "
 176             << setw(3) << maxConcave
 177             << " degrees                     : "
 178             << nNewWrongFaces-nWrongFaces << endl;
 179
 180         nWrongFaces = nNewWrongFaces;
 181     }
 182
 183     if (minArea > -SMALL)
 184     {
 185         polyMeshGeometry::checkFaceArea
 186         (
 187             report,
 188             minArea,
 189             mesh,
 190             mesh.faceAreas(),
 191             checkFaces,
 192             &wrongFaces
 193         );
 194
 195         label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
 196
 197         Info<< "    faces with area < "
 198             << setw(5) << minArea
 199             << " m^2                            : "
 200             << nNewWrongFaces-nWrongFaces << endl;
 201
 202         nWrongFaces = nNewWrongFaces;
 203     }
 204
 205     if (maxIntSkew > 0 || maxBounSkew > 0)
 206     {
 207         polyMeshGeometry::checkFaceSkewness
 208         (
 209             report,
 210             maxIntSkew,
 211             maxBounSkew,
 212             mesh,
 213             mesh.cellCentres(),
 214             mesh.faceCentres(),
 215             mesh.faceAreas(),
 216             checkFaces,
 217             baffles,
 218             &wrongFaces
 219         );
 220
 221         label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
 222
 223         Info<< "    faces with skewness > "
 224             << setw(3) << maxIntSkew
 225             << " (internal) or " << setw(3) << maxBounSkew
 226             << " (boundary) : " << nNewWrongFaces-nWrongFaces << endl;
 227
 228         nWrongFaces = nNewWrongFaces;
 229     }
 230
 231     if (minWeight >= 0 && minWeight < 1)
 232     {
 233         polyMeshGeometry::checkFaceWeights
 234         (
 235             report,
 236             minWeight,
 237             mesh,
 238             mesh.cellCentres(),
 239             mesh.faceCentres(),
 240             mesh.faceAreas(),
 241             checkFaces,
 242             baffles,
 243             &wrongFaces
 244         );
 245
 246         label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
 247
 248         Info<< "    faces with interpolation weights (0..1)  < "
 249             << setw(5) << minWeight
 250             << "       : "
 251             << nNewWrongFaces-nWrongFaces << endl;
 252
 253         nWrongFaces = nNewWrongFaces;
 254     }
 255
 256     if (minVolRatio >= 0)
 257     {
 258         polyMeshGeometry::checkVolRatio
 259         (
 260             report,
 261             minVolRatio,
 262             mesh,
 263             mesh.cellVolumes(),
 264             checkFaces,
 265             baffles,
 266             &wrongFaces
 267         );
 268
 269         label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
 270
 271         Info<< "    faces with volume ratio of neighbour cells < "
 272             << setw(5) << minVolRatio
 273             << "     : "
 274             << nNewWrongFaces-nWrongFaces << endl;
 275
 276         nWrongFaces = nNewWrongFaces;
 277     }
 278
 279     if (minTwist > -1)
 280     {
 281         //Pout<< "Checking face twist: dot product of face normal "
 282         //    << "with face triangle normals" << endl;
 283         polyMeshGeometry::checkFaceTwist
 284         (
 285             report,
 286             minTwist,
 287             mesh,
 288             mesh.cellCentres(),
 289             mesh.faceAreas(),
 290             mesh.faceCentres(),
 291             mesh.points(),
 292             checkFaces,
 293             &wrongFaces
 294         );
 295
 296         label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
 297
 298         Info<< "    faces with face twist < "
 299             << setw(5) << minTwist
 300             << "                          : "
 301             << nNewWrongFaces-nWrongFaces << endl;
 302
 303         nWrongFaces = nNewWrongFaces;
 304     }
 305
 306     if (minTriangleTwist > -1)
 307     {
 308         //Pout<< "Checking triangle twist: dot product of consecutive triangle"
 309         //    << " normals resulting from face-centre decomposition" << endl;
 310         polyMeshGeometry::checkTriangleTwist
 311         (
 312             report,
 313             minTriangleTwist,
 314             mesh,
 315             mesh.faceAreas(),
 316             mesh.faceCentres(),
 317             mesh.points(),
 318             checkFaces,
 319             &wrongFaces
 320         );
 321
 322         label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
 323
 324         Info<< "    faces with triangle twist < "
 325             << setw(5) << minTriangleTwist
 326             << "                      : "
 327             << nNewWrongFaces-nWrongFaces << endl;
 328
 329         nWrongFaces = nNewWrongFaces;
 330     }
 331
 332     if (minDet > -1)
 333     {
 334         polyMeshGeometry::checkCellDeterminant
 335         (
 336             report,
 337             minDet,
 338             mesh,
 339             mesh.faceAreas(),
 340             checkFaces,
 341             polyMeshGeometry::affectedCells(mesh, checkFaces),
 342             &wrongFaces
 343         );
 344
 345         label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
 346
 347         Info<< "    faces on cells with determinant < "
 348             << setw(5) << minDet << "                : "
 349             << nNewWrongFaces-nWrongFaces << endl;
 350
 351         nWrongFaces = nNewWrongFaces;
 352     }
 353
 354     //Pout.setf(ios_base::right);
 355
 356     return nWrongFaces > 0;
 357 }
 358
 359
 360 bool Foam::motionSmoother::checkMesh
 361 (
 362     const bool report,
 363     const polyMesh& mesh,
 364     const dictionary& dict,
 365     labelHashSet& wrongFaces
 366 )
 367 {
 368     return checkMesh
 369     (
 370         report,
 371         mesh,
 372         dict,
 373         identity(mesh.nFaces()),
 374         wrongFaces
 375     );
 376 }
 377
 378 bool Foam::motionSmoother::checkMesh
 379 (
 380     const bool report,
 381     const dictionary& dict,
 382     const polyMeshGeometry& meshGeom,
 383     const labelList& checkFaces,
 384     labelHashSet& wrongFaces
 385 )
 386 {
 387     List<labelPair> emptyBaffles;
 388
 389     return checkMesh
 390     (
 391         report,
 392         dict,
 393         meshGeom,
 394         checkFaces,
 395         emptyBaffles,
 396         wrongFaces
 397      );
 398 }
 399
 400
 401 bool Foam::motionSmoother::checkMesh
 402 (
 403     const bool report,
 404     const dictionary& dict,
 405     const polyMeshGeometry& meshGeom,
 406     const labelList& checkFaces,
 407     const List<labelPair>& baffles,
 408     labelHashSet& wrongFaces
 409 )
 410 {
 411     const scalar maxNonOrtho
 412     (
 413         readScalar(dict.lookup("maxNonOrtho", true))
 414     );
 415     const scalar minVol
 416     (
 417         readScalar(dict.lookup("minVol", true))
 418     );
 419     const scalar maxConcave
 420     (
 421         readScalar(dict.lookup("maxConcave", true))
 422     );
 423     const scalar minArea
 424     (
 425         readScalar(dict.lookup("minArea", true))
 426     );
 427     const scalar maxIntSkew
 428     (
 429         readScalar(dict.lookup("maxInternalSkewness", true))
 430     );
 431     const scalar maxBounSkew
 432     (
 433         readScalar(dict.lookup("maxBoundarySkewness", true))
 434     );
 435     const scalar minWeight
 436     (
 437         readScalar(dict.lookup("minFaceWeight", true))
 438     );
 439     const scalar minVolRatio
 440     (
 441         readScalar(dict.lookup("minVolRatio", true))
 442     );
 443     const scalar minTwist
 444     (
 445         readScalar(dict.lookup("minTwist", true))
 446     );
 447     const scalar minTriangleTwist
 448     (
 449         readScalar(dict.lookup("minTriangleTwist", true))
 450     );
 451     const scalar minDet
 452     (
 453         readScalar(dict.lookup("minDeterminant", true))
 454     );
 455
 456     label nWrongFaces = 0;
 457
 458     Info<< "Checking faces in error :" << endl;
 459     //Pout.setf(ios_base::left);
 460
 461     if (maxNonOrtho < 180.0-SMALL)
 462     {
 463         meshGeom.checkFaceDotProduct
 464         (
 465             report,
 466             maxNonOrtho,
 467             checkFaces,
 468             baffles,
 469             &wrongFaces
 470         );
 471
 472         label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
 473
 474         Info<< "    non-orthogonality > "
 475             << setw(3) << maxNonOrtho
 476             << " degrees                        : "
 477             << nNewWrongFaces-nWrongFaces << endl;
 478
 479         nWrongFaces = nNewWrongFaces;
 480     }
 481
 482     if (minVol > -GREAT)
 483     {
 484         meshGeom.checkFacePyramids
 485         (
 486             report,
 487             minVol,
 488             meshGeom.mesh().points(),
 489             checkFaces,
 490             baffles,
 491             &wrongFaces
 492         );
 493
 494         label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
 495
 496         Info<< "    faces with face pyramid volume < "
 497             << setw(5) << minVol << "                 : "
 498             << nNewWrongFaces-nWrongFaces << endl;
 499
 500         nWrongFaces = nNewWrongFaces;
 501     }
 502
 503     if (maxConcave < 180.0-SMALL)
 504     {
 505         meshGeom.checkFaceAngles
 506         (
 507             report,
 508             maxConcave,
 509             meshGeom.mesh().points(),
 510             checkFaces,
 511             &wrongFaces
 512         );
 513
 514         label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
 515
 516         Info<< "    faces with concavity > "
 517             << setw(3) << maxConcave
 518             << " degrees                     : "
 519             << nNewWrongFaces-nWrongFaces << endl;
 520
 521         nWrongFaces = nNewWrongFaces;
 522     }
 523
 524     if (minArea > -SMALL)
 525     {
 526         meshGeom.checkFaceArea(report, minArea, checkFaces, &wrongFaces);
 527
 528         label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
 529
 530         Info<< "    faces with area < "
 531             << setw(5) << minArea
 532             << " m^2                            : "
 533             << nNewWrongFaces-nWrongFaces << endl;
 534
 535         nWrongFaces = nNewWrongFaces;
 536     }
 537
 538     if (maxIntSkew > 0 || maxBounSkew > 0)
 539     {
 540         meshGeom.checkFaceSkewness
 541         (
 542             report,
 543             maxIntSkew,
 544             maxBounSkew,
 545             checkFaces,
 546             baffles,
 547             &wrongFaces
 548         );
 549
 550         label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
 551
 552         Info<< "    faces with skewness > "
 553             << setw(3) << maxIntSkew
 554             << " (internal) or " << setw(3) << maxBounSkew
 555             << " (boundary) : " << nNewWrongFaces-nWrongFaces << endl;
 556
 557         nWrongFaces = nNewWrongFaces;
 558     }
 559
 560     if (minWeight >= 0 && minWeight < 1)
 561     {
 562         meshGeom.checkFaceWeights
 563         (
 564             report,
 565             minWeight,
 566             checkFaces,
 567             baffles,
 568             &wrongFaces
 569         );
 570
 571         label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
 572
 573         Info<< "    faces with interpolation weights (0..1)  < "
 574             << setw(5) << minWeight
 575             << "       : "
 576             << nNewWrongFaces-nWrongFaces << endl;
 577
 578         nWrongFaces = nNewWrongFaces;
 579     }
 580
 581     if (minVolRatio >= 0)
 582     {
 583         meshGeom.checkVolRatio
 584         (
 585             report,
 586             minVolRatio,
 587             checkFaces,
 588             baffles,
 589             &wrongFaces
 590         );
 591
 592         label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
 593
 594         Info<< "    faces with volume ratio of neighbour cells < "
 595             << setw(5) << minVolRatio
 596             << "     : "
 597             << nNewWrongFaces-nWrongFaces << endl;
 598
 599         nWrongFaces = nNewWrongFaces;
 600     }
 601
 602     if (minTwist > -1)
 603     {
 604         //Pout<< "Checking face twist: dot product of face normal "
 605         //    << "with face triangle normals" << endl;
 606         meshGeom.checkFaceTwist
 607         (
 608             report,
 609             minTwist,
 610             meshGeom.mesh().points(),
 611             checkFaces,
 612             &wrongFaces
 613         );
 614
 615         label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
 616
 617         Info<< "    faces with face twist < "
 618             << setw(5) << minTwist
 619             << "                          : "
 620             << nNewWrongFaces-nWrongFaces << endl;
 621
 622         nWrongFaces = nNewWrongFaces;
 623     }
 624
 625     if (minTriangleTwist > -1)
 626     {
 627         //Pout<< "Checking triangle twist: dot product of consecutive triangle"
 628         //    << " normals resulting from face-centre decomposition" << endl;
 629         meshGeom.checkTriangleTwist
 630         (
 631             report,
 632             minTriangleTwist,
 633             meshGeom.mesh().points(),
 634             checkFaces,
 635             &wrongFaces
 636         );
 637
 638         label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
 639
 640         Info<< "    faces with triangle twist < "
 641             << setw(5) << minTriangleTwist
 642             << "                      : "
 643             << nNewWrongFaces-nWrongFaces << endl;
 644
 645         nWrongFaces = nNewWrongFaces;
 646     }
 647
 648     if (minDet > -1)
 649     {
 650         meshGeom.checkCellDeterminant
 651         (
 652             report,
 653             minDet,
 654             checkFaces,
 655             meshGeom.affectedCells(meshGeom.mesh(), checkFaces),
 656             &wrongFaces
 657         );
 658
 659         label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
 660
 661         Info<< "    faces on cells with determinant < "
 662             << setw(5) << minDet << "                : "
 663             << nNewWrongFaces-nWrongFaces << endl;
 664
 665         nWrongFaces = nNewWrongFaces;
 666     }
 667
 668     //Pout.setf(ios_base::right);
 669
 670     return nWrongFaces > 0;
 671 }
 672
 673
 674 // ************************ vim: set sw=4 sts=4 et: ************************ //