Merge pull request #518 from blattms/master-refactor-for-cpgrid-support

refactor for cpgrid support to apply to master

CMakeLists_files.cmake

@@ -29,6 +29,7 @@
 # originally generated with the command:
 # find opm -name '*.c*' -printf '\t%p\n' | sort
 list (APPEND MAIN_SOURCE_FILES
+  opm/core/grid/GridHelpers.cpp
 	opm/core/grid/GridManager.cpp
 	opm/core/grid/grid.c
 	opm/core/grid/cart_grid.c
@@ -74,6 +75,7 @@ list (APPEND MAIN_SOURCE_FILES
 	opm/core/pressure/tpfa/compr_quant_general.c
 	opm/core/pressure/tpfa/compr_source.c
 	opm/core/pressure/tpfa/ifs_tpfa.c
+	opm/core/pressure/tpfa/TransTpfa.cpp
 	opm/core/pressure/tpfa/trans_tpfa.c
 	opm/core/pressure/legacy_well.c
 	opm/core/props/BlackoilPropertiesBasic.cpp
@@ -248,6 +250,7 @@ list (APPEND PUBLIC_HEADER_FILES
 	opm/core/grid/CellQuadrature.hpp
 	opm/core/grid/ColumnExtract.hpp
 	opm/core/grid/FaceQuadrature.hpp
+	opm/core/grid/GridHelpers.hpp
 	opm/core/grid/GridManager.hpp
 	opm/core/grid/cart_grid.h
 	opm/core/grid/cornerpoint_grid.h
@@ -299,10 +302,13 @@ list (APPEND PUBLIC_HEADER_FILES
 	opm/core/pressure/tpfa/compr_quant_general.h
 	opm/core/pressure/tpfa/compr_source.h
 	opm/core/pressure/tpfa/ifs_tpfa.h
+	opm/core/pressure/tpfa/TransTpfa.hpp
+	opm/core/pressure/tpfa/TransTpfa_impl.hpp
 	opm/core/pressure/tpfa/trans_tpfa.h
 	opm/core/props/BlackoilPhases.hpp
 	opm/core/props/BlackoilPropertiesBasic.hpp
 	opm/core/props/BlackoilPropertiesFromDeck.hpp
+	opm/core/props/BlackoilPropertiesFromDeck_impl.hpp
 	opm/core/props/BlackoilPropertiesInterface.hpp
 	opm/core/props/IncompPropertiesBasic.hpp
 	opm/core/props/IncompPropertiesFromDeck.hpp
@@ -395,6 +401,7 @@ list (APPEND PUBLIC_HEADER_FILES
 	opm/core/utility/have_boost_redef.hpp
 	opm/core/utility/linearInterpolation.hpp
 	opm/core/utility/miscUtilities.hpp
+	opm/core/utility/miscUtilities_impl.hpp
 	opm/core/utility/miscUtilitiesBlackoil.hpp
 	opm/core/utility/parameters/Parameter.hpp
 	opm/core/utility/parameters/ParameterGroup.hpp
@@ -413,4 +420,5 @@ list (APPEND PUBLIC_HEADER_FILES
 	opm/core/wells/WellCollection.hpp
 	opm/core/wells/WellsGroup.hpp
 	opm/core/wells/WellsManager.hpp
+	opm/core/wells/WellsManager_impl.hpp
 	)

opm/core/grid/GridHelpers.cpp

@@ -0,0 +1,81 @@
+#include "config.h"
+#include <opm/core/grid/GridHelpers.hpp>
+namespace Opm
+{
+namespace UgGridHelpers
+{
+int numCells(const UnstructuredGrid& grid)
+{
+    return grid.number_of_cells;
+}
+
+int numFaces(const UnstructuredGrid& grid)
+{
+    return grid.number_of_faces;
+}
+int dimensions(const UnstructuredGrid& grid)
+{
+    return grid.dimensions;
+}
+int numCellFaces(const UnstructuredGrid& grid)
+{
+    return grid.cell_facepos[grid.number_of_cells];
+}
+
+const int* globalCell(const UnstructuredGrid& grid)
+{
+    return grid.global_cell;
+}
+
+const int* cartDims(const UnstructuredGrid& grid)
+{
+    return grid.cartdims;
+}
+                    
+const double* beginCellCentroids(const UnstructuredGrid& grid)
+{
+    return grid.cell_centroids;
+}
+
+double cellCentroidCoordinate(const UnstructuredGrid& grid, int cell_index,
+                                 int coordinate)
+{
+    return grid.cell_centroids[grid.dimensions*cell_index+coordinate];
+}
+
+const double* beginFaceCentroids(const UnstructuredGrid& grid)
+{
+    return grid.face_centroids;
+}
+
+const double* faceCentroid(const UnstructuredGrid& grid, int face_index)
+{
+    return grid.face_centroids+face_index*grid.dimensions;
+}
+
+const double* faceNormal(const UnstructuredGrid& grid, int face_index)
+{
+    return grid.face_normals+face_index*grid.dimensions;
+}
+
+double faceArea(const UnstructuredGrid& grid, int face_index)
+{
+    return grid.face_areas[face_index];
+}
+
+SparseTableView cell2Faces(const UnstructuredGrid& grid)
+{
+    return SparseTableView(grid.cell_faces, grid.cell_facepos, numCells(grid));
+}
+
+double cellVolume(const UnstructuredGrid& grid, int cell_index)
+{
+    return grid.cell_volumes[cell_index];
+}
+
+FaceCellTraits<UnstructuredGrid>::Type faceCells(const UnstructuredGrid& grid)
+{
+    return FaceCellsProxy(grid);
+}
+}
+}

opm/core/grid/GridHelpers.hpp

@@ -0,0 +1,277 @@
+/*
+  Copyright 2014 Dr. Markus Blatt - HPC-Simulation-Software & Services
+  Copyright 2014 Statoil AS
+
+  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 <http://www.gnu.org/licenses/>.
+*/
+#ifndef OPM_CORE_GRIDHELPERS_HEADER_INCLUDED
+#define OPM_CORE_GRIDHELPERS_HEADER_INCLUDED
+#include <opm/core/grid.h>
+#include <boost/range/iterator_range.hpp>
+
+namespace Opm
+{
+namespace UgGridHelpers
+{
+
+/// \brief Allows viewing a sparse table consisting out of C-array 
+///
+/// This class can be used to convert two int array (like they are
+/// in UnstructuredGrid for representing the cell to faces mapping
+/// as a sparse table object.
+class SparseTableView
+{
+public:
+    class IntRange : public boost::iterator_range<const int*>
+    {
+    public:
+        typedef boost::iterator_range<const int*> BaseRowType;
+        typedef BaseRowType::size_type size_type;
+        typedef int value_type;
+        
+        IntRange(const int* start, const int* end)
+            : BaseRowType(start, end)
+        {}
+    };
+    /// \brief The type of the roww.
+    typedef boost::iterator_range<const int*> row_type;
+
+    /// \brief Creates a sparse table view
+    /// \param data The array with data of the table.
+    /// \param offset The offsets of the rows. Row i starts
+    ///               at offset[i] and ends a offset[i+1]
+    /// \param size   The number of entries/rows of the table
+    SparseTableView(int* data, int *offset, std::size_t size)
+        : data_(data), offset_(offset), size_(size)
+    {}
+
+    /// \brief Get a row of the the table.
+    /// \param row The row index.
+    /// \return The corresponding row.
+    row_type operator[](std::size_t row) const
+    {
+        assert(row<=size());
+        return row_type(data_ + offset_[row], data_ + offset_[row+1]);
+    }
+    
+    /// \brief Get the size of the table.
+    /// \return the number rows.
+    std::size_t size() const
+    {
+        return size_;
+    }
+    
+    /// \brief Get the number of non-zero entries.
+    std::size_t noEntries() const
+    {
+        return offset_[size_];
+    }
+
+private:
+    /// \brief The array with data of the table.
+    const int* data_;
+    /// \brief offset The offsets of the rows. 
+    ///
+    /// Row i starts at offset[i] and ends a offset[i+1]
+    const int* offset_;
+    /// \brief The size, i.e. the number of rows.
+    std::size_t size_;
+};
+
+/// \brief Get the number of cells of a grid.
+int numCells(const UnstructuredGrid& grid);
+
+/// \brief Get the number of faces of a grid.
+int numFaces(const UnstructuredGrid& grid);
+
+/// \brief Get the dimensions of a grid
+int dimensions(const UnstructuredGrid& grid);
+
+/// \brief Get the number of faces, where each face counts as many times as there are adjacent faces
+int numCellFaces(const UnstructuredGrid& grid);
+
+/// \brief Get the cartesion dimension of the underlying structured grid.
+const int* cartDims(const UnstructuredGrid& grid);
+
+/// \brief Get the local to global index mapping.
+///
+/// The global index is the index of the active cell
+/// in the underlying structured grid.
+const int* globalCell(const UnstructuredGrid& grid);
+
+/// \brief Traits of the cell centroids of a grid.
+///
+/// This class exports two types: IteratorType, the type of the iterator
+/// over the cell centroids, and the ValueTpe, the type of the cell centroid.
+/// \tpatam G The type of the grid.
+template<class G>
+struct CellCentroidTraits
+{
+};
+
+template<>
+struct CellCentroidTraits<UnstructuredGrid>
+{
+    typedef const double* IteratorType;
+    typedef const double* ValueType;
+};
+
+/// \brief Get an iterator over the cell centroids positioned at the first cell.
+///
+/// The return type needs to be usable with the functions increment, and
+/// getCoordinate.
+CellCentroidTraits<UnstructuredGrid>::IteratorType
+beginCellCentroids(const UnstructuredGrid& grid);
+
+/// \brief Get a coordinate of a specific cell centroid.
+/// \brief grid The grid.
+/// \brief cell_index The index of the specific cell.
+/// \breif coordinate The coordinate index.
+double cellCentroidCoordinate(const UnstructuredGrid& grid, int cell_index,
+                                 int coordinate);
+
+/// \brief Traits of the face centroids of a grid.
+///
+/// This class exports two types: IteratorType, the type of the iterator
+/// over the face centroids, and the ValueTpe, the type of the face centroid.
+/// \tpatam G The type of the grid.
+template<class G>
+struct FaceCentroidTraits
+{
+};
+
+template<>
+struct FaceCentroidTraits<UnstructuredGrid>
+{
+    typedef const double* IteratorType;
+    typedef const double* ValueType;
+};
+
+/// \brief Get an iterator over the face centroids positioned at the first cell.
+FaceCentroidTraits<UnstructuredGrid>::IteratorType
+beginFaceCentroids(const UnstructuredGrid& grid);
+
+/// \brief Get a coordinate of a specific face centroid.
+/// \param grid The grid.
+/// \param face_index The index of the specific face.
+/// \param coordinate The coordinate index.
+FaceCentroidTraits<UnstructuredGrid>::ValueType
+faceCentroid(const UnstructuredGrid& grid, int face_index);
+
+/// \brief Get the normal of a face.
+/// \param grid The grid that the face is part of.
+/// \param face_index The index of the face in the grid.
+const double* faceNormal(const UnstructuredGrid& grid, int face_index);
+
+/// \brief Get the area of a face
+/// \param grid The grid that the face is part of.
+/// \param face_index The index of the face in the grid.
+double faceArea(const UnstructuredGrid& grid, int face_index);
+
+/// \brief Maps the grid type to the associated type of the cell to faces mapping.
+///
+/// Provides a type named Type.
+/// \tparam T The type of the grid.
+template<class T>
+struct Cell2FacesTraits
+{
+};
+
+template<>
+struct Cell2FacesTraits<UnstructuredGrid>
+{
+    typedef SparseTableView Type;
+};
+
+/// \brief Get the cell to faces mapping of a grid.
+Cell2FacesTraits<UnstructuredGrid>::Type 
+cell2Faces(const UnstructuredGrid& grid);
+
+class FaceCellsProxy
+{
+public:
+    FaceCellsProxy(const UnstructuredGrid& grid)
+    : face_cells_(grid.face_cells)
+    {}
+    int operator()(int cell_index, int local_index)
+    {
+        return face_cells_[2*cell_index+local_index];
+    }
+private:
+    const int* face_cells_;
+};
+
+/// \brief Traits of the face to attached cell mappping of a grid.
+///
+/// Exports the type Type, the type of the mapping
+/// \tparam T The type of the grid
+template<class T>
+struct FaceCellTraits
+{};
+
+template<>
+struct FaceCellTraits<UnstructuredGrid>
+{
+    typedef FaceCellsProxy Type;
+};
+
+/// \brief Get the face to cell mapping of a grid.
+FaceCellTraits<UnstructuredGrid>::Type faceCells(const UnstructuredGrid& grid);
+
+/// \brief Increment an iterator over an array that reresents a dense row-major
+///  matrix with dims columns
+/// \param cc The iterator.
+/// \param i The nzumber of rows to increment
+/// \param dim The number of columns of the matrix.
+template<class T>
+T* increment(T* cc, int i, int dim)
+{
+    return cc+(i*dim);
+}
+/// \brief Increment an iterator over an array that reresents a dense row-major
+///  matrix with dims columns
+/// \param cc The iterator.
+/// \param i The nzumber of rows to increment
+template<class T>
+T increment(const T& t, int i, int)
+{
+    return t+i;
+}
+
+/// \brief Get the i-th corrdinate of a centroid.
+/// \param cc The array with the coordinates.
+/// \param i The index of the coordinate.
+/// \tparam T The type of the coordinate of the centroid.
+template<class T>
+double getCoordinate(T* cc, int i)
+{
+    return cc[i];
+}
+
+/// \brief Get the i-th corrdinate of an array.
+/// \param t The iterator over the centroids
+/// \brief i The index of the coordinate.
+/// \tparam T The type of the iterator representing the centroid.
+/// Its value_type has to provide an operator[] to access the coordinates.
+template<class T>
+double getCoordinate(T t, int i)
+{
+    return (*t)[i];
+}
+
+} // end namespace UGGridHelpers
+} // end namespace OPM
+#endif

opm/core/io/eclipse/EclipseWriter.cpp

@@ -184,26 +184,27 @@ void restrictToActiveCells_(std::vector<double> &data, const std::vector<int> &a
 
 // throw away the data for all non-active cells in an array. (this is
 // the variant of the function which takes an UnstructuredGrid object.)
-void restrictToActiveCells_(std::vector<double> &data, const UnstructuredGrid &grid)
+void restrictToActiveCells_(std::vector<double> &data, int number_of_cells,
+                            const int* global_cell)
 {
-    if (!grid.global_cell)
+    if (!global_cell)
         // if there is no active -> global mapping, all cells
         // are considered active
         return;
 
     // activate those cells that are actually there
-    for (int i = 0; i < grid.number_of_cells; ++i) {
+    for (int i = 0; i < number_of_cells; ++i) {
         // make sure that global cell indices are always at least as
         // large as the active one and that the global cell indices
         // are in increasing order. the latter might become
         // problematic if cells are extensively re-ordered, but that
         // does not seem to be the case so far
-        assert(grid.global_cell[i] >= i);
-        assert(i == 0 || grid.global_cell[i - 1] < grid.global_cell[i]);
+        assert(global_cell[i] >= i);
+        assert(i == 0 || global_cell[i - 1] < global_cell[i]);
 
-        data[i] = data[grid.global_cell[i]];
+        data[i] = data[global_cell[i]];
     }
-    data.resize(grid.number_of_cells);
+    data.resize(number_of_cells);
 }
 
 // convert the units of an array
@@ -232,23 +233,25 @@ void extractFromStripedData_(std::vector<double> &data,
 }
 
 // enclosure of the current grid in a Cartesian space
-int getCartesianSize_(const UnstructuredGrid& grid) {
-    const int nx = grid.cartdims[0];
-    const int ny = grid.cartdims[1];
-    const int nz = grid.cartdims[2];
+int getCartesianSize_(const int* cartdims) {
+    const int nx = cartdims[0];
+    const int ny = cartdims[1];
+    const int nz = cartdims[2];
     return nx * ny * nz;
 }
 
-void getActiveCells_(const UnstructuredGrid& grid,
+void getActiveCells_(int number_of_cells,
+                     const int* cartdims,
+                     const int* global_cell,
                      std::vector <int>& actnum)
 {
     // we must fill the Cartesian grid with flags
-    const int size = getCartesianSize_(grid);
+    const int size = getCartesianSize_(cartdims);
 
     // if we don't have a global_cells field, then assume that all
     // grid cells is active
-    if (!grid.global_cell) {
-        if (grid.number_of_cells != size) {
+    if (!global_cell) {
+        if (number_of_cells != size) {
             OPM_THROW (std::runtime_error,
                        "No ACTNUM map but grid size != Cartesian size");
         }
@@ -259,8 +262,8 @@ void getActiveCells_(const UnstructuredGrid& grid,
         actnum.assign (size, 0);
 
         // activate those cells that are actually there
-        for (int i = 0; i < grid.number_of_cells; ++i) {
-            actnum[grid.global_cell[i]] = 1;
+        for (int i = 0; i < number_of_cells; ++i) {
+            actnum[global_cell[i]] = 1;
         }
     }
 }
@@ -503,7 +506,9 @@ private:
 struct EclipseGrid : public EclipseHandle <ecl_grid_type> {
     /// Create a grid based on the keywords available in input file
     static EclipseGrid make (Opm::DeckConstPtr newParserDeck,
-                             const UnstructuredGrid& grid)
+                             int number_of_cells,
+                             const int* cart_dims,
+                             const int* global_cell)
     {
         if (newParserDeck->hasKeyword("DXV")) {
             // make sure that the DYV and DZV keywords are present if the
@@ -528,7 +533,7 @@ struct EclipseGrid : public EclipseHandle <ecl_grid_type> {
 
             // get the actually active cells, after processing
             std::vector <int> actnum;
-            getActiveCells_(grid, actnum);
+            getActiveCells_(number_of_cells, cart_dims, global_cell, actnum);
             EclipseKeyword<int> actnum_kw (ACTNUM_KW, actnum);
 
             EclipseKeyword<float> mapaxes_kw (MAPAXES_KW);
@@ -537,7 +542,7 @@ struct EclipseGrid : public EclipseHandle <ecl_grid_type> {
                 mapaxes_kw = std::move (EclipseKeyword<float> (MAPAXES_KW, mapaxesData));
             }
 
-            return EclipseGrid (g.dims, zcorn_kw, coord_kw, actnum_kw, mapaxes_kw);
+            return EclipseGrid (cart_dims, zcorn_kw, coord_kw, actnum_kw, mapaxes_kw);
         }
         else {
             OPM_THROW(std::runtime_error,
@@ -629,15 +634,18 @@ struct EclipseInit : public EclipseHandle <fortio_type> {
         return EclipseInit (initFileName, fmt_file);
     }
 
-    void writeHeader (const UnstructuredGrid& grid,
+    void writeHeader (int number_of_cells,
+                      const int* cart_dims,
+                      const int* global_cell,
                       const SimulatorTimer& timer,
                       Opm::DeckConstPtr newParserDeck,
                       const PhaseUsage uses)
     {
         auto dataField = getAllSiDoubles_(newParserDeck->getKeyword(PORO_KW));
-        restrictToActiveCells_(dataField, grid);
+        restrictToActiveCells_(dataField, number_of_cells, global_cell);
 
-        EclipseGrid eclGrid = EclipseGrid::make (newParserDeck, grid);
+        EclipseGrid eclGrid = EclipseGrid::make (newParserDeck, number_of_cells,
+                                                 cart_dims, global_cell);
 
         EclipseKeyword<float> poro (PORO_KW, dataField);
         ecl_init_file_fwrite_header (*this,
@@ -1046,11 +1054,14 @@ void EclipseWriter::writeInit(const SimulatorTimer &timer)
         return;
     }
     /* Grid files */
-    EclipseGrid eclGrid = EclipseGrid::make (newParserDeck_, *grid_);
+    EclipseGrid eclGrid = EclipseGrid::make (newParserDeck_, number_of_cells_,
+                                             cart_dims_, global_cell_);
     eclGrid.write (outputDir_, baseName_, /*stepIdx=*/0);
 
     EclipseInit fortio = EclipseInit::make (outputDir_, baseName_, /*stepIdx=*/0);
-    fortio.writeHeader (*grid_,
+    fortio.writeHeader (number_of_cells_,
+                        cart_dims_,
+                        global_cell_,
                         timer,
                         newParserDeck_,
                         uses_);
@@ -1178,12 +1189,34 @@ void EclipseWriter::writeTimeStep(
 
 #endif // HAVE_ERT
 
-EclipseWriter::EclipseWriter(const ParameterGroup& params,
+EclipseWriter::EclipseWriter (
+        const ParameterGroup& params,
+        Opm::DeckConstPtr newParserDeck,
+        int number_of_cells, const int* global_cell, const int* cart_dims,
+        int dimensions)
+    : newParserDeck_ (newParserDeck)
+    , number_of_cells_(number_of_cells)
+    , dimensions_(dimensions)
+    , cart_dims_(cart_dims)
+    , global_cell_(global_cell)
+    , uses_ (phaseUsageFromDeck (newParserDeck_)) {
+    init(params);
+}
+
+EclipseWriter::EclipseWriter (
+        const ParameterGroup& params,
         Opm::DeckConstPtr newParserDeck,
         std::shared_ptr<const UnstructuredGrid> grid)
-    : newParserDeck_(newParserDeck)
+    : newParserDeck_ (newParserDeck)
     , grid_(grid)
-    , uses_(phaseUsageFromDeck(newParserDeck))
+    , number_of_cells_(grid->number_of_cells)
+    , dimensions_(grid->dimensions)
+    , cart_dims_(grid->cartdims)
+    , global_cell_(grid->global_cell)
+    , uses_ (phaseUsageFromDeck (newParserDeck_)) {
+    init(params);
+}
+void EclipseWriter::init(const ParameterGroup& params)
 {
     // get the base name from the name of the deck
     using boost::filesystem::path;

opm/core/io/eclipse/EclipseWriter.hpp

@@ -62,6 +62,17 @@ public:
                   Opm::DeckConstPtr newParserDeck,
                   std::shared_ptr <const UnstructuredGrid> grid);
 
+    /*!
+     * \brief Sets the common attributes required to write eclipse
+     *        binary files using ERT.
+     */
+    EclipseWriter(const parameter::ParameterGroup& params,
+                  Opm::DeckConstPtr newParserDeck,
+                  int number_of_cells, const int* global_cell, const int* cart_dims,
+                  int dimension);
+
+
+                 
     /**
      * We need a destructor in the compilation unit to avoid the
      * EclipseSummary being a complete type here.
@@ -94,6 +105,10 @@ public:
 private:
     Opm::DeckConstPtr newParserDeck_;
     std::shared_ptr <const UnstructuredGrid> grid_;
+    int number_of_cells_;
+    int dimensions_;
+    const int* cart_dims_;
+    const int* global_cell_;
     bool enableOutput_;
     int outputInterval_;
     int outputTimeStepIdx_;
@@ -101,6 +116,8 @@ private:
     std::string baseName_;
     PhaseUsage uses_;           // active phases in the input deck
     std::shared_ptr <EclipseSummary> summary_;
+
+    void init(const parameter::ParameterGroup& params);
 };
 } // namespace Opm
 

opm/core/pressure/tpfa/TransTpfa.cpp

@@ -0,0 +1,10 @@
+#include "TransTpfa.hpp"
+#include <opm/core/grid.h>
+
+/* Ecplicitly initialize UnstructuredGrid versions */
+
+template void tpfa_htrans_compute(const UnstructuredGrid*, const double*, double*);
+
+template void tpfa_trans_compute(const UnstructuredGrid*, const double*, double*);
+
+template void tpfa_eff_trans_compute(const UnstructuredGrid*, const double*, const double*, double*);

opm/core/pressure/tpfa/TransTpfa.hpp

@@ -0,0 +1,105 @@
+/*
+  Copyright 2010 SINTEF ICT, Applied Mathematics.
+  Copyright 2014 Dr. Blatt - HPC-Simulation-Software & Services
+
+  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 <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_TRANSTPFA_HEADER_INCLUDED
+#define OPM_TRANSTPFA_HEADER_INCLUDED
+
+/**
+ * \file
+ * Routines to assist in the calculation of two-point transmissibilities.
+ */
+
+/**
+ * Calculate static, one-sided transmissibilities for use in the two-point flux
+ * approximation method.
+ *
+ * The one-sided transmissibilities are defined by the formula
+ * \f[
+ * t_i = \frac{\vec{n}_f \mathsf{K}_c \vec{c}_c}{\lVert \vec{c}_c \rVert^2}
+ * \f]
+ * in which @c i is the half-face index corresponding to the cell-face index
+ * pair <CODE>(c,f)</CODE> and \f$\vec{c}_{cf} = \Bar{x}_f - \Bar{x}_c\f$ is the
+ * centroid difference vector.
+ *
+ * @param[in]  G       Grid.
+ * @param[in]  perm    Permeability.  One symmetric, positive definite tensor
+ *                     per grid cell.
+ * @param[out] htrans  One-sided transmissibilities.  Array of size at least
+ *                     <CODE>G->cell_facepos[ G->number_of_cells  ]</CODE>.
+ */
+template<class Grid>
+void
+tpfa_htrans_compute(const Grid   *G      ,
+                    const double *perm  ,
+                    double       *htrans);
+
+/**
+ * Compute two-point transmissibilities from one-sided transmissibilities.
+ *
+ * The two-point transmissibilities are given by the simple, near-harmonic
+ * average (save a factor of two)
+ * \f[
+ * \mathsf{T}_f = \big(\frac{1}{t_1} + \frac{1}{t_2}\big)^{-1}
+ *              = \frac{t_1t_2}{t_1 + t_2}
+ * \f]
+ * in which \f$t_1\f$ and \f$t_2\f$ are the one-sided transmissibilities that
+ * connect the neighbouring cells of face @c f.
+ *
+ * @param[in]  G       Grid.
+ * @param[in]  htrans  One-sided transmissibilities as defined by function
+ *                     tpfa_htrans_compute().
+ * @param[out] trans   Interface, two-point transmissibilities.  Array of size
+ *                     at least <CODE>numFaces(G)</CODE>.
+ */
+template<class Grid>
+void
+tpfa_trans_compute(const Grid   *G      ,
+                   const double *htrans,
+                   double       *trans );
+
+/**
+ * Calculate effective two-point transmissibilities from one-sided, total
+ * mobility weighted, transmissibilities.
+ *
+ * Specifically, compute the following product
+ * \f[
+ * \mathsf{T}_f = \big(\frac{1}{\lambda_1t_1} + \frac{1}{\lambda_2t_2}\big)^{-1}
+ *              = \lambda_1\lambda_2 \frac{t_1t_2}{t_1 + t_2}
+ * \f]
+ * in which \f$t_1\f$ and \f$t_2\f$ are the one-sided, static transmissibility
+ * values connecting the cells of face @c f and \f$\lambda_1\f$ and
+ * \f$\lambda_2\f$ denote the total mobilities of the respective cells.
+ *
+ * @param[in]  G      Grid.
+ * @param[in]  totmob Total mobilities. One positive scalar value for each cell.
+ * @param[in]  htrans One-sided transmissibilities as defined by function
+ *                    tpfa_htrans_compute().
+ * @param[out] trans  Effective, two-point transmissibilities.  Array of size at
+ *                    least <CODE>numFaces(G)</CODE>.
+ */
+template<class Grid>
+void
+tpfa_eff_trans_compute(const Grid   *G      ,
+                       const double *totmob,
+                       const double *htrans,
+                       double       *trans );
+
+#include "TransTpfa_impl.hpp"
+#endif  /* OPM_TRANS_TPFA_HEADER_INCLUDED */

opm/core/pressure/tpfa/TransTpfa_impl.hpp

@@ -0,0 +1,177 @@
+#include <opm/core/linalg/blas_lapack.h>
+#include <opm/core/pressure/tpfa/trans_tpfa.h>
+#include <opm/core/grid/GridHelpers.hpp>
+
+#include <cmath>
+
+namespace Dune
+{
+class CpGrid;
+}
+
+namespace Opm
+{
+namespace UgGridHelpers
+{
+int dimensions(const Dune::CpGrid&);
+
+double faceArea(const Dune::CpGrid&, int);
+}
+}
+
+namespace
+{
+#ifdef HAVE_DUNE_CORNERPOINT
+const double* multiplyFaceNormalWithArea(const Dune::CpGrid& grid, int face_index, const double* in)
+{
+    int d=Opm::UgGridHelpers::dimensions(grid);
+    double* out=new double[d];
+    double area=Opm::UgGridHelpers::faceArea(grid, face_index);
+    
+    for(int i=0;i<d;++i)
+        out[i]=in[i]*area;
+    return out;
+}
+
+inline void maybeFreeFaceNormal(const Dune::CpGrid&, const double* array)
+{
+    delete[] array;
+}
+#endif  // HAVE_DUNE_CORNERPOINT
+
+inline const double* multiplyFaceNormalWithArea(const UnstructuredGrid&, int, const double* in)
+{
+    return in;
+}
+
+inline void maybeFreeFaceNormal(const UnstructuredGrid&, const double*)
+{}
+}
+
+/* ---------------------------------------------------------------------- */
+/* htrans <- sum(C(:,i) .* K(cellNo,:) .* N(:,j), 2) ./ sum(C.*C, 2) */
+/* ---------------------------------------------------------------------- */
+template<class Grid>
+void
+tpfa_htrans_compute(const Grid* G, const double *perm, double *htrans)
+/* ---------------------------------------------------------------------- */
+{
+    using namespace Opm::UgGridHelpers;
+    int    d, j;
+    double s, dist, denom;
+
+    double Kn[3];
+    typename CellCentroidTraits<Grid>::IteratorType cc = beginCellCentroids(*G);
+    typename Cell2FacesTraits<Grid>::Type c2f = cell2Faces(*G);
+    typename FaceCellTraits<Grid>::Type face_cells = faceCells(*G);
+    
+    const double *n;
+    const double *K;
+
+    MAT_SIZE_T nrows, ncols, ldA, incx, incy;
+    double a1, a2;
+
+    d = dimensions(*G);
+
+    nrows = ncols    = ldA = d;
+    incx  = incy     = 1      ;
+    a1    = 1.0;  a2 = 0.0    ;
+
+    for (int c =0, i = 0; c < numCells(*G); c++) {
+        K  = perm + (c * d * d);
+        
+        typedef typename Cell2FacesTraits<Grid>::Type::row_type FaceRow;
+        FaceRow faces = c2f[c];
+        
+        for(typename FaceRow::const_iterator f=faces.begin(), end=faces.end();
+            f!=end; ++f, ++i)
+        {
+            s = 2.0*(face_cells(*f, 0) == c) - 1.0;
+            n = faceNormal(*G, *f);
+            const double* nn=multiplyFaceNormalWithArea(*G, *f, n);
+            const double* fc = &(faceCentroid(*G, *f)[0]);
+            dgemv_("No Transpose", &nrows, &ncols,
+                   &a1, K, &ldA, nn, &incx, &a2, &Kn[0], &incy);
+            maybeFreeFaceNormal(*G, nn);
+            
+            htrans[i] = denom = 0.0;
+            for (j = 0; j < d; j++) {
+                dist = fc[j] - getCoordinate(cc, j);
+
+                htrans[i] += s * dist * Kn[j];
+                denom     +=     dist * dist;
+            }
+
+            assert (denom > 0);
+            htrans[i] /= denom;
+            htrans[i]  = std::abs(htrans[i]);
+        }
+        // Move to next cell centroid.
+        cc = increment(cc, 1, d);
+    }
+}
+
+
+/* ---------------------------------------------------------------------- */
+template<class Grid>
+void
+tpfa_trans_compute(const Grid* G, const double *htrans, double *trans)
+/* ---------------------------------------------------------------------- */
+{
+    using namespace Opm::UgGridHelpers;
+
+    for (int f = 0; f < numFaces(*G); f++) {
+        trans[f] = 0.0;
+    }
+
+    typename Cell2FacesTraits<Grid>::Type c2f = cell2Faces(*G);
+
+    for (int c = 0, i = 0; c < numCells(*G); c++) {
+        typedef typename Cell2FacesTraits<Grid>::Type::row_type FaceRow;
+        FaceRow faces = c2f[c];
+        
+        for(typename FaceRow::const_iterator f=faces.begin(), end=faces.end();
+            f!=end; ++f, ++i)
+        {
+            trans[*f] += 1.0 / htrans[i];
+        }
+    }
+
+    for (int f = 0; f < numFaces(*G); f++) {
+        trans[f] = 1.0 / trans[f];
+    }
+}
+
+
+/* ---------------------------------------------------------------------- */
+ template<class Grid>
+void
+tpfa_eff_trans_compute(const Grid*        G,
+                       const double *totmob,
+                       const double *htrans,
+                       double       *trans)
+/* ---------------------------------------------------------------------- */
+{
+    using namespace Opm::UgGridHelpers;
+
+    for (int f = 0; f < numFaces(*G); f++) {
+        trans[f] = 0.0;
+    }
+
+    typename Cell2FacesTraits<Grid>::Type c2f = cell2Faces(*G);
+
+    for (int c = 0, i = 0; c < numCells(*G); c++) {
+        typedef typename Cell2FacesTraits<Grid>::Type::row_type FaceRow;
+        FaceRow faces = c2f[c];
+        
+        for(typename FaceRow::const_iterator f=faces.begin(), end=faces.end();
+            f!=end; ++f, ++i)
+        {
+            trans[*f] += 1.0 / (totmob[c] * htrans[i]);
+        }
+    }
+
+             for (int f = 0; f < numFaces(*G); f++) {
+        trans[f] = 1.0 / trans[f];
+    }
+}

opm/core/pressure/tpfa/trans_tpfa.c

@@ -7,6 +7,11 @@
 #include <opm/core/linalg/blas_lapack.h>
 #include <opm/core/pressure/tpfa/trans_tpfa.h>
 
+
+#ifdef __cplusplus
+#include "TransTpfa.hpp"
+#endif
+
 /* ---------------------------------------------------------------------- */
 /* htrans <- sum(C(:,i) .* K(cellNo,:) .* N(:,j), 2) ./ sum(C.*C, 2) */
 /* ---------------------------------------------------------------------- */
@@ -14,6 +19,10 @@ void
 tpfa_htrans_compute(struct UnstructuredGrid *G, const double *perm, double *htrans)
 /* ---------------------------------------------------------------------- */
 {
+    #ifdef __cplusplus
+    return tpfa_htrans_compute<UnstructuredGrid>(G, totmob, htrans, trans);
+    #endif
+    
     int    c, d, f, i, j;
     double s, dist, denom;
 
@@ -65,6 +74,10 @@ void
 tpfa_trans_compute(struct UnstructuredGrid *G, const double *htrans, double *trans)
 /* ---------------------------------------------------------------------- */
 {
+    #ifdef __cplusplus
+    return tpfa_trans_compute<UnstructuredGrid>(G, totmob, htrans, trans);
+    #endif
+    
     int c, i, f;
 
     for (f = 0; f < G->number_of_faces; f++) {
@@ -93,6 +106,10 @@ tpfa_eff_trans_compute(struct UnstructuredGrid       *G,
                        double       *trans)
 /* ---------------------------------------------------------------------- */
 {
+    #ifdef __cplusplus
+    return tpfa_eff_trans_compute<UnstructuredGrid>(G, totmob, htrans, trans);
+    #endif
+    
     int c, i, f;
 
     for (f = 0; f < G->number_of_faces; f++) {

opm/core/props/BlackoilPropertiesFromDeck.cpp

@@ -27,38 +27,17 @@ namespace Opm
                                                            const UnstructuredGrid& grid,
                                                            bool init_rock)
     {
-       if (init_rock){
-            rock_.init(deck, grid);
+        init(deck, grid.number_of_cells, grid.global_cell, grid.cartdims, grid.cell_centroids,
+             grid.dimensions, init_rock);
     }
-        pvt_.init(deck, 0);
-        SaturationPropsFromDeck<SatFuncSimpleNonuniform>* ptr
-            = new SaturationPropsFromDeck<SatFuncSimpleNonuniform>();
-        satprops_.reset(ptr);
-        ptr->init(deck, grid, 0);
 
-        if (pvt_.numPhases() != satprops_->numPhases()) {
-            OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data ("
-                  << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ").");
-        }
-    }
 
     BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck,
                                                            const UnstructuredGrid& grid,
                                                            bool init_rock)
     {
-        if (init_rock){
-           rock_.init(newParserDeck, grid);
-        }
-        pvt_.init(newParserDeck, /*numSamples=*/0);
-        SaturationPropsFromDeck<SatFuncSimpleNonuniform>* ptr
-            = new SaturationPropsFromDeck<SatFuncSimpleNonuniform>();
-        satprops_.reset(ptr);
-        ptr->init(newParserDeck, grid, /*numSamples=*/0);
-
-        if (pvt_.numPhases() != satprops_->numPhases()) {
-            OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data ("
-                  << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ").");
-        }
+        init(newParserDeck, grid.number_of_cells, grid.global_cell, grid.cartdims,
+             grid.cell_centroids, grid.dimensions, init_rock);
     }
 
     BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
@@ -66,63 +45,8 @@ namespace Opm
                                                            const parameter::ParameterGroup& param,
                                                            bool init_rock)
     {
-        if(init_rock){
-            rock_.init(deck, grid);
-        }
-
-        const int pvt_samples = param.getDefault("pvt_tab_size", 0);
-        pvt_.init(deck, pvt_samples);
-
-        // Unfortunate lack of pointer smartness here...
-        const int sat_samples = param.getDefault("sat_tab_size", 0);
-        std::string threephase_model = param.getDefault<std::string>("threephase_model", "simple");
-        if (deck.hasField("ENDSCALE") && threephase_model != "gwseg") {
-            OPM_THROW(std::runtime_error, "Sorry, end point scaling currently available for the 'gwseg' model only.");
-        }
-        if (sat_samples > 1) {
-            if (threephase_model == "stone2") {
-                SaturationPropsFromDeck<SatFuncStone2Uniform>* ptr
-                    = new SaturationPropsFromDeck<SatFuncStone2Uniform>();
-                satprops_.reset(ptr);
-                ptr->init(deck, grid, sat_samples);
-            } else if (threephase_model == "simple") {
-                SaturationPropsFromDeck<SatFuncSimpleUniform>* ptr
-                    = new SaturationPropsFromDeck<SatFuncSimpleUniform>();
-                satprops_.reset(ptr);
-                ptr->init(deck, grid, sat_samples);
-            } else if (threephase_model == "gwseg") {
-                SaturationPropsFromDeck<SatFuncGwsegUniform>* ptr
-                    = new SaturationPropsFromDeck<SatFuncGwsegUniform>();
-                satprops_.reset(ptr);
-                ptr->init(deck, grid, sat_samples);
-            } else {
-                OPM_THROW(std::runtime_error, "Unknown threephase_model: " << threephase_model);
-            }
-        } else {
-            if (threephase_model == "stone2") {
-                SaturationPropsFromDeck<SatFuncStone2Nonuniform>* ptr
-                    = new SaturationPropsFromDeck<SatFuncStone2Nonuniform>();
-                satprops_.reset(ptr);
-                ptr->init(deck, grid, sat_samples);
-            } else if (threephase_model == "simple") {
-                SaturationPropsFromDeck<SatFuncSimpleNonuniform>* ptr
-                    = new SaturationPropsFromDeck<SatFuncSimpleNonuniform>();
-                satprops_.reset(ptr);
-                ptr->init(deck, grid, sat_samples);
-            } else if (threephase_model == "gwseg") {
-                SaturationPropsFromDeck<SatFuncGwsegNonuniform>* ptr
-                    = new SaturationPropsFromDeck<SatFuncGwsegNonuniform>();
-                satprops_.reset(ptr);
-                ptr->init(deck, grid, sat_samples);
-            } else {
-                OPM_THROW(std::runtime_error, "Unknown threephase_model: " << threephase_model);
-            }
-        }
-
-        if (pvt_.numPhases() != satprops_->numPhases()) {
-            OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data ("
-                  << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ").");
-        }
+        init(deck, grid.number_of_cells, grid.global_cell, grid.cartdims, grid.cell_centroids, 
+             grid.dimensions, param, init_rock);
     }
 
     BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck,
@@ -130,63 +54,8 @@ namespace Opm
                                                            const parameter::ParameterGroup& param,
                                                            bool init_rock)
     {
-        if(init_rock){
-            rock_.init(newParserDeck, grid);
-        }
-
-        const int pvt_samples = param.getDefault("pvt_tab_size", 0);
-        pvt_.init(newParserDeck, pvt_samples);
-
-        // Unfortunate lack of pointer smartness here...
-        const int sat_samples = param.getDefault("sat_tab_size", 0);
-        std::string threephase_model = param.getDefault<std::string>("threephase_model", "gwseg");
-        if (newParserDeck->hasKeyword("ENDSCALE") && threephase_model != "gwseg") {
-            OPM_THROW(std::runtime_error, "Sorry, end point scaling currently available for the 'gwseg' model only.");
-        }
-        if (sat_samples > 1) {
-            if (threephase_model == "stone2") {
-                SaturationPropsFromDeck<SatFuncStone2Uniform>* ptr
-                    = new SaturationPropsFromDeck<SatFuncStone2Uniform>();
-                satprops_.reset(ptr);
-                ptr->init(newParserDeck, grid, sat_samples);
-            } else if (threephase_model == "simple") {
-                SaturationPropsFromDeck<SatFuncSimpleUniform>* ptr
-                    = new SaturationPropsFromDeck<SatFuncSimpleUniform>();
-                satprops_.reset(ptr);
-                ptr->init(newParserDeck, grid, sat_samples);
-            } else if (threephase_model == "gwseg") {
-                SaturationPropsFromDeck<SatFuncGwsegUniform>* ptr
-                    = new SaturationPropsFromDeck<SatFuncGwsegUniform>();
-                satprops_.reset(ptr);
-                ptr->init(newParserDeck, grid, sat_samples);
-            } else {
-                OPM_THROW(std::runtime_error, "Unknown threephase_model: " << threephase_model);
-            }
-        } else {
-            if (threephase_model == "stone2") {
-                SaturationPropsFromDeck<SatFuncStone2Nonuniform>* ptr
-                    = new SaturationPropsFromDeck<SatFuncStone2Nonuniform>();
-                satprops_.reset(ptr);
-                ptr->init(newParserDeck, grid, sat_samples);
-            } else if (threephase_model == "simple") {
-                SaturationPropsFromDeck<SatFuncSimpleNonuniform>* ptr
-                    = new SaturationPropsFromDeck<SatFuncSimpleNonuniform>();
-                satprops_.reset(ptr);
-                ptr->init(newParserDeck, grid, sat_samples);
-            } else if (threephase_model == "gwseg") {
-                SaturationPropsFromDeck<SatFuncGwsegNonuniform>* ptr
-                    = new SaturationPropsFromDeck<SatFuncGwsegNonuniform>();
-                satprops_.reset(ptr);
-                ptr->init(newParserDeck, grid, sat_samples);
-            } else {
-                OPM_THROW(std::runtime_error, "Unknown threephase_model: " << threephase_model);
-            }
-        }
-
-        if (pvt_.numPhases() != satprops_->numPhases()) {
-            OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data ("
-                  << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ").");
-        }
+        init(newParserDeck, grid.number_of_cells, grid.global_cell, grid.cartdims, grid.cell_centroids, 
+             grid.dimensions, param, init_rock);
     }
 
     BlackoilPropertiesFromDeck::~BlackoilPropertiesFromDeck()

opm/core/props/BlackoilPropertiesFromDeck.hpp

@@ -90,6 +90,45 @@ namespace Opm
                                    const parameter::ParameterGroup& param,
                                    bool init_rock=true);
 
+        template<class T>
+        BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
+                                   int number_of_cells,
+                                   const int* global_cell,
+                                   const int* cart_dims,
+                                   T begin_cell_centroids,
+                                   int dimension,
+                                   bool init_rock=true);
+
+
+        template<class T>
+        BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
+                                   int number_of_cells,
+                                   const int* global_cell,
+                                   const int* cart_dims,
+                                   T begin_cell_centroids,
+                                   int dimension,
+                                   const parameter::ParameterGroup& param,
+                                   bool init_rock=true);
+
+        template<class T>
+        BlackoilPropertiesFromDeck(Opm::DeckConstPtr  newParserDeck,
+                                   int number_of_cells,
+                                   const int* global_cell,
+                                   const int* cart_dims,
+                                   T begin_cell_centroids,
+                                   int dimension,
+                                   bool init_rock=true);
+
+        template<class T>
+        BlackoilPropertiesFromDeck(Opm::DeckConstPtr  newParserDeck,
+                                   int number_of_cells,
+                                   const int* global_cell,
+                                   const int* cart_dims,
+                                   T begin_cell_centroids,
+                                   int dimension,
+                                   const parameter::ParameterGroup& param,
+                                   bool init_rock=true);
+
         /// Destructor.
         virtual ~BlackoilPropertiesFromDeck();
 
@@ -211,6 +250,41 @@ namespace Opm
                               double* smax) const;
 
     private:
+        template<class T>
+        void init(const EclipseGridParser& deck,
+                  int number_of_cells,
+                  const int* global_cell,
+                  const int* cart_dims,
+                  T begin_cell_centroids,
+                  int dimension,
+                  bool init_rock);
+        template<class T>
+        void init(const EclipseGridParser& deck,
+                  int number_of_cells,
+                  const int* global_cell,
+                  const int* cart_dims,
+                  T begin_cell_centroids,
+                  int dimension,
+                  const parameter::ParameterGroup& param,
+                  bool init_rock);
+
+        template<class T>
+        void init(Opm::DeckConstPtr  newParserDeck,
+                  int number_of_cells,
+                  const int* global_cell,
+                  const int* cart_dims,
+                  T begin_cell_centroids,
+                  int dimension,
+                  bool init_rock);
+        template<class T>
+        void init(Opm::DeckConstPtr  newParserDeck,
+                  int number_of_cells,
+                  const int* global_cell,
+                  const int* cart_dims,
+                  T begin_cell_centroids,
+                  int dimension,
+                  const parameter::ParameterGroup& param,
+                  bool init_rock);
         RockFromDeck rock_;
         BlackoilPvtProperties pvt_;
         std::unique_ptr<SaturationPropsInterface> satprops_;
@@ -224,5 +298,6 @@ namespace Opm
 
 } // namespace Opm
 
+#include "BlackoilPropertiesFromDeck_impl.hpp"
 
 #endif // OPM_BLACKOILPROPERTIESFROMDECK_HEADER_INCLUDED

opm/core/props/BlackoilPropertiesFromDeck_impl.hpp

@@ -0,0 +1,260 @@
+namespace Opm
+{
+
+
+    template<class T>
+    BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
+                                                           int number_of_cells,
+                                                           const int* global_cell,
+                                                           const int* cart_dims,
+                                                           T begin_cell_centroids,
+                                                           int dimension,
+                                                           bool init_rock)
+    {
+        init(deck, number_of_cells, global_cell, cart_dims, begin_cell_centroids, dimension,
+             init_rock);
+    }
+
+
+    template<class T>
+    BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
+                                                           int number_of_cells,
+                                                           const int* global_cell,
+                                                           const int* cart_dims,
+                                                           T begin_cell_centroids,
+                                                           int dimension,
+                                                           const parameter::ParameterGroup& param,
+                                                           bool init_rock)
+    {
+        init(deck, number_of_cells, global_cell, cart_dims, begin_cell_centroids, dimension,
+             param, init_rock);
+    }
+
+    template<class T>
+    BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(Opm::DeckConstPtr  newParserDeck,
+                                                           int number_of_cells,
+                                                           const int* global_cell,
+                                                           const int* cart_dims,
+                                                           T begin_cell_centroids,
+                                                           int dimension,
+                                                           bool init_rock)
+    {
+        init(newParserDeck, number_of_cells, global_cell, cart_dims, begin_cell_centroids, dimension,
+             init_rock);
+    }
+
+    template<class T>
+    BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(Opm::DeckConstPtr  newParserDeck,
+                                                           int number_of_cells,
+                                                           const int* global_cell,
+                                                           const int* cart_dims,
+                                                           T begin_cell_centroids,
+                                                           int dimension,
+                                                           const parameter::ParameterGroup& param,
+                                                           bool init_rock)
+    {
+        init(newParserDeck, number_of_cells, global_cell, cart_dims, begin_cell_centroids, dimension,
+             param, init_rock);
+    }
+
+    template<class T>
+    void BlackoilPropertiesFromDeck::init(const EclipseGridParser& deck,
+                                          int number_of_cells,
+                                          const int* global_cell,
+                                          const int* cart_dims,
+                                          T begin_cell_centroids,
+                                          int dimension,
+                                          bool init_rock)
+    {
+        if (init_rock){
+           rock_.init(deck, number_of_cells, global_cell, cart_dims);
+        }
+        pvt_.init(deck, 0);
+        SaturationPropsFromDeck<SatFuncSimpleNonuniform>* ptr
+            = new SaturationPropsFromDeck<SatFuncSimpleNonuniform>();
+        satprops_.reset(ptr);
+        ptr->init(deck, number_of_cells, global_cell, begin_cell_centroids, dimension, 200);
+
+        if (pvt_.numPhases() != satprops_->numPhases()) {
+            OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data ("
+                  << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ").");
+        }
+    }
+
+    template<class T>
+    void BlackoilPropertiesFromDeck::init(Opm::DeckConstPtr newParserDeck,
+                                          int number_of_cells,
+                                          const int* global_cell,
+                                          const int* cart_dims,
+                                          T begin_cell_centroids,
+                                          int dimension,
+                                          bool init_rock)
+    {
+        if (init_rock){
+           rock_.init(newParserDeck, number_of_cells, global_cell, cart_dims);
+        }
+        pvt_.init(newParserDeck, /*numSamples=*/0);
+        SaturationPropsFromDeck<SatFuncSimpleNonuniform>* ptr
+            = new SaturationPropsFromDeck<SatFuncSimpleNonuniform>();
+        satprops_.reset(ptr);
+        ptr->init(newParserDeck, number_of_cells, global_cell, begin_cell_centroids, dimension,
+                  /*numSamples=*/0);
+
+        if (pvt_.numPhases() != satprops_->numPhases()) {
+            OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data ("
+                  << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ").");
+        }
+    }
+
+    template<class T>
+    void BlackoilPropertiesFromDeck::init(const EclipseGridParser& deck,
+                                          int number_of_cells,
+                                          const int* global_cell,
+                                          const int* cart_dims,
+                                          T begin_cell_centroids,
+                                          int dimension,
+                                          const parameter::ParameterGroup& param,
+                                          bool init_rock)
+    {
+        if(init_rock){
+            rock_.init(deck, number_of_cells, global_cell, cart_dims);
+        }
+
+        const int pvt_samples = param.getDefault("pvt_tab_size", 0);
+        pvt_.init(deck, pvt_samples);
+
+        // Unfortunate lack of pointer smartness here...
+        const int sat_samples = param.getDefault("sat_tab_size", 0);
+        std::string threephase_model = param.getDefault<std::string>("threephase_model", "simple");
+        if (deck.hasField("ENDSCALE") && threephase_model != "simple") {
+            OPM_THROW(std::runtime_error, "Sorry, end point scaling currently available for the 'simple' model only.");
+        }
+        if (sat_samples > 1) {
+            if (threephase_model == "stone2") {
+                SaturationPropsFromDeck<SatFuncStone2Uniform>* ptr
+                    = new SaturationPropsFromDeck<SatFuncStone2Uniform>();
+                satprops_.reset(ptr);
+                ptr->init(deck, number_of_cells, global_cell, begin_cell_centroids, dimension,
+                          sat_samples);
+            } else if (threephase_model == "simple") {
+                SaturationPropsFromDeck<SatFuncSimpleUniform>* ptr
+                    = new SaturationPropsFromDeck<SatFuncSimpleUniform>();
+                satprops_.reset(ptr);
+                ptr->init(deck, number_of_cells, global_cell, begin_cell_centroids, dimension,
+                          sat_samples);
+            } else if (threephase_model == "gwseg") {
+                SaturationPropsFromDeck<SatFuncGwsegUniform>* ptr
+                    = new SaturationPropsFromDeck<SatFuncGwsegUniform>();
+                satprops_.reset(ptr);
+                ptr->init(deck, number_of_cells, global_cell, begin_cell_centroids, dimension,
+                          sat_samples);
+            } else {
+                OPM_THROW(std::runtime_error, "Unknown threephase_model: " << threephase_model);
+            }
+        } else {
+            if (threephase_model == "stone2") {
+                SaturationPropsFromDeck<SatFuncStone2Nonuniform>* ptr
+                    = new SaturationPropsFromDeck<SatFuncStone2Nonuniform>();
+                satprops_.reset(ptr);
+                ptr->init(deck, number_of_cells, global_cell, begin_cell_centroids, dimension,
+                          sat_samples);
+            } else if (threephase_model == "simple") {
+                SaturationPropsFromDeck<SatFuncSimpleNonuniform>* ptr
+                    = new SaturationPropsFromDeck<SatFuncSimpleNonuniform>();
+                satprops_.reset(ptr);
+                ptr->init(deck, number_of_cells, global_cell, begin_cell_centroids, dimension,
+                          sat_samples);
+            } else if (threephase_model == "gwseg") {
+                SaturationPropsFromDeck<SatFuncGwsegNonuniform>* ptr
+                    = new SaturationPropsFromDeck<SatFuncGwsegNonuniform>();
+                satprops_.reset(ptr);
+                ptr->init(deck, number_of_cells, global_cell, begin_cell_centroids, dimension,
+                          sat_samples);
+            } else {
+                OPM_THROW(std::runtime_error, "Unknown threephase_model: " << threephase_model);
+            }
+        }
+
+        if (pvt_.numPhases() != satprops_->numPhases()) {
+            OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data ("
+                  << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ").");
+        }
+    }
+
+    template<class T>
+    void BlackoilPropertiesFromDeck::init(Opm::DeckConstPtr newParserDeck,
+                                          int number_of_cells,
+                                          const int* global_cell,
+                                          const int* cart_dims,
+                                          T begin_cell_centroids,
+                                          int dimension,
+                                          const parameter::ParameterGroup& param,
+                                          bool init_rock)
+    {
+
+        if(init_rock){
+            rock_.init(newParserDeck, number_of_cells, global_cell, cart_dims);
+        }
+
+        const int pvt_samples = param.getDefault("pvt_tab_size", 200);
+        pvt_.init(newParserDeck, pvt_samples);
+
+        // Unfortunate lack of pointer smartness here...
+        const int sat_samples = param.getDefault("sat_tab_size", 200);
+        std::string threephase_model = param.getDefault<std::string>("threephase_model", "simple");
+        if (newParserDeck->hasKeyword("ENDSCALE") && threephase_model != "simple") {
+            OPM_THROW(std::runtime_error, "Sorry, end point scaling currently available for the 'simple' model only.");
+        }
+        if (sat_samples > 1) {
+            if (threephase_model == "stone2") {
+                SaturationPropsFromDeck<SatFuncStone2Uniform>* ptr
+                    = new SaturationPropsFromDeck<SatFuncStone2Uniform>();
+                satprops_.reset(ptr);
+                ptr->init(newParserDeck, number_of_cells, global_cell, begin_cell_centroids,
+                          dimension, sat_samples);
+            } else if (threephase_model == "simple") {
+                SaturationPropsFromDeck<SatFuncSimpleUniform>* ptr
+                    = new SaturationPropsFromDeck<SatFuncSimpleUniform>();
+                satprops_.reset(ptr);
+                ptr->init(newParserDeck, number_of_cells, global_cell, begin_cell_centroids,
+                          dimension, sat_samples);
+            } else if (threephase_model == "gwseg") {
+                SaturationPropsFromDeck<SatFuncGwsegUniform>* ptr
+                    = new SaturationPropsFromDeck<SatFuncGwsegUniform>();
+                satprops_.reset(ptr);
+                ptr->init(newParserDeck, number_of_cells, global_cell, begin_cell_centroids,
+                          dimension, sat_samples);
+            } else {
+                OPM_THROW(std::runtime_error, "Unknown threephase_model: " << threephase_model);
+            }
+        } else {
+            if (threephase_model == "stone2") {
+                SaturationPropsFromDeck<SatFuncStone2Nonuniform>* ptr
+                    = new SaturationPropsFromDeck<SatFuncStone2Nonuniform>();
+                satprops_.reset(ptr);
+                ptr->init(newParserDeck, number_of_cells, global_cell, begin_cell_centroids,
+                          dimension, sat_samples);
+            } else if (threephase_model == "simple") {
+                SaturationPropsFromDeck<SatFuncSimpleNonuniform>* ptr
+                    = new SaturationPropsFromDeck<SatFuncSimpleNonuniform>();
+                satprops_.reset(ptr);
+                ptr->init(newParserDeck, number_of_cells, global_cell, begin_cell_centroids,
+                          dimension, sat_samples);
+            } else if (threephase_model == "gwseg") {
+                SaturationPropsFromDeck<SatFuncGwsegNonuniform>* ptr
+                    = new SaturationPropsFromDeck<SatFuncGwsegNonuniform>();
+                satprops_.reset(ptr);
+                ptr->init(newParserDeck, number_of_cells, global_cell, begin_cell_centroids,
+                          dimension, sat_samples);
+            } else {
+                OPM_THROW(std::runtime_error, "Unknown threephase_model: " << threephase_model);
+            }
+        }
+
+        if (pvt_.numPhases() != satprops_->numPhases()) {
+            OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data ("
+                  << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ").");
+        }
+    }
+}
+

opm/core/props/rock/RockFromDeck.cpp

@@ -65,45 +65,58 @@ namespace Opm
     void RockFromDeck::init(const EclipseGridParser& deck,
                             const UnstructuredGrid& grid)
     {
-        assignPorosity(deck, grid);
-        permfield_valid_.assign(grid.number_of_cells, false);
+        init(deck, grid.number_of_cells, grid.global_cell, grid.cartdims);
+    }
+
+    /// Initialize from deck and cell mapping.
+    /// \param  deck            Deck input parser
+    /// \param  number_of_cells The number of cells in the grid.
+    /// \param  global_cell     The mapping fom local to global cell indices.
+    ///                         global_cell[i] is the corresponding global index of i.
+    /// \param  cart_dims       The size of the underlying cartesian grid.
+    void RockFromDeck::init(const EclipseGridParser& deck,
+                            int number_of_cells, const int* global_cell,
+                            const int* cart_dims)
+    {
+        assignPorosity(deck, number_of_cells, global_cell);
+        permfield_valid_.assign(number_of_cells, false);
         const double perm_threshold = 0.0; // Maybe turn into parameter?
-        assignPermeability(deck, grid, perm_threshold);
+        assignPermeability(deck, number_of_cells, global_cell, cart_dims, perm_threshold);
     }
 
     void RockFromDeck::init(Opm::DeckConstPtr newParserDeck,
-                            const UnstructuredGrid& grid)
+                            int number_of_cells, const int* global_cell,
+                            const int* cart_dims)
     {
-        assignPorosity(newParserDeck, grid);
-        permfield_valid_.assign(grid.number_of_cells, false);
+        assignPorosity(newParserDeck, number_of_cells, global_cell);
+        permfield_valid_.assign(number_of_cells, false);
         const double perm_threshold = 0.0; // Maybe turn into parameter?
-        assignPermeability(newParserDeck, grid, perm_threshold);
+        assignPermeability(newParserDeck, number_of_cells, global_cell, cart_dims,
+                           perm_threshold);
     }
 
 
     void RockFromDeck::assignPorosity(const EclipseGridParser& parser,
-                                      const UnstructuredGrid& grid)
+                                      int number_of_cells, const int* global_cell)
     {
-        porosity_.assign(grid.number_of_cells, 1.0);
-        const int* gc = grid.global_cell;
+        porosity_.assign(number_of_cells, 1.0);
         if (parser.hasField("PORO")) {
             const std::vector<double>& poro = parser.getFloatingPointValue("PORO");
             for (int c = 0; c < int(porosity_.size()); ++c) {
-                const int deck_pos = (gc == NULL) ? c : gc[c];
+                const int deck_pos = (global_cell == NULL) ? c : global_cell[c];
                 porosity_[c] = poro[deck_pos];
             }
         }
     }
 
     void RockFromDeck::assignPorosity(Opm::DeckConstPtr newParserDeck,
-                                      const UnstructuredGrid& grid)
+                                      int number_of_cells, const int* global_cell)
     {
-        porosity_.assign(grid.number_of_cells, 1.0);
-        const int* gc = grid.global_cell;
+        porosity_.assign(number_of_cells, 1.0);
         if (newParserDeck->hasKeyword("PORO")) {
             const std::vector<double>& poro = newParserDeck->getKeyword("PORO")->getSIDoubleData();
             for (int c = 0; c < int(porosity_.size()); ++c) {
-                const int deck_pos = (gc == NULL) ? c : gc[c];
+                const int deck_pos = (global_cell == NULL) ? c : global_cell[c];
                 assert(0 <= c && c < (int) porosity_.size());
                 assert(0 <= deck_pos && deck_pos < (int) poro.size());
                 porosity_[c] = poro[deck_pos];
@@ -113,16 +126,17 @@ namespace Opm
 
 
     void RockFromDeck::assignPermeability(const EclipseGridParser& parser,
-                                          const UnstructuredGrid& grid,
+                                          int number_of_cells,
+                                          const int* global_cell,
+                                          const int* cartdims,
                                           double perm_threshold)
     {
         const int dim              = 3;
-        const int num_global_cells = grid.cartdims[0]*grid.cartdims[1]*grid.cartdims[2];
-        const int nc = grid.number_of_cells;
+        const int num_global_cells = cartdims[0]*cartdims[1]*cartdims[2];
 
         assert(num_global_cells > 0);
 
-        permeability_.assign(dim * dim * nc, 0.0);
+        permeability_.assign(dim * dim * number_of_cells, 0.0);
 
         std::vector<const std::vector<double>*> tensor;
         tensor.reserve(10);
@@ -144,13 +158,12 @@ namespace Opm
         // chosen) default value...
         //
         if (tensor.size() > 1) {
-            const int* gc = grid.global_cell;
             int off = 0;
 
-            for (int c = 0; c < nc; ++c, off += dim*dim) {
+            for (int c = 0; c < number_of_cells; ++c, off += dim*dim) {
                 // SharedPermTensor K(dim, dim, &permeability_[off]);
                 int       kix  = 0;
-                const int glob = (gc == NULL) ? c : gc[c];
+                const int glob = (global_cell == NULL) ? c : global_cell[c];
 
                 for (int i = 0; i < dim; ++i) {
                     for (int j = 0; j < dim; ++j, ++kix) {
@@ -167,12 +180,14 @@ namespace Opm
     }
 
     void RockFromDeck::assignPermeability(Opm::DeckConstPtr newParserDeck,
-                                          const UnstructuredGrid& grid,
+                                          int number_of_cells,
+                                          const int* global_cell,
+                                          const int* cartdims,
                                           double perm_threshold)
     {
         const int dim              = 3;
-        const int num_global_cells = grid.cartdims[0]*grid.cartdims[1]*grid.cartdims[2];
-        const int nc = grid.number_of_cells;
+        const int num_global_cells = cartdims[0]*cartdims[1]*cartdims[2];
+        const int nc = number_of_cells;
 
         assert(num_global_cells > 0);
 
@@ -198,7 +213,7 @@ namespace Opm
         // chosen) default value...
         //
         if (tensor.size() > 1) {
-            const int* gc = grid.global_cell;
+            const int* gc = global_cell;
             int off = 0;
 
             for (int c = 0; c < nc; ++c, off += dim*dim) {

opm/core/props/rock/RockFromDeck.hpp

@@ -46,13 +46,24 @@ namespace Opm
         void init(const EclipseGridParser& deck,
                   const UnstructuredGrid& grid);
 
-        /// Initialize from deck and grid.
+        /// Initialize from deck and cell mapping.
+        /// \param  deck            Deck input parser
+        /// \param  number_of_cells The number of cells in the grid.
+        /// \param  global_cell     The mapping fom local to global cell indices.
+        ///                         global_cell[i] is the corresponding global index of i.
+        /// \param  cart_dims       The size of the underlying cartesian grid.
+        void init(const EclipseGridParser& deck,
+                  int number_of_cells, const int* global_cell,
+                  const int* cart_dims);
+         /// Initialize from deck and cell mapping.
         /// \param  newParserDeck            Deck produced by the opm-parser code
-        /// \param  grid          Grid to which property object applies, needed for the
-        ///                       mapping from cell indices (typically from a processed grid)
-        ///                       to logical cartesian indices consistent with the deck.
+        /// \param  number_of_cells The number of cells in the grid.
+        /// \param  global_cell     The mapping fom local to global cell indices.
+        ///                         global_cell[i] is the corresponding global index of i.
+        /// \param  cart_dims       The size of the underlying cartesian grid.
         void init(Opm::DeckConstPtr newParserDeck,
-                  const UnstructuredGrid& grid);
+                  int number_of_cells, const int* global_cell,
+                  const int* cart_dims);
 
         /// \return   D, the number of spatial dimensions. Always 3 for deck input.
         int numDimensions() const
@@ -82,14 +93,20 @@ namespace Opm
 
     private:
         void assignPorosity(const EclipseGridParser& parser,
-                            const UnstructuredGrid& grid);
+                            int number_of_cells,
+                            const int* global_cell);
         void assignPorosity(Opm::DeckConstPtr newParserDeck,
-                            const UnstructuredGrid& grid);
+                            int number_of_cells,
+                            const int* global_cell);
         void assignPermeability(const EclipseGridParser& parser,
-                                const UnstructuredGrid& grid,
+                                int number_of_cells,
+                                const int* global_cell,
+                                const int* cart_dims,
                                 const double perm_threshold);
         void assignPermeability(Opm::DeckConstPtr newParserDeck,
-                                const UnstructuredGrid& grid,
+                                int number_of_cells,
+                                const int* global_cell,
+                                const int* cart_dims,
                                 double perm_threshold);
 
         std::vector<double> porosity_;

opm/core/props/satfunc/SaturationPropsFromDeck.hpp

@@ -65,6 +65,28 @@ namespace Opm
 
         /// Initialize from deck and grid.
         /// \param[in]  deck            Deck input parser
+        /// \param[in]  number_of_cells The number of cells of the grid to which property
+        ///                             object applies, needed for the
+        ///                             mapping from cell indices (typically from a processed
+        ///                             grid) to logical cartesian indices consistent with the
+        ///                             deck.
+        /// \param[in]  global_cell     The mapping from local cell indices of the grid to
+        ///                             global cell indices used in the deck.
+        /// \param[in]  begin_cell_centroids Pointer to the first cell_centroid of the grid.
+        /// \param[in]  dimensions      The dimensions of the grid. 
+        /// \param[in]  samples         Number of uniform sample points for saturation tables.
+        /// \tparam     T               The iterator Type for the cell centroids.
+        /// NOTE: samples will only be used with the SatFuncSetUniform template argument.
+        template<class T>
+        void init(const EclipseGridParser& deck,
+                  int number_of_cells,
+                  const int* global_cell,
+                  const T& begin_cell_centroids,
+                  int dimensions,
+                  const int samples);
+        
+        /// Initialize from deck and grid.
+        /// \param[in]  newParserDeck     Deck input parser
         /// \param[in]  grid     Grid to which property object applies, needed for the
         ///                      mapping from cell indices (typically from a processed grid)
         ///                      to logical cartesian indices consistent with the deck.
@@ -74,6 +96,27 @@ namespace Opm
                   const UnstructuredGrid& grid,
                   const int samples);
 
+        /// Initialize from deck and grid.
+        /// \param[in]  newParserDeck     Deck input parser
+        /// \param[in]  number_of_cells The number of cells of the grid to which property
+        ///                             object applies, needed for the
+        ///                             mapping from cell indices (typically from a processed
+        ///                             grid) to logical cartesian indices consistent with the
+        ///                             deck.
+        /// \param[in]  global_cell     The mapping from local cell indices of the grid to
+        ///                             global cell indices used in the deck.
+        /// \param[in]  begin_cell_centroids Pointer to the first cell_centroid of the grid.
+        /// \param[in]  dimensions      The dimensions of the grid. 
+        /// \param[in]  samples  Number of uniform sample points for saturation tables.
+        /// NOTE: samples will only be used with the SatFuncSetUniform template argument.
+        template<class T>
+        void init(Opm::DeckConstPtr newParserDeck,
+                  int number_of_cells,
+                  const int* global_cell,
+                  const T& begin_cell_centroids,
+                  int dimensions,
+                  const int samples);
+
         /// \return   P, the number of phases.
         int numPhases() const;
 
@@ -138,20 +181,45 @@ namespace Opm
         typedef SatFuncSet Funcs;
 
         const Funcs& funcForCell(const int cell) const;
+
+        template<class T>
         void initEPS(const EclipseGridParser& deck,
-                     const UnstructuredGrid& grid);
+                     int number_of_cells,
+                     const int* global_cell,
+                     const T& begin_cell_centroids,
+                     int dimensions);
+        template<class T>
         void initEPSHyst(const EclipseGridParser& deck,
-                         const UnstructuredGrid& grid);
+                         int number_of_cells,
+                         const int* global_cell,
+                         const T& begin_cell_centroids,
+                         int dimensions);
+        template<class T>
         void initEPSKey(const EclipseGridParser& deck,
-                        const UnstructuredGrid& grid,
+                        int number_of_cells,
+                        const int* global_cell,
+                        const T& begin_cell_centroids,
+                        int dimensions,
                         const std::string& keyword,
                         std::vector<double>& scaleparam);
+        template<class T>
         void initEPS(Opm::DeckConstPtr newParserDeck,
-                     const UnstructuredGrid& grid);
+                     int number_of_cells,
+                     const int* global_cell,
+                     const T& begin_cell_centroids,
+                     int dimensions);
+        template<class T>
         void initEPSHyst(Opm::DeckConstPtr newParserDeck,
-                         const UnstructuredGrid& grid);
+                         int number_of_cells,
+                         const int* global_cell,
+                         const T& begin_cell_centroids,
+                         int dimensions);
+        template<class T>
         void initEPSKey(Opm::DeckConstPtr newParserDeck,
-                        const UnstructuredGrid& grid,
+                        int number_of_cells,
+                        const int* global_cell,
+                        const T& begin_cell_centroids,
+                        int dimensions,
                         const std::string& keyword,
                         std::vector<double>& scaleparam);
         void initEPSParam(const int cell, 

opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp

@@ -25,6 +25,7 @@
 #include <opm/core/utility/NonuniformTableLinear.hpp>
 #include <opm/core/props/phaseUsageFromDeck.hpp>
 #include <opm/core/grid.h>
+#include <opm/core/grid/GridHelpers.hpp>
 
 #include <opm/parser/eclipse/Utility/EndscaleWrapper.hpp>
 #include <opm/parser/eclipse/Utility/ScalecrsWrapper.hpp>
@@ -51,6 +52,20 @@ namespace Opm
     void SaturationPropsFromDeck<SatFuncSet>::init(const EclipseGridParser& deck,
                                                    const UnstructuredGrid& grid,
                                                    const int samples)
+    {
+        init(deck, grid.number_of_cells, grid.global_cell, grid.cell_centroids,
+             grid.dimensions, samples);
+    }
+
+    /// Initialize from deck.
+    template <class SatFuncSet>
+    template< class T>
+    void SaturationPropsFromDeck<SatFuncSet>::init(const EclipseGridParser& deck,
+                                                   int number_of_cells,
+                                                   const int* global_cell,
+                                                   const T& begin_cell_centroids,
+                                                   int dimensions,
+                                                   const int samples)
     {
         phase_usage_ = phaseUsageFromDeck(deck);
 
@@ -66,11 +81,9 @@ namespace Opm
         if (deck.hasField("SATNUM")) {
             const std::vector<int>& satnum = deck.getIntegerValue("SATNUM");
             satfuncs_expected = *std::max_element(satnum.begin(), satnum.end());
-            const int num_cells = grid.number_of_cells;
-            cell_to_func_.resize(num_cells);
-            const int* gc = grid.global_cell;
-            for (int cell = 0; cell < num_cells; ++cell) {
-                const int deck_pos = (gc == NULL) ? cell : gc[cell];
+            cell_to_func_.resize(number_of_cells);
+            for (int cell = 0; cell < number_of_cells; ++cell) {
+                const int deck_pos = (global_cell == NULL) ? cell : global_cell[cell];
                 cell_to_func_[cell] = satnum[deck_pos] - 1;
             }
         }
@@ -124,8 +137,7 @@ namespace Opm
                 }
             }
             do_eps_ = true;
-            
-            initEPS(deck, grid);
+            initEPS(deck, number_of_cells, global_cell, begin_cell_centroids, dimensions);
             
             // For now, a primitive detection of hysteresis. TODO: SATOPTS HYSTER/ and EHYSTR
             do_hyst_ = deck.hasField("ISWL") || deck.hasField("ISWU") || deck.hasField("ISWCR") || deck.hasField("ISGL") ||
@@ -137,7 +149,8 @@ namespace Opm
                     deck.hasField("IKRGR") || deck.hasField("IKRORW") || deck.hasField("IKRORG") ) {
                     OPM_THROW(std::runtime_error, "SaturationPropsFromDeck::init()   --  ENDSCALE: Currently hysteresis and relperm value scaling can not be combined.");
                 }
-                initEPSHyst(deck, grid);
+                initEPSHyst(deck, number_of_cells, global_cell, begin_cell_centroids,
+                            dimensions);
             }
 
             //OPM_THROW(std::runtime_error, "SaturationPropsFromDeck::init()   --  ENDSCALE: Under construction ...");
@@ -150,6 +163,21 @@ namespace Opm
     void SaturationPropsFromDeck<SatFuncSet>::init(Opm::DeckConstPtr newParserDeck,
                                                    const UnstructuredGrid& grid,
                                                    const int samples)
+    {
+        this->template init<SatFuncSet>(newParserDeck, grid.number_of_cells, 
+                                        grid.global_cell, grid.cell_centroids,
+                                        grid.dimensions, samples);
+    }
+
+    /// Initialize from deck.
+    template <class SatFuncSet>
+    template<class T>
+    void SaturationPropsFromDeck<SatFuncSet>::init(Opm::DeckConstPtr newParserDeck,
+                                                   int number_of_cells,
+                                                   const int* global_cell,
+                                                   const T& begin_cell_centroids,
+                                                   int dimensions,
+                                                   const int samples)
     {
         phase_usage_ = phaseUsageFromDeck(newParserDeck);
 
@@ -165,9 +193,9 @@ namespace Opm
         if (newParserDeck->hasKeyword("SATNUM")) {
             const std::vector<int>& satnum = newParserDeck->getKeyword("SATNUM")->getIntData();
             satfuncs_expected = *std::max_element(satnum.begin(), satnum.end());
-            const int num_cells = grid.number_of_cells;
+            const int num_cells = number_of_cells;
             cell_to_func_.resize(num_cells);
-            const int* gc = grid.global_cell;
+            const int* gc = global_cell;
             for (int cell = 0; cell < num_cells; ++cell) {
                 const int newParserDeck_pos = (gc == NULL) ? cell : gc[cell];
                 cell_to_func_[cell] = satnum[newParserDeck_pos] - 1;
@@ -225,7 +253,8 @@ namespace Opm
             }
             do_eps_ = true;
 
-            initEPS(newParserDeck, grid);
+            initEPS(newParserDeck, number_of_cells, global_cell, begin_cell_centroids,
+                    dimensions);
 
             // For now, a primitive detection of hysteresis. TODO: SATOPTS HYSTER/ and EHYSTR
             do_hyst_ =
@@ -257,7 +286,8 @@ namespace Opm
                               "Currently hysteresis and relperm value scaling "
                               "cannot be combined.");
                 }
-                initEPSHyst(newParserDeck, grid);
+                initEPSHyst(newParserDeck, number_of_cells, global_cell, begin_cell_centroids,
+                            dimensions);
             }
         }
     }
@@ -443,29 +473,48 @@ namespace Opm
 
     // Initialize saturation scaling parameters
     template <class SatFuncSet>
+    template <class T>
     void SaturationPropsFromDeck<SatFuncSet>::initEPS(const EclipseGridParser& deck,
-                                                      const UnstructuredGrid& grid)
+                                                      int number_of_cells,
+                                                      const int* global_cell,
+                                                      const T& begin_cell_centroid,
+                                                      int dimensions)
     {
         std::vector<double> swl, swcr, swu, sgl, sgcr, sgu, sowcr, sogcr;
         std::vector<double> krw, krg, kro, krwr, krgr, krorw, krorg;
         // Initialize saturation scaling parameter
-        initEPSKey(deck, grid, std::string("SWL"),   swl);
-        initEPSKey(deck, grid, std::string("SWU"),   swu);
-        initEPSKey(deck, grid, std::string("SWCR"),  swcr);
-        initEPSKey(deck, grid, std::string("SGL"),   sgl);
-        initEPSKey(deck, grid, std::string("SGU"),   sgu);
-        initEPSKey(deck, grid, std::string("SGCR"),  sgcr);
-        initEPSKey(deck, grid, std::string("SOWCR"), sowcr);
-        initEPSKey(deck, grid, std::string("SOGCR"), sogcr);
-        initEPSKey(deck, grid, std::string("KRW"),   krw);
-        initEPSKey(deck, grid, std::string("KRG"),   krg);
-        initEPSKey(deck, grid, std::string("KRO"),   kro);
-        initEPSKey(deck, grid, std::string("KRWR"),  krwr);
-        initEPSKey(deck, grid, std::string("KRGR"),  krgr);
-        initEPSKey(deck, grid, std::string("KRORW"), krorw);
-        initEPSKey(deck, grid, std::string("KRORG"), krorg);
+        initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("SWL"),   swl);
+        initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("SWU"),   swu);
+        initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("SWCR"),  swcr);
+        initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("SGL"),   sgl);
+        initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("SGU"),   sgu);
+        initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("SGCR"),  sgcr);
+        initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("SOWCR"), sowcr);
+        initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("SOGCR"), sogcr);
+        initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("KRW"),   krw);
+        initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("KRG"),   krg);
+        initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("KRO"),   kro);
+        initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("KRWR"),  krwr);
+        initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("KRGR"),  krgr);
+        initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("KRORW"), krorw);
+        initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("KRORG"), krorg);
 
