src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixOperations.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 Description
  26     lduMatrix member operations.
  27
  28 \*---------------------------------------------------------------------------*/
  29
  30 #include "lduMatrix.H"
  31
  32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  33
  34 void Foam::lduMatrix::sumDiag()
  35 {
  36     const scalarField& Lower = const_cast<const lduMatrix&>(*this).lower();
  37     const scalarField& Upper = const_cast<const lduMatrix&>(*this).upper();
  38     scalarField& Diag = diag();
  39
  40     const unallocLabelList& l = lduAddr().lowerAddr();
  41     const unallocLabelList& u = lduAddr().upperAddr();
  42
  43     for (register label face=0; face<l.size(); face++)
  44     {
  45         Diag[l[face]] += Lower[face];
  46         Diag[u[face]] += Upper[face];
  47     }
  48 }
  49
  50
  51 void Foam::lduMatrix::negSumDiag()
  52 {
  53     const scalarField& Lower = const_cast<const lduMatrix&>(*this).lower();
  54     const scalarField& Upper = const_cast<const lduMatrix&>(*this).upper();
  55     scalarField& Diag = diag();
  56
  57     const unallocLabelList& l = lduAddr().lowerAddr();
  58     const unallocLabelList& u = lduAddr().upperAddr();
  59
  60     for (register label face=0; face<l.size(); face++)
  61     {
  62         Diag[l[face]] -= Lower[face];
  63         Diag[u[face]] -= Upper[face];
  64     }
  65 }
  66
  67
  68 void Foam::lduMatrix::sumMagOffDiag
  69 (
  70     scalarField& sumOff
  71 ) const
  72 {
  73     const scalarField& Lower = const_cast<const lduMatrix&>(*this).lower();
  74     const scalarField& Upper = const_cast<const lduMatrix&>(*this).upper();
  75
  76     const unallocLabelList& l = lduAddr().lowerAddr();
  77     const unallocLabelList& u = lduAddr().upperAddr();
  78
  79     for (register label face = 0; face < l.size(); face++)
  80     {
  81         sumOff[u[face]] += mag(Lower[face]);
  82         sumOff[l[face]] += mag(Upper[face]);
  83     }
  84 }
  85
  86
  87 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  88
  89 void Foam::lduMatrix::operator=(const lduMatrix& A)
  90 {
  91     if (this == &A)
  92     {
  93         FatalError
  94             << "lduMatrix::operator=(const lduMatrix&) : "
  95             << "attempted assignment to self"
  96             << abort(FatalError);
  97     }
  98
  99     if (A.lowerPtr_)
 100     {
 101         lower() = A.lower();
 102     }
 103     else if (lowerPtr_)
 104     {
 105         delete lowerPtr_;
 106         lowerPtr_ = NULL;
 107     }
 108
 109     if (A.upperPtr_)
 110     {
 111         upper() = A.upper();
 112     }
 113     else if (upperPtr_)
 114     {
 115         delete upperPtr_;
 116         upperPtr_ = NULL;
 117     }
 118
 119     if (A.diagPtr_)
 120     {
 121         diag() = A.diag();
 122     }
 123 }
 124
 125
 126 void Foam::lduMatrix::negate()
 127 {
 128     if (lowerPtr_)
 129     {
 130         lowerPtr_->negate();
 131     }
 132
 133     if (upperPtr_)
 134     {
 135         upperPtr_->negate();
 136     }
 137
 138     if (diagPtr_)
 139     {
 140         diagPtr_->negate();
 141     }
 142 }
 143
 144
 145 void Foam::lduMatrix::operator+=(const lduMatrix& A)
 146 {
 147     if (A.diagPtr_)
 148     {
 149         diag() += A.diag();
 150     }
 151
 152     if (symmetric() && A.symmetric())
 153     {
 154         upper() += A.upper();
 155     }
 156     else if (symmetric() && A.asymmetric())
 157     {
 158         if (upperPtr_)
 159         {
 160             lower();
 161         }
 162         else
 163         {
 164             upper();
 165         }
 166
 167         upper() += A.upper();
 168         lower() += A.lower();
 169     }
 170     else if (asymmetric() && A.symmetric())
 171     {
 172         if (A.upperPtr_)
 173         {
 174             lower() += A.upper();
 175             upper() += A.upper();
 176         }
 177         else
 178         {
 179             lower() += A.lower();
 180             upper() += A.lower();
 181         }
 182
 183     }
 184     else if (asymmetric() && A.asymmetric())
 185     {
 186         lower() += A.lower();
 187         upper() += A.upper();
 188     }
 189     else if (diagonal())
 190     {
 191         if (A.upperPtr_)
 192         {
 193             upper() = A.upper();
 194         }
 195
 196         if (A.lowerPtr_)
 197         {
 198             lower() = A.lower();
 199         }
 200     }
 201     else if (A.diagonal())
 202     {
 203     }
 204     else
 205     {
 206         FatalErrorIn("lduMatrix::operator+=(const lduMatrix& A)")
 207             << "Unknown matrix type combination"
 208             << abort(FatalError);
 209     }
 210 }
 211
 212
 213 void Foam::lduMatrix::operator-=(const lduMatrix& A)
 214 {
 215     if (A.diagPtr_)
 216     {
 217         diag() -= A.diag();
 218     }
 219
 220     if (symmetric() && A.symmetric())
 221     {
 222         upper() -= A.upper();
 223     }
 224     else if (symmetric() && A.asymmetric())
 225     {
 226         if (upperPtr_)
 227         {
 228             lower();
 229         }
 230         else
 231         {
 232             upper();
 233         }
 234
 235         upper() -= A.upper();
 236         lower() -= A.lower();
 237     }
 238     else if (asymmetric() && A.symmetric())
 239     {
 240         if (A.upperPtr_)
 241         {
 242             lower() -= A.upper();
 243             upper() -= A.upper();
 244         }
 245         else
 246         {
 247             lower() -= A.lower();
 248             upper() -= A.lower();
 249         }
 250
 251     }
 252     else if (asymmetric() && A.asymmetric())
 253     {
 254         lower() -= A.lower();
 255         upper() -= A.upper();
 256     }
 257     else if (diagonal())
 258     {
 259         if (A.upperPtr_)
 260         {
 261             upper() = -A.upper();
 262         }
 263
 264         if (A.lowerPtr_)
 265         {
 266             lower() = -A.lower();
 267         }
 268     }
 269     else if (A.diagonal())
 270     {
 271     }
 272     else
 273     {
 274         FatalErrorIn("lduMatrix::operator-=(const lduMatrix& A)")
 275             << "Unknown matrix type combination"
 276             << abort(FatalError);
 277     }
 278 }
 279
 280
 281 void Foam::lduMatrix::operator*=(const scalarField& sf)
 282 {
 283     if (diagPtr_)
 284     {
 285         *diagPtr_ *= sf;
 286     }
 287
 288     if (upperPtr_)
 289     {
 290         scalarField& upper = *upperPtr_;
 291
 292         const unallocLabelList& l = lduAddr().lowerAddr();
 293
 294         for (register label face=0; face<upper.size(); face++)
 295         {
 296             upper[face] *= sf[l[face]];
 297         }
 298     }
 299
 300     if (lowerPtr_)
 301     {
 302         scalarField& lower = *lowerPtr_;
 303
 304         const unallocLabelList& u = lduAddr().upperAddr();
 305
 306         for (register label face=0; face<lower.size(); face++)
 307         {
 308             lower[face] *= sf[u[face]];
 309         }
 310     }
 311 }
 312
 313
 314 void Foam::lduMatrix::operator*=(scalar s)
 315 {
 316     if (diagPtr_)
 317     {
 318         *diagPtr_ *= s;
 319     }
 320
 321     if (upperPtr_)
 322     {
 323         *upperPtr_ *= s;
 324     }
 325
 326     if (lowerPtr_)
 327     {
 328         *lowerPtr_ *= s;
 329     }
 330 }
 331
 332
 333 Foam::tmp<Foam::scalarField > Foam::lduMatrix::H1() const
 334 {
 335     tmp<scalarField > tH1
 336     (
 337         new scalarField(lduAddr().size(), 0.0)
 338     );
 339
 340     if (lowerPtr_ || upperPtr_)
 341     {
 342         scalarField& H1_ = tH1();
 343
 344         scalar* __restrict__ H1Ptr = H1_.begin();
 345
 346         const label* __restrict__ uPtr = lduAddr().upperAddr().begin();
 347         const label* __restrict__ lPtr = lduAddr().lowerAddr().begin();
 348
 349         const scalar* __restrict__ lowerPtr = lower().begin();
 350         const scalar* __restrict__ upperPtr = upper().begin();
 351
 352         register const label nFaces = upper().size();
 353
 354         for (register label face=0; face<nFaces; face++)
 355         {
 356             H1Ptr[uPtr[face]] -= lowerPtr[face];
 357             H1Ptr[lPtr[face]] -= upperPtr[face];
 358         }
 359     }
 360
 361     return tH1;
 362 }
 363
 364
 365 // ************************************************************************* //