src/dynamicMesh/dynamicTopoFvMesh/dynamicTopoFvMeshReOrder.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "objectMap.H"
  27 #include "coupledInfo.H"
  28 #include "dynamicTopoFvMesh.H"
  29
  30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  31
  32 namespace Foam
  33 {
  34
  35 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  36
  37 // Reorder points after a topology change
  38 void dynamicTopoFvMesh::reOrderPoints
  39 (
  40     pointField& points,
  41     pointField& preMotionPoints,
  42     labelListList& pointZoneMap,
  43     bool threaded
  44 )
  45 {
  46     // *** Point renumbering *** //
  47     // If points were deleted during topology change, the numerical order ceases
  48     // to be continuous. Loop through all points and renumber sequentially.
  49
  50     if (debug)
  51     {
  52         if (threaded)
  53         {
  54             Info<< "Thread: " << self() << ": ";
  55         }
  56
  57         Info<< "ReOrdering points..." << flush;
  58     }
  59
  60     // Allocate for the mapping information
  61     pointMap_.setSize(nPoints_, -1);
  62
  63     label pointInOrder = 0;
  64
  65     addedPointRenumbering_.clear();
  66
  67     for (label pointI = 0; pointI < nOldPoints_; pointI++)
  68     {
  69         // Check if this is a deleted point
  70         if (reversePointMap_[pointI] == -1)
  71         {
  72             continue;
  73         }
  74
  75         // Update the point info
  76         points[pointInOrder] = points_[pointI];
  77         preMotionPoints[pointInOrder] = oldPoints_[pointI];
  78
  79         // Update maps
  80         pointMap_[pointInOrder]  = pointI;
  81         reversePointMap_[pointI] = pointInOrder;
  82
  83         // Update the counter
  84         pointInOrder++;
  85     }
  86
  87     for (label pointI = nOldPoints_; pointI < points_.size(); pointI++)
  88     {
  89         // Was this point removed after addition?
  90         if (deletedPoints_.found(pointI))
  91         {
  92             continue;
  93         }
  94
  95         // Update the point info
  96         points[pointInOrder] = points_[pointI];
  97         preMotionPoints[pointInOrder] = oldPoints_[pointI];
  98
  99         // Put inserted points in a seperate hashSet
 100         addedPointRenumbering_.insert(pointI, pointInOrder);
 101
 102         // Update the counter
 103         pointInOrder++;
 104     }
 105
 106     // Now that we're done preparing the point maps, unlock the point mutex
 107     if (threaded)
 108     {
 109         entityMutex_[0].unlock();
 110     }
 111
 112     // Final check to ensure everything went okay
 113     if (pointInOrder != nPoints_)
 114     {
 115         FatalErrorIn("dynamicTopoFvMesh::reOrderPoints()") << nl
 116             << " Algorithm did not visit every point in the mesh."
 117             << " Something's messed up." << nl
 118             << abort(FatalError);
 119     }
 120
 121     // Renumber all maps.
 122     forAll(pointsFromPoints_, indexI)
 123     {
 124         objectMap& thisMap = pointsFromPoints_[indexI];
 125
 126         if (thisMap.index() < nOldPoints_)
 127         {
 128             thisMap.index() = reversePointMap_[thisMap.index()];
 129         }
 130         else
 131         {
 132             thisMap.index() = addedPointRenumbering_[thisMap.index()];
 133         }
 134     }
 135
 136     // Prepare the pointZoneMap
 137     pointZoneMesh& pointZones = polyMesh::pointZones();
 138
 139     labelListList newPointZoneAddr(pointZones.size());
 140     labelList nPointsInZone(pointZones.size(), 0);
 141
 142     // Prepare zone maps.
 143     forAll(pointZones, pzI)
 144     {
 145         // Get the list of old points
 146         const labelList& oldAddr = pointZones[pzI];
 147
 148         label& curNPoints = nPointsInZone[pzI];
 149
 150         // First count the actual number of points in each zone.
 151         forAll(oldAddr, pointI)
 152         {
 153             // Was this zone point deleted? Don't count it.
 154             if (reversePointMap_[oldAddr[pointI]] != -1)
 155             {
 156                 curNPoints++;
 157             }
 158         }
 159
 160         // Check for added points as well
 161         forAllIter(Map<label>, addedPointZones_, pIter)
 162         {
 163             if (pIter() == pzI)
 164             {
 165                 curNPoints++;
 166             }
 167         }
 168
 169         label pIndex = 0;
 170         labelList& newAddr = newPointZoneAddr[pzI];
 171
 172         // Set the sizes first
 173         newAddr.setSize(curNPoints);
 174         pointZoneMap[pzI].setSize(curNPoints, -1);
 175
 176         // Add existing zone points which have been renumbered.
 177         forAll(oldAddr, pointI)
 178         {
 179             if (reversePointMap_[oldAddr[pointI]] != -1)
 180             {
 181                 pointZoneMap[pzI][pIndex] = oldAddr[pointI];
 182                 newAddr[pIndex] = reversePointMap_[oldAddr[pointI]];
 183                 pIndex++;
 184             }
 185         }
 186
 187         // Next, add the newly added zone points.
 188         forAllIter(Map<label>, addedPointZones_, pIter)
 189         {
 190             if (pIter() == pzI)
 191             {
 192                 newAddr[pIndex++] = addedPointRenumbering_[pIter.key()];
 193             }
 194         }
 195
 196         // Finally, assign addressing to this zone.
 197         pointZones[pzI] = newPointZoneAddr[pzI];
 198     }
 199
 200     // Reset all zones
 201     pointZones.updateMesh();
 202
 203     // Loop through all local coupling maps, and renumber points.
 204     forAll(patchCoupling_, patchI)
 205     {
 206         if (!patchCoupling_(patchI))
 207         {
 208             continue;
 209         }
 210
 211         const coupleMap& cMap = patchCoupling_[patchI].map();
 212
 213         // Obtain references
 214         Map<label>& mtsMap = cMap.entityMap(coupleMap::POINT);
 215
 216         Map<label> newMtsMap, newStmMap;
 217
 218         forAllIter(Map<label>, mtsMap, pIter)
 219         {
 220             label newMaster = -1, newSlave = -1;
 221
 222             if (pIter.key() < nOldPoints_)
 223             {
 224                 newMaster = reversePointMap_[pIter.key()];
 225             }
 226             else
 227             {
 228                 newMaster = addedPointRenumbering_[pIter.key()];
 229             }
 230
 231             if (pIter() < nOldPoints_)
 232             {
 233                 newSlave = reversePointMap_[pIter()];
 234             }
 235             else
 236             {
 237                 newSlave = addedPointRenumbering_[pIter()];
 238             }
 239
 240             // Update the map.
 241             newMtsMap.insert(newMaster, newSlave);
 242             newStmMap.insert(newSlave, newMaster);
 243         }
 244
 245         // Overwrite the old maps.
 246         cMap.transferMaps(coupleMap::POINT, newMtsMap, newStmMap);
 247     }
 248
 249     // Clear local point copies
 250     points_.clear();
 251     oldPoints_.clear();
 252
 253     // Set values
 254     points_ = points;
 255     oldPoints_ = preMotionPoints;
 256
 257     if (debug)
 258     {
 259         Info<< "Done." << endl;
 260     }
 261 }
 262
 263
 264 // Static equivalent for multi-threading
 265 void dynamicTopoFvMesh::reOrderPointsThread
 266 (
 267     void *argument
 268 )
 269 {
 270     // Recast the argument
 271     meshHandler *thread = static_cast<meshHandler*>(argument);
 272
 273     dynamicTopoFvMesh& mesh = thread->reference();
 274
 275     // Lock the point mutex first
 276     mesh.entityMutex(0).lock();
 277
 278     // Signal the calling thread
 279     thread->sendSignal(meshHandler::START);
 280
 281     // Recast the pointers for the reOrderPoints argument
 282     pointField& points =
 283     (
 284         *(static_cast<pointField*>(thread->operator()(0)))
 285     );
 286
 287     pointField& preMotionPoints =
 288     (
 289         *(static_cast<pointField*>(thread->operator()(1)))
 290     );
 291
 292     labelListList& pointZoneMap =
 293     (
 294         *(static_cast<labelListList*>(thread->operator()(2)))
 295     );
 296
 297     // Reorder the points
 298     mesh.reOrderPoints(points, preMotionPoints, pointZoneMap, true);
 299 }
 300
 301
 302 // Reorder edges after a topology change
 303 void dynamicTopoFvMesh::reOrderEdges
 304 (
 305     edgeList& edges,
 306     labelListList& edgeFaces,
 307     labelListList& faceEdges,
 308     bool threaded
 309 )
 310 {
 311     // *** Edge renumbering *** //
 312     // If edges were deleted during topology change, the numerical order ceases
 313     // to be continuous. Edges are added to respective internal/boundary patches
 314
 315     if (debug)
 316     {
 317         if (threaded)
 318         {
 319             Info<< "Thread: " << self() << ": ";
 320         }
 321
 322         Info<< "ReOrdering edges..." << flush;
 323     }
 324
 325     // Allocate for mapping information
 326     edgeMap_.setSize(nEdges_, -1);
 327
 328     label edgeInOrder = 0, allEdges = edges_.size();
 329     edgeList oldEdges(allEdges);
 330     labelListList oldEdgeFaces(allEdges);
 331
 332     addedEdgeRenumbering_.clear();
 333
 334     // Transfer old edge-based lists, and clear them
 335     forAll(edges_, edgeI)
 336     {
 337         oldEdges[edgeI] = edges_[edgeI];
 338         oldEdgeFaces[edgeI].transfer(edgeFaces_[edgeI]);
 339     }
 340
 341     edges_.setSize(nEdges_); edgeFaces_.setSize(nEdges_);
 342
 343     // Keep track of inserted boundary edge indices
 344     labelList boundaryPatchIndices(edgePatchStarts_);
 345
 346     // Loop through all edges and add them
 347     forAll(oldEdges, edgeI)
 348     {
 349         // Ensure that we're adding valid edges
 350         if (oldEdgeFaces[edgeI].empty())
 351         {
 352             continue;
 353         }
 354
 355         // Determine which patch this edge belongs to
 356         label patch = whichEdgePatch(edgeI);
 357
 358         // Obtain references
 359         edge& thisEdge = oldEdges[edgeI];
 360         labelList& thisEF = oldEdgeFaces[edgeI];
 361
 362         if (patch >= 0)
 363         {
 364             label bEdgeIndex = boundaryPatchIndices[patch]++;
 365
 366             // Update the maps
 367             if (edgeI < nOldEdges_)
 368             {
 369                 edgeMap_[bEdgeIndex] = edgeI;
 370                 reverseEdgeMap_[edgeI] = bEdgeIndex;
 371             }
 372             else
 373             {
 374                 edgeMap_[bEdgeIndex] = -1;
 375                 addedEdgeRenumbering_.insert(edgeI, bEdgeIndex);
 376             }
 377
 378             // Insert entities into local lists...
 379             edges_[bEdgeIndex] = thisEdge;
 380             edgeFaces_[bEdgeIndex] = thisEF;
 381
 382             // Insert entities into mesh-reset lists
 383             edges[bEdgeIndex] = thisEdge;
 384             edgeFaces[bEdgeIndex].transfer(thisEF);
 385         }
 386         else
 387         {
 388             // Renumber internal edges and add normally.
 389             if (edgeI < nOldEdges_)
 390             {
 391                 edgeMap_[edgeInOrder] = edgeI;
 392                 reverseEdgeMap_[edgeI] = edgeInOrder;
 393             }
 394             else
 395             {
 396                 addedEdgeRenumbering_.insert(edgeI, edgeInOrder);
 397             }
 398
 399             // Insert entities into local lists...
 400             edges_[edgeInOrder] = thisEdge;
 401             edgeFaces_[edgeInOrder] = thisEF;
 402
 403             // Insert entities into mesh-reset lists
 404             edges[edgeInOrder] = thisEdge;
 405             edgeFaces[edgeInOrder].transfer(thisEF);
 406
 407             edgeInOrder++;
 408         }
 409     }
 410
 411     // Now that we're done with edges, unlock it
 412     if (threaded)
 413     {
 414         entityMutex_[1].unlock();
 415     }
 416
 417     // Final check to ensure everything went okay
 418     if (edgeInOrder != nInternalEdges_)
 419     {
 420         FatalErrorIn("dynamicTopoFvMesh::reOrderEdges()") << nl
 421             << " Algorithm did not visit every internal edge in the mesh."
 422             << " Something's messed up." << nl
 423             << abort(FatalError);
 424     }
 425
 426     // Renumber all edges with updated point information
 427     label pIndex = -1;
 428
 429     if (threaded)
 430     {
 431         entityMutex_[0].lock();
 432     }
 433
 434     forAll(edges_, edgeI)
 435     {
 436         // Obtain references
 437         edge& thisEdge = edges_[edgeI];
 438         edge& thisREdge = edges[edgeI];
 439
 440         // Renumber edges
 441         if (thisEdge[0] < nOldPoints_)
 442         {
 443             pIndex = reversePointMap_[thisEdge[0]];
 444         }
 445         else
 446         {
 447             pIndex = addedPointRenumbering_[thisEdge[0]];
 448         }
 449
 450         thisEdge[0] = pIndex;
 451         thisREdge[0] = pIndex;
 452
 453         if (thisEdge[1] < nOldPoints_)
 454         {
 455             pIndex = reversePointMap_[thisEdge[1]];
 456         }
 457         else
 458         {
 459             pIndex = addedPointRenumbering_[thisEdge[1]];
 460         }
 461
 462         thisEdge[1] = pIndex;
 463         thisREdge[1] = pIndex;
 464     }
 465
 466     if (threaded)
 467     {
 468         entityMutex_[0].unlock();
 469     }
 470
 471     // Renumber all faceEdges
 472     label eIndex = -1;
 473
 474     if (threaded)
 475     {
 476         entityMutex_[2].lock();
 477     }
 478
 479     forAll(faceEdges_, faceI)
 480     {
 481         // Obtain references
 482         labelList& fEdges = faceEdges_[faceI];
 483         labelList& rfEdges = faceEdges[faceI];
 484
 485         forAll(fEdges, edgeI)
 486         {
 487             if (fEdges[edgeI] < nOldEdges_)
 488             {
 489                 eIndex = reverseEdgeMap_[fEdges[edgeI]];
 490             }
 491             else
 492             {
 493                 eIndex = addedEdgeRenumbering_[fEdges[edgeI]];
 494             }
 495
 496             fEdges[edgeI] = eIndex;
 497             rfEdges[edgeI] = eIndex;
 498         }
 499     }
 500
 501     // Renumber all edgeFaces
 502     label fIndex = -1;
 503
 504     forAll(edgeFaces_, edgeI)
 505     {
 506         // Obtain references
 507         labelList& eFaces = edgeFaces_[edgeI];
 508         labelList& reFaces = edgeFaces[edgeI];
 509
 510         // Renumber edgeFaces
 511         forAll(eFaces, faceI)
 512         {
 513             if (eFaces[faceI] < nOldFaces_)
 514             {
 515                 fIndex = reverseFaceMap_[eFaces[faceI]];
 516             }
 517             else
 518             {
 519                 fIndex = addedFaceRenumbering_[eFaces[faceI]];
 520             }
 521
 522             eFaces[faceI] = fIndex;
 523             reFaces[faceI] = fIndex;
 524         }
 525     }
 526
 527     if (threaded)
 528     {
 529         entityMutex_[2].unlock();
 530     }
 531
 532     if (is3D())
 533     {
 534         // Invert edges to obtain pointEdges
 535         pointEdges_ = invertManyToMany<edge, labelList>(nPoints_, edges_);
 536
 537         // Loop through all local coupling maps, and renumber edges.
 538         forAll(patchCoupling_, patchI)
 539         {
 540             if (!patchCoupling_(patchI))
 541             {
 542                 continue;
 543             }
 544
 545             // Obtain references
 546             const coupleMap& cMap = patchCoupling_[patchI].map();
 547             const Map<label>& mtsMap = cMap.entityMap(coupleMap::EDGE);
 548
 549             Map<label> newMtsMap, newStmMap;
 550
 551             forAllConstIter(Map<label>, mtsMap, mIter)
 552             {
 553                 label newMaster = -1, newSlave = -1;
 554
 555                 if (mIter.key() < nOldEdges_)
 556                 {
 557                     newMaster = reverseEdgeMap_[mIter.key()];
 558                 }
 559                 else
 560                 {
 561                     newMaster = addedEdgeRenumbering_[mIter.key()];
 562                 }
 563
 564                 if (mIter() < nOldEdges_)
 565                 {
 566                     newSlave = reverseEdgeMap_[mIter()];
 567                 }
 568                 else
 569                 {
 570                     newSlave = addedEdgeRenumbering_[mIter()];
 571                 }
 572
 573                 // Update the map.
 574                 newMtsMap.insert(newMaster, newSlave);
 575                 newStmMap.insert(newSlave, newMaster);
 576             }
 577
 578             // Overwrite the old maps.
 579             cMap.transferMaps(coupleMap::EDGE, newMtsMap, newStmMap);
 580         }
 581     }
 582
 583     if (debug)
 584     {
 585         Info<< "Done." << endl;
 586     }
 587 }
 588
 589
 590 // Static equivalent for multi-threading
 591 void dynamicTopoFvMesh::reOrderEdgesThread
 592 (
 593     void *argument
 594 )
 595 {
 596     // Recast the argument
 597     meshHandler *thread = static_cast<meshHandler*>(argument);
 598
 599     dynamicTopoFvMesh& mesh = thread->reference();
 600
 601     // Lock the edge mutex first
 602     mesh.entityMutex(1).lock();
 603
 604     // Signal the calling thread
 605     thread->sendSignal(meshHandler::START);
 606
 607     // Recast the pointers for the argument
 608     edgeList& edges =
 609     (
 610         *(static_cast<edgeList*>(thread->operator()(0)))
 611     );
 612
 613     labelListList& edgeFaces =
 614     (
 615         *(static_cast<labelListList*>(thread->operator()(1)))
 616     );
 617
 618     labelListList& faceEdges =
 619     (
 620         *(static_cast<labelListList*>(thread->operator()(2)))
 621     );
 622
 623     // Reorder the edges
 624     mesh.reOrderEdges(edges, edgeFaces, faceEdges, true);
 625
 626     // Signal the calling thread
 627     thread->sendSignal(meshHandler::STOP);
 628 }
 629
 630
 631 // Reorder faces in upper-triangular order after a topology change
 632 void dynamicTopoFvMesh::reOrderFaces
 633 (
 634     faceList& faces,
 635     labelList& owner,
 636     labelList& neighbour,
 637     labelListList& faceEdges,
 638     labelListList& faceZoneFaceMap,
 639     bool threaded
 640 )
 641 {
 642     // *** Face renumbering *** //
 643     // Faces have to be renumbered if any were added/deleted/modified
 644     // Boundary faces are added to respective patches.
 645     // Internal faces, however, have to be added in upper-triangular ordering;
 646     // i.e., in the increasing order of neighbours
 647
 648     if (debug)
 649     {
 650         if (threaded)
 651         {
 652             Info<< "Thread: " << self() << ": ";
 653         }
 654
 655         Info<< "ReOrdering faces..." << flush;
 656     }
 657
 658     // Allocate for mapping information
 659     faceMap_.setSize(nFaces_, -1);
 660
 661     label faceInOrder = 0, allFaces = faces_.size();
 662     faceList oldFaces(allFaces);
 663     labelList oldOwner(allFaces), oldNeighbour(allFaces), visited(allFaces,0);
 664     labelListList oldFaceEdges(allFaces);
 665
 666     addedFaceRenumbering_.clear();
 667
 668     // Track reverse renumbering for added faces.
 669     //  - Required during coupled patch re-ordering.
 670     Map<label> addedFaceReverseRenumbering;
 671
 672     // Make a copy of the old face-based lists, and clear them
 673     forAll(faces_, faceI)
 674     {
 675         oldFaces[faceI].transfer(faces_[faceI]);
 676         oldOwner[faceI] = owner_[faceI];
 677         oldNeighbour[faceI] = neighbour_[faceI];
 678         oldFaceEdges[faceI].transfer(faceEdges_[faceI]);
 679     }
 680
 681     // Renumber all faces with updated point information
 682     label pIndex = -1;
 683
 684     if (threaded)
 685     {
 686         entityMutex_[0].lock();
 687     }
 688
 689     forAll(oldFaces, faceI)
 690     {
 691         face& thisFace = oldFaces[faceI];
 692
 693         forAll(thisFace, pointI)
 694         {
 695             if (thisFace[pointI] < nOldPoints_)
 696             {
 697                 pIndex = reversePointMap_[thisFace[pointI]];
 698             }
 699             else
 700             {
 701                 pIndex = addedPointRenumbering_[thisFace[pointI]];
 702             }
 703
 704             thisFace[pointI] = pIndex;
 705         }
 706     }
 707
 708     if (threaded)
 709     {
 710         entityMutex_[0].unlock();
 711     }
 712
 713     // Wait for the cell mutex to become available
 714     if (threaded)
 715     {
 716         entityMutex_[3].lock();
 717     }
 718
 719     faces_.setSize(nFaces_);
 720     owner_.setSize(nFaces_);
 721     neighbour_.setSize(nFaces_);
 722     faceEdges_.setSize(nFaces_);
 723
 724     // Mark the internal faces with -2 so that they are inserted first
 725     forAll(cells_, cellI)
 726     {
 727         const cell& curFaces = cells_[cellI];
 728
 729         forAll(curFaces, faceI)
 730         {
 731             visited[curFaces[faceI]]--;
 732         }
 733     }
 734
 735     // Keep track of inserted boundary face indices
 736     labelList boundaryPatchIndices(patchStarts_);
 737
 738     // Handle boundaries first. If any coupled interfaces need to be
 739     // updated, they can be reshuffled after interior faces are done.
 740     for (label faceI = nOldInternalFaces_; faceI < allFaces; faceI++)
 741     {
 742         if (visited[faceI] == -1)
 743         {
 744             label patchID = whichPatch(faceI);
 745             label bFaceIndex = boundaryPatchIndices[patchID]++;
 746
 747             // Update the maps
 748             if (faceI < nOldFaces_)
 749             {
 750                 faceMap_[bFaceIndex] = faceI;
 751                 reverseFaceMap_[faceI] = bFaceIndex;
 752             }
 753             else
 754             {
 755                 faceMap_[bFaceIndex] = -1;
 756                 addedFaceRenumbering_.insert(faceI, bFaceIndex);
 757                 addedFaceReverseRenumbering.insert(bFaceIndex, faceI);
 758             }
 759
 760             // Renumber owner
 761             label ownerRenumber =
 762             (
 763                 oldOwner[faceI] < nOldCells_
 764               ? reverseCellMap_[oldOwner[faceI]]
 765               : addedCellRenumbering_[oldOwner[faceI]]
 766             );
 767
 768             // Insert entities into local lists...
 769             owner_[bFaceIndex] = ownerRenumber;
 770             neighbour_[bFaceIndex] = -1;
 771             faces_[bFaceIndex] = oldFaces[faceI];
 772             faceEdges_[bFaceIndex] = oldFaceEdges[faceI];
 773
 774             // Insert entities into mesh-reset lists...
 775             owner[bFaceIndex] = ownerRenumber;
 776             faces[bFaceIndex].transfer(oldFaces[faceI]);
 777             faceEdges[bFaceIndex].transfer(oldFaceEdges[faceI]);
 778
 779             // Mark this face as visited
 780             visited[faceI] = 0;
 781         }
 782     }
 783
 784     // Prepare centres and anchors for coupled interfaces
 785     List<pointField> centres(nPatches_), anchors(nPatches_);
 786
 787     // Now that points and boundary faces are up-to-date,
 788     // send and receive centres and anchor points for coupled patches
 789     initCoupledBoundaryOrdering
 790     (
 791         centres,
 792         anchors
 793     );
 794
 795     // Upper-triangular ordering of internal faces:
 796
 797     // Insertion cannot be done in one go as the faces need to be
 798     // added into the list in the increasing order of neighbour
 799     // cells.  Therefore, all neighbours will be detected first
 800     // and then added in the correct order.
 801     forAll(cells_, cellI)
 802     {
 803         // Record the neighbour cell
 804         const cell& curFaces = cells_[cellI];
 805
 806         labelList neiCells(curFaces.size(), -1);
 807
 808         label nNeighbours = 0;
 809
 810         forAll(curFaces, faceI)
 811         {
 812             if (visited[curFaces[faceI]] == -2)
 813             {
 814                 // Face is internal and gets reordered
 815                 label own =
 816                 (
 817                     oldOwner[curFaces[faceI]] < nOldCells_
 818                   ? reverseCellMap_[oldOwner[curFaces[faceI]]]
 819                   : addedCellRenumbering_[oldOwner[curFaces[faceI]]]
 820                 );
 821
 822                 label nei =
 823                 (
 824                     oldNeighbour[curFaces[faceI]] < nOldCells_
 825                   ? reverseCellMap_[oldNeighbour[curFaces[faceI]]]
 826                   : addedCellRenumbering_[oldNeighbour[curFaces[faceI]]]
 827                 );
 828
 829                 label smallerIndex = own < nei ? own : nei;
 830                 label largerIndex  = own > nei ? own : nei;
 831
 832                 if (cellI == smallerIndex)
 833                 {
 834                     neiCells[faceI] = largerIndex;
 835                     nNeighbours++;
 836                 }
 837             }
 838         }
 839
 840         // Add internal faces in the increasing order of neighbours
 841         for (label neiSearch = 0; neiSearch < nNeighbours; neiSearch++)
 842         {
 843             // Find the lowest neighbour which is still valid
 844             label nextNei = -1;
 845             label minNei = nCells_;
 846
 847             forAll(neiCells, ncI)
 848             {
 849                 if (neiCells[ncI] > -1 && neiCells[ncI] < minNei)
 850                 {
 851                     nextNei = ncI;
 852                     minNei = neiCells[ncI];
 853                 }
 854             }
 855
 856             if (nextNei > -1)
 857             {
 858                 // Face is internal and gets reordered
 859                 if (curFaces[nextNei] < nOldFaces_)
 860                 {
 861                     faceMap_[faceInOrder] = curFaces[nextNei];
 862                     reverseFaceMap_[curFaces[nextNei]] = faceInOrder;
 863                 }
 864                 else
 865                 {
 866                     addedFaceRenumbering_.insert
 867                     (
 868                         curFaces[nextNei],
 869                         faceInOrder
 870                     );
 871                 }
 872
 873                 // Renumber owner/neighbour
 874                 label oldOwn = oldOwner[curFaces[nextNei]];
 875                 label oldNei = oldNeighbour[curFaces[nextNei]];
 876
 877                 label ownerRenumber =
 878                 (
 879                     oldOwn < nOldCells_
 880                   ? reverseCellMap_[oldOwn] : addedCellRenumbering_[oldOwn]
 881                 );
 882
 883                 label neighbourRenumber =
 884                 (
 885                     oldNei < nOldCells_
 886                   ? reverseCellMap_[oldNei] : addedCellRenumbering_[oldNei]
 887                 );
 888
 889                 // Cell-reordering may cause flipped faces.. Correct them.
 890                 face& faceRenumber = oldFaces[curFaces[nextNei]];
 891
 892                 if (neighbourRenumber < ownerRenumber)
 893                 {
 894                     faceRenumber = faceRenumber.reverseFace();
 895
 896                     if (!bandWidthReduction_)
 897                     {
 898                         FatalErrorIn("dynamicTopoFvMesh::reOrderFaces()")
 899                             << nl
 900                             << " Found an improperly ordered face."
 901                             << " Something's messed up." << nl
 902                             << abort(FatalError);
 903                     }
 904
 905                     setFlip(curFaces[nextNei]);
 906                 }
 907
 908                 // Insert entities into local lists...
 909                 faces_[faceInOrder] = faceRenumber;
 910                 owner_[faceInOrder] = cellI;
 911                 neighbour_[faceInOrder] = minNei;
 912                 faceEdges_[faceInOrder] =
 913                 (
 914                     oldFaceEdges[curFaces[nextNei]]
 915                 );
 916
 917                 // Insert entities into mesh-reset lists
 918                 faces[faceInOrder].transfer(faceRenumber);
 919                 owner[faceInOrder] = cellI;
 920                 neighbour[faceInOrder] = minNei;
 921                 faceEdges[faceInOrder].transfer
 922                 (
 923                     oldFaceEdges[curFaces[nextNei]]
 924                 );
 925
 926                 // Stop the neighbour from being used again
 927                 neiCells[nextNei] = -1;
 928
 929                 // Mark this face as visited
 930                 visited[curFaces[nextNei]] = 0;
 931
 932                 faceInOrder++;
 933             }
 934             else
 935             {
 936                 FatalErrorIn("dynamicTopoFvMesh::reOrderFaces()") << nl
 937                     << "Error in internal face insertion" << nl
 938                     << abort(FatalError);
 939             }
 940         }
 941     }
 942
 943     // Prepare faceMaps and rotations for coupled interfaces
 944     labelListList patchMaps(nPatches_), rotations(nPatches_);
 945
 946     // Now compute patchMaps for coupled interfaces
 947     bool anyChange =
 948     (
 949         syncCoupledBoundaryOrdering
 950         (
 951             centres,
 952             anchors,
 953             patchMaps,
 954             rotations
 955         )
 956     );
 957
 958     if (anyChange)
 959     {
 960         forAll(patchMaps, pI)
 961         {
 962             const labelList& patchMap = patchMaps[pI];
 963
 964             if (patchMap.size())
 965             {
 966                 // Faces in this patch need to be shuffled
 967                 label pStart = patchStarts_[pI];
 968
 969                 forAll(patchMap, fI)
 970                 {
 971                     label oldPos = fI + pStart;
 972                     label newPos = patchMap[fI] + pStart;
 973
 974                     // Rotate the face in-place
 975                     const face& oldFace = faces_[oldPos];
 976
 977                     // Copy the old face, and rotate if necessary
 978                     face newFace(oldFace);
 979                     label nPos = rotations[pI][fI];
 980
 981                     if (nPos != 0)
 982                     {
 983                         forAll(oldFace, fpI)
 984                         {
 985                             label nfpI = (fpI + nPos) % oldFace.size();
 986
 987                             if (nfpI < 0)
 988                             {
 989                                 nfpI += oldFace.size();
 990                             }
 991
 992                             newFace[nfpI] = oldFace[fpI];
 993                         }
 994                     }
 995
 996                     // Map to new locations, using the
 997                     // 'old' lists as temporaries.
 998                     //  - All boundary neighbours are '-1',
 999                     //    so use that as a faceMap temporary
1000                     oldNeighbour[newPos] = faceMap_[oldPos];
1001
1002                     // Check if this face was added
1003                     if (faceMap_[oldPos] == -1)
1004                     {
1005                         // Fetch the index for the new face
1006                         label addedIndex = addedFaceReverseRenumbering[oldPos];
1007
1008                         // Renumber for the added index
1009                         addedFaceRenumbering_[addedIndex] = newPos;
1010                     }
1011                     else
1012                     {
1013                         // Fetch the index for the existing face
1014                         label existIndex = faceMap_[oldPos];
1015
1016                         // Renumber the reverse map
1017                         reverseFaceMap_[existIndex] = newPos;
1018                     }
1019
1020                     oldOwner[newPos] = owner_[oldPos];
1021                     oldFaces[newPos].transfer(newFace);
1022                     oldFaceEdges[newPos].transfer(faceEdges_[oldPos]);
1023                 }
1024
1025                 // Now copy / transfer back to original lists
1026                 forAll(patchMap, fI)
1027                 {
1028                     label index = (fI + pStart);
1029
1030                     faceMap_[index] = oldNeighbour[index];
1031
1032                     owner_[index] = oldOwner[index];
1033                     faces_[index] = oldFaces[index];
1034                     faceEdges_[index] = oldFaceEdges[index];
1035
1036                     owner[index] = oldOwner[index];
1037                     faces[index].transfer(oldFaces[index]);
1038                     faceEdges[index].transfer(oldFaceEdges[index]);
1039                 }
1040             }
1041         }
1042     }
1043
1044     // Now that we're done with faces, unlock it
1045     if (threaded)
1046     {
1047         entityMutex_[2].unlock();
1048     }
1049
1050     // Renumber all maps.
1051     forAll(facesFromPoints_, indexI)
1052     {
1053         objectMap& thisMap = facesFromPoints_[indexI];
1054
1055         if (thisMap.index() < nOldFaces_)
1056         {
1057             thisMap.index() = reverseFaceMap_[thisMap.index()];
1058         }
1059         else
1060         {
1061             thisMap.index() = addedFaceRenumbering_[thisMap.index()];
1062         }
1063     }
1064
1065     forAll(facesFromEdges_, indexI)
1066     {
1067         objectMap& thisMap = facesFromEdges_[indexI];
1068
1069         if (thisMap.index() < nOldFaces_)
1070         {
1071             thisMap.index() = reverseFaceMap_[thisMap.index()];
1072         }
1073         else
1074         {
1075             thisMap.index() = addedFaceRenumbering_[thisMap.index()];
1076         }
1077     }
1078
1079     forAll(facesFromFaces_, indexI)
1080     {
1081         objectMap& thisMap = facesFromFaces_[indexI];
1082
1083         if (thisMap.index() < nOldFaces_)
1084         {
1085             thisMap.index() = reverseFaceMap_[thisMap.index()];
1086         }
1087         else
1088         {
1089             thisMap.index() = addedFaceRenumbering_[thisMap.index()];
1090         }
1091     }
1092
1093     // Renumber all flipFaces
1094     labelHashSet flipFaces;
1095
1096     forAllIter(labelHashSet, flipFaces_, fIter)
1097     {
1098         if (fIter.key() < nOldFaces_)
1099         {
1100             flipFaces.insert(reverseFaceMap_[fIter.key()]);
1101         }
1102         else
1103         {
1104             // Added faces cannot be flipped.
1105             FatalErrorIn("dynamicTopoFvMesh::reOrderFaces()") << nl
1106                 << " Face: " << fIter.key()
1107                 << " is new, and shouldn't be flipped." << nl
1108                 << " nOldFaces: " << nOldFaces_
1109                 << abort(FatalError);
1110         }
1111     }
1112
1113     flipFaces_.transfer(flipFaces);
1114
1115     // Renumber all cells with updated face information
1116     forAll(cells_, cellI)
1117     {
1118         cell& cellFaces = cells_[cellI];
1119
1120         forAll(cellFaces, faceI)
1121         {
1122             if (cellFaces[faceI] < nOldFaces_)
1123             {
1124                 cellFaces[faceI] = reverseFaceMap_[cellFaces[faceI]];
1125             }
1126             else
1127             {
1128                 cellFaces[faceI] = addedFaceRenumbering_[cellFaces[faceI]];
1129             }
1130         }
1131     }
1132
1133     // Now that we're done with cells, unlock it
1134     if (threaded)
1135     {
1136         entityMutex_[3].unlock();
1137     }
1138
1139     // Final check to ensure everything went okay
1140     if (debug > 1)
1141     {
1142         if (sum(visited) != 0)
1143         {
1144             FatalErrorIn("dynamicTopoFvMesh::reOrderFaces()") << nl
1145                 << " Algorithm did not visit every face in the mesh."
1146                 << " Something's messed up." << nl
1147                 << abort(FatalError);
1148         }
1149     }
1150
1151     // Prepare the faceZoneMap
1152     faceZoneMesh& faceZones = polyMesh::faceZones();
1153
1154     labelListList newFaceZoneAddr(faceZones.size());
1155     boolListList faceZoneFlipMap(faceZones.size());
1156     labelList nFacesInZone(faceZones.size(), 0);
1157
1158     // Prepare zone maps.
1159     forAll(faceZones, fzI)
1160     {
1161         // Get the list of old faces
1162         const labelList& oldAddr = faceZones[fzI];
1163
1164         label& curNFaces = nFacesInZone[fzI];
1165
1166         // First count the actual number of faces in each zone.
1167         forAll(oldAddr, faceI)
1168         {
1169             // Was this zone face deleted? Don't count it.
1170             if (reverseFaceMap_[oldAddr[faceI]] != -1)
1171             {
1172                 curNFaces++;
1173             }
1174         }
1175
1176         // Check for added faces as well
1177         forAllIter(Map<label>, addedFaceZones_, fIter)
1178         {
1179             if (fIter() == fzI)
1180             {
1181                 curNFaces++;
1182             }
1183         }
1184
1185         label fIndex = 0;
1186         labelList& newAddr = newFaceZoneAddr[fzI];
1187
1188         // Set the sizes first
1189         newAddr.setSize(curNFaces);
1190         faceZoneFaceMap[fzI].setSize(curNFaces, -1);
1191         faceZoneFlipMap[fzI].setSize(curNFaces, false);
1192
1193         // Add existing zone faces which have been renumbered.
1194         forAll(oldAddr, faceI)
1195         {
1196             if (reverseFaceMap_[oldAddr[faceI]] != -1)
1197             {
1198                 faceZoneFaceMap[fzI][fIndex] = oldAddr[faceI];
1199                 newAddr[fIndex] = reverseFaceMap_[oldAddr[faceI]];
1200                 fIndex++;
1201             }
1202         }
1203
1204         // Next, add the newly added zone faces.
1205         forAllIter(Map<label>, addedFaceZones_, fIter)
1206         {
1207             if (fIter() == fzI)
1208             {
1209                 newAddr[fIndex++] = addedFaceRenumbering_[fIter.key()];
1210             }
1211         }
1212
1213         // Finally, reset addressing for this zone.
1214         faceZones[fzI].resetAddressing
1215         (
1216             newFaceZoneAddr[fzI],
1217             faceZoneFlipMap[fzI]
1218         );
1219     }
1220
1221     // Reset all zones
1222     faceZones.updateMesh();
1223
1224     // Loop through all local coupling maps, and renumber faces.
1225     forAll(patchCoupling_, patchI)
1226     {
1227         if (!patchCoupling_(patchI))
1228         {
1229             continue;
1230         }
1231
1232         // Obtain references
1233         const coupleMap& cMap = patchCoupling_[patchI].map();
1234         const Map<label>& mtsMap = cMap.entityMap(coupleMap::FACE);
1235
1236         Map<label> newMtsMap, newStmMap;
1237
1238         forAllConstIter(Map<label>, mtsMap, mIter)
1239         {
1240             label newMaster = -1, newSlave = -1;
1241
1242             if (mIter.key() < nOldFaces_)
1243             {
1244                 newMaster = reverseFaceMap_[mIter.key()];
1245             }
1246             else
1247             {
1248                 newMaster = addedFaceRenumbering_[mIter.key()];
1249             }
1250
1251             if (mIter() < nOldFaces_)
1252             {
1253                 newSlave = reverseFaceMap_[mIter()];
1254             }
1255             else
1256             {
1257                 newSlave = addedFaceRenumbering_[mIter()];
1258             }
1259
1260             // Update the map.
1261             newMtsMap.insert(newMaster, newSlave);
1262             newStmMap.insert(newSlave, newMaster);
1263         }
1264
1265         // Overwrite the old maps.
1266         cMap.transferMaps(coupleMap::FACE, newMtsMap, newStmMap);
1267     }
1268
1269     if (debug)
1270     {
1271         Info<< "Done." << endl;
1272     }
1273 }
1274
1275
1276 // Static equivalent for multi-threading
1277 void dynamicTopoFvMesh::reOrderFacesThread
1278 (
1279     void *argument
1280 )
1281 {
1282     // Recast the argument
1283     meshHandler *thread = static_cast<meshHandler*>(argument);
1284
1285     dynamicTopoFvMesh& mesh = thread->reference();
1286
1287     // Lock the face mutex first
1288     mesh.entityMutex(2).lock();
1289
1290     // Signal the calling thread
1291     thread->sendSignal(meshHandler::START);
1292
1293     // Recast the pointers for the argument
1294     faceList& faces =
1295     (
1296         *(static_cast<faceList*>(thread->operator()(0)))
1297     );
1298
1299     labelList& owner =
1300     (
1301         *(static_cast<labelList*>(thread->operator()(1)))
1302     );
1303
1304     labelList& neighbour =
1305     (
1306         *(static_cast<labelList*>(thread->operator()(2)))
1307     );
1308
1309     labelListList& faceEdges =
1310     (
1311         *(static_cast<labelListList*>(thread->operator()(3)))
1312     );
1313
1314     labelListList& faceZoneFaceMap =
1315     (
1316         *(static_cast<labelListList*>(thread->operator()(4)))
1317     );
1318
1319     // Reorder the faces
1320     mesh.reOrderFaces
1321     (
1322         faces,
1323         owner,
1324         neighbour,
1325         faceEdges,
1326         faceZoneFaceMap,
1327         true
1328     );
1329 }
1330
1331
1332 // Reorder & renumber cells with bandwidth reduction after a topology change
1333 void dynamicTopoFvMesh::reOrderCells
1334 (
1335     labelListList& cellZoneMap,
1336     bool threaded
1337 )
1338 {
1339     // *** Cell renumbering *** //
1340     // If cells were deleted during topology change, the numerical order ceases
1341     // to be continuous. Also, cells are always added at the end of the list by
1342     // virtue of the append method. Thus, cells would now have to be
1343     // reordered so that bandwidth is reduced and renumbered to be sequential.
1344
1345     if (debug)
1346     {
1347         if (threaded)
1348         {
1349             Info<< "Thread: " << self() << ": ";
1350         }
1351
1352         Info<< "ReOrdering cells..." << flush;
1353     }
1354
1355     // Allocate for mapping information
1356     cellMap_.setSize(nCells_, -1);
1357
1358     label cellInOrder = 0;
1359
1360     addedCellRenumbering_.clear();
1361
1362     // Make a copy of the old cell-based lists, and clear them
1363     label allCells = cells_.size();
1364
1365     cellList oldCells(allCells);
1366
1367     forAll(cells_, cellI)
1368     {
1369         oldCells[cellI].transfer(cells_[cellI]);
1370     }
1371
1372     cells_.setSize(nCells_);
1373
1374     if (bandWidthReduction_)
1375     {
1376         label currentCell;
1377         SLList<label> nextCell;
1378         labelList ncc(allCells, 0), visited(allCells, 0);
1379         labelListList cellCellAddr(allCells);
1380
1381         // Build a cell-cell addressing list
1382         forAll(owner_, faceI)
1383         {
1384             if ((neighbour_[faceI] > -1) && (owner_[faceI] > -1))
1385             {
1386                 ncc[owner_[faceI]]++;
1387                 ncc[neighbour_[faceI]]++;
1388             }
1389         }
1390
1391         forAll(cellCellAddr, cellI)
1392         {
1393             cellCellAddr[cellI].setSize(ncc[cellI]);
1394
1395             // Mark off deleted cells as "visited"
1396             if (ncc[cellI] == 0)
1397             {
1398                 visited[cellI] = 1;
1399             }
1400         }
1401
1402         ncc = 0;
1403
1404         forAll(owner_, faceI)
1405         {
1406             if ((owner_[faceI] > -1) && (neighbour_[faceI] > -1))
1407             {
1408                 cellCellAddr[owner_[faceI]][ncc[owner_[faceI]]++] =
1409                 (
1410                     neighbour_[faceI]
1411                 );
1412
1413                 cellCellAddr[neighbour_[faceI]][ncc[neighbour_[faceI]]++] =
1414                 (
1415                     owner_[faceI]
1416                 );
1417             }
1418         }
1419
1420         // Let's get to the "business bit" of the band-compression
1421         forAll(visited, cellI)
1422         {
1423             // Find the first cell that has not been visited yet
1424             if (visited[cellI] == 0)
1425             {
1426                 // Use this cell as a start
1427                 currentCell = cellI;
1428
1429                 nextCell.append(currentCell);
1430
1431                 // Loop through the nextCell list. Add the first cell
1432                 // into the cell order if it has not already been visited
1433                 // and ask for its neighbours. If the neighbour in question
1434                 // has not been visited, add it to the end of the nextCell list
1435                 while (nextCell.size() > 0)
1436                 {
1437                     currentCell = nextCell.removeHead();
1438
1439                     if (visited[currentCell] == 0)
1440                     {
1441                         // Mark as visited and update cell mapping info
1442                         visited[currentCell] = 1;
1443
1444                         if (currentCell < nOldCells_)
1445                         {
1446                             cellMap_[cellInOrder] = currentCell;
1447                             reverseCellMap_[currentCell] = cellInOrder;
1448                         }
1449                         else
1450                         {
1451                             addedCellRenumbering_.insert
1452                             (
1453                                 currentCell,
1454                                 cellInOrder
1455                             );
1456                         }
1457
1458                         // Insert entities into local lists...
1459                         cells_[cellInOrder].transfer(oldCells[currentCell]);
1460
1461                         cellInOrder++;
1462
1463                         // Find if the neighbours have been visited
1464                         const labelList& neighbours = cellCellAddr[currentCell];
1465
1466                         forAll(neighbours, nI)
1467                         {
1468                             if (visited[neighbours[nI]] == 0)
1469                             {
1470                                 // Not visited, add to the list
1471                                 nextCell.append(neighbours[nI]);
1472                             }
1473                         }
1474                     }
1475                 }
1476             }
1477         }
1478
1479         if (debug > 1)
1480         {
1481             if (sum(visited) != allCells)
1482             {
1483                 FatalErrorIn("dynamicTopoFvMesh::reOrderCells()") << nl
1484                     << " Algorithm did not visit every cell in the mesh."
1485                     << " Something's messed up." << nl
1486                     << abort(FatalError);
1487             }
1488         }
1489     }
1490     else
1491     {
1492         // No bandwidth reduction. Fill-in sequentially.
1493         for (label cellI = 0; cellI < nOldCells_; cellI++)
1494         {
1495             // Check if this is a deleted cell
1496             if (reverseCellMap_[cellI] == -1)
1497             {
1498                 continue;
1499             }
1500
1501             // Update the cell info
1502             cells_[cellInOrder].transfer(oldCells[cellI]);
1503
1504             // Update maps
1505             cellMap_[cellInOrder] = cellI;
1506             reverseCellMap_[cellI] = cellInOrder;
1507
1508             // Update the counter
1509             cellInOrder++;
1510         }
1511
1512         for (label cellI = nOldCells_; cellI < allCells; cellI++)
1513         {
1514             // Was this cell removed after addition?
1515             if (deletedCells_.found(cellI))
1516             {
1517                 continue;
1518             }
1519
1520             // Update the cell info
1521             cells_[cellInOrder].transfer(oldCells[cellI]);
1522
1523             // Put inserted cells in a seperate hashSet
1524             addedCellRenumbering_.insert(cellI, cellInOrder);
1525
1526             // Update the counter
1527             cellInOrder++;
1528         }
1529
1530         // Final check to ensure everything went okay
1531         if (cellInOrder != nCells_)
1532         {
1533             FatalErrorIn("dynamicTopoFvMesh::reOrderCells()") << nl
1534                 << " Algorithm did not visit every cell in the mesh."
1535                 << " Something's messed up." << nl
1536                 << abort(FatalError);
1537         }
1538     }
1539
1540     // Now that we're done preparing the cell maps, unlock the cell mutex
1541     if (threaded)
1542     {
1543         entityMutex_[3].unlock();
1544     }
1545
1546     // Renumber all maps.
1547     forAll(cellsFromPoints_, cellI)
1548     {
1549         objectMap& thisMap = cellsFromPoints_[cellI];
1550
1551         if (thisMap.index() < nOldCells_)
1552         {
1553             thisMap.index() = reverseCellMap_[thisMap.index()];
1554         }
1555         else
1556         {
1557             thisMap.index() = addedCellRenumbering_[thisMap.index()];
1558         }
1559     }
1560
1561     forAll(cellsFromEdges_, cellI)
1562     {
1563         objectMap& thisMap = cellsFromEdges_[cellI];
1564
1565         if (thisMap.index() < nOldCells_)
1566         {
1567             thisMap.index() = reverseCellMap_[thisMap.index()];
1568         }
1569         else
1570         {
1571             thisMap.index() = addedCellRenumbering_[thisMap.index()];
1572         }
1573     }
1574
1575     forAll(cellsFromFaces_, cellI)
1576     {
1577         objectMap& thisMap = cellsFromFaces_[cellI];
1578
1579         if (thisMap.index() < nOldCells_)
1580         {
1581             thisMap.index() = reverseCellMap_[thisMap.index()];
1582         }
1583         else
1584         {
1585             thisMap.index() = addedCellRenumbering_[thisMap.index()];
1586         }
1587     }
1588
1589     forAll(cellsFromCells_, cellI)
1590     {
1591         objectMap& thisMap = cellsFromCells_[cellI];
1592
1593         if (thisMap.index() < nOldCells_)
1594         {
1595             thisMap.index() = reverseCellMap_[thisMap.index()];
1596         }
1597         else
1598         {
1599             thisMap.index() = addedCellRenumbering_[thisMap.index()];
1600         }
1601     }
1602
1603     // Prepare the cellZoneMap
1604     cellZoneMesh& cellZones = polyMesh::cellZones();
1605
1606     labelListList newCellZoneAddr(cellZones.size());
1607     labelList nCellsInZone(cellZones.size(), 0);
1608
1609     // Prepare zone maps.
1610     forAll(cellZones, czI)
1611     {
1612         // Get the list of old cells
1613         const labelList& oldAddr = cellZones[czI];
1614
1615         label& curNCells = nCellsInZone[czI];
1616
1617         // First count the actual number of cells in each zone.
1618         forAll(oldAddr, cellI)
1619         {
1620             // Was this zone cell deleted? Don't count it.
1621             if (reverseCellMap_[oldAddr[cellI]] != -1)
1622             {
1623                 curNCells++;
1624             }
1625         }
1626
1627         // Check for added cells as well
1628         forAllIter(Map<label>, addedCellZones_, cIter)
1629         {
1630             if (cIter() == czI)
1631             {
1632                 curNCells++;
1633             }
1634         }
1635
1636         label cIndex = 0;
1637         labelList& newAddr = newCellZoneAddr[czI];
1638
1639         // Set the sizes first
1640         newAddr.setSize(curNCells);
1641         cellZoneMap[czI].setSize(curNCells, -1);
1642
1643         // Add existing zone cells which have been renumbered.
1644         forAll(oldAddr, cellI)
1645         {
1646             if (reverseCellMap_[oldAddr[cellI]] != -1)
1647             {
1648                 cellZoneMap[czI][cIndex] = oldAddr[cellI];
1649                 newAddr[cIndex] = reverseCellMap_[oldAddr[cellI]];
1650                 cIndex++;
1651             }
1652         }
1653
1654         // Next, add the newly added zone cells.
1655         forAllIter(Map<label>, addedCellZones_, cIter)
1656         {
1657             if (cIter() == czI)
1658             {
1659                 newAddr[cIndex++] = addedCellRenumbering_[cIter.key()];
1660             }
1661         }
1662
1663         // Finally, assign addressing to this zone.
1664         cellZones[czI] = newCellZoneAddr[czI];
1665     }
1666
1667     // Reset all zones
1668     cellZones.updateMesh();
1669
1670     if (debug)
1671     {
1672         Info<< "Done." << endl;
1673     }
1674 }
1675
1676
1677 // Static equivalent for multi-threading
1678 void dynamicTopoFvMesh::reOrderCellsThread
1679 (
1680     void *argument
1681 )
1682 {
1683     // Recast the argument
1684     meshHandler *thread = static_cast<meshHandler*>(argument);
1685
1686     dynamicTopoFvMesh& mesh = thread->reference();
1687
1688     // Lock the cell mutex first
1689     mesh.entityMutex(3).lock();
1690
1691     // Signal the calling thread
1692     thread->sendSignal(meshHandler::START);
1693
1694     // Recast the pointers for the argument
1695     labelListList& cellZoneMap =
1696     (
1697         *(static_cast<labelListList*>(thread->operator()(0)))
1698     );
1699
1700     // Reorder the cells
1701     mesh.reOrderCells(cellZoneMap, true);
1702 }
1703
1704
1705 // Reorder the faces in upper-triangular order, and generate mapping information
1706 void dynamicTopoFvMesh::reOrderMesh
1707 (
1708     pointField& points,
1709     pointField& preMotionPoints,
1710     edgeList& edges,
1711     faceList& faces,
1712     labelList& owner,
1713     labelList& neighbour,
1714     labelListList& faceEdges,
1715     labelListList& edgeFaces,
1716     labelListList& pointZoneMap,
1717     labelListList& faceZoneFaceMap,
1718     labelListList& cellZoneMap
1719 )
1720 {
1721     if (debug)
1722     {
1723         // Buffered print-out of mesh-stats
1724         if (Pstream::parRun())
1725         {
1726             if (Pstream::master())
1727             {
1728                 for (label i = Pstream::nProcs() - 1; i >= 1; i--)
1729                 {
1730                     // Notify slave
1731                     meshOps::pWrite(i, i);
1732
1733                     // Wait for response from slave
1734                     label j = 0;
1735                     meshOps::pRead(i, j);
1736                 }
1737             }
1738             else
1739             {
1740                 // Wait for master
1741                 label i = 0;
1742                 meshOps::pRead(Pstream::masterNo(), i);
1743             }
1744         }
1745
1746         Pout<< nl
1747             << "=================" << nl
1748             << " Mesh reOrdering " << nl
1749             << "=================" << nl
1750             << " Mesh Info [n]:" << nl
1751             << " Points: " << nOldPoints_ << nl
1752             << " Edges: " << nOldEdges_ << nl
1753             << " Faces: " << nOldFaces_ << nl
1754             << " Cells: " << nOldCells_ << nl
1755             << " Internal Edges: " << nOldInternalEdges_ << nl
1756             << " Internal Faces: " << nOldInternalFaces_ << nl
1757             << " nPatches: " << nPatches_ << nl;
1758
1759         if (debug > 1)
1760         {
1761             Pout<< " Patch Starts [Face]: " << oldPatchStarts_ << nl
1762                 << " Patch Sizes  [Face]: " << oldPatchSizes_ << nl
1763                 << " Patch Starts [Edge]: " << oldEdgePatchStarts_ << nl
1764                 << " Patch Sizes  [Edge]: " << oldEdgePatchSizes_ << nl;
1765         }
1766
1767         Pout<< "=================" << nl
1768             << " Mesh Info [n+1]:" << nl
1769             << " Points: " << nPoints_ << nl
1770             << " Edges: " << nEdges_ << nl
1771             << " Faces: " << nFaces_ << nl
1772             << " Cells: " << nCells_ << nl
1773             << " Internal Edges: " << nInternalEdges_ << nl
1774             << " Internal Faces: " << nInternalFaces_ << nl
1775             << " nPatches: " << nPatches_ << nl;
1776
1777         if (debug > 1)
1778         {
1779             Pout<< " Patch Starts [Face]: " << patchStarts_ << nl
1780                 << " Patch Sizes: [Face]: " << patchSizes_ << nl
1781                 << " Patch Starts [Edge]: " << edgePatchStarts_ << nl
1782                 << " Patch Sizes: [Edge]: " << edgePatchSizes_ << nl;
1783         }
1784
1785         Pout<< "=================" << endl;
1786
1787         // Notify master
1788         if (Pstream::parRun() && !Pstream::master())
1789         {
1790             meshOps::pWrite(Pstream::masterNo(), Pstream::myProcNo());
1791         }
1792
1793         bool checkPoint = false;
1794         reduce(checkPoint, andOp<bool>());
1795
1796         // Check connectivity structures for consistency
1797         // before entering the reOrdering phase.
1798         checkConnectivity();
1799     }
1800
1801     if (threader_->multiThreaded() && !Pstream::parRun())
1802     {
1803         // Initialize multi-threaded reOrdering
1804         threadedMeshReOrdering
1805         (
1806             points,
1807             preMotionPoints,
1808             edges,
1809             faces,
1810             owner,
1811             neighbour,
1812             faceEdges,
1813             edgeFaces,
1814             pointZoneMap,
1815             faceZoneFaceMap,
1816             cellZoneMap
1817         );
1818     }
1819     else
1820     {
1821         // Reorder the points
1822         reOrderPoints(points, preMotionPoints, pointZoneMap);
1823
1824         // Reorder the cells
1825         reOrderCells(cellZoneMap);
1826
1827         // Reorder the faces
1828         reOrderFaces(faces, owner, neighbour, faceEdges, faceZoneFaceMap);
1829
1830         // Reorder the edges
1831         reOrderEdges(edges, edgeFaces, faceEdges);
1832     }
1833 }
1834
1835
1836 // Invoke reOrdering with multiple threads
1837 void dynamicTopoFvMesh::threadedMeshReOrdering
1838 (
1839     pointField& points,
1840     pointField& preMotionPoints,
1841     edgeList& edges,
1842     faceList& faces,
1843     labelList& owner,
1844     labelList& neighbour,
1845     labelListList& faceEdges,
1846     labelListList& edgeFaces,
1847     labelListList& pointZoneMap,
1848     labelListList& faceZoneFaceMap,
1849     labelListList& cellZoneMap
1850 )
1851 {
1852     // For reOrdering, one handler for each reOrdering method
1853     PtrList<meshHandler> reOrderPtr(4);
1854
1855     // Initialize reOrdering handlers
1856     forAll(reOrderPtr, memberI)
1857     {
1858         reOrderPtr.set
1859         (
1860             memberI,
1861             new meshHandler(*this, threader())
1862         );
1863     }
1864
1865     // Set argument sizes for individual members
1866
1867     // Points take three arguments
1868     // (Two pointFields and one labelListList)
1869     reOrderPtr[0].setSize(3);
1870
1871     // Prepare pointers for point reOrdering
1872     reOrderPtr[0].set(0, &points);
1873     reOrderPtr[0].set(1, &preMotionPoints);
1874     reOrderPtr[0].set(2, &pointZoneMap);
1875
1876     // Edges take three arguments
1877     // (One edgeList and two labelListLists)
1878     reOrderPtr[1].setSize(3);
1879
1880     // Prepare pointers for edge reOrdering
1881     reOrderPtr[1].set(0, &edges);
1882     reOrderPtr[1].set(1, &edgeFaces);
1883     reOrderPtr[1].set(2, &faceEdges);
1884
1885     // Faces take five arguments
1886     // (One faceList, two labelLists, and two labelListLists)
1887     reOrderPtr[2].setSize(5);
1888
1889     // Prepare pointers for face reOrdering
1890     reOrderPtr[2].set(0, &faces);
1891     reOrderPtr[2].set(1, &owner);
1892     reOrderPtr[2].set(2, &neighbour);
1893     reOrderPtr[2].set(3, &faceEdges);
1894     reOrderPtr[2].set(4, &faceZoneFaceMap);
1895
1896     // Cells take one argument
1897     // (One labelListList)
1898     reOrderPtr[3].setSize(1);
1899
1900     // Prepare pointers for cell reOrdering
1901     reOrderPtr[3].set(0, &cellZoneMap);
1902
1903     // Set the thread scheduling sequence
1904     labelList reOrderSeq(4, -1);
1905
1906     // Points, cells, faces and edges (in that order)
1907     reOrderSeq[0] = 0;
1908     reOrderSeq[1] = 3;
1909     reOrderSeq[2] = 2;
1910     reOrderSeq[3] = 1;
1911
1912     // Lock all slave threads first
1913     forAll(reOrderSeq, i)
1914     {
1915         reOrderPtr[reOrderSeq[i]].lock(meshHandler::START);
1916         reOrderPtr[reOrderSeq[i]].unsetPredicate(meshHandler::START);
1917     }
1918
1919     // Submit points to the work queue
1920     threader_->addToWorkQueue
1921     (
1922         &reOrderPointsThread,
1923         &(reOrderPtr[0])
1924     );
1925
1926     // Wait for a signal from this thread before moving on.
1927     reOrderPtr[0].waitForSignal(meshHandler::START);
1928
1929     // Submit cells to the work queue
1930     threader_->addToWorkQueue
1931     (
1932         &reOrderCellsThread,
1933         &(reOrderPtr[3])
1934     );
1935
1936     // Wait for a signal from this thread before moving on.
1937     reOrderPtr[3].waitForSignal(meshHandler::START);
1938
1939     // Submit faces to the work queue
1940     threader_->addToWorkQueue
1941     (
1942         &reOrderFacesThread,
1943         &(reOrderPtr[2])
1944     );
1945
1946     // Wait for a signal from this thread before moving on.
1947     reOrderPtr[2].waitForSignal(meshHandler::START);
1948
1949     // Lock the edge stop-mutex
1950     reOrderPtr[1].lock(meshHandler::STOP);
1951     reOrderPtr[1].unsetPredicate(meshHandler::STOP);
1952
1953     // Submit edges to the work queue
1954     threader_->addToWorkQueue
1955     (
1956         &reOrderEdgesThread,
1957         &(reOrderPtr[1])
1958     );
1959
1960     // Wait for a signal from this thread before moving on.
1961     reOrderPtr[1].waitForSignal(meshHandler::START);
1962
1963     // Wait for edges to be reOrdered before moving on.
1964     reOrderPtr[1].waitForSignal(meshHandler::STOP);
1965 }
1966
1967 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1968
1969 } // End namespace Foam
1970
1971 // ************************************************************************* //