src/foam/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomerateLduAddressing.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "GAMGAgglomeration.H"
  27 #include "GAMGInterface.H"
  28 #include "FieldFields.H"
  29
  30 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  31
  32 void Foam::GAMGAgglomeration::agglomerateLduAddressing
  33 (
  34     const label fineLevelIndex
  35 )
  36 {
  37     const lduMesh& fineMesh = meshLevel(fineLevelIndex);
  38     const lduAddressing& fineMeshAddr = fineMesh.lduAddr();
  39
  40     const unallocLabelList& upperAddr = fineMeshAddr.upperAddr();
  41     const unallocLabelList& lowerAddr = fineMeshAddr.lowerAddr();
  42
  43     label nFineFaces = upperAddr.size();
  44
  45     // Get restriction map for current level
  46     const labelField& restrictMap = restrictAddressing(fineLevelIndex);
  47
  48     if (min(restrictMap) == -1)
  49     {
  50         FatalErrorIn("GAMGAgglomeration::agglomerateLduAddressing")
  51             << "min(restrictMap) == -1" << exit(FatalError);
  52     }
  53
  54     if (restrictMap.size() != fineMeshAddr.size())
  55     {
  56         FatalErrorIn
  57         (
  58             "GAMGAgglomeration::agglomerateLduAddressing"
  59             "(const label fineLevelIndex)"
  60         )   << "restrict map does not correspond to fine level. " << endl
  61             << " Sizes: restrictMap: " << restrictMap.size()
  62             << " nEqns: " << fineMeshAddr.size()
  63             << abort(FatalError);
  64     }
  65
  66
  67     // Get the number of coarse cells
  68     const label nCoarseCells = nCells_[fineLevelIndex];
  69
  70     // Storage for coarse cell neighbours and coefficients
  71
  72     // Guess initial maximum number of neighbours in coarse cell
  73     label maxNnbrs = 10;
  74
  75     // Number of faces for each coarse-cell
  76     labelList cCellnFaces(nCoarseCells, 0);
  77
  78     // Setup initial packed storage for coarse-cell faces
  79     labelList cCellFaces(maxNnbrs*nCoarseCells);
  80
  81     // Create face-restriction addressing
  82     faceRestrictAddressing_.set(fineLevelIndex, new labelList(nFineFaces));
  83     labelList& faceRestrictAddr = faceRestrictAddressing_[fineLevelIndex];
  84
  85     // Initial neighbour array (not in upper-triangle order)
  86     labelList initCoarseNeighb(nFineFaces);
  87
  88     // Counter for coarse faces
  89     label nCoarseFaces = 0;
  90
  91     // Loop through all fine faces
  92     forAll (upperAddr, fineFacei)
  93     {
  94         label rmUpperAddr = restrictMap[upperAddr[fineFacei]];
  95         label rmLowerAddr = restrictMap[lowerAddr[fineFacei]];
  96
  97         if (rmUpperAddr == rmLowerAddr)
  98         {
  99             // For each fine face inside of a coarse cell keep the address
 100             // of the cell corresponding to the face in the faceRestrictAddr
 101             // as a negative index
 102             faceRestrictAddr[fineFacei] = -(rmUpperAddr + 1);
 103         }
 104         else
 105         {
 106             // this face is a part of a coarse face
 107
 108             label cOwn = rmUpperAddr;
 109             label cNei = rmLowerAddr;
 110
 111             // get coarse owner and neighbour
 112             if (rmUpperAddr > rmLowerAddr)
 113             {
 114                 cOwn = rmLowerAddr;
 115                 cNei = rmUpperAddr;
 116             }
 117
 118             // check the neighbour to see if this face has already been found
 119             label* ccFaces = &cCellFaces[maxNnbrs*cOwn];
 120
 121             bool nbrFound = false;
 122             label& ccnFaces = cCellnFaces[cOwn];
 123
 124             for (int i=0; i<ccnFaces; i++)
 125             {
 126                 if (initCoarseNeighb[ccFaces[i]] == cNei)
 127                 {
 128                     nbrFound = true;
 129                     faceRestrictAddr[fineFacei] = ccFaces[i];
 130                     break;
 131                 }
 132             }
 133
 134             if (!nbrFound)
 135             {
 136                 if (ccnFaces >= maxNnbrs)
 137                 {
 138                     label oldMaxNnbrs = maxNnbrs;
 139                     maxNnbrs *= 2;
 140
 141                     cCellFaces.setSize(maxNnbrs*nCoarseCells);
 142
 143                     forAllReverse(cCellnFaces, i)
 144                     {
 145                         label* oldCcNbrs = &cCellFaces[oldMaxNnbrs*i];
 146                         label* newCcNbrs = &cCellFaces[maxNnbrs*i];
 147
 148                         for (int j=0; j<cCellnFaces[i]; j++)
 149                         {
 150                             newCcNbrs[j] = oldCcNbrs[j];
 151                         }
 152                     }
 153
 154                     ccFaces = &cCellFaces[maxNnbrs*cOwn];
 155                 }
 156
 157                 ccFaces[ccnFaces] = nCoarseFaces;
 158                 initCoarseNeighb[nCoarseFaces] = cNei;
 159                 faceRestrictAddr[fineFacei] = nCoarseFaces;
 160                 ccnFaces++;
 161
 162                 // new coarse face created
 163                 nCoarseFaces++;
 164             }
 165         }
 166     } // end for all fine faces
 167
 168
 169     // Renumber into upper-triangular order
 170
 171     // All coarse owner-neighbour storage
 172     labelList coarseOwner(nCoarseFaces);
 173     labelList coarseNeighbour(nCoarseFaces);
 174     labelList coarseFaceMap(nCoarseFaces);
 175
 176     label coarseFacei = 0;
 177
 178     forAll (cCellnFaces, cci)
 179     {
 180         label* cFaces = &cCellFaces[maxNnbrs*cci];
 181         label ccnFaces = cCellnFaces[cci];
 182
 183         for (int i=0; i<ccnFaces; i++)
 184         {
 185             coarseOwner[coarseFacei] = cci;
 186             coarseNeighbour[coarseFacei] = initCoarseNeighb[cFaces[i]];
 187             coarseFaceMap[cFaces[i]] = coarseFacei;
 188             coarseFacei++;
 189         }
 190     }
 191
 192     forAll(faceRestrictAddr, fineFacei)
 193     {
 194         if (faceRestrictAddr[fineFacei] >= 0)
 195         {
 196             faceRestrictAddr[fineFacei] =
 197                 coarseFaceMap[faceRestrictAddr[fineFacei]];
 198         }
 199     }
 200
 201     // Clear the temporary storage for the coarse cell data
 202     cCellnFaces.setSize(0);
 203     cCellFaces.setSize(0);
 204     initCoarseNeighb.setSize(0);
 205     coarseFaceMap.setSize(0);
 206
 207
 208     // Create coarse-level interfaces
 209
 210     // Get reference to fine-level interfaces
 211     const lduInterfacePtrsList& fineInterfaces =
 212         interfaceLevels_[fineLevelIndex];
 213
 214     // Initialise transfer of restrict addressing on the interface
 215     forAll (fineInterfaces, inti)
 216     {
 217         if (fineInterfaces.set(inti))
 218         {
 219             if (fineInterfaces[inti].coupled())
 220             {
 221                 fineInterfaces[inti].initInternalFieldTransfer
 222                 (
 223                     Pstream::blocking,
 224                     restrictMap
 225                 );
 226             }
 227         }
 228     }
 229
 230     // Store coefficients to avoid tangled communications
 231     // HJ, 1/Apr/2009
 232     FieldField<Field, label> fineInterfaceAddr(fineInterfaces.size());
 233
 234     forAll (fineInterfaces, inti)
 235     {
 236         if (fineInterfaces.set(inti))
 237         {
 238             if (fineInterfaces[inti].coupled())
 239             {
 240                 fineInterfaceAddr.set
 241                 (
 242                     inti,
 243                     new labelField
 244                     (
 245                         fineInterfaces[inti].internalFieldTransfer
 246                         (
 247                             Pstream::blocking,
 248                             restrictMap
 249                         )
 250                     )
 251                 );
 252             }
 253         }
 254     }
 255
 256     // Add the coarse level
 257
 258     // Set the coarse ldu addressing onto the list
 259     meshLevels_.set
 260     (
 261         fineLevelIndex,
 262         new lduPrimitiveMesh
 263         (
 264             nCoarseCells,
 265             coarseOwner,
 266             coarseNeighbour,
 267             true
 268         )
 269     );
 270
 271     // Create coarse-level interfaces
 272     interfaceLevels_.set
 273     (
 274         fineLevelIndex + 1,
 275         new lduInterfacePtrsList(fineInterfaces.size())
 276     );
 277
 278     lduInterfacePtrsList& coarseInterfaces =
 279         interfaceLevels_[fineLevelIndex + 1];
 280
 281     labelListList coarseInterfaceAddr(fineInterfaces.size());
 282
 283     forAll (fineInterfaces, inti)
 284     {
 285         if (fineInterfaces.set(inti))
 286         {
 287             if (fineInterfaces[inti].coupled())
 288             {
 289                 coarseInterfaces.set
 290                 (
 291                     inti,
 292                     GAMGInterface::New
 293                     (
 294                         meshLevels_[fineLevelIndex],
 295                         fineInterfaces[inti],
 296                         fineInterfaces[inti].interfaceInternalField
 297                         (
 298                             restrictMap
 299                         ),
 300                         fineInterfaceAddr[inti]
 301                     ).ptr()
 302                 );
 303
 304                 coarseInterfaceAddr[inti] = coarseInterfaces[inti].faceCells();
 305             }
 306         }
 307     }
 308
 309     meshLevels_[fineLevelIndex].addInterfaces
 310     (
 311         coarseInterfaces,
 312         coarseInterfaceAddr,
 313         fineMeshAddr.patchSchedule()
 314     );
 315 }
 316
 317
 318 // ************************************************************************* //