1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | foam-extend: Open Source CFD
4 \\ / O peration | Version: 3.2
5 \\ / A nd | Web: http://www.foam-extend.org
6 \\/ M anipulation | For copyright notice see file Copyright
7 -------------------------------------------------------------------------------
9 This file is part of foam-extend.
11 foam-extend is free software: you can redistribute it and/or modify it
12 under the terms of the GNU General Public License as published by the
13 Free Software Foundation, either version 3 of the License, or (at your
14 option) any later version.
16 foam-extend is distributed in the hope that it will be useful, but
17 WITHOUT ANY WARRANTY; without even the implied warranty of
18 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 General Public License for more details.
21 You should have received a copy of the GNU General Public License
22 along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
24 \*---------------------------------------------------------------------------*/
26 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 namespace blockMatrixTools
36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
38 template<class BlockType>
43 Field<BlockType>& blockX
48 blockX[i](dir) = x[i];
53 template<class BlockType>
58 Field<BlockType>& blockX
63 blockX[i](dir) += x[i];
68 template<class BlockType>
73 const Field<BlockType>& blockX
78 x[i] = blockX[i](dir);
83 template<class BlockType>
87 const fvScalarMatrix& m,
88 BlockLduMatrix<BlockType>& blockM,
89 Field<BlockType>& blockB
92 // Prepare the diagonal and source
94 scalarField diag = m.diag();
95 scalarField source = m.source();
97 // Add boundary source contribution
98 m.addBoundaryDiag(diag, 0);
99 m.addBoundarySource(source, false);
101 if (blockM.diag().activeType() == blockCoeffBase::UNALLOCATED)
103 blockM.diag().asScalar() = diag;
107 blockM.diag().activeType() == blockCoeffBase::SCALAR
108 || blockM.diag().activeType() == blockCoeffBase::LINEAR
111 typename CoeffField<BlockType>::linearTypeField& blockDiag =
112 blockM.diag().asLinear();
116 blockDiag[i](dir) = diag[i];
119 else if (blockM.diag().activeType() == blockCoeffBase::SQUARE)
121 typename CoeffField<BlockType>::squareTypeField& blockDiag =
122 blockM.diag().asSquare();
126 blockDiag[i](dir, dir) = diag[i];
130 blockInsert(dir, source, blockB);
134 template<class BlockType>
135 void insertUpperLower
138 const fvScalarMatrix& m,
139 BlockLduMatrix<BlockType>& blockM
144 // Matrix for insertion is diagonal-only: nothing to do
148 if (m.symmetric() && blockM.symmetric())
150 Info<< "Both m and blockM are symmetric: inserting only upper triangle"
155 // Either scalar or block matrix is asymmetric: insert lower triangle
156 const scalarField& lower = m.lower();
158 if (blockM.lower().activeType() == blockCoeffBase::UNALLOCATED)
160 blockM.lower().asScalar() = lower;
164 blockM.lower().activeType() == blockCoeffBase::SCALAR
165 || blockM.lower().activeType() == blockCoeffBase::LINEAR
168 typename CoeffField<BlockType>::linearTypeField& blockLower =
169 blockM.lower().asLinear();
173 blockLower[i](dir) = lower[i];
176 else if (blockM.lower().activeType() == blockCoeffBase::SQUARE)
178 typename CoeffField<BlockType>::squareTypeField& blockLower =
179 blockM.lower().asSquare();
183 blockLower[i](dir, dir) = lower[i];
190 const scalarField& upper = m.upper();
192 if (blockM.upper().activeType() == blockCoeffBase::UNALLOCATED)
194 blockM.upper().asScalar() = upper;
198 blockM.upper().activeType() == blockCoeffBase::SCALAR
199 || blockM.upper().activeType() == blockCoeffBase::LINEAR
202 typename CoeffField<BlockType>::linearTypeField& blockUpper =
203 blockM.upper().asLinear();
207 blockUpper[i](dir) = upper[i];
210 else if (blockM.upper().activeType() == blockCoeffBase::SQUARE)
212 typename CoeffField<BlockType>::squareTypeField& blockUpper =
213 blockM.upper().asSquare();
217 blockUpper[i](dir, dir) = upper[i];
225 "void insertUpperLower\n"
227 " const direction dir,\n"
228 " const fvScalarMatrix& m,\n"
229 " BlockLduMatrix<BlockType>& blockM\n"
231 ) << "Error in matrix insertion: problem with block structure"
232 << abort(FatalError);
235 // Insert block interface fields
236 forAll(blockM.interfaces(), patchI)
238 if (blockM.interfaces().set(patchI))
240 // Couple upper and lower
241 const scalarField& cUpper = m.boundaryCoeffs()[patchI];
242 const scalarField& cLower = m.internalCoeffs()[patchI];
246 blockM.coupleUpper()[patchI].activeType()
247 == blockCoeffBase::UNALLOCATED
250 blockM.coupleUpper()[patchI].asScalar() = cUpper;
251 blockM.coupleLower()[patchI].asScalar() = cLower;
255 blockM.coupleUpper()[patchI].activeType()
256 == blockCoeffBase::SCALAR
257 || blockM.coupleUpper()[patchI].activeType()
258 == blockCoeffBase::LINEAR
261 typename CoeffField<BlockType>::linearTypeField& blockUpper =
262 blockM.coupleUpper()[patchI].asLinear();
264 typename CoeffField<BlockType>::linearTypeField& blockLower =
265 blockM.coupleLower()[patchI].asLinear();
269 blockUpper[i](dir) = cUpper[i];
270 blockLower[i](dir) = cLower[i];
275 blockM.coupleUpper()[patchI].activeType()
276 == blockCoeffBase::SQUARE
279 typename CoeffField<BlockType>::squareTypeField& blockUpper =
280 blockM.coupleUpper()[patchI].asSquare();
282 typename CoeffField<BlockType>::squareTypeField& blockLower =
283 blockM.coupleLower()[patchI].asSquare();
287 blockUpper[i](dir, dir) = cUpper[i];
288 blockLower[i](dir, dir) = cLower[i];
296 template<class BlockType>
300 const fvScalarMatrix& m,
301 BlockLduMatrix<BlockType>& blockM,
302 Field<BlockType>& blockX,
303 Field<BlockType>& blockB
306 insertDiagSource(dir, m, blockM, blockB);
307 insertUpperLower(dir, m, blockM);
308 blockInsert(dir, m.psi(), blockX);
312 template<class BlockType>
313 void insertCouplingDiagSource
315 const direction dirI,
316 const direction dirJ,
317 const fvScalarMatrix& m,
318 BlockLduMatrix<BlockType>& blockM,
319 Field<BlockType>& blockB
322 // Prepare the diagonal and source
324 scalarField diag = m.diag();
325 scalarField source = m.source();
327 // Add boundary source contribution
328 m.addBoundaryDiag(diag, 0);
329 m.addBoundarySource(source, false);
332 // Add off-diagonal block coefficients
333 typename CoeffField<BlockType>::squareTypeField& blockDiag =
334 blockM.diag().asSquare();
336 // Set off-diagonal coefficient
339 blockDiag[i](dirI, dirJ) = diag[i];
342 blockAdd(dirI, source, blockB);
346 template<class BlockType>
347 void insertCouplingUpperLower
349 const direction dirI,
350 const direction dirJ,
351 const fvScalarMatrix& m,
352 BlockLduMatrix<BlockType>& blockM
357 // Matrix for insertion is diagonal-only: nothing to do
361 if (m.symmetric() && blockM.symmetric())
363 Info<< "Both m and blockM are symmetric: inserting only upper triangle"
368 // Either scalar or block matrix is asymmetric: insert lower triangle
369 const scalarField& lower = m.lower();
371 typename CoeffField<BlockType>::squareTypeField& blockLower =
372 blockM.lower().asSquare();
376 blockLower[i](dirI, dirJ) = lower[i];
382 const scalarField& upper = m.upper();
384 typename CoeffField<BlockType>::squareTypeField& blockUpper =
385 blockM.upper().asSquare();
389 blockUpper[i](dirI, dirJ) = upper[i];
396 "void insertCouplingUpperLower\n"
398 " const direction dirI,\n"
399 " const direction dirJ,\n"
400 " const fvScalarMatrix& m,\n"
401 " BlockLduMatrix<BlockType>& blockM\n"
403 ) << "Error in matrix insertion: problem with block structure"
404 << abort(FatalError);
407 // Insert block interface fields
408 forAll(blockM.interfaces(), patchI)
410 if (blockM.interfaces().set(patchI))
412 // Couple upper and lower
413 const scalarField& cUpper = m.boundaryCoeffs()[patchI];
414 const scalarField& cLower = m.internalCoeffs()[patchI];
416 typename CoeffField<BlockType>::squareTypeField& blockUpper =
417 blockM.coupleUpper()[patchI].asSquare();
419 typename CoeffField<BlockType>::squareTypeField& blockLower =
420 blockM.coupleLower()[patchI].asSquare();
424 blockUpper[i](dirI, dirJ) = cUpper[i];
425 blockLower[i](dirI, dirJ) = cLower[i];
432 template<class BlockType>
435 const direction dirI,
436 const direction dirJ,
437 const fvScalarMatrix& m,
438 BlockLduMatrix<BlockType>& blockM,
439 Field<BlockType>& blockX,
440 Field<BlockType>& blockB
443 insertCouplingDiagSource(dirI, dirJ, m, blockM, blockB);
444 insertCouplingUpperLower(dirI, dirJ, m, blockM);
448 template<class BlockType>
449 void updateSourceCoupling
451 BlockLduMatrix<BlockType>& blockM,
456 // Eliminate off-diagonal block coefficients from the square diagonal
457 // With this change, the segregated matrix can be assembled with complete
458 // source terms and linearisation can be provided independently.
459 // Once the linearisation coefficients are set (off-diagonal entries
460 // in the square block matrix, they are multiplied by the current value
461 // of the field and subtracted from the source term
463 if (blockM.diag().activeType() == blockCoeffBase::SQUARE)
465 typename CoeffField<BlockType>::squareTypeField& blockDiag =
466 blockM.diag().asSquare();
468 typename CoeffField<BlockType>::linearTypeField lf(blockDiag.size());
469 typename CoeffField<BlockType>::squareTypeField sf(blockDiag.size());
471 // Expand and contract
473 // Take out the diagonal entries from the square coefficient
474 contractLinear(lf, blockDiag);
476 // Expand the diagonal for full square, with zeroes in the off-diagonal
477 expandLinear(sf, lf);
479 // Subtract from the source the difference between the full block
480 // diagonal and the diagonal terms only
481 // Sign is the same as in the derivative
482 b += (blockDiag - sf) & x;
486 template<class blockType, class fieldType>
487 void insertSolutionVector
490 const Field<fieldType>& xSingle,
491 Field<blockType>& xBlock
494 // Get number of field components and local copy of location, for
495 // consistency with member functions where locations need to be reset.
496 const label nCmpts = pTraits<fieldType>::nComponents;
497 label localLoc = loc;
499 for (label cmptI = 0; cmptI < nCmpts; cmptI++)
501 scalarField xSingleCurr(xSingle.component(cmptI));
503 forAll (xSingleCurr, cellI)
505 xBlock[cellI](localLoc) = xSingleCurr[cellI];
513 template<class blockType, class fieldType>
514 void retrieveSolution
517 Field<fieldType>& xSingle,
518 const Field<blockType>& xBlock
521 const label nCmpts = pTraits<fieldType>::nComponents;
522 label localLoc = loc;
524 for (label cmptI = 0; cmptI < nCmpts; cmptI++)
526 scalarField xSingleCurr(xSingle.component(cmptI));
528 forAll (xSingleCurr, cellI)
530 xSingleCurr[cellI] = xBlock[cellI](localLoc);
533 xSingle.replace(cmptI, xSingleCurr);
540 template<class blockType, class matrixType>
541 void insertDiagSource
544 fvMatrix<matrixType>& matrix,
545 BlockLduMatrix<blockType>& A,
549 matrix.completeAssembly();
551 // Save a copy for different components
552 scalarField& diag = matrix.diag();
553 scalarField saveDiag(diag);
555 // Add source boundary contribution
556 Field<matrixType>& source = matrix.source();
557 matrix.addBoundarySource(source, false);
559 const label nCmpts = pTraits<matrixType>::nComponents;
560 label localLoc = loc;
564 // This is needed if the matrixType is <vector>, then you need to grab
565 // coeffs as linear. Consider doing a matrixType check also.
567 A.diag().activeType() != blockCoeffBase::SQUARE
570 typename CoeffField<blockType>::linearTypeField& blockDiag =
573 for (label cmptI = 0; cmptI < nCmpts; cmptI++)
575 matrix.addBoundaryDiag(diag, cmptI);
576 scalarField sourceCmpt(source.component(cmptI));
578 // FieldField<Field, scalar> bouCoeffsCmpt
580 // matrix.boundaryCoeffs().component(cmptI)
583 // Possible problem for coupled non-aligned boundaries.
585 // matrix.correctImplicitBoundarySource
594 blockDiag[cellI](localLoc) = diag[cellI];
595 b[cellI](localLoc) += sourceCmpt[cellI];
604 else if (A.diag().activeType() == blockCoeffBase::SQUARE)
606 typename CoeffField<blockType>::squareTypeField& blockDiag =
609 for (label cmptI = 0; cmptI < nCmpts; cmptI++)
611 matrix.addBoundaryDiag(diag, cmptI);
612 scalarField sourceCmpt(source.component(cmptI));
614 // FieldField<Field, scalar> bouCoeffsCmpt
616 // matrix.boundaryCoeffs().component(cmptI)
619 // matrix.correctImplicitBoundarySource
628 blockDiag[cellI](localLoc, localLoc) = diag[cellI];
629 b[cellI](localLoc) += sourceCmpt[cellI];
641 template<class blockType, class matrixType>
642 void insertUpperLower
645 const fvMatrix<matrixType>& matrix,
646 BlockLduMatrix<blockType>& A
649 if (matrix.diagonal())
651 // Matrix for insertion is diagonal-only: nothing to do
655 const label nCmpts = pTraits<matrixType>::nComponents;
656 label localLoc = loc;
658 if (matrix.hasUpper())
660 const scalarField& upper = matrix.upper();
662 if (A.upper().activeType() == blockCoeffBase::UNALLOCATED)
664 A.upper().asScalar() = upper;
668 A.upper().activeType() == blockCoeffBase::SCALAR
669 || A.upper().activeType() == blockCoeffBase::LINEAR
672 typename CoeffField<blockType>::linearTypeField& blockUpper =
673 A.upper().asLinear();
675 for (label cmptI = 0; cmptI < nCmpts; cmptI++)
677 forAll (upper, faceI)
679 blockUpper[faceI](localLoc) = upper[faceI];
685 else if (A.upper().activeType() == blockCoeffBase::SQUARE)
687 typename CoeffField<blockType>::squareTypeField& blockUpper =
688 A.upper().asSquare();
690 for (label cmptI = 0; cmptI < nCmpts; cmptI++)
692 forAll (upper, faceI)
694 blockUpper[faceI](localLoc, localLoc) = upper[faceI];
705 "void insertUpperLower\n"
707 " const location loc,\n"
708 " const fvMatrix<matrixType>& matrix,\n"
709 " BlockLduMatrix<blockType>& A\n"
711 ) << "Error in matrix insertion: problem with block structure."
712 << abort(FatalError);
715 if (matrix.symmetric() && A.symmetric())
717 Info<< "Both matrices are symmetric: inserting only upper triangle"
725 // Either scalar or block matrix is asymmetric: insert lower triangle
726 const scalarField& lower = matrix.lower();
728 if (A.lower().activeType() == blockCoeffBase::UNALLOCATED)
730 A.lower().asScalar() = lower;
734 A.lower().activeType() == blockCoeffBase::SCALAR
735 || A.lower().activeType() == blockCoeffBase::LINEAR
738 typename CoeffField<blockType>::linearTypeField& blockLower =
739 A.lower().asLinear();
741 for (label cmptI = 0; cmptI < nCmpts; cmptI++)
743 forAll (lower, faceI)
745 blockLower[faceI](localLoc) = lower[faceI];
751 else if (A.lower().activeType() == blockCoeffBase::SQUARE)
753 typename CoeffField<blockType>::squareTypeField& blockLower =
754 A.lower().asSquare();
756 for (label cmptI = 0; cmptI < nCmpts; cmptI++)
758 forAll (lower, faceI)
760 blockLower[faceI](localLoc, localLoc) = lower[faceI];
770 template<class blockType, class matrixType>
771 void updateCouplingCoeffs
774 const fvMatrix<matrixType>& matrix,
775 BlockLduMatrix<blockType>& A
778 const label nCmpts = pTraits<matrixType>::nComponents;
779 label localLoc = loc;
781 const GeometricField<matrixType, fvPatchField, volMesh>& psi = matrix.psi();
782 forAll(psi.boundaryField(), patchI)
784 const fvPatchField<matrixType>& pf = psi.boundaryField()[patchI];
785 const fvPatch& patch = pf.patch();
789 const Field<matrixType>& icp = matrix.internalCoeffs()[patchI];
790 const Field<matrixType>& bcp = matrix.boundaryCoeffs()[patchI];
792 if (A.coupleUpper()[patchI].activeType() != blockCoeffBase::SQUARE)
794 typename CoeffField<blockType>::linearTypeField& pcoupleUpper =
795 A.coupleUpper()[patchI].asLinear();
796 typename CoeffField<blockType>::linearTypeField& pcoupleLower =
797 A.coupleLower()[patchI].asLinear();
799 for (label cmptI = 0; cmptI < nCmpts; cmptI++)
801 scalarField icpCmpt = icp.component(cmptI);
802 scalarField bcpCmpt = bcp.component(cmptI);
806 pcoupleUpper[faceI](localLoc) = bcpCmpt[faceI];
807 pcoupleLower[faceI](localLoc) = icpCmpt[faceI];
817 A.coupleUpper()[patchI].activeType() == blockCoeffBase::SQUARE
820 typename CoeffField<blockType>::squareTypeField& pcoupleUpper =
821 A.coupleUpper()[patchI].asSquare();
822 typename CoeffField<blockType>::squareTypeField& pcoupleLower =
823 A.coupleLower()[patchI].asSquare();
825 for (label cmptI = 0; cmptI < nCmpts; cmptI++)
827 scalarField icpCmpt = icp.component(cmptI);
828 scalarField bcpCmpt = bcp.component(cmptI);
832 pcoupleUpper[faceI](localLoc, localLoc) =
834 pcoupleLower[faceI](localLoc, localLoc) =
848 template<class blockType, class matrixType>
852 fvMatrix<matrixType>& matrix,
853 BlockLduMatrix<blockType>& A,
858 insertDiagSource(loc, matrix, A, b);
859 insertUpperLower(loc, matrix, A);
860 updateCouplingCoeffs(loc, matrix, A);
861 insertSolutionVector(loc, matrix.psi().internalField(), x);
865 template<class blockType1, class fieldType, class blockType2>
870 const BlockLduSystem<blockType1, fieldType>& blockSystem,
871 BlockLduMatrix<blockType2>& A,
877 const label blockMatrixSize = pTraits<blockType1>::nComponents;
878 const label blockMatrixASize = pTraits<blockType2>::nComponents;
880 if (blockMatrixSize > blockMatrixASize)
886 " const label loc1,\n"
887 " const label loc2,\n"
888 " BlockLduMatrix<blockType1>& blockMatrix,\n"
889 " BlockLduMatrix<blockType2>& A\n"
891 ) << "Trying to insert a block matrix into smaller one."
892 << abort(FatalError);
899 "void insertCoupling\n"
901 " const label loc1,\n"
902 " const label loc2,\n"
903 " BlockLduMatrix<blockType1>& blockMatrix,\n"
904 " BlockLduMatrix<blockType2>& A\n"
906 ) << "Trying to insert coupling in the position where equation "
907 << "should be, since loc1 = loc2. Try using insertEquatiion "
908 << "member function."
909 << abort(FatalError);
913 const label nCmpts = pTraits<blockType1>::nComponents;
914 label localLoc1 = loc1;
915 label localLoc2 = loc2;
917 // Get references to ldu fields of blockMatrix always as linear
918 const typename CoeffField<blockType1>::linearTypeField& bmd =
919 blockSystem.diag().asLinear();
920 const typename CoeffField<blockType1>::linearTypeField& bmu =
921 blockSystem.upper().asLinear();
922 const typename CoeffField<blockType1>::linearTypeField& bml =
923 blockSystem.lower().asLinear();
925 // Get references to ldu fields of A matrix always as square
926 typename CoeffField<blockType2>::squareTypeField& blockDiag =
928 typename CoeffField<blockType2>::squareTypeField& blockUpper =
929 A.upper().asSquare();
930 typename CoeffField<blockType2>::squareTypeField& blockLower =
931 A.lower().asSquare();
933 // Insert blockMatrix that represents coupling into larger system matrix
934 for (label cmptI = 0; cmptI < nCmpts; cmptI++)
938 blockDiag[cellI](localLoc1, localLoc2) +=
939 bmd[cellI].component(cmptI);
944 blockUpper[faceI](localLoc1, localLoc2) +=
945 bmu[faceI].component(cmptI);
946 blockLower[faceI](localLoc1, localLoc2) +=
947 bml[faceI].component(cmptI);
962 template<class blockType1, class blockType2, class fieldType>
963 void insertBoundaryContributions
967 const BlockLduSystem<blockType1, fieldType>& blockSystem,
968 BlockLduMatrix<blockType2>& A,
969 Field<blockType2>& b,
973 // Need to get reference to fvMesh instead of lduMesh
974 const fvBoundaryMesh& bmesh = refCast<const fvMesh>(A.mesh()).boundary();
976 const label nCmpts = pTraits<blockType1>::nComponents;
977 const label nSrcCmpts = pTraits<fieldType>::nComponents;
978 label localLoc1 = loc1;
979 label localLoc2 = loc2;
981 const Field<fieldType>& source = blockSystem.source();
983 // Insert source from block system to rhs
984 for (label cmptI = 0; cmptI < nSrcCmpts; cmptI++)
986 scalarField sourceCmpt(source.component(cmptI));
990 b[cellI](localLoc1) += sourceCmpt[cellI];
1003 // Reset local locations for coupling contributions
1007 // Insert coupling contributions into block matrix
1008 forAll(bmesh, patchI)
1010 if (bmesh[patchI].coupled())
1012 typename CoeffField<blockType2>::squareTypeField& pcoupleUpper =
1013 A.coupleUpper()[patchI].asSquare();
1014 typename CoeffField<blockType2>::squareTypeField& pcoupleLower =
1015 A.coupleLower()[patchI].asSquare();
1017 const typename CoeffField<blockType1>::linearTypeField& bmcu =
1018 blockSystem.coupleUpper()[patchI].asLinear();
1019 const typename CoeffField<blockType1>::linearTypeField& bmcl =
1020 blockSystem.coupleLower()[patchI].asLinear();
1022 for (label cmptI = 0; cmptI < nCmpts; cmptI++)
1026 pcoupleUpper[faceI](localLoc1, localLoc2) +=
1027 bmcu[faceI].component(cmptI);
1028 pcoupleLower[faceI](localLoc1, localLoc2) +=
1029 bmcl[faceI].component(cmptI);
1042 // Reset local locations for other patches
1050 template<class blockType1, class blockType2, class fieldType>
1051 void insertBlockCoupling
1055 const BlockLduSystem<blockType1, fieldType>& blockSystem,
1056 BlockLduMatrix<blockType2>& A,
1057 Field<blockType2>& b,
1061 insertBlock(loc1, loc2, blockSystem, A, incFirst);
1062 insertBoundaryContributions(loc1, loc2, blockSystem, A, b, incFirst);
1066 template<class blockType>
1067 void insertEquationCoupling
1071 fvScalarMatrix& matrix,
1072 BlockLduMatrix<blockType>& A,
1081 "void insertEquationCoupling\n"
1083 " const label loc1,\n"
1084 " const label loc2,\n"
1085 " fvScalarMatrix& matrix\n"
1086 " BlockLduMatrix<blockType>& A\n"
1087 " Field<blockType>& b\n"
1089 ) << "Trying to insert coupling in the position where equation "
1090 << "should be since loc1 = loc2. Try using insertEquatiion "
1091 << "member function."
1092 << abort(FatalError);
1095 // Get references to fvScalarMatrix fields, updating boundary contributions
1096 scalarField& diag = matrix.D();
1097 scalarField& source = matrix.source();
1098 matrix.addBoundarySource(source, false);
1100 const scalarField& upper = matrix.upper();
1101 const scalarField& lower = matrix.lower();
1103 // Get references to ldu fields of A matrix always as square
1104 typename CoeffField<blockType>::squareTypeField& blockDiag =
1105 A.diag().asSquare();
1106 typename CoeffField<blockType>::squareTypeField& blockUpper =
1107 A.upper().asSquare();
1108 typename CoeffField<blockType>::squareTypeField& blockLower =
1109 A.lower().asSquare();
1113 blockDiag[cellI](loc1, loc2) += diag[cellI];
1114 b[cellI](loc1) += source[cellI];
1117 forAll(upper, faceI)
1119 blockUpper[faceI](loc1, loc2) += upper[faceI];
1120 blockLower[faceI](loc1, loc2) += lower[faceI];
1123 // Update coupling contributions
1124 const volScalarField& psi = matrix.psi();
1125 forAll(psi.boundaryField(), patchI)
1127 const fvPatchScalarField& pf = psi.boundaryField()[patchI];
1128 const fvPatch& patch = pf.patch();
1130 if (patch.coupled())
1132 const scalarField& icp = matrix.internalCoeffs()[patchI];
1133 const scalarField& bcp = matrix.boundaryCoeffs()[patchI];
1135 typename CoeffField<blockType>::squareTypeField& pcoupleUpper =
1136 A.coupleUpper()[patchI].asSquare();
1137 typename CoeffField<blockType>::squareTypeField& pcoupleLower =
1138 A.coupleLower()[patchI].asSquare();
1142 pcoupleUpper[faceI](loc1, loc2) = bcp[faceI];
1143 pcoupleLower[faceI](loc1, loc2) = icp[faceI];
1150 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1152 } // End namespace blockMatrixTools
1154 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1156 } // End namespace Foam
1158 // ************************************************************************* //