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

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "GAMGAgglomeration.H"
  28 #include "GAMGInterface.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     // Create coarse-level interfaces
 215     interfaceLevels_.set
 216     (
 217         fineLevelIndex + 1,
 218         new lduInterfacePtrsList(fineInterfaces.size())
 219     );
 220
 221     lduInterfacePtrsList& coarseInterfaces =
 222         interfaceLevels_[fineLevelIndex + 1];
 223
 224     labelListList coarseInterfaceAddr(fineInterfaces.size());
 225
 226     // Initialise transfer of restrict addressing on the interface
 227     forAll (fineInterfaces, inti)
 228     {
 229         if (fineInterfaces.set(inti))
 230         {
 231             fineInterfaces[inti].initInternalFieldTransfer
 232             (
 233                 Pstream::blocking,
 234                 restrictMap
 235             );
 236         }
 237     }
 238
 239     // Add the coarse level
 240     forAll (fineInterfaces, inti)
 241     {
 242         if (fineInterfaces.set(inti))
 243         {
 244             coarseInterfaces.set
 245             (
 246                 inti,
 247                 GAMGInterface::New
 248                 (
 249                     fineInterfaces[inti],
 250                     fineInterfaces[inti].interfaceInternalField(restrictMap),
 251                     fineInterfaces[inti].internalFieldTransfer
 252                     (
 253                         Pstream::blocking,
 254                         restrictMap
 255                     )
 256                 ).ptr()
 257             );
 258
 259             coarseInterfaceAddr[inti] = coarseInterfaces[inti].faceCells();
 260         }
 261     }
 262
 263
 264     // Set the coarse ldu addressing onto the list
 265     meshLevels_.set
 266     (
 267         fineLevelIndex,
 268         new lduPrimitiveMesh
 269         (
 270             nCoarseCells,
 271             coarseOwner,
 272             coarseNeighbour,
 273             coarseInterfaceAddr,
 274             coarseInterfaces,
 275             fineMeshAddr.patchSchedule(),
 276             true
 277         )
 278     );
 279 }
 280
 281
 282 // ************************************************************************* //