src/mesh/cfMesh/utilities/smoothers/geometry/meshOptimizer/optimizeMeshFV.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | cfMesh: A library for mesh generation
   4    \\    /   O peration     |
   5     \\  /    A nd           | Author: Franjo Juretic (franjo.juretic@c-fields.com)
   6      \\/     M anipulation  | Copyright (C) Creative Fields, Ltd.
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of cfMesh.
  10
  11     cfMesh 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     cfMesh 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 cfMesh.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25
  26 \*---------------------------------------------------------------------------*/
  27
  28 #include <stdexcept>
  29
  30 #include "demandDrivenData.H"
  31 #include "meshOptimizer.H"
  32 #include "polyMeshGenAddressing.H"
  33 #include "polyMeshGenChecks.H"
  34 #include "partTetMesh.H"
  35 #include "HashSet.H"
  36
  37 #include "tetMeshOptimisation.H"
  38 #include "boundaryLayerOptimisation.H"
  39 #include "refineBoundaryLayers.H"
  40 #include "meshSurfaceEngine.H"
  41
  42 //#define DEBUGSmooth
  43
  44 # ifdef DEBUGSmooth
  45 #include "helperFunctions.H"
  46 #include "polyMeshGenModifier.H"
  47 # endif
  48
  49 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  50
  51 namespace Foam
  52 {
  53
  54 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  55
  56 void meshOptimizer::untangleMeshFV
  57 (
  58     const label maxNumGlobalIterations,
  59     const label maxNumIterations,
  60     const label maxNumSurfaceIterations,
  61     const bool relaxedCheck
  62 )
  63 {
  64     Info << "Starting untangling the mesh" << endl;
  65
  66     # ifdef DEBUGSmooth
  67     partTetMesh tm(mesh_);
  68     forAll(tm.tets(), tetI)
  69         if( tm.tets()[tetI].mag(tm.points()) < 0.0 )
  70             Info << "Tet " << tetI << " is inverted!" << endl;
  71     polyMeshGen tetPolyMesh(mesh_.returnTime());
  72     tm.createPolyMesh(tetPolyMesh);
  73     polyMeshGenModifier(tetPolyMesh).removeUnusedVertices();
  74     forAll(tm.smoothVertex(), pI)
  75         if( !tm.smoothVertex()[pI] )
  76             Info << "Point " << pI << " cannot be moved!" << endl;
  77
  78     const VRWGraph& pTets = tm.pointTets();
  79     forAll(pTets, pointI)
  80     {
  81         const LongList<partTet>& tets = tm.tets();
  82         forAllRow(pTets, pointI, i)
  83             if( tets[pTets(pointI, i)].whichPosition(pointI) < 0 )
  84                 FatalError << "Wrong partTet" << abort(FatalError);
  85
  86         partTetMeshSimplex simplex(tm, pointI);
  87     }
  88
  89     boolList boundaryVertex(tetPolyMesh.points().size(), false);
  90     const labelList& neighbour = tetPolyMesh.neighbour();
  91     forAll(neighbour, faceI)
  92         if( neighbour[faceI] == -1 )
  93         {
  94             const face& f = tetPolyMesh.faces()[faceI];
  95
  96             forAll(f, pI)
  97                 boundaryVertex[f[pI]] = true;
  98         }
  99
 100     forAll(boundaryVertex, pI)
 101     {
 102         if( boundaryVertex[pI] && tm.smoothVertex()[pI] )
 103             FatalErrorIn
 104             (
 105                 "void meshOptimizer::untangleMeshFV()"
 106             ) << "Boundary vertex should not be moved!" << abort(FatalError);
 107     }
 108     # endif
 109
 110     label nBadFaces, nGlobalIter(0), nIter;
 111
 112     const faceListPMG& faces = mesh_.faces();
 113
 114     boolList changedFace(faces.size(), true);
 115
 116     //- check if any points in the tet mesh shall not move
 117     labelLongList lockedPoints;
 118     forAll(vertexLocation_, pointI)
 119     {
 120         if( vertexLocation_[pointI] & LOCKED )
 121             lockedPoints.append(pointI);
 122     }
 123
 124     labelHashSet badFaces;
 125
 126     do
 127     {
 128         nIter = 0;
 129
 130         label minNumBadFaces(10 * faces.size()), minIter(-1);
 131         do
 132         {
 133             if( !relaxedCheck )
 134             {
 135                 nBadFaces =
 136                     polyMeshGenChecks::findBadFaces
 137                     (
 138                         mesh_,
 139                         badFaces,
 140                         false,
 141                         &changedFace
 142                     );
 143             }
 144             else
 145             {
 146                 nBadFaces =
 147                     polyMeshGenChecks::findBadFacesRelaxed
 148                     (
 149                         mesh_,
 150                         badFaces,
 151                         false,
 152                         &changedFace
 153                     );
 154             }
 155
 156             Info << "Iteration " << nIter
 157                 << ". Number of bad faces is " << nBadFaces << endl;
 158
 159             //- perform optimisation
 160             if( nBadFaces == 0 )
 161                 break;
 162
 163             if( nBadFaces < minNumBadFaces )
 164             {
 165                 minNumBadFaces = nBadFaces;
 166                 minIter = nIter;
 167             }
 168
 169             //- create a tet mesh from the mesh and the labels of bad faces
 170             partTetMesh tetMesh
 171             (
 172                 mesh_,
 173                 lockedPoints,
 174                 badFaces,
 175                 (nGlobalIter / 2) + 1
 176             );
 177
 178             //- construct tetMeshOptimisation and improve positions of
 179             //- points in the tet mesh
 180             tetMeshOptimisation tmo(tetMesh);
 181
 182             tmo.optimiseUsingKnuppMetric();
 183
 184             tmo.optimiseUsingMeshUntangler();
 185
 186             tmo.optimiseUsingVolumeOptimizer();
 187
 188             //- update points in the mesh from the coordinates in the tet mesh
 189             tetMesh.updateOrigMesh(&changedFace);
 190
 191         } while( (nIter < minIter+5) && (++nIter < maxNumIterations) );
 192
 193         if( (nBadFaces == 0) || (++nGlobalIter >= maxNumGlobalIterations) )
 194             break;
 195
 196         // move boundary vertices
 197         nIter = 0;
 198
 199         while( nIter++ < maxNumSurfaceIterations );
 200         {
 201             if( !relaxedCheck )
 202             {
 203                 nBadFaces =
 204                     polyMeshGenChecks::findBadFaces
 205                     (
 206                         mesh_,
 207                         badFaces,
 208                         false,
 209                         &changedFace
 210                     );
 211             }
 212             else
 213             {
 214                 nBadFaces =
 215                     polyMeshGenChecks::findBadFacesRelaxed
 216                     (
 217                         mesh_,
 218                         badFaces,
 219                         false,
 220                         &changedFace
 221                     );
 222             }
 223
 224             Info << "Iteration " << nIter
 225                 << ". Number of bad faces is " << nBadFaces << endl;
 226
 227             //- perform optimisation
 228             if( nBadFaces == 0 )
 229             {
 230                 break;
 231             }
 232             else if( enforceConstraints_ )
 233             {
 234                 const label subsetId =
 235                     mesh_.addPointSubset(badPointsSubsetName_);
 236
 237                 forAllConstIter(labelHashSet, badFaces, it)
 238                 {
 239                     const face& f = faces[it.key()];
 240                     forAll(f, pI)
 241                         mesh_.addPointToSubset(subsetId, f[pI]);
 242                 }
 243
 244                 WarningIn
 245                 (
 246                     "void meshOptimizer::untangleMeshFV()"
 247                 ) << "Writing mesh with " << badPointsSubsetName_
 248                   << " subset. These points cannot be untangled"
 249                   << " without sacrificing geometry constraints. Exitting.."
 250                   << endl;
 251
 252                 returnReduce(1, sumOp<label>());
 253                 throw std::logic_error
 254                 (
 255                     "void meshOptimizer::untangleMeshFV()"
 256                     "Cannot untangle mesh!!"
 257                 );
 258             }
 259
 260             //- create tethrahedral mesh from the cells which shall be smoothed
 261             partTetMesh tetMesh(mesh_, lockedPoints, badFaces, 0);
 262
 263             //- contruct tetMeshOptimisation
 264             tetMeshOptimisation tmo(tetMesh);
 265
 266             if( nGlobalIter < 2 )
 267             {
 268                 //- the point stays in the plane determined by the point normal
 269                 tmo.optimiseBoundaryVolumeOptimizer(true);
 270             }
 271             else if( nGlobalIter < 5 )
 272             {
 273                 //- move points without any constraints on the movement
 274                 tmo.optimiseBoundarySurfaceLaplace();
 275             }
 276             else
 277             {
 278                 //- move boundary points without any constraints
 279                 tmo.optimiseBoundaryVolumeOptimizer(false);
 280             }
 281
 282             tetMesh.updateOrigMesh(&changedFace);
 283
 284         }
 285
 286     } while( nBadFaces );
 287
 288     if( nBadFaces != 0 )
 289     {
 290         label subsetId = mesh_.faceSubsetIndex("badFaces");
 291         if( subsetId >= 0 )
 292             mesh_.removeFaceSubset(subsetId);
 293         subsetId = mesh_.addFaceSubset("badFaces");
 294
 295         forAllConstIter(labelHashSet, badFaces, it)
 296             mesh_.addFaceToSubset(subsetId, it.key());
 297     }
 298
 299     Info << "Finished untangling the mesh" << endl;
 300 }
 301
 302 void meshOptimizer::optimizeBoundaryLayer(const bool addBufferLayer)
 303 {
 304     if( mesh_.returnTime().foundObject<IOdictionary>("meshDict") )
 305     {
 306         const dictionary& meshDict =
 307             mesh_.returnTime().lookupObject<IOdictionary>("meshDict");
 308
 309         bool smoothLayer(false);
 310
 311         if( meshDict.found("boundaryLayers") )
 312         {
 313             const dictionary& layersDict = meshDict.subDict("boundaryLayers");
 314
 315             if( layersDict.found("optimiseLayer") )
 316                 smoothLayer = readBool(layersDict.lookup("optimiseLayer"));
 317         }
 318
 319         if( !smoothLayer )
 320             return;
 321
 322         if( addBufferLayer )
 323         {
 324             //- create a buffer layer which will not be modified by the smoother
 325             refineBoundaryLayers refLayers(mesh_);
 326
 327             refineBoundaryLayers::readSettings(meshDict, refLayers);
 328
 329             refLayers.activateSpecialMode();
 330
 331             refLayers.refineLayers();
 332
 333             clearSurface();
 334             calculatePointLocations();
 335         }
 336
 337         Info << "Starting optimising boundary layer" << endl;
 338
 339         const meshSurfaceEngine& mse = meshSurface();
 340         const labelList& faceOwner = mse.faceOwners();
 341
 342         boundaryLayerOptimisation optimiser(mesh_, mse);
 343
 344         boundaryLayerOptimisation::readSettings(meshDict, optimiser);
 345
 346         optimiser.optimiseLayer();
 347
 348         //- check if the bnd layer is tangled somewhere
 349         labelLongList bndLayerCells;
 350         const boolList& baseFace = optimiser.isBaseFace();
 351
 352         # ifdef DEBUGSmooth
 353         const label blCellsId = mesh_.addCellSubset("blCells");
 354         # endif
 355
 356         forAll(baseFace, bfI)
 357         {
 358             if( baseFace[bfI] )
 359             {
 360                 bndLayerCells.append(faceOwner[bfI]);
 361
 362                 # ifdef DEBUGSmooth
 363                 mesh_.addCellToSubset(blCellsId, faceOwner[bfI]);
 364                 # endif
 365             }
 366         }
 367
 368         clearSurface();
 369         mesh_.clearAddressingData();
 370
 371         //- lock boundary layer points, faces and cells
 372         lockCells(bndLayerCells);
 373
 374         # ifdef DEBUGSmooth
 375         pointField origPoints(mesh_.points().size());
 376         forAll(origPoints, pI)
 377             origPoints[pI] = mesh_.points()[pI];
 378         # endif
 379
 380         //- optimize mesh quality
 381         optimizeMeshFV(5, 1, 50, 0);
 382
 383         //- untangle remaining faces and lock the boundary layer cells
 384         untangleMeshFV(2, 50, 0);
 385
 386         # ifdef DEBUGSmooth
 387         forAll(vertexLocation_, pI)
 388         {
 389             if( vertexLocation_[pI] & LOCKED )
 390             {
 391                 if( mag(origPoints[pI] - mesh_.points()[pI]) > SMALL )
 392                     FatalError << "Locked points were moved"
 393                                << abort(FatalError);
 394             }
 395         }
 396         # endif
 397
 398         //- unlock bnd layer points
 399         removeUserConstraints();
 400
 401         Info << "Finished optimising boundary layer" << endl;
 402     }
 403 }
 404
 405 void meshOptimizer::untangleBoundaryLayer()
 406 {
 407     bool untangleLayer(true);
 408     if( mesh_.returnTime().foundObject<IOdictionary>("meshDict") )
 409     {
 410         const dictionary& meshDict =
 411             mesh_.returnTime().lookupObject<IOdictionary>("meshDict");
 412
 413         if( meshDict.found("boundaryLayers") )
 414         {
 415             const dictionary& layersDict = meshDict.subDict("boundaryLayers");
 416
 417             if( layersDict.found("untangleLayers") )
 418             {
 419                 untangleLayer =
 420                     readBool(layersDict.lookup("untangleLayers"));
 421             }
 422         }
 423     }
 424
 425     if( !untangleLayer )
 426     {
 427         labelHashSet badFaces;
 428         polyMeshGenChecks::checkFacePyramids(mesh_, false, VSMALL, &badFaces);
 429
 430         const label nInvalidFaces =
 431             returnReduce(badFaces.size(), sumOp<label>());
 432
 433         if( nInvalidFaces != 0 )
 434         {
 435             const labelList& owner = mesh_.owner();
 436             const labelList& neighbour = mesh_.neighbour();
 437
 438             const label badBlCellsId =
 439                 mesh_.addCellSubset("invalidBoundaryLayerCells");
 440
 441             forAllConstIter(labelHashSet, badFaces, it)
 442             {
 443                 mesh_.addCellToSubset(badBlCellsId, owner[it.key()]);
 444
 445                 if( neighbour[it.key()] < 0 )
 446                     continue;
 447
 448                 mesh_.addCellToSubset(badBlCellsId, neighbour[it.key()]);
 449             }
 450
 451             returnReduce(1, sumOp<label>());
 452
 453             throw std::logic_error
 454             (
 455                 "void meshOptimizer::untangleBoundaryLayer()"
 456                 "Found invalid faces in the boundary layer."
 457                 " Cannot untangle mesh!!"
 458             );
 459         }
 460     }
 461     else
 462     {
 463         optimizeLowQualityFaces();
 464         removeUserConstraints();
 465         untangleMeshFV(2, 50, 1, true);
 466     }
 467 }
 468
 469 void meshOptimizer::optimizeLowQualityFaces(const label maxNumIterations)
 470 {
 471     label nBadFaces, nIter(0);
 472
 473     const faceListPMG& faces = mesh_.faces();
 474     boolList changedFace(faces.size(), true);
 475
 476     //- check if any points in the tet mesh shall not move
 477     labelLongList lockedPoints;
 478     forAll(vertexLocation_, pointI)
 479     {
 480         if( vertexLocation_[pointI] & LOCKED )
 481             lockedPoints.append(pointI);
 482     }
 483
 484     do
 485     {
 486         labelHashSet lowQualityFaces;
 487         nBadFaces =
 488             polyMeshGenChecks::findLowQualityFaces
 489             (
 490                 mesh_,
 491                 lowQualityFaces,
 492                 false,
 493                 &changedFace
 494             );
 495
 496         changedFace = false;
 497         forAllConstIter(labelHashSet, lowQualityFaces, it)
 498             changedFace[it.key()] = true;
 499
 500         Info << "Iteration " << nIter
 501             << ". Number of bad faces is " << nBadFaces << endl;
 502
 503         //- perform optimisation
 504         if( nBadFaces == 0 )
 505             break;
 506
 507         partTetMesh tetMesh(mesh_, lockedPoints, lowQualityFaces, 2);
 508
 509         //- construct tetMeshOptimisation and improve positions
 510         //- of points in the tet mesh
 511         tetMeshOptimisation tmo(tetMesh);
 512
 513         tmo.optimiseUsingVolumeOptimizer();
 514
 515         //- update points in the mesh from the new coordinates in the tet mesh
 516         tetMesh.updateOrigMesh(&changedFace);
 517
 518     } while( ++nIter < maxNumIterations );
 519 }
 520
 521 void meshOptimizer::optimizeMeshNearBoundaries
 522 (
 523     const label maxNumIterations,
 524     const label numLayersOfCells
 525 )
 526 {
 527     label nIter(0);
 528
 529     const faceListPMG& faces = mesh_.faces();
 530     boolList changedFace(faces.size(), true);
 531
 532     //- check if any points in the tet mesh shall not move
 533     labelLongList lockedPoints;
 534     forAll(vertexLocation_, pointI)
 535     {
 536         if( vertexLocation_[pointI] & LOCKED )
 537             lockedPoints.append(pointI);
 538     }
 539
 540     partTetMesh tetMesh(mesh_, lockedPoints, numLayersOfCells);
 541     tetMeshOptimisation tmo(tetMesh);
 542     Info << "Iteration:" << flush;
 543     do
 544     {
 545         tmo.optimiseUsingVolumeOptimizer(1);
 546
 547         tetMesh.updateOrigMesh(&changedFace);
 548
 549         Info << "." << flush;
 550
 551     } while( ++nIter < maxNumIterations );
 552
 553     Info << endl;
 554 }
 555
 556 void meshOptimizer::optimizeMeshFV
 557 (
 558     const label numLaplaceIterations,
 559     const label maxNumGlobalIterations,
 560     const label maxNumIterations,
 561     const label maxNumSurfaceIterations
 562 )
 563 {
 564     Info << "Starting smoothing the mesh" << endl;
 565
 566     laplaceSmoother lps(mesh_, vertexLocation_);
 567     lps.optimizeLaplacianPC(numLaplaceIterations);
 568
 569     untangleMeshFV
 570     (
 571         maxNumGlobalIterations,
 572         maxNumIterations,
 573         maxNumSurfaceIterations
 574     );
 575
 576     Info << "Finished smoothing the mesh" << endl;
 577 }
 578
 579 void meshOptimizer::optimizeMeshFVBestQuality
 580 (
 581     const label maxNumIterations,
 582     const scalar threshold
 583 )
 584 {
 585     label nBadFaces, nIter(0);
 586     label minIter(-1);
 587
 588     const faceListPMG& faces = mesh_.faces();
 589     boolList changedFace(faces.size(), true);
 590
 591     //- check if any points in the tet mesh shall not move
 592     labelLongList lockedPoints;
 593     forAll(vertexLocation_, pointI)
 594     {
 595         if( vertexLocation_[pointI] & LOCKED )
 596             lockedPoints.append(pointI);
 597     }
 598
 599     label minNumBadFaces(10 * faces.size());
 600     do
 601     {
 602         labelHashSet lowQualityFaces;
 603         nBadFaces =
 604             polyMeshGenChecks::findWorstQualityFaces
 605             (
 606                 mesh_,
 607                 lowQualityFaces,
 608                 false,
 609                 &changedFace,
 610                 threshold
 611             );
 612
 613         changedFace = false;
 614         forAllConstIter(labelHashSet, lowQualityFaces, it)
 615             changedFace[it.key()] = true;
 616
 617         Info << "Iteration " << nIter
 618             << ". Number of worst quality faces is " << nBadFaces << endl;
 619
 620         //- perform optimisation
 621         if( nBadFaces == 0 )
 622             break;
 623
 624         if( nBadFaces < minNumBadFaces )
 625         {
 626             minNumBadFaces = nBadFaces;
 627
 628             //- update the iteration number when the minimum is achieved
 629             minIter = nIter;
 630         }
 631
 632         partTetMesh tetMesh(mesh_, lockedPoints, lowQualityFaces, 2);
 633
 634         //- construct tetMeshOptimisation and improve positions
 635         //- of points in the tet mesh
 636         tetMeshOptimisation tmo(tetMesh);
 637
 638         tmo.optimiseUsingVolumeOptimizer(20);
 639
 640         //- update points in the mesh from the new coordinates in the tet mesh
 641         tetMesh.updateOrigMesh(&changedFace);
 642
 643     } while( (nIter < minIter+5) && (++nIter < maxNumIterations) );
 644 }
 645
 646 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 647
 648 } // End namespace Foam
 649
 650 // ************************************************************************* //