Skip to content
GitLab
Commit c13421b1 authored by Henry Weller's avatar Henry Weller
Browse files

Matrix: Added (i, j) addressing to allow support for addressing blocks of the matrix 

This change brings OpenFOAM into line with the standard matrix
addressing in other C++ libraries for better interoperability.
parent 767ffc3c
Hide whitespace changes
Inline Side-by-side
applications/test/Matrix/Test-Matrix.C
View file @ c13421b1
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -36,15 +36,15 @@ int main(int argc, char *argv[])
{
SquareMatrix<scalar> hmm(3);
hmm[0][0] = -3.0;
hmm[0][1] = 10.0;
hmm[0][2] = -4.0;
hmm[1][0] = 2.0;
hmm[1][1] = 3.0;
hmm[1][2] = 10.0;
hmm[2][0] = 2.0;
hmm[2][1] = 6.0;
hmm[2][2] = 1.0;
hmm(0, 0) = -3.0;
hmm(0, 1) = 10.0;
hmm(0, 2) = -4.0;
hmm(1, 0) = 2.0;
hmm(1, 1) = 3.0;
hmm(1, 2) = 10.0;
hmm(2, 0) = 2.0;
hmm(2, 1) = 6.0;
hmm(2, 2) = 1.0;
//Info<< hmm << endl << hmm - 2.0*(-hmm) << endl;
Info<< max(hmm) << endl;
......@@ -106,15 +106,15 @@ int main(int argc, char *argv[])
{
scalarSquareMatrix squareMatrix(3, 3, 0);
squareMatrix[0][0] = 4;
squareMatrix[0][1] = 12;
squareMatrix[0][2] = -16;
squareMatrix[1][0] = 12;
squareMatrix[1][1] = 37;
squareMatrix[1][2] = -43;
squareMatrix[2][0] = -16;
squareMatrix[2][1] = -43;
squareMatrix[2][2] = 98;
squareMatrix(0, 0) = 4;
squareMatrix(0, 1) = 12;
squareMatrix(0, 2) = -16;
squareMatrix(1, 0) = 12;
squareMatrix(1, 1) = 37;
squareMatrix(1, 2) = -43;
squareMatrix(2, 0) = -16;
squareMatrix(2, 1) = -43;
squareMatrix(2, 2) = 98;
const scalarSquareMatrix squareMatrixCopy = squareMatrix;
Info<< nl << "Square Matrix = " << squareMatrix << endl;
......@@ -131,6 +131,22 @@ int main(int argc, char *argv[])
Info<< "det = " << detDecomposed(squareMatrix, sign) << endl;
}
{
scalarSquareMatrix squareMatrix(3000, 3000, 1);
for(label i=0; i<squareMatrix.n(); i++)
{
squareMatrix(i, i) = 10;
}
scalarField rhs(squareMatrix.n(), 0);
rhs[0] = 1;
rhs[1] = 2;
rhs[2] = 3;
LUsolve(squareMatrix, rhs);
}
Info<< "\nEnd\n" << endl;
return 0;
......
applications/test/ODE/Test-ODE.C
View file @ c13421b1
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -75,25 +75,25 @@ public:
dfdx[2] = (2.0/sqr(x))*y[2];
dfdx[3] = (3.0/sqr(x))*y[3];
dfdy[0][0] = 0.0;
dfdy[0][1] = -1.0;
dfdy[0][2] = 0.0;
dfdy[0][3] = 0.0;
dfdy[1][0] = 1.0;
dfdy[1][1] = -1.0/x;
dfdy[1][2] = 0.0;
dfdy[1][3] = 0.0;
dfdy[2][0] = 0.0;
dfdy[2][1] = 1.0;
dfdy[2][2] = -2.0/x;
dfdy[2][3] = 0.0;
dfdy[3][0] = 0.0;
dfdy[3][1] = 0.0;
dfdy[3][2] = 1.0;
dfdy[3][3] = -3.0/x;
dfdy(0, 0) = 0.0;
dfdy(0, 1) = -1.0;
dfdy(0, 2) = 0.0;
dfdy(0, 3) = 0.0;
dfdy(1, 0) = 1.0;
dfdy(1, 1) = -1.0/x;
dfdy(1, 2) = 0.0;
dfdy(1, 3) = 0.0;
dfdy(2, 0) = 0.0;
dfdy(2, 1) = 1.0;
dfdy(2, 2) = -2.0/x;
dfdy(2, 3) = 0.0;
dfdy(3, 0) = 0.0;
dfdy(3, 1) = 0.0;
dfdy(3, 2) = 1.0;
dfdy(3, 3) = -3.0/x;
}
};
......
applications/test/simpleMatrix/Test-simpleMatrix.C
View file @ c13421b1
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -35,15 +35,15 @@ int main(int argc, char *argv[])
{
simpleMatrix<vector> hmm(3);
hmm[0][0] = -3.0;
hmm[0][1] = 10.0;
hmm[0][2] = -4.0;
hmm[1][0] = 2.0;
hmm[1][1] = 3.0;
hmm[1][2] = 10.0;
hmm[2][0] = 2.0;
hmm[2][1] = 6.0;
hmm[2][2] = 1.0;
hmm(0, 0) = -3.0;
hmm(0, 1) = 10.0;
hmm(0, 2) = -4.0;
hmm(1, 0) = 2.0;
hmm(1, 1) = 3.0;
hmm(1, 2) = 10.0;
hmm(2, 0) = 2.0;
hmm(2, 1) = 6.0;
hmm(2, 2) = 1.0;
hmm.source()[0] = vector(2.0, 1.0, 3.0);
hmm.source()[1] = vector(1.0, 4.0, 3.0);
......
applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
View file @ c13421b1
......@@ -312,11 +312,11 @@ triSurfacePointScalarField calcCurvature
scalar x = edgeVectors[i] & faceCoordSys[0];
scalar y = edgeVectors[i] & faceCoordSys[1];
T[0][0] += sqr(x);
T[1][0] += x*y;
T[1][1] += sqr(x) + sqr(y);
T[2][1] += x*y;
T[2][2] += sqr(y);
T(0, 0) += sqr(x);
T(1, 0) += x*y;
T(1, 1) += sqr(x) + sqr(y);
T(2, 1) += x*y;
T(2, 2) += sqr(y);
scalar dndx = normalDifferences[i] & faceCoordSys[0];
scalar dndy = normalDifferences[i] & faceCoordSys[1];
......
src/ODE/ODESolvers/EulerSI/EulerSI.C
View file @ c13421b1
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013-2015 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -67,10 +67,10 @@ Foam::scalar Foam::EulerSI::solve
{
for (label j=0; j<n_; j++)
{
a_[i][j] = -dfdy_[i][j];
a_(i, j) = -dfdy_(i, j);
}
a_[i][i] += 1.0/dx;
a_(i, i) += 1.0/dx;
}
LUDecompose(a_, pivotIndices_);
......
src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.C
View file @ c13421b1
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013-2015 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -81,10 +81,10 @@ Foam::scalar Foam::Rosenbrock12::solve
{
for (label j=0; j<n_; j++)
{
a_[i][j] = -dfdy_[i][j];
a_(i, j) = -dfdy_(i, j);
}
a_[i][i] += 1.0/(gamma*dx);
a_(i, i) += 1.0/(gamma*dx);
}
LUDecompose(a_, pivotIndices_);
......
src/ODE/ODESolvers/Rosenbrock23/Rosenbrock23.C
View file @ c13421b1
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013-2015 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -94,10 +94,10 @@ Foam::scalar Foam::Rosenbrock23::solve
{
for (label j=0; j<n_; j++)
{
a_[i][j] = -dfdy_[i][j];
a_(i, j) = -dfdy_(i, j);
}
a_[i][i] += 1.0/(gamma*dx);
a_(i, i) += 1.0/(gamma*dx);
}
LUDecompose(a_, pivotIndices_);
......
src/ODE/ODESolvers/Rosenbrock34/Rosenbrock34.C
View file @ c13421b1
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013-2015 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -139,10 +139,10 @@ Foam::scalar Foam::Rosenbrock34::solve
{
for (label j=0; j<n_; j++)
{
a_[i][j] = -dfdy_[i][j];
a_(i, j) = -dfdy_(i, j);
}
a_[i][i] += 1.0/(gamma*dx);
a_(i, i) += 1.0/(gamma*dx);
}
LUDecompose(a_, pivotIndices_);
......
src/ODE/ODESolvers/SIBS/SIMPR.C
View file @ c13421b1
......@@ -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
......@@ -46,9 +46,9 @@ void Foam::SIBS::SIMPR
{
for (label j=0; j<n_; j++)
{
a[i][j] = -h*dfdy[i][j];
a(i, j) = -h*dfdy(i, j);
}
++a[i][i];
++a(i, i);
}
labelList pivotIndices(n_);
......
src/ODE/ODESolvers/SIBS/polyExtrapolate.C
View file @ c13421b1
......@@ -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
......@@ -51,7 +51,7 @@ void Foam::SIBS::polyExtrapolate
{
for (label j=0; j<n; j++)
{
d[j][0] = yest[j];
d(j, 0) = yest[j];
}
}
else
......
src/ODE/ODESolvers/rodas23/rodas23.C
View file @ c13421b1
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013-2015 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -85,10 +85,10 @@ Foam::scalar Foam::rodas23::solve
{
for (label j=0; j<n_; j++)
{
a_[i][j] = -dfdy_[i][j];
a_(i, j) = -dfdy_(i, j);
}
a_[i][i] += 1.0/(gamma*dx);
a_(i, i) += 1.0/(gamma*dx);
}
LUDecompose(a_, pivotIndices_);
......
src/ODE/ODESolvers/rodas34/rodas34.C
View file @ c13421b1
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013-2015 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -108,10 +108,10 @@ Foam::scalar Foam::rodas34::solve
{
for (label j=0; j<n_; j++)
{
a_[i][j] = -dfdy_[i][j];
a_(i, j) = -dfdy_(i, j);
}
a_[i][i] += 1.0/(gamma*dx);
a_(i, i) += 1.0/(gamma*dx);
}
LUDecompose(a_, pivotIndices_);
......
src/ODE/ODESolvers/seulex/seulex.C
View file @ c13421b1
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013-2015 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -92,7 +92,7 @@ Foam::seulex::seulex(const ODESystem& ode, const dictionary& dict)
for (int l=0; l<k; l++)
{
scalar ratio = scalar(nSeq_[k])/nSeq_[l];
coeff_[k][l] = 1.0/(ratio - 1.0);
coeff_(k, l) = 1.0/(ratio - 1.0);
}
}
}
......@@ -117,10 +117,10 @@ bool Foam::seulex::seul
{
for (label j=0; j<n_; j++)
{
a_[i][j] = -dfdy_[i][j];
a_(i, j) = -dfdy_(i, j);
}
a_[i][i] += 1.0/dx;
a_(i, i) += 1.0/dx;
}
LUDecompose(a_, pivotIndices_);
......@@ -192,13 +192,13 @@ void Foam::seulex::extrapolate
for (label i=0; i<n_; i++)
{
table[j-1][i] =
table[j][i] + coeff_[k][j]*(table[j][i] - table[j-1][i]);
table(j, i) + coeff_(k, j)*(table(j, i) - table[j-1][i]);
}
}
for (int i=0; i<n_; i++)
{
y[i] = table[0][i] + coeff_[k][0]*(table[0][i] - y[i]);
y[i] = table(0, i) + coeff_(k, 0)*(table(0, i) - y[i]);
}
}
......@@ -288,7 +288,7 @@ void Foam::seulex::solve
forAll(scale_, i)
{
scale_[i] = absTol_[i] + relTol_[i]*mag(y0_[i]);
err += sqr((y[i] - table_[0][i])/scale_[i]);
err += sqr((y[i] - table_(0, i))/scale_[i]);
}
err = sqrt(err/n_);
if (err > 1.0/SMALL || (k > 1 && err >= errOld))
......
src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C
View file @ c13421b1
......@@ -42,7 +42,7 @@ Foam::DiagonalMatrix<Type>::DiagonalMatrix(const Matrix<Form, Type>& a)
{
forAll(*this, i)
{
this->operator[](i) = a[i][i];
this->operator[](i) = a(i, i);
}
}
......
src/OpenFOAM/matrices/Matrix/Matrix.C
View file @ c13421b1
......@@ -156,7 +156,7 @@ Form Foam::Matrix<Form, Type>::T() const
{
for (label j=0; j<n(); j++)
{
At[j][i] = A[i][j];
At(j, i) = A(i, j);
}
}
......@@ -238,7 +238,7 @@ const Type& Foam::max(const Matrix<Form, Type>& a)
<< abort(FatalError);
// Return in error to keep compiler happy
return a[0][0];
return a(0, 0);
}
}
......@@ -270,7 +270,7 @@ const Type& Foam::min(const Matrix<Form, Type>& a)
<< abort(FatalError);
// Return in error to keep compiler happy
return a[0][0];
return a(0, 0);
}
}
......@@ -411,7 +411,7 @@ Form Foam::operator*(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
{
for (label l = 0; l < b.m(); l++)
{
ab[i][j] += a[i][l]*b[l][j];
ab(i, j) += a(i, l)*b(l, j);
}
}
}
......
src/OpenFOAM/matrices/Matrix/MatrixI.H
View file @ c13421b1
......@@ -136,7 +136,7 @@ inline const Type& Foam::Matrix<Form, Type>::operator()
) const
{
checki(i);
checki(j);
checkj(j);
return v_[i*nCols_ + j];
}
......@@ -149,7 +149,7 @@ inline Type& Foam::Matrix<Form, Type>::operator()
)
{
checki(i);
checki(j);
checkj(j);
return v_[i*nCols_ + j];
}
......
src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C
View file @ c13421b1
......@@ -39,7 +39,7 @@ Foam::scalar Foam::detDecomposed
for (label i = 0; i < matrix.m(); ++i)
{
diagProduct *= matrix[i][i];
diagProduct *= matrix(i, i);
}
return sign*diagProduct;
......
src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.C
View file @ c13421b1
......@@ -37,7 +37,7 @@ Foam::SymmetricSquareMatrix<Type> Foam::invDecomposed
for (label i = 0; i < matrix.m(); ++i)
{
inv[i][i] = 1.0/matrix[i][i];
inv(i, i) = 1.0/matrix(i, i);
for (label j = 0; j < i; ++j)
{
......@@ -45,10 +45,10 @@ Foam::SymmetricSquareMatrix<Type> Foam::invDecomposed
for (label k = j; k < i; k++)
{
sum -= matrix[i][k]*inv[k][j];
sum -= matrix(i, k)*inv(k, j);
}
inv[i][j] = sum/matrix[i][i];
inv(i, j) = sum/matrix(i, i);
}
}
......@@ -77,7 +77,7 @@ Foam::scalar Foam::detDecomposed(const SymmetricSquareMatrix<Type>& matrix)
for (label i = 0; i < matrix.m(); ++i)
{
diagProduct *= matrix[i][i];
diagProduct *= matrix(i, i);
}
return sqr(diagProduct);
......
src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H
View file @ c13421b1
......@@ -76,13 +76,6 @@ public:
//- Clone
inline autoPtr<SymmetricSquareMatrix<Type>> clone() const;
//- Return subscript-checked row of Matrix.
inline Type& operator()(const label r, const label c);
//- Return subscript-checked row of constant Matrix.
inline const Type& operator()(const label r, const label c) const;
};
......
src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrixI.H
View file @ c13421b1
......@@ -94,40 +94,4 @@ Foam::SymmetricSquareMatrix<Type>::clone() const
}
template<class Type>
inline Type& Foam::SymmetricSquareMatrix<Type>::operator()
(
const label r,
const label c
)
{
if (r > c)
{
return this->operator[](r)[c];
}
else
{
return this->operator[](c)[r];
}
}
template<class Type>
inline const Type& Foam::SymmetricSquareMatrix<Type>::operator()
(
const label r,
const label c
) const
{
if (r > c)
{
return this->operator[](r)[c];
}
else
{