-        eps_transf_.resize(grid.number_of_cells);
+        eps_transf_.resize(number_of_cells);
 
         const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua];
         const int gpos = phase_usage_.phase_pos[BlackoilPhases::Vapour];
@@ -473,7 +522,7 @@ namespace Opm
         const bool oilGas = !phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour];
         const bool threephase = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour];
 
-        for (int cell = 0; cell < grid.number_of_cells; ++cell) {
+        for (int cell = 0; cell < number_of_cells; ++cell) {
             if (oilWater) {
                 // ### krw
                 initEPSParam(cell, eps_transf_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos],
@@ -507,29 +556,48 @@ namespace Opm
 
     // Initialize saturation scaling parameters
     template <class SatFuncSet>
+    template<class T>
     void SaturationPropsFromDeck<SatFuncSet>::initEPS(Opm::DeckConstPtr newParserDeck,
-                                                      const UnstructuredGrid& grid)
+                                                      int number_of_cells,
+                                                      const int* global_cell,
+                                                      const T& begin_cell_centroid,
+                                                      int dimensions)
     {
         std::vector<double> swl, swcr, swu, sgl, sgcr, sgu, sowcr, sogcr;
         std::vector<double> krw, krg, kro, krwr, krgr, krorw, krorg;
         // Initialize saturation scaling parameter
-        initEPSKey(newParserDeck, grid, std::string("SWL"),   swl);
-        initEPSKey(newParserDeck, grid, std::string("SWU"),   swu);
-        initEPSKey(newParserDeck, grid, std::string("SWCR"),  swcr);
-        initEPSKey(newParserDeck, grid, std::string("SGL"),   sgl);
-        initEPSKey(newParserDeck, grid, std::string("SGU"),   sgu);
-        initEPSKey(newParserDeck, grid, std::string("SGCR"),  sgcr);
-        initEPSKey(newParserDeck, grid, std::string("SOWCR"), sowcr);
-        initEPSKey(newParserDeck, grid, std::string("SOGCR"), sogcr);
-        initEPSKey(newParserDeck, grid, std::string("KRW"),   krw);
-        initEPSKey(newParserDeck, grid, std::string("KRG"),   krg);
-        initEPSKey(newParserDeck, grid, std::string("KRO"),   kro);
-        initEPSKey(newParserDeck, grid, std::string("KRWR"),  krwr);
-        initEPSKey(newParserDeck, grid, std::string("KRGR"),  krgr);
-        initEPSKey(newParserDeck, grid, std::string("KRORW"), krorw);
-        initEPSKey(newParserDeck, grid, std::string("KRORG"), krorg);
+        initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("SWL"),   swl);
+        initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("SWU"),   swu);
+        initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("SWCR"),  swcr);
+        initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("SGL"),   sgl);
+        initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("SGU"),   sgu);
+        initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("SGCR"),  sgcr);
+        initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("SOWCR"), sowcr);
+        initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("SOGCR"), sogcr);
+        initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("KRW"),   krw);
+        initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("KRG"),   krg);
+        initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("KRO"),   kro);
+        initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("KRWR"),  krwr);
+        initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("KRGR"),  krgr);
+        initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("KRORW"), krorw);
+        initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("KRORG"), krorg);
 
