From 5cb721159dbfd39215c1741e3f31f2d37acdf555 Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Mon, 7 Sep 2015 11:40:00 +0100
Subject: [PATCH] Bugfix: Work-around for a cyclic interface bug for tensors

---
 .../constraint/cyclic/cyclicFvPatchField.H    | 35 +++++++++++++++++--
 .../cyclicLduInterfaceField.C                 | 20 ++++++++++-
 2 files changed, 52 insertions(+), 3 deletions(-)

diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.H
index d6b164635..1c66c4bba 100644
--- a/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.H
+++ b/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.H
@@ -38,6 +38,7 @@ SourceFiles
 #include "coupledFvPatchField.H"
 #include "cyclicLduInterfaceField.H"
 #include "cyclicFvPatch.H"
+#include "transformField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -159,11 +160,41 @@ public:
             //- Transform neighbour field
             virtual void transformCoupleField
             (
-                scalarField& f,
+                scalarField& pnf,
                 const direction cmpt
             ) const
             {
-                cyclicLduInterfaceField::transformCoupleField(f, cmpt);
+                // See comments in cyclicLduInterfaceField
+                // HJ, 3/Sep/2015
+//                 cyclicLduInterfaceField::transformCoupleField(pnf, cmpt);
+                if (doTransform())
+                {
+                    label sizeby2 = pnf.size()/2;
+
+                    Type powFwdTransform = transform
+                    (
+                        this->forwardT()[0],
+                        pTraits<Type>::one
+                    );
+
+                    scalar forwardScale =
+                        pow(component(powFwdTransform, cmpt), rank());
+
+                    Type powRevTransform = transform
+                    (
+                        this->reverseT()[0],
+                        pTraits<Type>::one
+                    );
+
+                    scalar reverseScale =
+                        pow(component(powRevTransform, cmpt), rank());
+
+                    for (label facei = 0; facei < sizeby2; facei++)
+                    {
+                        pnf[facei] *= forwardScale;
+                        pnf[facei + sizeby2] *= reverseScale;
+                    }
+                }
             }
 
             //- Update result field based on interface functionality
diff --git a/src/foam/matrices/lduMatrix/lduAddressing/lduInterfaceFields/cyclicLduInterfaceField/cyclicLduInterfaceField.C b/src/foam/matrices/lduMatrix/lduAddressing/lduInterfaceFields/cyclicLduInterfaceField/cyclicLduInterfaceField.C
index 338ae6b95..440e87d9e 100644
--- a/src/foam/matrices/lduMatrix/lduAddressing/lduInterfaceFields/cyclicLduInterfaceField/cyclicLduInterfaceField.C
+++ b/src/foam/matrices/lduMatrix/lduAddressing/lduInterfaceFields/cyclicLduInterfaceField/cyclicLduInterfaceField.C
@@ -52,13 +52,31 @@ void Foam::cyclicLduInterfaceField::transformCoupleField
     {
         label sizeby2 = pnf.size()/2;
 
+        //HJ This is wrong
+        // For a vector, the implicit contribution is correct:
+        // power of diagonal of the transformation tensor for the vector
+        // For a tensor this is wrong: the diag of transformation tensor
+        // is a vector and I am asking for a component cmpt which corresponds
+        // to the component of the variable I am offering.  If the variable
+        // is a tensor, it has got 6 or 9 components, so the component
+        // access is out of range and rubbish anyway
+        // The forwardScale/reverseScale should be the power of the complete
+        // transformation matrix for a tensor.  Something like:
+
+//         Type powTransform = transform(pTraits<Type>::one);
+//         forwardScale = pow(powTransform.component(cmpt), rank());
+
+        // This needs to be moved under templating into cyclicFvPatchField
+        // because it requires access to Type
+        // HJ, 3/Sep/2015
+
         scalar forwardScale =
             pow(diag(forwardT()[0]).component(cmpt), rank());
 
         scalar reverseScale =
             pow(diag(reverseT()[0]).component(cmpt), rank());
 
-        for (label facei=0; facei<sizeby2; facei++)
+        for (label facei = 0; facei < sizeby2; facei++)
         {
             pnf[facei] *= forwardScale;
             pnf[facei + sizeby2] *= reverseScale;