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