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