From 318b71e206e9e635aa649938e9d1c5dfbc0b7f6b Mon Sep 17 00:00:00 2001 From: henry Date: Sun, 28 Sep 2008 22:50:57 +0100 Subject: [PATCH] Made Matrix generic base-class for SquareMatrix and RectangularMatrix by introducing the type it will become as a template argument. Brought everything into line with this change. --- src/ODE/ODE/ODE.H | 4 +- src/ODE/ODESolvers/KRR4/KRR4.C | 19 +-- src/ODE/ODESolvers/KRR4/KRR4.H | 6 +- src/ODE/ODESolvers/SIBS/SIBS.H | 10 +- src/ODE/ODESolvers/SIBS/SIMPR.C | 13 +- src/ODE/ODESolvers/SIBS/polyExtrapolate.C | 2 +- src/OpenFOAM/Make/files | 6 +- .../IOstreams/Pstreams/PstreamCommsStruct.C | 2 - .../matrices/DiagonalMatrix/DiagonalMatrix.C | 3 +- .../matrices/DiagonalMatrix/DiagonalMatrix.H | 5 +- .../matrices/LUscalarMatrix/LUscalarMatrix.C | 8 +- .../matrices/LUscalarMatrix/LUscalarMatrix.H | 10 +- src/OpenFOAM/matrices/Matrix/Matrix.C | 106 ++++++------ src/OpenFOAM/matrices/Matrix/Matrix.H | 104 ++++++------ src/OpenFOAM/matrices/Matrix/MatrixI.H | 46 +++--- src/OpenFOAM/matrices/Matrix/MatrixIO.C | 24 +-- .../RectangularMatrix/RectangularMatrix.H | 92 +++++++++++ .../RectangularMatrix/RectangularMatrixI.H | 70 ++++++++ .../matrices/SquareMatrix/SquareMatrix.H | 96 +++++++++++ .../matrices/SquareMatrix/SquareMatrixI.H | 75 +++++++++ .../SVD/SVD.C | 13 +- .../SVD/SVD.H | 22 +-- .../SVD/SVDI.H | 8 +- .../scalarMatrices.C} | 99 +++++------ .../matrices/scalarMatrices/scalarMatrices.H | 137 ++++++++++++++++ .../scalarMatricesTemplates.C} | 25 +-- .../matrices/scalarMatrix/scalarMatrix.H | 155 ------------------ .../matrices/simpleMatrix/simpleMatrix.C | 28 ++-- .../matrices/simpleMatrix/simpleMatrix.H | 72 +++----- .../schemes/quadraticFit/quadraticFitData.C | 14 +- .../chemistryModel/chemistryModel.C | 38 ++--- .../chemistryModel/chemistryModel.H | 3 +- 32 files changed, 791 insertions(+), 524 deletions(-) create mode 100644 src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H create mode 100644 src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H create mode 100644 src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H create mode 100644 src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H rename src/OpenFOAM/matrices/{scalarMatrix => scalarMatrices}/SVD/SVD.C (96%) rename src/OpenFOAM/matrices/{scalarMatrix => scalarMatrices}/SVD/SVD.H (86%) rename src/OpenFOAM/matrices/{scalarMatrix => scalarMatrices}/SVD/SVDI.H (88%) rename src/OpenFOAM/matrices/{scalarMatrix/scalarMatrix.C => scalarMatrices/scalarMatrices.C} (76%) create mode 100644 src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H rename src/OpenFOAM/matrices/{scalarMatrix/scalarMatrixTemplates.C => scalarMatrices/scalarMatricesTemplates.C} (90%) delete mode 100644 src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.H 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 57dc7f8dc7..cf4e331e25 100644 --- a/src/OpenFOAM/Make/files +++ b/src/OpenFOAM/Make/files @@ -177,9 +177,9 @@ dimensionedTypes/dimensionedTensor/dimensionedTensor.C matrices/solution/solution.C -scalarMatrix = matrices/scalarMatrix -$(scalarMatrix)/scalarMatrix.C -$(scalarMatrix)/SVD/SVD.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 index c023e59c5a..a31058eace 100644 --- a/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C +++ b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C @@ -29,7 +29,8 @@ License // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // template -Foam::DiagonalMatrix::DiagonalMatrix(const Matrix& a) +template +Foam::DiagonalMatrix::DiagonalMatrix(const Matrix& a) : List(min(a.n(), a.m())) { diff --git a/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H index 36138fc25d..ae8ceffc1c 100644 --- a/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H +++ b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.H @@ -46,7 +46,7 @@ namespace Foam // * * * * * * * * * * * * Class Forward declaration * * * * * * * * * * * // -template class Matrix; +template class Matrix; /*---------------------------------------------------------------------------*\ Class DiagonalMatrix Declaration @@ -62,7 +62,8 @@ public: // Constructors //- Construct from diagonal component of a Matrix - DiagonalMatrix(const Matrix&); + template + DiagonalMatrix(const Matrix&); //- Construct empty from size DiagonalMatrix(const label size); diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C index 5de8693f2f..f425accf72 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, 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, 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/Matrix/Matrix.C b/src/OpenFOAM/matrices/Matrix/Matrix.C index 4fbef71c2b..577cbcabcb 100644 --- a/src/OpenFOAM/matrices/Matrix/Matrix.C +++ b/src/OpenFOAM/matrices/Matrix/Matrix.C @@ -28,8 +28,8 @@ License // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // -template -void Foam::Matrix::allocate() +template +void Foam::Matrix::allocate() { if (n_ && m_) { @@ -46,8 +46,8 @@ void Foam::Matrix::allocate() // * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * * // -template -Foam::Matrix::~Matrix() +template +Foam::Matrix::~Matrix() { if (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 Type& 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 Type& 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); } @@ -117,8 +117,8 @@ Foam::Matrix::Matrix(const label n, const label m, const Type& a) } -template -Foam::Matrix::Matrix(const Matrix& a) +template +Foam::Matrix::Matrix(const Matrix& a) : n_(a.n_), m_(a.m_), @@ -139,8 +139,8 @@ Foam::Matrix::Matrix(const Matrix& a) } -template -void Foam::Matrix::clear() +template +void Foam::Matrix::clear() { if (v_) { @@ -153,8 +153,8 @@ void Foam::Matrix::clear() } -template -void Foam::Matrix::transfer(Matrix& a) +template +void Foam::Matrix::transfer(Matrix& a) { clear(); @@ -169,13 +169,11 @@ void Foam::Matrix::transfer(Matrix& a) } -// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // - -template -Foam::Matrix Foam::Matrix::T() const +template +Form Foam::Matrix::T() const { - const Matrix& A = *this; - Matrix At(m(), n()); + const Matrix& A = *this; + Form At(m(), n()); for (register label i=0; i Foam::Matrix::T() const } -template -void Foam::Matrix::operator=(const Type& t) +// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // + +template +void Foam::Matrix::operator=(const Type& t) { if (v_) { @@ -206,12 +206,12 @@ void Foam::Matrix::operator=(const Type& 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); } @@ -240,8 +240,8 @@ void Foam::Matrix::operator=(const Matrix& a) // * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * // -template -const Type& Foam::max(const Matrix& a) +template +const Type& Foam::max(const Matrix& a) { label nm = a.n_*a.m_; @@ -262,7 +262,7 @@ const Type& Foam::max(const Matrix& a) } else { - FatalErrorIn("max(const Matrix&)") + FatalErrorIn("max(const Matrix&)") << "matrix is empty" << abort(FatalError); @@ -272,8 +272,8 @@ const Type& Foam::max(const Matrix& a) } -template -const Type& Foam::min(const Matrix& a) +template +const Type& Foam::min(const Matrix& a) { label nm = a.n_*a.m_; @@ -294,7 +294,7 @@ const Type& Foam::min(const Matrix& a) } else { - FatalErrorIn("min(const Matrix&)") + FatalErrorIn("min(const Matrix&)") << "matrix is empty" << abort(FatalError); @@ -306,10 +306,10 @@ const Type& 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_) { @@ -327,14 +327,14 @@ Foam::Matrix 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&)" + "Matrix::operator+(const Matrix&, const Matrix&)" ) << "attempted add matrices with different number of rows: " << a.n_ << ", " << b.n_ << abort(FatalError); @@ -344,13 +344,13 @@ Foam::Matrix Foam::operator+(const Matrix& a, const Matrix& b) { FatalErrorIn ( - "Matrix::operator+(const Matrix&, const Matrix&)" + "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_); Type* abv = ab.v_[0]; const Type* av = a.v_[0]; @@ -366,14 +366,14 @@ Foam::Matrix 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&)" + "Matrix::operator-(const Matrix&, const Matrix&)" ) << "attempted add matrices with different number of rows: " << a.n_ << ", " << b.n_ << abort(FatalError); @@ -383,13 +383,13 @@ Foam::Matrix Foam::operator-(const Matrix& a, const Matrix& b) { FatalErrorIn ( - "Matrix::operator-(const Matrix&, const Matrix&)" + "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_); Type* abv = ab.v_[0]; const Type* av = a.v_[0]; @@ -405,17 +405,17 @@ Foam::Matrix 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_) { 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 Type& max(const Matrix&); -template const Type& min(const Matrix&); - -template Matrix operator-(const Matrix&); -template Matrix operator+ +template Istream& operator>> ( - const Matrix&, - const Matrix& -); -template Matrix operator- -( - const Matrix&, - const Matrix& -); -template Matrix operator* -( - const scalar, - const Matrix& + Istream&, + 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 @@ -112,13 +101,13 @@ public: 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 @@ -129,7 +118,7 @@ public: // Member functions //- Return a null Matrix - static const Matrix& null(); + static const Matrix& null(); // Access @@ -160,7 +149,11 @@ 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 @@ -171,52 +164,49 @@ public: //- Return subscript-checked element of constant Matrix. inline const Type* operator[](const label) const; - //- Return the transpose of the matrix - Matrix T() 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 Type&); - // Friend functions - - friend const Type& max(const Matrix&); - friend const Type& 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& - ); - - // 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/matrices/Matrix/MatrixI.H b/src/OpenFOAM/matrices/Matrix/MatrixI.H index 8323174d4f..8d758cf5b1 100644 --- a/src/OpenFOAM/matrices/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,67 +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_; } -template -inline Foam::label Foam::Matrix::m() const +template +inline Foam::label Foam::Matrix::m() const { return m_; } -template -inline Foam::label Foam::Matrix::size() const +template +inline Foam::label Foam::Matrix::size() const { return n_*m_; } -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); } } -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); } @@ -104,8 +104,8 @@ inline void Foam::Matrix::checkj(const label j) const // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // -template -inline Type* Foam::Matrix::operator[](const label i) +template +inline Type* Foam::Matrix::operator[](const label i) { # ifdef FULLDEBUG checki(i); @@ -114,8 +114,8 @@ inline Type* Foam::Matrix::operator[](const label i) } -template -inline const Type* 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/matrices/Matrix/MatrixIO.C b/src/OpenFOAM/matrices/Matrix/MatrixIO.C index 9b5a922f0d..15f0d5345e 100644 --- a/src/OpenFOAM/matrices/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()) { @@ -88,7 +88,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) is.fatalCheck ( - "operator>>(Istream&, Matrix&) : " + "operator>>(Istream&, Matrix&) : " "reading entry" ); } @@ -103,7 +103,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) is.fatalCheck ( - "operator>>(Istream&, Matrix&) : " + "operator>>(Istream&, Matrix&) : " "reading the single entry" ); @@ -128,7 +128,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix& M) 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,8 +146,8 @@ 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_; 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/SquareMatrix/SquareMatrix.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H new file mode 100644 index 0000000000..6239781d85 --- /dev/null +++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H @@ -0,0 +1,96 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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::SquareMatrix + +Description + 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 + SquareMatrixI.H + SquareMatrix.C + +\*---------------------------------------------------------------------------*/ + +#ifndef SquareMatrix_H +#define SquareMatrix_H + +#include "Matrix.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class Matrix Declaration +\*---------------------------------------------------------------------------*/ + +template +class SquareMatrix +: + public Matrix, Type> +{ + +public: + + // Constructors + + //- Null constructor. + inline SquareMatrix(); + + //- Construct given number of rows/columns. + inline SquareMatrix(const label n); + + //- Construct given number of rows and columns, + // It checks that m == n. + inline SquareMatrix(const label m, const label n); + + //- Construct with given number of rows/columns + // and value for all elements. + inline SquareMatrix(const label n, const Type&); + + //- Construct from Istream. + inline SquareMatrix(Istream&); + + //- Clone + inline autoPtr > clone() const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +# include "SquareMatrixI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H new file mode 100644 index 0000000000..a7371a6a2e --- /dev/null +++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H @@ -0,0 +1,75 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 n, const Type& t) +: + Matrix, Type>(n, t) +{} + +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/scalarMatrix/SVD/SVD.C b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C similarity index 96% rename from src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.C rename to src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C index a1355c393d..13ef29ed9b 100644 --- a/src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.C +++ b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C @@ -26,12 +26,12 @@ License #include "SVD.H" #include "scalarList.H" -#include "scalarMatrix.H" +#include "scalarMatrices.H" #include "ListOps.H" // * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * // -Foam::SVD::SVD(const Matrix& A, const scalar minCondition) +Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition) : U_(A), V_(A.m(), A.m()), @@ -294,8 +294,11 @@ Foam::SVD::SVD(const Matrix& A, const scalar minCondition) } if (its == 34) { - WarningIn("SVD::SVD(Matrix& A)") - << "no convergence in 35 SVD iterations" + WarningIn + ( + "SVD::SVD" + "(scalarRectangularMatrix& A, const scalar minCondition)" + ) << "no convergence in 35 SVD iterations" << endl; } @@ -375,7 +378,7 @@ Foam::SVD::SVD(const Matrix& A, const scalar minCondition) multiply(VSinvUt_, V_, inv(S_), U_.T()); // test SVD - /*Matrix SVDA(A.n(), A.m()); + /*scalarRectangularMatrix SVDA(A.n(), A.m()); multiply(SVDA, U_, S_, transpose(V_)); scalar maxDiff = 0; scalar diff = 0; diff --git a/src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.H b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.H similarity index 86% rename from src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.H rename to src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.H index b597027732..48e581198a 100644 --- a/src/OpenFOAM/matrices/scalarMatrix/SVD/SVD.H +++ b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.H @@ -23,12 +23,13 @@ License Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA Class - SVD + Foam::SVD Description Singular value decomposition of a rectangular matrix. SourceFiles + SVDI.H SVD.C \*---------------------------------------------------------------------------*/ @@ -36,8 +37,7 @@ SourceFiles #ifndef SVD_H #define SVD_H -#include "DiagonalMatrix.H" -#include "Matrix.H" +#include "scalarMatrices.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -56,16 +56,16 @@ class SVD // Private data //- Rectangular matrix with the same dimensions as the input - Matrix U_; + scalarRectangularMatrix U_; //- square matrix V - Matrix V_; + scalarRectangularMatrix V_; //- The singular values DiagonalMatrix S_; //- The matrix product V S^(-1) U^T - Matrix VSinvUt_; + scalarRectangularMatrix VSinvUt_; //- The number of zero singular values label nZeros_; @@ -88,22 +88,22 @@ public: // Constructors //- Construct from a rectangular Matrix - SVD(const Matrix& A, const scalar minCondition = 0); + SVD(const scalarRectangularMatrix& A, const scalar minCondition = 0); // Access functions //- Return U - inline const Matrix& U() const; + inline const scalarRectangularMatrix& U() const; //- Return the square matrix V - inline const Matrix& V() const; + inline const scalarRectangularMatrix& V() const; //- Return the singular values - inline const DiagonalMatrix& S() const; + inline const scalarDiagonalMatrix& S() const; //- Return VSinvUt (the pseudo inverse) - inline const Matrix& VSinvUt() const; + inline const scalarRectangularMatrix& VSinvUt() const; //- Return the number of zero singular values inline label nZeros() const; diff --git a/src/OpenFOAM/matrices/scalarMatrix/SVD/SVDI.H b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVDI.H similarity index 88% rename from src/OpenFOAM/matrices/scalarMatrix/SVD/SVDI.H rename to src/OpenFOAM/matrices/scalarMatrices/SVD/SVDI.H index 337647edfb..298aac1aeb 100644 --- a/src/OpenFOAM/matrices/scalarMatrix/SVD/SVDI.H +++ b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVDI.H @@ -35,22 +35,22 @@ inline const T Foam::SVD::sign(const T& a, const T& b) // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -inline const Foam::Matrix& Foam::SVD::U() const +inline const Foam::scalarRectangularMatrix& Foam::SVD::U() const { return U_; } -inline const Foam::Matrix& Foam::SVD::V() const +inline const Foam::scalarRectangularMatrix& Foam::SVD::V() const { return V_; } -inline const Foam::DiagonalMatrix& Foam::SVD::S() const +inline const Foam::scalarDiagonalMatrix& Foam::SVD::S() const { return S_; } -inline const Foam::Matrix& Foam::SVD::VSinvUt() const +inline const Foam::scalarRectangularMatrix& Foam::SVD::VSinvUt() const { return VSinvUt_; } diff --git a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.C b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C similarity index 76% rename from src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.C rename to src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C index b88e167b6f..8ddd01a82c 100644 --- a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.C +++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C @@ -24,39 +24,14 @@ License \*---------------------------------------------------------------------------*/ -#include "scalarMatrix.H" +#include "scalarMatrices.H" #include "SVD.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 +void Foam::LUDecompose ( - Matrix& matrix, + scalarSquareMatrix& matrix, labelList& pivotIndices ) { @@ -81,8 +56,8 @@ void Foam::scalarMatrix::LUDecompose { FatalErrorIn ( - "scalarMatrix::LUdecompose" - "(Matrix& matrix, labelList& rowIndices)" + "LUdecompose" + "(scalarSquareMatrix& matrix, labelList& rowIndices)" ) << "Singular matrix" << exit(FatalError); } @@ -164,9 +139,9 @@ void Foam::scalarMatrix::LUDecompose void Foam::multiply ( - Matrix& ans, // value changed in return - const Matrix& A, - const Matrix& B + scalarRectangularMatrix& ans, // value changed in return + const scalarRectangularMatrix& A, + const scalarRectangularMatrix& B ) { if (A.m() != B.n()) @@ -174,15 +149,15 @@ void Foam::multiply FatalErrorIn ( "multiply(" - "Matrix& answer " - "const Matrix& A, " - "const Matrix& B)" + "scalarRectangularMatrix& answer " + "const scalarRectangularMatrix& A, " + "const scalarRectangularMatrix& B)" ) << "A and B must have identical inner dimensions but A.m = " << A.m() << " and B.n = " << B.n() << abort(FatalError); } - ans = Matrix(A.n(), B.m(), scalar(0)); + ans = scalarRectangularMatrix(A.n(), B.m(), scalar(0)); for(register label i = 0; i < A.n(); i++) { @@ -199,10 +174,10 @@ void Foam::multiply void Foam::multiply ( - Matrix& ans, // value changed in return - const Matrix& A, - const Matrix& B, - const Matrix& C + scalarRectangularMatrix& ans, // value changed in return + const scalarRectangularMatrix& A, + const scalarRectangularMatrix& B, + const scalarRectangularMatrix& C ) { if (A.m() != B.n()) @@ -210,10 +185,10 @@ void Foam::multiply FatalErrorIn ( "multiply(" - "const Matrix& A, " - "const Matrix& B, " - "const Matrix& C, " - "Matrix& answer)" + "const scalarRectangularMatrix& A, " + "const scalarRectangularMatrix& B, " + "const scalarRectangularMatrix& C, " + "scalarRectangularMatrix& answer)" ) << "A and B must have identical inner dimensions but A.m = " << A.m() << " and B.n = " << B.n() << abort(FatalError); @@ -224,16 +199,16 @@ void Foam::multiply FatalErrorIn ( "multiply(" - "const Matrix& A, " - "const Matrix& B, " - "const Matrix& C, " - "Matrix& answer)" + "const scalarRectangularMatrix& A, " + "const scalarRectangularMatrix& B, " + "const scalarRectangularMatrix& C, " + "scalarRectangularMatrix& answer)" ) << "B and C must have identical inner dimensions but B.m = " << B.m() << " and C.n = " << C.n() << abort(FatalError); } - ans = Matrix(A.n(), C.m(), scalar(0)); + ans = scalarRectangularMatrix(A.n(), C.m(), scalar(0)); for(register label i = 0; i < A.n(); i++) { @@ -255,10 +230,10 @@ void Foam::multiply void Foam::multiply ( - Matrix& ans, // value changed in return - const Matrix& A, + scalarRectangularMatrix& ans, // value changed in return + const scalarRectangularMatrix& A, const DiagonalMatrix& B, - const Matrix& C + const scalarRectangularMatrix& C ) { if (A.m() != B.size()) @@ -266,10 +241,10 @@ void Foam::multiply FatalErrorIn ( "multiply(" - "const Matrix& A, " + "const scalarRectangularMatrix& A, " "const DiagonalMatrix& B, " - "const Matrix& C, " - "Matrix& answer)" + "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); @@ -280,16 +255,16 @@ void Foam::multiply FatalErrorIn ( "multiply(" - "const Matrix& A, " + "const scalarRectangularMatrix& A, " "const DiagonalMatrix& B, " - "const Matrix& C, " - "Matrix& answer)" + "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 = Matrix(A.n(), C.m(), scalar(0)); + ans = scalarRectangularMatrix(A.n(), C.m(), scalar(0)); for(register label i = 0; i < A.n(); i++) { @@ -304,9 +279,9 @@ void Foam::multiply } -Foam::Matrix Foam::SVDinv +Foam::RectangularMatrix Foam::SVDinv ( - const Matrix& A, + const scalarRectangularMatrix& A, scalar minCondition ) { 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 90% rename from src/OpenFOAM/matrices/scalarMatrix/scalarMatrixTemplates.C rename to src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C index fb48a9b26d..131856b3fc 100644 --- a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrixTemplates.C +++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C @@ -24,15 +24,15 @@ License \*---------------------------------------------------------------------------*/ -#include "scalarMatrix.H" +#include "scalarMatrices.H" #include "Swap.H" // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template -void Foam::scalarMatrix::solve +void Foam::solve ( - Matrix& tmpMatrix, + scalarSquareMatrix& tmpMatrix, Field& sourceSol ) { @@ -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); } @@ -102,18 +102,23 @@ void Foam::scalarMatrix::solve template -void Foam::scalarMatrix::solve(Field& psi, const Field& source) const +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 +void Foam::LUBacksubstitute ( - const Matrix& luMatrix, + const scalarSquareMatrix& luMatrix, const labelList& pivotIndices, Field& sourceSol ) @@ -160,9 +165,9 @@ void Foam::scalarMatrix::LUBacksubstitute template -void Foam::scalarMatrix::LUsolve +void Foam::LUsolve ( - Matrix& matrix, + scalarSquareMatrix& matrix, Field& sourceSol ) { diff --git a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.H b/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.H deleted file mode 100644 index 2d22e71d4e..0000000000 --- a/src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.H +++ /dev/null @@ -1,155 +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 - -Class - Foam::scalarMatrix - -Description - Foam::scalarMatrix - -SourceFiles - scalarMatrix.C - -\*---------------------------------------------------------------------------*/ - -#ifndef scalarMatrix_H -#define scalarMatrix_H - -#include "Matrix.H" -#include "DiagonalMatrix.H" -#include "scalarField.H" -#include "labelList.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -/*---------------------------------------------------------------------------*\ - Class scalarMatrix Declaration -\*---------------------------------------------------------------------------*/ - -class scalarMatrix -: - public Matrix -{ - -public: - - // Constructors - - //- Construct null - scalarMatrix(); - - //- Construct given size - scalarMatrix(const label); - - //- Construct from Matrix - scalarMatrix(const Matrix&); - - //- Construct from Istream - scalarMatrix(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); -}; - - -// Global functions - -void multiply -( - Matrix& answer, // value changed in return - const Matrix& A, - const Matrix& B -); - -void multiply -( - Matrix& answer, // value changed in return - const Matrix& A, - const Matrix& B, - const Matrix& C -); - -void multiply -( - Matrix& answer, // value changed in return - const Matrix& A, - const DiagonalMatrix& B, - const Matrix& C -); - -//- Return the inverse of matrix A using SVD -Matrix SVDinv(const Matrix& A, scalar minCondition = 0); - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#ifdef NoRepository -# include "scalarMatrixTemplates.C" -#endif - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C index 1265b87ee5..f0f80750a5 100644 --- a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C +++ b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C @@ -31,7 +31,7 @@ License template Foam::simpleMatrix::simpleMatrix(const label mSize) : - scalarMatrix(mSize), + scalarSquareMatrix(mSize), source_(mSize, pTraits::zero) {} @@ -39,11 +39,11 @@ Foam::simpleMatrix::simpleMatrix(const label mSize) template Foam::simpleMatrix::simpleMatrix ( - const scalarMatrix& matrix, + const scalarSquareMatrix& matrix, const Field& source ) : - scalarMatrix(matrix), + scalarSquareMatrix(matrix), source_(source) {} @@ -51,7 +51,7 @@ Foam::simpleMatrix::simpleMatrix template Foam::simpleMatrix::simpleMatrix(Istream& is) : - scalarMatrix(is), + scalarSquareMatrix(is), source_(is) {} @@ -61,10 +61,10 @@ Foam::simpleMatrix::simpleMatrix(Istream& is) template Foam::Field Foam::simpleMatrix::solve() const { - scalarMatrix tmpMatrix = *this; + scalarSquareMatrix tmpMatrix = *this; Field sourceSol = source_; - scalarMatrix::solve(tmpMatrix, sourceSol); + Foam::solve(tmpMatrix, sourceSol); return sourceSol; } @@ -73,10 +73,10 @@ Foam::Field Foam::simpleMatrix::solve() const template Foam::Field Foam::simpleMatrix::LUsolve() const { - scalarMatrix luMatrix = *this; + scalarSquareMatrix luMatrix = *this; Field sourceSol = source_; - scalarMatrix::LUsolve(luMatrix, sourceSol); + Foam::LUsolve(luMatrix, sourceSol); return sourceSol; } @@ -108,7 +108,7 @@ void Foam::simpleMatrix::operator=(const simpleMatrix& m) << abort(FatalError); } - scalarMatrix::operator=(m); + scalarSquareMatrix::operator=(m); source_ = m.source_; } @@ -124,8 +124,8 @@ Foam::simpleMatrix Foam::operator+ { return simpleMatrix ( - static_cast(m1) - + static_cast(m2), + static_cast(m1) + + static_cast(m2), m1.source_ + m2.source_ ); } @@ -140,8 +140,8 @@ Foam::simpleMatrix Foam::operator- { return simpleMatrix ( - static_cast(m1) - - static_cast(m2), + static_cast(m1) + - static_cast(m2), m1.source_ - m2.source_ ); } @@ -159,7 +159,7 @@ Foam::simpleMatrix Foam::operator*(const scalar s, const simpleMatrix 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 d7b1a027be..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" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -48,27 +48,6 @@ namespace Foam 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 Ostream& operator<< ( @@ -84,7 +63,7 @@ Ostream& operator<< template class simpleMatrix : - public scalarMatrix + public scalarSquareMatrix { // Private data @@ -99,7 +78,7 @@ public: simpleMatrix(const label); //- Construct from components - simpleMatrix(const scalarMatrix&, const Field&); + simpleMatrix(const scalarSquareMatrix&, const Field&); //- Construct from Istream simpleMatrix(Istream&); @@ -137,27 +116,6 @@ public: 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& - ); - - // Ostream Operator friend Ostream& operator<< @@ -168,6 +126,30 @@ public: }; +// 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/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C index fd0e55fcdd..0caebcd421 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C @@ -49,11 +49,11 @@ Foam::quadraticFitData::quadraticFitData : MeshObject(mesh), centralWeight_(cWeight), -#ifdef SPHERICAL_GEOMETRY +# ifdef SPHERICAL_GEOMETRY dim_(2), -#else +# else dim_(mesh.nGeometricD()), -#endif +# endif minSize_ ( dim_ == 1 ? 3 : @@ -204,7 +204,7 @@ Foam::label Foam::quadraticFitData::calcFit scalar scale = 0; // calculate the matrix of the polynomial components - Matrix B(C.size(), minSize_, scalar(0)); + scalarRectangularMatrix B(C.size(), minSize_, scalar(0)); for(label ip = 0; ip < C.size(); ip++) { @@ -212,11 +212,11 @@ Foam::label Foam::quadraticFitData::calcFit scalar px = (p - p0)&idir; scalar py = (p - p0)&jdir; -#ifndef SPHERICAL_GEOMETRY +# ifndef SPHERICAL_GEOMETRY scalar pz = (p - p0)&kdir; -#else +# else scalar pz = mag(p) - mag(p0); -#endif +# endif if (ip == 0) { diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.C index 3dd49e126e..f2da51c746 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.C +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.C @@ -114,7 +114,7 @@ Foam::scalarField Foam::chemistryModel::omega forAll(reactions_, i) { const reaction& R = reactions_[i]; - + scalar omegai = omega ( R, c, T, p, pf, cf, lRef, pr, cr, rRef @@ -164,13 +164,13 @@ Foam::scalar Foam::chemistryModel::omega pf = 1.0; pr = 1.0; - + label Nl = R.lhs().size(); label Nr = R.rhs().size(); - + label slRef = 0; lRef = R.lhs()[slRef].index; - + pf = kf; for(label s=1; s& dfdc + scalarSquareMatrix& dfdc ) const { scalar T = c[Ns_]; scalar p = c[Ns_ + 1]; - + scalarField c2(Ns(), 0.0); for(label i=0; i Foam::chemistryModel::tc() const scalar pi = thermo_.p()[celli]; scalarField c(Ns_); scalar cSum = 0.0; - + for(label i=0; i& dfdc + scalarSquareMatrix& dfdc ) const; //- Calculates the reaction rates