diff --git a/opm/core/eclipse/CornerpointChopper.hpp b/opm/core/eclipse/CornerpointChopper.hpp
new file mode 100644
index 00000000..62a6a05c
--- /dev/null
+++ b/opm/core/eclipse/CornerpointChopper.hpp
@@ -0,0 +1,420 @@
+/*
+ Copyright 2010 SINTEF ICT, Applied Mathematics.
+
+ This file is part of the Open Porous Media project (OPM).
+
+ OPM is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OPM is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OPM. If not, see .
+*/
+
+#ifndef OPM_CORNERPOINTCHOPPER_HEADER_INCLUDED
+#define OPM_CORNERPOINTCHOPPER_HEADER_INCLUDED
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+namespace Dune
+{
+
+ class CornerPointChopper
+ {
+ public:
+ CornerPointChopper(const std::string& file)
+ : parser_(file, false)
+ {
+ for (int dd = 0; dd < 3; ++dd) {
+ dims_[dd] = parser_.getSPECGRID().dimensions[dd];
+ }
+
+ int layersz = 8*dims_[0]*dims_[1];
+ const std::vector& ZCORN = parser_.getFloatingPointValue("ZCORN");
+ botmax_ = *std::max_element(ZCORN.begin(), ZCORN.begin() + layersz/2);
+ topmin_ = *std::min_element(ZCORN.begin() + dims_[2]*layersz - layersz/2,
+ ZCORN.begin() + dims_[2]*layersz);
+
+ abszmax_ = *std::max_element(ZCORN.begin(), ZCORN.end());
+ abszmin_ = *std::min_element(ZCORN.begin(), ZCORN.end());
+
+ std::cout << "Parsed grdecl file with dimensions ("
+ << dims_[0] << ", " << dims_[1] << ", " << dims_[2] << ")" << std::endl;
+ }
+
+
+
+
+ const int* dimensions() const
+ {
+ return dims_;
+ }
+
+
+
+
+ const int* newDimensions() const
+ {
+ return new_dims_;
+ }
+
+
+
+
+ const std::pair zLimits() const
+ {
+ return std::make_pair(botmax_, topmin_);
+ }
+
+ const std::pair abszLimits() const
+ {
+ return std::make_pair(abszmin_, abszmax_);
+ }
+
+
+ void verifyInscribedShoebox(int imin, int ilen, int imax,
+ int jmin, int jlen, int jmax,
+ double zmin, double zlen, double zmax)
+ {
+ if (imin < 0) {
+ std::cerr << "Error! imin < 0 (imin = " << imin << ")\n";
+ throw std::runtime_error("Inconsistent user input.");
+ }
+ if (ilen > dims_[0]) {
+ std::cerr << "Error! ilen larger than grid (ilen = " << ilen <<")\n";
+ throw std::runtime_error("Inconsistent user input.");
+ }
+ if (imax > dims_[0]) {
+ std::cerr << "Error! imax larger than input grid (imax = " << imax << ")\n";
+ throw std::runtime_error("Inconsistent user input.");
+ }
+ if (jmin < 0) {
+ std::cerr << "Error! jmin < 0 (jmin = " << jmin << ")\n";
+ throw std::runtime_error("Inconsistent user input.");
+ }
+ if (jlen > dims_[1]) {
+ std::cerr << "Error! jlen larger than grid (jlen = " << jlen <<")\n";
+ throw std::runtime_error("Inconsistent user input.");
+ }
+ if (jmax > dims_[1]) {
+ std::cerr << "Error! jmax larger than input grid (jmax = " << jmax << ")\n";
+ throw std::runtime_error("Inconsistent user input.");
+ }
+ if (zmin < abszmin_) {
+ std::cerr << "Error! zmin ("<< zmin << ") less than minimum ZCORN value ("<< abszmin_ << ")\n";
+ throw std::runtime_error("Inconsistent user input.");
+ }
+ if (zmax > abszmax_) {
+ std::cerr << "Error! zmax ("<< zmax << ") larger than maximal ZCORN value ("<< abszmax_ << ")\n";
+ throw std::runtime_error("Inconsistent user input.");
+ }
+ if (zlen > (abszmax_ - abszmin_)) {
+ std::cerr << "Error! zlen ("<< zlen <<") larger than maximal ZCORN (" << abszmax_ << ") minus minimal ZCORN ("<< abszmin_ <<")\n";
+ throw std::runtime_error("Inconsistent user input.");
+ }
+ }
+
+ void chop(int imin, int imax, int jmin, int jmax, double zmin, double zmax, bool resettoorigin=true)
+ {
+ new_dims_[0] = imax - imin;
+ new_dims_[1] = jmax - jmin;
+
+ // Filter the coord field
+ const std::vector& COORD = parser_.getFloatingPointValue("COORD");
+ int num_coord = COORD.size();
+ if (num_coord != 6*(dims_[0] + 1)*(dims_[1] + 1)) {
+ std::cerr << "Error! COORD size (" << COORD.size() << ") not consistent with SPECGRID\n";
+ throw std::runtime_error("Inconsistent COORD and SPECGRID.");
+ }
+ int num_new_coord = 6*(new_dims_[0] + 1)*(new_dims_[1] + 1);
+ double x_correction = COORD[6*((dims_[0] + 1)*jmin + imin)];
+ double y_correction = COORD[6*((dims_[0] + 1)*jmin + imin) + 1];
+ new_COORD_.resize(num_new_coord, 1e100);
+ for (int j = jmin; j < jmax + 1; ++j) {
+ for (int i = imin; i < imax + 1; ++i) {
+ int pos = (dims_[0] + 1)*j + i;
+ int new_pos = (new_dims_[0] + 1)*(j-jmin) + (i-imin);
+ // Copy all 6 coordinates for a pillar.
+ std::copy(COORD.begin() + 6*pos, COORD.begin() + 6*(pos + 1), new_COORD_.begin() + 6*new_pos);
+ if (resettoorigin) {
+ // Substract lowest x value from all X-coords, similarly for y, and truncate in z-direction
+ new_COORD_[6*new_pos] -= x_correction;
+ new_COORD_[6*new_pos + 1] -= y_correction;
+ new_COORD_[6*new_pos + 2] = 0;
+ new_COORD_[6*new_pos + 3] -= x_correction;
+ new_COORD_[6*new_pos + 4] -= y_correction;
+ new_COORD_[6*new_pos + 5] = zmax-zmin;
+ }
+ }
+ }
+
+ // Get the z limits, check if they must be changed to make a shoe-box.
+ // This means that zmin must be greater than or equal to the highest
+ // coordinate of the bottom surface, while zmax must be less than or
+ // equal to the lowest coordinate of the top surface.
+ int layersz = 8*dims_[0]*dims_[1];
+ const std::vector& ZCORN = parser_.getFloatingPointValue("ZCORN");
+ int num_zcorn = ZCORN.size();
+ if (num_zcorn != layersz*dims_[2]) {
+ std::cerr << "Error! ZCORN size (" << ZCORN.size() << ") not consistent with SPECGRID\n";
+ throw std::runtime_error("Inconsistent ZCORN and SPECGRID.");
+ }
+
+ zmin = std::max(zmin, botmax_);
+ zmax = std::min(zmax, topmin_);
+ if (zmin >= zmax) {
+ std::cerr << "Error: zmin >= zmax (zmin = " << zmin << ", zmax = " << zmax << ")\n";
+ throw std::runtime_error("zmin >= zmax");
+ }
+ std::cout << "Chopping subsample, i: (" << imin << "--" << imax << ") j: (" << jmin << "--" << jmax << ") z: (" << zmin << "--" << zmax << ")" << std::endl;
+
+ // We must find the maximum and minimum k value for the given z limits.
+ // First, find the first layer with a z-coordinate strictly above zmin.
+ int kmin = -1;
+ for (int k = 0; k < dims_[2]; ++k) {
+ double layer_max = *std::max_element(ZCORN.begin() + k*layersz, ZCORN.begin() + (k + 1)*layersz);
+ if (layer_max > zmin) {
+ kmin = k;
+ break;
+ }
+ }
+ // Then, find the last layer with a z-coordinate strictly below zmax.
+ int kmax = -1;
+ for (int k = dims_[2]; k > 0; --k) {
+ double layer_min = *std::min_element(ZCORN.begin() + (k - 1)*layersz, ZCORN.begin() + k*layersz);
+ if (layer_min < zmax) {
+ kmax = k;
+ break;
+ }
+ }
+ new_dims_[2] = kmax - kmin;
+
+ // Filter the ZCORN field, build mapping from new to old cells.
+ double z_origin_correction = 0.0;
+ if (resettoorigin) {
+ z_origin_correction = zmin;
+ }
+ new_ZCORN_.resize(8*new_dims_[0]*new_dims_[1]*new_dims_[2], 1e100);
+ new_to_old_cell_.resize(new_dims_[0]*new_dims_[1]*new_dims_[2], -1);
+ int cellcount = 0;
+ int delta[3] = { 1, 2*dims_[0], 4*dims_[0]*dims_[1] };
+ int new_delta[3] = { 1, 2*new_dims_[0], 4*new_dims_[0]*new_dims_[1] };
+ for (int k = kmin; k < kmax; ++k) {
+ for (int j = jmin; j < jmax; ++j) {
+ for (int i = imin; i < imax; ++i) {
+ new_to_old_cell_[cellcount++] = dims_[0]*dims_[1]*k + dims_[0]*j + i;
+ int old_ix = 2*(i*delta[0] + j*delta[1] + k*delta[2]);
+ int new_ix = 2*((i-imin)*new_delta[0] + (j-jmin)*new_delta[1] + (k-kmin)*new_delta[2]);
+ int old_indices[8] = { old_ix, old_ix + delta[0],
+ old_ix + delta[1], old_ix + delta[1] + delta[0],
+ old_ix + delta[2], old_ix + delta[2] + delta[0],
+ old_ix + delta[2] + delta[1], old_ix + delta[2] + delta[1] + delta[0] };
+ int new_indices[8] = { new_ix, new_ix + new_delta[0],
+ new_ix + new_delta[1], new_ix + new_delta[1] + new_delta[0],
+ new_ix + new_delta[2], new_ix + new_delta[2] + new_delta[0],
+ new_ix + new_delta[2] + new_delta[1], new_ix + new_delta[2] + new_delta[1] + new_delta[0] };
+ for (int cc = 0; cc < 8; ++cc) {
+ new_ZCORN_[new_indices[cc]] = std::min(zmax, std::max(zmin, ZCORN[old_indices[cc]])) - z_origin_correction;
+ }
+ }
+ }
+ }
+
+ filterIntegerField("ACTNUM", new_ACTNUM_);
+ filterDoubleField("PORO", new_PORO_);
+ filterDoubleField("PERMX", new_PERMX_);
+ filterDoubleField("PERMY", new_PERMY_);
+ filterDoubleField("PERMZ", new_PERMZ_);
+ filterIntegerField("SATNUM", new_SATNUM_);
+ }
+
+
+
+ /// Return a subparser with fields corresponding to the selected subset.
+ /// Note that the returned parser is NOT converted to SI, that must be done
+ /// by the user afterwards with the parser's convertToSI() method.
+ EclipseGridParser subparser()
+ {
+ if (parser_.hasField("FIELD") || parser_.hasField("LAB") || parser_.hasField("PVT-M")) {
+ THROW("CornerPointChopper::subparser() cannot handle any eclipse unit system other than METRIC.");
+ }
+
+ EclipseGridParser sp;
+ std::tr1::shared_ptr sg(new SPECGRID);
+ for (int dd = 0; dd < 3; ++dd) {
+ sg->dimensions[dd] = new_dims_[dd];
+ }
+ sp.setSpecialField("SPECGRID", sg);
+ sp.setFloatingPointField("COORD", new_COORD_);
+ sp.setFloatingPointField("ZCORN", new_ZCORN_);
+ if (!new_ACTNUM_.empty()) sp.setIntegerField("ACTNUM", new_ACTNUM_);
+ if (!new_PORO_.empty()) sp.setFloatingPointField("PORO", new_PORO_);
+ if (!new_PERMX_.empty()) sp.setFloatingPointField("PERMX", new_PERMX_);
+ if (!new_PERMY_.empty()) sp.setFloatingPointField("PERMY", new_PERMY_);
+ if (!new_PERMZ_.empty()) sp.setFloatingPointField("PERMZ", new_PERMZ_);
+ if (!new_SATNUM_.empty()) sp.setIntegerField("SATNUM", new_SATNUM_);
+ sp.computeUnits(); // Always METRIC, since that is default.
+ return sp;
+ }
+
+
+
+
+ void writeGrdecl(const std::string& filename)
+ {
+ // Output new versions of SPECGRID, COORD, ZCORN, ACTNUM, PERMX, PORO, SATNUM.
+ std::ofstream out(filename.c_str());
+ if (!out) {
+ std::cerr << "Could not open file " << filename << "\n";
+ throw std::runtime_error("Could not open output file.");
+ }
+ out << "SPECGRID\n" << new_dims_[0] << ' ' << new_dims_[1] << ' ' << new_dims_[2]
+ << " 1 F\n/\n\n";
+ out << "COORD\n";
+ int num_new_coord = new_COORD_.size();
+ for (int i = 0; i < num_new_coord/6; ++i) {
+ for (int j = 0; j < 6; ++j) {
+ out << " " << new_COORD_[6*i + j];
+ }
+ out << '\n';
+ }
+ out << "/\n\n";
+ out << "ZCORN\n";
+ int num_new_zcorn = new_ZCORN_.size();
+ assert(num_new_zcorn%8 == 0);
+ for (int i = 0; i < num_new_zcorn/8; ++i) {
+ for (int j = 0; j < 8; ++j) {
+ out << " " << new_ZCORN_[8*i + j];
+ }
+ out << '\n';
+ }
+ out << "/\n\n";
+
+ outputField(out, new_ACTNUM_, "ACTNUM");
+ outputField(out, new_PORO_, "PORO");
+ outputField(out, new_PERMX_, "PERMX");
+ outputField(out, new_PERMY_, "PERMY");
+ outputField(out, new_PERMZ_, "PERMZ");
+ outputField(out, new_SATNUM_, "SATNUM");
+ }
+
+ private:
+ EclipseGridParser parser_;
+ double botmax_;
+ double topmin_;
+ double abszmin_;
+ double abszmax_;
+ std::vector new_COORD_;
+ std::vector new_ZCORN_;
+ std::vector new_ACTNUM_;
+ std::vector new_PORO_;
+ std::vector new_PERMX_;
+ std::vector new_PERMY_;
+ std::vector new_PERMZ_;
+ std::vector new_SATNUM_;
+ int dims_[3];
+ int new_dims_[3];
+ std::vector new_to_old_cell_;
+
+
+ template
+ void outputField(std::ostream& os,
+ const std::vector& field,
+ const std::string& keyword)
+ {
+ if (field.empty()) return;
+
+ os << keyword << '\n';
+ int sz = field.size();
+ //int num_new_zcorn = new_ZCORN_.size();
+ //assert(sz%20 == 0);
+ const int nel_per_row = 20;
+ int num_full_rows=sz/nel_per_row;
+ int num_extra_entries=sz%nel_per_row;
+ for (int i = 0; i < num_full_rows; ++i) {
+ for (int j = 0; j < nel_per_row; ++j) {
+ os << " " << field[nel_per_row*i + j];
+ }
+ os << '\n';
+ }
+ for (int i = 0; i < num_extra_entries; ++i) {
+ os << " " << field[num_full_rows*nel_per_row+i];
+ }
+ os << "/\n\n";
+ }
+
+// os << keyword << '\n';
+// int sz = field.size();
+// T last = std::numeric_limits::max();
+// int repeats = 0;
+// for (int i = 0; i < sz; ++i) {
+// T val = field[i];
+// if (val == last) {
+// ++repeats;
+// } else {
+// if (repeats == 1) {
+// os << last << '\n';
+// } else if (repeats > 1) {
+// os << repeats << '*' << last << '\n';
+// }
+// last = val;
+// repeats = 1;
+// }
+// }
+// if (repeats == 1) {
+// os << last << '\n';
+// } else if (repeats > 1) {
+// os << repeats << '*' << last << '\n';
+// }
+// os << "/\n\n";
+// }
+
+
+ template
+ void filterField(const std::vector& field,
+ std::vector& output_field)
+ {
+ int sz = new_to_old_cell_.size();
+ output_field.resize(sz);
+ for (int i = 0; i < sz; ++i) {
+ output_field[i] = field[new_to_old_cell_[i]];
+ }
+ }
+
+ void filterDoubleField(const std::string& keyword, std::vector& output_field)
+ {
+ if (parser_.hasField(keyword)) {
+ const std::vector& field = parser_.getFloatingPointValue(keyword);
+ filterField(field, output_field);
+ }
+ }
+
+ void filterIntegerField(const std::string& keyword, std::vector& output_field)
+ {
+ if (parser_.hasField(keyword)) {
+ const std::vector& field = parser_.getIntegerValue(keyword);
+ filterField(field, output_field);
+ }
+ }
+
+ };
+
+}
+
+
+
+
+#endif // OPM_CORNERPOINTCHOPPER_HEADER_INCLUDED
diff --git a/opm/core/eclipse/EclipseGridInspector.cpp b/opm/core/eclipse/EclipseGridInspector.cpp
new file mode 100644
index 00000000..b62e4198
--- /dev/null
+++ b/opm/core/eclipse/EclipseGridInspector.cpp
@@ -0,0 +1,335 @@
+//===========================================================================
+//
+// File: EclipseGridInspector.C
+//
+// Created: Mon Jun 2 12:17:51 2008
+//
+// Author: Atgeirr F Rasmussen
+//
+// $Date$
+//
+// $Revision$
+//
+// Revision: $Id: EclipseGridInspector.C,v 1.2 2008/08/18 14:16:13 atgeirr Exp $
+//
+//===========================================================================
+
+/*
+ Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+ Copyright 2009, 2010 Statoil ASA.
+
+ This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+ OpenRS is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenRS is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenRS. If not, see .
+*/
+
+#if HAVE_CONFIG_H
+#include "config.h"
+#endif
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+namespace Dune
+{
+
+EclipseGridInspector::EclipseGridInspector(const EclipseGridParser& parser)
+ : parser_(parser)
+{
+ std::vector keywords;
+ keywords.push_back("COORD");
+ keywords.push_back("ZCORN");
+
+ if (!parser_.hasFields(keywords)) {
+ THROW("Needed field is missing in file");
+ }
+
+ if (parser_.hasField("SPECGRID")) {
+ const SPECGRID& sgr = parser.getSPECGRID();
+ logical_gridsize_[0] = sgr.dimensions[0];
+ logical_gridsize_[1] = sgr.dimensions[1];
+ logical_gridsize_[2] = sgr.dimensions[2];
+ } else if (parser_.hasField("DIMENS")) {
+ const std::vector& dim = parser.getIntegerValue("DIMENS");
+ logical_gridsize_[0] = dim[0];
+ logical_gridsize_[1] = dim[1];
+ logical_gridsize_[2] = dim[2];
+ } else {
+ THROW("Found neither SPECGRID nor DIMENS in file. At least one is needed.");
+ }
+
+}
+
+/**
+ Return the dip slopes for the cell relative to xy-plane in x- and y- direction.
+ Dip slope is average rise in positive x-direction over cell length in x-direction.
+ Similarly for y.
+
+ Current implementation is for vertical pillars, but is not difficult to fix.
+
+ @returns a std::pair with x-dip in first component and y-dip in second.
+*/
+std::pair EclipseGridInspector::cellDips(int i, int j, int k) const
+{
+ checkLogicalCoords(i, j, k);
+ const std::vector& pillc = parser_.getFloatingPointValue("COORD");
+ int num_pillars = (logical_gridsize_[0] + 1)*(logical_gridsize_[1] + 1);
+ if (6*num_pillars != int(pillc.size())) {
+ throw std::runtime_error("Wrong size of COORD field.");
+ }
+ const std::vector& z = parser_.getFloatingPointValue("ZCORN");
+ int num_cells = logical_gridsize_[0]*logical_gridsize_[1]*logical_gridsize_[2];
+ if (8*num_cells != int(z.size())) {
+ throw std::runtime_error("Wrong size of ZCORN field");
+ }
+
+ // Pick ZCORN-value for all 8 corners of the given cell
+ boost::array cellz = cellZvals(i, j, k);
+
+ // Compute rise in positive x-direction for all four edges (and then find mean)
+ // Current implementation is for regularly placed and vertical pillars!
+ int numxpill = logical_gridsize_[0] + 1;
+ int pix = i + j*numxpill;
+ double cell_xlength = pillc[6*(pix + 1)] - pillc[6*pix];
+ flush(std::cout);
+ double xrise[4] = { (cellz[1] - cellz[0])/cell_xlength, // LLL -> HLL
+ (cellz[3] - cellz[2])/cell_xlength, // LHL -> HHL
+ (cellz[5] - cellz[4])/cell_xlength, // LLH -> HLH
+ (cellz[7] - cellz[6])/cell_xlength}; // LHH -> HHH
+
+ double cell_ylength = pillc[6*(pix + numxpill) + 1] - pillc[6*pix + 1];
+ double yrise[4] = { (cellz[2] - cellz[0])/cell_ylength, // LLL -> LHL
+ (cellz[3] - cellz[1])/cell_ylength, // HLL -> HHL
+ (cellz[6] - cellz[4])/cell_ylength, // LLH -> LHH
+ (cellz[7] - cellz[5])/cell_ylength}; // HLH -> HHH
+
+
+ // Now ignore those edges that touch the global top or bottom surface
+ // of the entire grdecl model. This is to avoid bias, as these edges probably
+ // don't follow an overall dip for the model if it exists.
+ int x_edges = 4;
+ int y_edges = 4;
+ boost::array gridlimits = getGridLimits();
+ double zmin = gridlimits[4];
+ double zmax = gridlimits[5];
+ // LLL -> HLL
+ if ((cellz[1] == zmin) || (cellz[0] == zmin)) {
+ xrise[0] = 0; x_edges--;
+ }
+ // LHL -> HHL
+ if ((cellz[2] == zmin) || (cellz[3] == zmin)) {
+ xrise[1] = 0; x_edges--;
+ }
+ // LLH -> HLH
+ if ((cellz[4] == zmax) || (cellz[5] == zmax)) {
+ xrise[2] = 0; x_edges--;
+ }
+ // LHH -> HHH
+ if ((cellz[6] == zmax) || (cellz[7] == zmax)) {
+ xrise[3] = 0; x_edges--;
+ }
+ // LLL -> LHL
+ if ((cellz[0] == zmin) || (cellz[2] == zmin)) {
+ yrise[0] = 0; y_edges--;
+ }
+ // HLL -> HHL
+ if ((cellz[1] == zmin) || (cellz[3] == zmin)) {
+ yrise[1] = 0; y_edges--;
+ }
+ // LLH -> LHH
+ if ((cellz[6] == zmax) || (cellz[4] == zmax)) {
+ yrise[2] = 0; y_edges--;
+ }
+ // HLH -> HHH
+ if ((cellz[7] == zmax) || (cellz[5] == zmax)) {
+ yrise[3] = 0; y_edges--;
+ }
+
+ return std::make_pair( (xrise[0] + xrise[1] + xrise[2] + xrise[3])/x_edges,
+ (yrise[0] + yrise[1] + yrise[2] + yrise[3])/y_edges);
+}
+/**
+ Wrapper for cellDips(i, j, k).
+*/
+std::pair EclipseGridInspector::cellDips(int cell_idx) const
+{
+ boost::array idxs = cellIdxToLogicalCoords(cell_idx);
+ return cellDips(idxs[0], idxs[1], idxs[2]);
+}
+
+boost::array EclipseGridInspector::cellIdxToLogicalCoords(int cell_idx) const
+{
+
+ int i,j,k; // Position of cell in cell hierarchy
+ int horIdx = (cell_idx+1) - int(std::floor(((double)(cell_idx+1))/((double)(logical_gridsize_[0]*logical_gridsize_[1]))))*logical_gridsize_[0]*logical_gridsize_[1]; // index in the corresponding horizon
+ if (horIdx == 0) {
+ horIdx = logical_gridsize_[0]*logical_gridsize_[1];
+ }
+ i = horIdx - int(std::floor(((double)horIdx)/((double)logical_gridsize_[0])))*logical_gridsize_[0];
+ if (i == 0) {
+ i = logical_gridsize_[0];
+ }
+ j = (horIdx-i)/logical_gridsize_[0]+1;
+ k = ((cell_idx+1)-logical_gridsize_[0]*(j-1)-1)/(logical_gridsize_[0]*logical_gridsize_[1])+1;
+
+ boost::array a = {{i-1, j-1, k-1}};
+ return a; //boost::array {{i-1, j-1, k-1}};
+}
+
+double EclipseGridInspector::cellVolumeVerticalPillars(int i, int j, int k) const
+{
+ // Checking parameters and obtaining values from parser.
+ checkLogicalCoords(i, j, k);
+ const std::vector& pillc = parser_.getFloatingPointValue("COORD");
+ int num_pillars = (logical_gridsize_[0] + 1)*(logical_gridsize_[1] + 1);
+ if (6*num_pillars != int(pillc.size())) {
+ throw std::runtime_error("Wrong size of COORD field.");
+ }
+ const std::vector& z = parser_.getFloatingPointValue("ZCORN");
+ int num_cells = logical_gridsize_[0]*logical_gridsize_[1]*logical_gridsize_[2];
+ if (8*num_cells != int(z.size())) {
+ throw std::runtime_error("Wrong size of ZCORN field");
+ }
+
+ // Computing the base area as half the 2d cross product of the diagonals.
+ int numxpill = logical_gridsize_[0] + 1;
+ int pix = i + j*numxpill;
+ double px[4] = { pillc[6*pix],
+ pillc[6*(pix + 1)],
+ pillc[6*(pix + numxpill)],
+ pillc[6*(pix + numxpill + 1)] };
+ double py[4] = { pillc[6*pix + 1],
+ pillc[6*(pix + 1) + 1],
+ pillc[6*(pix + numxpill) + 1],
+ pillc[6*(pix + numxpill + 1) + 1] };
+ double diag1[2] = { px[3] - px[0], py[3] - py[0] };
+ double diag2[2] = { px[2] - px[1], py[2] - py[1] };
+ double area = 0.5*(diag1[0]*diag2[1] - diag1[1]*diag2[0]);
+
+ // Computing the average of the z-differences along each pillar.
+ int delta[3] = { 1,
+ 2*logical_gridsize_[0],
+ 4*logical_gridsize_[0]*logical_gridsize_[1] };
+ int ix = 2*(i*delta[0] + j*delta[1] + k*delta[2]);
+ double cellz[8] = { z[ix], z[ix + delta[0]],
+ z[ix + delta[1]], z[ix + delta[1] + delta[0]],
+ z[ix + delta[2]], z[ix + delta[2] + delta[0]],
+ z[ix + delta[2] + delta[1]], z[ix + delta[2] + delta[1] + delta[0]] };
+ double diffz[4] = { cellz[4] - cellz[0],
+ cellz[5] - cellz[1],
+ cellz[6] - cellz[2],
+ cellz[7] - cellz[3] };
+ double averzdiff = 0.25*std::accumulate(diffz, diffz + 4, 0.0);
+ return averzdiff*area;
+}
+
+
+double EclipseGridInspector::cellVolumeVerticalPillars(int cell_idx) const
+{
+ boost::array idxs = cellIdxToLogicalCoords(cell_idx);
+ return cellVolumeVerticalPillars(idxs[0], idxs[1], idxs[2]);
+}
+
+void EclipseGridInspector::checkLogicalCoords(int i, int j, int k) const
+{
+ if (i < 0 || i >= logical_gridsize_[0])
+ throw std::runtime_error("First coordinate out of bounds");
+ if (j < 0 || j >= logical_gridsize_[1])
+ throw std::runtime_error("Second coordinate out of bounds");
+ if (k < 0 || k >= logical_gridsize_[2])
+ throw std::runtime_error("Third coordinate out of bounds");
+}
+
+
+boost::array EclipseGridInspector::getGridLimits() const
+{
+ if (! (parser_.hasField("COORD") && parser_.hasField("ZCORN") && parser_.hasField("SPECGRID")) ) {
+ throw std::runtime_error("EclipseGridInspector: Grid does not have SPECGRID, COORD, and ZCORN, can't find dimensions.");
+ }
+
+ std::vector coord = parser_.getFloatingPointValue("COORD");
+ std::vector zcorn = parser_.getFloatingPointValue("ZCORN");
+
+ double xmin = +DBL_MAX;
+ double xmax = -DBL_MAX;
+ double ymin = +DBL_MAX;
+ double ymax = -DBL_MAX;
+
+
+ int pillars = (logical_gridsize_[0]+1) * (logical_gridsize_[1]+1);
+
+ for (int pillarindex = 0; pillarindex < pillars; ++pillarindex) {
+ if (coord[pillarindex * 6 + 0] > xmax)
+ xmax = coord[pillarindex * 6 + 0];
+ if (coord[pillarindex * 6 + 0] < xmin)
+ xmin = coord[pillarindex * 6 + 0];
+ if (coord[pillarindex * 6 + 1] > ymax)
+ ymax = coord[pillarindex * 6 + 1];
+ if (coord[pillarindex * 6 + 1] < ymin)
+ ymin = coord[pillarindex * 6 + 1];
+ if (coord[pillarindex * 6 + 3] > xmax)
+ xmax = coord[pillarindex * 6 + 3];
+ if (coord[pillarindex * 6 + 3] < xmin)
+ xmin = coord[pillarindex * 6 + 3];
+ if (coord[pillarindex * 6 + 4] > ymax)
+ ymax = coord[pillarindex * 6 + 4];
+ if (coord[pillarindex * 6 + 4] < ymin)
+ ymin = coord[pillarindex * 6 + 4];
+ }
+
+ boost::array gridlimits = {{ xmin, xmax, ymin, ymax,
+ *min_element(zcorn.begin(), zcorn.end()),
+ *max_element(zcorn.begin(), zcorn.end()) }};
+ return gridlimits;
+}
+
+
+
+boost::array EclipseGridInspector::gridSize() const
+{
+ boost::array retval = {{ logical_gridsize_[0],
+ logical_gridsize_[1],
+ logical_gridsize_[2] }};
+ return retval;
+}
+
+
+boost::array EclipseGridInspector::cellZvals(int i, int j, int k) const
+{
+ // Get the zcorn field.
+ const std::vector& z = parser_.getFloatingPointValue("ZCORN");
+ int num_cells = logical_gridsize_[0]*logical_gridsize_[1]*logical_gridsize_[2];
+ if (8*num_cells != int(z.size())) {
+ throw std::runtime_error("Wrong size of ZCORN field");
+ }
+
+ // Make the coordinate array.
+ int delta[3] = { 1,
+ 2*logical_gridsize_[0],
+ 4*logical_gridsize_[0]*logical_gridsize_[1] };
+ int ix = 2*(i*delta[0] + j*delta[1] + k*delta[2]);
+ boost::array cellz = {{ z[ix], z[ix + delta[0]],
+ z[ix + delta[1]], z[ix + delta[1] + delta[0]],
+ z[ix + delta[2]], z[ix + delta[2] + delta[0]],
+ z[ix + delta[2] + delta[1]], z[ix + delta[2] + delta[1] + delta[0]] }};
+ return cellz;
+}
+
+
+} // namespace Dune
diff --git a/opm/core/eclipse/EclipseGridInspector.hpp b/opm/core/eclipse/EclipseGridInspector.hpp
new file mode 100644
index 00000000..c314d28e
--- /dev/null
+++ b/opm/core/eclipse/EclipseGridInspector.hpp
@@ -0,0 +1,105 @@
+//===========================================================================
+//
+// File: EclipseGridInspector.h
+//
+// Created: Mon Jun 2 09:46:08 2008
+//
+// Author: Atgeirr F Rasmussen
+//
+// $Date$
+//
+// Revision: $Id: EclipseGridInspector.h,v 1.2 2008/08/18 14:16:12 atgeirr Exp $
+//
+//===========================================================================
+
+/*
+Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+Copyright 2009, 2010 Statoil ASA.
+
+This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+OpenRS is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+OpenRS is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with OpenRS. If not, see .
+*/
+
+#ifndef SINTEF_ECLIPSEGRIDINSPECTOR_HEADER
+#define SINTEF_ECLIPSEGRIDINSPECTOR_HEADER
+
+#include
+#include
+
+
+namespace Dune
+{
+
+/**
+ @brief A class for inspecting the contents of an eclipse file.
+
+ Given an EclipseGridParser that has successfully read en Eclipse
+ .grdecl-file, this class may be used to answer certain queries about
+ its contents.
+
+ @author Atgeirr F. Rasmussen
+ @date 2008/06/02 09:46:08
+*/
+
+class EclipseGridParser;
+
+class EclipseGridInspector
+{
+public:
+ /// Constructor taking a parser as argument.
+ /// The parser must already have read an Eclipse file.
+ EclipseGridInspector(const EclipseGridParser& parser);
+
+ /// Assuming that the pillars are vertical, compute the
+ /// volume of the cell given by logical coordinates (i, j, k).
+ double cellVolumeVerticalPillars(int i, int j, int k) const;
+
+ /// Assuming that the pillars are vertical, compute the
+ /// volume of the cell given by the cell index
+ double cellVolumeVerticalPillars(int cell_idx) const;
+
+ /// Compute the average dip in x- and y-direction of the
+ /// cell tops and bottoms relative to the xy-plane
+ std::pair cellDips(int i, int j, int k) const;
+ std::pair cellDips(int cell_idx) const;
+
+ // Convert global cell index to logical ijk-coordinates
+ boost::array cellIdxToLogicalCoords(int cell_idx) const;
+
+ /// Returns a vector with the outer limits of grid (in the grid's unit).
+ /// The vector contains [xmin, xmax, ymin, ymax, zmin, zmax], as
+ /// read from COORDS and ZCORN
+ boost::array getGridLimits() const;
+
+ /// Returns the extent of the logical cartesian grid
+ /// as number of cells in the (i, j, k) directions.
+ boost::array gridSize() const;
+
+ /// Returns the eight z-values associated with a given cell.
+ /// The ordering is such that i runs fastest. That is, with
+ /// L = low and H = high:
+ /// {LLL, HLL, LHL, HHL, LLH, HLH, LHH, HHH }.
+ boost::array cellZvals(int i, int j, int k) const;
+
+private:
+ const EclipseGridParser& parser_;
+ int logical_gridsize_[3];
+ void checkLogicalCoords(int i, int j, int k) const;
+};
+
+} // namespace Dune
+
+#endif // SINTEF_ECLIPSEGRIDINSPECTOR_HEADER
+
diff --git a/opm/core/eclipse/EclipseGridParser.cpp b/opm/core/eclipse/EclipseGridParser.cpp
new file mode 100644
index 00000000..199c140c
--- /dev/null
+++ b/opm/core/eclipse/EclipseGridParser.cpp
@@ -0,0 +1,571 @@
+//===========================================================================
+//
+// File: EclipseGridParser.C
+//
+// Created: Thu Dec 6 08:46:05 2007
+//
+// Author: Atgeirr F Rasmussen
+//
+// $Date$
+//
+// $Revision$
+//
+// Revision: $Id: EclipseGridParser.C,v 1.4 2008/08/18 14:16:14 atgeirr Exp $
+//
+//===========================================================================
+
+/*
+ Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+ Copyright 2009, 2010 Statoil ASA.
+
+ This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+ OpenRS is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenRS is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenRS. If not, see .
+*/
+#if HAVE_CONFIG_H
+#include "config.h"
+#endif
+#include
+#include
+#include
+#include
+#include
+#include
+#include "EclipseGridParser.hpp"
+#include "EclipseGridParserHelpers.hpp"
+#include "SpecialEclipseFields.hpp"
+#include
+#include
+#include
+
+using namespace std;
+
+//#define VERBOSE
+
+namespace Dune
+{
+
+// ---------- List of supported keywords ----------
+
+namespace EclipseKeywords
+{
+ string integer_fields[] =
+ { string("ACTNUM"),
+ string("SATNUM"),
+ string("EQLNUM"),
+ string("REGNUM"),
+ string("ROCKTYPE"),
+ string("DIMENS"),
+ string("REGDIMS"),
+ string("WELLDIMS"),
+ string("TABDIMS"),
+ string("FIPNUM"),
+ string("GRIDFILE")
+ };
+ const int num_integer_fields = sizeof(integer_fields) / sizeof(integer_fields[0]);
+
+ string floating_fields[] =
+ { string("COORD"), string("ZCORN"), string("PERMX"),
+ string("PERMY"), string("PERMZ"), string("PERMXX"),
+ string("PERMYY"), string("PERMZZ"), string("PERMXY"),
+ string("PERMYZ"), string("PERMZX"), string("PORO"),
+ string("BULKMOD"), string("YOUNGMOD"), string("LAMEMOD"),
+ string("SHEARMOD"), string("POISSONMOD"), string("PWAVEMOD"),
+ string("MULTPV"), string("PRESSURE"), string("SGAS"),
+ string("SWAT")
+ };
+ const int num_floating_fields = sizeof(floating_fields) / sizeof(floating_fields[0]);
+
+ string special_fields[] =
+ { string("SPECGRID"), string("FAULTS"), string("MULTFLT"),
+ string("TITLE"), string("START"), string("DATES"),
+ string("DENSITY"), string("PVDG"), string("PVDO"),
+ string("PVTG"), string("PVTO"), string("PVTW"),
+ string("SGOF"), string("SWOF"), string("ROCK"),
+ string("ROCKTAB"), string("WELSPECS"), string("COMPDAT"),
+ string("WCONINJE"), string("WCONPROD"), string("WELTARG"),
+ string("EQUIL"), string("PVCDO"),
+ // The following fields only have a dummy implementation
+ // that allows us to ignore them.
+ "SWFN",
+ "SOF2",
+ "TUNING"
+ };
+
+ const int num_special_fields = sizeof(special_fields) / sizeof(special_fields[0]);
+
+ string ignore_with_data[] =
+ { string("MAPUNITS"), string("MAPAXES"), string("GRIDUNIT"),
+ string("NTG"), string("REGDIMS"), string("WELLDIMS"),
+ string("NSTACK"), string("SATNUM"),
+ string("RPTRST"), string("ROIP"), string("RWIP"),
+ string("RWSAT"), string("RPR"), string("WBHP"),
+ string("WOIR"), string("TSTEP"), string("BOX"),
+ string("COORDSYS"), string("PBVD")
+ };
+ const int num_ignore_with_data = sizeof(ignore_with_data) / sizeof(ignore_with_data[0]);
+
+ string ignore_no_data[] =
+ { string("RUNSPEC"), string("WATER"), string("OIL"),
+ string("METRIC"), string("FMTIN"), string("FMTOUT"),
+ string("GRID"), string("INIT"), string("NOECHO"),
+ string("ECHO"), string("EDIT"), string("PROPS"),
+ string("REGIONS"), string("SOLUTION"), string("SUMMARY"),
+ string("FPR"), string("FOIP"), string("FWIP"),
+ string("RUNSUM"), string("EXCEL"), string("SCHEDULE"),
+ string("END"), string("ENDBOX"), string("CONTINUE"),
+ string("NONNC"), string("GAS"), string("DISGAS"),
+ string("FIELD")
+ };
+ const int num_ignore_no_data = sizeof(ignore_no_data) / sizeof(ignore_no_data[0]);
+
+ string include_keywords[] = { string("INCLUDE") };
+ const int num_include_keywords = sizeof(include_keywords) / sizeof(include_keywords[0]);
+
+
+} // namespace EclipseKeywords
+
+namespace {
+
+ enum FieldType {
+ Integer,
+ FloatingPoint,
+ SpecialField,
+ IgnoreWithData,
+ IgnoreNoData,
+ Include,
+ Unknown
+ };
+
+ inline FieldType classifyKeyword(const string& keyword)
+ {
+ using namespace EclipseKeywords;
+ if (count(integer_fields, integer_fields + num_integer_fields, keyword)) {
+ return Integer;
+ } else if (count(floating_fields, floating_fields + num_floating_fields, keyword)) {
+ return FloatingPoint;
+ } else if (count(special_fields, special_fields + num_special_fields, keyword)) {
+ return SpecialField;
+ } else if (count(ignore_with_data, ignore_with_data + num_ignore_with_data, keyword)) {
+ return IgnoreWithData;
+ } else if (count(ignore_no_data, ignore_no_data + num_ignore_no_data, keyword)) {
+ return IgnoreNoData;
+ } else if (count(include_keywords, include_keywords + num_include_keywords, keyword)) {
+ return Include;
+ } else {
+ return Unknown;
+ }
+ }
+
+} // anon namespace
+
+
+
+// ---------- Member functions ----------
+
+/// Default constructor.
+//---------------------------------------------------------------------------
+EclipseGridParser::EclipseGridParser()
+//---------------------------------------------------------------------------
+{
+}
+
+
+/// Constructor taking an eclipse filename.
+//---------------------------------------------------------------------------
+EclipseGridParser::EclipseGridParser(const string& filename, bool convert_to_SI)
+//---------------------------------------------------------------------------
+{
+ // Store directory of filename
+ boost::filesystem::path p(filename);
+ directory_ = p.parent_path().string();
+ ifstream is(filename.c_str());
+ if (!is) {
+ cerr << "Unable to open file " << filename << endl;
+ throw exception();
+ }
+ read(is, convert_to_SI);
+}
+
+
+/// Read the given stream, overwriting any previous data.
+//---------------------------------------------------------------------------
+void EclipseGridParser::read(istream& is, bool convert_to_SI)
+//---------------------------------------------------------------------------
+{
+ integer_field_map_.clear();
+ floating_field_map_.clear();
+ special_field_map_.clear();
+
+ readImpl(is);
+ computeUnits();
+ if (convert_to_SI) {
+ convertToSI();
+ }
+}
+
+//---------------------------------------------------------------------------
+void EclipseGridParser::readImpl(istream& is)
+//---------------------------------------------------------------------------
+{
+ if (!is) {
+ cerr << "Could not read given input stream." << endl;
+ throw exception();
+ }
+
+ // Make temporary maps that will at the end be swapped with the
+ // member maps
+ // NOTE: Above is no longer true, for easier implementation of
+ // the INCLUDE keyword. We lose the strong exception guarantee,
+ // though (of course retaining the basic guarantee).
+ map >& intmap = integer_field_map_;
+ map >& floatmap = floating_field_map_;
+ map >& specialmap = special_field_map_;
+
+ // Actually read the data
+ is >> ignoreWhitespace;
+ while (!is.eof()) {
+ string keyword = readKeyword(is);
+#ifdef VERBOSE
+ cout << "Keyword found: " << keyword << endl;
+#endif
+ FieldType type = classifyKeyword(keyword);
+ switch (type) {
+ case Integer:
+ readVectorData(is, intmap[keyword]);
+ break;
+ case FloatingPoint:
+ readVectorData(is, floatmap[keyword]);
+ break;
+ case SpecialField: {
+ std::tr1::shared_ptr sb_ptr = createSpecialField(is, keyword);
+ if (sb_ptr) {
+ specialmap[keyword] = sb_ptr;
+ } else {
+ THROW("Could not create field " << keyword);
+ }
+ break;
+ }
+ case IgnoreWithData: {
+ ignored_fields_.insert(keyword);
+ is >> ignoreSlashLine;
+#ifdef VERBOSE
+ cout << "(ignored)" << endl;
+#endif
+ break;
+ }
+ case IgnoreNoData: {
+ ignored_fields_.insert(keyword);
+ is >> ignoreLine;
+#ifdef VERBOSE
+ cout << "(ignored)" << endl;
+#endif
+ break;
+ }
+ case Include: {
+ string include_filename = readString(is);
+ if (!directory_.empty()) {
+ include_filename = directory_ + '/' + include_filename;
+ }
+ ifstream include_is(include_filename.c_str());
+ if (!include_is) {
+ THROW("Unable to open INCLUDEd file " << include_filename);
+ }
+ readImpl(include_is);
+ is >> ignoreSlashLine;
+ break;
+ }
+ case Unknown:
+ default:
+ cerr << "Keyword " << keyword << " not recognized." << endl;
+ throw exception();
+ }
+ is >> ignoreWhitespace;
+ }
+
+#define VERBOSE_LIST_FIELDS 0
+#if VERBOSE_LIST_FIELDS
+ std::cout << "\nInteger fields:\n";
+ for (std::map >::iterator
+ i = intmap.begin(); i != intmap.end(); ++i)
+ std::cout << '\t' << i->first << '\n';
+
+ std::cout << "\nFloat fields:\n";
+ for (std::map >::iterator
+ i = floatmap.begin(); i != floatmap.end(); ++i)
+ std::cout << '\t' << i->first << '\n';
+
+ std::cout << "\nSpecial fields:\n";
+ for (std::map >::iterator
+ i = specialmap.begin(); i != specialmap.end(); ++i)
+ std::cout << '\t' << i->first << '\n';
+#endif
+}
+
+
+
+//---------------------------------------------------------------------------
+void EclipseGridParser::convertToSI()
+//---------------------------------------------------------------------------
+{
+ // Convert all special fields.
+ typedef std::map >::iterator SpecialIt;
+ for (SpecialIt i = special_field_map_.begin(); i != special_field_map_.end(); ++i) {
+ i->second->convertToSI(units_);
+ }
+
+ // Convert all floating point fields.
+ typedef std::map >::iterator FloatIt;
+ for (FloatIt i = floating_field_map_.begin(); i != floating_field_map_.end(); ++i) {
+ const std::string& key = i->first;
+ std::vector& field = i->second;
+ // Find the right unit.
+ double unit = 1e100;
+ if (key == "COORD" || key == "ZCORN") {
+ unit = units_.length;
+ } else if (key == "PERMX" || key == "PERMY" || key == "PERMZ" ||
+ key == "PERMXX" || key == "PERMYY" || key == "PERMZZ" ||
+ key == "PERMXY" || key == "PERMYZ" || key == "PERMZX") {
+ unit = units_.permeability;
+ } else if (key == "PORO" || key == "BULKMOD" || key == "YOUNGMOD" ||
+ key == "LAMEMOD" || key == "SHEARMOD" || key == "POISSONMOD" ||
+ key == "PWAVEMOD" || key == "MULTPV" || key == "PWAVEMOD" ||
+ key == "SGAS" || key == "SWAT") {
+ unit = 1.0;
+ } else if (key == "PRESSURE") {
+ unit = units_.pressure;
+ } else {
+ THROW("Units for field " << key << " not specified. Cannon convert to SI.");
+ }
+ // Convert.
+ for (int i = 0; i < int(field.size()); ++i) {
+ field[i] = Dune::unit::convert::from(field[i], unit);
+ }
+ }
+
+ // Set all units to one.
+ units_.setToOne();
+}
+
+
+/// Returns true if the given keyword corresponds to a field that
+/// was found in the file.
+//---------------------------------------------------------------------------
+bool EclipseGridParser::hasField(const string& keyword) const
+//---------------------------------------------------------------------------
+{
+ string ukey = upcase(keyword);
+ return integer_field_map_.count(ukey) || floating_field_map_.count(ukey) ||
+ special_field_map_.count(ukey) || ignored_fields_.count(ukey);
+}
+
+
+/// Returns true if all the given keywords correspond to fields
+/// that were found in the file.
+//---------------------------------------------------------------------------
+bool EclipseGridParser::hasFields(const vector& keywords) const
+//---------------------------------------------------------------------------
+{
+ int num_keywords = keywords.size();
+ for (int i = 0; i < num_keywords; ++i) {
+ if (!hasField(keywords[i])) {
+ return false;
+ }
+ }
+ return true;
+}
+
+//---------------------------------------------------------------------------
+vector EclipseGridParser::fieldNames() const
+//---------------------------------------------------------------------------
+{
+ vector names;
+ names.reserve(integer_field_map_.size() +
+ floating_field_map_.size() +
+ special_field_map_.size() +
+ ignored_fields_.size());
+ {
+ map >::const_iterator it = integer_field_map_.begin();
+ for (; it != integer_field_map_.end(); ++it) {
+ names.push_back(it->first);
+ }
+ }
+ {
+ map >::const_iterator it = floating_field_map_.begin();
+ for (; it != floating_field_map_.end(); ++it) {
+ names.push_back(it->first);
+ }
+ }
+ {
+ map >::const_iterator it = special_field_map_.begin();
+ for (; it != special_field_map_.end(); ++it) {
+ names.push_back(it->first);
+ }
+ }
+ {
+ set::const_iterator it = ignored_fields_.begin();
+ for (; it != ignored_fields_.end(); ++it) {
+ names.push_back(*it);
+ }
+ }
+ return names;
+}
+
+//---------------------------------------------------------------------------
+const std::vector& EclipseGridParser::getIntegerValue(const std::string& keyword) const
+//---------------------------------------------------------------------------
+{
+ if (keyword == "SPECGRID") {
+ cerr << "\nERROR. Interface has changed!\n"
+ << "const vector& dim = parser.getIntegerValue(""SPECGRID"") is deprecated.\n"
+ << "Use:\n"
+ << "const SPECGRID& specgrid = parser.getSPECGRID();\n"
+ << "const vector& dim = specgrid.dimensions;\n\n";
+ throw exception();
+ }
+
+ map >::const_iterator it
+ = integer_field_map_.find(keyword);
+ if (it == integer_field_map_.end()) {
+ return empty_integer_field_;
+ } else {
+ return it->second;
+ }
+}
+
+//---------------------------------------------------------------------------
+const std::vector& EclipseGridParser::getFloatingPointValue(const std::string& keyword) const
+//---------------------------------------------------------------------------
+{
+ map >::const_iterator it
+ = floating_field_map_.find(keyword);
+ if (it == floating_field_map_.end()) {
+ return empty_floating_field_;
+ } else {
+ return it->second;
+ }
+}
+
+
+//---------------------------------------------------------------------------
+const std::tr1::shared_ptr EclipseGridParser::getSpecialValue(const std::string& keyword) const
+//---------------------------------------------------------------------------
+{
+ map >::const_iterator it = special_field_map_.find(keyword);
+ if (it == special_field_map_.end()) {
+ THROW("No such field: " << keyword);
+ } else {
+ return it->second;
+ }
+}
+
+//---------------------------------------------------------------------------
+std::tr1::shared_ptr
+EclipseGridParser::createSpecialField(std::istream& is,
+ const std::string& fieldname)
+//---------------------------------------------------------------------------
+{
+ string ukey = upcase(fieldname);
+ std::tr1::shared_ptr spec_ptr
+ = Factory::createObject(fieldname);
+ spec_ptr->read(is);
+ return spec_ptr;
+}
+
+//---------------------------------------------------------------------------
+void EclipseGridParser::setIntegerField(const std::string& keyword,
+ const std::vector& field)
+//---------------------------------------------------------------------------
+{
+ integer_field_map_[keyword] = field;
+}
+
+//---------------------------------------------------------------------------
+void EclipseGridParser::setFloatingPointField(const std::string& keyword,
+ const std::vector& field)
+//---------------------------------------------------------------------------
+{
+ floating_field_map_[keyword] = field;
+}
+
+//---------------------------------------------------------------------------
+void EclipseGridParser::setSpecialField(const std::string& keyword,
+ std::tr1::shared_ptr field)
+//---------------------------------------------------------------------------
+{
+ special_field_map_[keyword] = field;
+}
+
+//---------------------------------------------------------------------------
+const EclipseUnits& EclipseGridParser::units() const
+//---------------------------------------------------------------------------
+{
+ return units_;
+}
+
+//---------------------------------------------------------------------------
+void EclipseGridParser::computeUnits()
+//---------------------------------------------------------------------------
+{
+ // Decide unit family.
+ enum EclipseUnitFamily { Metric = 0, Field = 1, Lab = 2, Pvtm = 3 };
+ EclipseUnitFamily unit_family = Metric; // The default.
+ if (hasField("FIELD")) unit_family = Field;
+ if (hasField("LAB")) unit_family = Lab;
+ if (hasField("PVT-M")) unit_family = Pvtm;
+
+ // Set units.
+ using namespace prefix;
+ using namespace unit;
+ switch (unit_family) {
+ case Metric:
+ units_.length = meter;
+ units_.time = day;
+ units_.density = kilogram/cubic(meter);
+ units_.pressure = barsa;
+ units_.compressibility = 1.0/barsa;
+ units_.viscosity = centi*Poise;
+ units_.permeability = milli*darcy;
+ units_.liqvol_s = cubic(meter);
+ units_.liqvol_r = cubic(meter);
+ units_.gasvol_s = cubic(meter);
+ units_.gasvol_r = cubic(meter);
+ units_.transmissibility = centi*Poise * cubic(meter) / (day * barsa);
+ break;
+ case Field:
+ units_.length = feet;
+ units_.time = day;
+ units_.density = pound/cubic(feet);
+ units_.pressure = psia;
+ units_.compressibility = 1.0/psia;
+ units_.viscosity = centi*Poise;
+ units_.permeability = milli*darcy;
+ units_.liqvol_s = stb;
+ units_.liqvol_r = stb;
+ units_.gasvol_s = 1000*cubic(feet); // Prefix 'M' is 1000
+ units_.gasvol_r = stb;
+ units_.transmissibility = centi*Poise * stb / (day * psia);
+ break;
+ case Lab:
+ THROW("Unhandled unit family " << unit_family);
+ break;
+ case Pvtm:
+ THROW("Unhandled unit family " << unit_family);
+ break;
+ default:
+ THROW("Unknown unit family " << unit_family);
+ }
+}
+
+} // namespace Dune
diff --git a/opm/core/eclipse/EclipseGridParser.hpp b/opm/core/eclipse/EclipseGridParser.hpp
new file mode 100644
index 00000000..4a922864
--- /dev/null
+++ b/opm/core/eclipse/EclipseGridParser.hpp
@@ -0,0 +1,191 @@
+//===========================================================================
+//
+// File: EclipseGridParser.h
+//
+// Created: Wed Dec 5 17:05:13 2007
+//
+// Author: Atgeirr F Rasmussen
+//
+// $Date$
+//
+// Revision: $Id: EclipseGridParser.h,v 1.3 2008/08/18 14:16:13 atgeirr Exp $
+//
+//===========================================================================
+
+/*
+Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+Copyright 2009, 2010 Statoil ASA.
+
+This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+OpenRS is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+OpenRS is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with OpenRS. If not, see .
+*/
+
+#ifndef SINTEF_ECLIPSEGRIDPARSER_HEADER
+#define SINTEF_ECLIPSEGRIDPARSER_HEADER
+
+#include
+#include
+#include
+#include