42 const Field<Type2>& pf,
46 if (addr.size() != pf.size())
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())
86 <<
"sizes of addressing and field are different"
92 intf[addr[faceI]] -= pf[faceI];
102 const tmp<Field<Type2>>& tpf,
106 subtractFromInternalField(addr, tpf(), intf);
118 forAll(internalCoeffs_, patchI)
122 lduAddr().patchAddr(patchI),
123 internalCoeffs_[patchI].
component(solveCmpt),
133 forAll(internalCoeffs_, patchI)
137 lduAddr().patchAddr(patchI),
138 cmptAv(internalCoeffs_[patchI]),
152 forAll(psi_.boundaryField(), patchI)
154 const faPatchField<Type>& ptf = psi_.boundaryField()[patchI];
155 const Field<Type>& pbc = boundaryCoeffs_[patchI];
159 addToInternalField(lduAddr().patchAddr(patchI), pbc, source);
163 tmp<Field<Type>> tpnf = ptf.patchNeighbourField();
164 const Field<Type>& pnf = tpnf();
166 const labelUList& addr = lduAddr().patchAddr(patchI);
170 source[addr[facei]] +=
cmptMultiply(pbc[facei], pnf[facei]);
190 internalCoeffs_(
psi.mesh().boundary().size()),
191 boundaryCoeffs_(
psi.mesh().boundary().size()),
192 faceFluxCorrectionPtr_(
nullptr)
195 <<
"constructing faMatrix<Type> for field " << psi_.name()
206 psi.mesh().boundary()[patchI].size(),
216 psi.mesh().boundary()[patchI].size(),
226 label currentStatePsi = psiRef.eventNo();
228 psiRef.eventNo() = currentStatePsi;
238 dimensions_(
fam.dimensions_),
239 source_(
fam.source_),
240 internalCoeffs_(
fam.internalCoeffs_),
241 boundaryCoeffs_(
fam.boundaryCoeffs_),
242 faceFluxCorrectionPtr_(nullptr)
245 <<
"copying faMatrix<Type> for field " << psi_.name()
248 if (
fam.faceFluxCorrectionPtr_)
250 faceFluxCorrectionPtr_ =
new
253 *(
fam.faceFluxCorrectionPtr_)
270 internalCoeffs_(
psi.mesh().boundary().size()),
271 boundaryCoeffs_(
psi.mesh().boundary().size()),
272 faceFluxCorrectionPtr_(
nullptr)
275 <<
"constructing faMatrix<Type> for field " << psi_.name()
286 psi.mesh().boundary()[patchI].size(),
296 psi.mesh().boundary()[patchI].size(),
322 <<
"destroying faMatrix<Type> for field " << psi_.name()
341 for (
const label i : faceLabels)
343 this->eliminatedEqns().
insert(i);
355 >(psi_).internalField();
360 label facei = faceLabels[i];
363 source_[facei] =
values[i]*Diag[facei];
365 if (symmetric() || asymmetric())
373 if (
mesh.isInternalEdge(edgei))
377 if (facei == own[edgei])
386 upper()[edgei] = 0.0;
390 if (facei == own[edgei])
399 upper()[edgei] = 0.0;
400 lower()[edgei] = 0.0;
407 if (internalCoeffs_[patchi].size())
412 internalCoeffs_[patchi][patchEdgei] =
415 boundaryCoeffs_[patchi][patchEdgei] =
432 if (psi_.needReference())
436 source()[facei] +=
diag()[facei]*value;
459 sumMagOffDiag(sumOff);
462 forAll(psi_.boundaryField(), patchI)
468 const labelUList& pa = lduAddr().patchAddr(patchI);
473 const Field<Type>& pCoeffs = boundaryCoeffs_[patchI];
491 Type iCoeff0 = iCoeffs[
face];
509 forAll(psi_.boundaryField(), patchI)
515 const labelUList& pa = lduAddr().patchAddr(patchI);
529 S += (
D - D0)*psi_.internalField();
536 if (psi_.mesh().relaxEquation(psi_.name()))
538 relax(psi_.mesh().equationRelaxationFactor(psi_.name()));
543 <<
"Relaxation factor for field " << psi_.name()
544 <<
" not found. Relaxation will not be used." <<
endl;
553 addCmptAvBoundaryDiag(tdiag.
ref());
567 "A("+psi_.name()+
')',
572 dimensions_/psi_.dimensions()/
dimArea,
573 zeroGradientFaPatchScalarField::typeName
594 "H("+psi_.name()+
')',
600 zeroGradientFaPatchScalarField::typeName
606 for (
direction cmpt=0; cmpt<Type::nComponents; ++cmpt)
608 scalarField psiCmpt(psi_.primitiveField().component(cmpt));
611 addBoundaryDiag(boundaryDiagCmpt, cmpt);
612 boundaryDiagCmpt.negate();
613 addCmptAvBoundaryDiag(boundaryDiagCmpt);
632 if (!psi_.mesh().fluxRequired(psi_.name()))
635 <<
"flux requested but " << psi_.name()
636 <<
" not specified in the fluxRequired sub-dictionary of faSchemes"
647 "flux("+psi_.name()+
')',
658 for (
direction cmpt=0; cmpt<pTraits<Type>::nComponents; ++cmpt)
669 forAll(InternalContrib, patchI)
671 InternalContrib[patchI] =
674 InternalContrib[patchI],
675 psi_.boundaryField()[patchI].patchInternalField()
681 forAll(NeighbourContrib, patchI)
683 if (psi_.boundaryField()[patchI].coupled())
685 NeighbourContrib[patchI] =
688 NeighbourContrib[patchI],
689 psi_.boundaryField()[patchI].patchNeighbourField()
697 InternalContrib[patchI] - NeighbourContrib[patchI];
700 if (faceFluxCorrectionPtr_)
702 fieldFlux += *faceFluxCorrectionPtr_;
719 if (&psi_ != &(famv.psi_))
722 <<
"different fields"
727 source_ = famv.source_;
728 internalCoeffs_ = famv.internalCoeffs_;
729 boundaryCoeffs_ = famv.boundaryCoeffs_;
731 if (faceFluxCorrectionPtr_ && famv.faceFluxCorrectionPtr_)
733 *faceFluxCorrectionPtr_ = *famv.faceFluxCorrectionPtr_;
735 else if (famv.faceFluxCorrectionPtr_)
737 faceFluxCorrectionPtr_ =
739 (*famv.faceFluxCorrectionPtr_);
757 internalCoeffs_.negate();
758 boundaryCoeffs_.negate();
760 if (faceFluxCorrectionPtr_)
762 faceFluxCorrectionPtr_->negate();
772 dimensions_ += famv.dimensions_;
774 source_ += famv.source_;
775 internalCoeffs_ += famv.internalCoeffs_;
776 boundaryCoeffs_ += famv.boundaryCoeffs_;
778 if (faceFluxCorrectionPtr_ && famv.faceFluxCorrectionPtr_)
780 *faceFluxCorrectionPtr_ += *famv.faceFluxCorrectionPtr_;
782 else if (famv.faceFluxCorrectionPtr_)
784 faceFluxCorrectionPtr_ =
new
787 *famv.faceFluxCorrectionPtr_
806 dimensions_ -= famv.dimensions_;
808 source_ -= famv.source_;
809 internalCoeffs_ -= famv.internalCoeffs_;
810 boundaryCoeffs_ -= famv.boundaryCoeffs_;
812 if (faceFluxCorrectionPtr_ && famv.faceFluxCorrectionPtr_)
814 *faceFluxCorrectionPtr_ -= *famv.faceFluxCorrectionPtr_;
816 else if (famv.faceFluxCorrectionPtr_)
818 faceFluxCorrectionPtr_ =
820 (-*famv.faceFluxCorrectionPtr_);
840 source() -= su.mesh().S()*su.internalField();
862 source() += su.mesh().S()*su.internalField();
883 source() -= su.mesh().S()*su;
893 source() += su.mesh().S()*su;
903 dimensions_ *= vsf.dimensions();
905 source_ *= vsf.internalField();
907 forAll(boundaryCoeffs_, patchI)
909 const scalarField psf(vsf.boundaryField()[patchI].patchInternalField());
910 internalCoeffs_[patchI] *= psf;
911 boundaryCoeffs_[patchI] *= psf;
914 if (faceFluxCorrectionPtr_)
917 <<
"cannot scale a matrix containing a faceFluxCorrection"
940 dimensions_ *= ds.dimensions();
942 source_ *= ds.value();
943 internalCoeffs_ *= ds.value();
944 boundaryCoeffs_ *= ds.value();
946 if (faceFluxCorrectionPtr_)
948 *faceFluxCorrectionPtr_ *= ds.value();
963 if (&fam1.
psi() != &fam2.
psi())
966 <<
"incompatible fields for operation "
968 <<
"[" << fam1.
psi().name() <<
"] "
970 <<
" [" << fam2.
psi().name() <<
"]"
977 <<
"incompatible dimensions for operation "
998 <<
"incompatible dimensions for operation "
1000 <<
"[" <<
fam.psi().name() <<
fam.dimensions()/
dimArea <<
" ] "
1002 <<
" [" << vf.name() << vf.dimensions() <<
" ]"
1008 template<
class Type>
1019 <<
"incompatible dimensions for operation "
1021 <<
"[" <<
fam.psi().name() <<
fam.dimensions()/
dimArea <<
" ] "
1029 template<
class Type>
1032 faMatrix<Type>&
fam,
1033 Istream& solverControls
1036 return fam.solve(solverControls);
1040 template<
class Type>
1043 const tmp<faMatrix<Type>>& tfam,
1044 Istream& solverControls
1047 SolverPerformance<Type> solverPerf =
1048 const_cast<faMatrix<Type>&
>(tfam()).
solve(solverControls);
1055 template<
class Type>
1062 template<
class Type>
1065 SolverPerformance<Type> solverPerf =
1066 const_cast<faMatrix<Type>&
>(tfam()).
solve();
1075 template<
class Type>
1078 const faMatrix<Type>&
A,
1079 const faMatrix<Type>&
B
1083 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1089 template<
class Type>
1092 const tmp<faMatrix<Type>>& tA,
1093 const faMatrix<Type>&
B
1097 tmp<faMatrix<Type>> tC(tA.
ptr());
1103 template<
class Type>
1106 const faMatrix<Type>&
A,
1107 const tmp<faMatrix<Type>>& tB
1111 tmp<faMatrix<Type>> tC(tB.ptr());
1117 template<
class Type>
1120 const tmp<faMatrix<Type>>& tA,
1121 const tmp<faMatrix<Type>>& tB
1125 tmp<faMatrix<Type>> tC(tA.
ptr());
1132 template<
class Type>
1135 const faMatrix<Type>&
A
1138 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1144 template<
class Type>
1147 const tmp<faMatrix<Type>>& tA
1150 tmp<faMatrix<Type>> tC(tA.
ptr());
1156 template<
class Type>
1159 const faMatrix<Type>&
A,
1160 const faMatrix<Type>&
B
1164 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1170 template<
class Type>
1173 const tmp<faMatrix<Type>>& tA,
1174 const faMatrix<Type>&
B
1178 tmp<faMatrix<Type>> tC(tA.
ptr());
1184 template<
class Type>
1187 const faMatrix<Type>&
A,
1188 const tmp<faMatrix<Type>>& tB
1192 tmp<faMatrix<Type>> tC(tB.ptr());
1199 template<
class Type>
1202 const tmp<faMatrix<Type>>& tA,
1203 const tmp<faMatrix<Type>>& tB
1207 tmp<faMatrix<Type>> tC(tA.
ptr());
1214 template<
class Type>
1217 const faMatrix<Type>&
A,
1218 const faMatrix<Type>&
B
1226 template<
class Type>
1229 const tmp<faMatrix<Type>>& tA,
1230 const faMatrix<Type>&
B
1238 template<
class Type>
1241 const faMatrix<Type>&
A,
1242 const tmp<faMatrix<Type>>& tB
1250 template<
class Type>
1253 const tmp<faMatrix<Type>>& tA,
1254 const tmp<faMatrix<Type>>& tB
1262 template<
class Type>
1265 const faMatrix<Type>&
A,
1266 const GeometricField<Type, faPatchField, areaMesh>& su
1270 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1271 tC.ref().source() -= su.mesh().S()*su.internalField();
1276 template<
class Type>
1279 const tmp<faMatrix<Type>>& tA,
1280 const GeometricField<Type, faPatchField, areaMesh>& su
1284 tmp<faMatrix<Type>> tC(tA.
ptr());
1285 tC.ref().source() -= su.mesh().S()*su.internalField();
1290 template<
class Type>
1293 const faMatrix<Type>&
A,
1294 const tmp<GeometricField<Type, faPatchField, areaMesh>>& tsu
1298 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1299 tC.ref().source() -= tsu().mesh().S()*tsu().internalField();
1305 template<
class Type>
1308 const tmp<faMatrix<Type>>& tA,
1309 const tmp<GeometricField<Type, faPatchField, areaMesh>>& tsu
1313 tmp<faMatrix<Type>> tC(tA.
ptr());
1314 tC.ref().source() -= tsu().mesh().S()*tsu().internalField();
1320 template<
class Type>
1323 const GeometricField<Type, faPatchField, areaMesh>& su,
1324 const faMatrix<Type>&
A
1328 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1329 tC.ref().source() -= su.mesh().S()*su.internalField();
1334 template<
class Type>
1337 const GeometricField<Type, faPatchField, areaMesh>& su,
1338 const tmp<faMatrix<Type>>& tA
1342 tmp<faMatrix<Type>> tC(tA.
ptr());
1343 tC.ref().source() -= su.mesh().S()*su.internalField();
1348 template<
class Type>
1351 const tmp<GeometricField<Type, faPatchField, areaMesh>>& tsu,
1352 const faMatrix<Type>&
A
1356 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1357 tC.ref().source() -= tsu().mesh().S()*tsu().internalField();
1363 template<
class Type>
1366 const tmp<GeometricField<Type, faPatchField, areaMesh>>& tsu,
1367 const tmp<faMatrix<Type>>& tA
1371 tmp<faMatrix<Type>> tC(tA.
ptr());
1372 tC.ref().source() -= tsu().mesh().S()*tsu().internalField();
1378 template<
class Type>
1381 const faMatrix<Type>&
A,
1382 const GeometricField<Type, faPatchField, areaMesh>& su
1386 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1387 tC.ref().source() += su.mesh().S()*su.internalField();
1392 template<
class Type>
1395 const tmp<faMatrix<Type>>& tA,
1396 const GeometricField<Type, faPatchField, areaMesh>& su
1400 tmp<faMatrix<Type>> tC(tA.
ptr());
1401 tC.ref().source() += su.mesh().S()*su.internalField();
1406 template<
class Type>
1409 const faMatrix<Type>&
A,
1410 const tmp<GeometricField<Type, faPatchField, areaMesh>>& tsu
1414 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1415 tC.ref().source() += tsu().mesh().S()*tsu().internalField();
1421 template<
class Type>
1424 const tmp<faMatrix<Type>>& tA,
1425 const tmp<GeometricField<Type, faPatchField, areaMesh>>& tsu
1429 tmp<faMatrix<Type>> tC(tA.
ptr());
1430 tC.ref().source() += tsu().mesh().S()*tsu().internalField();
1436 template<
class Type>
1439 const GeometricField<Type, faPatchField, areaMesh>& su,
1440 const faMatrix<Type>&
A
1444 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1446 tC.ref().source() -= su.mesh().S()*su.internalField();
1451 template<
class Type>
1454 const GeometricField<Type, faPatchField, areaMesh>& su,
1455 const tmp<faMatrix<Type>>& tA
1459 tmp<faMatrix<Type>> tC(tA.
ptr());
1461 tC.ref().source() -= su.mesh().S()*su.internalField();
1466 template<
class Type>
1469 const tmp<GeometricField<Type, faPatchField, areaMesh>>& tsu,
1470 const faMatrix<Type>&
A
1474 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1476 tC.ref().source() -= tsu().mesh().S()*tsu().internalField();
1482 template<
class Type>
1485 const tmp<GeometricField<Type, faPatchField, areaMesh>>& tsu,
1486 const tmp<faMatrix<Type>>& tA
1490 tmp<faMatrix<Type>> tC(tA.
ptr());
1492 tC.ref().source() -= tsu().mesh().S()*tsu().internalField();
1498 template<
class Type>
1501 const faMatrix<Type>&
A,
1502 const dimensioned<Type>& su
1506 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1507 tC.ref().source() -= su.value()*
A.psi().mesh().S();
1512 template<
class Type>
1515 const tmp<faMatrix<Type>>& tA,
1516 const dimensioned<Type>& su
1520 tmp<faMatrix<Type>> tC(tA.
ptr());
1521 tC.ref().source() -= su.value()*tC().psi().mesh().S();
1526 template<
class Type>
1529 const dimensioned<Type>& su,
1530 const faMatrix<Type>&
A
1534 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1535 tC.ref().source() -= su.value()*
A.psi().mesh().S();
1540 template<
class Type>
1543 const dimensioned<Type>& su,
1544 const tmp<faMatrix<Type>>& tA
1548 tmp<faMatrix<Type>> tC(tA.
ptr());
1549 tC.ref().source() -= su.value()*tC().psi().mesh().S();
1554 template<
class Type>
1557 const faMatrix<Type>&
A,
1558 const dimensioned<Type>& su
1562 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1563 tC.ref().source() += su.value()*tC().psi().mesh().S();
1568 template<
class Type>
1571 const tmp<faMatrix<Type>>& tA,
1572 const dimensioned<Type>& su
1576 tmp<faMatrix<Type>> tC(tA.
ptr());
1577 tC.ref().source() += su.value()*tC().psi().mesh().S();
1582 template<
class Type>
1585 const dimensioned<Type>& su,
1586 const faMatrix<Type>&
A
1590 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1592 tC.ref().source() -= su.value()*
A.psi().mesh().S();
1597 template<
class Type>
1600 const dimensioned<Type>& su,
1601 const tmp<faMatrix<Type>>& tA
1605 tmp<faMatrix<Type>> tC(tA.
ptr());
1607 tC.ref().source() -= su.value()*tC().psi().mesh().S();
1612 template<
class Type>
1615 const faMatrix<Type>&
A,
1616 const GeometricField<Type, faPatchField, areaMesh>& su
1620 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1621 tC.ref().source() += su.mesh().S()*su.internalField();
1626 template<
class Type>
1629 const tmp<faMatrix<Type>>& tA,
1630 const GeometricField<Type, faPatchField, areaMesh>& su
1634 tmp<faMatrix<Type>> tC(tA.ptr());
1635 tC.ref().source() += su.mesh().S()*su.internalField();
1640 template<
class Type>
1643 const faMatrix<Type>&
A,
1644 const tmp<GeometricField<Type, faPatchField, areaMesh>>& tsu
1648 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1649 tC.ref().source() += tsu().mesh().S()*tsu().internalField();
1655 template<
class Type>
1658 const tmp<faMatrix<Type>>& tA,
1659 const tmp<GeometricField<Type, faPatchField, areaMesh>>& tsu
1663 tmp<faMatrix<Type>> tC(tA.ptr());
1664 tC.ref().source() += tsu().mesh().S()*tsu().internalField();
1670 template<
class Type>
1673 const faMatrix<Type>&
A,
1674 const dimensioned<Type>& su
1678 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1679 tC.ref().source() +=
A.psi().mesh().S()*su.value();
1684 template<
class Type>
1687 const tmp<faMatrix<Type>>& tA,
1688 const dimensioned<Type>& su
1692 tmp<faMatrix<Type>> tC(tA.ptr());
1693 tC.ref().source() += tC().psi().mesh().S()*su.value();
1698 template<
class Type>
1702 const faMatrix<Type>&
A
1705 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1711 template<
class Type>
1714 const tmp<areaScalarField>& tvsf,
1715 const faMatrix<Type>&
A
1718 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1724 template<
class Type>
1728 const tmp<faMatrix<Type>>& tA
1731 tmp<faMatrix<Type>> tC(tA.ptr());
1737 template<
class Type>
1740 const tmp<areaScalarField>& tvsf,
1741 const tmp<faMatrix<Type>>& tA
1744 tmp<faMatrix<Type>> tC(tA.ptr());
1750 template<
class Type>
1753 const dimensioned<scalar>& ds,
1754 const faMatrix<Type>&
A
1757 tmp<faMatrix<Type>> tC(
new faMatrix<Type>(
A));
1763 template<
class Type>
1766 const dimensioned<scalar>& ds,
1767 const tmp<faMatrix<Type>>& tA
1770 tmp<faMatrix<Type>> tC(tA.ptr());
1778 template<
class Type>
1781 os << static_cast<const lduMatrix&>(
fam) <<
nl
1782 <<
fam.dimensions_ <<
nl
1783 <<
fam.source_ <<
nl
1784 <<
fam.internalCoeffs_ <<
nl
1785 <<
fam.boundaryCoeffs_ <<
endl;