diff --git a/applications/test/Matrix/Test-Matrix.C b/applications/test/Matrix/Test-Matrix.C
index ad48481fa1504ec94b4df8224ad1a3be0bd749a3..95cfa4c1af5f90da66af33f657cf1d438cd3713a 100644
--- a/applications/test/Matrix/Test-Matrix.C
+++ b/applications/test/Matrix/Test-Matrix.C
@@ -25,6 +25,7 @@ License
#include "scalarMatrices.H"
#include "vector.H"
+#include "tensor.H"
#include "IFstream.H"
using namespace Foam;
@@ -50,7 +51,7 @@ int main(int argc, char *argv[])
Info<< max(hmm) << endl;
Info<< min(hmm) << endl;
- SquareMatrix<scalar> hmm2(3, 3, 1.0);
+ SquareMatrix<scalar> hmm2(3, I);
hmm = hmm2;
@@ -72,7 +73,14 @@ int main(int argc, char *argv[])
Info<< hmm5 << endl;
{
- scalarSymmetricSquareMatrix symmMatrix(3, 3, 0);
+ RectangularMatrix<scalar> rm1(5, 6, 3.1);
+ rm1(0, 1) = 4.5;
+ RectangularMatrix<scalar> rm1b(rm1.block(2, 2, 0, 0));
+ Info<< "rm1b = " << rm1b << endl;
+ }
+
+ {
+ scalarSymmetricSquareMatrix symmMatrix(3, Zero);
symmMatrix(0, 0) = 4;
symmMatrix(1, 0) = 12;
@@ -104,7 +112,7 @@ int main(int argc, char *argv[])
}
{
- scalarSquareMatrix squareMatrix(3, 3, 0);
+ scalarSquareMatrix squareMatrix(3, Zero);
squareMatrix(0, 0) = 4;
squareMatrix(0, 1) = 12;
diff --git a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
index ed00d1a6cf33deb59aeff113a76b69b092c1dd9e..d5f295a62b5517a8f6739787115bb5246f70cc9c 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
+++ b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
@@ -675,7 +675,6 @@ int main(int argc, char *argv[])
// Matrix sum in j(Fij) for each i (if enclosure sum = 1)
scalarSquareMatrix sumViewFactorPatch
(
- totalPatches,
totalPatches,
0.0
);
@@ -888,7 +887,7 @@ int main(int argc, char *argv[])
if (Pstream::master())
{
- scalarSquareMatrix Fmatrix(totalNCoarseFaces, totalNCoarseFaces, 0.0);
+ scalarSquareMatrix Fmatrix(totalNCoarseFaces, 0.0);
labelListList globalFaceFaces(visibleFaceFaces.size());
diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
index 510305f28cfa0a58be73a3eca105391e6995a186..2882e6af2af442d4ea8039b204e7ad0e6577e988 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
@@ -303,7 +303,7 @@ triSurfacePointScalarField calcCurvature
faceCoordSys.normalize();
// Construct the matrix to solve
- scalarSymmetricSquareMatrix T(3, 3, 0);
+ scalarSymmetricSquareMatrix T(3, 0);
scalarDiagonalMatrix Z(3, 0);
// Least Squares
diff --git a/src/ODE/ODESolvers/SIBS/SIBS.C b/src/ODE/ODESolvers/SIBS/SIBS.C
index 36acb8a6b1abc36f54ba485a36c74f89ec4831b5..a095a7c50738b4c9b066bdbcd8631bd3a72751ed 100644
--- a/src/ODE/ODESolvers/SIBS/SIBS.C
+++ b/src/ODE/ODESolvers/SIBS/SIBS.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -50,7 +50,7 @@ Foam::SIBS::SIBS(const ODESystem& ode, const dictionary& dict)
:
ODESolver(ode, dict),
a_(iMaxX_, 0.0),
- alpha_(kMaxX_, kMaxX_, 0.0),
+ alpha_(kMaxX_, 0.0),
d_p_(n_, kMaxX_, 0.0),
x_p_(kMaxX_, 0.0),
err_(kMaxX_, 0.0),
@@ -60,7 +60,7 @@ Foam::SIBS::SIBS(const ODESystem& ode, const dictionary& dict)
yErr_(n_, 0.0),
dydx0_(n_),
dfdx_(n_, 0.0),
- dfdy_(n_, n_, 0.0),
+ dfdy_(n_, 0.0),
first_(1),
epsOld_(-1.0)
{}
diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C
index 74439c15defe0140c9e4e72725583d1638b5de74..cf34d78a212365979368df03cd06145458fbb9c4 100644
--- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C
+++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C
@@ -35,7 +35,7 @@ Foam::scalar Foam::detDecomposed
const label sign
)
{
- scalar diagProduct = 1.0;
+ Type diagProduct = pTraits<Type>::one;
for (label i = 0; i < matrix.m(); ++i)
{
diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H
index 270891a9551dd2a7d04029218d45b2a099a7297b..39403194d55be3cdbc8b2deeb206b5a0c8534fbe 100644
--- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H
+++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H
@@ -38,6 +38,7 @@ SourceFiles
#define SquareMatrix_H
#include "Matrix.H"
+#include "Identity.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -64,14 +65,17 @@ public:
//- 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 given number of rows/columns
+ // initializing all elements to zero
+ inline SquareMatrix(const label n, const zero);
+
+ //- Construct given number of rows/columns
+ // Initializing to the identity matrix
+ inline SquareMatrix(const label n, const Identity<Type>);
//- 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&);
+ // initializing all elements to the given value
+ inline SquareMatrix(const label n, const Type&);
//- Construct from Istream.
inline SquareMatrix(Istream&);
diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H
index 74d7bdb95c2df48fad77496dc47c054d1bc36767..ffa595d51da4c50bc68eab9dcc2e9cc67b62cd75 100644
--- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H
+++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H
@@ -31,40 +31,51 @@ inline Foam::SquareMatrix<Type>::SquareMatrix()
Matrix<SquareMatrix<Type>, Type>()
{}
+
template<class Type>
inline Foam::SquareMatrix<Type>::SquareMatrix(const label n)
:
Matrix<SquareMatrix<Type>, Type>(n, n)
{}
+
template<class Type>
-inline Foam::SquareMatrix<Type>::SquareMatrix(const label m, const label n)
+inline Foam::SquareMatrix<Type>::SquareMatrix
+(
+ const label n,
+ const zero
+)
:
- Matrix<SquareMatrix<Type>, Type>(m, n)
+ Matrix<SquareMatrix<Type>, Type>(n, n, Zero)
+{}
+
+
+template<class Type>
+inline Foam::SquareMatrix<Type>::SquareMatrix
+(
+ const label n,
+ const Identity<Type>
+)
+:
+ Matrix<SquareMatrix<Type>, Type>(n, n, Zero)
{
- if (m != n)
+ for (label i = 0; i < n; ++i)
{
- FatalErrorInFunction
- << "m != n for constructing a square matrix" << exit(FatalError);
+ this->operator()(i, i) = I;
}
}
+
template<class Type>
inline Foam::SquareMatrix<Type>::SquareMatrix
(
- const label m,
const label n,
const Type& t
)
:
- Matrix<SquareMatrix<Type>, Type>(m, n, t)
-{
- if (m != n)
- {
- FatalErrorInFunction
- << "m != n for constructing a square matrix" << exit(FatalError);
- }
-}
+ Matrix<SquareMatrix<Type>, Type>(n, n, t)
+{}
+
template<class Type>
inline Foam::SquareMatrix<Type>::SquareMatrix(Istream& is)
@@ -72,6 +83,7 @@ inline Foam::SquareMatrix<Type>::SquareMatrix(Istream& is)
Matrix<SquareMatrix<Type>, Type>(is)
{}
+
template<class Type>
inline Foam::autoPtr<Foam::SquareMatrix<Type>>
Foam::SquareMatrix<Type>::clone() const
diff --git a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.C b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.C
index 28ad1d2745d0d75eb1b7fa1b0cb643b3b299c6ed..dcb31ef24136242204836ecf52be5ef602f4e5a0 100644
--- a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.C
+++ b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.C
@@ -33,15 +33,17 @@ Foam::SymmetricSquareMatrix<Type> Foam::invDecomposed
const SymmetricSquareMatrix<Type>& matrix
)
{
- SymmetricSquareMatrix<Type> inv(matrix.m(), matrix.m(), 0.0);
+ const label n = matrix.n();
- for (label i = 0; i < matrix.m(); ++i)
+ SymmetricSquareMatrix<Type> inv(n, Zero);
+
+ for (label i = 0; i < n; ++i)
{
inv(i, i) = 1.0/matrix(i, i);
for (label j = 0; j < i; ++j)
{
- scalar sum = 0.0;
+ Type sum = Zero;
for (label k = j; k < i; k++)
{
@@ -52,7 +54,20 @@ Foam::SymmetricSquareMatrix<Type> Foam::invDecomposed
}
}
- return inv.T()*inv;
+ SymmetricSquareMatrix<Type> result(n, Zero);
+
+ for (label k = 0; k < n; k++)
+ {
+ for (label i = 0; i <= k; i++)
+ {
+ for (label j = 0; j <= k; j++)
+ {
+ result(i, j) += inv(k, i)*inv(k, j);
+ }
+ }
+ }
+
+ return result;
}
@@ -63,7 +78,6 @@ Foam::SymmetricSquareMatrix<Type> Foam::inv
)
{
SymmetricSquareMatrix<Type> matrixTmp(matrix);
-
LUDecompose(matrixTmp);
return invDecomposed(matrixTmp);
@@ -71,9 +85,9 @@ Foam::SymmetricSquareMatrix<Type> Foam::inv
template<class Type>
-Foam::scalar Foam::detDecomposed(const SymmetricSquareMatrix<Type>& matrix)
+Type Foam::detDecomposed(const SymmetricSquareMatrix<Type>& matrix)
{
- scalar diagProduct = 1.0;
+ Type diagProduct = pTraits<Type>::one;
for (label i = 0; i < matrix.m(); ++i)
{
@@ -85,10 +99,9 @@ Foam::scalar Foam::detDecomposed(const SymmetricSquareMatrix<Type>& matrix)
template<class Type>
-Foam::scalar Foam::det(const SymmetricSquareMatrix<Type>& matrix)
+Type Foam::det(const SymmetricSquareMatrix<Type>& matrix)
{
- SymmetricSquareMatrix<Type> matrixTmp = matrix;
-
+ SymmetricSquareMatrix<Type> matrixTmp(matrix);
LUDecompose(matrixTmp);
return detDecomposed(matrixTmp);
diff --git a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H
index 193ee1602f26369b0ae445ee86e5cdd1ecb7f1c9..543d16b8482fe16846cd03a2e3dc12278a1ec367 100644
--- a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H
+++ b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H
@@ -38,6 +38,7 @@ SourceFiles
#define SymmetricSquareMatrix_H
#include "SquareMatrix.H"
+#include "Identity.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -64,12 +65,15 @@ public:
//- Construct given number of rows/columns.
inline SymmetricSquareMatrix(const label n);
- //- Construct with given number of rows/columns
- inline SymmetricSquareMatrix(const label m, const label n);
+ //- Construct given number of rows/columns, initializing to zero
+ inline SymmetricSquareMatrix(const label n, const zero);
+
+ //- Construct given number of rows/columns,
+ inline SymmetricSquareMatrix(const label n, const Identity<Type>);
//- Construct with given number of rows/columns
- // and value for all elements.
- inline SymmetricSquareMatrix(const label m, const label n, const Type&);
+ // initializing all elements to the given value
+ inline SymmetricSquareMatrix(const label n, const Type&);
//- Construct from Istream.
inline SymmetricSquareMatrix(Istream&);
@@ -91,11 +95,11 @@ SymmetricSquareMatrix<Type> inv(const SymmetricSquareMatrix<Type>&);
//- Return the LU decomposed SymmetricSquareMatrix det
template<class Type>
-scalar detDecomposed(const SymmetricSquareMatrix<Type>&);
+Type detDecomposed(const SymmetricSquareMatrix<Type>&);
//- Return the SymmetricSquareMatrix det
template<class Type>
-scalar det(const SymmetricSquareMatrix<Type>&);
+Type det(const SymmetricSquareMatrix<Type>&);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrixI.H b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrixI.H
index a363026aac50e74a74bc8e1e957f19504126255d..821daf199bd2dd38cd295dce8e8700f202bcc5e0 100644
--- a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrixI.H
+++ b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrixI.H
@@ -42,17 +42,26 @@ inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix(const label n)
template<class Type>
inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix
(
- const label m,
- const label n
+ const label n,
+ const zero
)
:
- Matrix<SymmetricSquareMatrix<Type>, Type>(m, n)
+ Matrix<SymmetricSquareMatrix<Type>, Type>(n, n, Zero)
+{}
+
+
+template<class Type>
+inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix
+(
+ const label n,
+ const Identity<Type>
+)
+:
+ Matrix<SymmetricSquareMatrix<Type>, Type>(n, n, Zero)
{
- if (m != n)
+ for (label i = 0; i < n; i++)
{
- FatalErrorInFunction
- << "m != n for constructing a symmetric square matrix"
- << exit(FatalError);
+ this->operator()(i, i) = pTraits<Type>::one;
}
}
@@ -60,20 +69,12 @@ inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix
template<class Type>
inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix
(
- const label m,
const label n,
const Type& t
)
:
- Matrix<SymmetricSquareMatrix<Type>, Type>(m, n, t)
-{
- if (m != n)
- {
- FatalErrorInFunction
- << "m != n for constructing a symmetric square matrix"
- << exit(FatalError);
- }
-}
+ Matrix<SymmetricSquareMatrix<Type>, Type>(n, n, t)
+{}
template<class Type>
diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixTemplates.C b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixTemplates.C
index 739251c35879478fb47ac4112357846dbafa6939..7ac02d189f35dc172ad90815fc8dc959c6b400ae 100644
--- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixTemplates.C
+++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixTemplates.C
@@ -35,7 +35,7 @@ Foam::tmp<Foam::Field<Type>> Foam::lduMatrix::H(const Field<Type>& psi) const
{
tmp<Field<Type>> tHpsi
(
- new Field<Type>(lduAddr().size(), pTraits<Type>::zero)
+ new Field<Type>(lduAddr().size(), Zero)
);
if (lowerPtr_ || upperPtr_)
diff --git a/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C b/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C
index c6acc739e27dac7dbcc7d01b852c1b1e75cd3225..52413a38ddad9fc4b246a3a3e239a26bd27e0d74 100644
--- a/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C
+++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C
@@ -87,7 +87,7 @@ void Foam::solve
// Back-substitution
for (label j=m-1; j>=0; j--)
{
- Type ntempvec = pTraits<Type>::zero;
+ Type ntempvec = Zero;
for (label k=j+1; k<m; k++)
{
diff --git a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C
index 960ba0d563a90dd70e171f0725804bf363e3a154..0f1c6995b981bd1e9fc86d80e5824e73f767e17c 100644
--- a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C
+++ b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C
@@ -43,7 +43,7 @@ Foam::simpleMatrix<Type>::simpleMatrix
const Type& sourceVal
)
:
- scalarSquareMatrix(mSize, mSize, coeffVal),
+ scalarSquareMatrix(mSize, coeffVal),
source_(mSize, sourceVal)
{}
diff --git a/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C b/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C
index 2c8a7e2ba10252526763e0bf7e18c3e652e8d23a..d70eccfa9c306f5a59f51aaa0657819921ef5a16 100644
--- a/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C
+++ b/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C
@@ -173,7 +173,7 @@ void Foam::radiation::viewFactor::initialise()
{
Fmatrix_.reset
(
- new scalarSquareMatrix(totalNCoarseFaces_, totalNCoarseFaces_, 0.0)
+ new scalarSquareMatrix(totalNCoarseFaces_, 0.0)
);
if (debug)
@@ -224,12 +224,7 @@ void Foam::radiation::viewFactor::initialise()
{
CLU_.reset
(
- new scalarSquareMatrix
- (
- totalNCoarseFaces_,
- totalNCoarseFaces_,
- 0.0
- )
+ new scalarSquareMatrix(totalNCoarseFaces_, 0.0)
);
pivotIndices_.setSize(CLU_().m());
@@ -529,7 +524,7 @@ void Foam::radiation::viewFactor::calculate()
// Variable emissivity
if (!constEmissivity_)
{
- scalarSquareMatrix C(totalNCoarseFaces_, totalNCoarseFaces_, 0.0);
+ scalarSquareMatrix C(totalNCoarseFaces_, 0.0);
for (label i=0; i<totalNCoarseFaces_; i++)
{