FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

fvMeshDistribute.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 \*----------------------------------------------------------------------------*/
00025 
00026 #include <dynamicMesh/CompactListList_dev.H>
00027 #include "fvMeshDistribute.H"
00028 #include <OpenFOAM/PstreamCombineReduceOps.H>
00029 #include <dynamicMesh/fvMeshAdder.H>
00030 #include <dynamicMesh/faceCoupleInfo.H>
00031 #include <finiteVolume/processorFvPatchField.H>
00032 #include <finiteVolume/processorFvsPatchField.H>
00033 #include <dynamicMesh/polyTopoChange.H>
00034 #include <dynamicMesh/removeCells.H>
00035 #include <dynamicMesh/polyModifyFace.H>
00036 #include <dynamicMesh/polyRemovePoint.H>
00037 #include <OpenFOAM/mergePoints.H>
00038 #include <OpenFOAM/mapDistributePolyMesh.H>
00039 #include <finiteVolume/surfaceFields.H>
00040 #include <OpenFOAM/syncTools.H>
00041 
00042 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
00043 
00044 defineTypeNameAndDebug(Foam::fvMeshDistribute, 0);
00045 
00046 
00047 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
00048 
00049 Foam::labelList Foam::fvMeshDistribute::select
00050 (
00051     const bool selectEqual,
00052     const labelList& values,
00053     const label value
00054 )
00055 {
00056     label n = 0;
00057 
00058     forAll(values, i)
00059     {
00060         if (selectEqual == (values[i] == value))
00061         {
00062             n++;
00063         }
00064     }
00065 
00066     labelList indices(n);
00067     n = 0;
00068 
00069     forAll(values, i)
00070     {
00071         if (selectEqual == (values[i] == value))
00072         {
00073             indices[n++] = i;
00074         }
00075     }
00076     return indices;
00077 }
00078 
00079 
00080 // Check all procs have same names and in exactly same order.
00081 void Foam::fvMeshDistribute::checkEqualWordList
00082 (
00083     const string& msg,
00084     const wordList& lst
00085 )
00086 {
00087     List<wordList> allNames(Pstream::nProcs());
00088     allNames[Pstream::myProcNo()] = lst;
00089     Pstream::gatherList(allNames);
00090     Pstream::scatterList(allNames);
00091 
00092     for (label procI = 1; procI < Pstream::nProcs(); procI++)
00093     {
00094         if (allNames[procI] != allNames[0])
00095         {
00096             FatalErrorIn("fvMeshDistribute::checkEqualWordList(..)")
00097                 << "When checking for equal " << msg.c_str() << " :" << endl
00098                 << "processor0 has:" << allNames[0] << endl
00099                 << "processor" << procI << " has:" << allNames[procI] << endl
00100                 << msg.c_str() << " need to be synchronised on all processors."
00101                 << exit(FatalError);
00102         }
00103     }
00104 }
00105 
00106 
00107 Foam::wordList Foam::fvMeshDistribute::mergeWordList(const wordList& procNames)
00108 {
00109     List<wordList> allNames(Pstream::nProcs());
00110     allNames[Pstream::myProcNo()] = procNames;
00111     Pstream::gatherList(allNames);
00112     Pstream::scatterList(allNames);
00113 
00114     HashSet<word> mergedNames;
00115     forAll(allNames, procI)
00116     {
00117         forAll(allNames[procI], i)
00118         {
00119             mergedNames.insert(allNames[procI][i]);
00120         }
00121     }
00122     return mergedNames.toc();
00123 }
00124 
00125 
00126 // Print some info on mesh.
00127 void Foam::fvMeshDistribute::printMeshInfo(const fvMesh& mesh)
00128 {
00129     Pout<< "Primitives:" << nl
00130         << "    points       :" << mesh.nPoints() << nl
00131         << "    internalFaces:" << mesh.nInternalFaces() << nl
00132         << "    faces        :" << mesh.nFaces() << nl
00133         << "    cells        :" << mesh.nCells() << nl;
00134 
00135     const fvBoundaryMesh& patches = mesh.boundary();
00136 
00137     Pout<< "Patches:" << endl;
00138     forAll(patches, patchI)
00139     {
00140         const polyPatch& pp = patches[patchI].patch();
00141 
00142         Pout<< "    " << patchI << " name:" << pp.name()
00143             << " size:" << pp.size()
00144             << " start:" << pp.start()
00145             << " type:" << pp.type()
00146             << endl;
00147     }
00148 
00149     if (mesh.pointZones().size())
00150     {
00151         Pout<< "PointZones:" << endl;
00152         forAll(mesh.pointZones(), zoneI)
00153         {
00154             const pointZone& pz = mesh.pointZones()[zoneI];
00155             Pout<< "    " << zoneI << " name:" << pz.name()
00156                 << " size:" << pz.size()
00157                 << endl;
00158         }
00159     }
00160     if (mesh.faceZones().size())
00161     {
00162         Pout<< "FaceZones:" << endl;
00163         forAll(mesh.faceZones(), zoneI)
00164         {
00165             const faceZone& fz = mesh.faceZones()[zoneI];
00166             Pout<< "    " << zoneI << " name:" << fz.name()
00167                 << " size:" << fz.size()
00168                 << endl;
00169         }
00170     }
00171     if (mesh.cellZones().size())
00172     {
00173         Pout<< "CellZones:" << endl;
00174         forAll(mesh.cellZones(), zoneI)
00175         {
00176             const cellZone& cz = mesh.cellZones()[zoneI];
00177             Pout<< "    " << zoneI << " name:" << cz.name()
00178                 << " size:" << cz.size()
00179                 << endl;
00180         }
00181     }
00182 }
00183 
00184 
00185 void Foam::fvMeshDistribute::printCoupleInfo
00186 (
00187     const primitiveMesh& mesh,
00188     const labelList& sourceFace,
00189     const labelList& sourceProc,
00190     const labelList& sourceNewProc
00191 )
00192 {
00193     Pout<< nl
00194         << "Current coupling info:"
00195         << endl;
00196 
00197     forAll(sourceFace, bFaceI)
00198     {
00199         label meshFaceI = mesh.nInternalFaces() + bFaceI;
00200 
00201         Pout<< "    meshFace:" << meshFaceI
00202             << " fc:" << mesh.faceCentres()[meshFaceI]
00203             << " connects to proc:" << sourceProc[bFaceI]
00204             << "/face:" << sourceFace[bFaceI]
00205             << " which will move to proc:" << sourceNewProc[bFaceI]
00206             << endl;
00207     }
00208 }
00209 
00210 
00211 // Finds (non-empty) patch that exposed internal and proc faces can be put into.
00212 Foam::label Foam::fvMeshDistribute::findNonEmptyPatch() const
00213 {
00214     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
00215 
00216     label nonEmptyPatchI = -1;
00217 
00218     forAllReverse(patches, patchI)
00219     {
00220         const polyPatch& pp = patches[patchI];
00221 
00222         if (!isA<emptyPolyPatch>(pp) && !pp.coupled())
00223         {
00224             nonEmptyPatchI = patchI;
00225             break;
00226         }
00227     }
00228 
00229     if (nonEmptyPatchI == -1)
00230     {
00231         FatalErrorIn("fvMeshDistribute::findNonEmptyPatch() const")
00232             << "Cannot find a patch which is neither of type empty nor"
00233             << " coupled in patches " << patches.names() << endl
00234             << "There has to be at least one such patch for"
00235             << " distribution to work" << abort(FatalError);
00236     }
00237 
00238     if (debug)
00239     {
00240         Pout<< "findNonEmptyPatch : using patch " << nonEmptyPatchI
00241             << " name:" << patches[nonEmptyPatchI].name()
00242             << " type:" << patches[nonEmptyPatchI].type()
00243             << " to put exposed faces into." << endl;
00244     }
00245 
00246 
00247     // Do additional test for processor patches intermingled with non-proc
00248     // patches.
00249     label procPatchI = -1;
00250 
00251     forAll(patches, patchI)
00252     {
00253         if (isA<processorPolyPatch>(patches[patchI]))
00254         {
00255             procPatchI = patchI;
00256         }
00257         else if (procPatchI != -1)
00258         {
00259             FatalErrorIn("fvMeshDistribute::findNonEmptyPatch() const")
00260                 << "Processor patches should be at end of patch list."
00261                 << endl
00262                 << "Have processor patch " << procPatchI
00263                 << " followed by non-processor patch " << patchI
00264                 << " in patches " << patches.names()
00265                 << abort(FatalError);
00266         }
00267     }
00268 
00269     return nonEmptyPatchI;
00270 }
00271 
00272 
00273 // Appends processorPolyPatch. Returns patchID.
00274 Foam::label Foam::fvMeshDistribute::addProcPatch
00275 (
00276     const word& patchName,
00277     const label nbrProc
00278 )
00279 {
00280     // Clear local fields and e.g. polyMesh globalMeshData.
00281     mesh_.clearOut();
00282 
00283 
00284     polyBoundaryMesh& polyPatches =
00285         const_cast<polyBoundaryMesh&>(mesh_.boundaryMesh());
00286     fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh_.boundary());
00287 
00288     if (polyPatches.findPatchID(patchName) != -1)
00289     {
00290         FatalErrorIn("fvMeshDistribute::addProcPatch(const word&, const label)")
00291             << "Cannot create patch " << patchName << " since already exists."
00292             << nl
00293             << "Current patch names:" << polyPatches.names()
00294             << exit(FatalError);
00295     }
00296 
00297 
00298 
00299     // Add the patch
00300     // ~~~~~~~~~~~~~
00301 
00302     label sz = polyPatches.size();
00303 
00304     // Add polyPatch
00305     polyPatches.setSize(sz+1);
00306     polyPatches.set
00307     (
00308         sz,
00309         new processorPolyPatch
00310         (
00311             patchName,
00312             0,              // size
00313             mesh_.nFaces(),
00314             sz,
00315             mesh_.boundaryMesh(),
00316             Pstream::myProcNo(),
00317             nbrProc
00318         )
00319     );
00320     fvPatches.setSize(sz+1);
00321     fvPatches.set
00322     (
00323         sz,
00324         fvPatch::New
00325         (
00326             polyPatches[sz],  // point to newly added polyPatch
00327             mesh_.boundary()
00328         )
00329     );
00330 
00331     return sz;
00332 }
00333 
00334 
00335 // Deletes last patch
00336 void Foam::fvMeshDistribute::deleteTrailingPatch()
00337 {
00338     // Clear local fields and e.g. polyMesh globalMeshData.
00339     mesh_.clearOut();
00340 
00341     polyBoundaryMesh& polyPatches =
00342         const_cast<polyBoundaryMesh&>(mesh_.boundaryMesh());
00343     fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh_.boundary());
00344 
00345     if (polyPatches.empty())
00346     {
00347         FatalErrorIn("fvMeshDistribute::deleteTrailingPatch(fvMesh&)")
00348             << "No patches in mesh"
00349             << abort(FatalError);
00350     }
00351 
00352     label sz = polyPatches.size();
00353 
00354     label nFaces = polyPatches[sz-1].size();
00355 
00356     // Remove last polyPatch
00357     if (debug)
00358     {
00359         Pout<< "deleteTrailingPatch : Removing patch " << sz-1
00360             << " : " << polyPatches[sz-1].name() << " size:" << nFaces << endl;
00361     }
00362 
00363     if (nFaces)
00364     {
00365         FatalErrorIn("fvMeshDistribute::deleteTrailingPatch()")
00366             << "There are still " << nFaces << " faces in patch to be deleted "
00367             << sz-1 << ' ' << polyPatches[sz-1].name()
00368             << abort(FatalError);
00369     }
00370 
00371     // Remove actual patch
00372     polyPatches.setSize(sz-1);
00373     fvPatches.setSize(sz-1);
00374 }
00375 
00376 
00377 // Delete all processor patches. Move any processor faces into the last
00378 // non-processor patch.
00379 Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::deleteProcPatches
00380 (
00381     const label destinationPatch
00382 )
00383 {
00384     // New patchID per boundary faces to be repatched. Is -1 (no change)
00385     // or new patchID
00386     labelList newPatchID(mesh_.nFaces() - mesh_.nInternalFaces(), -1);
00387 
00388     label nProcPatches = 0;
00389 
00390     forAll(mesh_.boundaryMesh(), patchI)
00391     {
00392         const polyPatch& pp = mesh_.boundaryMesh()[patchI];
00393 
00394         if (isA<processorPolyPatch>(pp))
00395         {
00396             if (debug)
00397             {
00398                 Pout<< "Moving all faces of patch " << pp.name()
00399                     << " into patch " << destinationPatch
00400                     << endl;
00401             }
00402 
00403             label offset = pp.start() - mesh_.nInternalFaces();
00404 
00405             forAll(pp, i)
00406             {
00407                 newPatchID[offset+i] = destinationPatch;
00408             }
00409 
00410             nProcPatches++;
00411         }
00412     }
00413 
00414     // Note: order of boundary faces been kept the same since the
00415     // destinationPatch is at the end and we have visited the patches in
00416     // incremental order.
00417     labelListList dummyFaceMaps;
00418     autoPtr<mapPolyMesh> map = repatch(newPatchID, dummyFaceMaps);
00419 
00420 
00421     // Delete (now empty) processor patches.
00422     forAllReverse(mesh_.boundaryMesh(), patchI)
00423     {
00424         const polyPatch& pp = mesh_.boundaryMesh()[patchI];
00425 
00426         if (isA<processorPolyPatch>(pp))
00427         {
00428             deleteTrailingPatch();
00429             deleteTrailingPatchFields<volScalarField>();
00430             deleteTrailingPatchFields<volVectorField>();
00431             deleteTrailingPatchFields<volSphericalTensorField>();
00432             deleteTrailingPatchFields<volSymmTensorField>();
00433             deleteTrailingPatchFields<volTensorField>();
00434 
00435             deleteTrailingPatchFields<surfaceScalarField>();
00436             deleteTrailingPatchFields<surfaceVectorField>();
00437             deleteTrailingPatchFields<surfaceSphericalTensorField>();
00438             deleteTrailingPatchFields<surfaceSymmTensorField>();
00439             deleteTrailingPatchFields<surfaceTensorField>();
00440         }
00441     }
00442 
00443     return map;
00444 }
00445 
00446 
00447 // Repatch the mesh.
00448 Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::repatch
00449 (
00450     const labelList& newPatchID,         // per boundary face -1 or new patchID
00451     labelListList& constructFaceMap
00452 )
00453 {
00454     polyTopoChange meshMod(mesh_);
00455 
00456     forAll(newPatchID, bFaceI)
00457     {
00458         if (newPatchID[bFaceI] != -1)
00459         {
00460             label faceI = mesh_.nInternalFaces() + bFaceI;
00461 
00462             label zoneID = mesh_.faceZones().whichZone(faceI);
00463             bool zoneFlip = false;
00464 
00465             if (zoneID >= 0)
00466             {
00467                 const faceZone& fZone = mesh_.faceZones()[zoneID];
00468                 zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
00469             }
00470 
00471             meshMod.setAction
00472             (
00473                 polyModifyFace
00474                 (
00475                     mesh_.faces()[faceI],       // modified face
00476                     faceI,                      // label of face
00477                     mesh_.faceOwner()[faceI],   // owner
00478                     -1,                         // neighbour
00479                     false,                      // face flip
00480                     newPatchID[bFaceI],         // patch for face
00481                     false,                      // remove from zone
00482                     zoneID,                     // zone for face
00483                     zoneFlip                    // face flip in zone
00484                 )
00485             );
00486         }
00487     }
00488 
00489 
00490     // Do mapping of fields from one patchField to the other ourselves since
00491     // is currently not supported by updateMesh.
00492 
00493     // Store boundary fields (we only do this for surfaceFields)
00494     PtrList<FieldField<fvsPatchField, scalar> > sFlds;
00495     saveBoundaryFields<scalar, surfaceMesh>(sFlds);
00496     PtrList<FieldField<fvsPatchField, vector> > vFlds;
00497     saveBoundaryFields<vector, surfaceMesh>(vFlds);
00498     PtrList<FieldField<fvsPatchField, sphericalTensor> > sptFlds;
00499     saveBoundaryFields<sphericalTensor, surfaceMesh>(sptFlds);
00500     PtrList<FieldField<fvsPatchField, symmTensor> > sytFlds;
00501     saveBoundaryFields<symmTensor, surfaceMesh>(sytFlds);
00502     PtrList<FieldField<fvsPatchField, tensor> > tFlds;
00503     saveBoundaryFields<tensor, surfaceMesh>(tFlds);
00504 
00505     // Change the mesh (no inflation). Note: parallel comms allowed.
00506     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
00507 
00508     // Update fields. No inflation, parallel sync.
00509     mesh_.updateMesh(map);
00510 
00511     // Map patch fields using stored boundary fields. Note: assumes order
00512     // of fields has not changed in object registry!
00513     mapBoundaryFields<scalar, surfaceMesh>(map, sFlds);
00514     mapBoundaryFields<vector, surfaceMesh>(map, vFlds);
00515     mapBoundaryFields<sphericalTensor, surfaceMesh>(map, sptFlds);
00516     mapBoundaryFields<symmTensor, surfaceMesh>(map, sytFlds);
00517     mapBoundaryFields<tensor, surfaceMesh>(map, tFlds);
00518 
00519 
00520     // Move mesh (since morphing does not do this)
00521     if (map().hasMotionPoints())
00522     {
00523         mesh_.movePoints(map().preMotionPoints());
00524     }
00525 
00526     // Adapt constructMaps.
00527 
00528     if (debug)
00529     {
00530         label index = findIndex(map().reverseFaceMap(), -1);
00531 
00532         if (index != -1)
00533         {
00534             FatalErrorIn
00535             (
00536                 "fvMeshDistribute::repatch(const labelList&, labelListList&)"
00537             )   << "reverseFaceMap contains -1 at index:"
00538                 << index << endl
00539                 << "This means that the repatch operation was not just"
00540                 << " a shuffle?" << abort(FatalError);
00541         }
00542     }
00543 
00544     forAll(constructFaceMap, procI)
00545     {
00546         inplaceRenumber(map().reverseFaceMap(), constructFaceMap[procI]);
00547     }
00548 
00549 
00550     return map;
00551 }
00552 
00553 
00554 // Detect shared points. Need processor patches to be present.
00555 // Background: when adding bits of mesh one can get points which
00556 // share the same position but are only detectable to be topologically
00557 // the same point when doing parallel analysis. This routine will
00558 // merge those points.
00559 Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::mergeSharedPoints
00560 (
00561     labelListList& constructPointMap
00562 )
00563 {
00564     // Find out which sets of points get merged and create a map from
00565     // mesh point to unique point.
00566     Map<label> pointToMaster
00567     (
00568         fvMeshAdder::findSharedPoints
00569         (
00570             mesh_,
00571             mergeTol_
00572         )
00573     );
00574 
00575     if (returnReduce(pointToMaster.size(), sumOp<label>()) == 0)
00576     {
00577         return autoPtr<mapPolyMesh>(NULL);
00578     }
00579 
00580     polyTopoChange meshMod(mesh_);
00581 
00582     fvMeshAdder::mergePoints(mesh_, pointToMaster, meshMod);
00583 
00584     // Change the mesh (no inflation). Note: parallel comms allowed.
00585     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
00586 
00587     // Update fields. No inflation, parallel sync.
00588     mesh_.updateMesh(map);
00589 
00590     // Adapt constructMaps for merged points.
00591     forAll(constructPointMap, procI)
00592     {
00593         labelList& constructMap = constructPointMap[procI];
00594 
00595         forAll(constructMap, i)
00596         {
00597             label oldPointI = constructMap[i];
00598 
00599             label newPointI = map().reversePointMap()[oldPointI];
00600 
00601             if (newPointI < -1)
00602             {
00603                 constructMap[i] = -newPointI-2;
00604             }
00605             else if (newPointI >= 0)
00606             {
00607                 constructMap[i] = newPointI;
00608             }
00609             else
00610             {
00611                 FatalErrorIn("fvMeshDistribute::mergeSharedPoints()")
00612                     << "Problem. oldPointI:" << oldPointI
00613                     << " newPointI:" << newPointI << abort(FatalError);
00614             }
00615         }
00616     }
00617     return map;
00618 }
00619 
00620 
00621 // Construct the local environment of all boundary faces.
00622 void Foam::fvMeshDistribute::getNeighbourData
00623 (
00624     const labelList& distribution,
00625     labelList& sourceFace,
00626     labelList& sourceProc,
00627     labelList& sourceNewProc
00628 ) const
00629 {
00630     label nBnd = mesh_.nFaces() - mesh_.nInternalFaces();
00631     sourceFace.setSize(nBnd);
00632     sourceProc.setSize(nBnd);
00633     sourceNewProc.setSize(nBnd);
00634 
00635     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
00636 
00637     // Get neighbouring meshFace labels and new processor of coupled boundaries.
00638     labelList nbrFaces(nBnd, -1);
00639     labelList nbrNewProc(nBnd, -1);
00640 
00641     forAll(patches, patchI)
00642     {
00643         const polyPatch& pp = patches[patchI];
00644 
00645         if (isA<processorPolyPatch>(pp))
00646         {
00647             label offset = pp.start() - mesh_.nInternalFaces();
00648 
00649             // Mesh labels of faces on this side
00650             forAll(pp, i)
00651             {
00652                 label bndI = offset + i;
00653                 nbrFaces[bndI] = pp.start()+i;
00654             }
00655 
00656             // Which processor they will end up on
00657             SubList<label>(nbrNewProc, pp.size(), offset).assign
00658             (
00659                 UIndirectList<label>(distribution, pp.faceCells())()
00660             );
00661         }
00662     }
00663 
00664 
00665     // Exchange the boundary data
00666     syncTools::swapBoundaryFaceList(mesh_, nbrFaces, false);
00667     syncTools::swapBoundaryFaceList(mesh_, nbrNewProc, false);
00668 
00669 
00670     forAll(patches, patchI)
00671     {
00672         const polyPatch& pp = patches[patchI];
00673         label offset = pp.start() - mesh_.nInternalFaces();
00674 
00675         if (isA<processorPolyPatch>(pp))
00676         {
00677             const processorPolyPatch& procPatch =
00678                 refCast<const processorPolyPatch>(pp);
00679 
00680             // Check which of the two faces we store.
00681 
00682             if (Pstream::myProcNo() < procPatch.neighbProcNo())
00683             {
00684                 // Use my local face labels
00685                 forAll(pp, i)
00686                 {
00687                     label bndI = offset + i;
00688                     sourceFace[bndI] = pp.start()+i;
00689                     sourceProc[bndI] = Pstream::myProcNo();
00690                     sourceNewProc[bndI] = nbrNewProc[bndI];
00691                 }
00692             }
00693             else
00694             {
00695                 // Use my neighbours face labels
00696                 forAll(pp, i)
00697                 {
00698                     label bndI = offset + i;
00699                     sourceFace[bndI] = nbrFaces[bndI];
00700                     sourceProc[bndI] = procPatch.neighbProcNo();
00701                     sourceNewProc[bndI] = nbrNewProc[bndI];
00702                 }
00703             }
00704         }
00705         else
00706         {
00707             // Normal (physical) boundary
00708             forAll(pp, i)
00709             {
00710                 label bndI = offset + i;
00711                 sourceFace[bndI] = patchI;
00712                 sourceProc[bndI] = -1;
00713                 sourceNewProc[bndI] = -1;
00714             }
00715         }
00716     }
00717 }
00718 
00719 
00720 // Subset the neighbourCell/neighbourProc fields
00721 void Foam::fvMeshDistribute::subsetBoundaryData
00722 (
00723     const fvMesh& mesh,
00724     const labelList& faceMap,
00725     const labelList& cellMap,
00726 
00727     const labelList& oldDistribution,
00728     const labelList& oldFaceOwner,
00729     const labelList& oldFaceNeighbour,
00730     const label oldInternalFaces,
00731 
00732     const labelList& sourceFace,
00733     const labelList& sourceProc,
00734     const labelList& sourceNewProc,
00735 
00736     labelList& subFace,
00737     labelList& subProc,
00738     labelList& subNewProc
00739 )
00740 {
00741     subFace.setSize(mesh.nFaces() - mesh.nInternalFaces());
00742     subProc.setSize(mesh.nFaces() - mesh.nInternalFaces());
00743     subNewProc.setSize(mesh.nFaces() - mesh.nInternalFaces());
00744 
00745     forAll(subFace, newBFaceI)
00746     {
00747         label newFaceI = newBFaceI + mesh.nInternalFaces();
00748 
00749         label oldFaceI = faceMap[newFaceI];
00750 
00751         // Was oldFaceI internal face? If so which side did we get.
00752         if (oldFaceI < oldInternalFaces)
00753         {
00754             subFace[newBFaceI] = oldFaceI;
00755             subProc[newBFaceI] = Pstream::myProcNo();
00756 
00757             label oldOwn = oldFaceOwner[oldFaceI];
00758             label oldNei = oldFaceNeighbour[oldFaceI];
00759 
00760             if (oldOwn == cellMap[mesh.faceOwner()[newFaceI]])
00761             {
00762                 // We kept the owner side. Where does the neighbour move to?
00763                 subNewProc[newBFaceI] = oldDistribution[oldNei];
00764             }
00765             else
00766             {
00767                 // We kept the neighbour side.
00768                 subNewProc[newBFaceI] = oldDistribution[oldOwn];
00769             }
00770         }
00771         else
00772         {
00773             // Was boundary face. Take over boundary information
00774             label oldBFaceI = oldFaceI - oldInternalFaces;
00775 
00776             subFace[newBFaceI] = sourceFace[oldBFaceI];
00777             subProc[newBFaceI] = sourceProc[oldBFaceI];
00778             subNewProc[newBFaceI] = sourceNewProc[oldBFaceI];
00779         }
00780     }
00781 }
00782 
00783 
00784 // Find cells on mesh whose faceID/procID match the neighbour cell/proc of
00785 // domainMesh. Store the matching face.
00786 void Foam::fvMeshDistribute::findCouples
00787 (
00788     const primitiveMesh& mesh,
00789     const labelList& sourceFace,
00790     const labelList& sourceProc,
00791 
00792     const label domain,
00793     const primitiveMesh& domainMesh,
00794     const labelList& domainFace,
00795     const labelList& domainProc,
00796 
00797     labelList& masterCoupledFaces,
00798     labelList& slaveCoupledFaces
00799 )
00800 {
00801     // Store domain neighbour as map so we can easily look for pair
00802     // with same face+proc.
00803     HashTable<label, labelPair, labelPair::Hash<> > map(domainFace.size());
00804 
00805     forAll(domainFace, bFaceI)
00806     {
00807         map.insert(labelPair(domainFace[bFaceI], domainProc[bFaceI]), bFaceI);
00808     }
00809 
00810 
00811     // Try to match mesh data.
00812 
00813     masterCoupledFaces.setSize(domainFace.size());
00814     slaveCoupledFaces.setSize(domainFace.size());
00815     label coupledI = 0;
00816 
00817     forAll(sourceFace, bFaceI)
00818     {
00819         if (sourceProc[bFaceI] != -1)
00820         {
00821             labelPair myData(sourceFace[bFaceI], sourceProc[bFaceI]);
00822 
00823             HashTable<label, labelPair, labelPair::Hash<> >::const_iterator
00824                 iter = map.find(myData);
00825 
00826             if (iter != map.end())
00827             {
00828                 label nbrBFaceI = iter();
00829 
00830                 masterCoupledFaces[coupledI] = mesh.nInternalFaces() + bFaceI;
00831                 slaveCoupledFaces[coupledI] =
00832                     domainMesh.nInternalFaces()
00833                   + nbrBFaceI;
00834 
00835                 coupledI++;
00836             }
00837         }
00838     }
00839 
00840     masterCoupledFaces.setSize(coupledI);
00841     slaveCoupledFaces.setSize(coupledI);
00842 
00843     if (debug)
00844     {
00845         Pout<< "findCouples : found " << coupledI
00846             << " faces that will be stitched" << nl << endl;
00847     }
00848 }
00849 
00850 
00851 // Map data on boundary faces to new mesh (resulting from adding two meshes)
00852 Foam::labelList Foam::fvMeshDistribute::mapBoundaryData
00853 (
00854     const primitiveMesh& mesh,      // mesh after adding
00855     const mapAddedPolyMesh& map,
00856     const labelList& boundaryData0, // mesh before adding
00857     const label nInternalFaces1,
00858     const labelList& boundaryData1  // added mesh
00859 )
00860 {
00861     labelList newBoundaryData(mesh.nFaces() - mesh.nInternalFaces());
00862 
00863     forAll(boundaryData0, oldBFaceI)
00864     {
00865         label newFaceI = map.oldFaceMap()[oldBFaceI + map.nOldInternalFaces()];
00866 
00867         // Face still exists (is necessary?) and still boundary face
00868         if (newFaceI >= 0 && newFaceI >= mesh.nInternalFaces())
00869         {
00870             newBoundaryData[newFaceI - mesh.nInternalFaces()] =
00871                 boundaryData0[oldBFaceI];
00872         }
00873     }
00874 
00875     forAll(boundaryData1, addedBFaceI)
00876     {
00877         label newFaceI = map.addedFaceMap()[addedBFaceI + nInternalFaces1];
00878 
00879         if (newFaceI >= 0 && newFaceI >= mesh.nInternalFaces())
00880         {
00881             newBoundaryData[newFaceI - mesh.nInternalFaces()] =
00882                 boundaryData1[addedBFaceI];
00883         }
00884     }
00885 
00886     return newBoundaryData;
00887 }
00888 
00889 
00890 // Remove cells. Add all exposed faces to patch oldInternalPatchI
00891 Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::doRemoveCells
00892 (
00893     const labelList& cellsToRemove,
00894     const label oldInternalPatchI
00895 )
00896 {
00897     // Mesh change engine
00898     polyTopoChange meshMod(mesh_);
00899 
00900     // Cell removal topo engine. Do NOT synchronize parallel since
00901     // we are doing a local cell removal.
00902     removeCells cellRemover(mesh_, false);
00903 
00904     // Get all exposed faces
00905     labelList exposedFaces(cellRemover.getExposedFaces(cellsToRemove));
00906 
00907     // Insert the topo changes
00908     cellRemover.setRefinement
00909     (
00910         cellsToRemove,
00911         exposedFaces,
00912         labelList(exposedFaces.size(), oldInternalPatchI),  // patch for exposed
00913                                                             // faces.
00914         meshMod
00915     );
00916 
00917     // Change the mesh. No inflation. Note: no parallel comms allowed.
00918     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, false);
00919 
00920     // Update fields
00921     mesh_.updateMesh(map);
00922 
00923     // Move mesh (since morphing does not do this)
00924     if (map().hasMotionPoints())
00925     {
00926         mesh_.movePoints(map().preMotionPoints());
00927     }
00928 
00929     return map;
00930 }
00931 
00932 
00933 // Delete and add processor patches. Changes mesh and returns per neighbour proc
00934 // the processor patchID.
00935 void Foam::fvMeshDistribute::addProcPatches
00936 (
00937     const labelList& neighbourNewProc,   // processor that neighbour is on
00938     labelList& procPatchID
00939 )
00940 {
00941     // Now use the neighbourFace/Proc to repatch the mesh. These two lists
00942     // contain for all current boundary faces the global patchID (for non-proc
00943     // patch) or the processor.
00944 
00945     labelList procPatchSizes(Pstream::nProcs(), 0);
00946 
00947     forAll(neighbourNewProc, bFaceI)
00948     {
00949         if (neighbourNewProc[bFaceI] != -1)
00950         {
00951             procPatchSizes[neighbourNewProc[bFaceI]]++;
00952         }
00953     }
00954 
00955     // Per neighbour processor the label of the processor patch
00956     procPatchID.setSize(Pstream::nProcs());
00957 
00958     forAll(procPatchSizes, procI)
00959     {
00960         if (procPatchSizes[procI] > 0)
00961         {
00962             const word patchName =
00963                 "procBoundary"
00964               + name(Pstream::myProcNo())
00965               + "to"
00966               + name(procI);
00967 
00968 
00969             procPatchID[procI] = addProcPatch(patchName, procI);
00970             addPatchFields<volScalarField>
00971             (
00972                 processorFvPatchField<scalar>::typeName
00973             );
00974             addPatchFields<volVectorField>
00975             (
00976                 processorFvPatchField<vector>::typeName
00977             );
00978             addPatchFields<volSphericalTensorField>
00979             (
00980                 processorFvPatchField<sphericalTensor>::typeName
00981             );
00982             addPatchFields<volSymmTensorField>
00983             (
00984                 processorFvPatchField<symmTensor>::typeName
00985             );
00986             addPatchFields<volTensorField>
00987             (
00988                 processorFvPatchField<tensor>::typeName
00989             );
00990 
00991             addPatchFields<surfaceScalarField>
00992             (
00993                 processorFvPatchField<scalar>::typeName
00994             );
00995             addPatchFields<surfaceVectorField>
00996             (
00997                 processorFvPatchField<vector>::typeName
00998             );
00999             addPatchFields<surfaceSphericalTensorField>
01000             (
01001                 processorFvPatchField<sphericalTensor>::typeName
01002             );
01003             addPatchFields<surfaceSymmTensorField>
01004             (
01005                 processorFvPatchField<symmTensor>::typeName
01006             );
01007             addPatchFields<surfaceTensorField>
01008             (
01009                 processorFvPatchField<tensor>::typeName
01010             );
01011         }
01012         else
01013         {
01014             procPatchID[procI] = -1;
01015         }
01016     }
01017 }
01018 
01019 
01020 // Get boundary faces to be repatched. Is -1 or new patchID
01021 Foam::labelList Foam::fvMeshDistribute::getProcBoundaryPatch
01022 (
01023     const labelList& neighbourNewProc,  // new processor per boundary face
01024     const labelList& procPatchID        // patchID
01025 )
01026 {
01027     labelList patchIDs(neighbourNewProc);
01028 
01029     forAll(neighbourNewProc, bFaceI)
01030     {
01031         if (neighbourNewProc[bFaceI] != -1)
01032         {
01033             label nbrProc = neighbourNewProc[bFaceI];
01034 
01035             patchIDs[bFaceI] = procPatchID[nbrProc];
01036         }
01037         else
01038         {
01039             patchIDs[bFaceI] = -1;
01040         }
01041     }
01042     return patchIDs;
01043 }
01044 
01045 
01046 // Send mesh and coupling data.
01047 void Foam::fvMeshDistribute::sendMesh
01048 (
01049     const label domain,
01050     const fvMesh& mesh,
01051 
01052     const wordList& pointZoneNames,
01053     const wordList& faceZoneNames,
01054     const wordList& cellZoneNames,
01055 
01056     const labelList& sourceFace,
01057     const labelList& sourceProc,
01058     const labelList& sourceNewProc,
01059     OSstream& toDomain
01060 )
01061 {
01062     if (debug)
01063     {
01064         Pout<< "Sending to domain " << domain << nl
01065             << "    nPoints:" << mesh.nPoints() << nl
01066             << "    nFaces:" << mesh.nFaces() << nl
01067             << "    nCells:" << mesh.nCells() << nl
01068             << "    nPatches:" << mesh.boundaryMesh().size() << nl
01069             << endl;
01070     }
01071 
01072     // Assume sparse, possibly overlapping point zones. Get contents
01073     // in merged-zone indices.
01074     CompactListList_dev<label> zonePoints;
01075     {
01076         const pointZoneMesh& pointZones = mesh.pointZones();
01077 
01078         labelList rowSizes(pointZoneNames.size(), 0);
01079 
01080         forAll(pointZoneNames, nameI)
01081         {
01082             label myZoneID = pointZones.findZoneID(pointZoneNames[nameI]);
01083 
01084             if (myZoneID != -1)
01085             {
01086                 rowSizes[nameI] = pointZones[myZoneID].size();
01087             }
01088         }
01089         zonePoints.setSize(rowSizes);
01090 
01091         forAll(pointZoneNames, nameI)
01092         {
01093             label myZoneID = pointZones.findZoneID(pointZoneNames[nameI]);
01094 
01095             if (myZoneID != -1)
01096             {
01097                 zonePoints[nameI].assign(pointZones[myZoneID]);
01098             }
01099         }
01100     }
01101 
01102     // Assume sparse, possibly overlapping face zones
01103     CompactListList_dev<label> zoneFaces;
01104     CompactListList_dev<bool> zoneFaceFlip;
01105     {
01106         const faceZoneMesh& faceZones = mesh.faceZones();
01107 
01108         labelList rowSizes(faceZoneNames.size(), 0);
01109 
01110         forAll(faceZoneNames, nameI)
01111         {
01112             label myZoneID = faceZones.findZoneID(faceZoneNames[nameI]);
01113 
01114             if (myZoneID != -1)
01115             {
01116                 rowSizes[nameI] = faceZones[myZoneID].size();
01117             }
01118         }
01119 
01120         zoneFaces.setSize(rowSizes);
01121         zoneFaceFlip.setSize(rowSizes);
01122 
01123         forAll(faceZoneNames, nameI)
01124         {
01125             label myZoneID = faceZones.findZoneID(faceZoneNames[nameI]);
01126 
01127             if (myZoneID != -1)
01128             {
01129                 zoneFaces[nameI].assign(faceZones[myZoneID]);
01130                 zoneFaceFlip[nameI].assign(faceZones[myZoneID].flipMap());
01131             }
01132         }
01133     }
01134 
01135     // Assume sparse, possibly overlapping cell zones
01136     CompactListList_dev<label> zoneCells;
01137     {
01138         const cellZoneMesh& cellZones = mesh.cellZones();
01139 
01140         labelList rowSizes(cellZoneNames.size(), 0);
01141 
01142         forAll(cellZoneNames, nameI)
01143         {
01144             label myZoneID = cellZones.findZoneID(cellZoneNames[nameI]);
01145 
01146             if (myZoneID != -1)
01147             {
01148                 rowSizes[nameI] = cellZones[myZoneID].size();
01149             }
01150         }
01151 
01152         zoneCells.setSize(rowSizes);
01153 
01154         forAll(cellZoneNames, nameI)
01155         {
01156             label myZoneID = cellZones.findZoneID(cellZoneNames[nameI]);
01157 
01158             if (myZoneID != -1)
01159             {
01160                 zoneCells[nameI].assign(cellZones[myZoneID]);
01161             }
01162         }
01163     }
01165     //labelList cellZoneID;
01166     //if (hasCellZones)
01167     //{
01168     //    cellZoneID.setSize(mesh.nCells());
01169     //    cellZoneID = -1;
01170     //
01171     //    const cellZoneMesh& cellZones = mesh.cellZones();
01172     //
01173     //    forAll(cellZones, zoneI)
01174     //    {
01175     //        UIndirectList<label>(cellZoneID, cellZones[zoneI]) = zoneI;
01176     //    }
01177     //}
01178 
01179     // Send
01180     toDomain
01181         << mesh.points()
01182         << CompactListList_dev<label, face>(mesh.faces())
01183         << mesh.faceOwner()
01184         << mesh.faceNeighbour()
01185         << mesh.boundaryMesh()
01186 
01187         << zonePoints
01188         << zoneFaces
01189         << zoneFaceFlip
01190         << zoneCells
01191 
01192         << sourceFace
01193         << sourceProc
01194         << sourceNewProc;
01195 
01196 
01197     if (debug)
01198     {
01199         Pout<< "Started sending mesh to domain " << domain
01200             << endl;
01201     }
01202 }
01203 
01204 
01205 // Receive mesh. Opposite of sendMesh
01206 Foam::autoPtr<Foam::fvMesh> Foam::fvMeshDistribute::receiveMesh
01207 (
01208     const label domain,
01209     const wordList& pointZoneNames,
01210     const wordList& faceZoneNames,
01211     const wordList& cellZoneNames,
01212     const Time& runTime,
01213     labelList& domainSourceFace,
01214     labelList& domainSourceProc,
01215     labelList& domainSourceNewProc,
01216     ISstream& fromNbr
01217 )
01218 {
01219     pointField domainPoints(fromNbr);
01220     faceList domainFaces = CompactListList_dev<label, face>(fromNbr)();
01221     labelList domainAllOwner(fromNbr);
01222     labelList domainAllNeighbour(fromNbr);
01223     PtrList<entry> patchEntries(fromNbr);
01224 
01225     CompactListList_dev<label> zonePoints(fromNbr);
01226     CompactListList_dev<label> zoneFaces(fromNbr);
01227     CompactListList_dev<bool> zoneFaceFlip(fromNbr);
01228     CompactListList_dev<label> zoneCells(fromNbr);
01229 
01230     fromNbr
01231         >> domainSourceFace
01232         >> domainSourceProc
01233         >> domainSourceNewProc;
01234 
01235     // Construct fvMesh
01236     autoPtr<fvMesh> domainMeshPtr
01237     (
01238         new fvMesh
01239         (
01240             IOobject
01241             (
01242                 fvMesh::defaultRegion,
01243                 runTime.timeName(),
01244                 runTime,
01245                 IOobject::NO_READ
01246             ),
01247             xferMove(domainPoints),
01248             xferMove(domainFaces),
01249             xferMove(domainAllOwner),
01250             xferMove(domainAllNeighbour),
01251             false                   // no parallel comms
01252         )
01253     );
01254     fvMesh& domainMesh = domainMeshPtr();
01255 
01256     List<polyPatch*> patches(patchEntries.size());
01257 
01258     forAll(patchEntries, patchI)
01259     {
01260         patches[patchI] = polyPatch::New
01261         (
01262             patchEntries[patchI].keyword(),
01263             patchEntries[patchI].dict(),
01264             patchI,
01265             domainMesh.boundaryMesh()
01266         ).ptr();
01267     }
01268     // Add patches; no parallel comms
01269     domainMesh.addFvPatches(patches, false);
01270 
01271     // Construct zones
01272     List<pointZone*> pZonePtrs(pointZoneNames.size());
01273     forAll(pZonePtrs, i)
01274     {
01275         pZonePtrs[i] = new pointZone
01276         (
01277             pointZoneNames[i],
01278             zonePoints[i],
01279             i,
01280             domainMesh.pointZones()
01281         );
01282     }
01283 
01284     List<faceZone*> fZonePtrs(faceZoneNames.size());
01285     forAll(fZonePtrs, i)
01286     {
01287         fZonePtrs[i] = new faceZone
01288         (
01289             faceZoneNames[i],
01290             zoneFaces[i],
01291             zoneFaceFlip[i],
01292             i,
01293             domainMesh.faceZones()
01294         );
01295     }
01296 
01297     List<cellZone*> cZonePtrs(cellZoneNames.size());
01298     forAll(cZonePtrs, i)
01299     {
01300         cZonePtrs[i] = new cellZone
01301         (
01302             cellZoneNames[i],
01303             zoneCells[i],
01304             i,
01305             domainMesh.cellZones()
01306         );
01307     }
01308     domainMesh.addZones(pZonePtrs, fZonePtrs, cZonePtrs);
01309 
01310     return domainMeshPtr;
01311 }
01312 
01313 
01314 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
01315 
01316 // Construct from components
01317 Foam::fvMeshDistribute::fvMeshDistribute(fvMesh& mesh, const scalar mergeTol)
01318 :
01319     mesh_(mesh),
01320     mergeTol_(mergeTol)
01321 {}
01322 
01323 
01324 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
01325 
01326 Foam::labelList Foam::fvMeshDistribute::countCells
01327 (
01328     const labelList& distribution
01329 )
01330 {
01331     labelList nCells(Pstream::nProcs(), 0);
01332     forAll(distribution, cellI)
01333     {
01334         label newProc = distribution[cellI];
01335 
01336         if (newProc < 0 || newProc >= Pstream::nProcs())
01337         {
01338             FatalErrorIn("fvMeshDistribute::distribute(const labelList&)")
01339                 << "Distribution should be in range 0.." << Pstream::nProcs()-1
01340                 << endl
01341                 << "At index " << cellI << " distribution:" << newProc
01342                 << abort(FatalError);
01343         }
01344         nCells[newProc]++;
01345     }
01346     return nCells;
01347 }
01348 
01349 
01350 Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
01351 (
01352     const labelList& distribution
01353 )
01354 {
01355     // Some checks on distribution
01356     if (distribution.size() != mesh_.nCells())
01357     {
01358         FatalErrorIn("fvMeshDistribute::distribute(const labelList&)")
01359             << "Size of distribution:"
01360             << distribution.size() << " mesh nCells:" << mesh_.nCells()
01361             << abort(FatalError);
01362     }
01363 
01364 
01365     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
01366 
01367     // Check all processors have same non-proc patches in same order.
01368     if (patches.checkParallelSync(true))
01369     {
01370         FatalErrorIn("fvMeshDistribute::distribute(const labelList&)")
01371             << "This application requires all non-processor patches"
01372             << " to be present in the same order on all patches" << nl
01373             << "followed by the processor patches (which of course are unique)."
01374             << nl
01375             << "Local patches:" << mesh_.boundaryMesh().names()
01376             << abort(FatalError);
01377     }
01378 
01379     // Save some data for mapping later on
01380     const label nOldPoints(mesh_.nPoints());
01381     const label nOldFaces(mesh_.nFaces());
01382     const label nOldCells(mesh_.nCells());
01383     labelList oldPatchStarts(patches.size());
01384     labelList oldPatchNMeshPoints(patches.size());
01385     forAll(patches, patchI)
01386     {
01387         oldPatchStarts[patchI] = patches[patchI].start();
01388         oldPatchNMeshPoints[patchI] = patches[patchI].nPoints();
01389     }
01390 
01391 
01392 
01393     // Short circuit trivial case.
01394     if (!Pstream::parRun())
01395     {
01396         // Collect all maps and return
01397         return autoPtr<mapDistributePolyMesh>
01398         (
01399             new mapDistributePolyMesh
01400             (
01401                 mesh_,
01402 
01403                 nOldPoints,
01404                 nOldFaces,
01405                 nOldCells,
01406                 oldPatchStarts.xfer(),
01407                 oldPatchNMeshPoints.xfer(),
01408 
01409                 labelListList(1, identity(mesh_.nPoints())).xfer(),//subPointMap
01410                 labelListList(1, identity(mesh_.nFaces())).xfer(), //subFaceMap
01411                 labelListList(1, identity(mesh_.nCells())).xfer(), //subCellMap
01412                 labelListList(1, identity(patches.size())).xfer(), //subPatchMap
01413 
01414                 labelListList(1, identity(mesh_.nPoints())).xfer(),//pointMap
01415                 labelListList(1, identity(mesh_.nFaces())).xfer(), //faceMap
01416                 labelListList(1, identity(mesh_.nCells())).xfer(), //cellMap
01417                 labelListList(1, identity(patches.size())).xfer()  //patchMap
01418             )
01419         );
01420     }
01421 
01422 
01423     // Collect any zone names
01424     const wordList pointZoneNames(mergeWordList(mesh_.pointZones().names()));
01425     const wordList faceZoneNames(mergeWordList(mesh_.faceZones().names()));
01426     const wordList cellZoneNames(mergeWordList(mesh_.cellZones().names()));
01427 
01428 
01429 
01430     // Local environment of all boundary faces
01431     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
01432 
01433     // A face is uniquely defined by
01434     //  - proc
01435     //  - local face no
01436     //
01437     // To glue the parts of meshes which can get sent from anywhere we
01438     // need to know on boundary faces what the above tuple on both sides is.
01439     // So we need to maintain:
01440     //  - original face
01441     //  - original processor id (= trivial)
01442     // For coupled boundaries (where the faces are 'duplicate') we take the
01443     // lowest numbered processor as the data to store.
01444     //
01445     // Additionally to create the procboundaries we need to know where the owner
01446     // cell on the other side now is: newNeighbourProc.
01447     //
01448 
01449     // physical boundary:
01450     //     sourceProc = -1
01451     //     sourceNewProc = -1
01452     //     sourceFace = patchID
01453     // coupled boundary:
01454     //     sourceProc = proc
01455     //     sourceNewProc = distribution[cell on proc]
01456     //     sourceFace = face
01457     labelList sourceFace;
01458     labelList sourceProc;
01459     labelList sourceNewProc;
01460     getNeighbourData(distribution, sourceFace, sourceProc, sourceNewProc);
01461 
01462 
01463     // Remove meshPhi. Since this would otherwise dissappear anyway
01464     // during topo changes and we have to guarantee that all the fields
01465     // can be sent.
01466     mesh_.clearOut();
01467     mesh_.resetMotion();
01468 
01469     // Get data to send. Make sure is synchronised
01470     const wordList volScalars(mesh_.names(volScalarField::typeName));
01471     checkEqualWordList("volScalarFields", volScalars);
01472     const wordList volVectors(mesh_.names(volVectorField::typeName));
01473     checkEqualWordList("volVectorFields", volVectors);
01474     const wordList volSphereTensors
01475     (
01476         mesh_.names(volSphericalTensorField::typeName)
01477     );
01478     checkEqualWordList("volSphericalTensorFields", volSphereTensors);
01479     const wordList volSymmTensors(mesh_.names(volSymmTensorField::typeName));
01480     checkEqualWordList("volSymmTensorFields", volSymmTensors);
01481     const wordList volTensors(mesh_.names(volTensorField::typeName));
01482     checkEqualWordList("volTensorField", volTensors);
01483 
01484     const wordList surfScalars(mesh_.names(surfaceScalarField::typeName));
01485     checkEqualWordList("surfaceScalarFields", surfScalars);
01486     const wordList surfVectors(mesh_.names(surfaceVectorField::typeName));
01487     checkEqualWordList("surfaceVectorFields", surfVectors);
01488     const wordList surfSphereTensors
01489     (
01490         mesh_.names(surfaceSphericalTensorField::typeName)
01491     );
01492     checkEqualWordList("surfaceSphericalTensorFields", surfSphereTensors);
01493     const wordList surfSymmTensors
01494     (
01495         mesh_.names(surfaceSymmTensorField::typeName)
01496     );
01497     checkEqualWordList("surfaceSymmTensorFields", surfSymmTensors);
01498     const wordList surfTensors(mesh_.names(surfaceTensorField::typeName));
01499     checkEqualWordList("surfaceTensorFields", surfTensors);
01500 
01501 
01502 
01503 
01504     // Find patch to temporarily put exposed and processor faces into.
01505     label oldInternalPatchI = findNonEmptyPatch();
01506 
01507 
01508 
01509     // Delete processor patches, starting from the back. Move all faces into
01510     // oldInternalPatchI.
01511     labelList repatchFaceMap;
01512     {
01513         autoPtr<mapPolyMesh> repatchMap = deleteProcPatches(oldInternalPatchI);
01514 
01515         // Store face map (only face ordering that changed)
01516         repatchFaceMap = repatchMap().faceMap();
01517 
01518         // Reorder all boundary face data (sourceProc, sourceFace etc.)
01519         labelList bFaceMap
01520         (
01521             SubList<label>
01522             (
01523                 repatchMap().reverseFaceMap(),
01524                 mesh_.nFaces() - mesh_.nInternalFaces(),
01525                 mesh_.nInternalFaces()
01526             )
01527           - mesh_.nInternalFaces()
01528         );
01529 
01530         inplaceReorder(bFaceMap, sourceFace);
01531         inplaceReorder(bFaceMap, sourceProc);
01532         inplaceReorder(bFaceMap, sourceNewProc);
01533     }
01534 
01535 
01536 
01537     // Print a bit.
01538     if (debug)
01539     {
01540         Pout<< nl << "MESH WITH PROC PATCHES DELETED:" << endl;
01541         printMeshInfo(mesh_);
01542         printFieldInfo<volScalarField>(mesh_);
01543         printFieldInfo<volVectorField>(mesh_);
01544         printFieldInfo<volSphericalTensorField>(mesh_);
01545         printFieldInfo<volSymmTensorField>(mesh_);
01546         printFieldInfo<volTensorField>(mesh_);
01547         printFieldInfo<surfaceScalarField>(mesh_);
01548         printFieldInfo<surfaceVectorField>(mesh_);
01549         printFieldInfo<surfaceSphericalTensorField>(mesh_);
01550         printFieldInfo<surfaceSymmTensorField>(mesh_);
01551         printFieldInfo<surfaceTensorField>(mesh_);
01552         Pout<< nl << endl;
01553     }
01554 
01555 
01556 
01557     // Maps from subsetted mesh (that is sent) back to original maps
01558     labelListList subCellMap(Pstream::nProcs());
01559     labelListList subFaceMap(Pstream::nProcs());
01560     labelListList subPointMap(Pstream::nProcs());
01561     labelListList subPatchMap(Pstream::nProcs());
01562     // Maps from subsetted mesh to reconstructed mesh
01563     labelListList constructCellMap(Pstream::nProcs());
01564     labelListList constructFaceMap(Pstream::nProcs());
01565     labelListList constructPointMap(Pstream::nProcs());
01566     labelListList constructPatchMap(Pstream::nProcs());
01567 
01568 
01569 
01570 
01571     // Find out schedule
01572     // ~~~~~~~~~~~~~~~~~
01573 
01574     labelListList nSendCells(Pstream::nProcs());
01575     nSendCells[Pstream::myProcNo()] = countCells(distribution);
01576     Pstream::gatherList(nSendCells);
01577     Pstream::scatterList(nSendCells);
01578 
01579 
01580     // Allocate buffers
01581     PtrList<OStringStream> sendStr(Pstream::nProcs());
01582     forAll(sendStr, i)
01583     {
01584         sendStr.set(i, new OStringStream(IOstream::BINARY));
01585     }
01586     //PstreamBuffers pBufs(Pstream::nonBlocking);
01587 
01588 
01589     // What to send to neighbouring domains
01590     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
01591 
01592     forAll(nSendCells[Pstream::myProcNo()], recvProc)
01593     {
01594         if
01595         (
01596             recvProc != Pstream::myProcNo()
01597          && nSendCells[Pstream::myProcNo()][recvProc] > 0
01598         )
01599         {
01600             // Send to recvProc
01601 
01602             if (debug)
01603             {
01604                 Pout<< nl
01605                     << "SUBSETTING FOR DOMAIN " << recvProc
01606                     << " cells to send:"
01607                     << nSendCells[Pstream::myProcNo()][recvProc]
01608                     << nl << endl;
01609             }
01610 
01611             // Pstream for sending mesh and fields
01612             //OPstream str(Pstream::blocking, recvProc);
01613             //UOPstream str(recvProc, pBufs);
01614 
01615             // Mesh subsetting engine
01616             fvMeshSubset subsetter(mesh_);
01617 
01618             // Subset the cells of the current domain.
01619             subsetter.setLargeCellSubset
01620             (
01621                 distribution,
01622                 recvProc,
01623                 oldInternalPatchI,  // oldInternalFaces patch
01624                 false               // no parallel sync
01625             );
01626 
01627             subCellMap[recvProc] = subsetter.cellMap();
01628             subFaceMap[recvProc] = renumber
01629             (
01630                 repatchFaceMap,
01631                 subsetter.faceMap()
01632             );
01633             subPointMap[recvProc] = subsetter.pointMap();
01634             subPatchMap[recvProc] = subsetter.patchMap();
01635 
01636 
01637             // Subset the boundary fields (owner/neighbour/processor)
01638             labelList procSourceFace;
01639             labelList procSourceProc;
01640             labelList procSourceNewProc;
01641 
01642             subsetBoundaryData
01643             (
01644                 subsetter.subMesh(),
01645                 subsetter.faceMap(),        // from subMesh to mesh
01646                 subsetter.cellMap(),        //      ,,      ,,
01647 
01648                 distribution,               // old mesh distribution
01649                 mesh_.faceOwner(),          // old owner
01650                 mesh_.faceNeighbour(),
01651                 mesh_.nInternalFaces(),
01652 
01653                 sourceFace,
01654                 sourceProc,
01655                 sourceNewProc,
01656 
01657                 procSourceFace,
01658                 procSourceProc,
01659                 procSourceNewProc
01660             );
01661 
01662 
01663 
01664             // Send to neighbour
01665             sendMesh
01666             (
01667                 recvProc,
01668                 subsetter.subMesh(),
01669 
01670                 pointZoneNames,
01671                 faceZoneNames,
01672                 cellZoneNames,
01673 
01674                 procSourceFace,
01675                 procSourceProc,
01676                 procSourceNewProc,
01677                 sendStr[recvProc]
01678             );
01679             sendFields<volScalarField>
01680             (
01681                 recvProc,
01682                 volScalars,
01683                 subsetter,
01684                 sendStr[recvProc]
01685             );
01686             sendFields<volVectorField>
01687             (
01688                 recvProc,
01689                 volVectors,
01690                 subsetter,
01691                 sendStr[recvProc]
01692             );
01693             sendFields<volSphericalTensorField>
01694             (
01695                 recvProc,
01696                 volSphereTensors,
01697                 subsetter,
01698                 sendStr[recvProc]
01699             );
01700             sendFields<volSymmTensorField>
01701             (
01702                 recvProc,
01703                 volSymmTensors,
01704                 subsetter,
01705                 sendStr[recvProc]
01706             );
01707             sendFields<volTensorField>
01708             (
01709                 recvProc,
01710                 volTensors,
01711                 subsetter,
01712                 sendStr[recvProc]
01713             );
01714 
01715             sendFields<surfaceScalarField>
01716             (
01717                 recvProc,
01718                 surfScalars,
01719                 subsetter,
01720                 sendStr[recvProc]
01721             );
01722             sendFields<surfaceVectorField>
01723             (
01724                 recvProc,
01725                 surfVectors,
01726                 subsetter,
01727                 sendStr[recvProc]
01728             );
01729             sendFields<surfaceSphericalTensorField>
01730             (
01731                 recvProc,
01732                 surfSphereTensors,
01733                 subsetter,
01734                 sendStr[recvProc]
01735             );
01736             sendFields<surfaceSymmTensorField>
01737             (
01738                 recvProc,
01739                 surfSymmTensors,
01740                 subsetter,
01741                 sendStr[recvProc]
01742             );
01743             sendFields<surfaceTensorField>
01744             (
01745                 recvProc,
01746                 surfTensors,
01747                 subsetter,
01748                 sendStr[recvProc]
01749             );
01750         }
01751     }
01752 
01753 
01754     // Start sending&receiving from buffers
01755     //pBufs.finishedSends();
01756 
01757     // get the data.
01758     PtrList<IStringStream> recvStr(Pstream::nProcs());
01759     {
01760         List<List<char> > sendBufs(sendStr.size());
01761         forAll(sendStr, procI)
01762         {
01763             string contents = sendStr[procI].str();
01764             const char* ptr = contents.data();
01765 
01766             sendBufs[procI].setSize(contents.size());
01767             forAll(sendBufs[procI], i)
01768             {
01769                 sendBufs[procI][i] = *ptr++;
01770             }
01771             // Clear OStringStream
01772             sendStr.set(procI, NULL);
01773         }
01774 
01775         // Transfer sendBufs into recvBufs
01776         List<List<char> > recvBufs(Pstream::nProcs());
01777         labelListList sizes(Pstream::nProcs());
01778         exchange<List<char>, char>(sendBufs, recvBufs, sizes);
01779 
01780         forAll(recvStr, procI)
01781         {
01782             string contents(recvBufs[procI].begin(), recvBufs[procI].size());
01783             recvStr.set
01784             (
01785                 procI,
01786                 new IStringStream(contents, IOstream::BINARY)
01787             );
01788         }
01789     }
01790 
01791 
01792     // Subset the part that stays
01793     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
01794 
01795     {
01796         // Save old mesh maps before changing mesh
01797         const labelList oldFaceOwner(mesh_.faceOwner());
01798         const labelList oldFaceNeighbour(mesh_.faceNeighbour());
01799         const label oldInternalFaces = mesh_.nInternalFaces();
01800 
01801         // Remove cells.
01802         autoPtr<mapPolyMesh> subMap
01803         (
01804             doRemoveCells
01805             (
01806                 select(false, distribution, Pstream::myProcNo()),
01807                 oldInternalPatchI
01808             )
01809         );
01810 
01811         // Addressing from subsetted mesh
01812         subCellMap[Pstream::myProcNo()] = subMap().cellMap();
01813         subFaceMap[Pstream::myProcNo()] = renumber
01814         (
01815             repatchFaceMap,
01816             subMap().faceMap()
01817         );
01818         subPointMap[Pstream::myProcNo()] = subMap().pointMap();
01819         subPatchMap[Pstream::myProcNo()] = identity(patches.size());
01820 
01821         // Initialize all addressing into current mesh
01822         constructCellMap[Pstream::myProcNo()] = identity(mesh_.nCells());
01823         constructFaceMap[Pstream::myProcNo()] = identity(mesh_.nFaces());
01824         constructPointMap[Pstream::myProcNo()] = identity(mesh_.nPoints());
01825         constructPatchMap[Pstream::myProcNo()] = identity(patches.size());
01826 
01827         // Subset the mesh data: neighbourCell/neighbourProc
01828         // fields
01829         labelList domainSourceFace;
01830         labelList domainSourceProc;
01831         labelList domainSourceNewProc;
01832 
01833         subsetBoundaryData
01834         (
01835             mesh_,                          // new mesh
01836             subMap().faceMap(),             // from new to original mesh
01837             subMap().cellMap(),
01838 
01839             distribution,                   // distribution before subsetting
01840             oldFaceOwner,                   // owner before subsetting
01841             oldFaceNeighbour,               // neighbour        ,,
01842             oldInternalFaces,               // nInternalFaces   ,,
01843 
01844             sourceFace,
01845             sourceProc,
01846             sourceNewProc,
01847 
01848             domainSourceFace,
01849             domainSourceProc,
01850             domainSourceNewProc
01851         );
01852 
01853         sourceFace.transfer(domainSourceFace);
01854         sourceProc.transfer(domainSourceProc);
01855         sourceNewProc.transfer(domainSourceNewProc);
01856     }
01857 
01858 
01859     // Print a bit.
01860     if (debug)
01861     {
01862         Pout<< nl << "STARTING MESH:" << endl;
01863         printMeshInfo(mesh_);
01864         printFieldInfo<volScalarField>(mesh_);
01865         printFieldInfo<volVectorField>(mesh_);
01866         printFieldInfo<volSphericalTensorField>(mesh_);
01867         printFieldInfo<volSymmTensorField>(mesh_);
01868         printFieldInfo<volTensorField>(mesh_);
01869         printFieldInfo<surfaceScalarField>(mesh_);
01870         printFieldInfo<surfaceVectorField>(mesh_);
01871         printFieldInfo<surfaceSphericalTensorField>(mesh_);
01872         printFieldInfo<surfaceSymmTensorField>(mesh_);
01873         printFieldInfo<surfaceTensorField>(mesh_);
01874         Pout<< nl << endl;
01875     }
01876 
01877 
01878 
01879     // Receive and add what was sent
01880     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
01881 
01882     forAll(nSendCells, sendProc)
01883     {
01884         // Did processor sendProc send anything to me?
01885         if
01886         (
01887             sendProc != Pstream::myProcNo()
01888          && nSendCells[sendProc][Pstream::myProcNo()] > 0
01889         )
01890         {
01891             if (debug)
01892             {
01893                 Pout<< nl
01894                     << "RECEIVING FROM DOMAIN " << sendProc
01895                     << " cells to receive:"
01896                     << nSendCells[sendProc][Pstream::myProcNo()]
01897                     << nl << endl;
01898             }
01899 
01900 
01901             // Pstream for receiving mesh and fields
01902             //UIPstream str(sendProc, pBufs);
01903 
01904 
01905             // Receive from sendProc
01906             labelList domainSourceFace;
01907             labelList domainSourceProc;
01908             labelList domainSourceNewProc;
01909 
01910             autoPtr<fvMesh> domainMeshPtr;
01911             PtrList<volScalarField> vsf;
01912             PtrList<volVectorField> vvf;
01913             PtrList<volSphericalTensorField> vsptf;
01914             PtrList<volSymmTensorField> vsytf;
01915             PtrList<volTensorField> vtf;
01916             PtrList<surfaceScalarField> ssf;
01917             PtrList<surfaceVectorField> svf;
01918             PtrList<surfaceSphericalTensorField> ssptf;
01919             PtrList<surfaceSymmTensorField> ssytf;
01920             PtrList<surfaceTensorField> stf;
01921 
01922             // Opposite of sendMesh
01923             {
01924                 domainMeshPtr = receiveMesh
01925                 (
01926                     sendProc,
01927                     pointZoneNames,
01928                     faceZoneNames,
01929                     cellZoneNames,
01930 
01931                     const_cast<Time&>(mesh_.time()),
01932                     domainSourceFace,
01933                     domainSourceProc,
01934                     domainSourceNewProc,
01935                     recvStr[sendProc]
01936                 );
01937                 fvMesh& domainMesh = domainMeshPtr();
01938 
01939                 // Receive fields. Read as single dictionary because
01940                 // of problems reading consecutive fields from single stream.
01941                 dictionary fieldDicts(recvStr[sendProc]);
01942 
01943                 receiveFields<volScalarField>
01944                 (
01945                     sendProc,
01946                     volScalars,
01947                     domainMesh,
01948                     vsf,
01949                     fieldDicts.subDict(volScalarField::typeName)
01950                 );
01951                 receiveFields<volVectorField>
01952                 (
01953                     sendProc,
01954                     volVectors,
01955                     domainMesh,
01956                     vvf,
01957                     fieldDicts.subDict(volVectorField::typeName)
01958                 );
01959                 receiveFields<volSphericalTensorField>
01960                 (
01961                     sendProc,
01962                     volSphereTensors,
01963                     domainMesh,
01964                     vsptf,
01965                     fieldDicts.subDict(volSphericalTensorField::typeName)
01966                 );
01967                 receiveFields<volSymmTensorField>
01968                 (
01969                     sendProc,
01970                     volSymmTensors,
01971                     domainMesh,
01972                     vsytf,
01973                     fieldDicts.subDict(volSymmTensorField::typeName)
01974                 );
01975                 receiveFields<volTensorField>
01976                 (
01977                     sendProc,
01978                     volTensors,
01979                     domainMesh,
01980                     vtf,
01981                     fieldDicts.subDict(volTensorField::typeName)
01982                 );
01983 
01984                 receiveFields<surfaceScalarField>
01985                 (
01986                     sendProc,
01987                     surfScalars,
01988                     domainMesh,
01989                     ssf,
01990                     fieldDicts.subDict(surfaceScalarField::typeName)
01991                 );
01992                 receiveFields<surfaceVectorField>
01993                 (
01994                     sendProc,
01995                     surfVectors,
01996                     domainMesh,
01997                     svf,
01998                     fieldDicts.subDict(surfaceVectorField::typeName)
01999                 );
02000                 receiveFields<surfaceSphericalTensorField>
02001                 (
02002                     sendProc,
02003                     surfSphereTensors,
02004                     domainMesh,
02005                     ssptf,
02006                     fieldDicts.subDict(surfaceSphericalTensorField::typeName)
02007                 );
02008                 receiveFields<surfaceSymmTensorField>
02009                 (
02010                     sendProc,
02011                     surfSymmTensors,
02012                     domainMesh,
02013                     ssytf,
02014                     fieldDicts.subDict(surfaceSymmTensorField::typeName)
02015                 );
02016                 receiveFields<surfaceTensorField>
02017                 (
02018                     sendProc,
02019                     surfTensors,
02020                     domainMesh,
02021                     stf,
02022                     fieldDicts.subDict(surfaceTensorField::typeName)
02023                 );
02024             }
02025             const fvMesh& domainMesh = domainMeshPtr();
02026 
02027 
02028             constructCellMap[sendProc] = identity(domainMesh.nCells());
02029             constructFaceMap[sendProc] = identity(domainMesh.nFaces());
02030             constructPointMap[sendProc] = identity(domainMesh.nPoints());
02031             constructPatchMap[sendProc] =
02032                 identity(domainMesh.boundaryMesh().size());
02033 
02034 
02035             // Print a bit.
02036             if (debug)
02037             {
02038                 Pout<< nl << "RECEIVED MESH FROM:" << sendProc << endl;
02039                 printMeshInfo(domainMesh);
02040                 printFieldInfo<volScalarField>(domainMesh);
02041                 printFieldInfo<volVectorField>(domainMesh);
02042                 printFieldInfo<volSphericalTensorField>(domainMesh);
02043                 printFieldInfo<volSymmTensorField>(domainMesh);
02044                 printFieldInfo<volTensorField>(domainMesh);
02045                 printFieldInfo<surfaceScalarField>(domainMesh);
02046                 printFieldInfo<surfaceVectorField>(domainMesh);
02047                 printFieldInfo<surfaceSphericalTensorField>(domainMesh);
02048                 printFieldInfo<surfaceSymmTensorField>(domainMesh);
02049                 printFieldInfo<surfaceTensorField>(domainMesh);
02050             }
02051 
02052 
02053             // Now this mesh we received (from sendProc) needs to be merged
02054             // with the current mesh. On the current mesh we have for all
02055             // boundaryfaces the original face and processor. See if we can
02056             // match these up to the received domainSourceFace and
02057             // domainSourceProc.
02058             labelList masterCoupledFaces;
02059             labelList slaveCoupledFaces;
02060             findCouples
02061             (
02062                 mesh_,
02063 
02064                 sourceFace,
02065                 sourceProc,
02066 
02067                 sendProc,
02068                 domainMesh,
02069                 domainSourceFace,
02070                 domainSourceProc,
02071 
02072                 masterCoupledFaces,
02073                 slaveCoupledFaces
02074             );
02075 
02076             // Generate additional coupling info (points, edges) from
02077             // faces-that-match
02078             faceCoupleInfo couples
02079             (
02080                 mesh_,
02081                 masterCoupledFaces,
02082                 domainMesh,
02083                 slaveCoupledFaces,
02084                 mergeTol_,              // merge tolerance
02085                 true,                   // faces align
02086                 true,                   // couples are ordered already
02087                 false
02088             );
02089 
02090 
02091             // Add domainMesh to mesh
02092             // ~~~~~~~~~~~~~~~~~~~~~~
02093 
02094             autoPtr<mapAddedPolyMesh> map = fvMeshAdder::add
02095             (
02096                 mesh_,
02097                 domainMesh,
02098                 couples,
02099                 false           // no parallel comms
02100             );
02101 
02102             // Update mesh data: sourceFace,sourceProc for added
02103             // mesh.
02104 
02105             sourceFace =
02106                 mapBoundaryData
02107                 (
02108                     mesh_,
02109                     map(),
02110                     sourceFace,
02111                     domainMesh.nInternalFaces(),
02112                     domainSourceFace
02113                 );
02114             sourceProc =
02115                 mapBoundaryData
02116                 (
02117                     mesh_,
02118                     map(),
02119                     sourceProc,
02120                     domainMesh.nInternalFaces(),
02121                     domainSourceProc
02122                 );
02123             sourceNewProc =
02124                 mapBoundaryData
02125                 (
02126                     mesh_,
02127                     map(),
02128                     sourceNewProc,
02129                     domainMesh.nInternalFaces(),
02130                     domainSourceNewProc
02131                 );
02132 
02133             // Update all addressing so xxProcAddressing points to correct item
02134             // in masterMesh.
02135             const labelList& oldCellMap = map().oldCellMap();
02136             const labelList& oldFaceMap = map().oldFaceMap();
02137             const labelList& oldPointMap = map().oldPointMap();
02138             const labelList& oldPatchMap = map().oldPatchMap();
02139 
02140             forAll(constructPatchMap, procI)
02141             {
02142                 if (procI != sendProc && constructPatchMap[procI].size())
02143                 {
02144                     // Processor already in mesh (either myProcNo or received)
02145                     inplaceRenumber(oldCellMap, constructCellMap[procI]);
02146                     inplaceRenumber(oldFaceMap, constructFaceMap[procI]);
02147                     inplaceRenumber(oldPointMap, constructPointMap[procI]);
02148                     inplaceRenumber(oldPatchMap, constructPatchMap[procI]);
02149                 }
02150             }
02151 
02152             // Added processor
02153             inplaceRenumber(map().addedCellMap(), constructCellMap[sendProc]);
02154             inplaceRenumber(map().addedFaceMap(), constructFaceMap[sendProc]);
02155             inplaceRenumber(map().addedPointMap(), constructPointMap[sendProc]);
02156             inplaceRenumber(map().addedPatchMap(), constructPatchMap[sendProc]);
02157 
02158             if (debug)
02159             {
02160                 Pout<< nl << "MERGED MESH FROM:" << sendProc << endl;
02161                 printMeshInfo(mesh_);
02162                 printFieldInfo<volScalarField>(mesh_);
02163                 printFieldInfo<volVectorField>(mesh_);
02164                 printFieldInfo<volSphericalTensorField>(mesh_);
02165                 printFieldInfo<volSymmTensorField>(mesh_);
02166                 printFieldInfo<volTensorField>(mesh_);
02167                 printFieldInfo<surfaceScalarField>(mesh_);
02168                 printFieldInfo<surfaceVectorField>(mesh_);
02169                 printFieldInfo<surfaceSphericalTensorField>(mesh_);
02170                 printFieldInfo<surfaceSymmTensorField>(mesh_);
02171                 printFieldInfo<surfaceTensorField>(mesh_);
02172                 Pout<< nl << endl;
02173             }
02174         }
02175     }
02176 
02177 
02178     // Print a bit.
02179     if (debug)
02180     {
02181         Pout<< nl << "REDISTRIBUTED MESH:" << endl;
02182         printMeshInfo(mesh_);
02183         printFieldInfo<volScalarField>(mesh_);
02184         printFieldInfo<volVectorField>(mesh_);
02185         printFieldInfo<volSphericalTensorField>(mesh_);
02186         printFieldInfo<volSymmTensorField>(mesh_);
02187         printFieldInfo<volTensorField>(mesh_);
02188         printFieldInfo<surfaceScalarField>(mesh_);
02189         printFieldInfo<surfaceVectorField>(mesh_);
02190         printFieldInfo<surfaceSphericalTensorField>(mesh_);
02191         printFieldInfo<surfaceSymmTensorField>(mesh_);
02192         printFieldInfo<surfaceTensorField>(mesh_);
02193         Pout<< nl << endl;
02194     }
02195 
02196 
02197 
02198     // Add processorPatches
02199     // ~~~~~~~~~~~~~~~~~~~~
02200 
02201     // Per neighbour processor the patchID to it (or -1).
02202     labelList procPatchID;
02203 
02204     // Add processor patches.
02205     addProcPatches(sourceNewProc, procPatchID);
02206 
02207     // Put faces into correct patch. Note that we now have proper
02208     // processorPolyPatches again so repatching will take care of coupled face
02209     // ordering.
02210 
02211     // Get boundary faces to be repatched. Is -1 or new patchID
02212     labelList newPatchID
02213     (
02214         getProcBoundaryPatch
02215         (
02216             sourceNewProc,
02217             procPatchID
02218         )
02219     );
02220 
02221     // Change patches. Since this might change ordering of coupled faces
02222     // we also need to adapt our constructMaps.
02223     repatch(newPatchID, constructFaceMap);
02224 
02225     // See if any geometrically shared points need to be merged. Note: does
02226     // parallel comms.
02227     mergeSharedPoints(constructPointMap);
02228 
02229     // Bit of hack: processorFvPatchField does not get reset since created
02230     // from nothing so explicitly reset.
02231     initPatchFields<volScalarField, processorFvPatchField<scalar> >
02232     (
02233         pTraits<scalar>::zero
02234     );
02235     initPatchFields<volVectorField, processorFvPatchField<vector> >
02236     (
02237         pTraits<vector>::zero
02238     );
02239     initPatchFields
02240     <
02241         volSphericalTensorField,
02242         processorFvPatchField<sphericalTensor>
02243     >
02244     (
02245         pTraits<sphericalTensor>::zero
02246     );
02247     initPatchFields<volSymmTensorField, processorFvPatchField<symmTensor> >
02248     (
02249         pTraits<symmTensor>::zero
02250     );
02251     initPatchFields<volTensorField, processorFvPatchField<tensor> >
02252     (
02253         pTraits<tensor>::zero
02254     );
02255     initPatchFields<surfaceScalarField, processorFvsPatchField<scalar> >
02256     (
02257         pTraits<scalar>::zero
02258     );
02259     initPatchFields<surfaceVectorField, processorFvsPatchField<vector> >
02260     (
02261         pTraits<vector>::zero
02262     );
02263     initPatchFields
02264     <
02265         surfaceSphericalTensorField,
02266         processorFvsPatchField<sphericalTensor>
02267     >
02268     (
02269         pTraits<sphericalTensor>::zero
02270     );
02271     initPatchFields
02272     <
02273         surfaceSymmTensorField,
02274         processorFvsPatchField<symmTensor>
02275     >
02276     (
02277         pTraits<symmTensor>::zero
02278     );
02279     initPatchFields<surfaceTensorField, processorFvsPatchField<tensor> >
02280     (
02281         pTraits<tensor>::zero
02282     );
02283 
02284 
02285     mesh_.setInstance(mesh_.time().timeName());
02286 
02287 
02288     // Print a bit
02289     if (debug)
02290     {
02291         Pout<< nl << "FINAL MESH:" << endl;
02292         printMeshInfo(mesh_);
02293         printFieldInfo<volScalarField>(mesh_);
02294         printFieldInfo<volVectorField>(mesh_);
02295         printFieldInfo<volSphericalTensorField>(mesh_);
02296         printFieldInfo<volSymmTensorField>(mesh_);
02297         printFieldInfo<volTensorField>(mesh_);
02298         printFieldInfo<surfaceScalarField>(mesh_);
02299         printFieldInfo<surfaceVectorField>(mesh_);
02300         printFieldInfo<surfaceSphericalTensorField>(mesh_);
02301         printFieldInfo<surfaceSymmTensorField>(mesh_);
02302         printFieldInfo<surfaceTensorField>(mesh_);
02303         Pout<< nl << endl;
02304     }
02305 
02306     // Collect all maps and return
02307     return autoPtr<mapDistributePolyMesh>
02308     (
02309         new mapDistributePolyMesh
02310         (
02311             mesh_,
02312 
02313             nOldPoints,
02314             nOldFaces,
02315             nOldCells,
02316             oldPatchStarts,
02317             oldPatchNMeshPoints,
02318 
02319             subPointMap,
02320             subFaceMap,
02321             subCellMap,
02322             subPatchMap,
02323 
02324             constructPointMap,
02325             constructFaceMap,
02326             constructCellMap,
02327             constructPatchMap,
02328             true                // reuse storage
02329         )
02330     );
02331 }
02332 
02333 
02334 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines