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(),
325 <<
"Destroying faMatrix<Type> for field " << psi_.name() <<
endl;
334 template<
template<
class>
class ListType>
338 const ListType<Type>&
values
343 for (
const label i : faceLabels)
345 this->eliminatedEqns().insert(i);
364 const label facei = faceLabels[i];
365 const Type& value =
values[i];
368 source_[facei] = value*Diag[facei];
370 if (symmetric() || asymmetric())
372 for (
const label edgei : edges[facei])
374 if (
mesh.isInternalEdge(edgei))
378 if (facei == own[edgei])
380 source_[nei[edgei]] -=
upper()[edgei]*value;
384 source_[own[edgei]] -=
upper()[edgei]*value;
387 upper()[edgei] = 0.0;
391 if (facei == own[edgei])
393 source_[nei[edgei]] -=
lower()[edgei]*value;
397 source_[own[edgei]] -=
upper()[edgei]*value;
400 upper()[edgei] = 0.0;
401 lower()[edgei] = 0.0;
408 if (internalCoeffs_[patchi].size())
410 const label patchEdgei =
413 internalCoeffs_[patchi][patchEdgei] =
Zero;
414 boundaryCoeffs_[patchi][patchEdgei] =
Zero;
441 this->setValuesFromList(faceLabels,
values);
452 this->setValuesFromList(faceLabels,
values);
461 const bool forceReference
464 if ((forceReference || psi_.needReference()) && facei >= 0)
468 source()[facei] +=
diag()[facei]*value;
480 const bool forceReference
483 if (forceReference || psi_.needReference())
487 const label
faceId = faceLabels[facei];
503 const bool forceReference
506 if (forceReference || psi_.needReference())
510 const label
faceId = faceLabels[facei];
537 sumMagOffDiag(sumOff);
540 forAll(psi_.boundaryField(), patchI)
546 const labelUList& pa = lduAddr().patchAddr(patchI);
551 const Field<Type>& pCoeffs = boundaryCoeffs_[patchI];
569 Type iCoeff0 = iCoeffs[
face];
587 forAll(psi_.boundaryField(), patchI)
593 const labelUList& pa = lduAddr().patchAddr(patchI);
607 S += (
D - D0)*psi_.internalField();
614 if (psi_.mesh().relaxEquation(psi_.name()))
616 relax(psi_.mesh().equationRelaxationFactor(psi_.name()));
621 <<
"Relaxation factor for field " << psi_.name()
622 <<
" not found. Relaxation will not be used." <<
endl;
631 addCmptAvBoundaryDiag(tdiag.
ref());
645 "A("+psi_.name()+
')',
650 dimensions_/psi_.dimensions()/
dimArea,
651 zeroGradientFaPatchScalarField::typeName
672 "H("+psi_.name()+
')',
678 zeroGradientFaPatchScalarField::typeName
684 for (
direction cmpt=0; cmpt<Type::nComponents; ++cmpt)
686 scalarField psiCmpt(psi_.primitiveField().component(cmpt));
689 addBoundaryDiag(boundaryDiagCmpt, cmpt);
690 boundaryDiagCmpt.negate();
691 addCmptAvBoundaryDiag(boundaryDiagCmpt);
710 if (!psi_.mesh().fluxRequired(psi_.name()))
713 <<
"flux requested but " << psi_.name()
714 <<
" not specified in the fluxRequired sub-dictionary of faSchemes"
725 "flux("+psi_.name()+
')',
736 for (
direction cmpt=0; cmpt<pTraits<Type>::nComponents; ++cmpt)
747 forAll(InternalContrib, patchI)
749 InternalContrib[patchI] =
752 InternalContrib[patchI],
753 psi_.boundaryField()[patchI].patchInternalField()
759 forAll(NeighbourContrib, patchI)
761 if (psi_.boundaryField()[patchI].coupled())
763 NeighbourContrib[patchI] =
766 NeighbourContrib[patchI],
767 psi_.boundaryField()[patchI].patchNeighbourField()
775 InternalContrib[patchI] - NeighbourContrib[patchI];
778 if (faceFluxCorrectionPtr_)
780 fieldFlux += *faceFluxCorrectionPtr_;
797 if (&psi_ != &(famv.psi_))
800 <<
"different fields"
805 source_ = famv.source_;
806 internalCoeffs_ = famv.internalCoeffs_;
807 boundaryCoeffs_ = famv.boundaryCoeffs_;
809 if (faceFluxCorrectionPtr_ && famv.faceFluxCorrectionPtr_)
811 *faceFluxCorrectionPtr_ = *famv.faceFluxCorrectionPtr_;
813 else if (famv.faceFluxCorrectionPtr_)
815 faceFluxCorrectionPtr_ =
817 (*famv.faceFluxCorrectionPtr_);
835 internalCoeffs_.negate();
836 boundaryCoeffs_.negate();
838 if (faceFluxCorrectionPtr_)
840 faceFluxCorrectionPtr_->negate();
850 dimensions_ += famv.dimensions_;
852 source_ += famv.source_;
853 internalCoeffs_ += famv.internalCoeffs_;
854 boundaryCoeffs_ += famv.boundaryCoeffs_;
856 if (faceFluxCorrectionPtr_ && famv.faceFluxCorrectionPtr_)
858 *faceFluxCorrectionPtr_ += *famv.faceFluxCorrectionPtr_;
860 else if (famv.faceFluxCorrectionPtr_)
862 faceFluxCorrectionPtr_ =
new
865 *famv.faceFluxCorrectionPtr_
884 dimensions_ -= famv.dimensions_;
886 source_ -= famv.source_;
887 internalCoeffs_ -= famv.internalCoeffs_;
888 boundaryCoeffs_ -= famv.boundaryCoeffs_;
890 if (faceFluxCorrectionPtr_ && famv.faceFluxCorrectionPtr_)
892 *faceFluxCorrectionPtr_ -= *famv.faceFluxCorrectionPtr_;
894 else if (famv.faceFluxCorrectionPtr_)
896 faceFluxCorrectionPtr_ =
898 (-*famv.faceFluxCorrectionPtr_);
918 source() -= su.mesh().S()*su.internalField();
940 source() += su.mesh().S()*su.internalField();
961 source() -= su.mesh().S()*su;
971 source() += su.mesh().S()*su;
981 dimensions_ *= vsf.dimensions();
983 source_ *= vsf.internalField();
985 forAll(boundaryCoeffs_, patchI)
987 const scalarField psf(vsf.boundaryField()[patchI].patchInternalField());
988 internalCoeffs_[patchI] *= psf;
989 boundaryCoeffs_[patchI] *= psf;
992 if (faceFluxCorrectionPtr_)
995 <<
"cannot scale a matrix containing a faceFluxCorrection"
1001 template<
class Type>
1012 template<
class Type>
1018 dimensions_ *= ds.dimensions();
1020 source_ *= ds.value();
1021 internalCoeffs_ *= ds.value();
1022 boundaryCoeffs_ *= ds.value();
1024 if (faceFluxCorrectionPtr_)
1026 *faceFluxCorrectionPtr_ *= ds.value();
1033 template<
class Type>
1041 if (&fam1.
psi() != &fam2.
psi())
1044 <<
"incompatible fields for operation "
1046 <<
"[" << fam1.
psi().name() <<
"] "
1048 <<
" [" << fam2.
psi().name() <<
"]"
1055 <<
"incompatible dimensions for operation "
1065 template<
class Type>
1076 <<
"incompatible dimensions for operation "
1078 <<
"[" <<
fam.psi().name() <<
fam.dimensions()/
dimArea <<
" ] "
1080 <<
" [" << vf.name() << vf.dimensions() <<
" ]"
1086 template<
class Type>
1097 <<
"incompatible dimensions for operation "
1099 <<
"[" <<
fam.psi().name() <<
fam.dimensions()/
dimArea <<
" ] "
1107 template<
class Type>
1110 faMatrix<Type>&
fam,
1111 Istream& solverControls
1114 return fam.solve(solverControls);
1118 template<
class Type>
1121 const tmp<faMatrix<Type>>& tfam,
1122 Istream& solverControls
1125 SolverPerformance<Type> solverPerf =
1126 const_cast<faMatrix<Type>&
>(tfam()).
solve(solverControls);
1133 template<
class Type>
1140 template<
class Type>
1143 SolverPerformance<Type> solverPerf =
1144 const_cast<faMatrix<Type>&
>(tfam()).
solve();
1153 template<
class Type>
1156 const faMatrix<Type>&
A,
1157 const faMatrix<Type>&
B
1161 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1167 template<
class Type>
1170 const tmp<faMatrix<Type>>& tA,
1171 const faMatrix<Type>&
B
1175 tmp<faMatrix<Type>> tC(tA.
ptr());
1181 template<
class Type>
1184 const faMatrix<Type>&
A,
1185 const tmp<faMatrix<Type>>& tB
1189 tmp<faMatrix<Type>> tC(tB.ptr());
1195 template<
class Type>
1198 const tmp<faMatrix<Type>>& tA,
1199 const tmp<faMatrix<Type>>& tB
1203 tmp<faMatrix<Type>> tC(tA.
ptr());
1210 template<
class Type>
1213 const faMatrix<Type>&
A
1216 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1222 template<
class Type>
1225 const tmp<faMatrix<Type>>& tA
1228 tmp<faMatrix<Type>> tC(tA.
ptr());
1234 template<
class Type>
1237 const faMatrix<Type>&
A,
1238 const faMatrix<Type>&
B
1242 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1248 template<
class Type>
1251 const tmp<faMatrix<Type>>& tA,
1252 const faMatrix<Type>&
B
1256 tmp<faMatrix<Type>> tC(tA.
ptr());
1262 template<
class Type>
1265 const faMatrix<Type>&
A,
1266 const tmp<faMatrix<Type>>& tB
1270 tmp<faMatrix<Type>> tC(tB.ptr());
1277 template<
class Type>
1280 const tmp<faMatrix<Type>>& tA,
1281 const tmp<faMatrix<Type>>& tB
1285 tmp<faMatrix<Type>> tC(tA.
ptr());
1292 template<
class Type>
1295 const faMatrix<Type>&
A,
1296 const faMatrix<Type>&
B
1304 template<
class Type>
1307 const tmp<faMatrix<Type>>& tA,
1308 const faMatrix<Type>&
B
1316 template<
class Type>
1319 const faMatrix<Type>&
A,
1320 const tmp<faMatrix<Type>>& tB
1328 template<
class Type>
1331 const tmp<faMatrix<Type>>& tA,
1332 const tmp<faMatrix<Type>>& tB
1340 template<
class Type>
1343 const faMatrix<Type>&
A,
1344 const GeometricField<Type, faPatchField, areaMesh>& su
1348 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1349 tC.ref().source() -= su.mesh().S()*su.internalField();
1354 template<
class Type>
1357 const tmp<faMatrix<Type>>& tA,
1358 const GeometricField<Type, faPatchField, areaMesh>& su
1362 tmp<faMatrix<Type>> tC(tA.
ptr());
1363 tC.ref().source() -= su.mesh().S()*su.internalField();
1368 template<
class Type>
1371 const faMatrix<Type>&
A,
1372 const tmp<GeometricField<Type, faPatchField, areaMesh>>& tsu
1376 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1377 tC.ref().source() -= tsu().mesh().S()*tsu().internalField();
1383 template<
class Type>
1386 const tmp<faMatrix<Type>>& tA,
1387 const tmp<GeometricField<Type, faPatchField, areaMesh>>& tsu
1391 tmp<faMatrix<Type>> tC(tA.
ptr());
1392 tC.ref().source() -= tsu().mesh().S()*tsu().internalField();
1398 template<
class Type>
1401 const GeometricField<Type, faPatchField, areaMesh>& su,
1402 const faMatrix<Type>&
A
1406 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1407 tC.ref().source() -= su.mesh().S()*su.internalField();
1412 template<
class Type>
1415 const GeometricField<Type, faPatchField, areaMesh>& su,
1416 const tmp<faMatrix<Type>>& tA
1420 tmp<faMatrix<Type>> tC(tA.
ptr());
1421 tC.ref().source() -= su.mesh().S()*su.internalField();
1426 template<
class Type>
1429 const tmp<GeometricField<Type, faPatchField, areaMesh>>& tsu,
1430 const faMatrix<Type>&
A
1434 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1435 tC.ref().source() -= tsu().mesh().S()*tsu().internalField();
1441 template<
class Type>
1444 const tmp<GeometricField<Type, faPatchField, areaMesh>>& tsu,
1445 const tmp<faMatrix<Type>>& tA
1449 tmp<faMatrix<Type>> tC(tA.
ptr());
1450 tC.ref().source() -= tsu().mesh().S()*tsu().internalField();
1456 template<
class Type>
1459 const faMatrix<Type>&
A,
1460 const GeometricField<Type, faPatchField, areaMesh>& su
1464 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1465 tC.ref().source() += su.mesh().S()*su.internalField();
1470 template<
class Type>
1473 const tmp<faMatrix<Type>>& tA,
1474 const GeometricField<Type, faPatchField, areaMesh>& su
1478 tmp<faMatrix<Type>> tC(tA.
ptr());
1479 tC.ref().source() += su.mesh().S()*su.internalField();
1484 template<
class Type>
1487 const faMatrix<Type>&
A,
1488 const tmp<GeometricField<Type, faPatchField, areaMesh>>& tsu
1492 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1493 tC.ref().source() += tsu().mesh().S()*tsu().internalField();
1499 template<
class Type>
1502 const tmp<faMatrix<Type>>& tA,
1503 const tmp<GeometricField<Type, faPatchField, areaMesh>>& tsu
1507 tmp<faMatrix<Type>> tC(tA.
ptr());
1508 tC.ref().source() += tsu().mesh().S()*tsu().internalField();
1514 template<
class Type>
1517 const GeometricField<Type, faPatchField, areaMesh>& su,
1518 const faMatrix<Type>&
A
1522 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1524 tC.ref().source() -= su.mesh().S()*su.internalField();
1529 template<
class Type>
1532 const GeometricField<Type, faPatchField, areaMesh>& su,
1533 const tmp<faMatrix<Type>>& tA
1537 tmp<faMatrix<Type>> tC(tA.
ptr());
1539 tC.ref().source() -= su.mesh().S()*su.internalField();
1544 template<
class Type>
1547 const tmp<GeometricField<Type, faPatchField, areaMesh>>& tsu,
1548 const faMatrix<Type>&
A
1552 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1554 tC.ref().source() -= tsu().mesh().S()*tsu().internalField();
1560 template<
class Type>
1563 const tmp<GeometricField<Type, faPatchField, areaMesh>>& tsu,
1564 const tmp<faMatrix<Type>>& tA
1568 tmp<faMatrix<Type>> tC(tA.
ptr());
1570 tC.ref().source() -= tsu().mesh().S()*tsu().internalField();
1576 template<
class Type>
1579 const faMatrix<Type>&
A,
1580 const dimensioned<Type>& su
1584 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1585 tC.ref().source() -= su.value()*
A.psi().mesh().S();
1590 template<
class Type>
1593 const tmp<faMatrix<Type>>& tA,
1594 const dimensioned<Type>& su
1598 tmp<faMatrix<Type>> tC(tA.
ptr());
1599 tC.ref().source() -= su.value()*tC().psi().mesh().S();
1604 template<
class Type>
1607 const dimensioned<Type>& su,
1608 const faMatrix<Type>&
A
1612 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1613 tC.ref().source() -= su.value()*
A.psi().mesh().S();
1618 template<
class Type>
1621 const dimensioned<Type>& su,
1622 const tmp<faMatrix<Type>>& tA
1626 tmp<faMatrix<Type>> tC(tA.
ptr());
1627 tC.ref().source() -= su.value()*tC().psi().mesh().S();
1632 template<
class Type>
1635 const faMatrix<Type>&
A,
1636 const dimensioned<Type>& su
1640 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1641 tC.ref().source() += su.value()*tC().psi().mesh().S();
1646 template<
class Type>
1649 const tmp<faMatrix<Type>>& tA,
1650 const dimensioned<Type>& su
1654 tmp<faMatrix<Type>> tC(tA.
ptr());
1655 tC.ref().source() += su.value()*tC().psi().mesh().S();
1660 template<
class Type>
1663 const dimensioned<Type>& su,
1664 const faMatrix<Type>&
A
1668 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1670 tC.ref().source() -= su.value()*
A.psi().mesh().S();
1675 template<
class Type>
1678 const dimensioned<Type>& su,
1679 const tmp<faMatrix<Type>>& tA
1683 tmp<faMatrix<Type>> tC(tA.
ptr());
1685 tC.ref().source() -= su.value()*tC().psi().mesh().S();
1690 template<
class Type>
1693 const faMatrix<Type>&
A,
1694 const GeometricField<Type, faPatchField, areaMesh>& su
1698 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1699 tC.ref().source() += su.mesh().S()*su.internalField();
1704 template<
class Type>
1707 const tmp<faMatrix<Type>>& tA,
1708 const GeometricField<Type, faPatchField, areaMesh>& su
1712 tmp<faMatrix<Type>> tC(tA.ptr());
1713 tC.ref().source() += su.mesh().S()*su.internalField();
1718 template<
class Type>
1721 const faMatrix<Type>&
A,
1722 const tmp<GeometricField<Type, faPatchField, areaMesh>>& tsu
1726 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1727 tC.ref().source() += tsu().mesh().S()*tsu().internalField();
1733 template<
class Type>
1736 const tmp<faMatrix<Type>>& tA,
1737 const tmp<GeometricField<Type, faPatchField, areaMesh>>& tsu
1741 tmp<faMatrix<Type>> tC(tA.ptr());
1742 tC.ref().source() += tsu().mesh().S()*tsu().internalField();
1748 template<
class Type>
1751 const faMatrix<Type>&
A,
1752 const dimensioned<Type>& su
1756 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1757 tC.ref().source() +=
A.psi().mesh().S()*su.value();
1762 template<
class Type>
1765 const tmp<faMatrix<Type>>& tA,
1766 const dimensioned<Type>& su
1770 tmp<faMatrix<Type>> tC(tA.ptr());
1771 tC.ref().source() += tC().psi().mesh().S()*su.value();
1776 template<
class Type>
1780 const faMatrix<Type>&
A
1783 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1789 template<
class Type>
1792 const tmp<areaScalarField>& tvsf,
1793 const faMatrix<Type>&
A
1796 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1802 template<
class Type>
1806 const tmp<faMatrix<Type>>& tA
1809 tmp<faMatrix<Type>> tC(tA.ptr());
1815 template<
class Type>
1818 const tmp<areaScalarField>& tvsf,
1819 const tmp<faMatrix<Type>>& tA
1822 tmp<faMatrix<Type>> tC(tA.ptr());
1828 template<
class Type>
1831 const dimensioned<scalar>& ds,
1832 const faMatrix<Type>&
A
1835 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1841 template<
class Type>
1844 const dimensioned<scalar>& ds,
1845 const tmp<faMatrix<Type>>& tA
1848 tmp<faMatrix<Type>> tC(tA.ptr());
1856 template<
class Type>
1859 os << static_cast<const lduMatrix&>(
fam) <<
nl
1860 <<
fam.dimensions_ <<
nl
1861 <<
fam.source_ <<
nl
1862 <<
fam.internalCoeffs_ <<
nl
1863 <<
fam.boundaryCoeffs_ <<
endl;