src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.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 "autoRefineDriver.H"
  27 #include "meshRefinement.H"
  28 #include "fvMesh.H"
  29 #include "Time.H"
  30 #include "cellSet.H"
  31 #include "syncTools.H"
  32 #include "refinementParameters.H"
  33 #include "refinementSurfaces.H"
  34 #include "refinementFeatures.H"
  35 #include "shellSurfaces.H"
  36 #include "mapDistributePolyMesh.H"
  37 #include "unitConversion.H"
  38
  39 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  40
  41 namespace Foam
  42 {
  43
  44 defineTypeNameAndDebug(autoRefineDriver, 0);
  45
  46 } // End namespace Foam
  47
  48
  49 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  50
  51
  52 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  53
  54 // Construct from components
  55 Foam::autoRefineDriver::autoRefineDriver
  56 (
  57     meshRefinement& meshRefiner,
  58     decompositionMethod& decomposer,
  59     fvMeshDistribute& distributor,
  60     const labelList& globalToPatch
  61 )
  62 :
  63     meshRefiner_(meshRefiner),
  64     decomposer_(decomposer),
  65     distributor_(distributor),
  66     globalToPatch_(globalToPatch)
  67 {}
  68
  69
  70 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  71
  72 Foam::label Foam::autoRefineDriver::featureEdgeRefine
  73 (
  74     const refinementParameters& refineParams,
  75     const label maxIter,
  76     const label minRefine
  77 )
  78 {
  79     const fvMesh& mesh = meshRefiner_.mesh();
  80
  81     label iter = 0;
  82
  83     if (meshRefiner_.features().size() && maxIter > 0)
  84     {
  85         for (; iter < maxIter; iter++)
  86         {
  87             Info<< nl
  88                 << "Feature refinement iteration " << iter << nl
  89                 << "------------------------------" << nl
  90                 << endl;
  91
  92             labelList candidateCells
  93             (
  94                 meshRefiner_.refineCandidates
  95                 (
  96                     refineParams.keepPoints()[0],    // For now only use one.
  97                     refineParams.curvature(),
  98
  99                     true,               // featureRefinement
 100                     false,              // internalRefinement
 101                     false,              // surfaceRefinement
 102                     false,              // curvatureRefinement
 103                     refineParams.maxGlobalCells(),
 104                     refineParams.maxLocalCells()
 105                 )
 106             );
 107             labelList cellsToRefine
 108             (
 109                 meshRefiner_.meshCutter().consistentRefinement
 110                 (
 111                     candidateCells,
 112                     true
 113                 )
 114             );
 115             Info<< "Determined cells to refine in = "
 116                 << mesh.time().cpuTimeIncrement() << " s" << endl;
 117
 118
 119
 120             label nCellsToRefine = cellsToRefine.size();
 121             reduce(nCellsToRefine, sumOp<label>());
 122
 123             Info<< "Selected for feature refinement : " << nCellsToRefine
 124                 << " cells (out of " << mesh.globalData().nTotalCells()
 125                 << ')' << endl;
 126
 127             if (nCellsToRefine <= minRefine)
 128             {
 129                 Info<< "Stopping refining since too few cells selected."
 130                     << nl << endl;
 131                 break;
 132             }
 133
 134
 135             if (debug > 0)
 136             {
 137                 const_cast<Time&>(mesh.time())++;
 138             }
 139
 140
 141             if
 142             (
 143                 returnReduce
 144                 (
 145                     (mesh.nCells() >= refineParams.maxLocalCells()),
 146                     orOp<bool>()
 147                 )
 148             )
 149             {
 150                 meshRefiner_.balanceAndRefine
 151                 (
 152                     "feature refinement iteration " + name(iter),
 153                     decomposer_,
 154                     distributor_,
 155                     cellsToRefine,
 156                     refineParams.maxLoadUnbalance()
 157                 );
 158             }
 159             else
 160             {
 161                 meshRefiner_.refineAndBalance
 162                 (
 163                     "feature refinement iteration " + name(iter),
 164                     decomposer_,
 165                     distributor_,
 166                     cellsToRefine,
 167                     refineParams.maxLoadUnbalance()
 168                 );
 169             }
 170         }
 171     }
 172     return iter;
 173 }
 174
 175
 176 Foam::label Foam::autoRefineDriver::surfaceOnlyRefine
 177 (
 178     const refinementParameters& refineParams,
 179     const label maxIter
 180 )
 181 {
 182     const fvMesh& mesh = meshRefiner_.mesh();
 183
 184     // Determine the maximum refinement level over all surfaces. This
 185     // determines the minumum number of surface refinement iterations.
 186     label overallMaxLevel = max(meshRefiner_.surfaces().maxLevel());
 187
 188     label iter;
 189     for (iter = 0; iter < maxIter; iter++)
 190     {
 191         Info<< nl
 192             << "Surface refinement iteration " << iter << nl
 193             << "------------------------------" << nl
 194             << endl;
 195
 196
 197         // Determine cells to refine
 198         // ~~~~~~~~~~~~~~~~~~~~~~~~~
 199         // Only look at surface intersections (minLevel and surface curvature),
 200         // do not do internal refinement (refinementShells)
 201
 202         labelList candidateCells
 203         (
 204             meshRefiner_.refineCandidates
 205             (
 206                 refineParams.keepPoints()[0],
 207                 refineParams.curvature(),
 208
 209                 false,              // featureRefinement
 210                 false,              // internalRefinement
 211                 true,               // surfaceRefinement
 212                 true,               // curvatureRefinement
 213                 refineParams.maxGlobalCells(),
 214                 refineParams.maxLocalCells()
 215             )
 216         );
 217         labelList cellsToRefine
 218         (
 219             meshRefiner_.meshCutter().consistentRefinement
 220             (
 221                 candidateCells,
 222                 true
 223             )
 224         );
 225         Info<< "Determined cells to refine in = "
 226             << mesh.time().cpuTimeIncrement() << " s" << endl;
 227
 228
 229         label nCellsToRefine = cellsToRefine.size();
 230         reduce(nCellsToRefine, sumOp<label>());
 231
 232         Info<< "Selected for refinement : " << nCellsToRefine
 233             << " cells (out of " << mesh.globalData().nTotalCells()
 234             << ')' << endl;
 235
 236         // Stop when no cells to refine or have done minimum nessecary
 237         // iterations and not enough cells to refine.
 238         if
 239         (
 240             nCellsToRefine == 0
 241          || (
 242                 iter >= overallMaxLevel
 243              && nCellsToRefine <= refineParams.minRefineCells()
 244             )
 245         )
 246         {
 247             Info<< "Stopping refining since too few cells selected."
 248                 << nl << endl;
 249             break;
 250         }
 251
 252
 253         if (debug)
 254         {
 255             const_cast<Time&>(mesh.time())++;
 256         }
 257
 258
 259         if
 260         (
 261             returnReduce
 262             (
 263                 (mesh.nCells() >= refineParams.maxLocalCells()),
 264                 orOp<bool>()
 265             )
 266         )
 267         {
 268             meshRefiner_.balanceAndRefine
 269             (
 270                 "surface refinement iteration " + name(iter),
 271                 decomposer_,
 272                 distributor_,
 273                 cellsToRefine,
 274                 refineParams.maxLoadUnbalance()
 275             );
 276         }
 277         else
 278         {
 279             meshRefiner_.refineAndBalance
 280             (
 281                 "surface refinement iteration " + name(iter),
 282                 decomposer_,
 283                 distributor_,
 284                 cellsToRefine,
 285                 refineParams.maxLoadUnbalance()
 286             );
 287         }
 288     }
 289     return iter;
 290 }
 291
 292
 293 void Foam::autoRefineDriver::removeInsideCells
 294 (
 295     const refinementParameters& refineParams,
 296     const label nBufferLayers
 297 )
 298 {
 299     Info<< nl
 300         << "Removing mesh beyond surface intersections" << nl
 301         << "------------------------------------------" << nl
 302         << endl;
 303
 304     const fvMesh& mesh = meshRefiner_.mesh();
 305
 306     if (debug)
 307     {
 308        const_cast<Time&>(mesh.time())++;
 309     }
 310
 311     meshRefiner_.splitMesh
 312     (
 313         nBufferLayers,                  // nBufferLayers
 314         globalToPatch_,
 315         refineParams.keepPoints()[0]
 316     );
 317
 318     if (debug)
 319     {
 320         Pout<< "Writing subsetted mesh to time "
 321             << meshRefiner_.timeName() << '.' << endl;
 322         meshRefiner_.write(debug, mesh.time().path()/meshRefiner_.timeName());
 323         Pout<< "Dumped mesh in = "
 324             << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
 325     }
 326 }
 327
 328
 329 Foam::label Foam::autoRefineDriver::shellRefine
 330 (
 331     const refinementParameters& refineParams,
 332     const label maxIter
 333 )
 334 {
 335     const fvMesh& mesh = meshRefiner_.mesh();
 336
 337     // Mark current boundary faces with 0. Have meshRefiner maintain them.
 338     meshRefiner_.userFaceData().setSize(1);
 339
 340     // mark list to remove any refined faces
 341     meshRefiner_.userFaceData()[0].first() = meshRefinement::REMOVE;
 342     meshRefiner_.userFaceData()[0].second() = createWithValues<labelList>
 343     (
 344         mesh.nFaces(),
 345         -1,
 346         meshRefiner_.intersectedFaces(),
 347         0
 348     );
 349
 350     // Determine the maximum refinement level over all volume refinement
 351     // regions. This determines the minumum number of shell refinement
 352     // iterations.
 353     label overallMaxShellLevel = meshRefiner_.shells().maxLevel();
 354
 355     label iter;
 356     for (iter = 0; iter < maxIter; iter++)
 357     {
 358         Info<< nl
 359             << "Shell refinement iteration " << iter << nl
 360             << "----------------------------" << nl
 361             << endl;
 362
 363         labelList candidateCells
 364         (
 365             meshRefiner_.refineCandidates
 366             (
 367                 refineParams.keepPoints()[0],
 368                 refineParams.curvature(),
 369
 370                 false,              // featureRefinement
 371                 true,               // internalRefinement
 372                 false,              // surfaceRefinement
 373                 false,              // curvatureRefinement
 374                 refineParams.maxGlobalCells(),
 375                 refineParams.maxLocalCells()
 376             )
 377         );
 378
 379         if (debug)
 380         {
 381             Pout<< "Dumping " << candidateCells.size()
 382                 << " cells to cellSet candidateCellsFromShells." << endl;
 383
 384             cellSet c(mesh, "candidateCellsFromShells", candidateCells);
 385             c.instance() = meshRefiner_.timeName();
 386             c.write();
 387         }
 388
 389         // Problem choosing starting faces for bufferlayers (bFaces)
 390         //  - we can't use the current intersected boundary faces
 391         //    (intersectedFaces) since this grows indefinitely
 392         //  - if we use 0 faces we don't satisfy bufferLayers from the
 393         //    surface.
 394         //  - possibly we want to have bFaces only the initial set of faces
 395         //    and maintain the list while doing the refinement.
 396         labelList bFaces
 397         (
 398             findIndices(meshRefiner_.userFaceData()[0].second(), 0)
 399         );
 400
 401         //Info<< "Collected boundary faces : "
 402         //    << returnReduce(bFaces.size(), sumOp<label>()) << endl;
 403
 404         labelList cellsToRefine;
 405
 406         if (refineParams.nBufferLayers() <= 2)
 407         {
 408             cellsToRefine = meshRefiner_.meshCutter().consistentSlowRefinement
 409             (
 410                 refineParams.nBufferLayers(),
 411                 candidateCells,                     // cells to refine
 412                 bFaces,                             // faces for nBufferLayers
 413                 1,                                  // point difference
 414                 meshRefiner_.intersectedPoints()    // points to check
 415             );
 416         }
 417         else
 418         {
 419             cellsToRefine = meshRefiner_.meshCutter().consistentSlowRefinement2
 420             (
 421                 refineParams.nBufferLayers(),
 422                 candidateCells,                 // cells to refine
 423                 bFaces                          // faces for nBufferLayers
 424             );
 425         }
 426
 427         Info<< "Determined cells to refine in = "
 428             << mesh.time().cpuTimeIncrement() << " s" << endl;
 429
 430
 431         label nCellsToRefine = cellsToRefine.size();
 432         reduce(nCellsToRefine, sumOp<label>());
 433
 434         Info<< "Selected for internal refinement : " << nCellsToRefine
 435             << " cells (out of " << mesh.globalData().nTotalCells()
 436             << ')' << endl;
 437
 438         // Stop when no cells to refine or have done minimum nessecary
 439         // iterations and not enough cells to refine.
 440         if
 441         (
 442             nCellsToRefine == 0
 443          || (
 444                 iter >= overallMaxShellLevel
 445              && nCellsToRefine <= refineParams.minRefineCells()
 446             )
 447         )
 448         {
 449             Info<< "Stopping refining since too few cells selected."
 450                 << nl << endl;
 451             break;
 452         }
 453
 454
 455         if (debug)
 456         {
 457             const_cast<Time&>(mesh.time())++;
 458         }
 459
 460         if
 461         (
 462             returnReduce
 463             (
 464                 (mesh.nCells() >= refineParams.maxLocalCells()),
 465                 orOp<bool>()
 466             )
 467         )
 468         {
 469             meshRefiner_.balanceAndRefine
 470             (
 471                 "shell refinement iteration " + name(iter),
 472                 decomposer_,
 473                 distributor_,
 474                 cellsToRefine,
 475                 refineParams.maxLoadUnbalance()
 476             );
 477         }
 478         else
 479         {
 480             meshRefiner_.refineAndBalance
 481             (
 482                 "shell refinement iteration " + name(iter),
 483                 decomposer_,
 484                 distributor_,
 485                 cellsToRefine,
 486                 refineParams.maxLoadUnbalance()
 487             );
 488         }
 489     }
 490     meshRefiner_.userFaceData().clear();
 491
 492     return iter;
 493 }
 494
 495
 496 void Foam::autoRefineDriver::baffleAndSplitMesh
 497 (
 498     const refinementParameters& refineParams,
 499     const bool handleSnapProblems,
 500     const dictionary& motionDict
 501 )
 502 {
 503     Info<< nl
 504         << "Splitting mesh at surface intersections" << nl
 505         << "---------------------------------------" << nl
 506         << endl;
 507
 508     const fvMesh& mesh = meshRefiner_.mesh();
 509
 510     // Introduce baffles at surface intersections. Note:
 511     // meshRefiment::surfaceIndex() will
 512     // be like boundary face from now on so not coupled anymore.
 513     meshRefiner_.baffleAndSplitMesh
 514     (
 515         handleSnapProblems,             // detect&remove potential snap problem
 516         false,                          // perpendicular edge connected cells
 517         scalarField(0),                 // per region perpendicular angle
 518         !handleSnapProblems,            // merge free standing baffles?
 519         motionDict,
 520         const_cast<Time&>(mesh.time()),
 521         globalToPatch_,
 522         refineParams.keepPoints()[0]
 523     );
 524 }
 525
 526
 527 void Foam::autoRefineDriver::zonify
 528 (
 529     const refinementParameters& refineParams
 530 )
 531 {
 532     // Mesh is at its finest. Do zoning
 533     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 534     // This puts all faces with intersection across a zoneable surface
 535     // into that surface's faceZone. All cells inside faceZone get given the
 536     // same cellZone.
 537
 538     if (meshRefiner_.surfaces().getNamedSurfaces().size())
 539     {
 540         Info<< nl
 541             << "Introducing zones for interfaces" << nl
 542             << "--------------------------------" << nl
 543             << endl;
 544
 545         const fvMesh& mesh = meshRefiner_.mesh();
 546
 547         if (debug)
 548         {
 549             const_cast<Time&>(mesh.time())++;
 550         }
 551
 552         meshRefiner_.zonify
 553         (
 554             refineParams.keepPoints()[0],
 555             refineParams.allowFreeStandingZoneFaces()
 556         );
 557
 558         if (debug)
 559         {
 560             Pout<< "Writing zoned mesh to time "
 561                 << meshRefiner_.timeName() << '.' << endl;
 562             meshRefiner_.write
 563             (
 564                 debug,
 565                 mesh.time().path()/meshRefiner_.timeName()
 566             );
 567         }
 568
 569         // Check that all faces are synced
 570         meshRefinement::checkCoupledFaceZones(mesh);
 571     }
 572 }
 573
 574
 575 void Foam::autoRefineDriver::splitAndMergeBaffles
 576 (
 577     const refinementParameters& refineParams,
 578     const bool handleSnapProblems,
 579     const dictionary& motionDict
 580 )
 581 {
 582     Info<< nl
 583         << "Handling cells with snap problems" << nl
 584         << "---------------------------------" << nl
 585         << endl;
 586
 587     const fvMesh& mesh = meshRefiner_.mesh();
 588
 589     // Introduce baffles and split mesh
 590     if (debug)
 591     {
 592         const_cast<Time&>(mesh.time())++;
 593     }
 594
 595     const scalarField& perpAngle = meshRefiner_.surfaces().perpendicularAngle();
 596
 597     meshRefiner_.baffleAndSplitMesh
 598     (
 599         handleSnapProblems,
 600         handleSnapProblems,                 // remove perp edge connected cells
 601         perpAngle,                          // perp angle
 602         false,                              // merge free standing baffles?
 603         //true,                               // merge free standing baffles?
 604         motionDict,
 605         const_cast<Time&>(mesh.time()),
 606         globalToPatch_,
 607         refineParams.keepPoints()[0]
 608     );
 609
 610     if (debug)
 611     {
 612         const_cast<Time&>(mesh.time())++;
 613     }
 614
 615     // Duplicate points on baffles that are on more than one cell
 616     // region. This will help snapping pull them to separate surfaces.
 617     meshRefiner_.dupNonManifoldPoints();
 618
 619
 620     // Merge all baffles that are still remaining after duplicating points.
 621     List<labelPair> couples
 622     (
 623         meshRefiner_.getDuplicateFaces   // get all baffles
 624         (
 625             identity(mesh.nFaces()-mesh.nInternalFaces())
 626           + mesh.nInternalFaces()
 627         )
 628     );
 629
 630     label nCouples = returnReduce(couples.size(), sumOp<label>());
 631
 632     Info<< "Detected unsplittable baffles : "
 633         << nCouples << endl;
 634
 635     if (nCouples > 0)
 636     {
 637         // Actually merge baffles. Note: not exactly parallellized. Should
 638         // convert baffle faces into processor faces if they resulted
 639         // from them.
 640         meshRefiner_.mergeBaffles(couples);
 641
 642         if (debug)
 643         {
 644             // Debug:test all is still synced across proc patches
 645             meshRefiner_.checkData();
 646         }
 647         Info<< "Merged free-standing baffles in = "
 648             << mesh.time().cpuTimeIncrement() << " s." << endl;
 649     }
 650
 651     if (debug)
 652     {
 653         Pout<< "Writing handleProblemCells mesh to time "
 654             << meshRefiner_.timeName() << '.' << endl;
 655         meshRefiner_.write(debug, mesh.time().path()/meshRefiner_.timeName());
 656     }
 657 }
 658
 659
 660 void Foam::autoRefineDriver::mergePatchFaces
 661 (
 662     const refinementParameters& refineParams,
 663     const dictionary& motionDict
 664 )
 665 {
 666     const fvMesh& mesh = meshRefiner_.mesh();
 667
 668     Info<< nl
 669         << "Merge refined boundary faces" << nl
 670         << "----------------------------" << nl
 671         << endl;
 672
 673     if (debug)
 674     {
 675         const_cast<Time&>(mesh.time())++;
 676     }
 677
 678     meshRefiner_.mergePatchFacesUndo
 679     (
 680         Foam::cos(degToRad(45.0)),
 681         Foam::cos(degToRad(45.0)),
 682         meshRefiner_.meshedPatches(),
 683         motionDict
 684     );
 685
 686     if (debug)
 687     {
 688         meshRefiner_.checkData();
 689     }
 690
 691     meshRefiner_.mergeEdgesUndo(Foam::cos(degToRad(45.0)), motionDict);
 692
 693     if (debug)
 694     {
 695         meshRefiner_.checkData();
 696     }
 697 }
 698
 699
 700 void Foam::autoRefineDriver::doRefine
 701 (
 702     const dictionary& refineDict,
 703     const refinementParameters& refineParams,
 704     const bool prepareForSnapping,
 705     const dictionary& motionDict
 706 )
 707 {
 708     Info<< nl
 709         << "Refinement phase" << nl
 710         << "----------------" << nl
 711         << endl;
 712
 713     const fvMesh& mesh = meshRefiner_.mesh();
 714
 715     // Check that all the keep points are inside the mesh.
 716     refineParams.findCells(mesh);
 717
 718     // Refine around feature edges
 719     featureEdgeRefine
 720     (
 721         refineParams,
 722         100,    // maxIter
 723         0       // min cells to refine
 724     );
 725
 726     // Refine based on surface
 727     surfaceOnlyRefine
 728     (
 729         refineParams,
 730         100     // maxIter
 731     );
 732
 733     // Remove cells (a certain distance) beyond surface intersections
 734     removeInsideCells
 735     (
 736         refineParams,
 737         1       // nBufferLayers
 738     );
 739
 740     // Internal mesh refinement
 741     shellRefine
 742     (
 743         refineParams,
 744         100    // maxIter
 745     );
 746
 747     // Introduce baffles at surface intersections. Remove sections unreachable
 748     // from keepPoint.
 749     baffleAndSplitMesh(refineParams, prepareForSnapping, motionDict);
 750
 751     // Mesh is at its finest. Do optional zoning.
 752     zonify(refineParams);
 753
 754     // Pull baffles apart
 755     splitAndMergeBaffles(refineParams, prepareForSnapping, motionDict);
 756
 757     // Do something about cells with refined faces on the boundary
 758     if (prepareForSnapping)
 759     {
 760         mergePatchFaces(refineParams, motionDict);
 761     }
 762
 763
 764     if (Pstream::parRun())
 765     {
 766         Info<< nl
 767             << "Doing final balancing" << nl
 768             << "---------------------" << nl
 769             << endl;
 770
 771         //if (debug)
 772         //{
 773         //    const_cast<Time&>(mesh.time())++;
 774         //}
 775
 776         // Do final balancing. Keep zoned faces on one processor since the
 777         // snap phase will convert them to baffles and this only works for
 778         // internal faces.
 779         meshRefiner_.balance
 780         (
 781             true,
 782             false,
 783             scalarField(mesh.nCells(), 1), // dummy weights
 784             decomposer_,
 785             distributor_
 786         );
 787     }
 788 }
 789
 790
 791 // ************************************************************************* //