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