From 5a13369b6d6c7ef6f335353e1f681b7af3cf35c0 Mon Sep 17 00:00:00 2001
From: Joakim Hove <joakim.hove@gmail.com>
Date: Tue, 10 Jun 2014 19:52:22 +0200
Subject: [PATCH] Added EclipseGrid::getCellVolume() and getCellCenter()

---
 .../eclipse/EclipseState/Grid/EclipseGrid.cpp | 44 ++++++++++++++++++-
 .../eclipse/EclipseState/Grid/EclipseGrid.hpp |  7 +++
 .../Grid/tests/EclipseGridTests.cpp           | 34 +++++++++++++-
 3 files changed, 83 insertions(+), 2 deletions(-)

diff --git a/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.cpp b/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.cpp
index d49316869..c517149cd 100644
--- a/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.cpp
+++ b/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.cpp
@@ -19,6 +19,7 @@
 
 
 #include <iostream>
+#include <tuple>
 
 #include <boost/lexical_cast.hpp>
 
@@ -73,8 +74,49 @@ namespace Opm {
     size_t EclipseGrid::getCartesianSize( ) const {
         return static_cast<size_t>( ecl_grid_get_global_size( m_grid.get() ));
     }
-
     
+    void EclipseGrid::assertGlobalIndex(size_t globalIndex) const {
+        if (globalIndex >= getCartesianSize())
+            throw std::invalid_argument("input index above valid range");
+    }
+
+    void EclipseGrid::assertIJK(size_t i , size_t j , size_t k) const {
+        if (i >= getNX() || j >= getNY() || k >= getNZ())
+            throw std::invalid_argument("input index above valid range");
+    }
+
+
+    double EclipseGrid::getCellVolume(size_t globalIndex) const {
+        assertGlobalIndex( globalIndex );
+        return ecl_grid_get_cell_volume1( m_grid.get() , static_cast<int>(globalIndex));
+    }
+
+
+    double EclipseGrid::getCellVolume(size_t i , size_t j , size_t k) const {
+        assertIJK(i,j,k);
+        return ecl_grid_get_cell_volume3( m_grid.get() , static_cast<int>(i),static_cast<int>(j),static_cast<int>(k));
+    }
+
+    std::tuple<double,double,double> EclipseGrid::getCellCenter(size_t globalIndex) const {
+        assertGlobalIndex( globalIndex );
+        {
+            double x,y,z;
+            ecl_grid_get_xyz1( m_grid.get() , static_cast<int>(globalIndex) , &x , &y , &z);
+            return std::tuple<double,double,double> {x,y,z};
+        }
+    }
+
+
+    std::tuple<double,double,double> EclipseGrid::getCellCenter(size_t i,size_t j, size_t k) const {
+        assertIJK(i,j,k);
+        {
+            double x,y,z;
+            ecl_grid_get_xyz3( m_grid.get() , static_cast<int>(i),static_cast<int>(j),static_cast<int>(k), &x , &y , &z);
+            return std::tuple<double,double,double> {x,y,z};
+        }
+    }
+
+
     bool EclipseGrid::hasCornerPointKeywords(std::shared_ptr<const GRIDSection> gridSection) {
         if (gridSection->hasKeyword("ZCORN") && gridSection->hasKeyword("COORD"))
             return true;
diff --git a/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp b/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp
index b594ed949..59b834e33 100644
--- a/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp
+++ b/opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp
@@ -43,6 +43,13 @@ namespace Opm {
         size_t  getNZ( ) const;
         size_t  getCartesianSize( ) const;
 
+        void assertGlobalIndex(size_t globalIndex) const;
+        void assertIJK(size_t i , size_t j , size_t k) const;
+        std::tuple<double,double,double> getCellCenter(size_t i,size_t j, size_t k) const;
+        std::tuple<double,double,double> getCellCenter(size_t globalIndex) const;
+        double getCellVolume(size_t globalIndex) const;
+        double getCellVolume(size_t i , size_t j , size_t k) const;
+
         void exportMAPAXES( std::vector<double>& mapaxes) const;
         void exportCOORD( std::vector<double>& coord) const;
         void exportZCORN( std::vector<double>& zcorn) const;
diff --git a/opm/parser/eclipse/EclipseState/Grid/tests/EclipseGridTests.cpp b/opm/parser/eclipse/EclipseState/Grid/tests/EclipseGridTests.cpp
index af2603a22..dca409f81 100644
--- a/opm/parser/eclipse/EclipseState/Grid/tests/EclipseGridTests.cpp
+++ b/opm/parser/eclipse/EclipseState/Grid/tests/EclipseGridTests.cpp
@@ -170,7 +170,39 @@ BOOST_AUTO_TEST_CASE(DEPTHZ_EQUAL_TOPS) {
     std::shared_ptr<Opm::EclipseGrid> grid2(new Opm::EclipseGrid( runspecSection , gridSection2 ));
 
     BOOST_CHECK( grid1->equal( *(grid2.get()) ));
-    
+
+    {
+        BOOST_CHECK_THROW( grid1->getCellVolume(1000) , std::invalid_argument);
+        BOOST_CHECK_THROW( grid1->getCellVolume(10,0,0) , std::invalid_argument);
+        BOOST_CHECK_THROW( grid1->getCellVolume(0,10,0) , std::invalid_argument);
+        BOOST_CHECK_THROW( grid1->getCellVolume(0,0,10) , std::invalid_argument);
+        
+        for (size_t g=0; g < 1000; g++)
+            BOOST_CHECK_CLOSE( grid1->getCellVolume(g) , 0.25*0.25*0.25 , 0.001);
+        
+        
+        for (size_t k= 0; k < 10; k++)
+            for (size_t j= 0; j < 10; j++)  
+                for (size_t i= 0; i < 10; i++)
+                    BOOST_CHECK_CLOSE( grid1->getCellVolume(i,j,k) , 0.25*0.25*0.25 , 0.001 );
+    }
+    {
+        BOOST_CHECK_THROW( grid1->getCellCenter(1000) , std::invalid_argument);
+        BOOST_CHECK_THROW( grid1->getCellCenter(10,0,0) , std::invalid_argument);
+        BOOST_CHECK_THROW( grid1->getCellCenter(0,10,0) , std::invalid_argument);
+        BOOST_CHECK_THROW( grid1->getCellCenter(0,0,10) , std::invalid_argument);
+
+        for (size_t k= 0; k < 10; k++)
+            for (size_t j= 0; j < 10; j++)  
+                for (size_t i= 0; i < 10; i++) {
+                    auto pos = grid1->getCellCenter(i,j,k);
+                    
+                    BOOST_CHECK_CLOSE( std::get<0>(pos) , i*0.25 + 0.125, 0.001);
+                    BOOST_CHECK_CLOSE( std::get<1>(pos) , j*0.25 + 0.125, 0.001);
+                    BOOST_CHECK_CLOSE( std::get<2>(pos) , k*0.25 + 0.125 + 0.25, 0.001);
+                    
+                }
+    }
 }