src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixOperations.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 \*---------------------------------------------------------------------------*/
  25
  26 #include "BlockLduMatrix.H"
  27 #include "tmp.H"
  28
  29 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  30
  31 template<class Type>
  32 void Foam::BlockLduMatrix<Type>::sumDiag()
  33 {
  34     typedef typename TypeCoeffField::scalarTypeField scalarTypeField;
  35     typedef typename TypeCoeffField::linearTypeField linearTypeField;
  36     typedef typename TypeCoeffField::squareTypeField squareTypeField;
  37
  38     TypeCoeffField& Diag = this->diag();
  39
  40     const unallocLabelList& l = lduAddr().lowerAddr();
  41     const unallocLabelList& u = lduAddr().upperAddr();
  42
  43     if (this->symmetric())
  44     {
  45         // Symmetric matrix: re-use upper transpose for lower coefficients
  46
  47         const TypeCoeffField& Upper =
  48             const_cast<const BlockLduMatrix<Type>&>(*this).upper();
  49
  50         if
  51         (
  52             Upper.activeType() == blockCoeffBase::SQUARE
  53          || Diag.activeType() == blockCoeffBase::SQUARE
  54         )
  55         {
  56             const squareTypeField& activeUpper = Upper.asSquare();
  57             squareTypeField& activeDiag = Diag.asSquare();
  58
  59             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
  60             {
  61                 activeDiag[l[coeffI]] += activeUpper[coeffI].T();
  62                 activeDiag[u[coeffI]] += activeUpper[coeffI];
  63             }
  64         }
  65         else if
  66         (
  67             Upper.activeType() == blockCoeffBase::LINEAR
  68          || Diag.activeType() == blockCoeffBase::LINEAR
  69         )
  70         {
  71             const linearTypeField& activeUpper = Upper.asLinear();
  72             linearTypeField& activeDiag = Diag.asLinear();
  73
  74             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
  75             {
  76                 activeDiag[l[coeffI]] += activeUpper[coeffI];
  77                 activeDiag[u[coeffI]] += activeUpper[coeffI];
  78             }
  79         }
  80         else if
  81         (
  82             Upper.activeType() == blockCoeffBase::SCALAR
  83          || Diag.activeType() == blockCoeffBase::SCALAR
  84         )
  85         {
  86             const scalarTypeField& activeUpper = Upper.asScalar();
  87             scalarTypeField& activeDiag = Diag.asScalar();
  88
  89             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
  90             {
  91                 activeDiag[l[coeffI]] += activeUpper[coeffI];
  92                 activeDiag[u[coeffI]] += activeUpper[coeffI];
  93             }
  94         }
  95     }
  96     else if (this->asymmetric())
  97     {
  98         // Full asymmetric matrix
  99
 100         const TypeCoeffField& Lower =
 101             const_cast<const BlockLduMatrix<Type>&>(*this).lower();
 102
 103         const TypeCoeffField& Upper =
 104             const_cast<const BlockLduMatrix<Type>&>(*this).upper();
 105
 106         if
 107         (
 108             Lower.activeType() == blockCoeffBase::SQUARE
 109          || Upper.activeType() == blockCoeffBase::SQUARE
 110          || Diag.activeType() == blockCoeffBase::SQUARE
 111         )
 112         {
 113             const squareTypeField& activeLower = Lower.asSquare();
 114             const squareTypeField& activeUpper = Upper.asSquare();
 115             squareTypeField& activeDiag = Diag.asSquare();
 116
 117             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
 118             {
 119                 activeDiag[l[coeffI]] += activeLower[coeffI];
 120                 activeDiag[u[coeffI]] += activeUpper[coeffI];
 121             }
 122         }
 123         else if
 124         (
 125             Lower.activeType() == blockCoeffBase::LINEAR
 126          || Upper.activeType() == blockCoeffBase::LINEAR
 127          || Diag.activeType() == blockCoeffBase::LINEAR
 128         )
 129         {
 130             const linearTypeField& activeLower = Lower.asLinear();
 131             const linearTypeField& activeUpper = Upper.asLinear();
 132             linearTypeField& activeDiag = Diag.asLinear();
 133
 134             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
 135             {
 136                 activeDiag[l[coeffI]] += activeLower[coeffI];
 137                 activeDiag[u[coeffI]] += activeUpper[coeffI];
 138             }
 139         }
 140         else if
 141         (
 142             Lower.activeType() == blockCoeffBase::SCALAR
 143          || Upper.activeType() == blockCoeffBase::SCALAR
 144          || Diag.activeType() == blockCoeffBase::SCALAR
 145         )
 146         {
 147             const scalarTypeField& activeLower = Lower.asScalar();
 148             const scalarTypeField& activeUpper = Upper.asScalar();
 149             scalarTypeField& activeDiag = Diag.asScalar();
 150
 151             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
 152             {
 153                 activeDiag[l[coeffI]] += activeLower[coeffI];
 154                 activeDiag[u[coeffI]] += activeUpper[coeffI];
 155             }
 156         }
 157     }
 158     else
 159     {
 160         FatalErrorIn("void BlockLduMatrix<Type>::sumDiag()")
 161             << "No off-diagonal available"
 162             << abort(FatalError);
 163     }
 164 }
 165
 166
 167 template<class Type>
 168 void Foam::BlockLduMatrix<Type>::negSumDiag()
 169 {
 170     typedef typename TypeCoeffField::scalarTypeField scalarTypeField;
 171     typedef typename TypeCoeffField::linearTypeField linearTypeField;
 172     typedef typename TypeCoeffField::squareTypeField squareTypeField;
 173
 174     TypeCoeffField& Diag = this->diag();
 175
 176     const unallocLabelList& l = lduAddr().lowerAddr();
 177     const unallocLabelList& u = lduAddr().upperAddr();
 178
 179     if (this->symmetric())
 180     {
 181         // Symmetric matrix: re-use upper transpose for lower coefficients
 182
 183         const TypeCoeffField& Upper =
 184             const_cast<const BlockLduMatrix<Type>&>(*this).upper();
 185
 186         if
 187         (
 188             Upper.activeType() == blockCoeffBase::SQUARE
 189          || Diag.activeType() == blockCoeffBase::SQUARE
 190         )
 191         {
 192             const squareTypeField& activeUpper = Upper.asSquare();
 193             squareTypeField& activeDiag = Diag.asSquare();
 194
 195             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
 196             {
 197                 activeDiag[l[coeffI]] -= activeUpper[coeffI].T();
 198                 activeDiag[u[coeffI]] -= activeUpper[coeffI];
 199             }
 200         }
 201         else if
 202         (
 203             Upper.activeType() == blockCoeffBase::LINEAR
 204          || Diag.activeType() == blockCoeffBase::LINEAR
 205         )
 206         {
 207             const linearTypeField& activeUpper = Upper.asLinear();
 208             linearTypeField& activeDiag = Diag.asLinear();
 209
 210             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
 211             {
 212                 activeDiag[l[coeffI]] -= activeUpper[coeffI];
 213                 activeDiag[u[coeffI]] -= activeUpper[coeffI];
 214             }
 215         }
 216         else if
 217         (
 218             Upper.activeType() == blockCoeffBase::SCALAR
 219          || Diag.activeType() == blockCoeffBase::SCALAR
 220         )
 221         {
 222             const scalarTypeField& activeUpper = Upper.asScalar();
 223             scalarTypeField& activeDiag = Diag.asScalar();
 224
 225             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
 226             {
 227                 activeDiag[l[coeffI]] -= activeUpper[coeffI];
 228                 activeDiag[u[coeffI]] -= activeUpper[coeffI];
 229             }
 230         }
 231     }
 232     else if (this->asymmetric())
 233     {
 234         // Full asymmetric matrix
 235
 236         const TypeCoeffField& Lower =
 237             const_cast<const BlockLduMatrix<Type>&>(*this).lower();
 238
 239         const TypeCoeffField& Upper =
 240             const_cast<const BlockLduMatrix<Type>&>(*this).upper();
 241
 242         if
 243         (
 244             Lower.activeType() == blockCoeffBase::SQUARE
 245          || Upper.activeType() == blockCoeffBase::SQUARE
 246          || Diag.activeType() == blockCoeffBase::SQUARE
 247         )
 248         {
 249             const squareTypeField& activeLower = Lower.asSquare();
 250             const squareTypeField& activeUpper = Upper.asSquare();
 251             squareTypeField& activeDiag = Diag.asSquare();
 252
 253             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
 254             {
 255                 activeDiag[l[coeffI]] -= activeLower[coeffI];
 256                 activeDiag[u[coeffI]] -= activeUpper[coeffI];
 257             }
 258         }
 259         else if
 260         (
 261             Lower.activeType() == blockCoeffBase::LINEAR
 262          || Upper.activeType() == blockCoeffBase::LINEAR
 263          || Diag.activeType() == blockCoeffBase::LINEAR
 264         )
 265         {
 266             const linearTypeField& activeLower = Lower.asLinear();
 267             const linearTypeField& activeUpper = Upper.asLinear();
 268             linearTypeField& activeDiag = Diag.asLinear();
 269
 270             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
 271             {
 272                 activeDiag[l[coeffI]] -= activeLower[coeffI];
 273                 activeDiag[u[coeffI]] -= activeUpper[coeffI];
 274             }
 275         }
 276         else if
 277         (
 278             Lower.activeType() == blockCoeffBase::SCALAR
 279          || Upper.activeType() == blockCoeffBase::SCALAR
 280          || Diag.activeType() == blockCoeffBase::SCALAR
 281         )
 282         {
 283             const scalarTypeField& activeLower = Lower.asScalar();
 284             const scalarTypeField& activeUpper = Upper.asScalar();
 285             scalarTypeField& activeDiag = Diag.asScalar();
 286
 287             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
 288             {
 289                 activeDiag[l[coeffI]] -= activeLower[coeffI];
 290                 activeDiag[u[coeffI]] -= activeUpper[coeffI];
 291             }
 292         }
 293     }
 294     else
 295     {
 296         FatalErrorIn("void BlockLduMatrix<Type>::negSumDiag()")
 297             << "No off-diagonal available"
 298             << abort(FatalError);
 299     }
 300 }
 301
 302
 303 template<class Type>
 304 void Foam::BlockLduMatrix<Type>::check() const
 305 {
 306     typedef typename TypeCoeffField::scalarTypeField scalarTypeField;
 307     typedef typename TypeCoeffField::linearTypeField linearTypeField;
 308     typedef typename TypeCoeffField::squareTypeField squareTypeField;
 309
 310     // Copy the diagonal
 311     TypeCoeffField DiagCopy(this->diag().size());
 312 //     TypeCoeffField DiagCopy = this->diag();
 313
 314     const unallocLabelList& l = lduAddr().lowerAddr();
 315     const unallocLabelList& u = lduAddr().upperAddr();
 316
 317     if (this->symmetric())
 318     {
 319         // Symmetric matrix: re-use upper transpose for lower coefficients
 320
 321         const TypeCoeffField& Upper = this->upper();
 322
 323         if
 324         (
 325             Upper.activeType() == blockCoeffBase::SQUARE
 326          || DiagCopy.activeType() == blockCoeffBase::SQUARE
 327         )
 328         {
 329             const squareTypeField& activeUpper = Upper.asSquare();
 330             squareTypeField& activeDiagCopy = DiagCopy.asSquare();
 331
 332             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
 333             {
 334                 activeDiagCopy[l[coeffI]] += activeUpper[coeffI].T();
 335                 activeDiagCopy[u[coeffI]] += activeUpper[coeffI];
 336             }
 337
 338             Info<< "void BlockLduMatrix<Type>::check() const : "
 339                 << "Symmetric matrix: raw matrix difference: "
 340                 << sum(mag(activeDiagCopy))
 341                 << " scaled: "
 342                 << sum(mag(activeDiagCopy))/sum(mag(this->diag().asSquare()))
 343                 << endl;
 344         }
 345         else if
 346         (
 347             Upper.activeType() == blockCoeffBase::LINEAR
 348          || DiagCopy.activeType() == blockCoeffBase::LINEAR
 349         )
 350         {
 351             const linearTypeField& activeUpper = Upper.asLinear();
 352             linearTypeField& activeDiagCopy = DiagCopy.asLinear();
 353
 354             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
 355             {
 356                 activeDiagCopy[l[coeffI]] += activeUpper[coeffI];
 357                 activeDiagCopy[u[coeffI]] += activeUpper[coeffI];
 358             }
 359
 360             Info<< "void BlockLduMatrix<Type>::check() const : "
 361                 << "Symmetric matrix: raw matrix difference: "
 362                 << sum(mag(activeDiagCopy))
 363                 << " scaled: "
 364                 << sum(mag(activeDiagCopy))/sum(mag(this->diag().asLinear()))
 365                 << endl;
 366         }
 367         else if
 368         (
 369             Upper.activeType() == blockCoeffBase::SCALAR
 370          || DiagCopy.activeType() == blockCoeffBase::SCALAR
 371         )
 372         {
 373             const scalarTypeField& activeUpper = Upper.asScalar();
 374             scalarTypeField& activeDiagCopy = DiagCopy.asScalar();
 375
 376             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
 377             {
 378                 activeDiagCopy[l[coeffI]] += activeUpper[coeffI];
 379                 activeDiagCopy[u[coeffI]] += activeUpper[coeffI];
 380             }
 381
 382             Info<< "void BlockLduMatrix<Type>::check() const : "
 383                 << "Symmetric matrix: raw matrix difference: "
 384                 << sum(mag(activeDiagCopy))
 385                 << " scaled: "
 386                 << sum(mag(activeDiagCopy))/sum(mag(this->diag().asScalar()))
 387                 << endl;
 388         }
 389     }
 390     else if (this->asymmetric())
 391     {
 392         // Full asymmetric matrix
 393
 394         const TypeCoeffField& Lower = this->lower();
 395         const TypeCoeffField& Upper = this->upper();
 396
 397         if
 398         (
 399             Lower.activeType() == blockCoeffBase::SQUARE
 400          || Upper.activeType() == blockCoeffBase::SQUARE
 401          || DiagCopy.activeType() == blockCoeffBase::SQUARE
 402         )
 403         {
 404             const squareTypeField& activeLower = Lower.asSquare();
 405             const squareTypeField& activeUpper = Upper.asSquare();
 406             squareTypeField& activeDiagCopy = DiagCopy.asSquare();
 407
 408             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
 409             {
 410                 activeDiagCopy[l[coeffI]] += activeLower[coeffI];
 411                 activeDiagCopy[u[coeffI]] += activeUpper[coeffI];
 412             }
 413
 414             Info<< "void BlockLduMatrix<Type>::check() const : "
 415                 << "Asymmetric matrix: raw matrix difference: "
 416                 << sum(mag(activeDiagCopy))
 417                 << " scaled: "
 418                 << sum(mag(activeDiagCopy))/sum(mag(this->diag().asSquare()))
 419                 << endl;
 420         }
 421         else if
 422         (
 423             Lower.activeType() == blockCoeffBase::LINEAR
 424          || Upper.activeType() == blockCoeffBase::LINEAR
 425          || DiagCopy.activeType() == blockCoeffBase::LINEAR
 426         )
 427         {
 428             const linearTypeField& activeLower = Lower.asLinear();
 429             const linearTypeField& activeUpper = Upper.asLinear();
 430             linearTypeField& activeDiagCopy = DiagCopy.asLinear();
 431
 432             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
 433             {
 434                 activeDiagCopy[l[coeffI]] += activeLower[coeffI];
 435                 activeDiagCopy[u[coeffI]] += activeUpper[coeffI];
 436             }
 437
 438             Info<< "void BlockLduMatrix<Type>::check() const : "
 439                 << "Asymmetric matrix: raw matrix difference: "
 440                 << sum(mag(activeDiagCopy))
 441                 << " scaled: "
 442                 << sum(mag(activeDiagCopy))/sum(mag(this->diag().asLinear()))
 443                 << endl;
 444         }
 445         else if
 446         (
 447             Lower.activeType() == blockCoeffBase::SCALAR
 448          || Upper.activeType() == blockCoeffBase::SCALAR
 449          || DiagCopy.activeType() == blockCoeffBase::SCALAR
 450         )
 451         {
 452             const scalarTypeField& activeLower = Lower.asScalar();
 453             const scalarTypeField& activeUpper = Upper.asScalar();
 454             scalarTypeField& activeDiagCopy = DiagCopy.asScalar();
 455
 456             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
 457             {
 458                 activeDiagCopy[l[coeffI]] += activeLower[coeffI];
 459                 activeDiagCopy[u[coeffI]] += activeUpper[coeffI];
 460             }
 461
 462             Info<< "void BlockLduMatrix<Type>::check() const : "
 463                 << "Asymmetric matrix: raw matrix difference: "
 464                 << sum(mag(activeDiagCopy))
 465                 << " scaled: "
 466                 << sum(mag(activeDiagCopy))/sum(mag(this->diag().asScalar()))
 467                 << endl;
 468         }
 469     }
 470     else
 471     {
 472         Info<< "void BlockLduMatrix<Type>::check() const : "
 473             << "Diagonal matrix" << endl;
 474     }
 475 }
 476
 477
 478 template<class Type>
 479 void Foam::BlockLduMatrix<Type>::relax
 480 (
 481     const Field<Type>& x,
 482     Field<Type>& b,
 483     const scalar alpha
 484 )
 485 {
 486     typedef typename TypeCoeffField::scalarTypeField scalarTypeField;
 487     typedef typename TypeCoeffField::linearTypeField linearTypeField;
 488     typedef typename TypeCoeffField::squareTypeField squareTypeField;
 489
 490     //HJ Missing code: add coupling coefficients to under-relaxation
 491     //   HJ, 21/Feb/2008
 492
 493     if (alpha <= 0)
 494     {
 495         return;
 496     }
 497
 498     TypeCoeffField& Diag = this->diag();
 499
 500     // Create multiplication function object
 501     typename BlockCoeff<Type>::multiply mult;
 502
 503     const unallocLabelList& l = lduAddr().lowerAddr();
 504     const unallocLabelList& u = lduAddr().upperAddr();
 505
 506     if (this->symmetric())
 507     {
 508         // Symmetric matrix: re-use upper transpose for lower coefficients
 509
 510         const TypeCoeffField& Upper =
 511             const_cast<const BlockLduMatrix<Type>&>(*this).upper();
 512
 513         if
 514         (
 515             Upper.activeType() == blockCoeffBase::SQUARE
 516          || Diag.activeType() == blockCoeffBase::SQUARE
 517         )
 518         {
 519             const squareTypeField& activeUpper = Upper.asSquare();
 520             squareTypeField& activeDiag = Diag.asSquare();
 521
 522             // Make a copy of diagonal before relaxation
 523             squareTypeField activeDiagOld = activeDiag;
 524
 525             // There are options for under-relaxing the block diagonal
 526             // coefficient.  Currently, the complete diagonal block is
 527             // under-relaxed.  There's no checking on the off-diag sum
 528
 529             squareTypeField sumOff
 530             (
 531                 activeDiag.size(),
 532                 pTraits<typename TypeCoeffField::squareType>::zero
 533             );
 534
 535             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
 536             {
 537                 sumOff[u[coeffI]] += cmptMag(activeUpper[coeffI].T());
 538                 sumOff[l[coeffI]] += cmptMag(activeUpper[coeffI]);
 539             }
 540
 541             // Reconsider under-relaxation of square blocks.
 542             // HJ, 23/Sep/2011 (2 places)
 543             activeDiag = max(activeDiag, sumOff);
 544             activeDiag *= 1.0/alpha;
 545
 546             // Add the relaxation contribution to b
 547             forAll (b, i)
 548             {
 549                 b[i] += mult(activeDiag[i] - activeDiagOld[i], x[i]);
 550             }
 551         }
 552         else if
 553         (
 554             Upper.activeType() == blockCoeffBase::LINEAR
 555          || Diag.activeType() == blockCoeffBase::LINEAR
 556         )
 557         {
 558             const linearTypeField& activeUpper = Upper.asLinear();
 559             linearTypeField& activeDiag = Diag.asLinear();
 560
 561             // Make a copy of diagonal before relaxation
 562             linearTypeField activeDiagOld = activeDiag;
 563
 564             linearTypeField sumOff
 565             (
 566                 activeDiag.size(),
 567                 pTraits<typename TypeCoeffField::linearType>::zero
 568             );
 569
 570             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
 571             {
 572                 sumOff[u[coeffI]] += cmptMag(activeUpper[coeffI]);
 573                 sumOff[l[coeffI]] += cmptMag(activeUpper[coeffI]);
 574             }
 575
 576             activeDiag = max(activeDiag, sumOff);
 577             activeDiag *= 1.0/alpha;
 578
 579             // Add the relaxation contribution to b
 580             forAll (b, i)
 581             {
 582                 b[i] += mult(activeDiag[i] - activeDiagOld[i], x[i]);
 583             }
 584         }
 585         else if
 586         (
 587             Upper.activeType() == blockCoeffBase::SCALAR
 588          || Diag.activeType() == blockCoeffBase::SCALAR
 589         )
 590         {
 591             const scalarTypeField& activeUpper = Upper.asScalar();
 592             scalarTypeField& activeDiag = Diag.asScalar();
 593
 594             // Make a copy of diagonal before relaxation
 595             scalarTypeField activeDiagOld = activeDiag;
 596
 597             scalarTypeField sumOff
 598             (
 599                 activeDiag.size(),
 600                 pTraits<typename TypeCoeffField::scalarType>::zero
 601             );
 602
 603             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
 604             {
 605                 sumOff[u[coeffI]] += mag(activeUpper[coeffI]);
 606                 sumOff[l[coeffI]] += mag(activeUpper[coeffI]);
 607             }
 608
 609             activeDiag = max(activeDiag, sumOff);
 610             activeDiag *= 1.0/alpha;
 611
 612             // Add the relaxation contribution to b
 613             forAll (b, i)
 614             {
 615                 b[i] += (activeDiag[i] - activeDiagOld[i])*x[i];
 616             }
 617         }
 618     }
 619     else if (this->asymmetric())
 620     {
 621         // Full asymmetric matrix
 622
 623         const TypeCoeffField& Lower =
 624             const_cast<const BlockLduMatrix<Type>&>(*this).lower();
 625
 626         const TypeCoeffField& Upper =
 627             const_cast<const BlockLduMatrix<Type>&>(*this).upper();
 628
 629         if
 630         (
 631             Lower.activeType() == blockCoeffBase::SQUARE
 632          || Upper.activeType() == blockCoeffBase::SQUARE
 633          || Diag.activeType() == blockCoeffBase::SQUARE
 634         )
 635         {
 636             const squareTypeField& activeLower = Lower.asSquare();
 637             const squareTypeField& activeUpper = Upper.asSquare();
 638             squareTypeField& activeDiag = Diag.asSquare();
 639
 640             // Make a copy of diagonal before relaxation
 641             squareTypeField activeDiagOld = activeDiag;
 642
 643             // There are options for under-relaxing the block diagonal
 644             // coefficient.  Currently, the complete diagonal block is
 645             // under-relaxed.  There's no checking on the off-diag sum
 646
 647             squareTypeField sumOff
 648             (
 649                 activeDiag.size(),
 650                 pTraits<typename TypeCoeffField::squareType>::zero
 651             );
 652
 653             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
 654             {
 655                 sumOff[u[coeffI]] += cmptMag(activeLower[coeffI]);
 656                 sumOff[l[coeffI]] += cmptMag(activeUpper[coeffI]);
 657             }
 658
 659             // Reconsider under-relaxation of square blocks.
 660             // HJ, 23/Sep/2011 (2 places)
 661             activeDiag = max(activeDiag, sumOff);
 662             activeDiag *= 1.0/alpha;
 663
 664             // Add the relaxation contribution to b
 665             forAll (b, i)
 666             {
 667                 b[i] += mult(activeDiag[i] - activeDiagOld[i], x[i]);
 668             }
 669         }
 670         else if
 671         (
 672             Lower.activeType() == blockCoeffBase::LINEAR
 673          || Upper.activeType() == blockCoeffBase::LINEAR
 674          || Diag.activeType() == blockCoeffBase::LINEAR
 675         )
 676         {
 677             const linearTypeField& activeLower = Lower.asLinear();
 678             const linearTypeField& activeUpper = Upper.asLinear();
 679             linearTypeField& activeDiag = Diag.asLinear();
 680
 681             // Make a copy of diagonal before relaxation
 682             linearTypeField activeDiagOld = activeDiag;
 683
 684             linearTypeField sumOff
 685             (
 686                 activeDiag.size(),
 687                 pTraits<typename TypeCoeffField::linearType>::zero
 688             );
 689
 690             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
 691             {
 692                 sumOff[u[coeffI]] += cmptMag(activeLower[coeffI]);
 693                 sumOff[l[coeffI]] += cmptMag(activeUpper[coeffI]);
 694             }
 695
 696             activeDiag = max(activeDiag, sumOff);
 697             activeDiag *= 1.0/alpha;
 698
 699             // Add the relaxation contribution to b
 700             forAll (b, i)
 701             {
 702                 b[i] += mult(activeDiag[i] - activeDiagOld[i], x[i]);
 703             }
 704         }
 705         else if
 706         (
 707             Lower.activeType() == blockCoeffBase::SCALAR
 708          || Upper.activeType() == blockCoeffBase::SCALAR
 709          || Diag.activeType() == blockCoeffBase::SCALAR
 710         )
 711         {
 712             const scalarTypeField& activeLower = Lower.asScalar();
 713             const scalarTypeField& activeUpper = Upper.asScalar();
 714             scalarTypeField& activeDiag = Diag.asScalar();
 715
 716             // Make a copy of diagonal before relaxation
 717             scalarTypeField activeDiagOld = activeDiag;
 718
 719             scalarTypeField sumOff
 720             (
 721                 activeDiag.size(),
 722                 pTraits<typename TypeCoeffField::scalarType>::zero
 723             );
 724
 725             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
 726             {
 727                 sumOff[u[coeffI]] += mag(activeLower[coeffI]);
 728                 sumOff[l[coeffI]] += mag(activeUpper[coeffI]);
 729             }
 730
 731             activeDiag = max(activeDiag, sumOff);
 732             activeDiag *= 1.0/alpha;
 733
 734             // Add the relaxation contribution to b
 735             forAll (b, i)
 736             {
 737                 b[i] += activeDiag[i] - activeDiagOld[i]*x[i];
 738             }
 739         }
 740     }
 741 }
 742
 743
 744 template<class Type>
 745 void Foam::BlockLduMatrix<Type>::setValue
 746 (
 747     const label eqnIndex,
 748     const Type& value
 749 )
 750 {
 751     BlockConstraint<Type> cp(eqnIndex, value);
 752
 753     if (!fixedEqns_.found(eqnIndex))
 754     {
 755         fixedEqns_.insert(eqnIndex, cp);
 756     }
 757     else
 758     {
 759         WarningIn
 760         (
 761             "void BlockLduMatrix<Type>::setValue(const label eqnIndex, "
 762             "const Type& value)"
 763         )   << "Adding constraint on an already constrained point."
 764             << "  Point: " << eqnIndex << endl;
 765
 766         fixedEqns_[eqnIndex].combine(cp);
 767     }
 768 }
 769
 770
 771 template<class Type>
 772 Foam::tmp<Foam::Field<Type> > Foam::BlockLduMatrix<Type>::residual
 773 (
 774     const Field<Type>& x
 775 ) const
 776 {
 777     Field<Type> Ax(x.size());
 778     Amul(Ax, x);
 779     return -Ax;
 780 }
 781
 782
 783 template<class Type>
 784 typename Foam::tmp<Foam::Field<Type> > Foam::BlockLduMatrix<Type>::residual
 785 (
 786     const Field<Type>& x,
 787     const Field<Type>& b
 788 ) const
 789 {
 790     return b + residual(x);
 791 }
 792
 793
 794 template<class Type>
 795 void Foam::BlockLduMatrix<Type>::negate()
 796 {
 797     if (lowerPtr_)
 798     {
 799         lowerPtr_->negate();
 800     }
 801
 802     if (upperPtr_)
 803     {
 804         upperPtr_->negate();
 805     }
 806
 807     if (diagPtr_)
 808     {
 809         diagPtr_->negate();
 810     }
 811
 812     // Do coupling coefficients
 813     coupleUpper_.negate();
 814     coupleLower_.negate();
 815 }
 816
 817
 818 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 819
 820
 821 template<class Type>
 822 void Foam::BlockLduMatrix<Type>::operator=(const BlockLduMatrix<Type>& A)
 823 {
 824     if (this == &A)
 825     {
 826         FatalErrorIn
 827         (
 828             "void BlockLduMatrix<Type>::operator="
 829             "(const BlockLduMatrix<Type>& A)"
 830         )   << "attempted assignment to self"
 831             << abort(FatalError);
 832     }
 833
 834     if (A.lowerPtr_)
 835     {
 836         lower() = A.lower();
 837     }
 838     else if (lowerPtr_)
 839     {
 840         delete lowerPtr_;
 841         lowerPtr_ = NULL;
 842     }
 843
 844     if (A.upperPtr_)
 845     {
 846         upper() = A.upper();
 847     }
 848     else if (upperPtr_)
 849     {
 850         delete upperPtr_;
 851         upperPtr_ = NULL;
 852     }
 853
 854     if (A.diagPtr_)
 855     {
 856         diag() = A.diag();
 857     }
 858
 859     // Copy interface data
 860     interfaces_ = A.interfaces_;
 861     coupleUpper_ = A.coupleUpper_;
 862     coupleLower_ = A.coupleLower_;
 863
 864     // Copy constraints
 865     fixedEqns_ = A.fixedEqns_;
 866 }
 867
 868
 869 template<class Type>
 870 void Foam::BlockLduMatrix<Type>::operator+=(const BlockLduMatrix<Type>& A)
 871 {
 872     // Do diagonal first
 873     if (A.thereIsDiag())
 874     {
 875         diag() += A.diag();
 876     }
 877
 878     if (A.symmetric())
 879     {
 880         upper() += A.upper();
 881
 882         if (this->asymmetric())
 883         {
 884             lower() += A.upper().transpose();
 885         }
 886     }
 887
 888     if (A.asymmetric())
 889     {
 890         upper() += A.upper();
 891         lower() += A.lower();
 892     }
 893
 894     // Interface data
 895     coupleUpper_ += A.coupleUpper_;
 896     coupleLower_ += A.coupleLower_;
 897 }
 898
 899
 900 template<class Type>
 901 void Foam::BlockLduMatrix<Type>::operator-=(const BlockLduMatrix<Type>& A)
 902 {
 903     // Do diagonal first
 904     if (A.thereIsDiag())
 905     {
 906         diag() -= A.diag();
 907     }
 908
 909     if (A.symmetric())
 910     {
 911         upper() -= A.upper();
 912
 913         if (this->asymmetric())
 914         {
 915             lower() -= A.upper().transpose();
 916         }
 917     }
 918
 919     if (A.asymmetric())
 920     {
 921         upper() -= A.upper();
 922         lower() -= A.lower();
 923     }
 924
 925     // Interface data
 926     coupleUpper_ -= A.coupleUpper_;
 927     coupleLower_ -= A.coupleLower_;
 928 }
 929
 930
 931 template<class Type>
 932 void Foam::BlockLduMatrix<Type>::operator*=(const scalarField& sf)
 933 {
 934     typedef typename TypeCoeffField::scalarTypeField scalarTypeField;
 935     typedef typename TypeCoeffField::linearTypeField linearTypeField;
 936     typedef typename TypeCoeffField::squareTypeField squareTypeField;
 937
 938     //HJ Missing code: add coupling coefficients op*=
 939     //   HJ, 21/Feb/2008
 940     // IC - Complicated because we have to scale across the interfaces
 941     // We may need to include this functionality in lduInterfaceField
 942     // as additional initInterfaceScale and scaleInterface functions
 943
 944     if (diagPtr_)
 945     {
 946         *diagPtr_ *= sf;
 947     }
 948
 949     if (upperPtr_)
 950     {
 951         TypeCoeffField& Upper = *upperPtr_;
 952
 953         const unallocLabelList& l = lduAddr().lowerAddr();
 954
 955         if (Upper.activeType() == blockCoeffBase::SCALAR)
 956         {
 957             scalarTypeField& activeUpper = Upper.asScalar();
 958
 959             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
 960             {
 961                 activeUpper[coeffI] *= sf[l[coeffI]];
 962             }
 963         }
 964         else if (Upper.activeType() == blockCoeffBase::LINEAR)
 965         {
 966             linearTypeField& activeUpper = Upper.asLinear();
 967
 968             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
 969             {
 970                 activeUpper[coeffI] *= sf[l[coeffI]];
 971             }
 972         }
 973         else if (Upper.activeType() == blockCoeffBase::SQUARE)
 974         {
 975             squareTypeField& activeUpper = Upper.asSquare();
 976
 977             for (register label coeffI = 0; coeffI < l.size(); coeffI++)
 978             {
 979                 activeUpper[coeffI] *= sf[l[coeffI]];
 980             }
 981         }
 982     }
 983
 984     if (lowerPtr_)
 985     {
 986         TypeCoeffField& Lower = *lowerPtr_;
 987
 988         const unallocLabelList& u = lduAddr().upperAddr();
 989
 990         if (Lower.activeType() == blockCoeffBase::SCALAR)
 991         {
 992             scalarTypeField& activeLower = Lower.asScalar();
 993
 994             for (register label coeffI = 0; coeffI < u.size(); coeffI++)
 995             {
 996                 activeLower[coeffI] *= sf[u[coeffI]];
 997             }
 998         }
 999         else if (Lower.activeType() == blockCoeffBase::LINEAR)
1000         {
1001             linearTypeField& activeLower = Lower.asLinear();
1002
1003             for (register label coeffI = 0; coeffI < u.size(); coeffI++)
1004             {
1005                 activeLower[coeffI] *= sf[u[coeffI]];
1006             }
1007         }
1008         else if (Lower.activeType() == blockCoeffBase::SQUARE)
1009         {
1010             squareTypeField& activeLower = Lower.asSquare();
1011
1012             for (register label coeffI = 0; coeffI < u.size(); coeffI++)
1013             {
1014                 activeLower[coeffI] *= sf[u[coeffI]];
1015             }
1016         }
1017     }
1018 }
1019
1020
1021 template<class Type>
1022 void Foam::BlockLduMatrix<Type>::operator*=(const scalar s)
1023 {
1024     if (diagPtr_)
1025     {
1026         *diagPtr_ *= s;
1027     }
1028
1029     if (upperPtr_)
1030     {
1031         *upperPtr_ *= s;
1032     }
1033
1034     if (lowerPtr_)
1035     {
1036         *lowerPtr_ *= s;
1037     }
1038
1039     // Interface data
1040     coupleUpper_ *= s;
1041     coupleLower_ *= s;
1042 }
1043
1044
1045 // ************************************************************************* //