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