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