39 const Field<Type2>& pf,
43 if (addr.size() != pf.size())
47 "fvMatrix<Type>::addToInternalField(const unallocLabelList&, "
48 "const Field&, Field&)"
49 ) <<
"sizes of addressing and field are different"
55 intf[addr[faceI]] += pf[faceI];
65 const tmp<Field<Type2> >& tpf,
69 addToInternalField(addr, tpf(), intf);
79 const Field<Type2>& pf,
83 if (addr.size() != pf.size())
87 "fvMatrix<Type>::addToInternalField(const unallocLabelList&, "
88 "const Field&, Field&)"
89 ) <<
"sizes of addressing and field are different"
95 intf[addr[faceI]] -= pf[faceI];
101 template<
class Type2>
105 const tmp<Field<Type2> >& tpf,
109 subtractFromInternalField(addr, tpf(), intf);
121 forAll(internalCoeffs_, patchI)
125 lduAddr().patchAddr(patchI),
126 internalCoeffs_[patchI].
component(solveCmpt),
136 forAll(internalCoeffs_, patchI)
140 lduAddr().patchAddr(patchI),
141 cmptAv(internalCoeffs_[patchI]),
155 forAll(psi_.boundaryField(), patchI)
157 const fvPatchField<Type>& ptf = psi_.boundaryField()[patchI];
158 const Field<Type>& pbc = boundaryCoeffs_[patchI];
162 addToInternalField(lduAddr().patchAddr(patchI), pbc, source);
166 tmp<Field<Type> > tpnf = ptf.patchNeighbourField();
167 const Field<Type>& pnf = tpnf();
173 source[addr[facei]] +=
cmptMultiply(pbc[facei], pnf[facei]);
195 faceFluxCorrectionPtr_(NULL)
199 Info<<
"fvMatrix<Type>(GeometricField<Type, fvPatchField, volMesh>&,"
200 " const dimensionSet&) : "
201 "constructing fvMatrix<Type> for field " << psi_.name()
229 psi_.boundaryField().updateCoeffs();
239 dimensions_(fvm.dimensions_),
240 source_(fvm.source_),
241 internalCoeffs_(fvm.internalCoeffs_),
242 boundaryCoeffs_(fvm.boundaryCoeffs_),
243 faceFluxCorrectionPtr_(NULL)
247 Info<<
"fvMatrix<Type>::fvMatrix(const fvMatrix<Type>&) : "
248 <<
"copying fvMatrix<Type> for field " << psi_.
name()
252 if (fvm.faceFluxCorrectionPtr_)
254 faceFluxCorrectionPtr_ =
new
257 *(fvm.faceFluxCorrectionPtr_)
263 #ifdef ConstructFromTmp
270 const_cast<
fvMatrix<Type>&>(tfvm()),
274 dimensions_(tfvm().dimensions_),
277 const_cast<
fvMatrix<Type>&>(tfvm()).source_,
282 const_cast<
fvMatrix<Type>&>(tfvm()).internalCoeffs_,
287 const_cast<
fvMatrix<Type>&>(tfvm()).boundaryCoeffs_,
290 faceFluxCorrectionPtr_(NULL)
294 Info<<
"fvMatrix<Type>::fvMatrix(const tmp<fvMatrix<Type> >&) : "
295 <<
"copying fvMatrix<Type> for field " << psi_.
name()
299 if (tfvm().faceFluxCorrectionPtr_)
303 faceFluxCorrectionPtr_ = tfvm().faceFluxCorrectionPtr_;
304 tfvm().faceFluxCorrectionPtr_ = NULL;
308 faceFluxCorrectionPtr_ =
new
309 GeometricField<Type, fvsPatchField, surfaceMesh>
311 *(tfvm().faceFluxCorrectionPtr_)
334 faceFluxCorrectionPtr_(NULL)
338 Info<<
"fvMatrix<Type>"
339 "(GeometricField<Type, fvPatchField, volMesh>&, Istream&) : "
340 "constructing fvMatrix<Type> for field " << psi_.name()
376 Info<<
"fvMatrix<Type>::~fvMatrix<Type>() : "
377 <<
"destroying fvMatrix<Type> for field " << psi_.name()
381 if (faceFluxCorrectionPtr_)
383 delete faceFluxCorrectionPtr_;
409 label celli = cellLabels[i];
411 psi_[celli] = values[i];
412 source_[celli] = values[i]*Diag[celli];
414 if (symmetric() || asymmetric())
416 const cell& c = cells[celli];
426 if (celli == own[facei])
428 source_[nei[facei]] -= upper()[facei]*values[i];
432 source_[own[facei]] -= upper()[facei]*values[i];
435 upper()[facei] = 0.0;
439 if (celli == own[facei])
441 source_[nei[facei]] -= lower()[facei]*values[i];
445 source_[own[facei]] -= upper()[facei]*values[i];
448 upper()[facei] = 0.0;
449 lower()[facei] = 0.0;
456 if (internalCoeffs_[patchi].size())
461 internalCoeffs_[
patchi][patchFacei] =
464 boundaryCoeffs_[
patchi][patchFacei] =
479 const bool forceReference
482 if ((forceReference || psi_.needReference()) && celli >= 0)
484 source()[celli] +=
diag()[celli]*value;
506 sumMagOffDiag(sumOff);
509 forAll(psi_.boundaryField(), patchI)
520 const Field<Type>& pCoeffs = boundaryCoeffs_[patchI];
538 Type iCoeff0 = iCoeffs[
face];
540 sumOff[pa[face]] -=
cmptMin(iCoeffs[face]);
541 iCoeffs[face] /= alpha;
556 forAll(psi_.boundaryField(), patchI)
576 S += (D - D0)*psi_.internalField();
583 if (psi_.mesh().relax(psi_.name()))
585 relax(psi_.mesh().relaxationFactor(psi_.name()));
594 GeometricBoundaryField& bFields
599 bFields[patchI].manipulateMatrix(*
this);
608 addCmptAvBoundaryDiag(tdiag());
618 forAll(psi_.boundaryField(), patchI)
626 lduAddr().patchAddr(patchI),
627 internalCoeffs_[patchI],
646 "A("+psi_.name()+
')',
653 dimensions_/psi_.dimensions()/
dimVol,
654 zeroGradientFvPatchScalarField::typeName
658 tAphi().internalField() = D()/psi_.mesh().V();
659 tAphi().correctBoundaryConditions();
675 "H("+psi_.name()+
')',
683 zeroGradientFvPatchScalarField::typeName
689 for (
direction cmpt=0; cmpt<Type::nComponents; cmpt++)
694 addBoundaryDiag(boundaryDiagCmpt, cmpt);
695 boundaryDiagCmpt.
negate();
696 addCmptAvBoundaryDiag(boundaryDiagCmpt);
707 typename Type::labelType validComponents
711 psi_.mesh().solutionD(),
716 for(
direction cmpt=0; cmpt<Type::nComponents; cmpt++)
718 if (validComponents[cmpt] == -1)
748 dimensions_/(
dimVol*psi_.dimensions()),
749 zeroGradientFvPatchScalarField::typeName
773 H1_.internalField() /= psi_.mesh().V();
774 H1_.correctBoundaryConditions();
785 if (!psi_.mesh().fluxRequired(psi_.name()))
788 <<
"flux requested but " << psi_.name()
789 <<
" not specified in the fluxRequired sub-dictionary"
801 "flux("+psi_.name()+
')',
813 for (
direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
824 forAll(InternalContrib, patchI)
826 InternalContrib[patchI] =
829 InternalContrib[patchI],
830 psi_.boundaryField()[patchI].patchInternalField()
836 forAll(NeighbourContrib, patchI)
838 if (psi_.boundaryField()[patchI].coupled())
840 NeighbourContrib[patchI] =
843 NeighbourContrib[patchI],
844 psi_.boundaryField()[patchI].patchNeighbourField()
852 InternalContrib[patchI] - NeighbourContrib[patchI];
855 if (faceFluxCorrectionPtr_)
857 fieldFlux += *faceFluxCorrectionPtr_;
871 FatalErrorIn(
"fvMatrix<Type>::operator=(const fvMatrix<Type>&)")
872 <<
"attempted assignment to self"
876 if (&psi_ != &(fvmv.psi_))
878 FatalErrorIn(
"fvMatrix<Type>::operator=(const fvMatrix<Type>&)")
879 <<
"different fields"
884 source_ = fvmv.source_;
885 internalCoeffs_ = fvmv.internalCoeffs_;
886 boundaryCoeffs_ = fvmv.boundaryCoeffs_;
888 if (faceFluxCorrectionPtr_ && fvmv.faceFluxCorrectionPtr_)
890 *faceFluxCorrectionPtr_ = *fvmv.faceFluxCorrectionPtr_;
892 else if (fvmv.faceFluxCorrectionPtr_)
894 faceFluxCorrectionPtr_ =
896 (*fvmv.faceFluxCorrectionPtr_);
914 internalCoeffs_.negate();
915 boundaryCoeffs_.negate();
917 if (faceFluxCorrectionPtr_)
919 faceFluxCorrectionPtr_->negate();
929 dimensions_ += fvmv.dimensions_;
931 source_ += fvmv.source_;
932 internalCoeffs_ += fvmv.internalCoeffs_;
933 boundaryCoeffs_ += fvmv.boundaryCoeffs_;
935 if (faceFluxCorrectionPtr_ && fvmv.faceFluxCorrectionPtr_)
937 *faceFluxCorrectionPtr_ += *fvmv.faceFluxCorrectionPtr_;
939 else if (fvmv.faceFluxCorrectionPtr_)
941 faceFluxCorrectionPtr_ =
new
944 *fvmv.faceFluxCorrectionPtr_
963 dimensions_ -= fvmv.dimensions_;
965 source_ -= fvmv.source_;
966 internalCoeffs_ -= fvmv.internalCoeffs_;
967 boundaryCoeffs_ -= fvmv.boundaryCoeffs_;
969 if (faceFluxCorrectionPtr_ && fvmv.faceFluxCorrectionPtr_)
971 *faceFluxCorrectionPtr_ -= *fvmv.faceFluxCorrectionPtr_;
973 else if (fvmv.faceFluxCorrectionPtr_)
975 faceFluxCorrectionPtr_ =
977 (-*fvmv.faceFluxCorrectionPtr_);
997 source() -= su.mesh().V()*su.field();
1001 template<
class Type>
1002 void Foam::fvMatrix<Type>::operator+=
1012 template<
class Type>
1013 void Foam::fvMatrix<Type>::operator+=
1023 template<
class Type>
1024 void Foam::fvMatrix<Type>::operator-=
1030 source() += su.mesh().V()*su.field();
1034 template<
class Type>
1035 void Foam::fvMatrix<Type>::operator-=
1045 template<
class Type>
1046 void Foam::fvMatrix<Type>::operator-=
1056 template<
class Type>
1057 void Foam::fvMatrix<Type>::operator+=
1066 template<
class Type>
1067 void Foam::fvMatrix<Type>::operator-=
1076 template<
class Type>
1077 void Foam::fvMatrix<Type>::operator+=
1084 template<
class Type>
1085 void Foam::fvMatrix<Type>::operator-=
1092 template<
class Type>
1093 void Foam::fvMatrix<Type>::operator*=
1100 source_ *= dsf.field();
1102 forAll(boundaryCoeffs_, patchI)
1106 dsf.mesh().boundary()[patchI].patchInternalField(dsf.field())
1109 internalCoeffs_[patchI] *= psf;
1110 boundaryCoeffs_[patchI] *= psf;
1113 if (faceFluxCorrectionPtr_)
1117 "fvMatrix<Type>::operator*="
1118 "(const DimensionedField<scalar, volMesh>&)"
1119 ) <<
"cannot scale a matrix containing a faceFluxCorrection"
1125 template<
class Type>
1126 void Foam::fvMatrix<Type>::operator*=
1136 template<
class Type>
1137 void Foam::fvMatrix<Type>::operator*=
1147 template<
class Type>
1148 void Foam::fvMatrix<Type>::operator*=
1153 dimensions_ *= ds.dimensions();
1155 source_ *= ds.value();
1156 internalCoeffs_ *= ds.value();
1157 boundaryCoeffs_ *= ds.value();
1159 if (faceFluxCorrectionPtr_)
1161 *faceFluxCorrectionPtr_ *= ds.value();
1168 template<
class Type>
1176 if (&fvm1.
psi() != &fvm2.
psi())
1180 "checkMethod(const fvMatrix<Type>&, const fvMatrix<Type>&)"
1181 ) <<
"incompatible fields for operation "
1183 <<
"[" << fvm1.
psi().
name() <<
"] "
1185 <<
" [" << fvm2.
psi().
name() <<
"]"
1193 "checkMethod(const fvMatrix<Type>&, const fvMatrix<Type>&)"
1194 ) <<
"incompatible dimensions for operation "
1204 template<
class Type>
1216 "checkMethod(const fvMatrix<Type>&, const GeometricField<Type, "
1217 "fvPatchField, volMesh>&)"
1218 ) <<
"incompatible dimensions for operation "
1228 template<
class Type>
1240 "checkMethod(const fvMatrix<Type>&, const dimensioned<Type>&)"
1241 ) <<
"incompatible dimensions for operation "
1251 template<
class Type>
1258 return fvm.
solve(solverControls);
1261 template<
class Type>
1277 template<
class Type>
1283 template<
class Type>
1295 template<
class Type>
1309 tAcorr().faceFluxCorrectionPtr() = (-A.
flux()).ptr();
1316 template<
class Type>
1329 (A.hasUpper() || A.hasLower())
1330 && A.psi().mesh().fluxRequired(A.psi().name())
1333 tAcorr().faceFluxCorrectionPtr() = (-A.flux()).ptr();
1342 template<
class Type>
1353 template<
class Type>
1364 template<
class Type>
1375 template<
class Type>
1386 template<
class Type>
1395 tC().source() += su.mesh().V()*su.field();
1399 template<
class Type>
1408 tC().source() += tsu().mesh().V()*tsu().field();
1413 template<
class Type>
1422 tC().source() += tsu().mesh().V()*tsu().internalField();
1427 template<
class Type>
1436 tC().source() += su.mesh().V()*su.field();
1440 template<
class Type>
1449 tC().source() += tsu().mesh().V()*tsu().field();
1454 template<
class Type>
1463 tC().source() += tsu().mesh().V()*tsu().internalField();
1468 template<
class Type>
1477 tC().source() += A.
psi().
mesh().
V()*su.value();
1481 template<
class Type>
1490 tC().source() += tC().psi().mesh().V()*su.value();
1494 template<
class Type>
1505 template<
class Type>
1516 template<
class Type>
1527 template<
class Type>
1539 template<
class Type>
1552 template<
class Type>
1565 template<
class Type>
1578 template<
class Type>
1592 template<
class Type>
1601 tC().source() -= su.mesh().V()*su.field();
1605 template<
class Type>
1614 tC().source() -= tsu().mesh().V()*tsu().field();
1619 template<
class Type>
1628 tC().source() -= tsu().mesh().V()*tsu().internalField();
1633 template<
class Type>
1642 tC().source() -= su.mesh().V()*su.field();
1646 template<
class Type>
1655 tC().source() -= tsu().mesh().V()*tsu().field();
1660 template<
class Type>
1669 tC().source() -= tsu().mesh().V()*tsu().internalField();
1674 template<
class Type>
1687 template<
class Type>
1696 tC().source() -= tsu().mesh().V()*tsu().field();
1701 template<
class Type>
1710 tC().source() -= tsu().mesh().V()*tsu().internalField();
1715 template<
class Type>
1728 template<
class Type>
1737 tC().source() -= tsu().mesh().V()*tsu().field();
1742 template<
class Type>
1751 tC().source() -= tsu().mesh().V()*tsu().internalField();
1757 template<
class Type>
1770 template<
class Type>
1783 template<
class Type>
1797 template<
class Type>
1811 template<
class Type>
1824 template<
class Type>
1833 tC().source() += tsu().mesh().V()*tsu().field();
1838 template<
class Type>
1847 tC().source() += tsu().mesh().V()*tsu().internalField();
1852 template<
class Type>
1865 template<
class Type>
1874 tC().source() += tsu().mesh().V()*tsu().field();
1879 template<
class Type>
1888 tC().source() += tsu().mesh().V()*tsu().internalField();
1893 template<
class Type>
1907 template<
class Type>
1917 tC().source() -= tsu().mesh().V()*tsu().field();
1922 template<
class Type>
1932 tC().source() -= tsu().mesh().V()*tsu().internalField();
1937 template<
class Type>
1951 template<
class Type>
1961 tC().source() -= tsu().mesh().V()*tsu().field();
1966 template<
class Type>
1976 tC().source() -= tsu().mesh().V()*tsu().internalField();
1981 template<
class Type>
1990 tC().source() -= su.value()*A.
psi().
mesh().
V();
1994 template<
class Type>
2003 tC().source() -= su.value()*tC().psi().
mesh().
V();
2007 template<
class Type>
2020 template<
class Type>
2029 tC().source() -= su.
value()*tC().psi().mesh().V();
2033 template<
class Type>
2042 tC().source() += su.
value()*tC().psi().mesh().V();
2046 template<
class Type>
2055 tC().source() += su.
value()*tC().psi().mesh().V();
2059 template<
class Type>
2073 template<
class Type>
2083 tC().source() -= su.
value()*tC().psi().mesh().V();
2088 template<
class Type>
2100 template<
class Type>
2112 template<
class Type>
2124 template<
class Type>
2136 template<
class Type>
2148 template<
class Type>
2160 template<
class Type>
2172 template<
class Type>
2185 template<
class Type>
2207 zeroGradientFvPatchScalarField::typeName
2213 for (
direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
2217 M.addBoundaryDiag(boundaryDiagCmpt, cmpt);
2230 template<
class Type>
2243 template<
class Type>
2256 template<
class Type>
2269 template<
class Type>
2283 template<
class Type>
2300 template<
class Type>
2303 os << static_cast<const lduMatrix&>(fvm) <<
nl
2304 << fvm.dimensions_ <<
nl
2305 << fvm.source_ <<
nl
2306 << fvm.internalCoeffs_ <<
nl
2307 << fvm.boundaryCoeffs_ <<
endl;
2309 os.
check(
"Ostream& operator<<(Ostream&, fvMatrix<Type>&");