From f13764db1479e1a764b614fe281689c53e874c39 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 17 Mar 2016 11:51:50 +0000 Subject: [PATCH] Bugfix: Removed unused files --- .../blockMatrixTools/blockMatrixTools.C | 1158 ----------------- .../blockMatrixTools/blockMatrixTools.H | 288 ---- 2 files changed, 1446 deletions(-) delete mode 100644 src/finiteVolume/blockMatrixTools/blockMatrixTools.C delete mode 100644 src/finiteVolume/blockMatrixTools/blockMatrixTools.H diff --git a/src/finiteVolume/blockMatrixTools/blockMatrixTools.C b/src/finiteVolume/blockMatrixTools/blockMatrixTools.C deleted file mode 100644 index 92d6b8c28..000000000 --- a/src/finiteVolume/blockMatrixTools/blockMatrixTools.C +++ /dev/null @@ -1,1158 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 3.2 - \\ / A nd | Web: http://www.foam-extend.org - \\/ M anipulation | For copyright notice see file Copyright -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation, either version 3 of the License, or (at your - option) any later version. - - foam-extend is distributed in the hope that it will be useful, but - WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - General Public License for more details. - - You should have received a copy of the GNU General Public License - along with foam-extend. If not, see . - -\*---------------------------------------------------------------------------*/ - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace blockMatrixTools -{ - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -template -void blockInsert -( - const direction dir, - const scalarField& x, - Field& blockX -) -{ - forAll (x, i) - { - blockX[i](dir) = x[i]; - } -} - - -template -void blockAdd -( - const direction dir, - const scalarField& x, - Field& blockX -) -{ - forAll (x, i) - { - blockX[i](dir) += x[i]; - } -} - - -template -void blockRetrieve -( - const direction dir, - scalarField& x, - const Field& blockX -) -{ - forAll (x, i) - { - x[i] = blockX[i](dir); - } -} - - -template -void insertDiagSource -( - const direction dir, - const fvScalarMatrix& m, - BlockLduMatrix& blockM, - Field& blockB -) -{ - // Prepare the diagonal and source - - scalarField diag = m.diag(); - scalarField source = m.source(); - - // Add boundary source contribution - m.addBoundaryDiag(diag, 0); - m.addBoundarySource(source, false); - - if (blockM.diag().activeType() == blockCoeffBase::UNALLOCATED) - { - blockM.diag().asScalar() = diag; - } - else if - ( - blockM.diag().activeType() == blockCoeffBase::SCALAR - || blockM.diag().activeType() == blockCoeffBase::LINEAR - ) - { - typename CoeffField::linearTypeField& blockDiag = - blockM.diag().asLinear(); - - forAll (diag, i) - { - blockDiag[i](dir) = diag[i]; - } - } - else if (blockM.diag().activeType() == blockCoeffBase::SQUARE) - { - typename CoeffField::squareTypeField& blockDiag = - blockM.diag().asSquare(); - - forAll (diag, i) - { - blockDiag[i](dir, dir) = diag[i]; - } - } - - blockInsert(dir, source, blockB); -} - - -template -void insertUpperLower -( - const direction dir, - const fvScalarMatrix& m, - BlockLduMatrix& blockM -) -{ - if (m.diagonal()) - { - // Matrix for insertion is diagonal-only: nothing to do - return; - } - - if (m.symmetric() && blockM.symmetric()) - { - Info<< "Both m and blockM are symmetric: inserting only upper triangle" - << endl; - } - else - { - // Either scalar or block matrix is asymmetric: insert lower triangle - const scalarField& lower = m.lower(); - - if (blockM.lower().activeType() == blockCoeffBase::UNALLOCATED) - { - blockM.lower().asScalar() = lower; - } - else if - ( - blockM.lower().activeType() == blockCoeffBase::SCALAR - || blockM.lower().activeType() == blockCoeffBase::LINEAR - ) - { - typename CoeffField::linearTypeField& blockLower = - blockM.lower().asLinear(); - - forAll (lower, i) - { - blockLower[i](dir) = lower[i]; - } - } - else if (blockM.lower().activeType() == blockCoeffBase::SQUARE) - { - typename CoeffField::squareTypeField& blockLower = - blockM.lower().asSquare(); - - forAll (lower, i) - { - blockLower[i](dir, dir) = lower[i]; - } - } - } - - if (m.hasUpper()) - { - const scalarField& upper = m.upper(); - - if (blockM.upper().activeType() == blockCoeffBase::UNALLOCATED) - { - blockM.upper().asScalar() = upper; - } - else if - ( - blockM.upper().activeType() == blockCoeffBase::SCALAR - || blockM.upper().activeType() == blockCoeffBase::LINEAR - ) - { - typename CoeffField::linearTypeField& blockUpper = - blockM.upper().asLinear(); - - forAll (upper, i) - { - blockUpper[i](dir) = upper[i]; - } - } - else if (blockM.upper().activeType() == blockCoeffBase::SQUARE) - { - typename CoeffField::squareTypeField& blockUpper = - blockM.upper().asSquare(); - - forAll (upper, i) - { - blockUpper[i](dir, dir) = upper[i]; - } - } - } - else - { - FatalErrorIn - ( - "void insertUpperLower\n" - "(\n" - " const direction dir,\n" - " const fvScalarMatrix& m,\n" - " BlockLduMatrix& blockM\n" - ")" - ) << "Error in matrix insertion: problem with block structure" - << abort(FatalError); - } - - // Insert block interface fields - forAll(blockM.interfaces(), patchI) - { - if (blockM.interfaces().set(patchI)) - { - // Couple upper and lower - const scalarField& cUpper = m.boundaryCoeffs()[patchI]; - const scalarField& cLower = m.internalCoeffs()[patchI]; - - if - ( - blockM.coupleUpper()[patchI].activeType() - == blockCoeffBase::UNALLOCATED - ) - { - blockM.coupleUpper()[patchI].asScalar() = cUpper; - blockM.coupleLower()[patchI].asScalar() = cLower; - } - else if - ( - blockM.coupleUpper()[patchI].activeType() - == blockCoeffBase::SCALAR - || blockM.coupleUpper()[patchI].activeType() - == blockCoeffBase::LINEAR - ) - { - typename CoeffField::linearTypeField& blockUpper = - blockM.coupleUpper()[patchI].asLinear(); - - typename CoeffField::linearTypeField& blockLower = - blockM.coupleLower()[patchI].asLinear(); - - forAll (cUpper, i) - { - blockUpper[i](dir) = cUpper[i]; - blockLower[i](dir) = cLower[i]; - } - } - else if - ( - blockM.coupleUpper()[patchI].activeType() - == blockCoeffBase::SQUARE - ) - { - typename CoeffField::squareTypeField& blockUpper = - blockM.coupleUpper()[patchI].asSquare(); - - typename CoeffField::squareTypeField& blockLower = - blockM.coupleLower()[patchI].asSquare(); - - forAll (cUpper, i) - { - blockUpper[i](dir, dir) = cUpper[i]; - blockLower[i](dir, dir) = cLower[i]; - } - } - } - } -} - - -template -void insertEquation -( - const direction dir, - const fvScalarMatrix& m, - BlockLduMatrix& blockM, - Field& blockX, - Field& blockB -) -{ - insertDiagSource(dir, m, blockM, blockB); - insertUpperLower(dir, m, blockM); - blockInsert(dir, m.psi(), blockX); -} - - -template -void insertCouplingDiagSource -( - const direction dirI, - const direction dirJ, - const fvScalarMatrix& m, - BlockLduMatrix& blockM, - Field& blockB -) -{ - // Prepare the diagonal and source - - scalarField diag = m.diag(); - scalarField source = m.source(); - - // Add boundary source contribution - m.addBoundaryDiag(diag, 0); - m.addBoundarySource(source, false); - - - // Add off-diagonal block coefficients - typename CoeffField::squareTypeField& blockDiag = - blockM.diag().asSquare(); - - // Set off-diagonal coefficient - forAll (diag, i) - { - blockDiag[i](dirI, dirJ) = diag[i]; - } - - blockAdd(dirI, source, blockB); -} - - -template -void insertCouplingUpperLower -( - const direction dirI, - const direction dirJ, - const fvScalarMatrix& m, - BlockLduMatrix& blockM -) -{ - if (m.diagonal()) - { - // Matrix for insertion is diagonal-only: nothing to do - return; - } - - if (m.symmetric() && blockM.symmetric()) - { - Info<< "Both m and blockM are symmetric: inserting only upper triangle" - << endl; - } - else - { - // Either scalar or block matrix is asymmetric: insert lower triangle - const scalarField& lower = m.lower(); - - typename CoeffField::squareTypeField& blockLower = - blockM.lower().asSquare(); - - forAll (lower, i) - { - blockLower[i](dirI, dirJ) = lower[i]; - } - } - - if (m.hasUpper()) - { - const scalarField& upper = m.upper(); - - typename CoeffField::squareTypeField& blockUpper = - blockM.upper().asSquare(); - - forAll (upper, i) - { - blockUpper[i](dirI, dirJ) = upper[i]; - } - } - else - { - FatalErrorIn - ( - "void insertCouplingUpperLower\n" - "(\n" - " const direction dirI,\n" - " const direction dirJ,\n" - " const fvScalarMatrix& m,\n" - " BlockLduMatrix& blockM\n" - ")" - ) << "Error in matrix insertion: problem with block structure" - << abort(FatalError); - } - - // Insert block interface fields - forAll(blockM.interfaces(), patchI) - { - if (blockM.interfaces().set(patchI)) - { - // Couple upper and lower - const scalarField& cUpper = m.boundaryCoeffs()[patchI]; - const scalarField& cLower = m.internalCoeffs()[patchI]; - - typename CoeffField::squareTypeField& blockUpper = - blockM.coupleUpper()[patchI].asSquare(); - - typename CoeffField::squareTypeField& blockLower = - blockM.coupleLower()[patchI].asSquare(); - - forAll (cUpper, i) - { - blockUpper[i](dirI, dirJ) = cUpper[i]; - blockLower[i](dirI, dirJ) = cLower[i]; - } - } - } -} - - -template -void insertCoupling -( - const direction dirI, - const direction dirJ, - const fvScalarMatrix& m, - BlockLduMatrix& blockM, - Field& blockX, - Field& blockB -) -{ - insertCouplingDiagSource(dirI, dirJ, m, blockM, blockB); - insertCouplingUpperLower(dirI, dirJ, m, blockM); -} - - -template -void updateSourceCoupling -( - BlockLduMatrix& blockM, - Field& x, - Field& b -) -{ - // Eliminate off-diagonal block coefficients from the square diagonal - // With this change, the segregated matrix can be assembled with complete - // source terms and linearisation can be provided independently. - // Once the linearisation coefficients are set (off-diagonal entries - // in the square block matrix, they are multiplied by the current value - // of the field and subtracted from the source term - - if (blockM.diag().activeType() == blockCoeffBase::SQUARE) - { - typename CoeffField::squareTypeField& blockDiag = - blockM.diag().asSquare(); - - typename CoeffField::linearTypeField lf(blockDiag.size()); - typename CoeffField::squareTypeField sf(blockDiag.size()); - - // Expand and contract - - // Take out the diagonal entries from the square coefficient - contractLinear(lf, blockDiag); - - // Expand the diagonal for full square, with zeroes in the off-diagonal - expandLinear(sf, lf); - - // Subtract from the source the difference between the full block - // diagonal and the diagonal terms only - // Sign is the same as in the derivative - b += (blockDiag - sf) & x; - } -} - -template -void insertSolutionVector -( - const label loc, - const Field& xSingle, - Field& xBlock -) -{ - // Get number of field components and local copy of location, for - // consistency with member functions where locations need to be reset. - const label nCmpts = pTraits::nComponents; - label localLoc = loc; - - for (label cmptI = 0; cmptI < nCmpts; cmptI++) - { - scalarField xSingleCurr(xSingle.component(cmptI)); - - forAll (xSingleCurr, cellI) - { - xBlock[cellI](localLoc) = xSingleCurr[cellI]; - } - - localLoc++; - } -} - - -template -void retrieveSolution -( - const label loc, - Field& xSingle, - const Field& xBlock -) -{ - const label nCmpts = pTraits::nComponents; - label localLoc = loc; - - for (label cmptI = 0; cmptI < nCmpts; cmptI++) - { - scalarField xSingleCurr(xSingle.component(cmptI)); - - forAll (xSingleCurr, cellI) - { - xSingleCurr[cellI] = xBlock[cellI](localLoc); - } - - xSingle.replace(cmptI, xSingleCurr); - - localLoc++; - } -} - - -template -void insertDiagSource -( - const label loc, - fvMatrix& matrix, - BlockLduMatrix& A, - Field& b -) -{ - matrix.completeAssembly(); - - // Save a copy for different components - scalarField& diag = matrix.diag(); - scalarField saveDiag(diag); - - // Add source boundary contribution - Field& source = matrix.source(); - matrix.addBoundarySource(source, false); - - const label nCmpts = pTraits::nComponents; - label localLoc = loc; - - if - ( - // This is needed if the matrixType is , then you need to grab - // coeffs as linear. Consider doing a matrixType check also. - // VV, 17/March/2014 - A.diag().activeType() != blockCoeffBase::SQUARE - ) - { - typename CoeffField::linearTypeField& blockDiag = - A.diag().asLinear(); - - for (label cmptI = 0; cmptI < nCmpts; cmptI++) - { - matrix.addBoundaryDiag(diag, cmptI); - scalarField sourceCmpt(source.component(cmptI)); - -// FieldField bouCoeffsCmpt -// ( -// matrix.boundaryCoeffs().component(cmptI) -// ); - - // Possible problem for coupled non-aligned boundaries. - // VV, 14/May/2014. -// matrix.correctImplicitBoundarySource -// ( -// bouCoeffsCmpt, -// sourceCmpt, -// cmptI -// ); - - forAll (diag, cellI) - { - blockDiag[cellI](localLoc) = diag[cellI]; - b[cellI](localLoc) += sourceCmpt[cellI]; - } - - localLoc++; - - // Reset diagonal - diag = saveDiag; - } - } - else if (A.diag().activeType() == blockCoeffBase::SQUARE) - { - typename CoeffField::squareTypeField& blockDiag = - A.diag().asSquare(); - - for (label cmptI = 0; cmptI < nCmpts; cmptI++) - { - matrix.addBoundaryDiag(diag, cmptI); - scalarField sourceCmpt(source.component(cmptI)); - -// FieldField bouCoeffsCmpt -// ( -// matrix.boundaryCoeffs().component(cmptI) -// ); - -// matrix.correctImplicitBoundarySource -// ( -// bouCoeffsCmpt, -// sourceCmpt, -// cmptI -// ); - - forAll (diag, cellI) - { - blockDiag[cellI](localLoc, localLoc) = diag[cellI]; - b[cellI](localLoc) += sourceCmpt[cellI]; - } - - localLoc++; - - // Reset diagonal - diag = saveDiag; - } - } -} - - -template -void insertUpperLower -( - const label loc, - const fvMatrix& matrix, - BlockLduMatrix& A -) -{ - if (matrix.diagonal()) - { - // Matrix for insertion is diagonal-only: nothing to do - return; - } - - const label nCmpts = pTraits::nComponents; - label localLoc = loc; - - if (matrix.hasUpper()) - { - const scalarField& upper = matrix.upper(); - - if (A.upper().activeType() == blockCoeffBase::UNALLOCATED) - { - A.upper().asScalar() = upper; - } - else if - ( - A.upper().activeType() == blockCoeffBase::SCALAR - || A.upper().activeType() == blockCoeffBase::LINEAR - ) - { - typename CoeffField::linearTypeField& blockUpper = - A.upper().asLinear(); - - for (label cmptI = 0; cmptI < nCmpts; cmptI++) - { - forAll (upper, faceI) - { - blockUpper[faceI](localLoc) = upper[faceI]; - } - - localLoc++; - } - } - else if (A.upper().activeType() == blockCoeffBase::SQUARE) - { - typename CoeffField::squareTypeField& blockUpper = - A.upper().asSquare(); - - for (label cmptI = 0; cmptI < nCmpts; cmptI++) - { - forAll (upper, faceI) - { - blockUpper[faceI](localLoc, localLoc) = upper[faceI]; - } - - localLoc++; - } - } - } - else - { - FatalErrorIn - ( - "void insertUpperLower\n" - "(\n" - " const location loc,\n" - " const fvMatrix& matrix,\n" - " BlockLduMatrix& A\n" - ")" - ) << "Error in matrix insertion: problem with block structure." - << abort(FatalError); - } - - if (matrix.symmetric() && A.symmetric()) - { - Info<< "Both matrices are symmetric: inserting only upper triangle" - << endl; - } - else - { - // Reset localLoc - localLoc = loc; - - // Either scalar or block matrix is asymmetric: insert lower triangle - const scalarField& lower = matrix.lower(); - - if (A.lower().activeType() == blockCoeffBase::UNALLOCATED) - { - A.lower().asScalar() = lower; - } - else if - ( - A.lower().activeType() == blockCoeffBase::SCALAR - || A.lower().activeType() == blockCoeffBase::LINEAR - ) - { - typename CoeffField::linearTypeField& blockLower = - A.lower().asLinear(); - - for (label cmptI = 0; cmptI < nCmpts; cmptI++) - { - forAll (lower, faceI) - { - blockLower[faceI](localLoc) = lower[faceI]; - } - - localLoc++; - } - } - else if (A.lower().activeType() == blockCoeffBase::SQUARE) - { - typename CoeffField::squareTypeField& blockLower = - A.lower().asSquare(); - - for (label cmptI = 0; cmptI < nCmpts; cmptI++) - { - forAll (lower, faceI) - { - blockLower[faceI](localLoc, localLoc) = lower[faceI]; - } - - localLoc++; - } - } - } -} - - -template -void updateCouplingCoeffs -( - const label loc, - const fvMatrix& matrix, - BlockLduMatrix& A -) -{ - const label nCmpts = pTraits::nComponents; - label localLoc = loc; - - const GeometricField& psi = matrix.psi(); - forAll(psi.boundaryField(), patchI) - { - const fvPatchField& pf = psi.boundaryField()[patchI]; - const fvPatch& patch = pf.patch(); - - if (patch.coupled()) - { - const Field& icp = matrix.internalCoeffs()[patchI]; - const Field& bcp = matrix.boundaryCoeffs()[patchI]; - - if (A.coupleUpper()[patchI].activeType() != blockCoeffBase::SQUARE) - { - typename CoeffField::linearTypeField& pcoupleUpper = - A.coupleUpper()[patchI].asLinear(); - typename CoeffField::linearTypeField& pcoupleLower = - A.coupleLower()[patchI].asLinear(); - - for (label cmptI = 0; cmptI < nCmpts; cmptI++) - { - scalarField icpCmpt = icp.component(cmptI); - scalarField bcpCmpt = bcp.component(cmptI); - - forAll(pf, faceI) - { - pcoupleUpper[faceI](localLoc) = bcpCmpt[faceI]; - pcoupleLower[faceI](localLoc) = icpCmpt[faceI]; - } - - localLoc++; - } - - localLoc = loc; - } - else if - ( - A.coupleUpper()[patchI].activeType() == blockCoeffBase::SQUARE - ) - { - typename CoeffField::squareTypeField& pcoupleUpper = - A.coupleUpper()[patchI].asSquare(); - typename CoeffField::squareTypeField& pcoupleLower = - A.coupleLower()[patchI].asSquare(); - - for (label cmptI = 0; cmptI < nCmpts; cmptI++) - { - scalarField icpCmpt = icp.component(cmptI); - scalarField bcpCmpt = bcp.component(cmptI); - - forAll(pf, faceI) - { - pcoupleUpper[faceI](localLoc, localLoc) = - bcpCmpt[faceI]; - pcoupleLower[faceI](localLoc, localLoc) = - icpCmpt[faceI]; - } - - localLoc++; - } - - localLoc = loc; - } - } - } -} - - -template -void insertEquation -( - const label loc, - fvMatrix& matrix, - BlockLduMatrix& A, - Field& x, - Field& b -) -{ - insertDiagSource(loc, matrix, A, b); - insertUpperLower(loc, matrix, A); - updateCouplingCoeffs(loc, matrix, A); - insertSolutionVector(loc, matrix.psi().internalField(), x); -} - - -template -void insertBlock -( - const label loc1, - const label loc2, - const BlockLduSystem& blockSystem, - BlockLduMatrix& A, - const bool incFirst -) -{ - // Sanity checks - { - const label blockMatrixSize = pTraits::nComponents; - const label blockMatrixASize = pTraits::nComponents; - - if (blockMatrixSize > blockMatrixASize) - { - FatalErrorIn - ( - "void insertBlock\n" - "(\n" - " const label loc1,\n" - " const label loc2,\n" - " BlockLduMatrix& blockMatrix,\n" - " BlockLduMatrix& A\n" - ")" - ) << "Trying to insert a block matrix into smaller one." - << abort(FatalError); - } - - if (loc1 == loc2) - { - FatalErrorIn - ( - "void insertCoupling\n" - "(\n" - " const label loc1,\n" - " const label loc2,\n" - " BlockLduMatrix& blockMatrix,\n" - " BlockLduMatrix& A\n" - ")" - ) << "Trying to insert coupling in the position where equation " - << "should be, since loc1 = loc2. Try using insertEquatiion " - << "member function." - << abort(FatalError); - } - } - - const label nCmpts = pTraits::nComponents; - label localLoc1 = loc1; - label localLoc2 = loc2; - - // Get references to ldu fields of blockMatrix always as linear - const typename CoeffField::linearTypeField& bmd = - blockSystem.diag().asLinear(); - const typename CoeffField::linearTypeField& bmu = - blockSystem.upper().asLinear(); - const typename CoeffField::linearTypeField& bml = - blockSystem.lower().asLinear(); - - // Get references to ldu fields of A matrix always as square - typename CoeffField::squareTypeField& blockDiag = - A.diag().asSquare(); - typename CoeffField::squareTypeField& blockUpper = - A.upper().asSquare(); - typename CoeffField::squareTypeField& blockLower = - A.lower().asSquare(); - - // Insert blockMatrix that represents coupling into larger system matrix - for (label cmptI = 0; cmptI < nCmpts; cmptI++) - { - forAll(bmd, cellI) - { - blockDiag[cellI](localLoc1, localLoc2) += - bmd[cellI].component(cmptI); - } - - forAll(bmu, faceI) - { - blockUpper[faceI](localLoc1, localLoc2) += - bmu[faceI].component(cmptI); - blockLower[faceI](localLoc1, localLoc2) += - bml[faceI].component(cmptI); - } - - if (incFirst) - { - localLoc1++; - } - else - { - localLoc2++; - } - } -} - - -template -void insertBoundaryContributions -( - const label loc1, - const label loc2, - const BlockLduSystem& blockSystem, - BlockLduMatrix& A, - Field& b, - const bool incFirst -) -{ - // Need to get reference to fvMesh instead of lduMesh - const fvBoundaryMesh& bmesh = refCast(A.mesh()).boundary(); - - const label nCmpts = pTraits::nComponents; - const label nSrcCmpts = pTraits::nComponents; - label localLoc1 = loc1; - label localLoc2 = loc2; - - const Field& source = blockSystem.source(); - - // Insert source from block system to rhs - for (label cmptI = 0; cmptI < nSrcCmpts; cmptI++) - { - scalarField sourceCmpt(source.component(cmptI)); - - forAll(b, cellI) - { - b[cellI](localLoc1) += sourceCmpt[cellI]; - } - - if (incFirst) - { - localLoc1++; - } - else - { - localLoc2++; - } - } - - // Reset local locations for coupling contributions - localLoc1 = loc1; - localLoc2 = loc2; - - // Insert coupling contributions into block matrix - forAll(bmesh, patchI) - { - if (bmesh[patchI].coupled()) - { - typename CoeffField::squareTypeField& pcoupleUpper = - A.coupleUpper()[patchI].asSquare(); - typename CoeffField::squareTypeField& pcoupleLower = - A.coupleLower()[patchI].asSquare(); - - const typename CoeffField::linearTypeField& bmcu = - blockSystem.coupleUpper()[patchI].asLinear(); - const typename CoeffField::linearTypeField& bmcl = - blockSystem.coupleLower()[patchI].asLinear(); - - for (label cmptI = 0; cmptI < nCmpts; cmptI++) - { - forAll(bmcu, faceI) - { - pcoupleUpper[faceI](localLoc1, localLoc2) += - bmcu[faceI].component(cmptI); - pcoupleLower[faceI](localLoc1, localLoc2) += - bmcl[faceI].component(cmptI); - } - - if (incFirst) - { - localLoc1++; - } - else - { - localLoc2++; - } - } - - // Reset local locations for other patches - localLoc1 = loc1; - localLoc2 = loc2; - } - } -} - - -template -void insertBlockCoupling -( - const label loc1, - const label loc2, - const BlockLduSystem& blockSystem, - BlockLduMatrix& A, - Field& b, - const bool incFirst -) -{ - insertBlock(loc1, loc2, blockSystem, A, incFirst); - insertBoundaryContributions(loc1, loc2, blockSystem, A, b, incFirst); -} - - -template -void insertEquationCoupling -( - const label loc1, - const label loc2, - fvScalarMatrix& matrix, - BlockLduMatrix& A, - Field& b -) -{ - // Sanity check - if (loc1 == loc2) - { - FatalErrorIn - ( - "void insertEquationCoupling\n" - "(\n" - " const label loc1,\n" - " const label loc2,\n" - " fvScalarMatrix& matrix\n" - " BlockLduMatrix& A\n" - " Field& b\n" - ")" - ) << "Trying to insert coupling in the position where equation " - << "should be since loc1 = loc2. Try using insertEquatiion " - << "member function." - << abort(FatalError); - } - - // Get references to fvScalarMatrix fields, updating boundary contributions - scalarField& diag = matrix.D(); - scalarField& source = matrix.source(); - matrix.addBoundarySource(source, false); - - const scalarField& upper = matrix.upper(); - const scalarField& lower = matrix.lower(); - - // Get references to ldu fields of A matrix always as square - typename CoeffField::squareTypeField& blockDiag = - A.diag().asSquare(); - typename CoeffField::squareTypeField& blockUpper = - A.upper().asSquare(); - typename CoeffField::squareTypeField& blockLower = - A.lower().asSquare(); - - forAll(diag, cellI) - { - blockDiag[cellI](loc1, loc2) += diag[cellI]; - b[cellI](loc1) += source[cellI]; - } - - forAll(upper, faceI) - { - blockUpper[faceI](loc1, loc2) += upper[faceI]; - blockLower[faceI](loc1, loc2) += lower[faceI]; - } - - // Update coupling contributions - const volScalarField& psi = matrix.psi(); - forAll(psi.boundaryField(), patchI) - { - const fvPatchScalarField& pf = psi.boundaryField()[patchI]; - const fvPatch& patch = pf.patch(); - - if (patch.coupled()) - { - const scalarField& icp = matrix.internalCoeffs()[patchI]; - const scalarField& bcp = matrix.boundaryCoeffs()[patchI]; - - typename CoeffField::squareTypeField& pcoupleUpper = - A.coupleUpper()[patchI].asSquare(); - typename CoeffField::squareTypeField& pcoupleLower = - A.coupleLower()[patchI].asSquare(); - - forAll(pf, faceI) - { - pcoupleUpper[faceI](loc1, loc2) = bcp[faceI]; - pcoupleLower[faceI](loc1, loc2) = icp[faceI]; - } - } - } -} - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace blockMatrixTools - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// ************************************************************************* // diff --git a/src/finiteVolume/blockMatrixTools/blockMatrixTools.H b/src/finiteVolume/blockMatrixTools/blockMatrixTools.H deleted file mode 100644 index ab91e3fd9..000000000 --- a/src/finiteVolume/blockMatrixTools/blockMatrixTools.H +++ /dev/null @@ -1,288 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 3.2 - \\ / A nd | Web: http://www.foam-extend.org - \\/ M anipulation | For copyright notice see file Copyright -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation, either version 3 of the License, or (at your - option) any later version. - - foam-extend is distributed in the hope that it will be useful, but - WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - General Public License for more details. - - You should have received a copy of the GNU General Public License - along with foam-extend. If not, see . - -InNamespace - Foam:: - -Description - Block matrix insertion and retrieval tools. - Note: these will be obsolete since all of the functions are now member - functions in fvBlockMatrix class. All top level coupled solvers should - be rewritten using fvBlockMatrix - then we can get rid of this. - VV, 20/July/2014. - -SourceFiles - blockMatrixTools.C - -\*---------------------------------------------------------------------------*/ - - -#ifndef blockMatrixTools_H -#define blockMatrixTools_H - -#include "BlockLduSystem.H" -#include "fvMatrices.H" -#include "surfaceFieldsFwd.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -/*---------------------------------------------------------------------------*\ - Namespace blockMatrixTools functions Declaration -\*---------------------------------------------------------------------------*/ - -namespace blockMatrixTools -{ - // Field operations - - //- Insert field into block field - template - void blockInsert - ( - const direction dir, - const scalarField& x, - Field& blockX - ); - - //- Add field into block field - template - void blockAdd - ( - const direction dir, - const scalarField& x, - Field& blockX - ); - - //- Retrieve field from block field - template - void blockRetrieve - ( - const direction dir, - scalarField& x, - const Field& blockX - ); - - - // Diagonal block operations - - //- Insert matrix diagonal and source into the block system - template - void insertDiagSource - ( - const direction dir, - const fvScalarMatrix& m, - BlockLduMatrix& blockM, - Field& blockB - ); - - // Insert upper and lower coefficients matrix into the block system - template - void insertUpperLower - ( - const direction dir, - const fvScalarMatrix& m, - BlockLduMatrix& blockM - ); - - // Insert matrix into the block system - template - void insertEquation - ( - const direction dir, - const fvScalarMatrix& m, - BlockLduMatrix& blockM, - Field& blockX, - Field& blockB - ); - - - // Coupling block operations - - //- Insert coupling matrix diagonal and source into the block system - template - void insertCouplingDiagSource - ( - const direction dirI, - const direction dirJ, - const fvScalarMatrix& m, - BlockLduMatrix& blockM, - Field& blockB - ); - - // Insert coupling matrix into the block system - template - void insertCouplingUpperLower - ( - const direction dirI, - const direction dirJ, - const fvScalarMatrix& m, - BlockLduMatrix& blockM - ); - - // Insert coupling matrix into the block system - template - void insertCoupling - ( - const direction dir, - const fvScalarMatrix& m, - BlockLduMatrix& blockM, - Field& blockX, - Field& blockB - ); - - // Update coupling of block system - // Subtracts the block-coefficient coupling as specified by the user - // from the source, leaving the implicit update given by linearisation - template - void updateSourceCoupling - ( - BlockLduMatrix& blockM, - Field& x, - Field& b - ); - - //- Tools specific to pU coupled block matrix - - //- Insert field into block field to be solved for - template - void insertSolutionVector - ( - const label loc, - const Field& xSingle, - Field& xBlock - ); - - //- Retrieve solved field from block field - template - void retrieveSolution - ( - const label loc, - Field& xSingle, - const Field& xBlock - ); - - //- Insert matrix diagonal and source into the block system - template - void insertDiagSource - ( - const label loc, - fvMatrix& matrix, - BlockLduMatrix& A, - Field& b - ); - - // Insert upper and lower part into the block system - template - void insertUpperLower - ( - const label loc, - const fvMatrix& matrix, - BlockLduMatrix& A - ); - - // Update coupling coefficients in the block matrix - template - void updateCouplingCoeffs - ( - const label loc, - const fvMatrix& matrix, - BlockLduMatrix& A - ); - - // Insert matrix into the block system - template - void insertEquation - ( - const label loc, - fvMatrix& matrix, - BlockLduMatrix& A, - Field& x, - Field& b - ); - - // Insert diagonal, lower and upper of blockMatrix into A - template - void insertBlock - ( - const label loc1, - const label loc2, - const BlockLduSystem& blockMatrix, - BlockLduMatrix& A, - const bool incFirst - ); - - // Insert source and coupling coeffs of block system into A and b - template - void insertBoundaryContributions - ( - const label loc1, - const label loc2, - const BlockLduSystem& blockSystem, - BlockLduMatrix& A, - Field& b, - const bool incFirst - ); - - // Insert existing block system (obtained by implicit grad/div operator) - // into block system. - template - void insertBlockCoupling - ( - const label loc1, - const label loc2, - const BlockLduSystem& blockSystem, - BlockLduMatrix& A, - Field& b, - const bool incFirst - ); - - // Insert scalar equation coupling. Not tested: VV, 9/May/2014 - template - void insertEquationCoupling - ( - const label loc1, - const label loc2, - fvScalarMatrix& matrix, - BlockLduMatrix& A, - Field b - ); -} - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#ifdef NoRepository -# include "blockMatrixTools.C" -#endif - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* //