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         pointZone_.resize(pointZone_.size() + points.size()/100);
2215
2216         // Precalc offset zones
2217         labelList newZoneID(points.size(), -1);
2218
2219         forAll(pointZones, zoneI)
2220         {
2221             const labelList& pointLabels = pointZones[zoneI];
2222
2223             forAll(pointLabels, j)
2224             {
2225                 newZoneID[pointLabels[j]] = pointZoneMap[zoneI];
2226             }
2227         }
2228
2229         // Add points in mesh order
2230         for (label pointI = 0; pointI < mesh.nPoints(); pointI++)
2231         {
2232             addPoint
2233             (
2234                 points[pointI],
2235                 pointI,
2236                 newZoneID[pointI],
2237                 true
2238             );
2239         }
2240     }
2241
2242     // Add cells
2243     {
2244         const cellZoneMesh& cellZones = mesh.cellZones();
2245
2246         // Resize
2247
2248         // Note: polyMesh does not allow retired cells anymore. So allCells
2249         // always equals nCells
2250         label nAllCells = mesh.nCells();
2251
2252         cellMap_.setCapacity(cellMap_.size() + nAllCells);
2253         cellFromPoint_.resize(cellFromPoint_.size() + nAllCells/100);
2254         cellFromEdge_.resize(cellFromEdge_.size() + nAllCells/100);
2255         cellFromFace_.resize(cellFromFace_.size() + nAllCells/100);
2256         cellZone_.setCapacity(cellZone_.size() + nAllCells);
2257
2258
2259         // Precalc offset zones
2260         labelList newZoneID(nAllCells, -1);
2261
2262         forAll(cellZones, zoneI)
2263         {
2264             const labelList& cellLabels = cellZones[zoneI];
2265
2266             forAll(cellLabels, j)
2267             {
2268                 label cellI = cellLabels[j];
2269
2270                 if (newZoneID[cellI] != -1)
2271                 {
2272                     WarningIn
2273                     (
2274                         "polyTopoChange::addMesh"
2275                         "(const polyMesh&, const labelList&,"
2276                         "const labelList&, const labelList&,"
2277                         "const labelList&)"
2278                     )   << "Cell:" << cellI
2279                         << " centre:" << mesh.cellCentres()[cellI]
2280                         << " is in two zones:"
2281                         << cellZones[newZoneID[cellI]].name()
2282                         << " and " << cellZones[zoneI].name() << endl
2283                         << "    This is not supported."
2284                         << " Continuing with first zone only." << endl;
2285                 }
2286                 else
2287                 {
2288                     newZoneID[cellI] = cellZoneMap[zoneI];
2289                 }
2290             }
2291         }
2292
2293         // Add cells in mesh order
2294         for (label cellI = 0; cellI < nAllCells; cellI++)
2295         {
2296             // Add cell from cell
2297             addCell(-1, -1, -1, cellI, newZoneID[cellI]);
2298         }
2299     }
2300
2301     // Add faces
2302     {
2303         const polyBoundaryMesh& patches = mesh.boundaryMesh();
2304         const faceList& faces = mesh.faces();
2305         const labelList& faceOwner = mesh.faceOwner();
2306         const labelList& faceNeighbour = mesh.faceNeighbour();
2307         const faceZoneMesh& faceZones = mesh.faceZones();
2308
2309         // Resize
2310         label nAllFaces = mesh.faces().size();
2311
2312         faces_.setCapacity(faces_.size() + nAllFaces);
2313         region_.setCapacity(region_.size() + nAllFaces);
2314         faceOwner_.setCapacity(faceOwner_.size() + nAllFaces);
2315         faceNeighbour_.setCapacity(faceNeighbour_.size() + nAllFaces);
2316         faceMap_.setCapacity(faceMap_.size() + nAllFaces);
2317         faceFromPoint_.resize(faceFromPoint_.size() + nAllFaces/100);
2318         faceFromEdge_.resize(faceFromEdge_.size() + nAllFaces/100);
2319         flipFaceFlux_.setCapacity(faces_.size() + nAllFaces);
2320         faceZone_.resize(faceZone_.size() + nAllFaces/100);
2321         faceZoneFlip_.setCapacity(faces_.size() + nAllFaces);
2322
2323
2324         // Precalc offset zones
2325         labelList newZoneID(nAllFaces, -1);
2326         boolList zoneFlip(nAllFaces, false);
2327
2328         forAll(faceZones, zoneI)
2329         {
2330             const labelList& faceLabels = faceZones[zoneI];
2331             const boolList& flipMap = faceZones[zoneI].flipMap();
2332
2333             forAll(faceLabels, j)
2334             {
2335                 newZoneID[faceLabels[j]] = faceZoneMap[zoneI];
2336                 zoneFlip[faceLabels[j]] = flipMap[j];
2337             }
2338         }
2339
2340         // Add faces in mesh order
2341
2342         // 1. Internal faces
2343         for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
2344         {
2345             addFace
2346             (
2347                 faces[faceI],
2348                 faceOwner[faceI],
2349                 faceNeighbour[faceI],
2350                 -1,                         // masterPointID
2351                 -1,                         // masterEdgeID
2352                 faceI,                      // masterFaceID
2353                 false,                      // flipFaceFlux
2354                 -1,                         // patchID
2355                 newZoneID[faceI],           // zoneID
2356                 zoneFlip[faceI]             // zoneFlip
2357             );
2358         }
2359
2360         // 2. Patch faces
2361         forAll(patches, patchI)
2362         {
2363             const polyPatch& pp = patches[patchI];
2364
2365             if (pp.start() != faces_.size())
2366             {
2367                 FatalErrorIn
2368                 (
2369                     "polyTopoChange::polyTopoChange"
2370                     "(const polyMesh& mesh, const bool strict)"
2371                 )   << "Problem : "
2372                     << "Patch " << pp.name() << " starts at " << pp.start()
2373                     << endl
2374                     << "Current face counter at " << faces_.size() << endl
2375                     << "Are patches in incremental order?"
2376                     << abort(FatalError);
2377             }
2378             forAll(pp, patchFaceI)
2379             {
2380                 label faceI = pp.start() + patchFaceI;
2381
2382                 addFace
2383                 (
2384                     faces[faceI],
2385                     faceOwner[faceI],
2386                     -1,                         // neighbour
2387                     -1,                         // masterPointID
2388                     -1,                         // masterEdgeID
2389                     faceI,                      // masterFaceID
2390                     false,                      // flipFaceFlux
2391                     patchMap[patchI],           // patchID
2392                     newZoneID[faceI],           // zoneID
2393                     zoneFlip[faceI]             // zoneFlip
2394                 );
2395             }
2396         }
2397     }
2398 }
2399
2400
2401 Foam::label Foam::polyTopoChange::setAction(const topoAction& action)
2402 {
2403     if (isType<polyAddPoint>(action))
2404     {
2405         const polyAddPoint& pap = refCast<const polyAddPoint>(action);
2406
2407         return addPoint
2408         (
2409             pap.newPoint(),
2410             pap.masterPointID(),
2411             pap.zoneID(),
2412             pap.inCell()
2413         );
2414     }
2415     else if (isType<polyModifyPoint>(action))
2416     {
2417         const polyModifyPoint& pmp = refCast<const polyModifyPoint>(action);
2418
2419         modifyPoint
2420         (
2421             pmp.pointID(),
2422             pmp.newPoint(),
2423             pmp.zoneID(),
2424             pmp.inCell()
2425         );
2426
2427         return -1;
2428     }
2429     else if (isType<polyRemovePoint>(action))
2430     {
2431         const polyRemovePoint& prp = refCast<const polyRemovePoint>(action);
2432
2433         removePoint(prp.pointID(), prp.mergePointID());
2434
2435         return -1;
2436     }
2437     else if (isType<polyAddFace>(action))
2438     {
2439         const polyAddFace& paf = refCast<const polyAddFace>(action);
2440
2441         return addFace
2442         (
2443             paf.newFace(),
2444             paf.owner(),
2445             paf.neighbour(),
2446             paf.masterPointID(),
2447             paf.masterEdgeID(),
2448             paf.masterFaceID(),
2449             paf.flipFaceFlux(),
2450             paf.patchID(),
2451             paf.zoneID(),
2452             paf.zoneFlip()
2453         );
2454     }
2455     else if (isType<polyModifyFace>(action))
2456     {
2457         const polyModifyFace& pmf = refCast<const polyModifyFace>(action);
2458
2459         modifyFace
2460         (
2461             pmf.newFace(),
2462             pmf.faceID(),
2463             pmf.owner(),
2464             pmf.neighbour(),
2465             pmf.flipFaceFlux(),
2466             pmf.patchID(),
2467             pmf.zoneID(),
2468             pmf.zoneFlip()
2469         );
2470
2471         return -1;
2472     }
2473     else if (isType<polyRemoveFace>(action))
2474     {
2475         const polyRemoveFace& prf = refCast<const polyRemoveFace>(action);
2476
2477         removeFace(prf.faceID(), prf.mergeFaceID());
2478
2479         return -1;
2480     }
2481     else if (isType<polyAddCell>(action))
2482     {
2483         const polyAddCell& pac = refCast<const polyAddCell>(action);
2484
2485         return addCell
2486         (
2487             pac.masterPointID(),
2488             pac.masterEdgeID(),
2489             pac.masterFaceID(),
2490             pac.masterCellID(),
2491             pac.zoneID()
2492         );
2493     }
2494     else if (isType<polyModifyCell>(action))
2495     {
2496         const polyModifyCell& pmc = refCast<const polyModifyCell>(action);
2497
2498         if (pmc.removeFromZone())
2499         {
2500             modifyCell(pmc.cellID(), -1);
2501         }
2502         else
2503         {
2504             modifyCell(pmc.cellID(), pmc.zoneID());
2505         }
2506
2507         return -1;
2508     }
2509     else if (isType<polyRemoveCell>(action))
2510     {
2511         const polyRemoveCell& prc = refCast<const polyRemoveCell>(action);
2512
2513         removeCell(prc.cellID(), prc.mergeCellID());
2514
2515         return -1;
2516     }
2517     else
2518     {
2519         FatalErrorIn
2520         (
2521             "label polyTopoChange::setAction(const topoAction& action)"
2522         )   << "Unknown type of topoChange: " << action.type()
2523             << abort(FatalError);
2524
2525         // Dummy return to keep compiler happy
2526         return -1;
2527     }
2528 }
2529
2530
2531 Foam::label Foam::polyTopoChange::addPoint
2532 (
2533     const point& pt,
2534     const label masterPointID,
2535     const label zoneID,
2536     const bool inCell
2537 )
2538 {
2539     label pointI = points_.size();
2540
2541     points_.append(pt);
2542     pointMap_.append(masterPointID);
2543     reversePointMap_.append(pointI);
2544
2545     if (zoneID >= 0)
2546     {
2547         pointZone_.insert(pointI, zoneID);
2548     }
2549
2550     if (!inCell)
2551     {
2552         retiredPoints_.insert(pointI);
2553     }
2554
2555     return pointI;
2556 }
2557
2558
2559 void Foam::polyTopoChange::modifyPoint
2560 (
2561     const label pointI,
2562     const point& pt,
2563     const label newZoneID,
2564     const bool inCell
2565 )
2566 {
2567     if (pointI < 0 || pointI >= points_.size())
2568     {
2569         FatalErrorIn
2570         (
2571             "polyTopoChange::modifyPoint(const label, const point&)"
2572         )   << "illegal point label " << pointI << endl
2573             << "Valid point labels are 0 .. " << points_.size()-1
2574             << abort(FatalError);
2575     }
2576     if (pointRemoved(pointI) || pointMap_[pointI] == -1)
2577     {
2578         FatalErrorIn
2579         (
2580             "polyTopoChange::modifyPoint(const label, const point&)"
2581         )   << "point " << pointI << " already marked for removal"
2582             << abort(FatalError);
2583     }
2584     points_[pointI] = pt;
2585
2586     Map<label>::iterator pointFnd = pointZone_.find(pointI);
2587
2588     if (pointFnd != pointZone_.end())
2589     {
2590         if (newZoneID >= 0)
2591         {
2592             pointFnd() = newZoneID;
2593         }
2594         else
2595         {
2596             pointZone_.erase(pointFnd);
2597         }
2598     }
2599     else if (newZoneID >= 0)
2600     {
2601         pointZone_.insert(pointI, newZoneID);
2602     }
2603
2604     if (inCell)
2605     {
2606         retiredPoints_.erase(pointI);
2607     }
2608     else
2609     {
2610         retiredPoints_.insert(pointI);
2611     }
2612 }
2613
2614
2615 void Foam::polyTopoChange::movePoints(const pointField& newPoints)
2616 {
2617     if (newPoints.size() != points_.size())
2618     {
2619         FatalErrorIn("polyTopoChange::movePoints(const pointField&)")
2620             << "illegal pointField size." << endl
2621             << "Size:" << newPoints.size() << endl
2622             << "Points in mesh:" << points_.size()
2623             << abort(FatalError);
2624     }
2625
2626     forAll(points_, pointI)
2627     {
2628         points_[pointI] = newPoints[pointI];
2629     }
2630 }
2631
2632
2633 void Foam::polyTopoChange::removePoint
2634 (
2635     const label pointI,
2636     const label mergePointI
2637 )
2638 {
2639     if (pointI < 0 || pointI >= points_.size())
2640     {
2641         FatalErrorIn("polyTopoChange::removePoint(const label, const label)")
2642             << "illegal point label " << pointI << endl
2643             << "Valid point labels are 0 .. " << points_.size()-1
2644             << abort(FatalError);
2645     }
2646
2647     if
2648     (
2649         strict_
2650      && (pointRemoved(pointI) || pointMap_[pointI] == -1)
2651     )
2652     {
2653         FatalErrorIn("polyTopoChange::removePoint(const label, const label)")
2654             << "point " << pointI << " already marked for removal" << nl
2655             << "Point:" << points_[pointI] << " pointMap:" << pointMap_[pointI]
2656             << abort(FatalError);
2657     }
2658
2659     if (pointI == mergePointI)
2660     {
2661         FatalErrorIn("polyTopoChange::removePoint(const label, const label)")
2662             << "Cannot remove/merge point " << pointI << " onto itself."
2663             << abort(FatalError);
2664     }
2665
2666     points_[pointI] = greatPoint;
2667     pointMap_[pointI] = -1;
2668     if (mergePointI >= 0)
2669     {
2670         reversePointMap_[pointI] = -mergePointI-2;
2671     }
2672     else
2673     {
2674         reversePointMap_[pointI] = -1;
2675     }
2676     pointZone_.erase(pointI);
2677     retiredPoints_.erase(pointI);
2678 }
2679
2680
2681 Foam::label Foam::polyTopoChange::addFace
2682 (
2683     const face& f,
2684     const label own,
2685     const label nei,
2686     const label masterPointID,
2687     const label masterEdgeID,
2688     const label masterFaceID,
2689     const bool flipFaceFlux,
2690     const label patchID,
2691     const label zoneID,
2692     const bool zoneFlip
2693 )
2694 {
2695     // Check validity
2696     if (debug)
2697     {
2698         checkFace(f, -1, own, nei, patchID, zoneID);
2699     }
2700
2701     label faceI = faces_.size();
2702
2703     faces_.append(f);
2704     region_.append(patchID);
2705     faceOwner_.append(own);
2706     faceNeighbour_.append(nei);
2707
2708     if (masterPointID >= 0)
2709     {
2710         faceMap_.append(-1);
2711         faceFromPoint_.insert(faceI, masterPointID);
2712     }
2713     else if (masterEdgeID >= 0)
2714     {
2715         faceMap_.append(-1);
2716         faceFromEdge_.insert(faceI, masterEdgeID);
2717     }
2718     else if (masterFaceID >= 0)
2719     {
2720         faceMap_.append(masterFaceID);
2721     }
2722     else
2723     {
2724         // Allow inflate-from-nothing?
2725         //FatalErrorIn("polyTopoChange::addFace")
2726         //    << "Need to specify a master point, edge or face"
2727         //    << "face:" << f << " own:" << own << " nei:" << nei
2728         //    << abort(FatalError);
2729         faceMap_.append(-1);
2730     }
2731     reverseFaceMap_.append(faceI);
2732
2733     flipFaceFlux_[faceI] = (flipFaceFlux ? 1 : 0);
2734
2735     if (zoneID >= 0)
2736     {
2737         faceZone_.insert(faceI, zoneID);
2738     }
2739     faceZoneFlip_[faceI] = (zoneFlip ? 1 : 0);
2740
2741     return faceI;
2742 }
2743
2744
2745 void Foam::polyTopoChange::modifyFace
2746 (
2747     const face& f,
2748     const label faceI,
2749     const label own,
2750     const label nei,
2751     const bool flipFaceFlux,
2752     const label patchID,
2753     const label zoneID,
2754     const bool zoneFlip
2755 )
2756 {
2757     // Check validity
2758     if (debug)
2759     {
2760         checkFace(f, faceI, own, nei, patchID, zoneID);
2761     }
2762
2763     faces_[faceI] = f;
2764     faceOwner_[faceI] = own;
2765     faceNeighbour_[faceI] = nei;
2766     region_[faceI] = patchID;
2767
2768     flipFaceFlux_[faceI] = (flipFaceFlux ? 1 : 0);
2769
2770     Map<label>::iterator faceFnd = faceZone_.find(faceI);
2771
2772     if (faceFnd != faceZone_.end())
2773     {
2774         if (zoneID >= 0)
2775         {
2776             faceFnd() = zoneID;
2777         }
2778         else
2779         {
2780             faceZone_.erase(faceFnd);
2781         }
2782     }
2783     else if (zoneID >= 0)
2784     {
2785         faceZone_.insert(faceI, zoneID);
2786     }
2787     faceZoneFlip_[faceI] = (zoneFlip ? 1 : 0);
2788 }
2789
2790
2791 void Foam::polyTopoChange::removeFace(const label faceI, const label mergeFaceI)
2792 {
2793     if (faceI < 0 || faceI >= faces_.size())
2794     {
2795         FatalErrorIn("polyTopoChange::removeFace(const label, const label)")
2796             << "illegal face label " << faceI << endl
2797             << "Valid face labels are 0 .. " << faces_.size()-1
2798             << abort(FatalError);
2799     }
2800
2801     if
2802     (
2803         strict_
2804      && (faceRemoved(faceI) || faceMap_[faceI] == -1)
2805     )
2806     {
2807         FatalErrorIn("polyTopoChange::removeFace(const label, const label)")
2808             << "face " << faceI
2809             << " already marked for removal"
2810             << abort(FatalError);
2811     }
2812
2813     faces_[faceI].setSize(0);
2814     region_[faceI] = -1;
2815     faceOwner_[faceI] = -1;
2816     faceNeighbour_[faceI] = -1;
2817     faceMap_[faceI] = -1;
2818     if (mergeFaceI >= 0)
2819     {
2820         reverseFaceMap_[faceI] = -mergeFaceI-2;
2821     }
2822     else
2823     {
2824         reverseFaceMap_[faceI] = -1;
2825     }
2826     faceFromEdge_.erase(faceI);
2827     faceFromPoint_.erase(faceI);
2828     flipFaceFlux_[faceI] = 0;
2829     faceZone_.erase(faceI);
2830     faceZoneFlip_[faceI] = 0;
2831 }
2832
2833
2834 Foam::label Foam::polyTopoChange::addCell
2835 (
2836     const label masterPointID,
2837     const label masterEdgeID,
2838     const label masterFaceID,
2839     const label masterCellID,
2840     const label zoneID
2841 )
2842 {
2843     label cellI = cellMap_.size();
2844
2845     if (masterPointID >= 0)
2846     {
2847         cellMap_.append(-1);
2848         cellFromPoint_.insert(cellI, masterPointID);
2849     }
2850     else if (masterEdgeID >= 0)
2851     {
2852         cellMap_.append(-1);
2853         cellFromEdge_.insert(cellI, masterEdgeID);
2854     }
2855     else if (masterFaceID >= 0)
2856     {
2857         cellMap_.append(-1);
2858         cellFromFace_.insert(cellI, masterFaceID);
2859     }
2860     else
2861     {
2862         cellMap_.append(masterCellID);
2863     }
2864     reverseCellMap_.append(cellI);
2865     cellZone_.append(zoneID);
2866
2867     return cellI;
2868 }
2869
2870
2871 void Foam::polyTopoChange::modifyCell
2872 (
2873     const label cellI,
2874     const label zoneID
2875 )
2876 {
2877     cellZone_[cellI] = zoneID;
2878 }
2879
2880
2881 void Foam::polyTopoChange::removeCell(const label cellI, const label mergeCellI)
2882 {
2883     if (cellI < 0 || cellI >= cellMap_.size())
2884     {
2885         FatalErrorIn("polyTopoChange::removeCell(const label, const label)")
2886             << "illegal cell label " << cellI << endl
2887             << "Valid cell labels are 0 .. " << cellMap_.size()-1
2888             << abort(FatalError);
2889     }
2890
2891     if (strict_ && cellMap_[cellI] == -2)
2892     {
2893         FatalErrorIn("polyTopoChange::removeCell(const label, const label)")
2894             << "cell " << cellI
2895             << " already marked for removal"
2896             << abort(FatalError);
2897     }
2898
2899     cellMap_[cellI] = -2;
2900     if (mergeCellI >= 0)
2901     {
2902         reverseCellMap_[cellI] = -mergeCellI-2;
2903     }
2904     else
2905     {
2906         reverseCellMap_[cellI] = -1;
2907     }
2908     cellFromPoint_.erase(cellI);
2909     cellFromEdge_.erase(cellI);
2910     cellFromFace_.erase(cellI);
2911     cellZone_[cellI] = -1;
2912 }
2913
2914
2915 Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::changeMesh
2916 (
2917     polyMesh& mesh,
2918     const bool inflate,
2919     const bool syncParallel,
2920     const bool orderCells,
2921     const bool orderPoints
2922 )
2923 {
2924     if (debug)
2925     {
2926         Pout<< "polyTopoChange::changeMesh"
2927             << "(polyMesh&, const bool, const bool, const bool, const bool)"
2928             << endl;
2929     }
2930
2931     if (debug)
2932     {
2933         Pout<< "Old mesh:" << nl;
2934         writeMeshStats(mesh, Pout);
2935     }
2936
2937     // new mesh points
2938     pointField newPoints;
2939     // number of internal points
2940     label nInternalPoints;
2941     // patch slicing
2942     labelList patchSizes;
2943     labelList patchStarts;
2944     // inflate maps
2945     List<objectMap> pointsFromPoints;
2946     List<objectMap> facesFromPoints;
2947     List<objectMap> facesFromEdges;
2948     List<objectMap> facesFromFaces;
2949     List<objectMap> cellsFromPoints;
2950     List<objectMap> cellsFromEdges;
2951     List<objectMap> cellsFromFaces;
2952     List<objectMap> cellsFromCells;
2953     // old mesh info
2954     List<Map<label> > oldPatchMeshPointMaps;
2955     labelList oldPatchNMeshPoints;
2956     labelList oldPatchStarts;
2957     List<Map<label> > oldFaceZoneMeshPointMaps;
2958
2959     // Compact, reorder patch faces and calculate mesh/patch maps.
2960     compactAndReorder
2961     (
2962         mesh,
2963         syncParallel,
2964         orderCells,
2965         orderPoints,
2966
2967         nInternalPoints,
2968         newPoints,
2969         patchSizes,
2970         patchStarts,
2971         pointsFromPoints,
2972         facesFromPoints,
2973         facesFromEdges,
2974         facesFromFaces,
2975         cellsFromPoints,
2976         cellsFromEdges,
2977         cellsFromFaces,
2978         cellsFromCells,
2979         oldPatchMeshPointMaps,
2980         oldPatchNMeshPoints,
2981         oldPatchStarts,
2982         oldFaceZoneMeshPointMaps
2983     );
2984
2985     const label nOldPoints(mesh.nPoints());
2986     const label nOldFaces(mesh.nFaces());
2987     const label nOldCells(mesh.nCells());
2988
2989
2990     // Change the mesh
2991     // ~~~~~~~~~~~~~~~
2992     // This will invalidate any addressing so better make sure you have
2993     // all the information you need!!!
2994
2995     if (inflate)
2996     {
2997         // Keep (renumbered) mesh points, store new points in map for inflation
2998         // (appended points (i.e. from nowhere) get value zero)
2999         pointField renumberedMeshPoints(newPoints.size());
3000
3001         forAll(pointMap_, newPointI)
3002         {
3003             label oldPointI = pointMap_[newPointI];
3004
3005             if (oldPointI >= 0)
3006             {
3007                 renumberedMeshPoints[newPointI] = mesh.points()[oldPointI];
3008             }
3009             else
3010             {
3011                 renumberedMeshPoints[newPointI] = vector::zero;
3012             }
3013         }
3014
3015         mesh.resetPrimitives
3016         (
3017             xferMove(renumberedMeshPoints),
3018             faces_.xfer(),
3019             faceOwner_.xfer(),
3020             faceNeighbour_.xfer(),
3021             patchSizes,
3022             patchStarts,
3023             syncParallel
3024         );
3025
3026         mesh.changing(true);
3027     }
3028     else
3029     {
3030         // Set new points.
3031         mesh.resetPrimitives
3032         (
3033             xferMove(newPoints),
3034             faces_.xfer(),
3035             faceOwner_.xfer(),
3036             faceNeighbour_.xfer(),
3037             patchSizes,
3038             patchStarts,
3039             syncParallel
3040         );
3041         mesh.changing(true);
3042     }
3043
3044     // Clear out primitives
3045     {
3046         retiredPoints_.clearStorage();
3047         region_.clearStorage();
3048     }
3049
3050
3051     if (debug)
3052     {
3053         // Some stats on changes
3054         label nAdd, nInflate, nMerge, nRemove;
3055         countMap(pointMap_, reversePointMap_, nAdd, nInflate, nMerge, nRemove);
3056         Pout<< "Points:"
3057             << "  added(from point):" << nAdd
3058             << "  added(from nothing):" << nInflate
3059             << "  merged(into other point):" << nMerge
3060             << "  removed:" << nRemove
3061             << nl;
3062
3063         countMap(faceMap_, reverseFaceMap_, nAdd, nInflate, nMerge, nRemove);
3064         Pout<< "Faces:"
3065             << "  added(from face):" << nAdd
3066             << "  added(inflated):" << nInflate
3067             << "  merged(into other face):" << nMerge
3068             << "  removed:" << nRemove
3069             << nl;
3070
3071         countMap(cellMap_, reverseCellMap_, nAdd, nInflate, nMerge, nRemove);
3072         Pout<< "Cells:"
3073             << "  added(from cell):" << nAdd
3074             << "  added(inflated):" << nInflate
3075             << "  merged(into other cell):" << nMerge
3076             << "  removed:" << nRemove
3077             << nl
3078             << endl;
3079     }
3080
3081     if (debug)
3082     {
3083         Pout<< "New mesh:" << nl;
3084         writeMeshStats(mesh, Pout);
3085     }
3086
3087
3088     // Zones
3089     // ~~~~~
3090
3091     // Inverse of point/face/cell zone addressing.
3092     // For every preserved point/face/cells in zone give the old position.
3093     // For added points, the index is set to -1
3094     labelListList pointZoneMap(mesh.pointZones().size());
3095     labelListList faceZoneFaceMap(mesh.faceZones().size());
3096     labelListList cellZoneMap(mesh.cellZones().size());
3097
3098     resetZones(mesh, mesh, pointZoneMap, faceZoneFaceMap, cellZoneMap);
3099
3100     // Clear zone info
3101     {
3102         pointZone_.clearStorage();
3103         faceZone_.clearStorage();
3104         faceZoneFlip_.clearStorage();
3105         cellZone_.clearStorage();
3106     }
3107
3108
3109     // Patch point renumbering
3110     // For every preserved point on a patch give the old position.
3111     // For added points, the index is set to -1
3112     labelListList patchPointMap(mesh.boundaryMesh().size());
3113     calcPatchPointMap
3114     (
3115         oldPatchMeshPointMaps,
3116         mesh.boundaryMesh(),
3117         patchPointMap
3118     );
3119
3120     // Create the face zone mesh point renumbering
3121     labelListList faceZonePointMap(mesh.faceZones().size());
3122     calcFaceZonePointMap(mesh, oldFaceZoneMeshPointMaps, faceZonePointMap);
3123
3124     labelHashSet flipFaceFluxSet(getSetIndices(flipFaceFlux_));
3125
3126     return autoPtr<mapPolyMesh>
3127     (
3128         new mapPolyMesh
3129         (
3130             mesh,
3131             nOldPoints,
3132             nOldFaces,
3133             nOldCells,
3134
3135             pointMap_,
3136             pointsFromPoints,
3137
3138             faceMap_,
3139             facesFromPoints,
3140             facesFromEdges,
3141             facesFromFaces,
3142
3143             cellMap_,
3144             cellsFromPoints,
3145             cellsFromEdges,
3146             cellsFromFaces,
3147             cellsFromCells,
3148
3149             reversePointMap_,
3150             reverseFaceMap_,
3151             reverseCellMap_,
3152
3153             flipFaceFluxSet,
3154
3155             patchPointMap,
3156
3157             pointZoneMap,
3158
3159             faceZonePointMap,
3160             faceZoneFaceMap,
3161             cellZoneMap,
3162
3163             newPoints,          // if empty signals no inflation.
3164             oldPatchStarts,
3165             oldPatchNMeshPoints,
3166             true                // steal storage.
3167         )
3168     );
3169
3170     // At this point all member DynamicList (pointMap_, cellMap_ etc.) will
3171     // be invalid.
3172 }
3173
3174
3175 Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::makeMesh
3176 (
3177     autoPtr<fvMesh>& newMeshPtr,
3178     const IOobject& io,
3179     const fvMesh& mesh,
3180     const bool syncParallel,
3181     const bool orderCells,
3182     const bool orderPoints
3183 )
3184 {
3185     if (debug)
3186     {
3187         Pout<< "polyTopoChange::changeMesh"
3188             << "(autoPtr<fvMesh>&, const IOobject&, const fvMesh&"
3189             << ", const bool, const bool, const bool)"
3190             << endl;
3191     }
3192
3193     if (debug)
3194     {
3195         Pout<< "Old mesh:" << nl;
3196         writeMeshStats(mesh, Pout);
3197     }
3198
3199     // new mesh points
3200     pointField newPoints;
3201     // number of internal points
3202     label nInternalPoints;
3203     // patch slicing
3204     labelList patchSizes;
3205     labelList patchStarts;
3206     // inflate maps
3207     List<objectMap> pointsFromPoints;
3208     List<objectMap> facesFromPoints;
3209     List<objectMap> facesFromEdges;
3210     List<objectMap> facesFromFaces;
3211     List<objectMap> cellsFromPoints;
3212     List<objectMap> cellsFromEdges;
3213     List<objectMap> cellsFromFaces;
3214     List<objectMap> cellsFromCells;
3215
3216     // old mesh info
3217     List<Map<label> > oldPatchMeshPointMaps;
3218     labelList oldPatchNMeshPoints;
3219     labelList oldPatchStarts;
3220     List<Map<label> > oldFaceZoneMeshPointMaps;
3221
3222     // Compact, reorder patch faces and calculate mesh/patch maps.
3223     compactAndReorder
3224     (
3225         mesh,
3226         syncParallel,
3227         orderCells,
3228         orderPoints,
3229
3230         nInternalPoints,
3231         newPoints,
3232         patchSizes,
3233         patchStarts,
3234         pointsFromPoints,
3235         facesFromPoints,
3236         facesFromEdges,
3237         facesFromFaces,
3238         cellsFromPoints,
3239         cellsFromEdges,
3240         cellsFromFaces,
3241         cellsFromCells,
3242         oldPatchMeshPointMaps,
3243         oldPatchNMeshPoints,
3244         oldPatchStarts,
3245         oldFaceZoneMeshPointMaps
3246     );
3247
3248     const label nOldPoints(mesh.nPoints());
3249     const label nOldFaces(mesh.nFaces());
3250     const label nOldCells(mesh.nCells());
3251
3252
3253     // Create the mesh
3254     // ~~~~~~~~~~~~~~~
3255
3256     newMeshPtr.reset
3257     (
3258         new fvMesh
3259         (
3260             io,
3261             xferMove(newPoints),
3262             faces_.xfer(),
3263             faceOwner_.xfer(),
3264             faceNeighbour_.xfer()
3265         )
3266     );
3267     fvMesh& newMesh = newMeshPtr();
3268
3269     // Clear out primitives
3270     {
3271         retiredPoints_.clearStorage();
3272         region_.clearStorage();
3273     }
3274
3275
3276     if (debug)
3277     {
3278         // Some stats on changes
3279         label nAdd, nInflate, nMerge, nRemove;
3280         countMap(pointMap_, reversePointMap_, nAdd, nInflate, nMerge, nRemove);
3281         Pout<< "Points:"
3282             << "  added(from point):" << nAdd
3283             << "  added(from nothing):" << nInflate
3284             << "  merged(into other point):" << nMerge
3285             << "  removed:" << nRemove
3286             << nl;
3287
3288         countMap(faceMap_, reverseFaceMap_, nAdd, nInflate, nMerge, nRemove);
3289         Pout<< "Faces:"
3290             << "  added(from face):" << nAdd
3291             << "  added(inflated):" << nInflate
3292             << "  merged(into other face):" << nMerge
3293             << "  removed:" << nRemove
3294             << nl;
3295
3296         countMap(cellMap_, reverseCellMap_, nAdd, nInflate, nMerge, nRemove);
3297         Pout<< "Cells:"
3298             << "  added(from cell):" << nAdd
3299             << "  added(inflated):" << nInflate
3300             << "  merged(into other cell):" << nMerge
3301             << "  removed:" << nRemove
3302             << nl
3303             << endl;
3304     }
3305
3306
3307     {
3308         const polyBoundaryMesh& oldPatches = mesh.boundaryMesh();
3309
3310         List<polyPatch*> newBoundary(oldPatches.size());
3311
3312         forAll(oldPatches, patchI)
3313         {
3314             newBoundary[patchI] = oldPatches[patchI].clone
3315             (
3316                 newMesh.boundaryMesh(),
3317                 patchI,
3318                 patchSizes[patchI],
3319                 patchStarts[patchI]
3320             ).ptr();
3321         }
3322         newMesh.addFvPatches(newBoundary);
3323     }
3324
3325
3326     // Zones
3327     // ~~~~~
3328
3329     // Start off from empty zones.
3330     const pointZoneMesh& oldPointZones = mesh.pointZones();
3331     List<pointZone*> pZonePtrs(oldPointZones.size());
3332     {
3333         forAll(oldPointZones, i)
3334         {
3335             pZonePtrs[i] = new pointZone
3336             (
3337                 oldPointZones[i].name(),
3338                 labelList(0),
3339                 i,
3340                 newMesh.pointZones()
3341             );
3342         }
3343     }
3344
3345     const faceZoneMesh& oldFaceZones = mesh.faceZones();
3346     List<faceZone*> fZonePtrs(oldFaceZones.size());
3347     {
3348         forAll(oldFaceZones, i)
3349         {
3350             fZonePtrs[i] = new faceZone
3351             (
3352                 oldFaceZones[i].name(),
3353                 labelList(0),
3354                 boolList(0),
3355                 i,
3356                 newMesh.faceZones()
3357             );
3358         }
3359     }
3360
3361     const cellZoneMesh& oldCellZones = mesh.cellZones();
3362     List<cellZone*> cZonePtrs(oldCellZones.size());
3363     {
3364         forAll(oldCellZones, i)
3365         {
3366             cZonePtrs[i] = new cellZone
3367             (
3368                 oldCellZones[i].name(),
3369                 labelList(0),
3370                 i,
3371                 newMesh.cellZones()
3372             );
3373         }
3374     }
3375
3376     newMesh.addZones(pZonePtrs, fZonePtrs, cZonePtrs);
3377
3378     // Inverse of point/face/cell zone addressing.
3379     // For every preserved point/face/cells in zone give the old position.
3380     // For added points, the index is set to -1
3381     labelListList pointZoneMap(mesh.pointZones().size());
3382     labelListList faceZoneFaceMap(mesh.faceZones().size());
3383     labelListList cellZoneMap(mesh.cellZones().size());
3384
3385     resetZones(mesh, newMesh, pointZoneMap, faceZoneFaceMap, cellZoneMap);
3386
3387     // Clear zone info
3388     {
3389         pointZone_.clearStorage();
3390         faceZone_.clearStorage();
3391         faceZoneFlip_.clearStorage();
3392         cellZone_.clearStorage();
3393     }
3394
3395     // Patch point renumbering
3396     // For every preserved point on a patch give the old position.
3397     // For added points, the index is set to -1
3398     labelListList patchPointMap(newMesh.boundaryMesh().size());
3399     calcPatchPointMap
3400     (
3401         oldPatchMeshPointMaps,
3402         newMesh.boundaryMesh(),
3403         patchPointMap
3404     );
3405
3406     // Create the face zone mesh point renumbering
3407     labelListList faceZonePointMap(newMesh.faceZones().size());
3408     calcFaceZonePointMap(newMesh, oldFaceZoneMeshPointMaps, faceZonePointMap);
3409
3410     if (debug)
3411     {
3412         Pout<< "New mesh:" << nl;
3413         writeMeshStats(mesh, Pout);
3414     }
3415
3416     labelHashSet flipFaceFluxSet(getSetIndices(flipFaceFlux_));
3417
3418     return autoPtr<mapPolyMesh>
3419     (
3420         new mapPolyMesh
3421         (
3422             newMesh,
3423             nOldPoints,
3424             nOldFaces,
3425             nOldCells,
3426
3427             pointMap_,
3428             pointsFromPoints,
3429
3430             faceMap_,
3431             facesFromPoints,
3432             facesFromEdges,
3433             facesFromFaces,
3434
3435             cellMap_,
3436             cellsFromPoints,
3437             cellsFromEdges,
3438             cellsFromFaces,
3439             cellsFromCells,
3440
3441             reversePointMap_,
3442             reverseFaceMap_,
3443             reverseCellMap_,
3444
3445             flipFaceFluxSet,
3446
3447             patchPointMap,
3448
3449             pointZoneMap,
3450
3451             faceZonePointMap,
3452             faceZoneFaceMap,
3453             cellZoneMap,
3454
3455             newPoints,          // if empty signals no inflation.
3456             oldPatchStarts,
3457             oldPatchNMeshPoints,
3458             true                // steal storage.
3459         )
3460     );
3461
3462     // At this point all member DynamicList (pointMap_, cellMap_ etc.) will
3463     // be invalid.
3464 }
3465
3466
3467 // ************************************************************************* //