src/foam/matrices/lduMatrix/solvers/GAMG/GAMGSolverAgglomerateMatrix.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 "GAMGSolver.H"
  27 #include "GAMGInterfaceField.H"
  28
  29 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  30
  31 void Foam::GAMGSolver::agglomerateMatrix(const label fineLevelIndex)
  32 {
  33     // Get fine matrix
  34     const lduMatrix& fineMatrix = matrixLevel(fineLevelIndex);
  35
  36     // Set the coarse level matrix
  37     matrixLevels_.set
  38     (
  39         fineLevelIndex,
  40         new lduMatrix(agglomeration_.meshLevel(fineLevelIndex + 1))
  41     );
  42     lduMatrix& coarseMatrix = matrixLevels_[fineLevelIndex];
  43
  44     // Get face restriction map for current level
  45     const labelList& faceRestrictAddr =
  46         agglomeration_.faceRestrictAddressing(fineLevelIndex);
  47
  48     // Coarse matrix diagonal initialised by restricting the fine mesh diagonal
  49     scalarField& coarseDiag = coarseMatrix.diag();
  50     agglomeration_.restrictField
  51     (
  52         coarseDiag,
  53         fineMatrix.diag(),
  54         fineLevelIndex
  55     );
  56
  57     // Get reference to fine-level interfaces
  58     const lduInterfaceFieldPtrsList& fineInterfaces =
  59         interfaceLevel(fineLevelIndex);
  60
  61     // Get reference to fine-level boundary coefficients
  62     const FieldField<Field, scalar>& fineInterfaceBouCoeffs =
  63         coupleBouCoeffsLevel(fineLevelIndex);
  64
  65     // Get reference to fine-level internal coefficients
  66     const FieldField<Field, scalar>& fineInterfaceIntCoeffs =
  67         coupleIntCoeffsLevel(fineLevelIndex);
  68
  69
  70     // Create coarse-level interfaces
  71     interfaceLevels_.set
  72     (
  73         fineLevelIndex,
  74         new lduInterfaceFieldPtrsList(fineInterfaces.size())
  75     );
  76
  77     lduInterfaceFieldPtrsList& coarseInterfaces =
  78         interfaceLevels_[fineLevelIndex];
  79
  80     // Set coarse-level boundary coefficients
  81     coupleLevelsBouCoeffs_.set
  82     (
  83         fineLevelIndex,
  84         new FieldField<Field, scalar>(fineInterfaces.size())
  85     );
  86     FieldField<Field, scalar>& coarseInterfaceBouCoeffs =
  87         coupleLevelsBouCoeffs_[fineLevelIndex];
  88
  89     // Set coarse-level internal coefficients
  90     coupleLevelsIntCoeffs_.set
  91     (
  92         fineLevelIndex,
  93         new FieldField<Field, scalar>(fineInterfaces.size())
  94     );
  95     FieldField<Field, scalar>& coarseInterfaceIntCoeffs =
  96         coupleLevelsIntCoeffs_[fineLevelIndex];
  97
  98     // Add the coarse level
  99     forAll (fineInterfaces, inti)
 100     {
 101         if (fineInterfaces.set(inti))
 102         {
 103             const GAMGInterface& coarseInterface =
 104                 refCast<const GAMGInterface>
 105                 (
 106                     agglomeration_.interfaceLevel(fineLevelIndex + 1)[inti]
 107                 );
 108
 109             coarseInterfaces.set
 110             (
 111                 inti,
 112                 GAMGInterfaceField::New
 113                 (
 114                     coarseInterface,
 115                     fineInterfaces[inti]
 116                 ).ptr()
 117             );
 118
 119             coarseInterfaceBouCoeffs.set
 120             (
 121                 inti,
 122                 coarseInterface.agglomerateCoeffs(fineInterfaceBouCoeffs[inti])
 123             );
 124
 125             coarseInterfaceIntCoeffs.set
 126             (
 127                 inti,
 128                 coarseInterface.agglomerateCoeffs(fineInterfaceIntCoeffs[inti])
 129             );
 130         }
 131     }
 132
 133
 134     // Check if matrix is assymetric and if so agglomerate both upper and lower
 135     // coefficients ...
 136     if (fineMatrix.hasLower())
 137     {
 138         // Get off-diagonal matrix coefficients
 139         const scalarField& fineUpper = fineMatrix.upper();
 140         const scalarField& fineLower = fineMatrix.lower();
 141
 142         // Coarse matrix upper coefficients
 143         scalarField& coarseUpper = coarseMatrix.upper();
 144         scalarField& coarseLower = coarseMatrix.lower();
 145
 146         const labelList& restrictAddr =
 147             agglomeration_.restrictAddressing(fineLevelIndex);
 148
 149         const unallocLabelList& l = fineMatrix.lduAddr().lowerAddr();
 150         const unallocLabelList& cl = coarseMatrix.lduAddr().lowerAddr();
 151         const unallocLabelList& cu = coarseMatrix.lduAddr().upperAddr();
 152
 153         forAll(faceRestrictAddr, fineFacei)
 154         {
 155             label cFace = faceRestrictAddr[fineFacei];
 156
 157             if (cFace >= 0)
 158             {
 159                 // Check the orientation of the fine-face relative to the
 160                 // coarse face it is being agglomerated into
 161                 if (cl[cFace] == restrictAddr[l[fineFacei]])
 162                 {
 163                     coarseUpper[cFace] += fineUpper[fineFacei];
 164                     coarseLower[cFace] += fineLower[fineFacei];
 165                 }
 166                 else if (cu[cFace] == restrictAddr[l[fineFacei]])
 167                 {
 168                     coarseUpper[cFace] += fineLower[fineFacei];
 169                     coarseLower[cFace] += fineUpper[fineFacei];
 170                 }
 171                 else
 172                 {
 173                     FatalErrorIn
 174                     (
 175                         "GAMGSolver::agglomerateMatrix(const label)"
 176                     )   << "Inconsistent addressing between "
 177                            "fine and coarse grids"
 178                         << exit(FatalError);
 179                 }
 180             }
 181             else
 182             {
 183                 // Add the fine face coefficients into the diagonal.
 184                 coarseDiag[-1 - cFace] +=
 185                     fineUpper[fineFacei] + fineLower[fineFacei];
 186             }
 187         }
 188     }
 189     else // ... Otherwise it is symmetric so agglomerate just the upper
 190     {
 191         // Get off-diagonal matrix coefficients
 192         const scalarField& fineUpper = fineMatrix.upper();
 193
 194         // Coarse matrix upper coefficients
 195         scalarField& coarseUpper = coarseMatrix.upper();
 196
 197         forAll(faceRestrictAddr, fineFacei)
 198         {
 199             label cFace = faceRestrictAddr[fineFacei];
 200
 201             if (cFace >= 0)
 202             {
 203                 coarseUpper[cFace] += fineUpper[fineFacei];
 204             }
 205             else
 206             {
 207                 // Add the fine face coefficient into the diagonal.
 208                 coarseDiag[-1 - cFace] += 2*fineUpper[fineFacei];
 209             }
 210         }
 211     }
 212 }
 213
 214
 215 // ************************************************************************* //