src/dynamicMesh/dynamicFvMesh/fvMeshDistribute/fvMeshDistribute.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 "CompactListList_dev.H"
  28 #include "fvMeshDistribute.H"
  29 #include "PstreamCombineReduceOps.H"
  30 #include "fvMeshAdder.H"
  31 #include "faceCoupleInfo.H"
  32 #include "processorFvPatchField.H"
  33 #include "processorFvsPatchField.H"
  34 #include "directTopoChange.H"
  35 #include "removeCells.H"
  36 #include "polyModifyFace.H"
  37 #include "polyRemovePoint.H"
  38 #include "mergePoints.H"
  39 #include "mapDistributePolyMesh.H"
  40 #include "surfaceFields.H"
  41 #include "syncTools.H"
  42
  43 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  44
  45 defineTypeNameAndDebug(Foam::fvMeshDistribute, 0);
  46
  47
  48 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  49
  50 Foam::labelList Foam::fvMeshDistribute::select
  51 (
  52     const bool selectEqual,
  53     const labelList& values,
  54     const label value
  55 )
  56 {
  57     label n = 0;
  58
  59     forAll(values, i)
  60     {
  61         if (selectEqual == (values[i] == value))
  62         {
  63             n++;
  64         }
  65     }
  66
  67     labelList indices(n);
  68     n = 0;
  69
  70     forAll(values, i)
  71     {
  72         if (selectEqual == (values[i] == value))
  73         {
  74             indices[n++] = i;
  75         }
  76     }
  77     return indices;
  78 }
  79
  80
  81 // Check all procs have same names and in exactly same order.
  82 void Foam::fvMeshDistribute::checkEqualWordList
  83 (
  84     const string& msg,
  85     const wordList& lst
  86 )
  87 {
  88     List<wordList> allNames(Pstream::nProcs());
  89     allNames[Pstream::myProcNo()] = lst;
  90     Pstream::gatherList(allNames);
  91     Pstream::scatterList(allNames);
  92
  93     for (label procI = 1; procI < Pstream::nProcs(); procI++)
  94     {
  95         if (allNames[procI] != allNames[0])
  96         {
  97             FatalErrorIn("fvMeshDistribute::checkEqualWordList(..)")
  98                 << "When checking for equal " << msg.c_str() << " :" << endl
  99                 << "processor0 has:" << allNames[0] << endl
 100                 << "processor" << procI << " has:" << allNames[procI] << endl
 101                 << msg.c_str() << " need to be synchronised on all processors."
 102                 << exit(FatalError);
 103         }
 104     }
 105 }
 106
 107
 108 Foam::wordList Foam::fvMeshDistribute::mergeWordList(const wordList& procNames)
 109 {
 110     List<wordList> allNames(Pstream::nProcs());
 111     allNames[Pstream::myProcNo()] = procNames;
 112     Pstream::gatherList(allNames);
 113     Pstream::scatterList(allNames);
 114
 115     HashSet<word> mergedNames;
 116     forAll(allNames, procI)
 117     {
 118         forAll(allNames[procI], i)
 119         {
 120             mergedNames.insert(allNames[procI][i]);
 121         }
 122     }
 123     return mergedNames.toc();
 124 }
 125
 126
 127 // Print some info on mesh.
 128 void Foam::fvMeshDistribute::printMeshInfo(const fvMesh& mesh)
 129 {
 130     Pout<< "Primitives:" << nl
 131         << "    points       :" << mesh.nPoints() << nl
 132         << "    internalFaces:" << mesh.nInternalFaces() << nl
 133         << "    faces        :" << mesh.nFaces() << nl
 134         << "    cells        :" << mesh.nCells() << nl;
 135
 136     const fvBoundaryMesh& patches = mesh.boundary();
 137
 138     Pout<< "Patches:" << endl;
 139     forAll(patches, patchI)
 140     {
 141         const polyPatch& pp = patches[patchI].patch();
 142
 143         Pout<< "    " << patchI << " name:" << pp.name()
 144             << " size:" << pp.size()
 145             << " start:" << pp.start()
 146             << " type:" << pp.type()
 147             << endl;
 148     }
 149
 150     if (mesh.pointZones().size())
 151     {
 152         Pout<< "PointZones:" << endl;
 153         forAll(mesh.pointZones(), zoneI)
 154         {
 155             const pointZone& pz = mesh.pointZones()[zoneI];
 156             Pout<< "    " << zoneI << " name:" << pz.name()
 157                 << " size:" << pz.size()
 158                 << endl;
 159         }
 160     }
 161     if (mesh.faceZones().size())
 162     {
 163         Pout<< "FaceZones:" << endl;
 164         forAll(mesh.faceZones(), zoneI)
 165         {
 166             const faceZone& fz = mesh.faceZones()[zoneI];
 167             Pout<< "    " << zoneI << " name:" << fz.name()
 168                 << " size:" << fz.size()
 169                 << endl;
 170         }
 171     }
 172     if (mesh.cellZones().size())
 173     {
 174         Pout<< "CellZones:" << endl;
 175         forAll(mesh.cellZones(), zoneI)
 176         {
 177             const cellZone& cz = mesh.cellZones()[zoneI];
 178             Pout<< "    " << zoneI << " name:" << cz.name()
 179                 << " size:" << cz.size()
 180                 << endl;
 181         }
 182     }
 183 }
 184
 185
 186 void Foam::fvMeshDistribute::printCoupleInfo
 187 (
 188     const primitiveMesh& mesh,
 189     const labelList& sourceFace,
 190     const labelList& sourceProc,
 191     const labelList& sourceNewProc
 192 )
 193 {
 194     Pout<< nl
 195         << "Current coupling info:"
 196         << endl;
 197
 198     forAll(sourceFace, bFaceI)
 199     {
 200         label meshFaceI = mesh.nInternalFaces() + bFaceI;
 201
 202         Pout<< "    meshFace:" << meshFaceI
 203             << " fc:" << mesh.faceCentres()[meshFaceI]
 204             << " connects to proc:" << sourceProc[bFaceI]
 205             << "/face:" << sourceFace[bFaceI]
 206             << " which will move to proc:" << sourceNewProc[bFaceI]
 207             << endl;
 208     }
 209 }
 210
 211
 212 // Finds (non-empty) patch that exposed internal and proc faces can be put into
 213 Foam::label Foam::fvMeshDistribute::findNonEmptyPatch() const
 214 {
 215     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 216
 217     label nonEmptyPatchI = -1;
 218
 219     forAllReverse(patches, patchI)
 220     {
 221         const polyPatch& pp = patches[patchI];
 222
 223         if (!isA<emptyPolyPatch>(pp) && !pp.coupled())
 224         {
 225             nonEmptyPatchI = patchI;
 226             break;
 227         }
 228     }
 229
 230     if (nonEmptyPatchI == -1)
 231     {
 232         FatalErrorIn("fvMeshDistribute::findNonEmptyPatch() const")
 233             << "Cannot find a patch which is neither of type empty nor"
 234             << " coupled in patches " << patches.names() << endl
 235             << "There has to be at least one such patch for"
 236             << " distribution to work" << abort(FatalError);
 237     }
 238
 239     if (debug)
 240     {
 241         Pout<< "findNonEmptyPatch : using patch " << nonEmptyPatchI
 242             << " name:" << patches[nonEmptyPatchI].name()
 243             << " type:" << patches[nonEmptyPatchI].type()
 244             << " to put exposed faces into." << endl;
 245     }
 246
 247
 248     // Do additional test for processor patches intermingled with non-proc
 249     // patches.
 250     label procPatchI = -1;
 251
 252     forAll(patches, patchI)
 253     {
 254         if (isA<processorPolyPatch>(patches[patchI]))
 255         {
 256             procPatchI = patchI;
 257         }
 258         else if (procPatchI != -1)
 259         {
 260             FatalErrorIn("fvMeshDistribute::findNonEmptyPatch() const")
 261                 << "Processor patches should be at end of patch list."
 262                 << endl
 263                 << "Have processor patch " << procPatchI
 264                 << " followed by non-processor patch " << patchI
 265                 << " in patches " << patches.names()
 266                 << abort(FatalError);
 267         }
 268     }
 269
 270     return nonEmptyPatchI;
 271 }
 272
 273
 274 // Appends processorPolyPatch. Returns patchID.
 275 Foam::label Foam::fvMeshDistribute::addProcPatch
 276 (
 277     const word& patchName,
 278     const label nbrProc
 279 )
 280 {
 281     // Clear local fields and e.g. polyMesh globalMeshData.
 282     mesh_.clearOut();
 283
 284
 285     polyBoundaryMesh& polyPatches =
 286         const_cast<polyBoundaryMesh&>(mesh_.boundaryMesh());
 287     fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh_.boundary());
 288
 289     if (polyPatches.findPatchID(patchName) != -1)
 290     {
 291         FatalErrorIn
 292         (
 293             "fvMeshDistribute::addProcPatch(const word&, const label)"
 294         )   << "Cannot create patch " << patchName
 295             << " since it already exists." << nl
 296             << "Current patch names:" << polyPatches.names()
 297             << exit(FatalError);
 298     }
 299
 300
 301
 302     // Add the patch
 303     // ~~~~~~~~~~~~~
 304
 305     label sz = polyPatches.size();
 306
 307     // Add polyPatch
 308     polyPatches.setSize(sz+1);
 309     polyPatches.set
 310     (
 311         sz,
 312         new processorPolyPatch
 313         (
 314             patchName,
 315             0,              // size
 316             mesh_.nFaces(),
 317             sz,
 318             mesh_.boundaryMesh(),
 319             Pstream::myProcNo(),
 320             nbrProc
 321         )
 322     );
 323     fvPatches.setSize(sz+1);
 324     fvPatches.set
 325     (
 326         sz,
 327         fvPatch::New
 328         (
 329             polyPatches[sz],  // point to newly added polyPatch
 330             mesh_.boundary()
 331         )
 332     );
 333
 334     return sz;
 335 }
 336
 337
 338 // Deletes last patch
 339 void Foam::fvMeshDistribute::deleteTrailingPatch()
 340 {
 341     // Clear local fields and e.g. polyMesh globalMeshData.
 342     mesh_.clearOut();
 343
 344     polyBoundaryMesh& polyPatches =
 345         const_cast<polyBoundaryMesh&>(mesh_.boundaryMesh());
 346     fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh_.boundary());
 347
 348     if (polyPatches.empty())
 349     {
 350         FatalErrorIn("fvMeshDistribute::deleteTrailingPatch(fvMesh&)")
 351             << "No patches in mesh"
 352             << abort(FatalError);
 353     }
 354
 355     label sz = polyPatches.size();
 356
 357     label nFaces = polyPatches[sz-1].size();
 358
 359     // Remove last polyPatch
 360     if (debug)
 361     {
 362         Pout<< "deleteTrailingPatch : Removing patch " << sz-1
 363             << " : " << polyPatches[sz-1].name() << " size:" << nFaces << endl;
 364     }
 365
 366     if (nFaces)
 367     {
 368         FatalErrorIn("fvMeshDistribute::deleteTrailingPatch()")
 369             << "There are still " << nFaces << " faces in patch to be deleted "
 370             << sz-1 << ' ' << polyPatches[sz-1].name()
 371             << abort(FatalError);
 372     }
 373
 374     // Remove actual patch
 375     polyPatches.setSize(sz-1);
 376     fvPatches.setSize(sz-1);
 377 }
 378
 379
 380 // Delete all processor patches. Move any processor faces into the last
 381 // non-processor patch.
 382 Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::deleteProcPatches
 383 (
 384     const label destinationPatch
 385 )
 386 {
 387     // New patchID per boundary faces to be repatched. Is -1 (no change)
 388     // or new patchID
 389     labelList newPatchID(mesh_.nFaces() - mesh_.nInternalFaces(), -1);
 390
 391     label nProcPatches = 0;
 392
 393     forAll(mesh_.boundaryMesh(), patchI)
 394     {
 395         const polyPatch& pp = mesh_.boundaryMesh()[patchI];
 396
 397         if (isA<processorPolyPatch>(pp))
 398         {
 399             if (debug)
 400             {
 401                 Pout<< "Moving all faces of patch " << pp.name()
 402                     << " into patch " << destinationPatch
 403                     << endl;
 404             }
 405
 406             label offset = pp.start() - mesh_.nInternalFaces();
 407
 408             forAll(pp, i)
 409             {
 410                 newPatchID[offset+i] = destinationPatch;
 411             }
 412
 413             nProcPatches++;
 414         }
 415     }
 416
 417     // Note: order of boundary faces been kept the same since the
 418     // destinationPatch is at the end and we have visited the patches in
 419     // incremental order.
 420     labelListList dummyFaceMaps;
 421     autoPtr<mapPolyMesh> map = repatch(newPatchID, dummyFaceMaps);
 422
 423
 424     // Delete (now empty) processor patches.
 425     forAllReverse(mesh_.boundaryMesh(), patchI)
 426     {
 427         const polyPatch& pp = mesh_.boundaryMesh()[patchI];
 428
 429         if (isA<processorPolyPatch>(pp))
 430         {
 431             deleteTrailingPatch();
 432             deleteTrailingPatchFields<volScalarField>();
 433             deleteTrailingPatchFields<volVectorField>();
 434             deleteTrailingPatchFields<volSphericalTensorField>();
 435             deleteTrailingPatchFields<volSymmTensorField>();
 436             deleteTrailingPatchFields<volTensorField>();
 437
 438             deleteTrailingPatchFields<surfaceScalarField>();
 439             deleteTrailingPatchFields<surfaceVectorField>();
 440             deleteTrailingPatchFields<surfaceSphericalTensorField>();
 441             deleteTrailingPatchFields<surfaceSymmTensorField>();
 442             deleteTrailingPatchFields<surfaceTensorField>();
 443         }
 444     }
 445
 446     return map;
 447 }
 448
 449
 450 // Repatch the mesh.
 451 Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::repatch
 452 (
 453     const labelList& newPatchID,         // per boundary face -1 or new patchID
 454     labelListList& constructFaceMap
 455 )
 456 {
 457     directTopoChange meshMod(mesh_);
 458
 459     forAll(newPatchID, bFaceI)
 460     {
 461         if (newPatchID[bFaceI] != -1)
 462         {
 463             label faceI = mesh_.nInternalFaces() + bFaceI;
 464
 465             label zoneID = mesh_.faceZones().whichZone(faceI);
 466             bool zoneFlip = false;
 467
 468             if (zoneID >= 0)
 469             {
 470                 const faceZone& fZone = mesh_.faceZones()[zoneID];
 471                 zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
 472             }
 473
 474             meshMod.setAction
 475             (
 476                 polyModifyFace
 477                 (
 478                     mesh_.faces()[faceI],       // modified face
 479                     faceI,                      // label of face
 480                     mesh_.faceOwner()[faceI],   // owner
 481                     -1,                         // neighbour
 482                     false,                      // face flip
 483                     newPatchID[bFaceI],         // patch for face
 484                     false,                      // remove from zone
 485                     zoneID,                     // zone for face
 486                     zoneFlip                    // face flip in zone
 487                 )
 488             );
 489         }
 490     }
 491
 492
 493     // Do mapping of fields from one patchField to the other ourselves since
 494     // is currently not supported by updateMesh.
 495
 496     // Store boundary fields (we only do this for surfaceFields)
 497     PtrList<FieldField<fvsPatchField, scalar> > sFlds;
 498     saveBoundaryFields<scalar, surfaceMesh>(sFlds);
 499     PtrList<FieldField<fvsPatchField, vector> > vFlds;
 500     saveBoundaryFields<vector, surfaceMesh>(vFlds);
 501     PtrList<FieldField<fvsPatchField, sphericalTensor> > sptFlds;
 502     saveBoundaryFields<sphericalTensor, surfaceMesh>(sptFlds);
 503     PtrList<FieldField<fvsPatchField, symmTensor> > sytFlds;
 504     saveBoundaryFields<symmTensor, surfaceMesh>(sytFlds);
 505     PtrList<FieldField<fvsPatchField, tensor> > tFlds;
 506     saveBoundaryFields<tensor, surfaceMesh>(tFlds);
 507
 508     // Change the mesh (no inflation). Note: parallel comms allowed.
 509     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
 510
 511     // Update fields. No inflation, parallel sync.
 512     mesh_.updateMesh(map);
 513
 514     // Map patch fields using stored boundary fields. Note: assumes order
 515     // of fields has not changed in object registry!
 516     mapBoundaryFields<scalar, surfaceMesh>(map, sFlds);
 517     mapBoundaryFields<vector, surfaceMesh>(map, vFlds);
 518     mapBoundaryFields<sphericalTensor, surfaceMesh>(map, sptFlds);
 519     mapBoundaryFields<symmTensor, surfaceMesh>(map, sytFlds);
 520     mapBoundaryFields<tensor, surfaceMesh>(map, tFlds);
 521
 522
 523     // Move mesh (since morphing does not do this)
 524     if (map().hasMotionPoints())
 525     {
 526         mesh_.movePoints(map().preMotionPoints());
 527     }
 528
 529     // Adapt constructMaps.
 530
 531     if (debug)
 532     {
 533         label index = findIndex(map().reverseFaceMap(), -1);
 534
 535         if (index != -1)
 536         {
 537             FatalErrorIn
 538             (
 539                 "fvMeshDistribute::repatch(const labelList&, labelListList&)"
 540             )   << "reverseFaceMap contains -1 at index:"
 541                 << index << endl
 542                 << "This means that the repatch operation was not just"
 543                 << " a shuffle?" << abort(FatalError);
 544         }
 545     }
 546
 547     forAll(constructFaceMap, procI)
 548     {
 549         inplaceRenumber(map().reverseFaceMap(), constructFaceMap[procI]);
 550     }
 551
 552
 553     return map;
 554 }
 555
 556
 557 // Detect shared points. Need processor patches to be present.
 558 // Background: when adding bits of mesh one can get points which
 559 // share the same position but are only detectable to be topologically
 560 // the same point when doing parallel analysis. This routine will
 561 // merge those points.
 562 Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::mergeSharedPoints
 563 (
 564     labelListList& constructPointMap
 565 )
 566 {
 567     // Find out which sets of points get merged and create a map from
 568     // mesh point to unique point.
 569     Map<label> pointToMaster
 570     (
 571         fvMeshAdder::findSharedPoints
 572         (
 573             mesh_,
 574             mergeTol_
 575         )
 576     );
 577
 578     if (returnReduce(pointToMaster.size(), sumOp<label>()) == 0)
 579     {
 580         return autoPtr<mapPolyMesh>(NULL);
 581     }
 582
 583     directTopoChange meshMod(mesh_);
 584
 585     fvMeshAdder::mergePoints(mesh_, pointToMaster, meshMod);
 586
 587     // Change the mesh (no inflation). Note: parallel comms allowed.
 588     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
 589
 590     // Update fields. No inflation, parallel sync.
 591     mesh_.updateMesh(map);
 592
 593     // Adapt constructMaps for merged points.
 594     forAll(constructPointMap, procI)
 595     {
 596         labelList& constructMap = constructPointMap[procI];
 597
 598         forAll(constructMap, i)
 599         {
 600             label oldPointI = constructMap[i];
 601
 602             label newPointI = map().reversePointMap()[oldPointI];
 603
 604             if (newPointI < -1)
 605             {
 606                 constructMap[i] = -newPointI-2;
 607             }
 608             else if (newPointI >= 0)
 609             {
 610                 constructMap[i] = newPointI;
 611             }
 612             else
 613             {
 614                 FatalErrorIn("fvMeshDistribute::mergeSharedPoints()")
 615                     << "Problem. oldPointI:" << oldPointI
 616                     << " newPointI:" << newPointI << abort(FatalError);
 617             }
 618         }
 619     }
 620     return map;
 621 }
 622
 623
 624 // Construct the local environment of all boundary faces.
 625 void Foam::fvMeshDistribute::getNeighbourData
 626 (
 627     const labelList& distribution,
 628     labelList& sourceFace,
 629     labelList& sourceProc,
 630     labelList& sourceNewProc
 631 ) const
 632 {
 633     label nBnd = mesh_.nFaces() - mesh_.nInternalFaces();
 634     sourceFace.setSize(nBnd);
 635     sourceProc.setSize(nBnd);
 636     sourceNewProc.setSize(nBnd);
 637
 638     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 639
 640     // Get neighbouring meshFace labels and new processor of coupled boundaries
 641     labelList nbrFaces(nBnd, -1);
 642     labelList nbrNewProc(nBnd, -1);
 643
 644     forAll(patches, patchI)
 645     {
 646         const polyPatch& pp = patches[patchI];
 647
 648         if (isA<processorPolyPatch>(pp))
 649         {
 650             label offset = pp.start() - mesh_.nInternalFaces();
 651
 652             // Mesh labels of faces on this side
 653             forAll(pp, i)
 654             {
 655                 label bndI = offset + i;
 656                 nbrFaces[bndI] = pp.start()+i;
 657             }
 658
 659             // Which processor they will end up on
 660             SubList<label>(nbrNewProc, pp.size(), offset).assign
 661             (
 662                 UIndirectList<label>(distribution, pp.faceCells())()
 663             );
 664         }
 665     }
 666
 667
 668     // Exchange the boundary data
 669     syncTools::swapBoundaryFaceList(mesh_, nbrFaces, false);
 670     syncTools::swapBoundaryFaceList(mesh_, nbrNewProc, false);
 671
 672
 673     forAll(patches, patchI)
 674     {
 675         const polyPatch& pp = patches[patchI];
 676         label offset = pp.start() - mesh_.nInternalFaces();
 677
 678         if (isA<processorPolyPatch>(pp))
 679         {
 680             const processorPolyPatch& procPatch =
 681                 refCast<const processorPolyPatch>(pp);
 682
 683             // Check which of the two faces we store.
 684
 685             if (Pstream::myProcNo() < procPatch.neighbProcNo())
 686             {
 687                 // Use my local face labels
 688                 forAll(pp, i)
 689                 {
 690                     label bndI = offset + i;
 691                     sourceFace[bndI] = pp.start()+i;
 692                     sourceProc[bndI] = Pstream::myProcNo();
 693                     sourceNewProc[bndI] = nbrNewProc[bndI];
 694                 }
 695             }
 696             else
 697             {
 698                 // Use my neighbours face labels
 699                 forAll(pp, i)
 700                 {
 701                     label bndI = offset + i;
 702                     sourceFace[bndI] = nbrFaces[bndI];
 703                     sourceProc[bndI] = procPatch.neighbProcNo();
 704                     sourceNewProc[bndI] = nbrNewProc[bndI];
 705                 }
 706             }
 707         }
 708         else
 709         {
 710             // Normal (physical) boundary
 711             forAll(pp, i)
 712             {
 713                 label bndI = offset + i;
 714                 sourceFace[bndI] = patchI;
 715                 sourceProc[bndI] = -1;
 716                 sourceNewProc[bndI] = -1;
 717             }
 718         }
 719     }
 720 }
 721
 722
 723 // Subset the neighbourCell/neighbourProc fields
 724 void Foam::fvMeshDistribute::subsetBoundaryData
 725 (
 726     const fvMesh& mesh,
 727     const labelList& faceMap,
 728     const labelList& cellMap,
 729
 730     const labelList& oldDistribution,
 731     const labelList& oldFaceOwner,
 732     const labelList& oldFaceNeighbour,
 733     const label oldInternalFaces,
 734
 735     const labelList& sourceFace,
 736     const labelList& sourceProc,
 737     const labelList& sourceNewProc,
 738
 739     labelList& subFace,
 740     labelList& subProc,
 741     labelList& subNewProc
 742 )
 743 {
 744     subFace.setSize(mesh.nFaces() - mesh.nInternalFaces());
 745     subProc.setSize(mesh.nFaces() - mesh.nInternalFaces());
 746     subNewProc.setSize(mesh.nFaces() - mesh.nInternalFaces());
 747
 748     forAll(subFace, newBFaceI)
 749     {
 750         label newFaceI = newBFaceI + mesh.nInternalFaces();
 751
 752         label oldFaceI = faceMap[newFaceI];
 753
 754         // Was oldFaceI internal face? If so which side did we get.
 755         if (oldFaceI < oldInternalFaces)
 756         {
 757             subFace[newBFaceI] = oldFaceI;
 758             subProc[newBFaceI] = Pstream::myProcNo();
 759
 760             label oldOwn = oldFaceOwner[oldFaceI];
 761             label oldNei = oldFaceNeighbour[oldFaceI];
 762
 763             if (oldOwn == cellMap[mesh.faceOwner()[newFaceI]])
 764             {
 765                 // We kept the owner side. Where does the neighbour move to?
 766                 subNewProc[newBFaceI] = oldDistribution[oldNei];
 767             }
 768             else
 769             {
 770                 // We kept the neighbour side.
 771                 subNewProc[newBFaceI] = oldDistribution[oldOwn];
 772             }
 773         }
 774         else
 775         {
 776             // Was boundary face. Take over boundary information
 777             label oldBFaceI = oldFaceI - oldInternalFaces;
 778
 779             subFace[newBFaceI] = sourceFace[oldBFaceI];
 780             subProc[newBFaceI] = sourceProc[oldBFaceI];
 781             subNewProc[newBFaceI] = sourceNewProc[oldBFaceI];
 782         }
 783     }
 784 }
 785
 786
 787 // Find cells on mesh whose faceID/procID match the neighbour cell/proc of
 788 // domainMesh. Store the matching face.
 789 void Foam::fvMeshDistribute::findCouples
 790 (
 791     const primitiveMesh& mesh,
 792     const labelList& sourceFace,
 793     const labelList& sourceProc,
 794
 795     const label domain,
 796     const primitiveMesh& domainMesh,
 797     const labelList& domainFace,
 798     const labelList& domainProc,
 799
 800     labelList& masterCoupledFaces,
 801     labelList& slaveCoupledFaces
 802 )
 803 {
 804     // Store domain neighbour as map so we can easily look for pair
 805     // with same face+proc.
 806     HashTable<label, labelPair, labelPair::Hash<> > map(domainFace.size());
 807
 808     forAll(domainFace, bFaceI)
 809     {
 810         map.insert(labelPair(domainFace[bFaceI], domainProc[bFaceI]), bFaceI);
 811     }
 812
 813
 814     // Try to match mesh data.
 815
 816     masterCoupledFaces.setSize(domainFace.size());
 817     slaveCoupledFaces.setSize(domainFace.size());
 818     label coupledI = 0;
 819
 820     forAll(sourceFace, bFaceI)
 821     {
 822         if (sourceProc[bFaceI] != -1)
 823         {
 824             labelPair myData(sourceFace[bFaceI], sourceProc[bFaceI]);
 825
 826             HashTable<label, labelPair, labelPair::Hash<> >::const_iterator
 827                 iter = map.find(myData);
 828
 829             if (iter != map.end())
 830             {
 831                 label nbrBFaceI = iter();
 832
 833                 masterCoupledFaces[coupledI] = mesh.nInternalFaces() + bFaceI;
 834                 slaveCoupledFaces[coupledI] =
 835                     domainMesh.nInternalFaces()
 836                   + nbrBFaceI;
 837
 838                 coupledI++;
 839             }
 840         }
 841     }
 842
 843     masterCoupledFaces.setSize(coupledI);
 844     slaveCoupledFaces.setSize(coupledI);
 845
 846     if (debug)
 847     {
 848         Pout<< "findCouples : found " << coupledI
 849             << " faces that will be stitched" << nl << endl;
 850     }
 851 }
 852
 853
 854 // Map data on boundary faces to new mesh (resulting from adding two meshes)
 855 Foam::labelList Foam::fvMeshDistribute::mapBoundaryData
 856 (
 857     const primitiveMesh& mesh,      // mesh after adding
 858     const mapAddedPolyMesh& map,
 859     const labelList& boundaryData0, // mesh before adding
 860     const label nInternalFaces1,
 861     const labelList& boundaryData1  // added mesh
 862 )
 863 {
 864     labelList newBoundaryData(mesh.nFaces() - mesh.nInternalFaces());
 865
 866     forAll(boundaryData0, oldBFaceI)
 867     {
 868         label newFaceI = map.oldFaceMap()[oldBFaceI + map.nOldInternalFaces()];
 869
 870         // Face still exists (is necessary?) and still boundary face
 871         if (newFaceI >= 0 && newFaceI >= mesh.nInternalFaces())
 872         {
 873             newBoundaryData[newFaceI - mesh.nInternalFaces()] =
 874                 boundaryData0[oldBFaceI];
 875         }
 876     }
 877
 878     forAll(boundaryData1, addedBFaceI)
 879     {
 880         label newFaceI = map.addedFaceMap()[addedBFaceI + nInternalFaces1];
 881
 882         if (newFaceI >= 0 && newFaceI >= mesh.nInternalFaces())
 883         {
 884             newBoundaryData[newFaceI - mesh.nInternalFaces()] =
 885                 boundaryData1[addedBFaceI];
 886         }
 887     }
 888
 889     return newBoundaryData;
 890 }
 891
 892
 893 // Remove cells. Add all exposed faces to patch oldInternalPatchI
 894 Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::doRemoveCells
 895 (
 896     const labelList& cellsToRemove,
 897     const label oldInternalPatchI
 898 )
 899 {
 900     // Mesh change engine
 901     directTopoChange meshMod(mesh_);
 902
 903     // Cell removal topo engine. Do NOT synchronize parallel since
 904     // we are doing a local cell removal.
 905     removeCells cellRemover(mesh_, false);
 906
 907     // Get all exposed faces
 908     labelList exposedFaces(cellRemover.getExposedFaces(cellsToRemove));
 909
 910     // Insert the topo changes
 911     cellRemover.setRefinement
 912     (
 913         cellsToRemove,
 914         exposedFaces,
 915         labelList(exposedFaces.size(), oldInternalPatchI),  // patch for exposed
 916                                                             // faces.
 917         meshMod
 918     );
 919
 920     // Change the mesh. No inflation. Note: no parallel comms allowed.
 921     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, false);
 922
 923     // Update fields
 924     mesh_.updateMesh(map);
 925
 926     // Move mesh (since morphing does not do this)
 927     if (map().hasMotionPoints())
 928     {
 929         mesh_.movePoints(map().preMotionPoints());
 930     }
 931
 932     return map;
 933 }
 934
 935
 936 // Delete and add processor patches. Changes mesh and returns per
 937 // neighbour processor the processor patchID.
 938 void Foam::fvMeshDistribute::addProcPatches
 939 (
 940     const labelList& neighbourNewProc,   // processor that neighbour is on
 941     labelList& procPatchID
 942 )
 943 {
 944     // Now use the neighbourFace/Proc to repatch the mesh. These two lists
 945     // contain for all current boundary faces the global patchID (for non-proc
 946     // patch) or the processor.
 947
 948     labelList procPatchSizes(Pstream::nProcs(), 0);
 949
 950     forAll(neighbourNewProc, bFaceI)
 951     {
 952         if (neighbourNewProc[bFaceI] != -1)
 953         {
 954             procPatchSizes[neighbourNewProc[bFaceI]]++;
 955         }
 956     }
 957
 958     // Per neighbour processor the label of the processor patch
 959     procPatchID.setSize(Pstream::nProcs());
 960
 961     forAll(procPatchSizes, procI)
 962     {
 963         if (procPatchSizes[procI] > 0)
 964         {
 965             const word patchName =
 966                 "procBoundary"
 967               + name(Pstream::myProcNo())
 968               + "to"
 969               + name(procI);
 970
 971
 972             procPatchID[procI] = addProcPatch(patchName, procI);
 973             addPatchFields<volScalarField>
 974             (
 975                 processorFvPatchField<scalar>::typeName
 976             );
 977             addPatchFields<volVectorField>
 978             (
 979                 processorFvPatchField<vector>::typeName
 980             );
 981             addPatchFields<volSphericalTensorField>
 982             (
 983                 processorFvPatchField<sphericalTensor>::typeName
 984             );
 985             addPatchFields<volSymmTensorField>
 986             (
 987                 processorFvPatchField<symmTensor>::typeName
 988             );
 989             addPatchFields<volTensorField>
 990             (
 991                 processorFvPatchField<tensor>::typeName
 992             );
 993
 994             addPatchFields<surfaceScalarField>
 995             (
 996                 processorFvPatchField<scalar>::typeName
 997             );
 998             addPatchFields<surfaceVectorField>
 999             (
1000                 processorFvPatchField<vector>::typeName
1001             );
1002             addPatchFields<surfaceSphericalTensorField>
1003             (
1004                 processorFvPatchField<sphericalTensor>::typeName
1005             );
1006             addPatchFields<surfaceSymmTensorField>
1007             (
1008                 processorFvPatchField<symmTensor>::typeName
1009             );
1010             addPatchFields<surfaceTensorField>
1011             (
1012                 processorFvPatchField<tensor>::typeName
1013             );
1014         }
1015         else
1016         {
1017             procPatchID[procI] = -1;
1018         }
1019     }
1020 }
1021
1022
1023 // Get boundary faces to be repatched. Is -1 or new patchID
1024 Foam::labelList Foam::fvMeshDistribute::getProcBoundaryPatch
1025 (
1026     const labelList& neighbourNewProc,  // new processor per boundary face
1027     const labelList& procPatchID        // patchID
1028 )
1029 {
1030     labelList patchIDs(neighbourNewProc);
1031
1032     forAll(neighbourNewProc, bFaceI)
1033     {
1034         if (neighbourNewProc[bFaceI] != -1)
1035         {
1036             label nbrProc = neighbourNewProc[bFaceI];
1037
1038             patchIDs[bFaceI] = procPatchID[nbrProc];
1039         }
1040         else
1041         {
1042             patchIDs[bFaceI] = -1;
1043         }
1044     }
1045     return patchIDs;
1046 }
1047
1048
1049 // Send mesh and coupling data.
1050 void Foam::fvMeshDistribute::sendMesh
1051 (
1052     const label domain,
1053     const fvMesh& mesh,
1054
1055     const wordList& pointZoneNames,
1056     const wordList& faceZoneNames,
1057     const wordList& cellZoneNames,
1058
1059     const labelList& sourceFace,
1060     const labelList& sourceProc,
1061     const labelList& sourceNewProc,
1062     OSstream& toDomain
1063 )
1064 {
1065     if (debug)
1066     {
1067         Pout<< "Sending to domain " << domain << nl
1068             << "    nPoints:" << mesh.nPoints() << nl
1069             << "    nFaces:" << mesh.nFaces() << nl
1070             << "    nCells:" << mesh.nCells() << nl
1071             << "    nPatches:" << mesh.boundaryMesh().size() << nl
1072             << endl;
1073     }
1074
1075     // Assume sparse, possibly overlapping point zones. Get contents
1076     // in merged-zone indices.
1077     CompactListList_dev<label> zonePoints;
1078     {
1079         const pointZoneMesh& pointZones = mesh.pointZones();
1080
1081         labelList rowSizes(pointZoneNames.size(), 0);
1082
1083         forAll(pointZoneNames, nameI)
1084         {
1085             label myZoneID = pointZones.findZoneID(pointZoneNames[nameI]);
1086
1087             if (myZoneID != -1)
1088             {
1089                 rowSizes[nameI] = pointZones[myZoneID].size();
1090             }
1091         }
1092         zonePoints.setSize(rowSizes);
1093
1094         forAll(pointZoneNames, nameI)
1095         {
1096             label myZoneID = pointZones.findZoneID(pointZoneNames[nameI]);
1097
1098             if (myZoneID != -1)
1099             {
1100                 zonePoints[nameI].assign(pointZones[myZoneID]);
1101             }
1102         }
1103     }
1104
1105     // Assume sparse, possibly overlapping face zones
1106     CompactListList_dev<label> zoneFaces;
1107     CompactListList_dev<bool> zoneFaceFlip;
1108     {
1109         const faceZoneMesh& faceZones = mesh.faceZones();
1110
1111         labelList rowSizes(faceZoneNames.size(), 0);
1112
1113         forAll(faceZoneNames, nameI)
1114         {
1115             label myZoneID = faceZones.findZoneID(faceZoneNames[nameI]);
1116
1117             if (myZoneID != -1)
1118             {
1119                 rowSizes[nameI] = faceZones[myZoneID].size();
1120             }
1121         }
1122
1123         zoneFaces.setSize(rowSizes);
1124         zoneFaceFlip.setSize(rowSizes);
1125
1126         forAll(faceZoneNames, nameI)
1127         {
1128             label myZoneID = faceZones.findZoneID(faceZoneNames[nameI]);
1129
1130             if (myZoneID != -1)
1131             {
1132                 zoneFaces[nameI].assign(faceZones[myZoneID]);
1133                 zoneFaceFlip[nameI].assign(faceZones[myZoneID].flipMap());
1134             }
1135         }
1136     }
1137
1138     // Assume sparse, possibly overlapping cell zones
1139     CompactListList_dev<label> zoneCells;
1140     {
1141         const cellZoneMesh& cellZones = mesh.cellZones();
1142
1143         labelList rowSizes(cellZoneNames.size(), 0);
1144
1145         forAll(cellZoneNames, nameI)
1146         {
1147             label myZoneID = cellZones.findZoneID(cellZoneNames[nameI]);
1148
1149             if (myZoneID != -1)
1150             {
1151                 rowSizes[nameI] = cellZones[myZoneID].size();
1152             }
1153         }
1154
1155         zoneCells.setSize(rowSizes);
1156
1157         forAll(cellZoneNames, nameI)
1158         {
1159             label myZoneID = cellZones.findZoneID(cellZoneNames[nameI]);
1160
1161             if (myZoneID != -1)
1162             {
1163                 zoneCells[nameI].assign(cellZones[myZoneID]);
1164             }
1165         }
1166     }
1167     ////- Assume full cell zones
1168     //labelList cellZoneID;
1169     //if (hasCellZones)
1170     //{
1171     //    cellZoneID.setSize(mesh.nCells());
1172     //    cellZoneID = -1;
1173     //
1174     //    const cellZoneMesh& cellZones = mesh.cellZones();
1175     //
1176     //    forAll(cellZones, zoneI)
1177     //    {
1178     //        UIndirectList<label>(cellZoneID, cellZones[zoneI]) = zoneI;
1179     //    }
1180     //}
1181
1182     // Send
1183     toDomain
1184         << mesh.points()
1185         << CompactListList_dev<label, face>(mesh.faces())
1186         << mesh.faceOwner()
1187         << mesh.faceNeighbour()
1188         << mesh.boundaryMesh()
1189
1190         << zonePoints
1191         << zoneFaces
1192         << zoneFaceFlip
1193         << zoneCells
1194
1195         << sourceFace
1196         << sourceProc
1197         << sourceNewProc;
1198
1199
1200     if (debug)
1201     {
1202         Pout<< "Started sending mesh to domain " << domain
1203             << endl;
1204     }
1205 }
1206
1207
1208 // Receive mesh. Opposite of sendMesh
1209 Foam::autoPtr<Foam::fvMesh> Foam::fvMeshDistribute::receiveMesh
1210 (
1211     const label domain,
1212     const wordList& pointZoneNames,
1213     const wordList& faceZoneNames,
1214     const wordList& cellZoneNames,
1215     const Time& runTime,
1216     labelList& domainSourceFace,
1217     labelList& domainSourceProc,
1218     labelList& domainSourceNewProc,
1219     ISstream& fromNbr
1220 )
1221 {
1222     pointField domainPoints(fromNbr);
1223     faceList domainFaces = CompactListList_dev<label, face>(fromNbr)();
1224     labelList domainAllOwner(fromNbr);
1225     labelList domainAllNeighbour(fromNbr);
1226     PtrList<entry> patchEntries(fromNbr);
1227
1228     CompactListList_dev<label> zonePoints(fromNbr);
1229     CompactListList_dev<label> zoneFaces(fromNbr);
1230     CompactListList_dev<bool> zoneFaceFlip(fromNbr);
1231     CompactListList_dev<label> zoneCells(fromNbr);
1232
1233     fromNbr
1234         >> domainSourceFace
1235         >> domainSourceProc
1236         >> domainSourceNewProc;
1237
1238     // Construct fvMesh
1239     autoPtr<fvMesh> domainMeshPtr
1240     (
1241         new fvMesh
1242         (
1243             IOobject
1244             (
1245                 fvMesh::defaultRegion,
1246                 runTime.timeName(),
1247                 runTime,
1248                 IOobject::NO_READ
1249             ),
1250             xferMove(domainPoints),
1251             xferMove(domainFaces),
1252             xferMove(domainAllOwner),
1253             xferMove(domainAllNeighbour),
1254             false                   // no parallel comms
1255         )
1256     );
1257     fvMesh& domainMesh = domainMeshPtr();
1258
1259     List<polyPatch*> patches(patchEntries.size());
1260
1261     forAll(patchEntries, patchI)
1262     {
1263         patches[patchI] = polyPatch::New
1264         (
1265             patchEntries[patchI].keyword(),
1266             patchEntries[patchI].dict(),
1267             patchI,
1268             domainMesh.boundaryMesh()
1269         ).ptr();
1270     }
1271     // Add patches; no parallel comms
1272     domainMesh.addFvPatches(patches, false);
1273
1274     // Construct zones
1275     List<pointZone*> pZonePtrs(pointZoneNames.size());
1276     forAll(pZonePtrs, i)
1277     {
1278         pZonePtrs[i] = new pointZone
1279         (
1280             pointZoneNames[i],
1281             zonePoints[i],
1282             i,
1283             domainMesh.pointZones()
1284         );
1285     }
1286
1287     List<faceZone*> fZonePtrs(faceZoneNames.size());
1288     forAll(fZonePtrs, i)
1289     {
1290         fZonePtrs[i] = new faceZone
1291         (
1292             faceZoneNames[i],
1293             zoneFaces[i],
1294             zoneFaceFlip[i],
1295             i,
1296             domainMesh.faceZones()
1297         );
1298     }
1299
1300     List<cellZone*> cZonePtrs(cellZoneNames.size());
1301     forAll(cZonePtrs, i)
1302     {
1303         cZonePtrs[i] = new cellZone
1304         (
1305             cellZoneNames[i],
1306             zoneCells[i],
1307             i,
1308             domainMesh.cellZones()
1309         );
1310     }
1311     domainMesh.addZones(pZonePtrs, fZonePtrs, cZonePtrs);
1312
1313     return domainMeshPtr;
1314 }
1315
1316
1317 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
1318
1319 // Construct from components
1320 Foam::fvMeshDistribute::fvMeshDistribute(fvMesh& mesh, const scalar mergeTol)
1321 :
1322     mesh_(mesh),
1323     mergeTol_(mergeTol)
1324 {}
1325
1326
1327 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
1328
1329 Foam::labelList Foam::fvMeshDistribute::countCells
1330 (
1331     const labelList& distribution
1332 )
1333 {
1334     labelList nCells(Pstream::nProcs(), 0);
1335     forAll(distribution, cellI)
1336     {
1337         label newProc = distribution[cellI];
1338
1339         if (newProc < 0 || newProc >= Pstream::nProcs())
1340         {
1341             FatalErrorIn("fvMeshDistribute::distribute(const labelList&)")
1342                 << "Distribution should be in range 0.." << Pstream::nProcs()-1
1343                 << endl
1344                 << "At index " << cellI << " distribution:" << newProc
1345                 << abort(FatalError);
1346         }
1347         nCells[newProc]++;
1348     }
1349     return nCells;
1350 }
1351
1352
1353 Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
1354 (
1355     const labelList& distribution
1356 )
1357 {
1358     // Some checks on distribution
1359     if (distribution.size() != mesh_.nCells())
1360     {
1361         FatalErrorIn("fvMeshDistribute::distribute(const labelList&)")
1362             << "Size of distribution:"
1363             << distribution.size() << " mesh nCells:" << mesh_.nCells()
1364             << abort(FatalError);
1365     }
1366
1367
1368     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1369
1370     // Check all processors have same non-proc patches in same order.
1371     if (patches.checkParallelSync(true))
1372     {
1373         FatalErrorIn("fvMeshDistribute::distribute(const labelList&)")
1374             << "This application requires all non-processor patches"
1375             << " to be present in the same order on all patches" << nl
1376             << "followed by the processor patches (which are unique)."
1377             << nl
1378             << "Local patches:" << mesh_.boundaryMesh().names()
1379             << abort(FatalError);
1380     }
1381
1382     // Save some data for mapping later on
1383     const label nOldPoints(mesh_.nPoints());
1384     const label nOldFaces(mesh_.nFaces());
1385     const label nOldCells(mesh_.nCells());
1386     labelList oldPatchStarts(patches.size());
1387     labelList oldPatchNMeshPoints(patches.size());
1388     forAll(patches, patchI)
1389     {
1390         oldPatchStarts[patchI] = patches[patchI].start();
1391         oldPatchNMeshPoints[patchI] = patches[patchI].nPoints();
1392     }
1393
1394
1395
1396     // Short circuit trivial case.
1397     if (!Pstream::parRun())
1398     {
1399         // Collect all maps and return
1400         return autoPtr<mapDistributePolyMesh>
1401         (
1402             new mapDistributePolyMesh
1403             (
1404                 mesh_,
1405
1406                 nOldPoints,
1407                 nOldFaces,
1408                 nOldCells,
1409                 oldPatchStarts.xfer(),
1410                 oldPatchNMeshPoints.xfer(),
1411
1412                 labelListList(1, identity(mesh_.nPoints())).xfer(),//subPointMap
1413                 labelListList(1, identity(mesh_.nFaces())).xfer(), //subFaceMap
1414                 labelListList(1, identity(mesh_.nCells())).xfer(), //subCellMap
1415                 labelListList(1, identity(patches.size())).xfer(), //subPatchMap
1416
1417                 labelListList(1, identity(mesh_.nPoints())).xfer(),//pointMap
1418                 labelListList(1, identity(mesh_.nFaces())).xfer(), //faceMap
1419                 labelListList(1, identity(mesh_.nCells())).xfer(), //cellMap
1420                 labelListList(1, identity(patches.size())).xfer()  //patchMap
1421             )
1422         );
1423     }
1424
1425
1426     // Collect any zone names
1427     const wordList pointZoneNames(mergeWordList(mesh_.pointZones().names()));
1428     const wordList faceZoneNames(mergeWordList(mesh_.faceZones().names()));
1429     const wordList cellZoneNames(mergeWordList(mesh_.cellZones().names()));
1430
1431
1432
1433     // Local environment of all boundary faces
1434     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1435
1436     // A face is uniquely defined by
1437     //  - proc
1438     //  - local face no
1439     //
1440     // To glue the parts of meshes which can get sent from anywhere we
1441     // need to know on boundary faces what the above tuple on both sides is.
1442     // So we need to maintain:
1443     //  - original face
1444     //  - original processor id (= trivial)
1445     // For coupled boundaries (where the faces are 'duplicate') we take the
1446     // lowest numbered processor as the data to store.
1447     //
1448     // Additionally to create the procboundaries we need to know where the owner
1449     // cell on the other side now is: newNeighbourProc.
1450     //
1451
1452     // physical boundary:
1453     //     sourceProc = -1
1454     //     sourceNewProc = -1
1455     //     sourceFace = patchID
1456     // coupled boundary:
1457     //     sourceProc = proc
1458     //     sourceNewProc = distribution[cell on proc]
1459     //     sourceFace = face
1460     labelList sourceFace;
1461     labelList sourceProc;
1462     labelList sourceNewProc;
1463     getNeighbourData(distribution, sourceFace, sourceProc, sourceNewProc);
1464
1465
1466     // Remove meshPhi. Since this would otherwise dissappear anyway
1467     // during topo changes and we have to guarantee that all the fields
1468     // can be sent.
1469     mesh_.clearOut();
1470     mesh_.resetMotion();
1471
1472     // Get data to send. Make sure is synchronised
1473     const wordList volScalars(mesh_.names(volScalarField::typeName));
1474     checkEqualWordList("volScalarFields", volScalars);
1475     const wordList volVectors(mesh_.names(volVectorField::typeName));
1476     checkEqualWordList("volVectorFields", volVectors);
1477     const wordList volSphereTensors
1478     (
1479         mesh_.names(volSphericalTensorField::typeName)
1480     );
1481     checkEqualWordList("volSphericalTensorFields", volSphereTensors);
1482     const wordList volSymmTensors(mesh_.names(volSymmTensorField::typeName));
1483     checkEqualWordList("volSymmTensorFields", volSymmTensors);
1484     const wordList volTensors(mesh_.names(volTensorField::typeName));
1485     checkEqualWordList("volTensorField", volTensors);
1486
1487     const wordList surfScalars(mesh_.names(surfaceScalarField::typeName));
1488     checkEqualWordList("surfaceScalarFields", surfScalars);
1489     const wordList surfVectors(mesh_.names(surfaceVectorField::typeName));
1490     checkEqualWordList("surfaceVectorFields", surfVectors);
1491     const wordList surfSphereTensors
1492     (
1493         mesh_.names(surfaceSphericalTensorField::typeName)
1494     );
1495     checkEqualWordList("surfaceSphericalTensorFields", surfSphereTensors);
1496     const wordList surfSymmTensors
1497     (
1498         mesh_.names(surfaceSymmTensorField::typeName)
1499     );
1500     checkEqualWordList("surfaceSymmTensorFields", surfSymmTensors);
1501     const wordList surfTensors(mesh_.names(surfaceTensorField::typeName));
1502     checkEqualWordList("surfaceTensorFields", surfTensors);
1503
1504
1505
1506
1507     // Find patch to temporarily put exposed and processor faces into.
1508     label oldInternalPatchI = findNonEmptyPatch();
1509
1510
1511
1512     // Delete processor patches, starting from the back. Move all faces into
1513     // oldInternalPatchI.
1514     labelList repatchFaceMap;
1515     {
1516         autoPtr<mapPolyMesh> repatchMap = deleteProcPatches(oldInternalPatchI);
1517
1518         // Store face map (only face ordering that changed)
1519         repatchFaceMap = repatchMap().faceMap();
1520
1521         // Reorder all boundary face data (sourceProc, sourceFace etc.)
1522         labelList bFaceMap
1523         (
1524             SubList<label>
1525             (
1526                 repatchMap().reverseFaceMap(),
1527                 mesh_.nFaces() - mesh_.nInternalFaces(),
1528                 mesh_.nInternalFaces()
1529             )
1530           - mesh_.nInternalFaces()
1531         );
1532
1533         inplaceReorder(bFaceMap, sourceFace);
1534         inplaceReorder(bFaceMap, sourceProc);
1535         inplaceReorder(bFaceMap, sourceNewProc);
1536     }
1537
1538
1539
1540     // Print a bit.
1541     if (debug)
1542     {
1543         Pout<< nl << "MESH WITH PROC PATCHES DELETED:" << endl;
1544         printMeshInfo(mesh_);
1545         printFieldInfo<volScalarField>(mesh_);
1546         printFieldInfo<volVectorField>(mesh_);
1547         printFieldInfo<volSphericalTensorField>(mesh_);
1548         printFieldInfo<volSymmTensorField>(mesh_);
1549         printFieldInfo<volTensorField>(mesh_);
1550         printFieldInfo<surfaceScalarField>(mesh_);
1551         printFieldInfo<surfaceVectorField>(mesh_);
1552         printFieldInfo<surfaceSphericalTensorField>(mesh_);
1553         printFieldInfo<surfaceSymmTensorField>(mesh_);
1554         printFieldInfo<surfaceTensorField>(mesh_);
1555         Pout<< nl << endl;
1556     }
1557
1558
1559
1560     // Maps from subsetted mesh (that is sent) back to original maps
1561     labelListList subCellMap(Pstream::nProcs());
1562     labelListList subFaceMap(Pstream::nProcs());
1563     labelListList subPointMap(Pstream::nProcs());
1564     labelListList subPatchMap(Pstream::nProcs());
1565     // Maps from subsetted mesh to reconstructed mesh
1566     labelListList constructCellMap(Pstream::nProcs());
1567     labelListList constructFaceMap(Pstream::nProcs());
1568     labelListList constructPointMap(Pstream::nProcs());
1569     labelListList constructPatchMap(Pstream::nProcs());
1570
1571
1572
1573
1574     // Find out schedule
1575     // ~~~~~~~~~~~~~~~~~
1576
1577     labelListList nSendCells(Pstream::nProcs());
1578     nSendCells[Pstream::myProcNo()] = countCells(distribution);
1579     Pstream::gatherList(nSendCells);
1580     Pstream::scatterList(nSendCells);
1581
1582
1583     // Allocate buffers
1584     PtrList<OStringStream> sendStr(Pstream::nProcs());
1585     forAll(sendStr, i)
1586     {
1587         sendStr.set(i, new OStringStream(IOstream::BINARY));
1588     }
1589     //PstreamBuffers pBufs(Pstream::nonBlocking);
1590
1591
1592     // What to send to neighbouring domains
1593     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1594
1595     forAll(nSendCells[Pstream::myProcNo()], recvProc)
1596     {
1597         if
1598         (
1599             recvProc != Pstream::myProcNo()
1600          && nSendCells[Pstream::myProcNo()][recvProc] > 0
1601         )
1602         {
1603             // Send to recvProc
1604
1605             if (debug)
1606             {
1607                 Pout<< nl
1608                     << "SUBSETTING FOR DOMAIN " << recvProc
1609                     << " cells to send:"
1610                     << nSendCells[Pstream::myProcNo()][recvProc]
1611                     << nl << endl;
1612             }
1613
1614             // Pstream for sending mesh and fields
1615             //OPstream str(Pstream::blocking, recvProc);
1616             //UOPstream str(recvProc, pBufs);
1617
1618             // Mesh subsetting engine
1619             fvMeshSubset subsetter
1620             (
1621                 IOobject
1622                 (
1623                     "set",
1624                     mesh_.time().timeName(),
1625                     mesh_,
1626                     IOobject::NO_READ,
1627                     IOobject::NO_WRITE
1628                 ),
1629                 mesh_
1630             );
1631
1632             // Subset the cells of the current domain.
1633             subsetter.setLargeCellSubset
1634             (
1635                 distribution,
1636                 recvProc,
1637                 oldInternalPatchI,  // oldInternalFaces patch
1638                 false               // no parallel sync
1639             );
1640
1641             subCellMap[recvProc] = subsetter.cellMap();
1642             subFaceMap[recvProc] = renumber
1643             (
1644                 repatchFaceMap,
1645                 subsetter.faceMap()
1646             );
1647             subPointMap[recvProc] = subsetter.pointMap();
1648             subPatchMap[recvProc] = subsetter.patchMap();
1649
1650
1651             // Subset the boundary fields (owner/neighbour/processor)
1652             labelList procSourceFace;
1653             labelList procSourceProc;
1654             labelList procSourceNewProc;
1655
1656             subsetBoundaryData
1657             (
1658                 subsetter.subMesh(),
1659                 subsetter.faceMap(),        // from subMesh to mesh
1660                 subsetter.cellMap(),        //      ,,      ,,
1661
1662                 distribution,               // old mesh distribution
1663                 mesh_.faceOwner(),          // old owner
1664                 mesh_.faceNeighbour(),
1665                 mesh_.nInternalFaces(),
1666
1667                 sourceFace,
1668                 sourceProc,
1669                 sourceNewProc,
1670
1671                 procSourceFace,
1672                 procSourceProc,
1673                 procSourceNewProc
1674             );
1675
1676
1677
1678             // Send to neighbour
1679             sendMesh
1680             (
1681                 recvProc,
1682                 subsetter.subMesh(),
1683
1684                 pointZoneNames,
1685                 faceZoneNames,
1686                 cellZoneNames,
1687
1688                 procSourceFace,
1689                 procSourceProc,
1690                 procSourceNewProc,
1691                 sendStr[recvProc]
1692             );
1693             sendFields<volScalarField>
1694             (
1695                 recvProc,
1696                 volScalars,
1697                 subsetter,
1698                 sendStr[recvProc]
1699             );
1700             sendFields<volVectorField>
1701             (
1702                 recvProc,
1703                 volVectors,
1704                 subsetter,
1705                 sendStr[recvProc]
1706             );
1707             sendFields<volSphericalTensorField>
1708             (
1709                 recvProc,
1710                 volSphereTensors,
1711                 subsetter,
1712                 sendStr[recvProc]
1713             );
1714             sendFields<volSymmTensorField>
1715             (
1716                 recvProc,
1717                 volSymmTensors,
1718                 subsetter,
1719                 sendStr[recvProc]
1720             );
1721             sendFields<volTensorField>
1722             (
1723                 recvProc,
1724                 volTensors,
1725                 subsetter,
1726                 sendStr[recvProc]
1727             );
1728
1729             sendFields<surfaceScalarField>
1730             (
1731                 recvProc,
1732                 surfScalars,
1733                 subsetter,
1734                 sendStr[recvProc]
1735             );
1736             sendFields<surfaceVectorField>
1737             (
1738                 recvProc,
1739                 surfVectors,
1740                 subsetter,
1741                 sendStr[recvProc]
1742             );
1743             sendFields<surfaceSphericalTensorField>
1744             (
1745                 recvProc,
1746                 surfSphereTensors,
1747                 subsetter,
1748                 sendStr[recvProc]
1749             );
1750             sendFields<surfaceSymmTensorField>
1751             (
1752                 recvProc,
1753                 surfSymmTensors,
1754                 subsetter,
1755                 sendStr[recvProc]
1756             );
1757             sendFields<surfaceTensorField>
1758             (
1759                 recvProc,
1760                 surfTensors,
1761                 subsetter,
1762                 sendStr[recvProc]
1763             );
1764         }
1765     }
1766
1767
1768     // Start sending&receiving from buffers
1769     //pBufs.finishedSends();
1770
1771     // get the data.
1772     PtrList<IStringStream> recvStr(Pstream::nProcs());
1773     {
1774         List<List<char> > sendBufs(sendStr.size());
1775         forAll(sendStr, procI)
1776         {
1777             string contents = sendStr[procI].str();
1778             const char* ptr = contents.data();
1779
1780             sendBufs[procI].setSize(contents.size());
1781             forAll(sendBufs[procI], i)
1782             {
1783                 sendBufs[procI][i] = *ptr++;
1784             }
1785             // Clear OStringStream
1786             sendStr.set(procI, NULL);
1787         }
1788
1789         // Transfer sendBufs into recvBufs
1790         List<List<char> > recvBufs(Pstream::nProcs());
1791         labelListList sizes(Pstream::nProcs());
1792         exchange<List<char>, char>(sendBufs, recvBufs, sizes);
1793
1794         forAll(recvStr, procI)
1795         {
1796             string contents(recvBufs[procI].begin(), recvBufs[procI].size());
1797             recvStr.set
1798             (
1799                 procI,
1800                 new IStringStream(contents, IOstream::BINARY)
1801             );
1802         }
1803     }
1804
1805
1806     // Subset the part that stays
1807     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
1808
1809     {
1810         // Save old mesh maps before changing mesh
1811         const labelList oldFaceOwner(mesh_.faceOwner());
1812         const labelList oldFaceNeighbour(mesh_.faceNeighbour());
1813         const label oldInternalFaces = mesh_.nInternalFaces();
1814
1815         // Remove cells.
1816         autoPtr<mapPolyMesh> subMap
1817         (
1818             doRemoveCells
1819             (
1820                 select(false, distribution, Pstream::myProcNo()),
1821                 oldInternalPatchI
1822             )
1823         );
1824
1825         // Addressing from subsetted mesh
1826         subCellMap[Pstream::myProcNo()] = subMap().cellMap();
1827         subFaceMap[Pstream::myProcNo()] = renumber
1828         (
1829             repatchFaceMap,
1830             subMap().faceMap()
1831         );
1832         subPointMap[Pstream::myProcNo()] = subMap().pointMap();
1833         subPatchMap[Pstream::myProcNo()] = identity(patches.size());
1834
1835         // Initialize all addressing into current mesh
1836         constructCellMap[Pstream::myProcNo()] = identity(mesh_.nCells());
1837         constructFaceMap[Pstream::myProcNo()] = identity(mesh_.nFaces());
1838         constructPointMap[Pstream::myProcNo()] = identity(mesh_.nPoints());
1839         constructPatchMap[Pstream::myProcNo()] = identity(patches.size());
1840
1841         // Subset the mesh data: neighbourCell/neighbourProc
1842         // fields
1843         labelList domainSourceFace;
1844         labelList domainSourceProc;
1845         labelList domainSourceNewProc;
1846
1847         subsetBoundaryData
1848         (
1849             mesh_,                          // new mesh
1850             subMap().faceMap(),             // from new to original mesh
1851             subMap().cellMap(),
1852
1853             distribution,                   // distribution before subsetting
1854             oldFaceOwner,                   // owner before subsetting
1855             oldFaceNeighbour,               // neighbour        ,,
1856             oldInternalFaces,               // nInternalFaces   ,,
1857
1858             sourceFace,
1859             sourceProc,
1860             sourceNewProc,
1861
1862             domainSourceFace,
1863             domainSourceProc,
1864             domainSourceNewProc
1865         );
1866
1867         sourceFace.transfer(domainSourceFace);
1868         sourceProc.transfer(domainSourceProc);
1869         sourceNewProc.transfer(domainSourceNewProc);
1870     }
1871
1872
1873     // Print a bit.
1874     if (debug)
1875     {
1876         Pout<< nl << "STARTING MESH:" << endl;
1877         printMeshInfo(mesh_);
1878         printFieldInfo<volScalarField>(mesh_);
1879         printFieldInfo<volVectorField>(mesh_);
1880         printFieldInfo<volSphericalTensorField>(mesh_);
1881         printFieldInfo<volSymmTensorField>(mesh_);
1882         printFieldInfo<volTensorField>(mesh_);
1883         printFieldInfo<surfaceScalarField>(mesh_);
1884         printFieldInfo<surfaceVectorField>(mesh_);
1885         printFieldInfo<surfaceSphericalTensorField>(mesh_);
1886         printFieldInfo<surfaceSymmTensorField>(mesh_);
1887         printFieldInfo<surfaceTensorField>(mesh_);
1888         Pout<< nl << endl;
1889     }
1890
1891
1892
1893     // Receive and add what was sent
1894     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1895
1896     forAll(nSendCells, sendProc)
1897     {
1898         // Did processor sendProc send anything to me?
1899         if
1900         (
1901             sendProc != Pstream::myProcNo()
1902          && nSendCells[sendProc][Pstream::myProcNo()] > 0
1903         )
1904         {
1905             if (debug)
1906             {
1907                 Pout<< nl
1908                     << "RECEIVING FROM DOMAIN " << sendProc
1909                     << " cells to receive:"
1910                     << nSendCells[sendProc][Pstream::myProcNo()]
1911                     << nl << endl;
1912             }
1913
1914
1915             // Pstream for receiving mesh and fields
1916             //UIPstream str(sendProc, pBufs);
1917
1918
1919             // Receive from sendProc
1920             labelList domainSourceFace;
1921             labelList domainSourceProc;
1922             labelList domainSourceNewProc;
1923
1924             autoPtr<fvMesh> domainMeshPtr;
1925             PtrList<volScalarField> vsf;
1926             PtrList<volVectorField> vvf;
1927             PtrList<volSphericalTensorField> vsptf;
1928             PtrList<volSymmTensorField> vsytf;
1929             PtrList<volTensorField> vtf;
1930             PtrList<surfaceScalarField> ssf;
1931             PtrList<surfaceVectorField> svf;
1932             PtrList<surfaceSphericalTensorField> ssptf;
1933             PtrList<surfaceSymmTensorField> ssytf;
1934             PtrList<surfaceTensorField> stf;
1935
1936             // Opposite of sendMesh
1937             {
1938                 domainMeshPtr = receiveMesh
1939                 (
1940                     sendProc,
1941                     pointZoneNames,
1942                     faceZoneNames,
1943                     cellZoneNames,
1944
1945                     const_cast<Time&>(mesh_.time()),
1946                     domainSourceFace,
1947                     domainSourceProc,
1948                     domainSourceNewProc,
1949                     recvStr[sendProc]
1950                 );
1951                 fvMesh& domainMesh = domainMeshPtr();
1952
1953                 // Receive fields. Read as single dictionary because
1954                 // of problems reading consecutive fields from single stream.
1955                 dictionary fieldDicts(recvStr[sendProc]);
1956
1957                 receiveFields<volScalarField>
1958                 (
1959                     sendProc,
1960                     volScalars,
1961                     domainMesh,
1962                     vsf,
1963                     fieldDicts.subDict(volScalarField::typeName)
1964                 );
1965                 receiveFields<volVectorField>
1966                 (
1967                     sendProc,
1968                     volVectors,
1969                     domainMesh,
1970                     vvf,
1971                     fieldDicts.subDict(volVectorField::typeName)
1972                 );
1973                 receiveFields<volSphericalTensorField>
1974                 (
1975                     sendProc,
1976                     volSphereTensors,
1977                     domainMesh,
1978                     vsptf,
1979                     fieldDicts.subDict(volSphericalTensorField::typeName)
1980                 );
1981                 receiveFields<volSymmTensorField>
1982                 (
1983                     sendProc,
1984                     volSymmTensors,
1985                     domainMesh,
1986                     vsytf,
1987                     fieldDicts.subDict(volSymmTensorField::typeName)
1988                 );
1989                 receiveFields<volTensorField>
1990                 (
1991                     sendProc,
1992                     volTensors,
1993                     domainMesh,
1994                     vtf,
1995                     fieldDicts.subDict(volTensorField::typeName)
1996                 );
1997
1998                 receiveFields<surfaceScalarField>
1999                 (
2000                     sendProc,
2001                     surfScalars,
2002                     domainMesh,
2003                     ssf,
2004                     fieldDicts.subDict(surfaceScalarField::typeName)
2005                 );
2006                 receiveFields<surfaceVectorField>
2007                 (
2008                     sendProc,
2009                     surfVectors,
2010                     domainMesh,
2011                     svf,
2012                     fieldDicts.subDict(surfaceVectorField::typeName)
2013                 );
2014                 receiveFields<surfaceSphericalTensorField>
2015                 (
2016                     sendProc,
2017                     surfSphereTensors,
2018                     domainMesh,
2019                     ssptf,
2020                     fieldDicts.subDict(surfaceSphericalTensorField::typeName)
2021                 );
2022                 receiveFields<surfaceSymmTensorField>
2023                 (
2024                     sendProc,
2025                     surfSymmTensors,
2026                     domainMesh,
2027                     ssytf,
2028                     fieldDicts.subDict(surfaceSymmTensorField::typeName)
2029                 );
2030                 receiveFields<surfaceTensorField>
2031                 (
2032                     sendProc,
2033                     surfTensors,
2034                     domainMesh,
2035                     stf,
2036                     fieldDicts.subDict(surfaceTensorField::typeName)
2037                 );
2038             }
2039             const fvMesh& domainMesh = domainMeshPtr();
2040
2041
2042             constructCellMap[sendProc] = identity(domainMesh.nCells());
2043             constructFaceMap[sendProc] = identity(domainMesh.nFaces());
2044             constructPointMap[sendProc] = identity(domainMesh.nPoints());
2045             constructPatchMap[sendProc] =
2046                 identity(domainMesh.boundaryMesh().size());
2047
2048
2049             // Print a bit.
2050             if (debug)
2051             {
2052                 Pout<< nl << "RECEIVED MESH FROM:" << sendProc << endl;
2053                 printMeshInfo(domainMesh);
2054                 printFieldInfo<volScalarField>(domainMesh);
2055                 printFieldInfo<volVectorField>(domainMesh);
2056                 printFieldInfo<volSphericalTensorField>(domainMesh);
2057                 printFieldInfo<volSymmTensorField>(domainMesh);
2058                 printFieldInfo<volTensorField>(domainMesh);
2059                 printFieldInfo<surfaceScalarField>(domainMesh);
2060                 printFieldInfo<surfaceVectorField>(domainMesh);
2061                 printFieldInfo<surfaceSphericalTensorField>(domainMesh);
2062                 printFieldInfo<surfaceSymmTensorField>(domainMesh);
2063                 printFieldInfo<surfaceTensorField>(domainMesh);
2064             }
2065
2066
2067             // Now this mesh we received (from sendProc) needs to be merged
2068             // with the current mesh. On the current mesh we have for all
2069             // boundaryfaces the original face and processor. See if we can
2070             // match these up to the received domainSourceFace and
2071             // domainSourceProc.
2072             labelList masterCoupledFaces;
2073             labelList slaveCoupledFaces;
2074             findCouples
2075             (
2076                 mesh_,
2077
2078                 sourceFace,
2079                 sourceProc,
2080
2081                 sendProc,
2082                 domainMesh,
2083                 domainSourceFace,
2084                 domainSourceProc,
2085
2086                 masterCoupledFaces,
2087                 slaveCoupledFaces
2088             );
2089
2090             // Generate additional coupling info (points, edges) from
2091             // faces-that-match
2092             faceCoupleInfo couples
2093             (
2094                 mesh_,
2095                 masterCoupledFaces,
2096                 domainMesh,
2097                 slaveCoupledFaces,
2098                 mergeTol_,              // merge tolerance
2099                 true,                   // faces align
2100                 true,                   // couples are ordered already
2101                 false
2102             );
2103
2104
2105             // Add domainMesh to mesh
2106             // ~~~~~~~~~~~~~~~~~~~~~~
2107
2108             autoPtr<mapAddedPolyMesh> map = fvMeshAdder::add
2109             (
2110                 mesh_,
2111                 domainMesh,
2112                 couples,
2113                 false           // no parallel comms
2114             );
2115
2116             // Update mesh data: sourceFace,sourceProc for added
2117             // mesh.
2118
2119             sourceFace =
2120                 mapBoundaryData
2121                 (
2122                     mesh_,
2123                     map(),
2124                     sourceFace,
2125                     domainMesh.nInternalFaces(),
2126                     domainSourceFace
2127                 );
2128             sourceProc =
2129                 mapBoundaryData
2130                 (
2131                     mesh_,
2132                     map(),
2133                     sourceProc,
2134                     domainMesh.nInternalFaces(),
2135                     domainSourceProc
2136                 );
2137             sourceNewProc =
2138                 mapBoundaryData
2139                 (
2140                     mesh_,
2141                     map(),
2142                     sourceNewProc,
2143                     domainMesh.nInternalFaces(),
2144                     domainSourceNewProc
2145                 );
2146
2147             // Update all addressing so xxProcAddressing points to correct item
2148             // in masterMesh.
2149             const labelList& oldCellMap = map().oldCellMap();
2150             const labelList& oldFaceMap = map().oldFaceMap();
2151             const labelList& oldPointMap = map().oldPointMap();
2152             const labelList& oldPatchMap = map().oldPatchMap();
2153
2154             forAll(constructPatchMap, procI)
2155             {
2156                 if (procI != sendProc && constructPatchMap[procI].size())
2157                 {
2158                     // Processor already in mesh (either myProcNo or received)
2159                     inplaceRenumber(oldCellMap, constructCellMap[procI]);
2160                     inplaceRenumber(oldFaceMap, constructFaceMap[procI]);
2161                     inplaceRenumber(oldPointMap, constructPointMap[procI]);
2162                     inplaceRenumber(oldPatchMap, constructPatchMap[procI]);
2163                 }
2164             }
2165
2166             // Added processor
2167             inplaceRenumber(map().addedCellMap(), constructCellMap[sendProc]);
2168             inplaceRenumber(map().addedFaceMap(), constructFaceMap[sendProc]);
2169             inplaceRenumber(map().addedPointMap(), constructPointMap[sendProc]);
2170             inplaceRenumber(map().addedPatchMap(), constructPatchMap[sendProc]);
2171
2172             if (debug)
2173             {
2174                 Pout<< nl << "MERGED MESH FROM:" << sendProc << endl;
2175                 printMeshInfo(mesh_);
2176                 printFieldInfo<volScalarField>(mesh_);
2177                 printFieldInfo<volVectorField>(mesh_);
2178                 printFieldInfo<volSphericalTensorField>(mesh_);
2179                 printFieldInfo<volSymmTensorField>(mesh_);
2180                 printFieldInfo<volTensorField>(mesh_);
2181                 printFieldInfo<surfaceScalarField>(mesh_);
2182                 printFieldInfo<surfaceVectorField>(mesh_);
2183                 printFieldInfo<surfaceSphericalTensorField>(mesh_);
2184                 printFieldInfo<surfaceSymmTensorField>(mesh_);
2185                 printFieldInfo<surfaceTensorField>(mesh_);
2186                 Pout<< nl << endl;
2187             }
2188         }
2189     }
2190
2191
2192     // Print a bit.
2193     if (debug)
2194     {
2195         Pout<< nl << "REDISTRIBUTED MESH:" << endl;
2196         printMeshInfo(mesh_);
2197         printFieldInfo<volScalarField>(mesh_);
2198         printFieldInfo<volVectorField>(mesh_);
2199         printFieldInfo<volSphericalTensorField>(mesh_);
2200         printFieldInfo<volSymmTensorField>(mesh_);
2201         printFieldInfo<volTensorField>(mesh_);
2202         printFieldInfo<surfaceScalarField>(mesh_);
2203         printFieldInfo<surfaceVectorField>(mesh_);
2204         printFieldInfo<surfaceSphericalTensorField>(mesh_);
2205         printFieldInfo<surfaceSymmTensorField>(mesh_);
2206         printFieldInfo<surfaceTensorField>(mesh_);
2207         Pout<< nl << endl;
2208     }
2209
2210
2211
2212     // Add processorPatches
2213     // ~~~~~~~~~~~~~~~~~~~~
2214
2215     // Per neighbour processor the patchID to it (or -1).
2216     labelList procPatchID;
2217
2218     // Add processor patches.
2219     addProcPatches(sourceNewProc, procPatchID);
2220
2221     // Put faces into correct patch. Note that we now have proper
2222     // processorPolyPatches again so repatching will take care of coupled face
2223     // ordering.
2224
2225     // Get boundary faces to be repatched. Is -1 or new patchID
2226     labelList newPatchID
2227     (
2228         getProcBoundaryPatch
2229         (
2230             sourceNewProc,
2231             procPatchID
2232         )
2233     );
2234
2235     // Change patches. Since this might change ordering of coupled faces
2236     // we also need to adapt our constructMaps.
2237     repatch(newPatchID, constructFaceMap);
2238
2239     // See if any geometrically shared points need to be merged. Note: does
2240     // parallel comms.
2241     mergeSharedPoints(constructPointMap);
2242
2243     // Bit of hack: processorFvPatchField does not get reset since created
2244     // from nothing so explicitly reset.
2245     initPatchFields<volScalarField, processorFvPatchField<scalar> >
2246     (
2247         pTraits<scalar>::zero
2248     );
2249     initPatchFields<volVectorField, processorFvPatchField<vector> >
2250     (
2251         pTraits<vector>::zero
2252     );
2253     initPatchFields
2254     <
2255         volSphericalTensorField,
2256         processorFvPatchField<sphericalTensor>
2257     >
2258     (
2259         pTraits<sphericalTensor>::zero
2260     );
2261     initPatchFields<volSymmTensorField, processorFvPatchField<symmTensor> >
2262     (
2263         pTraits<symmTensor>::zero
2264     );
2265     initPatchFields<volTensorField, processorFvPatchField<tensor> >
2266     (
2267         pTraits<tensor>::zero
2268     );
2269     initPatchFields<surfaceScalarField, processorFvsPatchField<scalar> >
2270     (
2271         pTraits<scalar>::zero
2272     );
2273     initPatchFields<surfaceVectorField, processorFvsPatchField<vector> >
2274     (
2275         pTraits<vector>::zero
2276     );
2277     initPatchFields
2278     <
2279         surfaceSphericalTensorField,
2280         processorFvsPatchField<sphericalTensor>
2281     >
2282     (
2283         pTraits<sphericalTensor>::zero
2284     );
2285     initPatchFields
2286     <
2287         surfaceSymmTensorField,
2288         processorFvsPatchField<symmTensor>
2289     >
2290     (
2291         pTraits<symmTensor>::zero
2292     );
2293     initPatchFields<surfaceTensorField, processorFvsPatchField<tensor> >
2294     (
2295         pTraits<tensor>::zero
2296     );
2297
2298
2299     mesh_.setInstance(mesh_.time().timeName());
2300
2301
2302     // Print a bit
2303     if (debug)
2304     {
2305         Pout<< nl << "FINAL MESH:" << endl;
2306         printMeshInfo(mesh_);
2307         printFieldInfo<volScalarField>(mesh_);
2308         printFieldInfo<volVectorField>(mesh_);
2309         printFieldInfo<volSphericalTensorField>(mesh_);
2310         printFieldInfo<volSymmTensorField>(mesh_);
2311         printFieldInfo<volTensorField>(mesh_);
2312         printFieldInfo<surfaceScalarField>(mesh_);
2313         printFieldInfo<surfaceVectorField>(mesh_);
2314         printFieldInfo<surfaceSphericalTensorField>(mesh_);
2315         printFieldInfo<surfaceSymmTensorField>(mesh_);
2316         printFieldInfo<surfaceTensorField>(mesh_);
2317         Pout<< nl << endl;
2318     }
2319
2320     // Collect all maps and return
2321     return autoPtr<mapDistributePolyMesh>
2322     (
2323         new mapDistributePolyMesh
2324         (
2325             mesh_,
2326
2327             nOldPoints,
2328             nOldFaces,
2329             nOldCells,
2330             oldPatchStarts,
2331             oldPatchNMeshPoints,
2332
2333             subPointMap,
2334             subFaceMap,
2335             subCellMap,
2336             subPatchMap,
2337
2338             constructPointMap,
2339             constructFaceMap,
2340             constructCellMap,
2341             constructPatchMap,
2342             true                // reuse storage
2343         )
2344     );
2345 }
2346
2347
2348 // ************************************************************************* //