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