From 7a5c10e48da23c841dc424a0a755f8ce43ed38a7 Mon Sep 17 00:00:00 2001
From: Henrik Rusche <h.rusche@wikki-gmbh.de>
Date: Fri, 18 Mar 2011 12:25:39 +0100
Subject: [PATCH] Bug fix: harmonic interpolation @ coupled BC / optimisation

---
 .../regionCouple/regionCoupleFvPatchField.C   | 20 +++----
 .../schemes/harmonic/harmonic.H               | 52 +++++++++++++++----
 2 files changed, 52 insertions(+), 20 deletions(-)

diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/regionCouple/regionCoupleFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/constraint/regionCouple/regionCoupleFvPatchField.C
index e306ccf23..5ff86f070 100644
--- a/src/finiteVolume/fields/fvPatchFields/constraint/regionCouple/regionCoupleFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/constraint/regionCouple/regionCoupleFvPatchField.C
@@ -180,25 +180,27 @@ void regionCoupleFvPatchField<Type>::evaluate
     Field<Type> fOwn = this->patchInternalField();
     Field<Type> fNei = this->patchNeighbourField();
 
-    scalarField mOwn = mag(fOwn);
-    scalarField mNei = mag(fNei);
-    scalarField mean = 2*(mOwn*mNei)/(mOwn + mNei + SMALL);
+    scalarField magFOwn = mag(fOwn);
+    scalarField magFNei = mag(fNei);
 
-    scalarField weights(fOwn.size(), 0.5);
-
-    scalar den;
+    // Calculate internal weights using field magnitude
+    scalarField weights(fOwn.size());
 
     forAll (weights, faceI)
     {
-        den = mOwn[faceI] - mNei[faceI];
+        scalar mOwn = magFOwn[faceI];
+        scalar mNei = magFNei[faceI];
+
+        scalar den = mOwn - mNei;
 
         if (mag(den) > SMALL)
         {
-            weights[faceI] = (mean[faceI] - mNei[faceI])/den;
+            scalar mean = 2.0*mOwn*mNei/(mOwn + mNei);
+            weights[faceI] = (mean - mNei)/den;
         }
         else
         {
-            // Use 0.5 weights
+            weights[faceI] = 0.5;
         }
     }
 
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/harmonic/harmonic.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/harmonic/harmonic.H
index b7d2e5df5..e4c201a3a 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/harmonic/harmonic.H
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/harmonic/harmonic.H
@@ -131,33 +131,63 @@ public:
 
             scalarField magPhi = mag(phi);
 
-            // Initialise weights to 0.5 for uniform field (den = 0)
             scalarField& wIn = w.internalField();
-            wIn = 0.5;
 
             // Calculate internal weights using field magnitude
-            scalar mOwn, mNei, den, mean;
-
             forAll (owner, faceI)
             {
-                mOwn = magPhi[owner[faceI]];
-                mNei = magPhi[neighbour[faceI]];
+                scalar mOwn = magPhi[owner[faceI]];
+                scalar mNei = magPhi[neighbour[faceI]];
 
-                mean = 2*(mOwn*mNei)/(mOwn + mNei + SMALL);
-                den = mOwn - mNei;
+                scalar den = mOwn - mNei;
 
                 if (mag(den) > SMALL)
                 {
+                    scalar mean = 2.0*mOwn*mNei/(mOwn + mNei);
                     wIn[faceI] = (mean - mNei)/den;
                 }
                 else
                 {
-                    // Use 0.5 weights
+                    wIn[faceI] = 0.5;
                 }
             }
 
-            // Boundary weights are 1
-            w.boundaryField() = 1;
+            forAll (phi.boundaryField(), pi)
+            {
+                fvsPatchScalarField& wp = w.boundaryField()[pi];
+
+                const fvPatchField<Type>& patchPhi = phi.boundaryField()[pi];
+
+                if (patchPhi.coupled())
+                {
+                    scalarField magPhiOwn = mag(patchPhi.patchInternalField());
+                    scalarField magPhiNei = mag(patchPhi.patchNeighbourField());
+
+                    // Calculate internal weights using field magnitude
+                    forAll (patchPhi, faceI)
+                    {
+                        scalar mOwn = magPhiOwn[faceI];
+                        scalar mNei = magPhiNei[faceI];
+
+                        scalar den = mOwn - mNei;
+
+                        if (mag(den) > SMALL)
+                        {
+                             scalar mean = 2.0*mOwn*mNei/(mOwn + mNei);
+                             wp[faceI] = (mean - mNei)/den;
+                        }
+                        else
+                        {
+                             wp[faceI] = 0.5;
+                        }
+                    }
+                }
+                else
+                {
+                    // Boundary weights for uncoupled patches are 1
+                    wp = 1;
+                }
+            }
 
             return tw;
         }