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