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