-        eps_transf_.resize(grid.number_of_cells);
+        eps_transf_.resize(number_of_cells);
 
         const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua];
         const int gpos = phase_usage_.phase_pos[BlackoilPhases::Vapour];
@@ -537,7 +605,7 @@ namespace Opm
         const bool oilGas = !phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour];
         const bool threephase = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour];
 
-        for (int cell = 0; cell < grid.number_of_cells; ++cell) {
+        for (int cell = 0; cell < number_of_cells; ++cell) {
             if (oilWater) {
                 // ### krw
                 initEPSParam(cell, eps_transf_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos],
@@ -571,31 +639,50 @@ namespace Opm
 
     // Initialize hysteresis saturation scaling parameters
     template <class SatFuncSet>
+    template<class T>
     void SaturationPropsFromDeck<SatFuncSet>::initEPSHyst(const EclipseGridParser& deck,
-                                                          const UnstructuredGrid& grid)
+                                                      int number_of_cells,
+                                                      const int* global_cell,
+                                                      const T& begin_cell_centroid,
+                                                      int dimensions)
     {
         std::vector<double> iswl, iswcr, iswu, isgl, isgcr, isgu, isowcr, isogcr;
         std::vector<double> ikrw, ikrg, ikro, ikrwr, ikrgr, ikrorw, ikrorg;
         // Initialize hysteresis saturation scaling parameters
-        initEPSKey(deck, grid, std::string("ISWL"),   iswl);
-        initEPSKey(deck, grid, std::string("ISWU"),   iswu);
-        initEPSKey(deck, grid, std::string("ISWCR"),  iswcr);
-        initEPSKey(deck, grid, std::string("ISGL"),   isgl);
-        initEPSKey(deck, grid, std::string("ISGU"),   isgu);
-        initEPSKey(deck, grid, std::string("ISGCR"),  isgcr);
-        initEPSKey(deck, grid, std::string("ISOWCR"), isowcr);
-        initEPSKey(deck, grid, std::string("ISOGCR"), isogcr);
-        initEPSKey(deck, grid, std::string("IKRW"),   ikrw);
-        initEPSKey(deck, grid, std::string("IKRG"),   ikrg);
-        initEPSKey(deck, grid, std::string("IKRO"),   ikro);
-        initEPSKey(deck, grid, std::string("IKRWR"),  ikrwr);
-        initEPSKey(deck, grid, std::string("IKRGR"),  ikrgr);
-        initEPSKey(deck, grid, std::string("IKRORW"), ikrorw);
-        initEPSKey(deck, grid, std::string("IKRORG"), ikrorg);
+        initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("ISWL"),   iswl);
+        initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("ISWU"),   iswu);
+        initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("ISWCR"),  iswcr);
+        initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("ISGL"),   isgl);
+        initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("ISGU"),   isgu);
+        initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("ISGCR"),  isgcr);
+        initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("ISOWCR"), isowcr);
+        initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("ISOGCR"), isogcr);
+        initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("IKRW"),   ikrw);
+        initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("IKRG"),   ikrg);
+        initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("IKRO"),   ikro);
+        initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("IKRWR"),  ikrwr);
+        initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("IKRGR"),  ikrgr);
+        initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("IKRORW"), ikrorw);
+        initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("IKRORG"), ikrorg);
         
 
-        eps_transf_hyst_.resize(grid.number_of_cells);
-        sat_hyst_.resize(grid.number_of_cells);
+        eps_transf_hyst_.resize(number_of_cells);
+        sat_hyst_.resize(number_of_cells);
 
         const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua];
         const int gpos = phase_usage_.phase_pos[BlackoilPhases::Vapour];
