diff --git a/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructor.H b/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructor.H index 201ce66811..71ac334133 100644 --- a/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructor.H +++ b/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructor.H @@ -151,18 +151,20 @@ public: const IOobject& fieldIoObject ); - //- Reconstruct and write all volume fields + //- Reconstruct and write all/selected volume fields template void reconstructFvVolumeFields ( - const IOobjectList& objects + const IOobjectList& objects, + const HashSet& selectedFields ); - //- Reconstruct and write all volume fields + //- Reconstruct and write all/selected volume fields template void reconstructFvSurfaceFields ( - const IOobjectList& objects + const IOobjectList& objects, + const HashSet& selectedFields ); }; diff --git a/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructorReconstructFields.C b/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructorReconstructFields.C index 06f4f6a634..2ad0735ccf 100644 --- a/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructorReconstructFields.C +++ b/applications/utilities/parallelProcessing/reconstructPar/fvFieldReconstructorReconstructFields.C @@ -131,7 +131,7 @@ Foam::fvFieldReconstructor::reconstructFvVolumeField forAll(cp, faceI) { // Subtract one to take into account offsets for - // face direction. + // face direction. reverseAddressing[faceI] = cp[faceI] - 1 - curPatchStart; } @@ -151,7 +151,7 @@ Foam::fvFieldReconstructor::reconstructFvVolumeField forAll(cp, faceI) { // Subtract one to take into account offsets for - // face direction. + // face direction. label curF = cp[faceI] - 1; // Is the face on the boundary? @@ -282,7 +282,7 @@ Foam::fvFieldReconstructor::reconstructFvSurfaceField // It is necessary to create a copy of the addressing array to // take care of the face direction offset trick. - // + // { labelList curAddr(faceProcAddressing_[procI]); @@ -342,7 +342,7 @@ Foam::fvFieldReconstructor::reconstructFvSurfaceField forAll(cp, faceI) { // Subtract one to take into account offsets for - // face direction. + // face direction. reverseAddressing[faceI] = cp[faceI] - 1 - curPatchStart; } @@ -452,11 +452,12 @@ Foam::fvFieldReconstructor::reconstructFvSurfaceField } -// Reconstruct and write all volume fields +// Reconstruct and write all/selected volume fields template void Foam::fvFieldReconstructor::reconstructFvVolumeFields ( - const IOobjectList& objects + const IOobjectList& objects, + const HashSet& selectedFields ) { const word& fieldClassName = @@ -468,27 +469,29 @@ void Foam::fvFieldReconstructor::reconstructFvVolumeFields { Info<< " Reconstructing " << fieldClassName << "s\n" << endl; - for - ( - IOobjectList::const_iterator fieldIter = fields.begin(); - fieldIter != fields.end(); - ++fieldIter - ) + forAllConstIter(IOobjectList, fields, fieldIter) { - Info<< " " << fieldIter()->name() << endl; + if + ( + !selectedFields.size() + || selectedFields.found(fieldIter()->name()) + ) + { + Info<< " " << fieldIter()->name() << endl; - reconstructFvVolumeField(*fieldIter())().write(); + reconstructFvVolumeField(*fieldIter())().write(); + } } - Info<< endl; } } -// Reconstruct and write all surface fields +// Reconstruct and write all/selected surface fields template void Foam::fvFieldReconstructor::reconstructFvSurfaceFields ( - const IOobjectList& objects + const IOobjectList& objects, + const HashSet& selectedFields ) { const word& fieldClassName = @@ -500,18 +503,19 @@ void Foam::fvFieldReconstructor::reconstructFvSurfaceFields { Info<< " Reconstructing " << fieldClassName << "s\n" << endl; - for - ( - IOobjectList::const_iterator fieldIter = fields.begin(); - fieldIter != fields.end(); - ++fieldIter - ) + forAllConstIter(IOobjectList, fields, fieldIter) { - Info<< " " << fieldIter()->name() << endl; + if + ( + !selectedFields.size() + || selectedFields.found(fieldIter()->name()) + ) + { + Info<< " " << fieldIter()->name() << endl; - reconstructFvSurfaceField(*fieldIter())().write(); + reconstructFvSurfaceField(*fieldIter())().write(); + } } - Info<< endl; } } diff --git a/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C b/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C index db14efe3e5..50818f4db9 100644 --- a/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C +++ b/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C @@ -48,9 +48,17 @@ int main(int argc, char *argv[]) argList::noParallel(); timeSelector::addOptions(); # include "addRegionOption.H" + argList::validOptions.insert("fields", "\"(list of fields)\""); + # include "setRootCase.H" # include "createTime.H" + HashSet selectedFields; + if (args.options().found("fields")) + { + IStringStream(args.options()["fields"])() >> selectedFields; + } + // determine the processor count directly label nProcs = 0; while (dir(args.path()/(word("processor") + name(nProcs)))) @@ -184,13 +192,37 @@ int main(int argc, char *argv[]) procMeshes.boundaryProcAddressing() ); - fvReconstructor.reconstructFvVolumeFields(objects); - fvReconstructor.reconstructFvVolumeFields(objects); - fvReconstructor.reconstructFvVolumeFields(objects); - fvReconstructor.reconstructFvVolumeFields(objects); - fvReconstructor.reconstructFvVolumeFields(objects); + fvReconstructor.reconstructFvVolumeFields + ( + objects, + selectedFields + ); + fvReconstructor.reconstructFvVolumeFields + ( + objects, + selectedFields + ); + fvReconstructor.reconstructFvVolumeFields + ( + objects, + selectedFields + ); + fvReconstructor.reconstructFvVolumeFields + ( + objects, + selectedFields + ); + fvReconstructor.reconstructFvVolumeFields + ( + objects, + selectedFields + ); - fvReconstructor.reconstructFvSurfaceFields(objects); + fvReconstructor.reconstructFvSurfaceFields + ( + objects, + selectedFields + ); } else { diff --git a/applications/utilities/postProcessing/miscellaneous/execFlowFunctionObjects/execFlowFunctionObjects.C b/applications/utilities/postProcessing/miscellaneous/execFlowFunctionObjects/execFlowFunctionObjects.C index b848e330ea..3f7a99040b 100644 --- a/applications/utilities/postProcessing/miscellaneous/execFlowFunctionObjects/execFlowFunctionObjects.C +++ b/applications/utilities/postProcessing/miscellaneous/execFlowFunctionObjects/execFlowFunctionObjects.C @@ -139,7 +139,7 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh) IOobject LESPropertiesHeader ( - "RASProperties", + "LESProperties", runTime.constant(), mesh, IOobject::MUST_READ, diff --git a/src/ODE/ODE/ODE.H b/src/ODE/ODE/ODE.H index 6b3c69011c..fe57cb91f4 100644 --- a/src/ODE/ODE/ODE.H +++ b/src/ODE/ODE/ODE.H @@ -34,7 +34,7 @@ Description #define ODE_H #include "scalarField.H" -#include "Matrix.H" +#include "scalarMatrices.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -80,7 +80,7 @@ public: const scalar x, const scalarField& y, scalarField& dfdx, - Matrix& dfdy + scalarSquareMatrix& dfdy ) const = 0; }; diff --git a/src/ODE/ODESolvers/KRR4/KRR4.C b/src/ODE/ODESolvers/KRR4/KRR4.C index fa00919a42..ec3db04574 100644 --- a/src/ODE/ODESolvers/KRR4/KRR4.C +++ b/src/ODE/ODESolvers/KRR4/KRR4.C @@ -25,7 +25,6 @@ License \*---------------------------------------------------------------------------*/ #include "KRR4.H" -#include "simpleMatrix.H" #include "addToRunTimeSelectionTable.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -36,11 +35,11 @@ namespace Foam { addToRunTimeSelectionTable(ODESolver, KRR4, ODE); -const scalar +const scalar KRR4::safety = 0.9, KRR4::grow = 1.5, KRR4::pgrow = -0.25, KRR4::shrink = 0.5, KRR4::pshrink = (-1.0/3.0), KRR4::errcon = 0.1296; -const scalar +const scalar KRR4::gamma = 1.0/2.0, KRR4::a21 = 2.0, KRR4::a31 = 48.0/25.0, KRR4::a32 = 6.0/25.0, KRR4::c21 = -8.0, KRR4::c31 = 372.0/25.0, KRR4::c32 = 12.0/5.0, @@ -81,8 +80,8 @@ void Foam::KRR4::solve const ODE& ode, scalar& x, scalarField& y, - scalarField& dydx, - const scalar eps, + scalarField& dydx, + const scalar eps, const scalarField& yScale, const scalar hTry, scalar& hDid, @@ -109,14 +108,14 @@ void Foam::KRR4::solve a_[i][i] += 1.0/(gamma*h); } - simpleMatrix::LUDecompose(a_, pivotIndices_); + LUDecompose(a_, pivotIndices_); for (register label i=0; i::LUBacksubstitute(a_, pivotIndices_, g1_); + LUBacksubstitute(a_, pivotIndices_, g1_); for (register label i=0; i::LUBacksubstitute(a_, pivotIndices_, g2_); + LUBacksubstitute(a_, pivotIndices_, g2_); for (register label i=0; i::LUBacksubstitute(a_, pivotIndices_, g3_); + LUBacksubstitute(a_, pivotIndices_, g3_); for (register label i=0; i::LUBacksubstitute(a_, pivotIndices_, g4_); + LUBacksubstitute(a_, pivotIndices_, g4_); for (register label i=0; i dfdy_; - mutable Matrix a_; + mutable scalarSquareMatrix dfdy_; + mutable scalarSquareMatrix a_; mutable labelList pivotIndices_; static const int maxtry = 40; static const scalar safety, grow, pgrow, shrink, pshrink, errcon; - static const scalar + static const scalar gamma, a21, a31, a32, c21, c31, c32, c41, c42, c43, diff --git a/src/ODE/ODESolvers/SIBS/SIBS.H b/src/ODE/ODESolvers/SIBS/SIBS.H index 75383bc38a..164b92d299 100644 --- a/src/ODE/ODESolvers/SIBS/SIBS.H +++ b/src/ODE/ODESolvers/SIBS/SIBS.H @@ -60,8 +60,8 @@ class SIBS static const scalar safe1, safe2, redMax, redMin, scaleMX; mutable scalarField a_; - mutable Matrix alpha_; - mutable Matrix d_p_; + mutable scalarSquareMatrix alpha_; + mutable scalarSquareMatrix d_p_; mutable scalarField x_p_; mutable scalarField err_; @@ -69,7 +69,7 @@ class SIBS mutable scalarField ySeq_; mutable scalarField yErr_; mutable scalarField dfdx_; - mutable Matrix dfdy_; + mutable scalarSquareMatrix dfdy_; mutable label first_, kMax_, kOpt_; mutable scalar epsOld_, xNew_; @@ -82,7 +82,7 @@ void SIMPR const scalarField& y, const scalarField& dydx, const scalarField& dfdx, - const Matrix& dfdy, + const scalarSquareMatrix& dfdy, const scalar deltaX, const label nSteps, scalarField& yEnd @@ -97,7 +97,7 @@ void polyExtrapolate scalarField& yz, scalarField& dy, scalarField& x_p, - Matrix& d_p + scalarSquareMatrix& d_p ) const; diff --git a/src/ODE/ODESolvers/SIBS/SIMPR.C b/src/ODE/ODESolvers/SIBS/SIMPR.C index f540f26d11..671f1acec9 100644 --- a/src/ODE/ODESolvers/SIBS/SIMPR.C +++ b/src/ODE/ODESolvers/SIBS/SIMPR.C @@ -25,7 +25,6 @@ License \*---------------------------------------------------------------------------*/ #include "SIBS.H" -#include "simpleMatrix.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -36,7 +35,7 @@ void Foam::SIBS::SIMPR const scalarField& y, const scalarField& dydx, const scalarField& dfdx, - const Matrix& dfdy, + const scalarSquareMatrix& dfdy, const scalar deltaX, const label nSteps, scalarField& yEnd @@ -44,7 +43,7 @@ void Foam::SIBS::SIMPR { scalar h = deltaX/nSteps; - Matrix a(n_, n_); + scalarSquareMatrix a(n_); for (register label i=0; i::LUDecompose(a, pivotIndices); + LUDecompose(a, pivotIndices); for (register label i=0; i::LUBacksubstitute(a, pivotIndices, yEnd); + LUBacksubstitute(a, pivotIndices, yEnd); scalarField del(yEnd); scalarField ytemp(n_); @@ -83,7 +82,7 @@ void Foam::SIBS::SIMPR yEnd[i] = h*yEnd[i] - del[i]; } - simpleMatrix::LUBacksubstitute(a, pivotIndices, yEnd); + LUBacksubstitute(a, pivotIndices, yEnd); for (register label i=0; i::LUBacksubstitute(a, pivotIndices, yEnd); + LUBacksubstitute(a, pivotIndices, yEnd); for (register label i=0; i& d + scalarSquareMatrix& d ) const { label n = yz.size(); diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files index 7b137dedc0..cf4e331e25 100644 --- a/src/OpenFOAM/Make/files +++ b/src/OpenFOAM/Make/files @@ -177,8 +177,9 @@ dimensionedTypes/dimensionedTensor/dimensionedTensor.C matrices/solution/solution.C -scalarMatrix = matrices/scalarMatrix -$(scalarMatrix)/scalarMatrix.C +scalarMatrices = matrices/scalarMatrices +$(scalarMatrices)/scalarMatrices.C +$(scalarMatrices)/SVD/SVD.C LUscalarMatrix = matrices/LUscalarMatrix $(LUscalarMatrix)/LUscalarMatrix.C diff --git a/src/OpenFOAM/db/IOstreams/Pstreams/PstreamCommsStruct.C b/src/OpenFOAM/db/IOstreams/Pstreams/PstreamCommsStruct.C index 6b078aa25d..7a1a680643 100644 --- a/src/OpenFOAM/db/IOstreams/Pstreams/PstreamCommsStruct.C +++ b/src/OpenFOAM/db/IOstreams/Pstreams/PstreamCommsStruct.C @@ -22,8 +22,6 @@ License along with OpenFOAM; if not, write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA -Description - \*---------------------------------------------------------------------------*/ #include "Pstream.H" diff --git a/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C new file mode 100644 index 0000000000..a31058eace --- /dev/null +++ b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C @@ -0,0 +1,100 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM 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 2 of the License, or (at your + option) any later version. + + OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "DiagonalMatrix.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template +template +Foam::DiagonalMatrix::DiagonalMatrix(const Matrix& a) +: + List(min(a.n(), a.m())) +{ + forAll(*this, i) + { + this->operator[](i) = a[i][i]; + } +} + + +template +Foam::DiagonalMatrix::DiagonalMatrix(const label size) +: + List(size) +{} + + +template +Foam::DiagonalMatrix::DiagonalMatrix(const label size, const Type& val) +: + List(size, val) +{} + + +template +Foam::DiagonalMatrix& Foam::DiagonalMatrix::invert() +{ + forAll(*this, i) + { + Type x = this->operator[](i); + if (mag(x) < VSMALL) + { + this->operator[](i) = Type(0); + } + else + { + this->operator[](i) = Type(1)/x; + } + } + + return this; +} + + +template +Foam::DiagonalMatrix Foam::inv(const DiagonalMatrix& A) +{ + DiagonalMatrix Ainv = A; + + forAll(A, i) + { + Type x = A[i]; + if (mag(x) < VSMALL) + { + Ainv[i] = Type(0); + } + else + { + Ainv[i] = Type(1)/x; + } + } + + return Ainv; +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H new file mode 100644 index 0000000000..ae8ceffc1c --- /dev/null +++ b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H @@ -0,0 +1,103 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM 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 2 of the License, or (at your + option) any later version. + + OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::DiagonalMatrix + +Description + DiagonalMatrix is a 2D diagonal matrix of objects + of type Type, size nxn + +SourceFiles + DiagonalMatrix.C + +\*---------------------------------------------------------------------------*/ + +#ifndef DiagonalMatrix_H +#define DiagonalMatrix_H + +#include "List.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * Class Forward declaration * * * * * * * * * * * // + +template class Matrix; + +/*---------------------------------------------------------------------------*\ + Class DiagonalMatrix Declaration +\*---------------------------------------------------------------------------*/ + +template +class DiagonalMatrix +: + public List +{ +public: + + // Constructors + + //- Construct from diagonal component of a Matrix + template + DiagonalMatrix(const Matrix&); + + //- Construct empty from size + DiagonalMatrix(const label size); + + //- Construct from size and a value + DiagonalMatrix(const label, const Type&); + + + // Member functions + + //- Invert the diaganol matrix and return itself + DiagonalMatrix& invert(); +}; + + +// Global functions + +//- Return the diagonal Matrix inverse +template +DiagonalMatrix inv(const DiagonalMatrix&); + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "DiagonalMatrix.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C index 5de8693f2f..3e8fc3fe0f 100644 --- a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C +++ b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C @@ -31,9 +31,9 @@ License // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::LUscalarMatrix::LUscalarMatrix(const Matrix& matrix) +Foam::LUscalarMatrix::LUscalarMatrix(const scalarSquareMatrix& matrix) : - scalarMatrix(matrix), + scalarSquareMatrix(matrix), pivotIndices_(n()) { LUDecompose(*this, pivotIndices_); @@ -101,7 +101,7 @@ Foam::LUscalarMatrix::LUscalarMatrix nCells += lduMatrices[i].size(); } - Matrix m(nCells, nCells, 0.0); + scalarSquareMatrix m(nCells, nCells, 0.0); transfer(m); convert(lduMatrices); } @@ -109,7 +109,7 @@ Foam::LUscalarMatrix::LUscalarMatrix else { label nCells = ldum.lduAddr().size(); - Matrix m(nCells, nCells, 0.0); + scalarSquareMatrix m(nCells, nCells, 0.0); transfer(m); convert(ldum, interfaceCoeffs, interfaces); } diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H index d5e79faf19..2255cee236 100644 --- a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H +++ b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H @@ -36,7 +36,7 @@ SourceFiles #ifndef LUscalarMatrix_H #define LUscalarMatrix_H -#include "scalarMatrix.H" +#include "scalarMatrices.H" #include "labelList.H" #include "FieldField.H" #include "lduInterfaceFieldPtrsList.H" @@ -55,7 +55,7 @@ class procLduMatrix; class LUscalarMatrix : - public scalarMatrix + public scalarSquareMatrix { // Private data @@ -78,7 +78,7 @@ class LUscalarMatrix void convert(const PtrList& lduMatrices); - //- Print the ratio of the mag-sum of the off-diagonal coefficients + //- Print the ratio of the mag-sum of the off-diagonal coefficients // to the mag-diagonal void printDiagonalDominance() const; @@ -87,8 +87,8 @@ public: // Constructors - //- Construct from Matrix and perform LU decomposition - LUscalarMatrix(const Matrix&); + //- Construct from scalarSquareMatrix and perform LU decomposition + LUscalarMatrix(const scalarSquareMatrix&); //- Construct from lduMatrix and perform LU decomposition LUscalarMatrix diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C index cbb59e2c86..ac0bc17207 100644 --- a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C +++ b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C @@ -28,16 +28,16 @@ License // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -template -void Foam::LUscalarMatrix::solve(Field& sourceSol) const +template +void Foam::LUscalarMatrix::solve(Field& sourceSol) const { if (Pstream::parRun()) { - Field completeSourceSol(n()); + Field completeSourceSol(n()); if (Pstream::master()) { - typename Field::subField + typename Field::subField ( completeSourceSol, sourceSol.size() @@ -58,7 +58,7 @@ void Foam::LUscalarMatrix::solve(Field& sourceSol) const ( &(completeSourceSol[procOffsets_[slave]]) ), - (procOffsets_[slave + 1] - procOffsets_[slave])*sizeof(T) + (procOffsets_[slave + 1] - procOffsets_[slave])*sizeof(Type) ); } } @@ -77,7 +77,7 @@ void Foam::LUscalarMatrix::solve(Field& sourceSol) const { LUBacksubstitute(*this, pivotIndices_, completeSourceSol); - sourceSol = typename Field::subField + sourceSol = typename Field::subField ( completeSourceSol, sourceSol.size() @@ -98,7 +98,7 @@ void Foam::LUscalarMatrix::solve(Field& sourceSol) const ( &(completeSourceSol[procOffsets_[slave]]) ), - (procOffsets_[slave + 1] - procOffsets_[slave])*sizeof(T) + (procOffsets_[slave + 1] - procOffsets_[slave])*sizeof(Type) ); } } diff --git a/src/OpenFOAM/containers/Matrix/Matrix.C b/src/OpenFOAM/matrices/Matrix/Matrix.C similarity index 59% rename from src/OpenFOAM/containers/Matrix/Matrix.C rename to src/OpenFOAM/matrices/Matrix/Matrix.C index 733e04034e..577cbcabcb 100644 --- a/src/OpenFOAM/containers/Matrix/Matrix.C +++ b/src/OpenFOAM/matrices/Matrix/Matrix.C @@ -28,13 +28,13 @@ License // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // -template -void Foam::Matrix::allocate() +template +void Foam::Matrix::allocate() { if (n_ && m_) { - v_ = new T*[n_]; - v_[0] = new T[n_*m_]; + v_ = new Type*[n_]; + v_[0] = new Type[n_*m_]; for (register label i=1; i::allocate() // * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * * // -template -Foam::Matrix::~Matrix() +template +Foam::Matrix::~Matrix() { if (v_) { - delete[] (v_[0]); + delete[] (v_[0]); delete[] v_; } } @@ -59,16 +59,16 @@ Foam::Matrix::~Matrix() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -template -const Foam::Matrix& Foam::Matrix::null() +template +const Foam::Matrix& Foam::Matrix::null() { - Matrix* nullPtr = reinterpret_cast*>(NULL); + Matrix* nullPtr = reinterpret_cast*>(NULL); return *nullPtr; } -template -Foam::Matrix::Matrix(const label n, const label m) +template +Foam::Matrix::Matrix(const label n, const label m) : n_(n), m_(m), @@ -76,7 +76,7 @@ Foam::Matrix::Matrix(const label n, const label m) { if (n_ < 0 || m_ < 0) { - FatalErrorIn("Matrix::Matrix(const label n, const label m)") + FatalErrorIn("Matrix::Matrix(const label n, const label m)") << "bad n, m " << n_ << ", " << m_ << abort(FatalError); } @@ -85,8 +85,8 @@ Foam::Matrix::Matrix(const label n, const label m) } -template -Foam::Matrix::Matrix(const label n, const label m, const T& a) +template +Foam::Matrix::Matrix(const label n, const label m, const Type& a) : n_(n), m_(m), @@ -96,7 +96,7 @@ Foam::Matrix::Matrix(const label n, const label m, const T& a) { FatalErrorIn ( - "Matrix::Matrix(const label n, const label m, const T&)" + "Matrix::Matrix(const label n, const label m, const T&)" ) << "bad n, m " << n_ << ", " << m_ << abort(FatalError); } @@ -105,7 +105,7 @@ Foam::Matrix::Matrix(const label n, const label m, const T& a) if (v_) { - T* v = v_[0]; + Type* v = v_[0]; label nm = n_*m_; @@ -117,8 +117,8 @@ Foam::Matrix::Matrix(const label n, const label m, const T& a) } -template -Foam::Matrix::Matrix(const Matrix& a) +template +Foam::Matrix::Matrix(const Matrix& a) : n_(a.n_), m_(a.m_), @@ -127,8 +127,8 @@ Foam::Matrix::Matrix(const Matrix& a) if (a.v_) { allocate(); - T* v = v_[0]; - const T* av = a.v_[0]; + Type* v = v_[0]; + const Type* av = a.v_[0]; label nm = n_*m_; for (register label i=0; i::Matrix(const Matrix& a) } } - -template -void Foam::Matrix::clear() + +template +void Foam::Matrix::clear() { if (v_) { - delete[] (v_[0]); + delete[] (v_[0]); delete[] v_; } n_ = 0; @@ -153,8 +153,8 @@ void Foam::Matrix::clear() } -template -void Foam::Matrix::transfer(Matrix& a) +template +void Foam::Matrix::transfer(Matrix& a) { clear(); @@ -169,14 +169,32 @@ void Foam::Matrix::transfer(Matrix& a) } +template +Form Foam::Matrix::T() const +{ + const Matrix& A = *this; + Form At(m(), n()); + + for (register label i=0; i -void Foam::Matrix::operator=(const T& t) +template +void Foam::Matrix::operator=(const Type& t) { if (v_) { - T* v = v_[0]; + Type* v = v_[0]; label nm = n_*m_; for (register label i=0; i::operator=(const T& t) // Assignment operator. Takes linear time. -template -void Foam::Matrix::operator=(const Matrix& a) +template +void Foam::Matrix::operator=(const Matrix& a) { if (this == &a) { - FatalErrorIn("Matrix::operator=(const Matrix&)") + FatalErrorIn("Matrix::operator=(const Matrix&)") << "attempted assignment to self" << abort(FatalError); } @@ -204,12 +222,12 @@ void Foam::Matrix::operator=(const Matrix& a) n_ = a.n_; m_ = a.m_; allocate(); - } + } if (v_) { - T* v = v_[0]; - const T* av = a.v_[0]; + Type* v = v_[0]; + const Type* av = a.v_[0]; label nm = n_*m_; for (register label i=0; i::operator=(const Matrix& a) // * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * // -template -const T& Foam::max(const Matrix& a) +template +const Type& Foam::max(const Matrix& a) { label nm = a.n_*a.m_; if (nm) { label curMaxI = 0; - const T* v = a.v_[0]; + const Type* v = a.v_[0]; for (register label i=1; i& a) } else { - FatalErrorIn("max(const Matrix&)") + FatalErrorIn("max(const Matrix&)") << "matrix is empty" << abort(FatalError); @@ -254,15 +272,15 @@ const T& Foam::max(const Matrix& a) } -template -const T& Foam::min(const Matrix& a) +template +const Type& Foam::min(const Matrix& a) { label nm = a.n_*a.m_; if (nm) { label curMinI = 0; - const T* v = a.v_[0]; + const Type* v = a.v_[0]; for (register label i=1; i& a) } else { - FatalErrorIn("min(const Matrix&)") + FatalErrorIn("min(const Matrix&)") << "matrix is empty" << abort(FatalError); @@ -288,15 +306,15 @@ const T& Foam::min(const Matrix& a) // * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * // -template -Foam::Matrix Foam::operator-(const Matrix& a) +template +Form Foam::operator-(const Matrix& a) { - Matrix na(a.n_, a.m_); + Form na(a.n_, a.m_); if (a.n_ && a.m_) { - T* nav = na.v_[0]; - const T* av = a.v_[0]; + Type* nav = na.v_[0]; + const Type* av = a.v_[0]; label nm = a.n_*a.m_; for (register label i=0; i Foam::operator-(const Matrix& a) } -template -Foam::Matrix Foam::operator+(const Matrix& a, const Matrix& b) +template +Form Foam::operator+(const Matrix& a, const Matrix& b) { if (a.n_ != b.n_) { - FatalErrorIn("Matrix::operator+(const Matrix&, const Matrix&)") - << "attempted add matrices with different number of rows: " + FatalErrorIn + ( + "Matrix::operator+(const Matrix&, const Matrix&)" + ) << "attempted add matrices with different number of rows: " << a.n_ << ", " << b.n_ << abort(FatalError); } if (a.m_ != b.m_) { - FatalErrorIn("Matrix::operator+(const Matrix&, const Matrix&)") - << "attempted add matrices with different number of columns: " + FatalErrorIn + ( + "Matrix::operator+(const Matrix&, const Matrix&)" + ) << "attempted add matrices with different number of columns: " << a.m_ << ", " << b.m_ << abort(FatalError); } - Matrix ab(a.n_, a.m_); + Form ab(a.n_, a.m_); - T* abv = ab.v_[0]; - const T* av = a.v_[0]; - const T* bv = b.v_[0]; + Type* abv = ab.v_[0]; + const Type* av = a.v_[0]; + const Type* bv = b.v_[0]; label nm = a.n_*a.m_;; for (register label i=0; i Foam::operator+(const Matrix& a, const Matrix& b) } -template -Foam::Matrix Foam::operator-(const Matrix& a, const Matrix& b) +template +Form Foam::operator-(const Matrix& a, const Matrix& b) { if (a.n_ != b.n_) { - FatalErrorIn("Matrix::operator-(const Matrix&, const Matrix&)") - << "attempted add matrices with different number of rows: " + FatalErrorIn + ( + "Matrix::operator-(const Matrix&, const Matrix&)" + ) << "attempted add matrices with different number of rows: " << a.n_ << ", " << b.n_ << abort(FatalError); } if (a.m_ != b.m_) { - FatalErrorIn("Matrix::operator-(const Matrix&, const Matrix&)") - << "attempted add matrices with different number of columns: " + FatalErrorIn + ( + "Matrix::operator-(const Matrix&, const Matrix&)" + ) << "attempted add matrices with different number of columns: " << a.m_ << ", " << b.m_ << abort(FatalError); } - Matrix ab(a.n_, a.m_); + Form ab(a.n_, a.m_); - T* abv = ab.v_[0]; - const T* av = a.v_[0]; - const T* bv = b.v_[0]; + Type* abv = ab.v_[0]; + const Type* av = a.v_[0]; + const Type* bv = b.v_[0]; label nm = a.n_*a.m_;; for (register label i=0; i Foam::operator-(const Matrix& a, const Matrix& b) } -template -Foam::Matrix Foam::operator*(const scalar s, const Matrix& a) +template +Form Foam::operator*(const scalar s, const Matrix& a) { - Matrix sa(a.n_, a.m_); + Form sa(a.n_, a.m_); if (a.n_ && a.m_) { - T* sav = sa.v_[0]; - const T* av = a.v_[0]; + Type* sav = sa.v_[0]; + const Type* av = a.v_[0]; - label nm = a.n_*a.m_;; + label nm = a.n_*a.m_; for (register label i=0; i class Matrix; +template class Matrix; -template const T& max(const Matrix&); -template const T& min(const Matrix&); +template Istream& operator>> +( + Istream&, + Matrix& +); -template Matrix operator-(const Matrix&); -template Matrix operator+(const Matrix&, const Matrix&); -template Matrix operator-(const Matrix&, const Matrix&); -template Matrix operator*(const scalar, const Matrix&); - -template Istream& operator>>(Istream&, Matrix&); -template Ostream& operator<<(Ostream&, const Matrix&); +template Ostream& operator<< +( + Ostream&, + const Matrix& +); /*---------------------------------------------------------------------------*\ Class Matrix Declaration \*---------------------------------------------------------------------------*/ -template +template class Matrix { // Private data @@ -78,7 +79,7 @@ class Matrix label n_, m_; //- Row pointers - T** __restrict__ v_; + Type** __restrict__ v_; //- Allocate the storage for the row-pointers and the data // and set the row pointers @@ -97,16 +98,16 @@ public: //- Construct with given number of rows and columns // and value for all elements. - Matrix(const label n, const label m, const T&); + Matrix(const label n, const label m, const Type&); //- Copy constructor. - Matrix(const Matrix&); + Matrix(const Matrix&); //- Construct from Istream. Matrix(Istream&); //- Clone - inline autoPtr > clone() const; + inline autoPtr > clone() const; // Destructor @@ -117,7 +118,7 @@ public: // Member functions //- Return a null Matrix - static const Matrix& null(); + static const Matrix& null(); // Access @@ -148,48 +149,64 @@ public: //- Transfer the contents of the argument Matrix into this Matrix // and annull the argument Matrix. - void transfer(Matrix&); + void transfer(Matrix&); + + + //- Return the transpose of the matrix + Form T() const; // Member operators //- Return subscript-checked element of Matrix. - inline T* operator[](const label); + inline Type* operator[](const label); //- Return subscript-checked element of constant Matrix. - inline const T* operator[](const label) const; + inline const Type* operator[](const label) const; //- Assignment operator. Takes linear time. - void operator=(const Matrix&); + void operator=(const Matrix&); //- Assignment of all entries to the given value - void operator=(const T&); - - - // Friend functions - - friend const T& max(const Matrix&); - friend const T& min(const Matrix&); - - - // Friend operators - - friend Matrix operator-(const Matrix&); - friend Matrix operator+(const Matrix&, const Matrix&); - friend Matrix operator-(const Matrix&, const Matrix&); - friend Matrix operator*(const scalar, const Matrix&); + void operator=(const Type&); // IOstream operators //- Read Matrix from Istream, discarding contents of existing Matrix. - friend Istream& operator>> (Istream&, Matrix&); + friend Istream& operator>> (Istream&, Matrix&); // Write Matrix to Ostream. - friend Ostream& operator<< (Ostream&, const Matrix&); + friend Ostream& operator<< (Ostream&, const Matrix&); }; +// Global functions and operators + +template const Type& max(const Matrix&); +template const Type& min(const Matrix&); + +template Form operator-(const Matrix&); + +template Form operator+ +( + const Matrix&, + const Matrix& +); + +template Form operator- +( + const Matrix&, + const Matrix& +); + +template Form operator* +( + const scalar, + const Matrix& +); + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam diff --git a/src/OpenFOAM/containers/Matrix/MatrixI.H b/src/OpenFOAM/matrices/Matrix/MatrixI.H similarity index 66% rename from src/OpenFOAM/containers/Matrix/MatrixI.H rename to src/OpenFOAM/matrices/Matrix/MatrixI.H index d0e3fe8fdf..8d758cf5b1 100644 --- a/src/OpenFOAM/containers/Matrix/MatrixI.H +++ b/src/OpenFOAM/matrices/Matrix/MatrixI.H @@ -26,8 +26,8 @@ License // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -template -inline Foam::Matrix::Matrix() +template +inline Foam::Matrix::Matrix() : n_(0), m_(0), @@ -35,71 +35,67 @@ inline Foam::Matrix::Matrix() {} -template -inline Foam::autoPtr > Foam::Matrix::clone() const +template +inline Foam::autoPtr > Foam::Matrix::clone() const { - return autoPtr >(new Matrix(*this)); + return autoPtr >(new Matrix(*this)); } // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // //- Return the number of rows -template -inline Foam::label Foam::Matrix::n() const +template +inline Foam::label Foam::Matrix::n() const { return n_; } -//- Return the number of columns -template -inline Foam::label Foam::Matrix::m() const +template +inline Foam::label Foam::Matrix::m() const { return m_; } -//- Return the number of columns -template -inline Foam::label Foam::Matrix::size() const +template +inline Foam::label Foam::Matrix::size() const { return n_*m_; } -// Check index i is within valid range (0 ... n-1). -template -inline void Foam::Matrix::checki(const label i) const +template +inline void Foam::Matrix::checki(const label i) const { if (!n_) { - FatalErrorIn("Matrix::checki(const label)") + FatalErrorIn("Matrix::checki(const label)") << "attempt to access element from zero sized row" << abort(FatalError); } else if (i<0 || i>=n_) { - FatalErrorIn("Matrix::checki(const label)") + FatalErrorIn("Matrix::checki(const label)") << "index " << i << " out of range 0 ... " << n_-1 << abort(FatalError); } } -// Check index j is within valid range (0 ... n-1). -template -inline void Foam::Matrix::checkj(const label j) const +template +inline void Foam::Matrix::checkj(const label j) const { if (!m_) { - FatalErrorIn("Matrix::checkj(const label)") + FatalErrorIn("Matrix::checkj(const label)") << "attempt to access element from zero sized column" << abort(FatalError); } else if (j<0 || j>=m_) { - FatalErrorIn("Matrix::checkj(const label)") + FatalErrorIn("Matrix::checkj(const label)") << "index " << j << " out of range 0 ... " << m_-1 << abort(FatalError); } @@ -108,9 +104,8 @@ inline void Foam::Matrix::checkj(const label j) const // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // -// Return subscript-checked element access -template -inline T* Foam::Matrix::operator[](const label i) +template +inline Type* Foam::Matrix::operator[](const label i) { # ifdef FULLDEBUG checki(i); @@ -119,9 +114,8 @@ inline T* Foam::Matrix::operator[](const label i) } -// Return subscript-checked const element access -template -inline const T* Foam::Matrix::operator[](const label i) const +template +inline const Type* Foam::Matrix::operator[](const label i) const { # ifdef FULLDEBUG checki(i); diff --git a/src/OpenFOAM/containers/Matrix/MatrixIO.C b/src/OpenFOAM/matrices/Matrix/MatrixIO.C similarity index 83% rename from src/OpenFOAM/containers/Matrix/MatrixIO.C rename to src/OpenFOAM/matrices/Matrix/MatrixIO.C index 710e42876e..15f0d5345e 100644 --- a/src/OpenFOAM/containers/Matrix/MatrixIO.C +++ b/src/OpenFOAM/matrices/Matrix/MatrixIO.C @@ -32,8 +32,8 @@ License // * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * // -template -Foam::Matrix::Matrix(Istream& is) +template +Foam::Matrix::Matrix(Istream& is) : n_(0), m_(0), @@ -43,17 +43,17 @@ Foam::Matrix::Matrix(Istream& is) } -template -Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) +template +Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) { // Anull matrix M.clear(); - is.fatalCheck("operator>>(Istream&, Matrix&)"); + is.fatalCheck("operator>>(Istream&, Matrix&)"); token firstToken(is); - is.fatalCheck("operator>>(Istream&, Matrix&) : reading first token"); + is.fatalCheck("operator>>(Istream&, Matrix&) : reading first token"); if (firstToken.isLabel()) { @@ -63,7 +63,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) label nm = M.n_*M.m_; // Read list contents depending on data format - if (is.format() == IOstream::ASCII || !contiguous()) + if (is.format() == IOstream::ASCII || !contiguous()) { // Read beginning of contents char listDelimiter = is.readBeginList("Matrix"); @@ -71,7 +71,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) if (nm) { M.allocate(); - T* v = M.v_[0]; + Type* v = M.v_[0]; if (listDelimiter == token::BEGIN_LIST) { @@ -88,7 +88,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) is.fatalCheck ( - "operator>>(Istream&, Matrix&) : " + "operator>>(Istream&, Matrix&) : " "reading entry" ); } @@ -98,12 +98,12 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) } else { - T element; + Type element; is >> element; is.fatalCheck ( - "operator>>(Istream&, Matrix&) : " + "operator>>(Istream&, Matrix&) : " "reading the single entry" ); @@ -122,13 +122,13 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) if (nm) { M.allocate(); - T* v = M.v_[0]; + Type* v = M.v_[0]; - is.read(reinterpret_cast(v), nm*sizeof(T)); + is.read(reinterpret_cast(v), nm*sizeof(Type)); is.fatalCheck ( - "operator>>(Istream&, Matrix&) : " + "operator>>(Istream&, Matrix&) : " "reading the binary block" ); } @@ -136,7 +136,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) } else { - FatalIOErrorIn("operator>>(Istream&, Matrix&)", is) + FatalIOErrorIn("operator>>(Istream&, Matrix&)", is) << "incorrect first token, expected , found " << firstToken.info() << exit(FatalIOError); @@ -146,23 +146,23 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) } -template -Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix& M) +template +Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix& M) { label nm = M.n_*M.m_; os << M.n() << token::SPACE << M.m(); // Write list contents depending on data format - if (os.format() == IOstream::ASCII || !contiguous()) + if (os.format() == IOstream::ASCII || !contiguous()) { if (nm) { bool uniform = false; - const T* v = M.v_[0]; + const Type* v = M.v_[0]; - if (nm > 1 && contiguous()) + if (nm > 1 && contiguous()) { uniform = true; @@ -187,7 +187,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix& M) // Write end of contents delimiter os << token::END_BLOCK; } - else if (nm < 10 && contiguous()) + else if (nm < 10 && contiguous()) { // Write size of list and start contents delimiter os << token::BEGIN_LIST; @@ -246,7 +246,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix& M) { if (nm) { - os.write(reinterpret_cast(M.v_[0]), nm*sizeof(T)); + os.write(reinterpret_cast(M.v_[0]), nm*sizeof(Type)); } } diff --git a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H new file mode 100644 index 0000000000..e5e6b16777 --- /dev/null +++ b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H @@ -0,0 +1,92 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM 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 2 of the License, or (at your + option) any later version. + + OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::RectangularMatrix + +Description + A templated 2D rectangular matrix of objects of \, where the n x n matrix + dimension is known and used for subscript bounds checking, etc. + +SourceFiles + RectangularMatrixI.H + RectangularMatrix.C + +\*---------------------------------------------------------------------------*/ + +#ifndef RectangularMatrix_H +#define RectangularMatrix_H + +#include "Matrix.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class Matrix Declaration +\*---------------------------------------------------------------------------*/ + +template +class RectangularMatrix +: + public Matrix, Type> +{ + +public: + + // Constructors + + //- Null constructor. + inline RectangularMatrix(); + + //- Construct given number of rows and columns, + inline RectangularMatrix(const label m, const label n); + + //- Construct with given number of rows and columns + // and value for all elements. + inline RectangularMatrix(const label m, const label n, const Type&); + + //- Construct from Istream. + inline RectangularMatrix(Istream&); + + //- Clone + inline autoPtr > clone() const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +# include "RectangularMatrixI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H new file mode 100644 index 0000000000..1cacd2f11a --- /dev/null +++ b/src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H @@ -0,0 +1,70 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM 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 2 of the License, or (at your + option) any later version. + + OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +inline Foam::RectangularMatrix::RectangularMatrix() +: + Matrix, Type>() +{} + +template +inline Foam::RectangularMatrix::RectangularMatrix +( + const label m, + const label n +) +: + Matrix, Type>(m, n) +{} + +template +inline Foam::RectangularMatrix::RectangularMatrix +( + const label m, + const label n, + const Type& t +) +: + Matrix, Type>(m, n, t) +{} + +template +inline Foam::RectangularMatrix::RectangularMatrix(Istream& is) +: + Matrix, Type>(is) +{} + +template +inline Foam::autoPtr > +Foam::RectangularMatrix::clone() const +{ + return autoPtr >(new RectangularMatrix(*this)); +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H similarity index 52% rename from src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.H rename to src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H index f2b0764f1a..27160dc239 100644 --- a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.H +++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H @@ -23,22 +23,22 @@ License Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA Class - Foam::scalarMatrix + Foam::SquareMatrix Description - Foam::scalarMatrix + A templated 2D square matrix of objects of \, where the n x n matrix + dimension is known and used for subscript bounds checking, etc. SourceFiles - scalarMatrix.C + SquareMatrixI.H + SquareMatrix.C \*---------------------------------------------------------------------------*/ -#ifndef scalarMatrix_H -#define scalarMatrix_H +#ifndef SquareMatrix_H +#define SquareMatrix_H #include "Matrix.H" -#include "scalarField.H" -#include "labelList.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -46,65 +46,39 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class scalarMatrix Declaration + Class Matrix Declaration \*---------------------------------------------------------------------------*/ -class scalarMatrix +template +class SquareMatrix : - public Matrix + public Matrix, Type> { public: // Constructors - //- Construct null - scalarMatrix(); + //- Null constructor. + inline SquareMatrix(); - //- Construct given size - scalarMatrix(const label); + //- Construct given number of rows/columns. + inline SquareMatrix(const label n); - //- Construct from Matrix - scalarMatrix(const Matrix&); + //- Construct given number of rows and columns, + // It checks that m == n. + inline SquareMatrix(const label m, const label n); - //- Construct from Istream - scalarMatrix(Istream&); + //- Construct with given number of rows and rows + // and value for all elements. + // It checks that m == n. + inline SquareMatrix(const label m, const label n, const Type&); + //- Construct from Istream. + inline SquareMatrix(Istream&); - // Member Functions - - //- Solve the matrix using Gaussian elimination with pivoting, - // returning the solution in the source - template - static void solve(Matrix& matrix, Field& source); - - //- Solve the matrix using Gaussian elimination with pivoting - // and return the solution - template - void solve(Field& psi, const Field& source) const; - - - //- LU decompose the matrix with pivoting - static void LUDecompose - ( - Matrix& matrix, - labelList& pivotIndices - ); - - //- LU back-substitution with given source, returning the solution - // in the source - template - static void LUBacksubstitute - ( - const Matrix& luMmatrix, - const labelList& pivotIndices, - Field& source - ); - - //- Solve the matrix using LU decomposition with pivoting - // returning the LU form of the matrix and the solution in the source - template - static void LUsolve(Matrix& matrix, Field& source); + //- Clone + inline autoPtr > clone() const; }; @@ -114,9 +88,7 @@ public: // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -#ifdef NoRepository -# include "scalarMatrixTemplates.C" -#endif +# include "SquareMatrixI.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H new file mode 100644 index 0000000000..d3ee1cd6bf --- /dev/null +++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H @@ -0,0 +1,89 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM 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 2 of the License, or (at your + option) any later version. + + OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +inline Foam::SquareMatrix::SquareMatrix() +: + Matrix, Type>() +{} + +template +inline Foam::SquareMatrix::SquareMatrix(const label n) +: + Matrix, Type>(n, n) +{} + +template +inline Foam::SquareMatrix::SquareMatrix(const label m, const label n) +: + Matrix, Type>(m, n) +{ + if (m != n) + { + FatalErrorIn + ( + "SquareMatrix::SquareMatrix(const label m, const label n)" + ) << "m != n for constructing a square matrix" << exit(FatalError); + } +} + +template +inline Foam::SquareMatrix::SquareMatrix +( + const label m, + const label n, + const Type& t +) +: + Matrix, Type>(m, n, t) +{ + if (m != n) + { + FatalErrorIn + ( + "SquareMatrix::SquareMatrix" + "(const label m, const label n, const Type&)" + ) << "m != n for constructing a square matrix" << exit(FatalError); + } +} + +template +inline Foam::SquareMatrix::SquareMatrix(Istream& is) +: + Matrix, Type>(is) +{} + +template +inline Foam::autoPtr > +Foam::SquareMatrix::clone() const +{ + return autoPtr >(new SquareMatrix(*this)); +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C new file mode 100644 index 0000000000..13ef29ed9b --- /dev/null +++ b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C @@ -0,0 +1,403 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2005 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM 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 2 of the License, or (at your + option) any later version. + + OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "SVD.H" +#include "scalarList.H" +#include "scalarMatrices.H" +#include "ListOps.H" + +// * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * // + +Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition) +: + U_(A), + V_(A.m(), A.m()), + S_(A.m()), + VSinvUt_(A.m(), A.n()), + nZeros_(0) +{ + // SVDcomp to find U_, V_ and S_ - the singular values + + const label Um = U_.m(); + const label Un = U_.n(); + + scalarList rv1(Um); + scalar g = 0; + scalar scale = 0; + scalar s = 0; + scalar anorm = 0; + label l = 0; + + for (label i = 0; i < Um; i++) + { + l = i+2; + rv1[i] = scale*g; + g = s = scale = 0; + + if (i < Un) + { + for (label k = i; k < Un; k++) + { + scale += mag(U_[k][i]); + } + + if (scale != 0) + { + for (label k = i; k < Un; k++) + { + U_[k][i] /= scale; + s += U_[k][i]*U_[k][i]; + } + + scalar f = U_[i][i]; + g = -sign(Foam::sqrt(s), f); + scalar h = f*g - s; + U_[i][i] = f - g; + + for (label j = l-1; j < Um; j++) + { + s = 0; + for (label k = i; k < Un; k++) + { + s += U_[k][i]*U_[k][j]; + } + + f = s/h; + for (label k = i; k < A.n(); k++) + { + U_[k][j] += f*U_[k][i]; + } + } + + for (label k = i; k < Un; k++) + { + U_[k][i] *= scale; + } + } + } + + S_[i] = scale*g; + + g = s = scale = 0; + + if (i+1 <= Un && i != Um) + { + for (label k = l-1; k < Um; k++) + { + scale += mag(U_[i][k]); + } + + if (scale != 0) + { + for (label k=l-1; k < Um; k++) + { + U_[i][k] /= scale; + s += U_[i][k]*U_[i][k]; + } + + scalar f = U_[i][l-1]; + g = -sign(Foam::sqrt(s),f); + scalar h = f*g - s; + U_[i][l-1] = f - g; + + for (label k = l-1; k < Um; k++) + { + rv1[k] = U_[i][k]/h; + } + + for (label j = l-1; j < Un; j++) + { + s = 0; + for (label k = l-1; k < Um; k++) + { + s += U_[j][k]*U_[i][k]; + } + + for (label k = l-1; k < Um; k++) + { + U_[j][k] += s*rv1[k]; + } + } + for (label k = l-1; k < Um; k++) + { + U_[i][k] *= scale; + } + } + } + + anorm = max(anorm, mag(S_[i]) + mag(rv1[i])); + } + + for (label i = Um-1; i >= 0; i--) + { + if (i < Um-1) + { + if (g != 0) + { + for (label j = l; j < Um; j++) + { + V_[j][i] = (U_[i][j]/U_[i][l])/g; + } + + for (label j=l; j < Um; j++) + { + s = 0; + for (label k = l; k < Um; k++) + { + s += U_[i][k]*V_[k][j]; + } + + for (label k = l; k < Um; k++) + { + V_[k][j] += s*V_[k][i]; + } + } + } + + for (label j = l; j < Um;j++) + { + V_[i][j] = V_[j][i] = 0.0; + } + } + + V_[i][i] = 1; + g = rv1[i]; + l = i; + } + + for (label i = min(Um, Un) - 1; i >= 0; i--) + { + l = i+1; + g = S_[i]; + + for (label j = l; j < Um; j++) + { + U_[i][j] = 0.0; + } + + if (g != 0) + { + g = 1.0/g; + + for (label j = l; j < Um; j++) + { + s = 0; + for (label k = l; k < Un; k++) + { + s += U_[k][i]*U_[k][j]; + } + + scalar f = (s/U_[i][i])*g; + + for (label k = i; k < Un; k++) + { + U_[k][j] += f*U_[k][i]; + } + } + + for (label j = i; j < Un; j++) + { + U_[j][i] *= g; + } + } + else + { + for (label j = i; j < Un; j++) + { + U_[j][i] = 0.0; + } + } + + ++U_[i][i]; + } + + for (label k = Um-1; k >= 0; k--) + { + for (label its = 0; its < 35; its++) + { + bool flag = true; + + label nm; + for (l = k; l >= 0; l--) + { + nm = l-1; + if (mag(rv1[l]) + anorm == anorm) + { + flag = false; + break; + } + if (mag(S_[nm]) + anorm == anorm) break; + } + + if (flag) + { + scalar c = 0.0; + s = 1.0; + for (label i = l-1; i < k+1; i++) + { + scalar f = s*rv1[i]; + rv1[i] = c*rv1[i]; + + if (mag(f) + anorm == anorm) break; + + g = S_[i]; + scalar h = sqrtSumSqr(f, g); + S_[i] = h; + h = 1.0/h; + c = g*h; + s = -f*h; + + for (label j = 0; j < Un; j++) + { + scalar y = U_[j][nm]; + scalar z = U_[j][i]; + U_[j][nm] = y*c + z*s; + U_[j][i] = z*c - y*s; + } + } + } + + scalar z = S_[k]; + + if (l == k) + { + if (z < 0.0) + { + S_[k] = -z; + + for (label j = 0; j < Um; j++) V_[j][k] = -V_[j][k]; + } + break; + } + if (its == 34) + { + WarningIn + ( + "SVD::SVD" + "(scalarRectangularMatrix& A, const scalar minCondition)" + ) << "no convergence in 35 SVD iterations" + << endl; + } + + scalar x = S_[l]; + nm = k-1; + scalar y = S_[nm]; + g = rv1[nm]; + scalar h = rv1[k]; + scalar f = ((y - z)*(y + z) + (g - h)*(g + h))/(2.0*h*y); + g = sqrtSumSqr(f, 1.0); + f = ((x - z)*(x + z) + h*((y/(f + sign(g, f))) - h))/x; + scalar c = 1.0; + s = 1.0; + + for (label j = l; j <= nm; j++) + { + label i = j + 1; + g = rv1[i]; + y = S_[i]; + h = s*g; + g = c*g; + scalar z = sqrtSumSqr(f, h); + rv1[j] = z; + c = f/z; + s = h/z; + f = x*c + g*s; + g = g*c - x*s; + h = y*s; + y *= c; + + for (label jj = 0; jj < Um; jj++) + { + x = V_[jj][j]; + z = V_[jj][i]; + V_[jj][j] = x*c + z*s; + V_[jj][i] = z*c - x*s; + } + + z = sqrtSumSqr(f, h); + S_[j] = z; + if (z) + { + z = 1.0/z; + c = f*z; + s = h*z; + } + f = c*g + s*y; + x = c*y - s*g; + + for (label jj=0; jj < Un; jj++) + { + y = U_[jj][j]; + z = U_[jj][i]; + U_[jj][j] = y*c + z*s; + U_[jj][i] = z*c - y*s; + } + } + rv1[l] = 0.0; + rv1[k] = f; + S_[k] = x; + } + } + + // zero singular values that are less than minCondition*maxS + const scalar minS = minCondition*S_[findMax(S_)]; + for (label i = 0; i < S_.size(); i++) + { + if (S_[i] <= minS) + { + //Info << "Removing " << S_[i] << " < " << minS << endl; + S_[i] = 0; + nZeros_++; + } + } + + // now multiply out to find the pseudo inverse of A, VSinvUt_ + multiply(VSinvUt_, V_, inv(S_), U_.T()); + + // test SVD + /*scalarRectangularMatrix SVDA(A.n(), A.m()); + multiply(SVDA, U_, S_, transpose(V_)); + scalar maxDiff = 0; + scalar diff = 0; + for(label i = 0; i < A.n(); i++) + { + for(label j = 0; j < A.m(); j++) + { + diff = mag(A[i][j] - SVDA[i][j]); + if (diff > maxDiff) maxDiff = diff; + } + } + Info << "Maximum discrepancy between A and svd(A) = " << maxDiff << endl; + + if (maxDiff > 4) + { + Info << "singular values " << S_ << endl; + } + */ +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.H b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.H new file mode 100644 index 0000000000..48e581198a --- /dev/null +++ b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.H @@ -0,0 +1,128 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2005 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM 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 2 of the License, or (at your + option) any later version. + + OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::SVD + +Description + Singular value decomposition of a rectangular matrix. + +SourceFiles + SVDI.H + SVD.C + +\*---------------------------------------------------------------------------*/ + +#ifndef SVD_H +#define SVD_H + +#include "scalarMatrices.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + +/*---------------------------------------------------------------------------*\ + Class SVD Declaration +\*---------------------------------------------------------------------------*/ + +class SVD +{ + // Private data + + //- Rectangular matrix with the same dimensions as the input + scalarRectangularMatrix U_; + + //- square matrix V + scalarRectangularMatrix V_; + + //- The singular values + DiagonalMatrix S_; + + //- The matrix product V S^(-1) U^T + scalarRectangularMatrix VSinvUt_; + + //- The number of zero singular values + label nZeros_; + + + // Private Member Functions + + //- Disallow default bitwise copy construct + SVD(const SVD&); + + //- Disallow default bitwise assignment + void operator=(const SVD&); + + template + inline const T sign(const T& a, const T& b); + + +public: + + // Constructors + + //- Construct from a rectangular Matrix + SVD(const scalarRectangularMatrix& A, const scalar minCondition = 0); + + + // Access functions + + //- Return U + inline const scalarRectangularMatrix& U() const; + + //- Return the square matrix V + inline const scalarRectangularMatrix& V() const; + + //- Return the singular values + inline const scalarDiagonalMatrix& S() const; + + //- Return VSinvUt (the pseudo inverse) + inline const scalarRectangularMatrix& VSinvUt() const; + + //- Return the number of zero singular values + inline label nZeros() const; + + //- Return the minimum non-zero singular value + inline scalar minNonZeroS() const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "SVDI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/scalarMatrices/SVD/SVDI.H b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVDI.H new file mode 100644 index 0000000000..298aac1aeb --- /dev/null +++ b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVDI.H @@ -0,0 +1,75 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM 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 2 of the License, or (at your + option) any later version. + + OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +// * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * // + +template +inline const T Foam::SVD::sign(const T& a, const T& b) +{ + return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +inline const Foam::scalarRectangularMatrix& Foam::SVD::U() const +{ + return U_; +} + +inline const Foam::scalarRectangularMatrix& Foam::SVD::V() const +{ + return V_; +} + +inline const Foam::scalarDiagonalMatrix& Foam::SVD::S() const +{ + return S_; +} + +inline const Foam::scalarRectangularMatrix& Foam::SVD::VSinvUt() const +{ + return VSinvUt_; +} + +inline Foam::label Foam::SVD::nZeros() const +{ + return nZeros_; +} + +inline Foam::scalar Foam::SVD::minNonZeroS() const +{ + scalar minS = S_[0]; + for(label i = 1; i < S_.size(); i++) + { + scalar s = S_[i]; + if (s > VSMALL && s < minS) minS = s; + } + return minS; +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C new file mode 100644 index 0000000000..8ddd01a82c --- /dev/null +++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C @@ -0,0 +1,293 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM 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 2 of the License, or (at your + option) any later version. + + OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "scalarMatrices.H" +#include "SVD.H" + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::LUDecompose +( + scalarSquareMatrix& matrix, + labelList& pivotIndices +) +{ + label n = matrix.n(); + scalar vv[n]; + + for (register label i=0; i largestCoeff) + { + largestCoeff = temp; + } + } + + if (largestCoeff == 0.0) + { + FatalErrorIn + ( + "LUdecompose" + "(scalarSquareMatrix& matrix, labelList& rowIndices)" + ) << "Singular matrix" << exit(FatalError); + } + + vv[i] = 1.0/largestCoeff; + } + + for (register label j=0; j= largestCoeff) + { + largestCoeff = temp; + iMax = i; + } + } + + pivotIndices[j] = iMax; + + if (j != iMax) + { + scalar* __restrict__ matrixiMax = matrix[iMax]; + + for (register label k=0; k& B, " + "const scalarRectangularMatrix& C, " + "scalarRectangularMatrix& answer)" + ) << "A and B must have identical inner dimensions but A.m = " + << A.m() << " and B.n = " << B.size() + << abort(FatalError); + } + + if (B.size() != C.n()) + { + FatalErrorIn + ( + "multiply(" + "const scalarRectangularMatrix& A, " + "const DiagonalMatrix& B, " + "const scalarRectangularMatrix& C, " + "scalarRectangularMatrix& answer)" + ) << "B and C must have identical inner dimensions but B.m = " + << B.size() << " and C.n = " << C.n() + << abort(FatalError); + } + + ans = scalarRectangularMatrix(A.n(), C.m(), scalar(0)); + + for(register label i = 0; i < A.n(); i++) + { + for(register label g = 0; g < C.m(); g++) + { + for(register label l = 0; l < C.n(); l++) + { + ans[i][g] += C[l][g] * A[i][l]*B[l]; + } + } + } +} + + +Foam::RectangularMatrix Foam::SVDinv +( + const scalarRectangularMatrix& A, + scalar minCondition +) +{ + SVD svd(A, minCondition); + return svd.VSinvUt(); +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H new file mode 100644 index 0000000000..17d8727516 --- /dev/null +++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H @@ -0,0 +1,137 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM 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 2 of the License, or (at your + option) any later version. + + OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + scalarMatrices + +Description + Scalar matrices + +SourceFiles + scalarMatrices.C + scalarMatricesTemplates.C + +\*---------------------------------------------------------------------------*/ + +#ifndef scalarMatrices_H +#define scalarMatrices_H + +#include "RectangularMatrix.H" +#include "SquareMatrix.H" +#include "DiagonalMatrix.H" +#include "scalarField.H" +#include "labelList.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +typedef RectangularMatrix scalarRectangularMatrix; +typedef SquareMatrix scalarSquareMatrix; +typedef DiagonalMatrix scalarDiagonalMatrix; + +//- Solve the matrix using Gaussian elimination with pivoting, +// returning the solution in the source +template +void solve(scalarSquareMatrix& matrix, Field& source); + +//- Solve the matrix using Gaussian elimination with pivoting +// and return the solution +template +void solve +( + Field& psi, + const scalarSquareMatrix& matrix, + const Field& source +); + +//- LU decompose the matrix with pivoting +void LUDecompose +( + scalarSquareMatrix& matrix, + labelList& pivotIndices +); + +//- LU back-substitution with given source, returning the solution +// in the source +template +void LUBacksubstitute +( + const scalarSquareMatrix& luMmatrix, + const labelList& pivotIndices, + Field& source +); + +//- Solve the matrix using LU decomposition with pivoting +// returning the LU form of the matrix and the solution in the source +template +void LUsolve(scalarSquareMatrix& matrix, Field& source); + +void multiply +( + scalarRectangularMatrix& answer, // value changed in return + const scalarRectangularMatrix& A, + const scalarRectangularMatrix& B +); + +void multiply +( + scalarRectangularMatrix& answer, // value changed in return + const scalarRectangularMatrix& A, + const scalarRectangularMatrix& B, + const scalarRectangularMatrix& C +); + +void multiply +( + scalarRectangularMatrix& answer, // value changed in return + const scalarRectangularMatrix& A, + const DiagonalMatrix& B, + const scalarRectangularMatrix& C +); + +//- Return the inverse of matrix A using SVD +scalarRectangularMatrix SVDinv +( + const scalarRectangularMatrix& A, + scalar minCondition = 0 +); + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "scalarMatricesTemplates.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrixTemplates.C b/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C similarity index 83% rename from src/OpenFOAM/matrices/scalarMatrix/scalarMatrixTemplates.C rename to src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C index b729e32fcc..131856b3fc 100644 --- a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrixTemplates.C +++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C @@ -24,16 +24,16 @@ License \*---------------------------------------------------------------------------*/ -#include "scalarMatrix.H" +#include "scalarMatrices.H" #include "Swap.H" // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -template -void Foam::scalarMatrix::solve +template +void Foam::solve ( - Matrix& tmpMatrix, - Field& sourceSol + scalarSquareMatrix& tmpMatrix, + Field& sourceSol ) { label n = tmpMatrix.n(); @@ -68,7 +68,7 @@ void Foam::scalarMatrix::solve // Check that the system of equations isn't singular if (mag(tmpMatrix[i][i]) < 1e-20) { - FatalErrorIn("scalarMatrix::solve()") + FatalErrorIn("solve(scalarSquareMatrix&, Field& sourceSol)") << "Singular Matrix" << exit(FatalError); } @@ -89,7 +89,7 @@ void Foam::scalarMatrix::solve // Back-substitution for (register label j=n-1; j>=0; j--) { - T ntempvec = pTraits::zero; + Type ntempvec = pTraits::zero; for (register label k=j+1; k -void Foam::scalarMatrix::solve(Field& psi, const Field& source) const +template +void Foam::solve +( + Field& psi, + const scalarSquareMatrix& matrix, + const Field& source +) { - Matrix tmpMatrix = *this; + scalarSquareMatrix tmpMatrix = matrix; psi = source; solve(tmpMatrix, psi); } -template -void Foam::scalarMatrix::LUBacksubstitute +template +void Foam::LUBacksubstitute ( - const Matrix& luMatrix, + const scalarSquareMatrix& luMatrix, const labelList& pivotIndices, - Field& sourceSol + Field& sourceSol ) { label n = luMatrix.n(); @@ -125,7 +130,7 @@ void Foam::scalarMatrix::LUBacksubstitute for (register label i=0; i::zero) + else if (sum != pTraits::zero) { ii = i+1; } @@ -146,11 +151,11 @@ void Foam::scalarMatrix::LUBacksubstitute for (register label i=n-1; i>=0; i--) { - T sum = sourceSol[i]; + Type sum = sourceSol[i]; const scalar* __restrict__ luMatrixi = luMatrix[i]; for (register label j=i+1; j -void Foam::scalarMatrix::LUsolve +template +void Foam::LUsolve ( - Matrix& matrix, - Field& sourceSol + scalarSquareMatrix& matrix, + Field& sourceSol ) { labelList pivotIndices(matrix.n()); diff --git a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.C b/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.C deleted file mode 100644 index 2ff9e39f18..0000000000 --- a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.C +++ /dev/null @@ -1,161 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM 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 2 of the License, or (at your - option) any later version. - - OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -\*---------------------------------------------------------------------------*/ - -#include "scalarMatrix.H" - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -Foam::scalarMatrix::scalarMatrix() -{} - - -Foam::scalarMatrix::scalarMatrix(const label mSize) -: - Matrix(mSize, mSize, 0.0) -{} - - -Foam::scalarMatrix::scalarMatrix(const Matrix& matrix) -: - Matrix(matrix) -{} - - -Foam::scalarMatrix::scalarMatrix(Istream& is) -: - Matrix(is) -{} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -void Foam::scalarMatrix::LUDecompose -( - Matrix& matrix, - labelList& pivotIndices -) -{ - label n = matrix.n(); - scalar vv[n]; - - for (register label i=0; i largestCoeff) - { - largestCoeff = temp; - } - } - - if (largestCoeff == 0.0) - { - FatalErrorIn - ( - "scalarMatrix::LUdecompose" - "(Matrix& matrix, labelList& rowIndices)" - ) << "Singular matrix" << exit(FatalError); - } - - vv[i] = 1.0/largestCoeff; - } - - for (register label j=0; j= largestCoeff) - { - largestCoeff = temp; - iMax = i; - } - } - - pivotIndices[j] = iMax; - - if (j != iMax) - { - scalar* __restrict__ matrixiMax = matrix[iMax]; - - for (register label k=0; k -Foam::simpleMatrix::simpleMatrix(const label mSize) +template +Foam::simpleMatrix::simpleMatrix(const label mSize) : - scalarMatrix(mSize), - source_(mSize, pTraits::zero) + scalarSquareMatrix(mSize), + source_(mSize, pTraits::zero) {} -template -Foam::simpleMatrix::simpleMatrix +template +Foam::simpleMatrix::simpleMatrix ( - const scalarMatrix& matrix, - const Field& source + const scalarSquareMatrix& matrix, + const Field& source ) : - scalarMatrix(matrix), + scalarSquareMatrix(matrix), source_(source) {} -template -Foam::simpleMatrix::simpleMatrix(Istream& is) +template +Foam::simpleMatrix::simpleMatrix(Istream& is) : - scalarMatrix(is), + scalarSquareMatrix(is), source_(is) {} // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -template -Foam::Field Foam::simpleMatrix::solve() const +template +Foam::Field Foam::simpleMatrix::solve() const { - scalarMatrix tmpMatrix = *this; - Field sourceSol = source_; + scalarSquareMatrix tmpMatrix = *this; + Field sourceSol = source_; - scalarMatrix::solve(tmpMatrix, sourceSol); + Foam::solve(tmpMatrix, sourceSol); return sourceSol; } -template -Foam::Field Foam::simpleMatrix::LUsolve() const +template +Foam::Field Foam::simpleMatrix::LUsolve() const { - scalarMatrix luMatrix = *this; - Field sourceSol = source_; + scalarSquareMatrix luMatrix = *this; + Field sourceSol = source_; - scalarMatrix::LUsolve(luMatrix, sourceSol); + Foam::LUsolve(luMatrix, sourceSol); return sourceSol; } @@ -84,82 +84,82 @@ Foam::Field Foam::simpleMatrix::LUsolve() const // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // -template -void Foam::simpleMatrix::operator=(const simpleMatrix& m) +template +void Foam::simpleMatrix::operator=(const simpleMatrix& m) { if (this == &m) { - FatalErrorIn("simpleMatrix::operator=(const simpleMatrix&)") + FatalErrorIn("simpleMatrix::operator=(const simpleMatrix&)") << "Attempted assignment to self" << abort(FatalError); } if (n() != m.n()) { - FatalErrorIn("simpleMatrix::operator=(const simpleMatrix&)") + FatalErrorIn("simpleMatrix::operator=(const simpleMatrix&)") << "Different size matrices" << abort(FatalError); } if (source_.size() != m.source_.size()) { - FatalErrorIn("simpleMatrix::operator=(const simpleMatrix&)") + FatalErrorIn("simpleMatrix::operator=(const simpleMatrix&)") << "Different size source vectors" << abort(FatalError); } - scalarMatrix::operator=(m); + scalarSquareMatrix::operator=(m); source_ = m.source_; } // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * // -template -Foam::simpleMatrix Foam::operator+ +template +Foam::simpleMatrix Foam::operator+ ( - const simpleMatrix& m1, - const simpleMatrix& m2 + const simpleMatrix& m1, + const simpleMatrix& m2 ) { - return simpleMatrix + return simpleMatrix ( - static_cast(m1) - + static_cast(m2), + static_cast(m1) + + static_cast(m2), m1.source_ + m2.source_ ); } -template -Foam::simpleMatrix Foam::operator- +template +Foam::simpleMatrix Foam::operator- ( - const simpleMatrix& m1, - const simpleMatrix& m2 + const simpleMatrix& m1, + const simpleMatrix& m2 ) { - return simpleMatrix + return simpleMatrix ( - static_cast(m1) - - static_cast(m2), + static_cast(m1) + - static_cast(m2), m1.source_ - m2.source_ ); } -template -Foam::simpleMatrix Foam::operator*(const scalar s, const simpleMatrix& m) +template +Foam::simpleMatrix Foam::operator*(const scalar s, const simpleMatrix& m) { - return simpleMatrix(s*m.matrix_, s*m.source_); + return simpleMatrix(s*m.matrix_, s*m.source_); } // * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * // -template -Foam::Ostream& Foam::operator<<(Ostream& os, const simpleMatrix& m) +template +Foam::Ostream& Foam::operator<<(Ostream& os, const simpleMatrix& m) { - os << static_cast(m) << nl << m.source_; + os << static_cast(m) << nl << m.source_; return os; } diff --git a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H index d1c411452b..a972c787c8 100644 --- a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H +++ b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H @@ -36,7 +36,7 @@ SourceFiles #ifndef simpleMatrix_H #define simpleMatrix_H -#include "scalarMatrix.H" +#include "scalarMatrices.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -45,35 +45,14 @@ namespace Foam // Forward declaration of friend functions and operators -template +template class simpleMatrix; -template -simpleMatrix operator+ -( - const simpleMatrix&, - const simpleMatrix& -); - -template -simpleMatrix operator- -( - const simpleMatrix&, - const simpleMatrix& -); - -template -simpleMatrix operator* -( - const scalar, - const simpleMatrix& -); - -template +template Ostream& operator<< ( Ostream&, - const simpleMatrix& + const simpleMatrix& ); @@ -81,14 +60,14 @@ Ostream& operator<< Class simpleMatrix Declaration \*---------------------------------------------------------------------------*/ -template +template class simpleMatrix : - public scalarMatrix + public scalarSquareMatrix { // Private data - Field source_; + Field source_; public: @@ -99,25 +78,25 @@ public: simpleMatrix(const label); //- Construct from components - simpleMatrix(const scalarMatrix&, const Field&); + simpleMatrix(const scalarSquareMatrix&, const Field&); //- Construct from Istream simpleMatrix(Istream&); //- Construct as copy - simpleMatrix(const simpleMatrix&); + simpleMatrix(const simpleMatrix&); // Member Functions // Access - Field& source() + Field& source() { return source_; } - const Field& source() const + const Field& source() const { return source_; } @@ -125,49 +104,52 @@ public: //- Solve the matrix using Gaussian elimination with pivoting // and return the solution - Field solve() const; + Field solve() const; //- Solve the matrix using LU decomposition with pivoting // and return the solution - Field LUsolve() const; + Field LUsolve() const; // Member Operators - void operator=(const simpleMatrix&); - - - // Friend Operators - - friend simpleMatrix operator+ - ( - const simpleMatrix&, - const simpleMatrix& - ); - - friend simpleMatrix operator- - ( - const simpleMatrix&, - const simpleMatrix& - ); - - friend simpleMatrix operator* - ( - const scalar, - const simpleMatrix& - ); + void operator=(const simpleMatrix&); // Ostream Operator - friend Ostream& operator<< + friend Ostream& operator<< ( Ostream&, - const simpleMatrix& + const simpleMatrix& ); }; +// Global operators + +template +simpleMatrix operator+ +( + const simpleMatrix&, + const simpleMatrix& +); + +template +simpleMatrix operator- +( + const simpleMatrix&, + const simpleMatrix& +); + +template +simpleMatrix operator* +( + const scalar, + const simpleMatrix& +); + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files index 0c3fae1887..25137c58c7 100644 --- a/src/finiteVolume/Make/files +++ b/src/finiteVolume/Make/files @@ -38,6 +38,9 @@ fvMeshMapper = fvMesh/fvMeshMapper $(fvMeshMapper)/fvPatchMapper.C $(fvMeshMapper)/fvSurfaceMapper.C +extendedStencil = fvMesh/extendedStencil +$(extendedStencil)/extendedStencil.C + fvPatchFields = fields/fvPatchFields $(fvPatchFields)/fvPatchField/fvPatchFields.C @@ -165,6 +168,8 @@ $(schemes)/harmonic/harmonic.C $(schemes)/localBlended/localBlended.C $(schemes)/localMax/localMax.C $(schemes)/localMin/localMin.C +$(schemes)/quadraticFit/quadraticFit.C +$(schemes)/quadraticFit/quadraticFitData.C limitedSchemes = $(surfaceInterpolation)/limitedSchemes $(limitedSchemes)/limitedSurfaceInterpolationScheme/limitedSurfaceInterpolationSchemes.C diff --git a/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.C b/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.C new file mode 100644 index 0000000000..bb4c6c1534 --- /dev/null +++ b/src/finiteVolume/fvMesh/extendedStencil/extendedStencil.C @@ -0,0 +1,525 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM 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 2 of the License, or (at your + option) any later version. + + OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "extendedStencil.H" +#include "globalIndex.H" +#include "syncTools.H" +#include "SortableList.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +// Calculates per face a list of global cell/face indices. +void Foam::extendedStencil::calcFaceStencils +( + const polyMesh& mesh, + const globalIndex& globalNumbering +) +{ + const polyBoundaryMesh& patches = mesh.boundaryMesh(); + const label nBnd = mesh.nFaces()-mesh.nInternalFaces(); + const labelList& own = mesh.faceOwner(); + const labelList& nei = mesh.faceNeighbour(); + + + // Determine neighbouring global cell or boundary face + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + labelList neiGlobal(nBnd); + forAll(patches, patchI) + { + const polyPatch& pp = patches[patchI]; + label faceI = pp.start(); + + if (pp.coupled()) + { + // For coupled faces get the cell on the other side + forAll(pp, i) + { + label bFaceI = faceI-mesh.nInternalFaces(); + neiGlobal[bFaceI] = globalNumbering.toGlobal(own[faceI]); + faceI++; + } + } + else if (isA(pp)) + { + forAll(pp, i) + { + label bFaceI = faceI-mesh.nInternalFaces(); + neiGlobal[bFaceI] = -1; + faceI++; + } + } + else + { + // For noncoupled faces get the boundary face. + forAll(pp, i) + { + label bFaceI = faceI-mesh.nInternalFaces(); + neiGlobal[bFaceI] = + globalNumbering.toGlobal(mesh.nCells()+bFaceI); + faceI++; + } + } + } + syncTools::swapBoundaryFaceList(mesh, neiGlobal, false); + + + // Determine cellCells in global numbering + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + labelListList globalCellCells(mesh.nCells()); + forAll(globalCellCells, cellI) + { + const cell& cFaces = mesh.cells()[cellI]; + + labelList& cCells = globalCellCells[cellI]; + + cCells.setSize(cFaces.size()); + + // Collect neighbouring cells/faces + label nNbr = 0; + forAll(cFaces, i) + { + label faceI = cFaces[i]; + + if (mesh.isInternalFace(faceI)) + { + label nbrCellI = own[faceI]; + if (nbrCellI == cellI) + { + nbrCellI = nei[faceI]; + } + cCells[nNbr++] = globalNumbering.toGlobal(nbrCellI); + } + else + { + label nbrCellI = neiGlobal[faceI-mesh.nInternalFaces()]; + if (nbrCellI != -1) + { + cCells[nNbr++] = nbrCellI; + } + } + } + cCells.setSize(nNbr); + } + + + // Determine neighbouring global cell Cells + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + labelListList neiGlobalCellCells(nBnd); + for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++) + { + neiGlobalCellCells[faceI-mesh.nInternalFaces()] = + globalCellCells[own[faceI]]; + } + syncTools::swapBoundaryFaceList(mesh, neiGlobalCellCells, false); + + + + // Construct stencil in global numbering + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + stencil_.setSize(mesh.nFaces()); + + labelHashSet faceStencil; + + for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++) + { + faceStencil.clear(); + label globalOwn = globalNumbering.toGlobal(own[faceI]); + faceStencil.insert(globalOwn); + const labelList& ownCCells = globalCellCells[own[faceI]]; + forAll(ownCCells, i) + { + faceStencil.insert(ownCCells[i]); + } + + label globalNei = globalNumbering.toGlobal(nei[faceI]); + faceStencil.insert(globalNei); + const labelList& neiCCells = globalCellCells[nei[faceI]]; + forAll(neiCCells, i) + { + faceStencil.insert(neiCCells[i]); + } + + // Guarantee owner first, neighbour second. + stencil_[faceI].setSize(faceStencil.size()); + label n = 0; + stencil_[faceI][n++] = globalOwn; + stencil_[faceI][n++] = globalNei; + forAllConstIter(labelHashSet, faceStencil, iter) + { + if (iter.key() != globalOwn && iter.key() != globalNei) + { + stencil_[faceI][n++] = iter.key(); + } + } + //Pout<< "internalface:" << faceI << " toc:" << faceStencil.toc() + // << " stencil:" << stencil_[faceI] << endl; + } + forAll(patches, patchI) + { + const polyPatch& pp = patches[patchI]; + label faceI = pp.start(); + + if (pp.coupled()) + { + forAll(pp, i) + { + faceStencil.clear(); + label globalOwn = globalNumbering.toGlobal(own[faceI]); + faceStencil.insert(globalOwn); + const labelList& ownCCells = globalCellCells[own[faceI]]; + forAll(ownCCells, i) + { + faceStencil.insert(ownCCells[i]); + } + // Get the coupled cell + label globalNei = neiGlobal[faceI-mesh.nInternalFaces()]; + faceStencil.insert(globalNei); + // And the neighbours of the coupled cell + const labelList& neiCCells = + neiGlobalCellCells[faceI-mesh.nInternalFaces()]; + forAll(neiCCells, i) + { + faceStencil.insert(neiCCells[i]); + } + + // Guarantee owner first, neighbour second. + stencil_[faceI].setSize(faceStencil.size()); + label n = 0; + stencil_[faceI][n++] = globalOwn; + stencil_[faceI][n++] = globalNei; + forAllConstIter(labelHashSet, faceStencil, iter) + { + if (iter.key() != globalOwn && iter.key() != globalNei) + { + stencil_[faceI][n++] = iter.key(); + } + } + + //Pout<< "coupledface:" << faceI + // << " toc:" << faceStencil.toc() + // << " stencil:" << stencil_[faceI] << endl; + + faceI++; + } + } + else if (!isA(pp)) + { + forAll(pp, i) + { + faceStencil.clear(); + label globalOwn = globalNumbering.toGlobal(own[faceI]); + faceStencil.insert(globalOwn); + const labelList& ownCCells = globalCellCells[own[faceI]]; + forAll(ownCCells, i) + { + faceStencil.insert(ownCCells[i]); + } + + + // Guarantee owner first, neighbour second. + stencil_[faceI].setSize(faceStencil.size()); + label n = 0; + stencil_[faceI][n++] = globalOwn; + forAllConstIter(labelHashSet, faceStencil, iter) + { + if (iter.key() != globalOwn) + { + stencil_[faceI][n++] = iter.key(); + } + } + + //Pout<< "boundaryface:" << faceI + // << " toc:" << faceStencil.toc() + // << " stencil:" << stencil_[faceI] << endl; + + faceI++; + } + } + } +} + + +// Calculates extended stencil. This is per face +// - owner +// - cellCells of owner +// - neighbour +// - cellCells of neighbour +// It comes in two parts: +// - a map which collects/distributes all necessary data in a compact array +// - the stencil (a labelList per face) which is a set of indices into this +// compact array. +// The compact array is laid out as follows: +// - first data for current processor (Pstream::myProcNo()) +// - all cells +// - all boundary faces +// - then per processor +// - all used cells and boundary faces +void Foam::extendedStencil::calcExtendedFaceStencil(const polyMesh& mesh) +{ + const label nBnd = mesh.nFaces()-mesh.nInternalFaces(); + + // Global numbering for cells and boundary faces + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + globalIndex globalNumbering(mesh.nCells()+nBnd); + + + // Calculate stencil in global cell indices + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + calcFaceStencils(mesh, globalNumbering); + + + // Convert stencil to schedule + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + // We now know what information we need from other processors. This needs + // to be converted into what information I need to send as well + // (mapDistribute) + + + // 1. Construct per processor compact addressing of the global cells + // needed. The ones from the local processor are not included since + // these are always all needed. + List > globalToProc(Pstream::nProcs()); + { + const labelList& procPatchMap = mesh.globalData().procPatchMap(); + const polyBoundaryMesh& patches = mesh.boundaryMesh(); + + // Presize with (as estimate) size of patch to neighbour. + forAll(procPatchMap, procI) + { + if (procPatchMap[procI] != -1) + { + globalToProc[procI].resize + ( + patches[procPatchMap[procI]].size() + ); + } + } + + // Collect all (non-local) globalcells/faces needed. + forAll(stencil_, faceI) + { + const labelList& stencilCells = stencil_[faceI]; + + forAll(stencilCells, i) + { + label globalCellI = stencilCells[i]; + label procI = globalNumbering.whichProcID(stencilCells[i]); + + if (procI != Pstream::myProcNo()) + { + label nCompact = globalToProc[procI].size(); + globalToProc[procI].insert(globalCellI, nCompact); + } + } + } + // Sort global cells needed (not really necessary) + forAll(globalToProc, procI) + { + if (procI != Pstream::myProcNo()) + { + Map