src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverAgglomerateMatrix.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-2010 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
  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 "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 finer mesh diagonal
  49     scalarField& coarseDiag = coarseMatrix.diag();
  50     agglomeration_.restrictField(coarseDiag, fineMatrix.diag(), fineLevelIndex);
  51
  52     // Get reference to fine-level interfaces
  53     const lduInterfaceFieldPtrsList& fineInterfaces =
  54         interfaceLevel(fineLevelIndex);
  55
  56     // Get reference to fine-level boundary coefficients
  57     const FieldField<Field, scalar>& fineInterfaceBouCoeffs =
  58         interfaceBouCoeffsLevel(fineLevelIndex);
  59
  60     // Get reference to fine-level internal coefficients
  61     const FieldField<Field, scalar>& fineInterfaceIntCoeffs =
  62         interfaceIntCoeffsLevel(fineLevelIndex);
  63
  64
  65     // Create coarse-level interfaces
  66     interfaceLevels_.set
  67     (
  68         fineLevelIndex,
  69         new lduInterfaceFieldPtrsList(fineInterfaces.size())
  70     );
  71
  72     lduInterfaceFieldPtrsList& coarseInterfaces =
  73         interfaceLevels_[fineLevelIndex];
  74
  75     // Set coarse-level boundary coefficients
  76     interfaceLevelsBouCoeffs_.set
  77     (
  78         fineLevelIndex,
  79         new FieldField<Field, scalar>(fineInterfaces.size())
  80     );
  81     FieldField<Field, scalar>& coarseInterfaceBouCoeffs =
  82         interfaceLevelsBouCoeffs_[fineLevelIndex];
  83
  84     // Set coarse-level internal coefficients
  85     interfaceLevelsIntCoeffs_.set
  86     (
  87         fineLevelIndex,
  88         new FieldField<Field, scalar>(fineInterfaces.size())
  89     );
  90     FieldField<Field, scalar>& coarseInterfaceIntCoeffs =
  91         interfaceLevelsIntCoeffs_[fineLevelIndex];
  92
  93     // Add the coarse level
  94     forAll(fineInterfaces, inti)
  95     {
  96         if (fineInterfaces.set(inti))
  97         {
  98             const GAMGInterface& coarseInterface =
  99                 refCast<const GAMGInterface>
 100                 (
 101                     agglomeration_.interfaceLevel(fineLevelIndex + 1)[inti]
 102                 );
 103
 104             coarseInterfaces.set
 105             (
 106                 inti,
 107                 GAMGInterfaceField::New
 108                 (
 109                     coarseInterface,
 110                     fineInterfaces[inti]
 111                 ).ptr()
 112             );
 113
 114             coarseInterfaceBouCoeffs.set
 115             (
 116                 inti,
 117                 coarseInterface.agglomerateCoeffs(fineInterfaceBouCoeffs[inti])
 118             );
 119
 120             coarseInterfaceIntCoeffs.set
 121             (
 122                 inti,
 123                 coarseInterface.agglomerateCoeffs(fineInterfaceIntCoeffs[inti])
 124             );
 125         }
 126     }
 127
 128
 129     // Check if matrix is assymetric and if so agglomerate both upper and lower
 130     // coefficients ...
 131     if (fineMatrix.hasLower())
 132     {
 133         // Get off-diagonal matrix coefficients
 134         const scalarField& fineUpper = fineMatrix.upper();
 135         const scalarField& fineLower = fineMatrix.lower();
 136
 137         // Coarse matrix upper coefficients
 138         scalarField& coarseUpper = coarseMatrix.upper();
 139         scalarField& coarseLower = coarseMatrix.lower();
 140
 141         const labelList& restrictAddr =
 142             agglomeration_.restrictAddressing(fineLevelIndex);
 143
 144         const labelUList& l = fineMatrix.lduAddr().lowerAddr();
 145         const labelUList& cl = coarseMatrix.lduAddr().lowerAddr();
 146         const labelUList& cu = coarseMatrix.lduAddr().upperAddr();
 147
 148         forAll(faceRestrictAddr, fineFacei)
 149         {
 150             label cFace = faceRestrictAddr[fineFacei];
 151
 152             if (cFace >= 0)
 153             {
 154                 // Check the orientation of the fine-face relative to the
 155                 // coarse face it is being agglomerated into
 156                 if (cl[cFace] == restrictAddr[l[fineFacei]])
 157                 {
 158                     coarseUpper[cFace] += fineUpper[fineFacei];
 159                     coarseLower[cFace] += fineLower[fineFacei];
 160                 }
 161                 else if (cu[cFace] == restrictAddr[l[fineFacei]])
 162                 {
 163                     coarseUpper[cFace] += fineLower[fineFacei];
 164                     coarseLower[cFace] += fineUpper[fineFacei];
 165                 }
 166                 else
 167                 {
 168                     FatalErrorIn
 169                     (
 170                         "GAMGSolver::agglomerateMatrix(const label)"
 171                     )   << "Inconsistent addressing between "
 172                            "fine and coarse grids"
 173                         << exit(FatalError);
 174                 }
 175             }
 176             else
 177             {
 178                 // Add the fine face coefficients into the diagonal.
 179                 coarseDiag[-1 - cFace] +=
 180                     fineUpper[fineFacei] + fineLower[fineFacei];
 181             }
 182         }
 183     }
 184     else // ... Otherwise it is symmetric so agglomerate just the upper
 185     {
 186         // Get off-diagonal matrix coefficients
 187         const scalarField& fineUpper = fineMatrix.upper();
 188
 189         // Coarse matrix upper coefficients
 190         scalarField& coarseUpper = coarseMatrix.upper();
 191
 192         forAll(faceRestrictAddr, fineFacei)
 193         {
 194             label cFace = faceRestrictAddr[fineFacei];
 195
 196             if (cFace >= 0)
 197             {
 198                 coarseUpper[cFace] += fineUpper[fineFacei];
 199             }
 200             else
 201             {
 202                 // Add the fine face coefficient into the diagonal.
 203                 coarseDiag[-1 - cFace] += 2*fineUpper[fineFacei];
 204             }
 205         }
 206     }
 207 }
 208
 209
 210 // ************************************************************************* //