applications/utilities/parallelProcessing/decomposePar/distributeCells.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "domainDecomposition.H"
  27 #include "decompositionMethod.H"
  28 #include "cpuTime.H"
  29 #include "cyclicPolyPatch.H"
  30 #include "cellSet.H"
  31 #include "regionSplit.H"
  32
  33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  34
  35 void domainDecomposition::distributeCells()
  36 {
  37     Info<< "\nCalculating distribution of cells" << endl;
  38
  39     cpuTime decompositionTime;
  40
  41
  42     // See if any faces need to have owner and neighbour on same processor
  43     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  44
  45     labelHashSet sameProcFaces;
  46
  47     if (decompositionDict_.found("preservePatches"))
  48     {
  49         wordList pNames(decompositionDict_.lookup("preservePatches"));
  50
  51         Info<< "Keeping owner of faces in patches " << pNames
  52             << " on same processor. This only makes sense for cyclics." << endl;
  53
  54         const polyBoundaryMesh& patches = boundaryMesh();
  55
  56         forAll(pNames, i)
  57         {
  58             label patchI = patches.findPatchID(pNames[i]);
  59
  60             if (patchI == -1)
  61             {
  62                 FatalErrorIn("domainDecomposition::distributeCells()")
  63                     << "Unknown preservePatch " << pNames[i]
  64                     << endl << "Valid patches are " << patches.names()
  65                     << exit(FatalError);
  66             }
  67
  68             const polyPatch& pp = patches[patchI];
  69
  70             forAll(pp, i)
  71             {
  72                 sameProcFaces.insert(pp.start() + i);
  73             }
  74         }
  75     }
  76     if (decompositionDict_.found("preserveFaceZones"))
  77     {
  78         wordList zNames(decompositionDict_.lookup("preserveFaceZones"));
  79
  80         Info<< "Keeping owner and neighbour of faces in zones " << zNames
  81             << " on same processor" << endl;
  82
  83         const faceZoneMesh& fZones = faceZones();
  84
  85         forAll(zNames, i)
  86         {
  87             label zoneI = fZones.findZoneID(zNames[i]);
  88
  89             if (zoneI == -1)
  90             {
  91                 FatalErrorIn("domainDecomposition::distributeCells()")
  92                     << "Unknown preserveFaceZone " << zNames[i]
  93                     << endl << "Valid faceZones are " << fZones.names()
  94                     << exit(FatalError);
  95             }
  96
  97             const faceZone& fz = fZones[zoneI];
  98
  99             forAll(fz, i)
 100             {
 101                 sameProcFaces.insert(fz[i]);
 102             }
 103         }
 104     }
 105
 106
 107     // Construct decomposition method and either do decomposition on
 108     // cell centres or on agglomeration
 109
 110
 111     autoPtr<decompositionMethod> decomposePtr = decompositionMethod::New
 112     (
 113         decompositionDict_,
 114         *this
 115     );
 116
 117     if (sameProcFaces.empty())
 118     {
 119         if (decompositionDict_.found("weightField"))
 120         {
 121             word weightName = decompositionDict_.lookup("weightField");
 122
 123             volScalarField weights
 124             (
 125                 IOobject
 126                 (
 127                     weightName,
 128                     time().timeName(),
 129                     *this,
 130                     IOobject::MUST_READ,
 131                     IOobject::NO_WRITE
 132                 ),
 133                 *this
 134             );
 135
 136             cellToProc_ = decomposePtr().decompose
 137             (
 138                 cellCentres(),
 139                 weights.internalField()
 140             );
 141         }
 142         else
 143         {
 144             cellToProc_ = decomposePtr().decompose(cellCentres());
 145         }
 146     }
 147     else
 148     {
 149         Info<< "Selected " << sameProcFaces.size()
 150             << " faces whose owner and neighbour cell should be kept on the"
 151             << " same processor" << endl;
 152
 153         if (decompositionDict_.found("weightField"))
 154         {
 155             WarningIn("void domainDecomposition::distributeCells()")
 156                 << "weightField not supported when agglomerated "
 157                 << "decomposition required" << endl;
 158         }
 159
 160         // Faces where owner and neighbour are not 'connected' (= all except
 161         // sameProcFaces)
 162         boolList blockedFace(nFaces(), true);
 163
 164         forAllConstIter(labelHashSet, sameProcFaces, iter)
 165         {
 166             blockedFace[iter.key()] = false;
 167         }
 168
 169         // Connect coupled boundary faces
 170         const polyBoundaryMesh& patches =  boundaryMesh();
 171
 172         forAll(patches, patchI)
 173         {
 174             const polyPatch& pp = patches[patchI];
 175
 176             if (pp.coupled())
 177             {
 178                 forAll(pp, i)
 179                 {
 180                     blockedFace[pp.start()+i] = false;
 181                 }
 182             }
 183         }
 184
 185         // Determine global regions, separated by blockedFaces
 186         regionSplit globalRegion(*this, blockedFace);
 187
 188
 189         // Determine region cell centres
 190         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 191
 192         // This just takes the first cell in the region. Otherwise the problem
 193         // is with cyclics - if we'd average the region centre might be
 194         // somewhere in the middle of the domain which might not be anywhere
 195         // near any of the cells.
 196
 197         const point greatPoint(GREAT, GREAT, GREAT);
 198
 199         pointField regionCentres(globalRegion.nRegions(), greatPoint);
 200
 201         forAll(globalRegion, cellI)
 202         {
 203             label regionI = globalRegion[cellI];
 204
 205             if (regionCentres[regionI] == greatPoint)
 206             {
 207                 regionCentres[regionI] = cellCentres()[cellI];
 208             }
 209         }
 210
 211         // Do decomposition on agglomeration
 212         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 213         cellToProc_ = decomposePtr().decompose(globalRegion, regionCentres);
 214     }
 215
 216     Info<< "\nFinished decomposition in "
 217         << decompositionTime.elapsedCpuTime()
 218         << " s" << endl;
 219 }
 220
 221
 222 // ************************************************************************* //