src/blockMatrix/BlockLduMatrix/BlockLduMatrixATmul.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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
  24
  25 Description
  26     Vector-matrix multiplication operations for a block matrix
  27
  28 \*---------------------------------------------------------------------------*/
  29
  30 #include "BlockLduMatrix.H"
  31
  32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  33
  34 template<class Type>
  35 void Foam::BlockLduMatrix<Type>::Amul
  36 (
  37     TypeField& Ax,
  38     const TypeField& x
  39 ) const
  40 {
  41     // Initialise the update of coupled interfaces
  42     initInterfaces(Ax, x);
  43
  44     AmulCore(Ax, x);
  45
  46     // Update coupled interfaces
  47     updateInterfaces(coupleUpper_, Ax, x);
  48 }
  49
  50
  51 template<class Type>
  52 void Foam::BlockLduMatrix<Type>::AmulCore
  53 (
  54     TypeField& Ax,
  55     const TypeField& x
  56 ) const
  57 {
  58     typedef typename TypeCoeffField::scalarTypeField scalarTypeField;
  59     typedef typename TypeCoeffField::linearTypeField linearTypeField;
  60     typedef typename TypeCoeffField::squareTypeField squareTypeField;
  61
  62     const unallocLabelList& u = lduAddr().upperAddr();
  63     const unallocLabelList& l = lduAddr().lowerAddr();
  64
  65     const TypeCoeffField& Diag = this->diag();
  66     const TypeCoeffField& Upper = this->upper();
  67
  68     // Diagonal multiplication, no indirection
  69     multiply(Ax, Diag, x);
  70
  71     // Create multiplication function object
  72     typename BlockCoeff<Type>::multiply mult;
  73
  74     // Lower multiplication
  75
  76     if (symmetric())
  77     {
  78         if (Upper.activeType() == blockCoeffBase::SCALAR)
  79         {
  80             const scalarTypeField& activeUpper = Upper.asScalar();
  81
  82             for (register label coeffI = 0; coeffI < u.size(); coeffI++)
  83             {
  84                 Ax[u[coeffI]] += mult(activeUpper[coeffI], x[l[coeffI]]);
  85             }
  86         }
  87         else if (Upper.activeType() == blockCoeffBase::LINEAR)
  88         {
  89             const linearTypeField& activeUpper = Upper.asLinear();
  90
  91             for (register label coeffI = 0; coeffI < u.size(); coeffI++)
  92             {
  93                 Ax[u[coeffI]] += mult(activeUpper[coeffI], x[l[coeffI]]);
  94             }
  95         }
  96         else if (Upper.activeType() == blockCoeffBase::SQUARE)
  97         {
  98             const squareTypeField& activeUpper = Upper.asSquare();
  99
 100             for (register label coeffI = 0; coeffI < u.size(); coeffI++)
 101             {
 102                 // Use transpose upper coefficient
 103                 Ax[u[coeffI]] +=
 104                     mult(activeUpper[coeffI].T(), x[l[coeffI]]);
 105             }
 106         }
 107     }
 108     else // Asymmetric matrix
 109     {
 110         const TypeCoeffField& Lower = this->lower();
 111
 112         if (Lower.activeType() == blockCoeffBase::SCALAR)
 113         {
 114             const scalarTypeField& activeLower = Lower.asScalar();
 115
 116             for (register label coeffI = 0; coeffI < u.size(); coeffI++)
 117             {
 118                 Ax[u[coeffI]] += mult(activeLower[coeffI], x[l[coeffI]]);
 119             }
 120         }
 121         else if (Lower.activeType() == blockCoeffBase::LINEAR)
 122         {
 123             const linearTypeField& activeLower = Lower.asLinear();
 124
 125             for (register label coeffI = 0; coeffI < u.size(); coeffI++)
 126             {
 127                 Ax[u[coeffI]] += mult(activeLower[coeffI], x[l[coeffI]]);
 128             }
 129         }
 130         else if (Lower.activeType() == blockCoeffBase::SQUARE)
 131         {
 132             const squareTypeField& activeLower = Lower.asSquare();
 133
 134             for (register label coeffI = 0; coeffI < u.size(); coeffI++)
 135             {
 136                 Ax[u[coeffI]] += mult(activeLower[coeffI], x[l[coeffI]]);
 137             }
 138         }
 139     }
 140
 141
 142     // Upper multiplication
 143
 144     if (Upper.activeType() == blockCoeffBase::SCALAR)
 145     {
 146         const scalarTypeField& activeUpper = Upper.asScalar();
 147
 148         for (register label coeffI = 0; coeffI < u.size(); coeffI++)
 149         {
 150             Ax[l[coeffI]] += mult(activeUpper[coeffI], x[u[coeffI]]);
 151         }
 152     }
 153     else if (Upper.activeType() == blockCoeffBase::LINEAR)
 154     {
 155         const linearTypeField& activeUpper = Upper.asLinear();
 156
 157         for (register label coeffI = 0; coeffI < u.size(); coeffI++)
 158         {
 159             Ax[l[coeffI]] += mult(activeUpper[coeffI], x[u[coeffI]]);
 160         }
 161     }
 162     else if (Upper.activeType() == blockCoeffBase::SQUARE)
 163     {
 164         const squareTypeField& activeUpper = Upper.asSquare();
 165
 166         for (register label coeffI = 0; coeffI < u.size(); coeffI++)
 167         {
 168             Ax[l[coeffI]] += mult(activeUpper[coeffI], x[u[coeffI]]);
 169         }
 170     }
 171 }
 172
 173
 174 template<class Type>
 175 void Foam::BlockLduMatrix<Type>::Tmul
 176 (
 177     TypeField& Ax,
 178     const TypeField& x
 179 ) const
 180 {
 181     // Initialise the update of coupled interfaces
 182     initInterfaces(Ax, x);
 183
 184     TmulCore(Ax, x);
 185
 186     // Update coupled interfaces
 187     updateInterfaces(coupleLower_, Ax, x);
 188 }
 189
 190
 191 template<class Type>
 192 void Foam::BlockLduMatrix<Type>::TmulCore
 193 (
 194     TypeField& Tx,
 195     const TypeField& x
 196 ) const
 197 {
 198     typedef typename TypeCoeffField::scalarTypeField scalarTypeField;
 199     typedef typename TypeCoeffField::linearTypeField linearTypeField;
 200     typedef typename TypeCoeffField::squareTypeField squareTypeField;
 201
 202     const unallocLabelList& u = lduAddr().upperAddr();
 203     const unallocLabelList& l = lduAddr().lowerAddr();
 204
 205     const TypeCoeffField& Diag = this->diag();
 206     const TypeCoeffField& Upper = this->upper();
 207
 208     // Diagonal multiplication, no indirection
 209     multiply(Tx, Diag, x);
 210
 211     // Create multiplication function object
 212     typename BlockCoeff<Type>::multiply mult;
 213
 214     // Upper multiplication
 215
 216     if (Upper.activeType() == blockCoeffBase::SCALAR)
 217     {
 218         const scalarTypeField& activeUpper = Upper.asScalar();
 219
 220         for (register label coeffI = 0; coeffI < u.size(); coeffI++)
 221         {
 222             Tx[u[coeffI]] += mult(activeUpper[coeffI], x[l[coeffI]]);
 223         }
 224     }
 225     else if (Upper.activeType() == blockCoeffBase::LINEAR)
 226     {
 227         const linearTypeField& activeUpper = Upper.asLinear();
 228
 229         for (register label coeffI = 0; coeffI < u.size(); coeffI++)
 230         {
 231             Tx[u[coeffI]] += mult(activeUpper[coeffI], x[l[coeffI]]);
 232         }
 233     }
 234     else if (Upper.activeType() == blockCoeffBase::SQUARE)
 235     {
 236         const squareTypeField& activeUpper = Upper.asSquare();
 237
 238         for (register label coeffI = 0; coeffI < u.size(); coeffI++)
 239         {
 240             Tx[u[coeffI]] += mult(activeUpper[coeffI], x[l[coeffI]]);
 241         }
 242     }
 243
 244     // Lower multiplication
 245
 246     if (symmetric())
 247     {
 248         if (Upper.activeType() == blockCoeffBase::SCALAR)
 249         {
 250             const scalarTypeField& activeUpper = Upper.asScalar();
 251
 252             for (register label coeffI = 0; coeffI < u.size(); coeffI++)
 253             {
 254                 Tx[l[coeffI]] += mult(activeUpper[coeffI], x[u[coeffI]]);
 255             }
 256         }
 257         else if (Upper.activeType() == blockCoeffBase::LINEAR)
 258         {
 259             const linearTypeField& activeUpper = Upper.asLinear();
 260
 261             for (register label coeffI = 0; coeffI < u.size(); coeffI++)
 262             {
 263                 Tx[l[coeffI]] += mult(activeUpper[coeffI], x[u[coeffI]]);
 264             }
 265         }
 266         else if (Upper.activeType() == blockCoeffBase::SQUARE)
 267         {
 268             const squareTypeField& activeUpper = Upper.asSquare();
 269
 270             for (register label coeffI = 0; coeffI < u.size(); coeffI++)
 271             {
 272                 // Use transpose upper coefficient
 273                 Tx[l[coeffI]] +=
 274                     mult(activeUpper[coeffI].T(), x[u[coeffI]]);
 275             }
 276         }
 277     }
 278     else // Asymmetric matrix
 279     {
 280         const TypeCoeffField& Lower = this->lower();
 281
 282         if (Lower.activeType() == blockCoeffBase::SCALAR)
 283         {
 284             const scalarTypeField& activeLower = Lower.asScalar();
 285
 286             for (register label coeffI = 0; coeffI < u.size(); coeffI++)
 287             {
 288                 Tx[l[coeffI]] += mult(activeLower[coeffI], x[u[coeffI]]);
 289             }
 290         }
 291         else if (Lower.activeType() == blockCoeffBase::LINEAR)
 292         {
 293             const linearTypeField& activeLower = Lower.asLinear();
 294
 295             for (register label coeffI = 0; coeffI < u.size(); coeffI++)
 296             {
 297                 Tx[l[coeffI]] += mult(activeLower[coeffI], x[u[coeffI]]);
 298             }
 299         }
 300         else if (Lower.activeType() == blockCoeffBase::SQUARE)
 301         {
 302             const squareTypeField& activeLower = Lower.asSquare();
 303
 304             for (register label coeffI = 0; coeffI < u.size(); coeffI++)
 305             {
 306                 Tx[l[coeffI]] += mult(activeLower[coeffI], x[u[coeffI]]);
 307             }
 308         }
 309     }
 310 }
 311
 312
 313 template<class Type>
 314 void Foam::BlockLduMatrix<Type>::segregateB
 315 (
 316     TypeField& sMul,
 317     const TypeField& x
 318 ) const
 319 {
 320     typedef typename TypeCoeffField::linearType linearType;
 321     typedef typename TypeCoeffField::squareType squareType;
 322
 323     typedef typename TypeCoeffField::linearTypeField linearTypeField;
 324     typedef typename TypeCoeffField::squareTypeField squareTypeField;
 325
 326     const unallocLabelList& u = lduAddr().upperAddr();
 327     const unallocLabelList& l = lduAddr().lowerAddr();
 328
 329     // Diagonal multiplication
 330     if (thereIsDiag())
 331     {
 332         if (diag().activeType() == blockCoeffBase::SQUARE)
 333         {
 334             const squareTypeField& activeDiag = this->diag().asSquare();
 335             linearTypeField lf(activeDiag.size());
 336             squareTypeField sf(activeDiag.size());
 337
 338             // Expand and contract
 339             contractLinear(lf, activeDiag);
 340             expandLinear(sf, lf);
 341
 342             sMul -= (activeDiag - sf) & x;
 343         }
 344     }
 345
 346     // Lower multiplication
 347
 348     if (thereIsLower())
 349     {
 350         if (lower().activeType() == blockCoeffBase::SQUARE)
 351         {
 352             const squareTypeField& activeLower = this->lower().asSquare();
 353
 354             // Auxiliary variables used in expand/contract
 355             linearType lt;
 356             squareType st;
 357
 358             for (register label coeffI = 0; coeffI < u.size(); coeffI++)
 359             {
 360                 contractLinear(lt, activeLower[coeffI]);
 361                 expandLinear(st, lt);
 362
 363                 sMul[u[coeffI]] -= (activeLower[coeffI] - st) & x[l[coeffI]];
 364             }
 365         }
 366     }
 367
 368     // Upper multiplication
 369
 370     if (thereIsUpper())
 371     {
 372         if (upper().activeType() == blockCoeffBase::SQUARE)
 373         {
 374             const squareTypeField& activeUpper = this->upper().asSquare();
 375
 376             // Auxiliary variables used in expand/contract
 377             linearType lt;
 378             squareType st;
 379
 380             for (register label coeffI = 0; coeffI < u.size(); coeffI++)
 381             {
 382                 contractLinear(lt, activeUpper[coeffI]);
 383                 expandLinear(st, lt);
 384
 385                 sMul[l[coeffI]] -= (activeUpper[coeffI] - st) & x[u[coeffI]];
 386             }
 387
 388             // If the matrix is symmetric, the lower triangular product
 389             // is also needed
 390             if (symmetric())
 391             {
 392                 for (register label coeffI = 0; coeffI < u.size(); coeffI++)
 393                 {
 394                     // Use transpose upper coefficient
 395                     contractLinear(lt, activeUpper[coeffI]);
 396                     expandLinear(st, lt);
 397                     sMul[u[coeffI]] -=
 398                         (activeUpper[coeffI].T() - st) & x[l[coeffI]];
 399                 }
 400             }
 401         }
 402     }
 403 }
 404
 405
 406 // ************************************************************************* //