@@ -603,7 +690,7 @@ namespace Opm
         const bool oilGas = !phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour];
         const bool threephase = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour];
 
-        for (int cell = 0; cell < grid.number_of_cells; ++cell) {
+        for (int cell = 0; cell < number_of_cells; ++cell) {
             if (oilWater) {
                 // ### krw
                 initEPSParam(cell, eps_transf_hyst_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos],
@@ -637,30 +724,49 @@ namespace Opm
 
     // Initialize hysteresis saturation scaling parameters
     template <class SatFuncSet>
+    template<class T>
     void SaturationPropsFromDeck<SatFuncSet>::initEPSHyst(Opm::DeckConstPtr newParserDeck,
-                                                          const UnstructuredGrid& grid)
+                                                      int number_of_cells,
+                                                      const int* global_cell,
+                                                      const T& begin_cell_centroid,
+                                                      int dimensions)
     {
         std::vector<double> iswl, iswcr, iswu, isgl, isgcr, isgu, isowcr, isogcr;
         std::vector<double> ikrw, ikrg, ikro, ikrwr, ikrgr, ikrorw, ikrorg;
         // Initialize hysteresis saturation scaling parameters
-        initEPSKey(newParserDeck, grid, std::string("ISWL"),   iswl);
-        initEPSKey(newParserDeck, grid, std::string("ISWU"),   iswu);
-        initEPSKey(newParserDeck, grid, std::string("ISWCR"),  iswcr);
-        initEPSKey(newParserDeck, grid, std::string("ISGL"),   isgl);
-        initEPSKey(newParserDeck, grid, std::string("ISGU"),   isgu);
-        initEPSKey(newParserDeck, grid, std::string("ISGCR"),  isgcr);
-        initEPSKey(newParserDeck, grid, std::string("ISOWCR"), isowcr);
-        initEPSKey(newParserDeck, grid, std::string("ISOGCR"), isogcr);
-        initEPSKey(newParserDeck, grid, std::string("IKRW"),   ikrw);
-        initEPSKey(newParserDeck, grid, std::string("IKRG"),   ikrg);
-        initEPSKey(newParserDeck, grid, std::string("IKRO"),   ikro);
-        initEPSKey(newParserDeck, grid, std::string("IKRWR"),  ikrwr);
-        initEPSKey(newParserDeck, grid, std::string("IKRGR"),  ikrgr);
-        initEPSKey(newParserDeck, grid, std::string("IKRORW"), ikrorw);
-        initEPSKey(newParserDeck, grid, std::string("IKRORG"), ikrorg);
+        initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("ISWL"),   iswl);
+        initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("ISWU"),   iswu);
+        initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("ISWCR"),  iswcr);
+        initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("ISGL"),   isgl);
+        initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("ISGU"),   isgu);
+        initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("ISGCR"),  isgcr);
+        initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("ISOWCR"), isowcr);
+        initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("ISOGCR"), isogcr);
+        initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("IKRW"),   ikrw);
+        initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("IKRG"),   ikrg);
+        initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("IKRO"),   ikro);
+        initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("IKRWR"),  ikrwr);
+        initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("IKRGR"),  ikrgr);
+        initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("IKRORW"), ikrorw);
+        initEPSKey(newParserDeck, number_of_cells, global_cell, begin_cell_centroid, dimensions,
+                   std::string("IKRORG"), ikrorg);
 
-        eps_transf_hyst_.resize(grid.number_of_cells);
-        sat_hyst_.resize(grid.number_of_cells);
+        eps_transf_hyst_.resize(number_of_cells);
+        sat_hyst_.resize(number_of_cells);
 
         const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua];
         const int gpos = phase_usage_.phase_pos[BlackoilPhases::Vapour];
@@ -668,7 +774,7 @@ namespace Opm
         const bool oilGas = !phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour];
         const bool threephase = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour];
 
