src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
   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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your 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, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "volFields.H"
  27 #include "surfaceFields.H"
  28 #include "calculatedFvPatchFields.H"
  29 #include "zeroGradientFvPatchFields.H"
  30 #include "coupledFvPatchFields.H"
  31
  32 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  33
  34 template<class Type>
  35 template<class Type2>
  36 void Foam::fvMatrix<Type>::addToInternalField
  37 (
  38     const unallocLabelList& addr,
  39     const Field<Type2>& pf,
  40     Field<Type2>& intf
  41 ) const
  42 {
  43     if (addr.size() != pf.size())
  44     {
  45         FatalErrorIn
  46         (
  47             "fvMatrix<Type>::addToInternalField(const unallocLabelList&, "
  48             "const Field&, Field&)"
  49         )   << "sizes of addressing and field are different"
  50             << abort(FatalError);
  51     }
  52
  53     forAll(addr, faceI)
  54     {
  55         intf[addr[faceI]] += pf[faceI];
  56     }
  57 }
  58
  59
  60 template<class Type>
  61 template<class Type2>
  62 void Foam::fvMatrix<Type>::addToInternalField
  63 (
  64     const unallocLabelList& addr,
  65     const tmp<Field<Type2> >& tpf,
  66     Field<Type2>& intf
  67 ) const
  68 {
  69     addToInternalField(addr, tpf(), intf);
  70     tpf.clear();
  71 }
  72
  73
  74 template<class Type>
  75 template<class Type2>
  76 void Foam::fvMatrix<Type>::subtractFromInternalField
  77 (
  78     const unallocLabelList& addr,
  79     const Field<Type2>& pf,
  80     Field<Type2>& intf
  81 ) const
  82 {
  83     if (addr.size() != pf.size())
  84     {
  85         FatalErrorIn
  86         (
  87             "fvMatrix<Type>::addToInternalField(const unallocLabelList&, "
  88             "const Field&, Field&)"
  89         )   << "sizes of addressing and field are different"
  90             << abort(FatalError);
  91     }
  92
  93     forAll(addr, faceI)
  94     {
  95         intf[addr[faceI]] -= pf[faceI];
  96     }
  97 }
  98
  99
 100 template<class Type>
 101 template<class Type2>
 102 void Foam::fvMatrix<Type>::subtractFromInternalField
 103 (
 104     const unallocLabelList& addr,
 105     const tmp<Field<Type2> >& tpf,
 106     Field<Type2>& intf
 107 ) const
 108 {
 109     subtractFromInternalField(addr, tpf(), intf);
 110     tpf.clear();
 111 }
 112
 113
 114 template<class Type>
 115 void Foam::fvMatrix<Type>::addBoundaryDiag
 116 (
 117     scalarField& diag,
 118     const direction solveCmpt
 119 ) const
 120 {
 121     forAll(internalCoeffs_, patchI)
 122     {
 123         addToInternalField
 124         (
 125             lduAddr().patchAddr(patchI),
 126             internalCoeffs_[patchI].component(solveCmpt),
 127             diag
 128         );
 129     }
 130 }
 131
 132
 133 template<class Type>
 134 void Foam::fvMatrix<Type>::addCmptAvBoundaryDiag(scalarField& diag) const
 135 {
 136     forAll(internalCoeffs_, patchI)
 137     {
 138         addToInternalField
 139         (
 140             lduAddr().patchAddr(patchI),
 141             cmptAv(internalCoeffs_[patchI]),
 142             diag
 143         );
 144     }
 145 }
 146
 147
 148 template<class Type>
 149 void Foam::fvMatrix<Type>::addBoundarySource
 150 (
 151     Field<Type>& source,
 152     const bool couples
 153 ) const
 154 {
 155     forAll(psi_.boundaryField(), patchI)
 156     {
 157         const fvPatchField<Type>& ptf = psi_.boundaryField()[patchI];
 158         const Field<Type>& pbc = boundaryCoeffs_[patchI];
 159
 160         if (!ptf.coupled())
 161         {
 162             addToInternalField(lduAddr().patchAddr(patchI), pbc, source);
 163         }
 164         else if (couples)
 165         {
 166             tmp<Field<Type> > tpnf = ptf.patchNeighbourField();
 167             const Field<Type>& pnf = tpnf();
 168
 169             const unallocLabelList& addr = lduAddr().patchAddr(patchI);
 170
 171             forAll(addr, facei)
 172             {
 173                 source[addr[facei]] += cmptMultiply(pbc[facei], pnf[facei]);
 174             }
 175         }
 176     }
 177 }
 178
 179
 180 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 181
 182 template<class Type>
 183 Foam::fvMatrix<Type>::fvMatrix
 184 (
 185     GeometricField<Type, fvPatchField, volMesh>& psi,
 186     const dimensionSet& ds
 187 )
 188 :
 189     lduMatrix(psi.mesh()),
 190     psi_(psi),
 191     dimensions_(ds),
 192     source_(psi.size(), pTraits<Type>::zero),
 193     internalCoeffs_(psi.mesh().boundary().size()),
 194     boundaryCoeffs_(psi.mesh().boundary().size()),
 195     faceFluxCorrectionPtr_(NULL)
 196 {
 197     if (debug)
 198     {
 199         Info<< "fvMatrix<Type>(GeometricField<Type, fvPatchField, volMesh>&,"
 200                " const dimensionSet&) : "
 201                "constructing fvMatrix<Type> for field " << psi_.name()
 202             << endl;
 203     }
 204
 205     // Initialise coupling coefficients
 206     forAll (psi.mesh().boundary(), patchI)
 207     {
 208         internalCoeffs_.set
 209         (
 210             patchI,
 211             new Field<Type>
 212             (
 213                 psi.mesh().boundary()[patchI].size(),
 214                 pTraits<Type>::zero
 215             )
 216         );
 217
 218         boundaryCoeffs_.set
 219         (
 220             patchI,
 221             new Field<Type>
 222             (
 223                 psi.mesh().boundary()[patchI].size(),
 224                 pTraits<Type>::zero
 225             )
 226         );
 227     }
 228
 229     psi_.boundaryField().updateCoeffs();
 230 }
 231
 232
 233 template<class Type>
 234 Foam::fvMatrix<Type>::fvMatrix(const fvMatrix<Type>& fvm)
 235 :
 236     refCount(),
 237     lduMatrix(fvm),
 238     psi_(fvm.psi_),
 239     dimensions_(fvm.dimensions_),
 240     source_(fvm.source_),
 241     internalCoeffs_(fvm.internalCoeffs_),
 242     boundaryCoeffs_(fvm.boundaryCoeffs_),
 243     faceFluxCorrectionPtr_(NULL)
 244 {
 245     if (debug)
 246     {
 247         Info<< "fvMatrix<Type>::fvMatrix(const fvMatrix<Type>&) : "
 248             << "copying fvMatrix<Type> for field " << psi_.name()
 249             << endl;
 250     }
 251
 252     if (fvm.faceFluxCorrectionPtr_)
 253     {
 254         faceFluxCorrectionPtr_ = new
 255         GeometricField<Type, fvsPatchField, surfaceMesh>
 256         (
 257             *(fvm.faceFluxCorrectionPtr_)
 258         );
 259     }
 260 }
 261
 262
 263 #ifdef ConstructFromTmp
 264 template<class Type>
 265 Foam::fvMatrix<Type>::fvMatrix(const tmp<fvMatrix<Type> >& tfvm)
 266 :
 267     refCount(),
 268     lduMatrix
 269     (
 270         const_cast<fvMatrix<Type>&>(tfvm()),
 271         tfvm.isTmp()
 272     ),
 273     psi_(tfvm().psi_),
 274     dimensions_(tfvm().dimensions_),
 275     source_
 276     (
 277         const_cast<fvMatrix<Type>&>(tfvm()).source_,
 278         tfvm.isTmp()
 279     ),
 280     internalCoeffs_
 281     (
 282         const_cast<fvMatrix<Type>&>(tfvm()).internalCoeffs_,
 283         tfvm.isTmp()
 284     ),
 285     boundaryCoeffs_
 286     (
 287         const_cast<fvMatrix<Type>&>(tfvm()).boundaryCoeffs_,
 288         tfvm.isTmp()
 289     ),
 290     faceFluxCorrectionPtr_(NULL)
 291 {
 292     if (debug)
 293     {
 294         Info<< "fvMatrix<Type>::fvMatrix(const tmp<fvMatrix<Type> >&) : "
 295             << "copying fvMatrix<Type> for field " << psi_.name()
 296             << endl;
 297     }
 298
 299     if (tfvm().faceFluxCorrectionPtr_)
 300     {
 301         if (tfvm.isTmp())
 302         {
 303             faceFluxCorrectionPtr_ = tfvm().faceFluxCorrectionPtr_;
 304             tfvm().faceFluxCorrectionPtr_ = NULL;
 305         }
 306         else
 307         {
 308             faceFluxCorrectionPtr_ = new
 309                 GeometricField<Type, fvsPatchField, surfaceMesh>
 310                 (
 311                     *(tfvm().faceFluxCorrectionPtr_)
 312                 );
 313         }
 314     }
 315
 316     tfvm.clear();
 317 }
 318 #endif
 319
 320
 321 template<class Type>
 322 Foam::fvMatrix<Type>::fvMatrix
 323 (
 324     GeometricField<Type, fvPatchField, volMesh>& psi,
 325     Istream& is
 326 )
 327 :
 328     lduMatrix(psi.mesh()),
 329     psi_(psi),
 330     dimensions_(is),
 331     source_(is),
 332     internalCoeffs_(psi.mesh().boundary().size()),
 333     boundaryCoeffs_(psi.mesh().boundary().size()),
 334     faceFluxCorrectionPtr_(NULL)
 335 {
 336     if (debug)
 337     {
 338         Info<< "fvMatrix<Type>"
 339                "(GeometricField<Type, fvPatchField, volMesh>&, Istream&) : "
 340                "constructing fvMatrix<Type> for field " << psi_.name()
 341             << endl;
 342     }
 343
 344     // Initialise coupling coefficients
 345     forAll (psi.mesh().boundary(), patchI)
 346     {
 347         internalCoeffs_.set
 348         (
 349             patchI,
 350             new Field<Type>
 351             (
 352                 psi.mesh().boundary()[patchI].size(),
 353                 pTraits<Type>::zero
 354             )
 355         );
 356
 357         boundaryCoeffs_.set
 358         (
 359             patchI,
 360             new Field<Type>
 361             (
 362                 psi.mesh().boundary()[patchI].size(),
 363                 pTraits<Type>::zero
 364             )
 365         );
 366     }
 367
 368 }
 369
 370
 371 template<class Type>
 372 Foam::fvMatrix<Type>::~fvMatrix()
 373 {
 374     if (debug)
 375     {
 376         Info<< "fvMatrix<Type>::~fvMatrix<Type>() : "
 377             << "destroying fvMatrix<Type> for field " << psi_.name()
 378             << endl;
 379     }
 380
 381     if (faceFluxCorrectionPtr_)
 382     {
 383         delete faceFluxCorrectionPtr_;
 384     }
 385 }
 386
 387
 388 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 389
 390 // Set solution in given cells and eliminate corresponding
 391 // equations from the matrix
 392 template<class Type>
 393 void Foam::fvMatrix<Type>::setValues
 394 (
 395     const labelList& cellLabels,
 396     const Field<Type>& values
 397 )
 398 {
 399     const fvMesh& mesh = psi_.mesh();
 400
 401     const cellList& cells = mesh.cells();
 402     const unallocLabelList& own = mesh.owner();
 403     const unallocLabelList& nei = mesh.neighbour();
 404
 405     scalarField& Diag = diag();
 406
 407     forAll(cellLabels, i)
 408     {
 409         label celli = cellLabels[i];
 410
 411         psi_[celli] = values[i];
 412         source_[celli] = values[i]*Diag[celli];
 413
 414         if (symmetric() || asymmetric())
 415         {
 416             const cell& c = cells[celli];
 417
 418             forAll(c, j)
 419             {
 420                 label facei = c[j];
 421
 422                 if (mesh.isInternalFace(facei))
 423                 {
 424                     if (symmetric())
 425                     {
 426                         if (celli == own[facei])
 427                         {
 428                             source_[nei[facei]] -= upper()[facei]*values[i];
 429                         }
 430                         else
 431                         {
 432                             source_[own[facei]] -= upper()[facei]*values[i];
 433                         }
 434
 435                         upper()[facei] = 0.0;
 436                     }
 437                     else
 438                     {
 439                         if (celli == own[facei])
 440                         {
 441                             source_[nei[facei]] -= lower()[facei]*values[i];
 442                         }
 443                         else
 444                         {
 445                             source_[own[facei]] -= upper()[facei]*values[i];
 446                         }
 447
 448                         upper()[facei] = 0.0;
 449                         lower()[facei] = 0.0;
 450                     }
 451                 }
 452                 else
 453                 {
 454                     label patchi = mesh.boundaryMesh().whichPatch(facei);
 455
 456                     if (internalCoeffs_[patchi].size())
 457                     {
 458                         label patchFacei =
 459                             mesh.boundaryMesh()[patchi].whichFace(facei);
 460
 461                         internalCoeffs_[patchi][patchFacei] =
 462                             pTraits<Type>::zero;
 463
 464                         boundaryCoeffs_[patchi][patchFacei] =
 465                             pTraits<Type>::zero;
 466                     }
 467                 }
 468             }
 469         }
 470     }
 471 }
 472
 473
 474 template<class Type>
 475 void Foam::fvMatrix<Type>::setReference
 476 (
 477     const label celli,
 478     const Type& value,
 479     const bool forceReference
 480 )
 481 {
 482     if ((forceReference || psi_.needReference()) && celli >= 0)
 483     {
 484         source()[celli] += diag()[celli]*value;
 485         diag()[celli] += diag()[celli];
 486     }
 487 }
 488
 489
 490 template<class Type>
 491 void Foam::fvMatrix<Type>::relax(const scalar alpha)
 492 {
 493     if (alpha <= 0)
 494     {
 495         return;
 496     }
 497
 498     Field<Type>& S = source();
 499     scalarField& D = diag();
 500
 501     // Store the current unrelaxed diagonal for use in updating the source
 502     scalarField D0(D);
 503
 504     // Calculate the sum-mag off-diagonal from the interior faces
 505     scalarField sumOff(D.size(), 0.0);
 506     sumMagOffDiag(sumOff);
 507
 508     // Handle the boundary contributions to the diagonal
 509     forAll(psi_.boundaryField(), patchI)
 510     {
 511         const fvPatchField<Type>& ptf = psi_.boundaryField()[patchI];
 512
 513         if (ptf.size())
 514         {
 515             const unallocLabelList& pa = lduAddr().patchAddr(patchI);
 516             Field<Type>& iCoeffs = internalCoeffs_[patchI];
 517
 518             if (ptf.coupled())
 519             {
 520                 const Field<Type>& pCoeffs = boundaryCoeffs_[patchI];
 521
 522                 // For coupled boundaries add the diagonal and
 523                 // off-diagonal contributions
 524                 forAll(pa, face)
 525                 {
 526                     D[pa[face]] += component(iCoeffs[face], 0);
 527                     sumOff[pa[face]] += mag(component(pCoeffs[face], 0));
 528                 }
 529             }
 530             else
 531             {
 532                 // For non-coupled boundaries subtract the diagonal
 533                 // contribution off-diagonal sum which avoids having to remove
 534                 // it from the diagonal later.
 535                 // Also add the source contribution from the relaxation
 536                 forAll(pa, face)
 537                 {
 538                     Type iCoeff0 = iCoeffs[face];
 539                     iCoeffs[face] = cmptMag(iCoeffs[face]);
 540                     sumOff[pa[face]] -= cmptMin(iCoeffs[face]);
 541                     iCoeffs[face] /= alpha;
 542                     S[pa[face]] +=
 543                         cmptMultiply(iCoeffs[face] - iCoeff0, psi_[pa[face]]);
 544                 }
 545             }
 546         }
 547     }
 548
 549     // Ensure the matrix is diagonally dominant...
 550     max(D, D, sumOff);
 551
 552     // ... then relax
 553     D /= alpha;
 554
 555     // Now remove the diagonal contribution from coupled boundaries
 556     forAll(psi_.boundaryField(), patchI)
 557     {
 558         const fvPatchField<Type>& ptf = psi_.boundaryField()[patchI];
 559
 560         if (ptf.size())
 561         {
 562             const unallocLabelList& pa = lduAddr().patchAddr(patchI);
 563             Field<Type>& iCoeffs = internalCoeffs_[patchI];
 564
 565             if (ptf.coupled())
 566             {
 567                 forAll(pa, face)
 568                 {
 569                     D[pa[face]] -= component(iCoeffs[face], 0);
 570                 }
 571             }
 572         }
 573     }
 574
 575     // Finally add the relaxation contribution to the source.
 576     S += (D - D0)*psi_.internalField();
 577 }
 578
 579
 580 template<class Type>
 581 void Foam::fvMatrix<Type>::relax()
 582 {
 583     if (psi_.mesh().relax(psi_.name()))
 584     {
 585         relax(psi_.mesh().relaxationFactor(psi_.name()));
 586     }
 587 }
 588
 589
 590 template<class Type>
 591 void Foam::fvMatrix<Type>::boundaryManipulate
 592 (
 593     typename GeometricField<Type, fvPatchField, volMesh>::
 594         GeometricBoundaryField& bFields
 595 )
 596 {
 597     forAll(bFields, patchI)
 598     {
 599         bFields[patchI].manipulateMatrix(*this);
 600     }
 601 }
 602
 603
 604 template<class Type>
 605 Foam::tmp<Foam::scalarField> Foam::fvMatrix<Type>::D() const
 606 {
 607     tmp<scalarField> tdiag(new scalarField(diag()));
 608     addCmptAvBoundaryDiag(tdiag());
 609     return tdiag;
 610 }
 611
 612
 613 template<class Type>
 614 Foam::tmp<Foam::Field<Type> > Foam::fvMatrix<Type>::DD() const
 615 {
 616     tmp<Field<Type> > tdiag(pTraits<Type>::one*diag());
 617
 618     forAll(psi_.boundaryField(), patchI)
 619     {
 620         const fvPatchField<Type>& ptf = psi_.boundaryField()[patchI];
 621
 622         if (!ptf.coupled() && ptf.size())
 623         {
 624             addToInternalField
 625             (
 626                 lduAddr().patchAddr(patchI),
 627                 internalCoeffs_[patchI],
 628                 tdiag()
 629             );
 630         }
 631     }
 632
 633     return tdiag;
 634 }
 635
 636
 637 template<class Type>
 638 Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Type>::A() const
 639 {
 640     tmp<volScalarField> tAphi
 641     (
 642         new volScalarField
 643         (
 644             IOobject
 645             (
 646                 "A("+psi_.name()+')',
 647                 psi_.instance(),
 648                 psi_.mesh(),
 649                 IOobject::NO_READ,
 650                 IOobject::NO_WRITE
 651             ),
 652             psi_.mesh(),
 653             dimensions_/psi_.dimensions()/dimVol,
 654             zeroGradientFvPatchScalarField::typeName
 655         )
 656     );
 657
 658     tAphi().internalField() = D()/psi_.mesh().V();
 659     tAphi().correctBoundaryConditions();
 660
 661     return tAphi;
 662 }
 663
 664
 665 template<class Type>
 666 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
 667 Foam::fvMatrix<Type>::H() const
 668 {
 669     tmp<GeometricField<Type, fvPatchField, volMesh> > tHphi
 670     (
 671         new GeometricField<Type, fvPatchField, volMesh>
 672         (
 673             IOobject
 674             (
 675                 "H("+psi_.name()+')',
 676                 psi_.instance(),
 677                 psi_.mesh(),
 678                 IOobject::NO_READ,
 679                 IOobject::NO_WRITE
 680             ),
 681             psi_.mesh(),
 682             dimensions_/dimVol,
 683             zeroGradientFvPatchScalarField::typeName
 684         )
 685     );
 686     GeometricField<Type, fvPatchField, volMesh>& Hphi = tHphi();
 687
 688     // Loop over field components
 689     for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
 690     {
 691         scalarField psiCmpt = psi_.internalField().component(cmpt);
 692
 693         scalarField boundaryDiagCmpt(psi_.size(), 0.0);
 694         addBoundaryDiag(boundaryDiagCmpt, cmpt);
 695         boundaryDiagCmpt.negate();
 696         addCmptAvBoundaryDiag(boundaryDiagCmpt);
 697
 698         Hphi.internalField().replace(cmpt, boundaryDiagCmpt*psiCmpt);
 699     }
 700
 701     Hphi.internalField() += lduMatrix::H(psi_.internalField()) + source_;
 702     addBoundarySource(Hphi.internalField());
 703
 704     Hphi.internalField() /= psi_.mesh().V();
 705     Hphi.correctBoundaryConditions();
 706
 707     typename Type::labelType validComponents
 708     (
 709         pow
 710         (
 711             psi_.mesh().solutionD(),
 712             pTraits<typename powProduct<Vector<label>, Type::rank>::type>::zero
 713         )
 714     );
 715
 716     for(direction cmpt=0; cmpt<Type::nComponents; cmpt++)
 717     {
 718         if (validComponents[cmpt] == -1)
 719         {
 720             Hphi.replace
 721             (
 722                 cmpt,
 723                 dimensionedScalar("0", Hphi.dimensions(), 0.0)
 724             );
 725         }
 726     }
 727
 728     return tHphi;
 729 }
 730
 731
 732 template<class Type>
 733 Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Type>::H1() const
 734 {
 735     tmp<volScalarField> tH1
 736     (
 737         new volScalarField
 738         (
 739             IOobject
 740             (
 741                 "H(1)",
 742                 psi_.instance(),
 743                 psi_.mesh(),
 744                 IOobject::NO_READ,
 745                 IOobject::NO_WRITE
 746             ),
 747             psi_.mesh(),
 748             dimensions_/(dimVol*psi_.dimensions()),
 749             zeroGradientFvPatchScalarField::typeName
 750         )
 751     );
 752     volScalarField& H1_ = tH1();
 753
 754     // Loop over field components
 755     /*
 756     for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
 757     {
 758         scalarField psiCmpt = psi_.internalField().component(cmpt);
 759
 760         scalarField boundaryDiagCmpt(psi_.size(), 0.0);
 761         addBoundaryDiag(boundaryDiagCmpt, cmpt);
 762         boundaryDiagCmpt.negate();
 763         addCmptAvBoundaryDiag(boundaryDiagCmpt);
 764
 765         H1_.internalField().replace(cmpt, boundaryDiagCmpt);
 766     }
 767
 768     H1_.internalField() += lduMatrix::H1();
 769     */
 770
 771     H1_.internalField() = lduMatrix::H1();
 772
 773     H1_.internalField() /= psi_.mesh().V();
 774     H1_.correctBoundaryConditions();
 775
 776     return tH1;
 777 }
 778
 779
 780 template<class Type>
 781 Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> >
 782 Foam::fvMatrix<Type>::
 783 flux() const
 784 {
 785     if (!psi_.mesh().fluxRequired(psi_.name()))
 786     {
 787         FatalErrorIn("fvMatrix<Type>::flux()")
 788             << "flux requested but " << psi_.name()
 789             << " not specified in the fluxRequired sub-dictionary"
 790                " of fvSchemes."
 791             << abort(FatalError);
 792     }
 793
 794     // construct GeometricField<Type, fvsPatchField, surfaceMesh>
 795     tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tfieldFlux
 796     (
 797         new GeometricField<Type, fvsPatchField, surfaceMesh>
 798         (
 799             IOobject
 800             (
 801                 "flux("+psi_.name()+')',
 802                 psi_.instance(),
 803                 psi_.mesh(),
 804                 IOobject::NO_READ,
 805                 IOobject::NO_WRITE
 806             ),
 807             psi_.mesh(),
 808             dimensions()
 809         )
 810     );
 811     GeometricField<Type, fvsPatchField, surfaceMesh>& fieldFlux = tfieldFlux();
 812
 813     for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
 814     {
 815         fieldFlux.internalField().replace
 816         (
 817             cmpt,
 818             lduMatrix::faceH(psi_.internalField().component(cmpt))
 819         );
 820     }
 821
 822     FieldField<Field, Type> InternalContrib = internalCoeffs_;
 823
 824     forAll(InternalContrib, patchI)
 825     {
 826         InternalContrib[patchI] =
 827             cmptMultiply
 828             (
 829                 InternalContrib[patchI],
 830                 psi_.boundaryField()[patchI].patchInternalField()
 831             );
 832     }
 833
 834     FieldField<Field, Type> NeighbourContrib = boundaryCoeffs_;
 835
 836     forAll(NeighbourContrib, patchI)
 837     {
 838         if (psi_.boundaryField()[patchI].coupled())
 839         {
 840             NeighbourContrib[patchI] =
 841                 cmptMultiply
 842                 (
 843                     NeighbourContrib[patchI],
 844                     psi_.boundaryField()[patchI].patchNeighbourField()
 845                 );
 846         }
 847     }
 848
 849     forAll(fieldFlux.boundaryField(), patchI)
 850     {
 851         fieldFlux.boundaryField()[patchI] =
 852             InternalContrib[patchI] - NeighbourContrib[patchI];
 853     }
 854
 855     if (faceFluxCorrectionPtr_)
 856     {
 857         fieldFlux += *faceFluxCorrectionPtr_;
 858     }
 859
 860     return tfieldFlux;
 861 }
 862
 863
 864 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 865
 866 template<class Type>
 867 void Foam::fvMatrix<Type>::operator=(const fvMatrix<Type>& fvmv)
 868 {
 869     if (this == &fvmv)
 870     {
 871         FatalErrorIn("fvMatrix<Type>::operator=(const fvMatrix<Type>&)")
 872             << "attempted assignment to self"
 873             << abort(FatalError);
 874     }
 875
 876     if (&psi_ != &(fvmv.psi_))
 877     {
 878         FatalErrorIn("fvMatrix<Type>::operator=(const fvMatrix<Type>&)")
 879             << "different fields"
 880             << abort(FatalError);
 881     }
 882
 883     lduMatrix::operator=(fvmv);
 884     source_ = fvmv.source_;
 885     internalCoeffs_ = fvmv.internalCoeffs_;
 886     boundaryCoeffs_ = fvmv.boundaryCoeffs_;
 887
 888     if (faceFluxCorrectionPtr_ && fvmv.faceFluxCorrectionPtr_)
 889     {
 890         *faceFluxCorrectionPtr_ = *fvmv.faceFluxCorrectionPtr_;
 891     }
 892     else if (fvmv.faceFluxCorrectionPtr_)
 893     {
 894         faceFluxCorrectionPtr_ =
 895             new GeometricField<Type, fvsPatchField, surfaceMesh>
 896         (*fvmv.faceFluxCorrectionPtr_);
 897     }
 898 }
 899
 900
 901 template<class Type>
 902 void Foam::fvMatrix<Type>::operator=(const tmp<fvMatrix<Type> >& tfvmv)
 903 {
 904     operator=(tfvmv());
 905     tfvmv.clear();
 906 }
 907
 908
 909 template<class Type>
 910 void Foam::fvMatrix<Type>::negate()
 911 {
 912     lduMatrix::negate();
 913     source_.negate();
 914     internalCoeffs_.negate();
 915     boundaryCoeffs_.negate();
 916
 917     if (faceFluxCorrectionPtr_)
 918     {
 919         faceFluxCorrectionPtr_->negate();
 920     }
 921 }
 922
 923
 924 template<class Type>
 925 void Foam::fvMatrix<Type>::operator+=(const fvMatrix<Type>& fvmv)
 926 {
 927     checkMethod(*this, fvmv, "+=");
 928
 929     dimensions_ += fvmv.dimensions_;
 930     lduMatrix::operator+=(fvmv);
 931     source_ += fvmv.source_;
 932     internalCoeffs_ += fvmv.internalCoeffs_;
 933     boundaryCoeffs_ += fvmv.boundaryCoeffs_;
 934
 935     if (faceFluxCorrectionPtr_ && fvmv.faceFluxCorrectionPtr_)
 936     {
 937         *faceFluxCorrectionPtr_ += *fvmv.faceFluxCorrectionPtr_;
 938     }
 939     else if (fvmv.faceFluxCorrectionPtr_)
 940     {
 941         faceFluxCorrectionPtr_ = new
 942         GeometricField<Type, fvsPatchField, surfaceMesh>
 943         (
 944             *fvmv.faceFluxCorrectionPtr_
 945         );
 946     }
 947 }
 948
 949
 950 template<class Type>
 951 void Foam::fvMatrix<Type>::operator+=(const tmp<fvMatrix<Type> >& tfvmv)
 952 {
 953     operator+=(tfvmv());
 954     tfvmv.clear();
 955 }
 956
 957
 958 template<class Type>
 959 void Foam::fvMatrix<Type>::operator-=(const fvMatrix<Type>& fvmv)
 960 {
 961     checkMethod(*this, fvmv, "+=");
 962
 963     dimensions_ -= fvmv.dimensions_;
 964     lduMatrix::operator-=(fvmv);
 965     source_ -= fvmv.source_;
 966     internalCoeffs_ -= fvmv.internalCoeffs_;
 967     boundaryCoeffs_ -= fvmv.boundaryCoeffs_;
 968
 969     if (faceFluxCorrectionPtr_ && fvmv.faceFluxCorrectionPtr_)
 970     {
 971         *faceFluxCorrectionPtr_ -= *fvmv.faceFluxCorrectionPtr_;
 972     }
 973     else if (fvmv.faceFluxCorrectionPtr_)
 974     {
 975         faceFluxCorrectionPtr_ =
 976             new GeometricField<Type, fvsPatchField, surfaceMesh>
 977         (-*fvmv.faceFluxCorrectionPtr_);
 978     }
 979 }
 980
 981
 982 template<class Type>
 983 void Foam::fvMatrix<Type>::operator-=(const tmp<fvMatrix<Type> >& tfvmv)
 984 {
 985     operator-=(tfvmv());
 986     tfvmv.clear();
 987 }
 988
 989
 990 template<class Type>
 991 void Foam::fvMatrix<Type>::operator+=
 992 (
 993     const DimensionedField<Type, volMesh>& su
 994 )
 995 {
 996     checkMethod(*this, su, "+=");
 997     source() -= su.mesh().V()*su.field();
 998 }
 999
