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
 788     reorder(oldToNew, faceZone_);
 789     faceZone_.setCapacity(newSize);
 790
 791     inplaceReorder(oldToNew, faceZoneFlip_);
 792     faceZoneFlip_.setCapacity(newSize);
 793 }
 794
 795
 796 // Compact all and orders points and faces:
 797 // - points into internal followed by external points
 798 // - internalfaces upper-triangular
 799 // - externalfaces after internal ones.
 800 void Foam::polyTopoChange::compact
 801 (
 802     const bool orderCells,
 803     const bool orderPoints,
 804     label& nInternalPoints,
 805     labelList& patchSizes,
 806     labelList& patchStarts
 807 )
 808 {
 809     points_.shrink();
 810     pointMap_.shrink();
 811     reversePointMap_.shrink();
 812     pointZone_.shrink();
 813
 814     faces_.shrink();
 815     region_.shrink();
 816     faceOwner_.shrink();
 817     faceNeighbour_.shrink();
 818     faceMap_.shrink();
 819     reverseFaceMap_.shrink();
 820     faceZone_.shrink();
 821
 822     cellMap_.shrink();
 823     reverseCellMap_.shrink();
 824     cellZone_.shrink();
 825
 826
 827     // Compact points
 828     label nActivePoints = 0;
 829     {
 830         labelList localPointMap(points_.size(), -1);
 831         label newPointI = 0;
 832
 833         if (!orderPoints)
 834         {
 835             nInternalPoints = -1;
 836
 837             forAll(points_, pointI)
 838             {
 839                 if (!pointRemoved(pointI) && !retiredPoints_.found(pointI))
 840                 {
 841                     localPointMap[pointI] = newPointI++;
 842                 }
 843             }
 844             nActivePoints = newPointI;
 845         }
 846         else
 847         {
 848             forAll(points_, pointI)
 849             {
 850                 if (!pointRemoved(pointI) && !retiredPoints_.found(pointI))
 851                 {
 852                     nActivePoints++;
 853                 }
 854             }
 855
 856             // Mark boundary points
 857             forAll(faceOwner_, faceI)
 858             {
 859                 if
 860                 (
 861                    !faceRemoved(faceI)
 862                  && faceOwner_[faceI] >= 0
 863                  && faceNeighbour_[faceI] < 0
 864                 )
 865                 {
 866                     // Valid boundary face
 867                     const face& f = faces_[faceI];
 868
 869                     forAll(f, fp)
 870                     {
 871                         label pointI = f[fp];
 872
 873                         if (localPointMap[pointI] == -1)
 874                         {
 875                             if
 876                             (
 877                                 pointRemoved(pointI)
 878                              || retiredPoints_.found(pointI)
 879                             )
 880                             {
 881                                 FatalErrorIn("polyTopoChange::compact(..)")
 882                                     << "Removed or retired point " << pointI
 883                                     << " in face " << f
 884                                     << " at position " << faceI << endl
 885                                     << "Probably face has not been adapted for"
 886                                     << " removed points." << abort(FatalError);
 887                             }
 888                             localPointMap[pointI] = newPointI++;
 889                         }
 890                     }
 891                 }
 892             }
 893
 894             label nBoundaryPoints = newPointI;
 895             nInternalPoints = nActivePoints - nBoundaryPoints;
 896
 897             // Move the boundary addressing up
 898             forAll(localPointMap, pointI)
 899             {
 900                 if (localPointMap[pointI] != -1)
 901                 {
 902                     localPointMap[pointI] += nInternalPoints;
 903                 }
 904             }
 905
 906             newPointI = 0;
 907
 908             // Mark internal points
 909             forAll(faceOwner_, faceI)
 910             {
 911                 if
 912                 (
 913                    !faceRemoved(faceI)
 914                  && faceOwner_[faceI] >= 0
 915                  && faceNeighbour_[faceI] >= 0
 916                 )
 917                 {
 918                     // Valid internal face
 919                     const face& f = faces_[faceI];
 920
 921                     forAll(f, fp)
 922                     {
 923                         label pointI = f[fp];
 924
 925                         if (localPointMap[pointI] == -1)
 926                         {
 927                             if
 928                             (
 929                                 pointRemoved(pointI)
 930                              || retiredPoints_.found(pointI)
 931                             )
 932                             {
 933                                 FatalErrorIn("polyTopoChange::compact(..)")
 934                                     << "Removed or retired point " << pointI
 935                                     << " in face " << f
 936                                     << " at position " << faceI << endl
 937                                     << "Probably face has not been adapted for"
 938                                     << " removed points." << abort(FatalError);
 939                             }
 940                             localPointMap[pointI] = newPointI++;
 941                         }
 942                     }
 943                 }
 944             }
 945
 946             if (newPointI != nInternalPoints)
 947             {
 948                 FatalErrorIn("polyTopoChange::compact(..)")
 949                     << "Problem." << abort(FatalError);
 950             }
 951             newPointI = nActivePoints;
 952         }
 953
 954         forAllConstIter(labelHashSet, retiredPoints_, iter)
 955         {
 956             localPointMap[iter.key()] = newPointI++;
 957         }
 958
 959
 960         if (debug)
 961         {
 962             Pout<< "Points : active:" << nActivePoints
 963                 << "  removed:" << points_.size()-newPointI << endl;
 964         }
 965
 966         reorder(localPointMap, points_);
 967         points_.setCapacity(newPointI);
 968
 969         // Update pointMaps
 970         reorder(localPointMap, pointMap_);
 971         pointMap_.setCapacity(newPointI);
 972         renumberReverseMap(localPointMap, reversePointMap_);
 973
 974         reorder(localPointMap, pointZone_);
 975         pointZone_.setCapacity(newPointI);
 976         renumber(localPointMap, retiredPoints_);
 977
 978         // Use map to relabel face vertices
 979         forAll(faces_, faceI)
 980         {
 981             face& f = faces_[faceI];
 982
 983             //labelList oldF(f);
 984             renumberCompact(localPointMap, f);
 985
 986             if (!faceRemoved(faceI) && f.size() < 3)
 987             {
 988                 FatalErrorIn("polyTopoChange::compact(..)")
 989                     << "Created illegal face " << f
 990                     //<< " from face " << oldF
 991                     << " at position:" << faceI
 992                     << " when filtering removed points"
 993                     << abort(FatalError);
 994             }
 995         }
 996     }
 997
 998
 999     // Compact faces.
1000     {
1001         labelList localFaceMap(faces_.size(), -1);
1002         label newFaceI = 0;
1003
1004         forAll(faces_, faceI)
1005         {
1006             if (!faceRemoved(faceI) && faceOwner_[faceI] >= 0)
1007             {
1008                 localFaceMap[faceI] = newFaceI++;
1009             }
1010         }
1011         nActiveFaces_ = newFaceI;
1012
1013         forAll(faces_, faceI)
1014         {
1015             if (!faceRemoved(faceI) && faceOwner_[faceI] < 0)
1016             {
1017                 // Retired face
1018                 localFaceMap[faceI] = newFaceI++;
1019             }
1020         }
1021
1022         if (debug)
1023         {
1024             Pout<< "Faces : active:" << nActiveFaces_
1025                 << "  removed:" << faces_.size()-newFaceI << endl;
1026         }
1027
1028         // Reorder faces.
1029         reorderCompactFaces(newFaceI, localFaceMap);
1030     }
1031
1032     // Compact cells.
1033     {
1034         labelList localCellMap;
1035         label newCellI;
1036
1037         if (orderCells)
1038         {
1039             // Construct cellCell addressing
1040             CompactListList<label> cellCells;
1041             makeCellCells(nActiveFaces_, cellCells);
1042
1043             // Cell ordering (based on bandCompression). Handles removed cells.
1044             newCellI = getCellOrder(cellCells, localCellMap);
1045         }
1046         else
1047         {
1048             // Compact out removed cells
1049             localCellMap.setSize(cellMap_.size());
1050             localCellMap = -1;
1051
1052             newCellI = 0;
1053             forAll(cellMap_, cellI)
1054             {
1055                 if (!cellRemoved(cellI))
1056                 {
1057                     localCellMap[cellI] = newCellI++;
1058                 }
1059             }
1060         }
1061
1062         if (debug)
1063         {
1064             Pout<< "Cells : active:" << newCellI
1065                 << "  removed:" << cellMap_.size()-newCellI << endl;
1066         }
1067
1068         // Renumber -if cells reordered or -if cells removed
1069         if (orderCells || (newCellI != cellMap_.size()))
1070         {
1071             reorder(localCellMap, cellMap_);
1072             cellMap_.setCapacity(newCellI);
1073             renumberReverseMap(localCellMap, reverseCellMap_);
1074
1075             reorder(localCellMap, cellZone_);
1076             cellZone_.setCapacity(newCellI);
1077
1078             renumberKey(localCellMap, cellFromPoint_);
1079             renumberKey(localCellMap, cellFromEdge_);
1080             renumberKey(localCellMap, cellFromFace_);
1081
1082             // Renumber owner/neighbour. Take into account if neighbour suddenly
1083             // gets lower cell than owner.
1084             forAll(faceOwner_, faceI)
1085             {
1086                 label own = faceOwner_[faceI];
1087                 label nei = faceNeighbour_[faceI];
1088
1089                 if (own >= 0)
1090                 {
1091                     // Update owner
1092                     faceOwner_[faceI] = localCellMap[own];
1093
1094                     if (nei >= 0)
1095                     {
1096                         // Update neighbour.
1097                         faceNeighbour_[faceI] = localCellMap[nei];
1098
1099                         // Check if face needs reversing.
1100                         if
1101                         (
1102                             faceNeighbour_[faceI] >= 0
1103                          && faceNeighbour_[faceI] < faceOwner_[faceI]
1104                         )
1105                         {
1106                             faces_[faceI] = faces_[faceI].reverseFace();
1107                             Swap(faceOwner_[faceI], faceNeighbour_[faceI]);
1108                             flipFaceFlux_[faceI] =
1109                             (
1110                                 flipFaceFlux_[faceI]
1111                               ? 0
1112                               : 1
1113                             );
1114                             faceZoneFlip_[faceI] =
1115                             (
1116                                 faceZoneFlip_[faceI]
1117                               ? 0
1118                               : 1
1119                             );
1120                         }
1121                     }
1122                 }
1123                 else if (nei >= 0)
1124                 {
1125                     // Update neighbour.
1126                     faceNeighbour_[faceI] = localCellMap[nei];
1127                 }
1128             }
1129         }
1130     }
1131
1132     // Reorder faces into upper-triangular and patch ordering
1133     {
1134         // Create cells (packed storage)
1135         labelList cellFaces;
1136         labelList cellFaceOffsets;
1137         makeCells(nActiveFaces_, cellFaces, cellFaceOffsets);
1138
1139         // Do upper triangular order and patch sorting
1140         labelList localFaceMap;
1141         getFaceOrder
1142         (
1143             nActiveFaces_,
1144             cellFaces,
1145             cellFaceOffsets,
1146
1147             localFaceMap,
1148             patchSizes,
1149             patchStarts
1150         );
1151
1152         // Reorder faces.
1153         reorderCompactFaces(localFaceMap.size(), localFaceMap);
1154     }
1155 }
1156
1157
1158 // Find faces to interpolate to create value for new face. Only used if
1159 // face was inflated from edge or point. Internal faces should only be
1160 // created from internal faces, external faces only from external faces
1161 // (and ideally the same patch)
1162 // Is bit problematic if there are no faces to select, i.e. in polyDualMesh
1163 // an internal face can be created from a boundary edge with no internal
1164 // faces connected to it.
1165 Foam::labelList Foam::polyTopoChange::selectFaces
1166 (
1167     const primitiveMesh& mesh,
1168     const labelList& faceLabels,
1169     const bool internalFacesOnly
1170 )
1171 {
1172     label nFaces = 0;
1173
1174     forAll(faceLabels, i)
1175     {
1176         label faceI = faceLabels[i];
1177
1178         if (internalFacesOnly == mesh.isInternalFace(faceI))
1179         {
1180             nFaces++;
1181         }
1182     }
1183
1184     labelList collectedFaces;
1185
1186     if (nFaces == 0)
1187     {
1188         // Did not find any faces of the correct type so just use any old
1189         // face.
1190         collectedFaces = faceLabels;
1191     }
1192     else
1193     {
1194         collectedFaces.setSize(nFaces);
1195
1196         nFaces = 0;
1197
1198         forAll(faceLabels, i)
1199         {
1200             label faceI = faceLabels[i];
1201
1202             if (internalFacesOnly == mesh.isInternalFace(faceI))
1203             {
1204                 collectedFaces[nFaces++] = faceI;
1205             }
1206         }
1207     }
1208
1209     return collectedFaces;
1210 }
1211
1212
1213 // Calculate pointMap per patch (so from patch point label to old patch point
1214 // label)
1215 void Foam::polyTopoChange::calcPatchPointMap
1216 (
1217     const List<Map<label> >& oldPatchMeshPointMaps,
1218     const polyBoundaryMesh& boundary,
1219     labelListList& patchPointMap
1220 ) const
1221 {
1222     patchPointMap.setSize(boundary.size());
1223
1224     forAll(boundary, patchI)
1225     {
1226         const labelList& meshPoints = boundary[patchI].meshPoints();
1227
1228         const Map<label>& oldMeshPointMap = oldPatchMeshPointMaps[patchI];
1229
1230         labelList& curPatchPointRnb = patchPointMap[patchI];
1231
1232         curPatchPointRnb.setSize(meshPoints.size());
1233
1234         forAll(meshPoints, i)
1235         {
1236             if (meshPoints[i] < pointMap_.size())
1237             {
1238                 // Check if old point was part of same patch
1239                 Map<label>::const_iterator ozmpmIter = oldMeshPointMap.find
1240                 (
1241                     pointMap_[meshPoints[i]]
1242                 );
1243
1244                 if (ozmpmIter != oldMeshPointMap.end())
1245                 {
1246                     curPatchPointRnb[i] = ozmpmIter();
1247                 }
1248                 else
1249                 {
1250                     curPatchPointRnb[i] = -1;
1251                 }
1252             }
1253             else
1254             {
1255                 curPatchPointRnb[i] = -1;
1256             }
1257         }
1258     }
1259 }
1260
1261
1262 void Foam::polyTopoChange::calcFaceInflationMaps
1263 (
1264     const polyMesh& mesh,
1265     List<objectMap>& facesFromPoints,
1266     List<objectMap>& facesFromEdges,
1267     List<objectMap>& facesFromFaces
1268 ) const
1269 {
1270     // Faces inflated from points
1271     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
1272
1273     facesFromPoints.setSize(faceFromPoint_.size());
1274
1275     if (faceFromPoint_.size())
1276     {
1277         label nFacesFromPoints = 0;
1278
1279         // Collect all still existing faces connected to this point.
1280         forAllConstIter(Map<label>, faceFromPoint_, iter)
1281         {
1282             label newFaceI = iter.key();
1283
1284             if (region_[newFaceI] == -1)
1285             {
1286                 // Get internal faces using point on old mesh
1287                 facesFromPoints[nFacesFromPoints++] = objectMap
1288                 (
1289                     newFaceI,
1290                     selectFaces
1291                     (
1292                         mesh,
1293                         mesh.pointFaces()[iter()],
1294                         true
1295                     )
1296                 );
1297             }
1298             else
1299             {
1300                 // Get patch faces using point on old mesh
1301                 facesFromPoints[nFacesFromPoints++] = objectMap
1302                 (
1303                     newFaceI,
1304                     selectFaces
1305                     (
1306                         mesh,
1307                         mesh.pointFaces()[iter()],
1308                         false
1309                     )
1310                 );
1311             }
1312         }
1313     }
1314
1315
1316     // Faces inflated from edges
1317     // ~~~~~~~~~~~~~~~~~~~~~~~~~
1318
1319     facesFromEdges.setSize(faceFromEdge_.size());
1320
1321     if (faceFromEdge_.size())
1322     {
1323         label nFacesFromEdges = 0;
1324
1325         // Collect all still existing faces connected to this edge.
1326         forAllConstIter(Map<label>, faceFromEdge_, iter)
1327         {
1328             label newFaceI = iter.key();
1329
1330             if (region_[newFaceI] == -1)
1331             {
1332                 // Get internal faces using edge on old mesh
1333                 facesFromEdges[nFacesFromEdges++] = objectMap
1334                 (
1335                     newFaceI,
1336                     selectFaces
1337                     (
1338                         mesh,
1339                         mesh.edgeFaces(iter()),
1340                         true
1341                     )
1342                 );
1343             }
1344             else
1345             {
1346                 // Get patch faces using edge on old mesh
1347                 facesFromEdges[nFacesFromEdges++] = objectMap
1348                 (
1349                     newFaceI,
1350                     selectFaces
1351                     (
1352                         mesh,
1353                         mesh.edgeFaces(iter()),
1354                         false
1355                     )
1356                 );
1357             }
1358         }
1359     }
1360
1361
1362     // Faces from face merging
1363     // ~~~~~~~~~~~~~~~~~~~~~~~
1364
1365     getMergeSets
1366     (
1367         reverseFaceMap_,
1368         faceMap_,
1369         facesFromFaces
1370     );
1371 }
1372
1373
1374 void Foam::polyTopoChange::calcCellInflationMaps
1375 (
1376     const polyMesh& mesh,
1377     List<objectMap>& cellsFromPoints,
1378     List<objectMap>& cellsFromEdges,
1379     List<objectMap>& cellsFromFaces,
1380     List<objectMap>& cellsFromCells
1381 ) const
1382 {
1383     cellsFromPoints.setSize(cellFromPoint_.size());
1384
1385     if (cellFromPoint_.size())
1386     {
1387         label nCellsFromPoints = 0;
1388
1389         // Collect all still existing faces connected to this point.
1390         forAllConstIter(Map<label>, cellFromPoint_, iter)
1391         {
1392             cellsFromPoints[nCellsFromPoints++] = objectMap
1393             (
1394                 iter.key(),
1395                 mesh.pointCells()[iter()]
1396             );
1397         }
1398     }
1399
1400
1401     cellsFromEdges.setSize(cellFromEdge_.size());
1402
1403     if (cellFromEdge_.size())
1404     {
1405         label nCellsFromEdges = 0;
1406
1407         // Collect all still existing faces connected to this point.
1408         forAllConstIter(Map<label>, cellFromEdge_, iter)
1409         {
1410             cellsFromEdges[nCellsFromEdges++] = objectMap
1411             (
1412                 iter.key(),
1413                 mesh.edgeCells()[iter()]
1414             );
1415         }
1416     }
1417
1418
1419     cellsFromFaces.setSize(cellFromFace_.size());
1420
1421     if (cellFromFace_.size())
1422     {
1423         label nCellsFromFaces = 0;
1424
1425         labelList cellsAroundFace(2, -1);
1426
1427         // Collect all still existing faces connected to this point.
1428         forAllConstIter(Map<label>, cellFromFace_, iter)
1429         {
1430             label oldFaceI = iter();
1431
1432             if (mesh.isInternalFace(oldFaceI))
1433             {
1434                 cellsAroundFace[0] = mesh.faceOwner()[oldFaceI];
1435                 cellsAroundFace[1] = mesh.faceNeighbour()[oldFaceI];
1436             }
1437             else
1438             {
1439                 cellsAroundFace[0] = mesh.faceOwner()[oldFaceI];
1440                 cellsAroundFace[1] = -1;
1441             }
1442
1443             cellsFromFaces[nCellsFromFaces++] = objectMap
1444             (
1445                 iter.key(),
1446                 cellsAroundFace
1447             );
1448         }
1449     }
1450
1451
1452     // Cells from cell merging
1453     // ~~~~~~~~~~~~~~~~~~~~~~~
1454
1455     getMergeSets
1456     (
1457         reverseCellMap_,
1458         cellMap_,
1459         cellsFromCells
1460     );
1461 }
1462
1463
1464 void Foam::polyTopoChange::resetZones
1465 (
1466     const polyMesh& mesh,
1467     polyMesh& newMesh,
1468     labelListList& pointZoneMap,
1469     labelListList& faceZoneFaceMap,
1470     labelListList& cellZoneMap
1471 ) const
1472 {
1473     // pointZones
1474     // ~~~~~~~~~~
1475
1476     pointZoneMap.setSize(mesh.pointZones().size());
1477     {
1478         const pointZoneMesh& pointZones = mesh.pointZones();
1479
1480         // Count points per zone
1481
1482         labelList nPoints(pointZones.size(), 0);
1483
1484         forAll(pointZone_, pointI)
1485         {
1486             label zoneI = pointZone_[pointI];
1487
1488             if (zoneI >= pointZones.size())
1489             {
1490                 FatalErrorIn
1491                 (
1492                     "resetZones(const polyMesh&, polyMesh&, labelListList&"
1493                     "labelListList&, labelListList&)"
1494                 )   << "Illegal zoneID " << zoneI << " for point "
1495                     << pointI << " coord " << mesh.points()[pointI]
1496                     << abort(FatalError);
1497             }
1498
1499             if (zoneI >= 0)
1500             {
1501                 nPoints[zoneI]++;
1502             }
1503         }
1504
1505         // Distribute points per zone
1506
1507         labelListList addressing(pointZones.size());
1508         forAll(addressing, zoneI)
1509         {
1510             addressing[zoneI].setSize(nPoints[zoneI]);
1511         }
1512         nPoints = 0;
1513
1514         forAll(pointZone_, pointI)
1515         {
1516             label zoneI = pointZone_[pointI];
1517             if (zoneI >= 0)
1518             {
1519                 addressing[zoneI][nPoints[zoneI]++] = pointI;
1520             }
1521         }
1522         // Sort the addressing
1523         forAll(addressing, zoneI)
1524         {
1525             stableSort(addressing[zoneI]);
1526         }
1527
1528         // So now we both have old zones and the new addressing.
1529         // Invert the addressing to get pointZoneMap.
1530         forAll(addressing, zoneI)
1531         {
1532             const pointZone& oldZone = pointZones[zoneI];
1533             const labelList& newZoneAddr = addressing[zoneI];
1534
1535             labelList& curPzRnb = pointZoneMap[zoneI];
1536             curPzRnb.setSize(newZoneAddr.size());
1537
1538             forAll(newZoneAddr, i)
1539             {
1540                 if (newZoneAddr[i] < pointMap_.size())
1541                 {
1542                     curPzRnb[i] = oldZone.whichPoint(pointMap_[newZoneAddr[i]]);
1543                 }
1544                 else
1545                 {
1546                     curPzRnb[i] = -1;
1547                 }
1548             }
1549         }
1550
1551         // Reset the addresing on the zone
1552         newMesh.pointZones().clearAddressing();
1553         forAll(newMesh.pointZones(), zoneI)
1554         {
1555             if (debug)
1556             {
1557                 Pout<< "pointZone:" << zoneI
1558                     << "  name:" << newMesh.pointZones()[zoneI].name()
1559                     << "  size:" << addressing[zoneI].size()
1560                     << endl;
1561             }
1562
1563             newMesh.pointZones()[zoneI] = addressing[zoneI];
1564         }
1565     }
1566
1567
1568     // faceZones
1569     // ~~~~~~~~~
1570
1571     faceZoneFaceMap.setSize(mesh.faceZones().size());
1572     {
1573         const faceZoneMesh& faceZones = mesh.faceZones();
1574
1575         labelList nFaces(faceZones.size(), 0);
1576
1577         forAll(faceZone_, faceI)
1578         {
1579             label zoneI = faceZone_[faceI];
1580
1581             if (zoneI >= faceZones.size())
1582             {
1583                 FatalErrorIn
1584                 (
1585                     "resetZones(const polyMesh&, polyMesh&, labelListList&"
1586                     "labelListList&, labelListList&)"
1587                 )   << "Illegal zoneID " << zoneI << " for face "
1588                     << faceI
1589                     << abort(FatalError);
1590             }
1591             if (zoneI >= 0)
1592             {
1593                 nFaces[zoneI]++;
1594             }
1595         }
1596
1597         labelListList addressing(faceZones.size());
1598         boolListList flipMode(faceZones.size());
1599
1600         forAll(addressing, zoneI)
1601         {
1602             addressing[zoneI].setSize(nFaces[zoneI]);
1603             flipMode[zoneI].setSize(nFaces[zoneI]);
1604         }
1605         nFaces = 0;
1606
1607         forAll(faceZone_, faceI)
1608         {
1609             label zoneI = faceZone_[faceI];
1610
1611             if (zoneI >= 0)
1612             {
1613                 label index = nFaces[zoneI]++;
1614                 addressing[zoneI][index] = faceI;
1615                 flipMode[zoneI][index] = faceZoneFlip_[faceI];
1616             }
1617         }
1618         // Sort the addressing
1619         forAll(addressing, zoneI)
1620         {
1621             labelList newToOld;
1622             sortedOrder(addressing[zoneI], newToOld);
1623             {
1624                 labelList newAddressing(addressing[zoneI].size());
1625                 forAll(newAddressing, i)
1626                 {
1627                     newAddressing[i] = addressing[zoneI][newToOld[i]];
1628                 }
1629                 addressing[zoneI].transfer(newAddressing);
1630             }
1631             {
1632                 boolList newFlipMode(flipMode[zoneI].size());
1633                 forAll(newFlipMode, i)
1634                 {
1635                     newFlipMode[i] = flipMode[zoneI][newToOld[i]];
1636                 }
1637                 flipMode[zoneI].transfer(newFlipMode);
1638             }
1639         }
1640
1641         // So now we both have old zones and the new addressing.
1642         // Invert the addressing to get faceZoneFaceMap.
1643         forAll(addressing, zoneI)
1644         {
1645             const faceZone& oldZone = faceZones[zoneI];
1646             const labelList& newZoneAddr = addressing[zoneI];
1647
1648             labelList& curFzFaceRnb = faceZoneFaceMap[zoneI];
1649
1650             curFzFaceRnb.setSize(newZoneAddr.size());
1651
1652             forAll(newZoneAddr, i)
1653             {
1654                 if (newZoneAddr[i] < faceMap_.size())
1655                 {
1656                     curFzFaceRnb[i] =
1657                         oldZone.whichFace(faceMap_[newZoneAddr[i]]);
1658                 }
1659                 else
1660                 {
1661                     curFzFaceRnb[i] = -1;
1662                 }
1663             }
1664         }
1665
1666
1667         // Reset the addresing on the zone
1668         newMesh.faceZones().clearAddressing();
1669         forAll(newMesh.faceZones(), zoneI)
1670         {
1671             if (debug)
1672             {
1673                 Pout<< "faceZone:" << zoneI
1674                     << "  name:" << newMesh.faceZones()[zoneI].name()
1675                     << "  size:" << addressing[zoneI].size()
1676                     << endl;
1677             }
1678
1679             newMesh.faceZones()[zoneI].resetAddressing
1680             (
1681                 addressing[zoneI],
1682                 flipMode[zoneI]
1683             );
1684         }
1685     }
1686
1687
1688     // cellZones
1689     // ~~~~~~~~~
1690
1691     cellZoneMap.setSize(mesh.cellZones().size());
1692     {
1693         const cellZoneMesh& cellZones = mesh.cellZones();
1694
1695         labelList nCells(cellZones.size(), 0);
1696
1697         forAll(cellZone_, cellI)
1698         {
1699             label zoneI = cellZone_[cellI];
1700
1701             if (zoneI >= cellZones.size())
1702             {
1703                 FatalErrorIn
1704                 (
1705                     "resetZones(const polyMesh&, polyMesh&, labelListList&"
1706                     "labelListList&, labelListList&)"
1707                 )   << "Illegal zoneID " << zoneI << " for cell "
1708                     << cellI << abort(FatalError);
1709             }
1710
1711             if (zoneI >= 0)
1712             {
1713                 nCells[zoneI]++;
1714             }
1715         }
1716
1717         labelListList addressing(cellZones.size());
1718         forAll(addressing, zoneI)
1719         {
1720             addressing[zoneI].setSize(nCells[zoneI]);
1721         }
1722         nCells = 0;
1723
1724         forAll(cellZone_, cellI)
1725         {
1726             label zoneI = cellZone_[cellI];
1727
1728             if (zoneI >= 0)
1729             {
1730                 addressing[zoneI][nCells[zoneI]++] = cellI;
1731             }
1732         }
1733         // Sort the addressing
1734         forAll(addressing, zoneI)
1735         {
1736             stableSort(addressing[zoneI]);
1737         }
1738
1739         // So now we both have old zones and the new addressing.
1740         // Invert the addressing to get cellZoneMap.
1741         forAll(addressing, zoneI)
1742         {
1743             const cellZone& oldZone = cellZones[zoneI];
1744             const labelList& newZoneAddr = addressing[zoneI];
1745
1746             labelList& curCellRnb = cellZoneMap[zoneI];
1747
1748             curCellRnb.setSize(newZoneAddr.size());
1749
1750             forAll(newZoneAddr, i)
1751             {
1752                 if (newZoneAddr[i] < cellMap_.size())
1753                 {
1754                     curCellRnb[i] =
1755                         oldZone.whichCell(cellMap_[newZoneAddr[i]]);
1756                 }
1757                 else
1758                 {
1759                     curCellRnb[i] = -1;
1760                 }
1761             }
1762         }
1763
1764         // Reset the addresing on the zone
1765         newMesh.cellZones().clearAddressing();
1766         forAll(newMesh.cellZones(), zoneI)
1767         {
1768             if (debug)
1769             {
1770                 Pout<< "cellZone:" << zoneI
1771                     << "  name:" << newMesh.cellZones()[zoneI].name()
1772                     << "  size:" << addressing[zoneI].size()
1773                     << endl;
1774             }
1775
1776             newMesh.cellZones()[zoneI] = addressing[zoneI];
1777         }
1778     }
1779 }
1780
1781
1782 void Foam::polyTopoChange::calcFaceZonePointMap
1783 (
1784     const polyMesh& mesh,
1785     const List<Map<label> >& oldFaceZoneMeshPointMaps,
1786     labelListList& faceZonePointMap
1787 ) const
1788 {
1789     const faceZoneMesh& faceZones = mesh.faceZones();
1790
1791     faceZonePointMap.setSize(faceZones.size());
1792
1793     forAll(faceZones, zoneI)
1794     {
1795         const faceZone& newZone = faceZones[zoneI];
1796
1797         const labelList& newZoneMeshPoints = newZone().meshPoints();
1798
1799         const Map<label>& oldZoneMeshPointMap = oldFaceZoneMeshPointMaps[zoneI];
1800
1801         labelList& curFzPointRnb = faceZonePointMap[zoneI];
1802
1803         curFzPointRnb.setSize(newZoneMeshPoints.size());
1804
1805         forAll(newZoneMeshPoints, pointI)
1806         {
1807             if (newZoneMeshPoints[pointI] < pointMap_.size())
1808             {
1809                 Map<label>::const_iterator ozmpmIter =
1810                     oldZoneMeshPointMap.find
1811                     (
1812                         pointMap_[newZoneMeshPoints[pointI]]
1813                     );
1814
1815                 if (ozmpmIter != oldZoneMeshPointMap.end())
1816                 {
1817                     curFzPointRnb[pointI] = ozmpmIter();
1818                 }
1819                 else
1820                 {
1821                     curFzPointRnb[pointI] = -1;
1822                 }
1823             }
1824             else
1825             {
1826                 curFzPointRnb[pointI] = -1;
1827             }
1828         }
1829     }
1830 }
1831
1832
1833 Foam::face Foam::polyTopoChange::rotateFace
1834 (
1835     const face& f,
1836     const label nPos
1837 )
1838 {
1839     face newF(f.size());
1840
1841     forAll(f, fp)
1842     {
1843         label fp1 = (fp + nPos) % f.size();
1844
1845         if (fp1 < 0)
1846         {
1847             fp1 += f.size();
1848         }
1849
1850         newF[fp1] = f[fp];
1851     }
1852
1853     return newF;
1854 }
1855
1856
1857 void Foam::polyTopoChange::reorderCoupledFaces
1858 (
1859     const bool syncParallel,
1860     const polyBoundaryMesh& boundary,
1861     const labelList& patchStarts,
1862     const labelList& patchSizes,
1863     const pointField& points
1864 )
1865 {
1866     // Mapping for faces (old to new). Extends over all mesh faces for
1867     // convenience (could be just the external faces)
1868     labelList oldToNew(identity(faces_.size()));
1869
1870     // Rotation on new faces.
1871     labelList rotation(faces_.size(), 0);
1872
1873     // Send ordering
1874     forAll(boundary, patchI)
1875     {
1876         if (syncParallel || !isA<processorPolyPatch>(boundary[patchI]))
1877         {
1878             boundary[patchI].initOrder
1879             (
1880                 primitivePatch
1881                 (
1882                     SubList<face>
1883                     (
1884                         faces_,
1885                         patchSizes[patchI],
1886                         patchStarts[patchI]
1887                     ),
1888                     points
1889                 )
1890             );
1891         }
1892     }
1893
1894     // Receive and calculate ordering
1895
1896     bool anyChanged = false;
1897
1898     forAll(boundary, patchI)
1899     {
1900         if (syncParallel || !isA<processorPolyPatch>(boundary[patchI]))
1901         {
1902             labelList patchFaceMap(patchSizes[patchI], -1);
1903             labelList patchFaceRotation(patchSizes[patchI], 0);
1904
1905             bool changed = boundary[patchI].order
1906             (
1907                 primitivePatch
1908                 (
1909                     SubList<face>
1910                     (
1911                         faces_,
1912                         patchSizes[patchI],
1913                         patchStarts[patchI]
1914                     ),
1915                     points
1916                 ),
1917                 patchFaceMap,
1918                 patchFaceRotation
1919             );
1920
1921             if (changed)
1922             {
1923                 // Merge patch face reordering into mesh face reordering table
1924                 label start = patchStarts[patchI];
1925
1926                 forAll(patchFaceMap, patchFaceI)
1927                 {
1928                     oldToNew[patchFaceI + start] =
1929                         start + patchFaceMap[patchFaceI];
1930                 }
1931
1932                 forAll(patchFaceRotation, patchFaceI)
1933                 {
1934                     rotation[patchFaceI + start] =
1935                         patchFaceRotation[patchFaceI];
1936                 }
1937
1938                 anyChanged = true;
1939             }
1940         }
1941     }
1942
1943     if (syncParallel)
1944     {
1945         reduce(anyChanged, orOp<bool>());
1946     }
1947
1948     if (anyChanged)
1949     {
1950         // Reorder faces according to oldToNew.
1951         reorderCompactFaces(oldToNew.size(), oldToNew);
1952
1953         // Rotate faces (rotation is already in new face indices).
1954         forAll(rotation, faceI)
1955         {
1956             if (rotation[faceI] != 0)
1957             {
1958                 faces_[faceI] = rotateFace(faces_[faceI], rotation[faceI]);
1959             }
1960         }
1961     }
1962 }
1963
1964
1965 void Foam::polyTopoChange::compactAndReorder
1966 (
1967     const polyMesh& mesh,
1968     const bool syncParallel,
1969     const bool orderCells,
1970     const bool orderPoints,
1971
1972     label& nInternalPoints,
1973     pointField& newPoints,
1974     labelList& patchSizes,
1975     labelList& patchStarts,
1976     List<objectMap>& pointsFromPoints,
1977     List<objectMap>& facesFromPoints,
1978     List<objectMap>& facesFromEdges,
1979     List<objectMap>& facesFromFaces,
1980     List<objectMap>& cellsFromPoints,
1981     List<objectMap>& cellsFromEdges,
1982     List<objectMap>& cellsFromFaces,
1983     List<objectMap>& cellsFromCells,
1984     List<Map<label> >& oldPatchMeshPointMaps,
1985     labelList& oldPatchNMeshPoints,
1986     labelList& oldPatchStarts,
1987     List<Map<label> >& oldFaceZoneMeshPointMaps
1988 )
1989 {
1990     if (mesh.boundaryMesh().size() != nPatches_)
1991     {
1992         FatalErrorIn("polyTopoChange::compactAndReorder(..)")
1993             << "polyTopoChange was constructed with a mesh with "
1994             << nPatches_ << " patches." << endl
1995             << "The mesh now provided has a different number of patches "
1996             << mesh.boundaryMesh().size()
1997             << " which is illegal" << endl
1998             << abort(FatalError);
1999     }
2000
2001     // Remove any holes from points/faces/cells and sort faces.
2002     // Sets nActiveFaces_.
2003     compact(orderCells, orderPoints, nInternalPoints, patchSizes, patchStarts);
2004
2005     // Transfer points to pointField. points_ are now cleared!
2006     // Only done since e.g. reorderCoupledFaces requires pointField.
2007     newPoints.transfer(points_);
2008
2009     // Reorder any coupled faces
2010     reorderCoupledFaces
2011     (
2012         syncParallel,
2013         mesh.boundaryMesh(),
2014         patchStarts,
2015         patchSizes,
2016         newPoints
2017     );
2018
2019
2020     // Calculate inflation/merging maps
2021     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2022     // These are for the new face(/point/cell) the old faces whose value
2023     // needs to be
2024     // averaged/summed to get the new value. These old faces come either from
2025     // merged old faces (face remove into other face),
2026     // the old edgeFaces (inflate from edge) or the old pointFaces (inflate
2027     // from point). As an additional complexity will use only internal faces
2028     // to create new value for internal face and vice versa only patch
2029     // faces to to create patch face value.
2030
2031     // For point only point merging
2032     getMergeSets
2033     (
2034         reversePointMap_,
2035         pointMap_,
2036         pointsFromPoints
2037     );
2038
2039     calcFaceInflationMaps
2040     (
2041         mesh,
2042         facesFromPoints,
2043         facesFromEdges,
2044         facesFromFaces
2045     );
2046
2047     calcCellInflationMaps
2048     (
2049         mesh,
2050         cellsFromPoints,
2051         cellsFromEdges,
2052         cellsFromFaces,
2053         cellsFromCells
2054     );
2055
2056     // Clear inflation info
2057     {
2058         faceFromPoint_.clearStorage();
2059         faceFromEdge_.clearStorage();
2060
2061         cellFromPoint_.clearStorage();
2062         cellFromEdge_.clearStorage();
2063         cellFromFace_.clearStorage();
2064     }
2065
2066
2067     const polyBoundaryMesh& boundary = mesh.boundaryMesh();
2068
2069     // Grab patch mesh point maps
2070     oldPatchMeshPointMaps.setSize(boundary.size());
2071     oldPatchNMeshPoints.setSize(boundary.size());
2072     oldPatchStarts.setSize(boundary.size());
2073
2074     forAll(boundary, patchI)
2075     {
2076         // Copy old face zone mesh point maps
2077         oldPatchMeshPointMaps[patchI] = boundary[patchI].meshPointMap();
2078         oldPatchNMeshPoints[patchI] = boundary[patchI].meshPoints().size();
2079         oldPatchStarts[patchI] = boundary[patchI].start();
2080     }
2081
2082     // Grab old face zone mesh point maps.
2083     // These need to be saved before resetting the mesh and are used
2084     // later on to calculate the faceZone pointMaps.
2085     oldFaceZoneMeshPointMaps.setSize(mesh.faceZones().size());
2086
2087     forAll(mesh.faceZones(), zoneI)
2088     {
2089         const faceZone& oldZone = mesh.faceZones()[zoneI];
2090
2091         oldFaceZoneMeshPointMaps[zoneI] = oldZone().meshPointMap();
2092     }
2093 }
2094
2095
2096 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
2097
2098 // Construct from components
2099 Foam::polyTopoChange::polyTopoChange(const label nPatches, const bool strict)
2100 :
2101     strict_(strict),
2102     nPatches_(nPatches),
2103     points_(0),
2104     pointMap_(0),
2105     reversePointMap_(0),
2106     pointZone_(0),
2107     retiredPoints_(0),
2108     faces_(0),
2109     region_(0),
2110     faceOwner_(0),
2111     faceNeighbour_(0),
2112     faceMap_(0),
2113     reverseFaceMap_(0),
2114     faceFromPoint_(0),
2115     faceFromEdge_(0),
2116     flipFaceFlux_(0),
2117     faceZone_(0),
2118     faceZoneFlip_(0),
2119     nActiveFaces_(0),
2120     cellMap_(0),
2121     reverseCellMap_(0),
2122     cellFromPoint_(0),
2123     cellFromEdge_(0),
2124     cellFromFace_(0),
2125     cellZone_(0)
2126 {}
2127
2128
2129 // Construct from components
2130 Foam::polyTopoChange::polyTopoChange
2131 (
2132     const polyMesh& mesh,
2133     const bool strict
2134 )
2135 :
2136     strict_(strict),
2137     nPatches_(0),
2138     points_(0),
2139     pointMap_(0),
2140     reversePointMap_(0),
2141     pointZone_(0),
2142     retiredPoints_(0),
2143     faces_(0),
2144     region_(0),
2145     faceOwner_(0),
2146     faceNeighbour_(0),
2147     faceMap_(0),
2148     reverseFaceMap_(0),
2149     faceFromPoint_(0),
2150     faceFromEdge_(0),
2151     flipFaceFlux_(0),
2152     faceZone_(0),
2153     faceZoneFlip_(0),
2154     nActiveFaces_(0),
2155     cellMap_(0),
2156     reverseCellMap_(0),
2157     cellFromPoint_(0),
2158     cellFromEdge_(0),
2159     cellFromFace_(0),
2160     cellZone_(0)
2161 {
2162     addMesh
2163     (
2164         mesh,
2165         identity(mesh.boundaryMesh().size()),
2166         identity(mesh.pointZones().size()),
2167         identity(mesh.faceZones().size()),
2168         identity(mesh.cellZones().size())
2169     );
2170 }
2171
2172
2173 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
2174
2175 void Foam::polyTopoChange::clear()
2176 {
2177     points_.clearStorage();
2178     pointMap_.clearStorage();
2179     reversePointMap_.clearStorage();
2180     pointZone_.clearStorage();
2181     retiredPoints_.clearStorage();
2182
2183     faces_.clearStorage();
2184     region_.clearStorage();
2185     faceOwner_.clearStorage();
2186     faceNeighbour_.clearStorage();
2187     faceMap_.clearStorage();
2188     reverseFaceMap_.clearStorage();
2189     faceFromPoint_.clearStorage();
2190     faceFromEdge_.clearStorage();
2191     flipFaceFlux_.clearStorage();
2192     faceZone_.clearStorage();
2193     faceZoneFlip_.clearStorage();
2194     nActiveFaces_ = 0;
2195
2196     cellMap_.clearStorage();
2197     reverseCellMap_.clearStorage();
2198     cellZone_.clearStorage();
2199     cellFromPoint_.clearStorage();
2200     cellFromEdge_.clearStorage();
2201     cellFromFace_.clearStorage();
2202 }
2203
2204
2205 void Foam::polyTopoChange::addMesh
2206 (
2207     const polyMesh& mesh,
2208     const labelList& patchMap,
2209     const labelList& pointZoneMap,
2210     const labelList& faceZoneMap,
2211     const labelList& cellZoneMap
2212 )
2213 {
2214     label maxRegion = nPatches_ - 1;
2215     forAll(patchMap, i)
2216     {
2217         maxRegion = max(maxRegion, patchMap[i]);
2218     }
2219     nPatches_ = maxRegion + 1;
2220
2221
2222     // Add points
2223     {
2224         const pointField& points = mesh.points();
2225         const pointZoneMesh& pointZones = mesh.pointZones();
2226
2227         // Extend
2228         points_.setCapacity(points_.size() + points.size());
2229         pointMap_.setCapacity(pointMap_.size() + points.size());
2230         reversePointMap_.setCapacity(reversePointMap_.size() + points.size());
2231         pointZone_.setCapacity(pointZone_.size() + points.size());
2232
2233         // Precalc offset zones
2234         labelList newZoneID(points.size(), -1);
2235
2236         forAll(pointZones, zoneI)
2237         {
2238             const labelList& pointLabels = pointZones[zoneI];
2239
2240             forAll(pointLabels, j)
2241             {
2242                 newZoneID[pointLabels[j]] = pointZoneMap[zoneI];
2243             }
2244         }
2245
2246         // Add points in mesh order
2247         for (label pointI = 0; pointI < mesh.nPoints(); pointI++)
2248         {
2249             addPoint
2250             (
2251                 points[pointI],
2252                 pointI,
2253                 newZoneID[pointI],
2254                 true
2255             );
2256         }
2257     }
2258
2259     // Add cells
2260     {
2261         const cellZoneMesh& cellZones = mesh.cellZones();
2262
2263         // Resize
2264
2265         // Note: polyMesh does not allow retired cells anymore. So allCells
2266         // always equals nCells
2267         label nAllCells = mesh.nCells();
2268
2269         cellMap_.setCapacity(cellMap_.size() + nAllCells);
2270         reverseCellMap_.setCapacity(reverseCellMap_.size() + nAllCells);
2271         cellFromPoint_.resize(cellFromPoint_.size() + nAllCells/100);
2272         cellFromEdge_.resize(cellFromEdge_.size() + nAllCells/100);
2273         cellFromFace_.resize(cellFromFace_.size() + nAllCells/100);
2274         cellZone_.setCapacity(cellZone_.size() + nAllCells);
2275
2276
2277         // Precalc offset zones
2278         labelList newZoneID(nAllCells, -1);
2279
2280         forAll(cellZones, zoneI)
2281         {
2282             const labelList& cellLabels = cellZones[zoneI];
2283
2284             forAll(cellLabels, j)
2285             {
2286                 label cellI = cellLabels[j];
2287
2288                 if (newZoneID[cellI] != -1)
2289                 {
2290                     WarningIn
2291                     (
2292                         "polyTopoChange::addMesh"
2293                         "(const polyMesh&, const labelList&,"
2294                         "const labelList&, const labelList&,"
2295                         "const labelList&)"
2296                     )   << "Cell:" << cellI
2297                         << " centre:" << mesh.cellCentres()[cellI]
2298                         << " is in two zones:"
2299                         << cellZones[newZoneID[cellI]].name()
2300                         << " and " << cellZones[zoneI].name() << endl
2301                         << "    This is not supported."
2302                         << " Continuing with first zone only." << endl;
2303                 }
2304                 else
2305                 {
2306                     newZoneID[cellI] = cellZoneMap[zoneI];
2307                 }
2308             }
2309         }
2310
2311         // Add cells in mesh order
2312         for (label cellI = 0; cellI < nAllCells; cellI++)
2313         {
2314             // Add cell from cell
2315             addCell(-1, -1, -1, cellI, newZoneID[cellI]);
2316         }
2317     }
2318
2319     // Add faces
2320     {
2321         const polyBoundaryMesh& patches = mesh.boundaryMesh();
2322         const faceList& faces = mesh.faces();
2323         const labelList& faceOwner = mesh.faceOwner();
2324         const labelList& faceNeighbour = mesh.faceNeighbour();
2325         const faceZoneMesh& faceZones = mesh.faceZones();
2326
2327         // Resize
2328         label nAllFaces = mesh.faces().size();
2329
2330         faces_.setCapacity(faces_.size() + nAllFaces);
2331         region_.setCapacity(region_.size() + nAllFaces);
2332         faceOwner_.setCapacity(faceOwner_.size() + nAllFaces);
2333         faceNeighbour_.setCapacity(faceNeighbour_.size() + nAllFaces);
2334         faceMap_.setCapacity(faceMap_.size() + nAllFaces);
2335         reverseFaceMap_.setCapacity(reverseFaceMap_.size() + nAllFaces);
2336         faceFromPoint_.resize(faceFromPoint_.size() + nAllFaces/100);
2337         faceFromEdge_.resize(faceFromEdge_.size() + nAllFaces/100);
2338         flipFaceFlux_.setCapacity(flipFaceFlux_.size() + nAllFaces);
2339         faceZone_.setCapacity(faceZone_.size() + nAllFaces);
2340         faceZoneFlip_.setCapacity(faceZoneFlip_.size() + nAllFaces);
2341
2342
2343         // Precalc offset zones
2344         labelList newZoneID(nAllFaces, -1);
2345         boolList zoneFlip(nAllFaces, false);
2346
2347         forAll(faceZones, zoneI)
2348         {
2349             const labelList& faceLabels = faceZones[zoneI];
2350             const boolList& flipMap = faceZones[zoneI].flipMap();
2351
2352             forAll(faceLabels, j)
2353             {
2354                 newZoneID[faceLabels[j]] = faceZoneMap[zoneI];
2355                 zoneFlip[faceLabels[j]] = flipMap[j];
2356             }
2357         }
2358
2359         // Add faces in mesh order
2360
2361         // 1. Internal faces
2362         for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
2363         {
2364             addFace
2365             (
2366                 faces[faceI],
2367                 faceOwner[faceI],
2368                 faceNeighbour[faceI],
2369                 -1,                         // masterPointID
2370                 -1,                         // masterEdgeID
2371                 faceI,                      // masterFaceID
2372                 false,                      // flipFaceFlux
2373                 -1,                         // patchID
2374                 newZoneID[faceI],           // zoneID
2375                 zoneFlip[faceI]             // zoneFlip
2376             );
2377         }
2378
2379         // 2. Patch faces
2380         forAll(patches, patchI)
2381         {
2382             const polyPatch& pp = patches[patchI];
2383
2384             if (pp.start() != faces_.size())
2385             {
2386                 FatalErrorIn
2387                 (
2388                     "polyTopoChange::polyTopoChange"
2389                     "(const polyMesh& mesh, const bool strict)"
2390                 )   << "Problem : "
2391                     << "Patch " << pp.name() << " starts at " << pp.start()
2392                     << endl
2393                     << "Current face counter at " << faces_.size() << endl
2394                     << "Are patches in incremental order?"
2395                     << abort(FatalError);
2396             }
2397             forAll(pp, patchFaceI)
2398             {
2399                 label faceI = pp.start() + patchFaceI;
2400
2401                 addFace
2402                 (
2403                     faces[faceI],
2404                     faceOwner[faceI],
2405                     -1,                         // neighbour
2406                     -1,                         // masterPointID
2407                     -1,                         // masterEdgeID
2408                     faceI,                      // masterFaceID
2409                     false,                      // flipFaceFlux
2410                     patchMap[patchI],           // patchID
2411                     newZoneID[faceI],           // zoneID
2412                     zoneFlip[faceI]             // zoneFlip
2413                 );
2414             }
2415         }
2416     }
2417 }
2418
2419
2420 void Foam::polyTopoChange::setCapacity
2421 (
2422     const label nPoints,
2423     const label nFaces,
2424     const label nCells
2425 )
2426 {
2427     points_.setCapacity(nPoints);
2428     pointMap_.setCapacity(nPoints);
2429     reversePointMap_.setCapacity(nPoints);
2430     pointZone_.setCapacity(nPoints);
2431
2432     faces_.setCapacity(nFaces);
2433     region_.setCapacity(nFaces);
2434     faceOwner_.setCapacity(nFaces);
2435     faceNeighbour_.setCapacity(nFaces);
2436     faceMap_.setCapacity(nFaces);
2437     reverseFaceMap_.setCapacity(nFaces);
2438     faceFromPoint_.resize(faceFromPoint_.size() + nFaces/100);
2439     faceFromEdge_.resize(faceFromEdge_.size() + nFaces/100);
2440     flipFaceFlux_.setCapacity(nFaces);
2441     faceZone_.setCapacity(nFaces);
2442     faceZoneFlip_.setCapacity(nFaces);
2443
2444     cellMap_.setCapacity(nCells);
2445     reverseCellMap_.setCapacity(nCells);
2446     cellFromPoint_.resize(cellFromPoint_.size() + nCells/100);
2447     cellFromEdge_.resize(cellFromEdge_.size() + nCells/100);
2448     cellFromFace_.resize(cellFromFace_.size() + nCells/100);
2449     cellZone_.setCapacity(nCells);
2450 }
2451
2452
2453 Foam::label Foam::polyTopoChange::setAction(const topoAction& action)
2454 {
2455     if (isType<polyAddPoint>(action))
2456     {
2457         const polyAddPoint& pap = refCast<const polyAddPoint>(action);
2458
2459         return addPoint
2460         (
2461             pap.newPoint(),
2462             pap.masterPointID(),
2463             pap.zoneID(),
2464             pap.inCell()
2465         );
2466     }
2467     else if (isType<polyModifyPoint>(action))
2468     {
2469         const polyModifyPoint& pmp = refCast<const polyModifyPoint>(action);
2470
2471         modifyPoint
2472         (
2473             pmp.pointID(),
2474             pmp.newPoint(),
2475             pmp.zoneID(),
2476             pmp.inCell()
2477         );
2478
2479         return -1;
2480     }
2481     else if (isType<polyRemovePoint>(action))
2482     {
2483         const polyRemovePoint& prp = refCast<const polyRemovePoint>(action);
2484
2485         removePoint(prp.pointID(), prp.mergePointID());
2486
2487         return -1;
2488     }
2489     else if (isType<polyAddFace>(action))
2490     {
2491         const polyAddFace& paf = refCast<const polyAddFace>(action);
2492
2493         return addFace
2494         (
2495             paf.newFace(),
2496             paf.owner(),
2497             paf.neighbour(),
2498             paf.masterPointID(),
2499             paf.masterEdgeID(),
2500             paf.masterFaceID(),
2501             paf.flipFaceFlux(),
2502             paf.patchID(),
2503             paf.zoneID(),
2504             paf.zoneFlip()
2505         );
2506     }
2507     else if (isType<polyModifyFace>(action))
2508     {
2509         const polyModifyFace& pmf = refCast<const polyModifyFace>(action);
2510
2511         modifyFace
2512         (
2513             pmf.newFace(),
2514             pmf.faceID(),
2515             pmf.owner(),
2516             pmf.neighbour(),
2517             pmf.flipFaceFlux(),
2518             pmf.patchID(),
2519             pmf.zoneID(),
2520             pmf.zoneFlip()
2521         );
2522
2523         return -1;
2524     }
2525     else if (isType<polyRemoveFace>(action))
2526     {
2527         const polyRemoveFace& prf = refCast<const polyRemoveFace>(action);
2528
2529         removeFace(prf.faceID(), prf.mergeFaceID());
2530
2531         return -1;
2532     }
2533     else if (isType<polyAddCell>(action))
2534     {
2535         const polyAddCell& pac = refCast<const polyAddCell>(action);
2536
2537         return addCell
2538         (
2539             pac.masterPointID(),
2540             pac.masterEdgeID(),
2541             pac.masterFaceID(),
2542             pac.masterCellID(),
2543             pac.zoneID()
2544         );
2545     }
2546     else if (isType<polyModifyCell>(action))
2547     {
2548         const polyModifyCell& pmc = refCast<const polyModifyCell>(action);
2549
2550         if (pmc.removeFromZone())
2551         {
2552             modifyCell(pmc.cellID(), -1);
2553         }
2554         else
2555         {
2556             modifyCell(pmc.cellID(), pmc.zoneID());
2557         }
2558
2559         return -1;
2560     }
2561     else if (isType<polyRemoveCell>(action))
2562     {
2563         const polyRemoveCell& prc = refCast<const polyRemoveCell>(action);
2564
2565         removeCell(prc.cellID(), prc.mergeCellID());
2566
2567         return -1;
2568     }
2569     else
2570     {
2571         FatalErrorIn
2572         (
2573             "label polyTopoChange::setAction(const topoAction& action)"
2574         )   << "Unknown type of topoChange: " << action.type()
2575             << abort(FatalError);
2576
2577         // Dummy return to keep compiler happy
2578         return -1;
2579     }
2580 }
2581
2582
2583 Foam::label Foam::polyTopoChange::addPoint
2584 (
2585     const point& pt,
2586     const label masterPointID,
2587     const label zoneID,
2588     const bool inCell
2589 )
2590 {
2591     label pointI = points_.size();
2592
2593     points_.append(pt);
2594     pointMap_.append(masterPointID);
2595     reversePointMap_.append(pointI);
2596     pointZone_.append(zoneID);
2597
2598     if (!inCell)
2599     {
2600         retiredPoints_.insert(pointI);
2601     }
2602
2603     return pointI;
2604 }
2605
2606
2607 void Foam::polyTopoChange::modifyPoint
2608 (
2609     const label pointI,
2610     const point& pt,
2611     const label newZoneID,
2612     const bool inCell
2613 )
2614 {
2615     if (pointI < 0 || pointI >= points_.size())
2616     {
2617         FatalErrorIn
2618         (
2619             "polyTopoChange::modifyPoint(const label, const point&)"
2620         )   << "illegal point label " << pointI << endl
2621             << "Valid point labels are 0 .. " << points_.size()-1
2622             << abort(FatalError);
2623     }
2624     if (pointRemoved(pointI) || pointMap_[pointI] == -1)
2625     {
2626         FatalErrorIn
2627         (
2628             "polyTopoChange::modifyPoint(const label, const point&)"
2629         )   << "point " << pointI << " already marked for removal"
2630             << abort(FatalError);
2631     }
2632     points_[pointI] = pt;
2633     pointZone_[pointI] = newZoneID;
2634
2635     if (inCell)
2636     {
2637         retiredPoints_.erase(pointI);
2638     }
2639     else
2640     {
2641         retiredPoints_.insert(pointI);
2642     }
2643 }
2644
2645
2646 void Foam::polyTopoChange::movePoints(const pointField& newPoints)
2647 {
2648     if (newPoints.size() != points_.size())
2649     {
2650         FatalErrorIn("polyTopoChange::movePoints(const pointField&)")
2651             << "illegal pointField size." << endl
2652             << "Size:" << newPoints.size() << endl
2653             << "Points in mesh:" << points_.size()
2654             << abort(FatalError);
2655     }
2656
2657     forAll(points_, pointI)
2658     {
2659         points_[pointI] = newPoints[pointI];
2660     }
2661 }
2662
2663
2664 void Foam::polyTopoChange::removePoint
2665 (
2666     const label pointI,
2667     const label mergePointI
2668 )
2669 {
2670     if (pointI < 0 || pointI >= points_.size())
2671     {
2672         FatalErrorIn("polyTopoChange::removePoint(const label, const label)")
2673             << "illegal point label " << pointI << endl
2674             << "Valid point labels are 0 .. " << points_.size()-1
2675             << abort(FatalError);
2676     }
2677
2678     if
2679     (
2680         strict_
2681      && (pointRemoved(pointI) || pointMap_[pointI] == -1)
2682     )
2683     {
2684         FatalErrorIn("polyTopoChange::removePoint(const label, const label)")
2685             << "point " << pointI << " already marked for removal" << nl
2686             << "Point:" << points_[pointI] << " pointMap:" << pointMap_[pointI]
2687             << abort(FatalError);
2688     }
2689
2690     if (pointI == mergePointI)
2691     {
2692         FatalErrorIn("polyTopoChange::removePoint(const label, const label)")
2693             << "Cannot remove/merge point " << pointI << " onto itself."
2694             << abort(FatalError);
2695     }
2696
2697     points_[pointI] = greatPoint;
2698     pointMap_[pointI] = -1;
2699     if (mergePointI >= 0)
2700     {
2701         reversePointMap_[pointI] = -mergePointI-2;
2702     }
2703     else
2704     {
2705         reversePointMap_[pointI] = -1;
2706     }
2707     pointZone_[pointI] = -1;
2708     retiredPoints_.erase(pointI);
2709 }
2710
2711
2712 Foam::label Foam::polyTopoChange::addFace
2713 (
2714     const face& f,
2715     const label own,
2716     const label nei,
2717     const label masterPointID,
2718     const label masterEdgeID,
2719     const label masterFaceID,
2720     const bool flipFaceFlux,
2721     const label patchID,
2722     const label zoneID,
2723     const bool zoneFlip
2724 )
2725 {
2726     // Check validity
2727     if (debug)
2728     {
2729         checkFace(f, -1, own, nei, patchID, zoneID);
2730     }
2731
2732     label faceI = faces_.size();
2733
2734     faces_.append(f);
2735     region_.append(patchID);
2736     faceOwner_.append(own);
2737     faceNeighbour_.append(nei);
2738
2739     if (masterPointID >= 0)
2740     {
2741         faceMap_.append(-1);
2742         faceFromPoint_.insert(faceI, masterPointID);
2743     }
2744     else if (masterEdgeID >= 0)
2745     {
2746         faceMap_.append(-1);
2747         faceFromEdge_.insert(faceI, masterEdgeID);
2748     }
2749     else if (masterFaceID >= 0)
2750     {
2751         faceMap_.append(masterFaceID);
2752     }
2753     else
2754     {
2755         // Allow inflate-from-nothing?
2756         //FatalErrorIn("polyTopoChange::addFace")
2757         //    << "Need to specify a master point, edge or face"
2758         //    << "face:" << f << " own:" << own << " nei:" << nei
2759         //    << abort(FatalError);
2760         faceMap_.append(-1);
2761     }
2762     reverseFaceMap_.append(faceI);
2763
2764     flipFaceFlux_[faceI] = (flipFaceFlux ? 1 : 0);
2765     faceZone_.append(zoneID);
2766     faceZoneFlip_[faceI] = (zoneFlip ? 1 : 0);
2767
2768     return faceI;
2769 }
2770
2771
2772 void Foam::polyTopoChange::modifyFace
2773 (
2774     const face& f,
2775     const label faceI,
2776     const label own,
2777     const label nei,
2778     const bool flipFaceFlux,
2779     const label patchID,
2780     const label zoneID,
2781     const bool zoneFlip
2782 )
2783 {
2784     // Check validity
2785     if (debug)
2786     {
2787         checkFace(f, faceI, own, nei, patchID, zoneID);
2788     }
2789
2790     faces_[faceI] = f;
2791     faceOwner_[faceI] = own;
2792     faceNeighbour_[faceI] = nei;
2793     region_[faceI] = patchID;
2794
2795     flipFaceFlux_[faceI] = (flipFaceFlux ? 1 : 0);
2796     faceZone_[faceI] = zoneID;
2797     faceZoneFlip_[faceI] = (zoneFlip ? 1 : 0);
2798 }
2799
2800
2801 void Foam::polyTopoChange::removeFace(const label faceI, const label mergeFaceI)
2802 {
2803     if (faceI < 0 || faceI >= faces_.size())
2804     {
2805         FatalErrorIn("polyTopoChange::removeFace(const label, const label)")
2806             << "illegal face label " << faceI << endl
2807             << "Valid face labels are 0 .. " << faces_.size()-1
2808             << abort(FatalError);
2809     }
2810
2811     if
2812     (
2813         strict_
2814      && (faceRemoved(faceI) || faceMap_[faceI] == -1)
2815     )
2816     {
2817         FatalErrorIn("polyTopoChange::removeFace(const label, const label)")
2818             << "face " << faceI
2819             << " already marked for removal"
2820             << abort(FatalError);
2821     }
2822
2823     faces_[faceI].setSize(0);
2824     region_[faceI] = -1;
2825     faceOwner_[faceI] = -1;
2826     faceNeighbour_[faceI] = -1;
2827     faceMap_[faceI] = -1;
2828     if (mergeFaceI >= 0)
2829     {
2830         reverseFaceMap_[faceI] = -mergeFaceI-2;
2831     }
2832     else
2833     {
2834         reverseFaceMap_[faceI] = -1;
2835     }
2836     faceFromEdge_.erase(faceI);
2837     faceFromPoint_.erase(faceI);
2838     flipFaceFlux_[faceI] = 0;
2839     faceZone_[faceI] = -1;
2840     faceZoneFlip_[faceI] = 0;
2841 }
2842
2843
2844 Foam::label Foam::polyTopoChange::addCell
2845 (
2846     const label masterPointID,
2847     const label masterEdgeID,
2848     const label masterFaceID,
2849     const label masterCellID,
2850     const label zoneID
2851 )
2852 {
2853     label cellI = cellMap_.size();
2854
2855     if (masterPointID >= 0)
2856     {
2857         cellMap_.append(-1);
2858         cellFromPoint_.insert(cellI, masterPointID);
2859     }
2860     else if (masterEdgeID >= 0)
2861     {
2862         cellMap_.append(-1);
2863         cellFromEdge_.insert(cellI, masterEdgeID);
2864     }
2865     else if (masterFaceID >= 0)
2866     {
2867         cellMap_.append(-1);
2868         cellFromFace_.insert(cellI, masterFaceID);
2869     }
2870     else
2871     {
2872         cellMap_.append(masterCellID);
2873     }
2874     reverseCellMap_.append(cellI);
2875     cellZone_.append(zoneID);
2876
2877     return cellI;
2878 }
2879
2880
2881 void Foam::polyTopoChange::modifyCell
2882 (
2883     const label cellI,
2884     const label zoneID
2885 )
2886 {
2887     cellZone_[cellI] = zoneID;
2888 }
2889
2890
2891 void Foam::polyTopoChange::removeCell(const label cellI, const label mergeCellI)
2892 {
2893     if (cellI < 0 || cellI >= cellMap_.size())
2894     {
2895         FatalErrorIn("polyTopoChange::removeCell(const label, const label)")
2896             << "illegal cell label " << cellI << endl
2897             << "Valid cell labels are 0 .. " << cellMap_.size()-1
2898             << abort(FatalError);
2899     }
2900
2901     if (strict_ && cellMap_[cellI] == -2)
2902     {
2903         FatalErrorIn("polyTopoChange::removeCell(const label, const label)")
2904             << "cell " << cellI
2905             << " already marked for removal"
2906             << abort(FatalError);
2907     }
2908
2909     cellMap_[cellI] = -2;
2910     if (mergeCellI >= 0)
2911     {
2912         reverseCellMap_[cellI] = -mergeCellI-2;
2913     }
2914     else
2915     {
2916         reverseCellMap_[cellI] = -1;
2917     }
2918     cellFromPoint_.erase(cellI);
2919     cellFromEdge_.erase(cellI);
2920     cellFromFace_.erase(cellI);
2921     cellZone_[cellI] = -1;
2922 }
2923
2924
2925 Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::changeMesh
2926 (
2927     polyMesh& mesh,
2928     const bool inflate,
2929     const bool syncParallel,
2930     const bool orderCells,
2931     const bool orderPoints
2932 )
2933 {
2934     if (debug)
2935     {
2936         Pout<< "polyTopoChange::changeMesh"
2937             << "(polyMesh&, const bool, const bool, const bool, const bool)"
2938             << endl;
2939     }
2940
2941     if (debug)
2942     {
2943         Pout<< "Old mesh:" << nl;
2944         writeMeshStats(mesh, Pout);
2945     }
2946
2947     // new mesh points
2948     pointField newPoints;
2949     // number of internal points
2950     label nInternalPoints;
2951     // patch slicing
2952     labelList patchSizes;
2953     labelList patchStarts;
2954     // inflate maps
2955     List<objectMap> pointsFromPoints;
2956     List<objectMap> facesFromPoints;
2957     List<objectMap> facesFromEdges;
2958     List<objectMap> facesFromFaces;
2959     List<objectMap> cellsFromPoints;
2960     List<objectMap> cellsFromEdges;
2961     List<objectMap> cellsFromFaces;
2962     List<objectMap> cellsFromCells;
2963     // old mesh info
2964     List<Map<label> > oldPatchMeshPointMaps;
2965     labelList oldPatchNMeshPoints;
2966     labelList oldPatchStarts;
2967     List<Map<label> > oldFaceZoneMeshPointMaps;
2968
2969     // Compact, reorder patch faces and calculate mesh/patch maps.
2970     compactAndReorder
2971     (
2972         mesh,
2973         syncParallel,
2974         orderCells,
2975         orderPoints,
2976
2977         nInternalPoints,
2978         newPoints,
2979         patchSizes,
2980         patchStarts,
2981         pointsFromPoints,
2982         facesFromPoints,
2983         facesFromEdges,
2984         facesFromFaces,
2985         cellsFromPoints,
2986         cellsFromEdges,
2987         cellsFromFaces,
2988         cellsFromCells,
2989         oldPatchMeshPointMaps,
2990         oldPatchNMeshPoints,
2991         oldPatchStarts,
2992         oldFaceZoneMeshPointMaps
2993     );
2994
2995     const label nOldPoints(mesh.nPoints());
2996     const label nOldFaces(mesh.nFaces());
2997     const label nOldCells(mesh.nCells());
2998
2999
3000     // Change the mesh
3001     // ~~~~~~~~~~~~~~~
3002     // This will invalidate any addressing so better make sure you have
3003     // all the information you need!!!
3004
3005     if (inflate)
3006     {
3007         // Keep (renumbered) mesh points, store new points in map for inflation
3008         // (appended points (i.e. from nowhere) get value zero)
3009         pointField renumberedMeshPoints(newPoints.size());
3010
3011         forAll(pointMap_, newPointI)
3012         {
3013             label oldPointI = pointMap_[newPointI];
3014
3015             if (oldPointI >= 0)
3016             {
3017                 renumberedMeshPoints[newPointI] = mesh.points()[oldPointI];
3018             }
3019             else
3020             {
3021                 renumberedMeshPoints[newPointI] = vector::zero;
3022             }
3023         }
3024
3025         mesh.resetPrimitives
3026         (
3027             xferMove(renumberedMeshPoints),
3028             faces_.xfer(),
3029             faceOwner_.xfer(),
3030             faceNeighbour_.xfer(),
3031             patchSizes,
3032             patchStarts,
3033             syncParallel
3034         );
3035
3036         mesh.changing(true);
3037     }
3038     else
3039     {
3040         // Set new points.
3041         mesh.resetPrimitives
3042         (
3043             xferMove(newPoints),
3044             faces_.xfer(),
3045             faceOwner_.xfer(),
3046             faceNeighbour_.xfer(),
3047             patchSizes,
3048             patchStarts,
3049             syncParallel
3050         );
3051         mesh.changing(true);
3052     }
3053
3054     // Clear out primitives
3055     {
3056         retiredPoints_.clearStorage();
3057         region_.clearStorage();
3058     }
3059
3060
3061     if (debug)
3062     {
3063         // Some stats on changes
3064         label nAdd, nInflate, nMerge, nRemove;
3065         countMap(pointMap_, reversePointMap_, nAdd, nInflate, nMerge, nRemove);
3066         Pout<< "Points:"
3067             << "  added(from point):" << nAdd
3068             << "  added(from nothing):" << nInflate
3069             << "  merged(into other point):" << nMerge
3070             << "  removed:" << nRemove
3071             << nl;
3072
3073         countMap(faceMap_, reverseFaceMap_, nAdd, nInflate, nMerge, nRemove);
3074         Pout<< "Faces:"
3075             << "  added(from face):" << nAdd
3076             << "  added(inflated):" << nInflate
3077             << "  merged(into other face):" << nMerge
3078             << "  removed:" << nRemove
3079             << nl;
3080
3081         countMap(cellMap_, reverseCellMap_, nAdd, nInflate, nMerge, nRemove);
3082         Pout<< "Cells:"
3083             << "  added(from cell):" << nAdd
3084             << "  added(inflated):" << nInflate
3085             << "  merged(into other cell):" << nMerge
3086             << "  removed:" << nRemove
3087             << nl
3088             << endl;
3089     }
3090
3091     if (debug)
3092     {
3093         Pout<< "New mesh:" << nl;
3094         writeMeshStats(mesh, Pout);
3095     }
3096
3097
3098     // Zones
3099     // ~~~~~
3100
3101     // Inverse of point/face/cell zone addressing.
3102     // For every preserved point/face/cells in zone give the old position.
3103     // For added points, the index is set to -1
3104     labelListList pointZoneMap(mesh.pointZones().size());
3105     labelListList faceZoneFaceMap(mesh.faceZones().size());
3106     labelListList cellZoneMap(mesh.cellZones().size());
3107
3108     resetZones(mesh, mesh, pointZoneMap, faceZoneFaceMap, cellZoneMap);
3109
3110     // Clear zone info
3111     {
3112         pointZone_.clearStorage();
3113         faceZone_.clearStorage();
3114         faceZoneFlip_.clearStorage();
3115         cellZone_.clearStorage();
3116     }
3117
3118
3119     // Patch point renumbering
3120     // For every preserved point on a patch give the old position.
3121     // For added points, the index is set to -1
3122     labelListList patchPointMap(mesh.boundaryMesh().size());
3123     calcPatchPointMap
3124     (
3125         oldPatchMeshPointMaps,
3126         mesh.boundaryMesh(),
3127         patchPointMap
3128     );
3129
3130     // Create the face zone mesh point renumbering
3131     labelListList faceZonePointMap(mesh.faceZones().size());
3132     calcFaceZonePointMap(mesh, oldFaceZoneMeshPointMaps, faceZonePointMap);
3133
3134     labelHashSet flipFaceFluxSet(getSetIndices(flipFaceFlux_));
3135
3136     return autoPtr<mapPolyMesh>
3137     (
3138         new mapPolyMesh
3139         (
3140             mesh,
3141             nOldPoints,
3142             nOldFaces,
3143             nOldCells,
3144
3145             pointMap_,
3146             pointsFromPoints,
3147
3148             faceMap_,
3149             facesFromPoints,
3150             facesFromEdges,
3151             facesFromFaces,
3152
3153             cellMap_,
3154             cellsFromPoints,
3155             cellsFromEdges,
3156             cellsFromFaces,
3157             cellsFromCells,
3158
3159             reversePointMap_,
3160             reverseFaceMap_,
3161             reverseCellMap_,
3162
3163             flipFaceFluxSet,
3164
3165             patchPointMap,
3166
3167             pointZoneMap,
3168
3169             faceZonePointMap,
3170             faceZoneFaceMap,
3171             cellZoneMap,
3172
3173             newPoints,          // if empty signals no inflation.
3174             oldPatchStarts,
3175             oldPatchNMeshPoints,
3176             true                // steal storage.
3177         )
3178     );
3179
3180     // At this point all member DynamicList (pointMap_, cellMap_ etc.) will
3181     // be invalid.
3182 }
3183
3184
3185 Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::makeMesh
3186 (
3187     autoPtr<fvMesh>& newMeshPtr,
3188     const IOobject& io,
3189     const fvMesh& mesh,
3190     const bool syncParallel,
3191     const bool orderCells,
3192     const bool orderPoints
3193 )
3194 {
3195     if (debug)
3196     {
3197         Pout<< "polyTopoChange::changeMesh"
3198             << "(autoPtr<fvMesh>&, const IOobject&, const fvMesh&"
3199             << ", const bool, const bool, const bool)"
3200             << endl;
3201     }
3202
3203     if (debug)
3204     {
3205         Pout<< "Old mesh:" << nl;
3206         writeMeshStats(mesh, Pout);
3207     }
3208
3209     // new mesh points
3210     pointField newPoints;
3211     // number of internal points
3212     label nInternalPoints;
3213     // patch slicing
3214     labelList patchSizes;
3215     labelList patchStarts;
3216     // inflate maps
3217     List<objectMap> pointsFromPoints;
3218     List<objectMap> facesFromPoints;
3219     List<objectMap> facesFromEdges;
3220     List<objectMap> facesFromFaces;
3221     List<objectMap> cellsFromPoints;
3222     List<objectMap> cellsFromEdges;
3223     List<objectMap> cellsFromFaces;
3224     List<objectMap> cellsFromCells;
3225
3226     // old mesh info
3227     List<Map<label> > oldPatchMeshPointMaps;
3228     labelList oldPatchNMeshPoints;
3229     labelList oldPatchStarts;
3230     List<Map<label> > oldFaceZoneMeshPointMaps;
3231
3232     // Compact, reorder patch faces and calculate mesh/patch maps.
3233     compactAndReorder
3234     (
3235         mesh,
3236         syncParallel,
3237         orderCells,
3238         orderPoints,
3239
3240         nInternalPoints,
3241         newPoints,
3242         patchSizes,
3243         patchStarts,
3244         pointsFromPoints,
3245         facesFromPoints,
3246         facesFromEdges,
3247         facesFromFaces,
3248         cellsFromPoints,
3249         cellsFromEdges,
3250         cellsFromFaces,
3251         cellsFromCells,
3252         oldPatchMeshPointMaps,
3253         oldPatchNMeshPoints,
3254         oldPatchStarts,
3255         oldFaceZoneMeshPointMaps
3256     );
3257
3258     const label nOldPoints(mesh.nPoints());
3259     const label nOldFaces(mesh.nFaces());
3260     const label nOldCells(mesh.nCells());
3261
3262
3263     // Create the mesh
3264     // ~~~~~~~~~~~~~~~
3265
3266     newMeshPtr.reset
3267     (
3268         new fvMesh
3269         (
3270             io,
3271             xferMove(newPoints),
3272             faces_.xfer(),
3273             faceOwner_.xfer(),
3274             faceNeighbour_.xfer()
3275         )
3276     );
3277     fvMesh& newMesh = newMeshPtr();
3278
3279     // Clear out primitives
3280     {
3281         retiredPoints_.clearStorage();
3282         region_.clearStorage();
3283     }
3284
3285
3286     if (debug)
3287     {
3288         // Some stats on changes
3289         label nAdd, nInflate, nMerge, nRemove;
3290         countMap(pointMap_, reversePointMap_, nAdd, nInflate, nMerge, nRemove);
3291         Pout<< "Points:"
3292             << "  added(from point):" << nAdd
3293             << "  added(from nothing):" << nInflate
3294             << "  merged(into other point):" << nMerge
3295             << "  removed:" << nRemove
3296             << nl;
3297
3298         countMap(faceMap_, reverseFaceMap_, nAdd, nInflate, nMerge, nRemove);
3299         Pout<< "Faces:"
3300             << "  added(from face):" << nAdd
3301             << "  added(inflated):" << nInflate
3302             << "  merged(into other face):" << nMerge
3303             << "  removed:" << nRemove
3304             << nl;
3305
3306         countMap(cellMap_, reverseCellMap_, nAdd, nInflate, nMerge, nRemove);
3307         Pout<< "Cells:"
3308             << "  added(from cell):" << nAdd
3309             << "  added(inflated):" << nInflate
3310             << "  merged(into other cell):" << nMerge
3311             << "  removed:" << nRemove
3312             << nl
3313             << endl;
3314     }
3315
3316
3317     {
3318         const polyBoundaryMesh& oldPatches = mesh.boundaryMesh();
3319
3320         List<polyPatch*> newBoundary(oldPatches.size());
3321
3322         forAll(oldPatches, patchI)
3323         {
3324             newBoundary[patchI] = oldPatches[patchI].clone
3325             (
3326                 newMesh.boundaryMesh(),
3327                 patchI,
3328                 patchSizes[patchI],
3329                 patchStarts[patchI]
3330             ).ptr();
3331         }
3332         newMesh.addFvPatches(newBoundary);
3333     }
3334
3335
3336     // Zones
3337     // ~~~~~
3338
3339     // Start off from empty zones.
3340     const pointZoneMesh& oldPointZones = mesh.pointZones();
3341     List<pointZone*> pZonePtrs(oldPointZones.size());
3342     {
3343         forAll(oldPointZones, i)
3344         {
3345             pZonePtrs[i] = new pointZone
3346             (
3347                 oldPointZones[i].name(),
3348                 labelList(0),
3349                 i,
3350                 newMesh.pointZones()
3351             );
3352         }
3353     }
3354
3355     const faceZoneMesh& oldFaceZones = mesh.faceZones();
3356     List<faceZone*> fZonePtrs(oldFaceZones.size());
3357     {
3358         forAll(oldFaceZones, i)
3359         {
3360             fZonePtrs[i] = new faceZone
3361             (
3362                 oldFaceZones[i].name(),
3363                 labelList(0),
3364                 boolList(0),
3365                 i,
3366                 newMesh.faceZones()
3367             );
3368         }
3369     }
3370
3371     const cellZoneMesh& oldCellZones = mesh.cellZones();
3372     List<cellZone*> cZonePtrs(oldCellZones.size());
3373     {
3374         forAll(oldCellZones, i)
3375         {
3376             cZonePtrs[i] = new cellZone
3377             (
3378                 oldCellZones[i].name(),
3379                 labelList(0),
3380                 i,
3381                 newMesh.cellZones()
3382             );
3383         }
3384     }
3385
3386     newMesh.addZones(pZonePtrs, fZonePtrs, cZonePtrs);
3387
3388     // Inverse of point/face/cell zone addressing.
3389     // For every preserved point/face/cells in zone give the old position.
3390     // For added points, the index is set to -1
3391     labelListList pointZoneMap(mesh.pointZones().size());
3392     labelListList faceZoneFaceMap(mesh.faceZones().size());
3393     labelListList cellZoneMap(mesh.cellZones().size());
3394
3395     resetZones(mesh, newMesh, pointZoneMap, faceZoneFaceMap, cellZoneMap);
3396
3397     // Clear zone info
3398     {
3399         pointZone_.clearStorage();
3400         faceZone_.clearStorage();
3401         faceZoneFlip_.clearStorage();
3402         cellZone_.clearStorage();
3403     }
3404
3405     // Patch point renumbering
3406     // For every preserved point on a patch give the old position.
3407     // For added points, the index is set to -1
3408     labelListList patchPointMap(newMesh.boundaryMesh().size());
3409     calcPatchPointMap
3410     (
3411         oldPatchMeshPointMaps,
3412         newMesh.boundaryMesh(),
3413         patchPointMap
3414     );
3415
3416     // Create the face zone mesh point renumbering
3417     labelListList faceZonePointMap(newMesh.faceZones().size());
3418     calcFaceZonePointMap(newMesh, oldFaceZoneMeshPointMaps, faceZonePointMap);
3419
3420     if (debug)
3421     {
3422         Pout<< "New mesh:" << nl;
3423         writeMeshStats(mesh, Pout);
3424     }
3425
3426     labelHashSet flipFaceFluxSet(getSetIndices(flipFaceFlux_));
3427
3428     return autoPtr<mapPolyMesh>
3429     (
3430         new mapPolyMesh
3431         (
3432             newMesh,
3433             nOldPoints,
3434             nOldFaces,
3435             nOldCells,
3436
3437             pointMap_,
3438             pointsFromPoints,
3439
3440             faceMap_,
3441             facesFromPoints,
3442             facesFromEdges,
3443             facesFromFaces,
3444
3445             cellMap_,
3446             cellsFromPoints,
3447             cellsFromEdges,
3448             cellsFromFaces,
3449             cellsFromCells,
3450
3451             reversePointMap_,
3452             reverseFaceMap_,
3453             reverseCellMap_,
3454
3455             flipFaceFluxSet,
3456
3457             patchPointMap,
3458
3459             pointZoneMap,
3460
3461             faceZonePointMap,
3462             faceZoneFaceMap,
3463             cellZoneMap,
3464
3465             newPoints,          // if empty signals no inflation.
3466             oldPatchStarts,
3467             oldPatchNMeshPoints,
3468             true                // steal storage.
3469         )
3470     );
3471
3472     // At this point all member DynamicList (pointMap_, cellMap_ etc.) will
3473     // be invalid.
3474 }
3475
3476
3477 // ************************************************************************* //