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>"
 340                "(GeometricField<Type, fvPatchField, volMesh>&, 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 void Foam::fvMatrix<Type>::boundaryManipulate
 593 (
 594     typename GeometricField<Type, fvPatchField, volMesh>::
 595         GeometricBoundaryField& bFields
 596 )
 597 {
 598     forAll(bFields, patchI)
 599     {
 600         bFields[patchI].manipulateMatrix(*this);
 601     }
 602 }
 603
 604
 605 template<class Type>
 606 Foam::tmp<Foam::scalarField> Foam::fvMatrix<Type>::D() const
 607 {
 608     tmp<scalarField> tdiag(new scalarField(diag()));
 609     addCmptAvBoundaryDiag(tdiag());
 610     return tdiag;
 611 }
 612
 613
 614 template<class Type>
 615 Foam::tmp<Foam::Field<Type> > Foam::fvMatrix<Type>::DD() const
 616 {
 617     tmp<Field<Type> > tdiag(pTraits<Type>::one*diag());
 618
 619     forAll(psi_.boundaryField(), patchI)
 620     {
 621         const fvPatchField<Type>& ptf = psi_.boundaryField()[patchI];
 622
 623         if (!ptf.coupled() && ptf.size())
 624         {
 625             addToInternalField
 626             (
 627                 lduAddr().patchAddr(patchI),
 628                 internalCoeffs_[patchI],
 629                 tdiag()
 630             );
 631         }
 632     }
 633
 634     return tdiag;
 635 }
 636
 637
 638 template<class Type>
 639 Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Type>::A() const
 640 {
 641     tmp<volScalarField> tAphi
 642     (
 643         new volScalarField
 644         (
 645             IOobject
 646             (
 647                 "A("+psi_.name()+')',
 648                 psi_.instance(),
 649                 psi_.mesh(),
 650                 IOobject::NO_READ,
 651                 IOobject::NO_WRITE
 652             ),
 653             psi_.mesh(),
 654             dimensions_/psi_.dimensions()/dimVol,
 655             zeroGradientFvPatchScalarField::typeName
 656         )
 657     );
 658
 659     tAphi().internalField() = D()/psi_.mesh().V();
 660     tAphi().correctBoundaryConditions();
 661
 662     return tAphi;
 663 }
 664
 665
 666 template<class Type>
 667 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
 668 Foam::fvMatrix<Type>::H() const
 669 {
 670     tmp<GeometricField<Type, fvPatchField, volMesh> > tHphi
 671     (
 672         new GeometricField<Type, fvPatchField, volMesh>
 673         (
 674             IOobject
 675             (
 676                 "H("+psi_.name()+')',
 677                 psi_.instance(),
 678                 psi_.mesh(),
 679                 IOobject::NO_READ,
 680                 IOobject::NO_WRITE
 681             ),
 682             psi_.mesh(),
 683             dimensions_/dimVol,
 684             zeroGradientFvPatchScalarField::typeName
 685         )
 686     );
 687     GeometricField<Type, fvPatchField, volMesh>& Hphi = tHphi();
 688
 689     // Loop over field components
 690     for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
 691     {
 692         scalarField psiCmpt = psi_.internalField().component(cmpt);
 693
 694         scalarField boundaryDiagCmpt(psi_.size(), 0.0);
 695         addBoundaryDiag(boundaryDiagCmpt, cmpt);
 696         boundaryDiagCmpt.negate();
 697         addCmptAvBoundaryDiag(boundaryDiagCmpt);
 698
 699         Hphi.internalField().replace(cmpt, boundaryDiagCmpt*psiCmpt);
 700     }
 701
 702     Hphi.internalField() += lduMatrix::H(psi_.internalField()) + source_;
 703     addBoundarySource(Hphi.internalField());
 704
 705     Hphi.internalField() /= psi_.mesh().V();
 706     Hphi.correctBoundaryConditions();
 707
 708     typename Type::labelType validComponents
 709     (
 710         pow
 711         (
 712             psi_.mesh().solutionD(),
 713             pTraits<typename powProduct<Vector<label>, Type::rank>::type>::zero
 714         )
 715     );
 716
 717     for(direction cmpt=0; cmpt<Type::nComponents; cmpt++)
 718     {
 719         if (validComponents[cmpt] == -1)
 720         {
 721             Hphi.replace
 722             (
 723                 cmpt,
 724                 dimensionedScalar("0", Hphi.dimensions(), 0.0)
 725             );
 726         }
 727     }
 728
 729     return tHphi;
 730 }
 731
 732
 733 template<class Type>
 734 Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Type>::H1() const
 735 {
 736     tmp<volScalarField> tH1
 737     (
 738         new volScalarField
 739         (
 740             IOobject
 741             (
 742                 "H(1)",
 743                 psi_.instance(),
 744                 psi_.mesh(),
 745                 IOobject::NO_READ,
 746                 IOobject::NO_WRITE
 747             ),
 748             psi_.mesh(),
 749             dimensions_/(dimVol*psi_.dimensions()),
 750             zeroGradientFvPatchScalarField::typeName
 751         )
 752     );
 753     volScalarField& H1_ = tH1();
 754
 755     // Loop over field components
 756     /*
 757     for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
 758     {
 759         scalarField psiCmpt = psi_.internalField().component(cmpt);
 760
 761         scalarField boundaryDiagCmpt(psi_.size(), 0.0);
 762         addBoundaryDiag(boundaryDiagCmpt, cmpt);
 763         boundaryDiagCmpt.negate();
 764         addCmptAvBoundaryDiag(boundaryDiagCmpt);
 765
 766         H1_.internalField().replace(cmpt, boundaryDiagCmpt);
 767     }
 768
 769     H1_.internalField() += lduMatrix::H1();
 770     */
 771
 772     H1_.internalField() = lduMatrix::H1();
 773
 774     H1_.internalField() /= psi_.mesh().V();
 775     H1_.correctBoundaryConditions();
 776
 777     return tH1;
 778 }
 779
 780
 781 template<class Type>
 782 Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> >
 783 Foam::fvMatrix<Type>::
 784 flux() const
 785 {
 786     if (!psi_.mesh().fluxRequired(psi_.name()))
 787     {
 788         FatalErrorIn("fvMatrix<Type>::flux()")
 789             << "flux requested but " << psi_.name()
 790             << " not specified in the fluxRequired sub-dictionary"
 791                " of fvSchemes."
 792             << abort(FatalError);
 793     }
 794
 795     // construct GeometricField<Type, fvsPatchField, surfaceMesh>
 796     tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tfieldFlux
 797     (
 798         new GeometricField<Type, fvsPatchField, surfaceMesh>
 799         (
 800             IOobject
 801             (
 802                 "flux("+psi_.name()+')',
 803                 psi_.instance(),
 804                 psi_.mesh(),
 805                 IOobject::NO_READ,
 806                 IOobject::NO_WRITE
 807             ),
 808             psi_.mesh(),
 809             dimensions()
 810         )
 811     );
 812     GeometricField<Type, fvsPatchField, surfaceMesh>& fieldFlux = tfieldFlux();
 813
 814     for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
 815     {
 816         fieldFlux.internalField().replace
 817         (
 818             cmpt,
 819             lduMatrix::faceH(psi_.internalField().component(cmpt))
 820         );
 821     }
 822
 823     FieldField<Field, Type> InternalContrib = internalCoeffs_;
 824
 825     forAll(InternalContrib, patchI)
 826     {
 827         InternalContrib[patchI] =
 828             cmptMultiply
 829             (
 830                 InternalContrib[patchI],
 831                 psi_.boundaryField()[patchI].patchInternalField()
 832             );
 833     }
 834
 835     FieldField<Field, Type> NeighbourContrib = boundaryCoeffs_;
 836
 837     forAll(NeighbourContrib, patchI)
 838     {
 839         if (psi_.boundaryField()[patchI].coupled())
 840         {
 841             NeighbourContrib[patchI] =
 842                 cmptMultiply
 843                 (
 844                     NeighbourContrib[patchI],
 845                     psi_.boundaryField()[patchI].patchNeighbourField()
 846                 );
 847         }
 848     }
 849
 850     forAll(fieldFlux.boundaryField(), patchI)
 851     {
 852         fieldFlux.boundaryField()[patchI] =
 853             InternalContrib[patchI] - NeighbourContrib[patchI];
 854     }
 855
 856     if (faceFluxCorrectionPtr_)
 857     {
 858         fieldFlux += *faceFluxCorrectionPtr_;
 859     }
 860
 861     return tfieldFlux;
 862 }
 863
 864
 865 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 866
 867 template<class Type>
 868 void Foam::fvMatrix<Type>::operator=(const fvMatrix<Type>& fvmv)
 869 {
 870     if (this == &fvmv)
 871     {
 872         FatalErrorIn("fvMatrix<Type>::operator=(const fvMatrix<Type>&)")
 873             << "attempted assignment to self"
 874             << abort(FatalError);
 875     }
 876
 877     if (&psi_ != &(fvmv.psi_))
 878     {
 879         FatalErrorIn("fvMatrix<Type>::operator=(const fvMatrix<Type>&)")
 880             << "different fields"
 881             << abort(FatalError);
 882     }
 883
 884     lduMatrix::operator=(fvmv);
 885     source_ = fvmv.source_;
 886     internalCoeffs_ = fvmv.internalCoeffs_;
 887     boundaryCoeffs_ = fvmv.boundaryCoeffs_;
 888
 889     if (faceFluxCorrectionPtr_ && fvmv.faceFluxCorrectionPtr_)
 890     {
 891         *faceFluxCorrectionPtr_ = *fvmv.faceFluxCorrectionPtr_;
 892     }
 893     else if (fvmv.faceFluxCorrectionPtr_)
 894     {
 895         faceFluxCorrectionPtr_ =
 896             new GeometricField<Type, fvsPatchField, surfaceMesh>
 897         (*fvmv.faceFluxCorrectionPtr_);
 898     }
 899 }
 900
 901
 902 template<class Type>
 903 void Foam::fvMatrix<Type>::operator=(const tmp<fvMatrix<Type> >& tfvmv)
 904 {
 905     operator=(tfvmv());
 906     tfvmv.clear();
 907 }
 908
 909
 910 template<class Type>
 911 void Foam::fvMatrix<Type>::negate()
 912 {
 913     lduMatrix::negate();
 914     source_.negate();
 915     internalCoeffs_.negate();
 916     boundaryCoeffs_.negate();
 917
 918     if (faceFluxCorrectionPtr_)
 919     {
 920         faceFluxCorrectionPtr_->negate();
 921     }
 922 }
 923
 924
 925 template<class Type>
 926 void Foam::fvMatrix<Type>::operator+=(const fvMatrix<Type>& fvmv)
 927 {
 928     checkMethod(*this, fvmv, "+=");
 929
 930     dimensions_ += fvmv.dimensions_;
 931     lduMatrix::operator+=(fvmv);
 932     source_ += fvmv.source_;
 933     internalCoeffs_ += fvmv.internalCoeffs_;
 934     boundaryCoeffs_ += fvmv.boundaryCoeffs_;
 935
 936     if (faceFluxCorrectionPtr_ && fvmv.faceFluxCorrectionPtr_)
 937     {
 938         *faceFluxCorrectionPtr_ += *fvmv.faceFluxCorrectionPtr_;
 939     }
 940     else if (fvmv.faceFluxCorrectionPtr_)
 941     {
 942         faceFluxCorrectionPtr_ = new
 943         GeometricField<Type, fvsPatchField, surfaceMesh>
 944         (
 945             *fvmv.faceFluxCorrectionPtr_
 946         );
 947     }
 948 }
 949
 950
 951 template<class Type>
 952 void Foam::fvMatrix<Type>::operator+=(const tmp<fvMatrix<Type> >& tfvmv)
 953 {
 954     operator+=(tfvmv());
 955     tfvmv.clear();
 956 }
 957
 958
 959 template<class Type>
 960 void Foam::fvMatrix<Type>::operator-=(const fvMatrix<Type>& fvmv)
 961 {
 962     checkMethod(*this, fvmv, "+=");
 963
 964     dimensions_ -= fvmv.dimensions_;
 965     lduMatrix::operator-=(fvmv);
 966     source_ -= fvmv.source_;
 967     internalCoeffs_ -= fvmv.internalCoeffs_;
 968     boundaryCoeffs_ -= fvmv.boundaryCoeffs_;
 969
 970     if (faceFluxCorrectionPtr_ && fvmv.faceFluxCorrectionPtr_)
 971     {
 972         *faceFluxCorrectionPtr_ -= *fvmv.faceFluxCorrectionPtr_;
 973     }
 974     else if (fvmv.faceFluxCorrectionPtr_)
 975     {
 976         faceFluxCorrectionPtr_ =
 977             new GeometricField<Type, fvsPatchField, surfaceMesh>
 978         (-*fvmv.faceFluxCorrectionPtr_);
 979     }
 980 }
 981
 982
 983 template<class Type>
 984 void Foam::fvMatrix<Type>::operator-=(const tmp<fvMatrix<Type> >& tfvmv)
 985 {
 986     operator-=(tfvmv());
 987     tfvmv.clear();
 988 }
 989
 990
 991 template<class Type>
 992 void Foam::fvMatrix<Type>::operator+=
 993 (
 994     const DimensionedField<Type, volMesh>& su
 995 )
 996 {
 997     checkMethod(*this, su, "+=");
 998     source() -= su.mesh().V()*su.field();
 999 }
