From ef9966b188d64cce9611e87bb72a136d662399f8 Mon Sep 17 00:00:00 2001
From: Tor Harald Sandve <tor.harald.sandve@iris.no>
Date: Thu, 19 Jan 2017 08:57:44 +0100
Subject: [PATCH] Start using face geometry computed the Ecl way

Face centers are computed using the cell corners. With this
implementation the face center seen from a cell may be different from the
face center seen from its neighbour.

Face normals with area lenghts are calculated using the face corners
directly not using a triangulation point in the center of the face.

This gives transmissibility almost equal to eclipse.
---
 ebos/ecltransmissibility.hh | 28 ++++++++++++++--------------
 1 file changed, 14 insertions(+), 14 deletions(-)

diff --git a/ebos/ecltransmissibility.hh b/ebos/ecltransmissibility.hh
index d15c40ba3..b5a68bd0c 100644
--- a/ebos/ecltransmissibility.hh
+++ b/ebos/ecltransmissibility.hh
@@ -184,21 +184,26 @@ public:
                 unsigned insideFaceIdx  = intersection.indexInInside();
                 unsigned outsideFaceIdx = intersection.indexInOutside();
 
+                int faceIdx = intersection.id();
+                DimVector faceCenterInside = gridManager_.grid().faceCenterEcl(insideElemIdx,insideFaceIdx);
+                DimVector faceCenterOutside = gridManager_.grid().faceCenterEcl(outsideElemIdx,outsideFaceIdx);
+                DimVector faceAreaNormal = gridManager_.grid().faceAreaNormalEcl(faceIdx);
+
                 Scalar halfTrans1;
                 Scalar halfTrans2;
 
                 computeHalfTrans_(halfTrans1,
-                                  intersection,
+                                  faceAreaNormal,
                                   insideFaceIdx,
-                                  distanceVector_(intersection,
+                                  distanceVector_(faceCenterInside,
                                                   intersection.indexInInside(),
                                                   insideElemIdx,
                                                   axisCentroids),
                                   permeability_[insideElemIdx]);
                 computeHalfTrans_(halfTrans2,
-                                  intersection,
+                                  faceAreaNormal,
                                   outsideFaceIdx,
-                                  distanceVector_(intersection,
+                                  distanceVector_(faceCenterOutside,
                                                   intersection.indexInOutside(),
                                                   outsideElemIdx,
                                                   axisCentroids),
@@ -308,36 +313,31 @@ private:
     }
 
     void computeHalfTrans_(Scalar& halfTrans,
-                           const Intersection& is,
+                           const DimVector& areaNormal,
                            unsigned faceIdx, // in the reference element that contains the intersection
                            const DimVector& distance,
                            const DimMatrix& perm) const
     {
-        Scalar isArea = is.geometry().volume();
-        DimVector n = is.centerUnitOuterNormal();
-        n *= isArea;
-
         unsigned dimIdx = faceIdx/2;
         assert(dimIdx < dimWorld);
         halfTrans = perm[dimIdx][dimIdx];
-        //halfTrans *= isArea;
 
         Scalar val = 0;
-        for (unsigned i = 0; i < n.size(); ++i)
-            val += n[i]*distance[i];
+        for (unsigned i = 0; i < areaNormal.size(); ++i)
+            val += areaNormal[i]*distance[i];
 
         halfTrans *= std::abs<Scalar>(val);
         halfTrans /= distance.two_norm2();
     }
 
-    DimVector distanceVector_(const Intersection& is,
+    DimVector distanceVector_(const DimVector& center,
                               unsigned faceIdx, // in the reference element that contains the intersection
                               unsigned elemIdx,
                               const std::array<std::vector<DimVector>, dimWorld>& axisCentroids) const
     {
         unsigned dimIdx = faceIdx/2;
         assert(dimIdx < dimWorld);
-        DimVector x = is.geometry().center();
+        DimVector x = center;
         x -= axisCentroids[dimIdx][elemIdx];
 
         return x;