src/lduSolvers/crMatrix/crMatrix.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-6 H. Jasak All rights reserved
   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 Class
  26     crMatrix
  27
  28 Description
  29     Sparse matrix in compressed row format
  30
  31 Author
  32     Hrvoje Jasak, Wikki Ltd.  All rights reserved
  33
  34 \*----------------------------------------------------------------------------*/
  35
  36 #include "crMatrix.H"
  37
  38 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  39
  40 // Construct from addressing
  41 Foam::crMatrix::crMatrix(const crAddressing& addr)
  42 :
  43     refCount(),
  44     crAddr_(addr),
  45     coeffs_(addr.nEntries())
  46 {}
  47
  48
  49 // Construct given size.  Column and coefficients set later
  50 Foam::crMatrix::crMatrix
  51 (
  52     const label nRows,
  53     const label nCols,
  54     const labelList& count
  55 )
  56 :
  57     refCount(),
  58     crAddr_(nRows, nCols, count),
  59     coeffs_(crAddr_.nEntries())
  60 {}
  61
  62
  63 // Construct from components of addressing
  64 Foam::crMatrix::crMatrix
  65 (
  66     const label nRows,
  67     const label nCols,
  68     const labelList& row,
  69     const labelList& col
  70 )
  71 :
  72     refCount(),
  73     crAddr_(nRows, nCols, row, col),
  74     coeffs_(crAddr_.nEntries())
  75 {}
  76
  77 // Construct as copy
  78 Foam::crMatrix::crMatrix(const crMatrix& m)
  79 :
  80     refCount(),
  81     crAddr_(m.crAddr_),
  82     coeffs_(m.coeffs_)
  83 {}
  84
  85
  86 // Construct as copy of tmp<crMatrix> deleting argument
  87 Foam::crMatrix::crMatrix(const tmp<crMatrix>& tm)
  88 :
  89     refCount(),
  90     crAddr_(tm().crAddr()),
  91     coeffs_(tm().coeffs())
  92 {
  93     tm.clear();
  94 }
  95
  96
  97 // Construct from Istream
  98 Foam::crMatrix::crMatrix(Istream& is)
  99 :
 100     refCount(),
 101     crAddr_(is),
 102     coeffs_(is)
 103 {}
 104
 105
 106 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 107
 108 // Return transpose
 109 Foam::tmp<Foam::crMatrix> Foam::crMatrix::T() const
 110 {
 111     // My addressing
 112     const label myNRows = crAddr().nRows();
 113     const labelList& myRow = crAddr().row();
 114     const labelList& myCol = crAddr().col();
 115     const scalarField& myCoeffs = coeffs();
 116
 117     // Create transpose
 118     tmp<crMatrix> ttranspose(new crMatrix(crAddr().T()));
 119     crMatrix& transpose = ttranspose();
 120
 121     scalarField& tCoeffs = transpose.coeffs();
 122
 123     // Reset coefficients to zero to treat degenerate matrices
 124     tCoeffs = 0;
 125
 126     // Transpose addressing
 127     const labelList& tRow = transpose.crAddr().row();
 128     const labelList& tCol = transpose.crAddr().col();
 129
 130     // Set transpose coefficients
 131
 132     label ja;
 133
 134     for (label ia = 0; ia < myNRows; ia++)
 135     {
 136         for (label jpa = myRow[ia]; jpa < myRow[ia + 1]; jpa++)
 137         {
 138             ja = myCol[jpa];
 139
 140             for (label jpb = tRow[ja]; jpb < tRow[ja + 1]; jpb++)
 141             {
 142                 if (tCol[jpb] == ia)
 143                 {
 144                     tCoeffs[jpb] = myCoeffs[jpa];
 145                     break;
 146                 }
 147             }
 148         }
 149     }
 150
 151     return ttranspose;
 152 }
 153
 154 // Calculate b += A*x
 155 void Foam::crMatrix::dotPlus(scalarField& b, const scalarField& x) const
 156 {
 157     const labelList& row = crAddr_.row();
 158     const labelList& col = crAddr_.col();
 159
 160     forAll (b, i)
 161     {
 162         scalar& bi = b[i];
 163
 164         for (label ip = row[i]; ip < row[i + 1];  ip++)
 165         {
 166             bi += coeffs_[ip]*x[col[ip]];
 167         }
 168     }
 169 }
 170
 171
 172 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 173
 174 void Foam::crMatrix::operator=(const crMatrix& rhs)
 175 {
 176     // Check for assignment to self
 177     if (this == &rhs)
 178     {
 179         FatalErrorIn("Foam::crMatrix::operator=(const Foam::crMatrix&)")
 180             << "Attempted assignment to self"
 181             << abort(FatalError);
 182     }
 183
 184     crAddr_ = rhs.crAddr_;
 185     coeffs_ = rhs.coeffs_;
 186 }
 187
 188
 189 void Foam::crMatrix::operator=(const tmp<crMatrix>& trhs)
 190 {
 191     operator=(trhs());
 192     trhs.clear();
 193 }
 194
 195
 196 // * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
 197
 198
 199 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 200
 201 Foam::Ostream& Foam::operator<<(Ostream& os, const crMatrix& m)
 202 {
 203     os << m.crAddr_ << m.coeffs_ << endl;
 204
 205     return os;
 206 }
 207
 208
 209 // ************************************************************************* //