src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C

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