applications/utilities/mesh/advanced/autoRefineMesh/autoRefineMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25     Utility to refine cells near to a surface.
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "argList.H"
  30 #include "objectRegistry.H"
  31 #include "foamTime.H"
  32 #include "polyMesh.H"
  33 #include "twoDPointCorrector.H"
  34 #include "OFstream.H"
  35 #include "multiDirRefinement.H"
  36
  37 #include "triSurface.H"
  38 #include "triSurfaceSearch.H"
  39
  40 #include "cellSet.H"
  41 #include "pointSet.H"
  42 #include "surfaceToCell.H"
  43 #include "surfaceToPoint.H"
  44 #include "cellToPoint.H"
  45 #include "pointToCell.H"
  46 #include "cellToCell.H"
  47 #include "surfaceSets.H"
  48 #include "directTopoChange.H"
  49 #include "mapPolyMesh.H"
  50 #include "labelIOList.H"
  51 #include "emptyPolyPatch.H"
  52 #include "removeCells.H"
  53 #include "meshSearch.H"
  54
  55 using namespace Foam;
  56
  57 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  58
  59
  60 // Max cos angle for edges to be considered aligned with axis.
  61 static const scalar edgeTol = 1E-3;
  62
  63
  64 void writeSet(const cellSet& cells, const string& msg)
  65 {
  66     Info<< "Writing " << msg << " (" << cells.size() << ") to cellSet "
  67         << cells.instance()/cells.local()/cells.name()
  68         << endl;
  69     cells.write();
  70 }
  71
  72
  73 direction getNormalDir(const twoDPointCorrector* correct2DPtr)
  74 {
  75     direction dir = 3;
  76
  77     if (correct2DPtr)
  78     {
  79         const vector& normal = correct2DPtr->planeNormal();
  80
  81         if (mag(normal & vector(1, 0, 0)) > 1-edgeTol)
  82         {
  83             dir = 0;
  84         }
  85         else if (mag(normal & vector(0, 1, 0)) > 1-edgeTol)
  86         {
  87             dir = 1;
  88         }
  89         else if (mag(normal & vector(0, 0, 1)) > 1-edgeTol)
  90         {
  91             dir = 2;
  92         }
  93     }
  94     return dir;
  95 }
  96
  97
  98
  99 // Calculate some edge statistics on mesh. Return min. edge length over all
 100 // directions but exclude component (0=x, 1=y, 2=z, other=none)
 101 scalar getEdgeStats(const primitiveMesh& mesh, const direction excludeCmpt)
 102 {
 103     label nX = 0;
 104     label nY = 0;
 105     label nZ = 0;
 106
 107     scalar minX = GREAT;
 108     scalar maxX = -GREAT;
 109     vector x(1, 0, 0);
 110
 111     scalar minY = GREAT;
 112     scalar maxY = -GREAT;
 113     vector y(0, 1, 0);
 114
 115     scalar minZ = GREAT;
 116     scalar maxZ = -GREAT;
 117     vector z(0, 0, 1);
 118
 119     scalar minOther = GREAT;
 120     scalar maxOther = -GREAT;
 121
 122     const edgeList& edges = mesh.edges();
 123
 124     forAll(edges, edgeI)
 125     {
 126         const edge& e = edges[edgeI];
 127
 128         vector eVec(e.vec(mesh.points()));
 129
 130         scalar eMag = mag(eVec);
 131
 132         eVec /= eMag;
 133
 134         if (mag(eVec & x) > 1-edgeTol)
 135         {
 136             minX = min(minX, eMag);
 137             maxX = max(maxX, eMag);
 138             nX++;
 139         }
 140         else if (mag(eVec & y) > 1-edgeTol)
 141         {
 142             minY = min(minY, eMag);
 143             maxY = max(maxY, eMag);
 144             nY++;
 145         }
 146         else if (mag(eVec & z) > 1-edgeTol)
 147         {
 148             minZ = min(minZ, eMag);
 149             maxZ = max(maxZ, eMag);
 150             nZ++;
 151         }
 152         else
 153         {
 154             minOther = min(minOther, eMag);
 155             maxOther = max(maxOther, eMag);
 156         }
 157     }
 158
 159     Info<< "Mesh bounding box:" << boundBox(mesh.points()) << nl << nl
 160         << "Mesh edge statistics:" << nl
 161         << "    x aligned :  number:" << nX << "\tminLen:" << minX
 162         << "\tmaxLen:" << maxX << nl
 163         << "    y aligned :  number:" << nY << "\tminLen:" << minY
 164         << "\tmaxLen:" << maxY << nl
 165         << "    z aligned :  number:" << nZ << "\tminLen:" << minZ
 166         << "\tmaxLen:" << maxZ << nl
 167         << "    other     :  number:" << mesh.nEdges() - nX - nY - nZ
 168         << "\tminLen:" << minOther
 169         << "\tmaxLen:" << maxOther << nl << endl;
 170
 171     if (excludeCmpt == 0)
 172     {
 173         return min(minY, min(minZ, minOther));
 174     }
 175     else if (excludeCmpt == 1)
 176     {
 177         return min(minX, min(minZ, minOther));
 178     }
 179     else if (excludeCmpt == 2)
 180     {
 181         return min(minX, min(minY, minOther));
 182     }
 183     else
 184     {
 185         return min(minX, min(minY, min(minZ, minOther)));
 186     }
 187 }
 188
 189
 190 // Adds empty patch if not yet there. Returns patchID.
 191 label addPatch(polyMesh& mesh, const word& patchName)
 192 {
 193     label patchI = mesh.boundaryMesh().findPatchID(patchName);
 194
 195     if (patchI == -1)
 196     {
 197         const polyBoundaryMesh& patches = mesh.boundaryMesh();
 198
 199         List<polyPatch*> newPatches(patches.size() + 1);
 200
 201         // Add empty patch as 0th entry (Note: only since subsetMesh wants this)
 202         patchI = 0;
 203
 204         newPatches[patchI] =
 205             new emptyPolyPatch
 206             (
 207                 Foam::word(patchName),
 208                 0,
 209                 mesh.nInternalFaces(),
 210                 patchI,
 211                 patches
 212             );
 213
 214         forAll(patches, i)
 215         {
 216             const polyPatch& pp = patches[i];
 217
 218             newPatches[i+1] =
 219                 pp.clone
 220                 (
 221                     patches,
 222                     i+1,
 223                     pp.size(),
 224                     pp.start()
 225                 ).ptr();
 226         }
 227
 228         mesh.removeBoundary();
 229         mesh.addPatches(newPatches);
 230
 231         Info<< "Created patch oldInternalFaces at " << patchI << endl;
 232     }
 233     else
 234     {
 235         Info<< "Reusing patch oldInternalFaces at " << patchI << endl;
 236     }
 237
 238
 239     return patchI;
 240 }
 241
 242
 243 // Take surface and select cells based on surface curvature.
 244 void selectCurvatureCells
 245 (
 246     const polyMesh& mesh,
 247     const fileName& surfName,
 248     const triSurfaceSearch& querySurf,
 249     const scalar nearDist,
 250     const scalar curvature,
 251     cellSet& cells
 252 )
 253 {
 254     // Use surfaceToCell to do actual calculation.
 255
 256     // Since we're adding make sure set is on disk.
 257     cells.write();
 258
 259     // Take centre of cell 0 as outside point since info not used.
 260
 261     surfaceToCell cutSource
 262     (
 263         mesh,
 264         surfName,
 265         querySurf.surface(),
 266         querySurf,
 267         pointField(1, mesh.cellCentres()[0]),
 268         false,
 269         false,
 270         false,
 271         nearDist,
 272         curvature
 273     );
 274
 275     cutSource.applyToSet(topoSetSource::ADD, cells);
 276 }
 277
 278
 279 // cutCells contains currently selected cells to be refined. Add neighbours
 280 // on the inside or outside to them.
 281 void addCutNeighbours
 282 (
 283     const polyMesh& mesh,
 284     const bool selectInside,
 285     const bool selectOutside,
 286     const labelHashSet& inside,
 287     const labelHashSet& outside,
 288     labelHashSet& cutCells
 289 )
 290 {
 291     // Pick up face neighbours of cutCells
 292
 293     labelHashSet addCutFaces(cutCells.size());
 294
 295     for
 296     (
 297         labelHashSet::const_iterator iter = cutCells.begin();
 298         iter != cutCells.end();
 299         ++iter
 300     )
 301     {
 302         label cellI = iter.key();
 303
 304         const labelList& cFaces = mesh.cells()[cellI];
 305
 306         forAll(cFaces, i)
 307         {
 308             label faceI = cFaces[i];
 309
 310             if (mesh.isInternalFace(faceI))
 311             {
 312                 label nbr = mesh.faceOwner()[faceI];
 313
 314                 if (nbr == cellI)
 315                 {
 316                     nbr = mesh.faceNeighbour()[faceI];
 317                 }
 318
 319                 if (selectInside && inside.found(nbr))
 320                 {
 321                     addCutFaces.insert(nbr);
 322                 }
 323                 else if (selectOutside && outside.found(nbr))
 324                 {
 325                     addCutFaces.insert(nbr);
 326                 }
 327             }
 328         }
 329     }
 330
 331     Info<< "    Selected an additional " << addCutFaces.size()
 332         << " neighbours of cutCells to refine" << endl;
 333
 334     for
 335     (
 336         labelHashSet::const_iterator iter = addCutFaces.begin();
 337         iter != addCutFaces.end();
 338         ++iter
 339     )
 340     {
 341         cutCells.insert(iter.key());
 342     }
 343 }
 344
 345
 346 // Return true if any cells had to be split to keep a difference between
 347 // neighbouring refinement levels < limitDiff.
 348 // Gets cells which will be refined (so increase the refinement level) and
 349 // updates it.
 350 bool limitRefinementLevel
 351 (
 352     const primitiveMesh& mesh,
 353     const label limitDiff,
 354     const labelHashSet& excludeCells,
 355     const labelList& refLevel,
 356     labelHashSet& cutCells
 357 )
 358 {
 359     // Do simple check on validity of refinement level.
 360     forAll(refLevel, cellI)
 361     {
 362         if (!excludeCells.found(cellI))
 363         {
 364             const labelList& cCells = mesh.cellCells()[cellI];
 365
 366             forAll(cCells, i)
 367             {
 368                 label nbr = cCells[i];
 369
 370                 if (!excludeCells.found(nbr))
 371                 {
 372                     if (refLevel[cellI] - refLevel[nbr] >= limitDiff)
 373                     {
 374                         FatalErrorIn("limitRefinementLevel")
 375                             << "Level difference between neighbouring cells "
 376                             << cellI << " and " << nbr
 377                             << " greater than or equal to " << limitDiff << endl
 378                             << "refLevels:" << refLevel[cellI] << ' '
 379                             <<  refLevel[nbr] << abort(FatalError);
 380                     }
 381                 }
 382             }
 383         }
 384     }
 385
 386
 387     labelHashSet addCutCells(cutCells.size());
 388
 389     for
 390     (
 391         labelHashSet::const_iterator iter = cutCells.begin();
 392         iter != cutCells.end();
 393         ++iter
 394     )
 395     {
 396         // cellI will be refined.
 397         label cellI = iter.key();
 398
 399         const labelList& cCells = mesh.cellCells()[cellI];
 400
 401         forAll(cCells, i)
 402         {
 403             label nbr = cCells[i];
 404
 405             if (!excludeCells.found(nbr) && !cutCells.found(nbr))
 406             {
 407                 if (refLevel[cellI] + 1 - refLevel[nbr] >= limitDiff)
 408                 {
 409                     addCutCells.insert(nbr);
 410                 }
 411             }
 412         }
 413     }
 414
 415     if (addCutCells.size())
 416     {
 417         // Add cells to cutCells.
 418
 419         Info<< "Added an additional " << addCutCells.size() << " cells"
 420             << " to satisfy 1:" << limitDiff << " refinement level"
 421             << endl;
 422         for
 423         (
 424             labelHashSet::const_iterator iter = addCutCells.begin();
 425             iter != addCutCells.end();
 426             ++iter
 427         )
 428         {
 429             cutCells.insert(iter.key());
 430         }
 431         return true;
 432     }
 433     else
 434     {
 435         Info<< "Added no additional cells"
 436             << " to satisfy 1:" << limitDiff << " refinement level"
 437             << endl;
 438
 439         return false;
 440     }
 441 }
 442
 443
 444 // Do all refinement (i.e. refCells) according to refineDict and update
 445 // refLevel afterwards for added cells
 446 void doRefinement
 447 (
 448     polyMesh& mesh,
 449     const dictionary& refineDict,
 450     const labelHashSet& refCells,
 451     labelList& refLevel
 452 )
 453 {
 454     label oldCells = mesh.nCells();
 455
 456     // Multi-iteration, multi-direction topology modifier.
 457     multiDirRefinement multiRef
 458     (
 459         mesh,
 460         refCells.toc(),
 461         refineDict
 462     );
 463
 464     //
 465     // Update refLevel for split cells
 466     //
 467
 468     refLevel.setSize(mesh.nCells());
 469
 470     for (label cellI = oldCells; cellI < mesh.nCells(); cellI++)
 471     {
 472         refLevel[cellI] = 0;
 473     }
 474
 475     const labelListList& addedCells = multiRef.addedCells();
 476
 477     forAll(addedCells, oldCellI)
 478     {
 479         const labelList& added = addedCells[oldCellI];
 480
 481         if (added.size())
 482         {
 483             // Give all cells resulting from split the refinement level
 484             // of the master.
 485             label masterLevel = ++refLevel[oldCellI];
 486
 487             forAll(added, i)
 488             {
 489                 refLevel[added[i]] = masterLevel;
 490             }
 491         }
 492     }
 493 }
 494
 495
 496 // Subset mesh and update refLevel and cutCells
 497 void subsetMesh
 498 (
 499     polyMesh& mesh,
 500     const label writeMesh,
 501     const label patchI,                 // patchID for exposed faces
 502     const labelHashSet& cellsToRemove,
 503     cellSet& cutCells,
 504     labelIOList& refLevel
 505 )
 506 {
 507     removeCells cellRemover(mesh);
 508
 509     labelList cellLabels(cellsToRemove.toc());
 510
 511     Info<< "Mesh has:" << mesh.nCells() << " cells."
 512         << " Removing:" << cellLabels.size() << " cells" << endl;
 513
 514     // exposed faces
 515     labelList exposedFaces(cellRemover.getExposedFaces(cellLabels));
 516
 517     directTopoChange meshMod(mesh);
 518     cellRemover.setRefinement
 519     (
 520         cellLabels,
 521         exposedFaces,
 522         labelList(exposedFaces.size(), patchI),
 523         meshMod
 524     );
 525
 526     // Do all changes
 527     Info<< "Morphing ..." << endl;
 528
 529     const Time& runTime = mesh.time();
 530
 531     autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh, false);
 532
 533     if (morphMap().hasMotionPoints())
 534     {
 535         mesh.movePoints(morphMap().preMotionPoints());
 536     }
 537
 538     // Update topology on cellRemover
 539     cellRemover.updateMesh(morphMap());
 540
 541     // Update refLevel for removed cells.
 542     const labelList& cellMap = morphMap().cellMap();
 543
 544     labelList newRefLevel(cellMap.size());
 545
 546     forAll(cellMap, i)
 547     {
 548         newRefLevel[i] = refLevel[cellMap[i]];
 549     }
 550
 551     // Transfer back to refLevel
 552     refLevel.transfer(newRefLevel);
 553
 554     if (writeMesh)
 555     {
 556         Info<< "Writing refined mesh to time " << runTime.timeName() << nl
 557             << endl;
 558
 559         IOstream::defaultPrecision(10);
 560         mesh.write();
 561         refLevel.write();
 562     }
 563
 564     // Update cutCells for removed cells.
 565     cutCells.updateMesh(morphMap());
 566 }
 567
 568
 569 // Divide the cells according to status compared to surface. Constructs sets:
 570 // - cutCells : all cells cut by surface
 571 // - inside   : all cells inside surface
 572 // - outside  :   ,,      outside ,,
 573 // and a combined set:
 574 // - selected : sum of above according to selectCut/Inside/Outside flags.
 575 void classifyCells
 576 (
 577     const polyMesh& mesh,
 578     const fileName& surfName,
 579     const triSurfaceSearch& querySurf,
 580     const pointField& outsidePts,
 581
 582     const bool selectCut,
 583     const bool selectInside,
 584     const bool selectOutside,
 585
 586     const label nCutLayers,
 587
 588     cellSet& inside,
 589     cellSet& outside,
 590     cellSet& cutCells,
 591     cellSet& selected
 592 )
 593 {
 594     // Cut faces with surface and classify cells
 595     surfaceSets::getSurfaceSets
 596     (
 597         mesh,
 598         surfName,
 599         querySurf.surface(),
 600         querySurf,
 601         outsidePts,
 602
 603         nCutLayers,
 604
 605         inside,
 606         outside,
 607         cutCells
 608     );
 609
 610     // Combine wanted parts into selected
 611     if (selectCut)
 612     {
 613         selected.addSet(cutCells);
 614     }
 615     if (selectInside)
 616     {
 617         selected.addSet(inside);
 618     }
 619     if (selectOutside)
 620     {
 621         selected.addSet(outside);
 622     }
 623
 624     Info<< "Determined cell status:" << endl
 625         << "    inside  :" << inside.size() << endl
 626         << "    outside :" << outside.size() << endl
 627         << "    cutCells:" << cutCells.size() << endl
 628         << "    selected:" << selected.size() << endl
 629         << endl;
 630
 631     writeSet(inside, "inside cells");
 632     writeSet(outside, "outside cells");
 633     writeSet(cutCells, "cut cells");
 634     writeSet(selected, "selected cells");
 635 }
 636
 637
 638 // Main program:
 639
 640 int main(int argc, char *argv[])
 641 {
 642     argList::noParallel();
 643
 644 #   include "setRootCase.H"
 645 #   include "createTime.H"
 646 #   include "createPolyMesh.H"
 647
 648     // If nessecary add oldInternalFaces patch
 649     label newPatchI = addPatch(mesh, "oldInternalFaces");
 650
 651
 652     //
 653     // Read motionProperties dictionary
 654     //
 655
 656     Info<< "Checking for motionProperties\n" << endl;
 657
 658     IOobject motionObj
 659     (
 660         "motionProperties",
 661         runTime.constant(),
 662         mesh,
 663         IOobject::MUST_READ,
 664         IOobject::NO_WRITE
 665     );
 666
 667     // corrector for mesh motion
 668     twoDPointCorrector* correct2DPtr = NULL;
 669
 670     if (motionObj.headerOk())
 671     {
 672         Info<< "Reading " << runTime.constant() / "motionProperties"
 673             << endl << endl;
 674
 675         IOdictionary motionProperties(motionObj);
 676
 677         Switch twoDMotion(motionProperties.lookup("twoDMotion"));
 678
 679         if (twoDMotion)
 680         {
 681             Info<< "Correcting for 2D motion" << endl << endl;
 682             correct2DPtr = new twoDPointCorrector(mesh);
 683         }
 684     }
 685
 686     IOdictionary refineDict
 687     (
 688         IOobject
 689         (
 690             "autoRefineMeshDict",
 691             runTime.system(),
 692             mesh,
 693             IOobject::MUST_READ,
 694             IOobject::NO_WRITE
 695         )
 696     );
 697
 698     fileName surfName(refineDict.lookup("surface"));
 699     surfName.expand();
 700     label nCutLayers(readLabel(refineDict.lookup("nCutLayers")));
 701     label cellLimit(readLabel(refineDict.lookup("cellLimit")));
 702     bool selectCut(readBool(refineDict.lookup("selectCut")));
 703     bool selectInside(readBool(refineDict.lookup("selectInside")));
 704     bool selectOutside(readBool(refineDict.lookup("selectOutside")));
 705     bool selectHanging(readBool(refineDict.lookup("selectHanging")));
 706
 707     scalar minEdgeLen(readScalar(refineDict.lookup("minEdgeLen")));
 708     scalar maxEdgeLen(readScalar(refineDict.lookup("maxEdgeLen")));
 709     scalar curvature(readScalar(refineDict.lookup("curvature")));
 710     scalar curvDist(readScalar(refineDict.lookup("curvatureDistance")));
 711     pointField outsidePts(refineDict.lookup("outsidePoints"));
 712     label refinementLimit(readLabel(refineDict.lookup("splitLevel")));
 713     bool writeMesh(readBool(refineDict.lookup("writeMesh")));
 714
 715     Info<< "Cells to be used for meshing (0=false, 1=true):" << nl
 716         << "    cells cut by surface            : " << selectCut << nl
 717         << "    cells inside of surface         : " << selectInside << nl
 718         << "    cells outside of surface        : " << selectOutside << nl
 719         << "    hanging cells                   : " << selectHanging << nl
 720         << endl;
 721
 722
 723     if (nCutLayers > 0 && selectInside)
 724     {
 725         WarningIn(args.executable())
 726             << "Illogical settings : Both nCutLayers>0 and selectInside true."
 727             << endl
 728             << "This would mean that inside cells get removed but should be"
 729             << " included in final mesh" << endl;
 730     }
 731
 732     // Surface.
 733     triSurface surf(surfName);
 734
 735     // Search engine on surface
 736     triSurfaceSearch querySurf(surf);
 737
 738     // Search engine on mesh. No face decomposition since mesh unwarped.
 739     meshSearch queryMesh(mesh, false);
 740
 741     // Check all 'outside' points
 742     forAll(outsidePts, outsideI)
 743     {
 744         const point& outsidePoint = outsidePts[outsideI];
 745
 746         if (queryMesh.findCell(outsidePoint, -1, false) == -1)
 747         {
 748             FatalErrorIn(args.executable())
 749                 << "outsidePoint " << outsidePoint
 750                 << " is not inside any cell"
 751                 << exit(FatalError);
 752         }
 753     }
 754
 755
 756
 757     // Current refinement level. Read if present.
 758     labelIOList refLevel
 759     (
 760         IOobject
 761         (
 762             "refinementLevel",
 763             runTime.timeName(),
 764             polyMesh::defaultRegion,
 765             mesh,
 766             IOobject::READ_IF_PRESENT,
 767             IOobject::AUTO_WRITE
 768         ),
 769         labelList(mesh.nCells(), 0)
 770     );
 771
 772     label maxLevel = max(refLevel);
 773
 774     if (maxLevel > 0)
 775     {
 776         Info<< "Read existing refinement level from file "
 777             << refLevel.objectPath() << nl
 778             << "   min level : " << min(refLevel) << nl
 779             << "   max level : " << maxLevel << nl
 780             << endl;
 781     }
 782     else
 783     {
 784         Info<< "Created zero refinement level in file "
 785             << refLevel.objectPath() << nl
 786             << endl;
 787     }
 788
 789
 790
 791
 792     // Print edge stats on original mesh. Leave out 2d-normal direction
 793     direction normalDir(getNormalDir(correct2DPtr));
 794     scalar meshMinEdgeLen = getEdgeStats(mesh, normalDir);
 795
 796     while (meshMinEdgeLen > minEdgeLen)
 797     {
 798         // Get inside/outside/cutCells cellSets.
 799         cellSet inside(mesh, "inside", mesh.nCells()/10);
 800         cellSet outside(mesh, "outside", mesh.nCells()/10);
 801         cellSet cutCells(mesh, "cutCells", mesh.nCells()/10);
 802         cellSet selected(mesh, "selected", mesh.nCells()/10);
 803
 804         classifyCells
 805         (
 806             mesh,
 807             surfName,
 808             querySurf,
 809             outsidePts,
 810
 811             selectCut,      // for determination of selected
 812             selectInside,   // for determination of selected
 813             selectOutside,  // for determination of selected
 814
 815             nCutLayers,     // How many layers of cutCells to include
 816
 817             inside,
 818             outside,
 819             cutCells,
 820             selected        // not used but determined anyway.
 821         );
 822
 823         Info<< "    Selected " << cutCells.name() << " with "
 824             << cutCells.size() << " cells" << endl;
 825
 826         if ((curvDist > 0) && (meshMinEdgeLen < maxEdgeLen))
 827         {
 828             // Done refining enough close to surface. Now switch to more
 829             // intelligent refinement where only cutCells on surface curvature
 830             // are refined.
 831             cellSet curveCells(mesh, "curveCells", mesh.nCells()/10);
 832
 833             selectCurvatureCells
 834             (
 835                 mesh,
 836                 surfName,
 837                 querySurf,
 838                 maxEdgeLen,
 839                 curvature,
 840                 curveCells
 841             );
 842
 843             Info<< "    Selected " << curveCells.name() << " with "
 844                 << curveCells.size() << " cells" << endl;
 845
 846             // Add neighbours to cutCells. This is if selectCut is not
 847             // true and so outside and/or inside cells get exposed there is
 848             // also refinement in them.
 849             if (!selectCut)
 850             {
 851                 addCutNeighbours
 852                 (
 853                     mesh,
 854                     selectInside,
 855                     selectOutside,
 856                     inside,
 857                     outside,
 858                     cutCells
 859                 );
 860             }
 861
 862             // Subset cutCells to only curveCells
 863             cutCells.subset(curveCells);
 864
 865             Info<< "    Removed cells not on surface curvature. Selected "
 866                 << cutCells.size() << endl;
 867         }
 868
 869
 870         if (nCutLayers > 0)
 871         {
 872             // Subset mesh to remove inside cells altogether. Updates cutCells,
 873             // refLevel.
 874             subsetMesh(mesh, writeMesh, newPatchI, inside, cutCells, refLevel);
 875         }
 876
 877
 878         // Added cells from 2:1 refinement level restriction.
 879         while
 880         (
 881             limitRefinementLevel
 882             (
 883                 mesh,
 884                 refinementLimit,
 885                 labelHashSet(),
 886                 refLevel,
 887                 cutCells
 888             )
 889         )
 890         {}
 891
 892
 893         Info<< "    Current cells           : " << mesh.nCells() << nl
 894             << "    Selected for refinement :" << cutCells.size() << nl
 895             << endl;
 896
 897         if (cutCells.empty())
 898         {
 899             Info<< "Stopping refining since 0 cells selected to be refined ..."
 900                 << nl << endl;
 901             break;
 902         }
 903
 904         if ((mesh.nCells() + 8*cutCells.size()) > cellLimit)
 905         {
 906             Info<< "Stopping refining since cell limit reached ..." << nl
 907                 << "Would refine from " << mesh.nCells()
 908                 << " to " << mesh.nCells() + 8*cutCells.size() << " cells."
 909                 << nl << endl;
 910             break;
 911         }
 912
 913         doRefinement
 914         (
 915             mesh,
 916             refineDict,
 917             cutCells,
 918             refLevel
 919         );
 920
 921         Info<< "    After refinement:" << mesh.nCells() << nl
 922             << endl;
 923
 924         if (writeMesh)
 925         {
 926             Info<< "    Writing refined mesh to time " << runTime.timeName()
 927                 << nl << endl;
 928
 929             IOstream::defaultPrecision(10);
 930             mesh.write();
 931             refLevel.write();
 932         }
 933
 934         // Update mesh edge stats.
 935         meshMinEdgeLen = getEdgeStats(mesh, normalDir);
 936     }
 937
 938
 939     if (selectHanging)
 940     {
 941         // Get inside/outside/cutCells cellSets.
 942         cellSet inside(mesh, "inside", mesh.nCells()/10);
 943         cellSet outside(mesh, "outside", mesh.nCells()/10);
 944         cellSet cutCells(mesh, "cutCells", mesh.nCells()/10);
 945         cellSet selected(mesh, "selected", mesh.nCells()/10);
 946
 947         classifyCells
 948         (
 949             mesh,
 950             surfName,
 951             querySurf,
 952             outsidePts,
 953
 954             selectCut,
 955             selectInside,
 956             selectOutside,
 957
 958             nCutLayers,
 959
 960             inside,
 961             outside,
 962             cutCells,
 963             selected
 964         );
 965
 966
 967         // Find any cells which have all their points on the outside of the
 968         // selected set and refine them
 969         labelHashSet hanging = surfaceSets::getHangingCells(mesh, selected);
 970
 971         Info<< "Detected " << hanging.size() << " hanging cells"
 972             << " (cells with all points on"
 973             << " outside of cellSet 'selected').\nThese would probably result"
 974             << " in flattened cells when snapping the mesh to the surface"
 975             << endl;
 976
 977         Info<< "Refining " << hanging.size() << " hanging cells" << nl
 978             << endl;
 979
 980         // Added cells from 2:1 refinement level restriction.
 981         while
 982         (
 983             limitRefinementLevel
 984             (
 985                 mesh,
 986                 refinementLimit,
 987                 labelHashSet(),
 988                 refLevel,
 989                 hanging
 990             )
 991         )
 992         {}
 993
 994         doRefinement(mesh, refineDict, hanging, refLevel);
 995
 996         Info<< "Writing refined mesh to time " << runTime.timeName() << nl
 997             << endl;
 998
 999         // Write final mesh
1000         IOstream::defaultPrecision(10);
1001         mesh.write();
1002         refLevel.write();
1003
1004     }
1005     else if (!writeMesh)
1006     {
1007         Info<< "Writing refined mesh to time " << runTime.timeName() << nl
1008             << endl;
1009
1010         // Write final mesh. (will have been written already if writeMesh=true)
1011         IOstream::defaultPrecision(10);
1012         mesh.write();
1013         refLevel.write();
1014     }
1015
1016
1017     Info << "End\n" << endl;
1018
1019     return 0;
1020 }
1021
1022
1023 // ************************************************************************* //