-        for (int cell = 0; cell < grid.number_of_cells; ++cell) {
+        for (int cell = 0; cell < number_of_cells; ++cell) {
             if (oilWater) {
                 // ### krw
                 initEPSParam(cell, eps_transf_hyst_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos],
@@ -702,8 +808,12 @@ namespace Opm
 
     // Initialize saturation scaling parameter
     template <class SatFuncSet>
+    template<class T>
     void SaturationPropsFromDeck<SatFuncSet>::initEPSKey(const EclipseGridParser& deck,
-                                                         const UnstructuredGrid& grid,
+                                                         int number_of_cells,
+                                                         const int* global_cell,
+                                                         const T& begin_cell_centroid,
+                                                         int dimensions,
                                                          const std::string& keyword,
                                                          std::vector<double>& scaleparam)
     {
@@ -724,57 +834,57 @@ namespace Opm
             if (keyword == std::string("SWL") || keyword == std::string("ISWL") ) {
                 if (useAqua && (useKeyword || deck.getENPTVD().mask_[0])) {
                     itab = 1;
-                    scaleparam.resize(grid.number_of_cells);
-                    for (int i=0; i<grid.number_of_cells; ++i)
+                    scaleparam.resize(number_of_cells);
+                    for (int i=0; i<number_of_cells; ++i)
                         scaleparam[i] = funcForCell(i).smin_[phase_pos_aqua];
                 }
             } else if (keyword == std::string("SWCR") || keyword == std::string("ISWCR") ) {
                 if (useAqua && (useKeyword || deck.getENPTVD().mask_[1])) {
                     itab = 2;
-                    scaleparam.resize(grid.number_of_cells);
-                    for (int i=0; i<grid.number_of_cells; ++i)
+                    scaleparam.resize(number_of_cells);
+                    for (int i=0; i<number_of_cells; ++i)
                         scaleparam[i] = funcForCell(i).swcr_;
                 }
             } else if (keyword == std::string("SWU") || keyword == std::string("ISWU") ) {
                 if (useAqua && (useKeyword || deck.getENPTVD().mask_[2])) {
                     itab = 3;
-                    scaleparam.resize(grid.number_of_cells);
-                    for (int i=0; i<grid.number_of_cells; ++i)
+                    scaleparam.resize(number_of_cells);
+                    for (int i=0; i<number_of_cells; ++i)
                         scaleparam[i] = funcForCell(i).smax_[phase_pos_aqua];
                 }
             } else if (keyword == std::string("SGL") || keyword == std::string("ISGL") ) {
                 if (useVapour && (useKeyword || deck.getENPTVD().mask_[3])) {
                     itab = 4;
-                    scaleparam.resize(grid.number_of_cells);
-                    for (int i=0; i<grid.number_of_cells; ++i)
+                    scaleparam.resize(number_of_cells);
+                    for (int i=0; i<number_of_cells; ++i)
                         scaleparam[i] = funcForCell(i).smin_[phase_pos_vapour];
                 }
             } else if (keyword == std::string("SGCR") || keyword == std::string("ISGCR") ) {
                 if (useVapour && (useKeyword || deck.getENPTVD().mask_[4])) {
                     itab = 5;
-                    scaleparam.resize(grid.number_of_cells);
-                    for (int i=0; i<grid.number_of_cells; ++i)
+                    scaleparam.resize(number_of_cells);
+                    for (int i=0; i<number_of_cells; ++i)
                         scaleparam[i] = funcForCell(i).sgcr_;
                 }
             } else if (keyword == std::string("SGU") || keyword == std::string("ISGU") ) {
                 if (useVapour && (useKeyword || deck.getENPTVD().mask_[5])) {
                     itab = 6;
-                    scaleparam.resize(grid.number_of_cells);
-                    for (int i=0; i<grid.number_of_cells; ++i)
+                    scaleparam.resize(number_of_cells);
+                    for (int i=0; i<number_of_cells; ++i)
                         scaleparam[i] = funcForCell(i).smax_[phase_pos_vapour];
                 }
             } else if (keyword == std::string("SOWCR") || keyword == std::string("ISOWCR") ) {
                 if (useAqua && (useKeyword || deck.getENPTVD().mask_[6])) {
                     itab = 7;
-                    scaleparam.resize(grid.number_of_cells);
-                    for (int i=0; i<grid.number_of_cells; ++i)
+                    scaleparam.resize(number_of_cells);
+                    for (int i=0; i<number_of_cells; ++i)
                         scaleparam[i] = funcForCell(i).sowcr_;
                 }
             } else if (keyword == std::string("SOGCR") || keyword == std::string("ISOGCR") ) {
                 if (useVapour && (useKeyword || deck.getENPTVD().mask_[7])) {
                     itab = 8;
-                    scaleparam.resize(grid.number_of_cells);
-                    for (int i=0; i<grid.number_of_cells; ++i)
+                    scaleparam.resize(number_of_cells);
+                    for (int i=0; i<number_of_cells; ++i)
                         scaleparam[i] = funcForCell(i).sogcr_;
                 }
             } else {
@@ -787,50 +897,50 @@ namespace Opm
             if (keyword == std::string("KRW") || keyword == std::string("IKRW") ) {
                 if (useAqua && (useKeyword || deck.getENKRVD().mask_[0])) {
                     itab = 1;
-                    scaleparam.resize(grid.number_of_cells);
-                    for (int i=0; i<grid.number_of_cells; ++i)
+                    scaleparam.resize(number_of_cells);
+                    for (int i=0; i<number_of_cells; ++i)
                         scaleparam[i] = funcForCell(i).krwmax_;
                 }
             } else if (keyword == std::string("KRG") || keyword == std::string("IKRG") ) {
                 if (useVapour && (useKeyword || deck.getENKRVD().mask_[1])) {
                     itab = 2;
-                    scaleparam.resize(grid.number_of_cells);
-                    for (int i=0; i<grid.number_of_cells; ++i)
+                    scaleparam.resize(number_of_cells);
+                    for (int i=0; i<number_of_cells; ++i)
                         scaleparam[i] = funcForCell(i).krgmax_;
                 }
             } else if (keyword == std::string("KRO") || keyword == std::string("IKRO") ) {
                 if (useLiquid && (useKeyword || deck.getENKRVD().mask_[2])) {
                     itab = 3;
-                    scaleparam.resize(grid.number_of_cells);
-                    for (int i=0; i<grid.number_of_cells; ++i)
+                    scaleparam.resize(number_of_cells);
+                    for (int i=0; i<number_of_cells; ++i)
                         scaleparam[i] = funcForCell(i).kromax_;
                 }
             } else if (keyword == std::string("KRWR") || keyword == std::string("IKRWR") ) {
                 if (useAqua && (useKeyword || deck.getENKRVD().mask_[3])) {
                     itab = 4;
-                    scaleparam.resize(grid.number_of_cells);
-                    for (int i=0; i<grid.number_of_cells; ++i)
+                    scaleparam.resize(number_of_cells);
+                    for (int i=0; i<number_of_cells; ++i)
                         scaleparam[i] = funcForCell(i).krwr_;
                 }
             } else if (keyword == std::string("KRGR") || keyword == std::string("IKRGR") ) {
                 if (useVapour && (useKeyword || deck.getENKRVD().mask_[4])) {
                     itab = 5;
-                    scaleparam.resize(grid.number_of_cells);
-                    for (int i=0; i<grid.number_of_cells; ++i)
+                    scaleparam.resize(number_of_cells);
+                    for (int i=0; i<number_of_cells; ++i)
                         scaleparam[i] = funcForCell(i).krgr_;
                 }
             } else if (keyword == std::string("KRORW") || keyword == std::string("IKRORW") ) {
                 if (useAqua && (useKeyword || deck.getENKRVD().mask_[5])) {
                     itab = 6;
-                    scaleparam.resize(grid.number_of_cells);
-                    for (int i=0; i<grid.number_of_cells; ++i)
+                    scaleparam.resize(number_of_cells);
+                    for (int i=0; i<number_of_cells; ++i)
                         scaleparam[i] = funcForCell(i).krorw_;
                 }
             } else if (keyword == std::string("KRORG") || keyword == std::string("IKRORG") ) {
                 if (useVapour && (useKeyword || deck.getENKRVD().mask_[6])) {
                     itab = 7;
-                    scaleparam.resize(grid.number_of_cells);
-                    for (int i=0; i<grid.number_of_cells; ++i)
+                    scaleparam.resize(number_of_cells);
+                    for (int i=0; i<number_of_cells; ++i)
                         scaleparam[i] = funcForCell(i).krorg_;
                 }
             } else {
@@ -846,10 +956,9 @@ namespace Opm
         } else if (useKeyword) {
             // Keyword values from deck
             std::cout << "--- Scaling parameter '" << keyword << "' assigned." << std::endl;
-            const int* gc = grid.global_cell;
             const std::vector<double>& val = deck.getFloatingPointValue(keyword);
             for (int c = 0; c < int(scaleparam.size()); ++c) {
-                const int deck_pos = (gc == NULL) ? c : gc[c];
+                const int deck_pos = (global_cell == NULL) ? c : global_cell[c];
                 scaleparam[c] = val[deck_pos];
             }
         } else {
@@ -858,14 +967,15 @@ namespace Opm
                 deck.getENPTVD().write(std::cout);
             else
                 deck.getENKRVD().write(std::cout);
-            const double* cc = grid.cell_centroids;
-            const int dim = grid.dimensions;
-            for (int cell = 0; cell < grid.number_of_cells; ++cell) {
+            for (int cell = 0; cell < number_of_cells; ++cell) {
                 int jtab = cell_to_func_.empty() ? 0 : cell_to_func_[cell];
                 if (table[itab][jtab][0] != -1.0) {
                     std::vector<double>& depth = table[0][jtab];
                     std::vector<double>& val = table[itab][jtab];
-                    double zc = cc[dim*cell+dim-1];
+                    double zc = UgGridHelpers
+                        ::getCoordinate(UgGridHelpers::increment(begin_cell_centroid, cell,
+                                                                 dimensions),
+                                              dimensions-1);
                     if (zc >= depth.front() && zc <= depth.back()) { //don't want extrap outside depth interval
                         scaleparam[cell] = linearInterpolation(depth, val, zc);
                     }
@@ -876,8 +986,12 @@ namespace Opm
 
     // Initialize saturation scaling parameter
     template <class SatFuncSet>
+    template<class T>
     void SaturationPropsFromDeck<SatFuncSet>::initEPSKey(Opm::DeckConstPtr newParserDeck,
-                                                         const UnstructuredGrid& grid,
+                                                         int number_of_cells,
+                                                         const int* global_cell,
+                                                         const T& begin_cell_centroid,
+                                                         int dimensions,
                                                          const std::string& keyword,
                                                          std::vector<double>& scaleparam)
     { 
@@ -899,57 +1013,57 @@ namespace Opm
             if (keyword == std::string("SWL") || keyword == std::string("ISWL") ) {
                 if (useAqua && (useKeyword || columnIsMasked_(newParserDeck, "ENPTVD", 0))) {
                     itab = 1;
-                    scaleparam.resize(grid.number_of_cells);
-                    for (int i=0; i<grid.number_of_cells; ++i)
+                    scaleparam.resize(number_of_cells);
+                    for (int i=0; i<number_of_cells; ++i)
                         scaleparam[i] = funcForCell(i).smin_[phase_pos_aqua];
                 }
             } else if (keyword == std::string("SWCR") || keyword == std::string("ISWCR") ) {
                 if (useAqua && (useKeyword || columnIsMasked_(newParserDeck, "ENPTVD", 1))) {
                     itab = 2;
-                    scaleparam.resize(grid.number_of_cells);
-                    for (int i=0; i<grid.number_of_cells; ++i)
+                    scaleparam.resize(number_of_cells);
+                    for (int i=0; i<number_of_cells; ++i)
                         scaleparam[i] = funcForCell(i).swcr_;
                 }
             } else if (keyword == std::string("SWU") || keyword == std::string("ISWU") ) {
                 if (useAqua && (useKeyword || columnIsMasked_(newParserDeck, "ENPTVD", 2))) {
                     itab = 3;
-                    scaleparam.resize(grid.number_of_cells);
-                    for (int i=0; i<grid.number_of_cells; ++i)
+                    scaleparam.resize(number_of_cells);
+                    for (int i=0; i<number_of_cells; ++i)
                         scaleparam[i] = funcForCell(i).smax_[phase_pos_aqua];
                 }
             } else if (keyword == std::string("SGL") || keyword == std::string("ISGL") ) {
                 if (useVapour && (useKeyword || columnIsMasked_(newParserDeck, "ENPTVD", 3))) {
                     itab = 4;
-                    scaleparam.resize(grid.number_of_cells);
-                    for (int i=0; i<grid.number_of_cells; ++i)
+                    scaleparam.resize(number_of_cells);
+                    for (int i=0; i<number_of_cells; ++i)
                         scaleparam[i] = funcForCell(i).smin_[phase_pos_vapour];
                 }
             } else if (keyword == std::string("SGCR") || keyword == std::string("ISGCR") ) {
                 if (useVapour && (useKeyword || columnIsMasked_(newParserDeck, "ENPTVD", 4))) {
                     itab = 5;
-                    scaleparam.resize(grid.number_of_cells);
-                    for (int i=0; i<grid.number_of_cells; ++i)
+                    scaleparam.resize(number_of_cells);
+                    for (int i=0; i<number_of_cells; ++i)
                         scaleparam[i] = funcForCell(i).sgcr_;
                 }
             } else if (keyword == std::string("SGU") || keyword == std::string("ISGU") ) {
                 if (useVapour && (useKeyword || columnIsMasked_(newParserDeck, "ENPTVD", 5))) {
                     itab = 6;
-                    scaleparam.resize(grid.number_of_cells);
-                    for (int i=0; i<grid.number_of_cells; ++i)
+                    scaleparam.resize(number_of_cells);
+                    for (int i=0; i<number_of_cells; ++i)
                         scaleparam[i] = funcForCell(i).smax_[phase_pos_vapour];
                 }
             } else if (keyword == std::string("SOWCR") || keyword == std::string("ISOWCR") ) {
                 if (useAqua && (useKeyword || columnIsMasked_(newParserDeck, "ENPTVD", 6))) {
                     itab = 7;
-                    scaleparam.resize(grid.number_of_cells);
-                    for (int i=0; i<grid.number_of_cells; ++i)
+                    scaleparam.resize(number_of_cells);
+                    for (int i=0; i<number_of_cells; ++i)
                         scaleparam[i] = funcForCell(i).sowcr_;
                 }
             } else if (keyword == std::string("SOGCR") || keyword == std::string("ISOGCR") ) {
                 if (useVapour && (useKeyword || columnIsMasked_(newParserDeck, "ENPTVD", 7))) {
                     itab = 8;
-                    scaleparam.resize(grid.number_of_cells);
-                    for (int i=0; i<grid.number_of_cells; ++i)
+                    scaleparam.resize(number_of_cells);
+                    for (int i=0; i<number_of_cells; ++i)
                         scaleparam[i] = funcForCell(i).sogcr_;
                 }
             } else {
@@ -970,49 +1084,49 @@ namespace Opm
             if (keyword == std::string("KRW") || keyword == std::string("IKRW") ) {
                 if (useAqua && (useKeyword || columnIsMasked_(newParserDeck, "ENKRVD", 0))) {
                     itab = 1;
-                    scaleparam.resize(grid.number_of_cells);
-                    for (int i=0; i<grid.number_of_cells; ++i)
+                    scaleparam.resize(number_of_cells);
+                    for (int i=0; i<number_of_cells; ++i)
                         scaleparam[i] = funcForCell(i).krwmax_;
                 }
             } else if (keyword == std::string("KRG") || keyword == std::string("IKRG") ) {
                 if (useVapour && (useKeyword || columnIsMasked_(newParserDeck, "ENKRVD", 1))) {
                     itab = 2;
-                    scaleparam.resize(grid.number_of_cells);
-                    for (int i=0; i<grid.number_of_cells; ++i)
+                    scaleparam.resize(number_of_cells);
+                    for (int i=0; i<number_of_cells; ++i)
                         scaleparam[i] = funcForCell(i).krgmax_;
                 }
                 if (useLiquid && (useKeyword || columnIsMasked_(newParserDeck, "ENKRVD", 2))) {
                     itab = 3;
-                    scaleparam.resize(grid.number_of_cells);
-                    for (int i=0; i<grid.number_of_cells; ++i)
+                    scaleparam.resize(number_of_cells);
+                    for (int i=0; i<number_of_cells; ++i)
                         scaleparam[i] = funcForCell(i).kromax_;
                 }
             } else if (keyword == std::string("KRWR") || keyword == std::string("IKRWR") ) {
                 if (useAqua && (useKeyword || columnIsMasked_(newParserDeck, "ENKRVD", 3))) {
                     itab = 4;
-                    scaleparam.resize(grid.number_of_cells);
-                    for (int i=0; i<grid.number_of_cells; ++i)
+                    scaleparam.resize(number_of_cells);
+                    for (int i=0; i<number_of_cells; ++i)
                         scaleparam[i] = funcForCell(i).krwr_;
                 }
             } else if (keyword == std::string("KRGR") || keyword == std::string("IKRGR") ) {
                 if (useVapour && (useKeyword || columnIsMasked_(newParserDeck, "ENKRVD", 4))) {
                     itab = 5;
-                    scaleparam.resize(grid.number_of_cells);
-                    for (int i=0; i<grid.number_of_cells; ++i)
+                    scaleparam.resize(number_of_cells);
+                    for (int i=0; i<number_of_cells; ++i)
                         scaleparam[i] = funcForCell(i).krgr_;
                 }
             } else if (keyword == std::string("KRORW") || keyword == std::string("IKRORW") ) {
                 if (useAqua && (useKeyword || columnIsMasked_(newParserDeck, "ENKRVD", 5))) {
                     itab = 6;
-                    scaleparam.resize(grid.number_of_cells);
-                    for (int i=0; i<grid.number_of_cells; ++i)
+                    scaleparam.resize(number_of_cells);
+                    for (int i=0; i<number_of_cells; ++i)
                         scaleparam[i] = funcForCell(i).krorw_;
                 }
             } else if (keyword == std::string("KRORG") || keyword == std::string("IKRORG") ) {
                 if (useVapour && (useKeyword || columnIsMasked_(newParserDeck, "ENKRVD", 6))) {
                     itab = 7;
-                    scaleparam.resize(grid.number_of_cells);
-                    for (int i=0; i<grid.number_of_cells; ++i)
+                    scaleparam.resize(number_of_cells);
+                    for (int i=0; i<number_of_cells; ++i)
                         scaleparam[i] = funcForCell(i).krorg_;
                 }
             } else {
@@ -1036,19 +1150,20 @@ namespace Opm
         } else if (useKeyword) {
             // Keyword values from deck
             std::cout << "--- Scaling parameter '" << keyword << "' assigned." << std::endl;
-            const int* gc = grid.global_cell;
+            const int* gc = global_cell;
             const std::vector<double>& val = newParserDeck->getKeyword(keyword)->getSIDoubleData();
             for (int c = 0; c < int(scaleparam.size()); ++c) {
                 const int deck_pos = (gc == NULL) ? c : gc[c];
                 scaleparam[c] = val[deck_pos];
             }
         } else {
-            const double* cc = grid.cell_centroids;
-            const int dim = grid.dimensions;
-            for (int cell = 0; cell < grid.number_of_cells; ++cell) {
+            const int dim = dimensions;
+            for (int cell = 0; cell < number_of_cells; ++cell) {
                 int jtab = cell_to_func_.empty() ? 0 : cell_to_func_[cell];
                 if (param_col[jtab][0] >= 0.0) {
-                    double zc = cc[dim*cell+dim-1];
+                    double zc = UgGridHelpers
+                        ::getCoordinate(UgGridHelpers::increment(begin_cell_centroid, cell, dim),
+                                       dim-1);
                     if (zc >= depth_col[jtab].front() && zc <= depth_col[jtab].back()) { //don't want extrap outside depth interval
                         scaleparam[cell] = linearInterpolation(depth_col[jtab], param_col[jtab], zc);
                     }

opm/core/simulator/BlackoilState.cpp

@@ -4,13 +4,18 @@
 using namespace Opm;
 
 void
-BlackoilState::init(const UnstructuredGrid& g, int num_phases) {
-   SimulatorState::init(g, num_phases);
-   gor_.resize(g.number_of_cells, 0.) ;
-   rv_.resize(g.number_of_cells,0.);
+BlackoilState::init(int number_of_cells, int number_of_phases, int num_phases) {
+    SimulatorState::init(number_of_cells, number_of_phases, num_phases);
+   gor_.resize(number_of_cells, 0.) ;
+   rv_.resize(number_of_cells,0.);
    // surfvol_ intentionally empty, left to initBlackoilSurfvol
 }
 
+void
+BlackoilState::init(const UnstructuredGrid& g, int num_phases)
+{
+    init(g.number_of_cells, g.number_of_faces, num_phases);
+}
 /// Set the first saturation to either its min or max value in
 /// the indicated cells. The second saturation value s2 is set
 /// to (1.0 - s1) for each cell. Any further saturation values

opm/core/simulator/BlackoilState.hpp

@@ -32,7 +32,9 @@ namespace Opm
     class BlackoilState : public SimulatorState
     {
     public:
-        virtual void init(const UnstructuredGrid& g, int num_phases);
+        virtual void init(const UnstructuredGrid& grid, int num_phases);
+        
+        virtual void init(int number_of_cells, int number_of_faces, int num_phases);
 
         /// Set the first saturation to either its min or max value in
         /// the indicated cells. The second saturation value s2 is set

opm/core/simulator/SimulatorState.cpp

@@ -40,13 +40,18 @@ SimulatorState::vectorApproxEqual(const std::vector<double>& v1,
 
 void
 SimulatorState::init(const UnstructuredGrid& g, int num_phases)
+{
+    init(g.number_of_cells, g.number_of_faces, num_phases);
+}
+void
+SimulatorState::init(int number_of_cells, int number_of_faces, int num_phases)
 {
     num_phases_ = num_phases;
-    press_.resize(g.number_of_cells, 0.0);
-    fpress_.resize(g.number_of_faces, 0.0);
-    flux_.resize(g.number_of_faces, 0.0);
-    sat_.resize(num_phases_ * g.number_of_cells, 0.0);
-    for (int cell = 0; cell < g.number_of_cells; ++cell) {
+    press_.resize(number_of_cells, 0.0);
+    fpress_.resize(number_of_faces, 0.0);
+    flux_.resize(number_of_faces, 0.0);
+    sat_.resize(num_phases_ * number_of_cells, 0.0);
+    for (int cell = 0; cell < number_of_cells; ++cell) {
         // Defaulting the second saturation to 1.0.
         // This will usually be oil in a water-oil case,
         // gas in an oil-gas case.

opm/core/simulator/SimulatorState.hpp

@@ -17,6 +17,8 @@ namespace Opm
 
         virtual void init(const UnstructuredGrid& g, int num_phases);
 
+        virtual void init(int number_of_cells, int number_of_phases, int num_phases);
+        
         enum ExtremalSat { MinSat, MaxSat };
 
     protected:

opm/core/simulator/initState.hpp

@@ -20,6 +20,8 @@
 #ifndef OPM_INITSTATE_HEADER_INCLUDED
 #define OPM_INITSTATE_HEADER_INCLUDED
 
+#include <opm/parser/eclipse/Deck/Deck.hpp> // DeckConstPtr
+
 struct UnstructuredGrid;
 
 namespace Opm
@@ -65,6 +67,44 @@ namespace Opm
                         const double gravity,
                         State& state);
 
+    /// Initialize a two-phase state from parameters.
+    /// The following parameters are accepted (defaults):
+    ///    - convection_testcase   (false) -- Water in the 'left' part of the grid.
+    ///    - ref_pressure          (100)   -- Initial pressure in bar for all cells
+    ///                                     (if convection_testcase is true),
+    ///                                     or pressure at woc depth.
+    ///    - segregation_testcase  (false) -- Water above the woc instead of below.
+    ///    - water_oil_contact     (none)  -- Depth of water-oil contact (woc).
+    ///    - init_saturation       (none)  -- Initial water saturation for all cells.
+    ///
+    /// If convection_testcase is true, the saturation is initialised
+    /// as indicated, and pressure is initialised to a constant value
+    /// ('ref_pressure').
+    /// If segregation_testcase is true, the saturation is initialised
+    /// as indicated, and pressure is initialised hydrostatically.
+    /// Otherwise we have 3 cases:
+    ///   1. If 'water_oil_contact' is given, saturation is initialised
+    ///      accordingly.
+    ///   2. If 'water_oil_contact' is not given, but 'init_saturation'
+    ///      is given, water saturation is set to that value everywhere.
+    ///   3. If neither are given, water saturation is set to minimum.
+    ///
+    /// In all three cases, pressure is initialised hydrostatically.
+    /// In case 2) and 3), the depth of the first cell is used as reference depth.
+    template <class FaceCells, class CCI, class FCI, class State>
+    void initStateBasic(int number_of_cells,
+                        const int* global_cell,
+                        const int* cartdims,
+                        int number_of_faces,
+                        FaceCells face_cells,
+                        FCI begin_face_centroids,
+                        CCI begin_cell_centroids,
+                        int dimensions,
+                        const IncompPropertiesInterface& props,
+                        const parameter::ParameterGroup& param,
+                        const double gravity,
+                        State& state);
+
     /// Initialize a blackoil state from parameters.
     /// The following parameters are accepted (defaults):
     ///    - convection_testcase   (false) -- Water in the 'left' part of the grid.
@@ -88,6 +128,36 @@ namespace Opm
                         const double gravity,
                         State& state);
 
+    /// Initialize a blackoil state from parameters.
+    /// The following parameters are accepted (defaults):
+    ///    - convection_testcase   (false) -- Water in the 'left' part of the grid.
+    ///    - ref_pressure          (100)   -- Initial pressure in bar for all cells
+    ///                                     (if convection_testcase is true),
+    ///                                     or pressure at woc depth.
+    ///    - water_oil_contact     (none)  -- Depth of water-oil contact (woc).
+    /// If convection_testcase is true, the saturation is initialised
+    /// as indicated, and pressure is initialised to a constant value
+    /// ('ref_pressure').
+    /// Otherwise we have 2 cases:
+    ///   1. If 'water_oil_contact' is given, saturation is initialised
+    ///      accordingly.
+    ///   2. Water saturation is set to minimum.
+    /// In both cases, pressure is initialised hydrostatically.
+    /// In case 2., the depth of the first cell is used as reference depth.
+    template <class FaceCells, class FCI, class CCI, class State>
+    void initStateBasic(int number_of_cells,
+                        const int* global_cell,
+                        const int* cartdims,
+                        int number_of_faces,
+                        FaceCells face_cells,
+                        FCI begin_face_centroids,
+                        CCI begin_cell_centroids,
+                        int dimensions,
+                        const BlackoilPropertiesInterface& props,
+                        const parameter::ParameterGroup& param,
+                        const double gravity,
+                        State& state);
+
     /// Initialize a two-phase state from input deck.
     /// If EQUIL is present:
     ///   - saturation is set according to the water-oil contact,
@@ -102,6 +172,26 @@ namespace Opm
                            const double gravity,
                            State& state);
 
+    /// Initialize a two-phase state from input deck.
+    /// If EQUIL is present:
+    ///   - saturation is set according to the water-oil contact,
+    ///   - pressure is set to hydrostatic equilibrium.
+    /// Otherwise:
+    ///   - saturation is set according to SWAT,
+    ///   - pressure is set according to PRESSURE.
+    template <class FaceCells, class FCI, class CCI, class Props, class State>
+    void initStateFromDeck(int number_of_cells,
+                           const int* global_cell,
+                           int number_of_faces,
+                           FaceCells face_cells,
+                           FCI begin_face_centroids,
+                           CCI begin_cell_centroids,
+                           int dimensions,
+                           const Props& props,
+                           const EclipseGridParser& deck,
+                           const double gravity,
+                           State& state);
+
     /// Initialize a two-phase water-oil blackoil state from input deck.
     /// If EQUIL is present:
     ///   - saturation is set according to the water-oil contact,
@@ -117,7 +207,40 @@ namespace Opm
                                    const double gravity,
                                    State& state);
 
+    /// Initialize a two-phase water-oil blackoil state from input deck.
+    /// If EQUIL is present:
+    ///   - saturation is set according to the water-oil contact,
+    ///   - pressure is set to hydrostatic equilibrium.
+    /// Otherwise:
+    ///   - saturation is set according to SWAT,
+    ///   - pressure is set according to PRESSURE.
+    /// In addition, this function sets surfacevol.
+    template <class FaceCells, class FCI, class CCI, class Props, class State>
+    void initBlackoilStateFromDeck(int number_of_cells,
+                                   const int* global_cell,
+                                   int number_of_faces,
+                                   FaceCells face_cells,
+                                   FCI begin_face_centroids,
+                                   CCI begin_cell_centroids,
+                                   int dimensions,
+                                   const Props& props,
+                                   const EclipseGridParser& deck,
+                                   const double gravity,
+                                   State& state);
 
+    /// Initialize a blackoil state from input deck.
+    template <class FaceCells, class FCI, class CCI, class Props, class State>
+    void initBlackoilStateFromDeck(int number_of_cells,
+                                   const int* global_cell,
+                                   int number_of_faces,
+                                   FaceCells face_cells,
+                                   FCI begin_face_centroids,
+                                   CCI begin_cell_centroids,
+                                   int dimensions,
+                                   const Props& props,
+                                   Opm::DeckConstPtr newParserDeck,
+                                   const double gravity,
+                                   State& state);
 } // namespace Opm
 
 #include <opm/core/simulator/initState_impl.hpp>

opm/core/simulator/initState_impl.hpp

@@ -23,6 +23,7 @@
 
 #include <opm/core/utility/parameters/ParameterGroup.hpp>
 #include <opm/core/grid.h>
+#include <opm/core/grid/GridHelpers.hpp>
 #include <opm/core/io/eclipse/EclipseGridParser.hpp>
 #include <opm/core/utility/ErrorMacros.hpp>
 #include <opm/core/utility/MonotCubicInterpolator.hpp>
@@ -49,17 +50,21 @@ namespace Opm
         // Find the cells that are below and above a depth.
         // TODO: add 'anitialiasing', obtaining a more precise split
         //       by f. ex. subdividing cells cut by the split depths.
-        void cellsBelowAbove(const UnstructuredGrid& grid,
+       template<class T>
+        void cellsBelowAbove(int number_of_cells,
+                             T begin_cell_centroids,
+                             int dimension,
                              const double depth,
                              std::vector<int>& below,
                              std::vector<int>& above)
         {
-            const int num_cells = grid.number_of_cells;
-            below.reserve(num_cells);
-            above.reserve(num_cells);
-            const int dim = grid.dimensions;
-            for (int c = 0; c < num_cells; ++c) {
-                const double z = grid.cell_centroids[dim*c + dim - 1];
+            below.reserve(number_of_cells);
+            above.reserve(number_of_cells);
+            for (int c = 0; c < number_of_cells; ++c) {
+                const double z = UgGridHelpers
+                    ::getCoordinate(UgGridHelpers::increment(begin_cell_centroids, c,
+                                                             dimension),
+                                    dimension-1);
                 if (z > depth) {
                     below.push_back(c);
                 } else {
@@ -76,8 +81,10 @@ namespace Opm
         // Initialize saturations so that there is water below woc,
         // and oil above.
         // If invert is true, water is instead above, oil below.
-        template <class Props, class State>
-        void initWaterOilContact(const UnstructuredGrid& grid,
+        template <class T, class Props, class State>
+        void initWaterOilContact(int number_of_cells,
+                                 T begin_cell_centroids,
+                                 int dimensions,
                                  const Props& props,
                                  const double woc,
                                  const WaterInit waterinit,
@@ -93,10 +100,12 @@ namespace Opm
             // }
             switch (waterinit) {
             case WaterBelow:
-                cellsBelowAbove(grid, woc, water, oil);
+                cellsBelowAbove(number_of_cells, begin_cell_centroids, dimensions,
+                                woc, water, oil);
                 break;
             case WaterAbove:
-                cellsBelowAbove(grid, woc, oil, water);
+                cellsBelowAbove(number_of_cells, begin_cell_centroids, dimensions,
+                                woc, oil, water);
             }
             // Set saturations.
             state.setFirstSat(oil, props, State::MinSat);
@@ -115,8 +124,10 @@ namespace Opm
         // Note that by manipulating the densities parameter,
         // it is possible to initialise with water on top instead of
         // on the bottom etc.
-        template <class State>
-        void initHydrostaticPressure(const UnstructuredGrid& grid,
+        template <class T, class State>
+        void initHydrostaticPressure(int number_of_cells,
+                                     T begin_cell_centroids,
+                                     int dimensions,
                                      const double* densities,
                                      const double woc,
                                      const double gravity,
@@ -125,14 +136,14 @@ namespace Opm
                                      State& state)
         {
             std::vector<double>& p = state.pressure();
-            const int num_cells = grid.number_of_cells;
-            const int dim = grid.dimensions;
             // Compute pressure at woc
             const double rho_datum = datum_z > woc ? densities[0] : densities[1];
             const double woc_p = datum_p + (woc - datum_z)*gravity*rho_datum;
-            for (int c = 0; c < num_cells; ++c) {
+            for (int c = 0; c < number_of_cells; ++c) {
                 // Compute pressure as delta from woc pressure.
-                const double z = grid.cell_centroids[dim*c + dim - 1];
+                const double z = UgGridHelpers::
+                    getCoordinate(UgGridHelpers::increment(begin_cell_centroids, c, dimensions),
+                                  dimensions-1);
                 const double rho = z > woc ? densities[0] : densities[1];
                 p[c] = woc_p + (z - woc)*gravity*rho;
             }
@@ -141,8 +152,10 @@ namespace Opm
 
         // Facade to initHydrostaticPressure() taking a property object,
         // for similarity to initHydrostaticPressure() for compressible fluids.
-        template <class State>
-        void initHydrostaticPressure(const UnstructuredGrid& grid,
+        template <class T, class State>
+        void initHydrostaticPressure(int number_of_cells,
+                                     T begin_cell_centroids,
+                                     int dimensions,
                                      const IncompPropertiesInterface& props,
                                      const double woc,
                                      const double gravity,
@@ -151,7 +164,8 @@ namespace Opm
                                      State& state)
         {
             const double* densities = props.density();
-            initHydrostaticPressure(grid, densities, woc, gravity, datum_z, datum_p, state);
+            initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions,
+                                    densities, woc, gravity, datum_z, datum_p, state);
         }
 
 
@@ -183,8 +197,10 @@ namespace Opm
         // where rho is the oil density above the woc, water density
         // below woc. Note that even if there is (immobile) water in
         // the oil zone, it does not contribute to the pressure there.
-        template <class State>
-        void initHydrostaticPressure(const UnstructuredGrid& grid,
+       template <class T, class State>
+        void initHydrostaticPressure(int number_of_cells,
+                                     T begin_cell_centroids,
+                                     int dimensions,
                                      const BlackoilPropertiesInterface& props,
                                      const double woc,
                                      const double gravity,
@@ -195,11 +211,11 @@ namespace Opm
             assert(props.numPhases() == 2);
 
             // Obtain max and min z for which we will need to compute p.
-            const int num_cells = grid.number_of_cells;
-            const int dim = grid.dimensions;
             double zlim[2] = { 1e100, -1e100 };
-            for (int c = 0; c < num_cells; ++c) {
-                const double z = grid.cell_centroids[dim*c + dim - 1];
+            for (int c = 0; c < number_of_cells; ++c) {
+                const double z = UgGridHelpers::
+                    getCoordinate(UgGridHelpers::increment(begin_cell_centroids, c, dimensions),
+                                  dimensions-1);;
                 zlim[0] = std::min(zlim[0], z);
                 zlim[1] = std::max(zlim[1], z);
             }
@@ -268,30 +284,42 @@ namespace Opm
 
             // Evaluate pressure at each cell centroid.
             std::vector<double>& p = state.pressure();
-            for (int c = 0; c < num_cells; ++c) {
-                const double z = grid.cell_centroids[dim*c + dim - 1];
+            for (int c = 0; c < number_of_cells; ++c) {
+                const double z = UgGridHelpers::
+                    getCoordinate(UgGridHelpers::increment(begin_cell_centroids, c, dimensions),
+                                  dimensions-1);
                 p[c] = press(z);
             }
         }
 
         // Initialize face pressures to distance-weighted average of adjacent cell pressures.
-        template <class State>
-        void initFacePressure(const UnstructuredGrid& grid,
+    template <class State, class FaceCells, class FCI, class CCI>
+        void initFacePressure(int dimensions,
+                              int number_of_faces,
+                              FaceCells face_cells,
+                              FCI begin_face_centroids,
+                              CCI begin_cell_centroids,
                               State& state)
         {
-            const int dim = grid.dimensions;
             const std::vector<double>& cp = state.pressure();
             std::vector<double>& fp = state.facepressure();
-            for (int f = 0; f < grid.number_of_faces; ++f) {
+            for (int f = 0; f < number_of_faces; ++f) {
                 double dist[2] = { 0.0, 0.0 };
                 double press[2] = { 0.0, 0.0 };
                 int bdy_idx = -1;
                 for (int j = 0; j < 2; ++j) {
-                    const int c = grid.face_cells[2*f + j];
+                    const int c = face_cells(f, j);
                     if (c >= 0) {
                         dist[j] = 0.0;
-                        for (int dd = 0; dd < dim; ++dd) {
-                            double diff = grid.face_centroids[dim*f + dd] - grid.cell_centroids[dim*c + dd];
+                        for (int dd = 0; dd < dimensions; ++dd) {
+                            double diff = UgGridHelpers::
+                                getCoordinate(UgGridHelpers::increment(begin_face_centroids, f,
+                                                                       dimensions),
+                                              dd)
+                                - UgGridHelpers::
+                                getCoordinate(UgGridHelpers::increment(begin_cell_centroids, c,
+                                                                       dimensions),
+                                              dd);
                             dist[j] += diff*diff;
                         }
                         dist[j] = std::sqrt(dist[j]);
@@ -320,11 +348,32 @@ namespace Opm
                         const double gravity,
                         State& state)
     {
+        initStateBasic(grid.number_of_cells, grid.global_cell, grid.cartdims,
+                       grid.number_of_faces, UgGridHelpers::faceCells(grid),
+                       grid.face_centroids, grid.cell_centroids,
+                       grid.dimensions, props, param,
+                       gravity, state);
+    }
+
+    template <class FaceCells, class CCI, class FCI, class State>
+    void initStateBasic(int number_of_cells,
+                        const int* global_cell,
+                        const int* cartdims,
+                        int number_of_faces,
+                        FaceCells face_cells,
+                        FCI begin_face_centroids,
+                        CCI begin_cell_centroids,
+                        int dimensions,
+                        const IncompPropertiesInterface& props,
+                        const parameter::ParameterGroup& param,
+                        const double gravity,
+                        State& state)
+{
         const int num_phases = props.numPhases();
         if (num_phases != 2) {
             OPM_THROW(std::runtime_error, "initStateTwophaseBasic(): currently handling only two-phase scenarios.");
         }
-        state.init(grid, num_phases);
+        state.init(number_of_cells, number_of_faces, num_phases);
         const int num_cells = props.numCells();
         // By default: initialise water saturation to minimum everywhere.
         std::vector<int> all_cells(num_cells);
@@ -338,11 +387,9 @@ namespace Opm
             // Initialise water saturation to max in the 'left' part.
             std::vector<int> left_cells;
             left_cells.reserve(num_cells/2);
-            const int *glob_cell = grid.global_cell;
-            const int* cd = grid.cartdims;
             for (int cell = 0; cell < num_cells; ++cell) {
-                const int gc = glob_cell == 0 ? cell : glob_cell[cell];
-                bool left = (gc % cd[0]) < cd[0]/2;
+                const int gc = global_cell == 0 ? cell : global_cell[cell];
+                bool left = (gc % cartdims[0]) < cartdims[0]/2;
                 if (left) {
                     left_cells.push_back(cell);
                 }
@@ -355,30 +402,34 @@ namespace Opm
             if (gravity == 0.0) {
                 std::cout << "**** Warning: running gravity segregation scenario, but gravity is zero." << std::endl;
             }
-            if (grid.cartdims[2] <= 1) {
+            if (cartdims[2] <= 1) {
                 std::cout << "**** Warning: running gravity segregation scenario, which expects nz > 1." << std::endl;
             }
             // Initialise water saturation to max *above* water-oil contact.
             const double woc = param.get<double>("water_oil_contact");
-            initWaterOilContact(grid, props, woc, WaterAbove, state);
+            initWaterOilContact(number_of_cells, begin_cell_centroids, dimensions,
+                                props, woc, WaterAbove, state);
             // Initialise pressure to hydrostatic state.
             const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa;
             double dens[2] = { props.density()[1], props.density()[0] };
-            initHydrostaticPressure(grid, dens, woc, gravity, woc, ref_p, state);
+            initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions,
+                                    dens, woc, gravity, woc, ref_p, state);
         } else if (param.has("water_oil_contact")) {
             // Warn against error-prone usage.
             if (gravity == 0.0) {
                 std::cout << "**** Warning: running gravity convection scenario, but gravity is zero." << std::endl;
             }
-            if (grid.cartdims[2] <= 1) {
+            if (cartdims[2] <= 1) {
                 std::cout << "**** Warning: running gravity convection scenario, which expects nz > 1." << std::endl;
             }
             // Initialise water saturation to max below water-oil contact.
             const double woc = param.get<double>("water_oil_contact");
-            initWaterOilContact(grid, props, woc, WaterBelow, state);
+            initWaterOilContact(number_of_cells, begin_cell_centroids, dimensions,
+                                props, woc, WaterBelow, state);
             // Initialise pressure to hydrostatic state.
             const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa;
-            initHydrostaticPressure(grid, props.density(), woc, gravity, woc, ref_p, state);
+            initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions,
+                                    props.density(), woc, gravity, woc, ref_p, state);
         } else if (param.has("init_saturation")) {
             // Initialise water saturation to init_saturation parameter.
             const double init_saturation = param.get<double>("init_saturation");
@@ -390,20 +441,23 @@ namespace Opm
             const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa;
             const double rho =  props.density()[0]*init_saturation + props.density()[1]*(1.0 - init_saturation);
             const double dens[2] = { rho, rho };
-            const double ref_z = grid.cell_centroids[0 + grid.dimensions - 1];
-            initHydrostaticPressure(grid, dens, ref_z, gravity, ref_z, ref_p, state);
+            const double ref_z = UgGridHelpers::getCoordinate(begin_cell_centroids, dimensions - 1);
+            initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions,
+                                    dens, ref_z, gravity, ref_z, ref_p, state);
         } else {
             // Use default: water saturation is minimum everywhere.
             // Initialise pressure to hydrostatic state.
             const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa;
             const double rho =  props.density()[1];
             const double dens[2] = { rho, rho };
-            const double ref_z = grid.cell_centroids[0 + grid.dimensions - 1];
-            initHydrostaticPressure(grid, dens, ref_z, gravity, ref_z, ref_p, state);
+            const double ref_z = UgGridHelpers::getCoordinate(begin_cell_centroids, dimensions - 1);
+            initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions,
+                                    dens, ref_z, gravity, ref_z, ref_p, state);
         }
 
         // Finally, init face pressures.
-        initFacePressure(grid, state);
+        initFacePressure(dimensions, number_of_faces, face_cells, begin_face_centroids,
+                         begin_cell_centroids, state);
     }
 
 
@@ -414,13 +468,33 @@ namespace Opm
                         const parameter::ParameterGroup& param,
                         const double gravity,
                         State& state)
+    {
+        initStateBasic(grid.number_of_cells, grid.global_cell, grid.cartdims,
+                       grid.number_of_faces, UgGridHelpers::faceCells(grid), 
+                       grid.face_centroids, grid.cell_centroids, grid.dimensions,
+                       props, param, gravity, state);
+    }
+
+    template <class FaceCells, class FCI, class CCI, class State>
+    void initStateBasic(int number_of_cells,
+                        const int* global_cell,
+                        const int* cartdims,
+                        int number_of_faces,
+                        FaceCells face_cells,
+                        FCI begin_face_centroids,
+                        CCI begin_cell_centroids,
+                        int dimensions,
+                        const BlackoilPropertiesInterface& props,
+                        const parameter::ParameterGroup& param,
+                        const double gravity,
+                        State& state)
     {
         // TODO: Refactor to exploit similarity with IncompProp* case.
         const int num_phases = props.numPhases();
         if (num_phases != 2) {
             OPM_THROW(std::runtime_error, "initStateTwophaseBasic(): currently handling only two-phase scenarios.");
         }
-        state.init(grid, num_phases);
+        state.init(number_of_cells, number_of_faces, num_phases);
         const int num_cells = props.numCells();
         // By default: initialise water saturation to minimum everywhere.
         std::vector<int> all_cells(num_cells);
@@ -433,11 +507,9 @@ namespace Opm
             // Initialise water saturation to max in the 'left' part.
             std::vector<int> left_cells;
             left_cells.reserve(num_cells/2);
-            const int *glob_cell = grid.global_cell;
-            const int* cd = grid.cartdims;
             for (int cell = 0; cell < num_cells; ++cell) {
-                const int gc = glob_cell == 0 ? cell : glob_cell[cell];
-                bool left = (gc % cd[0]) < cd[0]/2;
+                const int gc = global_cell == 0 ? cell : global_cell[cell];
+                bool left = (gc % cartdims[0]) < cartdims[0]/2;
                 if (left) {
                     left_cells.push_back(cell);
                 }
@@ -450,26 +522,30 @@ namespace Opm
             if (gravity == 0.0) {
                 std::cout << "**** Warning: running gravity convection scenario, but gravity is zero." << std::endl;
             }
-            if (grid.cartdims[2] <= 1) {
+            if (cartdims[2] <= 1) {
                 std::cout << "**** Warning: running gravity convection scenario, which expects nz > 1." << std::endl;
             }
             // Initialise water saturation to max below water-oil contact.
             const double woc = param.get<double>("water_oil_contact");
-            initWaterOilContact(grid, props, woc, WaterBelow, state);
+            initWaterOilContact(number_of_cells, begin_cell_centroids, dimensions,
+                                props, woc, WaterBelow, state);
             // Initialise pressure to hydrostatic state.
             const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa;
-            initHydrostaticPressure(grid, props, woc, gravity, woc, ref_p, state);
+            initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions,
+                                    props, woc, gravity, woc, ref_p, state);
         } else {
             // Use default: water saturation is minimum everywhere.
             // Initialise pressure to hydrostatic state.
             const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa;
-            const double ref_z = grid.cell_centroids[0 + grid.dimensions - 1];
+            const double ref_z = UgGridHelpers::getCoordinate(begin_cell_centroids, dimensions - 1);
             const double woc = -1e100;
-            initHydrostaticPressure(grid, props, woc, gravity, ref_z, ref_p, state);
+            initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions,
+                                    props, woc, gravity, ref_z, ref_p, state);
         }
 
         // Finally, init face pressures.
-        initFacePressure(grid, state);
+        initFacePressure(dimensions, number_of_faces, face_cells, begin_face_centroids,
+                         begin_cell_centroids, state);
     }
 
 
@@ -480,6 +556,25 @@ namespace Opm
                            Opm::DeckConstPtr newParserDeck,
                            const double gravity,
                            State& state)
+    {
+        initStateFromDeck(grid.number_of_cells, grid.global_cell,
+                          grid.number_of_faces, UgGridHelpers::faceCells(grid),
+                          grid.face_centroids, grid.cell_centroids, grid.dimensions,
+                          props, newParserDeck, gravity, state);
+    }
+    /// Initialize a state from input deck.
+    template <class FaceCells, class FCI, class CCI, class Props, class State>
+    void initStateFromDeck(int number_of_cells,
+                           const int* global_cell,
+                           int number_of_faces,
+                           FaceCells face_cells,
+                           FCI begin_face_centroids,
+                           CCI begin_cell_centroids,
+                           int dimensions,
+                           const Props& props,
+                           Opm::DeckConstPtr newParserDeck,
+                           const double gravity,
+                           State& state)
     {
         const int num_phases = props.numPhases();
         const PhaseUsage pu = phaseUsageFromDeck(newParserDeck);
@@ -487,7 +582,7 @@ namespace Opm
             OPM_THROW(std::runtime_error, "initStateFromDeck():  user specified property object with " << num_phases << " phases, "
                   "found " << pu.num_phases << " phases in deck.");
         }
-        state.init(grid, num_phases);
+        state.init(number_of_cells, number_of_faces, num_phases);
         if (newParserDeck->hasKeyword("EQUIL") && newParserDeck->hasKeyword("PRESSURE")) {
             OPM_THROW(std::runtime_error, "initStateFromDeck(): The deck must either specify the initial "
                       "condition using the PRESSURE _or_ the EQUIL keyword (currently it has both)");
@@ -506,17 +601,19 @@ namespace Opm
                 OPM_THROW(std::runtime_error, "initStateFromDeck(): No region support yet.");
             }
             const double woc = equil.waterOilContactDepth(0);
-            initWaterOilContact(grid, props, woc, WaterBelow, state);
+            initWaterOilContact(number_of_cells, begin_cell_centroids, dimensions,
+                                props, woc, WaterBelow, state);
             // Set pressure depending on densities and depths.
             const double datum_z = equil.datumDepth(0);
             const double datum_p = equil.datumDepthPressure(0);
-            initHydrostaticPressure(grid, props, woc, gravity, datum_z, datum_p, state);
+            initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions,
+                                    props, woc, gravity, datum_z, datum_p, state);
         } else if (newParserDeck->hasKeyword("PRESSURE")) {
             // Set saturations from SWAT/SGAS, pressure from PRESSURE.
             std::vector<double>& s = state.saturation();
             std::vector<double>& p = state.pressure();
             const std::vector<double>& p_deck = newParserDeck->getKeyword("PRESSURE")->getSIDoubleData();
-            const int num_cells = grid.number_of_cells;
+            const int num_cells = number_of_cells;
             if (num_phases == 2) {
                 if (!pu.phase_used[BlackoilPhases::Aqua]) {
                     // oil-gas: we require SGAS
@@ -527,7 +624,7 @@ namespace Opm
                     const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
                     const int opos = pu.phase_pos[BlackoilPhases::Liquid];
                     for (int c = 0; c < num_cells; ++c) {
-                        int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
+                        int c_deck = (global_cell == NULL) ? c : global_cell[c];
                         s[2*c + gpos] = sg_deck[c_deck];
                         s[2*c + opos] = 1.0 - sg_deck[c_deck];
                         p[c] = p_deck[c_deck];
@@ -541,7 +638,7 @@ namespace Opm
                     const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
                     const int nwpos = (wpos + 1) % 2;
                     for (int c = 0; c < num_cells; ++c) {
-                        int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
+                        int c_deck = (global_cell == NULL) ? c : global_cell[c];
                         s[2*c + wpos] = sw_deck[c_deck];
                         s[2*c + nwpos] = 1.0 - sw_deck[c_deck];
                         p[c] = p_deck[c_deck];
@@ -558,7 +655,7 @@ namespace Opm
                 const std::vector<double>& sw_deck = newParserDeck->getKeyword("SWAT")->getSIDoubleData();
                 const std::vector<double>& sg_deck = newParserDeck->getKeyword("SGAS")->getSIDoubleData();
                 for (int c = 0; c < num_cells; ++c) {
-                    int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
+                    int c_deck = (global_cell == NULL) ? c : global_cell[c];
                     s[3*c + wpos] = sw_deck[c_deck];
                     s[3*c + opos] = 1.0 - (sw_deck[c_deck] + sg_deck[c_deck]);
                     s[3*c + gpos] = sg_deck[c_deck];
@@ -572,7 +669,8 @@ namespace Opm
         }
 
         // Finally, init face pressures.
-        initFacePressure(grid, state);
+        initFacePressure(dimensions, number_of_faces, face_cells, begin_face_centroids,
+                         begin_cell_centroids, state);
     }
 
     /// Initialize a state from input deck.
@@ -582,6 +680,26 @@ namespace Opm
                            const EclipseGridParser& deck,
                            const double gravity,
                            State& state)
+    {
+        initStateFromDeck(grid.number_of_cells, grid.global_cell,
+                          grid.number_of_faces, UgGridHelpers::faceCells(grid),
+                          grid.face_centroids, grid.cell_centroids,
+                          grid.dimensions, props, deck,
+                          gravity, state);
+    }
+
+    template <class FaceCells, class FCI, class CCI, class Props, class State>
+    void initStateFromDeck(int number_of_cells,
+                           const int* global_cell,
+                           int number_of_faces,
+                           FaceCells face_cells,
+                           FCI begin_face_centroids,
+                           CCI begin_cell_centroids,
+                           int dimensions,
+                           const Props& props,
+                           const EclipseGridParser& deck,
+                           const double gravity,
+                           State& state)
     {
         const int num_phases = props.numPhases();
         const PhaseUsage pu = phaseUsageFromDeck(deck);
@@ -593,7 +711,7 @@ namespace Opm
             OPM_THROW(std::runtime_error, "initStateFromDeck(): The deck must either specify the initial "
                       "condition using the PRESSURE _or_ the EQUIL keyword (currently it has both)");
         }
-        state.init(grid, num_phases);
+        state.init(number_of_cells, number_of_faces, num_phases);
         if (deck.hasField("EQUIL")) {
             if (num_phases != 2) {
                 OPM_THROW(std::runtime_error, "initStateFromDeck(): EQUIL-based init currently handling only two-phase scenarios.");
@@ -607,17 +725,18 @@ namespace Opm
                 OPM_THROW(std::runtime_error, "initStateFromDeck(): No region support yet.");
             }
             const double woc = equil.equil[0].water_oil_contact_depth_;
-            initWaterOilContact(grid, props, woc, WaterBelow, state);
+            initWaterOilContact(number_of_cells, begin_cell_centroids, dimensions,
+                                props, woc, WaterBelow, state);
             // Set pressure depending on densities and depths.
             const double datum_z = equil.equil[0].datum_depth_;
             const double datum_p = equil.equil[0].datum_depth_pressure_;
-            initHydrostaticPressure(grid, props, woc, gravity, datum_z, datum_p, state);
+            initHydrostaticPressure(number_of_cells, begin_cell_centroids, dimensions,
+                                    props, woc, gravity, datum_z, datum_p, state);
         } else if (deck.hasField("PRESSURE")) {
             // Set saturations from SWAT/SGAS, pressure from PRESSURE.
             std::vector<double>& s = state.saturation();
             std::vector<double>& p = state.pressure();
             const std::vector<double>& p_deck = deck.getFloatingPointValue("PRESSURE");
-            const int num_cells = grid.number_of_cells;
             if (num_phases == 2) {
                 if (!pu.phase_used[BlackoilPhases::Aqua]) {
                     // oil-gas: we require SGAS
@@ -627,8 +746,8 @@ namespace Opm
                     const std::vector<double>& sg_deck = deck.getFloatingPointValue("SGAS");
                     const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
                     const int opos = pu.phase_pos[BlackoilPhases::Liquid];
-                    for (int c = 0; c < num_cells; ++c) {
-                        int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
+                    for (int c = 0; c < number_of_cells; ++c) {
+                        int c_deck = (global_cell == NULL) ? c : global_cell[c];
                         s[2*c + gpos] = sg_deck[c_deck];
                         s[2*c + opos] = 1.0 - sg_deck[c_deck];
                         p[c] = p_deck[c_deck];
@@ -641,8 +760,8 @@ namespace Opm
                     const std::vector<double>& sw_deck = deck.getFloatingPointValue("SWAT");
                     const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
                     const int nwpos = (wpos + 1) % 2;
-                    for (int c = 0; c < num_cells; ++c) {
-                        int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
+                    for (int c = 0; c < number_of_cells; ++c) {
+                        int c_deck = (global_cell == NULL) ? c : global_cell[c];
                         s[2*c + wpos] = sw_deck[c_deck];
                         s[2*c + nwpos] = 1.0 - sw_deck[c_deck];
                         p[c] = p_deck[c_deck];
@@ -658,8 +777,8 @@ namespace Opm
                 const int opos = pu.phase_pos[BlackoilPhases::Liquid];
                 const std::vector<double>& sw_deck = deck.getFloatingPointValue("SWAT");
                 const std::vector<double>& sg_deck = deck.getFloatingPointValue("SGAS");
-                for (int c = 0; c < num_cells; ++c) {
-                    int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
+                for (int c = 0; c < number_of_cells; ++c) {
+                    int c_deck = (global_cell == NULL) ? c : global_cell[c];
                     s[3*c + wpos] = sw_deck[c_deck];
                     s[3*c + opos] = 1.0 - (sw_deck[c_deck] + sg_deck[c_deck]);
                     s[3*c + gpos] = sg_deck[c_deck];
@@ -673,7 +792,8 @@ namespace Opm
         }
 
         // Finally, init face pressures.
-        initFacePressure(grid, state);
+        initFacePressure(dimensions, number_of_faces, face_cells, begin_face_centroids,
+                         begin_cell_centroids, state);
     }
 
     /// Initialize surface volume from pressure and saturation by z = As.
@@ -683,10 +803,18 @@ namespace Opm
     void initBlackoilSurfvol(const UnstructuredGrid& grid,
                              const Props& props,
                              State& state)
+    {
+        initBlackoilSurfvol(grid.number_of_cells, props, state);
+    }
+
+    template <class Props, class State>
+    void initBlackoilSurfvol(int number_of_cells,
+                             const Props& props,
+                             State& state)
     {
         state.surfacevol() = state.saturation();
         const int np = props.numPhases();
-        const int nc = grid.number_of_cells;
+        const int nc = number_of_cells;
         std::vector<double> allA(nc*np*np);
         std::vector<int> allcells(nc);
         for (int c = 0; c < nc; ++c) {
@@ -720,6 +848,15 @@ namespace Opm
                              const Props& props,
                              State& state)
     {
+        initBlackoilSurfvolUsingRSorRV(grid.number_of_cells, props, state);
+    }
+
+    template <class Props, class State>
+    void initBlackoilSurfvolUsingRSorRV(int number_of_cells,
+                                        const Props& props,
+                                        State& state)
+    {
+        
         if (props.numPhases() != 3) {
             OPM_THROW(std::runtime_error, "initBlackoilSurfvol() is only supported in three-phase simulations.");
         }
@@ -731,27 +868,26 @@ namespace Opm
         const PhaseUsage pu = props.phaseUsage();
 
         const int np = props.numPhases();
-        const int nc = grid.number_of_cells;
 
-        std::vector<double> allA_a(nc*np*np);
-        std::vector<double> allA_l(nc*np*np);
-        std::vector<double> allA_v(nc*np*np);
+        std::vector<double> allA_a(number_of_cells*np*np);
+        std::vector<double> allA_l(number_of_cells*np*np);
+        std::vector<double> allA_v(number_of_cells*np*np);
 
-        std::vector<int> allcells(nc);
-        std::vector<double> z_init(nc*np);
+        std::vector<int> allcells(number_of_cells);
+        std::vector<double> z_init(number_of_cells*np);
         std::fill(z_init.begin(),z_init.end(),0.0);
 
-        for (int c = 0; c < nc; ++c) {
+        for (int c = 0; c < number_of_cells; ++c) {
             allcells[c] = c;
         }
 
-        std::vector<double> capPressures(nc*np);
-        props.capPress(nc,&state.saturation()[0],&allcells[0],&capPressures[0],NULL);
+        std::vector<double> capPressures(number_of_cells*np);
+        props.capPress(number_of_cells,&state.saturation()[0],&allcells[0],&capPressures[0],NULL);
 
-        std::vector<double> Pw(nc);
-        std::vector<double> Pg(nc);
+        std::vector<double> Pw(number_of_cells);
+        std::vector<double> Pg(number_of_cells);
 
-        for (int c = 0; c < nc; ++c){
+        for (int c = 0; c < number_of_cells; ++c){
             Pw[c] = state.pressure()[c] + capPressures[c*np + BlackoilPhases::Aqua];
             Pg[c] = state.pressure()[c] + capPressures[c*np + BlackoilPhases::Vapour];
         }
@@ -761,7 +897,7 @@ namespace Opm
 
         // Water phase
         if(pu.phase_used[BlackoilPhases::Aqua])
-           for (int c = 0; c <  nc ; ++c){
+           for (int c = 0; c <  number_of_cells ; ++c){
                for (int p = 0; p < np ; ++p){
                    if (p == BlackoilPhases::Aqua)
                        z_tmp = 1;
@@ -771,11 +907,11 @@ namespace Opm
                    z_init[c*np + p] = z_tmp;
                }
            }
-        props.matrix(nc, &state.pressure()[0], &z_init[0], &allcells[0], &allA_a[0], 0);
+        props.matrix(number_of_cells, &state.pressure()[0], &z_init[0], &allcells[0], &allA_a[0], 0);
 
         // Liquid phase
         if(pu.phase_used[BlackoilPhases::Liquid]){
-            for (int c = 0; c <  nc ; ++c){
+            for (int c = 0; c <  number_of_cells ; ++c){
                 for (int p = 0; p < np ; ++p){
                      if(p == BlackoilPhases::Vapour){
                          if(state.saturation()[np*c + p] > 0)
@@ -793,10 +929,10 @@ namespace Opm
                 }
             }
         }
-        props.matrix(nc, &state.pressure()[0], &z_init[0], &allcells[0], &allA_l[0], 0);
+        props.matrix(number_of_cells, &state.pressure()[0], &z_init[0], &allcells[0], &allA_l[0], 0);
 
         if(pu.phase_used[BlackoilPhases::Vapour]){
-            for (int c = 0; c <  nc ; ++c){
+            for (int c = 0; c <  number_of_cells ; ++c){
                 for (int p = 0; p < np ; ++p){
                      if(p == BlackoilPhases::Liquid){
                          if(state.saturation()[np*c + p] > 0)
@@ -814,9 +950,9 @@ namespace Opm
                 }
             }
         }
-        props.matrix(nc, &state.pressure()[0], &z_init[0], &allcells[0], &allA_v[0], 0);
+        props.matrix(number_of_cells, &state.pressure()[0], &z_init[0], &allcells[0], &allA_v[0], 0);
 
-        for (int c = 0; c < nc; ++c) {
+        for (int c = 0; c < number_of_cells; ++c) {
             // Using z = As
             double* z = &state.surfacevol()[c*np];
             const double* s = &state.saturation()[c*np];
@@ -848,24 +984,43 @@ namespace Opm
                                    const double gravity,
                                    State& state)
     {
-        initStateFromDeck(grid, props, deck, gravity, state);
+        initBlackoilStateFromDeck(grid.number_of_cells, grid.global_cell,
+                                  grid.number_of_faces, UgGridHelpers::faceCells(grid),
+                                  grid.face_centroids, grid.cell_centroids,grid.dimensions,
+                                  props, deck, gravity, state);
+    }
+
+    template <class FaceCells, class FCI, class CCI, class Props, class State>
+    void initBlackoilStateFromDeck(int number_of_cells,
+                                   const int* global_cell,
+                                   int number_of_faces,
+                                   FaceCells face_cells,
+                                   FCI begin_face_centroids,
+                                   CCI begin_cell_centroids,
+                                   int dimensions,
+                                   const Props& props,
+                                   const EclipseGridParser& deck,
+                                   const double gravity,
+                                   State& state)
+    {
+        initStateFromDeck(number_of_cells, global_cell, number_of_faces,
+                          face_cells, begin_face_centroids, begin_cell_centroids,
+                          dimensions, props, deck, gravity, state);
         if (deck.hasField("RS")) {
             const std::vector<double>& rs_deck = deck.getFloatingPointValue("RS");
-            const int num_cells = grid.number_of_cells;
-            for (int c = 0; c < num_cells; ++c) {
-                int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
+            for (int c = 0; c < number_of_cells; ++c) {
+                int c_deck = (global_cell == NULL) ? c : global_cell[c];
                 state.gasoilratio()[c] = rs_deck[c_deck];
             }
-            initBlackoilSurfvolUsingRSorRV(grid, props, state);
+            initBlackoilSurfvolUsingRSorRV(number_of_cells, props, state);
             computeSaturation(props,state);
         } else if (deck.hasField("RV")){
             const std::vector<double>& rv_deck = deck.getFloatingPointValue("RV");
-            const int num_cells = grid.number_of_cells;
-            for (int c = 0; c < num_cells; ++c) {
-                int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
+            for (int c = 0; c < number_of_cells; ++c) {
+                int c_deck = (global_cell == NULL) ? c : global_cell[c];
                 state.rv()[c] = rv_deck[c_deck];
             }
-            initBlackoilSurfvolUsingRSorRV(grid, props, state);
+            initBlackoilSurfvolUsingRSorRV(number_of_cells, props, state);
             computeSaturation(props,state);
         }
 
@@ -882,24 +1037,46 @@ namespace Opm
                                    const double gravity,
                                    State& state)
     {
-        initStateFromDeck(grid, props, newParserDeck, gravity, state);
+        initBlackoilStateFromDeck(grid.number_of_cells, grid.global_cell,
+                                  grid.number_of_faces, UgGridHelpers::faceCells(grid),
+                                  grid.face_centroids, grid.cell_centroids,grid.dimensions,
+                                  props, newParserDeck, gravity, state);
+    }
+
+    /// Initialize a blackoil state from input deck.
+    template <class FaceCells, class FCI, class CCI, class Props, class State>
+    void initBlackoilStateFromDeck(int number_of_cells,
+                                   const int* global_cell,
+                                   int number_of_faces,
+                                   FaceCells face_cells,
+                                   FCI begin_face_centroids,
+                                   CCI begin_cell_centroids,
+                                   int dimensions,
+                                   const Props& props,
+                                   Opm::DeckConstPtr newParserDeck,
+                                   const double gravity,
+                                   State& state)
+    {
+        initStateFromDeck(number_of_cells, global_cell, number_of_faces,
+                          face_cells, begin_face_centroids, begin_cell_centroids,
+                          dimensions, props, newParserDeck, gravity, state);
         if (newParserDeck->hasKeyword("RS")) {
             const std::vector<double>& rs_deck = newParserDeck->getKeyword("RS")->getSIDoubleData();
-            const int num_cells = grid.number_of_cells;
+            const int num_cells = number_of_cells;
             for (int c = 0; c < num_cells; ++c) {
-                int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
+                int c_deck = (global_cell == NULL) ? c : global_cell[c];
                 state.gasoilratio()[c] = rs_deck[c_deck];
             }
-            initBlackoilSurfvolUsingRSorRV(grid, props, state);
+            initBlackoilSurfvolUsingRSorRV(number_of_cells, props, state);
             computeSaturation(props,state);
         } else if (newParserDeck->hasKeyword("RV")){
             const std::vector<double>& rv_deck = newParserDeck->getKeyword("RV")->getSIDoubleData();
-            const int num_cells = grid.number_of_cells;
+            const int num_cells = number_of_cells;
             for (int c = 0; c < num_cells; ++c) {
-                int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
+                int c_deck = (global_cell == NULL) ? c : global_cell[c];
                 state.rv()[c] = rv_deck[c_deck];
             }
-            initBlackoilSurfvolUsingRSorRV(grid, props, state);
+            initBlackoilSurfvolUsingRSorRV(number_of_cells, props, state);
             computeSaturation(props,state);
         }
         else {

opm/core/utility/miscUtilities.cpp

@@ -45,15 +45,10 @@ namespace Opm
                            const double* porosity,
                            std::vector<double>& porevol)
     {
-        int num_cells = grid.number_of_cells;
-        porevol.resize(num_cells);
-        std::transform(porosity, porosity + num_cells,
-                       grid.cell_volumes,
-                       porevol.begin(),
-                       std::multiplies<double>());
+        computePorevolume(grid.number_of_cells, grid.cell_volumes,
+                          porosity, porevol);
     }
 
-
     /// @brief Computes pore volume of all cells in a grid, with rock compressibility effects.
     /// @param[in]  grid      a grid
     /// @param[in]  porosity  array of grid.number_of_cells porosity values
@@ -66,13 +61,11 @@ namespace Opm
                            const std::vector<double>& pressure,
                            std::vector<double>& porevol)
     {
-        int num_cells = grid.number_of_cells;
-        porevol.resize(num_cells);
-        for (int i = 0; i < num_cells; ++i) {
-            porevol[i] = porosity[i]*grid.cell_volumes[i]*rock_comp.poroMult(pressure[i]);
-        }
+        computePorevolume(grid.number_of_cells, grid.cell_volumes, porosity, rock_comp, pressure,
+                          porevol);
     }
 
+
     /// @brief Computes porosity of all cells in a grid, with rock compressibility effects.
     /// @param[in]  grid               a grid
     /// @param[in]  porosity_standard  array of grid.number_of_cells porosity values (at standard conditions)
@@ -401,24 +394,16 @@ namespace Opm
                               const std::vector<double>& face_flux,
                               std::vector<double>& cell_velocity)
     {
-        const int dim = grid.dimensions;
-        cell_velocity.clear();
-        cell_velocity.resize(grid.number_of_cells*dim, 0.0);
-        for (int face = 0; face < grid.number_of_faces; ++face) {
-            int c[2] = { grid.face_cells[2*face], grid.face_cells[2*face + 1] };
-            const double* fc = &grid.face_centroids[face*dim];
-            double flux = face_flux[face];
-            for (int i = 0; i < 2; ++i) {
-                if (c[i] >= 0) {
-                    const double* cc = &grid.cell_centroids[c[i]*dim];
-                    for (int d = 0; d < dim; ++d) {
-                        double v_contrib = fc[d] - cc[d];
-                        v_contrib *= flux/grid.cell_volumes[c[i]];
-                        cell_velocity[c[i]*dim + d] += (i == 0) ? v_contrib : -v_contrib;
-                    }
-                }
-            }
-        }
+        estimateCellVelocity(grid.number_of_cells,
+                             grid.number_of_faces,
+                             grid.face_centroids,
+                             UgGridHelpers::faceCells(grid),
+                             grid.cell_centroids,
+                             grid.cell_volumes,
+                             grid.dimensions,
+                             face_flux,
+                             cell_velocity);
+        
     }
 
     /// Extract a vector of water saturations from a vector of
@@ -490,44 +475,8 @@ namespace Opm
                     const double* densities, const double gravity, const bool per_grid_cell,
                     std::vector<double>& wdp)
     {
-        const int nw = wells.number_of_wells;
-        const size_t np = per_grid_cell ?
-            saturations.size()/grid.number_of_cells
-            : saturations.size()/wells.well_connpos[nw];
-        // Simple for now:
-        for (int i = 0; i < nw; i++) {
-            double depth_ref = wells.depth_ref[i];
-            for (int j = wells.well_connpos[i]; j < wells.well_connpos[i + 1]; j++) {
-                int cell = wells.well_cells[j];
-
-                // Is this correct wrt. depth_ref?
-                double cell_depth = grid.cell_centroids[3 * cell + 2];
-
-                double saturation_sum = 0.0;
-                for (size_t p = 0; p < np; p++) {
-                    if (!per_grid_cell) {
-                        saturation_sum += saturations[j * np + p];
-                    } else {
-                        saturation_sum += saturations[np * cell + p];
-                    }
-                }
-                if (saturation_sum == 0) {
-                    saturation_sum = 1.0;
-                }
-                double density = 0.0;
-                for (size_t p = 0; p < np; p++) {
-                    if (!per_grid_cell) {
-                        density += saturations[j * np + p] * densities[p] / saturation_sum;
-                    } else {
-                        // Is this a smart way of doing it?
-                        density += saturations[np * cell + p] * densities[p] / saturation_sum;
-                    }
-                }
-
-                // Is the sign correct?
-                wdp.push_back(density * (cell_depth - depth_ref) * gravity);
-            }
-        }
+        computeWDP(wells, grid.number_of_cells, grid.cell_centroids, saturations, densities,
+                   gravity, per_grid_cell, wdp);
     }
 
 

opm/core/utility/miscUtilities.hpp

@@ -41,6 +41,16 @@ namespace Opm
                            const double* porosity,
 			   std::vector<double>& porevol);
 
+    /// @brief Computes pore volume of all cells in a grid.
+    /// @param[in]  number_of_cells The number of cells of the grid.
+    /// @param[in]  begin_cell_volume Iterator to the volume of the first cell.
+    /// @param[in]  porosity  array of grid.number_of_cells porosity values
+    /// @param[out] porevol   the pore volume by cell.
+    template<class T>
+    void computePorevolume(int number_of_cells,
+                           T begin_cell_volume,
+                           const double* porosity,
+                           std::vector<double>& porevol);
 
     /// @brief Computes pore volume of all cells in a grid, with rock compressibility effects.
     /// @param[in]  grid      a grid
@@ -54,6 +64,21 @@ namespace Opm
                            const std::vector<double>& pressure,
                            std::vector<double>& porevol);
 
+    /// @brief Computes pore volume of all cells in a grid, with rock compressibility effects.
+    /// @param[in]  number_of_cells The number of cells of the grid.
+    /// @param[in]  Pointer to/ Iterator at the first cell volume.
+    /// @param[in]  porosity  array of grid.number_of_cells porosity values
+    /// @param[in]  rock_comp rock compressibility properties
+    /// @param[in]  pressure  pressure by cell
+    /// @param[out] porevol   the pore volume by cell.
+    template<class T>
+    void computePorevolume(int number_of_cells,
+                           T begin_cell_volume,
+                           const double* porosity,
+                           const RockCompressibility& rock_comp,
+                           const std::vector<double>& pressure,
+                           std::vector<double>& porevol);
+
     /// @brief Computes porosity of all cells in a grid, with rock compressibility effects.
     /// @param[in]  grid               a grid
     /// @param[in]  porosity_standard  array of grid.number_of_cells porosity values (at reference presure)
@@ -188,6 +213,26 @@ namespace Opm
 			      const std::vector<double>& face_flux,
 			      std::vector<double>& cell_velocity);
 
+    /// @brief Estimates a scalar cell velocity from face fluxes.
+    /// @param[in]  number_of_cells      The number of cells of the grid
+    /// @param[in]  number_of_faces      The number of cells of the grid
+    /// @param[in]  begin_face_centroids Iterator pointing to first face centroid.
+    /// @param[in]  face_cells           Mapping from faces to connected cells.
+    /// @param[in]  dimensions           The dimensions of the grid.
+    /// @param[in]  begin_cell_centroids Iterator pointing to first cell centroid.
+    /// @param[in]  face_flux            signed per-face fluxes
+    /// @param[out] cell_velocity        the estimated velocities.
+    template<class CC, class FC, class FC1, class CV>
+    void estimateCellVelocity(int number_of_cells,
+                              int number_of_faces,
+                              FC begin_face_centroids,
+                              FC1 face_cells,
+                              CC begin_cell_centroids,
+                              CV begin_cell_volumes,
+                              int dimension,
+                              const std::vector<double>& face_flux,
+                              std::vector<double>& cell_velocity);
+
     /// Extract a vector of water saturations from a vector of
     /// interleaved water and oil saturations.
     void toWaterSat(const std::vector<double>& sboth,
@@ -218,6 +263,24 @@ namespace Opm
                     const double* densities, const double gravity, const bool per_grid_cell,
                     std::vector<double>& wdp);
 
+    /// Computes the WDP for each well.
+    /// \param[in] wells        Wells that need their wdp calculated.
+    /// \param[in] number_of_cells The number of cells in the grid.
+    /// \param[in] begin_cell_centroids Pointer/Iterator to the first cell centroid. 
+    /// \param[in] saturations  A vector of weights for each cell for each phase
+    ///                         in the grid (or well, see per_grid_cell parameter). So for cell i,
+    ///                         saturations[i*densities.size() + p] should give the weight
+    ///                         of phase p in cell i.
+    /// \param[in] densities    Density for each phase.
+    /// \param[out] wdp         Will contain, for each well, the wdp of the well.
+    /// \param[in] per_grid_cell Whether or not the saturations are per grid cell or per
+    ///                          well cell.
+    template<class T>
+    void computeWDP(const Wells& wells, int number_of_cells, T begin_cell_centroids,
+                    const std::vector<double>& saturations,
+                    const double* densities, const double gravity, const bool per_grid_cell,
+                    std::vector<double>& wdp);
+
     /// Computes (sums) the flow rate for each well.
     /// \param[in] wells                The wells for which the flow rate should be computed.
     /// \param[in] flow_rates_per_cell  Flow rates per well cells. Should ordered the same way as
@@ -320,4 +383,5 @@ namespace Opm
 
 } // namespace Opm
 
+#include "miscUtilities_impl.hpp"
 #endif // OPM_MISCUTILITIES_HEADER_INCLUDED

opm/core/utility/miscUtilities_impl.hpp

@@ -0,0 +1,123 @@
+#include <opm/core/grid/GridHelpers.hpp>
+#include <opm/core/wells.h>
+#include <opm/core/props/rock/RockCompressibility.hpp>
+
+namespace Opm
+{
+    /// @brief Estimates a scalar cell velocity from face fluxes.
+    /// @param[in]  number_of_cells      The number of cells of the grid
+    /// @param[in]  begin_face_centroids Iterator pointing to first face centroid.
+    /// @param[in]  face_cells           Mapping from faces to connected cells.
+    /// @param[in]  dimensions           The dimensions of the grid.
+    /// @param[in]  begin_cell_centroids Iterator pointing to first cell centroid.
+    /// @param[in]  face_flux            signed per-face fluxes
+    /// @param[out] cell_velocity        the estimated velocities.
+    template<class CC, class FC, class FC1, class CV>
+    void estimateCellVelocity(int number_of_cells,
+                              int number_of_faces,
+                              FC begin_face_centroids,
+                              FC1 face_cells,
+                              CC begin_cell_centroids,
+                              CV begin_cell_volumes,
+                              int dimension,
+                              const std::vector<double>& face_flux,
+                              std::vector<double>& cell_velocity)
+    {
+        cell_velocity.clear();
+        cell_velocity.resize(number_of_cells*dimension, 0.0);
+        for (int face = 0; face < number_of_faces; ++face) {
+            int c[2] = { face_cells(face, 0), face_cells(face, 1) };
+            FC fc = UgGridHelpers::increment(begin_face_centroids, face, dimension);
+            double flux = face_flux[face];
+            for (int i = 0; i < 2; ++i) {
+                if (c[i] >= 0) {
+                    CC cc = UgGridHelpers::increment(begin_cell_centroids, c[i], dimension);
+                    for (int d = 0; d < dimension; ++d) {
+                        double v_contrib = UgGridHelpers::getCoordinate(fc, d) - UgGridHelpers::getCoordinate(cc, d);
+                        v_contrib *= flux/begin_cell_volumes[c[i]];
+                        cell_velocity[c[i]*dimension + d] += (i == 0) ? v_contrib : -v_contrib;
+                    }
+                }
+            }
+        }
+    }
+
+    template<class T>
+    void computePorevolume(int number_of_cells,
+                           T begin_cell_volume,
+                           const double* porosity,
+                           std::vector<double>& porevol)
+    {
+        porevol.resize(number_of_cells);
+        std::transform(porosity, porosity + number_of_cells,
+                       begin_cell_volume,
+                       porevol.begin(),
+                       std::multiplies<double>());
+    }
+
+    /// @brief Computes pore volume of all cells in a grid, with rock compressibility effects.
+    /// @param[in]  number_of_cells The number of cells of the grid.
+    /// @param[in]  porosity  array of grid.number_of_cells porosity values
+    /// @param[in]  rock_comp rock compressibility properties
+    /// @param[in]  pressure  pressure by cell
+    /// @param[out] porevol   the pore volume by cell.
+    template<class T>
+    void computePorevolume(int number_of_cells,
+                           T begin_cell_volumes,
+                           const double* porosity,
+                           const RockCompressibility& rock_comp,
+                           const std::vector<double>& pressure,
+                           std::vector<double>& porevol)
+    {
+        porevol.resize(number_of_cells);
+        for (int i = 0; i < number_of_cells; ++i) {
+            porevol[i] = porosity[i]*begin_cell_volumes[i]*rock_comp.poroMult(pressure[i]);
+        }
+    }
+
+    template<class T>
+    void computeWDP(const Wells& wells, int number_of_cells, T begin_cell_centroids, const std::vector<double>& saturations,
+                    const double* densities, const double gravity, const bool per_grid_cell,
+                    std::vector<double>& wdp)
+    {
+        const int nw = wells.number_of_wells;
+        const size_t np = per_grid_cell ?
+            saturations.size()/number_of_cells
+            : saturations.size()/wells.well_connpos[nw];
+        // Simple for now:
+        for (int i = 0; i < nw; i++) {
+            double depth_ref = wells.depth_ref[i];
+            for (int j = wells.well_connpos[i]; j < wells.well_connpos[i + 1]; j++) {
+                int cell = wells.well_cells[j];
+
+                // Is this correct wrt. depth_ref?
+                double cell_depth = UgGridHelpers
+                    ::getCoordinate(UgGridHelpers::increment(begin_cell_centroids, cell, 3), 2);
+
+                double saturation_sum = 0.0;
+                for (size_t p = 0; p < np; p++) {
+                    if (!per_grid_cell) {
+                        saturation_sum += saturations[j * np + p];
+                    } else {
+                        saturation_sum += saturations[np * cell + p];
+                    }
+                }
+                if (saturation_sum == 0) {
+                    saturation_sum = 1.0;
+                }
+                double density = 0.0;
+                for (size_t p = 0; p < np; p++) {
+                    if (!per_grid_cell) {
+                        density += saturations[j * np + p] * densities[p] / saturation_sum;
+                    } else {
+                        // Is this a smart way of doing it?
+                        density += saturations[np * cell + p] * densities[p] / saturation_sum;
+                    }
+                }
+
+                // Is the sign correct?
+                wdp.push_back(density * (cell_depth - depth_ref) * gravity);
+            }
+        }
+    }
+}

opm/core/wells/WellsManager.cpp

@@ -26,13 +26,11 @@
 #include <opm/core/wells.h>
 #include <opm/core/well_controls.h>
 #include <opm/core/utility/ErrorMacros.hpp>
-#include <opm/core/utility/Units.hpp>
 #include <opm/core/wells/WellCollection.hpp>
 #include <opm/core/props/phaseUsageFromDeck.hpp>
 
 #include <opm/parser/eclipse/EclipseState/Schedule/ScheduleEnums.hpp>
 
-#include <array>
 #include <algorithm>
 #include <cassert>
 #include <cmath>
@@ -44,17 +42,13 @@
 
 
 // Helper structs and functions for the implementation.
-namespace
+namespace WellsManagerDetail
 {
 
 
 
     namespace ProductionControl
     {
-        enum Mode { ORAT, WRAT, GRAT,
-                    LRAT, CRAT, RESV,
-                    BHP , THP , GRUP };
-
         namespace Details {
             std::map<std::string, Mode>
             init_mode_map() {
@@ -122,8 +116,6 @@ namespace
 
     namespace InjectionControl
     {
-        enum Mode { RATE, RESV, BHP,
-                    THP, GRUP };
 
         namespace Details {
             std::map<std::string, Mode>
@@ -177,28 +169,6 @@ namespace
 
     } // namespace InjectionControl
 
-
-    std::array<double, 3> getCubeDim(const UnstructuredGrid& grid, int cell)
-    {
-        using namespace std;
-        std::array<double, 3> cube;
-        int num_local_faces = grid.cell_facepos[cell + 1] - grid.cell_facepos[cell];
-        vector<double> x(num_local_faces);
-        vector<double> y(num_local_faces);
-        vector<double> z(num_local_faces);
-        for (int lf=0; lf<num_local_faces; ++ lf) {
-            int face = grid.cell_faces[grid.cell_facepos[cell] + lf];
-            const double* centroid = &grid.face_centroids[grid.dimensions*face];
-            x[lf] = centroid[0];
-            y[lf] = centroid[1];
-            z[lf] = centroid[2];
-        }
-        cube[0] = *max_element(x.begin(), x.end()) - *min_element(x.begin(), x.end());
-        cube[1] = *max_element(y.begin(), y.end()) - *min_element(y.begin(), y.end());
-        cube[2] = *max_element(z.begin(), z.end()) - *min_element(z.begin(), z.end());
-        return cube;
-    }
-
     // Use the Peaceman well model to compute well indices.
     // radius is the radius of the well.
     // cubical contains [dx, dy, dz] of the cell.
@@ -261,7 +231,6 @@ namespace Opm
     {
     }
 
-
     /// Construct from existing wells object.
     WellsManager::WellsManager(struct Wells* W)
         : w_(clone_wells(W))
@@ -275,83 +244,14 @@ namespace Opm
                                const double* permeability)
         : w_(0)
     {
-        if (grid.dimensions != 3) {
-            OPM_THROW(std::runtime_error, "We cannot initialize wells from a deck unless the corresponding grid is 3-dimensional.");
+        init(eclipseState, timeStep, UgGridHelpers::numCells(grid),
+             UgGridHelpers::globalCell(grid), UgGridHelpers::cartDims(grid), 
+             UgGridHelpers::dimensions(grid), UgGridHelpers::beginCellCentroids(grid),
+             UgGridHelpers::cell2Faces(grid), UgGridHelpers::beginFaceCentroids(grid),
+             permeability);
+
     }
 
-        if (eclipseState->getSchedule()->numWells() == 0) {
-            OPM_MESSAGE("No wells specified in Schedule section, initializing no wells");
-            return;
-        }
-
-        std::map<int,int> cartesian_to_compressed;
-        setupCompressedToCartesian(grid, cartesian_to_compressed);
-
-        // Obtain phase usage data.
-        PhaseUsage pu = phaseUsageFromDeck(eclipseState);
-
-        // These data structures will be filled in this constructor,
-        // then used to initialize the Wells struct.
-        std::vector<std::string> well_names;
-        std::vector<WellData> well_data;
-
-
-        // For easy lookup:
-        std::map<std::string, int> well_names_to_index;
-
-        ScheduleConstPtr schedule = eclipseState->getSchedule();
-        std::vector<WellConstPtr> wells = schedule->getWells(timeStep);
-
-        well_names.reserve(wells.size());
-        well_data.reserve(wells.size());
-
-        createWellsFromSpecs(wells, timeStep, grid, well_names, well_data, well_names_to_index, pu, cartesian_to_compressed, permeability);
-        setupWellControls(wells, timeStep, well_names, pu);
-
-        {
-            GroupTreeNodeConstPtr fieldNode = eclipseState->getSchedule()->getGroupTree(timeStep)->getNode("FIELD");
-            GroupConstPtr fieldGroup = eclipseState->getSchedule()->getGroup(fieldNode->name());
-            well_collection_.addField(fieldGroup, timeStep, pu);
-            addChildGroups(fieldNode, eclipseState->getSchedule(), timeStep, pu);
-        }
-
-        for (auto wellIter = wells.begin(); wellIter != wells.end(); ++wellIter ) {
-            well_collection_.addWell((*wellIter), timeStep, pu);
-        }
-
-        well_collection_.setWellsPointer(w_);
-        well_collection_.applyGroupControls();
-
-
-        setupGuideRates(wells, timeStep, well_data, well_names_to_index);
-
-        // Debug output.
-#define EXTRA_OUTPUT
-#ifdef EXTRA_OUTPUT
-        /*
-        std::cout << "\t WELL DATA" << std::endl;
-        for(int i = 0; i< num_wells; ++i) {
-            std::cout << i << ": " << well_data[i].type << "  "
-                      << well_data[i].control << "  " << well_data[i].target
-                      << std::endl;
-        }
-
-        std::cout << "\n\t PERF DATA" << std::endl;
-        for(int i=0; i< int(wellperf_data.size()); ++i) {
-            for(int j=0; j< int(wellperf_data[i].size()); ++j) {
-                std::cout << i << ": " << wellperf_data[i][j].cell << "  "
-                          << wellperf_data[i][j].well_index << std::endl;
-            }
-        }
-        */
-#endif
-    }
-
-
-
-
-
-
     /// Destructor.
     WellsManager::~WellsManager()
     {
@@ -405,141 +305,25 @@ namespace Opm
         well_collection_.applyExplicitReinjectionControls(well_reservoirrates_phase, well_surfacerates_phase);
     }
 
-    void WellsManager::setupCompressedToCartesian(const UnstructuredGrid& grid, std::map<int,int>& cartesian_to_compressed ) {
+    void WellsManager::setupCompressedToCartesian(const int* global_cell, int number_of_cells,
+                                                  std::map<int,int>& cartesian_to_compressed ) {
         // global_cell is a map from compressed cells to Cartesian grid cells.
         // We must make the inverse lookup.
-        const int* global_cell = grid.global_cell;
 
         if (global_cell) {
-            for (int i = 0; i < grid.number_of_cells; ++i) {
+            for (int i = 0; i < number_of_cells; ++i) {
                 cartesian_to_compressed.insert(std::make_pair(global_cell[i], i));
             }
         }
         else {
-            for (int i = 0; i < grid.number_of_cells; ++i) {
+            for (int i = 0; i < number_of_cells; ++i) {
                 cartesian_to_compressed.insert(std::make_pair(i, i));
             }
         }
 
     }
 
-    void WellsManager::createWellsFromSpecs(std::vector<WellConstPtr>& wells, size_t timeStep,
-                                            const UnstructuredGrid& grid,
-                                            std::vector<std::string>& well_names,
-                                            std::vector<WellData>& well_data,
-                                            std::map<std::string, int>& well_names_to_index,
-                                            const PhaseUsage& phaseUsage,
-                                            std::map<int,int> cartesian_to_compressed,
-                                            const double* permeability)
-    {
-        std::vector<std::vector<PerfData> > wellperf_data;
-        wellperf_data.resize(wells.size());
 
-        int well_index = 0;
-        for (auto wellIter= wells.begin(); wellIter != wells.end(); ++wellIter) {
-            WellConstPtr well = (*wellIter);
-            {   // WELSPECS handling
-                well_names_to_index[well->name()] = well_index;
-                well_names.push_back(well->name());
-                {
-                    WellData wd;
-                    // If negative (defaulted), set refdepth to a marker
-                    // value, will be changed after getting perforation
-                    // data to the centroid of the cell of the top well
-                    // perforation.
-                    wd.reference_bhp_depth = (well->getRefDepth() < 0.0) ? -1e100 : well->getRefDepth();
-                    wd.welspecsline = -1;
-                    if (well->isInjector( timeStep ))
-                        wd.type = INJECTOR;
-                    else
-                        wd.type = PRODUCER;
-                    well_data.push_back(wd);
-                }
-            }
-
-            {   // COMPDAT handling
-                CompletionSetConstPtr completionSet = well->getCompletions(timeStep);
-                for (size_t c=0; c<completionSet->size(); c++) {
-                    CompletionConstPtr completion = completionSet->get(c);
-                    int i = completion->getI();
-                    int j = completion->getJ();
-                    int k = completion->getK();
-
-                    const int* cpgdim = grid.cartdims;
-                    int cart_grid_indx = i + cpgdim[0]*(j + cpgdim[1]*k);
-                    std::map<int, int>::const_iterator cgit = cartesian_to_compressed.find(cart_grid_indx);
-                    if (cgit == cartesian_to_compressed.end()) {
-                        OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << i << ' ' << j << ' '
-                                  << k << " not found in grid (well = " << well->name() << ')');
-                    }
-                    int cell = cgit->second;
-                    PerfData pd;
-                    pd.cell = cell;
-                    if (completion->getCF() > 0.0) {
-                        pd.well_index = completion->getCF();
-                    } else {
-                        double radius = 0.5*completion->getDiameter();
-                        if (radius <= 0.0) {
-                            radius = 0.5*unit::feet;
-                            OPM_MESSAGE("**** Warning: Well bore internal radius set to " << radius);
-                        }
-                        std::array<double, 3> cubical = getCubeDim(grid, cell);
-                        const double* cell_perm = &permeability[grid.dimensions*grid.dimensions*cell];
-                        pd.well_index = computeWellIndex(radius, cubical, cell_perm, completion->getSkinFactor());
-                    }
-                    wellperf_data[well_index].push_back(pd);
-                }
-            }
-            well_index++;
-        }
-
-        // Set up reference depths that were defaulted. Count perfs.
-
-        const int num_wells = well_data.size();
-
-        int num_perfs = 0;
-        assert(grid.dimensions == 3);
-        for (int w = 0; w < num_wells; ++w) {
-            num_perfs += wellperf_data[w].size();
-            if (well_data[w].reference_bhp_depth < 0.0) {
-                // It was defaulted. Set reference depth to minimum perforation depth.
-                double min_depth = 1e100;
-                int num_wperfs = wellperf_data[w].size();
-                for (int perf = 0; perf < num_wperfs; ++perf) {
-                    double depth = grid.cell_centroids[3*wellperf_data[w][perf].cell + 2];
-                    min_depth = std::min(min_depth, depth);
-                }
-                well_data[w].reference_bhp_depth = min_depth;
-            }
-        }
-
-        // Create the well data structures.
-        w_ = create_wells(phaseUsage.num_phases, num_wells, num_perfs);
-        if (!w_) {
-            OPM_THROW(std::runtime_error, "Failed creating Wells struct.");
-        }
-
-
-        // Add wells.
-        for (int w = 0; w < num_wells; ++w) {
-            const int w_num_perf = wellperf_data[w].size();
-            std::vector<int> perf_cells(w_num_perf);
-            std::vector<double> perf_prodind(w_num_perf);
-            for (int perf = 0; perf < w_num_perf; ++perf) {
-                perf_cells[perf] = wellperf_data[w][perf].cell;
-                perf_prodind[perf] = wellperf_data[w][perf].well_index;
-            }
-            const double* comp_frac = NULL;
-            // We initialize all wells with a null component fraction,
-            // and must (for injection wells) overwrite it later.
-            int ok = add_well(well_data[w].type, well_data[w].reference_bhp_depth, w_num_perf,
-                              comp_frac, &perf_cells[0], &perf_prodind[0], well_names[w].c_str(), w_);
-            if (!ok) {
-                OPM_THROW(std::runtime_error, "Failed adding well " << well_names[w] << " to Wells data structure.");
-            }
-        }
-
-    }
 
     void WellsManager::setupWellControls(std::vector<WellConstPtr>& wells, size_t timeStep,
                                          std::vector<std::string>& well_names, const PhaseUsage& phaseUsage) {
@@ -558,7 +342,7 @@ namespace Opm
                                 
                 clear_well_controls(well_index, w_);
                 if (injectionProperties.hasInjectionControl(WellInjector::RATE)) {
-                    control_pos[InjectionControl::RATE] = well_controls_get_num(w_->ctrls[well_index]);
+                    control_pos[WellsManagerDetail::InjectionControl::RATE] = well_controls_get_num(w_->ctrls[well_index]);
                     double distr[3] = { 0.0, 0.0, 0.0 };
                     WellInjector::TypeEnum injectorType = injectionProperties.injectorType;
 
@@ -578,7 +362,7 @@ namespace Opm
                 }
 
                 if (ok && injectionProperties.hasInjectionControl(WellInjector::RESV)) {
-                    control_pos[InjectionControl::RESV] = well_controls_get_num(w_->ctrls[well_index]);
+                    control_pos[WellsManagerDetail::InjectionControl::RESV] = well_controls_get_num(w_->ctrls[well_index]);
                     double distr[3] = { 0.0, 0.0, 0.0 };
                     WellInjector::TypeEnum injectorType = injectionProperties.injectorType;
 
@@ -598,7 +382,8 @@ namespace Opm
                 }
 
                 if (ok && injectionProperties.hasInjectionControl(WellInjector::BHP)) {
-                    control_pos[InjectionControl::BHP] = well_controls_get_num(w_->ctrls[well_index]);
+                    control_pos[WellsManagerDetail::InjectionControl::BHP] = well_controls_get_num(w_->ctrls[well_index]);
+                    control_pos[WellsManagerDetail::InjectionControl::BHP] = well_controls_get_num(w_->ctrls[well_index]);
                     ok = append_well_controls(BHP,
                                               injectionProperties.BHPLimit, 
                                               NULL,
@@ -617,9 +402,9 @@ namespace Opm
 
 
                 {
-                    InjectionControl::Mode mode = InjectionControl::mode( injectionProperties.controlMode );
+                    WellsManagerDetail::InjectionControl::Mode mode = WellsManagerDetail::InjectionControl::mode( injectionProperties.controlMode );
                     int cpos = control_pos[mode];
-                    if (cpos == -1 && mode != InjectionControl::GRUP) {
+                    if (cpos == -1 && mode != WellsManagerDetail::InjectionControl::GRUP) {
                         OPM_THROW(std::runtime_error, "Control not specified in well " << well_names[well_index]);
                     }
 
@@ -670,7 +455,7 @@ namespace Opm
                         OPM_THROW(std::runtime_error, "Oil phase not active and ORAT control specified.");
                     }
 
-                    control_pos[ProductionControl::ORAT] = well_controls_get_num(w_->ctrls[well_index]);
+                    control_pos[WellsManagerDetail::ProductionControl::ORAT] = well_controls_get_num(w_->ctrls[well_index]);
                     double distr[3] = { 0.0, 0.0, 0.0 };
                     distr[phaseUsage.phase_pos[BlackoilPhases::Liquid]] = 1.0;
                     ok = append_well_controls(SURFACE_RATE,
@@ -684,7 +469,7 @@ namespace Opm
                     if (!phaseUsage.phase_used[BlackoilPhases::Aqua]) {
                         OPM_THROW(std::runtime_error, "Water phase not active and WRAT control specified.");
                     }
-                    control_pos[ProductionControl::WRAT] = well_controls_get_num(w_->ctrls[well_index]);
+                    control_pos[WellsManagerDetail::ProductionControl::WRAT] = well_controls_get_num(w_->ctrls[well_index]);
                     double distr[3] = { 0.0, 0.0, 0.0 };
                     distr[phaseUsage.phase_pos[BlackoilPhases::Aqua]] = 1.0;
                     ok = append_well_controls(SURFACE_RATE,
@@ -698,7 +483,7 @@ namespace Opm
                     if (!phaseUsage.phase_used[BlackoilPhases::Vapour]) {
                         OPM_THROW(std::runtime_error, "Gas phase not active and GRAT control specified.");
                     }
-                    control_pos[ProductionControl::GRAT] = well_controls_get_num(w_->ctrls[well_index]);
+                    control_pos[WellsManagerDetail::ProductionControl::GRAT] = well_controls_get_num(w_->ctrls[well_index]);
                     double distr[3] = { 0.0, 0.0, 0.0 };
                     distr[phaseUsage.phase_pos[BlackoilPhases::Vapour]] = 1.0;
                     ok = append_well_controls(SURFACE_RATE,
@@ -715,7 +500,7 @@ namespace Opm
                     if (!phaseUsage.phase_used[BlackoilPhases::Liquid]) {
                         OPM_THROW(std::runtime_error, "Oil phase not active and LRAT control specified.");
                     }
-                    control_pos[ProductionControl::LRAT] = well_controls_get_num(w_->ctrls[well_index]);
+                    control_pos[WellsManagerDetail::ProductionControl::LRAT] = well_controls_get_num(w_->ctrls[well_index]);
                     double distr[3] = { 0.0, 0.0, 0.0 };
                     distr[phaseUsage.phase_pos[BlackoilPhases::Aqua]] = 1.0;
                     distr[phaseUsage.phase_pos[BlackoilPhases::Liquid]] = 1.0;
@@ -727,7 +512,7 @@ namespace Opm
                 }
 
                 if (ok && productionProperties.hasProductionControl(WellProducer::RESV)) {
-                    control_pos[ProductionControl::RESV] = well_controls_get_num(w_->ctrls[well_index]);
+                    control_pos[WellsManagerDetail::ProductionControl::RESV] = well_controls_get_num(w_->ctrls[well_index]);
                     double distr[3] = { 1.0, 1.0, 1.0 };
                     ok = append_well_controls(RESERVOIR_RATE,
                                               -productionProperties.ResVRate , 
@@ -737,7 +522,7 @@ namespace Opm
                 }
 
                 if (ok && productionProperties.hasProductionControl(WellProducer::BHP)) {
-                    control_pos[ProductionControl::BHP] = well_controls_get_num(w_->ctrls[well_index]);
+                    control_pos[WellsManagerDetail::ProductionControl::BHP] = well_controls_get_num(w_->ctrls[well_index]);
                     ok = append_well_controls(BHP,
                                               productionProperties.BHPLimit , 
                                               NULL,
@@ -753,11 +538,15 @@ namespace Opm
                     OPM_THROW(std::runtime_error, "Failure occured appending controls for well " << well_names[well_index]);
                 }
 
-                ProductionControl::Mode mode = ProductionControl::mode(productionProperties.controlMode);
+                WellsManagerDetail::ProductionControl::Mode mode = WellsManagerDetail::ProductionControl::mode(productionProperties.controlMode);
                 int cpos = control_pos[mode];
+                if (cpos == -1 && mode != WellsManagerDetail::ProductionControl::GRUP) {
+                    OPM_THROW(std::runtime_error, "Control mode type " << mode << " not present in well " << well_names[well_index]);
+                }
+                // If it's shut, we complement the cpos
                 if (well->getStatus(timeStep) == WellCommon::SHUT) {
                     well_controls_shut_well( w_->ctrls[well_index] );
-                } else if (cpos == -1 && mode != ProductionControl::GRUP) {
+                } else if (cpos == -1 && mode != WellsManagerDetail::ProductionControl::GRUP) {
                     OPM_THROW(std::runtime_error, "Control mode type " << mode << " not present in well " << well_names[well_index]);
                 }
                 set_current_control(well_index, cpos, w_);

opm/core/wells/WellsManager.hpp

@@ -72,11 +72,22 @@ namespace Opm
         /// The permeability argument may be zero if the input contain
         /// well productivity indices, otherwise it must be given in
         /// order to approximate these by the Peaceman formula.
+        template<class CC, class F2C, class FC>
+        WellsManager(const Opm::EclipseStateConstPtr eclipseState,
+                     const size_t timeStep,
+                     int num_cells,
+                     const int* global_cell,
+                     const int* cart_dims,
+                     int dimensions,
+                     CC begin_cell_centroids,
+                     const F2C& f2c,
+                     FC begin_face_centroids,
+                     const double* permeability);
+
         WellsManager(const Opm::EclipseStateConstPtr eclipseState,
                      const size_t timeStep,
                      const UnstructuredGrid& grid,
                      const double* permeability);
-
         /// Destructor.
         ~WellsManager();
 
@@ -129,15 +140,30 @@ namespace Opm
 
 
     private:
+        template<class CC, class C2F, class FC>
+        void init(const Opm::EclipseStateConstPtr eclipseState,
+                  const size_t timeStep,
+                  int num_cells,
+                  const int* global_cell,
+                  const int* cart_dims,
+                  int dimensions,
+                  CC begin_cell_centroids,
+                  const C2F& cell_to_faces,
+                  FC begin_face_centroids,
+                  const double* permeability);
         // Disable copying and assignment.
         WellsManager(const WellsManager& other);
         WellsManager& operator=(const WellsManager& other);
-        static void setupCompressedToCartesian(const UnstructuredGrid& grid, std::map<int,int>& cartesian_to_compressed );
+        static void setupCompressedToCartesian(const int* global_cell, int number_of_cells, std::map<int,int>& cartesian_to_compressed );
         void setupWellControls(std::vector<WellConstPtr>& wells, size_t timeStep,
                                std::vector<std::string>& well_names, const PhaseUsage& phaseUsage);
-
+        template<class C2F, class CC, class FC>
         void createWellsFromSpecs( std::vector<WellConstPtr>& wells, size_t timeStep,
-                                   const UnstructuredGrid& grid,
+                                   const C2F& cell_to_faces, 
+                                   const int* cart_dims,
+                                   FC begin_face_centroids, 
+                                   CC begin_cell_centroids,
+                                   int dimensions,
                                    std::vector<std::string>& well_names,
                                    std::vector<WellData>& well_data,
                                    std::map<std::string, int> & well_names_to_index,
@@ -156,5 +182,5 @@ namespace Opm
 
 } // namespace Opm
 
-
+#include "WellsManager_impl.hpp"
 #endif // OPM_WELLSMANAGER_HEADER_INCLUDED

opm/core/wells/WellsManager_impl.hpp

@@ -0,0 +1,314 @@
+#include <opm/core/utility/Units.hpp>
+#include <opm/core/grid/GridHelpers.hpp>
+
+#include <array>
+
+namespace WellsManagerDetail
+{
+
+
+    namespace ProductionControl
+    {
+        enum Mode { ORAT, WRAT, GRAT,
+                    LRAT, CRAT, RESV,
+                    BHP , THP , GRUP };
+    /*
+        namespace Details {
+            std::map<std::string, Mode>
+            init_mode_map();
+        } // namespace Details
+    */
+        Mode mode(const std::string& control);
+
+
+        Mode mode(Opm::WellProducer::ControlModeEnum controlMode);
+    } // namespace ProductionControl
+
+
+    namespace InjectionControl
+    {
+        enum Mode { RATE, RESV, BHP,
+                    THP, GRUP };
+    /*
+        namespace Details {
+            std::map<std::string, Mode>
+            init_mode_map();
+        } // namespace Details
+    */
+        Mode mode(const std::string& control);
+
+        Mode mode(Opm::WellInjector::ControlModeEnum controlMode);
+
+    } // namespace InjectionControl
+
+double computeWellIndex(const double radius,
+                        const std::array<double, 3>& cubical,
+                        const double* cell_permeability,
+                        const double skin_factor);
+
+template<class C2F, class FC>
+std::array<double, 3> getCubeDim(const C2F& c2f, FC begin_face_centroids, int dimensions, 
+                                 int cell)
+{
+    using namespace std;
+    std::array<double, 3> cube;
+    //typedef Opm::UgGridHelpers::Cell2FacesTraits<UnstructuredGrid>::Type Cell2Faces;
+    //Cell2Faces c2f=Opm::UgGridHelpers::cell2Faces(grid);
+    typedef typename C2F::row_type FaceRow;
+    FaceRow faces=c2f[cell];
+    typename FaceRow::const_iterator face=faces.begin();
+    int num_local_faces = faces.end()-face;
+    vector<double> x(num_local_faces);
+    vector<double> y(num_local_faces);
+    vector<double> z(num_local_faces);
+    for (int lf=0; lf<num_local_faces; ++ lf, ++face) {
+        FC centroid = Opm::UgGridHelpers::increment(begin_face_centroids, *face, dimensions);
+        x[lf] = Opm::UgGridHelpers::getCoordinate(centroid, 0);
+        y[lf] = Opm::UgGridHelpers::getCoordinate(centroid, 1);
+        z[lf] = Opm::UgGridHelpers::getCoordinate(centroid, 2);
+    }
+    cube[0] = *max_element(x.begin(), x.end()) - *min_element(x.begin(), x.end());
+    cube[1] = *max_element(y.begin(), y.end()) - *min_element(y.begin(), y.end());
+    cube[2] = *max_element(z.begin(), z.end()) - *min_element(z.begin(), z.end());
+    return cube;
+}
+} // end namespace
+
+namespace Opm
+{
+template<class C2F, class CC, class FC>
+void WellsManager::createWellsFromSpecs(std::vector<WellConstPtr>& wells, size_t timeStep,
+                                        const C2F& c2f, 
+                                        const int* cart_dims,
+                                        FC begin_face_centroids, 
+                                        CC begin_cell_centroids,
+                                        int dimensions,
+                                        std::vector<std::string>& well_names,
+                                        std::vector<WellData>& well_data,
+                                        std::map<std::string, int>& well_names_to_index,
+                                        const PhaseUsage& phaseUsage,
+                                        std::map<int,int> cartesian_to_compressed,
+                                        const double* permeability)
+{
+    std::vector<std::vector<PerfData> > wellperf_data;
+    wellperf_data.resize(wells.size());
+
+    int well_index = 0;
+    for (auto wellIter= wells.begin(); wellIter != wells.end(); ++wellIter) {
+        WellConstPtr well = (*wellIter);
+        {   // WELSPECS handling
+            well_names_to_index[well->name()] = well_index;
+            well_names.push_back(well->name());
+            {
+                WellData wd;
+                // If negative (defaulted), set refdepth to a marker
+                // value, will be changed after getting perforation
+                // data to the centroid of the cell of the top well
+                // perforation.
+                wd.reference_bhp_depth = (well->getRefDepth() < 0.0) ? -1e100 : well->getRefDepth();
+                wd.welspecsline = -1;
+                if (well->isInjector( timeStep ))
+                    wd.type = INJECTOR;
+                else
+                    wd.type = PRODUCER;
+                well_data.push_back(wd);
+            }
+        }
+
+        {   // COMPDAT handling
+            CompletionSetConstPtr completionSet = well->getCompletions(timeStep);
+            for (size_t c=0; c<completionSet->size(); c++) {
+                CompletionConstPtr completion = completionSet->get(c);
+                int i = completion->getI();
+                int j = completion->getJ();
+                int k = completion->getK();
+
+                const int* cpgdim = cart_dims;
+                int cart_grid_indx = i + cpgdim[0]*(j + cpgdim[1]*k);
+                std::map<int, int>::const_iterator cgit = cartesian_to_compressed.find(cart_grid_indx);
+                if (cgit == cartesian_to_compressed.end()) {
+                        OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << i << ' ' << j << ' '
+                                  << k << " not found in grid (well = " << well->name() << ')');
+                }
+                int cell = cgit->second;
+                PerfData pd;
+                pd.cell = cell;
+                if (completion->getCF() > 0.0) {
+                    pd.well_index = completion->getCF();
+                } else {
+                    double radius = 0.5*completion->getDiameter();
+                    if (radius <= 0.0) {
+                        radius = 0.5*unit::feet;
+                        OPM_MESSAGE("**** Warning: Well bore internal radius set to " << radius);
+                    }
+                    std::array<double, 3> cubical = WellsManagerDetail::getCubeDim(c2f, begin_face_centroids,
+                                                                                   dimensions, cell);
+                    const double* cell_perm = &permeability[dimensions*dimensions*cell];
+                    pd.well_index = WellsManagerDetail::computeWellIndex(radius, cubical, cell_perm, completion->getSkinFactor());
+                }
+                wellperf_data[well_index].push_back(pd);
+            }
+        }
+        well_index++;
+    }
+
+    // Set up reference depths that were defaulted. Count perfs.
+
+    const int num_wells = well_data.size();
+
+    int num_perfs = 0;
+    assert(dimensions == 3);
+    for (int w = 0; w < num_wells; ++w) {
+        num_perfs += wellperf_data[w].size();
+        if (well_data[w].reference_bhp_depth < 0.0) {
+            // It was defaulted. Set reference depth to minimum perforation depth.
+            double min_depth = 1e100;
+            int num_wperfs = wellperf_data[w].size();
+            for (int perf = 0; perf < num_wperfs; ++perf) {
+                double depth = UgGridHelpers
+                    ::getCoordinate(UgGridHelpers::increment(begin_cell_centroids,
+                                                             wellperf_data[w][perf].cell,
+                                                             dimensions),
+                                    2);
+                min_depth = std::min(min_depth, depth);
+            }
+            well_data[w].reference_bhp_depth = min_depth;
+        }
+    }
+
+    // Create the well data structures.
+    w_ = create_wells(phaseUsage.num_phases, num_wells, num_perfs);
+    if (!w_) {
+        OPM_THROW(std::runtime_error, "Failed creating Wells struct.");
+    }
+
+
+    // Add wells.
+    for (int w = 0; w < num_wells; ++w) {
+        const int w_num_perf = wellperf_data[w].size();
+        std::vector<int> perf_cells(w_num_perf);
+        std::vector<double> perf_prodind(w_num_perf);
+        for (int perf = 0; perf < w_num_perf; ++perf) {
+            perf_cells[perf] = wellperf_data[w][perf].cell;
+            perf_prodind[perf] = wellperf_data[w][perf].well_index;
+        }
+        const double* comp_frac = NULL;
+        // We initialize all wells with a null component fraction,
+        // and must (for injection wells) overwrite it later.
+        int ok = add_well(well_data[w].type, well_data[w].reference_bhp_depth, w_num_perf,
+                          comp_frac, &perf_cells[0], &perf_prodind[0], well_names[w].c_str(), w_);
+        if (!ok) {
+            OPM_THROW(std::runtime_error, "Failed adding well " << well_names[w] << " to Wells data structure.");
+        }
+    }
+}
+
+template<class CC, class C2F, class FC>
+WellsManager::WellsManager(const Opm::EclipseStateConstPtr eclipseState,
+                           const size_t timeStep,
+                           int number_of_cells,
+                           const int* global_cell,
+                           const int* cart_dims,
+                           int dimensions,
+                           CC begin_cell_centroids,
+                           const C2F& cell_to_faces,
+                           FC begin_face_centroids,
+                           const double* permeability)
+    : w_(0)
+{
+    init(eclipseState, timeStep, number_of_cells, global_cell, cart_dims, dimensions,
+         begin_cell_centroids, cell_to_faces, begin_face_centroids, permeability);
+}
+
+/// Construct wells from deck.
+template<class CC, class C2F, class FC>
+void WellsManager::init(const Opm::EclipseStateConstPtr eclipseState,
+                        const size_t timeStep,
+                        int number_of_cells,
+                        const int* global_cell,
+                        const int* cart_dims,
+                        int dimensions,
+                        CC begin_cell_centroids,
+                        const C2F& cell_to_faces,
+                        FC begin_face_centroids,
+                        const double* permeability)
+{
+    if (dimensions != 3) {
+        OPM_THROW(std::runtime_error, "We cannot initialize wells from a deck unless the corresponding grid is 3-dimensional.");
+    }
+
+    if (eclipseState->getSchedule()->numWells() == 0) {
+        OPM_MESSAGE("No wells specified in Schedule section, initializing no wells");
+        return;
+    }
+
+    std::map<int,int> cartesian_to_compressed;
+    setupCompressedToCartesian(global_cell, 
+                               number_of_cells, cartesian_to_compressed);
+
+    // Obtain phase usage data.
+    PhaseUsage pu = phaseUsageFromDeck(eclipseState);
+
+    // These data structures will be filled in this constructor,
+    // then used to initialize the Wells struct.
+    std::vector<std::string> well_names;
+    std::vector<WellData> well_data;
+
+
+    // For easy lookup:
+    std::map<std::string, int> well_names_to_index;
+
+    ScheduleConstPtr schedule = eclipseState->getSchedule();
+    std::vector<WellConstPtr> wells = schedule->getWells(timeStep);
+
+    well_names.reserve(wells.size());
+    well_data.reserve(wells.size());
+
+    createWellsFromSpecs(wells, timeStep, cell_to_faces,
+                         cart_dims,
+                         begin_face_centroids,
+                         begin_cell_centroids,
+                         dimensions,
+                         well_names, well_data, well_names_to_index, pu, cartesian_to_compressed, permeability);
+    setupWellControls(wells, timeStep, well_names, pu);
+
+    {
+        GroupTreeNodeConstPtr fieldNode = eclipseState->getSchedule()->getGroupTree(timeStep)->getNode("FIELD");
+        GroupConstPtr fieldGroup = eclipseState->getSchedule()->getGroup(fieldNode->name());
+        well_collection_.addField(fieldGroup, timeStep, pu);
+        addChildGroups(fieldNode, eclipseState->getSchedule(), timeStep, pu);
+    }
+
+    for (auto wellIter = wells.begin(); wellIter != wells.end(); ++wellIter ) {
+        well_collection_.addWell((*wellIter), timeStep, pu);
+    }
+
+    well_collection_.setWellsPointer(w_);
+    well_collection_.applyGroupControls();
+
+
+    setupGuideRates(wells, timeStep, well_data, well_names_to_index);
+
+    // Debug output.
+#define EXTRA_OUTPUT
+#ifdef EXTRA_OUTPUT
+    /*
+      std::cout << "\t WELL DATA" << std::endl;
+      for(int i = 0; i< num_wells; ++i) {
+      std::cout << i << ": " << well_data[i].type << "  "
+      << well_data[i].control << "  " << well_data[i].target
+      << std::endl;
+      }
+
+      std::cout << "\n\t PERF DATA" << std::endl;
+      for(int i=0; i< int(wellperf_data.size()); ++i) {
+      for(int j=0; j< int(wellperf_data[i].size()); ++j) {
+      std::cout << i << ": " << wellperf_data[i][j].cell << "  "
+      << wellperf_data[i][j].well_index << std::endl;
+      }
+      }
+    */
+#endif
+}
+
+} // end namespace Opm