From aeb8807d890b2de37098da9f9b380e91fd01d108 Mon Sep 17 00:00:00 2001
From: Tor Harald Sandve <tor.harald.sandve@iris.no>
Date: Tue, 24 Nov 2015 10:01:43 +0100
Subject: [PATCH] Use cell centroid from EclipseGrid

The cell center in EclipseGrid is computed based on averaging the
corners of the cells as in Eclipse and not computing the centroid of the
cell as in the unstructured grid in opm.
---
 opm/autodiff/GeoProps.hpp | 26 ++++++++++++++++++++++----
 1 file changed, 22 insertions(+), 4 deletions(-)

diff --git a/opm/autodiff/GeoProps.hpp b/opm/autodiff/GeoProps.hpp
index 3e1722f2f..c10088a3a 100644
--- a/opm/autodiff/GeoProps.hpp
+++ b/opm/autodiff/GeoProps.hpp
@@ -125,7 +125,7 @@ namespace Opm
                 tpfa_htrans_compute(ug, props.permeability(), htrans.data());
             }
             else {
-                tpfa_loc_trans_compute_(grid,props.permeability(),htrans);
+                tpfa_loc_trans_compute_(grid,eclgrid, props.permeability(),htrans);
             }
 
             std::vector<double> mult;
@@ -191,6 +191,7 @@ namespace Opm
 
         template <class Grid>
         void tpfa_loc_trans_compute_(const Grid &grid,
+                                     Opm::EclipseGridConstPtr eclGrid,
                                      const double* perm,
                                      Vector &hTrans);
 
@@ -338,6 +339,7 @@ namespace Opm
 
     template <class GridType>
     inline void DerivedGeology::tpfa_loc_trans_compute_(const GridType& grid,
+                                                        Opm::EclipseGridConstPtr eclGrid,
                                                         const double* perm,
                                                         Vector& hTrans){
 
@@ -358,6 +360,7 @@ namespace Opm
             for(auto cellFaceIter = cellFacesRange.begin(), cellFaceEnd = cellFacesRange.end();
                 cellFaceIter != cellFaceEnd; ++cellFaceIter, ++cellFaceIdx)
             {
+
                 // The index of the face in the compressed grid
                 const int faceIdx = *cellFaceIter;
 
@@ -388,11 +391,26 @@ namespace Opm
                 const double* scaledFaceNormal = faceNormal;
 #endif
 
+                int cartesianCellIdx = AutoDiffGrid::globalCell(grid)[cellIdx];
+                auto cellCentroid = eclGrid->getCellCenter(cartesianCellIdx);
+
                 for (int indx = 0; indx < dim; ++indx) {
-                    const double Ci = Opm::UgGridHelpers::faceCentroid(grid, faceIdx)[indx] -
-                        Opm::UgGridHelpers::cellCentroidCoordinate(grid, cellIdx, indx);
+
+                    double Ci = Opm::UgGridHelpers::faceCentroid(grid, faceIdx)[indx];
+                    switch (indx) {
+                    case 0:
+                        Ci -= std::get<0>(cellCentroid);
+                        break;
+                    case 1:
+                        Ci -= std::get<1>(cellCentroid);
+                        break;
+                    case 2:
+                        Ci -= std::get<2>(cellCentroid);
+                        break;
+
+                    }
                     dist += Ci*Ci;
-                    cn += sgn * Ci * scaledFaceNormal[ indx ]; //Opm::UgGridHelpers::faceNormal(grid, faceIdx)[indx];
+                    cn += sgn * Ci * scaledFaceNormal[ indx ];
                 }
 
                 if (cn < 0){