From 297f63cce2eddaa453f3012a90edd5f1f8db7e2d Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Wed, 11 Jan 2017 11:55:37 +0000
Subject: [PATCH] Bugfix: consistency of finite area operators

---
 src/finiteArea/finiteArea/fac/facGrad.C  |  4 +--
 src/finiteArea/finiteArea/fac/facNDiv.C  | 38 ++++++++++++++----------
 src/finiteArea/finiteArea/fac/facNGrad.C | 24 ++++++++++-----
 src/finiteArea/finiteArea/fam/famDiv.C   |  4 ++-
 src/finiteArea/finiteArea/fam/famNDiv.C  |  2 ++
 5 files changed, 45 insertions(+), 27 deletions(-)

diff --git a/src/finiteArea/finiteArea/fac/facGrad.C b/src/finiteArea/finiteArea/fac/facGrad.C
index 568f4c2b6..d5b117638 100644
--- a/src/finiteArea/finiteArea/fac/facGrad.C
+++ b/src/finiteArea/finiteArea/fac/facGrad.C
@@ -63,7 +63,7 @@ grad
 
     GeometricField<GradType, faPatchField, areaMesh>& gGrad = tgGrad();
 
-    gGrad -= (gGrad & n)*n;
+    gGrad -= n*(n & gGrad);
     gGrad.correctBoundaryConditions();
 
     return tgGrad;
@@ -118,7 +118,7 @@ grad
 
     GeometricField<GradType, faPatchField, areaMesh>& gGrad = tgGrad();
 
-    gGrad -= (gGrad & n)*n;
+    gGrad -= n*(n & gGrad);
     gGrad.correctBoundaryConditions();
 
     return tgGrad;
diff --git a/src/finiteArea/finiteArea/fac/facNDiv.C b/src/finiteArea/finiteArea/fac/facNDiv.C
index 6a1953961..67553a50a 100644
--- a/src/finiteArea/finiteArea/fac/facNDiv.C
+++ b/src/finiteArea/finiteArea/fac/facNDiv.C
@@ -49,12 +49,13 @@ ndiv
     const GeometricField<Type, faePatchField, edgeMesh>& ssf
 )
 {
-    const areaVectorField &n = ssf.mesh().faceAreaNormals();
+    const areaVectorField& n = ssf.mesh().faceAreaNormals();
 
-    tmp<GeometricField<Type, faPatchField, areaMesh> > v = fac::edgeIntegrate(ssf);
+    tmp<GeometricField<Type, faPatchField, areaMesh> > v =
+        fac::edgeIntegrate(ssf);
 
-    //v.internalField() = transform(n*n, v.internalField());
-    v.internalField() = (v.internalField()&n)*n;
+    v.internalField() = n*(n & v.internalField());
+    v.correctBoundaryConditions();
 
     return v;
 }
@@ -87,7 +88,7 @@ ndiv
     const word& name
 )
 {
-    const areaVectorField &n = vf.mesh().faceAreaNormals();
+    const areaVectorField& n = vf.mesh().faceAreaNormals();
 
     tmp<GeometricField<Type, faPatchField, areaMesh> > tDiv
     (
@@ -96,11 +97,11 @@ ndiv
                 vf.mesh(), vf.mesh().schemesDict().divScheme(name)
             )().facDiv(vf)
     );
-
     GeometricField<Type, faPatchField, areaMesh>& Div = tDiv();
 
-    //Div.internalField() = transform(n*n, Div.internalField());
-    Div.internalField() = (Div.internalField()&n)*n;
+    Div.internalField() = n*(n & Div.internalField());
+    Div.correctBoundaryConditions();
+
     return tDiv;
 }
 
@@ -125,6 +126,7 @@ ndiv
         fac::ndiv(tvvf(), name)
     );
     tvvf.clear();
+
     return Div;
 }
 
@@ -178,22 +180,22 @@ ndiv
     const word& name
 )
 {
-    const areaVectorField &n = vf.mesh().faceAreaNormals();
+    const areaVectorField& n = vf.mesh().faceAreaNormals();
 
     tmp<GeometricField<Type, faPatchField, areaMesh> > tDiv
     (
         fa::convectionScheme<Type>::New
-            (
-                vf.mesh(),
-                flux,
-                vf.mesh().schemesDict().divScheme(name)
-            )().facDiv(flux, vf)
+        (
+            vf.mesh(),
+            flux,
+            vf.mesh().schemesDict().divScheme(name)
+        )().facDiv(flux, vf)
     );
 
     GeometricField<Type, faPatchField, areaMesh>& Div = tDiv();
 
-    //Div.internalField() = transform(n*n, Div.internalField());
-    Div.internalField() = (Div.internalField()&n)*n;
+    Div.internalField() = n*(n &Div.internalField());
+    Div.correctBoundaryConditions();
 
     return tDiv;
 
@@ -214,6 +216,7 @@ ndiv
         fac::ndiv(tflux(), vf, name)
     );
     tflux.clear();
+
     return Div;
 }
 