1000
1001 template<class Type>
1002 void Foam::fvMatrix<Type>::operator+=
1003 (
1004     const tmp<DimensionedField<Type, volMesh> >& tsu
1005 )
1006 {
1007     operator+=(tsu());
1008     tsu.clear();
1009 }
1010
1011
1012 template<class Type>
1013 void Foam::fvMatrix<Type>::operator+=
1014 (
1015     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
1016 )
1017 {
1018     operator+=(tsu());
1019     tsu.clear();
1020 }
1021
1022
1023 template<class Type>
1024 void Foam::fvMatrix<Type>::operator-=
1025 (
1026     const DimensionedField<Type, volMesh>& su
1027 )
1028 {
1029     checkMethod(*this, su, "-=");
1030     source() += su.mesh().V()*su.field();
1031 }
1032
1033
1034 template<class Type>
1035 void Foam::fvMatrix<Type>::operator-=
1036 (
1037     const tmp<DimensionedField<Type, volMesh> >& tsu
1038 )
1039 {
1040     operator-=(tsu());
1041     tsu.clear();
1042 }
1043
1044
1045 template<class Type>
1046 void Foam::fvMatrix<Type>::operator-=
1047 (
1048     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
1049 )
1050 {
1051     operator-=(tsu());
1052     tsu.clear();
1053 }
1054
1055
1056 template<class Type>
1057 void Foam::fvMatrix<Type>::operator+=
1058 (
1059     const dimensioned<Type>& su
1060 )
1061 {
1062     source() -= psi().mesh().V()*su;
1063 }
1064
1065
1066 template<class Type>
1067 void Foam::fvMatrix<Type>::operator-=
1068 (
1069     const dimensioned<Type>& su
1070 )
1071 {
1072     source() += psi().mesh().V()*su;
1073 }
1074
1075
1076 template<class Type>
1077 void Foam::fvMatrix<Type>::operator+=
1078 (
1079     const zeroField&
1080 )
1081 {}
1082
1083
1084 template<class Type>
1085 void Foam::fvMatrix<Type>::operator-=
1086 (
1087     const zeroField&
1088 )
1089 {}
1090
1091
1092 template<class Type>
1093 void Foam::fvMatrix<Type>::operator*=
1094 (
1095     const DimensionedField<scalar, volMesh>& dsf
1096 )
1097 {
1098     dimensions_ *= dsf.dimensions();
1099     lduMatrix::operator*=(dsf.field());
1100     source_ *= dsf.field();
1101
1102     forAll(boundaryCoeffs_, patchI)
1103     {
1104         scalarField psf
1105         (
1106             dsf.mesh().boundary()[patchI].patchInternalField(dsf.field())
1107         );
1108
1109         internalCoeffs_[patchI] *= psf;
1110         boundaryCoeffs_[patchI] *= psf;
1111     }
1112
1113     if (faceFluxCorrectionPtr_)
1114     {
1115         FatalErrorIn
1116         (
1117             "fvMatrix<Type>::operator*="
1118             "(const DimensionedField<scalar, volMesh>&)"
1119         )   << "cannot scale a matrix containing a faceFluxCorrection"
1120             << abort(FatalError);
1121     }
1122 }
1123
1124
1125 template<class Type>
1126 void Foam::fvMatrix<Type>::operator*=
1127 (
1128     const tmp<DimensionedField<scalar, volMesh> >& tdsf
1129 )
1130 {
1131     operator*=(tdsf());
1132     tdsf.clear();
1133 }
1134
1135
1136 template<class Type>
1137 void Foam::fvMatrix<Type>::operator*=
1138 (
1139     const tmp<volScalarField>& tvsf
1140 )
1141 {
1142     operator*=(tvsf());
1143     tvsf.clear();
1144 }
1145
1146
1147 template<class Type>
1148 void Foam::fvMatrix<Type>::operator*=
1149 (
1150     const dimensioned<scalar>& ds
1151 )
1152 {
1153     dimensions_ *= ds.dimensions();
1154     lduMatrix::operator*=(ds.value());
1155     source_ *= ds.value();
1156     internalCoeffs_ *= ds.value();
1157     boundaryCoeffs_ *= ds.value();
1158
1159     if (faceFluxCorrectionPtr_)
1160     {
1161         *faceFluxCorrectionPtr_ *= ds.value();
1162     }
1163 }
1164
1165
1166 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
1167
1168 template<class Type>
1169 void Foam::checkMethod
1170 (
1171     const fvMatrix<Type>& fvm1,
1172     const fvMatrix<Type>& fvm2,
1173     const char* op
1174 )
1175 {
1176     if (&fvm1.psi() != &fvm2.psi())
1177     {
1178         FatalErrorIn
1179         (
1180             "checkMethod(const fvMatrix<Type>&, const fvMatrix<Type>&)"
1181         )   << "incompatible fields for operation "
1182             << endl << "    "
1183             << "[" << fvm1.psi().name() << "] "
1184             << op
1185             << " [" << fvm2.psi().name() << "]"
1186             << abort(FatalError);
1187     }
1188
1189     if (dimensionSet::debug && fvm1.dimensions() != fvm2.dimensions())
1190     {
1191         FatalErrorIn
1192         (
1193             "checkMethod(const fvMatrix<Type>&, const fvMatrix<Type>&)"
1194         )   << "incompatible dimensions for operation "
1195             << endl << "    "
1196             << "[" << fvm1.psi().name() << fvm1.dimensions()/dimVolume << " ] "
1197             << op
1198             << " [" << fvm2.psi().name() << fvm2.dimensions()/dimVolume << " ]"
1199             << abort(FatalError);
1200     }
1201 }
1202
1203
1204 template<class Type>
1205 void Foam::checkMethod
1206 (
1207     const fvMatrix<Type>& fvm,
1208     const DimensionedField<Type, volMesh>& df,
1209     const char* op
1210 )
1211 {
1212     if (dimensionSet::debug && fvm.dimensions()/dimVolume != df.dimensions())
1213     {
1214         FatalErrorIn
1215         (
1216             "checkMethod(const fvMatrix<Type>&, const GeometricField<Type, "
1217             "fvPatchField, volMesh>&)"
1218         )   <<  "incompatible dimensions for operation "
1219             << endl << "    "
1220             << "[" << fvm.psi().name() << fvm.dimensions()/dimVolume << " ] "
1221             << op
1222             << " [" << df.name() << df.dimensions() << " ]"
1223             << abort(FatalError);
1224     }
1225 }
1226
1227
1228 template<class Type>
1229 void Foam::checkMethod
1230 (
1231     const fvMatrix<Type>& fvm,
1232     const dimensioned<Type>& dt,
1233     const char* op
1234 )
1235 {
1236     if (dimensionSet::debug && fvm.dimensions()/dimVolume != dt.dimensions())
1237     {
1238         FatalErrorIn
1239         (
1240             "checkMethod(const fvMatrix<Type>&, const dimensioned<Type>&)"
1241         )   << "incompatible dimensions for operation "
1242             << endl << "    "
1243             << "[" << fvm.psi().name() << fvm.dimensions()/dimVolume << " ] "
1244             << op
1245             << " [" << dt.name() << dt.dimensions() << " ]"
1246             << abort(FatalError);
1247     }
1248 }
1249
1250
1251 template<class Type>
1252 Foam::lduMatrix::solverPerformance Foam::solve
1253 (
1254     fvMatrix<Type>& fvm,
1255     const dictionary& solverControls
1256 )
1257 {
1258     return fvm.solve(solverControls);
1259 }
1260
1261 template<class Type>
1262 Foam::lduMatrix::solverPerformance Foam::solve
1263 (
1264     const tmp<fvMatrix<Type> >& tfvm,
1265     const dictionary& solverControls
1266 )
1267 {
1268     lduMatrix::solverPerformance solverPerf =
1269         const_cast<fvMatrix<Type>&>(tfvm()).solve(solverControls);
1270
1271     tfvm.clear();
1272
1273     return solverPerf;
1274 }
1275
1276
1277 template<class Type>
1278 Foam::lduMatrix::solverPerformance Foam::solve(fvMatrix<Type>& fvm)
1279 {
1280     return fvm.solve();
1281 }
1282
1283 template<class Type>
1284 Foam::lduMatrix::solverPerformance Foam::solve(const tmp<fvMatrix<Type> >& tfvm)
1285 {
1286     lduMatrix::solverPerformance solverPerf =
1287         const_cast<fvMatrix<Type>&>(tfvm()).solve();
1288
1289     tfvm.clear();
1290
1291     return solverPerf;
1292 }
1293
1294
1295 template<class Type>
1296 Foam::tmp<Foam::fvMatrix<Type> > Foam::correction
1297 (
1298     const fvMatrix<Type>& A
1299 )
1300 {
1301     tmp<Foam::fvMatrix<Type> > tAcorr = A - (A & A.psi());
1302
1303     if
1304     (
1305         (A.hasUpper() || A.hasLower())
1306      && A.psi().mesh().fluxRequired(A.psi().name())
1307     )
1308     {
1309         tAcorr().faceFluxCorrectionPtr() = (-A.flux()).ptr();
1310     }
1311
1312     return tAcorr;
1313 }
1314
1315
1316 template<class Type>
1317 Foam::tmp<Foam::fvMatrix<Type> > Foam::correction
1318 (
1319     const tmp<fvMatrix<Type> >& tA
1320 )
1321 {
1322     tmp<Foam::fvMatrix<Type> > tAcorr = tA - (tA() & tA().psi());
1323
1324     // Note the matrix coefficients are still that of matrix A
1325     const fvMatrix<Type>& A = tAcorr();
1326
1327     if
1328     (
1329         (A.hasUpper() || A.hasLower())
1330      && A.psi().mesh().fluxRequired(A.psi().name())
1331     )
1332     {
1333         tAcorr().faceFluxCorrectionPtr() = (-A.flux()).ptr();
1334     }
1335
1336     return tAcorr;
1337 }
1338
1339
1340 // * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
1341
1342 template<class Type>
1343 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1344 (
1345     const fvMatrix<Type>& A,
1346     const fvMatrix<Type>& B
1347 )
1348 {
1349     checkMethod(A, B, "==");
1350     return (A - B);
1351 }
1352
1353 template<class Type>
1354 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1355 (
1356     const tmp<fvMatrix<Type> >& tA,
1357     const fvMatrix<Type>& B
1358 )
1359 {
1360     checkMethod(tA(), B, "==");
1361     return (tA - B);
1362 }
1363
1364 template<class Type>
1365 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1366 (
1367     const fvMatrix<Type>& A,
1368     const tmp<fvMatrix<Type> >& tB
1369 )
1370 {
1371     checkMethod(A, tB(), "==");
1372     return (A - tB);
1373 }
1374
1375 template<class Type>
1376 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1377 (
1378     const tmp<fvMatrix<Type> >& tA,
1379     const tmp<fvMatrix<Type> >& tB
1380 )
1381 {
1382     checkMethod(tA(), tB(), "==");
1383     return (tA - tB);
1384 }
1385
1386 template<class Type>
1387 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1388 (
1389     const fvMatrix<Type>& A,
1390     const DimensionedField<Type, volMesh>& su
1391 )
1392 {
1393     checkMethod(A, su, "==");
1394     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1395     tC().source() += su.mesh().V()*su.field();
1396     return tC;
1397 }
1398
1399 template<class Type>
1400 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1401 (
1402     const fvMatrix<Type>& A,
1403     const tmp<DimensionedField<Type, volMesh> >& tsu
1404 )
1405 {
1406     checkMethod(A, tsu(), "==");
1407     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1408     tC().source() += tsu().mesh().V()*tsu().field();
1409     tsu.clear();
1410     return tC;
1411 }
1412
1413 template<class Type>
1414 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1415 (
1416     const fvMatrix<Type>& A,
1417     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
1418 )
1419 {
1420     checkMethod(A, tsu(), "==");
1421     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1422     tC().source() += tsu().mesh().V()*tsu().internalField();
1423     tsu.clear();
1424     return tC;
1425 }
1426
1427 template<class Type>
1428 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1429 (
1430     const tmp<fvMatrix<Type> >& tA,
1431     const DimensionedField<Type, volMesh>& su
1432 )
1433 {
1434     checkMethod(tA(), su, "==");
1435     tmp<fvMatrix<Type> > tC(tA.ptr());
1436     tC().source() += su.mesh().V()*su.field();
1437     return tC;
1438 }
1439
1440 template<class Type>
1441 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1442 (
1443     const tmp<fvMatrix<Type> >& tA,
1444     const tmp<DimensionedField<Type, volMesh> >& tsu
1445 )
1446 {
1447     checkMethod(tA(), tsu(), "==");
1448     tmp<fvMatrix<Type> > tC(tA.ptr());
1449     tC().source() += tsu().mesh().V()*tsu().field();
1450     tsu.clear();
1451     return tC;
1452 }
1453
1454 template<class Type>
1455 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1456 (
1457     const tmp<fvMatrix<Type> >& tA,
1458     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
1459 )
1460 {
1461     checkMethod(tA(), tsu(), "==");
1462     tmp<fvMatrix<Type> > tC(tA.ptr());
1463     tC().source() += tsu().mesh().V()*tsu().internalField();
1464     tsu.clear();
1465     return tC;
1466 }
1467
1468 template<class Type>
1469 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1470 (
1471     const fvMatrix<Type>& A,
1472     const dimensioned<Type>& su
1473 )
1474 {
1475     checkMethod(A, su, "==");
1476     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1477     tC().source() += A.psi().mesh().V()*su.value();
1478     return tC;
1479 }
1480
1481 template<class Type>
1482 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1483 (
1484     const tmp<fvMatrix<Type> >& tA,
1485     const dimensioned<Type>& su
1486 )
1487 {
1488     checkMethod(tA(), su, "==");
1489     tmp<fvMatrix<Type> > tC(tA.ptr());
1490     tC().source() += tC().psi().mesh().V()*su.value();
1491     return tC;
1492 }
1493
1494 template<class Type>
1495 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1496 (
1497     const fvMatrix<Type>& A,
1498     const zeroField&
1499 )
1500 {
1501     return A;
1502 }
1503
1504
1505 template<class Type>
1506 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1507 (
1508     const tmp<fvMatrix<Type> >& tA,
1509     const zeroField&
1510 )
1511 {
1512     return tA;
1513 }
1514
1515
1516 template<class Type>
1517 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1518 (
1519     const fvMatrix<Type>& A
1520 )
1521 {
1522     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1523     tC().negate();
1524     return tC;
1525 }
1526
1527 template<class Type>
1528 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1529 (
1530     const tmp<fvMatrix<Type> >& tA
1531 )
1532 {
1533     tmp<fvMatrix<Type> > tC(tA.ptr());
1534     tC().negate();
1535     return tC;
1536 }
1537
1538
1539 template<class Type>
1540 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1541 (
1542     const fvMatrix<Type>& A,
1543     const fvMatrix<Type>& B
1544 )
1545 {
1546     checkMethod(A, B, "+");
1547     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1548     tC() += B;
1549     return tC;
1550 }
1551
1552 template<class Type>
1553 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1554 (
1555     const tmp<fvMatrix<Type> >& tA,
1556     const fvMatrix<Type>& B
1557 )
1558 {
1559     checkMethod(tA(), B, "+");
1560     tmp<fvMatrix<Type> > tC(tA.ptr());
1561     tC() += B;
1562     return tC;
1563 }
1564
1565 template<class Type>
1566 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1567 (
1568     const fvMatrix<Type>& A,
1569     const tmp<fvMatrix<Type> >& tB
1570 )
1571 {
1572     checkMethod(A, tB(), "+");
1573     tmp<fvMatrix<Type> > tC(tB.ptr());
1574     tC() += A;
1575     return tC;
1576 }
1577
1578 template<class Type>
1579 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1580 (
1581     const tmp<fvMatrix<Type> >& tA,
1582     const tmp<fvMatrix<Type> >& tB
1583 )
1584 {
1585     checkMethod(tA(), tB(), "+");
1586     tmp<fvMatrix<Type> > tC(tA.ptr());
1587     tC() += tB();
1588     tB.clear();
1589     return tC;
1590 }
1591
1592 template<class Type>
1593 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1594 (
1595     const fvMatrix<Type>& A,
1596     const DimensionedField<Type, volMesh>& su
1597 )
1598 {
1599     checkMethod(A, su, "+");
1600     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1601     tC().source() -= su.mesh().V()*su.field();
1602     return tC;
1603 }
1604
1605 template<class Type>
1606 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1607 (
1608     const fvMatrix<Type>& A,
1609     const tmp<DimensionedField<Type, volMesh> >& tsu
1610 )
1611 {
1612     checkMethod(A, tsu(), "+");
1613     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1614     tC().source() -= tsu().mesh().V()*tsu().field();
1615     tsu.clear();
1616     return tC;
1617 }
1618
1619 template<class Type>
1620 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1621 (
1622     const fvMatrix<Type>& A,
1623     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
1624 )
1625 {
1626     checkMethod(A, tsu(), "+");
1627     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1628     tC().source() -= tsu().mesh().V()*tsu().internalField();
1629     tsu.clear();
1630     return tC;
1631 }
1632
1633 template<class Type>
1634 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1635 (
1636     const tmp<fvMatrix<Type> >& tA,
1637     const DimensionedField<Type, volMesh>& su
1638 )
1639 {
1640     checkMethod(tA(), su, "+");
1641     tmp<fvMatrix<Type> > tC(tA.ptr());
1642     tC().source() -= su.mesh().V()*su.field();
1643     return tC;
1644 }
1645
1646 template<class Type>
1647 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1648 (
1649     const tmp<fvMatrix<Type> >& tA,
1650     const tmp<DimensionedField<Type, volMesh> >& tsu
1651 )
1652 {
1653     checkMethod(tA(), tsu(), "+");
1654     tmp<fvMatrix<Type> > tC(tA.ptr());
1655     tC().source() -= tsu().mesh().V()*tsu().field();
1656     tsu.clear();
1657     return tC;
1658 }
1659
1660 template<class Type>
1661 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1662 (
1663     const tmp<fvMatrix<Type> >& tA,
1664     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
1665 )
1666 {
1667     checkMethod(tA(), tsu(), "+");
1668     tmp<fvMatrix<Type> > tC(tA.ptr());
1669     tC().source() -= tsu().mesh().V()*tsu().internalField();
1670     tsu.clear();
1671     return tC;
1672 }
1673
1674 template<class Type>
1675 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1676 (
1677     const DimensionedField<Type, volMesh>& su,
1678     const fvMatrix<Type>& A
1679 )
1680 {
1681     checkMethod(A, su, "+");
1682     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1683     tC().source() -= su.mesh().V()*su.field();
1684     return tC;
1685 }
1686
1687 template<class Type>
1688 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1689 (
1690     const tmp<DimensionedField<Type, volMesh> >& tsu,
1691     const fvMatrix<Type>& A
1692 )
1693 {
1694     checkMethod(A, tsu(), "+");
1695     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1696     tC().source() -= tsu().mesh().V()*tsu().field();
1697     tsu.clear();
1698     return tC;
1699 }
1700
1701 template<class Type>
1702 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1703 (
1704     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu,
1705     const fvMatrix<Type>& A
1706 )
1707 {
1708     checkMethod(A, tsu(), "+");
1709     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1710     tC().source() -= tsu().mesh().V()*tsu().internalField();
1711     tsu.clear();
1712     return tC;
1713 }
1714
1715 template<class Type>
1716 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1717 (
1718     const DimensionedField<Type, volMesh>& su,
1719     const tmp<fvMatrix<Type> >& tA
1720 )
1721 {
1722     checkMethod(tA(), su, "+");
1723     tmp<fvMatrix<Type> > tC(tA.ptr());
1724     tC().source() -= su.mesh().V()*su.field();
1725     return tC;
1726 }
1727
1728 template<class Type>
1729 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1730 (
1731     const tmp<DimensionedField<Type, volMesh> >& tsu,
1732     const tmp<fvMatrix<Type> >& tA
1733 )
1734 {
1735     checkMethod(tA(), tsu(), "+");
1736     tmp<fvMatrix<Type> > tC(tA.ptr());
1737     tC().source() -= tsu().mesh().V()*tsu().field();
1738     tsu.clear();
1739     return tC;
1740 }
1741
1742 template<class Type>
1743 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1744 (
1745     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu,
1746     const tmp<fvMatrix<Type> >& tA
1747 )
1748 {
1749     checkMethod(tA(), tsu(), "+");
1750     tmp<fvMatrix<Type> > tC(tA.ptr());
1751     tC().source() -= tsu().mesh().V()*tsu().internalField();
1752     tsu.clear();
1753     return tC;
1754 }
1755
1756
1757 template<class Type>
1758 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1759 (
1760     const fvMatrix<Type>& A,
1761     const fvMatrix<Type>& B
1762 )
1763 {
1764     checkMethod(A, B, "-");
1765     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1766     tC() -= B;
1767     return tC;
1768 }
1769
1770 template<class Type>
1771 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1772 (
1773     const tmp<fvMatrix<Type> >& tA,
1774     const fvMatrix<Type>& B
1775 )
1776 {
1777     checkMethod(tA(), B, "-");
1778     tmp<fvMatrix<Type> > tC(tA.ptr());
1779     tC() -= B;
1780     return tC;
1781 }
1782
1783 template<class Type>
1784 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1785 (
1786     const fvMatrix<Type>& A,
1787     const tmp<fvMatrix<Type> >& tB
1788 )
1789 {
1790     checkMethod(A, tB(), "-");
1791     tmp<fvMatrix<Type> > tC(tB.ptr());
1792     tC() -= A;
1793     tC().negate();
1794     return tC;
1795 }
1796
1797 template<class Type>
1798 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1799 (
1800     const tmp<fvMatrix<Type> >& tA,
1801     const tmp<fvMatrix<Type> >& tB
1802 )
1803 {
1804     checkMethod(tA(), tB(), "-");
1805     tmp<fvMatrix<Type> > tC(tA.ptr());
1806     tC() -= tB();
1807     tB.clear();
1808     return tC;
1809 }
1810
1811 template<class Type>
1812 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1813 (
1814     const fvMatrix<Type>& A,
1815     const DimensionedField<Type, volMesh>& su
1816 )
1817 {
1818     checkMethod(A, su, "-");
1819     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1820     tC().source() += su.mesh().V()*su.field();
1821     return tC;
1822 }
1823
1824 template<class Type>
1825 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1826 (
1827     const fvMatrix<Type>& A,
1828     const tmp<DimensionedField<Type, volMesh> >& tsu
1829 )
1830 {
1831     checkMethod(A, tsu(), "-");
1832     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1833     tC().source() += tsu().mesh().V()*tsu().field();
1834     tsu.clear();
1835     return tC;
1836 }
1837
1838 template<class Type>
1839 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1840 (
1841     const fvMatrix<Type>& A,
1842     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
1843 )
1844 {
1845     checkMethod(A, tsu(), "-");
1846     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1847     tC().source() += tsu().mesh().V()*tsu().internalField();
1848     tsu.clear();
1849     return tC;
1850 }
1851
1852 template<class Type>
1853 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1854 (
1855     const tmp<fvMatrix<Type> >& tA,
1856     const DimensionedField<Type, volMesh>& su
1857 )
1858 {
1859     checkMethod(tA(), su, "-");
1860     tmp<fvMatrix<Type> > tC(tA.ptr());
1861     tC().source() += su.mesh().V()*su.field();
1862     return tC;
1863 }
1864
1865 template<class Type>
1866 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1867 (
1868     const tmp<fvMatrix<Type> >& tA,
1869     const tmp<DimensionedField<Type, volMesh> >& tsu
1870 )
1871 {
1872     checkMethod(tA(), tsu(), "-");
1873     tmp<fvMatrix<Type> > tC(tA.ptr());
1874     tC().source() += tsu().mesh().V()*tsu().field();
1875     tsu.clear();
1876     return tC;
1877 }
1878
1879 template<class Type>
1880 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1881 (
1882     const tmp<fvMatrix<Type> >& tA,
1883     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
1884 )
1885 {
1886     checkMethod(tA(), tsu(), "-");
1887     tmp<fvMatrix<Type> > tC(tA.ptr());
1888     tC().source() += tsu().mesh().V()*tsu().internalField();
1889     tsu.clear();
1890     return tC;
1891 }
1892
1893 template<class Type>
1894 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1895 (
1896     const DimensionedField<Type, volMesh>& su,
1897     const fvMatrix<Type>& A
1898 )
1899 {
1900     checkMethod(A, su, "-");
1901     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1902     tC().negate();
1903     tC().source() -= su.mesh().V()*su.field();
1904     return tC;
1905 }
1906
1907 template<class Type>
1908 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1909 (
1910     const tmp<DimensionedField<Type, volMesh> >& tsu,
1911     const fvMatrix<Type>& A
1912 )
1913 {
1914     checkMethod(A, tsu(), "-");
1915     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1916     tC().negate();
1917     tC().source() -= tsu().mesh().V()*tsu().field();
1918     tsu.clear();
1919     return tC;
1920 }
1921
1922 template<class Type>
1923 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1924 (
1925     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu,
1926     const fvMatrix<Type>& A
1927 )
1928 {
1929     checkMethod(A, tsu(), "-");
1930     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1931     tC().negate();
1932     tC().source() -= tsu().mesh().V()*tsu().internalField();
1933     tsu.clear();
1934     return tC;
1935 }
1936
1937 template<class Type>
1938 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1939 (
1940     const DimensionedField<Type, volMesh>& su,
1941     const tmp<fvMatrix<Type> >& tA
1942 )
1943 {
1944     checkMethod(tA(), su, "-");
1945     tmp<fvMatrix<Type> > tC(tA.ptr());
1946     tC().negate();
1947     tC().source() -= su.mesh().V()*su.field();
1948     return tC;
1949 }
1950
1951 template<class Type>
1952 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1953 (
1954     const tmp<DimensionedField<Type, volMesh> >& tsu,
1955     const tmp<fvMatrix<Type> >& tA
1956 )
1957 {
1958     checkMethod(tA(), tsu(), "-");
1959     tmp<fvMatrix<Type> > tC(tA.ptr());
1960     tC().negate();
1961     tC().source() -= tsu().mesh().V()*tsu().field();
1962     tsu.clear();
1963     return tC;
1964 }
1965
1966 template<class Type>
1967 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1968 (
1969     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu,
1970     const tmp<fvMatrix<Type> >& tA
1971 )
1972 {
1973     checkMethod(tA(), tsu(), "-");
1974     tmp<fvMatrix<Type> > tC(tA.ptr());
1975     tC().negate();
1976     tC().source() -= tsu().mesh().V()*tsu().internalField();
1977     tsu.clear();
1978     return tC;
1979 }
1980
1981 template<class Type>
1982 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1983 (
1984     const fvMatrix<Type>& A,
1985     const dimensioned<Type>& su
1986 )
1987 {
1988     checkMethod(A, su, "+");
1989     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1990     tC().source() -= su.value()*A.psi().mesh().V();
1991     return tC;
1992 }
1993
1994 template<class Type>
1995 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1996 (
1997     const tmp<fvMatrix<Type> >& tA,
1998     const dimensioned<Type>& su
1999 )
2000 {
2001     checkMethod(tA(), su, "+");
2002     tmp<fvMatrix<Type> > tC(tA.ptr());
2003     tC().source() -= su.value()*tC().psi().mesh().V();
2004     return tC;
2005 }
2006
2007 template<class Type>
2008 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
2009 (
2010     const dimensioned<Type>& su,
2011     const fvMatrix<Type>& A
2012 )
2013 {
2014     checkMethod(A, su, "+");
2015     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
2016     tC().source() -= su.value()*A.psi().mesh().V();
2017     return tC;
2018 }
2019
2020 template<class Type>
2021 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
2022 (
2023     const dimensioned<Type>& su,
2024     const tmp<fvMatrix<Type> >& tA
2025 )
2026 {
2027     checkMethod(tA(), su, "+");
2028     tmp<fvMatrix<Type> > tC(tA.ptr());
2029     tC().source() -= su.value()*tC().psi().mesh().V();
2030     return tC;
2031 }
2032
2033 template<class Type>
2034 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
2035 (
2036     const fvMatrix<Type>& A,
2037     const dimensioned<Type>& su
2038 )
2039 {
2040     checkMethod(A, su, "-");
2041     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
2042     tC().source() += su.value()*tC().psi().mesh().V();
2043     return tC;
2044 }
2045
2046 template<class Type>
2047 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
2048 (
2049     const tmp<fvMatrix<Type> >& tA,
2050     const dimensioned<Type>& su
2051 )
2052 {
2053     checkMethod(tA(), su, "-");
2054     tmp<fvMatrix<Type> > tC(tA.ptr());
2055     tC().source() += su.value()*tC().psi().mesh().V();
2056     return tC;
2057 }
2058
2059 template<class Type>
2060 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
2061 (
2062     const dimensioned<Type>& su,
2063     const fvMatrix<Type>& A
2064 )
2065 {
2066     checkMethod(A, su, "-");
2067     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
2068     tC().negate();
2069     tC().source() -= su.value()*A.psi().mesh().V();
2070     return tC;
2071 }
2072
2073 template<class Type>
2074 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
2075 (
2076     const dimensioned<Type>& su,
2077     const tmp<fvMatrix<Type> >& tA
2078 )
2079 {
2080     checkMethod(tA(), su, "-");
2081     tmp<fvMatrix<Type> > tC(tA.ptr());
2082     tC().negate();
2083     tC().source() -= su.value()*tC().psi().mesh().V();
2084     return tC;
2085 }
2086
2087
2088 template<class Type>
2089 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator*
2090 (
2091     const DimensionedField<scalar, volMesh>& dsf,
2092     const fvMatrix<Type>& A
2093 )
2094 {
2095     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
2096     tC() *= dsf;
2097     return tC;
2098 }
2099
2100 template<class Type>
2101 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator*
2102 (
2103     const tmp< DimensionedField<scalar, volMesh> >& tdsf,
2104     const fvMatrix<Type>& A
2105 )
2106 {
2107     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
2108     tC() *= tdsf;
2109     return tC;
2110 }
2111
2112 template<class Type>
2113 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator*
2114 (
2115     const tmp<volScalarField>& tvsf,
2116     const fvMatrix<Type>& A
2117 )
2118 {
2119     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
2120     tC() *= tvsf;
2121     return tC;
2122 }
2123
2124 template<class Type>
2125 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator*
2126 (
2127     const DimensionedField<scalar, volMesh>& dsf,
2128     const tmp<fvMatrix<Type> >& tA
2129 )
2130 {
2131     tmp<fvMatrix<Type> > tC(tA.ptr());
2132     tC() *= dsf;
2133     return tC;
2134 }
2135
2136 template<class Type>
2137 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator*
2138 (
2139     const tmp<DimensionedField<scalar, volMesh> >& tdsf,
2140     const tmp<fvMatrix<Type> >& tA
2141 )
2142 {
2143     tmp<fvMatrix<Type> > tC(tA.ptr());
2144     tC() *= tdsf;
2145     return tC;
2146 }
2147
2148 template<class Type>
2149 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator*
2150 (
2151     const tmp<volScalarField>& tvsf,
2152     const tmp<fvMatrix<Type> >& tA
2153 )
2154 {
2155     tmp<fvMatrix<Type> > tC(tA.ptr());
2156     tC() *= tvsf;
2157     return tC;
2158 }
2159
2160 template<class Type>
2161 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator*
2162 (
2163     const dimensioned<scalar>& ds,
2164     const fvMatrix<Type>& A
2165 )
2166 {
2167     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
2168     tC() *= ds;
2169     return tC;
2170 }
2171
2172 template<class Type>
2173 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator*
2174 (
2175     const dimensioned<scalar>& ds,
2176     const tmp<fvMatrix<Type> >& tA
2177 )
2178 {
2179     tmp<fvMatrix<Type> > tC(tA.ptr());
2180     tC() *= ds;
2181     return tC;
2182 }
2183
2184
2185 template<class Type>
2186 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
2187 Foam::operator&
2188 (
2189     const fvMatrix<Type>& M,
2190     const DimensionedField<Type, volMesh>& psi
2191 )
2192 {
2193     tmp<GeometricField<Type, fvPatchField, volMesh> > tMphi
2194     (
2195         new GeometricField<Type, fvPatchField, volMesh>
2196         (
2197             IOobject
2198             (
2199                 "M&" + psi.name(),
2200                 psi.instance(),
2201                 psi.mesh(),
2202                 IOobject::NO_READ,
2203                 IOobject::NO_WRITE
2204             ),
2205             psi.mesh(),
2206             M.dimensions()/dimVol,
2207             zeroGradientFvPatchScalarField::typeName
2208         )
2209     );
2210     GeometricField<Type, fvPatchField, volMesh>& Mphi = tMphi();
2211
2212     // Loop over field components
2213     for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
2214     {
2215         scalarField psiCmpt = psi.field().component(cmpt);
2216         scalarField boundaryDiagCmpt(M.diag());
2217         M.addBoundaryDiag(boundaryDiagCmpt, cmpt);
2218         Mphi.internalField().replace(cmpt, -boundaryDiagCmpt*psiCmpt);
2219     }
2220
2221     Mphi.internalField() += M.lduMatrix::H(psi.field()) + M.source();
2222     M.addBoundarySource(Mphi.internalField());
2223
2224     Mphi.internalField() /= -psi.mesh().V();
2225     Mphi.correctBoundaryConditions();
2226
2227     return tMphi;
2228 }
2229
2230 template<class Type>
2231 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
2232 Foam::operator&
2233 (
2234     const fvMatrix<Type>& M,
2235     const tmp<DimensionedField<Type, volMesh> >& tpsi
2236 )
2237 {
2238     tmp<GeometricField<Type, fvPatchField, volMesh> > tMpsi = M & tpsi();
2239     tpsi.clear();
2240     return tMpsi;
2241 }
2242
2243 template<class Type>
2244 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
2245 Foam::operator&
2246 (
2247     const fvMatrix<Type>& M,
2248     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tpsi
2249 )
2250 {
2251     tmp<GeometricField<Type, fvPatchField, volMesh> > tMpsi = M & tpsi();
2252     tpsi.clear();
2253     return tMpsi;
2254 }
2255
2256 template<class Type>
2257 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
2258 Foam::operator&
2259 (
2260     const tmp<fvMatrix<Type> >& tM,
2261     const DimensionedField<Type, volMesh>& psi
2262 )
2263 {
2264     tmp<GeometricField<Type, fvPatchField, volMesh> > tMpsi = tM() & psi;
2265     tM.clear();
2266     return tMpsi;
2267 }
2268
2269 template<class Type>
2270 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
2271 Foam::operator&
2272 (
2273     const tmp<fvMatrix<Type> >& tM,
2274     const tmp<DimensionedField<Type, volMesh> >& tpsi
2275 )
2276 {
2277     tmp<GeometricField<Type, fvPatchField, volMesh> > tMpsi = tM() & tpsi();
2278     tM.clear();
2279     tpsi.clear();
2280     return tMpsi;
2281 }
2282
2283 template<class Type>
2284 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
2285 Foam::operator&
2286 (
2287     const tmp<fvMatrix<Type> >& tM,
2288     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tpsi
2289 )
2290 {
2291     tmp<GeometricField<Type, fvPatchField, volMesh> > tMpsi = tM() & tpsi();
2292     tM.clear();
2293     tpsi.clear();
2294     return tMpsi;
2295 }
2296
2297
2298 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
2299
2300 template<class Type>
2301 Foam::Ostream& Foam::operator<<(Ostream& os, const fvMatrix<Type>& fvm)
2302 {
2303     os  << static_cast<const lduMatrix&>(fvm) << nl
2304         << fvm.dimensions_ << nl
2305         << fvm.source_ << nl
2306         << fvm.internalCoeffs_ << nl
2307         << fvm.boundaryCoeffs_ << endl;
2308
2309     os.check("Ostream& operator<<(Ostream&, fvMatrix<Type>&");
2310
2311     return os;
2312 }
2313
2314
2315 // * * * * * * * * * * * * * * * * Solvers * * * * * * * * * * * * * * * * * //
2316
2317 #include "fvMatrixSolve.C"
2318
2319 // ************************************************************************* //