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