From 860fa2b28921360215a366bc84c75b305b3379ca Mon Sep 17 00:00:00 2001
From: Vuko Vukcevic <vuko.vukcevic@fsb.hr>
Date: Wed, 9 Sep 2015 15:59:41 +0200
Subject: [PATCH] Missing transpose in BlockLduMatrix TMulCore and new
 formulation of symmetry for block matrix

---
 src/foam/fields/CoeffField/CoeffField.H            |  2 ++
 .../blockLduMatrix/BlockLduMatrix/BlockLduMatrix.C | 14 +++++++++++++-
 .../BlockLduMatrix/BlockLduMatrixATmul.C           |  6 ++++--
 3 files changed, 19 insertions(+), 3 deletions(-)

diff --git a/src/foam/fields/CoeffField/CoeffField.H b/src/foam/fields/CoeffField/CoeffField.H
index a4ff4c7a5..1a7078cdd 100644
--- a/src/foam/fields/CoeffField/CoeffField.H
+++ b/src/foam/fields/CoeffField/CoeffField.H
@@ -257,6 +257,8 @@ public:
                 const CoeffField<Type>& f,
                 const labelList& addr
             );
+
+
     // Member operators
 
         void operator=(const CoeffField<Type>&);
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrix.C b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrix.C
index 86cf42339..b2a97e15d 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrix.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrix.C
@@ -295,7 +295,19 @@ bool Foam::BlockLduMatrix<Type>::symmetric() const
             << abort(FatalError);
     }
 
-    return (diagPtr_ && (!lowerPtr_ && upperPtr_));
+    // Assuming that the matrix is not symmetric if diag or upper are square
+    // type coefficients. VV and HJ, 9/Sep/2015.
+    bool activeSquare = false;
+    if
+    (
+        diagPtr_->activeType() == blockCoeffBase::SQUARE
+     || upperPtr_->activeType() == blockCoeffBase::SQUARE
+    )
+    {
+        activeSquare = true;
+    }
+
+    return (diagPtr_ && (!lowerPtr_ && upperPtr_) && !activeSquare);
 }
 
 
diff --git a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixATmul.C b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixATmul.C
index ac95f8c25..0cc879666 100644
--- a/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixATmul.C
+++ b/src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixATmul.C
@@ -237,7 +237,8 @@ void Foam::BlockLduMatrix<Type>::TmulCore
 
         for (register label coeffI = 0; coeffI < u.size(); coeffI++)
         {
-            Tx[u[coeffI]] += mult(activeUpper[coeffI], x[l[coeffI]]);
+            // Bug fix: Missing transpose. VV, 31/Aug/2015.
+            Tx[u[coeffI]] += mult(activeUpper[coeffI].T(), x[l[coeffI]]);
         }
     }
 
@@ -303,7 +304,8 @@ void Foam::BlockLduMatrix<Type>::TmulCore
 
             for (register label coeffI = 0; coeffI < u.size(); coeffI++)
             {
-                Tx[l[coeffI]] += mult(activeLower[coeffI], x[u[coeffI]]);
+                // Bug fix: Missing transpose. VV, 31/Aug/2015.
+                Tx[l[coeffI]] += mult(activeLower[coeffI].T(), x[u[coeffI]]);
             }
         }
     }