@@ -232,6 +235,7 @@ ndiv
         fac::ndiv(flux, tvf(), name)
     );
     tvf.clear();
+
     return Div;
 }
 
@@ -251,6 +255,7 @@ ndiv
     );
     tflux.clear();
     tvf.clear();
+
     return Div;
 }
 
@@ -318,6 +323,7 @@ ndiv
     );
     tflux.clear();
     tvf.clear();
+
     return Div;
 }
 
diff --git a/src/finiteArea/finiteArea/fac/facNGrad.C b/src/finiteArea/finiteArea/fac/facNGrad.C
index 9d4d96070..8c1df7809 100644
--- a/src/finiteArea/finiteArea/fac/facNGrad.C
+++ b/src/finiteArea/finiteArea/fac/facNGrad.C
@@ -56,18 +56,21 @@ ngrad
     const GeometricField<Type, faePatchField, edgeMesh>& ssf
 )
 {
-    const areaVectorField &n = ssf.mesh().faceAreaNormals();
+    const areaVectorField& n = ssf.mesh().faceAreaNormals();
     typedef typename outerProduct<vector,Type>::type GradType;
 
-    tmp<GeometricField<GradType, faPatchField, areaMesh> > tgGrad = fac::edgeIntegrate(ssf.mesh().Sf() * ssf);
+    tmp<GeometricField<GradType, faPatchField, areaMesh> > tgGrad =
+        fac::edgeIntegrate(ssf.mesh().Sf() * ssf);
 
-    GeometricField<GradType, faPatchField, areaMesh>& grad = tgGrad();
+    GeometricField<GradType, faPatchField, areaMesh>& gGrad = tgGrad();
 
-    grad = (grad&n)*n;
+    gGrad = n*(n & gGrad);
+    gGrad.correctBoundaryConditions();
 
     return tgGrad;
 }
 
+
 template<class Type>
 tmp
 <
@@ -87,6 +90,7 @@ ngrad
         fac::ngrad(tssf())
     );
     tssf.clear();
+
     return Grad;
 }
 
@@ -105,18 +109,20 @@ ngrad
     const word& name
 )
 {
-    const areaVectorField &n = vf.mesh().faceAreaNormals();
+    const areaVectorField& n = vf.mesh().faceAreaNormals();
     typedef typename outerProduct<vector,Type>::type GradType;
 
-    tmp<GeometricField<GradType, faPatchField, areaMesh> > tgGrad = fa::gradScheme<Type>::New
+    tmp<GeometricField<GradType, faPatchField, areaMesh> > tgGrad =
+        fa::gradScheme<Type>::New
         (
             vf.mesh(),
             vf.mesh().schemesDict().gradScheme(name)
         )().grad(vf);
 
-    GeometricField<GradType, faPatchField, areaMesh>& grad = tgGrad();
+    GeometricField<GradType, faPatchField, areaMesh>& gGrad = tgGrad();
 
-    grad = (grad&n)*n;
+    gGrad = n*(n & gGrad);
+    gGrad.correctBoundaryConditions();
 
     return tgGrad;
 }
@@ -147,6 +153,7 @@ ngrad
         fac::ngrad(tvf(), name)
     );
     tvf.clear();
+
     return tGrad;
 }
 
@@ -187,6 +194,7 @@ ngrad
         fac::ngrad(tvf())
     );
     tvf.clear();
+
     return Grad;
 }
 
diff --git a/src/finiteArea/finiteArea/fam/famDiv.C b/src/finiteArea/finiteArea/fam/famDiv.C
index 2296b224b..0b2af47a7 100644
--- a/src/finiteArea/finiteArea/fam/famDiv.C
+++ b/src/finiteArea/finiteArea/fam/famDiv.C
@@ -73,7 +73,7 @@ div
     );
 
     //HJ Check if the product is from left or right.  HJ, 6/Dec/2016
-    M -= (v & n)*n;
+    M -= n*(n & v);
 
     return tM;
 }
@@ -90,6 +90,7 @@ div
 {
     tmp<faMatrix<Type> > Div(fam::div(tflux(), vf, name));
     tflux.clear();
+
     return Div;
 }
 
@@ -115,6 +116,7 @@ div
 {
     tmp<faMatrix<Type> > Div(fam::div(tflux(), vf));
     tflux.clear();
+
     return Div;
 }
 
diff --git a/src/finiteArea/finiteArea/fam/famNDiv.C b/src/finiteArea/finiteArea/fam/famNDiv.C
index 33c55f3fa..f668bb626 100644
--- a/src/finiteArea/finiteArea/fam/famNDiv.C
+++ b/src/finiteArea/finiteArea/fam/famNDiv.C
@@ -68,6 +68,7 @@ ndiv
 {
     tmp<faMatrix<Type> > Div(fam::ndiv(tflux(), vf, name));
     tflux.clear();
+
     return Div;
 }
 
@@ -93,6 +94,7 @@ ndiv
 {
     tmp<faMatrix<Type> > Div(fam::ndiv(tflux(), vf));
     tflux.clear();
+
     return Div;
 }