src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/pairGAMGAgglomeration/pairGAMGAgglomerate.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 "pairGAMGAgglomeration.H"
  28 #include "lduAddressing.H"
  29
  30 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  31
  32 Foam::tmp<Foam::labelField> Foam::pairGAMGAgglomeration::agglomerate
  33 (
  34     label& nCoarseCells,
  35     const lduAddressing& fineMatrixAddressing,
  36     const scalarField& faceWeights
  37 )
  38 {
  39     const label nFineCells = fineMatrixAddressing.size();
  40
  41     const unallocLabelList& upperAddr = fineMatrixAddressing.upperAddr();
  42     const unallocLabelList& lowerAddr = fineMatrixAddressing.lowerAddr();
  43
  44     // For each cell calculate faces
  45     labelList cellFaces(upperAddr.size() + lowerAddr.size());
  46     labelList cellFaceOffsets(nFineCells + 1);
  47
  48     // memory management
  49     {
  50         labelList nNbrs(nFineCells, 0);
  51
  52         forAll (upperAddr, facei)
  53         {
  54             nNbrs[upperAddr[facei]]++;
  55         }
  56
  57         forAll (lowerAddr, facei)
  58         {
  59             nNbrs[lowerAddr[facei]]++;
  60         }
  61
  62         cellFaceOffsets[0] = 0;
  63         forAll (nNbrs, celli)
  64         {
  65             cellFaceOffsets[celli+1] = cellFaceOffsets[celli] + nNbrs[celli];
  66         }
  67
  68         // reset the whole list to use as counter
  69         nNbrs = 0;
  70
  71         forAll (upperAddr, facei)
  72         {
  73             cellFaces
  74             [
  75                 cellFaceOffsets[upperAddr[facei]] + nNbrs[upperAddr[facei]]
  76             ] = facei;
  77
  78             nNbrs[upperAddr[facei]]++;
  79         }
  80
  81         forAll (lowerAddr, facei)
  82         {
  83             cellFaces
  84             [
  85                 cellFaceOffsets[lowerAddr[facei]] + nNbrs[lowerAddr[facei]]
  86             ] = facei;
  87
  88             nNbrs[lowerAddr[facei]]++;
  89         }
  90     }
  91
  92
  93     // go through the faces and create clusters
  94
  95     tmp<labelField> tcoarseCellMap(new labelField(nFineCells, -1));
  96     labelField& coarseCellMap = tcoarseCellMap();
  97
  98     nCoarseCells = 0;
  99
 100     for (label celli=0; celli<nFineCells; celli++)
 101     {
 102         if (coarseCellMap[celli] < 0)
 103         {
 104             label matchFaceNo = -1;
 105             scalar maxFaceWeight = -GREAT;
 106
 107             // check faces to find ungrouped neighbour with largest face weight
 108             for
 109             (
 110                 label faceOs=cellFaceOffsets[celli];
 111                 faceOs<cellFaceOffsets[celli+1];
 112                 faceOs++
 113             )
 114             {
 115                 label facei = cellFaces[faceOs];
 116
 117                 // I don't know whether the current cell is owner or neighbour.
 118                 // Therefore I'll check both sides
 119                 if
 120                 (
 121                     coarseCellMap[upperAddr[facei]] < 0
 122                  && coarseCellMap[lowerAddr[facei]] < 0
 123                  && faceWeights[facei] > maxFaceWeight
 124                 )
 125                 {
 126                     // Match found. Pick up all the necessary data
 127                     matchFaceNo = facei;
 128                     maxFaceWeight = faceWeights[facei];
 129                 }
 130             }
 131
 132             if (matchFaceNo >= 0)
 133             {
 134                 // Make a new group
 135                 coarseCellMap[upperAddr[matchFaceNo]] = nCoarseCells;
 136                 coarseCellMap[lowerAddr[matchFaceNo]] = nCoarseCells;
 137                 nCoarseCells++;
 138             }
 139             else
 140             {
 141                 // No match. Find the best neighbouring cluster and
 142                 // put the cell there
 143                 label clusterMatchFaceNo = -1;
 144                 scalar clusterMaxFaceCoeff = -GREAT;
 145
 146                 for
 147                 (
 148                     label faceOs=cellFaceOffsets[celli];
 149                     faceOs<cellFaceOffsets[celli+1];
 150                     faceOs++
 151                 )
 152                 {
 153                     label facei = cellFaces[faceOs];
 154
 155                     if (faceWeights[facei] > clusterMaxFaceCoeff)
 156                     {
 157                         clusterMatchFaceNo = facei;
 158                         clusterMaxFaceCoeff = faceWeights[facei];
 159                     }
 160                 }
 161
 162                 if (clusterMatchFaceNo >= 0)
 163                 {
 164                     // Add the cell to the best cluster
 165                     coarseCellMap[celli] = max
 166                     (
 167                         coarseCellMap[upperAddr[clusterMatchFaceNo]],
 168                         coarseCellMap[lowerAddr[clusterMatchFaceNo]]
 169                     );
 170                 }
 171             }
 172         }
 173     }
 174
 175
 176     // Check that all cells are part of clusters,
 177     // if not create single-cell "clusters" for each
 178     for (label celli=0; celli<nFineCells; celli++)
 179     {
 180         if (coarseCellMap[celli] < 0)
 181         {
 182             coarseCellMap[celli] = nCoarseCells;
 183             nCoarseCells++;
 184         }
 185     }
 186
 187     // Reverese the map ordering to improve the next level of agglomeration
 188     // (doesn't always help and is sometimes detrimental)
 189     nCoarseCells--;
 190
 191     forAll (coarseCellMap, celli)
 192     {
 193         coarseCellMap[celli] = nCoarseCells - coarseCellMap[celli];
 194     }
 195
 196     nCoarseCells++;
 197
 198     return tcoarseCellMap;
 199 }
 200
 201
 202 void Foam::pairGAMGAgglomeration::agglomerate
 203 (
 204     const lduMesh& mesh,
 205     const scalarField& faceWeights
 206 )
 207 {
 208     // Get the finest-level interfaces from the mesh
 209     interfaceLevels_.set
 210     (
 211         0,
 212         new lduInterfacePtrsList(mesh.interfaces())
 213     );
 214
 215     // Start geometric agglomeration from the given faceWeights
 216     scalarField* faceWeightsPtr = const_cast<scalarField*>(&faceWeights);
 217
 218     // Agglomerate until the required number of cells in the coarsest level
 219     // is reached
 220
 221     label nPairLevels = 0;
 222     label nCreatedLevels = 0;
 223
 224     while (nCreatedLevels < maxLevels_ - 1)
 225     {
 226         label nCoarseCells = -1;
 227
 228         tmp<labelField> finalAgglomPtr = agglomerate
 229         (
 230             nCoarseCells,
 231             meshLevel(nCreatedLevels).lduAddr(),
 232             *faceWeightsPtr
 233         );
 234
 235         if (continueAgglomerating(nCoarseCells))
 236         {
 237             nCells_[nCreatedLevels] = nCoarseCells;
 238             restrictAddressing_.set(nCreatedLevels, finalAgglomPtr);
 239         }
 240         else
 241         {
 242             break;
 243         }
 244
 245         agglomerateLduAddressing(nCreatedLevels);
 246
 247         // Agglomerate the faceWeights field for the next level
 248         {
 249             scalarField* aggFaceWeightsPtr
 250             (
 251                 new scalarField
 252                 (
 253                     meshLevels_[nCreatedLevels].upperAddr().size(),
 254                     0.0
 255                 )
 256             );
 257
 258             restrictFaceField
 259             (
 260                 *aggFaceWeightsPtr,
 261                 *faceWeightsPtr,
 262                 nCreatedLevels
 263             );
 264
 265             if (nCreatedLevels)
 266             {
 267                 delete faceWeightsPtr;
 268             }
 269
 270             faceWeightsPtr = aggFaceWeightsPtr;
 271         }
 272
 273         if (nPairLevels % mergeLevels_)
 274         {
 275             combineLevels(nCreatedLevels);
 276         }
 277         else
 278         {
 279             nCreatedLevels++;
 280         }
 281
 282         nPairLevels++;
 283     }
 284
 285     // Shrink the storage of the levels to those created
 286     compactLevels(nCreatedLevels);
 287
 288     // Delete temporary geometry storage
 289     if (nCreatedLevels)
 290     {
 291         delete faceWeightsPtr;
 292     }
 293 }
 294
 295
 296 // ************************************************************************* //