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