src/finiteArea/faMatrices/faMatrix/faMatrix.C

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