src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C

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