Forward compatibility: flex
[foam-extend-3.2.git] / src / finiteVolume / blockMatrixTools / blockMatrixTools.C
blob92d6b8c28cc00170b5b7957748f5ec1276b18c4c
1 /*---------------------------------------------------------------------------*\
2   =========                 |
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 -------------------------------------------------------------------------------
8 License
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 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
28 namespace Foam
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 namespace blockMatrixTools
36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
38 template<class BlockType>
39 void blockInsert
41     const direction dir,
42     const scalarField& x,
43     Field<BlockType>& blockX
46     forAll (x, i)
47     {
48         blockX[i](dir) = x[i];
49     }
53 template<class BlockType>
54 void blockAdd
56     const direction dir,
57     const scalarField& x,
58     Field<BlockType>& blockX
61     forAll (x, i)
62     {
63         blockX[i](dir) += x[i];
64     }
68 template<class BlockType>
69 void blockRetrieve
71     const direction dir,
72     scalarField& x,
73     const Field<BlockType>& blockX
76     forAll (x, i)
77     {
78         x[i] = blockX[i](dir);
79     }
83 template<class BlockType>
84 void insertDiagSource
86     const direction dir,
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)
102     {
103         blockM.diag().asScalar() = diag;
104     }
105     else if
106     (
107         blockM.diag().activeType() == blockCoeffBase::SCALAR
108      || blockM.diag().activeType() == blockCoeffBase::LINEAR
109     )
110     {
111         typename CoeffField<BlockType>::linearTypeField& blockDiag =
112             blockM.diag().asLinear();
114         forAll (diag, i)
115         {
116             blockDiag[i](dir) = diag[i];
117         }
118     }
119     else if (blockM.diag().activeType() == blockCoeffBase::SQUARE)
120     {
121         typename CoeffField<BlockType>::squareTypeField& blockDiag =
122             blockM.diag().asSquare();
124         forAll (diag, i)
125         {
126             blockDiag[i](dir, dir) = diag[i];
127         }
128     }
130     blockInsert(dir, source, blockB);
134 template<class BlockType>
135 void insertUpperLower
137     const direction dir,
138     const fvScalarMatrix& m,
139     BlockLduMatrix<BlockType>& blockM
142     if (m.diagonal())
143     {
144         // Matrix for insertion is diagonal-only: nothing to do
145         return;
146     }
148     if (m.symmetric() && blockM.symmetric())
149     {
150         Info<< "Both m and blockM are symmetric: inserting only upper triangle"
151             << endl;
152     }
153     else
154     {
155         // Either scalar or block matrix is asymmetric: insert lower triangle
156         const scalarField& lower = m.lower();
158         if (blockM.lower().activeType() == blockCoeffBase::UNALLOCATED)
159         {
160             blockM.lower().asScalar() = lower;
161         }
162         else if
163         (
164             blockM.lower().activeType() == blockCoeffBase::SCALAR
165          || blockM.lower().activeType() == blockCoeffBase::LINEAR
166         )
167         {
168             typename CoeffField<BlockType>::linearTypeField& blockLower =
169                 blockM.lower().asLinear();
171             forAll (lower, i)
172             {
173                 blockLower[i](dir) = lower[i];
174             }
175         }
176         else if (blockM.lower().activeType() == blockCoeffBase::SQUARE)
177         {
178             typename CoeffField<BlockType>::squareTypeField& blockLower =
179                 blockM.lower().asSquare();
181             forAll (lower, i)
182             {
183                 blockLower[i](dir, dir) = lower[i];
184             }
185         }
186     }
188     if (m.hasUpper())
189     {
190         const scalarField& upper = m.upper();
192         if (blockM.upper().activeType() == blockCoeffBase::UNALLOCATED)
193         {
194             blockM.upper().asScalar() = upper;
195         }
196         else if
197         (
198             blockM.upper().activeType() == blockCoeffBase::SCALAR
199          || blockM.upper().activeType() == blockCoeffBase::LINEAR
200         )
201         {
202             typename CoeffField<BlockType>::linearTypeField& blockUpper =
203                 blockM.upper().asLinear();
205             forAll (upper, i)
206             {
207                 blockUpper[i](dir) = upper[i];
208             }
209         }
210         else if (blockM.upper().activeType() == blockCoeffBase::SQUARE)
211         {
212             typename CoeffField<BlockType>::squareTypeField& blockUpper =
213                 blockM.upper().asSquare();
215             forAll (upper, i)
216             {
217                 blockUpper[i](dir, dir) = upper[i];
218             }
219         }
220     }
221     else
222     {
223         FatalErrorIn
224         (
225             "void insertUpperLower\n"
226             "(\n"
227             "    const direction dir,\n"
228             "    const fvScalarMatrix& m,\n"
229             "    BlockLduMatrix<BlockType>& blockM\n"
230             ")"
231         )   << "Error in matrix insertion: problem with block structure"
232             << abort(FatalError);
233     }
235     // Insert block interface fields
236     forAll(blockM.interfaces(), patchI)
237     {
238         if (blockM.interfaces().set(patchI))
239         {
240             // Couple upper and lower
241             const scalarField& cUpper = m.boundaryCoeffs()[patchI];
242             const scalarField& cLower = m.internalCoeffs()[patchI];
244             if
245             (
246                 blockM.coupleUpper()[patchI].activeType()
247              == blockCoeffBase::UNALLOCATED
248             )
249             {
250                 blockM.coupleUpper()[patchI].asScalar() = cUpper;
251                 blockM.coupleLower()[patchI].asScalar() = cLower;
252             }
253             else if
254             (
255                 blockM.coupleUpper()[patchI].activeType()
256              == blockCoeffBase::SCALAR
257              || blockM.coupleUpper()[patchI].activeType()
258              == blockCoeffBase::LINEAR
259             )
260             {
261                 typename CoeffField<BlockType>::linearTypeField& blockUpper =
262                     blockM.coupleUpper()[patchI].asLinear();
264                 typename CoeffField<BlockType>::linearTypeField& blockLower =
265                     blockM.coupleLower()[patchI].asLinear();
267                 forAll (cUpper, i)
268                 {
269                     blockUpper[i](dir) = cUpper[i];
270                     blockLower[i](dir) = cLower[i];
271                 }
272             }
273             else if
274             (
275                 blockM.coupleUpper()[patchI].activeType()
276              == blockCoeffBase::SQUARE
277             )
278             {
279                 typename CoeffField<BlockType>::squareTypeField& blockUpper =
280                     blockM.coupleUpper()[patchI].asSquare();
282                 typename CoeffField<BlockType>::squareTypeField& blockLower =
283                     blockM.coupleLower()[patchI].asSquare();
285                 forAll (cUpper, i)
286                 {
287                     blockUpper[i](dir, dir) = cUpper[i];
288                     blockLower[i](dir, dir) = cLower[i];
289                 }
290             }
291         }
292     }
296 template<class BlockType>
297 void insertEquation
299     const direction dir,
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
337     forAll (diag, i)
338     {
339         blockDiag[i](dirI, dirJ) = diag[i];
340     }
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
355     if (m.diagonal())
356     {
357         // Matrix for insertion is diagonal-only: nothing to do
358         return;
359     }
361     if (m.symmetric() && blockM.symmetric())
362     {
363         Info<< "Both m and blockM are symmetric: inserting only upper triangle"
364             << endl;
365     }
366     else
367     {
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();
374         forAll (lower, i)
375         {
376             blockLower[i](dirI, dirJ) = lower[i];
377         }
378     }
380     if (m.hasUpper())
381     {
382         const scalarField& upper = m.upper();
384         typename CoeffField<BlockType>::squareTypeField& blockUpper =
385             blockM.upper().asSquare();
387         forAll (upper, i)
388         {
389             blockUpper[i](dirI, dirJ) = upper[i];
390         }
391     }
392     else
393     {
394         FatalErrorIn
395         (
396             "void insertCouplingUpperLower\n"
397             "(\n"
398             "    const direction dirI,\n"
399             "    const direction dirJ,\n"
400             "    const fvScalarMatrix& m,\n"
401             "    BlockLduMatrix<BlockType>& blockM\n"
402             ")"
403         )   << "Error in matrix insertion: problem with block structure"
404             << abort(FatalError);
405     }
407     // Insert block interface fields
408     forAll(blockM.interfaces(), patchI)
409     {
410         if (blockM.interfaces().set(patchI))
411         {
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();
422             forAll (cUpper, i)
423             {
424                 blockUpper[i](dirI, dirJ) = cUpper[i];
425                 blockLower[i](dirI, dirJ) = cLower[i];
426             }
427         }
428     }
432 template<class BlockType>
433 void insertCoupling
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,
452     Field<BlockType>& x,
453     Field<BlockType>& b
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)
464     {
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;
483     }
486 template<class blockType, class fieldType>
487 void insertSolutionVector
489     const label loc,
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++)
500     {
501         scalarField xSingleCurr(xSingle.component(cmptI));
503         forAll (xSingleCurr, cellI)
504         {
505             xBlock[cellI](localLoc) = xSingleCurr[cellI];
506         }
508         localLoc++;
509     }
513 template<class blockType, class fieldType>
514 void retrieveSolution
516     const label loc,
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++)
525     {
526         scalarField xSingleCurr(xSingle.component(cmptI));
528         forAll (xSingleCurr, cellI)
529         {
530             xSingleCurr[cellI] = xBlock[cellI](localLoc);
531         }
533         xSingle.replace(cmptI, xSingleCurr);
535         localLoc++;
536     }
540 template<class blockType, class matrixType>
541 void insertDiagSource
543     const label loc,
544     fvMatrix<matrixType>& matrix,
545     BlockLduMatrix<blockType>& A,
546     Field<blockType>& b
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;
562     if
563     (
564         // This is needed if the matrixType is <vector>, then you need to grab
565         // coeffs as linear. Consider doing a matrixType check also.
566         // VV, 17/March/2014
567         A.diag().activeType() != blockCoeffBase::SQUARE
568     )
569     {
570         typename CoeffField<blockType>::linearTypeField& blockDiag =
571             A.diag().asLinear();
573         for (label cmptI = 0; cmptI < nCmpts; cmptI++)
574         {
575             matrix.addBoundaryDiag(diag, cmptI);
576             scalarField sourceCmpt(source.component(cmptI));
578 //            FieldField<Field, scalar> bouCoeffsCmpt
579 //            (
580 //                matrix.boundaryCoeffs().component(cmptI)
581 //            );
583             // Possible problem for coupled non-aligned boundaries.
584             // VV, 14/May/2014.
585 //            matrix.correctImplicitBoundarySource
586 //            (
587 //                bouCoeffsCmpt,
588 //                sourceCmpt,
589 //                cmptI
590 //            );
592             forAll (diag, cellI)
593             {
594                 blockDiag[cellI](localLoc) = diag[cellI];
595                 b[cellI](localLoc) += sourceCmpt[cellI];
596             }
598             localLoc++;
600             // Reset diagonal
601             diag = saveDiag;
602         }
603     }
604     else if (A.diag().activeType() == blockCoeffBase::SQUARE)
605     {
606         typename CoeffField<blockType>::squareTypeField& blockDiag =
607             A.diag().asSquare();
609         for (label cmptI = 0; cmptI < nCmpts; cmptI++)
610         {
611             matrix.addBoundaryDiag(diag, cmptI);
612             scalarField sourceCmpt(source.component(cmptI));
614 //            FieldField<Field, scalar> bouCoeffsCmpt
615 //            (
616 //                matrix.boundaryCoeffs().component(cmptI)
617 //            );
619 //            matrix.correctImplicitBoundarySource
620 //            (
621 //                bouCoeffsCmpt,
622 //                sourceCmpt,
623 //                cmptI
624 //            );
626             forAll (diag, cellI)
627             {
628                 blockDiag[cellI](localLoc, localLoc) = diag[cellI];
629                 b[cellI](localLoc) += sourceCmpt[cellI];
630             }
632             localLoc++;
634             // Reset diagonal
635             diag = saveDiag;
636         }
637     }
641 template<class blockType, class matrixType>
642 void insertUpperLower
644     const label loc,
645     const fvMatrix<matrixType>& matrix,
646     BlockLduMatrix<blockType>& A
649     if (matrix.diagonal())
650     {
651         // Matrix for insertion is diagonal-only: nothing to do
652         return;
653     }
655     const label nCmpts = pTraits<matrixType>::nComponents;
656     label localLoc = loc;
658     if (matrix.hasUpper())
659     {
660         const scalarField& upper = matrix.upper();
662         if (A.upper().activeType() == blockCoeffBase::UNALLOCATED)
663         {
664             A.upper().asScalar() = upper;
665         }
666         else if
667         (
668             A.upper().activeType() == blockCoeffBase::SCALAR
669          || A.upper().activeType() == blockCoeffBase::LINEAR
670         )
671         {
672             typename CoeffField<blockType>::linearTypeField& blockUpper =
673                 A.upper().asLinear();
675             for (label cmptI = 0; cmptI < nCmpts; cmptI++)
676             {
677                 forAll (upper, faceI)
678                 {
679                     blockUpper[faceI](localLoc) = upper[faceI];
680                 }
682                 localLoc++;
683             }
684         }
685         else if (A.upper().activeType() == blockCoeffBase::SQUARE)
686         {
687             typename CoeffField<blockType>::squareTypeField& blockUpper =
688                 A.upper().asSquare();
690             for (label cmptI = 0; cmptI < nCmpts; cmptI++)
691             {
692                 forAll (upper, faceI)
693                 {
694                     blockUpper[faceI](localLoc, localLoc) = upper[faceI];
695                 }
697                 localLoc++;
698             }
699         }
700     }
701     else
702     {
703         FatalErrorIn
704         (
705             "void insertUpperLower\n"
706             "(\n"
707             "    const location loc,\n"
708             "    const fvMatrix<matrixType>& matrix,\n"
709             "    BlockLduMatrix<blockType>& A\n"
710             ")"
711         )   << "Error in matrix insertion: problem with block structure."
712             << abort(FatalError);
713     }
715     if (matrix.symmetric() && A.symmetric())
716     {
717         Info<< "Both matrices are symmetric: inserting only upper triangle"
718             << endl;
719     }
720     else
721     {
722         // Reset localLoc
723         localLoc = loc;
725         // Either scalar or block matrix is asymmetric: insert lower triangle
726         const scalarField& lower = matrix.lower();
728         if (A.lower().activeType() == blockCoeffBase::UNALLOCATED)
729         {
730             A.lower().asScalar() = lower;
731         }
732         else if
733         (
734             A.lower().activeType() == blockCoeffBase::SCALAR
735          || A.lower().activeType() == blockCoeffBase::LINEAR
736         )
737         {
738             typename CoeffField<blockType>::linearTypeField& blockLower =
739                 A.lower().asLinear();
741             for (label cmptI = 0; cmptI < nCmpts; cmptI++)
742             {
743                 forAll (lower, faceI)
744                 {
745                     blockLower[faceI](localLoc) = lower[faceI];
746                 }
748                 localLoc++;
749             }
750         }
751         else if (A.lower().activeType() == blockCoeffBase::SQUARE)
752         {
753             typename CoeffField<blockType>::squareTypeField& blockLower =
754                 A.lower().asSquare();
756             for (label cmptI = 0; cmptI < nCmpts; cmptI++)
757             {
758                 forAll (lower, faceI)
759                 {
760                     blockLower[faceI](localLoc, localLoc) = lower[faceI];
761                 }
763                 localLoc++;
764             }
765         }
766     }
770 template<class blockType, class matrixType>
771 void updateCouplingCoeffs
773     const label loc,
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)
783     {
784         const fvPatchField<matrixType>& pf = psi.boundaryField()[patchI];
785         const fvPatch& patch = pf.patch();
787         if (patch.coupled())
788         {
789             const Field<matrixType>& icp = matrix.internalCoeffs()[patchI];
790             const Field<matrixType>& bcp = matrix.boundaryCoeffs()[patchI];
792             if (A.coupleUpper()[patchI].activeType() != blockCoeffBase::SQUARE)
793             {
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++)
800                 {
801                     scalarField icpCmpt = icp.component(cmptI);
802                     scalarField bcpCmpt = bcp.component(cmptI);
804                     forAll(pf, faceI)
805                     {
806                         pcoupleUpper[faceI](localLoc) = bcpCmpt[faceI];
807                         pcoupleLower[faceI](localLoc) = icpCmpt[faceI];
808                     }
810                     localLoc++;
811                 }
813                 localLoc = loc;
814             }
815             else if
816             (
817                 A.coupleUpper()[patchI].activeType() == blockCoeffBase::SQUARE
818             )
819             {
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++)
826                 {
827                     scalarField icpCmpt = icp.component(cmptI);
828                     scalarField bcpCmpt = bcp.component(cmptI);
830                     forAll(pf, faceI)
831                     {
832                         pcoupleUpper[faceI](localLoc, localLoc) =
833                             bcpCmpt[faceI];
834                         pcoupleLower[faceI](localLoc, localLoc) =
835                             icpCmpt[faceI];
836                     }
838                     localLoc++;
839                 }
841                 localLoc = loc;
842             }
843         }
844     }
848 template<class blockType, class matrixType>
849 void insertEquation
851     const label loc,
852     fvMatrix<matrixType>& matrix,
853     BlockLduMatrix<blockType>& A,
854     Field<blockType>& x,
855     Field<blockType>& b
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>
866 void insertBlock
868     const label loc1,
869     const label loc2,
870     const BlockLduSystem<blockType1, fieldType>& blockSystem,
871     BlockLduMatrix<blockType2>& A,
872     const bool incFirst
875     // Sanity checks
876     {
877         const label blockMatrixSize = pTraits<blockType1>::nComponents;
878         const label blockMatrixASize = pTraits<blockType2>::nComponents;
880         if (blockMatrixSize > blockMatrixASize)
881         {
882             FatalErrorIn
883             (
884                 "void insertBlock\n"
885                 "(\n"
886                 "    const label loc1,\n"
887                 "    const label loc2,\n"
888                 "    BlockLduMatrix<blockType1>& blockMatrix,\n"
889                 "    BlockLduMatrix<blockType2>& A\n"
890                 ")"
891             )   << "Trying to insert a block matrix into smaller one."
892                 << abort(FatalError);
893         }
895         if (loc1 == loc2)
896         {
897             FatalErrorIn
898             (
899                 "void insertCoupling\n"
900                 "(\n"
901                 "    const label loc1,\n"
902                 "    const label loc2,\n"
903                 "    BlockLduMatrix<blockType1>& blockMatrix,\n"
904                 "    BlockLduMatrix<blockType2>& A\n"
905                 ")"
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);
910         }
911     }
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 =
927          A.diag().asSquare();
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++)
935     {
936         forAll(bmd, cellI)
937         {
938             blockDiag[cellI](localLoc1, localLoc2) +=
939                 bmd[cellI].component(cmptI);
940         }
942         forAll(bmu, faceI)
943         {
944             blockUpper[faceI](localLoc1, localLoc2) +=
945                 bmu[faceI].component(cmptI);
946             blockLower[faceI](localLoc1, localLoc2) +=
947                 bml[faceI].component(cmptI);
948         }
950         if (incFirst)
951         {
952             localLoc1++;
953         }
954         else
955         {
956             localLoc2++;
957         }
958     }
962 template<class blockType1, class blockType2, class fieldType>
963 void insertBoundaryContributions
965     const label loc1,
966     const label loc2,
967     const BlockLduSystem<blockType1, fieldType>& blockSystem,
968     BlockLduMatrix<blockType2>& A,
969     Field<blockType2>& b,
970     const bool incFirst
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++)
985     {
986         scalarField sourceCmpt(source.component(cmptI));
988         forAll(b, cellI)
989         {
990             b[cellI](localLoc1) += sourceCmpt[cellI];
991         }
993         if (incFirst)
994         {
995             localLoc1++;
996         }
997         else
998         {
999             localLoc2++;
1000         }
1001     }
1003     // Reset local locations for coupling contributions
1004     localLoc1 = loc1;
1005     localLoc2 = loc2;
1007     // Insert coupling contributions into block matrix
1008     forAll(bmesh, patchI)
1009     {
1010         if (bmesh[patchI].coupled())
1011         {
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++)
1023                 {
1024                     forAll(bmcu, faceI)
1025                     {
1026                         pcoupleUpper[faceI](localLoc1, localLoc2) +=
1027                               bmcu[faceI].component(cmptI);
1028                         pcoupleLower[faceI](localLoc1, localLoc2) +=
1029                               bmcl[faceI].component(cmptI);
1030                     }
1032                     if (incFirst)
1033                     {
1034                         localLoc1++;
1035                     }
1036                     else
1037                     {
1038                         localLoc2++;
1039                     }
1040                 }
1042             // Reset local locations for other patches
1043             localLoc1 = loc1;
1044             localLoc2 = loc2;
1045         }
1046     }
1050 template<class blockType1, class blockType2, class fieldType>
1051 void insertBlockCoupling
1053     const label loc1,
1054     const label loc2,
1055     const BlockLduSystem<blockType1, fieldType>& blockSystem,
1056     BlockLduMatrix<blockType2>& A,
1057     Field<blockType2>& b,
1058     const bool incFirst
1061     insertBlock(loc1, loc2, blockSystem, A, incFirst);
1062     insertBoundaryContributions(loc1, loc2, blockSystem, A, b, incFirst);
1066 template<class blockType>
1067 void insertEquationCoupling
1069     const label loc1,
1070     const label loc2,
1071     fvScalarMatrix& matrix,
1072     BlockLduMatrix<blockType>& A,
1073     Field<blockType>& b
1076     // Sanity check
1077     if (loc1 == loc2)
1078     {
1079         FatalErrorIn
1080         (
1081             "void insertEquationCoupling\n"
1082             "(\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"
1088             ")"
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);
1093     }
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();
1111     forAll(diag, cellI)
1112     {
1113         blockDiag[cellI](loc1, loc2) += diag[cellI];
1114         b[cellI](loc1) += source[cellI];
1115     }
1117     forAll(upper, faceI)
1118     {
1119         blockUpper[faceI](loc1, loc2) += upper[faceI];
1120         blockLower[faceI](loc1, loc2) += lower[faceI];
1121     }
1123     // Update coupling contributions
1124     const volScalarField& psi = matrix.psi();
1125     forAll(psi.boundaryField(), patchI)
1126     {
1127         const fvPatchScalarField& pf = psi.boundaryField()[patchI];
1128         const fvPatch& patch = pf.patch();
1130         if (patch.coupled())
1131         {
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();
1140              forAll(pf, faceI)
1141              {
1142                 pcoupleUpper[faceI](loc1, loc2) = bcp[faceI];
1143                 pcoupleLower[faceI](loc1, loc2) = icp[faceI];
1144              }
1145         }
1146     }
1150 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1152 } // End namespace blockMatrixTools
1154 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1156 } // End namespace Foam
1158 // ************************************************************************* //