src/finiteArea/faMatrices/faMatrix/faMatrix.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 Description
  26     Finite-Area matrix
  27
  28 \*---------------------------------------------------------------------------*/
  29
  30 #include "areaFields.H"
  31 #include "edgeFields.H"
  32 #include "calculatedFaPatchFields.H"
  33 #include "zeroGradientFaPatchFields.H"
  34 #include "coupledFaPatchFields.H"
  35
  36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  37
  38 namespace Foam
  39 {
  40
  41 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  42
  43 template<class Type>
  44 template<class Type2>
  45 void faMatrix<Type>::addToInternalField
  46 (
  47     const unallocLabelList& addr,
  48     const Field<Type2>& pf,
  49     Field<Type2>& intf
  50 ) const
  51 {
  52     if (addr.size() != pf.size())
  53     {
  54         FatalErrorIn
  55         (
  56             "faMatrix<Type>::addToInternalField(const unallocLabelList&, "
  57             "const Field&, Field&)"
  58         )   << "sizes of addressing and field are different"
  59             << abort(FatalError);
  60     }
  61
  62     forAll(addr, faceI)
  63     {
  64         intf[addr[faceI]] += pf[faceI];
  65     }
  66 }
  67
  68
  69 template<class Type>
  70 template<class Type2>
  71 void faMatrix<Type>::addToInternalField
  72 (
  73     const unallocLabelList& addr,
  74     const tmp<Field<Type2> >& tpf,
  75     Field<Type2>& intf
  76 ) const
  77 {
  78     addToInternalField(addr, tpf(), intf);
  79     tpf.clear();
  80 }
  81
  82
  83 template<class Type>
  84 template<class Type2>
  85 void faMatrix<Type>::subtractFromInternalField
  86 (
  87     const unallocLabelList& addr,
  88     const Field<Type2>& pf,
  89     Field<Type2>& intf
  90 ) const
  91 {
  92     if (addr.size() != pf.size())
  93     {
  94         FatalErrorIn
  95         (
  96             "faMatrix<Type>::addToInternalField(const unallocLabelList&, "
  97             "const Field&, Field&)"
  98         )   << "sizes of addressing and field are different"
  99             << abort(FatalError);
 100     }
 101
 102     forAll(addr, faceI)
 103     {
 104         intf[addr[faceI]] -= pf[faceI];
 105     }
 106 }
 107
 108
 109 template<class Type>
 110 template<class Type2>
 111 void faMatrix<Type>::subtractFromInternalField
 112 (
 113     const unallocLabelList& addr,
 114     const tmp<Field<Type2> >& tpf,
 115     Field<Type2>& intf
 116 ) const
 117 {
 118     subtractFromInternalField(addr, tpf(), intf);
 119     tpf.clear();
 120 }
 121
 122
 123 template<class Type>
 124 void faMatrix<Type>::addBoundaryDiag
 125 (
 126     scalarField& diag,
 127     const direction solveCmpt
 128 ) const
 129 {
 130     forAll(internalCoeffs_, patchI)
 131     {
 132         addToInternalField
 133         (
 134             lduAddr().patchAddr(patchI),
 135             internalCoeffs_[patchI].component(solveCmpt),
 136             diag
 137         );
 138     }
 139 }
 140
 141
 142 template<class Type>
 143 void faMatrix<Type>::addCmptAvBoundaryDiag(scalarField& diag) const
 144 {
 145     forAll(internalCoeffs_, patchI)
 146     {
 147         addToInternalField
 148         (
 149             lduAddr().patchAddr(patchI),
 150             cmptAv(internalCoeffs_[patchI]),
 151             diag
 152         );
 153     }
 154 }
 155
 156
 157 template<class Type>
 158 void faMatrix<Type>::addBoundarySource
 159 (
 160     Field<Type>& source,
 161     const bool couples
 162 ) const
 163 {
 164     forAll(psi_.boundaryField(), patchI)
 165     {
 166         const faPatchField<Type>& ptf = psi_.boundaryField()[patchI];
 167         const Field<Type>& pbc = boundaryCoeffs_[patchI];
 168
 169         if (!ptf.coupled())
 170         {
 171             addToInternalField(lduAddr().patchAddr(patchI), pbc, source);
 172         }
 173         else if (couples)
 174         {
 175             tmp<Field<Type> > tpnf = ptf.patchNeighbourField();
 176             const Field<Type>& pnf = tpnf();
 177
 178             const unallocLabelList& addr = lduAddr().patchAddr(patchI);
 179
 180             forAll(addr, facei)
 181             {
 182                 source[addr[facei]] += cmptMultiply(pbc[facei], pnf[facei]);
 183             }
 184         }
 185     }
 186 }
 187
 188
 189 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * //
 190
 191 template<class Type>
 192 faMatrix<Type>::faMatrix
 193 (
 194     GeometricField<Type, faPatchField, areaMesh>& psi,
 195     const dimensionSet& ds
 196 )
 197 :
 198     lduMatrix(psi.mesh()),
 199     psi_(psi),
 200     dimensions_(ds),
 201     source_(psi.size(), pTraits<Type>::zero),
 202     internalCoeffs_(psi.mesh().boundary().size()),
 203     boundaryCoeffs_(psi.mesh().boundary().size()),
 204     faceFluxCorrectionPtr_(NULL),
 205     solvingComponent(0)
 206 {
 207     if (debug)
 208     {
 209         Info<< "faMatrix<Type>(GeometricField<Type, faPatchField, areaMesh>&,"
 210                " const dimensionSet&) : "
 211                "constructing faMatrix<Type> for field " << psi_.name()
 212             << endl;
 213     }
 214
 215     // Initialise coupling coefficients
 216     forAll (psi.mesh().boundary(), patchI)
 217     {
 218         internalCoeffs_.set
 219         (
 220             patchI,
 221             new Field<Type>
 222             (
 223                 psi.mesh().boundary()[patchI].size(),
 224                 pTraits<Type>::zero
 225             )
 226         );
 227
 228         boundaryCoeffs_.set
 229         (
 230             patchI,
 231             new Field<Type>
 232             (
 233                 psi.mesh().boundary()[patchI].size(),
 234                 pTraits<Type>::zero
 235             )
 236         );
 237     }
 238
 239     psi_.boundaryField().updateCoeffs();
 240 }
 241
 242
 243 template<class Type>
 244 faMatrix<Type>::faMatrix(const faMatrix<Type>& fam)
 245 :
 246     refCount(),
 247     lduMatrix(fam),
 248     psi_(fam.psi_),
 249     dimensions_(fam.dimensions_),
 250     source_(fam.source_),
 251     internalCoeffs_(fam.internalCoeffs_),
 252     boundaryCoeffs_(fam.boundaryCoeffs_),
 253     faceFluxCorrectionPtr_(NULL),
 254     solvingComponent(fam.solvingComponent)
 255 {
 256     if (debug)
 257     {
 258         Info<< "faMatrix<Type>::faMatrix(const faMatrix<Type>&) : "
 259             << "copying faMatrix<Type> for field " << psi_.name()
 260             << endl;
 261     }
 262
 263     if (fam.faceFluxCorrectionPtr_)
 264     {
 265         faceFluxCorrectionPtr_ = new
 266         GeometricField<Type, faePatchField, edgeMesh>
 267         (
 268             *(fam.faceFluxCorrectionPtr_)
 269         );
 270     }
 271 }
 272
 273
 274 template<class Type>
 275 faMatrix<Type>::faMatrix
 276 (
 277     GeometricField<Type, faPatchField, areaMesh>& psi,
 278     Istream& is
 279 )
 280 :
 281     lduMatrix(psi.mesh()),
 282     psi_(psi),
 283     dimensions_(is),
 284     source_(is),
 285     internalCoeffs_(psi.mesh().boundary().size()),
 286     boundaryCoeffs_(psi.mesh().boundary().size()),
 287     faceFluxCorrectionPtr_(NULL),
 288     solvingComponent(0)
 289 {
 290     if (debug)
 291     {
 292         Info<< "faMatrix<Type>(GeometricField<Type, faPatchField, areaMesh>&,"
 293                " Istream&) : "
 294                "constructing faMatrix<Type> for field " << psi_.name()
 295             << endl;
 296     }
 297
 298     // Initialise coupling coefficients
 299     forAll (psi.mesh().boundary(), patchI)
 300     {
 301         internalCoeffs_.set
 302         (
 303             patchI,
 304             new Field<Type>
 305             (
 306                 psi.mesh().boundary()[patchI].size(),
 307                 pTraits<Type>::zero
 308             )
 309         );
 310
 311         boundaryCoeffs_.set
 312         (
 313             patchI,
 314             new Field<Type>
 315             (
 316                 psi.mesh().boundary()[patchI].size(),
 317                 pTraits<Type>::zero
 318             )
 319         );
 320     }
 321
 322 }
 323
 324
 325 template<class Type>
 326 faMatrix<Type>::~faMatrix()
 327 {
 328     if (debug)
 329     {
 330         Info<< "faMatrix<Type>::~faMatrix<Type>() : "
 331             << "destroying faMatrix<Type> for field " << psi_.name()
 332             << endl;
 333     }
 334
 335     if (faceFluxCorrectionPtr_)
 336     {
 337         delete faceFluxCorrectionPtr_;
 338     }
 339 }
 340
 341
 342 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 343
 344 // Set solution in given faces and eliminate corresponding
 345 // equations from the matrix
 346 template<class Type>
 347 void faMatrix<Type>::setValues
 348 (
 349     const labelList& faceLabels,
 350     const Field<Type>& values
 351 )
 352 {
 353     const faMesh& mesh = psi_.mesh();
 354
 355     const faceList& faces = mesh.faces();
 356     const unallocLabelList& own = mesh.owner();
 357     const unallocLabelList& nei = mesh.neighbour();
 358
 359     scalarField& Diag = diag();
 360
 361     forAll(faceLabels, i)
 362     {
 363         label facei = faceLabels[i];
 364
 365         psi_[facei] = values[i];
 366         source_[facei] = values[i]*Diag[facei];
 367
 368         if (symmetric() || asymmetric())
 369         {
 370             const face& c = faces[facei];
 371
 372             forAll(c, j)
 373             {
 374                 label edgei = c[j];
 375
 376                 if (mesh.isInternalEdge(edgei))
 377                 {
 378                     if (symmetric())
 379                     {
 380                         if (facei == own[edgei])
 381                         {
 382                             source_[nei[edgei]] -= upper()[edgei]*values[i];
 383                         }
 384                         else
 385                         {
 386                             source_[own[edgei]] -= upper()[edgei]*values[i];
 387                         }
 388
 389                         upper()[edgei] = 0.0;
 390                     }
 391                     else
 392                     {
 393                         if (facei == own[edgei])
 394                         {
 395                             source_[nei[edgei]] -= lower()[edgei]*values[i];
 396                         }
 397                         else
 398                         {
 399                             source_[own[edgei]] -= upper()[edgei]*values[i];
 400                         }
 401
 402                         upper()[edgei] = 0.0;
 403                         lower()[edgei] = 0.0;
 404                     }
 405                 }
 406                 else
 407                 {
 408                     label patchi = mesh.boundary().whichPatch(edgei);
 409
 410                     if (internalCoeffs_[patchi].size())
 411                     {
 412                         label patchEdgei =
 413                             mesh.boundary()[patchi].whichEdge(edgei);
 414
 415                         internalCoeffs_[patchi][patchEdgei] =
 416                             pTraits<Type>::zero;
 417
 418                         boundaryCoeffs_[patchi][patchEdgei] =
 419                             pTraits<Type>::zero;
 420                     }
 421                 }
 422             }
 423         }
 424     }
 425 }
 426
 427
 428 // Set reference level for solution
 429 template<class Type>
 430 void faMatrix<Type>::setReference
 431 (
 432     const label facei,
 433     const Type& value
 434 )
 435 {
 436     if (psi_.needReference())
 437     {
 438         if (Pstream::master())
 439         {
 440             source()[facei] += diag()[facei]*value;
 441             diag()[facei] += diag()[facei];
 442         }
 443     }
 444 }
 445
 446
 447 template<class Type>
 448 void faMatrix<Type>::relax(const scalar alpha)
 449 {
 450     if (alpha <= 0)
 451     {
 452         return;
 453     }
 454
 455     Field<Type>& S = source();
 456     scalarField& D = diag();
 457
 458     // Store the current unrelaxed diagonal for use in updating the source
 459     scalarField D0(D);
 460
 461     // Calculate the sum-mag off-diagonal from the interior faces
 462     scalarField sumOff(D.size(), 0.0);
 463     sumMagOffDiag(sumOff);
 464
 465     // Handle the boundary contributions to the diagonal
 466     forAll(psi_.boundaryField(), patchI)
 467     {
 468         const faPatchField<Type>& ptf = psi_.boundaryField()[patchI];
 469
 470         if (ptf.size())
 471         {
 472             const unallocLabelList& pa = lduAddr().patchAddr(patchI);
 473             Field<Type>& iCoeffs = internalCoeffs_[patchI];
 474
 475             if (ptf.coupled())
 476             {
 477                 const Field<Type>& pCoeffs = boundaryCoeffs_[patchI];
 478
 479                 // For coupled boundaries add the diagonal and
 480                 // off-diagonal contributions
 481                 forAll(pa, face)
 482                 {
 483                     D[pa[face]] += component(iCoeffs[face], 0);
 484                     sumOff[pa[face]] += mag(component(pCoeffs[face], 0));
 485                 }
 486             }
 487             else
 488             {
 489                 // For non-coupled boundaries subtract the diagonal
 490                 // contribution off-diagonal sum which avoids having to remove
 491                 // it from the diagonal later.
 492                 // Also add the source contribution from the relaxation
 493                 forAll(pa, face)
 494                 {
 495                     Type iCoeff0 = iCoeffs[face];
 496                     iCoeffs[face] = cmptMag(iCoeffs[face]);
 497                     sumOff[pa[face]] -= cmptMin(iCoeffs[face]);
 498                     iCoeffs[face] /= alpha;
 499                     S[pa[face]] +=
 500                         cmptMultiply(iCoeffs[face] - iCoeff0, psi_[pa[face]]);
 501                 }
 502             }
 503         }
 504     }
 505
 506     // Ensure the matrix is diagonally dominant...
 507     max(D, D, sumOff);
 508
 509     // ... then relax
 510     D /= alpha;
 511
 512     // Now remove the diagonal contribution from coupled boundaries
 513     forAll(psi_.boundaryField(), patchI)
 514     {
 515         const faPatchField<Type>& ptf = psi_.boundaryField()[patchI];
 516
 517         if (ptf.size())
 518         {
 519             const unallocLabelList& pa = lduAddr().patchAddr(patchI);
 520             Field<Type>& iCoeffs = internalCoeffs_[patchI];
 521
 522             if (ptf.coupled())
 523             {
 524                 forAll(pa, face)
 525                 {
 526                     D[pa[face]] -= component(iCoeffs[face], 0);
 527                 }
 528             }
 529         }
 530     }
 531
 532     // Finally add the relaxation contribution to the source.
 533     S += (D - D0)*psi_.internalField();
 534 }
 535
 536
 537 template<class Type>
 538 void faMatrix<Type>::relax()
 539 {
 540     scalar alpha = 0;
 541
 542     if (psi_.mesh().relax(psi_.name()))
 543     {
 544         alpha = psi_.mesh().relaxationFactor(psi_.name());
 545     }
 546
 547     if (alpha > 0)
 548     {
 549         relax(alpha);
 550     }
 551 }
 552
 553
 554 template<class Type>
 555 tmp<scalarField> faMatrix<Type>::D() const
 556 {
 557     tmp<scalarField> tdiag(new scalarField(diag()));
 558     addCmptAvBoundaryDiag(tdiag());
 559     return tdiag;
 560 }
 561
 562
 563 template<class Type>
 564 tmp<areaScalarField> faMatrix<Type>::A() const
 565 {
 566     tmp<areaScalarField> tAphi
 567     (
 568         new areaScalarField
 569         (
 570             IOobject
 571             (
 572                 "A("+psi_.name()+')',
 573                 psi_.instance(),
 574                 psi_.db()
 575             ),
 576             psi_.mesh(),
 577             dimensions_/psi_.dimensions()/dimArea,
 578             zeroGradientFaPatchScalarField::typeName
 579         )
 580     );
 581
 582     tAphi().internalField() = D()/psi_.mesh().S();
 583     tAphi().correctBoundaryConditions();
 584
 585     return tAphi;
 586 }
 587
 588
 589 template<class Type>
 590 tmp<GeometricField<Type, faPatchField, areaMesh> > faMatrix<Type>::H() const
 591 {
 592     tmp<GeometricField<Type, faPatchField, areaMesh> > tHphi
 593     (
 594         new GeometricField<Type, faPatchField, areaMesh>
 595         (
 596             IOobject
 597             (
 598                 "H("+psi_.name()+')',
 599                 psi_.instance(),
 600                 psi_.db()
 601             ),
 602             psi_.mesh(),
 603             dimensions_/dimArea,
 604             zeroGradientFaPatchScalarField::typeName
 605         )
 606     );
 607
 608     // Loop over field components
 609     for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
 610     {
 611         scalarField psiCmpt = psi_.internalField().component(cmpt);
 612
 613         scalarField boundaryDiagCmpt(psi_.size(), 0.0);
 614         addBoundaryDiag(boundaryDiagCmpt, cmpt);
 615         boundaryDiagCmpt.negate();
 616         addCmptAvBoundaryDiag(boundaryDiagCmpt);
 617
 618         tHphi().internalField().replace(cmpt, boundaryDiagCmpt*psiCmpt);
 619     }
 620
 621     tHphi().internalField() += lduMatrix::H(psi_.internalField()) + source_;
 622     addBoundarySource(tHphi().internalField());
 623
 624     tHphi().internalField() /= psi_.mesh().S();
 625     tHphi().correctBoundaryConditions();
 626
 627     return tHphi;
 628 }
 629
 630
 631 template<class Type>
 632 tmp<GeometricField<Type, faePatchField, edgeMesh> > faMatrix<Type>::
 633 flux() const
 634 {
 635     if (!psi_.mesh().fluxRequired(psi_.name()))
 636     {
 637         FatalErrorIn("faMatrix<Type>::flux()")
 638             << "flux requested but " << psi_.name()
 639             << " not specified in the fluxRequired sub-dictionary of faSchemes."
 640             << abort(FatalError);
 641     }
 642
 643     // construct GeometricField<Type, faePatchField, edgeMesh>
 644     tmp<GeometricField<Type, faePatchField, edgeMesh> > tfieldFlux
 645     (
 646         new GeometricField<Type, faePatchField, edgeMesh>
 647         (
 648             IOobject
 649             (
 650                 "flux("+psi_.name()+')',
 651                 psi_.instance(),
 652                 psi_.db()
 653             ),
 654             psi_.mesh(),
 655             dimensions()
 656         )
 657     );
 658     GeometricField<Type, faePatchField, edgeMesh>& fieldFlux = tfieldFlux();
 659
 660     for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
 661     {
 662         fieldFlux.internalField().replace
 663         (
 664             cmpt,
 665             lduMatrix::faceH(psi_.internalField().component(cmpt))
 666         );
 667     }
 668
 669     FieldField<Field, Type> InternalContrib = internalCoeffs_;
 670
 671     forAll(InternalContrib, patchI)
 672     {
 673         InternalContrib[patchI] =
 674             cmptMultiply
 675             (
 676                 InternalContrib[patchI],
 677                 psi_.boundaryField()[patchI].patchInternalField()
 678             );
 679     }
 680
 681     FieldField<Field, Type> NeighbourContrib = boundaryCoeffs_;
 682
 683     forAll(NeighbourContrib, patchI)
 684     {
 685         if (psi_.boundaryField()[patchI].coupled())
 686         {
 687             NeighbourContrib[patchI] =
 688                 cmptMultiply
 689                 (
 690                     NeighbourContrib[patchI],
 691                     psi_.boundaryField()[patchI].patchNeighbourField()
 692                 );
 693         }
 694     }
 695
 696     forAll(fieldFlux.boundaryField(), patchI)
 697     {
 698         fieldFlux.boundaryField()[patchI] =
 699             InternalContrib[patchI] - NeighbourContrib[patchI];
 700     }
 701
 702     if (faceFluxCorrectionPtr_)
 703     {
 704         fieldFlux += *faceFluxCorrectionPtr_;
 705     }
 706
 707     return tfieldFlux;
 708 }
 709
 710
 711 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 712
 713 template<class Type>
 714 void faMatrix<Type>::operator=(const faMatrix<Type>& famv)
 715 {
 716     if (this == &famv)
 717     {
 718         FatalErrorIn("faMatrix<Type>::operator=(const faMatrix<Type>&)")
 719             << "attempted to assignment to self"
 720             << abort(FatalError);
 721     }
 722
 723     if (&psi_ != &(famv.psi_))
 724     {
 725         FatalErrorIn("faMatrix<Type>::operator=(const faMatrix<Type>&)")
 726             << "different fields"
 727             << abort(FatalError);
 728     }
 729
 730     lduMatrix::operator=(famv);
 731     source_ = famv.source_;
 732     internalCoeffs_ = famv.internalCoeffs_;
 733     boundaryCoeffs_ = famv.boundaryCoeffs_;
 734
 735     if (faceFluxCorrectionPtr_ && famv.faceFluxCorrectionPtr_)
 736     {
 737         *faceFluxCorrectionPtr_ = *famv.faceFluxCorrectionPtr_;
 738     }
 739     else if (famv.faceFluxCorrectionPtr_)
 740     {
 741         faceFluxCorrectionPtr_ =
 742             new GeometricField<Type, faePatchField, edgeMesh>
 743         (*famv.faceFluxCorrectionPtr_);
 744     }
 745 }
 746
 747
 748 template<class Type>
 749 void faMatrix<Type>::operator=(const tmp<faMatrix<Type> >& tfamv)
 750 {
 751     operator=(tfamv());
 752     tfamv.clear();
 753 }
 754
 755
 756 template<class Type>
 757 void faMatrix<Type>::negate()
 758 {
 759     lduMatrix::negate();
 760     source_.negate();
 761     internalCoeffs_.negate();
 762     boundaryCoeffs_.negate();
 763
 764     if (faceFluxCorrectionPtr_)
 765     {
 766         faceFluxCorrectionPtr_->negate();
 767     }
 768 }
 769
 770
 771 template<class Type>
 772 void faMatrix<Type>::operator+=(const faMatrix<Type>& famv)
 773 {
 774     checkMethod(*this, famv, "+=");
 775
 776     dimensions_ += famv.dimensions_;
 777     lduMatrix::operator+=(famv);
 778     source_ += famv.source_;
 779     internalCoeffs_ += famv.internalCoeffs_;
 780     boundaryCoeffs_ += famv.boundaryCoeffs_;
 781
 782     if (faceFluxCorrectionPtr_ && famv.faceFluxCorrectionPtr_)
 783     {
 784         *faceFluxCorrectionPtr_ += *famv.faceFluxCorrectionPtr_;
 785     }
 786     else if (famv.faceFluxCorrectionPtr_)
 787     {
 788         faceFluxCorrectionPtr_ = new
 789         GeometricField<Type, faePatchField, edgeMesh>
 790         (
 791             *famv.faceFluxCorrectionPtr_
 792         );
 793     }
 794 }
 795
 796
 797 template<class Type>
 798 void faMatrix<Type>::operator+=(const tmp<faMatrix<Type> >& tfamv)
 799 {
 800     operator+=(tfamv());
 801     tfamv.clear();
 802 }
 803
 804
 805 template<class Type>
 806 void faMatrix<Type>::operator-=(const faMatrix<Type>& famv)
 807 {
 808     checkMethod(*this, famv, "+=");
 809
 810     dimensions_ -= famv.dimensions_;
 811     lduMatrix::operator-=(famv);
 812     source_ -= famv.source_;
 813     internalCoeffs_ -= famv.internalCoeffs_;
 814     boundaryCoeffs_ -= famv.boundaryCoeffs_;
 815
 816     if (faceFluxCorrectionPtr_ && famv.faceFluxCorrectionPtr_)
 817     {
 818         *faceFluxCorrectionPtr_ -= *famv.faceFluxCorrectionPtr_;
 819     }
 820     else if (famv.faceFluxCorrectionPtr_)
 821     {
 822         faceFluxCorrectionPtr_ =
 823             new GeometricField<Type, faePatchField, edgeMesh>
 824         (-*famv.faceFluxCorrectionPtr_);
 825     }
 826 }
 827
 828
 829 template<class Type>
 830 void faMatrix<Type>::operator-=(const tmp<faMatrix<Type> >& tfamv)
 831 {
 832     operator-=(tfamv());
 833     tfamv.clear();
 834 }
 835
 836
 837 template<class Type>
 838 void faMatrix<Type>::operator+=
 839 (
 840     const GeometricField<Type, faPatchField, areaMesh>& su
 841 )
 842 {
 843     checkMethod(*this, su, "+=");
 844     source() -= su.mesh().S()*su.internalField();
 845 }
 846
 847
 848 template<class Type>
 849 void faMatrix<Type>::operator+=
 850 (
 851     const tmp<GeometricField<Type, faPatchField, areaMesh> >& tsu
 852 )
 853 {
 854     operator+=(tsu());
 855     tsu.clear();
 856 }
 857
 858
 859 template<class Type>
 860 void faMatrix<Type>::operator-=
 861 (
 862     const GeometricField<Type, faPatchField, areaMesh>& su
 863 )
 864 {
 865     checkMethod(*this, su, "-=");
 866     source() += su.mesh().S()*su.internalField();
 867 }
 868
 869
 870 template<class Type>
 871 void faMatrix<Type>::operator-=
 872 (
 873     const tmp<GeometricField<Type, faPatchField, areaMesh> >& tsu
 874 )
 875 {
 876     operator-=(tsu());
 877     tsu.clear();
 878 }
 879
 880
 881 template<class Type>
 882 void faMatrix<Type>::operator+=
 883 (
 884     const dimensioned<Type>& su
 885 )
 886 {
 887     source() -= su.mesh().S()*su;
 888 }
 889
 890
 891 template<class Type>
 892 void faMatrix<Type>::operator-=
 893 (
 894     const dimensioned<Type>& su
 895 )
 896 {
 897     source() += su.mesh().S()*su;
 898 }
 899
 900
 901 template<class Type>
 902 void faMatrix<Type>::operator*=
 903 (
 904     const areaScalarField& vsf
 905 )
 906 {
 907     dimensions_ *= vsf.dimensions();
 908     lduMatrix::operator*=(vsf.internalField());
 909     source_ *= vsf.internalField();
 910
 911     forAll(boundaryCoeffs_, patchI)
 912     {
 913         scalarField psf = vsf.boundaryField()[patchI].patchInternalField();
 914         internalCoeffs_[patchI] *= psf;
 915         boundaryCoeffs_[patchI] *= psf;
 916     }
 917
 918     if (faceFluxCorrectionPtr_)
 919     {
 920         FatalErrorIn("faMatrix<Type>::operator*=(const tmp<areaScalarField>&)")
 921             << "cannot scale a matrix containing a faceFluxCorrection"
 922             << abort(FatalError);
 923     }
 924 }
 925
 926
 927 template<class Type>
 928 void faMatrix<Type>::operator*=
 929 (
 930     const tmp<areaScalarField>& tvsf
 931 )
 932 {
 933     operator*=(tvsf());
 934     tvsf.clear();
 935 }
 936
 937
 938 template<class Type>
 939 void faMatrix<Type>::operator*=
 940 (
 941     const dimensioned<scalar>& ds
 942 )
 943 {
 944     dimensions_ *= ds.dimensions();
 945     lduMatrix::operator*=(ds.value());
 946     source_ *= ds.value();
 947     internalCoeffs_ *= ds.value();
 948     boundaryCoeffs_ *= ds.value();
 949
 950     if (faceFluxCorrectionPtr_)
 951     {
 952         *faceFluxCorrectionPtr_ *= ds.value();
 953     }
 954 }
 955
 956
 957 // * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
 958
 959 template<class Type>
 960 void checkMethod
 961 (
 962     const faMatrix<Type>& fam1,
 963     const faMatrix<Type>& fam2,
 964     const char* op
 965 )
 966 {
 967     if (&fam1.psi() != &fam2.psi())
 968     {
 969         FatalErrorIn
 970         (
 971             "checkMethod(const faMatrix<Type>&, const faMatrix<Type>&)"
 972         )   << "incompatible fields for operation "
 973             << endl << "    "
 974             << "[" << fam1.psi().name() << "] "
 975             << op
 976             << " [" << fam2.psi().name() << "]"
 977             << abort(FatalError);
 978     }
 979
 980     if (dimensionSet::debug && fam1.dimensions() != fam2.dimensions())
 981     {
 982         FatalErrorIn
 983         (
 984             "checkMethod(const faMatrix<Type>&, const faMatrix<Type>&)"
 985         )   << "incompatible dimensions for operation "
 986             << endl << "    "
 987             << "[" << fam1.psi().name() << fam1.dimensions()/dimArea << " ] "
 988             << op
 989             << " [" << fam2.psi().name() << fam2.dimensions()/dimArea << " ]"
 990             << abort(FatalError);
 991     }
 992 }
 993
 994
 995 template<class Type>
 996 void checkMethod
 997 (
 998     const faMatrix<Type>& fam,
 999     const GeometricField<Type, faPatchField, areaMesh>& vf,
1000     const char* op
1001 )
1002 {
1003     if (dimensionSet::debug && fam.dimensions()/dimArea != vf.dimensions())
1004     {
1005         FatalErrorIn
1006         (
1007             "checkMethod(const faMatrix<Type>&, const GeometricField<Type, "
1008             "faPatchField, areaMesh>&)"
1009         )   <<  "incompatible dimensions for operation "
1010             << endl << "    "
1011             << "[" << fam.psi().name() << fam.dimensions()/dimArea << " ] "
1012             << op
1013             << " [" << vf.name() << vf.dimensions() << " ]"
1014             << abort(FatalError);
1015     }
1016 }
1017
1018
1019 template<class Type>
1020 void checkMethod
1021 (
1022     const faMatrix<Type>& fam,
1023     const dimensioned<Type>& dt,
1024     const char* op
1025 )
1026 {
1027     if (dimensionSet::debug && fam.dimensions()/dimArea != dt.dimensions())
1028     {
1029         FatalErrorIn
1030         (
1031             "checkMethod(const faMatrix<Type>&, const dimensioned<Type>&)"
1032         )   << "incompatible dimensions for operation "
1033             << endl << "    "
1034             << "[" << fam.psi().name() << fam.dimensions()/dimArea << " ] "
1035             << op
1036             << " [" << dt.name() << dt.dimensions() << " ]"
1037             << abort(FatalError);
1038     }
1039 }
1040
1041
1042 template<class Type>
1043 lduSolverPerformance solve
1044 (
1045     faMatrix<Type>& fam,
1046     Istream& solverControls
1047 )
1048 {
1049     return fam.solve(solverControls);
1050 }
1051
1052 template<class Type>
1053 lduSolverPerformance solve
1054 (
1055     const tmp<faMatrix<Type> >& tfam,
1056     Istream& solverControls
1057 )
1058 {
1059     lduSolverPerformance solverPerf =
1060         const_cast<faMatrix<Type>&>(tfam()).solve(solverControls);
1061
1062     tfam.clear();
1063     return solverPerf;
1064 }
1065
1066
1067 template<class Type>
1068 lduSolverPerformance solve(faMatrix<Type>& fam)
1069 {
1070     return fam.solve();
1071 }
1072
1073 template<class Type>
1074 lduSolverPerformance solve(const tmp<faMatrix<Type> >& tfam)
1075 {
1076     lduSolverPerformance solverPerf =
1077         const_cast<faMatrix<Type>&>(tfam()).solve();
1078
1079     tfam.clear();
1080     return solverPerf;
1081 }
1082
1083
1084 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
1085
1086 template<class Type>
1087 tmp<faMatrix<Type> > operator+
1088 (
1089     const faMatrix<Type>& A,
1090     const faMatrix<Type>& B
1091 )
1092 {
1093     checkMethod(A, B, "+");
1094     tmp<faMatrix<Type> > tC(new faMatrix<Type>(A));
1095     tC() += B;
1096     return tC;
1097 }
1098
1099
1100 template<class Type>
1101 tmp<faMatrix<Type> > operator+
1102 (
1103     const tmp<faMatrix<Type> >& tA,
1104     const faMatrix<Type>& B
1105 )
1106 {
1107     checkMethod(tA(), B, "+");
1108     tmp<faMatrix<Type> > tC(tA.ptr());
1109     tC() += B;
1110     return tC;
1111 }
1112
1113
1114 template<class Type>
1115 tmp<faMatrix<Type> > operator+
1116 (
1117     const faMatrix<Type>& A,
1118     const tmp<faMatrix<Type> >& tB
1119 )
1120 {
1121     checkMethod(A, tB(), "+");
1122     tmp<faMatrix<Type> > tC(tB.ptr());
1123     tC() += A;
1124     return tC;
1125 }
1126
1127
1128 template<class Type>
1129 tmp<faMatrix<Type> > operator+
1130 (
1131     const tmp<faMatrix<Type> >& tA,
1132     const tmp<faMatrix<Type> >& tB
1133 )
1134 {
1135     checkMethod(tA(), tB(), "+");
1136     tmp<faMatrix<Type> > tC(tA.ptr());
1137     tC() += tB();
1138     tB.clear();
1139     return tC;
1140 }
1141
1142
1143 template<class Type>
1144 tmp<faMatrix<Type> > operator-
1145 (
1146     const faMatrix<Type>& A
1147 )
1148 {
1149     tmp<faMatrix<Type> > tC(new faMatrix<Type>(A));
1150     tC().negate();
1151     return tC;
1152 }
1153
1154
1155 template<class Type>
1156 tmp<faMatrix<Type> > operator-
1157 (
1158     const tmp<faMatrix<Type> >& tA
1159 )
1160 {
1161     tmp<faMatrix<Type> > tC(tA.ptr());
1162     tC().negate();
1163     return tC;
1164 }
1165
1166
1167 template<class Type>
1168 tmp<faMatrix<Type> > operator-
1169 (
1170     const faMatrix<Type>& A,
1171     const faMatrix<Type>& B
1172 )
1173 {
1174     checkMethod(A, B, "-");
1175     tmp<faMatrix<Type> > tC(new faMatrix<Type>(A));
1176     tC() -= B;
1177     return tC;
1178 }
1179
1180
1181 template<class Type>
1182 tmp<faMatrix<Type> > operator-
1183 (
1184     const tmp<faMatrix<Type> >& tA,
1185     const faMatrix<Type>& B
1186 )
1187 {
1188     checkMethod(tA(), B, "-");
1189     tmp<faMatrix<Type> > tC(tA.ptr());
1190     tC() -= B;
1191     return tC;
1192 }
1193
1194
1195 template<class Type>
1196 tmp<faMatrix<Type> > operator-
1197 (
1198     const faMatrix<Type>& A,
1199     const tmp<faMatrix<Type> >& tB
1200 )
1201 {
1202     checkMethod(A, tB(), "-");
1203     tmp<faMatrix<Type> > tC(tB.ptr());
1204     tC() -= A;
1205     tC().negate();
1206     return tC;
1207 }
1208
1209
1210 template<class Type>
1211 tmp<faMatrix<Type> > operator-
1212 (
1213     const tmp<faMatrix<Type> >& tA,
1214     const tmp<faMatrix<Type> >& tB
1215 )
1216 {
1217     checkMethod(tA(), tB(), "-");
1218     tmp<faMatrix<Type> > tC(tA.ptr());
1219     tC() -= tB();
1220     tB.clear();
1221     return tC;
1222 }
1223
1224
1225 template<class Type>
1226 tmp<faMatrix<Type> > operator==
1227 (
1228     const faMatrix<Type>& A,
1229     const faMatrix<Type>& B
1230 )
1231 {
1232     checkMethod(A, B, "==");
1233     return (A - B);
1234 }
1235
1236
1237 template<class Type>
1238 tmp<faMatrix<Type> > operator==
1239 (
1240     const tmp<faMatrix<Type> >& tA,
1241     const faMatrix<Type>& B
1242 )
1243 {
1244     checkMethod(tA(), B, "==");
1245     return (tA - B);
1246 }
1247
1248
1249 template<class Type>
1250 tmp<faMatrix<Type> > operator==
1251 (
1252     const faMatrix<Type>& A,
1253     const tmp<faMatrix<Type> >& tB
1254 )
1255 {
1256     checkMethod(A, tB(), "==");
1257     return (A - tB);
1258 }
1259
1260
1261 template<class Type>
1262 tmp<faMatrix<Type> > operator==
1263 (
1264     const tmp<faMatrix<Type> >& tA,
1265     const tmp<faMatrix<Type> >& tB
1266 )
1267 {
1268     checkMethod(tA(), tB(), "==");
1269     return (tA - tB);
1270 }
1271
1272
1273 template<class Type>
1274 tmp<faMatrix<Type> > operator+
1275 (
1276     const faMatrix<Type>& A,
1277     const GeometricField<Type, faPatchField, areaMesh>& su
1278 )
1279 {
1280     checkMethod(A, su, "+");
1281     tmp<faMatrix<Type> > tC(new faMatrix<Type>(A));
1282     tC().source() -= su.mesh().S()*su.internalField();
1283     return tC;
1284 }
1285
1286 template<class Type>
1287 tmp<faMatrix<Type> > operator+
1288 (
1289     const tmp<faMatrix<Type> >& tA,
1290     const GeometricField<Type, faPatchField, areaMesh>& su
1291 )
1292 {
1293     checkMethod(tA(), su, "+");
1294     tmp<faMatrix<Type> > tC(tA.ptr());
1295     tC().source() -= su.mesh().S()*su.internalField();
1296     return tC;
1297 }
1298
1299 template<class Type>
1300 tmp<faMatrix<Type> > operator+
1301 (
1302     const faMatrix<Type>& A,
1303     const tmp<GeometricField<Type, faPatchField, areaMesh> >& tsu
1304 )
1305 {
1306     checkMethod(A, tsu(), "+");
1307     tmp<faMatrix<Type> > tC(new faMatrix<Type>(A));
1308     tC().source() -= tsu().mesh().S()*tsu().internalField();
1309     tsu.clear();
1310     return tC;
1311 }
1312
1313
1314 template<class Type>
1315 tmp<faMatrix<Type> > operator+
1316 (
1317     const tmp<faMatrix<Type> >& tA,
1318     const tmp<GeometricField<Type, faPatchField, areaMesh> >& tsu
1319 )
1320 {
1321     checkMethod(tA(), tsu(), "+");
1322     tmp<faMatrix<Type> > tC(tA.ptr());
1323     tC().source() -= tsu().mesh().S()*tsu().internalField();
1324     tsu.clear();
1325     return tC;
1326 }
1327
1328 template<class Type>
1329 tmp<faMatrix<Type> > operator+
1330 (
1331     const GeometricField<Type, faPatchField, areaMesh>& su,
1332     const faMatrix<Type>& A
1333 )
1334 {
1335     checkMethod(A, su, "+");
1336     tmp<faMatrix<Type> > tC(new faMatrix<Type>(A));
1337     tC().source() -= su.mesh().S()*su.internalField();
1338     return tC;
1339 }
1340
1341 template<class Type>
1342 tmp<faMatrix<Type> > operator+
1343 (
1344     const GeometricField<Type, faPatchField, areaMesh>& su,
1345     const tmp<faMatrix<Type> >& tA
1346 )
1347 {
1348     checkMethod(tA(), su, "+");
1349     tmp<faMatrix<Type> > tC(tA.ptr());
1350     tC().source() -= su.mesh().S()*su.internalField();
1351     return tC;
1352 }
1353
1354 template<class Type>
1355 tmp<faMatrix<Type> > operator+
1356 (
1357     const tmp<GeometricField<Type, faPatchField, areaMesh> >& tsu,
1358     const faMatrix<Type>& A
1359 )
1360 {
1361     checkMethod(A, tsu(), "+");
1362     tmp<faMatrix<Type> > tC(new faMatrix<Type>(A));
1363     tC().source() -= tsu().mesh().S()*tsu().internalField();
1364     tsu.clear();
1365     return tC;
1366 }
1367
1368 template<class Type>
1369 tmp<faMatrix<Type> > operator+
1370 (
1371     const tmp<GeometricField<Type, faPatchField, areaMesh> >& tsu,
1372     const tmp<faMatrix<Type> >& tA
1373 )
1374 {
1375     checkMethod(tA(), tsu(), "+");
1376     tmp<faMatrix<Type> > tC(tA.ptr());
1377     tC().source() -= tsu().mesh().S()*tsu().internalField();
1378     tsu.clear();
1379     return tC;
1380 }
1381
1382
1383 template<class Type>
1384 tmp<faMatrix<Type> > operator-
1385 (
1386     const faMatrix<Type>& A,
1387     const GeometricField<Type, faPatchField, areaMesh>& su
1388 )
1389 {
1390     checkMethod(A, su, "-");
1391     tmp<faMatrix<Type> > tC(new faMatrix<Type>(A));
1392     tC().source() += su.mesh().S()*su.internalField();
1393     return tC;
1394 }
1395
1396 template<class Type>
1397 tmp<faMatrix<Type> > operator-
1398 (
1399     const tmp<faMatrix<Type> >& tA,
1400     const GeometricField<Type, faPatchField, areaMesh>& su
1401 )
1402 {
1403     checkMethod(tA(), su, "-");
1404     tmp<faMatrix<Type> > tC(tA.ptr());
1405     tC().source() += su.mesh().S()*su.internalField();
1406     return tC;
1407 }
1408
1409 template<class Type>
1410 tmp<faMatrix<Type> > operator-
1411 (
1412     const faMatrix<Type>& A,
1413     const tmp<GeometricField<Type, faPatchField, areaMesh> >& tsu
1414 )
1415 {
1416     checkMethod(A, tsu(), "-");
1417     tmp<faMatrix<Type> > tC(new faMatrix<Type>(A));
1418     tC().source() += tsu().mesh().S()*tsu().internalField();
1419     tsu.clear();
1420     return tC;
1421 }
1422
1423 template<class Type>
1424 tmp<faMatrix<Type> > operator-
1425 (
1426     const tmp<faMatrix<Type> >& tA,
1427     const tmp<GeometricField<Type, faPatchField, areaMesh> >& tsu
1428 )
1429 {
1430     checkMethod(tA(), tsu(), "-");
1431     tmp<faMatrix<Type> > tC(tA.ptr());
1432     tC().source() += tsu().mesh().S()*tsu().internalField();
1433     tsu.clear();
1434     return tC;
1435 }
1436
1437
1438 template<class Type>
1439 tmp<faMatrix<Type> > operator-
1440 (
1441     const GeometricField<Type, faPatchField, areaMesh>& su,
1442     const faMatrix<Type>& A
1443 )
1444 {
1445     checkMethod(A, su, "-");
1446     tmp<faMatrix<Type> > tC(new faMatrix<Type>(A));
1447     tC().negate();
1448     tC().source() -= su.mesh().S()*su.internalField();
1449     return tC;
1450 }
1451
1452
1453 template<class Type>
1454 tmp<faMatrix<Type> > operator-
1455 (
1456     const GeometricField<Type, faPatchField, areaMesh>& su,
1457     const tmp<faMatrix<Type> >& tA
1458 )
1459 {
1460     checkMethod(tA(), su, "-");
1461     tmp<faMatrix<Type> > tC(tA.ptr());
1462     tC().negate();
1463     tC().source() -= su.mesh().S()*su.internalField();
1464     return tC;
1465 }
1466
1467 template<class Type>
1468 tmp<faMatrix<Type> > operator-
1469 (
1470     const tmp<GeometricField<Type, faPatchField, areaMesh> >& tsu,
1471     const faMatrix<Type>& A
1472 )
1473 {
1474     checkMethod(A, tsu(), "-");
1475     tmp<faMatrix<Type> > tC(new faMatrix<Type>(A));
1476     tC().negate();
1477     tC().source() -= tsu().mesh().S()*tsu().internalField();
1478     tsu.clear();
1479     return tC;
1480 }
1481
1482
1483 template<class Type>
1484 tmp<faMatrix<Type> > operator-
1485 (
1486     const tmp<GeometricField<Type, faPatchField, areaMesh> >& tsu,
1487     const tmp<faMatrix<Type> >& tA
1488 )
1489 {
1490     checkMethod(tA(), tsu(), "-");
1491     tmp<faMatrix<Type> > tC(tA.ptr());
1492     tC().negate();
1493     tC().source() -= tsu().mesh().S()*tsu().internalField();
1494     tsu.clear();
1495     return tC;
1496 }
1497
1498
1499 template<class Type>
1500 tmp<faMatrix<Type> > operator+
1501 (
1502     const faMatrix<Type>& A,
1503     const dimensioned<Type>& su
1504 )
1505 {
1506     checkMethod(A, su, "+");
1507     tmp<faMatrix<Type> > tC(new faMatrix<Type>(A));
1508     tC().source() -= su.value()*A.psi().mesh().S();
1509     return tC;
1510 }
1511
1512
1513 template<class Type>
1514 tmp<faMatrix<Type> > operator+
1515 (
1516     const tmp<faMatrix<Type> >& tA,
1517     const dimensioned<Type>& su
1518 )
1519 {
1520     checkMethod(tA(), su, "+");
1521     tmp<faMatrix<Type> > tC(tA.ptr());
1522     tC().source() -= su.value()*tC().psi().mesh().S();
1523     return tC;
1524 }
1525
1526
1527 template<class Type>
1528 tmp<faMatrix<Type> > operator+
1529 (
1530     const dimensioned<Type>& su,
1531     const faMatrix<Type>& A
1532 )
1533 {
1534     checkMethod(A, su, "+");
1535     tmp<faMatrix<Type> > tC(new faMatrix<Type>(A));
1536     tC().source() -= su.value()*A.psi().mesh().S();
1537     return tC;
1538 }
1539
1540
1541 template<class Type>
1542 tmp<faMatrix<Type> > operator+
1543 (
1544     const dimensioned<Type>& su,
1545     const tmp<faMatrix<Type> >& tA
1546 )
1547 {
1548     checkMethod(tA(), su, "+");
1549     tmp<faMatrix<Type> > tC(tA.ptr());
1550     tC().source() -= su.value()*tC().psi().mesh().S();
1551     return tC;
1552 }
1553
1554
1555 template<class Type>
1556 tmp<faMatrix<Type> > operator-
1557 (
1558     const faMatrix<Type>& A,
1559     const dimensioned<Type>& su
1560 )
1561 {
1562     checkMethod(A, su, "-");
1563     tmp<faMatrix<Type> > tC(new faMatrix<Type>(A));
1564     tC().source() += su.value()*tC().psi().mesh().S();
1565     return tC;
1566 }
1567
1568
1569 template<class Type>
1570 tmp<faMatrix<Type> > operator-
1571 (
1572     const tmp<faMatrix<Type> >& tA,
1573     const dimensioned<Type>& su
1574 )
1575 {
1576     checkMethod(tA(), su, "-");
1577     tmp<faMatrix<Type> > tC(tA.ptr());
1578     tC().source() += su.value()*tC().psi().mesh().S();
1579     return tC;
1580 }
1581
1582
1583 template<class Type>
1584 tmp<faMatrix<Type> > operator-
1585 (
1586     const dimensioned<Type>& su,
1587     const faMatrix<Type>& A
1588 )
1589 {
1590     checkMethod(A, su, "-");
1591     tmp<faMatrix<Type> > tC(new faMatrix<Type>(A));
1592     tC().negate();
1593     tC().source() -= su.value()*A.psi().mesh().S();
1594     return tC;
1595 }
1596
1597
1598 template<class Type>
1599 tmp<faMatrix<Type> > operator-
1600 (
1601     const dimensioned<Type>& su,
1602     const tmp<faMatrix<Type> >& tA
1603 )
1604 {
1605     checkMethod(tA(), su, "-");
1606     tmp<faMatrix<Type> > tC(tA.ptr());
1607     tC().negate();
1608     tC().source() -= su.value()*tC().psi().mesh().S();
1609     return tC;
1610 }
1611
1612
1613 template<class Type>
1614 tmp<faMatrix<Type> > operator==
1615 (
1616     const faMatrix<Type>& A,
1617     const GeometricField<Type, faPatchField, areaMesh>& su
1618 )
1619 {
1620     checkMethod(A, su, "==");
1621     tmp<faMatrix<Type> > tC(new faMatrix<Type>(A));
1622     tC().source() += su.mesh().S()*su.internalField();
1623     return tC;
1624 }
1625
1626 template<class Type>
1627 tmp<faMatrix<Type> > operator==
1628 (
1629     const tmp<faMatrix<Type> >& tA,
1630     const GeometricField<Type, faPatchField, areaMesh>& su
1631 )
1632 {
1633     checkMethod(tA(), su, "==");
1634     tmp<faMatrix<Type> > tC(tA.ptr());
1635     tC().source() += su.mesh().S()*su.internalField();
1636     return tC;
1637 }
1638
1639 template<class Type>
1640 tmp<faMatrix<Type> > operator==
1641 (
1642     const faMatrix<Type>& A,
1643     const tmp<GeometricField<Type, faPatchField, areaMesh> >& tsu
1644 )
1645 {
1646     checkMethod(A, tsu(), "==");
1647     tmp<faMatrix<Type> > tC(new faMatrix<Type>(A));
1648     tC().source() += tsu().mesh().S()*tsu().internalField();
1649     tsu.clear();
1650     return tC;
1651 }
1652
1653 template<class Type>
1654 tmp<faMatrix<Type> > operator==
1655 (
1656     const tmp<faMatrix<Type> >& tA,
1657     const tmp<GeometricField<Type, faPatchField, areaMesh> >& tsu
1658 )
1659 {
1660     checkMethod(tA(), tsu(), "==");
1661     tmp<faMatrix<Type> > tC(tA.ptr());
1662     tC().source() += tsu().mesh().S()*tsu().internalField();
1663     tsu.clear();
1664     return tC;
1665 }
1666
1667
1668 template<class Type>
1669 tmp<faMatrix<Type> > operator==
1670 (
1671     const faMatrix<Type>& A,
1672     const dimensioned<Type>& su
1673 )
1674 {
1675     checkMethod(A, su, "==");
1676     tmp<faMatrix<Type> > tC(new faMatrix<Type>(A));
1677     tC().source() += A.psi().mesh().S()*su.value();
1678     return tC;
1679 }
1680
1681
1682 template<class Type>
1683 tmp<faMatrix<Type> > operator==
1684 (
1685     const tmp<faMatrix<Type> >& tA,
1686     const dimensioned<Type>& su
1687 )
1688 {
1689     checkMethod(tA(), su, "==");
1690     tmp<faMatrix<Type> > tC(tA.ptr());
1691     tC().source() += tC().psi().mesh().S()*su.value();
1692     return tC;
1693 }
1694
1695
1696 template<class Type>
1697 tmp<faMatrix<Type> > operator*
1698 (
1699     const areaScalarField& vsf,
1700     const faMatrix<Type>& A
1701 )
1702 {
1703     tmp<faMatrix<Type> > tC(new faMatrix<Type>(A));
1704     tC() *= vsf;
1705     return tC;
1706 }
1707
1708 template<class Type>
1709 tmp<faMatrix<Type> > operator*
1710 (
1711     const tmp<areaScalarField>& tvsf,
1712     const faMatrix<Type>& A
1713 )
1714 {
1715     tmp<faMatrix<Type> > tC(new faMatrix<Type>(A));
1716     tC() *= tvsf;
1717     return tC;
1718 }
1719
1720 template<class Type>
1721 tmp<faMatrix<Type> > operator*
1722 (
1723     const areaScalarField& vsf,
1724     const tmp<faMatrix<Type> >& tA
1725 )
1726 {
1727     tmp<faMatrix<Type> > tC(tA.ptr());
1728     tC() *= vsf;
1729     return tC;
1730 }
1731
1732 template<class Type>
1733 tmp<faMatrix<Type> > operator*
1734 (
1735     const tmp<areaScalarField>& tvsf,
1736     const tmp<faMatrix<Type> >& tA
1737 )
1738 {
1739     tmp<faMatrix<Type> > tC(tA.ptr());
1740     tC() *= tvsf;
1741     return tC;
1742 }
1743
1744
1745 template<class Type>
1746 tmp<faMatrix<Type> > operator*
1747 (
1748     const dimensioned<scalar>& ds,
1749     const faMatrix<Type>& A
1750 )
1751 {
1752     tmp<faMatrix<Type> > tC(new faMatrix<Type>(A));
1753     tC() *= ds;
1754     return tC;
1755 }
1756
1757
1758 template<class Type>
1759 tmp<faMatrix<Type> > operator*
1760 (
1761     const dimensioned<scalar>& ds,
1762     const tmp<faMatrix<Type> >& tA
1763 )
1764 {
1765     tmp<faMatrix<Type> > tC(tA.ptr());
1766     tC() *= ds;
1767     return tC;
1768 }
1769
1770
1771 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
1772
1773 template<class Type>
1774 Ostream& operator<<(Ostream& os, const faMatrix<Type>& fam)
1775 {
1776     os  << static_cast<const lduMatrix&>(fam) << nl
1777         << fam.dimensions_ << nl
1778         << fam.source_ << nl
1779         << fam.internalCoeffs_ << nl
1780         << fam.boundaryCoeffs_ << endl;
1781
1782     os.check("Ostream& operator<<(Ostream&, faMatrix<Type>&");
1783
1784     return os;
1785 }
1786
1787
1788 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1789
1790 } // End namespace Foam
1791
1792 // * * * * * * * * * * * * * * * * Solvers * * * * * * * * * * * * * * * * * //
1793
1794 #include "faMatrixSolve.C"
1795
1796 // ************************************************************************* //