src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.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 "polyTopoChange.H"
  28 #include "SortableList.H"
  29 #include "polyMesh.H"
  30 #include "polyAddPoint.H"
  31 #include "polyModifyPoint.H"
  32 #include "polyRemovePoint.H"
  33 #include "polyAddFace.H"
  34 #include "polyModifyFace.H"
  35 #include "polyRemoveFace.H"
  36 #include "polyAddCell.H"
  37 #include "polyModifyCell.H"
  38 #include "polyRemoveCell.H"
  39 #include "objectMap.H"
  40 #include "processorPolyPatch.H"
  41 #include "fvMesh.H"
  42
  43 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  44
  45 namespace Foam
  46 {
  47     defineTypeNameAndDebug(polyTopoChange, 0);
  48 }
  49
  50
  51 const Foam::point Foam::polyTopoChange::greatPoint
  52 (
  53     GREAT,
  54     GREAT,
  55     GREAT
  56 );
  57
  58
  59 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  60
  61 // Renumber with special handling for merged items (marked with <-1)
  62 void Foam::polyTopoChange::renumberReverseMap
  63 (
  64     const labelList& map,
  65     DynamicList<label>& elems
  66 )
  67 {
  68     forAll(elems, elemI)
  69     {
  70         label val = elems[elemI];
  71
  72         if (val >= 0)
  73         {
  74             elems[elemI] = map[val];
  75         }
  76         else if (val < -1)
  77         {
  78             label mergedVal = -val-2;
  79             elems[elemI] = -map[mergedVal]-2;
  80         }
  81     }
  82 }
  83
  84
  85 void Foam::polyTopoChange::renumber
  86 (
  87     const labelList& map,
  88     labelHashSet& elems
  89 )
  90 {
  91     labelHashSet newElems(elems.size());
  92
  93     forAllConstIter(labelHashSet, elems, iter)
  94     {
  95         label newElem = map[iter.key()];
  96
  97         if (newElem >= 0)
  98         {
  99             newElems.insert(newElem);
 100         }
 101     }
 102
 103     elems.transfer(newElems);
 104 }
 105
 106
 107 // Renumber and remove -1 elements.
 108 void Foam::polyTopoChange::renumberCompact
 109 (
 110     const labelList& map,
 111     labelList& elems
 112 )
 113 {
 114     label newElemI = 0;
 115
 116     forAll(elems, elemI)
 117     {
 118         label newVal = map[elems[elemI]];
 119
 120         if (newVal != -1)
 121         {
 122             elems[newElemI++] = newVal;
 123         }
 124     }
 125     elems.setSize(newElemI);
 126 }
 127
 128
 129 void Foam::polyTopoChange::countMap
 130 (
 131     const labelList& map,
 132     const labelList& reverseMap,
 133     label& nAdd,
 134     label& nInflate,
 135     label& nMerge,
 136     label& nRemove
 137 )
 138 {
 139     nAdd = 0;
 140     nInflate = 0;
 141     nMerge = 0;
 142     nRemove = 0;
 143
 144     forAll(map, newCellI)
 145     {
 146         label oldCellI = map[newCellI];
 147
 148         if (oldCellI >= 0)
 149         {
 150             if (reverseMap[oldCellI] == newCellI)
 151             {
 152                 // unchanged
 153             }
 154             else
 155             {
 156                 // Added (from another cell v.s. inflated from face/point)
 157                 nAdd++;
 158             }
 159         }
 160         else if (oldCellI == -1)
 161         {
 162             // Created from nothing
 163             nInflate++;
 164         }
 165         else
 166         {
 167             FatalErrorIn("countMap") << "old:" << oldCellI
 168                 << " new:" << newCellI << abort(FatalError);
 169         }
 170     }
 171
 172     forAll(reverseMap, oldCellI)
 173     {
 174         label newCellI = reverseMap[oldCellI];
 175
 176         if (newCellI >= 0)
 177         {
 178             // unchanged
 179         }
 180         else if (newCellI == -1)
 181         {
 182             // removed
 183             nRemove++;
 184         }
 185         else
 186         {
 187             // merged into -newCellI-2
 188             nMerge++;
 189         }
 190     }
 191 }
 192
 193
 194 Foam::labelHashSet Foam::polyTopoChange::getSetIndices
 195 (
 196     const PackedBoolList& lst
 197 )
 198 {
 199     labelHashSet values(lst.count());
 200     forAll(lst, i)
 201     {
 202         if (lst[i])
 203         {
 204             values.insert(i);
 205         }
 206     }
 207     return values;
 208 }
 209
 210
 211 void Foam::polyTopoChange::writeMeshStats(const polyMesh& mesh, Ostream& os)
 212 {
 213     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 214
 215     labelList patchSizes(patches.size());
 216     labelList patchStarts(patches.size());
 217     forAll(patches, patchI)
 218     {
 219         patchSizes[patchI] = patches[patchI].size();
 220         patchStarts[patchI] = patches[patchI].start();
 221     }
 222
 223     os  << "    Points      : " << mesh.nPoints() << nl
 224         << "    Faces       : " << mesh.nFaces() << nl
 225         << "    Cells       : " << mesh.nCells() << nl
 226         << "    PatchSizes  : " << patchSizes << nl
 227         << "    PatchStarts : " << patchStarts << nl
 228         << endl;
 229 }
 230
 231
 232 void Foam::polyTopoChange::getMergeSets
 233 (
 234     const labelList& reverseCellMap,
 235     const labelList& cellMap,
 236     List<objectMap>& cellsFromCells
 237 )
 238 {
 239     // Per new cell the number of old cells that have been merged into it
 240     labelList nMerged(cellMap.size(), 1);
 241
 242     forAll(reverseCellMap, oldCellI)
 243     {
 244         label newCellI = reverseCellMap[oldCellI];
 245
 246         if (newCellI < -1)
 247         {
 248             label mergeCellI = -newCellI-2;
 249
 250             nMerged[mergeCellI]++;
 251         }
 252     }
 253
 254     // From merged cell to set index
 255     labelList cellToMergeSet(cellMap.size(), -1);
 256
 257     label nSets = 0;
 258
 259     forAll(nMerged, cellI)
 260     {
 261         if (nMerged[cellI] > 1)
 262         {
 263             cellToMergeSet[cellI] = nSets++;
 264         }
 265     }
 266
 267     // Collect cell labels.
 268     // Each objectMap will have
 269     // - index : new mesh cell label
 270     // - masterObjects : list of old cells that have been merged. Element 0
 271     //                   will be the original destination cell label.
 272
 273     cellsFromCells.setSize(nSets);
 274
 275     forAll(reverseCellMap, oldCellI)
 276     {
 277         label newCellI = reverseCellMap[oldCellI];
 278
 279         if (newCellI < -1)
 280         {
 281             label mergeCellI = -newCellI-2;
 282
 283             // oldCellI was merged into mergeCellI
 284
 285             label setI = cellToMergeSet[mergeCellI];
 286
 287             objectMap& mergeSet = cellsFromCells[setI];
 288
 289             if (mergeSet.masterObjects().empty())
 290             {
 291                 // First occurrence of master cell mergeCellI
 292
 293                 mergeSet.index() = mergeCellI;
 294                 mergeSet.masterObjects().setSize(nMerged[mergeCellI]);
 295
 296                 // old master label
 297                 mergeSet.masterObjects()[0] = cellMap[mergeCellI];
 298
 299                 // old slave label
 300                 mergeSet.masterObjects()[1] = oldCellI;
 301
 302                 nMerged[mergeCellI] = 2;
 303             }
 304             else
 305             {
 306                 mergeSet.masterObjects()[nMerged[mergeCellI]++] = oldCellI;
 307             }
 308         }
 309     }
 310 }
 311
 312
 313 void Foam::polyTopoChange::checkFace
 314 (
 315     const face& f,
 316     const label faceI,
 317     const label own,
 318     const label nei,
 319     const label patchI,
 320     const label zoneI
 321 ) const
 322 {
 323     if (nei == -1)
 324     {
 325         if (own == -1 && zoneI != -1)
 326         {
 327             // retired face
 328         }
 329         else if (patchI == -1 || patchI >= nPatches_)
 330         {
 331             FatalErrorIn
 332             (
 333                 "polyTopoChange::checkFace(const face&, const label"
 334                 ", const label, const label, const label)"
 335             )   << "Face has no neighbour (so external) but does not have"
 336                 << " a valid patch" << nl
 337                 << "f:" << f
 338                 << " faceI(-1 if added face):" << faceI
 339                 << " own:" << own << " nei:" << nei
 340                 << " patchI:" << patchI << abort(FatalError);
 341         }
 342     }
 343     else
 344     {
 345         if (patchI != -1)
 346         {
 347             FatalErrorIn
 348             (
 349                 "polyTopoChange::checkFace(const face&, const label"
 350                 ", const label, const label, const label)"
 351             )   << "Cannot both have valid patchI and neighbour" << nl
 352                 << "f:" << f
 353                 << " faceI(-1 if added face):" << faceI
 354                 << " own:" << own << " nei:" << nei
 355                 << " patchI:" << patchI << abort(FatalError);
 356         }
 357
 358         if (nei <= own)
 359         {
 360             FatalErrorIn
 361             (
 362                 "polyTopoChange::checkFace(const face&, const label"
 363                 ", const label, const label, const label)"
 364             )   << "Owner cell label should be less than neighbour cell label"
 365                 << nl
 366                 << "f:" << f
 367                 << " faceI(-1 if added face):" << faceI
 368                 << " own:" << own << " nei:" << nei
 369                 << " patchI:" << patchI << abort(FatalError);
 370         }
 371     }
 372
 373     if (f.size() < 3 || findIndex(f, -1) != -1)
 374     {
 375         FatalErrorIn
 376         (
 377             "polyTopoChange::checkFace(const face&, const label"
 378             ", const label, const label, const label)"
 379         )   << "Illegal vertices in face"
 380             << nl
 381             << "f:" << f
 382             << " faceI(-1 if added face):" << faceI
 383             << " own:" << own << " nei:" << nei
 384             << " patchI:" << patchI << abort(FatalError);
 385     }
 386     if (faceI >= 0 && faceI < faces_.size() && faceRemoved(faceI))
 387     {
 388         FatalErrorIn
 389         (
 390             "polyTopoChange::checkFace(const face&, const label"
 391             ", const label, const label, const label)"
 392         )   << "Face already marked for removal"
 393             << nl
 394             << "f:" << f
 395             << " faceI(-1 if added face):" << faceI
 396             << " own:" << own << " nei:" << nei
 397             << " patchI:" << patchI << abort(FatalError);
 398     }
 399     forAll(f, fp)
 400     {
 401         if (f[fp] < points_.size() && pointRemoved(f[fp]))
 402         {
 403             FatalErrorIn
 404             (
 405                 "polyTopoChange::checkFace(const face&, const label"
 406                 ", const label, const label, const label)"
 407             )   << "Face uses removed vertices"
 408                 << nl
 409                 << "f:" << f
 410                 << " faceI(-1 if added face):" << faceI
 411                 << " own:" << own << " nei:" << nei
 412                 << " patchI:" << patchI << abort(FatalError);
 413         }
 414     }
 415 }
 416
 417
 418 void Foam::polyTopoChange::makeCells
 419 (
 420     const label nActiveFaces,
 421     labelList& cellFaces,
 422     labelList& cellFaceOffsets
 423 ) const
 424 {
 425     cellFaces.setSize(2*nActiveFaces);
 426     cellFaceOffsets.setSize(cellMap_.size() + 1);
 427
 428     // Faces per cell
 429     labelList nNbrs(cellMap_.size(), 0);
 430
 431     // 1. Count faces per cell
 432
 433     for (label faceI = 0; faceI < nActiveFaces; faceI++)
 434     {
 435         nNbrs[faceOwner_[faceI]]++;
 436     }
 437     for (label faceI = 0; faceI < nActiveFaces; faceI++)
 438     {
 439         if (faceNeighbour_[faceI] >= 0)
 440         {
 441             nNbrs[faceNeighbour_[faceI]]++;
 442         }
 443     }
 444
 445     // 2. Calculate offsets
 446
 447     cellFaceOffsets[0] = 0;
 448     forAll(nNbrs, cellI)
 449     {
 450         cellFaceOffsets[cellI+1] = cellFaceOffsets[cellI] + nNbrs[cellI];
 451     }
 452
 453     // 3. Fill faces per cell
 454
 455     // reset the whole list to use as counter
 456     nNbrs = 0;
 457
 458     for (label faceI = 0; faceI < nActiveFaces; faceI++)
 459     {
 460         label cellI = faceOwner_[faceI];
 461
 462         cellFaces[cellFaceOffsets[cellI] + nNbrs[cellI]++] = faceI;
 463     }
 464
 465     for (label faceI = 0; faceI < nActiveFaces; faceI++)
 466     {
 467         label cellI = faceNeighbour_[faceI];
 468
 469         if (cellI >= 0)
 470         {
 471             cellFaces[cellFaceOffsets[cellI] + nNbrs[cellI]++] = faceI;
 472         }
 473     }
 474
 475     // Last offset points to beyond end of cellFaces.
 476     cellFaces.setSize(cellFaceOffsets[cellMap_.size()]);
 477 }
 478
 479
 480 // Create cell-cell addressing. Called after compaction (but before ordering)
 481 // of faces
 482 void Foam::polyTopoChange::makeCellCells
 483 (
 484     const label nActiveFaces,
 485     CompactListList<label>& cellCells
 486 ) const
 487 {
 488     // Neighbours per cell
 489     labelList nNbrs(cellMap_.size(), 0);
 490
 491     // Overall number of cellCells
 492     label nCellCells = 0;
 493
 494     // 1. Count neighbours (through internal faces) per cell
 495
 496     for (label faceI = 0; faceI < nActiveFaces; faceI++)
 497     {
 498         if (faceNeighbour_[faceI] >= 0)
 499         {
 500             nNbrs[faceOwner_[faceI]]++;
 501             nNbrs[faceNeighbour_[faceI]]++;
 502             nCellCells += 2;
 503         }
 504     }
 505
 506     cellCells.setSize(cellMap_.size(), nCellCells);
 507
 508     // 2. Calculate offsets
 509
 510     labelList& offsets = cellCells.offsets();
 511
 512     label sumSize = 0;
 513     forAll(nNbrs, cellI)
 514     {
 515         sumSize += nNbrs[cellI];
 516         offsets[cellI] = sumSize;
 517     }
 518
 519     // 3. Fill faces per cell
 520
 521     // reset the whole list to use as counter
 522     nNbrs = 0;
 523
 524     for (label faceI = 0; faceI < nActiveFaces; faceI++)
 525     {
 526         label nei = faceNeighbour_[faceI];
 527
 528         if (nei >= 0)
 529         {
 530             label own = faceOwner_[faceI];
 531             cellCells.m()[cellCells.index(own, nNbrs[own]++)] = nei;
 532             cellCells.m()[cellCells.index(nei, nNbrs[nei]++)] = own;
 533         }
 534     }
 535 }
 536
 537
 538 // Cell ordering (based on bandCompression).
 539 // Handles removed cells. Returns number of remaining cells.
 540 Foam::label Foam::polyTopoChange::getCellOrder
 541 (
 542     const CompactListList<label>& cellCellAddressing,
 543     labelList& oldToNew
 544 ) const
 545 {
 546     const labelList& offsets = cellCellAddressing.offsets();
 547
 548     labelList newOrder(cellCellAddressing.size());
 549
 550     // Fifo buffer for string of cells
 551     SLList<label> nextCell;
 552
 553     // Whether cell has been done already
 554     PackedBoolList visited(cellCellAddressing.size());
 555
 556     label cellInOrder = 0;
 557
 558
 559     // loop over the cells
 560     forAll (visited, cellI)
 561     {
 562         // find the first non-removed cell that has not been visited yet
 563         if (!cellRemoved(cellI) && visited.get(cellI) == 0)
 564         {
 565             // use this cell as a start
 566             nextCell.append(cellI);
 567
 568             // loop through the nextCell list. Add the first cell into the
 569             // cell order if it has not already been visited and ask for its
 570             // neighbours. If the neighbour in question has not been visited,
 571             // add it to the end of the nextCell list
 572
 573             do
 574             {
 575                 label currentCell = nextCell.removeHead();
 576
 577                 if (visited.get(currentCell) == 0)
 578                 {
 579                     visited.set(currentCell, 1);
 580
 581                     // add into cellOrder
 582                     newOrder[cellInOrder] = currentCell;
 583                     cellInOrder++;
 584
 585                     // find if the neighbours have been visited
 586                     label i0 = (currentCell == 0 ? 0 : offsets[currentCell-1]);
 587                     label i1 = offsets[currentCell];
 588
 589                     for (label i = i0; i < i1; i++)
 590                     {
 591                         label nbr = cellCellAddressing.m()[i];
 592
 593                         if (!cellRemoved(nbr) && visited.get(nbr) == 0)
 594                         {
 595                             // not visited, add to the list
 596                             nextCell.append(nbr);
 597                         }
 598                     }
 599                 }
 600             }
 601             while (nextCell.size());
 602         }
 603     }
 604
 605
 606     // Now we have new-to-old in newOrder.
 607     newOrder.setSize(cellInOrder);
 608
 609     // Invert to get old-to-new. Make sure removed (i.e. unmapped) cells are -1.
 610     oldToNew = invert(cellCellAddressing.size(), newOrder);
 611
 612     return cellInOrder;
 613 }
 614
 615
 616 // Determine order for faces:
 617 // - upper-triangular order for internal faces
 618 // - external faces after internal faces and in patch order.
 619 void Foam::polyTopoChange::getFaceOrder
 620 (
 621     const label nActiveFaces,
 622     const labelList& cellFaces,
 623     const labelList& cellFaceOffsets,
 624
 625     labelList& oldToNew,
 626     labelList& patchSizes,
 627     labelList& patchStarts
 628 ) const
 629 {
 630     oldToNew.setSize(faceOwner_.size());
 631     oldToNew = -1;
 632
 633     // First unassigned face
 634     label newFaceI = 0;
 635
 636     forAll(cellMap_, cellI)
 637     {
 638         label startOfCell = cellFaceOffsets[cellI];
 639         label nFaces = cellFaceOffsets[cellI+1] - startOfCell;
 640
 641         // Neighbouring cells
 642         SortableList<label> nbr(nFaces);
 643
 644         for (label i = 0; i < nFaces; i++)
 645         {
 646             label faceI = cellFaces[startOfCell + i];
 647
 648             label nbrCellI = faceNeighbour_[faceI];
 649
 650             if (faceI >= nActiveFaces)
 651             {
 652                 // Retired face.
 653                 nbr[i] = -1;
 654             }
 655             else if (nbrCellI != -1)
 656             {
 657                 // Internal face. Get cell on other side.
 658                 if (nbrCellI == cellI)
 659                 {
 660                     nbrCellI = faceOwner_[faceI];
 661                 }
 662
 663                 if (cellI < nbrCellI)
 664                 {
 665                     // CellI is master
 666                     nbr[i] = nbrCellI;
 667                 }
 668                 else
 669                 {
 670                     // nbrCell is master. Let it handle this face.
 671                     nbr[i] = -1;
 672                 }
 673             }
 674             else
 675             {
 676                 // External face. Do later.
 677                 nbr[i] = -1;
 678             }
 679         }
 680
 681         nbr.sort();
 682
 683         forAll(nbr, i)
 684         {
 685             if (nbr[i] != -1)
 686             {
 687                 oldToNew[cellFaces[startOfCell + nbr.indices()[i]]] =
 688                     newFaceI++;
 689             }
 690         }
 691     }
 692
 693
 694     // Pick up all patch faces in patch face order.
 695     patchStarts.setSize(nPatches_);
 696     patchStarts = 0;
 697     patchSizes.setSize(nPatches_);
 698     patchSizes = 0;
 699
 700     patchStarts[0] = newFaceI;
 701
 702     for (label faceI = 0; faceI < nActiveFaces; faceI++)
 703     {
 704         if (region_[faceI] >= 0)
 705         {
 706             patchSizes[region_[faceI]]++;
 707         }
 708     }
 709
 710     label faceI = patchStarts[0];
 711
 712     forAll(patchStarts, patchI)
 713     {
 714         patchStarts[patchI] = faceI;
 715         faceI += patchSizes[patchI];
 716     }
 717
 718     //if (debug)
 719     //{
 720     //    Pout<< "patchSizes:" << patchSizes << nl
 721     //        << "patchStarts:" << patchStarts << endl;
 722     //}
 723
 724     labelList workPatchStarts(patchStarts);
 725
 726     for (label faceI = 0; faceI < nActiveFaces; faceI++)
 727     {
 728         if (region_[faceI] >= 0)
 729         {
 730             oldToNew[faceI] = workPatchStarts[region_[faceI]]++;
 731         }
 732     }
 733
 734     // Retired faces.
 735     for (label faceI = nActiveFaces; faceI < oldToNew.size(); faceI++)
 736     {
 737         oldToNew[faceI] = faceI;
 738     }
 739
 740     // Check done all faces.
 741     forAll(oldToNew, faceI)
 742     {
 743         if (oldToNew[faceI] == -1)
 744         {
 745             FatalErrorIn
 746             (
 747                 "polyTopoChange::getFaceOrder"
 748                 "(const label, const labelList&, const labelList&)"
 749                 " const"
 750             )   << "Did not determine new position"
 751                 << " for face " << faceI
 752                 << abort(FatalError);
 753         }
 754     }
 755 }
 756
 757
 758 // Reorder and compact faces according to map.
 759 void Foam::polyTopoChange::reorderCompactFaces
 760 (
 761     const label newSize,
 762     const labelList& oldToNew
 763 )
 764 {
 765     reorder(oldToNew, faces_);
 766     faces_.setCapacity(newSize);
 767
 768     reorder(oldToNew, region_);
 769     region_.setCapacity(newSize);
 770
 771     reorder(oldToNew, faceOwner_);
 772     faceOwner_.setCapacity(newSize);
 773
 774     reorder(oldToNew, faceNeighbour_);
 775     faceNeighbour_.setCapacity(newSize);
 776
 777     // Update faceMaps.
 778     reorder(oldToNew, faceMap_);
 779     faceMap_.setCapacity(newSize);
 780
 781     renumberReverseMap(oldToNew, reverseFaceMap_);
 782
 783     renumberKey(oldToNew, faceFromPoint_);
 784     renumberKey(oldToNew, faceFromEdge_);
 785     inplaceReorder(oldToNew, flipFaceFlux_);
 786     flipFaceFlux_.setCapacity(newSize);
 787     renumberKey(oldToNew, faceZone_);
 788     inplaceReorder(oldToNew, faceZoneFlip_);
 789     faceZoneFlip_.setCapacity(newSize);
 790 }
 791
 792
 793 // Compact all and orders points and faces:
 794 // - points into internal followed by external points
 795 // - internalfaces upper-triangular
 796 // - externalfaces after internal ones.
 797 void Foam::polyTopoChange::compact
 798 (
 799     const bool orderCells,
 800     const bool orderPoints,
 801     label& nInternalPoints,
 802     labelList& patchSizes,
 803     labelList& patchStarts
 804 )
 805 {
 806     points_.shrink();
 807     pointMap_.shrink();
 808     reversePointMap_.shrink();
 809
 810     faces_.shrink();
 811     region_.shrink();
 812     faceOwner_.shrink();
 813     faceNeighbour_.shrink();
 814     faceMap_.shrink();
 815     reverseFaceMap_.shrink();
 816
 817     cellMap_.shrink();
 818     reverseCellMap_.shrink();
 819     cellZone_.shrink();
 820
 821
 822     // Compact points
 823     label nActivePoints = 0;
 824     {
 825         labelList localPointMap(points_.size(), -1);
 826         label newPointI = 0;
 827
 828         if (!orderPoints)
 829         {
 830             nInternalPoints = -1;
 831
 832             forAll(points_, pointI)
 833             {
 834                 if (!pointRemoved(pointI) && !retiredPoints_.found(pointI))
 835                 {
 836                     localPointMap[pointI] = newPointI++;
 837                 }
 838             }
 839             nActivePoints = newPointI;
 840         }
 841         else
 842         {
 843             forAll(points_, pointI)
 844             {
 845                 if (!pointRemoved(pointI) && !retiredPoints_.found(pointI))
 846                 {
 847                     nActivePoints++;
 848                 }
 849             }
 850
 851             // Mark boundary points
 852             forAll(faceOwner_, faceI)
 853             {
 854                 if
 855                 (
 856                    !faceRemoved(faceI)
 857                  && faceOwner_[faceI] >= 0
 858                  && faceNeighbour_[faceI] < 0
 859                 )
 860                 {
 861                     // Valid boundary face
 862                     const face& f = faces_[faceI];
 863
 864                     forAll(f, fp)
 865                     {
 866                         label pointI = f[fp];
 867
 868                         if (localPointMap[pointI] == -1)
 869                         {
 870                             if
 871                             (
 872                                 pointRemoved(pointI)
 873                              || retiredPoints_.found(pointI)
 874                             )
 875                             {
 876                                 FatalErrorIn("polyTopoChange::compact(..)")
 877                                     << "Removed or retired point " << pointI
 878                                     << " in face " << f
 879                                     << " at position " << faceI << endl
 880                                     << "Probably face has not been adapted for"
 881                                     << " removed points." << abort(FatalError);
 882                             }
 883                             localPointMap[pointI] = newPointI++;
 884                         }
 885                     }
 886                 }
 887             }
 888
 889             label nBoundaryPoints = newPointI;
 890             nInternalPoints = nActivePoints - nBoundaryPoints;
 891
 892             // Move the boundary addressing up
 893             forAll(localPointMap, pointI)
 894             {
 895                 if (localPointMap[pointI] != -1)
 896                 {
 897                     localPointMap[pointI] += nInternalPoints;
 898                 }
 899             }
 900
 901             newPointI = 0;
 902
 903             // Mark internal points
 904             forAll(faceOwner_, faceI)
 905             {
 906                 if
 907                 (
 908                    !faceRemoved(faceI)
 909                  && faceOwner_[faceI] >= 0
 910                  && faceNeighbour_[faceI] >= 0
 911                 )
 912                 {
 913                     // Valid internal face
 914                     const face& f = faces_[faceI];
 915
 916                     forAll(f, fp)
 917                     {
 918                         label pointI = f[fp];
 919
 920                         if (localPointMap[pointI] == -1)
 921                         {
 922                             if
 923                             (
 924                                 pointRemoved(pointI)
 925                              || retiredPoints_.found(pointI)
 926                             )
 927                             {
 928                                 FatalErrorIn("polyTopoChange::compact(..)")
 929                                     << "Removed or retired point " << pointI
 930                                     << " in face " << f
 931                                     << " at position " << faceI << endl
 932                                     << "Probably face has not been adapted for"
 933                                     << " removed points." << abort(FatalError);
 934                             }
 935                             localPointMap[pointI] = newPointI++;
 936                         }
 937                     }
 938                 }
 939             }
 940
 941             if (newPointI != nInternalPoints)
 942             {
 943                 FatalErrorIn("polyTopoChange::compact(..)")
 944                     << "Problem." << abort(FatalError);
 945             }
 946             newPointI = nActivePoints;
 947         }
 948
 949         forAllConstIter(labelHashSet, retiredPoints_, iter)
 950         {
 951             localPointMap[iter.key()] = newPointI++;
 952         }
 953
 954
 955         if (debug)
 956         {
 957             Pout<< "Points : active:" << nActivePoints
 958                 << "  removed:" << points_.size()-newPointI << endl;
 959         }
 960
 961         reorder(localPointMap, points_);
 962         points_.setCapacity(newPointI);
 963
 964         // Update pointMaps
 965         reorder(localPointMap, pointMap_);
 966         pointMap_.setCapacity(newPointI);
 967         renumberReverseMap(localPointMap, reversePointMap_);
 968
 969         renumberKey(localPointMap, pointZone_);
 970         renumber(localPointMap, retiredPoints_);
 971
 972         // Use map to relabel face vertices
 973         forAll(faces_, faceI)
 974         {
 975             face& f = faces_[faceI];
 976
 977             //labelList oldF(f);
 978             renumberCompact(localPointMap, f);
 979
 980             if (!faceRemoved(faceI) && f.size() < 3)
 981             {
 982                 FatalErrorIn("polyTopoChange::compact(..)")
 983                     << "Created illegal face " << f
 984                     //<< " from face " << oldF
 985                     << " at position:" << faceI
 986                     << " when filtering removed points"
 987                     << abort(FatalError);
 988             }
 989         }
 990     }
 991
 992
 993     // Compact faces.
 994     {
 995         labelList localFaceMap(faces_.size(), -1);
 996         label newFaceI = 0;
 997
 998         forAll(faces_, faceI)
 999         {
1000             if (!faceRemoved(faceI) && faceOwner_[faceI] >= 0)
1001             {
1002                 localFaceMap[faceI] = newFaceI++;
1003             }
1004         }
1005         nActiveFaces_ = newFaceI;
1006
1007         forAll(faces_, faceI)
1008         {
1009             if (!faceRemoved(faceI) && faceOwner_[faceI] < 0)
1010             {
1011                 // Retired face
1012                 localFaceMap[faceI] = newFaceI++;
1013             }
1014         }
1015
1016         if (debug)
1017         {
1018             Pout<< "Faces : active:" << nActiveFaces_
1019                 << "  removed:" << faces_.size()-newFaceI << endl;
1020         }
1021
1022         // Reorder faces.
1023         reorderCompactFaces(newFaceI, localFaceMap);
1024     }
1025
1026     // Compact cells.
1027     {
1028         labelList localCellMap;
1029         label newCellI;
1030
1031         if (orderCells)
1032         {
1033             // Construct cellCell addressing
1034             CompactListList<label> cellCells;
1035             makeCellCells(nActiveFaces_, cellCells);
1036
1037             // Cell ordering (based on bandCompression). Handles removed cells.
1038             newCellI = getCellOrder(cellCells, localCellMap);
1039         }
1040         else
1041         {
1042             // Compact out removed cells
1043             localCellMap.setSize(cellMap_.size());
1044             localCellMap = -1;
1045
1046             newCellI = 0;
1047             forAll(cellMap_, cellI)
1048             {
1049                 if (!cellRemoved(cellI))
1050                 {
1051                     localCellMap[cellI] = newCellI++;
1052                 }
1053             }
1054         }
1055
1056         if (debug)
1057         {
1058             Pout<< "Cells : active:" << newCellI
1059                 << "  removed:" << cellMap_.size()-newCellI << endl;
1060         }
1061
1062         // Renumber -if cells reordered or -if cells removed
1063         if (orderCells || (newCellI != cellMap_.size()))
1064         {
1065             reorder(localCellMap, cellMap_);
1066             cellMap_.setCapacity(newCellI);
1067             renumberReverseMap(localCellMap, reverseCellMap_);
1068
1069             reorder(localCellMap, cellZone_);
1070             cellZone_.setCapacity(newCellI);
1071
1072             renumberKey(localCellMap, cellFromPoint_);
1073             renumberKey(localCellMap, cellFromEdge_);
1074             renumberKey(localCellMap, cellFromFace_);
1075
1076             // Renumber owner/neighbour. Take into account if neighbour suddenly
1077             // gets lower cell than owner.
1078             forAll(faceOwner_, faceI)
1079             {
1080                 label own = faceOwner_[faceI];
1081                 label nei = faceNeighbour_[faceI];
1082
1083                 if (own >= 0)
1084                 {
1085                     // Update owner
1086                     faceOwner_[faceI] = localCellMap[own];
1087
1088                     if (nei >= 0)
1089                     {
1090                         // Update neighbour.
1091                         faceNeighbour_[faceI] = localCellMap[nei];
1092
1093                         // Check if face needs reversing.
1094                         if
1095                         (
1096                             faceNeighbour_[faceI] >= 0
1097                          && faceNeighbour_[faceI] < faceOwner_[faceI]
1098                         )
1099                         {
1100                             faces_[faceI] = faces_[faceI].reverseFace();
1101                             Swap(faceOwner_[faceI], faceNeighbour_[faceI]);
1102                             flipFaceFlux_[faceI] =
1103                             (
1104                                 flipFaceFlux_[faceI]
1105                               ? 0
1106                               : 1
1107                             );
1108                             faceZoneFlip_[faceI] =
1109                             (
1110                                 faceZoneFlip_[faceI]
1111                               ? 0
1112                               : 1
1113                             );
1114                         }
1115                     }
1116                 }
1117                 else if (nei >= 0)
1118                 {
1119                     // Update neighbour.
1120                     faceNeighbour_[faceI] = localCellMap[nei];
1121                 }
1122             }
1123         }
1124     }
1125
1126     // Reorder faces into upper-triangular and patch ordering
1127     {
1128         // Create cells (packed storage)
1129         labelList cellFaces;
1130         labelList cellFaceOffsets;
1131         makeCells(nActiveFaces_, cellFaces, cellFaceOffsets);
1132
1133         // Do upper triangular order and patch sorting
1134         labelList localFaceMap;
1135         getFaceOrder
1136         (
1137             nActiveFaces_,
1138             cellFaces,
1139             cellFaceOffsets,
1140
1141             localFaceMap,
1142             patchSizes,
1143             patchStarts
1144         );
1145
1146         // Reorder faces.
1147         reorderCompactFaces(localFaceMap.size(), localFaceMap);
1148     }
1149 }
1150
1151
1152 // Find faces to interpolate to create value for new face. Only used if
1153 // face was inflated from edge or point. Internal faces should only be
1154 // created from internal faces, external faces only from external faces
1155 // (and ideally the same patch)
1156 // Is bit problematic if there are no faces to select, i.e. in polyDualMesh
1157 // an internal face can be created from a boundary edge with no internal
1158 // faces connected to it.
1159 Foam::labelList Foam::polyTopoChange::selectFaces
1160 (
1161     const primitiveMesh& mesh,
1162     const labelList& faceLabels,
1163     const bool internalFacesOnly
1164 )
1165 {
1166     label nFaces = 0;
1167
1168     forAll(faceLabels, i)
1169     {
1170         label faceI = faceLabels[i];
1171
1172         if (internalFacesOnly == mesh.isInternalFace(faceI))
1173         {
1174             nFaces++;
1175         }
1176     }
1177
1178     labelList collectedFaces;
1179
1180     if (nFaces == 0)
1181     {
1182         // Did not find any faces of the correct type so just use any old
1183         // face.
1184         collectedFaces = faceLabels;
1185     }
1186     else
1187     {
1188         collectedFaces.setSize(nFaces);
1189
1190         nFaces = 0;
1191
1192         forAll(faceLabels, i)
1193         {
1194             label faceI = faceLabels[i];
1195
1196             if (internalFacesOnly == mesh.isInternalFace(faceI))
1197             {
1198                 collectedFaces[nFaces++] = faceI;
1199             }
1200         }
1201     }
1202
1203     return collectedFaces;
1204 }
1205
1206
1207 // Calculate pointMap per patch (so from patch point label to old patch point
1208 // label)
1209 void Foam::polyTopoChange::calcPatchPointMap
1210 (
1211     const List<Map<label> >& oldPatchMeshPointMaps,
1212     const polyBoundaryMesh& boundary,
1213     labelListList& patchPointMap
1214 ) const
1215 {
1216     patchPointMap.setSize(boundary.size());
1217
1218     forAll(boundary, patchI)
1219     {
1220         const labelList& meshPoints = boundary[patchI].meshPoints();
1221
1222         const Map<label>& oldMeshPointMap = oldPatchMeshPointMaps[patchI];
1223
1224         labelList& curPatchPointRnb = patchPointMap[patchI];
1225
1226         curPatchPointRnb.setSize(meshPoints.size());
1227
1228         forAll(meshPoints, i)
1229         {
1230             if (meshPoints[i] < pointMap_.size())
1231             {
1232                 // Check if old point was part of same patch
1233                 Map<label>::const_iterator ozmpmIter = oldMeshPointMap.find
1234                 (
1235                     pointMap_[meshPoints[i]]
1236                 );
1237
1238                 if (ozmpmIter != oldMeshPointMap.end())
1239                 {
1240                     curPatchPointRnb[i] = ozmpmIter();
1241                 }
1242                 else
1243                 {
1244                     curPatchPointRnb[i] = -1;
1245                 }
1246             }
1247             else
1248             {
1249                 curPatchPointRnb[i] = -1;
1250             }
1251         }
1252     }
1253 }
1254
1255
1256 void Foam::polyTopoChange::calcFaceInflationMaps
1257 (
1258     const polyMesh& mesh,
1259     List<objectMap>& facesFromPoints,
1260     List<objectMap>& facesFromEdges,
1261     List<objectMap>& facesFromFaces
1262 ) const
1263 {
1264     // Faces inflated from points
1265     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
1266
1267     facesFromPoints.setSize(faceFromPoint_.size());
1268
1269     if (faceFromPoint_.size())
1270     {
1271         label nFacesFromPoints = 0;
1272
1273         // Collect all still existing faces connected to this point.
1274         forAllConstIter(Map<label>, faceFromPoint_, iter)
1275         {
1276             label newFaceI = iter.key();
1277
1278             if (region_[newFaceI] == -1)
1279             {
1280                 // Get internal faces using point on old mesh
1281                 facesFromPoints[nFacesFromPoints++] = objectMap
1282                 (
1283                     newFaceI,
1284                     selectFaces
1285                     (
1286                         mesh,
1287                         mesh.pointFaces()[iter()],
1288                         true
1289                     )
1290                 );
1291             }
1292             else
1293             {
1294                 // Get patch faces using point on old mesh
1295                 facesFromPoints[nFacesFromPoints++] = objectMap
1296                 (
1297                     newFaceI,
1298                     selectFaces
1299                     (
1300                         mesh,
1301                         mesh.pointFaces()[iter()],
1302                         false
1303                     )
1304                 );
1305             }
1306         }
1307     }
1308
1309
1310     // Faces inflated from edges
1311     // ~~~~~~~~~~~~~~~~~~~~~~~~~
1312
1313     facesFromEdges.setSize(faceFromEdge_.size());
1314
1315     if (faceFromEdge_.size())
1316     {
1317         label nFacesFromEdges = 0;
1318
1319         // Collect all still existing faces connected to this edge.
1320         forAllConstIter(Map<label>, faceFromEdge_, iter)
1321         {
1322             label newFaceI = iter.key();
1323
1324             if (region_[newFaceI] == -1)
1325             {
1326                 // Get internal faces using edge on old mesh
1327                 facesFromEdges[nFacesFromEdges++] = objectMap
1328                 (
1329                     newFaceI,
1330                     selectFaces
1331                     (
1332                         mesh,
1333                         mesh.edgeFaces(iter()),
1334                         true
1335                     )
1336                 );
1337             }
1338             else
1339             {
1340                 // Get patch faces using edge on old mesh
1341                 facesFromEdges[nFacesFromEdges++] = objectMap
1342                 (
1343                     newFaceI,
1344                     selectFaces
1345                     (
1346                         mesh,
1347                         mesh.edgeFaces(iter()),
1348                         false
1349                     )
1350                 );
1351             }
1352         }
1353     }
1354
1355
1356     // Faces from face merging
1357     // ~~~~~~~~~~~~~~~~~~~~~~~
1358
1359     getMergeSets
1360     (
1361         reverseFaceMap_,
1362         faceMap_,
1363         facesFromFaces
1364     );
1365 }
1366
1367
1368 void Foam::polyTopoChange::calcCellInflationMaps
1369 (
1370     const polyMesh& mesh,
1371     List<objectMap>& cellsFromPoints,
1372     List<objectMap>& cellsFromEdges,
1373     List<objectMap>& cellsFromFaces,
1374     List<objectMap>& cellsFromCells
1375 ) const
1376 {
1377     cellsFromPoints.setSize(cellFromPoint_.size());
1378
1379     if (cellFromPoint_.size())
1380     {
1381         label nCellsFromPoints = 0;
1382
1383         // Collect all still existing faces connected to this point.
1384         forAllConstIter(Map<label>, cellFromPoint_, iter)
1385         {
1386             cellsFromPoints[nCellsFromPoints++] = objectMap
1387             (
1388                 iter.key(),
1389                 mesh.pointCells()[iter()]
1390             );
1391         }
1392     }
1393
1394
1395     cellsFromEdges.setSize(cellFromEdge_.size());
1396
1397     if (cellFromEdge_.size())
1398     {
1399         label nCellsFromEdges = 0;
1400
1401         // Collect all still existing faces connected to this point.
1402         forAllConstIter(Map<label>, cellFromEdge_, iter)
1403         {
1404             cellsFromEdges[nCellsFromEdges++] = objectMap
1405             (
1406                 iter.key(),
1407                 mesh.edgeCells()[iter()]
1408             );
1409         }
1410     }
1411
1412
1413     cellsFromFaces.setSize(cellFromFace_.size());
1414
1415     if (cellFromFace_.size())
1416     {
1417         label nCellsFromFaces = 0;
1418
1419         labelList cellsAroundFace(2, -1);
1420
1421         // Collect all still existing faces connected to this point.
1422         forAllConstIter(Map<label>, cellFromFace_, iter)
1423         {
1424             label oldFaceI = iter();
1425
1426             if (mesh.isInternalFace(oldFaceI))
1427             {
1428                 cellsAroundFace[0] = mesh.faceOwner()[oldFaceI];
1429                 cellsAroundFace[1] = mesh.faceNeighbour()[oldFaceI];
1430             }
1431             else
1432             {
1433                 cellsAroundFace[0] = mesh.faceOwner()[oldFaceI];
1434                 cellsAroundFace[1] = -1;
1435             }
1436
1437             cellsFromFaces[nCellsFromFaces++] = objectMap
1438             (
1439                 iter.key(),
1440                 cellsAroundFace
1441             );
1442         }
1443     }
1444
1445
1446     // Cells from cell merging
1447     // ~~~~~~~~~~~~~~~~~~~~~~~
1448
1449     getMergeSets
1450     (
1451         reverseCellMap_,
1452         cellMap_,
1453         cellsFromCells
1454     );
1455 }
1456
1457
1458 void Foam::polyTopoChange::resetZones
1459 (
1460     const polyMesh& mesh,
1461     polyMesh& newMesh,
1462     labelListList& pointZoneMap,
1463     labelListList& faceZoneFaceMap,
1464     labelListList& cellZoneMap
1465 ) const
1466 {
1467     // pointZones
1468     // ~~~~~~~~~~
1469
1470     pointZoneMap.setSize(mesh.pointZones().size());
1471     {
1472         const pointZoneMesh& pointZones = mesh.pointZones();
1473
1474         // Count points per zone
1475
1476         labelList nPoints(pointZones.size(), 0);
1477
1478         forAllConstIter(Map<label>, pointZone_, iter)
1479         {
1480             label zoneI = iter();
1481
1482             if (zoneI < 0 || zoneI >= pointZones.size())
1483             {
1484                 FatalErrorIn
1485                 (
1486                     "resetZones(const polyMesh&, polyMesh&, labelListList&"
1487                     "labelListList&, labelListList&)"
1488                 )   << "Illegal zoneID " << zoneI << " for point "
1489                     << iter.key() << " coord " << mesh.points()[iter.key()]
1490                     << abort(FatalError);
1491             }
1492             nPoints[zoneI]++;
1493         }
1494
1495         // Distribute points per zone
1496
1497         labelListList addressing(pointZones.size());
1498         forAll(addressing, zoneI)
1499         {
1500             addressing[zoneI].setSize(nPoints[zoneI]);
1501         }
1502         nPoints = 0;
1503
1504         forAllConstIter(Map<label>, pointZone_, iter)
1505         {
1506             label zoneI = iter();
1507
1508             addressing[zoneI][nPoints[zoneI]++] = iter.key();
1509         }
1510         // Sort the addressing
1511         forAll(addressing, zoneI)
1512         {
1513             stableSort(addressing[zoneI]);
1514         }
1515
1516         // So now we both have old zones and the new addressing.
1517         // Invert the addressing to get pointZoneMap.
1518         forAll(addressing, zoneI)
1519         {
1520             const pointZone& oldZone = pointZones[zoneI];
1521             const labelList& newZoneAddr = addressing[zoneI];
1522
1523             labelList& curPzRnb = pointZoneMap[zoneI];
1524             curPzRnb.setSize(newZoneAddr.size());
1525
1526             forAll(newZoneAddr, i)
1527             {
1528                 if (newZoneAddr[i] < pointMap_.size())
1529                 {
1530                     curPzRnb[i] = oldZone.whichPoint(pointMap_[newZoneAddr[i]]);
1531                 }
1532                 else
1533                 {
1534                     curPzRnb[i] = -1;
1535                 }
1536             }
1537         }
1538
1539         // Reset the addresing on the zone
1540         newMesh.pointZones().clearAddressing();
1541         forAll(newMesh.pointZones(), zoneI)
1542         {
1543             if (debug)
1544             {
1545                 Pout<< "pointZone:" << zoneI
1546                     << "  name:" << newMesh.pointZones()[zoneI].name()
1547                     << "  size:" << addressing[zoneI].size()
1548                     << endl;
1549             }
1550
1551             newMesh.pointZones()[zoneI] = addressing[zoneI];
1552         }
1553     }
1554
1555
1556     // faceZones
1557     // ~~~~~~~~~
1558
1559     faceZoneFaceMap.setSize(mesh.faceZones().size());
1560     {
1561         const faceZoneMesh& faceZones = mesh.faceZones();
1562
1563         labelList nFaces(faceZones.size(), 0);
1564
1565         forAllConstIter(Map<label>, faceZone_, iter)
1566         {
1567             label zoneI = iter();
1568
1569             if (zoneI < 0 || zoneI >= faceZones.size())
1570             {
1571                 FatalErrorIn
1572                 (
1573                     "resetZones(const polyMesh&, polyMesh&, labelListList&"
1574                     "labelListList&, labelListList&)"
1575                 )   << "Illegal zoneID " << zoneI << " for face "
1576                     << iter.key()
1577                     << abort(FatalError);
1578             }
1579             nFaces[zoneI]++;
1580         }
1581
1582         labelListList addressing(faceZones.size());
1583         boolListList flipMode(faceZones.size());
1584
1585         forAll(addressing, zoneI)
1586         {
1587             addressing[zoneI].setSize(nFaces[zoneI]);
1588             flipMode[zoneI].setSize(nFaces[zoneI]);
1589         }
1590         nFaces = 0;
1591
1592         forAllConstIter(Map<label>, faceZone_, iter)
1593         {
1594             label zoneI = iter();
1595             label faceI = iter.key();
1596
1597             label index = nFaces[zoneI]++;
1598
1599             addressing[zoneI][index] = faceI;
1600             flipMode[zoneI][index] = faceZoneFlip_[faceI];
1601         }
1602         // Sort the addressing
1603         forAll(addressing, zoneI)
1604         {
1605             labelList newToOld;
1606             sortedOrder(addressing[zoneI], newToOld);
1607             {
1608                 labelList newAddressing(addressing[zoneI].size());
1609                 forAll(newAddressing, i)
1610                 {
1611                     newAddressing[i] = addressing[zoneI][newToOld[i]];
1612                 }
1613                 addressing[zoneI].transfer(newAddressing);
1614             }
1615             {
1616                 boolList newFlipMode(flipMode[zoneI].size());
1617                 forAll(newFlipMode, i)
1618                 {
1619                     newFlipMode[i] = flipMode[zoneI][newToOld[i]];
1620                 }
1621                 flipMode[zoneI].transfer(newFlipMode);
1622             }
1623         }
1624
1625         // So now we both have old zones and the new addressing.
1626         // Invert the addressing to get faceZoneFaceMap.
1627         forAll(addressing, zoneI)
1628         {
1629             const faceZone& oldZone = faceZones[zoneI];
1630             const labelList& newZoneAddr = addressing[zoneI];
1631
1632             labelList& curFzFaceRnb = faceZoneFaceMap[zoneI];
1633
1634             curFzFaceRnb.setSize(newZoneAddr.size());
1635
1636             forAll(newZoneAddr, i)
1637             {
1638                 if (newZoneAddr[i] < faceMap_.size())
1639                 {
1640                     curFzFaceRnb[i] =
1641                         oldZone.whichFace(faceMap_[newZoneAddr[i]]);
1642                 }
1643                 else
1644                 {
1645                     curFzFaceRnb[i] = -1;
1646                 }
1647             }
1648         }
1649
1650
1651         // Reset the addresing on the zone
1652         newMesh.faceZones().clearAddressing();
1653         forAll(newMesh.faceZones(), zoneI)
1654         {
1655             if (debug)
1656             {
1657                 Pout<< "faceZone:" << zoneI
1658                     << "  name:" << newMesh.faceZones()[zoneI].name()
1659                     << "  size:" << addressing[zoneI].size()
1660                     << endl;
1661             }
1662
1663             newMesh.faceZones()[zoneI].resetAddressing
1664             (
1665                 addressing[zoneI],
1666                 flipMode[zoneI]
1667             );
1668         }
1669     }
1670
1671
1672     // cellZones
1673     // ~~~~~~~~~
1674
1675     cellZoneMap.setSize(mesh.cellZones().size());
1676     {
1677         const cellZoneMesh& cellZones = mesh.cellZones();
1678
1679         labelList nCells(cellZones.size(), 0);
1680
1681         forAll(cellZone_, cellI)
1682         {
1683             label zoneI = cellZone_[cellI];
1684
1685             if (zoneI >= cellZones.size())
1686             {
1687                 FatalErrorIn
1688                 (
1689                     "resetZones(const polyMesh&, polyMesh&, labelListList&"
1690                     "labelListList&, labelListList&)"
1691                 )   << "Illegal zoneID " << zoneI << " for cell "
1692                     << cellI << abort(FatalError);
1693             }
1694
1695             if (zoneI >= 0)
1696             {
1697                 nCells[zoneI]++;
1698             }
1699         }
1700
1701         labelListList addressing(cellZones.size());
1702         forAll(addressing, zoneI)
1703         {
1704             addressing[zoneI].setSize(nCells[zoneI]);
1705         }
1706         nCells = 0;
1707
1708         forAll(cellZone_, cellI)
1709         {
1710             label zoneI = cellZone_[cellI];
1711
1712             if (zoneI >= 0)
1713             {
1714                 addressing[zoneI][nCells[zoneI]++] = cellI;
1715             }
1716         }
1717         // Sort the addressing
1718         forAll(addressing, zoneI)
1719         {
1720             stableSort(addressing[zoneI]);
1721         }
1722
1723         // So now we both have old zones and the new addressing.
1724         // Invert the addressing to get cellZoneMap.
1725         forAll(addressing, zoneI)
1726         {
1727             const cellZone& oldZone = cellZones[zoneI];
1728             const labelList& newZoneAddr = addressing[zoneI];
1729
1730             labelList& curCellRnb = cellZoneMap[zoneI];
1731
1732             curCellRnb.setSize(newZoneAddr.size());
1733
1734             forAll(newZoneAddr, i)
1735             {
1736                 if (newZoneAddr[i] < cellMap_.size())
1737                 {
1738                     curCellRnb[i] =
1739                         oldZone.whichCell(cellMap_[newZoneAddr[i]]);
1740                 }
1741                 else
1742                 {
1743                     curCellRnb[i] = -1;
1744                 }
1745             }
1746         }
1747
1748         // Reset the addresing on the zone
1749         newMesh.cellZones().clearAddressing();
1750         forAll(newMesh.cellZones(), zoneI)
1751         {
1752             if (debug)
1753             {
1754                 Pout<< "cellZone:" << zoneI
1755                     << "  name:" << newMesh.cellZones()[zoneI].name()
1756                     << "  size:" << addressing[zoneI].size()
1757                     << endl;
1758             }
1759
1760             newMesh.cellZones()[zoneI] = addressing[zoneI];
1761         }
1762     }
1763 }
1764
1765
1766 void Foam::polyTopoChange::calcFaceZonePointMap
1767 (
1768     const polyMesh& mesh,
1769     const List<Map<label> >& oldFaceZoneMeshPointMaps,
1770     labelListList& faceZonePointMap
1771 ) const
1772 {
1773     const faceZoneMesh& faceZones = mesh.faceZones();
1774
1775     faceZonePointMap.setSize(faceZones.size());
1776
1777     forAll(faceZones, zoneI)
1778     {
1779         const faceZone& newZone = faceZones[zoneI];
1780
1781         const labelList& newZoneMeshPoints = newZone().meshPoints();
1782
1783         const Map<label>& oldZoneMeshPointMap = oldFaceZoneMeshPointMaps[zoneI];
1784
1785         labelList& curFzPointRnb = faceZonePointMap[zoneI];
1786
1787         curFzPointRnb.setSize(newZoneMeshPoints.size());
1788
1789         forAll(newZoneMeshPoints, pointI)
1790         {
1791             if (newZoneMeshPoints[pointI] < pointMap_.size())
1792             {
1793                 Map<label>::const_iterator ozmpmIter =
1794                     oldZoneMeshPointMap.find
1795                     (
1796                         pointMap_[newZoneMeshPoints[pointI]]
1797                     );
1798
1799                 if (ozmpmIter != oldZoneMeshPointMap.end())
1800                 {
1801                     curFzPointRnb[pointI] = ozmpmIter();
1802                 }
1803                 else
1804                 {
1805                     curFzPointRnb[pointI] = -1;
1806                 }
1807             }
1808             else
1809             {
1810                 curFzPointRnb[pointI] = -1;
1811             }
1812         }
1813     }
1814 }
1815
1816
1817 Foam::face Foam::polyTopoChange::rotateFace
1818 (
1819     const face& f,
1820     const label nPos
1821 )
1822 {
1823     face newF(f.size());
1824
1825     forAll(f, fp)
1826     {
1827         label fp1 = (fp + nPos) % f.size();
1828
1829         if (fp1 < 0)
1830         {
1831             fp1 += f.size();
1832         }
1833
1834         newF[fp1] = f[fp];
1835     }
1836
1837     return newF;
1838 }
1839
1840
1841 void Foam::polyTopoChange::reorderCoupledFaces
1842 (
1843     const bool syncParallel,
1844     const polyBoundaryMesh& boundary,
1845     const labelList& patchStarts,
1846     const labelList& patchSizes,
1847     const pointField& points
1848 )
1849 {
1850     // Mapping for faces (old to new). Extends over all mesh faces for
1851     // convenience (could be just the external faces)
1852     labelList oldToNew(identity(faces_.size()));
1853
1854     // Rotation on new faces.
1855     labelList rotation(faces_.size(), 0);
1856
1857     // Send ordering
1858     forAll(boundary, patchI)
1859     {
1860         if (syncParallel || !isA<processorPolyPatch>(boundary[patchI]))
1861         {
1862             boundary[patchI].initOrder
1863             (
1864                 primitivePatch
1865                 (
1866                     SubList<face>
1867                     (
1868                         faces_,
1869                         patchSizes[patchI],
1870                         patchStarts[patchI]
1871                     ),
1872                     points
1873                 )
1874             );
1875         }
1876     }
1877
1878     // Receive and calculate ordering
1879
1880     bool anyChanged = false;
1881
1882     forAll(boundary, patchI)
1883     {
1884         if (syncParallel || !isA<processorPolyPatch>(boundary[patchI]))
1885         {
1886             labelList patchFaceMap(patchSizes[patchI], -1);
1887             labelList patchFaceRotation(patchSizes[patchI], 0);
1888
1889             bool changed = boundary[patchI].order
1890             (
1891                 primitivePatch
1892                 (
1893                     SubList<face>
1894                     (
1895                         faces_,
1896                         patchSizes[patchI],
1897                         patchStarts[patchI]
1898                     ),
1899                     points
1900                 ),
1901                 patchFaceMap,
1902                 patchFaceRotation
1903             );
1904
1905             if (changed)
1906             {
1907                 // Merge patch face reordering into mesh face reordering table
1908                 label start = patchStarts[patchI];
1909
1910                 forAll(patchFaceMap, patchFaceI)
1911                 {
1912                     oldToNew[patchFaceI + start] =
1913                         start + patchFaceMap[patchFaceI];
1914                 }
1915
1916                 forAll(patchFaceRotation, patchFaceI)
1917                 {
1918                     rotation[patchFaceI + start] =
1919                         patchFaceRotation[patchFaceI];
1920                 }
1921
1922                 anyChanged = true;
1923             }
1924         }
1925     }
1926
1927     if (syncParallel)
1928     {
1929         reduce(anyChanged, orOp<bool>());
1930     }
1931
1932     if (anyChanged)
1933     {
1934         // Reorder faces according to oldToNew.
1935         reorderCompactFaces(oldToNew.size(), oldToNew);
1936
1937         // Rotate faces (rotation is already in new face indices).
1938         forAll(rotation, faceI)
1939         {
1940             if (rotation[faceI] != 0)
1941             {
1942                 faces_[faceI] = rotateFace(faces_[faceI], rotation[faceI]);
1943             }
1944         }
1945     }
1946 }
1947
1948
1949 void Foam::polyTopoChange::compactAndReorder
1950 (
1951     const polyMesh& mesh,
1952     const bool syncParallel,
1953     const bool orderCells,
1954     const bool orderPoints,
1955
1956     label& nInternalPoints,
1957     pointField& newPoints,
1958     labelList& patchSizes,
1959     labelList& patchStarts,
1960     List<objectMap>& pointsFromPoints,
1961     List<objectMap>& facesFromPoints,
1962     List<objectMap>& facesFromEdges,
1963     List<objectMap>& facesFromFaces,
1964     List<objectMap>& cellsFromPoints,
1965     List<objectMap>& cellsFromEdges,
1966     List<objectMap>& cellsFromFaces,
1967     List<objectMap>& cellsFromCells,
1968     List<Map<label> >& oldPatchMeshPointMaps,
1969     labelList& oldPatchNMeshPoints,
1970     labelList& oldPatchStarts,
1971     List<Map<label> >& oldFaceZoneMeshPointMaps
1972 )
1973 {
1974     if (mesh.boundaryMesh().size() != nPatches_)
1975     {
1976         FatalErrorIn("polyTopoChange::compactAndReorder(..)")
1977             << "polyTopoChange was constructed with a mesh with "
1978             << nPatches_ << " patches." << endl
1979             << "The mesh now provided has a different number of patches "
1980             << mesh.boundaryMesh().size()
1981             << " which is illegal" << endl
1982             << abort(FatalError);
1983     }
1984
1985     // Remove any holes from points/faces/cells and sort faces.
1986     // Sets nActiveFaces_.
1987     compact(orderCells, orderPoints, nInternalPoints, patchSizes, patchStarts);
1988
1989     // Transfer points to pointField. points_ are now cleared!
1990     // Only done since e.g. reorderCoupledFaces requires pointField.
1991     newPoints.transfer(points_);
1992
1993     // Reorder any coupled faces
1994     reorderCoupledFaces
1995     (
1996         syncParallel,
1997         mesh.boundaryMesh(),
1998         patchStarts,
1999         patchSizes,
2000         newPoints
2001     );
2002
2003
2004     // Calculate inflation/merging maps
2005     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2006     // These are for the new face(/point/cell) the old faces whose value
2007     // needs to be
2008     // averaged/summed to get the new value. These old faces come either from
2009     // merged old faces (face remove into other face),
2010     // the old edgeFaces (inflate from edge) or the old pointFaces (inflate
2011     // from point). As an additional complexity will use only internal faces
2012     // to create new value for internal face and vice versa only patch
2013     // faces to to create patch face value.
2014
2015     // For point only point merging
2016     getMergeSets
2017     (
2018         reversePointMap_,
2019         pointMap_,
2020         pointsFromPoints
2021     );
2022
2023     calcFaceInflationMaps
2024     (
2025         mesh,
2026         facesFromPoints,
2027         facesFromEdges,
2028         facesFromFaces
2029     );
2030
2031     calcCellInflationMaps
2032     (
2033         mesh,
2034         cellsFromPoints,
2035         cellsFromEdges,
2036         cellsFromFaces,
2037         cellsFromCells
2038     );
2039
2040     // Clear inflation info
2041     {
2042         faceFromPoint_.clearStorage();
2043         faceFromEdge_.clearStorage();
2044
2045         cellFromPoint_.clearStorage();
2046         cellFromEdge_.clearStorage();
2047         cellFromFace_.clearStorage();
2048     }
2049
2050
2051     const polyBoundaryMesh& boundary = mesh.boundaryMesh();
2052
2053     // Grab patch mesh point maps
2054     oldPatchMeshPointMaps.setSize(boundary.size());
2055     oldPatchNMeshPoints.setSize(boundary.size());
2056     oldPatchStarts.setSize(boundary.size());
2057
2058     forAll(boundary, patchI)
2059     {
2060         // Copy old face zone mesh point maps
2061         oldPatchMeshPointMaps[patchI] = boundary[patchI].meshPointMap();
2062         oldPatchNMeshPoints[patchI] = boundary[patchI].meshPoints().size();
2063         oldPatchStarts[patchI] = boundary[patchI].start();
2064     }
2065
2066     // Grab old face zone mesh point maps.
2067     // These need to be saved before resetting the mesh and are used
2068     // later on to calculate the faceZone pointMaps.
2069     oldFaceZoneMeshPointMaps.setSize(mesh.faceZones().size());
2070
2071     forAll(mesh.faceZones(), zoneI)
2072     {
2073         const faceZone& oldZone = mesh.faceZones()[zoneI];
2074
2075         oldFaceZoneMeshPointMaps[zoneI] = oldZone().meshPointMap();
2076     }
2077 }
2078
2079
2080 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
2081
2082 // Construct from components
2083 Foam::polyTopoChange::polyTopoChange(const label nPatches, const bool strict)
2084 :
2085     strict_(strict),
2086     nPatches_(nPatches),
2087     points_(0),
2088     pointMap_(0),
2089     reversePointMap_(0),
2090     pointZone_(0),
2091     retiredPoints_(0),
2092     faces_(0),
2093     region_(0),
2094     faceOwner_(0),
2095     faceNeighbour_(0),
2096     faceMap_(0),
2097     reverseFaceMap_(0),
2098     faceFromPoint_(0),
2099     faceFromEdge_(0),
2100     flipFaceFlux_(0),
2101     faceZone_(0),
2102     faceZoneFlip_(0),
2103     nActiveFaces_(0),
2104     cellMap_(0),
2105     reverseCellMap_(0),
2106     cellFromPoint_(0),
2107     cellFromEdge_(0),
2108     cellFromFace_(0),
2109     cellZone_(0)
2110 {}
2111
2112
2113 // Construct from components
2114 Foam::polyTopoChange::polyTopoChange
2115 (
2116     const polyMesh& mesh,
2117     const bool strict
2118 )
2119 :
2120     strict_(strict),
2121     nPatches_(0),
2122     points_(0),
2123     pointMap_(0),
2124     reversePointMap_(0),
2125     pointZone_(0),
2126     retiredPoints_(0),
2127     faces_(0),
2128     region_(0),
2129     faceOwner_(0),
2130     faceNeighbour_(0),
2131     faceMap_(0),
2132     reverseFaceMap_(0),
2133     faceFromPoint_(0),
2134     faceFromEdge_(0),
2135     flipFaceFlux_(0),
2136     faceZone_(0),
2137     faceZoneFlip_(0),
2138     nActiveFaces_(0),
2139     cellMap_(0),
2140     reverseCellMap_(0),
2141     cellFromPoint_(0),
2142     cellFromEdge_(0),
2143     cellFromFace_(0),
2144     cellZone_(0)
2145 {
2146     addMesh
2147     (
2148         mesh,
2149         identity(mesh.boundaryMesh().size()),
2150         identity(mesh.pointZones().size()),
2151         identity(mesh.faceZones().size()),
2152         identity(mesh.cellZones().size())
2153     );
2154 }
2155
2156
2157 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
2158
2159 void Foam::polyTopoChange::clear()
2160 {
2161     points_.clearStorage();
2162     pointMap_.clearStorage();
2163     reversePointMap_.clearStorage();
2164     pointZone_.clearStorage();
2165     retiredPoints_.clearStorage();
2166
2167     faces_.clearStorage();
2168     region_.clearStorage();
2169     faceOwner_.clearStorage();
2170     faceNeighbour_.clearStorage();
2171     faceMap_.clearStorage();
2172     reverseFaceMap_.clearStorage();
2173     faceFromPoint_.clearStorage();
2174     faceFromEdge_.clearStorage();
2175     flipFaceFlux_.clearStorage();
2176     faceZone_.clearStorage();
2177     faceZoneFlip_.clearStorage();
2178     nActiveFaces_ = 0;
2179
2180     cellMap_.clearStorage();
2181     reverseCellMap_.clearStorage();
2182     cellZone_.clearStorage();
2183     cellFromPoint_.clearStorage();
2184     cellFromEdge_.clearStorage();
2185     cellFromFace_.clearStorage();
2186 }
2187
2188
2189 void Foam::polyTopoChange::addMesh
2190 (
2191     const polyMesh& mesh,
2192     const labelList& patchMap,
2193     const labelList& pointZoneMap,
2194     const labelList& faceZoneMap,
2195     const labelList& cellZoneMap
2196 )
2197 {
2198     label maxRegion = nPatches_ - 1;
2199     forAll(patchMap, i)
2200     {
2201         maxRegion = max(maxRegion, patchMap[i]);
2202     }
2203     nPatches_ = maxRegion + 1;
2204
2205
2206     // Add points
2207     {
2208         const pointField& points = mesh.points();
2209         const pointZoneMesh& pointZones = mesh.pointZones();
2210
2211         // Extend
2212         points_.setCapacity(points_.size() + points.size());
2213         pointMap_.setCapacity(pointMap_.size() + points.size());
2214         reversePointMap_.setCapacity(reversePointMap_.size() + points.size());
2215         pointZone_.resize(pointZone_.size() + points.size()/100);
2216
2217         // Precalc offset zones
2218         labelList newZoneID(points.size(), -1);
2219
2220         forAll(pointZones, zoneI)
2221         {
2222             const labelList& pointLabels = pointZones[zoneI];
2223
2224             forAll(pointLabels, j)
2225             {
2226                 newZoneID[pointLabels[j]] = pointZoneMap[zoneI];
2227             }
2228         }
2229
2230         // Add points in mesh order
2231         for (label pointI = 0; pointI < mesh.nPoints(); pointI++)
2232         {
2233             addPoint
2234             (
2235                 points[pointI],
2236                 pointI,
2237                 newZoneID[pointI],
2238                 true
2239             );
2240         }
2241     }
2242
2243     // Add cells
2244     {
2245         const cellZoneMesh& cellZones = mesh.cellZones();
2246
2247         // Resize
2248
2249         // Note: polyMesh does not allow retired cells anymore. So allCells
2250         // always equals nCells
2251         label nAllCells = mesh.nCells();
2252
2253         cellMap_.setCapacity(cellMap_.size() + nAllCells);
2254         reverseCellMap_.setCapacity(reverseCellMap_.size() + nAllCells);
2255         cellFromPoint_.resize(cellFromPoint_.size() + nAllCells/100);
2256         cellFromEdge_.resize(cellFromEdge_.size() + nAllCells/100);
2257         cellFromFace_.resize(cellFromFace_.size() + nAllCells/100);
2258         cellZone_.setCapacity(cellZone_.size() + nAllCells);
2259
2260
2261         // Precalc offset zones
2262         labelList newZoneID(nAllCells, -1);
2263
2264         forAll(cellZones, zoneI)
2265         {
2266             const labelList& cellLabels = cellZones[zoneI];
2267
2268             forAll(cellLabels, j)
2269             {
2270                 label cellI = cellLabels[j];
2271
2272                 if (newZoneID[cellI] != -1)
2273                 {
2274                     WarningIn
2275                     (
2276                         "polyTopoChange::addMesh"
2277                         "(const polyMesh&, const labelList&,"
2278                         "const labelList&, const labelList&,"
2279                         "const labelList&)"
2280                     )   << "Cell:" << cellI
2281                         << " centre:" << mesh.cellCentres()[cellI]
2282                         << " is in two zones:"
2283                         << cellZones[newZoneID[cellI]].name()
2284                         << " and " << cellZones[zoneI].name() << endl
2285                         << "    This is not supported."
2286                         << " Continuing with first zone only." << endl;
2287                 }
2288                 else
2289                 {
2290                     newZoneID[cellI] = cellZoneMap[zoneI];
2291                 }
2292             }
2293         }
2294
2295         // Add cells in mesh order
2296         for (label cellI = 0; cellI < nAllCells; cellI++)
2297         {
2298             // Add cell from cell
2299             addCell(-1, -1, -1, cellI, newZoneID[cellI]);
2300         }
2301     }
2302
2303     // Add faces
2304     {
2305         const polyBoundaryMesh& patches = mesh.boundaryMesh();
2306         const faceList& faces = mesh.faces();
2307         const labelList& faceOwner = mesh.faceOwner();
2308         const labelList& faceNeighbour = mesh.faceNeighbour();
2309         const faceZoneMesh& faceZones = mesh.faceZones();
2310
2311         // Resize
2312         label nAllFaces = mesh.faces().size();
2313
2314         faces_.setCapacity(faces_.size() + nAllFaces);
2315         region_.setCapacity(region_.size() + nAllFaces);
2316         faceOwner_.setCapacity(faceOwner_.size() + nAllFaces);
2317         faceNeighbour_.setCapacity(faceNeighbour_.size() + nAllFaces);
2318         faceMap_.setCapacity(faceMap_.size() + nAllFaces);
2319         reverseFaceMap_.setCapacity(reverseFaceMap_.size() + nAllFaces);
2320         faceFromPoint_.resize(faceFromPoint_.size() + nAllFaces/100);
2321         faceFromEdge_.resize(faceFromEdge_.size() + nAllFaces/100);
2322         flipFaceFlux_.setCapacity(faces_.size() + nAllFaces);
2323         faceZone_.resize(faceZone_.size() + nAllFaces/100);
2324         faceZoneFlip_.setCapacity(faces_.size() + nAllFaces);
2325
2326
2327         // Precalc offset zones
2328         labelList newZoneID(nAllFaces, -1);
2329         boolList zoneFlip(nAllFaces, false);
2330
2331         forAll(faceZones, zoneI)
2332         {
2333             const labelList& faceLabels = faceZones[zoneI];
2334             const boolList& flipMap = faceZones[zoneI].flipMap();
2335
2336             forAll(faceLabels, j)
2337             {
2338                 newZoneID[faceLabels[j]] = faceZoneMap[zoneI];
2339                 zoneFlip[faceLabels[j]] = flipMap[j];
2340             }
2341         }
2342
2343         // Add faces in mesh order
2344
2345         // 1. Internal faces
2346         for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
2347         {
2348             addFace
2349             (
2350                 faces[faceI],
2351                 faceOwner[faceI],
2352                 faceNeighbour[faceI],
2353                 -1,                         // masterPointID
2354                 -1,                         // masterEdgeID
2355                 faceI,                      // masterFaceID
2356                 false,                      // flipFaceFlux
2357                 -1,                         // patchID
2358                 newZoneID[faceI],           // zoneID
2359                 zoneFlip[faceI]             // zoneFlip
2360             );
2361         }
2362
2363         // 2. Patch faces
2364         forAll(patches, patchI)
2365         {
2366             const polyPatch& pp = patches[patchI];
2367
2368             if (pp.start() != faces_.size())
2369             {
2370                 FatalErrorIn
2371                 (
2372                     "polyTopoChange::polyTopoChange"
2373                     "(const polyMesh& mesh, const bool strict)"
2374                 )   << "Problem : "
2375                     << "Patch " << pp.name() << " starts at " << pp.start()
2376                     << endl
2377                     << "Current face counter at " << faces_.size() << endl
2378                     << "Are patches in incremental order?"
2379                     << abort(FatalError);
2380             }
2381             forAll(pp, patchFaceI)
2382             {
2383                 label faceI = pp.start() + patchFaceI;
2384
2385                 addFace
2386                 (
2387                     faces[faceI],
2388                     faceOwner[faceI],
2389                     -1,                         // neighbour
2390                     -1,                         // masterPointID
2391                     -1,                         // masterEdgeID
2392                     faceI,                      // masterFaceID
2393                     false,                      // flipFaceFlux
2394                     patchMap[patchI],           // patchID
2395                     newZoneID[faceI],           // zoneID
2396                     zoneFlip[faceI]             // zoneFlip
2397                 );
2398             }
2399         }
2400     }
2401 }
2402
2403
2404 void Foam::polyTopoChange::setCapacity
2405 (
2406     const label nPoints,
2407     const label nFaces,
2408     const label nCells
2409 )
2410 {
2411     points_.setCapacity(nPoints);
2412     pointMap_.setCapacity(nPoints);
2413     reversePointMap_.setCapacity(nPoints);
2414     pointZone_.resize(pointZone_.size() + nPoints/100);
2415
2416     faces_.setCapacity(nFaces);
2417     region_.setCapacity(nFaces);
2418     faceOwner_.setCapacity(nFaces);
2419     faceNeighbour_.setCapacity(nFaces);
2420     faceMap_.setCapacity(nFaces);
2421     reverseFaceMap_.setCapacity(nFaces);
2422     faceFromPoint_.resize(faceFromPoint_.size() + nFaces/100);
2423     faceFromEdge_.resize(faceFromEdge_.size() + nFaces/100);
2424     flipFaceFlux_.setCapacity(nFaces);
2425     faceZone_.resize(faceZone_.size() + nFaces/100);
2426     faceZoneFlip_.setCapacity(nFaces);
2427
2428     cellMap_.setCapacity(nCells);
2429     reverseCellMap_.setCapacity(nCells);
2430     cellFromPoint_.resize(cellFromPoint_.size() + nCells/100);
2431     cellFromEdge_.resize(cellFromEdge_.size() + nCells/100);
2432     cellFromFace_.resize(cellFromFace_.size() + nCells/100);
2433     cellZone_.setCapacity(nCells);
2434 }
2435
2436
2437 Foam::label Foam::polyTopoChange::setAction(const topoAction& action)
2438 {
2439     if (isType<polyAddPoint>(action))
2440     {
2441         const polyAddPoint& pap = refCast<const polyAddPoint>(action);
2442
2443         return addPoint
2444         (
2445             pap.newPoint(),
2446             pap.masterPointID(),
2447             pap.zoneID(),
2448             pap.inCell()
2449         );
2450     }
2451     else if (isType<polyModifyPoint>(action))
2452     {
2453         const polyModifyPoint& pmp = refCast<const polyModifyPoint>(action);
2454
2455         modifyPoint
2456         (
2457             pmp.pointID(),
2458             pmp.newPoint(),
2459             pmp.zoneID(),
2460             pmp.inCell()
2461         );
2462
2463         return -1;
2464     }
2465     else if (isType<polyRemovePoint>(action))
2466     {
2467         const polyRemovePoint& prp = refCast<const polyRemovePoint>(action);
2468
2469         removePoint(prp.pointID(), prp.mergePointID());
2470
2471         return -1;
2472     }
2473     else if (isType<polyAddFace>(action))
2474     {
2475         const polyAddFace& paf = refCast<const polyAddFace>(action);
2476
2477         return addFace
2478         (
2479             paf.newFace(),
2480             paf.owner(),
2481             paf.neighbour(),
2482             paf.masterPointID(),
2483             paf.masterEdgeID(),
2484             paf.masterFaceID(),
2485             paf.flipFaceFlux(),
2486             paf.patchID(),
2487             paf.zoneID(),
2488             paf.zoneFlip()
2489         );
2490     }
2491     else if (isType<polyModifyFace>(action))
2492     {
2493         const polyModifyFace& pmf = refCast<const polyModifyFace>(action);
2494
2495         modifyFace
2496         (
2497             pmf.newFace(),
2498             pmf.faceID(),
2499             pmf.owner(),
2500             pmf.neighbour(),
2501             pmf.flipFaceFlux(),
2502             pmf.patchID(),
2503             pmf.zoneID(),
2504             pmf.zoneFlip()
2505         );
2506
2507         return -1;
2508     }
2509     else if (isType<polyRemoveFace>(action))
2510     {
2511         const polyRemoveFace& prf = refCast<const polyRemoveFace>(action);
2512
2513         removeFace(prf.faceID(), prf.mergeFaceID());
2514
2515         return -1;
2516     }
2517     else if (isType<polyAddCell>(action))
2518     {
2519         const polyAddCell& pac = refCast<const polyAddCell>(action);
2520
2521         return addCell
2522         (
2523             pac.masterPointID(),
2524             pac.masterEdgeID(),
2525             pac.masterFaceID(),
2526             pac.masterCellID(),
2527             pac.zoneID()
2528         );
2529     }
2530     else if (isType<polyModifyCell>(action))
2531     {
2532         const polyModifyCell& pmc = refCast<const polyModifyCell>(action);
2533
2534         if (pmc.removeFromZone())
2535         {
2536             modifyCell(pmc.cellID(), -1);
2537         }
2538         else
2539         {
2540             modifyCell(pmc.cellID(), pmc.zoneID());
2541         }
2542
2543         return -1;
2544     }
2545     else if (isType<polyRemoveCell>(action))
2546     {
2547         const polyRemoveCell& prc = refCast<const polyRemoveCell>(action);
2548
2549         removeCell(prc.cellID(), prc.mergeCellID());
2550
2551         return -1;
2552     }
2553     else
2554     {
2555         FatalErrorIn
2556         (
2557             "label polyTopoChange::setAction(const topoAction& action)"
2558         )   << "Unknown type of topoChange: " << action.type()
2559             << abort(FatalError);
2560
2561         // Dummy return to keep compiler happy
2562         return -1;
2563     }
2564 }
2565
2566
2567 Foam::label Foam::polyTopoChange::addPoint
2568 (
2569     const point& pt,
2570     const label masterPointID,
2571     const label zoneID,
2572     const bool inCell
2573 )
2574 {
2575     label pointI = points_.size();
2576
2577     points_.append(pt);
2578     pointMap_.append(masterPointID);
2579     reversePointMap_.append(pointI);
2580
2581     if (zoneID >= 0)
2582     {
2583         pointZone_.insert(pointI, zoneID);
2584     }
2585
2586     if (!inCell)
2587     {
2588         retiredPoints_.insert(pointI);
2589     }
2590
2591     return pointI;
2592 }
2593
2594
2595 void Foam::polyTopoChange::modifyPoint
2596 (
2597     const label pointI,
2598     const point& pt,
2599     const label newZoneID,
2600     const bool inCell
2601 )
2602 {
2603     if (pointI < 0 || pointI >= points_.size())
2604     {
2605         FatalErrorIn
2606         (
2607             "polyTopoChange::modifyPoint(const label, const point&)"
2608         )   << "illegal point label " << pointI << endl
2609             << "Valid point labels are 0 .. " << points_.size()-1
2610             << abort(FatalError);
2611     }
2612     if (pointRemoved(pointI) || pointMap_[pointI] == -1)
2613     {
2614         FatalErrorIn
2615         (
2616             "polyTopoChange::modifyPoint(const label, const point&)"
2617         )   << "point " << pointI << " already marked for removal"
2618             << abort(FatalError);
2619     }
2620     points_[pointI] = pt;
2621
2622     Map<label>::iterator pointFnd = pointZone_.find(pointI);
2623
2624     if (pointFnd != pointZone_.end())
2625     {
2626         if (newZoneID >= 0)
2627         {
2628             pointFnd() = newZoneID;
2629         }
2630         else
2631         {
2632             pointZone_.erase(pointFnd);
2633         }
2634     }
2635     else if (newZoneID >= 0)
2636     {
2637         pointZone_.insert(pointI, newZoneID);
2638     }
2639
2640     if (inCell)
2641     {
2642         retiredPoints_.erase(pointI);
2643     }
2644     else
2645     {
2646         retiredPoints_.insert(pointI);
2647     }
2648 }
2649
2650
2651 void Foam::polyTopoChange::movePoints(const pointField& newPoints)
2652 {
2653     if (newPoints.size() != points_.size())
2654     {
2655         FatalErrorIn("polyTopoChange::movePoints(const pointField&)")
2656             << "illegal pointField size." << endl
2657             << "Size:" << newPoints.size() << endl
2658             << "Points in mesh:" << points_.size()
2659             << abort(FatalError);
2660     }
2661
2662     forAll(points_, pointI)
2663     {
2664         points_[pointI] = newPoints[pointI];
2665     }
2666 }
2667
2668
2669 void Foam::polyTopoChange::removePoint
2670 (
2671     const label pointI,
2672     const label mergePointI
2673 )
2674 {
2675     if (pointI < 0 || pointI >= points_.size())
2676     {
2677         FatalErrorIn("polyTopoChange::removePoint(const label, const label)")
2678             << "illegal point label " << pointI << endl
2679             << "Valid point labels are 0 .. " << points_.size()-1
2680             << abort(FatalError);
2681     }
2682
2683     if
2684     (
2685         strict_
2686      && (pointRemoved(pointI) || pointMap_[pointI] == -1)
2687     )
2688     {
2689         FatalErrorIn("polyTopoChange::removePoint(const label, const label)")
2690             << "point " << pointI << " already marked for removal" << nl
2691             << "Point:" << points_[pointI] << " pointMap:" << pointMap_[pointI]
2692             << abort(FatalError);
2693     }
2694
2695     if (pointI == mergePointI)
2696     {
2697         FatalErrorIn("polyTopoChange::removePoint(const label, const label)")
2698             << "Cannot remove/merge point " << pointI << " onto itself."
2699             << abort(FatalError);
2700     }
2701
2702     points_[pointI] = greatPoint;
2703     pointMap_[pointI] = -1;
2704     if (mergePointI >= 0)
2705     {
2706         reversePointMap_[pointI] = -mergePointI-2;
2707     }
2708     else
2709     {
2710         reversePointMap_[pointI] = -1;
2711     }
2712     pointZone_.erase(pointI);
2713     retiredPoints_.erase(pointI);
2714 }
2715
2716
2717 Foam::label Foam::polyTopoChange::addFace
2718 (
2719     const face& f,
2720     const label own,
2721     const label nei,
2722     const label masterPointID,
2723     const label masterEdgeID,
2724     const label masterFaceID,
2725     const bool flipFaceFlux,
2726     const label patchID,
2727     const label zoneID,
2728     const bool zoneFlip
2729 )
2730 {
2731     // Check validity
2732     if (debug)
2733     {
2734         checkFace(f, -1, own, nei, patchID, zoneID);
2735     }
2736
2737     label faceI = faces_.size();
2738
2739     faces_.append(f);
2740     region_.append(patchID);
2741     faceOwner_.append(own);
2742     faceNeighbour_.append(nei);
2743
2744     if (masterPointID >= 0)
2745     {
2746         faceMap_.append(-1);
2747         faceFromPoint_.insert(faceI, masterPointID);
2748     }
2749     else if (masterEdgeID >= 0)
2750     {
2751         faceMap_.append(-1);
2752         faceFromEdge_.insert(faceI, masterEdgeID);
2753     }
2754     else if (masterFaceID >= 0)
2755     {
2756         faceMap_.append(masterFaceID);
2757     }
2758     else
2759     {
2760         // Allow inflate-from-nothing?
2761         //FatalErrorIn("polyTopoChange::addFace")
2762         //    << "Need to specify a master point, edge or face"
2763         //    << "face:" << f << " own:" << own << " nei:" << nei
2764         //    << abort(FatalError);
2765         faceMap_.append(-1);
2766     }
2767     reverseFaceMap_.append(faceI);
2768
2769     flipFaceFlux_[faceI] = (flipFaceFlux ? 1 : 0);
2770
2771     if (zoneID >= 0)
2772     {
2773         faceZone_.insert(faceI, zoneID);
2774     }
2775     faceZoneFlip_[faceI] = (zoneFlip ? 1 : 0);
2776
2777     return faceI;
2778 }
2779
2780
2781 void Foam::polyTopoChange::modifyFace
2782 (
2783     const face& f,
2784     const label faceI,
2785     const label own,
2786     const label nei,
2787     const bool flipFaceFlux,
2788     const label patchID,
2789     const label zoneID,
2790     const bool zoneFlip
2791 )
2792 {
2793     // Check validity
2794     if (debug)
2795     {
2796         checkFace(f, faceI, own, nei, patchID, zoneID);
2797     }
2798
2799     faces_[faceI] = f;
2800     faceOwner_[faceI] = own;
2801     faceNeighbour_[faceI] = nei;
2802     region_[faceI] = patchID;
2803
2804     flipFaceFlux_[faceI] = (flipFaceFlux ? 1 : 0);
2805
2806     Map<label>::iterator faceFnd = faceZone_.find(faceI);
2807
2808     if (faceFnd != faceZone_.end())
2809     {
2810         if (zoneID >= 0)
2811         {
2812             faceFnd() = zoneID;
2813         }
2814         else
2815         {
2816             faceZone_.erase(faceFnd);
2817         }
2818     }
2819     else if (zoneID >= 0)
2820     {
2821         faceZone_.insert(faceI, zoneID);
2822     }
2823     faceZoneFlip_[faceI] = (zoneFlip ? 1 : 0);
2824 }
2825
2826
2827 void Foam::polyTopoChange::removeFace(const label faceI, const label mergeFaceI)
2828 {
2829     if (faceI < 0 || faceI >= faces_.size())
2830     {
2831         FatalErrorIn("polyTopoChange::removeFace(const label, const label)")
2832             << "illegal face label " << faceI << endl
2833             << "Valid face labels are 0 .. " << faces_.size()-1
2834             << abort(FatalError);
2835     }
2836
2837     if
2838     (
2839         strict_
2840      && (faceRemoved(faceI) || faceMap_[faceI] == -1)
2841     )
2842     {
2843         FatalErrorIn("polyTopoChange::removeFace(const label, const label)")
2844             << "face " << faceI
2845             << " already marked for removal"
2846             << abort(FatalError);
2847     }
2848
2849     faces_[faceI].setSize(0);
2850     region_[faceI] = -1;
2851     faceOwner_[faceI] = -1;
2852     faceNeighbour_[faceI] = -1;
2853     faceMap_[faceI] = -1;
2854     if (mergeFaceI >= 0)
2855     {
2856         reverseFaceMap_[faceI] = -mergeFaceI-2;
2857     }
2858     else
2859     {
2860         reverseFaceMap_[faceI] = -1;
2861     }
2862     faceFromEdge_.erase(faceI);
2863     faceFromPoint_.erase(faceI);
2864     flipFaceFlux_[faceI] = 0;
2865     faceZone_.erase(faceI);
2866     faceZoneFlip_[faceI] = 0;
2867 }
2868
2869
2870 Foam::label Foam::polyTopoChange::addCell
2871 (
2872     const label masterPointID,
2873     const label masterEdgeID,
2874     const label masterFaceID,
2875     const label masterCellID,
2876     const label zoneID
2877 )
2878 {
2879     label cellI = cellMap_.size();
2880
2881     if (masterPointID >= 0)
2882     {
2883         cellMap_.append(-1);
2884         cellFromPoint_.insert(cellI, masterPointID);
2885     }
2886     else if (masterEdgeID >= 0)
2887     {
2888         cellMap_.append(-1);
2889         cellFromEdge_.insert(cellI, masterEdgeID);
2890     }
2891     else if (masterFaceID >= 0)
2892     {
2893         cellMap_.append(-1);
2894         cellFromFace_.insert(cellI, masterFaceID);
2895     }
2896     else
2897     {
2898         cellMap_.append(masterCellID);
2899     }
2900     reverseCellMap_.append(cellI);
2901     cellZone_.append(zoneID);
2902
2903     return cellI;
2904 }
2905
2906
2907 void Foam::polyTopoChange::modifyCell
2908 (
2909     const label cellI,
2910     const label zoneID
2911 )
2912 {
2913     cellZone_[cellI] = zoneID;
2914 }
2915
2916
2917 void Foam::polyTopoChange::removeCell(const label cellI, const label mergeCellI)
2918 {
2919     if (cellI < 0 || cellI >= cellMap_.size())
2920     {
2921         FatalErrorIn("polyTopoChange::removeCell(const label, const label)")
2922             << "illegal cell label " << cellI << endl
2923             << "Valid cell labels are 0 .. " << cellMap_.size()-1
2924             << abort(FatalError);
2925     }
2926
2927     if (strict_ && cellMap_[cellI] == -2)
2928     {
2929         FatalErrorIn("polyTopoChange::removeCell(const label, const label)")
2930             << "cell " << cellI
2931             << " already marked for removal"
2932             << abort(FatalError);
2933     }
2934
2935     cellMap_[cellI] = -2;
2936     if (mergeCellI >= 0)
2937     {
2938         reverseCellMap_[cellI] = -mergeCellI-2;
2939     }
2940     else
2941     {
2942         reverseCellMap_[cellI] = -1;
2943     }
2944     cellFromPoint_.erase(cellI);
2945     cellFromEdge_.erase(cellI);
2946     cellFromFace_.erase(cellI);
2947     cellZone_[cellI] = -1;
2948 }
2949
2950
2951 Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::changeMesh
2952 (
2953     polyMesh& mesh,
2954     const bool inflate,
2955     const bool syncParallel,
2956     const bool orderCells,
2957     const bool orderPoints
2958 )
2959 {
2960     if (debug)
2961     {
2962         Pout<< "polyTopoChange::changeMesh"
2963             << "(polyMesh&, const bool, const bool, const bool, const bool)"
2964             << endl;
2965     }
2966
2967     if (debug)
2968     {
2969         Pout<< "Old mesh:" << nl;
2970         writeMeshStats(mesh, Pout);
2971     }
2972
2973     // new mesh points
2974     pointField newPoints;
2975     // number of internal points
2976     label nInternalPoints;
2977     // patch slicing
2978     labelList patchSizes;
2979     labelList patchStarts;
2980     // inflate maps
2981     List<objectMap> pointsFromPoints;
2982     List<objectMap> facesFromPoints;
2983     List<objectMap> facesFromEdges;
2984     List<objectMap> facesFromFaces;
2985     List<objectMap> cellsFromPoints;
2986     List<objectMap> cellsFromEdges;
2987     List<objectMap> cellsFromFaces;
2988     List<objectMap> cellsFromCells;
2989     // old mesh info
2990     List<Map<label> > oldPatchMeshPointMaps;
2991     labelList oldPatchNMeshPoints;
2992     labelList oldPatchStarts;
2993     List<Map<label> > oldFaceZoneMeshPointMaps;
2994
2995     // Compact, reorder patch faces and calculate mesh/patch maps.
2996     compactAndReorder
2997     (
2998         mesh,
2999         syncParallel,
3000         orderCells,
3001         orderPoints,
3002
3003         nInternalPoints,
3004         newPoints,
3005         patchSizes,
3006         patchStarts,
3007         pointsFromPoints,
3008         facesFromPoints,
3009         facesFromEdges,
3010         facesFromFaces,
3011         cellsFromPoints,
3012         cellsFromEdges,
3013         cellsFromFaces,
3014         cellsFromCells,
3015         oldPatchMeshPointMaps,
3016         oldPatchNMeshPoints,
3017         oldPatchStarts,
3018         oldFaceZoneMeshPointMaps
3019     );
3020
3021     const label nOldPoints(mesh.nPoints());
3022     const label nOldFaces(mesh.nFaces());
3023     const label nOldCells(mesh.nCells());
3024
3025
3026     // Change the mesh
3027     // ~~~~~~~~~~~~~~~
3028     // This will invalidate any addressing so better make sure you have
3029     // all the information you need!!!
3030
3031     if (inflate)
3032     {
3033         // Keep (renumbered) mesh points, store new points in map for inflation
3034         // (appended points (i.e. from nowhere) get value zero)
3035         pointField renumberedMeshPoints(newPoints.size());
3036
3037         forAll(pointMap_, newPointI)
3038         {
3039             label oldPointI = pointMap_[newPointI];
3040
3041             if (oldPointI >= 0)
3042             {
3043                 renumberedMeshPoints[newPointI] = mesh.points()[oldPointI];
3044             }
3045             else
3046             {
3047                 renumberedMeshPoints[newPointI] = vector::zero;
3048             }
3049         }
3050
3051         mesh.resetPrimitives
3052         (
3053             xferMove(renumberedMeshPoints),
3054             faces_.xfer(),
3055             faceOwner_.xfer(),
3056             faceNeighbour_.xfer(),
3057             patchSizes,
3058             patchStarts,
3059             syncParallel
3060         );
3061
3062         mesh.changing(true);
3063     }
3064     else
3065     {
3066         // Set new points.
3067         mesh.resetPrimitives
3068         (
3069             xferMove(newPoints),
3070             faces_.xfer(),
3071             faceOwner_.xfer(),
3072             faceNeighbour_.xfer(),
3073             patchSizes,
3074             patchStarts,
3075             syncParallel
3076         );
3077         mesh.changing(true);
3078     }
3079
3080     // Clear out primitives
3081     {
3082         retiredPoints_.clearStorage();
3083         region_.clearStorage();
3084     }
3085
3086
3087     if (debug)
3088     {
3089         // Some stats on changes
3090         label nAdd, nInflate, nMerge, nRemove;
3091         countMap(pointMap_, reversePointMap_, nAdd, nInflate, nMerge, nRemove);
3092         Pout<< "Points:"
3093             << "  added(from point):" << nAdd
3094             << "  added(from nothing):" << nInflate
3095             << "  merged(into other point):" << nMerge
3096             << "  removed:" << nRemove
3097             << nl;
3098
3099         countMap(faceMap_, reverseFaceMap_, nAdd, nInflate, nMerge, nRemove);
3100         Pout<< "Faces:"
3101             << "  added(from face):" << nAdd
3102             << "  added(inflated):" << nInflate
3103             << "  merged(into other face):" << nMerge
3104             << "  removed:" << nRemove
3105             << nl;
3106
3107         countMap(cellMap_, reverseCellMap_, nAdd, nInflate, nMerge, nRemove);
3108         Pout<< "Cells:"
3109             << "  added(from cell):" << nAdd
3110             << "  added(inflated):" << nInflate
3111             << "  merged(into other cell):" << nMerge
3112             << "  removed:" << nRemove
3113             << nl
3114             << endl;
3115     }
3116
3117     if (debug)
3118     {
3119         Pout<< "New mesh:" << nl;
3120         writeMeshStats(mesh, Pout);
3121     }
3122
3123
3124     // Zones
3125     // ~~~~~
3126
3127     // Inverse of point/face/cell zone addressing.
3128     // For every preserved point/face/cells in zone give the old position.
3129     // For added points, the index is set to -1
3130     labelListList pointZoneMap(mesh.pointZones().size());
3131     labelListList faceZoneFaceMap(mesh.faceZones().size());
3132     labelListList cellZoneMap(mesh.cellZones().size());
3133
3134     resetZones(mesh, mesh, pointZoneMap, faceZoneFaceMap, cellZoneMap);
3135
3136     // Clear zone info
3137     {
3138         pointZone_.clearStorage();
3139         faceZone_.clearStorage();
3140         faceZoneFlip_.clearStorage();
3141         cellZone_.clearStorage();
3142     }
3143
3144
3145     // Patch point renumbering
3146     // For every preserved point on a patch give the old position.
3147     // For added points, the index is set to -1
3148     labelListList patchPointMap(mesh.boundaryMesh().size());
3149     calcPatchPointMap
3150     (
3151         oldPatchMeshPointMaps,
3152         mesh.boundaryMesh(),
3153         patchPointMap
3154     );
3155
3156     // Create the face zone mesh point renumbering
3157     labelListList faceZonePointMap(mesh.faceZones().size());
3158     calcFaceZonePointMap(mesh, oldFaceZoneMeshPointMaps, faceZonePointMap);
3159
3160     labelHashSet flipFaceFluxSet(getSetIndices(flipFaceFlux_));
3161
3162     return autoPtr<mapPolyMesh>
3163     (
3164         new mapPolyMesh
3165         (
3166             mesh,
3167             nOldPoints,
3168             nOldFaces,
3169             nOldCells,
3170
3171             pointMap_,
3172             pointsFromPoints,
3173
3174             faceMap_,
3175             facesFromPoints,
3176             facesFromEdges,
3177             facesFromFaces,
3178
3179             cellMap_,
3180             cellsFromPoints,
3181             cellsFromEdges,
3182             cellsFromFaces,
3183             cellsFromCells,
3184
3185             reversePointMap_,
3186             reverseFaceMap_,
3187             reverseCellMap_,
3188
3189             flipFaceFluxSet,
3190
3191             patchPointMap,
3192
3193             pointZoneMap,
3194
3195             faceZonePointMap,
3196             faceZoneFaceMap,
3197             cellZoneMap,
3198
3199             newPoints,          // if empty signals no inflation.
3200             oldPatchStarts,
3201             oldPatchNMeshPoints,
3202             true                // steal storage.
3203         )
3204     );
3205
3206     // At this point all member DynamicList (pointMap_, cellMap_ etc.) will
3207     // be invalid.
3208 }
3209
3210
3211 Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::makeMesh
3212 (
3213     autoPtr<fvMesh>& newMeshPtr,
3214     const IOobject& io,
3215     const fvMesh& mesh,
3216     const bool syncParallel,
3217     const bool orderCells,
3218     const bool orderPoints
3219 )
3220 {
3221     if (debug)
3222     {
3223         Pout<< "polyTopoChange::changeMesh"
3224             << "(autoPtr<fvMesh>&, const IOobject&, const fvMesh&"
3225             << ", const bool, const bool, const bool)"
3226             << endl;
3227     }
3228
3229     if (debug)
3230     {
3231         Pout<< "Old mesh:" << nl;
3232         writeMeshStats(mesh, Pout);
3233     }
3234
3235     // new mesh points
3236     pointField newPoints;
3237     // number of internal points
3238     label nInternalPoints;
3239     // patch slicing
3240     labelList patchSizes;
3241     labelList patchStarts;
3242     // inflate maps
3243     List<objectMap> pointsFromPoints;
3244     List<objectMap> facesFromPoints;
3245     List<objectMap> facesFromEdges;
3246     List<objectMap> facesFromFaces;
3247     List<objectMap> cellsFromPoints;
3248     List<objectMap> cellsFromEdges;
3249     List<objectMap> cellsFromFaces;
3250     List<objectMap> cellsFromCells;
3251
3252     // old mesh info
3253     List<Map<label> > oldPatchMeshPointMaps;
3254     labelList oldPatchNMeshPoints;
3255     labelList oldPatchStarts;
3256     List<Map<label> > oldFaceZoneMeshPointMaps;
3257
3258     // Compact, reorder patch faces and calculate mesh/patch maps.
3259     compactAndReorder
3260     (
3261         mesh,
3262         syncParallel,
3263         orderCells,
3264         orderPoints,
3265
3266         nInternalPoints,
3267         newPoints,
3268         patchSizes,
3269         patchStarts,
3270         pointsFromPoints,
3271         facesFromPoints,
3272         facesFromEdges,
3273         facesFromFaces,
3274         cellsFromPoints,
3275         cellsFromEdges,
3276         cellsFromFaces,
3277         cellsFromCells,
3278         oldPatchMeshPointMaps,
3279         oldPatchNMeshPoints,
3280         oldPatchStarts,
3281         oldFaceZoneMeshPointMaps
3282     );
3283
3284     const label nOldPoints(mesh.nPoints());
3285     const label nOldFaces(mesh.nFaces());
3286     const label nOldCells(mesh.nCells());
3287
3288
3289     // Create the mesh
3290     // ~~~~~~~~~~~~~~~
3291
3292     newMeshPtr.reset
3293     (
3294         new fvMesh
3295         (
3296             io,
3297             xferMove(newPoints),
3298             faces_.xfer(),
3299             faceOwner_.xfer(),
3300             faceNeighbour_.xfer()
3301         )
3302     );
3303     fvMesh& newMesh = newMeshPtr();
3304
3305     // Clear out primitives
3306     {
3307         retiredPoints_.clearStorage();
3308         region_.clearStorage();
3309     }
3310
3311
3312     if (debug)
3313     {
3314         // Some stats on changes
3315         label nAdd, nInflate, nMerge, nRemove;
3316         countMap(pointMap_, reversePointMap_, nAdd, nInflate, nMerge, nRemove);
3317         Pout<< "Points:"
3318             << "  added(from point):" << nAdd
3319             << "  added(from nothing):" << nInflate
3320             << "  merged(into other point):" << nMerge
3321             << "  removed:" << nRemove
3322             << nl;
3323
3324         countMap(faceMap_, reverseFaceMap_, nAdd, nInflate, nMerge, nRemove);
3325         Pout<< "Faces:"
3326             << "  added(from face):" << nAdd
3327             << "  added(inflated):" << nInflate
3328             << "  merged(into other face):" << nMerge
3329             << "  removed:" << nRemove
3330             << nl;
3331
3332         countMap(cellMap_, reverseCellMap_, nAdd, nInflate, nMerge, nRemove);
3333         Pout<< "Cells:"
3334             << "  added(from cell):" << nAdd
3335             << "  added(inflated):" << nInflate
3336             << "  merged(into other cell):" << nMerge
3337             << "  removed:" << nRemove
3338             << nl
3339             << endl;
3340     }
3341
3342
3343     {
3344         const polyBoundaryMesh& oldPatches = mesh.boundaryMesh();
3345
3346         List<polyPatch*> newBoundary(oldPatches.size());
3347
3348         forAll(oldPatches, patchI)
3349         {
3350             newBoundary[patchI] = oldPatches[patchI].clone
3351             (
3352                 newMesh.boundaryMesh(),
3353                 patchI,
3354                 patchSizes[patchI],
3355                 patchStarts[patchI]
3356             ).ptr();
3357         }
3358         newMesh.addFvPatches(newBoundary);
3359     }
3360
3361
3362     // Zones
3363     // ~~~~~
3364
3365     // Start off from empty zones.
3366     const pointZoneMesh& oldPointZones = mesh.pointZones();
3367     List<pointZone*> pZonePtrs(oldPointZones.size());
3368     {
3369         forAll(oldPointZones, i)
3370         {
3371             pZonePtrs[i] = new pointZone
3372             (
3373                 oldPointZones[i].name(),
3374                 labelList(0),
3375                 i,
3376                 newMesh.pointZones()
3377             );
3378         }
3379     }
3380
3381     const faceZoneMesh& oldFaceZones = mesh.faceZones();
3382     List<faceZone*> fZonePtrs(oldFaceZones.size());
3383     {
3384         forAll(oldFaceZones, i)
3385         {
3386             fZonePtrs[i] = new faceZone
3387             (
3388                 oldFaceZones[i].name(),
3389                 labelList(0),
3390                 boolList(0),
3391                 i,
3392                 newMesh.faceZones()
3393             );
3394         }
3395     }
3396
3397     const cellZoneMesh& oldCellZones = mesh.cellZones();
3398     List<cellZone*> cZonePtrs(oldCellZones.size());
3399     {
3400         forAll(oldCellZones, i)
3401         {
3402             cZonePtrs[i] = new cellZone
3403             (
3404                 oldCellZones[i].name(),
3405                 labelList(0),
3406                 i,
3407                 newMesh.cellZones()
3408             );
3409         }
3410     }
3411
3412     newMesh.addZones(pZonePtrs, fZonePtrs, cZonePtrs);
3413
3414     // Inverse of point/face/cell zone addressing.
3415     // For every preserved point/face/cells in zone give the old position.
3416     // For added points, the index is set to -1
3417     labelListList pointZoneMap(mesh.pointZones().size());
3418     labelListList faceZoneFaceMap(mesh.faceZones().size());
3419     labelListList cellZoneMap(mesh.cellZones().size());
3420
3421     resetZones(mesh, newMesh, pointZoneMap, faceZoneFaceMap, cellZoneMap);
3422
3423     // Clear zone info
3424     {
3425         pointZone_.clearStorage();
3426         faceZone_.clearStorage();
3427         faceZoneFlip_.clearStorage();
3428         cellZone_.clearStorage();
3429     }
3430
3431     // Patch point renumbering
3432     // For every preserved point on a patch give the old position.
3433     // For added points, the index is set to -1
3434     labelListList patchPointMap(newMesh.boundaryMesh().size());
3435     calcPatchPointMap
3436     (
3437         oldPatchMeshPointMaps,
3438         newMesh.boundaryMesh(),
3439         patchPointMap
3440     );
3441
3442     // Create the face zone mesh point renumbering
3443     labelListList faceZonePointMap(newMesh.faceZones().size());
3444     calcFaceZonePointMap(newMesh, oldFaceZoneMeshPointMaps, faceZonePointMap);
3445
3446     if (debug)
3447     {
3448         Pout<< "New mesh:" << nl;
3449         writeMeshStats(mesh, Pout);
3450     }
3451
3452     labelHashSet flipFaceFluxSet(getSetIndices(flipFaceFlux_));
3453
3454     return autoPtr<mapPolyMesh>
3455     (
3456         new mapPolyMesh
3457         (
3458             newMesh,
3459             nOldPoints,
3460             nOldFaces,
3461             nOldCells,
3462
3463             pointMap_,
3464             pointsFromPoints,
3465
3466             faceMap_,
3467             facesFromPoints,
3468             facesFromEdges,
3469             facesFromFaces,
3470
3471             cellMap_,
3472             cellsFromPoints,
3473             cellsFromEdges,
3474             cellsFromFaces,
3475             cellsFromCells,
3476
3477             reversePointMap_,
3478             reverseFaceMap_,
3479             reverseCellMap_,
3480
3481             flipFaceFluxSet,
3482
3483             patchPointMap,
3484
3485             pointZoneMap,
3486
3487             faceZonePointMap,
3488             faceZoneFaceMap,
3489             cellZoneMap,
3490
3491             newPoints,          // if empty signals no inflation.
3492             oldPatchStarts,
3493             oldPatchNMeshPoints,
3494             true                // steal storage.
3495         )
3496     );
3497
3498     // At this point all member DynamicList (pointMap_, cellMap_ etc.) will
3499     // be invalid.
3500 }
3501
3502
3503 // ************************************************************************* //