48 if (addr.
size() != pf.size())
51 <<
"addressing (" << addr.
size()
52 <<
") and field (" << pf.size() <<
") are different sizes" <<
endl
58 intf[addr[faceI]] += pf[faceI];
72 addToInternalField(addr, tpf(), intf);
86 if (addr.
size() != pf.size())
89 <<
"addressing (" << addr.
size()
90 <<
") and field (" << pf.size() <<
") are different sizes" <<
endl
96 intf[addr[faceI]] -= pf[faceI];
102 template<
class Type2>
110 subtractFromInternalField(addr, tpf(), intf);
122 forAll(internalCoeffs_, patchI)
126 lduAddr().patchAddr(patchI),
127 internalCoeffs_[patchI].
component(solveCmpt),
137 forAll(internalCoeffs_, patchI)
141 lduAddr().patchAddr(patchI),
142 cmptAv(internalCoeffs_[patchI]),
156 forAll(psi_.boundaryField(), patchI)
163 addToInternalField(lduAddr().patchAddr(patchI), pbc, source);
170 const labelUList& addr = lduAddr().patchAddr(patchI);
174 source[addr[facei]] +=
cmptMultiply(pbc[facei], pnf[facei]);
194 internalCoeffs_(
psi.mesh().boundary().size()),
195 boundaryCoeffs_(
psi.mesh().boundary().size()),
196 faceFluxCorrectionPtr_(
nullptr)
199 <<
"constructing faMatrix<Type> for field " << psi_.name()
210 psi.mesh().boundary()[patchI].size(),
220 psi.mesh().boundary()[patchI].size(),
230 const label currentStatePsi = psiRef.eventNo();
232 psiRef.eventNo() = currentStatePsi;
242 dimensions_(
fam.dimensions_),
243 source_(
fam.source_),
244 internalCoeffs_(
fam.internalCoeffs_),
245 boundaryCoeffs_(
fam.boundaryCoeffs_),
246 faceFluxCorrectionPtr_(nullptr)
249 <<
"copying faMatrix<Type> for field " << psi_.name()
252 if (
fam.faceFluxCorrectionPtr_)
254 faceFluxCorrectionPtr_ =
new
257 *(
fam.faceFluxCorrectionPtr_)
274 internalCoeffs_(
psi.mesh().boundary().size()),
275 boundaryCoeffs_(
psi.mesh().boundary().size()),
276 faceFluxCorrectionPtr_(
nullptr)
279 <<
"constructing faMatrix<Type> for field " << psi_.name()
290 psi.mesh().boundary()[patchI].size(),
300 psi.mesh().boundary()[patchI].size(),
322 <<
"Destroying faMatrix<Type> for field " << psi_.name() <<
endl;
331 template<
template<
class>
class ListType>
335 const ListType<Type>&
values
340 for (
const label i : faceLabels)
342 this->eliminatedEqns().insert(i);
361 const label facei = faceLabels[i];
362 const Type& value =
values[i];
365 source_[facei] = value*Diag[facei];
367 if (symmetric() || asymmetric())
369 for (
const label edgei : edges[facei])
371 if (
mesh.isInternalEdge(edgei))
375 if (facei == own[edgei])
377 source_[nei[edgei]] -=
upper()[edgei]*value;
381 source_[own[edgei]] -=
upper()[edgei]*value;
384 upper()[edgei] = 0.0;
388 if (facei == own[edgei])
390 source_[nei[edgei]] -=
lower()[edgei]*value;
394 source_[own[edgei]] -=
upper()[edgei]*value;
397 upper()[edgei] = 0.0;
398 lower()[edgei] = 0.0;
405 if (internalCoeffs_[patchi].size())
407 const label patchEdgei =
410 internalCoeffs_[patchi][patchEdgei] =
Zero;
411 boundaryCoeffs_[patchi][patchEdgei] =
Zero;
438 this->setValuesFromList(faceLabels,
values);
449 this->setValuesFromList(faceLabels,
values);
458 const bool forceReference
461 if ((forceReference || psi_.needReference()) && facei >= 0)
465 source()[facei] +=
diag()[facei]*value;
477 const bool forceReference
480 if (forceReference || psi_.needReference())
484 const label
faceId = faceLabels[facei];
500 const bool forceReference
503 if (forceReference || psi_.needReference())
507 const label
faceId = faceLabels[facei];
534 sumMagOffDiag(sumOff);
537 forAll(psi_.boundaryField(), patchI)
543 const labelUList& pa = lduAddr().patchAddr(patchI);
548 const Field<Type>& pCoeffs = boundaryCoeffs_[patchI];
566 Type iCoeff0 = iCoeffs[
face];
584 forAll(psi_.boundaryField(), patchI)
590 const labelUList& pa = lduAddr().patchAddr(patchI);
604 S += (
D - D0)*psi_.internalField();
611 if (psi_.mesh().relaxEquation(psi_.name()))
613 relax(psi_.mesh().equationRelaxationFactor(psi_.name()));
618 <<
"Relaxation factor for field " << psi_.name()
619 <<
" not found. Relaxation will not be used." <<
endl;
628 addCmptAvBoundaryDiag(tdiag.
ref());
642 "A("+psi_.name()+
')',
647 dimensions_/psi_.dimensions()/
dimArea,
648 zeroGradientFaPatchScalarField::typeName
669 "H("+psi_.name()+
')',
675 zeroGradientFaPatchScalarField::typeName
681 for (
direction cmpt=0; cmpt<Type::nComponents; ++cmpt)
683 scalarField psiCmpt(psi_.primitiveField().component(cmpt));
686 addBoundaryDiag(boundaryDiagCmpt, cmpt);
687 boundaryDiagCmpt.negate();
688 addCmptAvBoundaryDiag(boundaryDiagCmpt);
707 if (!psi_.mesh().fluxRequired(psi_.name()))
710 <<
"flux requested but " << psi_.name()
711 <<
" not specified in the fluxRequired sub-dictionary of faSchemes"
722 "flux("+psi_.name()+
')',
733 for (
direction cmpt=0; cmpt<pTraits<Type>::nComponents; ++cmpt)
744 forAll(InternalContrib, patchI)
746 InternalContrib[patchI] =
749 InternalContrib[patchI],
750 psi_.boundaryField()[patchI].patchInternalField()
756 forAll(NeighbourContrib, patchI)
758 if (psi_.boundaryField()[patchI].coupled())
760 NeighbourContrib[patchI] =
763 NeighbourContrib[patchI],
764 psi_.boundaryField()[patchI].patchNeighbourField()
772 InternalContrib[patchI] - NeighbourContrib[patchI];
775 if (faceFluxCorrectionPtr_)
777 fieldFlux += *faceFluxCorrectionPtr_;
794 if (&psi_ != &(famv.psi_))
797 <<
"different fields"
802 source_ = famv.source_;
803 internalCoeffs_ = famv.internalCoeffs_;
804 boundaryCoeffs_ = famv.boundaryCoeffs_;
806 if (faceFluxCorrectionPtr_ && famv.faceFluxCorrectionPtr_)
808 *faceFluxCorrectionPtr_ = *famv.faceFluxCorrectionPtr_;
810 else if (famv.faceFluxCorrectionPtr_)
812 faceFluxCorrectionPtr_ =
814 (*famv.faceFluxCorrectionPtr_);
832 internalCoeffs_.negate();
833 boundaryCoeffs_.negate();
835 if (faceFluxCorrectionPtr_)
837 faceFluxCorrectionPtr_->negate();
847 dimensions_ += famv.dimensions_;
849 source_ += famv.source_;
850 internalCoeffs_ += famv.internalCoeffs_;
851 boundaryCoeffs_ += famv.boundaryCoeffs_;
853 if (faceFluxCorrectionPtr_ && famv.faceFluxCorrectionPtr_)
855 *faceFluxCorrectionPtr_ += *famv.faceFluxCorrectionPtr_;
857 else if (famv.faceFluxCorrectionPtr_)
859 faceFluxCorrectionPtr_ =
new
862 *famv.faceFluxCorrectionPtr_
881 dimensions_ -= famv.dimensions_;
883 source_ -= famv.source_;
884 internalCoeffs_ -= famv.internalCoeffs_;
885 boundaryCoeffs_ -= famv.boundaryCoeffs_;
887 if (faceFluxCorrectionPtr_ && famv.faceFluxCorrectionPtr_)
889 *faceFluxCorrectionPtr_ -= *famv.faceFluxCorrectionPtr_;
891 else if (famv.faceFluxCorrectionPtr_)
893 faceFluxCorrectionPtr_ =
895 (-*famv.faceFluxCorrectionPtr_);
915 source() -= su.mesh().S()*su.internalField();
937 source() += su.mesh().S()*su.internalField();
958 source() -= su.mesh().S()*su;
968 source() += su.mesh().S()*su;
978 dimensions_ *= vsf.dimensions();
980 source_ *= vsf.internalField();
982 forAll(boundaryCoeffs_, patchI)
984 const scalarField psf(vsf.boundaryField()[patchI].patchInternalField());
985 internalCoeffs_[patchI] *= psf;
986 boundaryCoeffs_[patchI] *= psf;
989 if (faceFluxCorrectionPtr_)
992 <<
"cannot scale a matrix containing a faceFluxCorrection"
1009 template<
class Type>
1015 dimensions_ *= ds.dimensions();
1017 source_ *= ds.value();
1018 internalCoeffs_ *= ds.value();
1019 boundaryCoeffs_ *= ds.value();
1021 if (faceFluxCorrectionPtr_)
1023 *faceFluxCorrectionPtr_ *= ds.value();
1030 template<
class Type>
1038 if (&fam1.
psi() != &fam2.
psi())
1041 <<
"incompatible fields for operation "
1043 <<
"[" << fam1.
psi().name() <<
"] "
1045 <<
" [" << fam2.
psi().name() <<
"]"
1056 <<
"incompatible dimensions for operation "
1066 template<
class Type>
1077 &&
fam.dimensions()/
dimArea != vf.dimensions()
1081 <<
"incompatible dimensions for operation "
1083 <<
"[" <<
fam.psi().name() <<
fam.dimensions()/
dimArea <<
" ] "
1085 <<
" [" << vf.name() << vf.dimensions() <<
" ]"
1091 template<
class Type>
1106 <<
"incompatible dimensions for operation "
1108 <<
"[" <<
fam.psi().name() <<
fam.dimensions()/
dimArea <<
" ] "
1116 template<
class Type>
1119 faMatrix<Type>&
fam,
1120 Istream& solverControls
1123 return fam.solve(solverControls);
1127 template<
class Type>
1130 const tmp<faMatrix<Type>>& tfam,
1131 Istream& solverControls
1134 SolverPerformance<Type> solverPerf =
1135 const_cast<faMatrix<Type>&
>(tfam()).
solve(solverControls);
1142 template<
class Type>
1149 template<
class Type>
1152 SolverPerformance<Type> solverPerf =
1153 const_cast<faMatrix<Type>&
>(tfam()).
solve();
1162 template<
class Type>
1165 const faMatrix<Type>&
A,
1166 const faMatrix<Type>&
B
1170 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1176 template<
class Type>
1179 const tmp<faMatrix<Type>>& tA,
1180 const faMatrix<Type>&
B
1184 tmp<faMatrix<Type>> tC(tA.
ptr());
1190 template<
class Type>
1193 const faMatrix<Type>&
A,
1194 const tmp<faMatrix<Type>>& tB
1198 tmp<faMatrix<Type>> tC(tB.ptr());
1204 template<
class Type>
1207 const tmp<faMatrix<Type>>& tA,
1208 const tmp<faMatrix<Type>>& tB
1212 tmp<faMatrix<Type>> tC(tA.
ptr());
1219 template<
class Type>
1222 const faMatrix<Type>&
A
1225 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1231 template<
class Type>
1234 const tmp<faMatrix<Type>>& tA
1237 tmp<faMatrix<Type>> tC(tA.
ptr());
1243 template<
class Type>
1246 const faMatrix<Type>&
A,
1247 const faMatrix<Type>&
B
1251 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1257 template<
class Type>
1260 const tmp<faMatrix<Type>>& tA,
1261 const faMatrix<Type>&
B
1265 tmp<faMatrix<Type>> tC(tA.
ptr());
1271 template<
class Type>
1274 const faMatrix<Type>&
A,
1275 const tmp<faMatrix<Type>>& tB
1279 tmp<faMatrix<Type>> tC(tB.ptr());
1286 template<
class Type>
1289 const tmp<faMatrix<Type>>& tA,
1290 const tmp<faMatrix<Type>>& tB
1294 tmp<faMatrix<Type>> tC(tA.
ptr());
1301 template<
class Type>
1304 const faMatrix<Type>&
A,
1305 const faMatrix<Type>&
B
1313 template<
class Type>
1316 const tmp<faMatrix<Type>>& tA,
1317 const faMatrix<Type>&
B
1325 template<
class Type>
1328 const faMatrix<Type>&
A,
1329 const tmp<faMatrix<Type>>& tB
1337 template<
class Type>
1340 const tmp<faMatrix<Type>>& tA,
1341 const tmp<faMatrix<Type>>& tB
1349 template<
class Type>
1352 const faMatrix<Type>&
A,
1353 const GeometricField<Type, faPatchField, areaMesh>& su
1357 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1358 tC.ref().source() -= su.mesh().S()*su.internalField();
1363 template<
class Type>
1366 const tmp<faMatrix<Type>>& tA,
1367 const GeometricField<Type, faPatchField, areaMesh>& su
1371 tmp<faMatrix<Type>> tC(tA.
ptr());
1372 tC.ref().source() -= su.mesh().S()*su.internalField();
1377 template<
class Type>
1380 const faMatrix<Type>&
A,
1381 const tmp<GeometricField<Type, faPatchField, areaMesh>>& tsu
1385 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1386 tC.ref().source() -= tsu().mesh().S()*tsu().internalField();
1392 template<
class Type>
1395 const tmp<faMatrix<Type>>& tA,
1396 const tmp<GeometricField<Type, faPatchField, areaMesh>>& tsu
1400 tmp<faMatrix<Type>> tC(tA.
ptr());
1401 tC.ref().source() -= tsu().mesh().S()*tsu().internalField();
1407 template<
class Type>
1410 const GeometricField<Type, faPatchField, areaMesh>& su,
1411 const faMatrix<Type>&
A
1415 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1416 tC.ref().source() -= su.mesh().S()*su.internalField();
1421 template<
class Type>
1424 const GeometricField<Type, faPatchField, areaMesh>& su,
1425 const tmp<faMatrix<Type>>& tA
1429 tmp<faMatrix<Type>> tC(tA.
ptr());
1430 tC.ref().source() -= su.mesh().S()*su.internalField();
1435 template<
class Type>
1438 const tmp<GeometricField<Type, faPatchField, areaMesh>>& tsu,
1439 const faMatrix<Type>&
A
1443 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1444 tC.ref().source() -= tsu().mesh().S()*tsu().internalField();
1450 template<
class Type>
1453 const tmp<GeometricField<Type, faPatchField, areaMesh>>& tsu,
1454 const tmp<faMatrix<Type>>& tA
1458 tmp<faMatrix<Type>> tC(tA.
ptr());
1459 tC.ref().source() -= tsu().mesh().S()*tsu().internalField();
1465 template<
class Type>
1468 const faMatrix<Type>&
A,
1469 const GeometricField<Type, faPatchField, areaMesh>& su
1473 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1474 tC.ref().source() += su.mesh().S()*su.internalField();
1479 template<
class Type>
1482 const tmp<faMatrix<Type>>& tA,
1483 const GeometricField<Type, faPatchField, areaMesh>& su
1487 tmp<faMatrix<Type>> tC(tA.
ptr());
1488 tC.ref().source() += su.mesh().S()*su.internalField();
1493 template<
class Type>
1496 const faMatrix<Type>&
A,
1497 const tmp<GeometricField<Type, faPatchField, areaMesh>>& tsu
1501 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1502 tC.ref().source() += tsu().mesh().S()*tsu().internalField();
1508 template<
class Type>
1511 const tmp<faMatrix<Type>>& tA,
1512 const tmp<GeometricField<Type, faPatchField, areaMesh>>& tsu
1516 tmp<faMatrix<Type>> tC(tA.
ptr());
1517 tC.ref().source() += tsu().mesh().S()*tsu().internalField();
1523 template<
class Type>
1526 const GeometricField<Type, faPatchField, areaMesh>& su,
1527 const faMatrix<Type>&
A
1531 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1533 tC.ref().source() -= su.mesh().S()*su.internalField();
1538 template<
class Type>
1541 const GeometricField<Type, faPatchField, areaMesh>& su,
1542 const tmp<faMatrix<Type>>& tA
1546 tmp<faMatrix<Type>> tC(tA.
ptr());
1548 tC.ref().source() -= su.mesh().S()*su.internalField();
1553 template<
class Type>
1556 const tmp<GeometricField<Type, faPatchField, areaMesh>>& tsu,
1557 const faMatrix<Type>&
A
1561 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1563 tC.ref().source() -= tsu().mesh().S()*tsu().internalField();
1569 template<
class Type>
1572 const tmp<GeometricField<Type, faPatchField, areaMesh>>& tsu,
1573 const tmp<faMatrix<Type>>& tA
1577 tmp<faMatrix<Type>> tC(tA.
ptr());
1579 tC.ref().source() -= tsu().mesh().S()*tsu().internalField();
1585 template<
class Type>
1588 const faMatrix<Type>&
A,
1589 const dimensioned<Type>& su
1593 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1594 tC.ref().source() -= su.value()*
A.psi().mesh().S();
1599 template<
class Type>
1602 const tmp<faMatrix<Type>>& tA,
1603 const dimensioned<Type>& su
1607 tmp<faMatrix<Type>> tC(tA.
ptr());
1608 tC.ref().source() -= su.value()*tC().psi().mesh().S();
1613 template<
class Type>
1616 const dimensioned<Type>& su,
1617 const faMatrix<Type>&
A
1621 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1622 tC.ref().source() -= su.value()*
A.psi().mesh().S();
1627 template<
class Type>
1630 const dimensioned<Type>& su,
1631 const tmp<faMatrix<Type>>& tA
1635 tmp<faMatrix<Type>> tC(tA.
ptr());
1636 tC.ref().source() -= su.value()*tC().psi().mesh().S();
1641 template<
class Type>
1644 const faMatrix<Type>&
A,
1645 const dimensioned<Type>& su
1649 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1650 tC.ref().source() += su.value()*tC().psi().mesh().S();
1655 template<
class Type>
1658 const tmp<faMatrix<Type>>& tA,
1659 const dimensioned<Type>& su
1663 tmp<faMatrix<Type>> tC(tA.
ptr());
1664 tC.ref().source() += su.value()*tC().psi().mesh().S();
1669 template<
class Type>
1672 const dimensioned<Type>& su,
1673 const faMatrix<Type>&
A
1677 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1679 tC.ref().source() -= su.value()*
A.psi().mesh().S();
1684 template<
class Type>
1687 const dimensioned<Type>& su,
1688 const tmp<faMatrix<Type>>& tA
1692 tmp<faMatrix<Type>> tC(tA.
ptr());
1694 tC.ref().source() -= su.value()*tC().psi().mesh().S();
1699 template<
class Type>
1702 const faMatrix<Type>&
A,
1703 const GeometricField<Type, faPatchField, areaMesh>& su
1707 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1708 tC.ref().source() += su.mesh().S()*su.internalField();
1713 template<
class Type>
1716 const tmp<faMatrix<Type>>& tA,
1717 const GeometricField<Type, faPatchField, areaMesh>& su
1721 tmp<faMatrix<Type>> tC(tA.ptr());
1722 tC.ref().source() += su.mesh().S()*su.internalField();
1727 template<
class Type>
1730 const faMatrix<Type>&
A,
1731 const tmp<GeometricField<Type, faPatchField, areaMesh>>& tsu
1735 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1736 tC.ref().source() += tsu().mesh().S()*tsu().internalField();
1742 template<
class Type>
1745 const tmp<faMatrix<Type>>& tA,
1746 const tmp<GeometricField<Type, faPatchField, areaMesh>>& tsu
1750 tmp<faMatrix<Type>> tC(tA.ptr());
1751 tC.ref().source() += tsu().mesh().S()*tsu().internalField();
1757 template<
class Type>
1760 const faMatrix<Type>&
A,
1761 const dimensioned<Type>& su
1765 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1766 tC.ref().source() +=
A.psi().mesh().S()*su.value();
1771 template<
class Type>
1774 const tmp<faMatrix<Type>>& tA,
1775 const dimensioned<Type>& su
1779 tmp<faMatrix<Type>> tC(tA.ptr());
1780 tC.ref().source() += tC().psi().mesh().S()*su.value();
1785 template<
class Type>
1789 const faMatrix<Type>&
A
1792 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1798 template<
class Type>
1801 const tmp<areaScalarField>& tvsf,
1802 const faMatrix<Type>&
A
1805 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1811 template<
class Type>
1815 const tmp<faMatrix<Type>>& tA
1818 tmp<faMatrix<Type>> tC(tA.ptr());
1824 template<
class Type>
1827 const tmp<areaScalarField>& tvsf,
1828 const tmp<faMatrix<Type>>& tA
1831 tmp<faMatrix<Type>> tC(tA.ptr());
1837 template<
class Type>
1840 const dimensioned<scalar>& ds,
1841 const faMatrix<Type>&
A
1844 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1850 template<
class Type>
1853 const dimensioned<scalar>& ds,
1854 const tmp<faMatrix<Type>>& tA
1857 tmp<faMatrix<Type>> tC(tA.ptr());
1865 template<
class Type>
1868 os << static_cast<const lduMatrix&>(
fam) <<
nl
1869 <<
fam.dimensions_ <<
nl
1870 <<
fam.source_ <<
nl
1871 <<
fam.internalCoeffs_ <<
nl
1872 <<
fam.boundaryCoeffs_ <<
endl;