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-2008 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() > 0)
 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() > 0)
 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() > 0)
 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.size() == 0)
 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 != 0)
 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 a shuffle?"
 590                 << 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     bool merged = pointToMaster.size() > 0;
 626
 627     if (!returnReduce(merged, orOp<bool>()))
 628     {
 629         return autoPtr<mapPolyMesh>(NULL);
 630     }
 631
 632     polyTopoChange meshMod(mesh_);
 633
 634     fvMeshAdder::mergePoints(mesh_, pointToMaster, meshMod);
 635
 636     // Change the mesh (no inflation). Note: parallel comms allowed.
 637     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
 638
 639     // Update fields. No inflation, parallel sync.
 640     mesh_.updateMesh(map);
 641
 642     // Move mesh (since morphing does not do this)
 643     if (map().hasMotionPoints())
 644     {
 645         mesh_.movePoints(map().preMotionPoints());
 646     }
 647
 648     // Adapt constructMaps for merged points.
 649     // 1.4.1: use reversePointMap < -1 feature.
 650     forAll(constructPointMap, procI)
 651     {
 652         labelList& constructMap = constructPointMap[procI];
 653
 654         forAll(constructMap, i)
 655         {
 656             label oldPointI = constructMap[i];
 657
 658             // See if merged into other point
 659             Map<label>::const_iterator iter = pointToMaster.find(oldPointI);
 660
 661             if (iter != pointToMaster.end())
 662             {
 663                 oldPointI = iter();
 664             }
 665
 666             constructMap[i] = map().reversePointMap()[oldPointI];
 667         }
 668     }
 669     return map;
 670 }
 671
 672
 673 // Construct the local environment of all boundary faces.
 674 void Foam::fvMeshDistribute::getNeighbourData
 675 (
 676     const labelList& distribution,
 677     labelList& sourceFace,
 678     labelList& sourceProc,
 679     labelList& sourceNewProc
 680 ) const
 681 {
 682     sourceFace.setSize(mesh_.nFaces() - mesh_.nInternalFaces());
 683     sourceProc.setSize(sourceFace.size());
 684     sourceNewProc.setSize(sourceFace.size());
 685
 686     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 687
 688     // Send meshFace labels of processor patches and destination processor.
 689     forAll(patches, patchI)
 690     {
 691         const polyPatch& pp = patches[patchI];
 692
 693         if (isA<processorPolyPatch>(pp))
 694         {
 695             const processorPolyPatch& procPatch =
 696                 refCast<const processorPolyPatch>(pp);
 697
 698             // Labels of faces on this side
 699             labelList meshFaceLabels(pp.size());
 700             forAll(meshFaceLabels, i)
 701             {
 702                 meshFaceLabels[i] = pp.start()+i;
 703             }
 704
 705             // Which processor they will end up on
 706             const labelList newProc
 707             (
 708                 IndirectList<label>(distribution, pp.faceCells())
 709             );
 710
 711             OPstream toNeighbour(Pstream::blocking, procPatch.neighbProcNo());
 712
 713             toNeighbour << meshFaceLabels << newProc;
 714         }
 715     }
 716
 717     // Receive meshFace labels and destination processors of processor faces.
 718     forAll(patches, patchI)
 719     {
 720         const polyPatch& pp = patches[patchI];
 721
 722         label offset = pp.start() - mesh_.nInternalFaces();
 723
 724         if (isA<processorPolyPatch>(pp))
 725         {
 726             const processorPolyPatch& procPatch =
 727                 refCast<const processorPolyPatch>(pp);
 728
 729             // Receive the data
 730             IPstream fromNeighbour(Pstream::blocking, procPatch.neighbProcNo());
 731             labelList nbrFaces(fromNeighbour);
 732             labelList nbrNewProc(fromNeighbour);
 733
 734             // Check which of the two faces we store.
 735
 736             if (Pstream::myProcNo() < procPatch.neighbProcNo())
 737             {
 738                 // Use my local face labels
 739                 forAll(pp, i)
 740                 {
 741                     sourceFace[offset + i] = pp.start()+i;
 742                     sourceProc[offset + i] = Pstream::myProcNo();
 743                     sourceNewProc[offset + i] = nbrNewProc[i];
 744                 }
 745             }
 746             else
 747             {
 748                 // Use my neighbours face labels
 749                 forAll(pp, i)
 750                 {
 751                     sourceFace[offset + i] = nbrFaces[i];
 752                     sourceProc[offset + i] = procPatch.neighbProcNo();
 753                     sourceNewProc[offset + i] = nbrNewProc[i];
 754                 }
 755             }
 756         }
 757         else
 758         {
 759             // Normal (physical) boundary
 760             forAll(pp, i)
 761             {
 762                 sourceFace[offset + i] = patchI;
 763                 sourceProc[offset + i] = -1;
 764                 sourceNewProc[offset + i] = -1;
 765             }
 766         }
 767     }
 768 }
 769
 770
 771 // Subset the neighbourCell/neighbourProc fields
 772 void Foam::fvMeshDistribute::subsetBoundaryData
 773 (
 774     const fvMesh& mesh,
 775     const labelList& faceMap,
 776     const labelList& cellMap,
 777
 778     const labelList& oldDistribution,
 779     const labelList& oldFaceOwner,
 780     const labelList& oldFaceNeighbour,
 781     const label oldInternalFaces,
 782
 783     const labelList& sourceFace,
 784     const labelList& sourceProc,
 785     const labelList& sourceNewProc,
 786
 787     labelList& subFace,
 788     labelList& subProc,
 789     labelList& subNewProc
 790 )
 791 {
 792     subFace.setSize(mesh.nFaces() - mesh.nInternalFaces());
 793     subProc.setSize(mesh.nFaces() - mesh.nInternalFaces());
 794     subNewProc.setSize(mesh.nFaces() - mesh.nInternalFaces());
 795
 796     forAll(subFace, newBFaceI)
 797     {
 798         label newFaceI = newBFaceI + mesh.nInternalFaces();
 799
 800         label oldFaceI = faceMap[newFaceI];
 801
 802         // Was oldFaceI internal face? If so which side did we get.
 803         if (oldFaceI < oldInternalFaces)
 804         {
 805             subFace[newBFaceI] = oldFaceI;
 806             subProc[newBFaceI] = Pstream::myProcNo();
 807
 808             label oldOwn = oldFaceOwner[oldFaceI];
 809             label oldNei = oldFaceNeighbour[oldFaceI];
 810
 811             if (oldOwn == cellMap[mesh.faceOwner()[newFaceI]])
 812             {
 813                 // We kept the owner side. Where does the neighbour move to?
 814                 subNewProc[newBFaceI] = oldDistribution[oldNei];
 815             }
 816             else
 817             {
 818                 // We kept the neighbour side.
 819                 subNewProc[newBFaceI] = oldDistribution[oldOwn];
 820             }
 821         }
 822         else
 823         {
 824             // Was boundary face. Take over boundary information
 825             label oldBFaceI = oldFaceI - oldInternalFaces;
 826
 827             subFace[newBFaceI] = sourceFace[oldBFaceI];
 828             subProc[newBFaceI] = sourceProc[oldBFaceI];
 829             subNewProc[newBFaceI] = sourceNewProc[oldBFaceI];
 830         }
 831     }
 832 }
 833
 834
 835 // Find cells on mesh whose faceID/procID match the neighbour cell/proc of
 836 // domainMesh. Store the matching face.
 837 void Foam::fvMeshDistribute::findCouples
 838 (
 839     const primitiveMesh& mesh,
 840     const labelList& sourceFace,
 841     const labelList& sourceProc,
 842
 843     const label domain,
 844     const primitiveMesh& domainMesh,
 845     const labelList& domainFace,
 846     const labelList& domainProc,
 847
 848     labelList& masterCoupledFaces,
 849     labelList& slaveCoupledFaces
 850 )
 851 {
 852     // Store domain neighbour as map so we can easily look for pair
 853     // with same face+proc.
 854     HashTable<label, labelPair, labelPairHash> map(domainFace.size());
 855
 856     forAll(domainFace, bFaceI)
 857     {
 858         map.insert(labelPair(domainFace[bFaceI], domainProc[bFaceI]), bFaceI);
 859     }
 860
 861
 862     // Try to match mesh data.
 863
 864     masterCoupledFaces.setSize(domainFace.size());
 865     slaveCoupledFaces.setSize(domainFace.size());
 866     label coupledI = 0;
 867
 868     forAll(sourceFace, bFaceI)
 869     {
 870         if (sourceProc[bFaceI] != -1)
 871         {
 872             labelPair myData(sourceFace[bFaceI], sourceProc[bFaceI]);
 873
 874             HashTable<label, labelPair, labelPairHash>::const_iterator iter =
 875                 map.find(myData);
 876
 877             if (iter != map.end())
 878             {
 879                 label nbrBFaceI = iter();
 880
 881                 masterCoupledFaces[coupledI] = mesh.nInternalFaces() + bFaceI;
 882                 slaveCoupledFaces[coupledI] =
 883                     domainMesh.nInternalFaces()
 884                   + nbrBFaceI;
 885
 886                 coupledI++;
 887             }
 888         }
 889     }
 890
 891     masterCoupledFaces.setSize(coupledI);
 892     slaveCoupledFaces.setSize(coupledI);
 893
 894     if (debug)
 895     {
 896         Pout<< "findCouples : found " << coupledI
 897             << " faces that will be stitched" << nl << endl;
 898     }
 899 }
 900
 901
 902 // Map data on boundary faces to new mesh (resulting from adding two meshes)
 903 Foam::labelList Foam::fvMeshDistribute::mapBoundaryData
 904 (
 905     const primitiveMesh& mesh,      // mesh after adding
 906     const mapAddedPolyMesh& map,
 907     const labelList& boundaryData0, // mesh before adding
 908     const label nInternalFaces1,
 909     const labelList& boundaryData1  // added mesh
 910 )
 911 {
 912     labelList newBoundaryData(mesh.nFaces() - mesh.nInternalFaces());
 913
 914     forAll(boundaryData0, oldBFaceI)
 915     {
 916         label newFaceI = map.oldFaceMap()[oldBFaceI + map.nOldInternalFaces()];
 917
 918         // Face still exists (is nessecary?) and still boundary face
 919         if (newFaceI >= 0 && newFaceI >= mesh.nInternalFaces())
 920         {
 921             newBoundaryData[newFaceI - mesh.nInternalFaces()] =
 922                 boundaryData0[oldBFaceI];
 923         }
 924     }
 925
 926     forAll(boundaryData1, addedBFaceI)
 927     {
 928         label newFaceI = map.addedFaceMap()[addedBFaceI + nInternalFaces1];
 929
 930         if (newFaceI >= 0 && newFaceI >= mesh.nInternalFaces())
 931         {
 932             newBoundaryData[newFaceI - mesh.nInternalFaces()] =
 933                 boundaryData1[addedBFaceI];
 934         }
 935     }
 936
 937     return newBoundaryData;
 938 }
 939
 940
 941 // Remove cells. Add all exposed faces to patch oldInternalPatchI
 942 Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::doRemoveCells
 943 (
 944     const labelList& cellsToRemove,
 945     const label oldInternalPatchI
 946 )
 947 {
 948     // Mesh change engine
 949     polyTopoChange meshMod(mesh_);
 950
 951     // Cell removal topo engine. Do NOT synchronize parallel since
 952     // we are doing a local cell removal.
 953     removeCells cellRemover(mesh_, false);
 954
 955     // Get all exposed faces
 956     labelList exposedFaces(cellRemover.getExposedFaces(cellsToRemove));
 957
 958     // Insert the topo changes
 959     cellRemover.setRefinement
 960     (
 961         cellsToRemove,
 962         exposedFaces,
 963         labelList(exposedFaces.size(), oldInternalPatchI),  // patch for exposed
 964                                                             // faces.
 965         meshMod
 966     );
 967
 968     // Change the mesh. No inflation. Note: no parallel comms allowed.
 969     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, false);
 970
 971     // Update fields
 972     mesh_.updateMesh(map);
 973
 974     // Move mesh (since morphing does not do this)
 975     if (map().hasMotionPoints())
 976     {
 977         mesh_.movePoints(map().preMotionPoints());
 978     }
 979
 980     return map;
 981 }
 982
 983
 984 // Delete and add processor patches. Changes mesh and returns per neighbour proc
 985 // the processor patchID.
 986 void Foam::fvMeshDistribute::addProcPatches
 987 (
 988     const labelList& neighbourNewProc,   // processor that neighbour is on
 989     labelList& procPatchID
 990 )
 991 {
 992     // Now use the neighbourFace/Proc to repatch the mesh. These two lists
 993     // contain for all current boundary faces the global patchID (for non-proc
 994     // patch) or the processor.
 995
 996     labelList procPatchSizes(Pstream::nProcs(), 0);
 997
 998     forAll(neighbourNewProc, bFaceI)
 999     {
1000         if (neighbourNewProc[bFaceI] != -1)
1001         {
1002             procPatchSizes[neighbourNewProc[bFaceI]]++;
1003         }
1004     }
1005
1006     // Per neighbour processor the label of the processor patch
1007     procPatchID.setSize(Pstream::nProcs());
1008
1009     forAll(procPatchSizes, procI)
1010     {
1011         if (procPatchSizes[procI] > 0)
1012         {
1013             const word patchName =
1014                 "procBoundary"
1015               + name(Pstream::myProcNo())
1016               + "to"
1017               + name(procI);
1018
1019
1020             procPatchID[procI] = addProcPatch(patchName, procI);
1021             addPatchFields<volScalarField>
1022             (
1023                 processorFvPatchField<scalar>::typeName
1024             );
1025             addPatchFields<volVectorField>
1026             (
1027                 processorFvPatchField<vector>::typeName
1028             );
1029             addPatchFields<volSphericalTensorField>
1030             (
1031                 processorFvPatchField<sphericalTensor>::typeName
1032             );
1033             addPatchFields<volSymmTensorField>
1034             (
1035                 processorFvPatchField<symmTensor>::typeName
1036             );
1037             addPatchFields<volTensorField>
1038             (
1039                 processorFvPatchField<tensor>::typeName
1040             );
1041
1042             addPatchFields<surfaceScalarField>
1043             (
1044                 processorFvPatchField<scalar>::typeName
1045             );
1046             addPatchFields<surfaceVectorField>
1047             (
1048                 processorFvPatchField<vector>::typeName
1049             );
1050             addPatchFields<surfaceSphericalTensorField>
1051             (
1052                 processorFvPatchField<sphericalTensor>::typeName
1053             );
1054             addPatchFields<surfaceSymmTensorField>
1055             (
1056                 processorFvPatchField<symmTensor>::typeName
1057             );
1058             addPatchFields<surfaceTensorField>
1059             (
1060                 processorFvPatchField<tensor>::typeName
1061             );
1062         }
1063         else
1064         {
1065             procPatchID[procI] = -1;
1066         }
1067     }
1068 }
1069
1070
1071 // Get boundary faces to be repatched. Is -1 or new patchID
1072 Foam::labelList Foam::fvMeshDistribute::getProcBoundaryPatch
1073 (
1074     const labelList& neighbourNewProc,  // new processor per boundary face
1075     const labelList& procPatchID        // patchID
1076 )
1077 {
1078     labelList patchIDs(neighbourNewProc);
1079
1080     forAll(neighbourNewProc, bFaceI)
1081     {
1082         if (neighbourNewProc[bFaceI] != -1)
1083         {
1084             label nbrProc = neighbourNewProc[bFaceI];
1085
1086             patchIDs[bFaceI] = procPatchID[nbrProc];
1087         }
1088         else
1089         {
1090             patchIDs[bFaceI] = -1;
1091         }
1092     }
1093     return patchIDs;
1094 }
1095
1096
1097 // Send mesh and coupling data.
1098 void Foam::fvMeshDistribute::sendMesh
1099 (
1100     const label domain,
1101     const fvMesh& mesh,
1102
1103     const wordList& pointZoneNames,
1104     const wordList& faceZoneNames,
1105     const wordList& cellZoneNames,
1106
1107     const labelList& sourceFace,
1108     const labelList& sourceProc,
1109     const labelList& sourceNewProc
1110 )
1111 {
1112     if (debug)
1113     {
1114         Pout<< "Sending to domain " << domain << nl
1115             << "    nPoints:" << mesh.nPoints() << nl
1116             << "    nFaces:" << mesh.nFaces() << nl
1117             << "    nCells:" << mesh.nCells() << nl
1118             << "    nPatches:" << mesh.boundaryMesh().size() << nl
1119             << endl;
1120     }
1121
1122     // Assume sparse point zones. Get contents in merged-zone indices.
1123     labelListList zonePoints(pointZoneNames.size());
1124     {
1125         const pointZoneMesh& pointZones = mesh.pointZones();
1126
1127         forAll(pointZoneNames, nameI)
1128         {
1129             label myZoneID = pointZones.findZoneID(pointZoneNames[nameI]);
1130
1131             if (myZoneID != -1)
1132             {
1133                 zonePoints[nameI] = pointZones[myZoneID];
1134             }
1135         }
1136     }
1137
1138     // Assume sparse face zones
1139     labelListList zoneFaces(faceZoneNames.size());
1140     boolListList zoneFaceFlip(faceZoneNames.size());
1141     {
1142         const faceZoneMesh& faceZones = mesh.faceZones();
1143
1144         forAll(faceZoneNames, nameI)
1145         {
1146             label myZoneID = faceZones.findZoneID(faceZoneNames[nameI]);
1147
1148             if (myZoneID != -1)
1149             {
1150                 zoneFaces[nameI] = faceZones[myZoneID];
1151                 zoneFaceFlip[nameI] = faceZones[myZoneID].flipMap();
1152             }
1153         }
1154     }
1155
1156     // Assume sparse cell zones
1157     labelListList zoneCells(cellZoneNames.size());
1158     {
1159         const cellZoneMesh& cellZones = mesh.cellZones();
1160
1161         forAll(cellZoneNames, nameI)
1162         {
1163             label myZoneID = cellZones.findZoneID(cellZoneNames[nameI]);
1164
1165             if (myZoneID != -1)
1166             {
1167                 zoneCells[nameI] = cellZones[myZoneID];
1168             }
1169         }
1170     }
1171     ////- Assume full cell zones
1172     //labelList cellZoneID;
1173     //if (hasCellZones)
1174     //{
1175     //    cellZoneID.setSize(mesh.nCells());;
1176     //    cellZoneID = -1;
1177     //
1178     //    const cellZoneMesh& cellZones = mesh.cellZones();
1179     //
1180     //    forAll(cellZones, zoneI)
1181     //    {
1182     //        IndirectList<label>(cellZoneID, cellZones[zoneI]) = zoneI;
1183     //    }
1184     //}
1185
1186     // Send
1187     OPstream toDomain(Pstream::blocking, domain);
1188     toDomain
1189         << mesh.points()
1190         << mesh.faces()
1191         << mesh.faceOwner()
1192         << mesh.faceNeighbour()
1193         << mesh.boundaryMesh()
1194
1195         << zonePoints
1196         << zoneFaces
1197         << zoneFaceFlip
1198         << zoneCells
1199
1200         << sourceFace
1201         << sourceProc
1202         << sourceNewProc;
1203 }
1204
1205
1206 // Receive mesh. Opposite of sendMesh
1207 Foam::autoPtr<Foam::fvMesh> Foam::fvMeshDistribute::receiveMesh
1208 (
1209     const label domain,
1210     const wordList& pointZoneNames,
1211     const wordList& faceZoneNames,
1212     const wordList& cellZoneNames,
1213     const Time& runTime,
1214     labelList& domainSourceFace,
1215     labelList& domainSourceProc,
1216     labelList& domainSourceNewProc
1217 )
1218 {
1219     IPstream fromNbr(Pstream::blocking, domain);
1220
1221     pointField domainPoints(fromNbr);
1222     faceList domainFaces(fromNbr);
1223     labelList domainAllOwner(fromNbr);
1224     labelList domainAllNeighbour(fromNbr);
1225     PtrList<entry> patchEntries(fromNbr);
1226
1227     labelListList zonePoints(fromNbr);
1228     labelListList zoneFaces(fromNbr);
1229     boolListList zoneFaceFlip(fromNbr);
1230     labelListList zoneCells(fromNbr);
1231
1232     fromNbr
1233         >> domainSourceFace
1234         >> domainSourceProc
1235         >> domainSourceNewProc;
1236
1237     // Construct fvMesh
1238     autoPtr<fvMesh> domainMeshPtr
1239     (
1240         new fvMesh
1241         (
1242             IOobject
1243             (
1244                 fvMesh::defaultRegion,
1245                 runTime.timeName(),
1246                 runTime,
1247                 IOobject::NO_READ
1248             ),
1249             domainPoints,
1250             domainFaces,
1251             domainAllOwner,
1252             domainAllNeighbour,
1253             false                   // no parallel comms
1254         )
1255     );
1256     fvMesh& domainMesh = domainMeshPtr();
1257
1258     List<polyPatch*> patches(patchEntries.size());
1259
1260     forAll(patchEntries, patchI)
1261     {
1262         patches[patchI] = polyPatch::New
1263         (
1264             patchEntries[patchI].keyword(),
1265             patchEntries[patchI].dict(),
1266             patchI,
1267             domainMesh.boundaryMesh()
1268         ).ptr();
1269     }
1270     // Add patches; no parallel comms
1271     domainMesh.addFvPatches(patches, false);
1272
1273
1274     // Construct zones
1275     List<pointZone*> pZonePtrs(pointZoneNames.size());
1276     forAll(pZonePtrs, i)
1277     {
1278         pZonePtrs[i] = new pointZone
1279         (
1280             pointZoneNames[i],
1281             zonePoints[i],
1282             i,
1283             domainMesh.pointZones()
1284         );
1285     }
1286
1287     List<faceZone*> fZonePtrs(faceZoneNames.size());
1288     forAll(fZonePtrs, i)
1289     {
1290         fZonePtrs[i] = new faceZone
1291         (
1292             faceZoneNames[i],
1293             zoneFaces[i],
1294             zoneFaceFlip[i],
1295             i,
1296             domainMesh.faceZones()
1297         );
1298     }
1299
1300     List<cellZone*> cZonePtrs(cellZoneNames.size());
1301     forAll(cZonePtrs, i)
1302     {
1303         cZonePtrs[i] = new cellZone
1304         (
1305             cellZoneNames[i],
1306             zoneCells[i],
1307             i,
1308             domainMesh.cellZones()
1309         );
1310     }
1311     domainMesh.addZones(pZonePtrs, fZonePtrs, cZonePtrs);
1312
1313     return domainMeshPtr;
1314 }
1315
1316
1317 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
1318
1319 // Construct from components
1320 Foam::fvMeshDistribute::fvMeshDistribute(fvMesh& mesh, const scalar mergeTol)
1321 :
1322     mesh_(mesh),
1323     mergeTol_(mergeTol)
1324 {}
1325
1326
1327 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
1328
1329 Foam::labelList Foam::fvMeshDistribute::countCells
1330 (
1331     const labelList& distribution
1332 )
1333 {
1334     labelList nCells(Pstream::nProcs(), 0);
1335     forAll(distribution, cellI)
1336     {
1337         label newProc = distribution[cellI];
1338
1339         if (newProc < 0 || newProc >= Pstream::nProcs())
1340         {
1341             FatalErrorIn("fvMeshDistribute::distribute(const labelList&)")
1342                 << "Distribution should be in range 0.." << Pstream::nProcs()-1
1343                 << endl
1344                 << "At index " << cellI << " distribution:" << newProc
1345                 << abort(FatalError);
1346         }
1347         nCells[newProc]++;
1348     }
1349     return nCells;
1350 }
1351
1352
1353 Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
1354 (
1355     const labelList& distribution
1356 )
1357 {
1358     // Some checks on distribution
1359     if (distribution.size() != mesh_.nCells())
1360     {
1361         FatalErrorIn("fvMeshDistribute::distribute(const labelList&)")
1362             << "Size of distribution:"
1363             << distribution.size() << " mesh nCells:" << mesh_.nCells()
1364             << abort(FatalError);
1365     }
1366
1367
1368     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1369
1370     // Check all processors have same non-proc patches in same order.
1371     if (patches.checkParallelSync(true))
1372     {
1373         FatalErrorIn("fvMeshDistribute::distribute(const labelList&)")
1374             << "This application requires all non-processor patches"
1375             << " to be present in the same order on all patches" << nl
1376             << "followed by the processor patches (which of course are unique)."
1377             << nl
1378             << "Local patches:" << mesh_.boundaryMesh().names()
1379             << abort(FatalError);
1380     }
1381
1382     // Save some data for mapping later on
1383     const label nOldPoints(mesh_.nPoints());
1384     const label nOldFaces(mesh_.nFaces());
1385     const label nOldCells(mesh_.nCells());
1386     labelList oldPatchStarts(patches.size());
1387     labelList oldPatchNMeshPoints(patches.size());
1388     forAll(patches, patchI)
1389     {
1390         oldPatchStarts[patchI] = patches[patchI].start();
1391         oldPatchNMeshPoints[patchI] = patches[patchI].nPoints();
1392     }
1393
1394
1395
1396     // Short circuit trivial case.
1397     if (!Pstream::parRun())
1398     {
1399         // Collect all maps and return
1400         return autoPtr<mapDistributePolyMesh>
1401         (
1402             new mapDistributePolyMesh
1403             (
1404                 mesh_,
1405
1406                 nOldPoints,
1407                 nOldFaces,
1408                 nOldCells,
1409                 oldPatchStarts,
1410                 oldPatchNMeshPoints,
1411
1412                 labelListList(1, identity(mesh_.nPoints())),//subPointMap
1413                 labelListList(1, identity(mesh_.nFaces())), //subFaceMap
1414                 labelListList(1, identity(mesh_.nCells())), //subCellMap
1415                 labelListList(1, identity(patches.size())), //subPatchMap
1416
1417                 labelListList(1, identity(mesh_.nPoints())),//constructPointMap
1418                 labelListList(1, identity(mesh_.nFaces())), //constructFaceMap
1419                 labelListList(1, identity(mesh_.nCells())), //constructCellMap
1420                 labelListList(1, identity(patches.size()))  //constructPatchMap
1421             )
1422         );
1423     }
1424
1425
1426     // Collect any zone names
1427     const wordList pointZoneNames(mergeWordList(mesh_.pointZones().names()));
1428     const wordList faceZoneNames(mergeWordList(mesh_.faceZones().names()));
1429     const wordList cellZoneNames(mergeWordList(mesh_.cellZones().names()));
1430
1431
1432     // Local environment of all boundary faces
1433     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1434
1435     // A face is uniquely defined by
1436     //  - proc
1437     //  - local face no
1438     //
1439     // To glue the parts of meshes which can get sent from anywhere we
1440     // need to know on boundary faces what the above tuple on both sides is.
1441     // So we need to maintain:
1442     //  - original face
1443     //  - original processor id (= trivial)
1444     // For coupled boundaries (where the faces are 'duplicate') we take the
1445     // lowest numbered processor as the data to store.
1446     //
1447     // Additionally to create the procboundaries we need to know where the owner
1448     // cell on the other side now is: newNeighbourProc.
1449     //
1450
1451     // physical boundary:
1452     //     sourceProc = -1
1453     //     sourceNewProc = -1
1454     //     sourceFace = patchID
1455     // coupled boundary:
1456     //     sourceProc = proc
1457     //     sourceNewProc = distribution[cell on proc]
1458     //     sourceFace = face
1459     labelList sourceFace;
1460     labelList sourceProc;
1461     labelList sourceNewProc;
1462     getNeighbourData(distribution, sourceFace, sourceProc, sourceNewProc);
1463
1464
1465     // Remove meshPhi. Since this would otherwise dissappear anyway
1466     // during topo changes and we have to guarantee that all the fields
1467     // can be sent.
1468     mesh_.clearOut();
1469     mesh_.resetMotion();
1470
1471     const wordList volScalars(mesh_.names(volScalarField::typeName));
1472     checkEqualWordList(volScalars);
1473     const wordList volVectors(mesh_.names(volVectorField::typeName));
1474     checkEqualWordList(volVectors);
1475     const wordList volSphereTensors
1476     (
1477         mesh_.names(volSphericalTensorField::typeName)
1478     );
1479     checkEqualWordList(volSphereTensors);
1480     const wordList volSymmTensors(mesh_.names(volSymmTensorField::typeName));
1481     checkEqualWordList(volSymmTensors);
1482     const wordList volTensors(mesh_.names(volTensorField::typeName));
1483     checkEqualWordList(volTensors);
1484
1485     const wordList surfScalars(mesh_.names(surfaceScalarField::typeName));
1486     checkEqualWordList(surfScalars);
1487     const wordList surfVectors(mesh_.names(surfaceVectorField::typeName));
1488     checkEqualWordList(surfVectors);
1489     const wordList surfSphereTensors
1490     (
1491         mesh_.names(surfaceSphericalTensorField::typeName)
1492     );
1493     checkEqualWordList(surfSphereTensors);
1494     const wordList surfSymmTensors
1495     (
1496         mesh_.names(surfaceSymmTensorField::typeName)
1497     );
1498     checkEqualWordList(surfSymmTensors);
1499     const wordList surfTensors(mesh_.names(surfaceTensorField::typeName));
1500     checkEqualWordList(surfTensors);
1501
1502
1503     // Find patch to temporarily put exposed and processor faces into.
1504     label oldInternalPatchI = findNonEmptyPatch();
1505
1506
1507
1508     // Delete processor patches, starting from the back. Move all faces into
1509     // oldInternalPatchI.
1510     labelList repatchFaceMap;
1511     {
1512         autoPtr<mapPolyMesh> repatchMap = deleteProcPatches(oldInternalPatchI);
1513
1514         // Store face map (only face ordering that changed)
1515         repatchFaceMap = repatchMap().faceMap();
1516
1517         // Reorder all boundary face data (sourceProc, sourceFace etc.)
1518         labelList bFaceMap
1519         (
1520             SubList<label>
1521             (
1522                 repatchMap().reverseFaceMap(),
1523                 mesh_.nFaces() - mesh_.nInternalFaces(),
1524                 mesh_.nInternalFaces()
1525             )
1526           - mesh_.nInternalFaces()
1527         );
1528
1529         inplaceReorder(bFaceMap, sourceFace);
1530         inplaceReorder(bFaceMap, sourceProc);
1531         inplaceReorder(bFaceMap, sourceNewProc);
1532     }
1533
1534
1535
1536     // Print a bit.
1537     if (debug)
1538     {
1539         Pout<< nl << "MESH WITH PROC PATCHES DELETED:" << endl;
1540         printMeshInfo(mesh_);
1541         printFieldInfo<volScalarField>(mesh_);
1542         printFieldInfo<volVectorField>(mesh_);
1543         printFieldInfo<volSphericalTensorField>(mesh_);
1544         printFieldInfo<volSymmTensorField>(mesh_);
1545         printFieldInfo<volTensorField>(mesh_);
1546         printFieldInfo<surfaceScalarField>(mesh_);
1547         printFieldInfo<surfaceVectorField>(mesh_);
1548         printFieldInfo<surfaceSphericalTensorField>(mesh_);
1549         printFieldInfo<surfaceSymmTensorField>(mesh_);
1550         printFieldInfo<surfaceTensorField>(mesh_);
1551         Pout<< nl << endl;
1552     }
1553
1554
1555
1556     // Maps from subsetted mesh (that is sent) back to original maps
1557     labelListList subCellMap(Pstream::nProcs());
1558     labelListList subFaceMap(Pstream::nProcs());
1559     labelListList subPointMap(Pstream::nProcs());
1560     labelListList subPatchMap(Pstream::nProcs());
1561     // Maps from subsetted mesh to reconstructed mesh
1562     labelListList constructCellMap(Pstream::nProcs());
1563     labelListList constructFaceMap(Pstream::nProcs());
1564     labelListList constructPointMap(Pstream::nProcs());
1565     labelListList constructPatchMap(Pstream::nProcs());
1566
1567
1568
1569
1570     // Find out schedule
1571     // ~~~~~~~~~~~~~~~~~
1572
1573     labelListList nSendCells(Pstream::nProcs());
1574     nSendCells[Pstream::myProcNo()] = countCells(distribution);
1575     Pstream::gatherList(nSendCells);
1576     Pstream::scatterList(nSendCells);
1577
1578
1579
1580     // What to send to neighbouring domains
1581     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1582
1583     forAll(nSendCells[Pstream::myProcNo()], recvProc)
1584     {
1585         if
1586         (
1587             recvProc != Pstream::myProcNo()
1588          && nSendCells[Pstream::myProcNo()][recvProc] > 0
1589         )
1590         {
1591             // Send to recvProc
1592
1593             if (debug)
1594             {
1595                 Pout<< nl
1596                     << "SUBSETTING FOR DOMAIN " << recvProc
1597                     << " cells to send:"
1598                     << nSendCells[Pstream::myProcNo()][recvProc]
1599                     << nl << endl;
1600             }
1601
1602             // Mesh subsetting engine
1603             fvMeshSubset subsetter(mesh_);
1604
1605             // Subset the cells of the current domain.
1606             subsetter.setLargeCellSubset
1607             (
1608                 distribution,
1609                 recvProc,
1610                 oldInternalPatchI,  // oldInternalFaces patch
1611                 false               // no parallel sync
1612             );
1613
1614             subCellMap[recvProc] = subsetter.cellMap();
1615             subFaceMap[recvProc] = renumber
1616             (
1617                 repatchFaceMap,
1618                 subsetter.faceMap()
1619             );
1620             subPointMap[recvProc] = subsetter.pointMap();
1621             subPatchMap[recvProc] = subsetter.patchMap();
1622
1623
1624             // Subset the boundary fields (owner/neighbour/processor)
1625             labelList procSourceFace;
1626             labelList procSourceProc;
1627             labelList procSourceNewProc;
1628
1629             subsetBoundaryData
1630             (
1631                 subsetter.subMesh(),
1632                 subsetter.faceMap(),        // from subMesh to mesh
1633                 subsetter.cellMap(),        //      ,,      ,,
1634
1635                 distribution,               // old mesh distribution
1636                 mesh_.faceOwner(),          // old owner
1637                 mesh_.faceNeighbour(),
1638                 mesh_.nInternalFaces(),
1639
1640                 sourceFace,
1641                 sourceProc,
1642                 sourceNewProc,
1643
1644                 procSourceFace,
1645                 procSourceProc,
1646                 procSourceNewProc
1647             );
1648
1649             // Send to neighbour
1650             sendMesh
1651             (
1652                 recvProc,
1653                 subsetter.subMesh(),
1654
1655                 pointZoneNames,
1656                 faceZoneNames,
1657                 cellZoneNames,
1658
1659                 procSourceFace,
1660                 procSourceProc,
1661                 procSourceNewProc
1662             );
1663             sendFields<volScalarField>(recvProc, volScalars, subsetter);
1664             sendFields<volVectorField>(recvProc, volVectors, subsetter);
1665             sendFields<volSphericalTensorField>
1666             (
1667                 recvProc,
1668                 volSphereTensors,
1669                 subsetter
1670             );
1671             sendFields<volSymmTensorField>
1672             (
1673                 recvProc,
1674                 volSymmTensors,
1675                 subsetter
1676             );
1677             sendFields<volTensorField>(recvProc, volTensors, subsetter);
1678
1679             sendFields<surfaceScalarField>(recvProc, surfScalars, subsetter);
1680             sendFields<surfaceVectorField>(recvProc, surfVectors, subsetter);
1681             sendFields<surfaceSphericalTensorField>
1682             (
1683                 recvProc,
1684                 surfSphereTensors,
1685                 subsetter
1686             );
1687             sendFields<surfaceSymmTensorField>
1688             (
1689                 recvProc,
1690                 surfSymmTensors,
1691                 subsetter
1692             );
1693             sendFields<surfaceTensorField>(recvProc, surfTensors, subsetter);
1694         }
1695     }
1696
1697
1698
1699     // Subset the part that stays
1700     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
1701
1702     {
1703         // Save old mesh maps before changing mesh
1704         const labelList oldFaceOwner(mesh_.faceOwner());
1705         const labelList oldFaceNeighbour(mesh_.faceNeighbour());
1706         const label oldInternalFaces = mesh_.nInternalFaces();
1707
1708         // Remove cells.
1709         autoPtr<mapPolyMesh> subMap
1710         (
1711             doRemoveCells
1712             (
1713                 select(false, distribution, Pstream::myProcNo()),
1714                 oldInternalPatchI
1715             )
1716         );
1717
1718         // Addressing from subsetted mesh
1719         subCellMap[Pstream::myProcNo()] = subMap().cellMap();
1720         subFaceMap[Pstream::myProcNo()] = renumber
1721         (
1722             repatchFaceMap,
1723             subMap().faceMap()
1724         );
1725         subPointMap[Pstream::myProcNo()] = subMap().pointMap();
1726         subPatchMap[Pstream::myProcNo()] = identity(patches.size());
1727
1728         // Initialize all addressing into current mesh
1729         constructCellMap[Pstream::myProcNo()] = identity(mesh_.nCells());
1730         constructFaceMap[Pstream::myProcNo()] = identity(mesh_.nFaces());
1731         constructPointMap[Pstream::myProcNo()] = identity(mesh_.nPoints());
1732         constructPatchMap[Pstream::myProcNo()] = identity(patches.size());
1733
1734         // Subset the mesh data: neighbourCell/neighbourProc
1735         // fields
1736         labelList domainSourceFace;
1737         labelList domainSourceProc;
1738         labelList domainSourceNewProc;
1739
1740         subsetBoundaryData
1741         (
1742             mesh_,                          // new mesh
1743             subMap().faceMap(),             // from new to original mesh
1744             subMap().cellMap(),
1745
1746             distribution,                   // distribution before subsetting
1747             oldFaceOwner,                   // owner before subsetting
1748             oldFaceNeighbour,               // neighbour        ,,
1749             oldInternalFaces,               // nInternalFaces   ,,
1750
1751             sourceFace,
1752             sourceProc,
1753             sourceNewProc,
1754
1755             domainSourceFace,
1756             domainSourceProc,
1757             domainSourceNewProc
1758         );
1759
1760         sourceFace.transfer(domainSourceFace);
1761         sourceProc.transfer(domainSourceProc);
1762         sourceNewProc.transfer(domainSourceNewProc);
1763     }
1764
1765
1766     // Print a bit.
1767     if (debug)
1768     {
1769         Pout<< nl << "STARTING MESH:" << endl;
1770         printMeshInfo(mesh_);
1771         printFieldInfo<volScalarField>(mesh_);
1772         printFieldInfo<volVectorField>(mesh_);
1773         printFieldInfo<volSphericalTensorField>(mesh_);
1774         printFieldInfo<volSymmTensorField>(mesh_);
1775         printFieldInfo<volTensorField>(mesh_);
1776         printFieldInfo<surfaceScalarField>(mesh_);
1777         printFieldInfo<surfaceVectorField>(mesh_);
1778         printFieldInfo<surfaceSphericalTensorField>(mesh_);
1779         printFieldInfo<surfaceSymmTensorField>(mesh_);
1780         printFieldInfo<surfaceTensorField>(mesh_);
1781         Pout<< nl << endl;
1782     }
1783
1784
1785
1786     // Receive and add what was sent
1787     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1788
1789     forAll(nSendCells, sendProc)
1790     {
1791         // Did processor sendProc send anything to me?
1792         if
1793         (
1794             sendProc != Pstream::myProcNo()
1795          && nSendCells[sendProc][Pstream::myProcNo()] > 0
1796         )
1797         {
1798             if (debug)
1799             {
1800                 Pout<< nl
1801                     << "RECEIVING FROM DOMAIN " << sendProc
1802                     << " cells to receive:"
1803                     << nSendCells[sendProc][Pstream::myProcNo()]
1804                     << nl << endl;
1805             }
1806
1807             // Receive from sendProc
1808             labelList domainSourceFace;
1809             labelList domainSourceProc;
1810             labelList domainSourceNewProc;
1811
1812             // Opposite of sendMesh
1813             autoPtr<fvMesh> domainMeshPtr = receiveMesh
1814             (
1815                 sendProc,
1816                 pointZoneNames,
1817                 faceZoneNames,
1818                 cellZoneNames,
1819
1820                 const_cast<Time&>(mesh_.time()),
1821                 domainSourceFace,
1822                 domainSourceProc,
1823                 domainSourceNewProc
1824             );
1825             fvMesh& domainMesh = domainMeshPtr();
1826
1827             // Receive fields
1828             PtrList<volScalarField> vsf;
1829             receiveFields<volScalarField>
1830             (
1831                 sendProc,
1832                 volScalars,
1833                 domainMesh,
1834                 vsf
1835             );
1836
1837             PtrList<volVectorField> vvf;
1838             receiveFields<volVectorField>
1839             (
1840                 sendProc,
1841                 volVectors,
1842                 domainMesh,
1843                 vvf
1844             );
1845             PtrList<volSphericalTensorField> vsptf;
1846             receiveFields<volSphericalTensorField>
1847             (
1848                 sendProc,
1849                 volSphereTensors,
1850                 domainMesh,
1851                 vsptf
1852             );
1853             PtrList<volSymmTensorField> vsytf;
1854             receiveFields<volSymmTensorField>
1855             (
1856                 sendProc,
1857                 volSymmTensors,
1858                 domainMesh,
1859                 vsytf
1860             );
1861             PtrList<volTensorField> vtf;
1862             receiveFields<volTensorField>
1863             (
1864                 sendProc,
1865                 volTensors,
1866                 domainMesh,
1867                 vtf
1868             );
1869
1870             PtrList<surfaceScalarField> ssf;
1871             receiveFields<surfaceScalarField>
1872             (
1873                 sendProc,
1874                 surfScalars,
1875                 domainMesh,
1876                 ssf
1877             );
1878             PtrList<surfaceVectorField> svf;
1879             receiveFields<surfaceVectorField>
1880             (
1881                 sendProc,
1882                 surfVectors,
1883                 domainMesh,
1884                 svf
1885             );
1886             PtrList<surfaceSphericalTensorField> ssptf;
1887             receiveFields<surfaceSphericalTensorField>
1888             (
1889                 sendProc,
1890                 surfSphereTensors,
1891                 domainMesh,
1892                 ssptf
1893             );
1894             PtrList<surfaceSymmTensorField> ssytf;
1895             receiveFields<surfaceSymmTensorField>
1896             (
1897                 sendProc,
1898                 surfSymmTensors,
1899                 domainMesh,
1900                 ssytf
1901             );
1902             PtrList<surfaceTensorField> stf;
1903             receiveFields<surfaceTensorField>
1904             (
1905                 sendProc,
1906                 surfTensors,
1907                 domainMesh,
1908                 stf
1909             );
1910
1911
1912             constructCellMap[sendProc] = identity(domainMesh.nCells());
1913             constructFaceMap[sendProc] = identity(domainMesh.nFaces());
1914             constructPointMap[sendProc] = identity(domainMesh.nPoints());
1915             constructPatchMap[sendProc] =
1916                 identity(domainMesh.boundaryMesh().size());
1917
1918
1919             // Print a bit.
1920             if (debug)
1921             {
1922                 Pout<< nl << "RECEIVED MESH FROM:" << sendProc << endl;
1923                 printMeshInfo(domainMesh);
1924                 printFieldInfo<volScalarField>(domainMesh);
1925                 printFieldInfo<volVectorField>(domainMesh);
1926                 printFieldInfo<volSphericalTensorField>(domainMesh);
1927                 printFieldInfo<volSymmTensorField>(domainMesh);
1928                 printFieldInfo<volTensorField>(domainMesh);
1929                 printFieldInfo<surfaceScalarField>(domainMesh);
1930                 printFieldInfo<surfaceVectorField>(domainMesh);
1931                 printFieldInfo<surfaceSphericalTensorField>(domainMesh);
1932                 printFieldInfo<surfaceSymmTensorField>(domainMesh);
1933                 printFieldInfo<surfaceTensorField>(domainMesh);
1934             }
1935
1936
1937             // Now this mesh we received (from sendProc) needs to be merged
1938             // with the current mesh. On the current mesh we have for all
1939             // boundaryfaces the original face and processor. See if we can
1940             // match these up to the received domainSourceFace and
1941             // domainSourceProc.
1942             labelList masterCoupledFaces;
1943             labelList slaveCoupledFaces;
1944             findCouples
1945             (
1946                 mesh_,
1947
1948                 sourceFace,
1949                 sourceProc,
1950
1951                 sendProc,
1952                 domainMesh,
1953                 domainSourceFace,
1954                 domainSourceProc,
1955
1956                 masterCoupledFaces,
1957                 slaveCoupledFaces
1958             );
1959
1960             // Generate additional coupling info (points, edges) from
1961             // faces-that-match
1962             faceCoupleInfo couples
1963             (
1964                 mesh_,
1965                 masterCoupledFaces,
1966                 domainMesh,
1967                 slaveCoupledFaces,
1968                 mergeTol_,              // merge tolerance
1969                 true,                   // faces align
1970                 true,                   // couples are ordered already
1971                 false
1972             );
1973
1974
1975             // Add domainMesh to mesh
1976             // ~~~~~~~~~~~~~~~~~~~~~~
1977
1978             autoPtr<mapAddedPolyMesh> map = fvMeshAdder::add
1979             (
1980                 mesh_,
1981                 domainMesh,
1982                 couples,
1983                 false           // no parallel comms
1984             );
1985
1986             // Update mesh data: sourceFace,sourceProc for added
1987             // mesh.
1988
1989             sourceFace =
1990                 mapBoundaryData
1991                 (
1992                     mesh_,
1993                     map(),
1994                     sourceFace,
1995                     domainMesh.nInternalFaces(),
1996                     domainSourceFace
1997                 );
1998             sourceProc =
1999                 mapBoundaryData
2000                 (
2001                     mesh_,
2002                     map(),
2003                     sourceProc,
2004                     domainMesh.nInternalFaces(),
2005                     domainSourceProc
2006                 );
2007             sourceNewProc =
2008                 mapBoundaryData
2009                 (
2010                     mesh_,
2011                     map(),
2012                     sourceNewProc,
2013                     domainMesh.nInternalFaces(),
2014                     domainSourceNewProc
2015                 );
2016
2017             // Update all addressing so xxProcAddressing points to correct item
2018             // in masterMesh.
2019             const labelList& oldCellMap = map().oldCellMap();
2020             const labelList& oldFaceMap = map().oldFaceMap();
2021             const labelList& oldPointMap = map().oldPointMap();
2022             const labelList& oldPatchMap = map().oldPatchMap();
2023
2024             forAll(constructPatchMap, procI)
2025             {
2026                 if (procI != sendProc && constructPatchMap[procI].size() > 0)
2027                 {
2028                     // Processor already in mesh (either myProcNo or received)
2029                     inplaceRenumber(oldCellMap, constructCellMap[procI]);
2030                     inplaceRenumber(oldFaceMap, constructFaceMap[procI]);
2031                     inplaceRenumber(oldPointMap, constructPointMap[procI]);
2032                     inplaceRenumber(oldPatchMap, constructPatchMap[procI]);
2033                 }
2034             }
2035
2036             // Added processor
2037             inplaceRenumber(map().addedCellMap(), constructCellMap[sendProc]);
2038             inplaceRenumber(map().addedFaceMap(), constructFaceMap[sendProc]);
2039             inplaceRenumber(map().addedPointMap(), constructPointMap[sendProc]);
2040             inplaceRenumber(map().addedPatchMap(), constructPatchMap[sendProc]);
2041
2042             if (debug)
2043             {
2044                 Pout<< nl << "MERGED MESH FROM:" << sendProc << endl;
2045                 printMeshInfo(mesh_);
2046                 printFieldInfo<volScalarField>(mesh_);
2047                 printFieldInfo<volVectorField>(mesh_);
2048                 printFieldInfo<volSphericalTensorField>(mesh_);
2049                 printFieldInfo<volSymmTensorField>(mesh_);
2050                 printFieldInfo<volTensorField>(mesh_);
2051                 printFieldInfo<surfaceScalarField>(mesh_);
2052                 printFieldInfo<surfaceVectorField>(mesh_);
2053                 printFieldInfo<surfaceSphericalTensorField>(mesh_);
2054                 printFieldInfo<surfaceSymmTensorField>(mesh_);
2055                 printFieldInfo<surfaceTensorField>(mesh_);
2056                 Pout<< nl << endl;
2057             }
2058         }
2059     }
2060
2061
2062     // Print a bit.
2063     if (debug)
2064     {
2065         Pout<< nl << "REDISTRIBUTED MESH:" << endl;
2066         printMeshInfo(mesh_);
2067         printFieldInfo<volScalarField>(mesh_);
2068         printFieldInfo<volVectorField>(mesh_);
2069         printFieldInfo<volSphericalTensorField>(mesh_);
2070         printFieldInfo<volSymmTensorField>(mesh_);
2071         printFieldInfo<volTensorField>(mesh_);
2072         printFieldInfo<surfaceScalarField>(mesh_);
2073         printFieldInfo<surfaceVectorField>(mesh_);
2074         printFieldInfo<surfaceSphericalTensorField>(mesh_);
2075         printFieldInfo<surfaceSymmTensorField>(mesh_);
2076         printFieldInfo<surfaceTensorField>(mesh_);
2077         Pout<< nl << endl;
2078     }
2079
2080
2081
2082     // Add processorPatches
2083     // ~~~~~~~~~~~~~~~~~~~~
2084
2085     // Per neighbour processor the patchID to it (or -1).
2086     labelList procPatchID;
2087
2088     // Add processor patches.
2089     addProcPatches(sourceNewProc, procPatchID);
2090
2091     // Put faces into correct patch. Note that we now have proper
2092     // processorPolyPatches again so repatching will take care of coupled face
2093     // ordering.
2094
2095     // Get boundary faces to be repatched. Is -1 or new patchID
2096     labelList newPatchID
2097     (
2098         getProcBoundaryPatch
2099         (
2100             sourceNewProc,
2101             procPatchID
2102         )
2103     );
2104
2105     // Change patches. Since this might change ordering of coupled faces
2106     // we also need to adapt our constructMaps.
2107     repatch(newPatchID, constructFaceMap);
2108
2109     // See if any geometrically shared points need to be merged. Note: does
2110     // parallel comms.
2111     mergeSharedPoints(constructPointMap);
2112
2113     // Bit of hack: processorFvPatchField does not get reset since created
2114     // from nothing so explicitly reset.
2115     initPatchFields<volScalarField>
2116     (
2117         processorFvPatchField<scalar>::typeName,
2118         pTraits<scalar>::zero
2119     );
2120     initPatchFields<volVectorField>
2121     (
2122         processorFvPatchField<vector>::typeName,
2123         pTraits<vector>::zero
2124     );
2125     initPatchFields<volSphericalTensorField>
2126     (
2127         processorFvPatchField<sphericalTensor>::typeName,
2128         pTraits<sphericalTensor>::zero
2129     );
2130     initPatchFields<volSymmTensorField>
2131     (
2132         processorFvPatchField<symmTensor>::typeName,
2133         pTraits<symmTensor>::zero
2134     );
2135     initPatchFields<volTensorField>
2136     (
2137         processorFvPatchField<tensor>::typeName,
2138         pTraits<tensor>::zero
2139     );
2140     initPatchFields<surfaceScalarField>
2141     (
2142         processorFvPatchField<scalar>::typeName,
2143         pTraits<scalar>::zero
2144     );
2145     initPatchFields<surfaceVectorField>
2146     (
2147         processorFvPatchField<vector>::typeName,
2148         pTraits<vector>::zero
2149     );
2150     initPatchFields<surfaceSphericalTensorField>
2151     (
2152         processorFvPatchField<sphericalTensor>::typeName,
2153         pTraits<sphericalTensor>::zero
2154     );
2155     initPatchFields<surfaceSymmTensorField>
2156     (
2157         processorFvPatchField<symmTensor>::typeName,
2158         pTraits<symmTensor>::zero
2159     );
2160     initPatchFields<surfaceTensorField>
2161     (
2162         processorFvPatchField<tensor>::typeName,
2163         pTraits<tensor>::zero
2164     );
2165
2166
2167     mesh_.setInstance(mesh_.time().timeName());
2168
2169
2170     // Print a bit
2171     if (debug)
2172     {
2173         Pout<< nl << "FINAL MESH:" << endl;
2174         printMeshInfo(mesh_);
2175         printFieldInfo<volScalarField>(mesh_);
2176         printFieldInfo<volVectorField>(mesh_);
2177         printFieldInfo<volSphericalTensorField>(mesh_);
2178         printFieldInfo<volSymmTensorField>(mesh_);
2179         printFieldInfo<volTensorField>(mesh_);
2180         printFieldInfo<surfaceScalarField>(mesh_);
2181         printFieldInfo<surfaceVectorField>(mesh_);
2182         printFieldInfo<surfaceSphericalTensorField>(mesh_);
2183         printFieldInfo<surfaceSymmTensorField>(mesh_);
2184         printFieldInfo<surfaceTensorField>(mesh_);
2185         Pout<< nl << endl;
2186     }
2187
2188     // Collect all maps and return
2189     return autoPtr<mapDistributePolyMesh>
2190     (
2191         new mapDistributePolyMesh
2192         (
2193             mesh_,
2194
2195             nOldPoints,
2196             nOldFaces,
2197             nOldCells,
2198             oldPatchStarts,
2199             oldPatchNMeshPoints,
2200
2201             subPointMap,
2202             subFaceMap,
2203             subCellMap,
2204             subPatchMap,
2205
2206             constructPointMap,
2207             constructFaceMap,
2208             constructCellMap,
2209             constructPatchMap,
2210             true                // reuse storage
2211         )
2212     );
2213 }
2214
2215
2216 // ************************************************************************* //