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

regionSplit.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 "regionSplit.H"
00027 #include <OpenFOAM/cyclicPolyPatch.H>
00028 #include <OpenFOAM/processorPolyPatch.H>
00029 #include <OpenFOAM/globalIndex.H>
00030 #include <OpenFOAM/syncTools.H>
00031 
00032 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
00033 
00034 namespace Foam
00035 {
00036 
00037 defineTypeNameAndDebug(regionSplit, 0);
00038 
00039 }
00040 
00041 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
00042 
00043 // Handle (non-processor) coupled faces.
00044 void Foam::regionSplit::transferCoupledFaceRegion
00045 (
00046     const label faceI,
00047     const label otherFaceI,
00048 
00049     labelList& faceRegion,
00050     DynamicList<label>& newChangedFaces
00051 ) const
00052 {
00053     if (faceRegion[faceI] >= 0)
00054     {
00055         if (faceRegion[otherFaceI] == -1)
00056         {
00057             faceRegion[otherFaceI] = faceRegion[faceI];
00058             newChangedFaces.append(otherFaceI);
00059         }
00060         else if (faceRegion[otherFaceI] == -2)
00061         {
00062             // otherFaceI blocked but faceI is not. Is illegal for coupled
00063             // faces, not for explicit connections.
00064         }
00065         else if (faceRegion[otherFaceI] != faceRegion[faceI])
00066         {
00067             FatalErrorIn
00068             (
00069                   "regionSplit::transferCoupledFaceRegion"
00070                   "(const label, const label, labelList&, labelList&) const"
00071               )   << "Problem : coupled face " << faceI
00072                   << " on patch " << mesh_.boundaryMesh().whichPatch(faceI)
00073                   << " has region " << faceRegion[faceI]
00074                   << " but coupled face " << otherFaceI
00075                   << " has region " << faceRegion[otherFaceI]
00076                   << endl
00077                   << "Is your blocked faces specification"
00078                   << " synchronized across coupled boundaries?"
00079                   << abort(FatalError);
00080         }
00081     }
00082     else if (faceRegion[faceI] == -1)
00083     {
00084         if (faceRegion[otherFaceI] >= 0)
00085         {
00086             faceRegion[faceI] = faceRegion[otherFaceI];
00087             newChangedFaces.append(faceI);
00088         }
00089         else if (faceRegion[otherFaceI] == -2)
00090         {
00091             // otherFaceI blocked but faceI is not. Is illegal for coupled
00092             // faces, not for explicit connections.
00093         }
00094     }
00095 }
00096 
00097 
00098 void Foam::regionSplit::fillSeedMask
00099 (
00100     const List<labelPair>& explicitConnections,
00101     labelList& cellRegion,
00102     labelList& faceRegion,
00103     const label seedCellID,
00104     const label markValue
00105 ) const
00106 {
00107     // Do seed cell
00108     cellRegion[seedCellID] = markValue;
00109 
00110 
00111     // Collect faces on seed cell
00112     const cell& cFaces = mesh_.cells()[seedCellID];
00113 
00114     label nFaces = 0;
00115 
00116     labelList changedFaces(cFaces.size());
00117 
00118     forAll(cFaces, i)
00119     {
00120         label faceI = cFaces[i];
00121 
00122         if (faceRegion[faceI] == -1)
00123         {
00124             faceRegion[faceI] = markValue;
00125             changedFaces[nFaces++] = faceI;
00126         }
00127     }
00128     changedFaces.setSize(nFaces);
00129 
00130 
00131     // Loop over changed faces. MeshWave in small.
00132 
00133     while (changedFaces.size())
00134     {
00135         //if (debug)
00136         //{
00137         //    Pout<< "regionSplit::fillSeedMask : changedFaces:"
00138         //        << changedFaces.size() << endl;
00139         //}
00140 
00141         DynamicList<label> changedCells(changedFaces.size());
00142 
00143         forAll(changedFaces, i)
00144         {
00145             label faceI = changedFaces[i];
00146 
00147             label own = mesh_.faceOwner()[faceI];
00148 
00149             if (cellRegion[own] == -1)
00150             {
00151                 cellRegion[own] = markValue;
00152                 changedCells.append(own);
00153             }
00154 
00155             if (mesh_.isInternalFace(faceI))
00156             {
00157                 label nei = mesh_.faceNeighbour()[faceI];
00158 
00159                 if (cellRegion[nei] == -1)
00160                 {
00161                     cellRegion[nei] = markValue;
00162                     changedCells.append(nei);
00163                 }
00164             }
00165         }
00166 
00167 
00168         //if (debug)
00169         //{
00170         //    Pout<< "regionSplit::fillSeedMask : changedCells:"
00171         //        << changedCells.size() << endl;
00172         //}
00173 
00174         // Loop over changedCells and collect faces
00175         DynamicList<label> newChangedFaces(changedCells.size());
00176 
00177         forAll(changedCells, i)
00178         {
00179             label cellI = changedCells[i];
00180 
00181             const cell& cFaces = mesh_.cells()[cellI];
00182 
00183             forAll(cFaces, cFaceI)
00184             {
00185                 label faceI = cFaces[cFaceI];
00186 
00187                 if (faceRegion[faceI] == -1)
00188                 {
00189                     faceRegion[faceI] = markValue;
00190                     newChangedFaces.append(faceI);
00191                 }
00192             }
00193         }
00194 
00195 
00196         //if (debug)
00197         //{
00198         //    Pout<< "regionSplit::fillSeedMask : changedFaces before sync:"
00199         //        << changedFaces.size() << endl;
00200         //}
00201 
00202 
00203         // Check for changes to any locally coupled face.
00204         // Global connections are done later.
00205 
00206         const polyBoundaryMesh& patches = mesh_.boundaryMesh();
00207 
00208         forAll(patches, patchI)
00209         {
00210             const polyPatch& pp = patches[patchI];
00211 
00212             if (isA<cyclicPolyPatch>(pp))
00213             {
00214                 label faceI = pp.start();
00215 
00216                 label halfSz = pp.size()/2;
00217 
00218                 for (label i = 0; i < halfSz; i++)
00219                 {
00220                     label otherFaceI = refCast<const cyclicPolyPatch>(pp)
00221                         .transformGlobalFace(faceI);
00222 
00223                     transferCoupledFaceRegion
00224                     (
00225                         faceI,
00226                         otherFaceI,
00227                         faceRegion,
00228                         newChangedFaces
00229                     );
00230 
00231                     faceI++;
00232                 }
00233             }
00234         }
00235         forAll(explicitConnections, i)
00236         {
00237             transferCoupledFaceRegion
00238             (
00239                 explicitConnections[i][0],
00240                 explicitConnections[i][1],
00241                 faceRegion,
00242                 newChangedFaces
00243             );
00244         }
00245 
00246         //if (debug)
00247         //{
00248         //    Pout<< "regionSplit::fillSeedMask : changedFaces after sync:"
00249         //        << newChangedFaces.size() << endl;
00250         //}
00251 
00252         changedFaces.transfer(newChangedFaces);
00253     }
00254 }
00255 
00256 
00257 Foam::label Foam::regionSplit::calcRegionSplit
00258 (
00259     const boolList& blockedFace,
00260     const List<labelPair>& explicitConnections,
00261 
00262     labelList& cellRegion
00263 ) const
00264 {
00265     if (debug)
00266     {
00267         if (blockedFace.size())
00268         {
00269             // Check that blockedFace is synced.
00270             boolList syncBlockedFace(blockedFace);
00271             syncTools::swapFaceList(mesh_, syncBlockedFace, false);
00272 
00273             forAll(syncBlockedFace, faceI)
00274             {
00275                 if (syncBlockedFace[faceI] != blockedFace[faceI])
00276                 {
00277                     FatalErrorIn
00278                     (
00279                         "regionSplit::calcRegionSplit(..)"
00280                     )   << "Face " << faceI << " not synchronised. My value:"
00281                         << blockedFace[faceI] << "  coupled value:"
00282                         << syncBlockedFace[faceI]
00283                         << abort(FatalError);
00284                 }
00285             }
00286         }
00287     }
00288 
00289     // Region per face.
00290     // -1 unassigned
00291     // -2 blocked
00292     labelList faceRegion(mesh_.nFaces(), -1);
00293 
00294     if (blockedFace.size())
00295     {
00296         forAll(blockedFace, faceI)
00297         {
00298             if (blockedFace[faceI])
00299             {
00300                 faceRegion[faceI] = -2;
00301             }
00302         }
00303     }
00304 
00305 
00306     // Assign local regions
00307     // ~~~~~~~~~~~~~~~~~~~~
00308 
00309     // Start with region 0
00310     label nRegions = 0;
00311 
00312     label unsetCellI = 0;
00313 
00314     do
00315     {
00316         // Find first unset cell
00317 
00318         for (; unsetCellI < mesh_.nCells(); unsetCellI++)
00319         {
00320             if (cellRegion[unsetCellI] == -1)
00321             {
00322                 break;
00323             }
00324         }
00325 
00326         if (unsetCellI >= mesh_.nCells())
00327         {
00328             break;
00329         }
00330 
00331         fillSeedMask
00332         (
00333             explicitConnections,
00334             cellRegion,
00335             faceRegion,
00336             unsetCellI,
00337             nRegions
00338         );
00339 
00340         // Current unsetCell has now been handled. Go to next region.
00341         nRegions++;
00342         unsetCellI++;
00343     }
00344     while(true);
00345 
00346 
00347     if (debug)
00348     {
00349         forAll(cellRegion, cellI)
00350         {
00351             if (cellRegion[cellI] < 0)
00352             {
00353                 FatalErrorIn("regionSplit::calcRegionSplit(..)")
00354                     << "cell:" << cellI << " region:" << cellRegion[cellI]
00355                     << abort(FatalError);
00356             }
00357         }
00358 
00359         forAll(faceRegion, faceI)
00360         {
00361             if (faceRegion[faceI] == -1)
00362             {
00363                 FatalErrorIn("regionSplit::calcRegionSplit(..)")
00364                     << "face:" << faceI << " region:" << faceRegion[faceI]
00365                     << abort(FatalError);
00366             }
00367         }
00368     }
00369 
00370 
00371 
00372     // Assign global regions
00373     // ~~~~~~~~~~~~~~~~~~~~~
00374     // Offset local regions to create unique global regions.
00375 
00376     globalIndex globalRegions(nRegions);
00377 
00378 
00379     // Merge global regions
00380     // ~~~~~~~~~~~~~~~~~~~~
00381     // Regions across non-blocked proc patches get merged.
00382     // This will set merged global regions to be the min of both.
00383     // (this will create gaps in the global region list so they will get
00384     // merged later on)
00385 
00386     // Map from global to merged global
00387     labelList mergedGlobal(identity(globalRegions.size()));
00388 
00389 
00390     // See if any regions get merged. Only nessecary for parallel
00391     while (Pstream::parRun())
00392     {
00393         if (debug)
00394         {
00395             Pout<< nl << "-- Starting Iteration --" << endl;
00396         }
00397 
00398         const polyBoundaryMesh& patches = mesh_.boundaryMesh();
00399 
00400         // Send global regions across (or -2 if blocked face)
00401         forAll(patches, patchI)
00402         {
00403             const polyPatch& pp = patches[patchI];
00404 
00405             if (isA<processorPolyPatch>(pp))
00406             {
00407                 labelList myGlobalRegions(pp.size());
00408 
00409                 label faceI = pp.start();
00410 
00411                 forAll(pp, i)
00412                 {
00413                     if (faceRegion[faceI] < 0)
00414                     {
00415                         myGlobalRegions[i] = faceRegion[faceI];
00416                     }
00417                     else
00418                     {
00419                         myGlobalRegions[i] = mergedGlobal
00420                         [globalRegions.toGlobal(faceRegion[faceI])];
00421                     }
00422 
00423                     faceI++;
00424                 }
00425 
00426                 OPstream toProcNbr
00427                 (
00428                     Pstream::blocking,
00429                     refCast<const processorPolyPatch>(pp).neighbProcNo()
00430                 );
00431 
00432                 toProcNbr << myGlobalRegions;
00433             }
00434         }
00435 
00436 
00437         // Receive global regions
00438 
00439         label nMerged = 0;
00440 
00441         forAll(patches, patchI)
00442         {
00443             const polyPatch& pp = patches[patchI];
00444 
00445             if (isA<processorPolyPatch>(pp))
00446             {
00447                 const processorPolyPatch& procPp =
00448                     refCast<const processorPolyPatch>(pp);
00449 
00450                 IPstream fromProcNbr(Pstream::blocking, procPp.neighbProcNo());
00451 
00452                 labelList nbrRegions(fromProcNbr);
00453 
00454 
00455                 // Compare with my regions to see which get merged.
00456 
00457                 label faceI = pp.start();
00458 
00459                 forAll(pp, i)
00460                 {
00461                     if
00462                     (
00463                         faceRegion[faceI] < 0
00464                      || nbrRegions[i] < 0
00465                     )
00466                     {
00467                         if (faceRegion[faceI] != nbrRegions[i])
00468                         {
00469                             FatalErrorIn("regionSplit::calcRegionSplit(..)")
00470                                 << "On patch:" << pp.name()
00471                                 << " face:" << faceI
00472                                 << " my local region:" << faceRegion[faceI]
00473                                 << " neighbouring region:"
00474                                 << nbrRegions[i] << nl
00475                                 << "Maybe your blockedFaces are not"
00476                                 << " synchronized across coupled faces?"
00477                                 << abort(FatalError);
00478                         }
00479                     }
00480                     else
00481                     {
00482                         label uncompactGlobal =
00483                             globalRegions.toGlobal(faceRegion[faceI]);
00484 
00485                         label myGlobal = mergedGlobal[uncompactGlobal];
00486 
00487                         if (myGlobal != nbrRegions[i])
00488                         {
00489                             label minRegion = min(myGlobal, nbrRegions[i]);
00490 
00491                             if (debug)
00492                             {
00493                                 Pout<< "Merging region " << myGlobal
00494                                     << " (on proc " << Pstream::myProcNo()
00495                                     << ") and region " << nbrRegions[i]
00496                                     << " (on proc " << procPp.neighbProcNo()
00497                                     << ") into region " << minRegion << endl;
00498                             }
00499 
00500                             mergedGlobal[uncompactGlobal] = minRegion;
00501                             mergedGlobal[myGlobal] = minRegion;
00502                             mergedGlobal[nbrRegions[i]] = minRegion;
00503 
00504                             nMerged++;
00505                         }
00506                     }
00507 
00508                     faceI++;
00509                 }
00510             }
00511         }
00512 
00513 
00514         reduce(nMerged, sumOp<label>());
00515 
00516         if (debug)
00517         {
00518             Pout<< "nMerged:" << nMerged << endl;
00519         }
00520 
00521         if (nMerged == 0)
00522         {
00523             break;
00524         }
00525 
00526         // Merge the compacted regions.
00527         Pstream::listCombineGather(mergedGlobal, minEqOp<label>());
00528         Pstream::listCombineScatter(mergedGlobal);
00529     }
00530 
00531 
00532     // Compact global regions
00533     // ~~~~~~~~~~~~~~~~~~~~~~
00534 
00535     // All procs will have the same global mergedGlobal region.
00536     // There might be gaps in it however so compact.
00537 
00538     labelList mergedToCompacted(globalRegions.size(), -1);
00539 
00540     label compactI = 0;
00541 
00542     forAll(mergedGlobal, i)
00543     {
00544         label merged = mergedGlobal[i];
00545 
00546         if (mergedToCompacted[merged] == -1)
00547         {
00548             mergedToCompacted[merged] = compactI++;
00549         }
00550     }
00551 
00552     if (debug)
00553     {
00554         Pout<< "Compacted down to " << compactI << " regions." << endl;
00555     }
00556 
00557     // Renumber cellRegion to be global regions
00558     forAll(cellRegion, cellI)
00559     {
00560         label region = cellRegion[cellI];
00561 
00562         if (region >= 0)
00563         {
00564             label merged = mergedGlobal[globalRegions.toGlobal(region)];
00565 
00566             cellRegion[cellI] = mergedToCompacted[merged];
00567         }
00568     }
00569 
00570     return compactI;
00571 }
00572 
00573 
00574 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
00575 
00576 Foam::regionSplit::regionSplit(const polyMesh& mesh)
00577 :
00578     labelList(mesh.nCells(), -1),
00579     mesh_(mesh),
00580     nRegions_(calcRegionSplit(boolList(0, false), List<labelPair>(0), *this))
00581 {}
00582 
00583 
00584 Foam::regionSplit::regionSplit
00585 (
00586     const polyMesh& mesh,
00587     const boolList& blockedFace
00588 )
00589 :
00590     labelList(mesh.nCells(), -1),
00591     mesh_(mesh),
00592     nRegions_(calcRegionSplit(blockedFace, List<labelPair>(0), *this))
00593 {}
00594 
00595 
00596 Foam::regionSplit::regionSplit
00597 (
00598     const polyMesh& mesh,
00599     const boolList& blockedFace,
00600     const List<labelPair>& explicitConnections
00601 )
00602 :
00603     labelList(mesh.nCells(), -1),
00604     mesh_(mesh),
00605     nRegions_(calcRegionSplit(blockedFace, explicitConnections, *this))
00606 {}
00607 
00608 
00609 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00610 
00611 
00612 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines