From e4f345d1249924bd36a83e0fe0e237b080722af3 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Thu, 19 Jan 2017 08:27:01 +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. --- opm/autodiff/GeoProps.hpp | 23 ++++++----------------- 1 file changed, 6 insertions(+), 17 deletions(-) diff --git a/opm/autodiff/GeoProps.hpp b/opm/autodiff/GeoProps.hpp index d769306e5..d312a23b3 100644 --- a/opm/autodiff/GeoProps.hpp +++ b/opm/autodiff/GeoProps.hpp @@ -615,27 +615,15 @@ namespace Opm double sgn = 2.0 * (faceCells(faceIdx, 0) == cellIdx) - 1; const int dim = Opm::UgGridHelpers::dimensions(grid); - const double* faceNormal = Opm::UgGridHelpers::faceNormal(grid, faceIdx); -#if HAVE_OPM_GRID - assert( dim <= 3 ); - Dune::FieldVector< double, 3 > scaledFaceNormal( 0 ); - for (int indx = 0; indx < dim; ++indx) { - scaledFaceNormal[ indx ] = faceNormal[ indx ]; - } - // compute unit normal incase the normal is already scaled - scaledFaceNormal /= scaledFaceNormal.two_norm(); - // compute proper normal scaled with face area - scaledFaceNormal *= Opm::UgGridHelpers::faceArea(grid, faceIdx); -#else - const double* scaledFaceNormal = faceNormal; -#endif - int cartesianCellIdx = AutoDiffGrid::globalCell(grid)[cellIdx]; auto cellCenter = eclGrid.getCellCenter(cartesianCellIdx); + const auto& faceCenter = Opm::UgGridHelpers::faceCenterEcl(grid, cellIdx, faceTag); + const auto& faceAreaNormalEcl = Opm::UgGridHelpers::faceAreaNormalEcl(grid, faceIdx); + for (int indx = 0; indx < dim; ++indx) { - const double Ci = Opm::UgGridHelpers::faceCentroid(grid, faceIdx)[indx] - cellCenter[indx]; + const double Ci = faceCenter[indx] - cellCenter[indx]; dist += Ci*Ci; - cn += sgn * Ci * scaledFaceNormal[ indx ]; + cn += sgn * Ci * faceAreaNormalEcl[ indx ]; } if (cn < 0){ @@ -656,6 +644,7 @@ namespace Opm cn = -cn; } hTrans[cellFaceIdx] = perm[cellIdx*dim*dim + d] * cn / dist; + } }