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