fvMatrix.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
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"
30 #include "coupledFvPatchFields.H"
31 #include "UIndirectList.H"
32 #include "UCompactListList.H"
33 #include "fvmDdt.H"
34 
35 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
36 
37 template<class Type>
38 template<class Type2>
40 (
41  const labelUList& addr,
42  const Field<Type2>& pf,
43  Field<Type2>& intf
44 ) const
45 {
46  if (addr.size() != pf.size())
47  {
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 labelUList& 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 labelUList& addr,
79  const Field<Type2>& pf,
80  Field<Type2>& intf
81 ) const
82 {
83  if (addr.size() != pf.size())
84  {
86  << "sizes of addressing and field are different"
87  << abort(FatalError);
88  }
89 
90  forAll(addr, facei)
91  {
92  intf[addr[facei]] -= pf[facei];
93  }
94 }
95 
96 
97 template<class Type>
98 template<class Type2>
100 (
101  const labelUList& addr,
102  const tmp<Field<Type2>>& tpf,
103  Field<Type2>& intf
104 ) const
105 {
106  subtractFromInternalField(addr, tpf(), intf);
107  tpf.clear();
108 }
109 
110 
111 template<class Type>
113 (
114  scalarField& diag,
115  const direction solveCmpt
116 ) const
117 {
118  forAll(internalCoeffs_, patchi)
119  {
120  addToInternalField
121  (
122  lduAddr().patchAddr(patchi),
123  internalCoeffs_[patchi].component(solveCmpt),
124  diag
125  );
126  }
127 }
128 
129 
130 template<class Type>
132 {
133  forAll(internalCoeffs_, patchi)
134  {
135  addToInternalField
136  (
137  lduAddr().patchAddr(patchi),
138  cmptAv(internalCoeffs_[patchi]),
139  diag
140  );
141  }
142 }
143 
144 
145 template<class Type>
147 (
148  Field<Type>& source,
149  const bool couples
150 ) const
151 {
152  forAll(psi_.boundaryField(), patchi)
153  {
154  const fvPatchField<Type>& ptf = psi_.boundaryField()[patchi];
155  const Field<Type>& pbc = boundaryCoeffs_[patchi];
156 
157  if (!ptf.coupled())
158  {
159  addToInternalField(lduAddr().patchAddr(patchi), pbc, source);
160  }
161  else if (couples)
162  {
163  const tmp<Field<Type>> tpnf = ptf.patchNeighbourField();
164  const Field<Type>& pnf = tpnf();
165 
166  const labelUList& addr = lduAddr().patchAddr(patchi);
167 
168  forAll(addr, facei)
169  {
170  source[addr[facei]] += cmptMultiply(pbc[facei], pnf[facei]);
171  }
172  }
173  }
174 }
175 
176 
177 template<class Type>
179 (
180  const label celli,
181  const Type& value
182 )
183 {
184  const fvMesh& mesh = psi_.mesh();
185 
186  const cellList& cells = mesh.cells();
187  const labelUList& own = mesh.owner();
188  const labelUList& nei = mesh.neighbour();
189 
190  scalarField& Diag = diag();
191 
192  Field<Type>& psi =
193  const_cast<VolField<Type>&>(psi_)
195 
196  psi[celli] = value;
197  source_[celli] = value*Diag[celli];
198 
199  if (symmetric() || asymmetric())
200  {
201  const cell& c = cells[celli];
202 
203  forAll(c, j)
204  {
205  const label facei = c[j];
206 
207  if (mesh.isInternalFace(facei))
208  {
209  if (symmetric())
210  {
211  if (celli == own[facei])
212  {
213  source_[nei[facei]] -= upper()[facei]*value;
214  }
215  else
216  {
217  source_[own[facei]] -= upper()[facei]*value;
218  }
219 
220  upper()[facei] = 0.0;
221  }
222  else
223  {
224  if (celli == own[facei])
225  {
226  source_[nei[facei]] -= lower()[facei]*value;
227  }
228  else
229  {
230  source_[own[facei]] -= upper()[facei]*value;
231  }
232 
233  upper()[facei] = 0.0;
234  lower()[facei] = 0.0;
235  }
236  }
237  else
238  {
239  const label bFacei = facei - mesh.nInternalFaces();
240 
241  const labelUList patches =
242  mesh.polyBFacePatches()[bFacei];
243  const labelUList patchFaces =
244  mesh.polyBFacePatchFaces()[bFacei];
245 
246  forAll(patches, i)
247  {
248  internalCoeffs_[patches[i]][patchFaces[i]] = Zero;
249  boundaryCoeffs_[patches[i]][patchFaces[i]] = Zero;
250  }
251  }
252  }
253  }
254 }
255 
256 
257 template<class Type>
259 (
260  const label celli,
261  const Type& value,
262  const scalar fraction,
263  const scalarField& ddtDiag
264 )
265 {
266  Field<Type>& psi =
267  const_cast<VolField<Type>&>(psi_)
269 
270  psi[celli] = (1 - fraction)*psi[celli] + fraction*value;
271 
272  const scalar coeff =
273  fraction/(1 - fraction)
274  *(diag()[celli] - ddtDiag[celli]);
275 
276  diag()[celli] += coeff;
277  source()[celli] += coeff*value;
278 }
279 
280 
281 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
282 
283 template<class Type>
285 (
286  const VolField<Type>& psi,
287  const dimensionSet& ds
288 )
289 :
290  lduMatrix(psi.mesh()),
291  psi_(psi),
292  dimensions_(ds),
293  source_(psi.size(), Zero),
294  internalCoeffs_(psi.mesh().boundary().size()),
295  boundaryCoeffs_(psi.mesh().boundary().size()),
296  faceFluxCorrectionPtr_(nullptr)
297 {
298  if (debug)
299  {
301  << "Constructing fvMatrix<Type> for field " << psi_.name() << endl;
302  }
303 
304  // Initialise coupling coefficients
305  forAll(psi.mesh().boundary(), patchi)
306  {
307  internalCoeffs_.set
308  (
309  patchi,
310  new Field<Type>
311  (
312  psi.mesh().boundary()[patchi].size(),
313  Zero
314  )
315  );
316 
317  boundaryCoeffs_.set
318  (
319  patchi,
320  new Field<Type>
321  (
322  psi.mesh().boundary()[patchi].size(),
323  Zero
324  )
325  );
326  }
327 
328  // Update the boundary coefficients of psi without changing its event No.
329  VolField<Type>& psiRef =
330  const_cast<VolField<Type>&>(psi_);
331 
332  label currentStatePsi = psiRef.eventNo();
333  psiRef.boundaryFieldRef().updateCoeffs();
334  psiRef.eventNo() = currentStatePsi;
335 }
336 
337 
338 template<class Type>
340 :
341  tmp<fvMatrix<Type>>::refCount(),
342  lduMatrix(fvm),
343  psi_(fvm.psi_),
344  dimensions_(fvm.dimensions_),
345  source_(fvm.source_),
346  internalCoeffs_(fvm.internalCoeffs_),
347  boundaryCoeffs_(fvm.boundaryCoeffs_),
348  faceFluxCorrectionPtr_(nullptr)
349 {
350  if (debug)
351  {
353  << "Copying fvMatrix<Type> for field " << psi_.name() << endl;
354  }
355 
356  if (fvm.faceFluxCorrectionPtr_)
357  {
358  faceFluxCorrectionPtr_ = new
360  (
361  *(fvm.faceFluxCorrectionPtr_)
362  );
363  }
364 }
365 
366 
367 template<class Type>
369 :
370  lduMatrix
371  (
372  const_cast<fvMatrix<Type>&>(tfvm()),
373  tfvm.isTmp()
374  ),
375  psi_(tfvm().psi_),
376  dimensions_(tfvm().dimensions_),
377  source_
378  (
379  const_cast<fvMatrix<Type>&>(tfvm()).source_,
380  tfvm.isTmp()
381  ),
382  internalCoeffs_
383  (
384  const_cast<fvMatrix<Type>&>(tfvm()).internalCoeffs_,
385  tfvm.isTmp()
386  ),
387  boundaryCoeffs_
388  (
389  const_cast<fvMatrix<Type>&>(tfvm()).boundaryCoeffs_,
390  tfvm.isTmp()
391  ),
392  faceFluxCorrectionPtr_(nullptr)
393 {
394  if (debug)
395  {
397  << "Copying fvMatrix<Type> for field " << psi_.name() << endl;
398  }
399 
400  if (tfvm().faceFluxCorrectionPtr_)
401  {
402  if (tfvm.isTmp())
403  {
404  faceFluxCorrectionPtr_ = tfvm().faceFluxCorrectionPtr_;
405  tfvm().faceFluxCorrectionPtr_ = nullptr;
406  }
407  else
408  {
409  faceFluxCorrectionPtr_ = new
411  (
412  *(tfvm().faceFluxCorrectionPtr_)
413  );
414  }
415  }
416 
417  tfvm.clear();
418 }
419 
420 
421 template<class Type>
423 (
424  const VolField<Type>& psi,
425  Istream& is
426 )
427 :
428  lduMatrix(psi.mesh()),
429  psi_(psi),
430  dimensions_(is),
431  source_(is),
432  internalCoeffs_(psi.mesh().boundary().size()),
433  boundaryCoeffs_(psi.mesh().boundary().size()),
434  faceFluxCorrectionPtr_(nullptr)
435 {
436  if (debug)
437  {
439  << "Constructing fvMatrix<Type> for field " << psi_.name() << endl;
440  }
441 
442  // Initialise coupling coefficients
443  forAll(psi.mesh().boundary(), patchi)
444  {
445  internalCoeffs_.set
446  (
447  patchi,
448  new Field<Type>
449  (
450  psi.mesh().boundary()[patchi].size(),
451  Zero
452  )
453  );
454 
455  boundaryCoeffs_.set
456  (
457  patchi,
458  new Field<Type>
459  (
460  psi.mesh().boundary()[patchi].size(),
461  Zero
462  )
463  );
464  }
465 
466 }
467 
468 
469 template<class Type>
471 {
472  return tmp<fvMatrix<Type>>
473  (
474  new fvMatrix<Type>(*this)
475  );
476 }
477 
478 
479 // * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * * * //
480 
481 template<class Type>
483 {
484  if (debug)
485  {
487  << "Destroying fvMatrix<Type> for field " << psi_.name() << endl;
488  }
489 
490  if (faceFluxCorrectionPtr_)
491  {
492  delete faceFluxCorrectionPtr_;
493  }
494 }
495 
496 
497 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
498 
499 template<class Type>
500 template<template<class> class ListType>
502 (
503  const labelUList& cellLabels,
504  const ListType<Type>& values
505 )
506 {
507  // Fix the values
508  forAll(cellLabels, i)
509  {
510  setValue(cellLabels[i], values[i]);
511  }
512 }
513 
514 
515 template<class Type>
516 template<template<class> class ListType>
518 (
519  const labelUList& cellLabels,
520  const ListType<Type>& values,
521  const scalarList& fractions,
522  const bool hasDdt
523 )
524 {
525  // Get the proportion of the diagonal associated with iterative update
526  scalarField ddtDiag(diag().size(), 0);
527  const scalar alpha = relaxationFactor();
528  if (alpha > 0)
529  {
530  ddtDiag += (1 - alpha)*diag();
531  }
532  if (hasDdt)
533  {
534  const fvMatrix<Type> ddtEqn(fvm::ddt(psi_));
535  if (ddtEqn.hasDiag())
536  {
537  ddtDiag += ddtEqn.diag();
538  }
539  }
540 
541  forAll(cellLabels, i)
542  {
543  if (- rootSmall < fractions[i] && fractions[i] < rootSmall)
544  {
545  // Do nothing
546  }
547  else if (1 - rootSmall < fractions[i] && fractions[i] < 1 + rootSmall)
548  {
549  // Fix the values
550  setValue(cellLabels[i], values[i]);
551  }
552  else
553  {
554  // Fractionally fix the values
555  setValue(cellLabels[i], values[i], fractions[i], ddtDiag);
556  }
557  }
558 }
559 
560 
561 template<class Type>
563 (
564  const label celli,
565  const Type& value,
566  const bool forceReference
567 )
568 {
569  if ((forceReference || psi_.needReference()) && celli >= 0)
570  {
571  source()[celli] += diag()[celli]*value;
572  diag()[celli] += diag()[celli];
573  }
574 }
575 
576 
577 template<class Type>
579 {
580  if
581  (
582  psi_.mesh().data::template lookupOrDefault<bool>
583  ("finalIteration", false)
584  && psi_.mesh().solution().relaxEquation(psi_.name() + "Final")
585  )
586  {
587  return psi_.mesh().solution().equationRelaxationFactor
588  (
589  psi_.name() + "Final"
590  );
591  }
592  else if (psi_.mesh().solution().relaxEquation(psi_.name()))
593  {
594  return psi_.mesh().solution().equationRelaxationFactor(psi_.name());
595  }
596  else
597  {
598  return 0;
599  }
600 }
601 
602 
603 template<class Type>
605 {
606  if (alpha <= 0)
607  {
608  return;
609  }
610 
611  if (debug)
612  {
614  << "Relaxing " << psi_.name() << " by " << alpha << endl;
615  }
616 
617  Field<Type>& S = source();
618  scalarField& D = diag();
619 
620  // Store the current unrelaxed diagonal for use in updating the source
621  scalarField D0(D);
622 
623  // Calculate the sum-mag off-diagonal from the interior faces
624  scalarField sumOff(D.size(), 0.0);
625  sumMagOffDiag(sumOff);
626 
627  // Handle the boundary contributions to the diagonal
628  forAll(psi_.boundaryField(), patchi)
629  {
630  const fvPatchField<Type>& ptf = psi_.boundaryField()[patchi];
631 
632  if (ptf.size())
633  {
634  const labelUList& pa = lduAddr().patchAddr(patchi);
635  Field<Type>& iCoeffs = internalCoeffs_[patchi];
636 
637  if (ptf.coupled())
638  {
639  const Field<Type>& pCoeffs = boundaryCoeffs_[patchi];
640 
641  // For coupled boundaries add the diagonal and
642  // off-diagonal contributions
643  forAll(pa, face)
644  {
645  D[pa[face]] += component(iCoeffs[face], 0);
646  sumOff[pa[face]] += mag(component(pCoeffs[face], 0));
647  }
648  }
649  else
650  {
651  // For non-coupled boundaries add the maximum magnitude diagonal
652  // contribution to ensure stability
653  forAll(pa, face)
654  {
655  D[pa[face]] += cmptMax(cmptMag(iCoeffs[face]));
656  }
657  }
658  }
659  }
660 
661 
662  if (debug)
663  {
664  // Calculate amount of non-dominance.
665  label nNon = 0;
666  scalar maxNon = 0.0;
667  scalar sumNon = 0.0;
668  forAll(D, celli)
669  {
670  scalar d = (sumOff[celli] - D[celli])/mag(D[celli]);
671 
672  if (d > 0)
673  {
674  nNon++;
675  maxNon = max(maxNon, d);
676  sumNon += d;
677  }
678  }
679 
680  reduce(nNon, sumOp<label>(), UPstream::msgType(), psi_.mesh().comm());
681  reduce
682  (
683  maxNon,
684  maxOp<scalar>(),
686  psi_.mesh().comm()
687  );
688  reduce
689  (
690  sumNon,
691  sumOp<scalar>(),
693  psi_.mesh().comm()
694  );
695  sumNon /= returnReduce
696  (
697  D.size(),
698  sumOp<label>(),
700  psi_.mesh().comm()
701  );
702 
704  << "Matrix dominance test for " << psi_.name() << nl
705  << " number of non-dominant cells : " << nNon << nl
706  << " maximum relative non-dominance : " << maxNon << nl
707  << " average relative non-dominance : " << sumNon << nl
708  << endl;
709  }
710 
711 
712  // Ensure the matrix is diagonally dominant...
713  // Assumes that the central coefficient is positive and ensures it is
714  forAll(D, celli)
715  {
716  D[celli] = max(mag(D[celli]), sumOff[celli]);
717  }
718 
719  // ... then relax
720  D /= alpha;
721 
722  // Now remove the diagonal contribution from coupled boundaries
723  forAll(psi_.boundaryField(), patchi)
724  {
725  const fvPatchField<Type>& ptf = psi_.boundaryField()[patchi];
726 
727  if (ptf.size())
728  {
729  const labelUList& pa = lduAddr().patchAddr(patchi);
730  Field<Type>& iCoeffs = internalCoeffs_[patchi];
731 
732  if (ptf.coupled())
733  {
734  forAll(pa, face)
735  {
736  D[pa[face]] -= component(iCoeffs[face], 0);
737  }
738  }
739  else
740  {
741  forAll(pa, face)
742  {
743  D[pa[face]] -= cmptMin(iCoeffs[face]);
744  }
745  }
746  }
747  }
748 
749  // Finally add the relaxation contribution to the source.
750  S += (D - D0)*psi_.primitiveField();
751 }
752 
753 
754 template<class Type>
756 {
757  relax(relaxationFactor());
758 }
759 
760 
761 template<class Type>
763 (
764  typename VolField<Type>::
765  Boundary& bFields
766 )
767 {
768  forAll(bFields, patchi)
769  {
770  bFields[patchi].manipulateMatrix(*this);
771  }
772 }
773 
774 
775 template<class Type>
777 {
778  tmp<scalarField> tdiag(new scalarField(diag()));
779  addCmptAvBoundaryDiag(tdiag.ref());
780  return tdiag;
781 }
782 
783 
784 template<class Type>
786 {
788 
789  forAll(psi_.boundaryField(), patchi)
790  {
791  const fvPatchField<Type>& ptf = psi_.boundaryField()[patchi];
792 
793  if (!ptf.coupled() && ptf.size())
794  {
795  addToInternalField
796  (
797  lduAddr().patchAddr(patchi),
798  internalCoeffs_[patchi],
799  tdiag.ref()
800  );
801  }
802  }
803 
804  return tdiag;
805 }
806 
807 
808 template<class Type>
810 {
811  tmp<volScalarField> tAphi
812  (
814  (
815  "A(" + psi_.name() + ')',
816  psi_.mesh(),
817  dimensions_/psi_.dimensions()/dimVol,
818  extrapolatedCalculatedFvPatchScalarField::typeName
819  )
820  );
821 
822  tAphi.ref().primitiveFieldRef() = D()/psi_.mesh().V();
823  tAphi.ref().correctBoundaryConditions();
824 
825  return tAphi;
826 }
827 
828 
829 template<class Type>
831 {
833  (
835  (
836  "Su(" +psi_.name() + ')',
837  psi_.mesh(),
838  dimensions_/dimVol,
839  -source()/psi_.mesh().V()
840  )
841  );
842 
843  return tSu;
844 }
845 
846 
847 template<class Type>
849 {
851  (
853  (
854  "Sp(" + psi_.name() + ')',
855  psi_.mesh(),
856  dimensions_/psi_.dimensions()/dimVol,
857  diag()/psi_.mesh().V()
858  )
859  );
860 
861  return tSp;
862 }
863 
864 
865 template<class Type>
868 {
869  tmp<VolField<Type>> tHphi
870  (
872  (
873  "H(" + psi_.name() + ')',
874  psi_.mesh(),
875  dimensions_/dimVol,
876  extrapolatedCalculatedFvPatchScalarField::typeName
877  )
878  );
879  VolField<Type>& Hphi = tHphi.ref();
880 
881  // Loop over field components
882  for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
883  {
884  scalarField psiCmpt(psi_.primitiveField().component(cmpt));
885 
886  scalarField boundaryDiagCmpt(psi_.size(), 0.0);
887  addBoundaryDiag(boundaryDiagCmpt, cmpt);
888  boundaryDiagCmpt.negate();
889  addCmptAvBoundaryDiag(boundaryDiagCmpt);
890 
891  Hphi.primitiveFieldRef().replace(cmpt, boundaryDiagCmpt*psiCmpt);
892  }
893 
894  Hphi.primitiveFieldRef() += lduMatrix::H(psi_.primitiveField()) + source_;
895  addBoundarySource(Hphi.primitiveFieldRef());
896 
897  Hphi.primitiveFieldRef() /= psi_.mesh().V();
899 
900  typename Type::labelType validComponents
901  (
902  psi_.mesh().template validComponents<Type>()
903  );
904 
905  for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
906  {
907  if (validComponents[cmpt] == -1)
908  {
909  Hphi.replace
910  (
911  cmpt,
912  dimensionedScalar(Hphi.dimensions(), 0)
913  );
914  }
915  }
916 
917  return tHphi;
918 }
919 
920 
921 template<class Type>
923 {
925  (
927  (
928  "H(1)",
929  psi_.mesh(),
930  dimensions_/(dimVol*psi_.dimensions()),
931  extrapolatedCalculatedFvPatchScalarField::typeName
932  )
933  );
934  volScalarField& H1_ = tH1.ref();
935 
937 
938  forAll(psi_.boundaryField(), patchi)
939  {
940  const fvPatchField<Type>& ptf = psi_.boundaryField()[patchi];
941 
942  if (ptf.coupled() && ptf.size())
943  {
944  addToInternalField
945  (
946  lduAddr().patchAddr(patchi),
947  boundaryCoeffs_[patchi].component(0),
948  H1_
949  );
950  }
951  }
952 
953  H1_.primitiveFieldRef() /= psi_.mesh().V();
955 
956  return tH1;
957 }
958 
959 
960 template<class Type>
963 flux() const
964 {
965  if (!psi_.mesh().schemes().fluxRequired(psi_.name()))
966  {
968  << "flux requested but " << psi_.name()
969  << " not specified in the fluxRequired sub-dictionary"
970  " of fvSchemes."
971  << abort(FatalError);
972  }
973 
974  // construct SurfaceField<Type>
975  tmp<SurfaceField<Type>> tfieldFlux
976  (
978  (
979  "flux(" + psi_.name() + ')',
980  psi_.mesh(),
981  dimensions()
982  )
983  );
984  SurfaceField<Type>& fieldFlux =
985  tfieldFlux.ref();
986 
987  for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
988  {
989  fieldFlux.primitiveFieldRef().replace
990  (
991  cmpt,
992  lduMatrix::faceH(psi_.primitiveField().component(cmpt))
993  );
994  }
995 
996  FieldField<Field, Type> InternalContrib = internalCoeffs_;
997 
998  forAll(InternalContrib, patchi)
999  {
1000  InternalContrib[patchi] =
1001  cmptMultiply
1002  (
1003  InternalContrib[patchi],
1004  psi_.boundaryField()[patchi].patchInternalField()
1005  );
1006  }
1007 
1008  FieldField<Field, Type> NeighbourContrib = boundaryCoeffs_;
1009 
1010  forAll(NeighbourContrib, patchi)
1011  {
1012  if (psi_.boundaryField()[patchi].coupled())
1013  {
1014  NeighbourContrib[patchi] =
1015  cmptMultiply
1016  (
1017  NeighbourContrib[patchi],
1018  psi_.boundaryField()[patchi].patchNeighbourField()
1019  );
1020  }
1021  }
1022 
1023  typename SurfaceField<Type>::
1024  Boundary& ffbf = fieldFlux.boundaryFieldRef();
1025 
1026  forAll(ffbf, patchi)
1027  {
1028  ffbf[patchi] = InternalContrib[patchi] - NeighbourContrib[patchi];
1029  }
1030 
1031  if (faceFluxCorrectionPtr_)
1032  {
1033  fieldFlux += *faceFluxCorrectionPtr_;
1034  }
1035 
1036  return tfieldFlux;
1037 }
1038 
1039 
1040 // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
1041 
1042 template<class Type>
1044 {
1045  if (this == &fvmv)
1046  {
1048  << "attempted assignment to self"
1049  << abort(FatalError);
1050  }
1051 
1052  if (&psi_ != &(fvmv.psi_))
1053  {
1055  << "different fields"
1056  << abort(FatalError);
1057  }
1058 
1059  dimensions_ = fvmv.dimensions_;
1060  lduMatrix::operator=(fvmv);
1061  source_ = fvmv.source_;
1062  internalCoeffs_ = fvmv.internalCoeffs_;
1063  boundaryCoeffs_ = fvmv.boundaryCoeffs_;
1064 
1065  if (faceFluxCorrectionPtr_ && fvmv.faceFluxCorrectionPtr_)
1066  {
1067  *faceFluxCorrectionPtr_ = *fvmv.faceFluxCorrectionPtr_;
1068  }
1069  else if (fvmv.faceFluxCorrectionPtr_)
1070  {
1071  faceFluxCorrectionPtr_ =
1072  new SurfaceField<Type>
1073  (*fvmv.faceFluxCorrectionPtr_);
1074  }
1075 }
1076 
1077 
1078 template<class Type>
1080 {
1081  operator=(tfvmv());
1082  tfvmv.clear();
1083 }
1084 
1085 
1086 template<class Type>
1088 {
1090  source_.negate();
1091  internalCoeffs_.negate();
1092  boundaryCoeffs_.negate();
1093 
1094  if (faceFluxCorrectionPtr_)
1095  {
1096  faceFluxCorrectionPtr_->negate();
1097  }
1098 }
1099 
1100 
1101 template<class Type>
1103 {
1104  checkMethod(*this, fvmv, "+=");
1105 
1106  dimensions_ += fvmv.dimensions_;
1107  lduMatrix::operator+=(fvmv);
1108  source_ += fvmv.source_;
1109  internalCoeffs_ += fvmv.internalCoeffs_;
1110  boundaryCoeffs_ += fvmv.boundaryCoeffs_;
1111 
1112  if (faceFluxCorrectionPtr_ && fvmv.faceFluxCorrectionPtr_)
1113  {
1114  *faceFluxCorrectionPtr_ += *fvmv.faceFluxCorrectionPtr_;
1115  }
1116  else if (fvmv.faceFluxCorrectionPtr_)
1117  {
1118  faceFluxCorrectionPtr_ = new
1120  (
1121  *fvmv.faceFluxCorrectionPtr_
1122  );
1123  }
1124 }
1125 
1126 
1127 template<class Type>
1129 {
1130  operator+=(tfvmv());
1131  tfvmv.clear();
1132 }
1133 
1134 
1135 template<class Type>
1137 {
1138  checkMethod(*this, fvmv, "-=");
1139 
1140  dimensions_ -= fvmv.dimensions_;
1141  lduMatrix::operator-=(fvmv);
1142  source_ -= fvmv.source_;
1143  internalCoeffs_ -= fvmv.internalCoeffs_;
1144  boundaryCoeffs_ -= fvmv.boundaryCoeffs_;
1145 
1146  if (faceFluxCorrectionPtr_ && fvmv.faceFluxCorrectionPtr_)
1147  {
1148  *faceFluxCorrectionPtr_ -= *fvmv.faceFluxCorrectionPtr_;
1149  }
1150  else if (fvmv.faceFluxCorrectionPtr_)
1151  {
1152  faceFluxCorrectionPtr_ =
1153  new SurfaceField<Type>
1154  (-*fvmv.faceFluxCorrectionPtr_);
1155  }
1156 }
1157 
1158 
1159 template<class Type>
1161 {
1162  operator-=(tfvmv());
1163  tfvmv.clear();
1164 }
1165 
1166 
1167 template<class Type>
1169 (
1171 )
1172 {
1173  checkMethod(*this, su, "+=");
1174  source() -= su.mesh().V()*su.field();
1175 }
1176 
1177 
1178 template<class Type>
1180 (
1182 )
1183 {
1184  operator+=(tsu());
1185  tsu.clear();
1186 }
1187 
1188 
1189 template<class Type>
1191 (
1192  const tmp<VolField<Type>>& tsu
1193 )
1194 {
1195  operator+=(tsu());
1196  tsu.clear();
1197 }
1198 
1199 
1200 template<class Type>
1202 (
1204 )
1205 {
1206  checkMethod(*this, su, "-=");
1207  source() += su.mesh().V()*su.field();
1208 }
1209 
1210 
1211 template<class Type>
1213 (
1215 )
1216 {
1217  operator-=(tsu());
1218  tsu.clear();
1219 }
1220 
1221 
1222 template<class Type>
1224 (
1225  const tmp<VolField<Type>>& tsu
1226 )
1227 {
1228  operator-=(tsu());
1229  tsu.clear();
1230 }
1231 
1232 
1233 template<class Type>
1235 (
1236  const dimensioned<Type>& su
1237 )
1238 {
1239  source() -= psi().mesh().V()*su;
1240 }
1241 
1242 
1243 template<class Type>
1245 (
1246  const dimensioned<Type>& su
1247 )
1248 {
1249  source() += psi().mesh().V()*su;
1250 }
1251 
1252 
1253 template<class Type>
1255 (
1256  const zero&
1257 )
1258 {}
1259 
1260 
1261 template<class Type>
1263 (
1264  const zero&
1265 )
1266 {}
1267 
1268 
1269 template<class Type>
1271 (
1272  const volScalarField::Internal& dsf
1273 )
1274 {
1275  dimensions_ *= dsf.dimensions();
1276  lduMatrix::operator*=(dsf.field());
1277  source_ *= dsf.field();
1278 
1279  forAll(boundaryCoeffs_, patchi)
1280  {
1281  scalarField pisf
1282  (
1283  dsf.mesh().boundary()[patchi].patchInternalField(dsf.field())
1284  );
1285 
1286  internalCoeffs_[patchi] *= pisf;
1287  boundaryCoeffs_[patchi] *= pisf;
1288  }
1289 
1290  if (faceFluxCorrectionPtr_)
1291  {
1293  << "cannot scale a matrix containing a faceFluxCorrection"
1294  << abort(FatalError);
1295  }
1296 }
1297 
1298 
1299 template<class Type>
1301 (
1302  const tmp<volScalarField::Internal>& tdsf
1303 )
1304 {
1305  operator*=(tdsf());
1306  tdsf.clear();
1307 }
1308 
1309 
1310 template<class Type>
1312 (
1313  const tmp<volScalarField>& tvsf
1314 )
1315 {
1316  operator*=(tvsf());
1317  tvsf.clear();
1318 }
1319 
1320 
1321 template<class Type>
1323 (
1324  const dimensioned<scalar>& ds
1325 )
1326 {
1327  dimensions_ *= ds.dimensions();
1328  lduMatrix::operator*=(ds.value());
1329  source_ *= ds.value();
1330  internalCoeffs_ *= ds.value();
1331  boundaryCoeffs_ *= ds.value();
1332 
1333  if (faceFluxCorrectionPtr_)
1334  {
1335  *faceFluxCorrectionPtr_ *= ds.value();
1336  }
1337 }
1338 
1339 
1340 template<class Type>
1342 (
1343  const volScalarField::Internal& dsf
1344 )
1345 {
1346  dimensions_ /= dsf.dimensions();
1347  lduMatrix::operator/=(dsf.field());
1348  source_ /= dsf.field();
1349 
1350  forAll(boundaryCoeffs_, patchi)
1351  {
1352  scalarField pisf
1353  (
1354  dsf.mesh().boundary()[patchi].patchInternalField(dsf.field())
1355  );
1356 
1357  internalCoeffs_[patchi] /= pisf;
1358  boundaryCoeffs_[patchi] /= pisf;
1359  }
1360 
1361  if (faceFluxCorrectionPtr_)
1362  {
1364  << "cannot scale a matrix containing a faceFluxCorrection"
1365  << abort(FatalError);
1366  }
1367 }
1368 
1369 
1370 template<class Type>
1372 (
1373  const tmp<volScalarField::Internal>& tdsf
1374 )
1375 {
1376  operator/=(tdsf());
1377  tdsf.clear();
1378 }
1379 
1380 
1381 template<class Type>
1383 (
1384  const tmp<volScalarField>& tvsf
1385 )
1386 {
1387  operator/=(tvsf());
1388  tvsf.clear();
1389 }
1390 
1391 
1392 template<class Type>
1394 (
1395  const dimensioned<scalar>& ds
1396 )
1397 {
1398  dimensions_ /= ds.dimensions();
1399  lduMatrix::operator/=(ds.value());
1400  source_ /= ds.value();
1401  internalCoeffs_ /= ds.value();
1402  boundaryCoeffs_ /= ds.value();
1403 
1404  if (faceFluxCorrectionPtr_)
1405  {
1406  *faceFluxCorrectionPtr_ /= ds.value();
1407  }
1408 }
1409 
1410 
1411 // * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
1412 
1413 template<class Type>
1415 (
1416  const fvMatrix<Type>& fvm1,
1417  const fvMatrix<Type>& fvm2,
1418  const char* op
1419 )
1420 {
1421  if (&fvm1.psi() != &fvm2.psi())
1422  {
1424  << "incompatible fields for operation "
1425  << endl << " "
1426  << "[" << fvm1.psi().name() << "] "
1427  << op
1428  << " [" << fvm2.psi().name() << "]"
1429  << abort(FatalError);
1430  }
1431 
1432  if (dimensionSet::debug && fvm1.dimensions() != fvm2.dimensions())
1433  {
1435  << "incompatible dimensions for operation "
1436  << endl << " "
1437  << "[" << fvm1.psi().name() << fvm1.dimensions()/dimVolume << " ] "
1438  << op
1439  << " [" << fvm2.psi().name() << fvm2.dimensions()/dimVolume << " ]"
1440  << abort(FatalError);
1441  }
1442 }
1443 
1444 
1445 template<class Type>
1447 (
1448  const fvMatrix<Type>& fvm,
1450  const char* op
1451 )
1452 {
1453  if (dimensionSet::debug && fvm.dimensions()/dimVolume != df.dimensions())
1454  {
1456  << endl << " "
1457  << "[" << fvm.psi().name() << fvm.dimensions()/dimVolume << " ] "
1458  << op
1459  << " [" << df.name() << df.dimensions() << " ]"
1460  << abort(FatalError);
1461  }
1462 }
1463 
1464 
1465 template<class Type>
1467 (
1468  const fvMatrix<Type>& fvm,
1469  const dimensioned<Type>& dt,
1470  const char* op
1471 )
1472 {
1473  if (dimensionSet::debug && fvm.dimensions()/dimVolume != dt.dimensions())
1474  {
1476  << "incompatible dimensions for operation "
1477  << endl << " "
1478  << "[" << fvm.psi().name() << fvm.dimensions()/dimVolume << " ] "
1479  << op
1480  << " [" << dt.name() << dt.dimensions() << " ]"
1481  << abort(FatalError);
1482  }
1483 }
1484 
1485 
1486 template<class Type>
1488 (
1489  const fvMatrix<Type>& A
1490 )
1491 {
1492  tmp<Foam::fvMatrix<Type>> tAcorr = A - (A & A.psi());
1493 
1494  // Delete the faceFluxCorrection from the correction matrix
1495  // as it does not have a clear meaning or purpose
1496  deleteDemandDrivenData(tAcorr.ref().faceFluxCorrectionPtr());
1497 
1498  return tAcorr;
1499 }
1500 
1501 
1502 template<class Type>
1504 (
1505  const tmp<fvMatrix<Type>>& tA
1506 )
1507 {
1508  tmp<Foam::fvMatrix<Type>> tAcorr = tA - (tA() & tA().psi());
1509 
1510  // Delete the faceFluxCorrection from the correction matrix
1511  // as it does not have a clear meaning or purpose
1512  deleteDemandDrivenData(tAcorr.ref().faceFluxCorrectionPtr());
1513 
1514  return tAcorr;
1515 }
1516 
1517 
1518 // * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * //
1519 
1520 template<class Type>
1521 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
1522 (
1523  const fvMatrix<Type>& A,
1524  const fvMatrix<Type>& B
1525 )
1526 {
1527  checkMethod(A, B, "==");
1528  return (A - B);
1529 }
1530 
1531 template<class Type>
1532 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
1533 (
1534  const tmp<fvMatrix<Type>>& tA,
1535  const fvMatrix<Type>& B
1536 )
1537 {
1538  checkMethod(tA(), B, "==");
1539  return (tA - B);
1540 }
1541 
1542 template<class Type>
1543 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
1544 (
1545  const fvMatrix<Type>& A,
1546  const tmp<fvMatrix<Type>>& tB
1547 )
1548 {
1549  checkMethod(A, tB(), "==");
1550  return (A - tB);
1551 }
1552 
1553 template<class Type>
1554 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
1555 (
1556  const tmp<fvMatrix<Type>>& tA,
1557  const tmp<fvMatrix<Type>>& tB
1558 )
1559 {
1560  checkMethod(tA(), tB(), "==");
1561  return (tA - tB);
1562 }
1563 
1564 template<class Type>
1565 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
1566 (
1567  const fvMatrix<Type>& A,
1568  const DimensionedField<Type, volMesh>& su
1569 )
1570 {
1571  checkMethod(A, su, "==");
1572  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
1573  tC.ref().source() += su.mesh().V()*su.field();
1574  return tC;
1575 }
1576 
1577 template<class Type>
1578 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
1579 (
1580  const fvMatrix<Type>& A,
1581  const tmp<DimensionedField<Type, volMesh>>& tsu
1582 )
1583 {
1584  checkMethod(A, tsu(), "==");
1585  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
1586  tC.ref().source() += tsu().mesh().V()*tsu().field();
1587  tsu.clear();
1588  return tC;
1589 }
1590 
1591 template<class Type>
1592 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
1593 (
1594  const fvMatrix<Type>& A,
1595  const tmp<VolField<Type>>& tsu
1596 )
1597 {
1598  checkMethod(A, tsu(), "==");
1599  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
1600  tC.ref().source() += tsu().mesh().V()*tsu().primitiveField();
1601  tsu.clear();
1602  return tC;
1603 }
1604 
1605 template<class Type>
1606 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
1607 (
1608  const tmp<fvMatrix<Type>>& tA,
1609  const DimensionedField<Type, volMesh>& su
1610 )
1611 {
1612  checkMethod(tA(), su, "==");
1613  tmp<fvMatrix<Type>> tC(tA.ptr());
1614  tC.ref().source() += su.mesh().V()*su.field();
1615  return tC;
1616 }
1617 
1618 template<class Type>
1619 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
1620 (
1621  const tmp<fvMatrix<Type>>& tA,
1622  const tmp<DimensionedField<Type, volMesh>>& tsu
1623 )
1624 {
1625  checkMethod(tA(), tsu(), "==");
1626  tmp<fvMatrix<Type>> tC(tA.ptr());
1627  tC.ref().source() += tsu().mesh().V()*tsu().field();
1628  tsu.clear();
1629  return tC;
1630 }
1631 
1632 template<class Type>
1633 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
1634 (
1635  const tmp<fvMatrix<Type>>& tA,
1636  const tmp<VolField<Type>>& tsu
1637 )
1638 {
1639  checkMethod(tA(), tsu(), "==");
1640  tmp<fvMatrix<Type>> tC(tA.ptr());
1641  tC.ref().source() += tsu().mesh().V()*tsu().primitiveField();
1642  tsu.clear();
1643  return tC;
1644 }
1645 
1646 template<class Type>
1647 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
1648 (
1649  const fvMatrix<Type>& A,
1650  const dimensioned<Type>& su
1651 )
1652 {
1653  checkMethod(A, su, "==");
1654  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
1655  tC.ref().source() += A.psi().mesh().V()*su.value();
1656  return tC;
1657 }
1658 
1659 template<class Type>
1660 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
1661 (
1662  const tmp<fvMatrix<Type>>& tA,
1663  const dimensioned<Type>& su
1664 )
1665 {
1666  checkMethod(tA(), su, "==");
1667  tmp<fvMatrix<Type>> tC(tA.ptr());
1668  tC.ref().source() += tC().psi().mesh().V()*su.value();
1669  return tC;
1670 }
1671 
1672 template<class Type>
1673 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
1674 (
1675  const fvMatrix<Type>& A,
1676  const zero&
1677 )
1678 {
1679  return A;
1680 }
1681 
1682 
1683 template<class Type>
1684 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
1685 (
1686  const tmp<fvMatrix<Type>>& tA,
1687  const zero&
1688 )
1689 {
1690  return tA;
1691 }
1692 
1693 
1694 template<class Type>
1695 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
1696 (
1697  const fvMatrix<Type>& A
1698 )
1699 {
1700  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
1701  tC.ref().negate();
1702  return tC;
1703 }
1704 
1705 template<class Type>
1706 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
1707 (
1708  const tmp<fvMatrix<Type>>& tA
1709 )
1710 {
1711  tmp<fvMatrix<Type>> tC(tA.ptr());
1712  tC.ref().negate();
1713  return tC;
1714 }
1715 
1716 
1717 template<class Type>
1718 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
1719 (
1720  const fvMatrix<Type>& A,
1721  const fvMatrix<Type>& B
1722 )
1723 {
1724  checkMethod(A, B, "+");
1725  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
1726  tC.ref() += B;
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 fvMatrix<Type>& B
1735 )
1736 {
1737  checkMethod(tA(), B, "+");
1738  tmp<fvMatrix<Type>> tC(tA.ptr());
1739  tC.ref() += B;
1740  return tC;
1741 }
1742 
1743 template<class Type>
1744 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
1745 (
1746  const fvMatrix<Type>& A,
1747  const tmp<fvMatrix<Type>>& tB
1748 )
1749 {
1750  checkMethod(A, tB(), "+");
1751  tmp<fvMatrix<Type>> tC(tB.ptr());
1752  tC.ref() += A;
1753  return tC;
1754 }
1755 
1756 template<class Type>
1757 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
1758 (
1759  const tmp<fvMatrix<Type>>& tA,
1760  const tmp<fvMatrix<Type>>& tB
1761 )
1762 {
1763  checkMethod(tA(), tB(), "+");
1764  tmp<fvMatrix<Type>> tC(tA.ptr());
1765  tC.ref() += tB();
1766  tB.clear();
1767  return tC;
1768 }
1769 
1770 template<class Type>
1771 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
1772 (
1773  const fvMatrix<Type>& A,
1774  const DimensionedField<Type, volMesh>& su
1775 )
1776 {
1777  checkMethod(A, su, "+");
1778  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
1779  tC.ref().source() -= su.mesh().V()*su.field();
1780  return tC;
1781 }
1782 
1783 template<class Type>
1784 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
1785 (
1786  const fvMatrix<Type>& A,
1787  const tmp<DimensionedField<Type, volMesh>>& tsu
1788 )
1789 {
1790  checkMethod(A, tsu(), "+");
1791  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
1792  tC.ref().source() -= tsu().mesh().V()*tsu().field();
1793  tsu.clear();
1794  return tC;
1795 }
1796 
1797 template<class Type>
1798 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
1799 (
1800  const fvMatrix<Type>& A,
1801  const tmp<VolField<Type>>& tsu
1802 )
1803 {
1804  checkMethod(A, tsu(), "+");
1805  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
1806  tC.ref().source() -= tsu().mesh().V()*tsu().primitiveField();
1807  tsu.clear();
1808  return tC;
1809 }
1810 
1811 template<class Type>
1812 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
1813 (
1814  const tmp<fvMatrix<Type>>& tA,
1815  const DimensionedField<Type, volMesh>& su
1816 )
1817 {
1818  checkMethod(tA(), su, "+");
1819  tmp<fvMatrix<Type>> tC(tA.ptr());
1820  tC.ref().source() -= su.mesh().V()*su.field();
1821  return tC;
1822 }
1823 
1824 template<class Type>
1825 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
1826 (
1827  const tmp<fvMatrix<Type>>& tA,
1828  const tmp<DimensionedField<Type, volMesh>>& tsu
1829 )
1830 {
1831  checkMethod(tA(), tsu(), "+");
1832  tmp<fvMatrix<Type>> tC(tA.ptr());
1833  tC.ref().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+
1840 (
1841  const tmp<fvMatrix<Type>>& tA,
1842  const tmp<VolField<Type>>& tsu
1843 )
1844 {
1845  checkMethod(tA(), tsu(), "+");
1846  tmp<fvMatrix<Type>> tC(tA.ptr());
1847  tC.ref().source() -= tsu().mesh().V()*tsu().primitiveField();
1848  tsu.clear();
1849  return tC;
1850 }
1851 
1852 template<class Type>
1853 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
1854 (
1855  const DimensionedField<Type, volMesh>& su,
1856  const fvMatrix<Type>& A
1857 )
1858 {
1859  checkMethod(A, su, "+");
1860  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
1861  tC.ref().source() -= su.mesh().V()*su.field();
1862  return tC;
1863 }
1864 
1865 template<class Type>
1866 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
1867 (
1868  const tmp<DimensionedField<Type, volMesh>>& tsu,
1869  const fvMatrix<Type>& A
1870 )
1871 {
1872  checkMethod(A, tsu(), "+");
1873  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
1874  tC.ref().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+
1881 (
1882  const tmp<VolField<Type>>& tsu,
1883  const fvMatrix<Type>& A
1884 )
1885 {
1886  checkMethod(A, tsu(), "+");
1887  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
1888  tC.ref().source() -= tsu().mesh().V()*tsu().primitiveField();
1889  tsu.clear();
1890  return tC;
1891 }
1892 
1893 template<class Type>
1894 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
1895 (
1896  const DimensionedField<Type, volMesh>& su,
1897  const tmp<fvMatrix<Type>>& tA
1898 )
1899 {
1900  checkMethod(tA(), su, "+");
1901  tmp<fvMatrix<Type>> tC(tA.ptr());
1902  tC.ref().source() -= su.mesh().V()*su.field();
1903  return tC;
1904 }
1905 
1906 template<class Type>
1907 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
1908 (
1909  const tmp<DimensionedField<Type, volMesh>>& tsu,
1910  const tmp<fvMatrix<Type>>& tA
1911 )
1912 {
1913  checkMethod(tA(), tsu(), "+");
1914  tmp<fvMatrix<Type>> tC(tA.ptr());
1915  tC.ref().source() -= tsu().mesh().V()*tsu().field();
1916  tsu.clear();
1917  return tC;
1918 }
1919 
1920 template<class Type>
1921 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
1922 (
1923  const tmp<VolField<Type>>& tsu,
1924  const tmp<fvMatrix<Type>>& tA
1925 )
1926 {
1927  checkMethod(tA(), tsu(), "+");
1928  tmp<fvMatrix<Type>> tC(tA.ptr());
1929  tC.ref().source() -= tsu().mesh().V()*tsu().primitiveField();
1930  tsu.clear();
1931  return tC;
1932 }
1933 
1934 
1935 template<class Type>
1936 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
1937 (
1938  const fvMatrix<Type>& A,
1939  const fvMatrix<Type>& B
1940 )
1941 {
1942  checkMethod(A, B, "-");
1943  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
1944  tC.ref() -= B;
1945  return tC;
1946 }
1947 
1948 template<class Type>
1949 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
1950 (
1951  const tmp<fvMatrix<Type>>& tA,
1952  const fvMatrix<Type>& B
1953 )
1954 {
1955  checkMethod(tA(), B, "-");
1956  tmp<fvMatrix<Type>> tC(tA.ptr());
1957  tC.ref() -= B;
1958  return tC;
1959 }
1960 
1961 template<class Type>
1962 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
1963 (
1964  const fvMatrix<Type>& A,
1965  const tmp<fvMatrix<Type>>& tB
1966 )
1967 {
1968  checkMethod(A, tB(), "-");
1969  tmp<fvMatrix<Type>> tC(tB.ptr());
1970  tC.ref() -= A;
1971  tC.ref().negate();
1972  return tC;
1973 }
1974 
1975 template<class Type>
1976 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
1977 (
1978  const tmp<fvMatrix<Type>>& tA,
1979  const tmp<fvMatrix<Type>>& tB
1980 )
1981 {
1982  checkMethod(tA(), tB(), "-");
1983  tmp<fvMatrix<Type>> tC(tA.ptr());
1984  tC.ref() -= tB();
1985  tB.clear();
1986  return tC;
1987 }
1988 
1989 template<class Type>
1990 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
1991 (
1992  const fvMatrix<Type>& A,
1993  const DimensionedField<Type, volMesh>& su
1994 )
1995 {
1996  checkMethod(A, su, "-");
1997  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
1998  tC.ref().source() += su.mesh().V()*su.field();
1999  return tC;
2000 }
2001 
2002 template<class Type>
2003 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
2004 (
2005  const fvMatrix<Type>& A,
2006  const tmp<DimensionedField<Type, volMesh>>& tsu
2007 )
2008 {
2009  checkMethod(A, tsu(), "-");
2010  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
2011  tC.ref().source() += tsu().mesh().V()*tsu().field();
2012  tsu.clear();
2013  return tC;
2014 }
2015 
2016 template<class Type>
2017 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
2018 (
2019  const fvMatrix<Type>& A,
2020  const tmp<VolField<Type>>& tsu
2021 )
2022 {
2023  checkMethod(A, tsu(), "-");
2024  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
2025  tC.ref().source() += tsu().mesh().V()*tsu().primitiveField();
2026  tsu.clear();
2027  return tC;
2028 }
2029 
2030 template<class Type>
2031 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
2032 (
2033  const tmp<fvMatrix<Type>>& tA,
2034  const DimensionedField<Type, volMesh>& su
2035 )
2036 {
2037  checkMethod(tA(), su, "-");
2038  tmp<fvMatrix<Type>> tC(tA.ptr());
2039  tC.ref().source() += su.mesh().V()*su.field();
2040  return tC;
2041 }
2042 
2043 template<class Type>
2044 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
2045 (
2046  const tmp<fvMatrix<Type>>& tA,
2047  const tmp<DimensionedField<Type, volMesh>>& tsu
2048 )
2049 {
2050  checkMethod(tA(), tsu(), "-");
2051  tmp<fvMatrix<Type>> tC(tA.ptr());
2052  tC.ref().source() += tsu().mesh().V()*tsu().field();
2053  tsu.clear();
2054  return tC;
2055 }
2056 
2057 template<class Type>
2058 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
2059 (
2060  const tmp<fvMatrix<Type>>& tA,
2061  const tmp<VolField<Type>>& tsu
2062 )
2063 {
2064  checkMethod(tA(), tsu(), "-");
2065  tmp<fvMatrix<Type>> tC(tA.ptr());
2066  tC.ref().source() += tsu().mesh().V()*tsu().primitiveField();
2067  tsu.clear();
2068  return tC;
2069 }
2070 
2071 template<class Type>
2072 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
2073 (
2074  const DimensionedField<Type, volMesh>& su,
2075  const fvMatrix<Type>& A
2076 )
2077 {
2078  checkMethod(A, su, "-");
2079  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
2080  tC.ref().negate();
2081  tC.ref().source() -= su.mesh().V()*su.field();
2082  return tC;
2083 }
2084 
2085 template<class Type>
2086 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
2087 (
2088  const tmp<DimensionedField<Type, volMesh>>& tsu,
2089  const fvMatrix<Type>& A
2090 )
2091 {
2092  checkMethod(A, tsu(), "-");
2093  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
2094  tC.ref().negate();
2095  tC.ref().source() -= tsu().mesh().V()*tsu().field();
2096  tsu.clear();
2097  return tC;
2098 }
2099 
2100 template<class Type>
2101 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
2102 (
2103  const tmp<VolField<Type>>& tsu,
2104  const fvMatrix<Type>& A
2105 )
2106 {
2107  checkMethod(A, tsu(), "-");
2108  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
2109  tC.ref().negate();
2110  tC.ref().source() -= tsu().mesh().V()*tsu().primitiveField();
2111  tsu.clear();
2112  return tC;
2113 }
2114 
2115 template<class Type>
2116 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
2117 (
2118  const DimensionedField<Type, volMesh>& su,
2119  const tmp<fvMatrix<Type>>& tA
2120 )
2121 {
2122  checkMethod(tA(), su, "-");
2123  tmp<fvMatrix<Type>> tC(tA.ptr());
2124  tC.ref().negate();
2125  tC.ref().source() -= su.mesh().V()*su.field();
2126  return tC;
2127 }
2128 
2129 template<class Type>
2130 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
2131 (
2132  const tmp<DimensionedField<Type, volMesh>>& tsu,
2133  const tmp<fvMatrix<Type>>& tA
2134 )
2135 {
2136  checkMethod(tA(), tsu(), "-");
2137  tmp<fvMatrix<Type>> tC(tA.ptr());
2138  tC.ref().negate();
2139  tC.ref().source() -= tsu().mesh().V()*tsu().field();
2140  tsu.clear();
2141  return tC;
2142 }
2143 
2144 template<class Type>
2145 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
2146 (
2147  const tmp<VolField<Type>>& tsu,
2148  const tmp<fvMatrix<Type>>& tA
2149 )
2150 {
2151  checkMethod(tA(), tsu(), "-");
2152  tmp<fvMatrix<Type>> tC(tA.ptr());
2153  tC.ref().negate();
2154  tC.ref().source() -= tsu().mesh().V()*tsu().primitiveField();
2155  tsu.clear();
2156  return tC;
2157 }
2158 
2159 template<class Type>
2160 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
2161 (
2162  const fvMatrix<Type>& A,
2163  const dimensioned<Type>& su
2164 )
2165 {
2166  checkMethod(A, su, "+");
2167  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
2168  tC.ref().source() -= su.value()*A.psi().mesh().V();
2169  return tC;
2170 }
2171 
2172 template<class Type>
2173 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
2174 (
2175  const tmp<fvMatrix<Type>>& tA,
2176  const dimensioned<Type>& su
2177 )
2178 {
2179  checkMethod(tA(), su, "+");
2180  tmp<fvMatrix<Type>> tC(tA.ptr());
2181  tC.ref().source() -= su.value()*tC().psi().mesh().V();
2182  return tC;
2183 }
2184 
2185 template<class Type>
2186 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
2187 (
2188  const dimensioned<Type>& su,
2189  const fvMatrix<Type>& A
2190 )
2191 {
2192  checkMethod(A, su, "+");
2193  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
2194  tC.ref().source() -= su.value()*A.psi().mesh().V();
2195  return tC;
2196 }
2197 
2198 template<class Type>
2199 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
2200 (
2201  const dimensioned<Type>& su,
2202  const tmp<fvMatrix<Type>>& tA
2203 )
2204 {
2205  checkMethod(tA(), su, "+");
2206  tmp<fvMatrix<Type>> tC(tA.ptr());
2207  tC.ref().source() -= su.value()*tC().psi().mesh().V();
2208  return tC;
2209 }
2210 
2211 template<class Type>
2212 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
2213 (
2214  const fvMatrix<Type>& A,
2215  const dimensioned<Type>& su
2216 )
2217 {
2218  checkMethod(A, su, "-");
2219  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
2220  tC.ref().source() += su.value()*tC().psi().mesh().V();
2221  return tC;
2222 }
2223 
2224 template<class Type>
2225 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
2226 (
2227  const tmp<fvMatrix<Type>>& tA,
2228  const dimensioned<Type>& su
2229 )
2230 {
2231  checkMethod(tA(), su, "-");
2232  tmp<fvMatrix<Type>> tC(tA.ptr());
2233  tC.ref().source() += su.value()*tC().psi().mesh().V();
2234  return tC;
2235 }
2236 
2237 template<class Type>
2238 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
2239 (
2240  const dimensioned<Type>& su,
2241  const fvMatrix<Type>& A
2242 )
2243 {
2244  checkMethod(A, su, "-");
2245  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
2246  tC.ref().negate();
2247  tC.ref().source() -= su.value()*A.psi().mesh().V();
2248  return tC;
2249 }
2250 
2251 template<class Type>
2252 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
2253 (
2254  const dimensioned<Type>& su,
2255  const tmp<fvMatrix<Type>>& tA
2256 )
2257 {
2258  checkMethod(tA(), su, "-");
2259  tmp<fvMatrix<Type>> tC(tA.ptr());
2260  tC.ref().negate();
2261  tC.ref().source() -= su.value()*tC().psi().mesh().V();
2262  return tC;
2263 }
2264 
2265 
2266 template<class Type>
2267 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator*
2268 (
2269  const volScalarField::Internal& dsf,
2270  const fvMatrix<Type>& A
2271 )
2272 {
2273  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
2274  tC.ref() *= dsf;
2275  return tC;
2276 }
2277 
2278 template<class Type>
2279 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator*
2280 (
2281  const tmp<volScalarField::Internal>& tdsf,
2282  const fvMatrix<Type>& A
2283 )
2284 {
2285  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
2286  tC.ref() *= tdsf;
2287  return tC;
2288 }
2289 
2290 template<class Type>
2291 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator*
2292 (
2293  const tmp<volScalarField>& tvsf,
2294  const fvMatrix<Type>& A
2295 )
2296 {
2297  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
2298  tC.ref() *= tvsf;
2299  return tC;
2300 }
2301 
2302 template<class Type>
2303 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator*
2304 (
2305  const volScalarField::Internal& dsf,
2306  const tmp<fvMatrix<Type>>& tA
2307 )
2308 {
2309  tmp<fvMatrix<Type>> tC(tA.ptr());
2310  tC.ref() *= dsf;
2311  return tC;
2312 }
2313 
2314 template<class Type>
2315 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator*
2316 (
2317  const tmp<volScalarField::Internal>& tdsf,
2318  const tmp<fvMatrix<Type>>& tA
2319 )
2320 {
2321  tmp<fvMatrix<Type>> tC(tA.ptr());
2322  tC.ref() *= tdsf;
2323  return tC;
2324 }
2325 
2326 template<class Type>
2327 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator*
2328 (
2329  const tmp<volScalarField>& tvsf,
2330  const tmp<fvMatrix<Type>>& tA
2331 )
2332 {
2333  tmp<fvMatrix<Type>> tC(tA.ptr());
2334  tC.ref() *= tvsf;
2335  return tC;
2336 }
2337 
2338 template<class Type>
2339 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator*
2340 (
2341  const dimensioned<scalar>& ds,
2342  const fvMatrix<Type>& A
2343 )
2344 {
2345  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
2346  tC.ref() *= ds;
2347  return tC;
2348 }
2349 
2350 template<class Type>
2351 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator*
2352 (
2353  const dimensioned<scalar>& ds,
2354  const tmp<fvMatrix<Type>>& tA
2355 )
2356 {
2357  tmp<fvMatrix<Type>> tC(tA.ptr());
2358  tC.ref() *= ds;
2359  return tC;
2360 }
2361 
2362 
2363 template<class Type>
2364 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator/
2365 (
2366  const fvMatrix<Type>& A,
2367  const volScalarField::Internal& dsf
2368 )
2369 {
2370  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
2371  tC.ref() /= dsf;
2372  return tC;
2373 }
2374 
2375 template<class Type>
2376 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator/
2377 (
2378  const fvMatrix<Type>& A,
2379  const tmp<volScalarField::Internal>& tdsf
2380 )
2381 {
2382  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
2383  tC.ref() /= tdsf;
2384  return tC;
2385 }
2386 
2387 template<class Type>
2388 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator/
2389 (
2390  const fvMatrix<Type>& A,
2391  const tmp<volScalarField>& tvsf
2392 )
2393 {
2394  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
2395  tC.ref() /= tvsf;
2396  return tC;
2397 }
2398 
2399 template<class Type>
2400 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator/
2401 (
2402  const tmp<fvMatrix<Type>>& tA,
2403  const volScalarField::Internal& dsf
2404 )
2405 {
2406  tmp<fvMatrix<Type>> tC(tA.ptr());
2407  tC.ref() /= dsf;
2408  return tC;
2409 }
2410 
2411 template<class Type>
2412 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator/
2413 (
2414  const tmp<fvMatrix<Type>>& tA,
2415  const tmp<volScalarField::Internal>& tdsf
2416 )
2417 {
2418  tmp<fvMatrix<Type>> tC(tA.ptr());
2419  tC.ref() /= tdsf;
2420  return tC;
2421 }
2422 
2423 template<class Type>
2424 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator/
2425 (
2426  const tmp<fvMatrix<Type>>& tA,
2427  const tmp<volScalarField>& tvsf
2428 )
2429 {
2430  tmp<fvMatrix<Type>> tC(tA.ptr());
2431  tC.ref() /= tvsf;
2432  return tC;
2433 }
2434 
2435 template<class Type>
2436 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator/
2437 (
2438  const fvMatrix<Type>& A,
2439  const dimensioned<scalar>& ds
2440 )
2441 {
2442  tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
2443  tC.ref() /= ds;
2444  return tC;
2445 }
2446 
2447 template<class Type>
2448 Foam::tmp<Foam::fvMatrix<Type>> Foam::operator/
2449 (
2450  const tmp<fvMatrix<Type>>& tA,
2451  const dimensioned<scalar>& ds
2452 )
2453 {
2454  tmp<fvMatrix<Type>> tC(tA.ptr());
2455  tC.ref() /= ds;
2456  return tC;
2457 }
2458 
2459 
2460 template<class Type>
2462 Foam::operator&
2463 (
2464  const fvMatrix<Type>& M,
2465  const DimensionedField<Type, volMesh>& psi
2466 )
2467 {
2468  tmp<VolField<Type>> tMphi
2469  (
2471  (
2472  "M&" + psi.name(),
2473  psi.mesh(),
2474  M.dimensions()/dimVol,
2475  extrapolatedCalculatedFvPatchScalarField::typeName
2476  )
2477  );
2478  VolField<Type>& Mphi = tMphi.ref();
2479 
2480  // Loop over field components
2481  if (M.hasDiag())
2482  {
2483  for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
2484  {
2485  scalarField psiCmpt(psi.field().component(cmpt));
2486  scalarField boundaryDiagCmpt(M.diag());
2487  M.addBoundaryDiag(boundaryDiagCmpt, cmpt);
2488  Mphi.primitiveFieldRef().replace(cmpt, -boundaryDiagCmpt*psiCmpt);
2489  }
2490  }
2491  else
2492  {
2493  Mphi.primitiveFieldRef() = Zero;
2494  }
2495 
2496  Mphi.primitiveFieldRef() += M.lduMatrix::H(psi.field()) + M.source();
2497  M.addBoundarySource(Mphi.primitiveFieldRef());
2498 
2499  Mphi.primitiveFieldRef() /= -psi.mesh().V();
2500  Mphi.correctBoundaryConditions();
2501 
2502  return tMphi;
2503 }
2504 
2505 template<class Type>
2507 Foam::operator&
2508 (
2509  const fvMatrix<Type>& M,
2510  const tmp<DimensionedField<Type, volMesh>>& tpsi
2511 )
2512 {
2513  tmp<VolField<Type>> tMpsi = M & tpsi();
2514  tpsi.clear();
2515  return tMpsi;
2516 }
2517 
2518 template<class Type>
2520 Foam::operator&
2521 (
2522  const fvMatrix<Type>& M,
2523  const tmp<VolField<Type>>& tpsi
2524 )
2525 {
2526  tmp<VolField<Type>> tMpsi = M & tpsi();
2527  tpsi.clear();
2528  return tMpsi;
2529 }
2530 
2531 template<class Type>
2533 Foam::operator&
2534 (
2535  const tmp<fvMatrix<Type>>& tM,
2536  const DimensionedField<Type, volMesh>& psi
2537 )
2538 {
2539  tmp<VolField<Type>> tMpsi = tM() & psi;
2540  tM.clear();
2541  return tMpsi;
2542 }
2543 
2544 template<class Type>
2546 Foam::operator&
2547 (
2548  const tmp<fvMatrix<Type>>& tM,
2549  const tmp<DimensionedField<Type, volMesh>>& tpsi
2550 )
2551 {
2552  tmp<VolField<Type>> tMpsi = tM() & tpsi();
2553  tM.clear();
2554  tpsi.clear();
2555  return tMpsi;
2556 }
2557 
2558 template<class Type>
2560 Foam::operator&
2561 (
2562  const tmp<fvMatrix<Type>>& tM,
2563  const tmp<VolField<Type>>& tpsi
2564 )
2565 {
2566  tmp<VolField<Type>> tMpsi = tM() & tpsi();
2567  tM.clear();
2568  tpsi.clear();
2569  return tMpsi;
2570 }
2571 
2572 
2573 // * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * //
2574 
2575 template<class Type>
2576 Foam::Ostream& Foam::operator<<(Ostream& os, const fvMatrix<Type>& fvm)
2577 {
2578  os << static_cast<const lduMatrix&>(fvm) << nl
2579  << fvm.dimensions_ << nl
2580  << fvm.source_ << nl
2581  << fvm.internalCoeffs_ << nl
2582  << fvm.boundaryCoeffs_ << endl;
2583 
2584  os.check("Ostream& operator<<(Ostream&, fvMatrix<Type>&");
2585 
2586  return os;
2587 }
2588 
2589 
2590 // * * * * * * * * * * * * * * * * Solvers * * * * * * * * * * * * * * * * * //
2591 
2592 #include "fvMatrixSolve.C"
2593 
2594 // ************************************************************************* //
static const Foam::dimensionedScalar A("A", Foam::dimPressure, 611.21)
static const Foam::dimensionedScalar B("B", Foam::dimless, 18.678)
static const Foam::dimensionedScalar D("D", Foam::dimTemperature, 257.14)
EaEqn relax()
#define M(I)
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
const dimensionSet & dimensions() const
Return dimensions.
const Mesh & mesh() const
Return mesh.
const Field< Type > & field() const
static tmp< DimensionedField< Type, GeoMesh > > New(const word &name, const Mesh &mesh, const dimensionSet &, const Field< Type > &)
Return a temporary field constructed from name, mesh,.
Generic field type.
Definition: FieldField.H:77
Pre-declare SubField and related Field type.
Definition: Field.H:82
void replace(const direction, const UList< cmptType > &)
Replace a component field of the field.
Definition: Field.C:467
void negate()
Negate this field.
Definition: Field.C:446
void updateCoeffs()
Update the boundary condition coefficients.
Generic GeometricField class.
Internal::FieldType & primitiveFieldRef()
Return a reference to the internal field.
DimensionedField< Type, GeoMesh > Internal
Type of the internal field from which this GeometricField is derived.
static tmp< GeometricField< Type, PatchField, GeoMesh > > New(const word &name, const Internal &, const PtrList< PatchField< Type >> &)
Return a temporary field constructed from name,.
void replace(const direction, const GeometricField< cmptType, PatchField, GeoMesh > &)
void correctBoundaryConditions()
Correct boundary field.
Boundary & boundaryFieldRef()
Return a reference to the boundary field.
const word & name() const
Return name.
Definition: IOobject.H:310
virtual bool check(const char *operation) const
Check IOstream status for given operation.
Definition: IOstream.C:92
An Istream is an abstract base class for all input systems (streams, files, token lists etc)....
Definition: Istream.H:60
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
Definition: Ostream.H:57
label size() const
Return the number of elements in the UList.
Definition: UListI.H:311
static int & msgType()
Message tag of standard messages.
Definition: UPstream.H:476
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:60
Dimension set for the base types.
Definition: dimensionSet.H:122
Generic dimensioned Type class.
const dimensionSet & dimensions() const
Return const reference to dimensions.
const word & name() const
Return const reference to name.
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:76
A special matrix type and solver, designed for finite volume solutions of scalar equations....
Definition: fvMatrix.H:118
void operator-=(const fvMatrix< Type > &)
Definition: fvMatrix.C:1136
tmp< Field< Type > > DD() const
Return the matrix Type diagonal.
Definition: fvMatrix.C:785
tmp< volScalarField > A() const
Return the central coefficient.
Definition: fvMatrix.C:809
void addToInternalField(const labelUList &addr, const Field< Type2 > &pf, Field< Type2 > &intf) const
Add patch contribution to internal field.
Definition: fvMatrix.C:40
void setValue(const label celli, const Type &value)
Set solution in the given cell to the specified value.
Definition: fvMatrix.C:179
void boundaryManipulate(typename VolField< Type >::Boundary &values)
Manipulate based on a boundary field.
Definition: fvMatrix.C:763
void relax()
Relax matrix (for steady-state solution).
Definition: fvMatrix.C:755
tmp< volScalarField::Internal > Su() const
Return the explicit source.
Definition: fvMatrix.C:830
tmp< volScalarField > H1() const
Return H(1)
Definition: fvMatrix.C:922
tmp< VolField< Type > > H() const
Return the H operation source.
Definition: fvMatrix.C:867
void operator+=(const fvMatrix< Type > &)
Definition: fvMatrix.C:1102
virtual ~fvMatrix()
Destructor.
Definition: fvMatrix.C:482
VolField< Type > & psi()
Definition: fvMatrix.H:289
tmp< SurfaceField< Type > > flux() const
Return the face-flux field from the matrix.
Definition: fvMatrix.C:963
void addCmptAvBoundaryDiag(scalarField &diag) const
Definition: fvMatrix.C:131
void addBoundarySource(Field< Type > &source, const bool couples=true) const
Definition: fvMatrix.C:147
fvMatrix(const VolField< Type > &, const dimensionSet &)
Construct given a field to solve for.
Definition: fvMatrix.C:285
tmp< volScalarField::Internal > Sp() const
Return the implicit source.
Definition: fvMatrix.C:848
void setReference(const label celli, const Type &value, const bool forceReference=false)
Set reference level for solution.
Definition: fvMatrix.C:563
scalar relaxationFactor() const
Return the relaxation factor for this equation.
Definition: fvMatrix.C:578
void addBoundaryDiag(scalarField &diag, const direction cmpt) const
Definition: fvMatrix.C:113
void negate()
Definition: fvMatrix.C:1087
void setValues(const labelUList &cells, const ListType< Type > &values)
Set solution in given cells to the specified values.
Definition: fvMatrix.C:502
tmp< fvMatrix< Type > > clone() const
Clone.
Definition: fvMatrix.C:470
tmp< scalarField > D() const
Return the matrix scalar diagonal.
Definition: fvMatrix.C:776
void subtractFromInternalField(const labelUList &addr, const Field< Type2 > &pf, Field< Type2 > &intf) const
Subtract patch contribution from internal field.
Definition: fvMatrix.C:77
const dimensionSet & dimensions() const
Definition: fvMatrix.H:302
void operator=(const fvMatrix< Type > &)
Definition: fvMatrix.C:1043
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:101
const labelUList & owner() const
Internal face owner.
Definition: fvMesh.H:451
const UCompactListList< label > & polyBFacePatches() const
Return poly-bFace-patch addressing.
Definition: fvMesh.C:951
const UCompactListList< label > & polyBFacePatchFaces() const
Return poly-bFace-patch-face addressing.
Definition: fvMesh.C:1027
const labelUList & neighbour() const
Internal face neighbour.
Definition: fvMesh.H:457
Abstract base class with a fat-interface to all derived classes covering all possible ways in which t...
Definition: fvPatchField.H:87
virtual bool coupled() const
Return true if this patch field is coupled.
Definition: fvPatchField.H:327
virtual tmp< Field< Type > > patchNeighbourField(const Pstream::commsTypes commsType=Pstream::commsTypes::blocking) const
Return patchField on the opposite patch of a coupled patch.
Definition: fvPatchField.H:436
lduMatrix is a general matrix class in which the coefficients are stored as three arrays,...
Definition: lduMatrix.H:80
tmp< scalarField > H1() const
tmp< Field< Type > > faceH(const Field< Type > &) const
void operator=(const lduMatrix &)
tmp< Field< Type > > H(const Field< Type > &) const
void operator*=(const scalarField &)
scalarField & diag()
Definition: lduMatrix.C:186
void operator/=(const scalarField &)
void operator+=(const lduMatrix &)
void operator-=(const lduMatrix &)
bool hasDiag() const
Definition: lduMatrix.H:580
Traits class for primitives.
Definition: pTraits.H:53
label nInternalFaces() const
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
const cellList & cells() const
Reference counter for various OpenFOAM components.
Definition: refCount.H:50
label eventNo() const
Event number at last update.
Definition: regIOobjectI.H:89
A class for managing temporary objects.
Definition: tmp.H:55
T * ptr() const
Return tmp pointer for reuse.
Definition: tmpI.H:205
void clear() const
If object pointer points to valid object:
Definition: tmpI.H:237
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
A class representing the concept of 0 used to avoid unnecessary manipulations for objects that are kn...
Definition: zero.H:50
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
Calculate the matrix for the first temporal derivative.
label patchi
const cellShapeList & cells
const fvPatchList & patches
const volScalarField & psi
volScalarField alpha(IOobject("alpha", runTime.name(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
#define InfoInFunction
Report an information message using Foam::Info.
const dimensionedScalar c
Speed of light in a vacuum.
tmp< fvMatrix< Type > > S(const Pair< tmp< volScalarField::Internal >> &, const VolField< Type > &)
tmp< fvMatrix< Type > > ddt(const VolField< Type > &vf)
Definition: fvmDdt.C:46
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
static const zero Zero
Definition: zero.H:97
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
void component(FieldField< Field, typename FieldField< Field, Type >::cmptType > &sf, const FieldField< Field, Type > &f, const direction d)
errorManip< error > abort(error &err)
Definition: errorManip.H:131
void cmptMin(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
void checkMethod(const fvMatrix< Type > &, const fvMatrix< Type > &, const char *)
Definition: fvMatrix.C:1415
void deleteDemandDrivenData(DataPtr &dataPtr)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
tmp< DimensionedField< typename DimensionedField< Type, GeoMesh >::cmptType, GeoMesh >> cmptAv(const DimensionedField< Type, GeoMesh > &df)
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
const dimensionSet dimVol
void diag(pointPatchField< vector > &, const pointPatchField< tensor > &)
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
const dimensionSet dimVolume
dimensioned< scalar > mag(const dimensioned< Type > &)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
tmp< fvMatrix< Type > > correction(const fvMatrix< Type > &)
Return the correction form of the given matrix.
dimensioned< Type > cmptMultiply(const dimensioned< Type > &, const dimensioned< Type > &)
dimensionSet cmptMag(const dimensionSet &)
Definition: dimensionSet.C:275
void cmptMax(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
error FatalError
Ostream & operator<<(Ostream &, const ensightPart &)
static const char nl
Definition: Ostream.H:260
uint8_t direction
Definition: direction.H:45
dictionary fractions(initialConditions.subDict("fractions"))
faceListList boundary(nPatches)
Foam::surfaceFields.
conserve primitiveFieldRef()+