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