src/dynamicMesh/meshCut/meshModifiers/undoableMeshCutter/undoableMeshCutter.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 "undoableMeshCutter.H"
  28 #include "polyMesh.H"
  29 #include "polyTopoChange.H"
  30 #include "DynamicList.H"
  31 #include "meshCutter.H"
  32 #include "cellCuts.H"
  33 #include "splitCell.H"
  34 #include "mapPolyMesh.H"
  35 #include "mathematicalConstants.H"
  36 #include "meshTools.H"
  37
  38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  39
  40 namespace Foam
  41 {
  42
  43 defineTypeNameAndDebug(undoableMeshCutter, 0);
  44
  45 }
  46
  47 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  48
  49 // For debugging
  50 void Foam::undoableMeshCutter::printCellRefTree
  51 (
  52     Ostream& os,
  53     const word& indent,
  54     const splitCell* splitCellPtr
  55 ) const
  56 {
  57     if (splitCellPtr)
  58     {
  59         os << indent << splitCellPtr->cellLabel() << endl;
  60
  61         word subIndent = indent + "--";
  62
  63         printCellRefTree(os, subIndent, splitCellPtr->master());
  64
  65         printCellRefTree(os, subIndent, splitCellPtr->slave());
  66     }
  67 }
  68
  69
  70 // For debugging
  71 void Foam::undoableMeshCutter::printRefTree(Ostream& os) const
  72 {
  73     for
  74     (
  75         Map<splitCell*>::const_iterator iter = liveSplitCells_.begin();
  76         iter != liveSplitCells_.end();
  77         ++iter
  78     )
  79     {
  80         const splitCell* splitPtr = iter();
  81
  82         // Walk to top (master path only)
  83         while (splitPtr->parent())
  84         {
  85             if (!splitPtr->isMaster())
  86             {
  87                 splitPtr = NULL;
  88
  89                 break;
  90             }
  91             else
  92             {
  93                 splitPtr = splitPtr->parent();
  94             }
  95         }
  96
  97         // If we have reached top along master path start printing.
  98         if (splitPtr)
  99         {
 100             // Print from top down
 101             printCellRefTree(os, word(""), splitPtr);
 102         }
 103     }
 104 }
 105
 106
 107 // Update all (cell) labels on splitCell structure.
 108 // Do in two passes to prevent allocation if nothing changed.
 109 void Foam::undoableMeshCutter::updateLabels
 110 (
 111     const labelList& map,
 112     Map<splitCell*>& liveSplitCells
 113 )
 114 {
 115     // Pass1 : check if changed
 116
 117     bool changed = false;
 118
 119     forAllConstIter(Map<splitCell*>,liveSplitCells, iter)
 120     {
 121         const splitCell* splitPtr = iter();
 122
 123         if (!splitPtr)
 124         {
 125             FatalErrorIn
 126             (
 127                 "undoableMeshCutter::updateLabels"
 128                 "(const labelList&, Map<splitCell*>&)"
 129             )   << "Problem: null pointer on liveSplitCells list"
 130                 << abort(FatalError);
 131         }
 132
 133         label cellI = splitPtr->cellLabel();
 134
 135         if (cellI != map[cellI])
 136         {
 137             changed = true;
 138
 139             break;
 140         }
 141     }
 142
 143
 144     // Pass2: relabel
 145
 146     if (changed)
 147     {
 148         // Build new liveSplitCells
 149         // since new labels (= keys in Map) might clash with existing ones.
 150         Map<splitCell*> newLiveSplitCells(2*liveSplitCells.size());
 151
 152         forAllIter(Map<splitCell*>, liveSplitCells, iter)
 153         {
 154             splitCell* splitPtr = iter();
 155
 156             label cellI = splitPtr->cellLabel();
 157
 158             label newCellI = map[cellI];
 159
 160             if (debug && (cellI != newCellI))
 161             {
 162                 Pout<< "undoableMeshCutter::updateLabels :"
 163                     << " Updating live (split)cell from " << cellI
 164                     << " to " << newCellI << endl;
 165             }
 166
 167             if (newCellI >= 0)
 168             {
 169                 // Update splitCell. Can do inplace since only one cellI will
 170                 // refer to this structure.
 171                 splitPtr->cellLabel() = newCellI;
 172
 173                 // Update liveSplitCells
 174                 newLiveSplitCells.insert(newCellI, splitPtr);
 175             }
 176         }
 177         liveSplitCells = newLiveSplitCells;
 178     }
 179 }
 180
 181
 182 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 183
 184 // Construct from components
 185 Foam::undoableMeshCutter::undoableMeshCutter
 186 (
 187     const polyMesh& mesh,
 188     const bool undoable
 189 )
 190 :
 191     meshCutter(mesh),
 192     undoable_(undoable),
 193     liveSplitCells_(mesh.nCells()/100 + 100),
 194     faceRemover_
 195     (
 196         mesh,
 197         Foam::cos(30./180. * mathematicalConstant::pi)
 198     )
 199 {}
 200
 201
 202 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 203
 204 Foam::undoableMeshCutter::~undoableMeshCutter()
 205 {
 206     // Clean split cell tree.
 207
 208     forAllIter(Map<splitCell*>, liveSplitCells_, iter)
 209     {
 210         splitCell* splitPtr = iter();
 211
 212         while (splitPtr)
 213         {
 214             splitCell* parentPtr = splitPtr->parent();
 215
 216             // Sever ties with parent. Also of other side of refinement since
 217             // we are handling rest of tree so other side will not have to.
 218             if (parentPtr)
 219             {
 220                 splitCell* otherSidePtr = splitPtr->getOther();
 221
 222                 otherSidePtr->parent() = NULL;
 223
 224                 splitPtr->parent() = NULL;
 225             }
 226
 227             // Delete splitCell (updates pointer on parent to itself)
 228             delete splitPtr;
 229
 230             splitPtr = parentPtr;
 231         }
 232     }
 233 }
 234
 235
 236 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 237
 238 void Foam::undoableMeshCutter::setRefinement
 239 (
 240     const cellCuts& cuts,
 241     polyTopoChange& meshMod
 242 )
 243 {
 244     // Insert commands to actually cut cells
 245     meshCutter::setRefinement(cuts, meshMod);
 246
 247     if (undoable_)
 248     {
 249         // Use cells cut in this iteration to update splitCell tree.
 250         forAllConstIter(Map<label>, addedCells(), iter)
 251         {
 252             label cellI = iter.key();
 253
 254             label addedCellI = iter();
 255
 256
 257             // Newly created split cell. (cellI ->  cellI + addedCellI)
 258
 259             // Check if cellI already part of split.
 260             Map<splitCell*>::iterator findCell =
 261                 liveSplitCells_.find(cellI);
 262
 263             if (findCell == liveSplitCells_.end())
 264             {
 265                 // CellI not yet split. It cannot be unlive split cell
 266                 // since that would be illegal to split in the first
 267                 // place.
 268
 269                 // Create 0th level. Null parent to denote this.
 270                 splitCell* parentPtr = new splitCell(cellI, NULL);
 271
 272                 splitCell* masterPtr = new splitCell(cellI, parentPtr);
 273
 274                 splitCell* slavePtr = new splitCell(addedCellI, parentPtr);
 275
 276                 // Store newly created cells on parent together with face
 277                 // that splits them
 278                 parentPtr->master() = masterPtr;
 279                 parentPtr->slave() = slavePtr;
 280
 281                 // Insert master and slave into live splitcell list
 282
 283                 if (liveSplitCells_.found(addedCellI))
 284                 {
 285                     FatalErrorIn("undoableMeshCutter::setRefinement")
 286                         << "problem addedCell:" << addedCellI
 287                         << abort(FatalError);
 288                 }
 289
 290                 liveSplitCells_.insert(cellI, masterPtr);
 291                 liveSplitCells_.insert(addedCellI, slavePtr);
 292             }
 293             else
 294             {
 295                 // Cell that was split has been split again.
 296                 splitCell* parentPtr = findCell();
 297
 298                 // It is no longer live
 299                 liveSplitCells_.erase(findCell);
 300
 301                 splitCell* masterPtr = new splitCell(cellI, parentPtr);
 302
 303                 splitCell* slavePtr = new splitCell(addedCellI, parentPtr);
 304
 305                 // Store newly created cells on parent together with face
 306                 // that splits them
 307                 parentPtr->master() = masterPtr;
 308                 parentPtr->slave() = slavePtr;
 309
 310                 // Insert master and slave into live splitcell list
 311
 312                 if (liveSplitCells_.found(addedCellI))
 313                 {
 314                     FatalErrorIn("undoableMeshCutter::setRefinement")
 315                         << "problem addedCell:" << addedCellI
 316                         << abort(FatalError);
 317                 }
 318
 319                 liveSplitCells_.insert(cellI, masterPtr);
 320                 liveSplitCells_.insert(addedCellI, slavePtr);
 321             }
 322         }
 323
 324         if (debug & 2)
 325         {
 326             Pout<< "** After refinement: liveSplitCells_:" << endl;
 327
 328             printRefTree(Pout);
 329         }
 330     }
 331 }
 332
 333
 334 void Foam::undoableMeshCutter::updateMesh(const mapPolyMesh& morphMap)
 335 {
 336     // Update mesh cutter for new labels.
 337     meshCutter::updateMesh(morphMap);
 338
 339     // No need to update cell walker for new labels since does not store any.
 340
 341     // Update faceRemover for new labels
 342     faceRemover_.updateMesh(morphMap);
 343
 344     if (undoable_)
 345     {
 346         // Update all live split cells for mesh mapper.
 347         updateLabels(morphMap.reverseCellMap(), liveSplitCells_);
 348     }
 349 }
 350
 351
 352 Foam::labelList Foam::undoableMeshCutter::getSplitFaces() const
 353 {
 354     if (!undoable_)
 355     {
 356         FatalErrorIn("undoableMeshCutter::getSplitFaces()")
 357             << "Only call if constructed with unrefinement capability"
 358             << abort(FatalError);
 359     }
 360
 361     DynamicList<label> liveSplitFaces(liveSplitCells_.size());
 362
 363     forAllConstIter(Map<splitCell*>, liveSplitCells_, iter)
 364     {
 365         const splitCell* splitPtr = iter();
 366
 367         if (!splitPtr->parent())
 368         {
 369             FatalErrorIn("undoableMeshCutter::getSplitFaces()")
 370                 << "Live split cell without parent" << endl
 371                 << "splitCell:" << splitPtr->cellLabel()
 372                 << abort(FatalError);
 373         }
 374
 375         // Check if not top of refinement and whether it is the master side
 376         if (splitPtr->isMaster())
 377         {
 378             splitCell* slavePtr = splitPtr->getOther();
 379
 380             if
 381             (
 382                 liveSplitCells_.found(slavePtr->cellLabel())
 383              && splitPtr->isUnrefined()
 384              && slavePtr->isUnrefined()
 385             )
 386             {
 387                 // Both master and slave are live and are not refined.
 388                 // Find common face.
 389
 390                 label cellI = splitPtr->cellLabel();
 391
 392                 label slaveCellI = slavePtr->cellLabel();
 393
 394                 label commonFaceI =
 395                     meshTools::getSharedFace
 396                     (
 397                         mesh(),
 398                         cellI,
 399                         slaveCellI
 400                     );
 401
 402                 liveSplitFaces.append(commonFaceI);
 403             }
 404         }
 405     }
 406
 407     return liveSplitFaces.shrink();
 408 }
 409
 410
 411 Foam::Map<Foam::label> Foam::undoableMeshCutter::getAddedCells() const
 412 {
 413     // (code copied from getSplitFaces)
 414
 415     if (!undoable_)
 416     {
 417         FatalErrorIn("undoableMeshCutter::getAddedCells()")
 418             << "Only call if constructed with unrefinement capability"
 419             << abort(FatalError);
 420     }
 421
 422     Map<label> addedCells(liveSplitCells_.size());
 423
 424     forAllConstIter(Map<splitCell*>, liveSplitCells_, iter)
 425     {
 426         const splitCell* splitPtr = iter();
 427
 428         if (!splitPtr->parent())
 429         {
 430             FatalErrorIn("undoableMeshCutter::getAddedCells()")
 431                 << "Live split cell without parent" << endl
 432                 << "splitCell:" << splitPtr->cellLabel()
 433                 << abort(FatalError);
 434         }
 435
 436         // Check if not top of refinement and whether it is the master side
 437         if (splitPtr->isMaster())
 438         {
 439             splitCell* slavePtr = splitPtr->getOther();
 440
 441             if
 442             (
 443                 liveSplitCells_.found(slavePtr->cellLabel())
 444              && splitPtr->isUnrefined()
 445              && slavePtr->isUnrefined()
 446             )
 447             {
 448                 // Both master and slave are live and are not refined.
 449                 addedCells.insert(splitPtr->cellLabel(), slavePtr->cellLabel());
 450             }
 451         }
 452     }
 453     return addedCells;
 454 }
 455
 456
 457 Foam::labelList Foam::undoableMeshCutter::removeSplitFaces
 458 (
 459     const labelList& splitFaces,
 460     polyTopoChange& meshMod
 461 )
 462 {
 463     if (!undoable_)
 464     {
 465         FatalErrorIn("undoableMeshCutter::removeSplitFaces(const labelList&)")
 466             << "Only call if constructed with unrefinement capability"
 467             << abort(FatalError);
 468     }
 469
 470     // Check with faceRemover what faces will get removed. Note that this can
 471     // be more (but never less) than splitFaces provided.
 472     labelList cellRegion;
 473     labelList cellRegionMaster;
 474     labelList facesToRemove;
 475
 476     faceRemover().compatibleRemoves
 477     (
 478         splitFaces,         // pierced faces
 479         cellRegion,         // per cell -1 or region it is merged into
 480         cellRegionMaster,   // per region the master cell
 481         facesToRemove       // new faces to be removed.
 482     );
 483
 484     if (facesToRemove.size() != splitFaces.size())
 485     {
 486         Pout<< "cellRegion:" << cellRegion << endl;
 487         Pout<< "cellRegionMaster:" << cellRegionMaster << endl;
 488
 489         FatalErrorIn
 490         (
 491             "undoableMeshCutter::removeSplitFaces(const labelList&)"
 492         )   << "Faces to remove:" << splitFaces << endl
 493             << "to be removed:" << facesToRemove
 494             << abort(FatalError);
 495     }
 496
 497
 498     // Every face removed will result in neighbour and owner being merged
 499     // into owner.
 500     forAll(facesToRemove, facesToRemoveI)
 501     {
 502         label faceI = facesToRemove[facesToRemoveI];
 503
 504         if (!mesh().isInternalFace(faceI))
 505         {
 506             FatalErrorIn
 507             (
 508                 "undoableMeshCutter::removeSplitFaces(const labelList&)"
 509             )   << "Trying to remove face that is not internal"
 510                 << abort(FatalError);
 511         }
 512
 513         label own = mesh().faceOwner()[faceI];
 514
 515         label nbr = mesh().faceNeighbour()[faceI];
 516
 517         Map<splitCell*>::iterator ownFind = liveSplitCells_.find(own);
 518
 519         Map<splitCell*>::iterator nbrFind = liveSplitCells_.find(nbr);
 520
 521         if
 522         (
 523             (ownFind == liveSplitCells_.end())
 524          || (nbrFind == liveSplitCells_.end())
 525         )
 526         {
 527             // Can happen because of removeFaces adding extra faces to
 528             // original splitFaces
 529         }
 530         else
 531         {
 532             // Face is original splitFace.
 533
 534             splitCell* ownPtr = ownFind();
 535
 536             splitCell* nbrPtr = nbrFind();
 537
 538             splitCell* parentPtr = ownPtr->parent();
 539
 540             // Various checks on sanity.
 541
 542             if (debug)
 543             {
 544                 Pout<< "Updating for removed splitFace " << faceI
 545                     << " own:" << own <<  " nbr:" << nbr
 546                     << " ownPtr:" << ownPtr->cellLabel()
 547                     << " nbrPtr:" << nbrPtr->cellLabel()
 548                     << endl;
 549             }
 550             if (!parentPtr)
 551             {
 552                 FatalErrorIn
 553                 (
 554                     "undoableMeshCutter::removeSplitFaces(const labelList&)"
 555                 )   << "No parent for owner " << ownPtr->cellLabel()
 556                     << abort(FatalError);
 557             }
 558
 559             if (!nbrPtr->parent())
 560             {
 561                 FatalErrorIn
 562                 (
 563                     "undoableMeshCutter::removeSplitFaces(const labelList&)"
 564                 )   << "No parent for neighbour " << nbrPtr->cellLabel()
 565                     << abort(FatalError);
 566             }
 567
 568             if (parentPtr != nbrPtr->parent())
 569             {
 570                 FatalErrorIn
 571                 (
 572                     "undoableMeshCutter::removeSplitFaces(const labelList&)"
 573                 )   << "Owner and neighbour liveSplitCell entries do not have"
 574                     << " same parent. faceI:" << faceI << "  owner:" << own
 575                     << "  ownparent:" << parentPtr->cellLabel()
 576                     << " neighbour:" << nbr
 577                     << "  nbrparent:" << nbrPtr->parent()->cellLabel()
 578                     << abort(FatalError);
 579             }
 580
 581             if
 582             (
 583                 !ownPtr->isUnrefined()
 584              || !nbrPtr->isUnrefined()
 585              || parentPtr->isUnrefined()
 586             )
 587             {
 588                 // Live owner and neighbour are refined themselves.
 589                 FatalErrorIn
 590                 (
 591                     "undoableMeshCutter::removeSplitFaces(const labelList&)"
 592                 )   << "Owner and neighbour liveSplitCell entries are"
 593                     << " refined themselves or the parent is not refined"
 594                     << endl
 595                     << "owner unrefined:" << ownPtr->isUnrefined()
 596                     << "  neighbour unrefined:" << nbrPtr->isUnrefined()
 597                     << "  master unrefined:" << parentPtr->isUnrefined()
 598                     << abort(FatalError);
 599             }
 600
 601             // Delete from liveSplitCell
 602             liveSplitCells_.erase(ownFind);
 603
 604             //!important: Redo search since ownFind entry deleted.
 605             liveSplitCells_.erase(liveSplitCells_.find(nbr));
 606
 607             // Delete entries themselves
 608             delete ownPtr;
 609             delete nbrPtr;
 610
 611             //
 612             // Update parent:
 613             //   - has parent itself: is part of split cell. Update cellLabel
 614             //     with merged cell one.
 615             //   - has no parent: is start of tree. Completely remove.
 616
 617             if (parentPtr->parent())
 618             {
 619                 // Update parent cell label to be new merged cell label
 620                 // (will be owner)
 621                 parentPtr->cellLabel() = own;
 622
 623                 // And insert into live cells (is ok since old entry with
 624                 // own as key has been removed above)
 625                 liveSplitCells_.insert(own, parentPtr);
 626             }
 627             else
 628             {
 629                 // No parent so is start of tree. No need to keep splitCell
 630                 // tree.
 631                 delete parentPtr;
 632             }
 633         }
 634     }
 635
 636     // Insert all commands to combine cells. Never fails so don't have to
 637     // test for success.
 638     faceRemover().setRefinement
 639     (
 640         facesToRemove,
 641         cellRegion,
 642         cellRegionMaster,
 643         meshMod
 644     );
 645
 646     return facesToRemove;
 647 }
 648
 649
 650 // ************************************************************************* //