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) 2011 OpenFOAM Foundation
   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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your 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, 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 "Time.H"
  31 #include "polyMesh.H"
  32 #include "twoDPointCorrector.H"
  33 #include "OFstream.H"
  34 #include "multiDirRefinement.H"
  35
  36 #include "triSurface.H"
  37 #include "triSurfaceSearch.H"
  38
  39 #include "cellSet.H"
  40 #include "pointSet.H"
  41 #include "surfaceToCell.H"
  42 #include "surfaceToPoint.H"
  43 #include "cellToPoint.H"
  44 #include "pointToCell.H"
  45 #include "cellToCell.H"
  46 #include "surfaceSets.H"
  47 #include "polyTopoChange.H"
  48 #include "polyTopoChanger.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     forAllConstIter(labelHashSet, cutCells, iter)
 296     {
 297         const label cellI = iter.key();
 298         const labelList& cFaces = mesh.cells()[cellI];
 299
 300         forAll(cFaces, i)
 301         {
 302             const label faceI = cFaces[i];
 303
 304             if (mesh.isInternalFace(faceI))
 305             {
 306                 label nbr = mesh.faceOwner()[faceI];
 307
 308                 if (nbr == cellI)
 309                 {
 310                     nbr = mesh.faceNeighbour()[faceI];
 311                 }
 312
 313                 if (selectInside && inside.found(nbr))
 314                 {
 315                     addCutFaces.insert(nbr);
 316                 }
 317                 else if (selectOutside && outside.found(nbr))
 318                 {
 319                     addCutFaces.insert(nbr);
 320                 }
 321             }
 322         }
 323     }
 324
 325     Info<< "    Selected an additional " << addCutFaces.size()
 326         << " neighbours of cutCells to refine" << endl;
 327
 328     forAllConstIter(labelHashSet, addCutFaces, iter)
 329     {
 330         cutCells.insert(iter.key());
 331     }
 332 }
 333
 334
 335 // Return true if any cells had to be split to keep a difference between
 336 // neighbouring refinement levels < limitDiff.
 337 // Gets cells which will be refined (so increase the refinement level) and
 338 // updates it.
 339 bool limitRefinementLevel
 340 (
 341     const primitiveMesh& mesh,
 342     const label limitDiff,
 343     const labelHashSet& excludeCells,
 344     const labelList& refLevel,
 345     labelHashSet& cutCells
 346 )
 347 {
 348     // Do simple check on validity of refinement level.
 349     forAll(refLevel, cellI)
 350     {
 351         if (!excludeCells.found(cellI))
 352         {
 353             const labelList& cCells = mesh.cellCells()[cellI];
 354
 355             forAll(cCells, i)
 356             {
 357                 label nbr = cCells[i];
 358
 359                 if (!excludeCells.found(nbr))
 360                 {
 361                     if (refLevel[cellI] - refLevel[nbr] >= limitDiff)
 362                     {
 363                         FatalErrorIn("limitRefinementLevel")
 364                             << "Level difference between neighbouring cells "
 365                             << cellI << " and " << nbr
 366                             << " greater than or equal to " << limitDiff << endl
 367                             << "refLevels:" << refLevel[cellI] << ' '
 368                             <<  refLevel[nbr] << abort(FatalError);
 369                     }
 370                 }
 371             }
 372         }
 373     }
 374
 375
 376     labelHashSet addCutCells(cutCells.size());
 377
 378     forAllConstIter(labelHashSet, cutCells, iter)
 379     {
 380         // cellI will be refined.
 381         const label cellI = iter.key();
 382         const labelList& cCells = mesh.cellCells()[cellI];
 383
 384         forAll(cCells, i)
 385         {
 386             const label nbr = cCells[i];
 387
 388             if (!excludeCells.found(nbr) && !cutCells.found(nbr))
 389             {
 390                 if (refLevel[cellI] + 1 - refLevel[nbr] >= limitDiff)
 391                 {
 392                     addCutCells.insert(nbr);
 393                 }
 394             }
 395         }
 396     }
 397
 398     if (addCutCells.size())
 399     {
 400         // Add cells to cutCells.
 401
 402         Info<< "Added an additional " << addCutCells.size() << " cells"
 403             << " to satisfy 1:" << limitDiff << " refinement level"
 404             << endl;
 405
 406         forAllConstIter(labelHashSet, addCutCells, iter)
 407         {
 408             cutCells.insert(iter.key());
 409         }
 410         return true;
 411     }
 412     else
 413     {
 414         Info<< "Added no additional cells"
 415             << " to satisfy 1:" << limitDiff << " refinement level"
 416             << endl;
 417
 418         return false;
 419     }
 420 }
 421
 422
 423 // Do all refinement (i.e. refCells) according to refineDict and update
 424 // refLevel afterwards for added cells
 425 void doRefinement
 426 (
 427     polyMesh& mesh,
 428     const dictionary& refineDict,
 429     const labelHashSet& refCells,
 430     labelList& refLevel
 431 )
 432 {
 433     label oldCells = mesh.nCells();
 434
 435     // Multi-iteration, multi-direction topology modifier.
 436     multiDirRefinement multiRef
 437     (
 438         mesh,
 439         refCells.toc(),
 440         refineDict
 441     );
 442
 443     //
 444     // Update refLevel for split cells
 445     //
 446
 447     refLevel.setSize(mesh.nCells());
 448
 449     for (label cellI = oldCells; cellI < mesh.nCells(); cellI++)
 450     {
 451         refLevel[cellI] = 0;
 452     }
 453
 454     const labelListList& addedCells = multiRef.addedCells();
 455
 456     forAll(addedCells, oldCellI)
 457     {
 458         const labelList& added = addedCells[oldCellI];
 459
 460         if (added.size())
 461         {
 462             // Give all cells resulting from split the refinement level
 463             // of the master.
 464             label masterLevel = ++refLevel[oldCellI];
 465
 466             forAll(added, i)
 467             {
 468                 refLevel[added[i]] = masterLevel;
 469             }
 470         }
 471     }
 472 }
 473
 474
 475 // Subset mesh and update refLevel and cutCells
 476 void subsetMesh
 477 (
 478     polyMesh& mesh,
 479     const label writeMesh,
 480     const label patchI,                 // patchID for exposed faces
 481     const labelHashSet& cellsToRemove,
 482     cellSet& cutCells,
 483     labelIOList& refLevel
 484 )
 485 {
 486     removeCells cellRemover(mesh);
 487
 488     labelList cellLabels(cellsToRemove.toc());
 489
 490     Info<< "Mesh has:" << mesh.nCells() << " cells."
 491         << " Removing:" << cellLabels.size() << " cells" << endl;
 492
 493     // exposed faces
 494     labelList exposedFaces(cellRemover.getExposedFaces(cellLabels));
 495
 496     polyTopoChange meshMod(mesh);
 497     cellRemover.setRefinement
 498     (
 499         cellLabels,
 500         exposedFaces,
 501         labelList(exposedFaces.size(), patchI),
 502         meshMod
 503     );
 504
 505     // Do all changes
 506     Info<< "Morphing ..." << endl;
 507
 508     const Time& runTime = mesh.time();
 509
 510     autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh, false);
 511
 512     if (morphMap().hasMotionPoints())
 513     {
 514         mesh.movePoints(morphMap().preMotionPoints());
 515     }
 516
 517     // Update topology on cellRemover
 518     cellRemover.updateMesh(morphMap());
 519
 520     // Update refLevel for removed cells.
 521     const labelList& cellMap = morphMap().cellMap();
 522
 523     labelList newRefLevel(cellMap.size());
 524
 525     forAll(cellMap, i)
 526     {
 527         newRefLevel[i] = refLevel[cellMap[i]];
 528     }
 529
 530     // Transfer back to refLevel
 531     refLevel.transfer(newRefLevel);
 532
 533     if (writeMesh)
 534     {
 535         Info<< "Writing refined mesh to time " << runTime.timeName() << nl
 536             << endl;
 537
 538         IOstream::defaultPrecision(10);
 539         mesh.write();
 540         refLevel.write();
 541     }
 542
 543     // Update cutCells for removed cells.
 544     cutCells.updateMesh(morphMap());
 545 }
 546
 547
 548 // Divide the cells according to status compared to surface. Constructs sets:
 549 // - cutCells : all cells cut by surface
 550 // - inside   : all cells inside surface
 551 // - outside  :   ,,      outside ,,
 552 // and a combined set:
 553 // - selected : sum of above according to selectCut/Inside/Outside flags.
 554 void classifyCells
 555 (
 556     const polyMesh& mesh,
 557     const fileName& surfName,
 558     const triSurfaceSearch& querySurf,
 559     const pointField& outsidePts,
 560
 561     const bool selectCut,
 562     const bool selectInside,
 563     const bool selectOutside,
 564
 565     const label nCutLayers,
 566
 567     cellSet& inside,
 568     cellSet& outside,
 569     cellSet& cutCells,
 570     cellSet& selected
 571 )
 572 {
 573     // Cut faces with surface and classify cells
 574     surfaceSets::getSurfaceSets
 575     (
 576         mesh,
 577         surfName,
 578         querySurf.surface(),
 579         querySurf,
 580         outsidePts,
 581
 582         nCutLayers,
 583
 584         inside,
 585         outside,
 586         cutCells
 587     );
 588
 589     // Combine wanted parts into selected
 590     if (selectCut)
 591     {
 592         selected.addSet(cutCells);
 593     }
 594     if (selectInside)
 595     {
 596         selected.addSet(inside);
 597     }
 598     if (selectOutside)
 599     {
 600         selected.addSet(outside);
 601     }
 602
 603     Info<< "Determined cell status:" << endl
 604         << "    inside  :" << inside.size() << endl
 605         << "    outside :" << outside.size() << endl
 606         << "    cutCells:" << cutCells.size() << endl
 607         << "    selected:" << selected.size() << endl
 608         << endl;
 609
 610     writeSet(inside, "inside cells");
 611     writeSet(outside, "outside cells");
 612     writeSet(cutCells, "cut cells");
 613     writeSet(selected, "selected cells");
 614 }
 615
 616
 617 // Main program:
 618
 619 int main(int argc, char *argv[])
 620 {
 621     argList::noParallel();
 622
 623 #   include "setRootCase.H"
 624 #   include "createTime.H"
 625 #   include "createPolyMesh.H"
 626
 627     // If nessecary add oldInternalFaces patch
 628     label newPatchI = addPatch(mesh, "oldInternalFaces");
 629
 630
 631     //
 632     // Read motionProperties dictionary
 633     //
 634
 635     Info<< "Checking for motionProperties\n" << endl;
 636
 637     IOobject motionObj
 638     (
 639         "motionProperties",
 640         runTime.constant(),
 641         mesh,
 642         IOobject::MUST_READ_IF_MODIFIED,
 643         IOobject::NO_WRITE
 644     );
 645
 646     // corrector for mesh motion
 647     twoDPointCorrector* correct2DPtr = NULL;
 648
 649     if (motionObj.headerOk())
 650     {
 651         Info<< "Reading " << runTime.constant() / "motionProperties"
 652             << endl << endl;
 653
 654         IOdictionary motionProperties(motionObj);
 655
 656         Switch twoDMotion(motionProperties.lookup("twoDMotion"));
 657
 658         if (twoDMotion)
 659         {
 660             Info<< "Correcting for 2D motion" << endl << endl;
 661             correct2DPtr = new twoDPointCorrector(mesh);
 662         }
 663     }
 664
 665     IOdictionary refineDict
 666     (
 667         IOobject
 668         (
 669             "autoRefineMeshDict",
 670             runTime.system(),
 671             mesh,
 672             IOobject::MUST_READ_IF_MODIFIED,
 673             IOobject::NO_WRITE
 674         )
 675     );
 676
 677     fileName surfName(refineDict.lookup("surface"));
 678     surfName.expand();
 679     label nCutLayers(readLabel(refineDict.lookup("nCutLayers")));
 680     label cellLimit(readLabel(refineDict.lookup("cellLimit")));
 681     bool selectCut(readBool(refineDict.lookup("selectCut")));
 682     bool selectInside(readBool(refineDict.lookup("selectInside")));
 683     bool selectOutside(readBool(refineDict.lookup("selectOutside")));
 684     bool selectHanging(readBool(refineDict.lookup("selectHanging")));
 685
 686     scalar minEdgeLen(readScalar(refineDict.lookup("minEdgeLen")));
 687     scalar maxEdgeLen(readScalar(refineDict.lookup("maxEdgeLen")));
 688     scalar curvature(readScalar(refineDict.lookup("curvature")));
 689     scalar curvDist(readScalar(refineDict.lookup("curvatureDistance")));
 690     pointField outsidePts(refineDict.lookup("outsidePoints"));
 691     label refinementLimit(readLabel(refineDict.lookup("splitLevel")));
 692     bool writeMesh(readBool(refineDict.lookup("writeMesh")));
 693
 694     Info<< "Cells to be used for meshing (0=false, 1=true):" << nl
 695         << "    cells cut by surface            : " << selectCut << nl
 696         << "    cells inside of surface         : " << selectInside << nl
 697         << "    cells outside of surface        : " << selectOutside << nl
 698         << "    hanging cells                   : " << selectHanging << nl
 699         << endl;
 700
 701
 702     if (nCutLayers > 0 && selectInside)
 703     {
 704         WarningIn(args.executable())
 705             << "Illogical settings : Both nCutLayers>0 and selectInside true."
 706             << endl
 707             << "This would mean that inside cells get removed but should be"
 708             << " included in final mesh" << endl;
 709     }
 710
 711     // Surface.
 712     triSurface surf(surfName);
 713
 714     // Search engine on surface
 715     triSurfaceSearch querySurf(surf);
 716
 717     // Search engine on mesh. No face decomposition since mesh unwarped.
 718     meshSearch queryMesh(mesh, false);
 719
 720     // Check all 'outside' points
 721     forAll(outsidePts, outsideI)
 722     {
 723         const point& outsidePoint = outsidePts[outsideI];
 724
 725         if (queryMesh.findCell(outsidePoint, -1, false) == -1)
 726         {
 727             FatalErrorIn(args.executable())
 728                 << "outsidePoint " << outsidePoint
 729                 << " is not inside any cell"
 730                 << exit(FatalError);
 731         }
 732     }
 733
 734
 735
 736     // Current refinement level. Read if present.
 737     labelIOList refLevel
 738     (
 739         IOobject
 740         (
 741             "refinementLevel",
 742             runTime.timeName(),
 743             polyMesh::defaultRegion,
 744             mesh,
 745             IOobject::READ_IF_PRESENT,
 746             IOobject::AUTO_WRITE
 747         ),
 748         labelList(mesh.nCells(), 0)
 749     );
 750
 751     label maxLevel = max(refLevel);
 752
 753     if (maxLevel > 0)
 754     {
 755         Info<< "Read existing refinement level from file "
 756             << refLevel.objectPath() << nl
 757             << "   min level : " << min(refLevel) << nl
 758             << "   max level : " << maxLevel << nl
 759             << endl;
 760     }
 761     else
 762     {
 763         Info<< "Created zero refinement level in file "
 764             << refLevel.objectPath() << nl
 765             << endl;
 766     }
 767
 768
 769
 770
 771     // Print edge stats on original mesh. Leave out 2d-normal direction
 772     direction normalDir(getNormalDir(correct2DPtr));
 773     scalar meshMinEdgeLen = getEdgeStats(mesh, normalDir);
 774
 775     while (meshMinEdgeLen > minEdgeLen)
 776     {
 777         // Get inside/outside/cutCells cellSets.
 778         cellSet inside(mesh, "inside", mesh.nCells()/10);
 779         cellSet outside(mesh, "outside", mesh.nCells()/10);
 780         cellSet cutCells(mesh, "cutCells", mesh.nCells()/10);
 781         cellSet selected(mesh, "selected", mesh.nCells()/10);
 782
 783         classifyCells
 784         (
 785             mesh,
 786             surfName,
 787             querySurf,
 788             outsidePts,
 789
 790             selectCut,      // for determination of selected
 791             selectInside,   // for determination of selected
 792             selectOutside,  // for determination of selected
 793
 794             nCutLayers,     // How many layers of cutCells to include
 795
 796             inside,
 797             outside,
 798             cutCells,
 799             selected        // not used but determined anyway.
 800         );
 801
 802         Info<< "    Selected " << cutCells.name() << " with "
 803             << cutCells.size() << " cells" << endl;
 804
 805         if ((curvDist > 0) && (meshMinEdgeLen < maxEdgeLen))
 806         {
 807             // Done refining enough close to surface. Now switch to more
 808             // intelligent refinement where only cutCells on surface curvature
 809             // are refined.
 810             cellSet curveCells(mesh, "curveCells", mesh.nCells()/10);
 811
 812             selectCurvatureCells
 813             (
 814                 mesh,
 815                 surfName,
 816                 querySurf,
 817                 maxEdgeLen,
 818                 curvature,
 819                 curveCells
 820             );
 821
 822             Info<< "    Selected " << curveCells.name() << " with "
 823                 << curveCells.size() << " cells" << endl;
 824
 825             // Add neighbours to cutCells. This is if selectCut is not
 826             // true and so outside and/or inside cells get exposed there is
 827             // also refinement in them.
 828             if (!selectCut)
 829             {
 830                 addCutNeighbours
 831                 (
 832                     mesh,
 833                     selectInside,
 834                     selectOutside,
 835                     inside,
 836                     outside,
 837                     cutCells
 838                 );
 839             }
 840
 841             // Subset cutCells to only curveCells
 842             cutCells.subset(curveCells);
 843
 844             Info<< "    Removed cells not on surface curvature. Selected "
 845                 << cutCells.size() << endl;
 846         }
 847
 848
 849         if (nCutLayers > 0)
 850         {
 851             // Subset mesh to remove inside cells altogether. Updates cutCells,
 852             // refLevel.
 853             subsetMesh(mesh, writeMesh, newPatchI, inside, cutCells, refLevel);
 854         }
 855
 856
 857         // Added cells from 2:1 refinement level restriction.
 858         while
 859         (
 860             limitRefinementLevel
 861             (
 862                 mesh,
 863                 refinementLimit,
 864                 labelHashSet(),
 865                 refLevel,
 866                 cutCells
 867             )
 868         )
 869         {}
 870
 871
 872         Info<< "    Current cells           : " << mesh.nCells() << nl
 873             << "    Selected for refinement :" << cutCells.size() << nl
 874             << endl;
 875
 876         if (cutCells.empty())
 877         {
 878             Info<< "Stopping refining since 0 cells selected to be refined ..."
 879                 << nl << endl;
 880             break;
 881         }
 882
 883         if ((mesh.nCells() + 8*cutCells.size()) > cellLimit)
 884         {
 885             Info<< "Stopping refining since cell limit reached ..." << nl
 886                 << "Would refine from " << mesh.nCells()
 887                 << " to " << mesh.nCells() + 8*cutCells.size() << " cells."
 888                 << nl << endl;
 889             break;
 890         }
 891
 892         doRefinement
 893         (
 894             mesh,
 895             refineDict,
 896             cutCells,
 897             refLevel
 898         );
 899
 900         Info<< "    After refinement:" << mesh.nCells() << nl
 901             << endl;
 902
 903         if (writeMesh)
 904         {
 905             Info<< "    Writing refined mesh to time " << runTime.timeName()
 906                 << nl << endl;
 907
 908             IOstream::defaultPrecision(10);
 909             mesh.write();
 910             refLevel.write();
 911         }
 912
 913         // Update mesh edge stats.
 914         meshMinEdgeLen = getEdgeStats(mesh, normalDir);
 915     }
 916
 917
 918     if (selectHanging)
 919     {
 920         // Get inside/outside/cutCells cellSets.
 921         cellSet inside(mesh, "inside", mesh.nCells()/10);
 922         cellSet outside(mesh, "outside", mesh.nCells()/10);
 923         cellSet cutCells(mesh, "cutCells", mesh.nCells()/10);
 924         cellSet selected(mesh, "selected", mesh.nCells()/10);
 925
 926         classifyCells
 927         (
 928             mesh,
 929             surfName,
 930             querySurf,
 931             outsidePts,
 932
 933             selectCut,
 934             selectInside,
 935             selectOutside,
 936
 937             nCutLayers,
 938
 939             inside,
 940             outside,
 941             cutCells,
 942             selected
 943         );
 944
 945
 946         // Find any cells which have all their points on the outside of the
 947         // selected set and refine them
 948         labelHashSet hanging = surfaceSets::getHangingCells(mesh, selected);
 949
 950         Info<< "Detected " << hanging.size() << " hanging cells"
 951             << " (cells with all points on"
 952             << " outside of cellSet 'selected').\nThese would probably result"
 953             << " in flattened cells when snapping the mesh to the surface"
 954             << endl;
 955
 956         Info<< "Refining " << hanging.size() << " hanging cells" << nl
 957             << endl;
 958
 959         // Added cells from 2:1 refinement level restriction.
 960         while
 961         (
 962             limitRefinementLevel
 963             (
 964                 mesh,
 965                 refinementLimit,
 966                 labelHashSet(),
 967                 refLevel,
 968                 hanging
 969             )
 970         )
 971         {}
 972
 973         doRefinement(mesh, refineDict, hanging, refLevel);
 974
 975         Info<< "Writing refined mesh to time " << runTime.timeName() << nl
 976             << endl;
 977
 978         // Write final mesh
 979         IOstream::defaultPrecision(10);
 980         mesh.write();
 981         refLevel.write();
 982
 983     }
 984     else if (!writeMesh)
 985     {
 986         Info<< "Writing refined mesh to time " << runTime.timeName() << nl
 987             << endl;
 988
 989         // Write final mesh. (will have been written already if writeMesh=true)
 990         IOstream::defaultPrecision(10);
 991         mesh.write();
 992         refLevel.write();
 993     }
 994
 995
 996     Info<< "End\n" << endl;
 997
 998     return 0;
 999 }
1000
1001
1002 // ************************************************************************* //