diff --git a/src/OpenFOAM/primitives/MatrixSpace/MatrixSpace.H b/src/OpenFOAM/primitives/MatrixSpace/MatrixSpace.H
new file mode 100644
index 0000000000..1d7bcec20d
--- /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 0000000000..e8df098a94
--- /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 cf92644c16..ce5d305be7 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 35bfdef0c4..5673703285 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)
 {}
 
 
@@ -106,57 +124,31 @@ inline Foam::Tensor<Cmpt>::Tensor
 }
 
 
+template<class Cmpt>
+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
+)
+:
+    Tensor::msType(block)
+{}
+
+
 template<class Cmpt>
 inline Foam::Tensor<Cmpt>::Tensor(Istream& is)
 :
-    VectorSpace<Tensor<Cmpt>, Cmpt, 9>(is)
+    Tensor::msType(is)
 {}
 
 
 // * * * * * * * * * * * * * * * 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]);
-}
-
-
-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 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 27c1bff3c1..1d47ae3a07 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 fe4a47be48..10f6964e7b 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 4214f32f0f..5992444b8b 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
 {};