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