1000
1001
1002 template<class Type>
1003 void Foam::fvMatrix<Type>::operator+=
1004 (
1005     const tmp<DimensionedField<Type, volMesh> >& tsu
1006 )
1007 {
1008     operator+=(tsu());
1009     tsu.clear();
1010 }
1011
1012
1013 template<class Type>
1014 void Foam::fvMatrix<Type>::operator+=
1015 (
1016     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
1017 )
1018 {
1019     operator+=(tsu());
1020     tsu.clear();
1021 }
1022
1023
1024 template<class Type>
1025 void Foam::fvMatrix<Type>::operator-=
1026 (
1027     const DimensionedField<Type, volMesh>& su
1028 )
1029 {
1030     checkMethod(*this, su, "-=");
1031     source() += su.mesh().V()*su.field();
1032 }
1033
1034
1035 template<class Type>
1036 void Foam::fvMatrix<Type>::operator-=
1037 (
1038     const tmp<DimensionedField<Type, volMesh> >& tsu
1039 )
1040 {
1041     operator-=(tsu());
1042     tsu.clear();
1043 }
1044
1045
1046 template<class Type>
1047 void Foam::fvMatrix<Type>::operator-=
1048 (
1049     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
1050 )
1051 {
1052     operator-=(tsu());
1053     tsu.clear();
1054 }
1055
1056
1057 template<class Type>
1058 void Foam::fvMatrix<Type>::operator+=
1059 (
1060     const dimensioned<Type>& su
1061 )
1062 {
1063     source() -= psi().mesh().V()*su;
1064 }
1065
1066
1067 template<class Type>
1068 void Foam::fvMatrix<Type>::operator-=
1069 (
1070     const dimensioned<Type>& su
1071 )
1072 {
1073     source() += psi().mesh().V()*su;
1074 }
1075
1076
1077 template<class Type>
1078 void Foam::fvMatrix<Type>::operator+=
1079 (
1080     const zeroField&
1081 )
1082 {}
1083
1084
1085 template<class Type>
1086 void Foam::fvMatrix<Type>::operator-=
1087 (
1088     const zeroField&
1089 )
1090 {}
1091
1092
1093 template<class Type>
1094 void Foam::fvMatrix<Type>::operator*=
1095 (
1096     const DimensionedField<scalar, volMesh>& dsf
1097 )
1098 {
1099     dimensions_ *= dsf.dimensions();
1100     lduMatrix::operator*=(dsf.field());
1101     source_ *= dsf.field();
1102
1103     forAll(boundaryCoeffs_, patchI)
1104     {
1105         scalarField psf
1106         (
1107             dsf.mesh().boundary()[patchI].patchInternalField(dsf.field())
1108         );
1109
1110         internalCoeffs_[patchI] *= psf;
1111         boundaryCoeffs_[patchI] *= psf;
1112     }
1113
1114     if (faceFluxCorrectionPtr_)
1115     {
1116         FatalErrorIn
1117         (
1118             "fvMatrix<Type>::operator*="
1119             "(const DimensionedField<scalar, volMesh>&)"
1120         )   << "cannot scale a matrix containing a faceFluxCorrection"
1121             << abort(FatalError);
1122     }
1123 }
1124
1125
1126 template<class Type>
1127 void Foam::fvMatrix<Type>::operator*=
1128 (
1129     const tmp<DimensionedField<scalar, volMesh> >& tdsf
1130 )
1131 {
1132     operator*=(tdsf());
1133     tdsf.clear();
1134 }
1135
1136
1137 template<class Type>
1138 void Foam::fvMatrix<Type>::operator*=
1139 (
1140     const tmp<volScalarField>& tvsf
1141 )
1142 {
1143     operator*=(tvsf());
1144     tvsf.clear();
1145 }
1146
1147
1148 template<class Type>
1149 void Foam::fvMatrix<Type>::operator*=
1150 (
1151     const dimensioned<scalar>& ds
1152 )
1153 {
1154     dimensions_ *= ds.dimensions();
1155     lduMatrix::operator*=(ds.value());
1156     source_ *= ds.value();
1157     internalCoeffs_ *= ds.value();
1158     boundaryCoeffs_ *= ds.value();
1159
1160     if (faceFluxCorrectionPtr_)
1161     {
1162         *faceFluxCorrectionPtr_ *= ds.value();
1163     }
1164 }
1165
1166
1167 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
1168
1169 template<class Type>
1170 void Foam::checkMethod
1171 (
1172     const fvMatrix<Type>& fvm1,
1173     const fvMatrix<Type>& fvm2,
1174     const char* op
1175 )
1176 {
1177     if (&fvm1.psi() != &fvm2.psi())
1178     {
1179         FatalErrorIn
1180         (
1181             "checkMethod(const fvMatrix<Type>&, const fvMatrix<Type>&)"
1182         )   << "incompatible fields for operation "
1183             << endl << "    "
1184             << "[" << fvm1.psi().name() << "] "
1185             << op
1186             << " [" << fvm2.psi().name() << "]"
1187             << abort(FatalError);
1188     }
1189
1190     if (dimensionSet::debug && fvm1.dimensions() != fvm2.dimensions())
1191     {
1192         FatalErrorIn
1193         (
1194             "checkMethod(const fvMatrix<Type>&, const fvMatrix<Type>&)"
1195         )   << "incompatible dimensions for operation "
1196             << endl << "    "
1197             << "[" << fvm1.psi().name() << fvm1.dimensions()/dimVolume << " ] "
1198             << op
1199             << " [" << fvm2.psi().name() << fvm2.dimensions()/dimVolume << " ]"
1200             << abort(FatalError);
1201     }
1202 }
1203
1204
1205 template<class Type>
1206 void Foam::checkMethod
1207 (
1208     const fvMatrix<Type>& fvm,
1209     const DimensionedField<Type, volMesh>& df,
1210     const char* op
1211 )
1212 {
1213     if (dimensionSet::debug && fvm.dimensions()/dimVolume != df.dimensions())
1214     {
1215         FatalErrorIn
1216         (
1217             "checkMethod(const fvMatrix<Type>&, const GeometricField<Type, "
1218             "fvPatchField, volMesh>&)"
1219         )   <<  "incompatible dimensions for operation "
1220             << endl << "    "
1221             << "[" << fvm.psi().name() << fvm.dimensions()/dimVolume << " ] "
1222             << op
1223             << " [" << df.name() << df.dimensions() << " ]"
1224             << abort(FatalError);
1225     }
1226 }
1227
1228
1229 template<class Type>
1230 void Foam::checkMethod
1231 (
1232     const fvMatrix<Type>& fvm,
1233     const dimensioned<Type>& dt,
1234     const char* op
1235 )
1236 {
1237     if (dimensionSet::debug && fvm.dimensions()/dimVolume != dt.dimensions())
1238     {
1239         FatalErrorIn
1240         (
1241             "checkMethod(const fvMatrix<Type>&, const dimensioned<Type>&)"
1242         )   << "incompatible dimensions for operation "
1243             << endl << "    "
1244             << "[" << fvm.psi().name() << fvm.dimensions()/dimVolume << " ] "
1245             << op
1246             << " [" << dt.name() << dt.dimensions() << " ]"
1247             << abort(FatalError);
1248     }
1249 }
1250
1251
1252 template<class Type>
1253 Foam::lduMatrix::solverPerformance Foam::solve
1254 (
1255     fvMatrix<Type>& fvm,
1256     const dictionary& solverControls
1257 )
1258 {
1259     return fvm.solve(solverControls);
1260 }
1261
1262 template<class Type>
1263 Foam::lduMatrix::solverPerformance Foam::solve
1264 (
1265     const tmp<fvMatrix<Type> >& tfvm,
1266     const dictionary& solverControls
1267 )
1268 {
1269     lduMatrix::solverPerformance solverPerf =
1270         const_cast<fvMatrix<Type>&>(tfvm()).solve(solverControls);
1271
1272     tfvm.clear();
1273
1274     return solverPerf;
1275 }
1276
1277
1278 template<class Type>
1279 Foam::lduMatrix::solverPerformance Foam::solve(fvMatrix<Type>& fvm)
1280 {
1281     return fvm.solve();
1282 }
1283
1284 template<class Type>
1285 Foam::lduMatrix::solverPerformance Foam::solve(const tmp<fvMatrix<Type> >& tfvm)
1286 {
1287     lduMatrix::solverPerformance solverPerf =
1288         const_cast<fvMatrix<Type>&>(tfvm()).solve();
1289
1290     tfvm.clear();
1291
1292     return solverPerf;
1293 }
1294
1295
1296 template<class Type>
1297 Foam::tmp<Foam::fvMatrix<Type> > Foam::correction
1298 (
1299     const fvMatrix<Type>& A
1300 )
1301 {
1302     tmp<Foam::fvMatrix<Type> > tAcorr = A - (A & A.psi());
1303
1304     if
1305     (
1306         (A.hasUpper() || A.hasLower())
1307      && A.psi().mesh().fluxRequired(A.psi().name())
1308     )
1309     {
1310         tAcorr().faceFluxCorrectionPtr() = (-A.flux()).ptr();
1311     }
1312
1313     return tAcorr;
1314 }
1315
1316
1317 template<class Type>
1318 Foam::tmp<Foam::fvMatrix<Type> > Foam::correction
1319 (
1320     const tmp<fvMatrix<Type> >& tA
1321 )
1322 {
1323     tmp<Foam::fvMatrix<Type> > tAcorr = tA - (tA() & tA().psi());
1324
1325     // Note the matrix coefficients are still that of matrix A
1326     const fvMatrix<Type>& A = tAcorr();
1327
1328     if
1329     (
1330         (A.hasUpper() || A.hasLower())
1331      && A.psi().mesh().fluxRequired(A.psi().name())
1332     )
1333     {
1334         tAcorr().faceFluxCorrectionPtr() = (-A.flux()).ptr();
1335     }
1336
1337     return tAcorr;
1338 }
1339
1340
1341 // * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
1342
1343 template<class Type>
1344 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1345 (
1346     const fvMatrix<Type>& A,
1347     const fvMatrix<Type>& B
1348 )
1349 {
1350     checkMethod(A, B, "==");
1351     return (A - B);
1352 }
1353
1354 template<class Type>
1355 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1356 (
1357     const tmp<fvMatrix<Type> >& tA,
1358     const fvMatrix<Type>& B
1359 )
1360 {
1361     checkMethod(tA(), B, "==");
1362     return (tA - B);
1363 }
1364
1365 template<class Type>
1366 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1367 (
1368     const fvMatrix<Type>& A,
1369     const tmp<fvMatrix<Type> >& tB
1370 )
1371 {
1372     checkMethod(A, tB(), "==");
1373     return (A - tB);
1374 }
1375
1376 template<class Type>
1377 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1378 (
1379     const tmp<fvMatrix<Type> >& tA,
1380     const tmp<fvMatrix<Type> >& tB
1381 )
1382 {
1383     checkMethod(tA(), tB(), "==");
1384     return (tA - tB);
1385 }
1386
1387 template<class Type>
1388 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1389 (
1390     const fvMatrix<Type>& A,
1391     const DimensionedField<Type, volMesh>& su
1392 )
1393 {
1394     checkMethod(A, su, "==");
1395     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1396     tC().source() += su.mesh().V()*su.field();
1397     return tC;
1398 }
1399
1400 template<class Type>
1401 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1402 (
1403     const fvMatrix<Type>& A,
1404     const tmp<DimensionedField<Type, volMesh> >& tsu
1405 )
1406 {
1407     checkMethod(A, tsu(), "==");
1408     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1409     tC().source() += tsu().mesh().V()*tsu().field();
1410     tsu.clear();
1411     return tC;
1412 }
1413
1414 template<class Type>
1415 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1416 (
1417     const fvMatrix<Type>& A,
1418     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
1419 )
1420 {
1421     checkMethod(A, tsu(), "==");
1422     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1423     tC().source() += tsu().mesh().V()*tsu().internalField();
1424     tsu.clear();
1425     return tC;
1426 }
1427
1428 template<class Type>
1429 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1430 (
1431     const tmp<fvMatrix<Type> >& tA,
1432     const DimensionedField<Type, volMesh>& su
1433 )
1434 {
1435     checkMethod(tA(), su, "==");
1436     tmp<fvMatrix<Type> > tC(tA.ptr());
1437     tC().source() += su.mesh().V()*su.field();
1438     return tC;
1439 }
1440
1441 template<class Type>
1442 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1443 (
1444     const tmp<fvMatrix<Type> >& tA,
1445     const tmp<DimensionedField<Type, volMesh> >& tsu
1446 )
1447 {
1448     checkMethod(tA(), tsu(), "==");
1449     tmp<fvMatrix<Type> > tC(tA.ptr());
1450     tC().source() += tsu().mesh().V()*tsu().field();
1451     tsu.clear();
1452     return tC;
1453 }
1454
1455 template<class Type>
1456 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1457 (
1458     const tmp<fvMatrix<Type> >& tA,
1459     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
1460 )
1461 {
1462     checkMethod(tA(), tsu(), "==");
1463     tmp<fvMatrix<Type> > tC(tA.ptr());
1464     tC().source() += tsu().mesh().V()*tsu().internalField();
1465     tsu.clear();
1466     return tC;
1467 }
1468
1469 template<class Type>
1470 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1471 (
1472     const fvMatrix<Type>& A,
1473     const dimensioned<Type>& su
1474 )
1475 {
1476     checkMethod(A, su, "==");
1477     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1478     tC().source() += A.psi().mesh().V()*su.value();
1479     return tC;
1480 }
1481
1482 template<class Type>
1483 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1484 (
1485     const tmp<fvMatrix<Type> >& tA,
1486     const dimensioned<Type>& su
1487 )
1488 {
1489     checkMethod(tA(), su, "==");
1490     tmp<fvMatrix<Type> > tC(tA.ptr());
1491     tC().source() += tC().psi().mesh().V()*su.value();
1492     return tC;
1493 }
1494
1495 template<class Type>
1496 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1497 (
1498     const fvMatrix<Type>& A,
1499     const zeroField&
1500 )
1501 {
1502     return A;
1503 }
1504
1505
1506 template<class Type>
1507 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1508 (
1509     const tmp<fvMatrix<Type> >& tA,
1510     const zeroField&
1511 )
1512 {
1513     return tA;
1514 }
1515
1516
1517 template<class Type>
1518 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1519 (
1520     const fvMatrix<Type>& A
1521 )
1522 {
1523     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1524     tC().negate();
1525     return tC;
1526 }
1527
1528 template<class Type>
1529 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1530 (
1531     const tmp<fvMatrix<Type> >& tA
1532 )
1533 {
1534     tmp<fvMatrix<Type> > tC(tA.ptr());
1535     tC().negate();
1536     return tC;
1537 }
1538
1539
1540 template<class Type>
1541 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1542 (
1543     const fvMatrix<Type>& A,
1544     const fvMatrix<Type>& B
1545 )
1546 {
1547     checkMethod(A, B, "+");
1548     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1549     tC() += B;
1550     return tC;
1551 }
1552
1553 template<class Type>
1554 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1555 (
1556     const tmp<fvMatrix<Type> >& tA,
1557     const fvMatrix<Type>& B
1558 )
1559 {
1560     checkMethod(tA(), B, "+");
1561     tmp<fvMatrix<Type> > tC(tA.ptr());
1562     tC() += B;
1563     return tC;
1564 }
1565
1566 template<class Type>
1567 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1568 (
1569     const fvMatrix<Type>& A,
1570     const tmp<fvMatrix<Type> >& tB
1571 )
1572 {
1573     checkMethod(A, tB(), "+");
1574     tmp<fvMatrix<Type> > tC(tB.ptr());
1575     tC() += A;
1576     return tC;
1577 }
1578
1579 template<class Type>
1580 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1581 (
1582     const tmp<fvMatrix<Type> >& tA,
1583     const tmp<fvMatrix<Type> >& tB
1584 )
1585 {
1586     checkMethod(tA(), tB(), "+");
1587     tmp<fvMatrix<Type> > tC(tA.ptr());
1588     tC() += tB();
1589     tB.clear();
1590     return tC;
1591 }
1592
1593 template<class Type>
1594 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1595 (
1596     const fvMatrix<Type>& A,
1597     const DimensionedField<Type, volMesh>& su
1598 )
1599 {
1600     checkMethod(A, su, "+");
1601     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1602     tC().source() -= su.mesh().V()*su.field();
1603     return tC;
1604 }
1605
1606 template<class Type>
1607 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1608 (
1609     const fvMatrix<Type>& A,
1610     const tmp<DimensionedField<Type, volMesh> >& tsu
1611 )
1612 {
1613     checkMethod(A, tsu(), "+");
1614     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1615     tC().source() -= tsu().mesh().V()*tsu().field();
1616     tsu.clear();
1617     return tC;
1618 }
1619
1620 template<class Type>
1621 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1622 (
1623     const fvMatrix<Type>& A,
1624     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
1625 )
1626 {
1627     checkMethod(A, tsu(), "+");
1628     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1629     tC().source() -= tsu().mesh().V()*tsu().internalField();
1630     tsu.clear();
1631     return tC;
1632 }
1633
1634 template<class Type>
1635 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1636 (
1637     const tmp<fvMatrix<Type> >& tA,
1638     const DimensionedField<Type, volMesh>& su
1639 )
1640 {
1641     checkMethod(tA(), su, "+");
1642     tmp<fvMatrix<Type> > tC(tA.ptr());
1643     tC().source() -= su.mesh().V()*su.field();
1644     return tC;
1645 }
1646
1647 template<class Type>
1648 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1649 (
1650     const tmp<fvMatrix<Type> >& tA,
1651     const tmp<DimensionedField<Type, volMesh> >& tsu
1652 )
1653 {
1654     checkMethod(tA(), tsu(), "+");
1655     tmp<fvMatrix<Type> > tC(tA.ptr());
1656     tC().source() -= tsu().mesh().V()*tsu().field();
1657     tsu.clear();
1658     return tC;
1659 }
1660
1661 template<class Type>
1662 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1663 (
1664     const tmp<fvMatrix<Type> >& tA,
1665     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
1666 )
1667 {
1668     checkMethod(tA(), tsu(), "+");
1669     tmp<fvMatrix<Type> > tC(tA.ptr());
1670     tC().source() -= tsu().mesh().V()*tsu().internalField();
1671     tsu.clear();
1672     return tC;
1673 }
1674
1675 template<class Type>
1676 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1677 (
1678     const DimensionedField<Type, volMesh>& su,
1679     const fvMatrix<Type>& A
1680 )
1681 {
1682     checkMethod(A, su, "+");
1683     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1684     tC().source() -= su.mesh().V()*su.field();
1685     return tC;
1686 }
1687
1688 template<class Type>
1689 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1690 (
1691     const tmp<DimensionedField<Type, volMesh> >& tsu,
1692     const fvMatrix<Type>& A
1693 )
1694 {
1695     checkMethod(A, tsu(), "+");
1696     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1697     tC().source() -= tsu().mesh().V()*tsu().field();
1698     tsu.clear();
1699     return tC;
1700 }
1701
1702 template<class Type>
1703 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1704 (
1705     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu,
1706     const fvMatrix<Type>& A
1707 )
1708 {
1709     checkMethod(A, tsu(), "+");
1710     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1711     tC().source() -= tsu().mesh().V()*tsu().internalField();
1712     tsu.clear();
1713     return tC;
1714 }
1715
1716 template<class Type>
1717 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1718 (
1719     const DimensionedField<Type, volMesh>& su,
1720     const tmp<fvMatrix<Type> >& tA
1721 )
1722 {
1723     checkMethod(tA(), su, "+");
1724     tmp<fvMatrix<Type> > tC(tA.ptr());
1725     tC().source() -= su.mesh().V()*su.field();
1726     return tC;
1727 }
1728
1729 template<class Type>
1730 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1731 (
1732     const tmp<DimensionedField<Type, volMesh> >& tsu,
1733     const tmp<fvMatrix<Type> >& tA
1734 )
1735 {
1736     checkMethod(tA(), tsu(), "+");
1737     tmp<fvMatrix<Type> > tC(tA.ptr());
1738     tC().source() -= tsu().mesh().V()*tsu().field();
1739     tsu.clear();
1740     return tC;
1741 }
1742
1743 template<class Type>
1744 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1745 (
1746     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu,
1747     const tmp<fvMatrix<Type> >& tA
1748 )
1749 {
1750     checkMethod(tA(), tsu(), "+");
1751     tmp<fvMatrix<Type> > tC(tA.ptr());
1752     tC().source() -= tsu().mesh().V()*tsu().internalField();
1753     tsu.clear();
1754     return tC;
1755 }
1756
1757
1758 template<class Type>
1759 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1760 (
1761     const fvMatrix<Type>& A,
1762     const fvMatrix<Type>& B
1763 )
1764 {
1765     checkMethod(A, B, "-");
1766     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1767     tC() -= B;
1768     return tC;
1769 }
1770
1771 template<class Type>
1772 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1773 (
1774     const tmp<fvMatrix<Type> >& tA,
1775     const fvMatrix<Type>& B
1776 )
1777 {
1778     checkMethod(tA(), B, "-");
1779     tmp<fvMatrix<Type> > tC(tA.ptr());
1780     tC() -= B;
1781     return tC;
1782 }
1783
1784 template<class Type>
1785 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1786 (
1787     const fvMatrix<Type>& A,
1788     const tmp<fvMatrix<Type> >& tB
1789 )
1790 {
1791     checkMethod(A, tB(), "-");
1792     tmp<fvMatrix<Type> > tC(tB.ptr());
1793     tC() -= A;
1794     tC().negate();
1795     return tC;
1796 }
1797
1798 template<class Type>
1799 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1800 (
1801     const tmp<fvMatrix<Type> >& tA,
1802     const tmp<fvMatrix<Type> >& tB
1803 )
1804 {
1805     checkMethod(tA(), tB(), "-");
1806     tmp<fvMatrix<Type> > tC(tA.ptr());
1807     tC() -= tB();
1808     tB.clear();
1809     return tC;
1810 }
1811
1812 template<class Type>
1813 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1814 (
1815     const fvMatrix<Type>& A,
1816     const DimensionedField<Type, volMesh>& su
1817 )
1818 {
1819     checkMethod(A, su, "-");
1820     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1821     tC().source() += su.mesh().V()*su.field();
1822     return tC;
1823 }
1824
1825 template<class Type>
1826 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1827 (
1828     const fvMatrix<Type>& A,
1829     const tmp<DimensionedField<Type, volMesh> >& tsu
1830 )
1831 {
1832     checkMethod(A, tsu(), "-");
1833     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1834     tC().source() += tsu().mesh().V()*tsu().field();
1835     tsu.clear();
1836     return tC;
1837 }
1838
1839 template<class Type>
1840 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1841 (
1842     const fvMatrix<Type>& A,
1843     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
1844 )
1845 {
1846     checkMethod(A, tsu(), "-");
1847     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1848     tC().source() += tsu().mesh().V()*tsu().internalField();
1849     tsu.clear();
1850     return tC;
1851 }
1852
1853 template<class Type>
1854 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1855 (
1856     const tmp<fvMatrix<Type> >& tA,
1857     const DimensionedField<Type, volMesh>& su
1858 )
1859 {
1860     checkMethod(tA(), su, "-");
1861     tmp<fvMatrix<Type> > tC(tA.ptr());
1862     tC().source() += su.mesh().V()*su.field();
1863     return tC;
1864 }
1865
1866 template<class Type>
1867 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1868 (
1869     const tmp<fvMatrix<Type> >& tA,
1870     const tmp<DimensionedField<Type, volMesh> >& tsu
1871 )
1872 {
1873     checkMethod(tA(), tsu(), "-");
1874     tmp<fvMatrix<Type> > tC(tA.ptr());
1875     tC().source() += tsu().mesh().V()*tsu().field();
1876     tsu.clear();
1877     return tC;
1878 }
1879
1880 template<class Type>
1881 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1882 (
1883     const tmp<fvMatrix<Type> >& tA,
1884     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
1885 )
1886 {
1887     checkMethod(tA(), tsu(), "-");
1888     tmp<fvMatrix<Type> > tC(tA.ptr());
1889     tC().source() += tsu().mesh().V()*tsu().internalField();
1890     tsu.clear();
1891     return tC;
1892 }
1893
1894 template<class Type>
1895 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1896 (
1897     const DimensionedField<Type, volMesh>& su,
1898     const fvMatrix<Type>& A
1899 )
1900 {
1901     checkMethod(A, su, "-");
1902     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1903     tC().negate();
1904     tC().source() -= su.mesh().V()*su.field();
1905     return tC;
1906 }
1907
1908 template<class Type>
1909 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1910 (
1911     const tmp<DimensionedField<Type, volMesh> >& tsu,
1912     const fvMatrix<Type>& A
1913 )
1914 {
1915     checkMethod(A, tsu(), "-");
1916     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1917     tC().negate();
1918     tC().source() -= tsu().mesh().V()*tsu().field();
1919     tsu.clear();
1920     return tC;
1921 }
1922
1923 template<class Type>
1924 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1925 (
1926     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu,
1927     const fvMatrix<Type>& A
1928 )
1929 {
1930     checkMethod(A, tsu(), "-");
1931     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1932     tC().negate();
1933     tC().source() -= tsu().mesh().V()*tsu().internalField();
1934     tsu.clear();
1935     return tC;
1936 }
1937
1938 template<class Type>
1939 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1940 (
1941     const DimensionedField<Type, volMesh>& su,
1942     const tmp<fvMatrix<Type> >& tA
1943 )
1944 {
1945     checkMethod(tA(), su, "-");
1946     tmp<fvMatrix<Type> > tC(tA.ptr());
1947     tC().negate();
1948     tC().source() -= su.mesh().V()*su.field();
1949     return tC;
1950 }
1951
1952 template<class Type>
1953 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1954 (
1955     const tmp<DimensionedField<Type, volMesh> >& tsu,
1956     const tmp<fvMatrix<Type> >& tA
1957 )
1958 {
1959     checkMethod(tA(), tsu(), "-");
1960     tmp<fvMatrix<Type> > tC(tA.ptr());
1961     tC().negate();
1962     tC().source() -= tsu().mesh().V()*tsu().field();
1963     tsu.clear();
1964     return tC;
1965 }
1966
1967 template<class Type>
1968 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1969 (
1970     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu,
1971     const tmp<fvMatrix<Type> >& tA
1972 )
1973 {
1974     checkMethod(tA(), tsu(), "-");
1975     tmp<fvMatrix<Type> > tC(tA.ptr());
1976     tC().negate();
1977     tC().source() -= tsu().mesh().V()*tsu().internalField();
1978     tsu.clear();
1979     return tC;
1980 }
1981
1982 template<class Type>
1983 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1984 (
1985     const fvMatrix<Type>& A,
1986     const dimensioned<Type>& su
1987 )
1988 {
1989     checkMethod(A, su, "+");
1990     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1991     tC().source() -= su.value()*A.psi().mesh().V();
1992     return tC;
1993 }
1994
1995 template<class Type>
1996 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1997 (
1998     const tmp<fvMatrix<Type> >& tA,
1999     const dimensioned<Type>& su
2000 )
2001 {
2002     checkMethod(tA(), su, "+");
2003     tmp<fvMatrix<Type> > tC(tA.ptr());
2004     tC().source() -= su.value()*tC().psi().mesh().V();
2005     return tC;
2006 }
2007
2008 template<class Type>
2009 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
2010 (
2011     const dimensioned<Type>& su,
2012     const fvMatrix<Type>& A
2013 )
2014 {
2015     checkMethod(A, su, "+");
2016     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
2017     tC().source() -= su.value()*A.psi().mesh().V();
2018     return tC;
2019 }
2020
2021 template<class Type>
2022 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
2023 (
2024     const dimensioned<Type>& su,
2025     const tmp<fvMatrix<Type> >& tA
2026 )
2027 {
2028     checkMethod(tA(), su, "+");
2029     tmp<fvMatrix<Type> > tC(tA.ptr());
2030     tC().source() -= su.value()*tC().psi().mesh().V();
2031     return tC;
2032 }
2033
2034 template<class Type>
2035 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
2036 (
2037     const fvMatrix<Type>& A,
2038     const dimensioned<Type>& su
2039 )
2040 {
2041     checkMethod(A, su, "-");
2042     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
2043     tC().source() += su.value()*tC().psi().mesh().V();
2044     return tC;
2045 }
2046
2047 template<class Type>
2048 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
2049 (
2050     const tmp<fvMatrix<Type> >& tA,
2051     const dimensioned<Type>& su
2052 )
2053 {
2054     checkMethod(tA(), su, "-");
2055     tmp<fvMatrix<Type> > tC(tA.ptr());
2056     tC().source() += su.value()*tC().psi().mesh().V();
2057     return tC;
2058 }
2059
2060 template<class Type>
2061 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
2062 (
2063     const dimensioned<Type>& su,
2064     const fvMatrix<Type>& A
2065 )
2066 {
2067     checkMethod(A, su, "-");
2068     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
2069     tC().negate();
2070     tC().source() -= su.value()*A.psi().mesh().V();
2071     return tC;
2072 }
2073
2074 template<class Type>
2075 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
2076 (
2077     const dimensioned<Type>& su,
2078     const tmp<fvMatrix<Type> >& tA
2079 )
2080 {
2081     checkMethod(tA(), su, "-");
2082     tmp<fvMatrix<Type> > tC(tA.ptr());
2083     tC().negate();
2084     tC().source() -= su.value()*tC().psi().mesh().V();
2085     return tC;
2086 }
2087
2088
2089 template<class Type>
2090 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator*
2091 (
2092     const DimensionedField<scalar, volMesh>& dsf,
2093     const fvMatrix<Type>& A
2094 )
2095 {
2096     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
2097     tC() *= dsf;
2098     return tC;
2099 }
2100
2101 template<class Type>
2102 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator*
2103 (
2104     const tmp< DimensionedField<scalar, volMesh> >& tdsf,
2105     const fvMatrix<Type>& A
2106 )
2107 {
2108     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
2109     tC() *= tdsf;
2110     return tC;
2111 }
2112
2113 template<class Type>
2114 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator*
2115 (
2116     const tmp<volScalarField>& tvsf,
2117     const fvMatrix<Type>& A
2118 )
2119 {
2120     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
2121     tC() *= tvsf;
2122     return tC;
2123 }
2124
2125 template<class Type>
2126 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator*
2127 (
2128     const DimensionedField<scalar, volMesh>& dsf,
2129     const tmp<fvMatrix<Type> >& tA
2130 )
2131 {
2132     tmp<fvMatrix<Type> > tC(tA.ptr());
2133     tC() *= dsf;
2134     return tC;
2135 }
2136
2137 template<class Type>
2138 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator*
2139 (
2140     const tmp<DimensionedField<scalar, volMesh> >& tdsf,
2141     const tmp<fvMatrix<Type> >& tA
2142 )
2143 {
2144     tmp<fvMatrix<Type> > tC(tA.ptr());
2145     tC() *= tdsf;
2146     return tC;
2147 }
2148
2149 template<class Type>
2150 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator*
2151 (
2152     const tmp<volScalarField>& tvsf,
2153     const tmp<fvMatrix<Type> >& tA
2154 )
2155 {
2156     tmp<fvMatrix<Type> > tC(tA.ptr());
2157     tC() *= tvsf;
2158     return tC;
2159 }
2160
2161 template<class Type>
2162 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator*
2163 (
2164     const dimensioned<scalar>& ds,
2165     const fvMatrix<Type>& A
2166 )
2167 {
2168     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
2169     tC() *= ds;
2170     return tC;
2171 }
2172
2173 template<class Type>
2174 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator*
2175 (
2176     const dimensioned<scalar>& ds,
2177     const tmp<fvMatrix<Type> >& tA
2178 )
2179 {
2180     tmp<fvMatrix<Type> > tC(tA.ptr());
2181     tC() *= ds;
2182     return tC;
2183 }
2184
2185
2186 template<class Type>
2187 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
2188 Foam::operator&
2189 (
2190     const fvMatrix<Type>& M,
2191     const DimensionedField<Type, volMesh>& psi
2192 )
2193 {
2194     tmp<GeometricField<Type, fvPatchField, volMesh> > tMphi
2195     (
2196         new GeometricField<Type, fvPatchField, volMesh>
2197         (
2198             IOobject
2199             (
2200                 "M&" + psi.name(),
2201                 psi.instance(),
2202                 psi.mesh(),
2203                 IOobject::NO_READ,
2204                 IOobject::NO_WRITE
2205             ),
2206             psi.mesh(),
2207             M.dimensions()/dimVol,
2208             zeroGradientFvPatchScalarField::typeName
2209         )
2210     );
2211     GeometricField<Type, fvPatchField, volMesh>& Mphi = tMphi();
2212
2213     // Loop over field components
2214     for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
2215     {
2216         scalarField psiCmpt = psi.field().component(cmpt);
2217         scalarField boundaryDiagCmpt(M.diag());
2218         M.addBoundaryDiag(boundaryDiagCmpt, cmpt);
2219         Mphi.internalField().replace(cmpt, -boundaryDiagCmpt*psiCmpt);
2220     }
2221
2222     Mphi.internalField() += M.lduMatrix::H(psi.field()) + M.source();
2223     M.addBoundarySource(Mphi.internalField());
2224
2225     Mphi.internalField() /= -psi.mesh().V();
2226     Mphi.correctBoundaryConditions();
2227
2228     return tMphi;
2229 }
2230
2231 template<class Type>
2232 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
2233 Foam::operator&
2234 (
2235     const fvMatrix<Type>& M,
2236     const tmp<DimensionedField<Type, volMesh> >& tpsi
2237 )
2238 {
2239     tmp<GeometricField<Type, fvPatchField, volMesh> > tMpsi = M & tpsi();
2240     tpsi.clear();
2241     return tMpsi;
2242 }
2243
2244 template<class Type>
2245 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
2246 Foam::operator&
2247 (
2248     const fvMatrix<Type>& M,
2249     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tpsi
2250 )
2251 {
2252     tmp<GeometricField<Type, fvPatchField, volMesh> > tMpsi = M & tpsi();
2253     tpsi.clear();
2254     return tMpsi;
2255 }
2256
2257 template<class Type>
2258 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
2259 Foam::operator&
2260 (
2261     const tmp<fvMatrix<Type> >& tM,
2262     const DimensionedField<Type, volMesh>& psi
2263 )
2264 {
2265     tmp<GeometricField<Type, fvPatchField, volMesh> > tMpsi = tM() & psi;
2266     tM.clear();
2267     return tMpsi;
2268 }
2269
2270 template<class Type>
2271 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
2272 Foam::operator&
2273 (
2274     const tmp<fvMatrix<Type> >& tM,
2275     const tmp<DimensionedField<Type, volMesh> >& tpsi
2276 )
2277 {
2278     tmp<GeometricField<Type, fvPatchField, volMesh> > tMpsi = tM() & tpsi();
2279     tM.clear();
2280     tpsi.clear();
2281     return tMpsi;
2282 }
2283
2284 template<class Type>
2285 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
2286 Foam::operator&
2287 (
2288     const tmp<fvMatrix<Type> >& tM,
2289     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tpsi
2290 )
2291 {
2292     tmp<GeometricField<Type, fvPatchField, volMesh> > tMpsi = tM() & tpsi();
2293     tM.clear();
2294     tpsi.clear();
2295     return tMpsi;
2296 }
2297
2298
2299 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
2300
2301 template<class Type>
2302 Foam::Ostream& Foam::operator<<(Ostream& os, const fvMatrix<Type>& fvm)
2303 {
2304     os  << static_cast<const lduMatrix&>(fvm) << nl
2305         << fvm.dimensions_ << nl
2306         << fvm.source_ << nl
2307         << fvm.internalCoeffs_ << nl
2308         << fvm.boundaryCoeffs_ << endl;
2309
2310     os.check("Ostream& operator<<(Ostream&, fvMatrix<Type>&");
2311
2312     return os;
2313 }
2314
2315
2316 // * * * * * * * * * * * * * * * * Solvers * * * * * * * * * * * * * * * * * //
2317
2318 #include "fvMatrixSolve.C"
2319
2320 // ************************************************************************* //