From 878866b16e432faa3febb0bfd0f3f82534b51d39 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Wed, 23 Mar 2016 12:50:34 +0000
Subject: [PATCH] MatrixBlock: Separate Matrix::Block into the separate class
MatrixBlock
This avoids serious problems with template parameter deduction when
manipulating blocks of different matrix types e.g. Square and
Rectangular.
---
src/OpenFOAM/matrices/Matrix/Matrix.C | 562 +++++++-----------
src/OpenFOAM/matrices/Matrix/Matrix.H | 214 ++-----
src/OpenFOAM/matrices/Matrix/MatrixI.H | 220 +------
src/OpenFOAM/matrices/Matrix/MatrixIO.C | 8 +-
.../matrices/MatrixBlock/MatrixBlock.C | 342 +++++++++++
.../matrices/MatrixBlock/MatrixBlock.H | 240 ++++++++
.../matrices/MatrixBlock/MatrixBlockI.H | 203 +++++++
7 files changed, 1060 insertions(+), 729 deletions(-)
create mode 100644 src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.C
create mode 100644 src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.H
create mode 100644 src/OpenFOAM/matrices/MatrixBlock/MatrixBlockI.H
diff --git a/src/OpenFOAM/matrices/Matrix/Matrix.C b/src/OpenFOAM/matrices/Matrix/Matrix.C
index 808e00b86c2..46508d8e3f6 100644
--- a/src/OpenFOAM/matrices/Matrix/Matrix.C
+++ b/src/OpenFOAM/matrices/Matrix/Matrix.C
@@ -30,7 +30,7 @@ License
template<class Form, class Type>
void Foam::Matrix<Form, Type>::allocate()
{
- if (nRows_ && nCols_)
+ if (mRows_ && nCols_)
{
v_ = new Type[size()];
}
@@ -42,14 +42,14 @@ void Foam::Matrix<Form, Type>::allocate()
template<class Form, class Type>
Foam::Matrix<Form, Type>::Matrix(const label m, const label n)
:
- nRows_(m),
+ mRows_(m),
nCols_(n),
v_(NULL)
{
- if (nRows_ < 0 || nCols_ < 0)
+ if (mRows_ < 0 || nCols_ < 0)
{
FatalErrorInFunction
- << "Bad m, n " << nRows_ << ", " << nCols_
+ << "Incorrect m, n " << mRows_ << ", " << nCols_
<< abort(FatalError);
}
@@ -60,14 +60,14 @@ Foam::Matrix<Form, Type>::Matrix(const label m, const label n)
template<class Form, class Type>
Foam::Matrix<Form, Type>::Matrix(const label m, const label n, const zero)
:
- nRows_(m),
+ mRows_(m),
nCols_(n),
v_(NULL)
{
- if (nRows_ < 0 || nCols_ < 0)
+ if (mRows_ < 0 || nCols_ < 0)
{
FatalErrorInFunction
- << "Bad m, n " << nRows_ << ", " << nCols_
+ << "Incorrect m, n " << mRows_ << ", " << nCols_
<< abort(FatalError);
}
@@ -85,16 +85,16 @@ Foam::Matrix<Form, Type>::Matrix(const label m, const label n, const zero)
template<class Form, class Type>
-Foam::Matrix<Form, Type>::Matrix(const label m, const label n, const Type& a)
+Foam::Matrix<Form, Type>::Matrix(const label m, const label n, const Type& s)
:
- nRows_(m),
+ mRows_(m),
nCols_(n),
v_(NULL)
{
- if (nRows_ < 0 || nCols_ < 0)
+ if (mRows_ < 0 || nCols_ < 0)
{
FatalErrorInFunction
- << "Bad m, n " << nRows_ << ", " << nCols_
+ << "Incorrect m, n " << mRows_ << ", " << nCols_
<< abort(FatalError);
}
@@ -105,45 +105,92 @@ Foam::Matrix<Form, Type>::Matrix(const label m, const label n, const Type& a)
const label mn = size();
for (label i=0; i<mn; i++)
{
- v_[i] = a;
+ v_[i] = s;
}
}
}
template<class Form, class Type>
-inline Foam::Matrix<Form, Type>::Matrix(const mType::Block& block)
+Foam::Matrix<Form, Type>::Matrix(const Matrix<Form, Type>& M)
:
- nRows_(block.m()),
- nCols_(block.n())
+ mRows_(M.mRows_),
+ nCols_(M.nCols_),
+ v_(NULL)
{
- allocate();
-
- for (label i=0; i<nRows_; i++)
+ if (M.v_)
{
- for (label j=0; j<nCols_; j++)
+ allocate();
+
+ const label mn = size();
+ for (label i=0; i<mn; i++)
{
- (*this)(i,j) = block(i,j);
+ v_[i] = M.v_[i];
}
}
}
template<class Form, class Type>
-Foam::Matrix<Form, Type>::Matrix(const Matrix<Form, Type>& a)
+template<class Form2>
+Foam::Matrix<Form, Type>::Matrix(const Matrix<Form2, Type>& M)
:
- nRows_(a.nRows_),
- nCols_(a.nCols_),
+ mRows_(M.m()),
+ nCols_(M.n()),
v_(NULL)
{
- if (a.v_)
+ if (M.v())
{
allocate();
const label mn = size();
for (label i=0; i<mn; i++)
{
- v_[i] = a.v_[i];
+ v_[i] = M.v()[i];
+ }
+ }
+}
+
+
+template<class Form, class Type>
+template<class MatrixType>
+inline Foam::Matrix<Form, Type>::Matrix
+(
+ const ConstMatrixBlock<MatrixType>& Mb
+)
+:
+ mRows_(Mb.m()),
+ nCols_(Mb.n())
+{
+ allocate();
+
+ for (label i=0; i<mRows_; i++)
+ {
+ for (label j=0; j<nCols_; j++)
+ {
+ (*this)(i,j) = Mb(i,j);
+ }
+ }
+}
+
+
+template<class Form, class Type>
+template<class MatrixType>
+inline Foam::Matrix<Form, Type>::Matrix
+(
+ const MatrixBlock<MatrixType>& Mb
+)
+:
+ mRows_(Mb.m()),
+ nCols_(Mb.n())
+{
+ allocate();
+
+ for (label i=0; i<mRows_; i++)
+ {
+ for (label j=0; j<nCols_; j++)
+ {
+ (*this)(i,j) = Mb(i,j);
}
}
}
@@ -172,24 +219,24 @@ void Foam::Matrix<Form, Type>::clear()
v_ = NULL;
}
- nRows_ = 0;
+ mRows_ = 0;
nCols_ = 0;
}
template<class Form, class Type>
-void Foam::Matrix<Form, Type>::transfer(Matrix<Form, Type>& a)
+void Foam::Matrix<Form, Type>::transfer(Matrix<Form, Type>& M)
{
clear();
- nRows_ = a.nRows_;
- a.nRows_ = 0;
+ mRows_ = M.mRows_;
+ M.mRows_ = 0;
- nCols_ = a.nCols_;
- a.nCols_ = 0;
+ nCols_ = M.nCols_;
+ M.nCols_ = 0;
- v_ = a.v_;
- a.v_ = NULL;
+ v_ = M.v_;
+ M.v_ = NULL;
}
@@ -198,7 +245,7 @@ void Foam::Matrix<Form, Type>::setSize(const label m, const label n)
{
mType newMatrix(m, n, Zero);
- label minM = min(m, nRows_);
+ label minM = min(m, mRows_);
label minN = min(n, nCols_);
for (label i=0; i<minM; i++)
@@ -231,65 +278,23 @@ Form Foam::Matrix<Form, Type>::T() const
}
-template<class Form, class Type>
-Foam::Matrix<Form, Type>::ConstBlock::operator Field<cmptType>() const
-{
- if (nCols_ != 1)
- {
- FatalErrorInFunction
- << "Number of columns " << nCols_ << " != 1"
- << abort(FatalError);
- }
-
- Field<cmptType> f(nRows_);
-
- forAll(f, i)
- {
- f[i] = operator()(i, 0);
- }
-
- return f;
-}
-
-
-template<class Form, class Type>
-Foam::Matrix<Form, Type>::Block::operator Field<cmptType>() const
-{
- if (nCols_ != 1)
- {
- FatalErrorInFunction
- << "Number of columns " << nCols_ << " != 1"
- << abort(FatalError);
- }
-
- Field<cmptType> f(nRows_);
-
- forAll(f, i)
- {
- f[i] = operator()(i, 0);
- }
-
- return f;
-}
-
-
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
template<class Form, class Type>
-void Foam::Matrix<Form, Type>::operator=(const Matrix<Form, Type>& a)
+void Foam::Matrix<Form, Type>::operator=(const Matrix<Form, Type>& M)
{
- if (this == &a)
+ if (this == &M)
{
FatalErrorInFunction
<< "Attempted assignment to self"
<< abort(FatalError);
}
- if (nRows_ != a.nRows_ || nCols_ != a.nCols_)
+ if (mRows_ != M.mRows_ || nCols_ != M.nCols_)
{
clear();
- nRows_ = a.nRows_;
- nCols_ = a.nCols_;
+ mRows_ = M.mRows_;
+ nCols_ = M.nCols_;
allocate();
}
@@ -298,106 +303,55 @@ void Foam::Matrix<Form, Type>::operator=(const Matrix<Form, Type>& a)
const label mn = size();
for (label i=0; i<mn; i++)
{
- v_[i] = a.v_[i];
+ v_[i] = M.v_[i];
}
}
}
template<class Form, class Type>
-void Foam::Matrix<Form, Type>::operator=(const mType::Block& block)
+template<class MatrixType>
+void Foam::Matrix<Form, Type>::operator=
+(
+ const ConstMatrixBlock<MatrixType>& Mb
+)
{
- for (label i=0; i<nRows_; i++)
+ for (label i=0; i<mRows_; i++)
{
for (label j=0; j<nCols_; j++)
{
- (*this)(i,j) = block(i,j);
+ (*this)(i,j) = Mb(i,j);
}
}
}
template<class Form, class Type>
-void Foam::Matrix<Form, Type>::Block::operator=(const Block& block)
-{
- if (this != &block)
- {
- if (nRows_ != block.m() || nCols_ != block.n())
- {
- FatalErrorInFunction
- << "Attempt to assign blocks of different sizes: "
- << nRows_ << "x" << nCols_ << " != "
- << block.m() << "x" << block.n()
- << abort(FatalError);
- }
-
- for (label i=0; i<nRows_; i++)
- {
- for (label j=0; j<nCols_; j++)
- {
- (*this)(i, j) = block(i, j);
- }
- }
- }
-}
-
-
-template<class Form, class Type>
-void Foam::Matrix<Form, Type>::Block::operator=(const ConstBlock& block)
-{
- if (this != &block)
- {
- if (nRows_ != block.m() || nCols_ != block.n())
- {
- FatalErrorInFunction
- << "Attempt to assign blocks of different sizes: "
- << nRows_ << "x" << nCols_ << " != "
- << block.m() << "x" << block.n()
- << abort(FatalError);
- }
-
- for (label i=0; i<nRows_; i++)
- {
- for (label j=0; j<nCols_; j++)
- {
- (*this)(i, j) = block(i, j);
- }
- }
- }
-}
-
-
-template<class Form, class Type>
-void Foam::Matrix<Form, Type>::Block::operator=(const mType& block)
+template<class MatrixType>
+void Foam::Matrix<Form, Type>::operator=
+(
+ const MatrixBlock<MatrixType>& Mb
+)
{
- if (nRows_ != block.m() || nCols_ != block.n())
- {
- FatalErrorInFunction
- << "Attempt to assign blocks of different sizes: "
- << nRows_ << "x" << nCols_ << " != "
- << block.m() << "x" << block.n()
- << abort(FatalError);
- }
-
- for (label i=0; i<nRows_; i++)
+ for (label i=0; i<mRows_; i++)
{
for (label j=0; j<nCols_; j++)
{
- (*this)(i, j) = block(i, j);
+ (*this)(i,j) = Mb(i,j);
}
}
}
template<class Form, class Type>
-void Foam::Matrix<Form, Type>::operator=(const Type& t)
+void Foam::Matrix<Form, Type>::operator=(const Type& s)
{
if (v_)
{
const label mn = size();
for (label i=0; i<mn; i++)
{
- v_[i] = t;
+ v_[i] = s;
}
}
}
@@ -417,155 +371,27 @@ void Foam::Matrix<Form, Type>::operator=(const zero)
}
-template<class Form, class Type>
-template
-<
- template<class, Foam::direction, Foam::direction> class MSBlock,
- class SubTensor,
- Foam::direction BRowStart,
- Foam::direction BColStart
->
-void Foam::Matrix<Form, Type>::Block::operator=
-(
- const MSBlock<SubTensor, BRowStart, BColStart>& block
-)
-{
- if (nRows_ != block.nRows || nCols_ != block.nCols)
- {
- FatalErrorInFunction
- << "Attempt to assign blocks of different sizes: "
- << nRows_ << "x" << nCols_ << " != "
- << block.nRows << "x" << block.nCols
- << abort(FatalError);
- }
-
- for (direction i=0; i<nRows_; ++i)
- {
- for (direction j=0; j<nCols_; ++j)
- {
- operator()(i, j) = block(i, j);
- }
- }
-}
-
-
-template<class Form, class Type>
-template
-<
- template<class, Foam::direction> class VSBlock,
- class SubVector,
- Foam::direction BStart
->
-void Foam::Matrix<Form, Type>::Block::operator=
-(
- const VSBlock<SubVector, BStart>& block
-)
-{
- if (nRows_ != block.nComponents || nCols_ != 1)
- {
- FatalErrorInFunction
- << "Attempt to assign blocks of different sizes: "
- << nRows_ << "x" << nCols_ << " != "
- << block.nComponents << "x" << 1
- << abort(FatalError);
- }
-
- for (direction i=0; i<nRows_; ++i)
- {
- operator()(i, 0) = block[i];
- }
-}
-
-
-template<class Form, class Type>
-template<class MSForm, Foam::direction Nrows, Foam::direction Ncols>
-void Foam::Matrix<Form, Type>::Block::operator=
-(
- const MatrixSpace<MSForm, cmptType, Nrows, Ncols>& ms
-)
-{
- if (nRows_ != Nrows || nCols_ != Ncols)
- {
- FatalErrorInFunction
- << "Attempt to assign blocks of different sizes: "
- << nRows_ << "x" << nCols_ << " != "
- << Nrows << "x" << Ncols
- << abort(FatalError);
- }
-
- for (label i=0; i<nRows_; i++)
- {
- for (label j=0; j<nCols_; j++)
- {
- (*this)(i, j) = ms(i, j);
- }
- }
-}
-
-
-template<class Form, class Type>
-template<class VSForm, Foam::direction Ncmpts>
-void Foam::Matrix<Form, Type>::Block::operator=
-(
- const VectorSpace<VSForm, cmptType, Ncmpts>& ms
-)
-{
- if (nRows_ != Ncmpts || nCols_ != 1)
- {
- FatalErrorInFunction
- << "Attempt to assign blocks of different sizes: "
- << nRows_ << "x" << nCols_ << " != "
- << Ncmpts << "x" << 1
- << abort(FatalError);
- }
-
- for (direction i=0; i<Ncmpts; ++i)
- {
- operator()(i, 0) = ms[i];
- }
-}
-
-
-template<class Form, class Type>
-void Foam::Matrix<Form, Type>::Block::operator=(const Field<cmptType>& f)
-{
- if (nRows_ != f.size() || nCols_ != 1)
- {
- FatalErrorInFunction
- << "Error: cannot assign blocks of different size (left is "
- << nRows_ << "x" << nCols_ << " != "
- << f.size() << "x" << 1
- << abort(FatalError);
- }
-
- forAll(f, i)
- {
- operator()(i, 0) = f[i];
- }
-}
-
-
// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
template<class Form, class Type>
-const Type& Foam::max(const Matrix<Form, Type>& a)
+const Type& Foam::max(const Matrix<Form, Type>& M)
{
- const label mn = a.size();
+ const label mn = M.size();
if (mn)
{
label curMaxI = 0;
- const Type* v = a.v();
+ const Type* Mv = M.v();
for (label i=1; i<mn; i++)
{
- if (v[i] > v[curMaxI])
+ if (Mv[i] > Mv[curMaxI])
{
curMaxI = i;
}
}
- return v[curMaxI];
+ return Mv[curMaxI];
}
else
{
@@ -574,30 +400,30 @@ const Type& Foam::max(const Matrix<Form, Type>& a)
<< abort(FatalError);
// Return in error to keep compiler happy
- return a(0, 0);
+ return M(0, 0);
}
}
template<class Form, class Type>
-const Type& Foam::min(const Matrix<Form, Type>& a)
+const Type& Foam::min(const Matrix<Form, Type>& M)
{
- const label mn = a.size();
+ const label mn = M.size();
if (mn)
{
label curMinI = 0;
- const Type* v = a.v();
+ const Type* Mv = M.v();
for (label i=1; i<mn; i++)
{
- if (v[i] < v[curMinI])
+ if (Mv[i] < Mv[curMinI])
{
curMinI = i;
}
}
- return v[curMinI];
+ return Mv[curMinI];
}
else
{
@@ -606,7 +432,7 @@ const Type& Foam::min(const Matrix<Form, Type>& a)
<< abort(FatalError);
// Return in error to keep compiler happy
- return a(0, 0);
+ return M(0, 0);
}
}
@@ -614,187 +440,197 @@ const Type& Foam::min(const Matrix<Form, Type>& a)
// * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * //
template<class Form, class Type>
-Form Foam::operator-(const Matrix<Form, Type>& a)
+Form Foam::operator-(const Matrix<Form, Type>& M)
{
- Form na(a.m(), a.n());
+ Form nM(M.m(), M.n());
- if (a.m() && a.n())
+ if (M.m() && M.n())
{
- Type* nav = na.v();
- const Type* av = a.v();
+ Type* nMv = nM.v();
+ const Type* Mv = M.v();
- const label mn = a.size();
+ const label mn = M.size();
for (label i=0; i<mn; i++)
{
- nav[i] = -av[i];
+ nMv[i] = -Mv[i];
}
}
- return na;
+ return nM;
}
template<class Form, class Type>
-Form Foam::operator+(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
+Form Foam::operator+(const Matrix<Form, Type>& A, const Matrix<Form, Type>& B)
{
- if (a.m() != b.m())
+ if (A.m() != B.m())
{
FatalErrorInFunction
- << "Attempt to add matrices with different number of rows: "
- << a.m() << ", " << b.m()
+ << "Attempt to add matrices with different numbers of rows: "
+ << A.m() << ", " << B.m()
<< abort(FatalError);
}
- if (a.n() != b.n())
+ if (A.n() != B.n())
{
FatalErrorInFunction
- << "Attempt to add matrices with different number of columns: "
- << a.n() << ", " << b.n()
+ << "Attempt to add matrices with different numbers of columns: "
+ << A.n() << ", " << B.n()
<< abort(FatalError);
}
- Form ab(a.m(), a.n());
+ Form AB(A.m(), A.n());
- Type* abv = ab.v();
- const Type* av = a.v();
- const Type* bv = b.v();
+ Type* ABv = AB.v();
+ const Type* Av = A.v();
+ const Type* Bv = B.v();
- const label mn = a.size();
+ const label mn = A.size();
for (label i=0; i<mn; i++)
{
- abv[i] = av[i] + bv[i];
+ ABv[i] = Av[i] + Bv[i];
}
- return ab;
+ return AB;
}
template<class Form, class Type>
-Form Foam::operator-(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
+Form Foam::operator-(const Matrix<Form, Type>& A, const Matrix<Form, Type>& B)
{
- if (a.m() != b.m())
+ if (A.m() != B.m())
{
FatalErrorInFunction
- << "Attempt to add matrices with different number of rows: "
- << a.m() << ", " << b.m()
+ << "Attempt to add matrices with different numbers of rows: "
+ << A.m() << ", " << B.m()
<< abort(FatalError);
}
- if (a.n() != b.n())
+ if (A.n() != B.n())
{
FatalErrorInFunction
- << "Attempt to add matrices with different number of columns: "
- << a.n() << ", " << b.n()
+ << "Attempt to add matrices with different numbers of columns: "
+ << A.n() << ", " << B.n()
<< abort(FatalError);
}
- Form ab(a.m(), a.n());
+ Form AB(A.m(), A.n());
- Type* abv = ab.v();
- const Type* av = a.v();
- const Type* bv = b.v();
+ Type* ABv = AB.v();
+ const Type* Av = A.v();
+ const Type* Bv = B.v();
- const label mn = a.size();
+ const label mn = A.size();
for (label i=0; i<mn; i++)
{
- abv[i] = av[i] - bv[i];
+ ABv[i] = Av[i] - Bv[i];
}
- return ab;
+ return AB;
}
template<class Form, class Type>
-Form Foam::operator*(const scalar s, const Matrix<Form, Type>& a)
+Form Foam::operator*(const scalar s, const Matrix<Form, Type>& M)
{
- Form sa(a.m(), a.n());
+ Form sM(M.m(), M.n());
- if (a.m() && a.n())
+ if (M.m() && M.n())
{
- Type* sav = sa.v();
- const Type* av = a.v();
+ Type* sMv = sM.v();
+ const Type* Mv = M.v();
- const label mn = a.size();
+ const label mn = M.size();
for (label i=0; i<mn; i++)
{
- sav[i] = s*av[i];
+ sMv[i] = s*Mv[i];
}
}
- return sa;
+ return sM;
}
template<class Form, class Type>
-Form Foam::operator*(const Matrix<Form, Type>& a, const scalar s)
+Form Foam::operator*(const Matrix<Form, Type>& M, const scalar s)
{
- Form sa(a.m(), a.n());
+ Form sM(M.m(), M.n());
- if (a.m() && a.n())
+ if (M.m() && M.n())
{
- Type* sav = sa.v();
- const Type* av = a.v();
+ Type* sMv = sM.v();
+ const Type* Mv = M.v();
- const label mn = a.size();
+ const label mn = M.size();
for (label i=0; i<mn; i++)
{
- sav[i] = av[i]*s;
+ sMv[i] = Mv[i]*s;
}
}
- return sa;
+ return sM;
}
template<class Form, class Type>
-Form Foam::operator/(const Matrix<Form, Type>& a, const scalar s)
+Form Foam::operator/(const Matrix<Form, Type>& M, const scalar s)
{
- Form sa(a.m(), a.n());
+ Form sM(M.m(), M.n());
- if (a.m() && a.n())
+ if (M.m() && M.n())
{
- Type* sav = sa.v();
- const Type* av = a.v();
+ Type* sMv = sM.v();
+ const Type* Mv = M.v();
- const label mn = a.size();
+ const label mn = M.size();
for (label i=0; i<mn; i++)
{
- sav[i] = av[i]/s;
+ sMv[i] = Mv[i]/s;
}
}
- return sa;
+ return sM;
}
-template<class Form, class Type>
-Form Foam::operator*(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
+template<class Form1, class Form2, class Type>
+typename Foam::typeOfInnerProduct<Type, Form1, Form2>::type
+Foam::operator*
+(
+ const Matrix<Form1, Type>& A,
+ const Matrix<Form2, Type>& B
+)
{
- if (a.n() != b.m())
+ if (A.n() != B.m())
{
FatalErrorInFunction
<< "Attempt to multiply incompatible matrices:" << nl
- << "Matrix A : " << a.m() << " rows, " << a.n() << " columns" << nl
- << "Matrix B : " << b.m() << " rows, " << b.n() << " columns" << nl
+ << "Matrix A : " << A.m() << " x " << A.n() << nl
+ << "Matrix B : " << B.m() << " x " << B.n() << nl
<< "In order to multiply matrices, columns of A must equal "
<< "rows of B"
<< abort(FatalError);
}
- Form ab(a.m(), b.n(), scalar(0));
+ typename typeOfInnerProduct<Type, Form1, Form2>::type AB
+ (
+ A.m(),
+ B.n(),
+ Zero
+ );
- for (label i=0; i<ab.m(); i++)
+ for (label i=0; i<AB.m(); i++)
{
- for (label j=0; j<ab.n(); j++)
+ for (label j=0; j<AB.n(); j++)
{
- for (label k=0; k<b.m(); k++)
+ for (label k=0; k<B.m(); k++)
{
- ab(i, j) += a(i, k)*b(k, j);
+ AB(i, j) += A(i, k)*B(k, j);
}
}
}
- return ab;
+ return AB;
}
@@ -809,7 +645,7 @@ inline Foam::tmp<Foam::Field<Type>> Foam::operator*
{
FatalErrorInFunction
<< "Attempt to multiply incompatible matrix and field:" << nl
- << "Matrix : " << M.m() << " rows, " << M.n() << " columns" << nl
+ << "Matrix : " << M.m() << " x " << M.n() << nl
<< "Field : " << f.size() << " rows" << nl
<< "In order to multiply a matrix M and field f, "
"columns of M must equal rows of f"
@@ -819,9 +655,9 @@ inline Foam::tmp<Foam::Field<Type>> Foam::operator*
tmp<Field<Type>> tMf(new Field<Type>(f.size(), Zero));
Field<Type>& Mf = tMf.ref();
- for (label i=0; i<M.m(); ++i)
+ for (label i=0; i<M.m(); i++)
{
- for (label j=0; j<M.n(); ++j)
+ for (label j=0; j<M.n(); j++)
{
Mf[i] += M(i, j)*f[j];
}
diff --git a/src/OpenFOAM/matrices/Matrix/Matrix.H b/src/OpenFOAM/matrices/Matrix/Matrix.H
index 5f905b465bc..9385d42d873 100644
--- a/src/OpenFOAM/matrices/Matrix/Matrix.H
+++ b/src/OpenFOAM/matrices/Matrix/Matrix.H
@@ -25,8 +25,7 @@ Class
Foam::Matrix
Description
- A templated 2D matrix of objects of \<T\>, where the n x m matrix
- dimensions are known and used for subscript bounds checking, etc.
+ A templated (m x n) matrix of objects of \<T\>.
SourceFiles
Matrix.C
@@ -42,7 +41,6 @@ SourceFiles
#include "label.H"
#include "uLabel.H"
#include "Field.H"
-#include "MatrixSpace.H"
#include "zero.H"
#include "autoPtr.H"
@@ -67,6 +65,12 @@ template<class Form, class Type> Ostream& operator<<
const Matrix<Form, Type>&
);
+template<class MatrixType>
+class ConstMatrixBlock;
+
+template<class MatrixType>
+class MatrixBlock;
+
/*---------------------------------------------------------------------------*\
Class Matrix Declaration
@@ -78,7 +82,7 @@ class Matrix
// Private data
//- Number of rows and columns in Matrix.
- label nRows_, nCols_;
+ label mRows_, nCols_;
//- Row pointers
Type* __restrict__ v_;
@@ -102,157 +106,6 @@ public:
inline static const mType& null();
- // Sub-Block Classes
-
- class ConstBlock
- {
- // Private data
-
- //- Const reference to the parent matrix
- const mType& matrix_;
-
- // Block size
- const label nRows_;
- const label nCols_;
-
- // Block location in parent matrix
- const label rowStart_;
- const label colStart_;
-
- public:
-
- typedef typename mType::cmptType cmptType;
-
- // Constructors
-
- //- Construct block for matrix, size and location
- inline ConstBlock
- (
- const mType& matrix,
- const label m,
- const label n,
- const label mStart,
- const label nStart
- );
-
-
- // Member Functions
-
- //- Return the number of rows in the block
- inline label m() const;
-
- //- Return the number of columns in the block
- inline label n() const;
-
- //- (i, j) const element access operator
- inline const Type& operator()
- (
- const label i,
- const label j
- ) const;
-
- //- Convert a column of a matrix to a Field
- operator Field<Type>() const;
- };
-
-
- class Block
- {
- // Private data
-
- //- Reference to the parent matrix
- mType& matrix_;
-
- // Block size
- const label nRows_;
- const label nCols_;
-
- // Block location in parent matrix
- const label rowStart_;
- const label colStart_;
-
- public:
-
- typedef typename mType::cmptType cmptType;
-
- // Constructors
-
- //- Construct block for matrix, size and location
- inline Block
- (
- mType& matrix,
- const label m,
- const label n,
- const label mStart,
- const label nStart
- );
-
-
- // Member Functions
-
- //- Return the number of rows in the block
- inline label m() const;
-
- //- Return the number of columns in the block
- inline label n() const;
-
- //- (i, j) const element access operator
- inline const Type& operator()
- (
- const label i,
- const label j
- ) const;
-
- //- (i, j) element access operator
- inline Type& operator()(const label i, const label j);
-
- //- Convert a column of a matrix to a Field
- operator Field<Type>() const;
-
-
- // Member operators
-
- //- Assignment to a compatible matrix
- void operator=(const mType&);
-
- //- Assignment to a compatible block
- void operator=(const Block&);
-
- //- Assignment to a compatible const block
- void operator=(const ConstBlock&);
-
- //- Assignment to a compatible MatrixSpace
- template<class MSForm, direction Nrows, direction Ncols>
- void operator=(const MatrixSpace<MSForm, Type, Nrows, Ncols>&);
-
- //- Assignment to a compatible MatrixSpace block
- template
- <
- template<class, direction, direction> class Block,
- class SubTensor,
- direction BRowStart,
- direction BColStart
- >
- void operator=(const Block<SubTensor, BRowStart, BColStart>&);
-
- //- Assignment to a compatible VectorSpace (column-vector)
- template<class VSForm, direction Ncmpts>
- void operator=(const VectorSpace<VSForm, Type, Ncmpts>&);
-
- //- Assignment to a compatible VectorSpace (column-vector) block
- template
- <
- template<class, direction> class Block,
- class SubVector,
- direction BStart
- >
- void operator=(const Block<SubVector, BStart>&);
-
- //- Assignment to a Field (column-vector)
- void operator=(const Field<Type>&);
- };
-
-
// Constructors
//- Null constructor.
@@ -272,8 +125,17 @@ public:
//- Copy constructor.
Matrix(const mType&);
+ //- Copy constructor from matrix of a different form
+ template<class Form2>
+ explicit Matrix(const Matrix<Form2, Type>&);
+
+ //- Construct from a block of another matrix
+ template<class MatrixType>
+ Matrix(const ConstMatrixBlock<MatrixType>&);
+
//- Construct from a block of another matrix
- Matrix(const mType::Block&);
+ template<class MatrixType>
+ Matrix(const MatrixBlock<MatrixType>&);
//- Construct from Istream.
Matrix(Istream&);
@@ -308,7 +170,7 @@ public:
// Block access
- inline ConstBlock block
+ inline ConstMatrixBlock<mType> block
(
const label m,
const label n,
@@ -317,19 +179,19 @@ public:
) const;
template<class VectorSpace>
- inline ConstBlock block
+ inline ConstMatrixBlock<mType> block
(
const label mStart,
const label nStart
) const;
- inline ConstBlock col
+ inline ConstMatrixBlock<mType> col
(
const label m,
const label rowStart
) const;
- inline ConstBlock col
+ inline ConstMatrixBlock<mType> col
(
const label m,
const label mStart,
@@ -337,7 +199,7 @@ public:
) const;
- inline Block block
+ inline MatrixBlock<mType> block
(
const label m,
const label n,
@@ -346,11 +208,19 @@ public:
);
template<class VectorSpace>
- inline Block block(const label mStart, const label nStart);
+ inline MatrixBlock<mType> block
+ (
+ const label mStart,
+ const label nStart
+ );
- inline Block col(const label m, const label rowStart);
+ inline MatrixBlock<mType> col
+ (
+ const label m,
+ const label rowStart
+ );
- inline Block col
+ inline MatrixBlock<mType> col
(
const label m,
const label mStart,
@@ -402,7 +272,12 @@ public:
void operator=(const mType&);
//- Assignment to a block of another matrix
- void operator=(const mType::Block&);
+ template<class MatrixType>
+ void operator=(const ConstMatrixBlock<MatrixType>&);
+
+ //- Assignment to a block of another matrix
+ template<class MatrixType>
+ void operator=(const MatrixBlock<MatrixType>&);
//- Assignment of all entries to zero
void operator=(const zero);
@@ -475,11 +350,12 @@ Form operator/
const scalar
);
-template<class Form, class Type>
-Form operator*
+template<class Form1, class Form2, class Type>
+typename typeOfInnerProduct<Type, Form1, Form2>::type
+operator*
(
- const Matrix<Form, Type>&,
- const Matrix<Form, Type>&
+ const Matrix<Form1, Type>& a,
+ const Matrix<Form2, Type>& b
);
template<class Form, class Type>
diff --git a/src/OpenFOAM/matrices/Matrix/MatrixI.H b/src/OpenFOAM/matrices/Matrix/MatrixI.H
index 73992f3285d..73992aac77b 100644
--- a/src/OpenFOAM/matrices/Matrix/MatrixI.H
+++ b/src/OpenFOAM/matrices/Matrix/MatrixI.H
@@ -23,12 +23,14 @@ License
\*---------------------------------------------------------------------------*/
+#include "MatrixBlock.H"
+
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Form, class Type>
inline Foam::Matrix<Form, Type>::Matrix()
:
- nRows_(0),
+ mRows_(0),
nCols_(0),
v_(NULL)
{}
@@ -42,68 +44,6 @@ clone() const
}
-template<class Form, class Type>
-Foam::Matrix<Form, Type>::ConstBlock::ConstBlock
-(
- const mType& matrix,
- const label m,
- const label n,
- const label mStart,
- const label nStart
-)
-:
- matrix_(matrix),
- nRows_(m),
- nCols_(n),
- rowStart_(mStart),
- colStart_(nStart)
-{
- #ifdef FULLDEBUG
- if
- (
- rowStart_ + nRows_ > matrix.m()
- || colStart_ + nCols_ > matrix.n()
- )
- {
- FatalErrorInFunction
- << "Block addresses outside matrix"
- << abort(FatalError);
- }
- #endif
-}
-
-
-template<class Form, class Type>
-Foam::Matrix<Form, Type>::Block::Block
-(
- mType& matrix,
- const label m,
- const label n,
- const label mStart,
- const label nStart
-)
-:
- matrix_(matrix),
- nRows_(m),
- nCols_(n),
- rowStart_(mStart),
- colStart_(nStart)
-{
- #ifdef FULLDEBUG
- if
- (
- rowStart_ + nRows_ > matrix.m()
- || colStart_ + nCols_ > matrix.n()
- )
- {
- FatalErrorInFunction
- << "Block addresses outside matrix"
- << abort(FatalError);
- }
- #endif
-}
-
-
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Form, class Type>
@@ -116,7 +56,7 @@ inline const Foam::Matrix<Form, Type>& Foam::Matrix<Form, Type>::null()
template<class Form, class Type>
inline Foam::label Foam::Matrix<Form, Type>::m() const
{
- return nRows_;
+ return mRows_;
}
@@ -130,7 +70,7 @@ inline Foam::label Foam::Matrix<Form, Type>::n() const
template<class Form, class Type>
inline Foam::label Foam::Matrix<Form, Type>::size() const
{
- return nRows_*nCols_;
+ return mRows_*nCols_;
}
@@ -138,16 +78,16 @@ template<class Form, class Type>
inline void Foam::Matrix<Form, Type>::checki(const label i) const
{
#ifdef FULLDEBUG
- if (!nRows_ || !nCols_)
+ if (!mRows_ || !nCols_)
{
FatalErrorInFunction
<< "Attempt to access element from empty matrix"
<< abort(FatalError);
}
- else if (i<0 || i>=nRows_)
+ else if (i<0 || i>=mRows_)
{
FatalErrorInFunction
- << "Index " << i << " out of range 0 ... " << nRows_-1
+ << "Index " << i << " out of range 0 ... " << mRows_-1
<< abort(FatalError);
}
#endif
@@ -158,7 +98,7 @@ template<class Form, class Type>
inline void Foam::Matrix<Form, Type>::checkj(const label j) const
{
#ifdef FULLDEBUG
- if (!nRows_ || !nCols_)
+ if (!mRows_ || !nCols_)
{
FatalErrorInFunction
<< "Attempt to access element from empty matrix"
@@ -174,34 +114,6 @@ inline void Foam::Matrix<Form, Type>::checkj(const label j) const
}
-template<class Form, class Type>
-inline Foam::label Foam::Matrix<Form, Type>::ConstBlock::m() const
-{
- return nRows_;
-}
-
-
-template<class Form, class Type>
-inline Foam::label Foam::Matrix<Form, Type>::ConstBlock::n() const
-{
- return nCols_;
-}
-
-
-template<class Form, class Type>
-inline Foam::label Foam::Matrix<Form, Type>::Block::m() const
-{
- return nRows_;
-}
-
-
-template<class Form, class Type>
-inline Foam::label Foam::Matrix<Form, Type>::Block::n() const
-{
- return nCols_;
-}
-
-
template<class Form, class Type>
inline const Type* Foam::Matrix<Form, Type>::v() const
{
@@ -217,7 +129,7 @@ inline Type* Foam::Matrix<Form, Type>::v()
template<class Form, class Type>
-inline typename Foam::Matrix<Form, Type>::ConstBlock
+inline Foam::ConstMatrixBlock<Foam::Matrix<Form, Type>>
Foam::Matrix<Form, Type>::block
(
const label m,
@@ -226,7 +138,7 @@ Foam::Matrix<Form, Type>::block
const label nStart
) const
{
- return ConstBlock
+ return ConstMatrixBlock<mType>
(
*this,
m,
@@ -239,17 +151,17 @@ Foam::Matrix<Form, Type>::block
template<class Form, class Type>
template<class VectorSpace>
-inline typename Foam::Matrix<Form, Type>::ConstBlock
+inline Foam::ConstMatrixBlock<Foam::Matrix<Form, Type>>
Foam::Matrix<Form, Type>::block
(
const label mStart,
const label nStart
) const
{
- return ConstBlock
+ return ConstMatrixBlock<mType>
(
*this,
- VectorSpace::nRows,
+ VectorSpace::mRows,
VectorSpace::nCols,
mStart,
nStart
@@ -258,14 +170,14 @@ Foam::Matrix<Form, Type>::block
template<class Form, class Type>
-inline typename Foam::Matrix<Form, Type>::ConstBlock
+inline Foam::ConstMatrixBlock<Foam::Matrix<Form, Type>>
Foam::Matrix<Form, Type>::col
(
const label m,
const label mStart
) const
{
- return ConstBlock
+ return ConstMatrixBlock<mType>
(
*this,
m,
@@ -277,7 +189,7 @@ Foam::Matrix<Form, Type>::col
template<class Form, class Type>
-inline typename Foam::Matrix<Form, Type>::ConstBlock
+inline Foam::ConstMatrixBlock<Foam::Matrix<Form, Type>>
Foam::Matrix<Form, Type>::col
(
const label m,
@@ -285,7 +197,7 @@ Foam::Matrix<Form, Type>::col
const label nStart
) const
{
- return ConstBlock
+ return ConstMatrixBlock<mType>
(
*this,
m,
@@ -297,7 +209,7 @@ Foam::Matrix<Form, Type>::col
template<class Form, class Type>
-inline typename Foam::Matrix<Form, Type>::Block
+inline Foam::MatrixBlock<Foam::Matrix<Form, Type>>
Foam::Matrix<Form, Type>::block
(
const label m,
@@ -306,7 +218,7 @@ Foam::Matrix<Form, Type>::block
const label nStart
)
{
- return Block
+ return MatrixBlock<mType>
(
*this,
m,
@@ -319,13 +231,13 @@ Foam::Matrix<Form, Type>::block
template<class Form, class Type>
template<class VectorSpace>
-inline typename Foam::Matrix<Form, Type>::Block
+inline Foam::MatrixBlock<Foam::Matrix<Form, Type>>
Foam::Matrix<Form, Type>::block(const label mStart, const label nStart)
{
- return Block
+ return MatrixBlock<mType>
(
*this,
- VectorSpace::nRows,
+ VectorSpace::mRows,
VectorSpace::nCols,
mStart,
nStart
@@ -334,10 +246,10 @@ Foam::Matrix<Form, Type>::block(const label mStart, const label nStart)
template<class Form, class Type>
-inline typename Foam::Matrix<Form, Type>::Block
+inline Foam::MatrixBlock<Foam::Matrix<Form, Type>>
Foam::Matrix<Form, Type>::col(const label m, const label mStart)
{
- return Block
+ return MatrixBlock<mType>
(
*this,
m,
@@ -349,7 +261,7 @@ Foam::Matrix<Form, Type>::col(const label m, const label mStart)
template<class Form, class Type>
-inline typename Foam::Matrix<Form, Type>::Block
+inline Foam::MatrixBlock<Foam::Matrix<Form, Type>>
Foam::Matrix<Form, Type>::col
(
const label m,
@@ -357,7 +269,7 @@ Foam::Matrix<Form, Type>::col
const label nStart
)
{
- return Block
+ return MatrixBlock<mType>
(
*this,
m,
@@ -370,84 +282,6 @@ Foam::Matrix<Form, Type>::col
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
-template<class Form, class Type>
-inline const Type& Foam::Matrix<Form, Type>::ConstBlock::operator()
-(
- const label i,
- const label j
-) const
-{
- #ifdef FULLDEBUG
- if (i<0 || i>=nRows_)
- {
- FatalErrorInFunction
- << "Index " << i << " out of range 0 ... " << nRows_-1
- << abort(FatalError);
- }
- if (j<0 || j>=nCols_)
- {
- FatalErrorInFunction
- << "Index " << j << " out of range 0 ... " << nCols_-1
- << abort(FatalError);
- }
- #endif
-
- return matrix_(i + rowStart_, j + colStart_);
-}
-
-
-template<class Form, class Type>
-inline const Type& Foam::Matrix<Form, Type>::Block::operator()
-(
- const label i,
- const label j
-) const
-{
- #ifdef FULLDEBUG
- if (i<0 || i>=nRows_)
- {
- FatalErrorInFunction
- << "Index " << i << " out of range 0 ... " << nRows_-1
- << abort(FatalError);
- }
- if (j<0 || j>=nCols_)
- {
- FatalErrorInFunction
- << "Index " << j << " out of range 0 ... " << nCols_-1
- << abort(FatalError);
- }
- #endif
-
- return matrix_(i + rowStart_, j + colStart_);
-}
-
-
-template<class Form, class Type>
-inline Type& Foam::Matrix<Form, Type>::Block::operator()
-(
- const label i,
- const label j
-)
-{
- #ifdef FULLDEBUG
- if (i<0 || i>=nRows_)
- {
- FatalErrorInFunction
- << "Index " << i << " out of range 0 ... " << nRows_-1
- << abort(FatalError);
- }
- if (j<0 || j>=nCols_)
- {
- FatalErrorInFunction
- << "Index " << j << " out of range 0 ... " << nCols_-1
- << abort(FatalError);
- }
- #endif
-
- return matrix_(i + rowStart_, j + colStart_);
-}
-
-
template<class Form, class Type>
inline const Type& Foam::Matrix<Form, Type>::operator()
(
diff --git a/src/OpenFOAM/matrices/Matrix/MatrixIO.C b/src/OpenFOAM/matrices/Matrix/MatrixIO.C
index 9be178c12f3..51943ada59f 100644
--- a/src/OpenFOAM/matrices/Matrix/MatrixIO.C
+++ b/src/OpenFOAM/matrices/Matrix/MatrixIO.C
@@ -34,7 +34,7 @@ License
template<class Form, class Type>
Foam::Matrix<Form, Type>::Matrix(Istream& is)
:
- nRows_(0),
+ mRows_(0),
nCols_(0),
v_(NULL)
{
@@ -59,10 +59,10 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<Form, Type>& M)
if (firstToken.isLabel())
{
- M.nRows_ = firstToken.labelToken();
+ M.mRows_ = firstToken.labelToken();
M.nCols_ = readLabel(is);
- label mn = M.nRows_*M.nCols_;
+ label mn = M.mRows_*M.nCols_;
// Read list contents depending on data format
if (is.format() == IOstream::ASCII || !contiguous<Type>())
@@ -151,7 +151,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<Form, Type>& M)
template<class Form, class Type>
Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<Form, Type>& M)
{
- label mn = M.nRows_*M.nCols_;
+ label mn = M.mRows_*M.nCols_;
os << M.m() << token::SPACE << M.n();
diff --git a/src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.C b/src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.C
new file mode 100644
index 00000000000..3aecbe8377d
--- /dev/null
+++ b/src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.C
@@ -0,0 +1,342 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2016 OpenFOAM Foundation
+ \\/ 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "MatrixBlock.H"
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+template<class MatrixType>
+Foam::ConstMatrixBlock<MatrixType>::operator Field<cmptType>() const
+{
+ if (nCols_ != 1)
+ {
+ FatalErrorInFunction
+ << "Number of columns " << nCols_ << " != 1"
+ << abort(FatalError);
+ }
+
+ Field<cmptType> f(mRows_);
+
+ forAll(f, i)
+ {
+ f[i] = operator()(i, 0);
+ }
+
+ return f;
+}
+
+
+template<class MatrixType>
+Foam::MatrixBlock<MatrixType>::operator Field<cmptType>() const
+{
+ if (nCols_ != 1)
+ {
+ FatalErrorInFunction
+ << "Number of columns " << nCols_ << " != 1"
+ << abort(FatalError);
+ }
+
+ Field<cmptType> f(mRows_);
+
+ forAll(f, i)
+ {
+ f[i] = operator()(i, 0);
+ }
+
+ return f;
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
+
+template<class MatrixType>
+template<class Form>
+void Foam::MatrixBlock<MatrixType>::operator=
+(
+ const Matrix<Form, cmptType>& Mb
+)
+{
+ if (mRows_ != Mb.m() || nCols_ != Mb.n())
+ {
+ FatalErrorInFunction
+ << "Attempt to assign blocks of different sizes: "
+ << mRows_ << "x" << nCols_ << " != "
+ << Mb.m() << "x" << Mb.n()
+ << abort(FatalError);
+ }
+
+ for (label i=0; i<mRows_; i++)
+ {
+ for (label j=0; j<nCols_; j++)
+ {
+ (*this)(i, j) = Mb(i, j);
+ }
+ }
+}
+
+
+template<class MatrixType>
+void Foam::MatrixBlock<MatrixType>::operator=
+(
+ const ConstMatrixBlock<MatrixType>& Mb
+)
+{
+ if (this != &Mb)
+ {
+ if (mRows_ != Mb.m() || nCols_ != Mb.n())
+ {
+ FatalErrorInFunction
+ << "Attempt to assign blocks of different sizes: "
+ << mRows_ << "x" << nCols_ << " != "
+ << Mb.m() << "x" << Mb.n()
+ << abort(FatalError);
+ }
+
+ for (label i=0; i<mRows_; i++)
+ {
+ for (label j=0; j<nCols_; j++)
+ {
+ (*this)(i, j) = Mb(i, j);
+ }
+ }
+ }
+}
+
+
+template<class MatrixType>
+void Foam::MatrixBlock<MatrixType>::operator=
+(
+ const MatrixBlock<MatrixType>& Mb
+)
+{
+ if (this != &Mb)
+ {
+ if (mRows_ != Mb.m() || nCols_ != Mb.n())
+ {
+ FatalErrorInFunction
+ << "Attempt to assign blocks of different sizes: "
+ << mRows_ << "x" << nCols_ << " != "
+ << Mb.m() << "x" << Mb.n()
+ << abort(FatalError);
+ }
+
+ for (label i=0; i<mRows_; i++)
+ {
+ for (label j=0; j<nCols_; j++)
+ {
+ (*this)(i, j) = Mb(i, j);
+ }
+ }
+ }
+}
+
+
+template<class MatrixType>
+template<class MatrixType2>
+void Foam::MatrixBlock<MatrixType>::operator=
+(
+ const ConstMatrixBlock<MatrixType2>& Mb
+)
+{
+ if (this != &Mb)
+ {
+ if (mRows_ != Mb.m() || nCols_ != Mb.n())
+ {
+ FatalErrorInFunction
+ << "Attempt to assign blocks of different sizes: "
+ << mRows_ << "x" << nCols_ << " != "
+ << Mb.m() << "x" << Mb.n()
+ << abort(FatalError);
+ }
+
+ for (label i=0; i<mRows_; i++)
+ {
+ for (label j=0; j<nCols_; j++)
+ {
+ (*this)(i, j) = Mb(i, j);
+ }
+ }
+ }
+}
+
+
+template<class MatrixType>
+template<class MatrixType2>
+void Foam::MatrixBlock<MatrixType>::operator=
+(
+ const MatrixBlock<MatrixType2>& Mb
+)
+{
+ if (this != &Mb)
+ {
+ if (mRows_ != Mb.m() || nCols_ != Mb.n())
+ {
+ FatalErrorInFunction
+ << "Attempt to assign blocks of different sizes: "
+ << mRows_ << "x" << nCols_ << " != "
+ << Mb.m() << "x" << Mb.n()
+ << abort(FatalError);
+ }
+
+ for (label i=0; i<mRows_; i++)
+ {
+ for (label j=0; j<nCols_; j++)
+ {
+ (*this)(i, j) = Mb(i, j);
+ }
+ }
+ }
+}
+
+
+template<class MatrixType>
+template
+<
+ template<class, Foam::direction, Foam::direction> class MSBlock,
+ class SubTensor,
+ Foam::direction BRowStart,
+ Foam::direction BColStart
+>
+void Foam::MatrixBlock<MatrixType>::operator=
+(
+ const MSBlock<SubTensor, BRowStart, BColStart>& Mb
+)
+{
+ if (mRows_ != Mb.mRows || nCols_ != Mb.nCols)
+ {
+ FatalErrorInFunction
+ << "Attempt to assign blocks of different sizes: "
+ << mRows_ << "x" << nCols_ << " != "
+ << Mb.mRows << "x" << Mb.nCols
+ << abort(FatalError);
+ }
+
+ for (direction i=0; i<mRows_; ++i)
+ {
+ for (direction j=0; j<nCols_; ++j)
+ {
+ operator()(i, j) = Mb(i, j);
+ }
+ }
+}
+
+
+template<class MatrixType>
+template
+<
+ template<class, Foam::direction> class VSBlock,
+ class SubVector,
+ Foam::direction BStart
+>
+void Foam::MatrixBlock<MatrixType>::operator=
+(
+ const VSBlock<SubVector, BStart>& Mb
+)
+{
+ if (mRows_ != Mb.nComponents || nCols_ != 1)
+ {
+ FatalErrorInFunction
+ << "Attempt to assign blocks of different sizes: "
+ << mRows_ << "x" << nCols_ << " != "
+ << Mb.nComponents << "x" << 1
+ << abort(FatalError);
+ }
+
+ for (direction i=0; i<mRows_; ++i)
+ {
+ operator()(i, 0) = Mb[i];
+ }
+}
+
+
+template<class MatrixType>
+template<class MSForm, Foam::direction Nrows, Foam::direction Ncols>
+void Foam::MatrixBlock<MatrixType>::operator=
+(
+ const MatrixSpace<MSForm, cmptType, Nrows, Ncols>& ms
+)
+{
+ if (mRows_ != Nrows || nCols_ != Ncols)
+ {
+ FatalErrorInFunction
+ << "Attempt to assign blocks of different sizes: "
+ << mRows_ << "x" << nCols_ << " != "
+ << Nrows << "x" << Ncols
+ << abort(FatalError);
+ }
+
+ for (label i=0; i<mRows_; i++)
+ {
+ for (label j=0; j<nCols_; j++)
+ {
+ (*this)(i, j) = ms(i, j);
+ }
+ }
+}
+
+
+template<class MatrixType>
+template<class VSForm, Foam::direction Ncmpts>
+void Foam::MatrixBlock<MatrixType>::operator=
+(
+ const VectorSpace<VSForm, cmptType, Ncmpts>& ms
+)
+{
+ if (mRows_ != Ncmpts || nCols_ != 1)
+ {
+ FatalErrorInFunction
+ << "Attempt to assign blocks of different sizes: "
+ << mRows_ << "x" << nCols_ << " != "
+ << Ncmpts << "x" << 1
+ << abort(FatalError);
+ }
+
+ for (direction i=0; i<Ncmpts; ++i)
+ {
+ operator()(i, 0) = ms[i];
+ }
+}
+
+
+template<class MatrixType>
+void Foam::MatrixBlock<MatrixType>::operator=(const Field<cmptType>& f)
+{
+ if (mRows_ != f.size() || nCols_ != 1)
+ {
+ FatalErrorInFunction
+ << "Error: cannot assign blocks of different size (left is "
+ << mRows_ << "x" << nCols_ << " != "
+ << f.size() << "x" << 1
+ << abort(FatalError);
+ }
+
+ forAll(f, i)
+ {
+ operator()(i, 0) = f[i];
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.H b/src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.H
new file mode 100644
index 00000000000..67ece198eed
--- /dev/null
+++ b/src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.H
@@ -0,0 +1,240 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2016 OpenFOAM Foundation
+ \\/ 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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+ Foam::MatrixBlock
+
+Description
+ A templated block of an (m x n) matrix of type \<MatrixType\>.
+
+ Foam::ConstMatrixBlock: block of a const matrix
+ Foam::MatrixBlock: block of a non-const matrix
+
+ The block may be assigned to a block of another matrix or to a VectorSpace
+ or MatrixSpace e.g. \c tensor. Conversion of a column block to a \c
+ Field<T> is also provide.
+
+SourceFiles
+ MatrixBlock.C
+ MatrixBlockI.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef MatrixBlock_H
+#define MatrixBlock_H
+
+#include "Matrix.H"
+#include "MatrixSpace.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class ConstMatrixBlock Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class MatrixType>
+class ConstMatrixBlock
+{
+ // Private data
+
+ //- Const reference to the parent matrix
+ const MatrixType& matrix_;
+
+ // Block size
+ const label mRows_;
+ const label nCols_;
+
+ // Block location in parent matrix
+ const label rowStart_;
+ const label colStart_;
+
+public:
+
+ typedef typename MatrixType::cmptType cmptType;
+
+ // Constructors
+
+ //- Construct block for matrix, size and location
+ inline ConstMatrixBlock
+ (
+ const MatrixType& matrix,
+ const label m,
+ const label n,
+ const label mStart,
+ const label nStart
+ );
+
+
+ // Member Functions
+
+ //- Return the number of rows in the block
+ inline label m() const;
+
+ //- Return the number of columns in the block
+ inline label n() const;
+
+ //- (i, j) const element access operator
+ inline const cmptType& operator()
+ (
+ const label i,
+ const label j
+ ) const;
+
+ //- Convert a column of a matrix to a Field
+ operator Field<cmptType>() const;
+};
+
+
+/*---------------------------------------------------------------------------*\
+ Class MatrixBlock Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class MatrixType>
+class MatrixBlock
+{
+ // Private data
+
+ //- Reference to the parent matrix
+ MatrixType& matrix_;
+
+ // Block size
+ const label mRows_;
+ const label nCols_;
+
+ // Block location in parent matrix
+ const label rowStart_;
+ const label colStart_;
+
+public:
+
+ typedef typename MatrixType::cmptType cmptType;
+
+ // Constructors
+
+ //- Construct block for matrix, size and location
+ inline MatrixBlock
+ (
+ MatrixType& matrix,
+ const label m,
+ const label n,
+ const label mStart,
+ const label nStart
+ );
+
+
+ // Member Functions
+
+ //- Return the number of rows in the block
+ inline label m() const;
+
+ //- Return the number of columns in the block
+ inline label n() const;
+
+ //- (i, j) const element access operator
+ inline const cmptType& operator()
+ (
+ const label i,
+ const label j
+ ) const;
+
+ //- (i, j) element access operator
+ inline cmptType& operator()(const label i, const label j);
+
+ //- Convert a column of a matrix to a Field
+ operator Field<cmptType>() const;
+
+
+ // Member operators
+
+ //- Assignment to a compatible matrix
+ template<class Form>
+ void operator=(const Matrix<Form, cmptType>&);
+
+ //- Assignment to a compatible const block
+ void operator=(const ConstMatrixBlock<MatrixType>&);
+
+ //- Assignment to a compatible block
+ void operator=(const MatrixBlock<MatrixType>&);
+
+ //- Assignment to a compatible const block
+ template<class MatrixType2>
+ void operator=(const ConstMatrixBlock<MatrixType2>&);
+
+ //- Assignment to a compatible block
+ template<class MatrixType2>
+ void operator=(const MatrixBlock<MatrixType2>&);
+
+ //- Assignment to a compatible MatrixSpace
+ template<class MSForm, direction Nrows, direction Ncols>
+ void operator=(const MatrixSpace<MSForm, cmptType, Nrows, Ncols>&);
+
+ //- Assignment to a compatible MatrixSpace block
+ template
+ <
+ template<class, direction, direction> class Block,
+ class SubTensor,
+ direction BRowStart,
+ direction BColStart
+ >
+ void operator=(const Block<SubTensor, BRowStart, BColStart>&);
+
+ //- Assignment to a compatible VectorSpace (column-vector)
+ template<class VSForm, direction Ncmpts>
+ void operator=(const VectorSpace<VSForm, cmptType, Ncmpts>&);
+
+ //- Assignment to a compatible VectorSpace (column-vector) block
+ template
+ <
+ template<class, direction> class Block,
+ class SubVector,
+ direction BStart
+ >
+ void operator=(const Block<SubVector, BStart>&);
+
+ //- Assignment to a Field (column-vector)
+ void operator=(const Field<cmptType>&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "MatrixBlockI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+ #include "MatrixBlock.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/MatrixBlock/MatrixBlockI.H b/src/OpenFOAM/matrices/MatrixBlock/MatrixBlockI.H
new file mode 100644
index 00000000000..e764dea0320
--- /dev/null
+++ b/src/OpenFOAM/matrices/MatrixBlock/MatrixBlockI.H
@@ -0,0 +1,203 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2016 OpenFOAM Foundation
+ \\/ 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+template<class MatrixType>
+Foam::ConstMatrixBlock<MatrixType>::ConstMatrixBlock
+(
+ const MatrixType& matrix,
+ const label m,
+ const label n,
+ const label mStart,
+ const label nStart
+)
+:
+ matrix_(matrix),
+ mRows_(m),
+ nCols_(n),
+ rowStart_(mStart),
+ colStart_(nStart)
+{
+ #ifdef FULLDEBUG
+ if
+ (
+ rowStart_ + mRows_ > matrix.m()
+ || colStart_ + nCols_ > matrix.n()
+ )
+ {
+ FatalErrorInFunction
+ << "Block addresses outside matrix"
+ << abort(FatalError);
+ }
+ #endif
+}
+
+
+template<class MatrixType>
+Foam::MatrixBlock<MatrixType>::MatrixBlock
+(
+ MatrixType& matrix,
+ const label m,
+ const label n,
+ const label mStart,
+ const label nStart
+)
+:
+ matrix_(matrix),
+ mRows_(m),
+ nCols_(n),
+ rowStart_(mStart),
+ colStart_(nStart)
+{
+ #ifdef FULLDEBUG
+ if
+ (
+ rowStart_ + mRows_ > matrix.m()
+ || colStart_ + nCols_ > matrix.n()
+ )
+ {
+ FatalErrorInFunction
+ << "Block addresses outside matrix"
+ << abort(FatalError);
+ }
+ #endif
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+template<class MatrixType>
+inline Foam::label Foam::ConstMatrixBlock<MatrixType>::m() const
+{
+ return mRows_;
+}
+
+
+template<class MatrixType>
+inline Foam::label Foam::ConstMatrixBlock<MatrixType>::n() const
+{
+ return nCols_;
+}
+
+
+template<class MatrixType>
+inline Foam::label Foam::MatrixBlock<MatrixType>::m() const
+{
+ return mRows_;
+}
+
+
+template<class MatrixType>
+inline Foam::label Foam::MatrixBlock<MatrixType>::n() const
+{
+ return nCols_;
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
+
+template<class MatrixType>
+inline const typename MatrixType::cmptType&
+Foam::ConstMatrixBlock<MatrixType>::operator()
+(
+ const label i,
+ const label j
+) const
+{
+ #ifdef FULLDEBUG
+ if (i<0 || i>=mRows_)
+ {
+ FatalErrorInFunction
+ << "Index " << i << " out of range 0 ... " << mRows_-1
+ << abort(FatalError);
+ }
+ if (j<0 || j>=nCols_)
+ {
+ FatalErrorInFunction
+ << "Index " << j << " out of range 0 ... " << nCols_-1
+ << abort(FatalError);
+ }
+ #endif
+
+ return matrix_(i + rowStart_, j + colStart_);
+}
+
+
+template<class MatrixType>
+inline const typename MatrixType::cmptType&
+Foam::MatrixBlock<MatrixType>::operator()
+(
+ const label i,
+ const label j
+) const
+{
+ #ifdef FULLDEBUG
+ if (i<0 || i>=mRows_)
+ {
+ FatalErrorInFunction
+ << "Index " << i << " out of range 0 ... " << mRows_-1
+ << abort(FatalError);
+ }
+ if (j<0 || j>=nCols_)
+ {
+ FatalErrorInFunction
+ << "Index " << j << " out of range 0 ... " << nCols_-1
+ << abort(FatalError);
+ }
+ #endif
+
+ return matrix_(i + rowStart_, j + colStart_);
+}
+
+
+template<class MatrixType>
+inline typename MatrixType::cmptType&
+Foam::MatrixBlock<MatrixType>::operator()
+(
+ const label i,
+ const label j
+)
+{
+ #ifdef FULLDEBUG
+ if (i<0 || i>=mRows_)
+ {
+ FatalErrorInFunction
+ << "Index " << i << " out of range 0 ... " << mRows_-1
+ << abort(FatalError);
+ }
+ if (j<0 || j>=nCols_)
+ {
+ FatalErrorInFunction
+ << "Index " << j << " out of range 0 ... " << nCols_-1
+ << abort(FatalError);
+ }
+ #endif
+
+ return matrix_(i + rowStart_, j + colStart_);
+}
+
+
+// ************************************************************************* //
--
GitLab