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