diff --git a/src/OpenFOAM/primitives/MatrixSpace/MatrixSpace.H b/src/OpenFOAM/primitives/MatrixSpace/MatrixSpace.H
new file mode 100644
index 0000000000000000000000000000000000000000..1d7bcec20ddfa37c0ef5ba23fb8ad592712b3c25
--- /dev/null
+++ b/src/OpenFOAM/primitives/MatrixSpace/MatrixSpace.H
@@ -0,0 +1,322 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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::MatrixSpace
+
+Description
+ Templated matrix space.
+
+ Template arguments are the Form the matrix space will be used to create,
+ the type of the elements and the number of rows and columns of the matrix.
+
+SourceFiles
+ MatrixSpaceI.H
+
+SeeAlso
+ Foam::VectorSpace
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef MatrixSpace_H
+#define MatrixSpace_H
+
+#include "VectorSpace.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class MatrixSpace Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Form, class Cmpt, direction Nrows, direction Ncols>
+class MatrixSpace
+:
+ public VectorSpace<Form, Cmpt, Nrows*Ncols>
+{
+
+public:
+
+ //- MatrixSpace type
+ typedef MatrixSpace<Form, Cmpt, Nrows, Ncols> msType;
+
+
+ // Member constants
+
+ static const direction nRows = Nrows;
+ static const direction nCols = Ncols;
+
+
+ // Static member functions
+
+ //- Return the number of rows
+ static direction m()
+ {
+ return Nrows;
+ }
+
+ //- Return the number of columns
+ static direction n()
+ {
+ return Ncols;
+ }
+
+ //- Return the identity matrix for square matrix spaces
+ inline static msType identity();
+
+
+ // Sub-Block Classes
+
+ //- Const sub-block type
+ template<class SubTensor, direction BRowStart, direction BColStart>
+ class ConstBlock
+ {
+ //- Reference to parent matrix
+ const msType& matrix_;
+
+ public:
+
+ static const direction nRows = SubTensor::nRows;
+ static const direction nCols = SubTensor::nCols;
+
+ //- Return the number of rows in the block
+ static direction m()
+ {
+ return nRows;
+ }
+
+ //- Return the number of columns in the block
+ static direction n()
+ {
+ return nCols;
+ }
+
+ //- Construct for the given matrix
+ inline ConstBlock(const msType& matrix);
+
+ //- (i, j) const element access operator
+ inline const Cmpt& operator()
+ (
+ const direction i,
+ const direction j
+ ) const;
+ };
+
+
+ //- Sub-block type
+ template
+ <
+ class SubTensor,
+ direction BRowStart,
+ direction BColStart
+ >
+ class Block
+ {
+ //- Reference to parent matrix
+ msType& matrix_;
+
+ public:
+
+ static const direction nRows = SubTensor::nRows;
+ static const direction nCols = SubTensor::nCols;
+
+ //- Return the number of rows in the block
+ static direction m()
+ {
+ return nRows;
+ }
+
+ //- Return the number of columns in the block
+ static direction n()
+ {
+ return nCols;
+ }
+
+ //- Construct for the given matrix
+ inline Block(msType& matrix);
+
+ //- Assignment to a matrix
+ template<class Form2>
+ inline void operator=
+ (
+ const MatrixSpace
+ <
+ Form2,
+ Cmpt,
+ SubTensor::nRows,
+ SubTensor::nCols
+ >& matrix
+ );
+
+ //- Assignment to a column vector
+ template<class VSForm>
+ inline void operator=
+ (
+ const VectorSpace<VSForm, Cmpt, SubTensor::nRows>& v
+ );
+
+ //- (i, j) const element access operator
+ inline const Cmpt& operator()
+ (
+ const direction i,
+ const direction j
+ ) const;
+
+ //- (i, j) element access operator
+ inline Cmpt& operator()(const direction i, const direction j);
+ };
+
+
+ // Constructors
+
+ //- Construct null
+ inline MatrixSpace();
+
+ //- Construct initialized to zero
+ inline explicit MatrixSpace(const Foam::zero);
+
+ //- Construct as copy of a VectorSpace with the same size
+ template<class Form2, class Cmpt2>
+ inline explicit MatrixSpace
+ (
+ const VectorSpace<Form2, Cmpt2, Nrows*Ncols>&
+ );
+
+ //- Construct from a block of another matrix space
+ template
+ <
+ template<class, direction, direction> class Block2,
+ direction BRowStart,
+ direction BColStart
+ >
+ inline MatrixSpace
+ (
+ const Block2<Form, BRowStart, BColStart>& block
+ );
+
+ //- Construct from Istream
+ MatrixSpace(Istream&);
+
+
+ // Member Functions
+
+ //- Fast const element access using compile-time addressing
+ template<direction Row, direction Col>
+ inline const Cmpt& elmt() const;
+
+ //- Fast element access using compile-time addressing
+ template<direction Row, direction Col>
+ inline Cmpt& elmt();
+
+ // Const element access functions for a 3x3
+ // Compile-time errors are generated for inappropriate use
+
+ inline const Cmpt& xx() const;
+ inline const Cmpt& xy() const;
+ inline const Cmpt& xz() const;
+ inline const Cmpt& yx() const;
+ inline const Cmpt& yy() const;
+ inline const Cmpt& yz() const;
+ inline const Cmpt& zx() const;
+ inline const Cmpt& zy() const;
+ inline const Cmpt& zz() const;
+
+ // Element access functions for a 3x3
+ // Compile-time errors are generated for inappropriate use
+
+ inline Cmpt& xx();
+ inline Cmpt& xy();
+ inline Cmpt& xz();
+ inline Cmpt& yx();
+ inline Cmpt& yy();
+ inline Cmpt& yz();
+ inline Cmpt& zx();
+ inline Cmpt& zy();
+ inline Cmpt& zz();
+
+ //- Return the transpose of the matrix
+ inline typename typeOfTranspose<Cmpt, Form>::type T() const;
+
+ //- Return a const sub-block corresponding to the specified type
+ // starting at the specified row and column
+ template<class SubTensor, direction BRowStart, direction BColStart>
+ inline ConstBlock<SubTensor, BRowStart, BColStart> block() const;
+
+ //- Return a sub-block corresponding to the specified type
+ // starting at the specified row and column
+ template<class SubTensor, direction BRowStart, direction BColStart>
+ inline Block<SubTensor, BRowStart, BColStart> block();
+
+ //- (i, j) const element access operator
+ inline const Cmpt& operator()
+ (
+ const direction& row,
+ const direction& col
+ ) const;
+
+ //- (i, j) element access operator
+ inline Cmpt& operator()(const direction& row, const direction& col);
+
+
+ // Member Operators
+
+ //- Assignment to zero
+ inline void operator=(const Foam::zero);
+
+ //- Assignment to a block of another matrix space
+ template
+ <
+ template<class, direction, direction> class Block2,
+ direction BRowStart,
+ direction BColStart
+ >
+ inline void operator=
+ (
+ const Block2<Form, BRowStart, BColStart>& block
+ );
+
+ //- Inner product with a compatible square matrix
+ template<class Form2>
+ inline void operator&=
+ (
+ const MatrixSpace<Form, Cmpt, Ncols, Ncols>& matrix
+ );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "MatrixSpaceI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/primitives/MatrixSpace/MatrixSpaceI.H b/src/OpenFOAM/primitives/MatrixSpace/MatrixSpaceI.H
new file mode 100644
index 0000000000000000000000000000000000000000..e8df098a9489e14c7bfa2a579ad330ba60ad64c9
--- /dev/null
+++ b/src/OpenFOAM/primitives/MatrixSpace/MatrixSpaceI.H
@@ -0,0 +1,609 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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 "StaticAssert.H"
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+inline Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::MatrixSpace()
+{}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+inline Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::MatrixSpace
+(
+ const Foam::zero z
+)
+:
+ MatrixSpace::vsType(z)
+{}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+template<class Form2, class Cmpt2>
+inline Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::MatrixSpace
+(
+ const VectorSpace<Form2, Cmpt2, Nrows*Ncols>& vs
+)
+:
+ MatrixSpace::vsType(vs)
+{}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+template
+<
+ template<class, Foam::direction, Foam::direction> class Block2,
+ Foam::direction BRowStart,
+ Foam::direction BColStart
+>
+inline Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::MatrixSpace
+(
+ const Block2<Form, BRowStart, BColStart>& block
+)
+{
+ for (direction i=0; i<Nrows; ++i)
+ {
+ for (direction j=0; j<Ncols; ++j)
+ {
+ operator()(i, j) = block(i, j);
+ }
+ }
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+inline Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::MatrixSpace(Istream& is)
+:
+ MatrixSpace::vsType(is)
+{}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+template<class SubTensor, Foam::direction BRowStart, Foam::direction BColStart>
+inline Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::
+ConstBlock<SubTensor, BRowStart, BColStart>::
+ConstBlock(const msType& matrix)
+:
+ matrix_(matrix)
+{
+ StaticAssert(msType::nRows >= BRowStart + nRows);
+ StaticAssert(msType::nCols >= BColStart + nCols);
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+template<class SubTensor, Foam::direction BRowStart, Foam::direction BColStart>
+inline Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::
+Block<SubTensor, BRowStart, BColStart>::
+Block(msType& matrix)
+:
+ matrix_(matrix)
+{
+ StaticAssert(msType::nRows >= BRowStart + nRows);
+ StaticAssert(msType::nCols >= BColStart + nCols);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+template<Foam::direction Row, Foam::direction Col>
+inline const Cmpt& Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::elmt() const
+{
+ StaticAssert(Row < Nrows && Col < Ncols);
+ return this->v_[Row*Ncols + Col];
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+template<Foam::direction Row, Foam::direction Col>
+inline Cmpt& Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::elmt()
+{
+ StaticAssert(Row < Nrows && Col < Ncols);
+ return this->v_[Row*Ncols + Col];
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+inline const Cmpt& Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::xx() const
+{
+ return elmt<0, 0>();
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+inline Cmpt& Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::xx()
+{
+ return elmt<0, 0>();
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+inline const Cmpt& Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::xy() const
+{
+ return elmt<0,1>();
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+inline Cmpt& Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::xy()
+{
+ return elmt<0,1>();
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+inline const Cmpt& Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::xz() const
+{
+ return elmt<0,2>();
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+inline Cmpt& Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::xz()
+{
+ return elmt<0,2>();
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+inline const Cmpt& Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::yx() const
+{
+ return elmt<1,0>();
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+inline Cmpt& Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::yx()
+{
+ return elmt<1,0>();
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+inline const Cmpt& Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::yy() const
+{
+ return elmt<1,1>();
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+inline Cmpt& Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::yy()
+{
+ return elmt<1,1>();
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+inline const Cmpt& Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::yz() const
+{
+ return elmt<1,2>();
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+inline Cmpt& Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::yz()
+{
+ return elmt<1,2>();
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+inline const Cmpt& Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::zx() const
+{
+ return elmt<2,0>();
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+inline Cmpt& Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::zx()
+{
+ return elmt<2,0>();
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+inline const Cmpt& Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::zy() const
+{
+ return elmt<2,1>();
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+inline Cmpt& Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::zy()
+{
+ return elmt<2,1>();
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+inline const Cmpt& Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::zz() const
+{
+ return elmt<2,2>();
+}
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+inline Cmpt& Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::zz()
+{
+ return elmt<2,2>();
+}
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+inline Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>
+Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::identity()
+{
+ StaticAssert(Nrows == Ncols);
+ msType result((Foam::zero()));
+
+ for (direction i=0; i<Ncols; ++i)
+ {
+ result(i, i) = 1;
+ }
+
+ return result;
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+inline typename Foam::typeOfTranspose<Cmpt, Form>::type
+Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::T() const
+{
+ typename typeOfTranspose<Cmpt, Form>::type result;
+
+ for (direction i=0; i<Nrows; ++i)
+ {
+ for (direction j=0; j<Ncols; ++j)
+ {
+ result(j,i) = operator()(i, j);
+ }
+ }
+
+ return result;
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+template
+<
+ class SubTensor,
+ Foam::direction BRowStart,
+ Foam::direction BColStart
+>
+inline typename Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::template
+ ConstBlock<SubTensor, BRowStart, BColStart>
+Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::block() const
+{
+ return *this;
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+template
+<
+ class SubTensor,
+ Foam::direction BRowStart,
+ Foam::direction BColStart
+>
+inline
+typename Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::template
+ Block<SubTensor, BRowStart, BColStart>
+Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::block()
+{
+ return *this;
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+inline const Cmpt& Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::operator()
+(
+ const direction& row,
+ const direction& col
+) const
+{
+ #ifdef FULLDEBUG
+ if (row < 0 || row > Nrows-1 || col < 0 || col > Ncols-1)
+ {
+ FatalErrorInFunction
+ << "indices out of range"
+ << abort(FatalError);
+ }
+ #endif
+
+ return this->v_[row*Ncols + col];
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+inline Cmpt& Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::operator()
+(
+ const direction& row,
+ const direction& col
+)
+{
+ #ifdef FULLDEBUG
+ if (row < 0 || row > Nrows-1 || col < 0 || col > Ncols-1)
+ {
+ FatalErrorInFunction
+ << "indices out of range"
+ << abort(FatalError);
+ }
+ #endif
+
+ return this->v_[row*Ncols + col];
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+template<class SubTensor, Foam::direction BRowStart, Foam::direction BColStart>
+inline const Cmpt&
+Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::
+ConstBlock<SubTensor, BRowStart, BColStart>::
+operator()(const direction i, const direction j) const
+{
+ return matrix_(BRowStart + i, BColStart + j);
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+template<class SubTensor, Foam::direction BRowStart, Foam::direction BColStart>
+inline const Cmpt&
+Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::
+Block<SubTensor, BRowStart, BColStart>::
+operator()(const direction i, const direction j) const
+{
+ return matrix_(BRowStart + i, BColStart + j);
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+template<class SubTensor, Foam::direction BRowStart, Foam::direction BColStart>
+inline Cmpt&
+Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::
+Block<SubTensor, BRowStart, BColStart>::
+operator()(const direction i, const direction j)
+{
+ return matrix_(BRowStart + i, BColStart + j);
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+inline void Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::operator=
+(
+ const Foam::zero z
+)
+{
+ MatrixSpace::vsType::operator=(z);
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+template<class Form2>
+inline void Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::operator&=
+(
+ const MatrixSpace<Form, Cmpt, Ncols, Ncols>& matrix
+)
+{
+ *this = *this & matrix;
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+template
+<
+ template<class, Foam::direction, Foam::direction> class Block2,
+ Foam::direction BRowStart,
+ Foam::direction BColStart
+>
+inline void Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::operator=
+(
+ const Block2<Form, BRowStart, BColStart>& block
+)
+{
+ for (direction i = 0; i < Nrows; ++i)
+ {
+ for (direction j = 0; j < Ncols; ++j)
+ {
+ operator()(i, j) = block(i, j);
+ }
+ }
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+template<class SubTensor, Foam::direction BRowStart, Foam::direction BColStart>
+template<class Form2>
+inline void
+Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::
+Block<SubTensor, BRowStart, BColStart>::
+operator=
+(
+ const MatrixSpace<Form2, Cmpt, SubTensor::nRows, SubTensor::nCols>& matrix
+)
+{
+ for (direction i=0; i<nRows; ++i)
+ {
+ for (direction j=0; j<nCols; ++j)
+ {
+ operator()(i,j) = matrix(i,j);
+ }
+ }
+}
+
+
+template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
+template<class SubTensor, Foam::direction BRowStart, Foam::direction BColStart>
+template<class VSForm>
+inline void
+Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::
+Block<SubTensor, BRowStart, BColStart>::
+operator=
+(
+ const VectorSpace<VSForm, Cmpt, SubTensor::nRows>& v
+)
+{
+ StaticAssert(nCols == 1);
+
+ for (direction i=0; i<SubTensor::nRows; ++i)
+ {
+ operator()(i,0) = v[i];
+ }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
+
+template<class Form, class Cmpt, direction Nrows, direction Ncols>
+inline typename typeOfTranspose<Cmpt, Form>::type T
+(
+ const MatrixSpace<Form, Cmpt, Ncols, Nrows>& matrix
+)
+{
+ return matrix.T();
+}
+
+
+template<class Form, class Cmpt, direction Ncmpts>
+inline typename typeOfTranspose<Cmpt, Form>::type T
+(
+ const VectorSpace<Form, Cmpt, Ncmpts>& v
+)
+{
+ typename typeOfTranspose<Cmpt, Form>::type result;
+
+ for (direction i=0; i<Ncmpts; ++i)
+ {
+ result[i] = v[i];
+ }
+
+ return result;
+}
+
+
+// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
+
+template
+<
+ class Form1,
+ class Form2,
+ class Cmpt,
+ direction Nrows1,
+ direction Ncols1,
+ direction Nrows2,
+ direction Ncols2
+>
+inline typename typeOfInnerProduct<Cmpt, Form1, Form2>::type operator&
+(
+ const MatrixSpace<Form1, Cmpt, Nrows1, Ncols1>& matrix1,
+ const MatrixSpace<Form2, Cmpt, Nrows2, Ncols2>& matrix2
+)
+{
+ StaticAssert(Ncols1 == Nrows2);
+
+ typename typeOfInnerProduct<Cmpt, Form1, Form2>::type result
+ (
+ (Foam::zero())
+ );
+
+ for (direction i=0; i<Nrows1; ++i)
+ {
+ for (direction j=0; j<Ncols2; ++j)
+ {
+ for (direction k=0; k<Nrows2; k++)
+ {
+ result(i, j) += matrix1(i, k)*matrix2(k, j);
+ }
+ }
+ }
+
+ return result;
+}
+
+
+template<class Form, class VSForm, class Cmpt, direction Nrows, direction Ncols>
+inline typename typeOfInnerProduct<Cmpt, Form, VSForm>::type operator&
+(
+ const MatrixSpace<Form, Cmpt, Nrows, Ncols>& matrix,
+ const VectorSpace<VSForm, Cmpt, Ncols>& v
+)
+{
+ typename typeOfInnerProduct<Cmpt, Form, VSForm>::type result
+ (
+ (Foam::zero())
+ );
+
+ for (direction i=0; i<Nrows; ++i)
+ {
+ for (direction j=0; j<Ncols; ++j)
+ {
+ result[i] += matrix(i, j)*v[j];
+ }
+ }
+
+ return result;
+}
+
+
+template
+<
+ class Form1,
+ class Form2,
+ class Cmpt,
+ direction Ncmpts1,
+ direction Ncmpts2
+>
+inline typename typeOfOuterProduct<Cmpt, Form1, Form2>::type operator*
+(
+ const VectorSpace<Form1, Cmpt, Ncmpts1>& v1,
+ const VectorSpace<Form2, Cmpt, Ncmpts2>& v2
+)
+{
+ typename typeOfOuterProduct<Cmpt, Form1, Form2>::type result;
+
+ for (direction i=0; i<Ncmpts1; ++i)
+ {
+ for (direction j=0; j<Ncmpts2; ++j)
+ {
+ result(i, j) = v1[i]*v2[j];
+ }
+ }
+
+ return result;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/primitives/Tensor/Tensor.H b/src/OpenFOAM/primitives/Tensor/Tensor.H
index cf92644c16d2b667ac61ffe15b2effe9d258b5e8..ce5d305be7ff07446663c7f7779af31cc49f9cb6 100644
--- a/src/OpenFOAM/primitives/Tensor/Tensor.H
+++ b/src/OpenFOAM/primitives/Tensor/Tensor.H
@@ -25,7 +25,7 @@ Class
Foam::Tensor
Description
- Templated 3D tensor derived from VectorSpace adding construction from
+ Templated 3D tensor derived from MatrixSpace adding construction from
9 components, element access using xx(), xy() etc. member functions and
the inner-product (dot-product) and outer-product of two Vectors
(tensor-product) operators.
@@ -33,11 +33,16 @@ Description
SourceFiles
TensorI.H
+SeeAlso
+ Foam::MatrixSpace
+ Foam::Vector
+
\*---------------------------------------------------------------------------*/
#ifndef Tensor_H
#define Tensor_H
+#include "MatrixSpace.H"
#include "Vector.H"
#include "SphericalTensor.H"
@@ -56,7 +61,7 @@ class SymmTensor;
template<class Cmpt>
class Tensor
:
- public VectorSpace<Tensor<Cmpt>, Cmpt, 9>
+ public MatrixSpace<Tensor<Cmpt>, Cmpt, 3, 3>
{
public:
@@ -85,6 +90,13 @@ public:
//- Construct null
inline Tensor();
+ //- Construct initialized to zero
+ inline explicit Tensor(const Foam::zero);
+
+ //- Construct given MatrixSpace of the same rank
+ template<class Cmpt2>
+ inline Tensor(const MatrixSpace<Tensor<Cmpt2>, Cmpt2, 3, 3>&);
+
//- Construct given VectorSpace of the same rank
template<class Cmpt2>
inline Tensor(const VectorSpace<Tensor<Cmpt2>, Cmpt2, 9>&);
@@ -114,13 +126,25 @@ public:
const Cmpt tzx, const Cmpt tzy, const Cmpt tzz
);
+ //- Construct from a block of another matrix space
+ template
+ <
+ template<class, direction, direction> class Block2,
+ direction BRowStart,
+ direction BColStart
+ >
+ Tensor
+ (
+ const Block2<Tensor<Cmpt>, BRowStart, BColStart>& block
+ );
+
//- Construct from Istream
- Tensor(Istream&);
+ inline Tensor(Istream&);
// Member Functions
- // Access
+ // Component access
inline const Cmpt& xx() const;
inline const Cmpt& xy() const;
@@ -142,7 +166,7 @@ public:
inline Cmpt& zy();
inline Cmpt& zz();
- // Access vector components.
+ // Row-vector access.
inline Vector<Cmpt> x() const;
inline Vector<Cmpt> y() const;
@@ -161,6 +185,10 @@ public:
//- Inner-product with a Tensor
inline void operator&=(const Tensor<Cmpt>&);
+ //- Assign to an equivalent vector space
+ template<class Cmpt2>
+ inline void operator=(const VectorSpace<Tensor<Cmpt2>, Cmpt2, 9>&);
+
//- Assign to a SphericalTensor
inline void operator=(const SphericalTensor<Cmpt>&);
@@ -181,6 +209,15 @@ public:
};
+template<class Cmpt>
+class typeOfTranspose<Cmpt, Tensor<Cmpt>>
+{
+public:
+
+ typedef Tensor<Cmpt> type;
+};
+
+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
diff --git a/src/OpenFOAM/primitives/Tensor/TensorI.H b/src/OpenFOAM/primitives/Tensor/TensorI.H
index 35bfdef0c469f823a546c9b796e9798f06f5daa6..56737032858015fcd5cc78bded6bf6d8ab516e07 100644
--- a/src/OpenFOAM/primitives/Tensor/TensorI.H
+++ b/src/OpenFOAM/primitives/Tensor/TensorI.H
@@ -32,6 +32,24 @@ inline Foam::Tensor<Cmpt>::Tensor()
{}
+template<class Cmpt>
+inline Foam::Tensor<Cmpt>::Tensor(const Foam::zero z)
+:
+ Tensor::msType(z)
+{}
+
+
+template<class Cmpt>
+template<class Cmpt2>
+inline Foam::Tensor<Cmpt>::Tensor
+(
+ const MatrixSpace<Tensor<Cmpt2>, Cmpt2, 3, 3>& vs
+)
+:
+ Tensor::msType(vs)
+{}
+
+
template<class Cmpt>
template<class Cmpt2>
inline Foam::Tensor<Cmpt>::Tensor
@@ -39,7 +57,7 @@ inline Foam::Tensor<Cmpt>::Tensor
const VectorSpace<Tensor<Cmpt2>, Cmpt2, 9>& vs
)
:
- VectorSpace<Tensor<Cmpt>, Cmpt, 9>(vs)
+ Tensor::msType(vs)
{}
@@ -107,55 +125,29 @@ inline Foam::Tensor<Cmpt>::Tensor
template<class Cmpt>
-inline Foam::Tensor<Cmpt>::Tensor(Istream& is)
+template
+<
+ template<class, Foam::direction, Foam::direction> class Block2,
+ Foam::direction BRowStart,
+ Foam::direction BColStart
+>
+inline Foam::Tensor<Cmpt>::Tensor
+(
+ const Block2<Tensor<Cmpt>, BRowStart, BColStart>& block
+)
:
- VectorSpace<Tensor<Cmpt>, Cmpt, 9>(is)
+ Tensor::msType(block)
{}
-// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-
-template<class Cmpt>
-inline Foam::Vector<Cmpt> Foam::Tensor<Cmpt>::x() const
-{
- return Vector<Cmpt>(this->v_[XX], this->v_[XY], this->v_[XZ]);
-}
-
-
template<class Cmpt>
-inline Foam::Vector<Cmpt> Foam::Tensor<Cmpt>::y() const
-{
- return Vector<Cmpt>(this->v_[YX], this->v_[YY], this->v_[YZ]);
-}
-
-
-template<class Cmpt>
-inline Foam::Vector<Cmpt> Foam::Tensor<Cmpt>::z() const
-{
- return Vector<Cmpt>(this->v_[ZX], this->v_[ZY], this->v_[ZZ]);
-}
+inline Foam::Tensor<Cmpt>::Tensor(Istream& is)
+:
+ Tensor::msType(is)
+{}
-template<class Cmpt>
-inline Foam::Vector<Cmpt> Foam::Tensor<Cmpt>::vectorComponent
-(
- const direction cmpt
-) const
-{
- switch (cmpt)
- {
- case 0:
- return x();
- break;
- case 1:
- return y();
- break;
- case 2:
- return z();
- break;
- }
-}
-
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Cmpt>
inline const Cmpt& Foam::Tensor<Cmpt>::xx() const
@@ -283,6 +275,48 @@ inline Cmpt& Foam::Tensor<Cmpt>::zz()
}
+template<class Cmpt>
+inline Foam::Vector<Cmpt> Foam::Tensor<Cmpt>::x() const
+{
+ return Vector<Cmpt>(this->v_[XX], this->v_[XY], this->v_[XZ]);
+}
+
+
+template<class Cmpt>
+inline Foam::Vector<Cmpt> Foam::Tensor<Cmpt>::y() const
+{
+ return Vector<Cmpt>(this->v_[YX], this->v_[YY], this->v_[YZ]);
+}
+
+
+template<class Cmpt>
+inline Foam::Vector<Cmpt> Foam::Tensor<Cmpt>::z() const
+{
+ return Vector<Cmpt>(this->v_[ZX], this->v_[ZY], this->v_[ZZ]);
+}
+
+
+template<class Cmpt>
+inline Foam::Vector<Cmpt> Foam::Tensor<Cmpt>::vectorComponent
+(
+ const direction cmpt
+) const
+{
+ switch (cmpt)
+ {
+ case 0:
+ return x();
+ break;
+ case 1:
+ return y();
+ break;
+ case 2:
+ return z();
+ break;
+ }
+}
+
+
template<class Cmpt>
inline Foam::Tensor<Cmpt> Foam::Tensor<Cmpt>::T() const
{
@@ -320,6 +354,17 @@ inline void Foam::Tensor<Cmpt>::operator&=(const Tensor<Cmpt>& t)
}
+template<class Cmpt>
+template<class Cmpt2>
+inline void Foam::Tensor<Cmpt>::operator=
+(
+ const VectorSpace<Tensor<Cmpt2>, Cmpt2, 9>& vs
+)
+{
+ VectorSpace<Tensor<Cmpt>, Cmpt, 9>::operator=(vs);
+}
+
+
template<class Cmpt>
inline void Foam::Tensor<Cmpt>::operator=(const SphericalTensor<Cmpt>& st)
{
diff --git a/src/OpenFOAM/primitives/VectorSpace/VectorSpace.H b/src/OpenFOAM/primitives/VectorSpace/VectorSpace.H
index 27c1bff3c1cd555208c932f146de793564a36007..1d47ae3a0788fa78a4ce1eacef24cd625b54c2cd 100644
--- a/src/OpenFOAM/primitives/VectorSpace/VectorSpace.H
+++ b/src/OpenFOAM/primitives/VectorSpace/VectorSpace.H
@@ -88,7 +88,7 @@ public:
typedef Cmpt cmptType;
- // Member constants
+ // Static constants
//- Dimensionality of space
static const direction dim = 3;
@@ -97,6 +97,13 @@ public:
static const direction nComponents = Ncmpts;
+ // VectorSpace currently defaults to a column-vector
+ // This will be removed when column-vector is introduced
+ // as a specialization
+ static const direction nRows = Ncmpts;
+ static const direction nCols = 1;
+
+
// Static data members
static const char* const typeName;
@@ -109,6 +116,38 @@ public:
static const Form rootMin;
+ // Sub-Block Classes
+
+ //- Const sub-block type
+ template
+ <
+ class SubVector,
+ direction BStart
+ >
+ class ConstBlock
+ {
+ const vsType& vs_;
+
+ public:
+
+ //- Number of components in this vector space
+ static const direction nComponents = SubVector::nComponents;
+
+ //- Construct for a given vector
+ inline ConstBlock(const vsType& vs);
+
+ //- [i] const element access operator
+ inline const Cmpt& operator[](const direction i) const;
+
+ //- (i, 0) const element access operator
+ inline const Cmpt& operator()
+ (
+ const direction i,
+ const direction
+ ) const;
+ };
+
+
// Constructors
//- Construct null
@@ -123,9 +162,9 @@ public:
//- Construct as copy
inline VectorSpace(const VectorSpace<Form, Cmpt, Ncmpts>&);
- //- Construct as copy of another VectorSpace type of the same rank
+ //- Construct as copy of a VectorSpace with the same size
template<class Form2, class Cmpt2>
- inline VectorSpace(const VectorSpace<Form2, Cmpt2, Ncmpts>&);
+ inline explicit VectorSpace(const VectorSpace<Form2, Cmpt2, Ncmpts>&);
// Member Functions
@@ -142,6 +181,9 @@ public:
//- Return a VectorSpace with all elements = s
inline static Form uniform(const Cmpt& s);
+ template<class SubVector, direction BStart>
+ inline const ConstBlock<SubVector, BStart> block() const;
+
// Member Operators
diff --git a/src/OpenFOAM/primitives/VectorSpace/VectorSpaceI.H b/src/OpenFOAM/primitives/VectorSpace/VectorSpaceI.H
index fe4a47be48ebbf6e60d164c4bbdd358f7deceb80..10f6964e7b2ff8376b38dba56a6f4e65c9c52347 100644
--- a/src/OpenFOAM/primitives/VectorSpace/VectorSpaceI.H
+++ b/src/OpenFOAM/primitives/VectorSpace/VectorSpaceI.H
@@ -27,6 +27,7 @@ License
#include "products.H"
#include "VectorSpaceOps.H"
#include "ops.H"
+#include "StaticAssert.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -68,6 +69,20 @@ inline VectorSpace<Form, Cmpt, Ncmpts>::VectorSpace
}
+template<class Form, class Cmpt, direction Ncmpts>
+template<class SubVector, direction BStart>
+inline
+VectorSpace<Form, Cmpt, Ncmpts>::ConstBlock<SubVector, BStart>::ConstBlock
+(
+ const vsType& vs
+)
+:
+ vs_(vs)
+{
+ StaticAssert(vsType::nComponents >= BStart + nComponents);
+}
+
+
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Form, class Cmpt, direction Ncmpts>
@@ -164,6 +179,16 @@ inline Form VectorSpace<Form, Cmpt, Ncmpts>::uniform(const Cmpt& s)
}
+template<class Form, class Cmpt, direction Ncmpts>
+template<class SubVector, direction BStart>
+inline const typename VectorSpace<Form, Cmpt, Ncmpts>::template
+ ConstBlock<SubVector, BStart>
+VectorSpace<Form, Cmpt, Ncmpts>::block() const
+{
+ return *this;
+}
+
+
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
template<class Form, class Cmpt, direction Ncmpts>
@@ -204,6 +229,58 @@ inline Cmpt& VectorSpace<Form, Cmpt, Ncmpts>::operator[]
}
+template<class Form, class Cmpt, direction Ncmpts>
+template<class SubVector, direction BStart>
+inline const Cmpt&
+VectorSpace<Form, Cmpt, Ncmpts>::
+ConstBlock<SubVector, BStart>::operator[]
+(
+ const direction i
+) const
+{
+ #ifdef FULLDEBUG
+ if (i >= Ncmpts)
+ {
+ FatalErrorInFunction
+ << "index out of range"
+ << abort(FatalError);
+ }
+ #endif
+
+ return vs_[BStart + i];
+}
+
+
+template<class Form, class Cmpt, direction Ncmpts>
+template<class SubVector, direction BStart>
+inline const Cmpt&
+VectorSpace<Form, Cmpt, Ncmpts>::
+ConstBlock<SubVector, BStart>::operator()
+(
+ const direction i,
+ const direction j
+) const
+{
+ #ifdef FULLDEBUG
+ if (i >= Ncmpts)
+ {
+ FatalErrorInFunction
+ << "index out of range"
+ << abort(FatalError);
+ }
+
+ if (j != 0)
+ {
+ FatalErrorInFunction
+ << "j != 0"
+ << abort(FatalError);
+ }
+ #endif
+
+ return vs_[BStart + i];
+}
+
+
template<class Form, class Cmpt, direction Ncmpts>
inline void VectorSpace<Form, Cmpt, Ncmpts>::operator=
(
diff --git a/src/OpenFOAM/primitives/VectorSpace/products.H b/src/OpenFOAM/primitives/VectorSpace/products.H
index 4214f32f0f27dce8549764d7e1451fcdc2e81d06..5992444b8b5ac8d2054b8d3b1f074fa352d39408 100644
--- a/src/OpenFOAM/primitives/VectorSpace/products.H
+++ b/src/OpenFOAM/primitives/VectorSpace/products.H
@@ -41,6 +41,24 @@ namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+//- Abstract template class to provide the form resulting from
+// the inner-product of two forms
+template<class Cmpt, class Form1, class Form2>
+class typeOfInnerProduct
+{};
+
+//- Abstract template class to provide the form resulting from
+// the outer-product of two forms
+template<class Cmpt, class Form1, class Form2>
+class typeOfOuterProduct
+{};
+
+//- Abstract template class to provide the transpose form of a form
+template<class Cmpt, class Form>
+class typeOfTranspose
+{};
+
+
template<class Cmpt, direction rank>
class typeOfRank
{};