SquareMatrix, RectangularMatrix: Updated block handling

Added 'typeOfInnerProduct' support to ensure the correct type is
returned from the matrix product operator.

src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrix.H

@@ -39,6 +39,7 @@ SourceFiles
 #define RectangularMatrix_H
 
 #include "Matrix.H"
+#include "SquareMatrix.H"
 #include "Identity.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -67,7 +68,12 @@ public:
         inline RectangularMatrix(const label m, const label n);
 
         //- Construct from a block of another matrix
-        inline RectangularMatrix(const typename RectangularMatrix::Block&);
+        template<class MatrixType>
+        inline RectangularMatrix(const ConstMatrixBlock<MatrixType>&);
+
+        //- Construct from a block of another matrix
+        template<class MatrixType>
+        inline RectangularMatrix(const MatrixBlock<MatrixType>&);
 
         //- Construct with given number of rows and columns
         //  initializing all elements to zero
@@ -81,6 +87,9 @@ public:
         //  and value for all elements.
         inline RectangularMatrix(const label m, const label n, const Type&);
 
+        //- Construct as copy of a square matrix
+        inline RectangularMatrix(const SquareMatrix<Type>&);
+
         //- Construct from Istream.
         inline RectangularMatrix(Istream&);
 
@@ -97,6 +106,31 @@ public:
 
 // Global functions and operators
 
+template<class Type>
+class typeOfInnerProduct<Type, RectangularMatrix<Type>, RectangularMatrix<Type>>
+{
+public:
+
+    typedef RectangularMatrix<Type> type;
+};
+
+template<class Type>
+class typeOfInnerProduct<Type, RectangularMatrix<Type>, SquareMatrix<Type>>
+{
+public:
+
+    typedef RectangularMatrix<Type> type;
+};
+
+template<class Type>
+class typeOfInnerProduct<Type, SquareMatrix<Type>, RectangularMatrix<Type>>
+{
+public:
+
+    typedef RectangularMatrix<Type> type;
+};
+
+
 template<class Type>
 RectangularMatrix<Type> outer(const Field<Type>& f1, const Field<Type>& f2);
 

src/OpenFOAM/matrices/RectangularMatrix/RectangularMatrixI.H

@@ -44,9 +44,21 @@ inline Foam::RectangularMatrix<Type>::RectangularMatrix
 
 
 template<class Type>
+template<class MatrixType>
 inline Foam::RectangularMatrix<Type>::RectangularMatrix
 (
-    const typename RectangularMatrix::Block& block
+    const ConstMatrixBlock<MatrixType>& block
+)
+:
+    Matrix<RectangularMatrix<Type>, Type>(block)
+{}
+
+
+template<class Type>
+template<class MatrixType>
+inline Foam::RectangularMatrix<Type>::RectangularMatrix
+(
+    const MatrixBlock<MatrixType>& block
 )
 :
     Matrix<RectangularMatrix<Type>, Type>(block)
