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