From ae6b76038fb7bdfc6deca95ede29c3835217f791 Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Wed, 14 Aug 2013 12:54:34 +0100
Subject: [PATCH] Block multiplication bug fix

---
 .../BlockLduMatrix/BlockLduMatrixATmul.C      | 69 ++-----------------
 1 file changed, 6 insertions(+), 63 deletions(-)

diff --git a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixATmul.C b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixATmul.C
index 661dfe335..2d8c79231 100644
--- a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixATmul.C
+++ b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixATmul.C
@@ -64,44 +64,15 @@ void Foam::BlockLduMatrix<Type>::AmulCore
     const unallocLabelList& u = lduAddr().upperAddr();
     const unallocLabelList& l = lduAddr().lowerAddr();
 
-
     const TypeCoeffField& Diag = this->diag();
     const TypeCoeffField& Upper = this->upper();
 
+    // Diagonal multiplication, no indirection
+    multiply(Ax, Diag, x);
+
     // Create multiplication function object
     typename BlockCoeff<Type>::multiply mult;
 
-    // AmulCore must be additive to account for initialisation step
-    // in ldu interfaces.  HJ, 6/Nov/2007
-    // Fixed by IC 19/Oct/2011
-    if (Diag.activeType() == blockCoeffBase::SCALAR)
-    {
-        const scalarTypeField& activeDiag = Diag.asScalar();
-
-        for (register label cellI = 0; cellI < u.size(); cellI++)
-        {
-            Ax[cellI] += mult(activeDiag[cellI], x[cellI]);
-        }
-    }
-    else if (Diag.activeType() == blockCoeffBase::LINEAR)
-    {
-        const linearTypeField& activeDiag = Diag.asLinear();
-
-        for (register label cellI = 0; cellI < u.size(); cellI++)
-        {
-            Ax[cellI] += mult(activeDiag[cellI], x[cellI]);
-        }
-    }
-    else if (Diag.activeType() == blockCoeffBase::SQUARE)
-    {
-        const squareTypeField& activeDiag = Diag.asSquare();
-
-        for (register label cellI = 0; cellI < u.size(); cellI++)
-        {
-            Ax[cellI] += mult(activeDiag[cellI], x[cellI]);
-        }
-    }
-
     // Lower multiplication
 
     if (symmetric())
@@ -238,40 +209,12 @@ void Foam::BlockLduMatrix<Type>::TmulCore
     const TypeCoeffField& Diag = this->diag();
     const TypeCoeffField& Upper = this->upper();
 
+    // Diagonal multiplication, no indirection
+    multiply(Tx, Diag, x);
+
     // Create multiplication function object
     typename BlockCoeff<Type>::multiply mult;
 
-    // AmulCore must be additive to account for initialisation step
-    // in ldu interfaces.  HJ, 6/Nov/2007
-    // Fixed by IC 19/Oct/2011
-    if (Diag.activeType() == blockCoeffBase::SCALAR)
-    {
-        const scalarTypeField& activeDiag = Diag.asScalar();
-
-        for (register label cellI = 0; cellI < u.size(); cellI++)
-        {
-            Tx[cellI] += mult(activeDiag[cellI], x[cellI]);
-        }
-    }
-    else if (Diag.activeType() == blockCoeffBase::LINEAR)
-    {
-        const linearTypeField& activeDiag = Diag.asLinear();
-
-        for (register label cellI = 0; cellI < u.size(); cellI++)
-        {
-            Tx[cellI] += mult(activeDiag[cellI], x[cellI]);
-        }
-    }
-    else if (Diag.activeType() == blockCoeffBase::SQUARE)
-    {
-        const squareTypeField& activeDiag = Diag.asSquare();
-
-        for (register label cellI = 0; cellI < u.size(); cellI++)
-        {
-            Tx[cellI] += mult(activeDiag[cellI].T(), x[cellI]);
-        }
-    }
-
     // Upper multiplication
 
     if (Upper.activeType() == blockCoeffBase::SCALAR)