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(const tmp<fvMatrix<Type> >& tfvm)
1272 {
1273     lduMatrix::solverPerformance solverPerf =
1274         const_cast<fvMatrix<Type>&>(tfvm()).solve();
1275
1276     tfvm.clear();
1277
1278     return solverPerf;
1279 }
1280
1281
1282 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
1283
1284 template<class Type>
1285 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1286 (
1287     const fvMatrix<Type>& A,
1288     const fvMatrix<Type>& B
1289 )
1290 {
1291     checkMethod(A, B, "==");
1292     return (A - B);
1293 }
1294
1295 template<class Type>
1296 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1297 (
1298     const tmp<fvMatrix<Type> >& tA,
1299     const fvMatrix<Type>& B
1300 )
1301 {
1302     checkMethod(tA(), B, "==");
1303     return (tA - B);
1304 }
1305
1306 template<class Type>
1307 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1308 (
1309     const fvMatrix<Type>& A,
1310     const tmp<fvMatrix<Type> >& tB
1311 )
1312 {
1313     checkMethod(A, tB(), "==");
1314     return (A - tB);
1315 }
1316
1317 template<class Type>
1318 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1319 (
1320     const tmp<fvMatrix<Type> >& tA,
1321     const tmp<fvMatrix<Type> >& tB
1322 )
1323 {
1324     checkMethod(tA(), tB(), "==");
1325     return (tA - tB);
1326 }
1327
1328 template<class Type>
1329 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1330 (
1331     const fvMatrix<Type>& A,
1332     const DimensionedField<Type, volMesh>& su
1333 )
1334 {
1335     checkMethod(A, su, "==");
1336     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1337     tC().source() += su.mesh().V()*su.field();
1338     return tC;
1339 }
1340
1341 template<class Type>
1342 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1343 (
1344     const fvMatrix<Type>& A,
1345     const tmp<DimensionedField<Type, volMesh> >& tsu
1346 )
1347 {
1348     checkMethod(A, tsu(), "==");
1349     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1350     tC().source() += tsu().mesh().V()*tsu().field();
1351     tsu.clear();
1352     return tC;
1353 }
1354
1355 template<class Type>
1356 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1357 (
1358     const fvMatrix<Type>& A,
1359     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
1360 )
1361 {
1362     checkMethod(A, tsu(), "==");
1363     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1364     tC().source() += tsu().mesh().V()*tsu().internalField();
1365     tsu.clear();
1366     return tC;
1367 }
1368
1369 template<class Type>
1370 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1371 (
1372     const tmp<fvMatrix<Type> >& tA,
1373     const DimensionedField<Type, volMesh>& su
1374 )
1375 {
1376     checkMethod(tA(), su, "==");
1377     tmp<fvMatrix<Type> > tC(tA.ptr());
1378     tC().source() += su.mesh().V()*su.field();
1379     return tC;
1380 }
1381
1382 template<class Type>
1383 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1384 (
1385     const tmp<fvMatrix<Type> >& tA,
1386     const tmp<DimensionedField<Type, volMesh> >& tsu
1387 )
1388 {
1389     checkMethod(tA(), tsu(), "==");
1390     tmp<fvMatrix<Type> > tC(tA.ptr());
1391     tC().source() += tsu().mesh().V()*tsu().field();
1392     tsu.clear();
1393     return tC;
1394 }
1395
1396 template<class Type>
1397 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1398 (
1399     const tmp<fvMatrix<Type> >& tA,
1400     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
1401 )
1402 {
1403     checkMethod(tA(), tsu(), "==");
1404     tmp<fvMatrix<Type> > tC(tA.ptr());
1405     tC().source() += tsu().mesh().V()*tsu().internalField();
1406     tsu.clear();
1407     return tC;
1408 }
1409
1410 template<class Type>
1411 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1412 (
1413     const fvMatrix<Type>& A,
1414     const dimensioned<Type>& su
1415 )
1416 {
1417     checkMethod(A, su, "==");
1418     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1419     tC().source() += A.psi().mesh().V()*su.value();
1420     return tC;
1421 }
1422
1423 template<class Type>
1424 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1425 (
1426     const tmp<fvMatrix<Type> >& tA,
1427     const dimensioned<Type>& su
1428 )
1429 {
1430     checkMethod(tA(), su, "==");
1431     tmp<fvMatrix<Type> > tC(tA.ptr());
1432     tC().source() += tC().psi().mesh().V()*su.value();
1433     return tC;
1434 }
1435
1436 template<class Type>
1437 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1438 (
1439     const fvMatrix<Type>& A,
1440     const zeroField&
1441 )
1442 {
1443     return A;
1444 }
1445
1446
1447 template<class Type>
1448 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1449 (
1450     const tmp<fvMatrix<Type> >& tA,
1451     const zeroField&
1452 )
1453 {
1454     return tA;
1455 }
1456
1457
1458 template<class Type>
1459 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1460 (
1461     const fvMatrix<Type>& A
1462 )
1463 {
1464     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1465     tC().negate();
1466     return tC;
1467 }
1468
1469 template<class Type>
1470 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1471 (
1472     const tmp<fvMatrix<Type> >& tA
1473 )
1474 {
1475     tmp<fvMatrix<Type> > tC(tA.ptr());
1476     tC().negate();
1477     return tC;
1478 }
1479
1480
1481 template<class Type>
1482 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1483 (
1484     const fvMatrix<Type>& A,
1485     const fvMatrix<Type>& B
1486 )
1487 {
1488     checkMethod(A, B, "+");
1489     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1490     tC() += B;
1491     return tC;
1492 }
1493
1494 template<class Type>
1495 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1496 (
1497     const tmp<fvMatrix<Type> >& tA,
1498     const fvMatrix<Type>& B
1499 )
1500 {
1501     checkMethod(tA(), B, "+");
1502     tmp<fvMatrix<Type> > tC(tA.ptr());
1503     tC() += B;
1504     return tC;
1505 }
1506
1507 template<class Type>
1508 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1509 (
1510     const fvMatrix<Type>& A,
1511     const tmp<fvMatrix<Type> >& tB
1512 )
1513 {
1514     checkMethod(A, tB(), "+");
1515     tmp<fvMatrix<Type> > tC(tB.ptr());
1516     tC() += A;
1517     return tC;
1518 }
1519
1520 template<class Type>
1521 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1522 (
1523     const tmp<fvMatrix<Type> >& tA,
1524     const tmp<fvMatrix<Type> >& tB
1525 )
1526 {
1527     checkMethod(tA(), tB(), "+");
1528     tmp<fvMatrix<Type> > tC(tA.ptr());
1529     tC() += tB();
1530     tB.clear();
1531     return tC;
1532 }
1533
1534 template<class Type>
1535 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1536 (
1537     const fvMatrix<Type>& A,
1538     const DimensionedField<Type, volMesh>& su
1539 )
1540 {
1541     checkMethod(A, su, "+");
1542     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1543     tC().source() -= su.mesh().V()*su.field();
1544     return tC;
1545 }
1546
1547 template<class Type>
1548 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1549 (
1550     const fvMatrix<Type>& A,
1551     const tmp<DimensionedField<Type, volMesh> >& tsu
1552 )
1553 {
1554     checkMethod(A, tsu(), "+");
1555     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1556     tC().source() -= tsu().mesh().V()*tsu().field();
1557     tsu.clear();
1558     return tC;
1559 }
1560
1561 template<class Type>
1562 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1563 (
1564     const fvMatrix<Type>& A,
1565     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
1566 )
1567 {
1568     checkMethod(A, tsu(), "+");
1569     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1570     tC().source() -= tsu().mesh().V()*tsu().internalField();
1571     tsu.clear();
1572     return tC;
1573 }
1574
1575 template<class Type>
1576 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1577 (
1578     const tmp<fvMatrix<Type> >& tA,
1579     const DimensionedField<Type, volMesh>& su
1580 )
1581 {
1582     checkMethod(tA(), su, "+");
1583     tmp<fvMatrix<Type> > tC(tA.ptr());
1584     tC().source() -= su.mesh().V()*su.field();
1585     return tC;
1586 }
1587
1588 template<class Type>
1589 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1590 (
1591     const tmp<fvMatrix<Type> >& tA,
1592     const tmp<DimensionedField<Type, volMesh> >& tsu
1593 )
1594 {
1595     checkMethod(tA(), tsu(), "+");
1596     tmp<fvMatrix<Type> > tC(tA.ptr());
1597     tC().source() -= tsu().mesh().V()*tsu().field();
1598     tsu.clear();
1599     return tC;
1600 }
1601
1602 template<class Type>
1603 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1604 (
1605     const tmp<fvMatrix<Type> >& tA,
1606     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
1607 )
1608 {
1609     checkMethod(tA(), tsu(), "+");
1610     tmp<fvMatrix<Type> > tC(tA.ptr());
1611     tC().source() -= tsu().mesh().V()*tsu().internalField();
1612     tsu.clear();
1613     return tC;
1614 }
1615
1616 template<class Type>
1617 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1618 (
1619     const DimensionedField<Type, volMesh>& su,
1620     const fvMatrix<Type>& A
1621 )
1622 {
1623     checkMethod(A, su, "+");
1624     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1625     tC().source() -= su.mesh().V()*su.field();
1626     return tC;
1627 }
1628
1629 template<class Type>
1630 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1631 (
1632     const tmp<DimensionedField<Type, volMesh> >& tsu,
1633     const fvMatrix<Type>& A
1634 )
1635 {
1636     checkMethod(A, tsu(), "+");
1637     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1638     tC().source() -= tsu().mesh().V()*tsu().field();
1639     tsu.clear();
1640     return tC;
1641 }
1642
1643 template<class Type>
1644 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1645 (
1646     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu,
1647     const fvMatrix<Type>& A
1648 )
1649 {
1650     checkMethod(A, tsu(), "+");
1651     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1652     tC().source() -= tsu().mesh().V()*tsu().internalField();
1653     tsu.clear();
1654     return tC;
1655 }
1656
1657 template<class Type>
1658 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1659 (
1660     const DimensionedField<Type, volMesh>& su,
1661     const tmp<fvMatrix<Type> >& tA
1662 )
1663 {
1664     checkMethod(tA(), su, "+");
1665     tmp<fvMatrix<Type> > tC(tA.ptr());
1666     tC().source() -= su.mesh().V()*su.field();
1667     return tC;
1668 }
1669
1670 template<class Type>
1671 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1672 (
1673     const tmp<DimensionedField<Type, volMesh> >& tsu,
1674     const tmp<fvMatrix<Type> >& tA
1675 )
1676 {
1677     checkMethod(tA(), tsu(), "+");
1678     tmp<fvMatrix<Type> > tC(tA.ptr());
1679     tC().source() -= tsu().mesh().V()*tsu().field();
1680     tsu.clear();
1681     return tC;
1682 }
1683
1684 template<class Type>
1685 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1686 (
1687     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu,
1688     const tmp<fvMatrix<Type> >& tA
1689 )
1690 {
1691     checkMethod(tA(), tsu(), "+");
1692     tmp<fvMatrix<Type> > tC(tA.ptr());
1693     tC().source() -= tsu().mesh().V()*tsu().internalField();
1694     tsu.clear();
1695     return tC;
1696 }
1697
1698
1699 template<class Type>
1700 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1701 (
1702     const fvMatrix<Type>& A,
1703     const fvMatrix<Type>& B
1704 )
1705 {
1706     checkMethod(A, B, "-");
1707     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1708     tC() -= B;
1709     return tC;
1710 }
1711
1712 template<class Type>
1713 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1714 (
1715     const tmp<fvMatrix<Type> >& tA,
1716     const fvMatrix<Type>& B
1717 )
1718 {
1719     checkMethod(tA(), B, "-");
1720     tmp<fvMatrix<Type> > tC(tA.ptr());
1721     tC() -= B;
1722     return tC;
1723 }
1724
1725 template<class Type>
1726 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1727 (
1728     const fvMatrix<Type>& A,
1729     const tmp<fvMatrix<Type> >& tB
1730 )
1731 {
1732     checkMethod(A, tB(), "-");
1733     tmp<fvMatrix<Type> > tC(tB.ptr());
1734     tC() -= A;
1735     tC().negate();
1736     return tC;
1737 }
1738
1739 template<class Type>
1740 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1741 (
1742     const tmp<fvMatrix<Type> >& tA,
1743     const tmp<fvMatrix<Type> >& tB
1744 )
1745 {
1746     checkMethod(tA(), tB(), "-");
1747     tmp<fvMatrix<Type> > tC(tA.ptr());
1748     tC() -= tB();
1749     tB.clear();
1750     return tC;
1751 }
1752
1753 template<class Type>
1754 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1755 (
1756     const fvMatrix<Type>& A,
1757     const DimensionedField<Type, volMesh>& su
1758 )
1759 {
1760     checkMethod(A, su, "-");
1761     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1762     tC().source() += su.mesh().V()*su.field();
1763     return tC;
1764 }
1765
1766 template<class Type>
1767 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1768 (
1769     const fvMatrix<Type>& A,
1770     const tmp<DimensionedField<Type, volMesh> >& tsu
1771 )
1772 {
1773     checkMethod(A, tsu(), "-");
1774     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1775     tC().source() += tsu().mesh().V()*tsu().field();
1776     tsu.clear();
1777     return tC;
1778 }
1779
1780 template<class Type>
1781 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1782 (
1783     const fvMatrix<Type>& A,
1784     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
1785 )
1786 {
1787     checkMethod(A, tsu(), "-");
1788     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1789     tC().source() += tsu().mesh().V()*tsu().internalField();
1790     tsu.clear();
1791     return tC;
1792 }
1793
1794 template<class Type>
1795 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1796 (
1797     const tmp<fvMatrix<Type> >& tA,
1798     const DimensionedField<Type, volMesh>& su
1799 )
1800 {
1801     checkMethod(tA(), su, "-");
1802     tmp<fvMatrix<Type> > tC(tA.ptr());
1803     tC().source() += su.mesh().V()*su.field();
1804     return tC;
1805 }
1806
1807 template<class Type>
1808 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1809 (
1810     const tmp<fvMatrix<Type> >& tA,
1811     const tmp<DimensionedField<Type, volMesh> >& tsu
1812 )
1813 {
1814     checkMethod(tA(), tsu(), "-");
1815     tmp<fvMatrix<Type> > tC(tA.ptr());
1816     tC().source() += tsu().mesh().V()*tsu().field();
1817     tsu.clear();
1818     return tC;
1819 }
1820
1821 template<class Type>
1822 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1823 (
1824     const tmp<fvMatrix<Type> >& tA,
1825     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu
1826 )
1827 {
1828     checkMethod(tA(), tsu(), "-");
1829     tmp<fvMatrix<Type> > tC(tA.ptr());
1830     tC().source() += tsu().mesh().V()*tsu().internalField();
1831     tsu.clear();
1832     return tC;
1833 }
1834
1835 template<class Type>
1836 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1837 (
1838     const DimensionedField<Type, volMesh>& su,
1839     const fvMatrix<Type>& A
1840 )
1841 {
1842     checkMethod(A, su, "-");
1843     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1844     tC().negate();
1845     tC().source() -= su.mesh().V()*su.field();
1846     return tC;
1847 }
1848
1849 template<class Type>
1850 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1851 (
1852     const tmp<DimensionedField<Type, volMesh> >& tsu,
1853     const fvMatrix<Type>& A
1854 )
1855 {
1856     checkMethod(A, tsu(), "-");
1857     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1858     tC().negate();
1859     tC().source() -= tsu().mesh().V()*tsu().field();
1860     tsu.clear();
1861     return tC;
1862 }
1863
1864 template<class Type>
1865 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1866 (
1867     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu,
1868     const fvMatrix<Type>& A
1869 )
1870 {
1871     checkMethod(A, tsu(), "-");
1872     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1873     tC().negate();
1874     tC().source() -= tsu().mesh().V()*tsu().internalField();
1875     tsu.clear();
1876     return tC;
1877 }
1878
1879 template<class Type>
1880 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1881 (
1882     const DimensionedField<Type, volMesh>& su,
1883     const tmp<fvMatrix<Type> >& tA
1884 )
1885 {
1886     checkMethod(tA(), su, "-");
1887     tmp<fvMatrix<Type> > tC(tA.ptr());
1888     tC().negate();
1889     tC().source() -= su.mesh().V()*su.field();
1890     return tC;
1891 }
1892
1893 template<class Type>
1894 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1895 (
1896     const tmp<DimensionedField<Type, volMesh> >& tsu,
1897     const tmp<fvMatrix<Type> >& tA
1898 )
1899 {
1900     checkMethod(tA(), tsu(), "-");
1901     tmp<fvMatrix<Type> > tC(tA.ptr());
1902     tC().negate();
1903     tC().source() -= tsu().mesh().V()*tsu().field();
1904     tsu.clear();
1905     return tC;
1906 }
1907
1908 template<class Type>
1909 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1910 (
1911     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tsu,
1912     const tmp<fvMatrix<Type> >& tA
1913 )
1914 {
1915     checkMethod(tA(), tsu(), "-");
1916     tmp<fvMatrix<Type> > tC(tA.ptr());
1917     tC().negate();
1918     tC().source() -= tsu().mesh().V()*tsu().internalField();
1919     tsu.clear();
1920     return tC;
1921 }
1922
1923 template<class Type>
1924 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1925 (
1926     const fvMatrix<Type>& A,
1927     const dimensioned<Type>& su
1928 )
1929 {
1930     checkMethod(A, su, "+");
1931     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1932     tC().source() -= su.value()*A.psi().mesh().V();
1933     return tC;
1934 }
1935
1936 template<class Type>
1937 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1938 (
1939     const tmp<fvMatrix<Type> >& tA,
1940     const dimensioned<Type>& su
1941 )
1942 {
1943     checkMethod(tA(), su, "+");
1944     tmp<fvMatrix<Type> > tC(tA.ptr());
1945     tC().source() -= su.value()*tC().psi().mesh().V();
1946     return tC;
1947 }
1948
1949 template<class Type>
1950 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1951 (
1952     const dimensioned<Type>& su,
1953     const fvMatrix<Type>& A
1954 )
1955 {
1956     checkMethod(A, su, "+");
1957     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1958     tC().source() -= su.value()*A.psi().mesh().V();
1959     return tC;
1960 }
1961
1962 template<class Type>
1963 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator+
1964 (
1965     const dimensioned<Type>& su,
1966     const tmp<fvMatrix<Type> >& tA
1967 )
1968 {
1969     checkMethod(tA(), su, "+");
1970     tmp<fvMatrix<Type> > tC(tA.ptr());
1971     tC().source() -= su.value()*tC().psi().mesh().V();
1972     return tC;
1973 }
1974
1975 template<class Type>
1976 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1977 (
1978     const fvMatrix<Type>& A,
1979     const dimensioned<Type>& su
1980 )
1981 {
1982     checkMethod(A, su, "-");
1983     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1984     tC().source() += su.value()*tC().psi().mesh().V();
1985     return tC;
1986 }
1987
1988 template<class Type>
1989 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1990 (
1991     const tmp<fvMatrix<Type> >& tA,
1992     const dimensioned<Type>& su
1993 )
1994 {
1995     checkMethod(tA(), su, "-");
1996     tmp<fvMatrix<Type> > tC(tA.ptr());
1997     tC().source() += su.value()*tC().psi().mesh().V();
1998     return tC;
1999 }
2000
2001 template<class Type>
2002 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
2003 (
2004     const dimensioned<Type>& su,
2005     const fvMatrix<Type>& A
2006 )
2007 {
2008     checkMethod(A, su, "-");
2009     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
2010     tC().negate();
2011     tC().source() -= su.value()*A.psi().mesh().V();
2012     return tC;
2013 }
2014
2015 template<class Type>
2016 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
2017 (
2018     const dimensioned<Type>& su,
2019     const tmp<fvMatrix<Type> >& tA
2020 )
2021 {
2022     checkMethod(tA(), su, "-");
2023     tmp<fvMatrix<Type> > tC(tA.ptr());
2024     tC().negate();
2025     tC().source() -= su.value()*tC().psi().mesh().V();
2026     return tC;
2027 }
2028
2029
2030 template<class Type>
2031 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator*
2032 (
2033     const DimensionedField<scalar, volMesh>& dsf,
2034     const fvMatrix<Type>& A
2035 )
2036 {
2037     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
2038     tC() *= dsf;
2039     return tC;
2040 }
2041
2042 template<class Type>
2043 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator*
2044 (
2045     const tmp< DimensionedField<scalar, volMesh> >& tdsf,
2046     const fvMatrix<Type>& A
2047 )
2048 {
2049     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
2050     tC() *= tdsf;
2051     return tC;
2052 }
2053
2054 template<class Type>
2055 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator*
2056 (
2057     const tmp<volScalarField>& tvsf,
2058     const fvMatrix<Type>& A
2059 )
2060 {
2061     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
2062     tC() *= tvsf;
2063     return tC;
2064 }
2065
2066 template<class Type>
2067 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator*
2068 (
2069     const DimensionedField<scalar, volMesh>& dsf,
2070     const tmp<fvMatrix<Type> >& tA
2071 )
2072 {
2073     tmp<fvMatrix<Type> > tC(tA.ptr());
2074     tC() *= dsf;
2075     return tC;
2076 }
2077
2078 template<class Type>
2079 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator*
2080 (
2081     const tmp<DimensionedField<scalar, volMesh> >& tdsf,
2082     const tmp<fvMatrix<Type> >& tA
2083 )
2084 {
2085     tmp<fvMatrix<Type> > tC(tA.ptr());
2086     tC() *= tdsf;
2087     return tC;
2088 }
2089
2090 template<class Type>
2091 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator*
2092 (
2093     const tmp<volScalarField>& tvsf,
2094     const tmp<fvMatrix<Type> >& tA
2095 )
2096 {
2097     tmp<fvMatrix<Type> > tC(tA.ptr());
2098     tC() *= tvsf;
2099     return tC;
2100 }
2101
2102 template<class Type>
2103 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator*
2104 (
2105     const dimensioned<scalar>& ds,
2106     const fvMatrix<Type>& A
2107 )
2108 {
2109     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
2110     tC() *= ds;
2111     return tC;
2112 }
2113
2114 template<class Type>
2115 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator*
2116 (
2117     const dimensioned<scalar>& ds,
2118     const tmp<fvMatrix<Type> >& tA
2119 )
2120 {
2121     tmp<fvMatrix<Type> > tC(tA.ptr());
2122     tC() *= ds;
2123     return tC;
2124 }
2125
2126
2127 template<class Type>
2128 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
2129 Foam::operator&
2130 (
2131     const fvMatrix<Type>& M,
2132     const DimensionedField<Type, volMesh>& psi
2133 )
2134 {
2135     tmp<GeometricField<Type, fvPatchField, volMesh> > tMphi
2136     (
2137         new GeometricField<Type, fvPatchField, volMesh>
2138         (
2139             IOobject
2140             (
2141                 "M&" + psi.name(),
2142                 psi.instance(),
2143                 psi.mesh(),
2144                 IOobject::NO_READ,
2145                 IOobject::NO_WRITE
2146             ),
2147             psi.mesh(),
2148             M.dimensions()/dimVol,
2149             zeroGradientFvPatchScalarField::typeName
2150         )
2151     );
2152     GeometricField<Type, fvPatchField, volMesh>& Mphi = tMphi();
2153
2154     // Loop over field components
2155     for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
2156     {
2157         scalarField psiCmpt = psi.field().component(cmpt);
2158         scalarField boundaryDiagCmpt(M.diag());
2159         M.addBoundaryDiag(boundaryDiagCmpt, cmpt);
2160         Mphi.internalField().replace(cmpt, -boundaryDiagCmpt*psiCmpt);
2161     }
2162
2163     Mphi.internalField() += M.lduMatrix::H(psi.field()) + M.source();
2164     M.addBoundarySource(Mphi.internalField());
2165
2166     Mphi.internalField() /= -psi.mesh().V();
2167     Mphi.correctBoundaryConditions();
2168
2169     return tMphi;
2170 }
2171
2172 template<class Type>
2173 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
2174 Foam::operator&
2175 (
2176     const fvMatrix<Type>& M,
2177     const tmp<DimensionedField<Type, volMesh> >& tpsi
2178 )
2179 {
2180     tmp<GeometricField<Type, fvPatchField, volMesh> > tMpsi = M & tpsi();
2181     tpsi.clear();
2182     return tMpsi;
2183 }
2184
2185 template<class Type>
2186 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
2187 Foam::operator&
2188 (
2189     const fvMatrix<Type>& M,
2190     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tpsi
2191 )
2192 {
2193     tmp<GeometricField<Type, fvPatchField, volMesh> > tMpsi = M & tpsi();
2194     tpsi.clear();
2195     return tMpsi;
2196 }
2197
2198 template<class Type>
2199 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
2200 Foam::operator&
2201 (
2202     const tmp<fvMatrix<Type> >& tM,
2203     const DimensionedField<Type, volMesh>& psi
2204 )
2205 {
2206     tmp<GeometricField<Type, fvPatchField, volMesh> > tMpsi = tM() & psi;
2207     tM.clear();
2208     return tMpsi;
2209 }
2210
2211 template<class Type>
2212 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
2213 Foam::operator&
2214 (
2215     const tmp<fvMatrix<Type> >& tM,
2216     const tmp<DimensionedField<Type, volMesh> >& tpsi
2217 )
2218 {
2219     tmp<GeometricField<Type, fvPatchField, volMesh> > tMpsi = tM() & tpsi();
2220     tM.clear();
2221     tpsi.clear();
2222     return tMpsi;
2223 }
2224
2225 template<class Type>
2226 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
2227 Foam::operator&
2228 (
2229     const tmp<fvMatrix<Type> >& tM,
2230     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tpsi
2231 )
2232 {
2233     tmp<GeometricField<Type, fvPatchField, volMesh> > tMpsi = tM() & tpsi();
2234     tM.clear();
2235     tpsi.clear();
2236     return tMpsi;
2237 }
2238
2239
2240 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
2241
2242 template<class Type>
2243 Foam::Ostream& Foam::operator<<(Ostream& os, const fvMatrix<Type>& fvm)
2244 {
2245     os  << static_cast<const lduMatrix&>(fvm) << nl
2246         << fvm.dimensions_ << nl
2247         << fvm.source_ << nl
2248         << fvm.internalCoeffs_ << nl
2249         << fvm.boundaryCoeffs_ << endl;
2250
2251     os.check("Ostream& operator<<(Ostream&, fvMatrix<Type>&");
2252
2253     return os;
2254 }
2255
2256
2257 // * * * * * * * * * * * * * * * * Solvers * * * * * * * * * * * * * * * * * //
2258
2259 #include "fvMatrixSolve.C"
2260
2261 // ************************************************************************* //