src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   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 InClass
  25     decompositionMethod
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "decompositionMethod.H"
  30 #include "globalIndex.H"
  31 #include "cyclicPolyPatch.H"
  32 #include "syncTools.H"
  33
  34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  35
  36 namespace Foam
  37 {
  38     defineTypeNameAndDebug(decompositionMethod, 0);
  39     defineRunTimeSelectionTable(decompositionMethod, dictionary);
  40 }
  41
  42 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  43
  44 Foam::autoPtr<Foam::decompositionMethod> Foam::decompositionMethod::New
  45 (
  46     const dictionary& decompositionDict
  47 )
  48 {
  49     const word methodType(decompositionDict.lookup("method"));
  50
  51     Info<< "Selecting decompositionMethod " << methodType << endl;
  52
  53     dictionaryConstructorTable::iterator cstrIter =
  54         dictionaryConstructorTablePtr_->find(methodType);
  55
  56     if (cstrIter == dictionaryConstructorTablePtr_->end())
  57     {
  58         FatalErrorIn
  59         (
  60             "decompositionMethod::New"
  61             "(const dictionary& decompositionDict)"
  62         )   << "Unknown decompositionMethod "
  63             << methodType << nl << nl
  64             << "Valid decompositionMethods are : " << endl
  65             << dictionaryConstructorTablePtr_->sortedToc()
  66             << exit(FatalError);
  67     }
  68
  69     return autoPtr<decompositionMethod>(cstrIter()(decompositionDict));
  70 }
  71
  72
  73 Foam::labelList Foam::decompositionMethod::decompose
  74 (
  75     const polyMesh& mesh,
  76     const pointField& points
  77 )
  78 {
  79     scalarField weights(points.size(), 1.0);
  80
  81     return decompose(mesh, points, weights);
  82 }
  83
  84
  85 Foam::labelList Foam::decompositionMethod::decompose
  86 (
  87     const polyMesh& mesh,
  88     const labelList& fineToCoarse,
  89     const pointField& coarsePoints,
  90     const scalarField& coarseWeights
  91 )
  92 {
  93     CompactListList<label> coarseCellCells;
  94     calcCellCells
  95     (
  96         mesh,
  97         fineToCoarse,
  98         coarsePoints.size(),
  99         coarseCellCells
 100     );
 101
 102     // Decompose based on agglomerated points
 103     labelList coarseDistribution
 104     (
 105         decompose
 106         (
 107             coarseCellCells(),
 108             coarsePoints,
 109             coarseWeights
 110         )
 111     );
 112
 113     // Rework back into decomposition for original mesh_
 114     labelList fineDistribution(fineToCoarse.size());
 115
 116     forAll(fineDistribution, i)
 117     {
 118         fineDistribution[i] = coarseDistribution[fineToCoarse[i]];
 119     }
 120
 121     return fineDistribution;
 122 }
 123
 124
 125 Foam::labelList Foam::decompositionMethod::decompose
 126 (
 127     const polyMesh& mesh,
 128     const labelList& fineToCoarse,
 129     const pointField& coarsePoints
 130 )
 131 {
 132     scalarField cWeights(coarsePoints.size(), 1.0);
 133
 134     return decompose
 135     (
 136         mesh,
 137         fineToCoarse,
 138         coarsePoints,
 139         cWeights
 140     );
 141 }
 142
 143
 144 Foam::labelList Foam::decompositionMethod::decompose
 145 (
 146     const labelListList& globalCellCells,
 147     const pointField& cc
 148 )
 149 {
 150     scalarField cWeights(cc.size(), 1.0);
 151
 152     return decompose(globalCellCells, cc, cWeights);
 153 }
 154
 155
 156 void Foam::decompositionMethod::calcCellCells
 157 (
 158     const polyMesh& mesh,
 159     const labelList& agglom,
 160     const label nCoarse,
 161     CompactListList<label>& cellCells
 162 )
 163 {
 164     const labelList& faceOwner = mesh.faceOwner();
 165     const labelList& faceNeighbour = mesh.faceNeighbour();
 166     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 167
 168
 169     // Create global cell numbers
 170     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
 171
 172     globalIndex globalAgglom(nCoarse);
 173
 174
 175     // Get agglomerate owner on other side of coupled faces
 176     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 177
 178     labelList globalNeighbour(mesh.nFaces()-mesh.nInternalFaces());
 179
 180     forAll(patches, patchI)
 181     {
 182         const polyPatch& pp = patches[patchI];
 183
 184         if (pp.coupled())
 185         {
 186             label faceI = pp.start();
 187             label bFaceI = pp.start() - mesh.nInternalFaces();
 188
 189             forAll(pp, i)
 190             {
 191                 globalNeighbour[bFaceI] = globalAgglom.toGlobal
 192                 (
 193                     agglom[faceOwner[faceI]]
 194                 );
 195
 196                 bFaceI++;
 197                 faceI++;
 198             }
 199         }
 200     }
 201
 202     // Get the cell on the other side of coupled patches
 203     syncTools::swapBoundaryFaceList(mesh, globalNeighbour);
 204
 205
 206     // Count number of faces (internal + coupled)
 207     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 208
 209     // Number of faces per coarse cell
 210     labelList nFacesPerCell(nCoarse, 0);
 211
 212     for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
 213     {
 214         label own = agglom[faceOwner[faceI]];
 215         label nei = agglom[faceNeighbour[faceI]];
 216
 217         nFacesPerCell[own]++;
 218         nFacesPerCell[nei]++;
 219     }
 220
 221     forAll(patches, patchI)
 222     {
 223         const polyPatch& pp = patches[patchI];
 224
 225         if (pp.coupled())
 226         {
 227             label faceI = pp.start();
 228             label bFaceI = pp.start()-mesh.nInternalFaces();
 229
 230             forAll(pp, i)
 231             {
 232                 label own = agglom[faceOwner[faceI]];
 233
 234                 label globalNei = globalNeighbour[bFaceI];
 235                 if
 236                 (
 237                    !globalAgglom.isLocal(globalNei)
 238                  || globalAgglom.toLocal(globalNei) != own
 239                 )
 240                 {
 241                     nFacesPerCell[own]++;
 242                 }
 243
 244                 faceI++;
 245                 bFaceI++;
 246             }
 247         }
 248     }
 249
 250
 251     // Fill in offset and data
 252     // ~~~~~~~~~~~~~~~~~~~~~~~
 253
 254     cellCells.setSize(nFacesPerCell);
 255
 256     nFacesPerCell = 0;
 257
 258     labelList& m = cellCells.m();
 259     const labelList& offsets = cellCells.offsets();
 260
 261     // For internal faces is just offsetted owner and neighbour
 262     for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
 263     {
 264         label own = agglom[faceOwner[faceI]];
 265         label nei = agglom[faceNeighbour[faceI]];
 266
 267         m[offsets[own] + nFacesPerCell[own]++] = globalAgglom.toGlobal(nei);
 268         m[offsets[nei] + nFacesPerCell[nei]++] = globalAgglom.toGlobal(own);
 269     }
 270
 271     // For boundary faces is offsetted coupled neighbour
 272     forAll(patches, patchI)
 273     {
 274         const polyPatch& pp = patches[patchI];
 275
 276         if (pp.coupled())
 277         {
 278             label faceI = pp.start();
 279             label bFaceI = pp.start()-mesh.nInternalFaces();
 280
 281             forAll(pp, i)
 282             {
 283                 label own = agglom[faceOwner[faceI]];
 284
 285                 label globalNei = globalNeighbour[bFaceI];
 286
 287                 if
 288                 (
 289                    !globalAgglom.isLocal(globalNei)
 290                  || globalAgglom.toLocal(globalNei) != own
 291                 )
 292                 {
 293                     m[offsets[own] + nFacesPerCell[own]++] = globalNei;
 294                 }
 295
 296                 faceI++;
 297                 bFaceI++;
 298             }
 299         }
 300     }
 301
 302
 303     // Check for duplicates connections between cells
 304     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 305     // Done as postprocessing step since we now have cellCells.
 306     label newIndex = 0;
 307     labelHashSet nbrCells;
 308     forAll(cellCells, cellI)
 309     {
 310         nbrCells.clear();
 311         nbrCells.insert(globalAgglom.toGlobal(cellI));
 312
 313         label startIndex = cellCells.offsets()[cellI];
 314         label endIndex = cellCells.offsets()[cellI+1];
 315
 316         for (label i = startIndex; i < endIndex; i++)
 317         {
 318             if (nbrCells.insert(cellCells.m()[i]))
 319             {
 320                 cellCells.m()[newIndex++] = cellCells.m()[i];
 321             }
 322         }
 323         cellCells.offsets()[cellI+1] = newIndex;
 324     }
 325
 326     cellCells.m().setSize(newIndex);
 327
 328
 329     //forAll(cellCells, cellI)
 330     //{
 331     //    const labelUList cCells = cellCells[cellI];
 332     //    Pout<< "Compacted: Coarse cell " << cellI << endl;
 333     //    forAll(cCells, i)
 334     //    {
 335     //        Pout<< "    nbr:" << cCells[i] << endl;
 336     //    }
 337     //}
 338 }
 339
 340
 341 // ************************************************************************* //