@@ -74,7 +86,7 @@ inline Foam::RectangularMatrix<Type>::RectangularMatrix
 :
     Matrix<RectangularMatrix<Type>, Type>(n, n, Zero)
 {
-    for (label i = 0; i < n; ++i)
+    for (label i = 0; i < n; i++)
     {
         this->operator()(i, i) = I;
     }
@@ -93,6 +105,16 @@ inline Foam::RectangularMatrix<Type>::RectangularMatrix
 {}
 
 
+template<class Type>
+inline Foam::RectangularMatrix<Type>::RectangularMatrix
+(
+    const SquareMatrix<Type>& SM
+)
+:
+    Matrix<RectangularMatrix<Type>, Type>(SM)
+{}
+
+
 template<class Type>
 inline Foam::RectangularMatrix<Type>::RectangularMatrix(Istream& is)
 :
@@ -136,9 +158,9 @@ inline Foam::RectangularMatrix<Type> outer
 {
     RectangularMatrix<Type> f1f2T(f1.size(), f2.size());
 
-    for (label i=0; i<f1f2T.m(); ++i)
+    for (label i=0; i<f1f2T.m(); i++)
     {
-        for (label j=0; j<f1f2T.n(); ++j)
+        for (label j=0; j<f1f2T.n(); j++)
         {
             f1f2T(i, j) = f1[i]*f2[j];
         }

src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C

@@ -37,7 +37,7 @@ Foam::scalar Foam::detDecomposed
 {
     Type diagProduct = pTraits<Type>::one;
 
-    for (label i = 0; i < matrix.m(); ++i)
+    for (label i = 0; i < matrix.m(); i++)
     {
         diagProduct *= matrix(i, i);
     }

src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H

@@ -45,6 +45,12 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of friend functions and operators
+
+template<class Type>
+class RectangularMatrix;
+
+
 /*---------------------------------------------------------------------------*\
                            Class Matrix Declaration
 \*---------------------------------------------------------------------------*/
@@ -65,10 +71,22 @@ public:
         //- Construct given number of rows/columns.
         inline SquareMatrix(const label n);
 
+        //- Construct from a block of another matrix
+        template<class MatrixType>
+        inline SquareMatrix(const ConstMatrixBlock<MatrixType>&);
+
+        //- Construct from a block of another matrix
+        template<class MatrixType>
+        inline SquareMatrix(const MatrixBlock<MatrixType>&);
+
         //- Construct given number of rows/columns
         //  initializing all elements to zero
         inline SquareMatrix(const label n, const zero);
 
+        //- Construct given number of rows and columns (checked to be equal)
+        //  initializing all elements to zero
+        inline SquareMatrix(const label m, const label n, const zero);
+
         //- Construct given number of rows/columns
         //  Initializing to the identity matrix
         inline SquareMatrix(const label n, const Identity<Type>);
@@ -77,15 +95,25 @@ public:
         //  initializing all elements to the given value
         inline SquareMatrix(const label n, const Type&);
 
+        //- Construct as copy of a RectangularMatrix
+        //  which is checked to be square
+        inline explicit SquareMatrix(const RectangularMatrix<Type>&);
+
         //- Construct from Istream.
         inline SquareMatrix(Istream&);
 
         //- Clone
         inline autoPtr<SquareMatrix<Type>> clone() const;
+
+
+    // Member operators
+
+        //- Assignment of all entries to zero
+        void operator=(const zero);
 };
 
 
-// Global functions
+// Global functions and operators
 
 //- Return the LU decomposed SquareMatrix det
 template<class Type>
@@ -99,6 +127,14 @@ scalar det(const SquareMatrix<Type>&);
 template<class Type>
 scalar det(SquareMatrix<Type>&);
 
+template<class Type>
+class typeOfInnerProduct<Type, SquareMatrix<Type>, SquareMatrix<Type>>
+{
+public:
+
+    typedef SquareMatrix<Type> type;
+};
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -106,7 +142,7 @@ scalar det(SquareMatrix<Type>&);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-    #include "SquareMatrixI.H"
+#include "SquareMatrixI.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H

@@ -39,6 +39,28 @@ inline Foam::SquareMatrix<Type>::SquareMatrix(const label n)
 {}
 
 
+template<class Type>
+template<class MatrixType>
+inline Foam::SquareMatrix<Type>::SquareMatrix
+(
+    const ConstMatrixBlock<MatrixType>& block
+)
+:
+    Matrix<SquareMatrix<Type>, Type>(block)
+{}
+
+
+template<class Type>
+template<class MatrixType>
+inline Foam::SquareMatrix<Type>::SquareMatrix
+(
+    const MatrixBlock<MatrixType>& block
+)
+:
+    Matrix<SquareMatrix<Type>, Type>(block)
+{}
+
+
 template<class Type>
 inline Foam::SquareMatrix<Type>::SquareMatrix
 (
@@ -50,6 +72,26 @@ inline Foam::SquareMatrix<Type>::SquareMatrix
 {}
 
 
+template<class Type>
+inline Foam::SquareMatrix<Type>::SquareMatrix
+(
+    const label m,
+    const label n,
+    const zero
+)
+:
+    Matrix<SquareMatrix<Type>, Type>(m, n, Zero)
+{
+    if (m != n)
+    {
+        FatalErrorInFunction
+            << "Attempt to construct a square matrix "
+            << m << " x " << n << nl
+            << abort(FatalError);
+    }
+}
+
+
 template<class Type>
 inline Foam::SquareMatrix<Type>::SquareMatrix
 (
@@ -59,7 +101,7 @@ inline Foam::SquareMatrix<Type>::SquareMatrix
 :
     Matrix<SquareMatrix<Type>, Type>(n, n, Zero)
 {
-    for (label i = 0; i < n; ++i)
+    for (label i = 0; i < n; i++)
     {
         this->operator()(i, i) = I;
     }
@@ -77,6 +119,24 @@ inline Foam::SquareMatrix<Type>::SquareMatrix
 {}
 
 
+template<class Type>
+inline Foam::SquareMatrix<Type>::SquareMatrix
+(
+    const RectangularMatrix<Type>& RM
+)
+:
+    Matrix<SquareMatrix<Type>, Type>(RM)
+{
+    if (this->m() != this->n())
+    {
+        FatalErrorInFunction
+            << "Attempt to construct a square matrix from a rectangular matrix "
+            << this->m() << " x " << this->n() << nl
+            << abort(FatalError);
+    }
+}
+
+
 template<class Type>
 inline Foam::SquareMatrix<Type>::SquareMatrix(Istream& is)
 :
@@ -92,4 +152,45 @@ Foam::SquareMatrix<Type>::clone() const
 }
 
 
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::SquareMatrix<Type>::operator=(const zero)
+{
+    Matrix<SquareMatrix<Type>, Type>::operator=(Zero);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+inline Foam::SquareMatrix<Type> symmOuter
+(
+    const Field<Type>& f1,
+    const Field<Type>& f2
+)
+{
+    SquareMatrix<Type> f1f2T(f1.size());
+
+    for (label i=0; i<f1f2T.m(); i++)
+    {
+        for (label j=0; j<f1f2T.n(); j++)
+        {
+            f1f2T(i, j) = f1[i]*f2[j];
+        }
+    }
+
+    return f1f2T;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
 // ************************************************************************* //

src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.C

@@ -37,11 +37,11 @@ Foam::SymmetricSquareMatrix<Type> Foam::invDecomposed
 
     SymmetricSquareMatrix<Type> inv(n, Zero);
 
-    for (label i = 0; i < n; ++i)
+    for (label i = 0; i < n; i++)
     {
         inv(i, i) = 1.0/matrix(i, i);
 
-        for (label j = 0; j < i; ++j)
+        for (label j = 0; j < i; j++)
         {
             Type sum = Zero;
 
@@ -89,7 +89,7 @@ Type Foam::detDecomposed(const SymmetricSquareMatrix<Type>& matrix)
 {
     Type diagProduct = pTraits<Type>::one;
 
-    for (label i = 0; i < matrix.m(); ++i)
+    for (label i = 0; i < matrix.m(); i++)
     {
         diagProduct *= matrix(i, i);
     }