src/foam/matrices/constraint/constraintTools.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25     A storage mechanism which allows setting of the fixed value and
  26     consequently recovering the equation for a single row of the matrix as
  27     well as the source. The equation is taken out of the matrix using a
  28     variant of compact matrix storage mechanism.
  29
  30 \*---------------------------------------------------------------------------*/
  31
  32 #include "constraint.H"
  33
  34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  35
  36 namespace Foam
  37 {
  38
  39 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  40
  41 template<class Type>
  42 template<template<class> class Matrix>
  43 void constraint<Type>::setMatrix
  44 (
  45     const Matrix<Type>& matrix
  46 )
  47 {
  48     if (matrixCoeffsSet_)
  49     {
  50         FatalErrorIn
  51         (
  52             "const scalarField& constraint<Type>::setMatrix"
  53         )   << "(const Matrix<Type>& matrix)"
  54             << "matrix coefficients already set"
  55             << abort(FatalError);
  56     }
  57
  58     matrixCoeffsSet_ = true;
  59
  60     if (matrix.hasDiag())
  61     {
  62         diagCoeff_ = matrix.diag()[rowID_];
  63     }
  64
  65     source_ = matrix.source()[rowID_];
  66
  67     const label startFaceOwn =
  68         matrix.psi().mesh().lduAddr().ownerStartAddr()[rowID_];
  69     const label endFaceOwn =
  70         matrix.psi().mesh().lduAddr().ownerStartAddr()[rowID_ + 1];
  71
  72     const label startFaceNbr =
  73         matrix.psi().mesh().lduAddr().losortStartAddr()[rowID_];
  74     const label endFaceNbr =
  75         matrix.psi().mesh().lduAddr().losortStartAddr()[rowID_ + 1];
  76
  77     const unallocLabelList& losort = matrix.psi().mesh().lduAddr().losortAddr();
  78
  79     if (matrix.hasUpper())
  80     {
  81         // get the upper coefficients
  82
  83         const scalarField& matrixUpper = matrix.upper();
  84
  85         // owner side
  86         upperCoeffsOwnerPtr_ = new scalarField(endFaceOwn - startFaceOwn);
  87
  88         scalarField& uOwn = *upperCoeffsOwnerPtr_;
  89
  90         label faceIndex = startFaceOwn;
  91
  92         forAll(uOwn, uOwnI)
  93         {
  94             uOwn[uOwnI] = matrixUpper[faceIndex];
  95
  96             faceIndex++;
  97         }
  98
  99         // neighbour side
 100         upperCoeffsNeighbourPtr_ = new scalarField(endFaceNbr - startFaceNbr);
 101
 102         scalarField& uNbr = *upperCoeffsNeighbourPtr_;
 103
 104         faceIndex = startFaceNbr;
 105
 106         forAll (uNbr, uNbrI)
 107         {
 108             uNbr[uNbrI] = matrixUpper[losort[faceIndex]];
 109
 110             faceIndex++;
 111         }
 112     }
 113
 114     if (matrix.hasLower())
 115     {
 116         // get the lower coefficients
 117
 118         const scalarField& matrixLower = matrix.lower();
 119
 120         // owner side
 121         lowerCoeffsOwnerPtr_ = new scalarField(endFaceOwn - startFaceOwn);
 122
 123         scalarField& lOwn = *lowerCoeffsOwnerPtr_;
 124
 125         label faceIndex = startFaceOwn;
 126
 127         forAll (lOwn, lOwnI)
 128         {
 129             lOwn[lOwnI] = matrixLower[faceIndex];
 130
 131             faceIndex++;
 132         }
 133
 134         // neighbour side
 135         lowerCoeffsNeighbourPtr_ = new scalarField(endFaceNbr - startFaceNbr);
 136
 137         scalarField& lNbr = *lowerCoeffsNeighbourPtr_;
 138
 139         faceIndex = startFaceNbr;
 140
 141         forAll (lNbr, lNbrI)
 142         {
 143             lNbr[lNbrI] = matrixLower[losort[faceIndex]];
 144
 145             faceIndex++;
 146         }
 147     }
 148 }
 149
 150
 151 template<class Type>
 152 template<template<class> class Matrix>
 153 void constraint<Type>::eliminateEquation
 154 (
 155     Matrix<Type>& matrix,
 156     const label rowID,
 157     const Type& value
 158 )
 159 {
 160     // Record equation as eliminated
 161     matrix.eliminatedEqns().insert(rowID);
 162
 163     Field<Type>& source = matrix.source();
 164
 165     const label startFaceOwn =
 166         matrix.psi().mesh().lduAddr().ownerStartAddr()[rowID];
 167     const label endFaceOwn =
 168         matrix.psi().mesh().lduAddr().ownerStartAddr()[rowID + 1];
 169
 170     const label startFaceNbr =
 171         matrix.psi().mesh().lduAddr().losortStartAddr()[rowID];
 172     const label endFaceNbr =
 173         matrix.psi().mesh().lduAddr().losortStartAddr()[rowID + 1];
 174
 175     const unallocLabelList& owner = matrix.psi().mesh().lduAddr().lowerAddr();
 176     const unallocLabelList& neighbour =
 177         matrix.psi().mesh().lduAddr().upperAddr();
 178     const unallocLabelList& losort =
 179         matrix.psi().mesh().lduAddr().losortAddr();
 180
 181     // My index =  rowID
 182     if (matrix.symmetric())
 183     {
 184         // get the coefficients
 185         scalarField* coeffs = NULL;
 186
 187         if (matrix.hasUpper())
 188         {
 189             coeffs = &matrix.upper();
 190         }
 191         else if (matrix.hasLower())
 192         {
 193             coeffs = &matrix.lower();
 194         }
 195
 196         scalarField& matrixCoeffs = *coeffs;
 197
 198         // owner side
 199         label faceIndex = startFaceOwn;
 200
 201         while (faceIndex < endFaceOwn)
 202         {
 203             // add contribution to source of the neighbour (I am the owner)
 204             source[neighbour[faceIndex]] -= matrixCoeffs[faceIndex]*value;
 205
 206             // I am the owner =  owner[faceIndex]
 207             // set off-diagonal coefficient to zero
 208             matrixCoeffs[faceIndex] = 0.0;
 209
 210             faceIndex++;
 211         }
 212
 213         // neighbour side
 214         faceIndex = startFaceNbr;
 215
 216         while (faceIndex < endFaceNbr)
 217         {
 218             // add contribution to source of owner (I am the neighbour)
 219             source[owner[losort[faceIndex]]] -=
 220                 matrixCoeffs[losort[faceIndex]]*value;
 221
 222             // I am the neighbour = neighbour[losort[faceIndex]]
 223             // set off-diagonal coefficient to zero
 224             matrixCoeffs[losort[faceIndex]] = 0.0;
 225
 226             faceIndex++;
 227         }
 228     }
 229     else if (matrix.asymmetric())
 230     {
 231         scalarField& matrixUpper = matrix.upper();
 232
 233         scalarField& matrixLower = matrix.lower();
 234
 235         // owner side
 236         label faceIndex = startFaceOwn;
 237
 238         while (faceIndex < endFaceOwn)
 239         {
 240             // add contribution to source of the neighbour (I am the owner)
 241             source[neighbour[faceIndex]] -= matrixLower[faceIndex]*value;
 242
 243             // set off-diagonal coefficient to zero
 244             matrixLower[faceIndex] = 0.0;
 245
 246             faceIndex++;
 247         }
 248
 249         // neighbour side
 250         faceIndex = startFaceNbr;
 251
 252         while (faceIndex < endFaceNbr)
 253         {
 254             // add contribution to source of owner (I am the neighbour)
 255             source[owner[losort[faceIndex]]] -=
 256                 matrixUpper[losort[faceIndex]]*value;
 257
 258             // set off-diagonal coefficient to zero
 259             matrixUpper[losort[faceIndex]] = 0.0;
 260
 261             faceIndex++;
 262         }
 263     }
 264 }
 265
 266
 267 template<class Type>
 268 template<template<class> class Matrix>
 269 void constraint<Type>::eliminateEquation(Matrix<Type>& matrix) const
 270 {
 271     eliminateEquation(matrix, rowID_, value_);
 272 }
 273
 274
 275 template<class Type>
 276 template<template<class> class Matrix>
 277 void constraint<Type>::eliminateEquation
 278 (
 279     Matrix<Type>& matrix,
 280     const direction d,
 281     scalarField& sourceCmpt
 282 ) const
 283 {
 284     // Record equation as eliminated
 285     matrix.eliminatedEqns().insert(rowID_);
 286
 287     const Type& fc = fixedComponents();
 288
 289     const scalar fcOfD = componentOfValue(fc, d);
 290
 291     if (fcOfD > SMALL)
 292     {
 293         const label startFaceOwn =
 294             matrix.psi().mesh().lduAddr().ownerStartAddr()[rowID_];
 295         const label endFaceOwn =
 296             matrix.psi().mesh().lduAddr().ownerStartAddr()[rowID_ + 1];
 297
 298         const label startFaceNbr =
 299             matrix.psi().mesh().lduAddr().losortStartAddr()[rowID_];
 300         const label endFaceNbr =
 301             matrix.psi().mesh().lduAddr().losortStartAddr()[rowID_ + 1];
 302
 303         const unallocLabelList& owner =
 304             matrix.psi().mesh().lduAddr().lowerAddr();
 305         const unallocLabelList& neighbour =
 306             matrix.psi().mesh().lduAddr().upperAddr();
 307         const unallocLabelList& losort =
 308             matrix.psi().mesh().lduAddr().losortAddr();
 309
 310         // My index =  rowID_
 311         if (matrix.symmetric())
 312         {
 313             // get the coefficients
 314             scalarField* coeffs = NULL;
 315
 316             if (matrix.hasUpper())
 317             {
 318                 coeffs = &matrix.upper();
 319             }
 320             else if (matrix.hasLower())
 321             {
 322                 coeffs = &matrix.lower();
 323             }
 324
 325             scalarField& matrixCoeffs = *coeffs;
 326
 327             // owner side
 328             label faceIndex = startFaceOwn;
 329
 330             while (faceIndex < endFaceOwn)
 331             {
 332                 // add contribution to source of the neighbour (I am the owner)
 333                 sourceCmpt[neighbour[faceIndex]] -=
 334                     fcOfD*matrixCoeffs[faceIndex]*componentOfValue(value(), d);
 335
 336                 // I am the owner =  owner[faceIndex]
 337                 // set off-diagonal coefficient to zero
 338                 matrixCoeffs[faceIndex] *= 1.0 - fcOfD;
 339
 340                 faceIndex++;
 341             }
 342
 343             // neighbour side
 344             faceIndex = startFaceNbr;
 345
 346             while (faceIndex < endFaceNbr)
 347             {
 348                 // add contribution to source of owner (I am the neighbour)
 349                 sourceCmpt[owner[losort[faceIndex]]] -=
 350                     fcOfD*matrixCoeffs[losort[faceIndex]]
 351                     *componentOfValue(value(), d);
 352
 353                 // I am the neighbour = neighbour[losort[faceIndex]]
 354                 // set off-diagonal coefficient to zero
 355                 matrixCoeffs[losort[faceIndex]] *= 1.0 - fcOfD;
 356
 357                 faceIndex++;
 358             }
 359         }
 360         else if (matrix.asymmetric())
 361         {
 362             scalarField& matrixUpper = matrix.upper();
 363
 364             scalarField& matrixLower = matrix.lower();
 365
 366             // owner side
 367             label faceIndex = startFaceOwn;
 368
 369             while (faceIndex < endFaceOwn)
 370             {
 371                 // add contribution to source of the neighbour (I am the owner)
 372                 sourceCmpt[neighbour[faceIndex]] -=
 373                     fcOfD*matrixLower[faceIndex]*componentOfValue(value(), d);
 374
 375                 // set off-diagonal coefficient to zero
 376                 matrixLower[faceIndex] *= 1.0 - fcOfD;
 377
 378                 faceIndex++;
 379             }
 380
 381             // neighbour side
 382             faceIndex = startFaceNbr;
 383
 384             while (faceIndex < endFaceNbr)
 385             {
 386                 // add contribution to source of owner (I am the neighbour)
 387                 sourceCmpt[owner[losort[faceIndex]]] -=
 388                     fcOfD*matrixUpper[losort[faceIndex]]
 389                     *componentOfValue(value(), d);
 390
 391                 // set off-diagonal coefficient to zero
 392                 matrixUpper[losort[faceIndex]] *= 1.0 - fcOfD;
 393
 394                 faceIndex++;
 395             }
 396         }
 397     }
 398 }
 399
 400
 401 template<class Type>
 402 template<template<class> class Matrix>
 403 void constraint<Type>::setSource
 404 (
 405     Matrix<Type>& matrix,
 406     const label rowID,
 407     const Type& value
 408 )
 409 {
 410     matrix.source()[rowID] = matrix.diag()[rowID]*value;
 411     matrix.psi()[rowID] = value;
 412 }
 413
 414
 415 template<class Type>
 416 template<template<class> class Matrix>
 417 void constraint<Type>::setSource(Matrix<Type>& matrix) const
 418 {
 419     setSource(matrix, rowID_, value_);
 420 }
 421
 422
 423 template<class Type>
 424 template<template<class> class Matrix>
 425 void constraint<Type>::setSourceDiag
 426 (
 427     Matrix<Type>& matrix,
 428     const direction d,
 429     scalarField& psiCmpt,
 430     scalarField& sourceCmpt
 431 ) const
 432 {
 433     const Type& fc = fixedComponents();
 434
 435     if (componentOfValue(fc, d) > SMALL)
 436     {
 437         sourceCmpt[rowID()] = componentOfValue(fc, d)*matrix.diag()[rowID()]
 438             *componentOfValue(value(), d);
 439
 440         // set the solution to the right value as well
 441         psiCmpt[rowID()] =
 442             componentOfValue(fc, d)*componentOfValue(value(), d);
 443     }
 444 }
 445
 446
 447 template<class Type>
 448 template<template<class> class Matrix>
 449 void constraint<Type>::reconstructMatrix
 450 (
 451     Matrix<Type>& matrix
 452 ) const
 453 {
 454     if (!matrixCoeffsSet_)
 455     {
 456         FatalErrorIn
 457         (
 458             "void constraint<Type>::reconstructMatrix("
 459             "Matrix<Type>& matrix)"
 460         )   << "matrix coefficients not set"
 461             << abort(FatalError);
 462     }
 463
 464     if (matrix.hasDiag())
 465     {
 466         matrix.diag()[rowID_] = diagCoeff_;
 467     }
 468
 469     const label startFaceOwn =
 470         matrix.psi().mesh().lduAddr().ownerStartAddr()[rowID_];
 471
 472     const label startFaceNbr =
 473         matrix.psi().mesh().lduAddr().losortStartAddr()[rowID_];
 474
 475     const unallocLabelList& losort = matrix.psi().mesh().lduAddr().losortAddr();
 476
 477     if (matrix.hasUpper())
 478     {
 479         // get the upper coefficients
 480
 481         scalarField& matrixUpper = matrix.upper();
 482
 483         // owner side
 484         const scalarField& uOwn = upperCoeffsOwner();
 485
 486         label faceIndex = startFaceOwn;
 487
 488         forAll(uOwn, uOwnI)
 489         {
 490             matrixUpper[faceIndex] = uOwn[uOwnI];
 491
 492             faceIndex++;
 493         }
 494
 495         // neighbour side
 496         const scalarField& uNbr = upperCoeffsNeighbour();
 497
 498         faceIndex = startFaceNbr;
 499
 500         forAll (uNbr, uNbrI)
 501         {
 502             matrixUpper[losort[faceIndex]] = uNbr[uNbrI];
 503
 504             faceIndex++;
 505         }
 506     }
 507
 508     if (matrix.hasLower())
 509     {
 510         // get the lower coefficients
 511
 512         scalarField& matrixLower = matrix.lower();
 513
 514         // owner side
 515         const scalarField& lOwn = lowerCoeffsOwner();
 516
 517         label faceIndex = startFaceOwn;
 518
 519         forAll (lOwn, lOwnI)
 520         {
 521             matrixLower[faceIndex] = lOwn[lOwnI];
 522
 523             faceIndex++;
 524         }
 525
 526         // neighbour side
 527         const scalarField& lNbr = lowerCoeffsNeighbour();
 528
 529         faceIndex = startFaceNbr;
 530
 531         forAll (lNbr, lNbrI)
 532         {
 533             matrixLower[losort[faceIndex]] = lNbr[lNbrI];
 534
 535             faceIndex++;
 536         }
 537     }
 538 }
 539
 540
 541 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 542
 543 } // End namespace Foam
 544
 545 // ************************************************************************* //