src/blockMatrix/BlockLduMatrix/BlockLduMatrix.H

   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 Class
  26     BlockLduMatrix
  27
  28 Description
  29     BlockLduMatrix is a general matrix class in which the coefficients are
  30     stored as three arrays, one for the upper triangle, one for the
  31     lower triangle and a third for the diagonal.  Addressing object must
  32     be supplied for the upper and lower triangles.
  33
  34 Author
  35     Hrvoje Jasak, Wikki Ltd.  All rights reserved.
  36
  37 SourceFiles
  38     BlockLduMatrix.C
  39     BlockLduMatrixOperations.C
  40     BlockLduMatrixUpdateInterfaces.C
  41     BlockLduMatrixATmul.C
  42     BlockLduMatrixHOps.C
  43     BlockLduMatrixDecouple.C
  44     BlockLduMatrixDecoupledHOps.C
  45
  46 \*---------------------------------------------------------------------------*/
  47
  48 #ifndef BlockLduMatrix_H
  49 #define BlockLduMatrix_H
  50
  51 #include "coeffFields.H"
  52 #include "lduMesh.H"
  53 #include "BlockLduInterface.H"
  54 #include "Map.H"
  55
  56 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  57
  58 namespace Foam
  59 {
  60
  61 // * * * * * * Forward declaration of template friend fuctions * * * * * * * //
  62
  63 template<class Type>
  64 class BlockLduMatrix;
  65
  66 template<class Type>
  67 Ostream& operator<<(Ostream&, const BlockLduMatrix<Type>&);
  68
  69 template<class Type>
  70 class BlockConstraint;
  71
  72
  73 /*---------------------------------------------------------------------------*\
  74                         Class BlockLduMatrix Declaration
  75 \*---------------------------------------------------------------------------*/
  76
  77 template<class Type>
  78 class BlockLduMatrix
  79 :
  80     public refCount
  81 {
  82 public:
  83
  84     // Public data types
  85
  86         typedef CoeffField<Type> TypeCoeffField;
  87         typedef Field<Type> TypeField;
  88         typedef BlockConstraint<Type> ConstraintType;
  89
  90
  91 private:
  92
  93     // Private data
  94
  95         // LDU mesh reference
  96         const lduMesh& lduMesh_;
  97
  98
  99         // Block matrix elements
 100
 101             //- Diagonal coefficients
 102             CoeffField<Type>* diagPtr_;
 103
 104             //- Upper triangle coefficients.  Also used for symmetric matrix
 105             CoeffField<Type>* upperPtr_;
 106
 107             //- Lower triangle coefficients
 108             CoeffField<Type> *lowerPtr_;
 109
 110
 111         // Coupling
 112
 113             //- List of coupled interfaces
 114             PtrList<BlockLduInterface<Type> > interfaces_;
 115
 116             //- Coupled interface coefficients, upper
 117             FieldField<CoeffField, Type> coupleUpper_;
 118
 119             //- Coupled interface coefficients, lower
 120             FieldField<CoeffField, Type> coupleLower_;
 121
 122
 123         // Constraints
 124
 125             //- Equation triangle map
 126             mutable Map<ConstraintType> fixedEqns_;
 127
 128
 129     // Private static data
 130
 131         //- Matrix constraint fill-in
 132         //  Equals to the estimated fraction of fixed nodes in the matrix
 133         static const label fixFillIn;
 134
 135
 136     // Private member functions
 137
 138         // Decoupled versions of nmatrix operations
 139
 140             //- Sum off-diagonal coefficients and add to diagonal,
 141             //  decoupled version
 142             void decoupledSumDiag();
 143
 144             //- Sum negative off-diagonal coefficients and add to diagonal,
 145             //  decoupled version
 146             void decoupledNegSumDiag();
 147
 148             //- Check matrix for diagonal dominance, decoupled version
 149             void decoupledCheck() const;
 150
 151             //- Relax matrix, decoupled version
 152             void decoupledRelax
 153             (
 154                 const TypeField& x,
 155                 TypeField& b,
 156                 const scalar alpha
 157             );
 158
 159             //- Matrix scaling with scalar field, decoupled version
 160             void decoupledMultEqOp(const scalarField& sf);
 161
 162             //- Matrix multiplication without coupled interface update,
 163             //  decoupled version
 164             void decoupledAmulCore
 165             (
 166                 TypeField& Ax,
 167                 const TypeField& x
 168             ) const;
 169
 170             //- Matrix transpose multiplication without coupled interface update
 171             //  decoupled version
 172             void decoupledTmulCore
 173             (
 174                 TypeField& Tx,
 175                 const TypeField& x
 176             ) const;
 177
 178             //- Return L-U vector-matrix multiplication in row-form,
 179             //  decoupled version
 180             tmp<TypeField> decoupledH(const TypeField& x) const;
 181
 182             //- Return L-U  vector-matrix multiplication in off-diagonal form,
 183             //  decoupled version
 184             tmp<TypeField> decoupledFaceH(const TypeField& x) const;
 185
 186
 187 protected:
 188
 189     // Access to constraints
 190
 191         //- Return constraint map
 192         const Map<ConstraintType>& fixedEqns() const
 193         {
 194             return fixedEqns_;
 195         }
 196
 197
 198         //- Return access constraint map
 199         Map<ConstraintType>& fixedEqns()
 200         {
 201             return fixedEqns_;
 202         }
 203
 204
 205 public:
 206
 207     //- Runtime type information
 208     TypeName("BlockLduMatrix");
 209
 210
 211     // Constructors
 212
 213         //- Construct given addressing
 214         explicit BlockLduMatrix(const lduMesh&);
 215
 216         //- Construct as copy
 217         BlockLduMatrix(const BlockLduMatrix<Type>&);
 218
 219
 220     // Destructor
 221
 222         virtual ~BlockLduMatrix();
 223
 224
 225     // Member functions
 226
 227         // Access to addressing
 228
 229             //- Return the LDU mesh from which the addressing is obtained
 230             const lduMesh& mesh() const
 231             {
 232                 return lduMesh_;
 233             }
 234
 235             //- Return the LDU addressing
 236             const lduAddressing& lduAddr() const
 237             {
 238                 return lduMesh_.lduAddr();
 239             }
 240
 241             //- Return the patch evaluation schedule
 242             const lduSchedule& patchSchedule() const
 243             {
 244                 return lduAddr().patchSchedule();
 245             }
 246
 247
 248         // Access to coefficients
 249
 250             //- Return access to diagonal coefficients
 251             TypeCoeffField& diag();
 252
 253             //- Return diagonal coefficients
 254             const TypeCoeffField& diag() const;
 255
 256             //- Return access to upper coefficients
 257             //  Also used for symmetric matrices
 258             TypeCoeffField& upper();
 259
 260             //- Return upper coefficients
 261             //  Also used for symmetric matrices
 262             const TypeCoeffField& upper() const;
 263
 264             //- Return access to lower coefficients
 265             TypeCoeffField& lower();
 266
 267             //- Return lower coefficients
 268             const TypeCoeffField& lower() const;
 269
 270             //- Return access to coupled interface coefficients, upper
 271             FieldField<CoeffField, Type>& coupleUpper()
 272             {
 273                 return coupleUpper_;
 274             }
 275
 276             //- Return coupled interface coefficients, upper
 277             const FieldField<CoeffField, Type>& coupleUpper() const
 278             {
 279                 return coupleUpper_;
 280             }
 281
 282             //- Return access to coupled interface coefficients, lower
 283             FieldField<CoeffField, Type>& coupleLower()
 284             {
 285                 return coupleLower_;
 286             }
 287
 288             //- Return coupled interface coefficients, lower
 289             const FieldField<CoeffField, Type>& coupleLower() const
 290             {
 291                 return coupleLower_;
 292             }
 293
 294
 295         // Matrix structure
 296
 297             //- Return true if there is a diagonal
 298             bool thereIsDiag() const
 299             {
 300                 return (diagPtr_);
 301             }
 302
 303             //- Return true if upper triangle is allocated
 304             bool thereIsUpper() const
 305             {
 306                 return (upperPtr_);
 307             }
 308
 309             //- Return true if lower triangle is allocated
 310             bool thereIsLower() const
 311             {
 312                 return (lowerPtr_);
 313             }
 314
 315             //- Return true if matrix is empty
 316             bool empty() const;
 317
 318             //- Return true if matrix is diagonal-only
 319             bool diagonal() const;
 320
 321             //- Return true if matrix is symmetric
 322             bool symmetric() const;
 323
 324             //- Return true if matrix is asymmetric
 325             bool asymmetric() const;
 326
 327             //- Return true if matrix is component-coupled
 328             bool componentCoupled() const;
 329
 330
 331         // Operations
 332
 333             //- Sum off-diagonal coefficients and add to diagonal
 334             void sumDiag();
 335
 336             //- Sum negative off-diagonal coefficients and add to diagonal
 337             void negSumDiag();
 338
 339             //- Check matrix for diagonal dominance
 340             void check() const;
 341
 342             //- Relax matrix
 343             void relax
 344             (
 345                 const TypeField& x,
 346                 TypeField& b,
 347                 const scalar alpha
 348             );
 349
 350             //- Matrix multiplication
 351             void Amul
 352             (
 353                 TypeField& Ax,
 354                 const TypeField& x
 355             ) const;
 356
 357             //- Matrix multiplication without coupled interface update
 358             void AmulCore
 359             (
 360                 TypeField& Ax,
 361                 const TypeField& x
 362             ) const;
 363
 364             //- Matrix transpose multiplication
 365             void Tmul
 366             (
 367                 TypeField& Ax,
 368                 const TypeField& x
 369             ) const;
 370
 371             //- Matrix transpose multiplication without
 372             //  coupled interface update
 373             void TmulCore
 374             (
 375                 TypeField& Ax,
 376                 const TypeField& x
 377             ) const;
 378
 379
 380             //- Return decoupled b
 381             void segregateB
 382             (
 383                 TypeField& sMul,
 384                 const TypeField& x
 385             ) const;
 386
 387
 388         // Coupled interface functionality
 389
 390             //- Initialise the update of coupled interfaces
 391             //  for Amul operations
 392             void initInterfaces
 393             (
 394                 TypeField& Ax,
 395                 const TypeField& x
 396             ) const;
 397
 398             //- Update coupled interfaces
 399             void updateInterfaces
 400             (
 401                 const FieldField<CoeffField, Type>& coeffs,
 402                 TypeField& Ax,
 403                 const TypeField& x
 404             ) const;
 405
 406
 407         // Constraint manipulation
 408
 409             //- Set constrained value in a prescribed point
 410             void setValue
 411             (
 412                 const label eqnIndex,
 413                 const Type& value
 414             );
 415
 416
 417         // Residual calculation
 418
 419             //- Calculate residual
 420             tmp<TypeField> residual
 421             (
 422                 const TypeField& x
 423             ) const;
 424
 425             tmp<TypeField> residual
 426             (
 427                 const TypeField& x,
 428                 const TypeField& b
 429             ) const;
 430
 431
 432         // H-operations
 433
 434             //- Return L-U vector-matrix multiplication in row-form
 435             tmp<TypeField> H(const TypeField&) const;
 436
 437             //- Return L-U  vector-matrix multiplication in off-diagonal form
 438             tmp<TypeField> faceH(const TypeField&) const;
 439
 440
 441     // Member operators
 442
 443         void operator=(const BlockLduMatrix<Type>&);
 444
 445         void negate();
 446
 447         void operator+=(const BlockLduMatrix<Type>&);
 448         void operator-=(const BlockLduMatrix<Type>&);
 449
 450         void operator*=(const scalarField&);
 451         void operator*=(const scalar);
 452
 453
 454     // Ostream operator
 455
 456         friend Ostream& operator<< <Type>
 457         (
 458             Ostream&,
 459             const BlockLduMatrix<Type>&
 460         );
 461 };
 462
 463
 464 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 465
 466 } // End namespace Foam
 467
 468 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 469
 470 #ifdef NoRepository
 471 #   include "BlockLduMatrix.C"
 472 #   include "BlockLduMatrixOperations.C"
 473 #   include "BlockLduMatrixUpdateInterfaces.C"
 474 #   include "BlockLduMatrixATmul.C"
 475 #   include "BlockLduMatrixHOps.C"
 476 #   include "BlockLduMatrixDecouple.C"
 477 #   include "BlockLduMatrixDecoupledHOps.C"
 478 #endif
 479
 480 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 481
 482 #endif
 483
 484 // ************************************************************************* //