Commit caf8776f authored by Henry Weller's avatar Henry Weller
Browse files

SquareMatrix, SymmetricSquareMatrix: Changed the constructor from size to require only n

This avoids the need to check that the m and n dimensions are the same.
parent 08226143
......@@ -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;
......
......@@ -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());
......
......@@ -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
......
......@@ -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)
{}
......
......@@ -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)
{
......
......@@ -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&);
......
......@@ -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
......
......@@ -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);
......
......@@ -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>&);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......
......@@ -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>
......
......@@ -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_)
......
......@@ -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++)
{
......
......@@ -43,7 +43,7 @@ Foam::simpleMatrix<Type>::simpleMatrix
const Type& sourceVal
)
:
scalarSquareMatrix(mSize, mSize, coeffVal),
scalarSquareMatrix(mSize, coeffVal),
source_(mSize, sourceVal)
{}
......
......@@ -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++)
{
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment