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