src/dynamicMesh/meshCut/meshModifiers/multiDirRefinement/multiDirRefinement.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 \*---------------------------------------------------------------------------*/
  25
  26 #include "multiDirRefinement.H"
  27 #include "polyMesh.H"
  28 #include "polyTopoChanger.H"
  29 #include "Time.H"
  30 #include "undoableMeshCutter.H"
  31 #include "hexCellLooper.H"
  32 #include "geomCellLooper.H"
  33 #include "topoSet.H"
  34 #include "directions.H"
  35 #include "hexRef8.H"
  36 #include "mapPolyMesh.H"
  37 #include "polyTopoChange.H"
  38 #include "ListOps.H"
  39 #include "cellModeller.H"
  40
  41 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  42
  43 defineTypeNameAndDebug(Foam::multiDirRefinement, 0);
  44
  45
  46 // * * * * * * * * * * * * * Private Statc Functions * * * * * * * * * * * * //
  47
  48 // Update refCells pattern for split cells. Note refCells is current
  49 // list of cells to refine (these should all have been refined)
  50 void Foam::multiDirRefinement::addCells
  51 (
  52     const Map<label>& splitMap,
  53     List<refineCell>& refCells
  54 )
  55 {
  56     label newRefI = refCells.size();
  57
  58     label oldSize = refCells.size();
  59
  60     refCells.setSize(newRefI + splitMap.size());
  61
  62     for (label refI = 0; refI < oldSize; refI++)
  63     {
  64         const refineCell& refCell = refCells[refI];
  65
  66         Map<label>::const_iterator iter = splitMap.find(refCell.cellNo());
  67
  68         if (iter == splitMap.end())
  69         {
  70             FatalErrorIn
  71             (
  72                 "multiDirRefinement::addCells(const Map<label>&"
  73                 ", List<refineCell>&)"
  74             )   << "Problem : cannot find added cell for cell "
  75                 << refCell.cellNo() << abort(FatalError);
  76         }
  77
  78         refCells[newRefI++] = refineCell(iter(), refCell.direction());
  79     }
  80 }
  81
  82
  83 // Update vectorField for all the new cells. Takes over value of
  84 // original cell.
  85 void Foam::multiDirRefinement::update
  86 (
  87     const Map<label>& splitMap,
  88     vectorField& field
  89 )
  90 {
  91     field.setSize(field.size() + splitMap.size());
  92
  93     forAllConstIter(Map<label>, splitMap, iter)
  94     {
  95         field[iter()] = field[iter.key()];
  96     }
  97 }
  98
  99
 100 // Add added cells to labelList
 101 void Foam::multiDirRefinement::addCells
 102 (
 103     const Map<label>& splitMap,
 104     labelList& labels
 105 )
 106 {
 107     label newCellI = labels.size();
 108
 109     labels.setSize(labels.size() + splitMap.size());
 110
 111     forAllConstIter(Map<label>, splitMap, iter)
 112     {
 113         labels[newCellI++] = iter();
 114     }
 115 }
 116
 117
 118 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 119
 120 void Foam::multiDirRefinement::addCells
 121 (
 122     const primitiveMesh& mesh,
 123     const Map<label>& splitMap
 124 )
 125 {
 126     // Construct inverse addressing: from new to original cell.
 127     labelList origCell(mesh.nCells(), -1);
 128
 129     forAll(addedCells_, cellI)
 130     {
 131         const labelList& added = addedCells_[cellI];
 132
 133         forAll(added, i)
 134         {
 135             label slave = added[i];
 136
 137             if (origCell[slave] == -1)
 138             {
 139                 origCell[slave] = cellI;
 140             }
 141             else if (origCell[slave] != cellI)
 142             {
 143                 FatalErrorIn
 144                 (
 145                     "multiDirRefinement::addCells(const primitiveMesh&"
 146                     ", const Map<label>&"
 147                 )   << "Added cell " << slave << " has two different masters:"
 148                     << origCell[slave] << " , " << cellI
 149                     << abort(FatalError);
 150             }
 151         }
 152     }
 153
 154
 155     forAllConstIter(Map<label>, splitMap, iter)
 156     {
 157         label masterI = iter.key();
 158         label newCellI = iter();
 159
 160         while (origCell[masterI] != -1 && origCell[masterI] != masterI)
 161         {
 162             masterI = origCell[masterI];
 163         }
 164
 165         if (masterI >= addedCells_.size())
 166         {
 167             FatalErrorIn
 168             (
 169                 "multiDirRefinement::addCells(const primitiveMesh&"
 170                 ", const Map<label>&"
 171             )   << "Map of added cells contains master cell " << masterI
 172                 << " which is not a valid cell number" << endl
 173                 << "This means that the mesh is not consistent with the"
 174                 << " done refinement" << endl
 175                 << "newCell:" << newCellI << abort(FatalError);
 176         }
 177
 178         labelList& added = addedCells_[masterI];
 179
 180         if (added.empty())
 181         {
 182             added.setSize(2);
 183             added[0] = masterI;
 184             added[1] = newCellI;
 185         }
 186         else if (findIndex(added, newCellI) == -1)
 187         {
 188             label sz = added.size();
 189             added.setSize(sz + 1);
 190             added[sz] = newCellI;
 191         }
 192     }
 193 }
 194
 195
 196 Foam::labelList Foam::multiDirRefinement::splitOffHex(const primitiveMesh& mesh)
 197 {
 198     const cellModel& hex = *(cellModeller::lookup("hex"));
 199
 200     const cellShapeList& cellShapes = mesh.cellShapes();
 201
 202     // Split cellLabels_ into two lists.
 203
 204     labelList nonHexLabels(cellLabels_.size());
 205     label nonHexI = 0;
 206
 207     labelList hexLabels(cellLabels_.size());
 208     label hexI = 0;
 209
 210     forAll(cellLabels_, i)
 211     {
 212         label cellI = cellLabels_[i];
 213
 214         if (cellShapes[cellI].model() == hex)
 215         {
 216             hexLabels[hexI++] = cellI;
 217         }
 218         else
 219         {
 220             nonHexLabels[nonHexI++] = cellI;
 221         }
 222     }
 223
 224     nonHexLabels.setSize(nonHexI);
 225
 226     cellLabels_.transfer(nonHexLabels);
 227
 228     hexLabels.setSize(hexI);
 229
 230     return hexLabels;
 231 }
 232
 233
 234 void Foam::multiDirRefinement::refineHex8
 235 (
 236     polyMesh& mesh,
 237     const labelList& hexCells,
 238     const bool writeMesh
 239 )
 240 {
 241     if (debug)
 242     {
 243         Pout<< "multiDirRefinement : Refining hexes " << hexCells.size()
 244             << endl;
 245     }
 246
 247     hexRef8 hexRefiner
 248     (
 249         mesh,
 250         labelList(mesh.nCells(), 0),    // cellLevel
 251         labelList(mesh.nPoints(), 0),   // pointLevel
 252         refinementHistory
 253         (
 254             IOobject
 255             (
 256                 "refinementHistory",
 257                 mesh.facesInstance(),
 258                 polyMesh::meshSubDir,
 259                 mesh,
 260                 IOobject::NO_READ,
 261                 IOobject::NO_WRITE,
 262                 false
 263             ),
 264             List<refinementHistory::splitCell8>(0),
 265             labelList(0)
 266         )                                   // refinement history
 267     );
 268
 269     polyTopoChange meshMod(mesh);
 270
 271     labelList consistentCells
 272     (
 273         hexRefiner.consistentRefinement
 274         (
 275             hexCells,
 276             true                  // buffer layer
 277         )
 278     );
 279
 280     // Check that consistentRefinement hasn't added cells
 281     {
 282         // Create count 1 for original cells
 283         Map<label> hexCellSet(2*hexCells.size());
 284         forAll(hexCells, i)
 285         {
 286             hexCellSet.insert(hexCells[i], 1);
 287         }
 288
 289         // Increment count
 290         forAll(consistentCells, i)
 291         {
 292             const label cellI = consistentCells[i];
 293
 294             Map<label>::iterator iter = hexCellSet.find(cellI);
 295
 296             if (iter == hexCellSet.end())
 297             {
 298                 FatalErrorIn
 299                 (
 300                     "multiDirRefinement::refineHex8"
 301                     "(polyMesh&, const labelList&, const bool)"
 302                 )   << "Resulting mesh would not satisfy 2:1 ratio"
 303                     << " when refining cell " << cellI << abort(FatalError);
 304             }
 305             else
 306             {
 307                 iter() = 2;
 308             }
 309         }
 310
 311         // Check if all been visited (should always be since
 312         // consistentRefinement set up to extend set.
 313         forAllConstIter(Map<label>, hexCellSet, iter)
 314         {
 315             if (iter() != 2)
 316             {
 317                 FatalErrorIn
 318                 (
 319                     "multiDirRefinement::refineHex8"
 320                     "(polyMesh&, const labelList&, const bool)"
 321                 )   << "Resulting mesh would not satisfy 2:1 ratio"
 322                     << " when refining cell " << iter.key()
 323                     << abort(FatalError);
 324             }
 325         }
 326     }
 327
 328
 329     hexRefiner.setRefinement(consistentCells, meshMod);
 330
 331     // Change mesh, no inflation
 332     autoPtr<mapPolyMesh> morphMapPtr = meshMod.changeMesh(mesh, false, true);
 333     const mapPolyMesh& morphMap = morphMapPtr();
 334
 335     if (morphMap.hasMotionPoints())
 336     {
 337         mesh.movePoints(morphMap.preMotionPoints());
 338     }
 339
 340     if (writeMesh)
 341     {
 342         mesh.write();
 343     }
 344
 345     if (debug)
 346     {
 347         Pout<< "multiDirRefinement : updated mesh at time "
 348             << mesh.time().timeName() << endl;
 349     }
 350
 351     hexRefiner.updateMesh(morphMap);
 352
 353     // Collect all cells originating from same old cell (original + 7 extra)
 354
 355     forAll(consistentCells, i)
 356     {
 357         addedCells_[consistentCells[i]].setSize(8);
 358     }
 359     labelList nAddedCells(addedCells_.size(), 0);
 360
 361     const labelList& cellMap = morphMap.cellMap();
 362
 363     forAll(cellMap, cellI)
 364     {
 365         const label oldCellI = cellMap[cellI];
 366
 367         if (addedCells_[oldCellI].size())
 368         {
 369             addedCells_[oldCellI][nAddedCells[oldCellI]++] = cellI;
 370         }
 371     }
 372 }
 373
 374
 375 void Foam::multiDirRefinement::refineAllDirs
 376 (
 377     polyMesh& mesh,
 378     List<vectorField>& cellDirections,
 379     const cellLooper& cellWalker,
 380     undoableMeshCutter& cutter,
 381     const bool writeMesh
 382 )
 383 {
 384     // Iterator
 385     refinementIterator refiner(mesh, cutter, cellWalker, writeMesh);
 386
 387     forAll(cellDirections, dirI)
 388     {
 389         if (debug)
 390         {
 391             Pout<< "multiDirRefinement : Refining " << cellLabels_.size()
 392                 << " cells in direction " << dirI << endl
 393                 << endl;
 394         }
 395
 396         const vectorField& dirField = cellDirections[dirI];
 397
 398         // Combine cell to be cut with direction to cut in.
 399         // If dirField is only one element use this for all cells.
 400
 401         List<refineCell> refCells(cellLabels_.size());
 402
 403         if (dirField.size() == 1)
 404         {
 405             // Uniform directions.
 406             if (debug)
 407             {
 408                 Pout<< "multiDirRefinement : Uniform refinement:"
 409                     << dirField[0] << endl;
 410             }
 411
 412             forAll(refCells, refI)
 413             {
 414                 label cellI = cellLabels_[refI];
 415
 416                 refCells[refI] = refineCell(cellI, dirField[0]);
 417             }
 418         }
 419         else
 420         {
 421             // Non uniform directions.
 422             forAll(refCells, refI)
 423             {
 424                 const label cellI = cellLabels_[refI];
 425
 426                 refCells[refI] = refineCell(cellI, dirField[cellI]);
 427             }
 428         }
 429
 430         // Do refine mesh (multiple iterations). Remember added cells.
 431         Map<label> splitMap = refiner.setRefinement(refCells);
 432
 433         // Update overall map of added cells
 434         addCells(mesh, splitMap);
 435
 436         // Add added cells to list of cells to refine in next iteration
 437         addCells(splitMap, cellLabels_);
 438
 439         // Update refinement direction for added cells.
 440         if (dirField.size() != 1)
 441         {
 442             forAll(cellDirections, i)
 443             {
 444                 update(splitMap, cellDirections[i]);
 445             }
 446         }
 447
 448         if (debug)
 449         {
 450             Pout<< "multiDirRefinement : Done refining direction " << dirI
 451                 << " resulting in " << cellLabels_.size() << " cells" << nl
 452                 << endl;
 453         }
 454     }
 455 }
 456
 457
 458 void Foam::multiDirRefinement::refineFromDict
 459 (
 460     polyMesh& mesh,
 461     List<vectorField>& cellDirections,
 462     const dictionary& dict,
 463     const bool writeMesh
 464 )
 465 {
 466     // How to walk cell circumference.
 467     Switch pureGeomCut(dict.lookup("geometricCut"));
 468
 469     autoPtr<cellLooper> cellWalker(NULL);
 470     if (pureGeomCut)
 471     {
 472         cellWalker.reset(new geomCellLooper(mesh));
 473     }
 474     else
 475     {
 476         cellWalker.reset(new hexCellLooper(mesh));
 477     }
 478
 479     // Construct undoable refinement topology modifier.
 480     //Note: undoability switched off.
 481     // Might want to reconsider if needs to be possible. But then can always
 482     // use other constructor.
 483     undoableMeshCutter cutter(mesh, false);
 484
 485     refineAllDirs(mesh, cellDirections, cellWalker(), cutter, writeMesh);
 486 }
 487
 488
 489
 490 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 491
 492 // Construct from dictionary
 493 Foam::multiDirRefinement::multiDirRefinement
 494 (
 495     polyMesh& mesh,
 496     const labelList& cellLabels,        // list of cells to refine
 497     const dictionary& dict
 498 )
 499 :
 500     cellLabels_(cellLabels),
 501     addedCells_(mesh.nCells())
 502 {
 503     Switch useHex(dict.lookup("useHexTopology"));
 504
 505     Switch writeMesh(dict.lookup("writeMesh"));
 506
 507     wordList dirNames(dict.lookup("directions"));
 508
 509     if (useHex && dirNames.size() == 3)
 510     {
 511         // Filter out hexes from cellLabels_
 512         labelList hexCells(splitOffHex(mesh));
 513
 514         refineHex8(mesh, hexCells, writeMesh);
 515     }
 516
 517     label nRemainingCells = cellLabels_.size();
 518
 519     reduce(nRemainingCells, sumOp<label>());
 520
 521     if (nRemainingCells > 0)
 522     {
 523         // Any cells to refine using meshCutter
 524
 525         // Determine directions for every cell. Can either be uniform
 526         // (size = 1) or non-uniform (one vector per cell)
 527         directions cellDirections(mesh, dict);
 528
 529         refineFromDict(mesh, cellDirections, dict, writeMesh);
 530     }
 531 }
 532
 533
 534 // Construct from directionary and directions to refine.
 535 Foam::multiDirRefinement::multiDirRefinement
 536 (
 537     polyMesh& mesh,
 538     const labelList& cellLabels,        // list of cells to refine
 539     const List<vectorField>& cellDirs,  // Explicitly provided directions
 540     const dictionary& dict
 541 )
 542 :
 543     cellLabels_(cellLabels),
 544     addedCells_(mesh.nCells())
 545 {
 546     Switch useHex(dict.lookup("useHexTopology"));
 547
 548     Switch writeMesh(dict.lookup("writeMesh"));
 549
 550     wordList dirNames(dict.lookup("directions"));
 551
 552     if (useHex && dirNames.size() == 3)
 553     {
 554         // Filter out hexes from cellLabels_
 555         labelList hexCells(splitOffHex(mesh));
 556
 557         refineHex8(mesh, hexCells, writeMesh);
 558     }
 559
 560     label nRemainingCells = cellLabels_.size();
 561
 562     reduce(nRemainingCells, sumOp<label>());
 563
 564     if (nRemainingCells > 0)
 565     {
 566         // Any cells to refine using meshCutter
 567
 568         // Working copy of cells to refine
 569         List<vectorField> cellDirections(cellDirs);
 570
 571         refineFromDict(mesh, cellDirections, dict, writeMesh);
 572     }
 573 }
 574
 575
 576 // Construct from components. Implies meshCutter since directions provided.
 577 Foam::multiDirRefinement::multiDirRefinement
 578 (
 579     polyMesh& mesh,
 580     undoableMeshCutter& cutter,     // actual mesh modifier
 581     const cellLooper& cellWalker,   // how to cut a single cell with
 582                                     // a plane
 583     const labelList& cellLabels,    // list of cells to refine
 584     const List<vectorField>& cellDirs,
 585     const bool writeMesh            // write intermediate meshes
 586 )
 587 :
 588     cellLabels_(cellLabels),
 589     addedCells_(mesh.nCells())
 590 {
 591     // Working copy of cells to refine
 592     List<vectorField> cellDirections(cellDirs);
 593
 594     refineAllDirs(mesh, cellDirections, cellWalker, cutter, writeMesh);
 595 }
 596
 597
 598 // ************************************************************************* //