Move 'common' directory into 'dune'.

dune/common/Average.hpp

@@ -0,0 +1,98 @@
+//===========================================================================
+//
+// File: Average.hpp<2>
+//
+// Created: Wed Jul  1 12:25:57 2009
+//
+// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
+//            Bård Skaflestad     <bard.skaflestad@sintef.no>
+//
+// $Date$
+//
+// $Revision$
+//
+//===========================================================================
+
+/*
+  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+  Copyright 2009, 2010 Statoil ASA.
+
+  This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+  OpenRS is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OpenRS is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPENRS_AVERAGE_HEADER
+#define OPENRS_AVERAGE_HEADER
+
+
+#include <boost/static_assert.hpp>
+#include <tr1/type_traits>
+#include <cmath>
+
+namespace Dune {
+    namespace utils {
+
+
+
+	/// Computes the arithmetic average:
+	/// a_A(t_1, t_2) = (t_1 + t_2)/2.
+	template <typename T, typename Tresult>
+	Tresult arithmeticAverage(const T& t1, const T& t2)
+	{
+	    // To avoid some user errors, we disallow taking averages of
+	    // integral values.
+	    BOOST_STATIC_ASSERT(std::tr1::is_integral<T>::value == false);
+	    Tresult retval(t1);
+	    retval += t2;
+	    retval *= 0.5;
+	    return retval;
+	}
+
+
+
+	/// Computes the geometric average:
+	/// a_G(t_1, t_2) = \sqrt{t_1 t_2}.
+	/// Since we use std::sqrt(), this function will only
+	/// compile when T is convertible to and from double.
+	template <typename T>
+	T geometricAverage(const T& t1, const T& t2)
+	{
+	    // To avoid some user errors, we disallow taking averages of
+	    // integral values.
+	    BOOST_STATIC_ASSERT(std::tr1::is_integral<T>::value == false);
+	    return std::sqrt(t1*t2);
+	}
+
+
+
+	/// Computes the harmonic average:
+	/// a_H(t_1, t_2) = \frac{2}{1/t_1 + 1/t_2} = \frac{2 t_1 t_2}{t_1 + t_2}.
+	template <typename T>
+	T harmonicAverage(const T& t1, const T& t2)
+	{
+	    // To avoid some user errors, we disallow taking averages of
+	    // integral values.
+	    BOOST_STATIC_ASSERT(std::tr1::is_integral<T>::value == false);
+	    return (2*t1*t2)/(t1 + t2);
+	}
+
+
+
+    } // namespace utils
+} // namespace Dune
+
+
+
+#endif // OPENRS_AVERAGE_HEADER

dune/common/CornerpointChopper.hpp

@@ -0,0 +1,399 @@
+/*
+  Copyright 2010 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_CORNERPOINTCHOPPER_HEADER_INCLUDED
+#define OPM_CORNERPOINTCHOPPER_HEADER_INCLUDED
+
+#include <dune/common/EclipseGridParser.hpp>
+#include <dune/common/param/ParameterGroup.hpp>
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <vector>
+#include <stdexcept>
+
+namespace Dune
+{
+
+    class CornerPointChopper
+    {
+    public:
+        CornerPointChopper(const std::string& file)
+            : parser_(file, false)
+        {
+            for (int dd = 0; dd < 3; ++dd) {
+                dims_[dd] = parser_.getSPECGRID().dimensions[dd];
+            }
+
+            int layersz = 8*dims_[0]*dims_[1];
+            const std::vector<double>& ZCORN = parser_.getFloatingPointValue("ZCORN");
+            botmax_ = *std::max_element(ZCORN.begin(), ZCORN.begin() + layersz/2);
+            topmin_ = *std::min_element(ZCORN.begin() + dims_[2]*layersz - layersz/2,
+                                        ZCORN.begin() + dims_[2]*layersz);
+
+            abszmax_ = *std::max_element(ZCORN.begin(), ZCORN.end());
+            abszmin_ = *std::min_element(ZCORN.begin(), ZCORN.end());
+
+            std::cout << "Parsed grdecl file with dimensions ("
+                      << dims_[0] << ", " << dims_[1] << ", " << dims_[2] << ")" << std::endl;
+        }
+
+
+
+
+        const int* dimensions() const
+        {
+            return dims_;
+        }
+
+
+
+
+        const int* newDimensions() const
+        {
+            return new_dims_;
+        }
+
+
+
+
+        const std::pair<double, double> zLimits() const
+        {
+            return std::make_pair(botmax_, topmin_);
+        }
+
+        const std::pair<double, double> abszLimits() const
+	{
+            return std::make_pair(abszmin_, abszmax_);
+        }
+
+
+	void verifyInscribedShoebox(int imin, int ilen, int imax,
+				    int jmin, int jlen, int jmax,
+				    double zmin, double zlen, double zmax) 
+	{
+	    if (imin < 0) {
+		std::cerr << "Error! imin < 0 (imin = " << imin << ")\n";
+		throw std::runtime_error("Inconsistent user input.");
+	    }
+	    if (ilen > dims_[0]) {
+		std::cerr << "Error! ilen larger than grid (ilen = " << ilen <<")\n";
+		throw std::runtime_error("Inconsistent user input.");
+	    }
+	    if (imax > dims_[0]) {
+		std::cerr << "Error! imax larger than input grid (imax = " << imax << ")\n";
+		throw std::runtime_error("Inconsistent user input.");
+	    }
+	    if (jmin < 0) {
+		std::cerr << "Error! jmin < 0 (jmin = " << jmin << ")\n";
+		throw std::runtime_error("Inconsistent user input.");
+	    }
+	    if (jlen > dims_[1]) {
+		std::cerr << "Error! jlen larger than grid (jlen = " << jlen <<")\n";
+		throw std::runtime_error("Inconsistent user input.");
+	    }
+	    if (jmax > dims_[1]) {
+		std::cerr << "Error! jmax larger than input grid (jmax = " << jmax << ")\n";
+		throw std::runtime_error("Inconsistent user input.");
+	    }
+	    if (zmin < abszmin_) {
+		std::cerr << "Error! zmin ("<< zmin << ") less than minimum ZCORN value ("<< abszmin_ << ")\n";
+		throw std::runtime_error("Inconsistent user input.");
+	    }
+	    if (zmax > abszmax_) {
+		std::cerr << "Error! zmax ("<< zmax << ") larger than maximal ZCORN value ("<< abszmax_ << ")\n";
+		throw std::runtime_error("Inconsistent user input.");
+	    }
+	    if (zlen > (abszmax_ - abszmin_)) {
+		std::cerr << "Error! zlen ("<< zlen <<") larger than maximal ZCORN (" << abszmax_ << ") minus minimal ZCORN ("<< abszmin_ <<")\n";
+		throw std::runtime_error("Inconsistent user input.");
+	    }
+	}
+
+        void chop(int imin, int imax, int jmin, int jmax, double zmin, double zmax, bool resettoorigin=true)
+        {
+            new_dims_[0] = imax - imin;
+            new_dims_[1] = jmax - jmin;
+
+            // Filter the coord field
+            const std::vector<double>& COORD = parser_.getFloatingPointValue("COORD");
+            int num_coord = COORD.size();
+            if (num_coord != 6*(dims_[0] + 1)*(dims_[1] + 1)) {
+                std::cerr << "Error! COORD size (" << COORD.size() << ") not consistent with SPECGRID\n";
+                throw std::runtime_error("Inconsistent COORD and SPECGRID.");
+            }
+            int num_new_coord = 6*(new_dims_[0] + 1)*(new_dims_[1] + 1);
+	    double x_correction = COORD[6*((dims_[0] + 1)*jmin + imin)];
+	    double y_correction = COORD[6*((dims_[0] + 1)*jmin + imin) + 1];
+            new_COORD_.resize(num_new_coord, 1e100);
+            for (int j = jmin; j < jmax + 1; ++j) {
+                for (int i = imin; i < imax + 1; ++i) {
+                    int pos = (dims_[0] + 1)*j + i;
+                    int new_pos = (new_dims_[0] + 1)*(j-jmin) + (i-imin);
+		    // Copy all 6 coordinates for a pillar.
+                    std::copy(COORD.begin() + 6*pos, COORD.begin() + 6*(pos + 1), new_COORD_.begin() + 6*new_pos);
+		    if (resettoorigin) {
+		        // Substract lowest x value from all X-coords, similarly for y, and truncate in z-direction
+		      new_COORD_[6*new_pos]     -= x_correction;
+		      new_COORD_[6*new_pos + 1] -= y_correction;
+		      new_COORD_[6*new_pos + 2]  = 0;
+		      new_COORD_[6*new_pos + 3] -= x_correction;
+		      new_COORD_[6*new_pos + 4] -= y_correction;
+		      new_COORD_[6*new_pos + 5]  = zmax-zmin;
+		    }
+                }
+            }
+
+            // Get the z limits, check if they must be changed to make a shoe-box.
+            // This means that zmin must be greater than or equal to the highest
+            // coordinate of the bottom surface, while zmax must be less than or
+            // equal to the lowest coordinate of the top surface.
+            int layersz = 8*dims_[0]*dims_[1];
+            const std::vector<double>& ZCORN = parser_.getFloatingPointValue("ZCORN");
+            int num_zcorn = ZCORN.size();
+            if (num_zcorn != layersz*dims_[2]) {
+                std::cerr << "Error! ZCORN size (" << ZCORN.size() << ") not consistent with SPECGRID\n";
+                throw std::runtime_error("Inconsistent ZCORN and SPECGRID.");
+            }
+
+            zmin = std::max(zmin, botmax_);
+            zmax = std::min(zmax, topmin_);
+            if (zmin >= zmax) {
+                std::cerr << "Error: zmin >= zmax (zmin = " << zmin << ", zmax = " << zmax << ")\n";
+                throw std::runtime_error("zmin >= zmax");
+            }
+            std::cout << "Chopping subsample,  i: (" << imin << "--" << imax << ")  j: (" << jmin << "--" << jmax << ")   z: (" << zmin << "--" << zmax << ")" <<  std::endl;
+
+            // We must find the maximum and minimum k value for the given z limits.
+            // First, find the first layer with a z-coordinate strictly above zmin.
+            int kmin = -1;
+            for (int k = 0; k < dims_[2]; ++k) {
+                double layer_max = *std::max_element(ZCORN.begin() + k*layersz, ZCORN.begin() + (k + 1)*layersz);
+                if (layer_max > zmin) {
+                    kmin = k;
+                    break;
+                }
+            }
+            // Then, find the last layer with a z-coordinate strictly below zmax.
+            int kmax = -1;
+            for (int k = dims_[2]; k > 0; --k) {
+                double layer_min = *std::min_element(ZCORN.begin() + (k - 1)*layersz, ZCORN.begin() + k*layersz);
+                if (layer_min < zmax) {
+                    kmax = k;
+                    break;
+                }
+            }
+            new_dims_[2] = kmax - kmin;
+
+            // Filter the ZCORN field, build mapping from new to old cells.
+	    double z_origin_correction = 0.0;
+	    if (resettoorigin) {
+		z_origin_correction = zmin;
+	    }
+            new_ZCORN_.resize(8*new_dims_[0]*new_dims_[1]*new_dims_[2], 1e100);
+            new_to_old_cell_.resize(new_dims_[0]*new_dims_[1]*new_dims_[2], -1);
+            int cellcount = 0;
+            int delta[3] = { 1, 2*dims_[0], 4*dims_[0]*dims_[1] };
+            int new_delta[3] = { 1, 2*new_dims_[0], 4*new_dims_[0]*new_dims_[1] };
+            for (int k = kmin; k < kmax; ++k) {
+                for (int j = jmin; j < jmax; ++j) {
+                    for (int i = imin; i < imax; ++i) {
+                        new_to_old_cell_[cellcount++] = dims_[0]*dims_[1]*k + dims_[0]*j + i;
+                        int old_ix = 2*(i*delta[0] + j*delta[1] + k*delta[2]);
+                        int new_ix = 2*((i-imin)*new_delta[0] + (j-jmin)*new_delta[1] + (k-kmin)*new_delta[2]);
+                        int old_indices[8] = { old_ix, old_ix + delta[0],
+                                               old_ix + delta[1], old_ix + delta[1] + delta[0],
+                                               old_ix + delta[2], old_ix + delta[2] + delta[0],
+                                               old_ix + delta[2] + delta[1], old_ix + delta[2] + delta[1] + delta[0] };
+                        int new_indices[8] = { new_ix, new_ix + new_delta[0],
+                                               new_ix + new_delta[1], new_ix + new_delta[1] + new_delta[0],
+                                               new_ix + new_delta[2], new_ix + new_delta[2] + new_delta[0],
+                                               new_ix + new_delta[2] + new_delta[1], new_ix + new_delta[2] + new_delta[1] + new_delta[0] };
+                        for (int cc = 0; cc < 8; ++cc) {
+                            new_ZCORN_[new_indices[cc]] = std::min(zmax, std::max(zmin, ZCORN[old_indices[cc]])) - z_origin_correction;
+                        }
+                    }
+                }
+            }
+            
+            filterIntegerField("ACTNUM", new_ACTNUM_);
+            filterDoubleField("PORO", new_PORO_);
+            filterDoubleField("PERMX", new_PERMX_);
+            filterDoubleField("PERMY", new_PERMY_);
+            filterDoubleField("PERMZ", new_PERMZ_);
+            filterIntegerField("SATNUM", new_SATNUM_);
+        }
+
+
+
+        /// Return a subparser with fields corresponding to the selected subset.
+        /// Note that the returned parser is NOT converted to SI, that must be done
+        /// by the user afterwards with the parser's convertToSI() method.
+        EclipseGridParser subparser()
+        {
+            if (parser_.hasField("FIELD") || parser_.hasField("LAB") || parser_.hasField("PVT-M")) {
+                THROW("CornerPointChopper::subparser() cannot handle any eclipse unit system other than METRIC.");
+            }
+
+            EclipseGridParser sp;
+            boost::shared_ptr<SPECGRID> sg(new SPECGRID);
+            for (int dd = 0; dd < 3; ++dd) {
+                sg->dimensions[dd] = new_dims_[dd];
+            }
+            sp.setSpecialField("SPECGRID", sg);
+            sp.setFloatingPointField("COORD", new_COORD_);
+            sp.setFloatingPointField("ZCORN", new_ZCORN_);
+            if (!new_ACTNUM_.empty()) sp.setIntegerField("ACTNUM", new_ACTNUM_);
+            if (!new_PORO_.empty()) sp.setFloatingPointField("PORO", new_PORO_);
+            if (!new_PERMX_.empty()) sp.setFloatingPointField("PERMX", new_PERMX_);
+            if (!new_PERMY_.empty()) sp.setFloatingPointField("PERMY", new_PERMY_);
+            if (!new_PERMZ_.empty()) sp.setFloatingPointField("PERMZ", new_PERMZ_);
+            if (!new_SATNUM_.empty()) sp.setIntegerField("SATNUM", new_SATNUM_);
+            sp.computeUnits(); // Always METRIC, since that is default.
+            return sp;
+        }
+
+
+
+
+        void writeGrdecl(const std::string& filename)
+        {
+            // Output new versions of SPECGRID, COORD, ZCORN, ACTNUM, PERMX, PORO, SATNUM.
+            std::ofstream out(filename.c_str());
+            if (!out) {
+                std::cerr << "Could not open file " << filename << "\n";
+                throw std::runtime_error("Could not open output file.");
+            }
+            out << "SPECGRID\n" << new_dims_[0] << ' ' << new_dims_[1] << ' ' << new_dims_[2]
+                << " 1 F\n/\n\n";
+            out << "COORD\n";
+            int num_new_coord = new_COORD_.size();
+            for (int i = 0; i < num_new_coord/6; ++i) {
+                for (int j = 0; j < 6; ++j) {
+                    out << "  " << new_COORD_[6*i + j];
+                }
+                out << '\n';
+            }
+            out << "/\n\n";
+            out << "ZCORN\n";
+            int num_new_zcorn = new_ZCORN_.size();
+            assert(num_new_zcorn%8 == 0);
+            for (int i = 0; i < num_new_zcorn/8; ++i) {
+                for (int j = 0; j < 8; ++j) {
+                    out << "  " << new_ZCORN_[8*i + j];
+                }
+                out << '\n';
+            }
+            out << "/\n\n";
+
+            outputField(out, new_ACTNUM_, "ACTNUM");
+            outputField(out, new_PORO_, "PORO");
+            outputField(out, new_PERMX_, "PERMX");
+            outputField(out, new_PERMY_, "PERMY");
+            outputField(out, new_PERMZ_, "PERMZ");
+            outputField(out, new_SATNUM_, "SATNUM");
+        }
+
+    private:
+        EclipseGridParser parser_;
+        double botmax_;
+        double topmin_;
+        double abszmin_;
+        double abszmax_;
+        std::vector<double> new_COORD_;
+        std::vector<double> new_ZCORN_;
+        std::vector<int> new_ACTNUM_;
+        std::vector<double> new_PORO_;
+        std::vector<double> new_PERMX_;
+        std::vector<double> new_PERMY_;
+        std::vector<double> new_PERMZ_;
+        std::vector<int> new_SATNUM_;
+        int dims_[3];
+        int new_dims_[3];
+        std::vector<int> new_to_old_cell_;
+
+
+        template <typename T>
+        void outputField(std::ostream& os,
+                         const std::vector<T>& field,
+                         const std::string& keyword)
+        {
+            if (field.empty()) return;
+            os << keyword << '\n';
+            int sz = field.size();
+            T last = std::numeric_limits<T>::max();
+            int repeats = 0;
+            for (int i = 0; i < sz; ++i) {
+                T val = field[i];
+                if (val == last) {
+                    ++repeats;
+                } else {
+                    if (repeats == 1) {
+                        os << last << '\n';
+                    } else if (repeats > 1) {
+                        os << repeats << '*' << last << '\n';
+                    }
+                    last = val;
+                    repeats = 1;
+                }
+            }
+            if (repeats == 1) {
+                os << last << '\n';
+            } else if (repeats > 1) {
+                os << repeats << '*' << last << '\n';
+            }
+            os << "/\n\n";
+        }
+
+
+        template <typename T>
+        void filterField(const std::vector<T>& field,
+                                std::vector<T>& output_field)
+        {
+            int sz = new_to_old_cell_.size();
+            output_field.resize(sz);
+            for (int i = 0; i < sz; ++i) {
+                output_field[i] = field[new_to_old_cell_[i]];
+            }
+        }
+
+        void filterDoubleField(const std::string& keyword, std::vector<double>& output_field)
+        {
+            if (parser_.hasField(keyword)) {
+                const std::vector<double>& field = parser_.getFloatingPointValue(keyword);
+                filterField(field, output_field);
+            }
+        }
+
+        void filterIntegerField(const std::string& keyword, std::vector<int>& output_field)
+        {
+            if (parser_.hasField(keyword)) {
+                const std::vector<int>& field = parser_.getIntegerValue(keyword);
+                filterField(field, output_field);
+            }
+        }
+
+    };
+
+}
+
+
+
+
+#endif // OPM_CORNERPOINTCHOPPER_HEADER_INCLUDED

dune/common/EclipseGridInspector.cpp

@@ -0,0 +1,335 @@
+//===========================================================================
+//
+// File: EclipseGridInspector.C
+//
+// Created: Mon Jun  2 12:17:51 2008
+//
+// Author: Atgeirr F Rasmussen <atgeirr@sintef.no>
+//
+// $Date$
+//
+// $Revision$
+//
+// Revision: $Id: EclipseGridInspector.C,v 1.2 2008/08/18 14:16:13 atgeirr Exp $
+//
+//===========================================================================
+
+/*
+  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+  Copyright 2009, 2010 Statoil ASA.
+
+  This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+  OpenRS is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OpenRS is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#if HAVE_CONFIG_H
+#include "config.h"
+#endif
+#include <stdexcept>
+#include <numeric>
+#include <cmath>
+#include <cfloat>
+#include <algorithm>
+#include "EclipseGridInspector.hpp"
+#include "EclipseGridParser.hpp"
+#include "SpecialEclipseFields.hpp"
+
+namespace Dune
+{
+
+EclipseGridInspector::EclipseGridInspector(const EclipseGridParser& parser)
+    : parser_(parser)
+{
+    std::vector<std::string> keywords;
+    keywords.push_back("COORD");
+    keywords.push_back("ZCORN");
+
+    if (!parser_.hasFields(keywords)) {
+	THROW("Needed field is missing in file");
+    }
+
+    if (parser_.hasField("SPECGRID")) {
+        const SPECGRID& sgr = parser.getSPECGRID();
+        logical_gridsize_[0] = sgr.dimensions[0];
+        logical_gridsize_[1] = sgr.dimensions[1];
+        logical_gridsize_[2] = sgr.dimensions[2];
+    } else if (parser_.hasField("DIMENS")) {
+        const std::vector<int>& dim = parser.getIntegerValue("DIMENS");
+        logical_gridsize_[0] = dim[0];
+        logical_gridsize_[1] = dim[1];
+        logical_gridsize_[2] = dim[2];
+    } else {
+        THROW("Found neither SPECGRID nor DIMENS in file. At least one is needed.");
+    }
+
+}
+
+/**
+   Return the dip slopes for the cell relative to xy-plane in x- and y- direction. 
+   Dip slope is average rise in positive x-direction over cell length in x-direction. 
+   Similarly for y.
+
+   Current implementation is for vertical pillars, but is not difficult to fix.
+
+   @returns a std::pair<double,double> with x-dip in first component and y-dip in second.
+*/  
+std::pair<double,double> EclipseGridInspector::cellDips(int i, int j, int k) const
+{
+    checkLogicalCoords(i, j, k);
+    const std::vector<double>& pillc = parser_.getFloatingPointValue("COORD");
+    int num_pillars = (logical_gridsize_[0] + 1)*(logical_gridsize_[1] + 1);
+        if (6*num_pillars != int(pillc.size())) {
+        throw std::runtime_error("Wrong size of COORD field.");
+    }
+    const std::vector<double>& z = parser_.getFloatingPointValue("ZCORN");
+    int num_cells = logical_gridsize_[0]*logical_gridsize_[1]*logical_gridsize_[2];
+    if (8*num_cells != int(z.size())) {
+        throw std::runtime_error("Wrong size of ZCORN field");
+    }
+
+    // Pick ZCORN-value for all 8 corners of the given cell
+    boost::array<double, 8> cellz = cellZvals(i, j, k);
+
+    // Compute rise in positive x-direction for all four edges (and then find mean)
+    // Current implementation is for regularly placed and vertical pillars!
+    int numxpill = logical_gridsize_[0] + 1;
+    int pix = i + j*numxpill;
+    double cell_xlength = pillc[6*(pix + 1)] - pillc[6*pix];
+    flush(std::cout);
+    double xrise[4] = { (cellz[1] - cellz[0])/cell_xlength,  // LLL -> HLL
+                        (cellz[3] - cellz[2])/cell_xlength,  // LHL -> HHL
+                        (cellz[5] - cellz[4])/cell_xlength,  // LLH -> HLH
+                        (cellz[7] - cellz[6])/cell_xlength}; // LHH -> HHH
+
+    double cell_ylength = pillc[6*(pix + numxpill) + 1] - pillc[6*pix + 1];
+    double yrise[4] = { (cellz[2] - cellz[0])/cell_ylength,  // LLL -> LHL
+                        (cellz[3] - cellz[1])/cell_ylength,  // HLL -> HHL
+                        (cellz[6] - cellz[4])/cell_ylength,  // LLH -> LHH
+                        (cellz[7] - cellz[5])/cell_ylength}; // HLH -> HHH
+
+
+    // Now ignore those edges that touch the global top or bottom surface
+    // of the entire grdecl model. This is to avoid bias, as these edges probably
+    // don't follow an overall dip for the model if it exists.
+    int x_edges = 4;
+    int y_edges = 4;
+    boost::array<double, 6> gridlimits = getGridLimits();
+    double zmin = gridlimits[4];
+    double zmax = gridlimits[5];
+    // LLL -> HLL
+    if ((cellz[1] == zmin) || (cellz[0] == zmin)) {
+	xrise[0] = 0; x_edges--;
+    }
+    // LHL -> HHL
+    if ((cellz[2] == zmin) || (cellz[3] == zmin)) {
+	xrise[1] = 0; x_edges--;
+    }
+    // LLH -> HLH
+    if ((cellz[4] == zmax) || (cellz[5] == zmax)) {
+	xrise[2] = 0; x_edges--;
+    }
+    // LHH -> HHH
+    if ((cellz[6] == zmax) || (cellz[7] == zmax)) {
+	xrise[3] = 0; x_edges--;
+    }
+    // LLL -> LHL
+    if ((cellz[0] == zmin) || (cellz[2] == zmin)) {
+	yrise[0] = 0; y_edges--;
+    }
+    // HLL -> HHL
+    if ((cellz[1] == zmin) || (cellz[3] == zmin)) {
+	yrise[1] = 0; y_edges--;
+    }
+    // LLH -> LHH
+    if ((cellz[6] == zmax) || (cellz[4] == zmax)) {
+	yrise[2] = 0; y_edges--;
+    }
+    // HLH -> HHH
+    if ((cellz[7] == zmax) || (cellz[5] == zmax)) {
+	yrise[3] = 0; y_edges--;
+    }
+
+    return std::make_pair( (xrise[0] + xrise[1] + xrise[2] + xrise[3])/x_edges,
+                          (yrise[0] + yrise[1] + yrise[2] + yrise[3])/y_edges);
+}
+/**
+  Wrapper for cellDips(i, j, k).
+*/
+std::pair<double,double> EclipseGridInspector::cellDips(int cell_idx) const
+{
+    boost::array<int, 3> idxs = cellIdxToLogicalCoords(cell_idx);
+    return cellDips(idxs[0], idxs[1], idxs[2]);
+}
+
+boost::array<int, 3> EclipseGridInspector::cellIdxToLogicalCoords(int cell_idx) const
+{
+    
+    int i,j,k; // Position of cell in cell hierarchy
+    int horIdx = (cell_idx+1) - int(std::floor(((double)(cell_idx+1))/((double)(logical_gridsize_[0]*logical_gridsize_[1]))))*logical_gridsize_[0]*logical_gridsize_[1]; // index in the corresponding horizon
+    if (horIdx == 0) { 
+        horIdx = logical_gridsize_[0]*logical_gridsize_[1]; 
+    }
+    i = horIdx - int(std::floor(((double)horIdx)/((double)logical_gridsize_[0])))*logical_gridsize_[0];
+    if (i == 0) { 
+        i = logical_gridsize_[0]; 
+    }
+    j = (horIdx-i)/logical_gridsize_[0]+1;
+    k = ((cell_idx+1)-logical_gridsize_[0]*(j-1)-1)/(logical_gridsize_[0]*logical_gridsize_[1])+1;
+
+    boost::array<int, 3> a = {{i-1, j-1, k-1}};
+    return a; //boost::array<int, 3> {{i-1, j-1, k-1}};
+}
+
+double EclipseGridInspector::cellVolumeVerticalPillars(int i, int j, int k) const
+{
+    // Checking parameters and obtaining values from parser.
+    checkLogicalCoords(i, j, k);
+    const std::vector<double>& pillc = parser_.getFloatingPointValue("COORD");
+    int num_pillars = (logical_gridsize_[0] + 1)*(logical_gridsize_[1] + 1);
+    if (6*num_pillars != int(pillc.size())) {
+	throw std::runtime_error("Wrong size of COORD field.");
+    }
+    const std::vector<double>& z = parser_.getFloatingPointValue("ZCORN");
+    int num_cells = logical_gridsize_[0]*logical_gridsize_[1]*logical_gridsize_[2];
+    if (8*num_cells != int(z.size())) {
+	throw std::runtime_error("Wrong size of ZCORN field");
+    }
+
+    // Computing the base area as half the 2d cross product of the diagonals.
+    int numxpill = logical_gridsize_[0] + 1;
+    int pix = i + j*numxpill;
+    double px[4] = { pillc[6*pix],
+		     pillc[6*(pix + 1)],
+		     pillc[6*(pix + numxpill)],
+		     pillc[6*(pix + numxpill + 1)] };
+    double py[4] = { pillc[6*pix + 1],
+		     pillc[6*(pix + 1) + 1],
+		     pillc[6*(pix + numxpill) + 1],
+		     pillc[6*(pix + numxpill + 1) + 1] };
+    double diag1[2] = { px[3] - px[0], py[3] - py[0] };
+    double diag2[2] = { px[2] - px[1], py[2] - py[1] };
+    double area = 0.5*(diag1[0]*diag2[1] - diag1[1]*diag2[0]);
+
+    // Computing the average of the z-differences along each pillar.
+    int delta[3] = { 1,
+		     2*logical_gridsize_[0],
+		     4*logical_gridsize_[0]*logical_gridsize_[1] };
+    int ix = 2*(i*delta[0] + j*delta[1] + k*delta[2]);
+    double cellz[8] = { z[ix], z[ix + delta[0]],
+			z[ix + delta[1]], z[ix + delta[1] + delta[0]],
+			z[ix + delta[2]], z[ix + delta[2] + delta[0]],
+			z[ix + delta[2] + delta[1]], z[ix + delta[2] + delta[1] + delta[0]] };
+    double diffz[4] = { cellz[4] - cellz[0],
+			cellz[5] - cellz[1],
+			cellz[6] - cellz[2],
+			cellz[7] - cellz[3] };
+    double averzdiff = 0.25*std::accumulate(diffz, diffz + 4, 0.0);
+    return averzdiff*area;
+}
+
+
+double EclipseGridInspector::cellVolumeVerticalPillars(int cell_idx) const
+{
+    boost::array<int, 3> idxs = cellIdxToLogicalCoords(cell_idx);
+    return cellVolumeVerticalPillars(idxs[0], idxs[1], idxs[2]);
+}
+
+void EclipseGridInspector::checkLogicalCoords(int i, int j, int k) const
+{
+    if (i < 0 || i >= logical_gridsize_[0])
+	throw std::runtime_error("First coordinate out of bounds");
+    if (j < 0 || j >= logical_gridsize_[1])
+	throw std::runtime_error("Second coordinate out of bounds");
+    if (k < 0 || k >= logical_gridsize_[2])
+	throw std::runtime_error("Third coordinate out of bounds");
+}
+
+
+boost::array<double, 6> EclipseGridInspector::getGridLimits() const
+{
+    if (! (parser_.hasField("COORD") && parser_.hasField("ZCORN") && parser_.hasField("SPECGRID")) ) {
+        throw std::runtime_error("EclipseGridInspector: Grid does not have SPECGRID, COORD, and ZCORN, can't find dimensions.");
+    }
+
+    std::vector<double> coord = parser_.getFloatingPointValue("COORD");
+    std::vector<double> zcorn = parser_.getFloatingPointValue("ZCORN");
+
+    double xmin = +DBL_MAX;
+    double xmax = -DBL_MAX;
+    double ymin = +DBL_MAX;
+    double ymax = -DBL_MAX;
+
+
+    int pillars = (logical_gridsize_[0]+1) * (logical_gridsize_[1]+1);
+
+    for (int pillarindex = 0; pillarindex < pillars; ++pillarindex) {
+        if        (coord[pillarindex * 6 + 0] > xmax)
+            xmax = coord[pillarindex * 6 + 0];
+        if        (coord[pillarindex * 6 + 0] < xmin)
+            xmin = coord[pillarindex * 6 + 0];
+        if        (coord[pillarindex * 6 + 1] > ymax)
+            ymax = coord[pillarindex * 6 + 1];
+        if        (coord[pillarindex * 6 + 1] < ymin)
+            ymin = coord[pillarindex * 6 + 1];
+        if        (coord[pillarindex * 6 + 3] > xmax)
+            xmax = coord[pillarindex * 6 + 3];
+        if        (coord[pillarindex * 6 + 3] < xmin)
+            xmin = coord[pillarindex * 6 + 3];
+        if        (coord[pillarindex * 6 + 4] > ymax)
+            ymax = coord[pillarindex * 6 + 4];
+        if        (coord[pillarindex * 6 + 4] < ymin)
+            ymin = coord[pillarindex * 6 + 4];
+    }
+
+    boost::array<double, 6> gridlimits = {{ xmin, xmax, ymin, ymax,
+                                            *min_element(zcorn.begin(), zcorn.end()),
+                                            *max_element(zcorn.begin(), zcorn.end()) }};
+    return gridlimits;
+}
+
+
+
+boost::array<int, 3> EclipseGridInspector::gridSize() const
+{
+    boost::array<int, 3> retval = {{ logical_gridsize_[0],
+				     logical_gridsize_[1],
+				     logical_gridsize_[2] }};
+    return retval;
+}
+
+
+boost::array<double, 8> EclipseGridInspector::cellZvals(int i, int j, int k) const
+{
+    // Get the zcorn field.
+    const std::vector<double>& z = parser_.getFloatingPointValue("ZCORN");
+    int num_cells = logical_gridsize_[0]*logical_gridsize_[1]*logical_gridsize_[2];
+    if (8*num_cells != int(z.size())) {
+	throw std::runtime_error("Wrong size of ZCORN field");
+    }
+
+    // Make the coordinate array.
+    int delta[3] = { 1,
+		     2*logical_gridsize_[0],
+		     4*logical_gridsize_[0]*logical_gridsize_[1] };
+    int ix = 2*(i*delta[0] + j*delta[1] + k*delta[2]);
+    boost::array<double, 8> cellz = {{ z[ix], z[ix + delta[0]],
+				       z[ix + delta[1]], z[ix + delta[1] + delta[0]],
+				       z[ix + delta[2]], z[ix + delta[2] + delta[0]],
+				       z[ix + delta[2] + delta[1]], z[ix + delta[2] + delta[1] + delta[0]] }};
+    return cellz;
+}
+
+
+} // namespace Dune

dune/common/EclipseGridInspector.hpp

@@ -0,0 +1,105 @@
+//===========================================================================
+//                                                                           
+// File: EclipseGridInspector.h                                              
+//                                                                           
+// Created: Mon Jun  2 09:46:08 2008                                         
+//                                                                           
+// Author: Atgeirr F Rasmussen <atgeirr@sintef.no>
+//
+// $Date$
+//                                                                           
+// Revision: $Id: EclipseGridInspector.h,v 1.2 2008/08/18 14:16:12 atgeirr Exp $
+//                                                                           
+//===========================================================================
+
+/*
+Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+Copyright 2009, 2010 Statoil ASA.
+
+This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+OpenRS is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+OpenRS is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef SINTEF_ECLIPSEGRIDINSPECTOR_HEADER
+#define SINTEF_ECLIPSEGRIDINSPECTOR_HEADER
+
+#include <vector>
+#include <boost/array.hpp>
+
+
+namespace Dune
+{
+
+/**
+   @brief A class for inspecting the contents of an eclipse file.
+
+   Given an EclipseGridParser that has successfully read en Eclipse
+   .grdecl-file, this class may be used to answer certain queries about
+   its contents.
+
+   @author Atgeirr F. Rasmussen <atgeirr@sintef.no>
+   @date 2008/06/02 09:46:08
+*/
+
+class EclipseGridParser;
+
+class EclipseGridInspector
+{
+public:
+    /// Constructor taking a parser as argument.
+    /// The parser must already have read an Eclipse file.
+    EclipseGridInspector(const EclipseGridParser& parser);
+
+    /// Assuming that the pillars are vertical, compute the
+    /// volume of the cell given by logical coordinates (i, j, k).
+    double cellVolumeVerticalPillars(int i, int j, int k) const;
+
+    /// Assuming that the pillars are vertical, compute the
+    /// volume of the cell given by the cell index
+    double cellVolumeVerticalPillars(int cell_idx) const;
+
+    /// Compute the average dip in x- and y-direction of the
+    /// cell tops and bottoms relative to the xy-plane
+    std::pair<double,double> cellDips(int i, int j, int k) const;
+    std::pair<double,double> cellDips(int cell_idx) const;
+
+    // Convert global cell index to logical ijk-coordinates
+    boost::array<int, 3> cellIdxToLogicalCoords(int cell_idx) const;
+    
+    /// Returns a vector with the outer limits of grid (in the grid's unit).
+    /// The vector contains [xmin, xmax, ymin, ymax, zmin, zmax], as 
+    /// read from COORDS and ZCORN
+    boost::array<double, 6> getGridLimits() const;
+
+    /// Returns the extent of the logical cartesian grid
+    /// as number of cells in the (i, j, k) directions.
+    boost::array<int, 3> gridSize() const;
+
+    /// Returns the eight z-values associated with a given cell.
+    /// The ordering is such that i runs fastest. That is, with
+    /// L = low and H = high:
+    /// {LLL, HLL, LHL, HHL, LLH, HLH, LHH, HHH }.
+    boost::array<double, 8> cellZvals(int i, int j, int k) const;
+
+private:
+    const EclipseGridParser& parser_;
+    int logical_gridsize_[3];
+    void checkLogicalCoords(int i, int j, int k) const;
+};
+
+} // namespace Dune
+
+#endif // SINTEF_ECLIPSEGRIDINSPECTOR_HEADER
+

dune/common/EclipseGridParser.cpp

@@ -0,0 +1,571 @@
+//===========================================================================
+//
+// File: EclipseGridParser.C
+//
+// Created: Thu Dec  6 08:46:05 2007
+//
+// Author: Atgeirr F Rasmussen <atgeirr@sintef.no>
+//
+// $Date$
+//
+// $Revision$
+//
+// Revision: $Id: EclipseGridParser.C,v 1.4 2008/08/18 14:16:14 atgeirr Exp $
+//
+//===========================================================================
+
+/*
+  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+  Copyright 2009, 2010 Statoil ASA.
+
+  This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+  OpenRS is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OpenRS is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+#if HAVE_CONFIG_H
+#include "config.h"
+#endif
+#include <iostream>
+#include <fstream>
+#include <exception>
+#include <algorithm>
+#include <limits>
+#include <cfloat>
+#include "EclipseGridParser.hpp"
+#include "EclipseGridParserHelpers.hpp"
+#include "SpecialEclipseFields.hpp"
+#include <dune/common/ErrorMacros.hpp>
+#include <boost/filesystem.hpp>
+#include <dune/common/Units.hpp>
+
+using namespace std;
+
+//#define VERBOSE
+
+namespace Dune
+{
+
+// ---------- List of supported keywords ----------
+
+namespace EclipseKeywords
+{
+    string integer_fields[] =
+	{ string("ACTNUM"),
+	  string("SATNUM"),
+	  string("EQLNUM"),
+	  string("REGNUM"),
+	  string("ROCKTYPE"),
+	  string("DIMENS"),
+	  string("REGDIMS"),
+	  string("WELLDIMS"),
+	  string("TABDIMS"),
+          string("FIPNUM"),
+          string("GRIDFILE")
+	};
+    const int num_integer_fields = sizeof(integer_fields) / sizeof(integer_fields[0]);
+
+    string floating_fields[] =
+        { string("COORD"),    string("ZCORN"),      string("PERMX"),
+          string("PERMY"),    string("PERMZ"),      string("PERMXX"),
+          string("PERMYY"),   string("PERMZZ"),     string("PERMXY"),
+          string("PERMYZ"),   string("PERMZX"),     string("PORO"),
+          string("BULKMOD"),  string("YOUNGMOD"),   string("LAMEMOD"),
+          string("SHEARMOD"), string("POISSONMOD"), string("PWAVEMOD"),
+          string("MULTPV"),   string("PRESSURE"),   string("SGAS"),
+	  string("SWAT")
+        };
+    const int num_floating_fields = sizeof(floating_fields) / sizeof(floating_fields[0]);
+
+    string special_fields[] =
+        { string("SPECGRID"), string("FAULTS"), string("MULTFLT"),
+          string("TITLE"),    string("START"),  string("DATES"),
+	  string("DENSITY"),  string("PVDG"),   string("PVDO"),
+	  string("PVTG"),     string("PVTO"),   string("PVTW"),
+	  string("SGOF"),     string("SWOF"),   string("ROCK"),
+	  string("ROCKTAB"),  string("WELSPECS"), string("COMPDAT"),
+	  string("WCONINJE"), string("WCONPROD"), string("WELTARG"),
+	  string("EQUIL"),    string("PVCDO"),
+	  // The following fields only have a dummy implementation
+	  // that allows us to ignore them.
+          "SWFN",
+          "SOF2",
+          "TUNING"
+        };
+
+    const int num_special_fields = sizeof(special_fields) / sizeof(special_fields[0]);
+
+    string ignore_with_data[] =
+	{ string("MAPUNITS"), string("MAPAXES"),  string("GRIDUNIT"),
+	  string("NTG"),      string("REGDIMS"),  string("WELLDIMS"),
+	  string("NSTACK"),   string("SATNUM"),
+	  string("RPTRST"),   string("ROIP"),     string("RWIP"),
+	  string("RWSAT"),    string("RPR"),      string("WBHP"),
+	  string("WOIR"),     string("TSTEP"),    string("BOX"),
+          string("COORDSYS"), string("PBVD")
+	};
+    const int num_ignore_with_data = sizeof(ignore_with_data) / sizeof(ignore_with_data[0]);
+
+    string ignore_no_data[] =
+	{ string("RUNSPEC"), string("WATER"),    string("OIL"),
+	  string("METRIC"),  string("FMTIN"),    string("FMTOUT"),
+	  string("GRID"),    string("INIT"),     string("NOECHO"),
+	  string("ECHO"),    string("EDIT"),     string("PROPS"),
+	  string("REGIONS"), string("SOLUTION"), string("SUMMARY"),
+	  string("FPR"),     string("FOIP"),     string("FWIP"),
+	  string("RUNSUM"),  string("EXCEL"),    string("SCHEDULE"),
+	  string("END"),     string("ENDBOX"),   string("CONTINUE"),
+          string("NONNC"),   string("GAS"),      string("DISGAS"),
+          string("FIELD")
+	};
+    const int num_ignore_no_data = sizeof(ignore_no_data) / sizeof(ignore_no_data[0]);
+
+    string include_keywords[] = { string("INCLUDE") };
+    const int num_include_keywords = sizeof(include_keywords) / sizeof(include_keywords[0]);
+
+
+} // namespace EclipseKeywords
+
+namespace {
+
+    enum FieldType {
+	Integer,
+	FloatingPoint,
+	SpecialField,
+	IgnoreWithData,
+	IgnoreNoData,
+        Include,
+	Unknown
+    };
+
+    inline FieldType classifyKeyword(const string& keyword)
+    {
+	using namespace EclipseKeywords;
+	if (count(integer_fields, integer_fields + num_integer_fields, keyword)) {
+	    return Integer;
+	} else if (count(floating_fields, floating_fields + num_floating_fields, keyword)) {
+	    return FloatingPoint;
+	} else if (count(special_fields, special_fields + num_special_fields, keyword)) {
+	    return SpecialField;
+	} else if (count(ignore_with_data, ignore_with_data + num_ignore_with_data, keyword)) {
+	    return IgnoreWithData;
+	} else if (count(ignore_no_data, ignore_no_data + num_ignore_no_data, keyword)) {
+	    return IgnoreNoData;
+	} else if (count(include_keywords, include_keywords + num_include_keywords, keyword)) {
+	    return Include;
+	} else {
+	    return Unknown;
+	}
+    }
+
+} // anon namespace
+
+
+
+// ---------- Member functions ----------
+
+/// Default constructor.
+//---------------------------------------------------------------------------
+EclipseGridParser::EclipseGridParser()
+//---------------------------------------------------------------------------
+{
+}
+
+
+/// Constructor taking an eclipse filename.
+//---------------------------------------------------------------------------
+EclipseGridParser::EclipseGridParser(const string& filename, bool convert_to_SI)
+//---------------------------------------------------------------------------
+{
+    // Store directory of filename
+    boost::filesystem::path p(filename);
+    directory_ = p.parent_path().string();
+    ifstream is(filename.c_str());
+    if (!is) {
+	cerr << "Unable to open file " << filename << endl;
+	throw exception();
+    }
+    read(is, convert_to_SI);
+}
+
+
+/// Read the given stream, overwriting any previous data.
+//---------------------------------------------------------------------------
+void EclipseGridParser::read(istream& is, bool convert_to_SI)
+//---------------------------------------------------------------------------
+{
+    integer_field_map_.clear();
+    floating_field_map_.clear();
+    special_field_map_.clear();
+
+    readImpl(is);
+    computeUnits();
+    if (convert_to_SI) {
+        convertToSI();
+    }
+}
+
+//---------------------------------------------------------------------------
+void EclipseGridParser::readImpl(istream& is)
+//---------------------------------------------------------------------------
+{
+    if (!is) {
+	cerr << "Could not read given input stream." << endl;
+	throw exception();
+    }
+
+    // Make temporary maps that will at the end be swapped with the
+    // member maps
+    // NOTE: Above is no longer true, for easier implementation of 
+    //       the INCLUDE keyword. We lose the strong exception guarantee,
+    //       though (of course retaining the basic guarantee).
+    map<string, vector<int> >& intmap = integer_field_map_;
+    map<string, vector<double> >& floatmap = floating_field_map_;
+    map<string, boost::shared_ptr<SpecialBase> >& specialmap = special_field_map_;
+
+    // Actually read the data
+    is >> ignoreWhitespace;
+    while (!is.eof()) {
+	string keyword = readKeyword(is);
+#ifdef VERBOSE
+	cout << "Keyword found: " << keyword << endl;
+#endif
+	FieldType type = classifyKeyword(keyword);
+	switch (type) {
+	case Integer:
+	    readVectorData(is, intmap[keyword]);
+	    break;
+	case FloatingPoint:
+	    readVectorData(is, floatmap[keyword]);
+	    break;
+	case SpecialField: {
+	    boost::shared_ptr<SpecialBase> sb_ptr = createSpecialField(is, keyword);
+	    if (sb_ptr) {
+		specialmap[keyword] = sb_ptr;
+	    } else {
+		THROW("Could not create field " << keyword);
+	    }
+	    break;
+	}
+	case IgnoreWithData: {
+            ignored_fields_.insert(keyword);
+            is >> ignoreSlashLine;
+#ifdef VERBOSE
+            cout << "(ignored)" << endl;
+#endif
+	    break;
+	}
+	case IgnoreNoData: {
+            ignored_fields_.insert(keyword);
+	    is >> ignoreLine;
+#ifdef VERBOSE
+            cout << "(ignored)" << endl;
+#endif
+	    break;
+	}
+        case Include: {
+            string include_filename = readString(is);
+            if (!directory_.empty()) {
+                include_filename = directory_ + '/' + include_filename;
+            }
+            ifstream include_is(include_filename.c_str());
+            if (!include_is) {
+                THROW("Unable to open INCLUDEd file " << include_filename);
+            }
+            readImpl(include_is);
+            is >> ignoreSlashLine;
+            break;
+        }
+	case Unknown:
+	default:
+	    cerr << "Keyword " << keyword << " not recognized." << endl;
+	    throw exception();
+	}
+	is >> ignoreWhitespace;
+    }
+
+#define VERBOSE_LIST_FIELDS 0
+#if VERBOSE_LIST_FIELDS
+    std::cout << "\nInteger fields:\n";
+    for (std::map<string, std::vector<int> >::iterator
+             i = intmap.begin(); i != intmap.end(); ++i)
+        std::cout << '\t' << i->first << '\n';
+
+    std::cout << "\nFloat fields:\n";
+    for (std::map<string, std::vector<double> >::iterator
+             i = floatmap.begin(); i != floatmap.end(); ++i)
+        std::cout << '\t' << i->first << '\n';
+
+    std::cout << "\nSpecial fields:\n";
+    for (std::map<string, boost::shared_ptr<SpecialBase> >::iterator
+              i = specialmap.begin(); i != specialmap.end(); ++i)
+        std::cout << '\t' << i->first << '\n';
+#endif
+}
+
+
+
+//---------------------------------------------------------------------------
+void EclipseGridParser::convertToSI()
+//---------------------------------------------------------------------------
+{
+    // Convert all special fields.
+    typedef std::map<string, boost::shared_ptr<SpecialBase> >::iterator SpecialIt;
+    for (SpecialIt i = special_field_map_.begin(); i != special_field_map_.end(); ++i) {
+        i->second->convertToSI(units_);
+    }
+
+    // Convert all floating point fields.
+    typedef std::map<string, std::vector<double> >::iterator FloatIt;
+    for (FloatIt i = floating_field_map_.begin(); i != floating_field_map_.end(); ++i) {
+        const std::string& key = i->first;
+        std::vector<double>& field = i->second;
+        // Find the right unit.
+        double unit = 1e100;
+        if (key == "COORD" || key == "ZCORN") {
+            unit = units_.length;
+        } else if (key == "PERMX"  || key == "PERMY"  || key == "PERMZ"  ||
+		   key == "PERMXX" || key == "PERMYY" || key == "PERMZZ" ||
+		   key == "PERMXY" || key == "PERMYZ" || key == "PERMZX") {
+            unit = units_.permeability;
+        } else if (key == "PORO"     || key == "BULKMOD"  || key == "YOUNGMOD" ||
+		   key == "LAMEMOD"  || key == "SHEARMOD" || key == "POISSONMOD" ||
+		   key == "PWAVEMOD" || key == "MULTPV"   || key == "PWAVEMOD" ||
+		   key == "SGAS"     || key == "SWAT") {
+            unit = 1.0;
+	} else if (key == "PRESSURE") {
+	    unit = units_.pressure;	    
+        } else {
+            THROW("Units for field " << key << " not specified. Cannon convert to SI.");
+        }
+        // Convert.
+        for (int i = 0; i < int(field.size()); ++i) {
+            field[i] = Dune::unit::convert::from(field[i], unit);
+        }
+    }
+
+    // Set all units to one.
+    units_.setToOne();
+}
+
+
+/// Returns true if the given keyword corresponds to a field that
+/// was found in the file.
+//---------------------------------------------------------------------------
+bool EclipseGridParser::hasField(const string& keyword) const
+//---------------------------------------------------------------------------
+{
+    string ukey = upcase(keyword);
+    return integer_field_map_.count(ukey) || floating_field_map_.count(ukey) ||
+	special_field_map_.count(ukey) || ignored_fields_.count(ukey);
+}
+
+
+/// Returns true if all the given keywords correspond to fields
+/// that were found in the file.
+//---------------------------------------------------------------------------
+bool EclipseGridParser::hasFields(const vector<string>& keywords) const
+//---------------------------------------------------------------------------
+{
+    int num_keywords = keywords.size();
+    for (int i = 0; i < num_keywords; ++i) {
+	if (!hasField(keywords[i])) {
+	    return false;
+	}
+    }
+    return true;
+}
+
+//---------------------------------------------------------------------------
+vector<string> EclipseGridParser::fieldNames() const
+//---------------------------------------------------------------------------
+{
+    vector<string> names;
+    names.reserve(integer_field_map_.size() +
+                  floating_field_map_.size() +
+                  special_field_map_.size() +
+                  ignored_fields_.size());
+    {
+	map<string, vector<int> >::const_iterator it = integer_field_map_.begin();
+	for (; it != integer_field_map_.end(); ++it) {
+	    names.push_back(it->first);
+	}
+    }
+    {
+	map<string, vector<double> >::const_iterator it = floating_field_map_.begin();
+	for (; it != floating_field_map_.end(); ++it) {
+	    names.push_back(it->first);
+	}
+    }
+    {
+	map<string, boost::shared_ptr<SpecialBase> >::const_iterator it = special_field_map_.begin();
+	for (; it != special_field_map_.end(); ++it) {
+	    names.push_back(it->first);
+	}
+    }
+    {
+	set<string>::const_iterator it = ignored_fields_.begin();
+	for (; it != ignored_fields_.end(); ++it) {
+	    names.push_back(*it);
+	}
+    }
+    return names;
+}
+
+//---------------------------------------------------------------------------
+const std::vector<int>& EclipseGridParser::getIntegerValue(const std::string& keyword) const
+//---------------------------------------------------------------------------
+{
+    if (keyword == "SPECGRID") {
+	cerr << "\nERROR. Interface has changed!\n"
+	     << "const vector<int>& dim = parser.getIntegerValue(""SPECGRID"") is deprecated.\n"
+	     << "Use:\n"
+	     << "const SPECGRID& specgrid = parser.getSPECGRID();\n"
+	     << "const vector<int>& dim = specgrid.dimensions;\n\n";
+	throw exception();
+    }
+
+    map<string, vector<int> >::const_iterator it
+	= integer_field_map_.find(keyword);
+    if (it == integer_field_map_.end()) {
+ 	return empty_integer_field_;
+    } else {
+	return it->second;
+    }
+}
+
+//---------------------------------------------------------------------------
+const std::vector<double>& EclipseGridParser::getFloatingPointValue(const std::string& keyword) const
+//---------------------------------------------------------------------------
+{
+    map<string, vector<double> >::const_iterator it
+	= floating_field_map_.find(keyword);
+    if (it == floating_field_map_.end()) {
+ 	return empty_floating_field_;
+    } else {
+	return it->second;
+    }
+}
+
+
+//---------------------------------------------------------------------------
+const boost::shared_ptr<SpecialBase> EclipseGridParser::getSpecialValue(const std::string& keyword) const
+//---------------------------------------------------------------------------
+{
+    map<string, boost::shared_ptr<SpecialBase> >::const_iterator it = special_field_map_.find(keyword);
+    if (it == special_field_map_.end()) {
+        THROW("No such field: " << keyword);
+    } else {
+	return it->second;
+    }
+}
+
+//---------------------------------------------------------------------------
+boost::shared_ptr<SpecialBase>
+EclipseGridParser::createSpecialField(std::istream& is,
+				      const std::string& fieldname)
+//---------------------------------------------------------------------------
+{
+    string ukey = upcase(fieldname);
+    boost::shared_ptr<SpecialBase> spec_ptr
+        = Factory<SpecialBase>::createObject(fieldname);
+    spec_ptr->read(is);
+    return spec_ptr;
+}
+
+//---------------------------------------------------------------------------
+void EclipseGridParser::setIntegerField(const std::string& keyword,
+                                        const std::vector<int>& field)
+//---------------------------------------------------------------------------
+{
+    integer_field_map_[keyword] = field;
+}
+
+//---------------------------------------------------------------------------
+void EclipseGridParser::setFloatingPointField(const std::string& keyword,
+                                              const std::vector<double>& field)
+//---------------------------------------------------------------------------
+{
+    floating_field_map_[keyword] = field;
+}
+
+//---------------------------------------------------------------------------
+void EclipseGridParser::setSpecialField(const std::string& keyword,
+                                        boost::shared_ptr<SpecialBase> field)
+//---------------------------------------------------------------------------
+{
+    special_field_map_[keyword] = field;
+}
+
+//---------------------------------------------------------------------------
+const EclipseUnits& EclipseGridParser::units() const
+//---------------------------------------------------------------------------
+{
+    return units_;
+}
+
+//---------------------------------------------------------------------------
+void EclipseGridParser::computeUnits()
+//---------------------------------------------------------------------------
+{
+    // Decide unit family.
+    enum EclipseUnitFamily { Metric = 0, Field = 1, Lab = 2, Pvtm = 3 };
+    EclipseUnitFamily unit_family = Metric; // The default.
+    if (hasField("FIELD")) unit_family = Field;
+    if (hasField("LAB")) unit_family = Lab;
+    if (hasField("PVT-M")) unit_family = Pvtm;
+
+    // Set units.
+    using namespace prefix;
+    using namespace unit;
+    switch (unit_family) {
+    case Metric:
+        units_.length = meter;
+        units_.time = day;
+        units_.density = kilogram/cubic(meter);
+        units_.pressure = barsa;
+        units_.compressibility = 1.0/barsa;
+        units_.viscosity = centi*Poise;
+        units_.permeability = milli*darcy;
+        units_.liqvol_s = cubic(meter);
+        units_.liqvol_r = cubic(meter);
+        units_.gasvol_s = cubic(meter);
+        units_.gasvol_r = cubic(meter);
+        units_.transmissibility = centi*Poise * cubic(meter) / (day * barsa);
+        break;
+    case Field:
+        units_.length = feet;
+        units_.time = day;
+        units_.density = pound/cubic(feet);
+        units_.pressure = psia;
+        units_.compressibility = 1.0/psia;
+        units_.viscosity = centi*Poise;
+        units_.permeability = milli*darcy;
+        units_.liqvol_s = stb;
+        units_.liqvol_r = stb;
+        units_.gasvol_s = 1000*cubic(feet);  // Prefix 'M' is 1000
+        units_.gasvol_r = stb;
+        units_.transmissibility = centi*Poise * stb / (day * psia);
+        break;
+    case Lab:
+        THROW("Unhandled unit family " << unit_family);
+        break;
+    case Pvtm:
+        THROW("Unhandled unit family " << unit_family);
+        break;
+    default:
+        THROW("Unknown unit family " << unit_family);
+    }
+}
+
+} // namespace Dune

dune/common/EclipseGridParser.hpp

@@ -0,0 +1,191 @@
+//===========================================================================
+//                                                                           
+// File: EclipseGridParser.h                                                 
+//                                                                           
+// Created: Wed Dec  5 17:05:13 2007                                         
+//                                                                           
+// Author: Atgeirr F Rasmussen <atgeirr@sintef.no>
+//
+// $Date$
+//                                                                           
+// Revision: $Id: EclipseGridParser.h,v 1.3 2008/08/18 14:16:13 atgeirr Exp $
+//                                                                           
+//===========================================================================
+
+/*
+Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+Copyright 2009, 2010 Statoil ASA.
+
+This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+OpenRS is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+OpenRS is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef SINTEF_ECLIPSEGRIDPARSER_HEADER
+#define SINTEF_ECLIPSEGRIDPARSER_HEADER
+
+#include <iosfwd>
+#include <string>
+#include <vector>
+#include <map>
+#include <set>
+#include <boost/shared_ptr.hpp>
+#include "SpecialEclipseFields.hpp"
+#include "EclipseUnits.hpp"
+#include <dune/common/Factory.hpp>
+
+namespace Dune
+{
+
+/**
+   @brief A class for reading and parsing all fields of an eclipse file.
+   
+   This object is constructed using an Eclipse .grdecl-file. All data
+   fields are extracted upon construction and written to vector data
+   structures, which can then be read out in O(1) time afterwards via
+   convenience functions.
+
+   There is also a convenience function to easily check which fields
+   were successfully parsed.
+
+   @author Atgeirr F. Rasmussen <atgeirr@sintef.no>
+   @date 2007/12/06 13:00:03
+
+*/
+
+class EclipseGridParser
+{
+public:
+    /// Default constructor.
+    EclipseGridParser();
+    /// Constructor taking an eclipse filename. Unless the second
+    /// argument 'convert_to_SI' is false, all fields will be
+    /// converted to SI units.
+    explicit EclipseGridParser(const std::string& filename, bool convert_to_SI = true);
+
+    /// Read the given stream, overwriting any previous data.  Unless
+    /// the second argument 'convert_to_SI' is false, all fields will
+    /// be converted to SI units.
+    void read(std::istream& is, bool convert_to_SI = true);
+
+    /// Convert all data to SI units, according to unit category as
+    /// specified in the eclipse file (METRIC, FIELD etc.).
+    /// It is unnecessary to call this if constructed with
+    /// 'convert_to_SI' equal to true, but it is not an error.
+    void convertToSI();
+
+    /// Returns true if the given keyword corresponds to a field that
+    /// was found in the file.
+    bool hasField(const std::string& keyword) const;
+    /// Returns true if all the given keywords correspond to fields
+    /// that were found in the file.
+    bool hasFields(const std::vector<std::string>& keywords) const;
+    /// The keywords/fields found in the file.
+    std::vector<std::string> fieldNames() const;
+
+    /// Returns a reference to a vector containing the values
+    /// corresponding to the given integer keyword.
+    const std::vector<int>& getIntegerValue(const std::string& keyword) const;
+
+    /// Returns a reference to a vector containing the values
+    /// corresponding to the given floating-point keyword.
+    const std::vector<double>& getFloatingPointValue(const std::string& keyword) const;
+
+    /// Returns a reference to a vector containing pointers to the values 
+    /// corresponding to the given keyword when the values are not only integers
+    /// or floats.
+    const boost::shared_ptr<SpecialBase> getSpecialValue(const std::string& keyword) const;
+
+    // This macro implements support for a special field keyword. It requires that a subclass
+    // of SpecialBase exists, that has the same name as the keyword.
+    // After using SPECIAL_FIELD(KEYWORD), the public member getKEYWORD will be available.
+#define SPECIAL_FIELD(keyword)                                                                   \
+private:                                                                                         \
+    struct X##keyword { X##keyword() { Factory<SpecialBase>::addCreator<keyword>(#keyword); } }; \
+    X##keyword x##keyword;                                                                       \
+public:                                                                                          \
+    const keyword& get##keyword() const                                                          \
+    { return dynamic_cast<const keyword&>(*getSpecialValue(#keyword)); }
+
+    // Support for special fields.
+    SPECIAL_FIELD(SPECGRID);
+    SPECIAL_FIELD(FAULTS);
+    SPECIAL_FIELD(MULTFLT);
+    SPECIAL_FIELD(TITLE);
+    SPECIAL_FIELD(START);
+    SPECIAL_FIELD(DATES);
+    SPECIAL_FIELD(DENSITY);
+    SPECIAL_FIELD(PVDG);
+    SPECIAL_FIELD(PVDO);
+    SPECIAL_FIELD(PVTG);
+    SPECIAL_FIELD(PVTO);
+    SPECIAL_FIELD(PVTW);
+    SPECIAL_FIELD(SGOF);
+    SPECIAL_FIELD(SWOF);
+    SPECIAL_FIELD(ROCK);
+    SPECIAL_FIELD(ROCKTAB);
+    SPECIAL_FIELD(WELSPECS);
+    SPECIAL_FIELD(COMPDAT);
+    SPECIAL_FIELD(WCONINJE);
+    SPECIAL_FIELD(WCONPROD);
+    SPECIAL_FIELD(WELTARG);
+    SPECIAL_FIELD(EQUIL);
+    SPECIAL_FIELD(PVCDO);
+
+    // The following fields only have a dummy implementation
+    // that allows us to ignore them.
+    SPECIAL_FIELD(SWFN);
+    SPECIAL_FIELD(SOF2);
+    SPECIAL_FIELD(TUNING);
+
+#undef SPECIAL_FIELD
+
+
+    /// Sets an integer field to have a particular value.
+    void setIntegerField(const std::string& keyword, const std::vector<int>& field);
+
+    /// Sets a floating point field to have a particular value.
+    void setFloatingPointField(const std::string& keyword, const std::vector<double>& field);
+
+    /// Sets a special field to have a particular value.
+    void setSpecialField(const std::string& keyword, boost::shared_ptr<SpecialBase> field);
+
+    /// Compute the units used by the deck, depending on the presence
+    /// of keywords such as METRIC, FIELD etc.  It is an error to call
+    /// this after conversion to SI has taken place.
+    void computeUnits();
+
+    /// The units specified by the eclipse file read.
+    const EclipseUnits& units() const;
+
+private:
+    boost::shared_ptr<SpecialBase> createSpecialField(std::istream& is, const std::string& fieldname);
+    void readImpl(std::istream& is);
+
+
+    std::string directory_;
+    std::map<std::string, std::vector<int> > integer_field_map_;
+    std::map<std::string, std::vector<double> > floating_field_map_;
+    std::map<std::string, boost::shared_ptr<SpecialBase> > special_field_map_;
+    std::set<std::string> ignored_fields_;
+    std::vector<int> empty_integer_field_;
+    std::vector<double> empty_floating_field_;
+    boost::shared_ptr<SpecialBase> empty_special_field_;
+    EclipseUnits units_;
+};
+
+
+} // namespace Dune
+
+#endif // SINTEF_ECLIPSEGRIDPARSER_HEADER

dune/common/EclipseGridParserHelpers.hpp

@@ -0,0 +1,469 @@
+//===========================================================================
+//
+// File: EclipseGridParserHelpers.hpp
+//
+// Created: Tue Dec 22 11:35:32 2009
+//
+// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
+//            Bård Skaflestad     <bard.skaflestad@sintef.no>
+//
+// $Date$
+//
+// $Revision$
+//
+//===========================================================================
+
+/*
+  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+  Copyright 2009, 2010 Statoil ASA.
+
+  This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+  OpenRS is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OpenRS is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPENRS_ECLIPSEGRIDPARSERHELPERS_HEADER
+#define OPENRS_ECLIPSEGRIDPARSERHELPERS_HEADER
+
+#include <limits>
+#include <string>
+#include <istream>
+#include <vector>
+#include <dune/common/ErrorMacros.hpp>
+#include "linInt.hpp"
+#include <boost/date_time/gregorian/gregorian.hpp>
+
+namespace Dune
+{
+
+namespace
+{
+
+    inline std::istream& ignoreLine(std::istream& is)
+    {
+	is.ignore(std::numeric_limits<int>::max(), '\n');
+	return is;
+    }
+
+    inline std::istream& ignoreSlashLine(std::istream& is)
+    {
+	is.ignore(std::numeric_limits<int>::max(), '/');
+	is.ignore(std::numeric_limits<int>::max(), '\n');
+	return is;
+    }
+
+    inline std::istream& ignoreWhitespace(std::istream& is)
+    {
+	// Getting the character type facet for is()
+	// We use the classic (i.e. C) locale.
+	const std::ctype<char>& ct = std::use_facet< std::ctype<char> >(std::locale::classic());
+	char c;
+	while (is.get(c)) {
+	    if (!ct.is(std::ctype_base::space, c)) {
+		is.putback(c);
+		break;
+	    }
+	}
+	return is;
+    }
+
+    inline std::string upcase(const std::string& s)
+    {
+	std::string us(s);
+	// Getting the character type facet for toupper().
+	// We use the classic (i.e. C) locale.
+	const std::ctype<char>& ct = std::use_facet< std::ctype<char> >(std::locale::classic());
+	for (int i = 0; i < int(s.size()); ++i) {
+	    us[i] = ct.toupper(s[i]);
+	}
+	return us;
+    }
+
+    inline std::string readKeyword(std::istream& is)
+    {
+	std::string keyword_candidate;
+	while (!is.eof()) {
+	    is >> keyword_candidate;
+	    if(keyword_candidate.find("--") == 0) {
+		is >> ignoreLine;  // This line is a comment
+	    } else {
+		return upcase(keyword_candidate);
+	    }
+	}
+	return "CONTINUE";  // Last line in included file is a comment
+    }
+
+    inline std::string readString(std::istream& is)
+    {
+	std::string string_candidate;
+	is >> string_candidate;
+        const char quote('\'');
+        int beg = string_candidate[0] == quote ? 1 : 0;
+        int len = string_candidate[0] == quote ? string_candidate.size() - 2 : string_candidate.size();
+        return string_candidate.substr(beg, len);
+    }
+
+    // Reads data until '/' or an error is encountered.
+    template<typename T>
+    inline void readVectorData(std::istream& is, std::vector<T>& data)
+    {
+	data.clear();
+	while (is) {
+	    T candidate;
+	    is >> candidate;
+	    if (is.rdstate() & std::ios::failbit) {
+		is.clear(is.rdstate() & ~std::ios::failbit);
+                is >> ignoreWhitespace;
+		char dummy;
+		is >> dummy;
+		if (dummy == '/') {
+                    is >> ignoreLine;
+		    break;
+		} else if (dummy == '-') {  // "comment test"
+		    is >> ignoreLine; // This line is a comment
+		} else {
+                    THROW("Encountered format error while reading data values. Value = " << dummy);
+		}
+	    } else {
+		if (is.peek() == int('*')) {
+		    is.ignore(); // ignore the '*'
+		    int multiplier = int(candidate);
+		    is >> candidate;
+		    data.insert(data.end(), multiplier, candidate);
+		} else {
+		    data.push_back(candidate);
+		}
+	    }
+	}
+	if (!is) {
+	    THROW("Encountered error while reading data values.");
+	}
+    }
+
+
+    // Reads data items of type T. Not more than 'max_values' items.
+    // Asterisks may be used to signify 'repeat counts'. 5*3.14 will
+    // insert 3.14 five times. Asterisk followed by a space is used to
+    // signify default values.  n* will default n consecutive quantities.
+    template<class Vec>
+    inline int readDefaultedVectorData(std::istream& is, Vec& data, int max_values)
+    {
+        ASSERT(int(data.size()) >= max_values);
+	const std::ctype<char>& ct = std::use_facet< std::ctype<char> >(std::locale::classic());
+        int num_values = 0;
+        while (is) {
+            typename Vec::value_type candidate;
+            is >> candidate;
+            if (is.rdstate() & std::ios::failbit) {
+                is.clear(is.rdstate() & ~std::ios::failbit);
+                std::string dummy;
+                is >> dummy;
+                if (dummy == "/") {
+                    is >> ignoreLine;	
+                    break;
+		} else if (dummy[0] == '-') {  // "comment test"
+                    is >> ignoreLine;   // This line is a comment
+                } else {
+                    THROW("Encountered format error while reading data values. Value = " << dummy);
+                }
+            } else {
+                if (is.peek() == int('*')) {
+                    is.ignore(); // ignore the '*'
+                    int multiplier = (int)candidate;
+                    if (ct.is(std::ctype_base::space, is.peek())) {
+                        num_values += multiplier;  // Use default value(s)
+                    } else {
+                        is >> candidate;         // Use candidate 'multipler' times
+                        for (int i=0; i<multiplier; ++i, ++num_values) {
+                            data[num_values] = candidate;
+                        }
+                    }
+                } else {
+                    data[num_values] = candidate;
+                    ++num_values;
+                }
+            }
+            if (num_values >= max_values) {
+		//is >> ignoreLine;
+                break;
+            }
+        }
+        if (!is) {
+            THROW("Encountered error while reading data values.");
+        }
+        return num_values;
+    }
+
+
+    // Keywords SGOF and SWOF. Reads data until '/' or an error is encountered.
+    // Default values represented by 1* is replaced by -1. Use linear interpolation
+    // outside this function to replace -1.
+    template<typename T>
+    inline void readRelPermTable(std::istream& is, std::vector<T>& data)
+    {
+	data.clear();
+	while (is) {
+	    T candidate;
+	    is >> candidate;
+	    if (is.rdstate() & std::ios::failbit) {
+		is.clear(is.rdstate() & ~std::ios::failbit);
+		std::string dummy;
+		is >> dummy;
+		if (dummy == "/") {
+                    is >> ignoreLine;
+		    break;
+		} else if (dummy[0] == '-') {  // "comment test"
+		    is >> ignoreLine; // This line is a comment
+		} else {
+                    THROW("Encountered format error while reading data values. Value = " << dummy);
+		}
+	    } else {
+		if (is.peek() == int('*')) {
+		    is.ignore(); // ignore the '*'
+		    ASSERT(int(candidate) == 1);
+		    data.push_back(-1); // Set new flag for interpolation.
+		} else {
+		    data.push_back(candidate);
+		}
+	    }
+	}
+	if (!is) {
+	    THROW("Encountered error while reading data values.");
+	}
+    }
+
+
+    // Returns month number 1-12. Returns 0 if illegal month name. 
+    inline int getMonthNumber(const std::string& month_name)
+    {
+        const int num_months = 12;
+        std::string months[num_months] = {"JAN", "FEB", "MAR", "APR", "MAY", "JUN",
+                                          "JLY", "AUG", "SEP", "OCT", "NOV", "DEC"};
+        if (month_name == "JUL") {
+            return 7;    // "JUL" is an acceptable alternative to 'JLY'
+        }
+        int m = 0;
+        for (int i=0; i<num_months; ++i) {
+            if (month_name == months[i]) {
+                m = i+1;
+                break;
+            }
+        }
+        return m;
+    }
+
+
+    inline boost::gregorian::date readDate(std::istream& is)
+    {
+        while (is.peek() == int('-')) {
+            is >> ignoreLine;   // This line is a comment
+        }
+        int day, year;
+        std::string month_name;
+        is >> day;
+        month_name = readString(is);
+        is >> year;
+        ignoreSlashLine(is);
+        int month = getMonthNumber(month_name);
+        return boost::gregorian::date(boost::gregorian::greg_year(year),
+                                      boost::gregorian::greg_month(month),
+                                      boost::gregorian::greg_day(day));
+    }
+
+    // Get first character after whitespace and comments, and decide
+    // next action.
+    // NB! Will fail for negative number in column one.
+    inline int next_action(std::istream& is)
+    {
+	int task(0);  // 0:continue  1:return  2:throw
+	const std::ctype<char>& ct =
+	    std::use_facet< std::ctype<char> >(std::locale::classic());
+
+	while (!is.eof()) {
+	    is >> ignoreWhitespace;
+	    char c;
+	    is.get(c);
+	    if (is.eof()) {
+		return task;
+	    }
+	    is.putback(c);
+	    if (ct.is(std::ctype_base::digit, c) || c== '.') {
+		task = 0;   // Decimal digit. Read more data.
+		break;
+	    }
+	    if (ct.is(std::ctype_base::alpha, c)) {
+		task = 1;   // Alphabetic char. Read next keyword.
+		break;
+	    }
+	    if (c == '-') {
+		is >> ignoreLine;  // This line is a comment
+	    } else {
+		task = 2;  // Read error. Unexpected character.
+		break;
+	    }
+	}
+	return task;
+    }
+
+    // Reads keywords PVTG, PVTO
+    typedef std::vector<std::vector<std::vector<double> > > table_t;
+    inline void readPvtTable(std::istream& is, table_t& pvt_table,
+			     const std::string& field_name)
+    {
+	const std::ctype<char>& ct =
+	    std::use_facet< std::ctype<char> >(std::locale::classic());
+	std::vector<double> record;
+	std::vector<std::vector<double> > table;	
+	while (!is.eof()) {
+	    record.clear();
+	    readVectorData(is, record);
+	    table.push_back(record);
+	    while (!is.eof()) {
+		is >> ignoreWhitespace;
+		char c;
+		is.get(c);
+                if (is.eof()) {
+                    // Reached end of file, we should have pushed
+                    // the last table, and emptied it. If not,
+                    // we have an error.
+                    if (!table.empty()) {
+                        THROW("Reached EOF while still building PVT table. Missing end-of-table (slash)?");
+                    }
+                    return;
+                }
+		is.putback(c);
+		if (ct.is(std::ctype_base::digit, c) || c== '.') {
+		    break;   // Decimal digit. Read more records.
+		}
+		if (ct.is(std::ctype_base::alpha, c)) {
+		    return;   // Alphabetic char. Read next keyword.
+		}
+		if (c == '-') {
+		    is >> ignoreLine;  // This line is a comment
+		    continue;
+		}
+		if (c == '/') {
+		    is >> ignoreLine;
+		    pvt_table.push_back(table);
+		    table.clear();
+		} else {
+		    std::ostringstream oss;
+		    oss << "Error reading " << field_name
+		       << ". Next character is " <<  (char)is.peek();
+		    THROW(oss.str());
+		}
+	    }
+	}
+    }
+
+    // Reads keywords PVDG, PVDO, ROCKTAB
+    inline void readPvdTable(std::istream& is, table_t& pvd_table,
+			     const std::string& field_name, int ncol)
+    {
+	std::vector<double> record;
+	std::vector<std::vector<double> > table(ncol);	
+	while (!is.eof()) {
+	    record.clear();
+	    readVectorData(is, record);
+	    const int rec_size = record.size()/ncol;
+	    for (int k=0; k<ncol; ++k) {
+		table[k].resize(rec_size);
+	    }
+	    for (int i=0, n=-1; i<rec_size; ++i) {
+		for (int k=0; k<ncol; ++k) {
+		    table[k][i] = record[++n];
+		}
+	    }
+	    pvd_table.push_back(table);
+
+	    int action = next_action(is); // 0:continue  1:return  2:throw
+	    if (action == 1) {
+		return;     // Alphabetic char. Read next keyword.
+	    } else if (action == 2) {
+		std::ostringstream oss;
+		oss << "Error reading " << field_name
+		    << ". Next character is " <<  (char)is.peek();
+		THROW(oss.str());
+	    }
+	}
+    }
+
+    // Replace default values -1 by linear interpolation
+    inline void insertDefaultValues(std::vector<std::vector<double> >& table, int ncol)
+    {
+	const int sz = table[0].size();
+	for (int k=1; k<ncol; ++k) {
+	    std::vector<int> indx;
+	    std::vector<double> x;
+	    for (int i=0; i<sz; ++i) {
+		if (table[k][i] == -1) {
+		    indx.push_back(i);
+		    x.push_back(table[0][i]);    
+		}
+	    }
+	    if (!indx.empty()) {
+		std::vector<double> xv, yv;
+		for (int i=0; i<sz; ++i) {
+		    if (table[k][i] != -1) {
+			xv.push_back(table[0][i]);		    
+			yv.push_back(table[k][i]);		    
+		    }
+		}
+		// Interpolate
+		for (int i=0; i<int(indx.size()); ++i) {
+		    table[k][indx[i]] = linearInterpolationExtrap(xv, yv, x[i]);
+		}
+	    }
+	}
+    }
+
+    // Reads keywords SGOF and SWOF
+    inline void readSGWOF(std::istream& is, table_t& relperm_table,
+			  const std::string& field_name, int ncol)
+    {
+	std::vector<double> record;
+	std::vector<std::vector<double> > table(ncol);	
+	while (!is.eof()) {
+	    record.clear();
+	    readRelPermTable(is, record);
+	    const int rec_size = record.size()/ncol;
+	    for (int k=0; k<ncol; ++k) {
+		table[k].resize(rec_size);
+	    }
+	    for (int i=0, n=-1; i<rec_size; ++i) {
+		for (int k=0; k<ncol; ++k) {
+		    table[k][i] = record[++n];
+		}
+	    }
+	    insertDefaultValues(table, ncol);
+	    relperm_table.push_back(table);
+
+	    int action = next_action(is); // 0:continue  1:return  2:throw
+	    if (action == 1) {
+		return;     // Alphabetic char. Read next keyword.
+	    } else if (action == 2) {
+		std::ostringstream oss;
+		oss << "Error reading " << field_name
+		    << ". Next character is " <<  (char)is.peek();
+		THROW(oss.str());
+	    }
+	}
+    }
+
+} // anon namespace
+
+} // namespace Dune
+
+
+#endif // OPENRS_ECLIPSEGRIDPARSERHELPERS_HEADER

dune/common/EclipseUnits.hpp

@@ -0,0 +1,61 @@
+/*
+  Copyright 2010 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_ECLIPSEUNITS_HEADER_INCLUDED
+#define OPM_ECLIPSEUNITS_HEADER_INCLUDED
+
+namespace Dune
+{
+    struct EclipseUnits
+    {
+        double length;
+        double time;
+        double density;
+        double pressure;
+        double compressibility;
+        double viscosity;
+        double permeability;
+        double liqvol_s;
+        double liqvol_r;
+        double gasvol_s;
+        double gasvol_r;
+        double transmissibility;
+
+        void setToOne()
+        {
+            length = 1.0;
+            time = 1.0;
+            density = 1.0;
+            pressure = 1.0;
+            compressibility = 1.0;
+            viscosity = 1.0;
+            permeability = 1.0;
+            liqvol_s = 1.0;
+            liqvol_r = 1.0;
+            gasvol_s = 1.0;
+            gasvol_r = 1.0;
+            transmissibility = 1.0;
+        }
+    };
+
+
+} // namespace Dune
+
+
+#endif // OPM_ECLIPSEUNITS_HEADER_INCLUDED

dune/common/ErrorMacros.hpp

@@ -0,0 +1,107 @@
+//===========================================================================
+//
+// File: ErrorMacros.hpp
+//
+// Created: Tue May 14 12:22:16 2002
+//
+// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
+//
+// $Date$
+//
+// $Revision$
+//
+//===========================================================================
+
+/*
+Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+Copyright 2009, 2010 Statoil ASA.
+
+This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+OpenRS is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+OpenRS is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPENRS_ERRORMACROS_HEADER
+#define OPENRS_ERRORMACROS_HEADER
+
+
+/// Error macros. In order to use some of them, you must also
+/// include <iostream> or <exception>, so they are included by
+/// this file. The compile defines NDEBUG and NVERBOSE control
+/// the behaviour of these macros.
+
+#ifndef NVERBOSE
+  #include <iostream>
+#endif
+#include <exception>
+
+/// Usage: REPORT;
+/// Usage: MESSAGE("Message string.");
+#ifdef NVERBOSE // Not verbose mode
+#  ifndef REPORT
+#    define REPORT
+#  endif
+#  ifndef MESSAGE
+#    define MESSAGE(x)
+#  endif
+#  ifndef MESSAGE_IF
+#    define MESSAGE_IF(cond, m)
+#  endif
+#else // Verbose mode
+#  ifndef REPORT
+#    define REPORT std::cerr << "\nIn file " << __FILE__ << ", line " << __LINE__ << std::endl
+#  endif
+#  ifndef MESSAGE
+#    define MESSAGE(x) std::cerr << "\nIn file " << __FILE__ << ", line " << __LINE__ << ": " << x << std::endl
+#  endif
+#  ifndef MESSAGE_IF
+#    define MESSAGE_IF(cond, m) do {if(cond) MESSAGE(m);} while(0)
+#  endif
+#endif
+
+
+/// Usage: THROW("Error message string.");
+#ifndef THROW
+#  define THROW(x) do { MESSAGE(x); throw std::exception(); } while(0)
+#endif
+
+#define ALWAYS_ERROR_IF(condition, message) do {if(condition){ THROW(message);}} while(0)
+
+/// Usage: ASSERT(condition)
+/// Usage: ASSERT2(condition, "Error message string.")
+/// Usage: DEBUG_ERROR_IF(condition, "Error message string.");
+#ifdef NDEBUG // Not in debug mode
+#  ifndef ASSERT
+#    define ASSERT(x)
+#  endif
+#  ifndef ASSERT2
+#    define ASSERT2(cond, x)
+#  endif
+#  ifndef DEBUG_ERROR_IF
+#    define DEBUG_ERROR_IF(cond, x)
+#  endif
+#else // Debug mode
+#  ifndef ASSERT
+#    define ASSERT(cond) if (!(cond)) THROW("Assertion \'" #cond "\' failed.")
+#  endif
+#  ifndef ASSERT2
+#    define ASSERT2(cond, x) do { if (!(cond)) THROW(x);} while(0)
+#  endif
+#  ifndef DEBUG_ERROR_IF
+#    define DEBUG_ERROR_IF(cond, x) do { if (cond) THROW(x); } while(0)
+#  endif
+#endif
+
+
+#endif // OPENRS_ERRORMACROS_HEADER

dune/common/Factory.hpp

@@ -0,0 +1,145 @@
+//===========================================================================
+//
+// File: Factory.hpp
+//
+// Created: Mon Nov  6 10:00:43 2000
+//
+// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
+//
+// $Date$
+//
+// $Revision$
+//
+//===========================================================================
+
+/*
+  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+  Copyright 2009, 2010 Statoil ASA.
+
+  This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+  OpenRS is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OpenRS is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPENRS_FACTORY_HEADER
+#define OPENRS_FACTORY_HEADER
+
+
+#include <map>
+#include <boost/shared_ptr.hpp>
+
+namespace Dune
+{
+    
+    /** This is an object factory for creating objects of some type
+     *  requested by the user, with a shared base class.  The user
+     *  need only interact with the factory through the static
+     *  template member addCreator() and the static member function
+     *  createObject().
+     */
+    template <class Base>
+    class Factory
+    {
+    public:
+        /// The type of pointer returned by createObject().
+        typedef boost::shared_ptr<Base> ProductPtr;
+
+	/// Creates a new object of the class associated with the given type string,
+        /// and returns a pointer to it.
+	/// \param type the type string of the class that the user wants to have
+	///             constructed.
+        /// \return (smart) pointer to the created object.
+	static ProductPtr createObject(const std::string& type)
+	{
+	    return instance().doCreateObject(type);
+	}
+
+	/// Add a creator to the Factory.
+        /// After the call, the user may obtain new objects of the Derived type by
+        /// calling createObject() with the given type string as an argument.
+        /// \tparam Derived the class we want to add a creator for, must inherit
+        ///                 the class template parameter Base.
+	/// \param type the type string with which we want the Factory to associate
+	///             the class Derived.
+        template <class Derived>
+	static void addCreator(const std::string& type)
+	{
+            instance().doAddCreator<Derived>(type);
+	}
+
+    private:
+        // The method that implements the singleton pattern,
+        // using the Meyers singleton technique.
+	static Factory& instance()
+	{
+            static Factory singleton;
+	    return singleton;
+	}
+
+	// Private constructor, to keep users from creating a Factory.
+	Factory()
+	{
+	}
+
+        // Abstract base class for Creators.
+        class Creator
+        {
+        public:
+            virtual ProductPtr create() = 0;
+            virtual ~Creator() {}
+        };
+
+        /// This is the concrete Creator subclass for generating Derived objects.
+        template <class Derived>
+        class ConcreteCreator : public Creator
+        {
+        public:
+            virtual ProductPtr create()
+            {
+                return ProductPtr(new Derived);
+            }
+        };
+
+        typedef boost::shared_ptr<Creator> CreatorPtr;
+        typedef std::map<std::string, CreatorPtr> CreatorMap;
+        // This map contains the whole factory, i.e. all the Creators.
+	CreatorMap string_to_creator_;
+
+        // Actually creates the product object.
+	ProductPtr doCreateObject(const std::string& type)
+	{
+	    typename CreatorMap::iterator it;
+	    it = string_to_creator_.find(type);
+	    if (it == string_to_creator_.end()) {
+		THROW("Creator type " << type
+		      << " is not registered in the factory.");
+	    }
+	    return it->second->create();
+	}
+
+        // Actually adds the creator.
+        template <class Derived>
+        void doAddCreator(const std::string& type)
+        {
+            boost::shared_ptr<Creator> c(new ConcreteCreator<Derived>);
+	    string_to_creator_[type] = c;
+        }
+    };
+
+
+} // namespace Dune
+
+#endif // OPENRS_FACTORY_HEADER
+
+

dune/common/Makefile.am

@@ -0,0 +1,38 @@
+# $Date$
+# $Revision$
+
+SUBDIRS = . param test
+
+common_HEADERS	= \
+	Average.hpp \
+	EclipseGridInspector.hpp \
+	EclipseGridParser.hpp \
+	EclipseGridParserHelpers.hpp \
+        ErrorMacros.hpp \
+	Factory.hpp \
+	linInt.hpp \
+        MonotCubicInterpolator.hpp \
+	RootFinders.hpp \
+	SparseTable.hpp \
+	SparseVector.hpp \
+	SpecialEclipseFields.hpp \
+	StopWatch.hpp \
+	Units.hpp
+
+commondir = $(includedir)/dune/common
+
+noinst_LTLIBRARIES = libcommon.la
+
+libcommon_la_SOURCES = \
+	EclipseGridInspector.cpp \
+	EclipseGridParser.cpp \
+	MonotCubicInterpolator.cpp \
+	StopWatch.cpp
+
+
+libcommon_la_CPPFLAGS = $(DUNE_CPPFLAGS) $(BOOST_CPPFLAGS)
+libcommon_la_LDFLAGS  = $(DUNE_LDFLAGS) $(BOOST_LDFLAGS) \
+                        $(BOOST_DATE_TIME_LIB)
+
+
+include $(top_srcdir)/am/global-rules

dune/common/MonotCubicInterpolator.cpp

@@ -0,0 +1,726 @@
+/*
+  MonotCubicInterpolator
+  Copyright (C) 2006 Statoil ASA
+  
+  This program 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 2
+  of the License, or (at your option) any later version.
+  
+  This program 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 this program; if not, write to the Free Software
+  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+*/
+
+/**
+   @file MonotCubicInterpolator.C
+   @brief Represents one dimensional function f with single valued argument x
+
+   Class to represent a one-dimensional function with single-valued
+   argument.  Cubic interpolation, preserving the monotonicity of the
+   discrete known function values
+   
+   @see MonotCubicInterpolator.h for further documentation.
+
+*/
+
+
+#include "MonotCubicInterpolator.hpp"
+#include <iostream>
+#include <iomanip>
+#include <fstream>
+#include <sstream>
+#include <string>
+#include <vector>
+#include <map>
+#include <cmath>
+
+using namespace std;
+
+/*
+
+  SOME DISCUSSION ON DATA STORAGE:
+
+  Internal data structure of points and values:
+
+  vector(s): 
+   - Easier coding 
+   - Faster vector operations when setting up derivatives.
+   - sorting slightly more complex.
+   - insertion of further values bad.
+
+  vector<double,double>
+   - easy sorting
+   - code complexity almost as for map.
+   - insertion of additional values bad
+
+  vector over struct datapoint { x, f, d} 
+   - nice code
+   - not as sortable, insertion is cumbersome.
+
+   ** This is used currently: ** 
+  map<double, double> one for (x,f) and one for (x,d)
+   - Naturally sorted on x-values (done by the map-construction)
+   - Slower to set up, awkward loop coding (?)
+   - easy to add more points.
+   - easier to just add code to linear interpolation code
+   - x-data is duplicated, but that memory waste is
+     unlikely to represent a serious issue.
+
+  map<double, <double, double> >
+   - naturally couples x-value, f-value and d-value
+   - even more unreadable(??) code?
+   - higher skills needed by programmer.
+
+  
+  MONOTONE CUBIC INTERPOLATION:
+
+  It is a local algorithm. Adding one point only incur recomputation
+  of values in a neighbourhood of the point (in the interval getting
+  divided).
+
+  NOTE: We do not currently make use of this local fact. Upon
+  insertion of a new data pair, everything is recomputed. Revisit 
+  this when needed.
+
+*/
+
+
+
+
+MonotCubicInterpolator::
+MonotCubicInterpolator(const vector<double> & x, const vector<double> & f) {
+  if (x.size() != f.size()) {
+    throw("Unable to constuct MonotCubicInterpolator from vectors.") ;
+  }
+  
+  // Add the contents of the input vectors to our map of values.
+  vector<double>::const_iterator posx, posf;
+  for (posx = x.begin(), posf = f.begin(); posx != x.end(); ++posx, ++posf) {
+    data[*posx] = *posf ;
+  }
+
+  computeInternalFunctionData();
+}
+
+
+
+bool 
+MonotCubicInterpolator::
+read(const std::string & datafilename, int xColumn, int fColumn) 
+{
+  data.clear() ;
+  ddata.clear() ;
+
+  ifstream datafile_fs(datafilename.c_str());
+  if (!datafile_fs) {
+    return false ;
+  }
+  
+  string linestring;
+  while (!datafile_fs.eof()) {
+    getline(datafile_fs, linestring);
+
+    // Replace commas with space:
+    string::size_type pos = 0;
+    while ( (pos = linestring.find(",", pos)) != string::npos ) {
+        // cout << "Found comma at position " << pos << endl;
+        linestring.replace(pos, 1, " ");
+        pos++;
+    }
+
+    stringstream strs(linestring);
+    int columnindex = 0;
+    std::vector<double> value; 
+    if (linestring.size() > 0 && linestring.at(0) != '#') {
+        while (!(strs.rdstate() & std::ios::failbit)) {
+            double tmp;
+            strs >> tmp;
+            value.push_back(tmp);
+            columnindex++;
+        }
+    }
+    if (columnindex >= (max(xColumn, fColumn))) {
+      data[value[xColumn-1]] =  value[fColumn-1] ;
+    }
+  }   
+  datafile_fs.close();
+  
+  if (data.size() == 0) {
+    return false ;
+  }
+  
+  computeInternalFunctionData();
+  return true ;
+}
+
+
+void 
+MonotCubicInterpolator::
+addPair(double newx, double newf) throw(const char*) {
+  if (isnan(newx) || isinf(newx) || isnan(newf) || isinf(newf)) {
+    throw("MonotCubicInterpolator: addPair() received inf/nan input.");
+  }
+  data[newx] = newf ;
+  
+  // In a critical application, we should only update the
+  // internal function data for the offended interval, 
+  // not for all function values over again.
+  computeInternalFunctionData();
+}
+
+   
+double 
+MonotCubicInterpolator::
+evaluate(double x) const throw(const char*){
+
+  if (isnan(x) || isinf(x)) {
+    throw("MonotCubicInterpolator: evaluate() received inf/nan input.");
+  }
+
+  // xf becomes the first (xdata,fdata) pair where xdata >= x
+  map<double,double>::const_iterator xf_iterator = data.lower_bound(x);
+  
+  // First check if we must extrapolate:
+  if (xf_iterator == data.begin()) {
+    if (data.begin()->first == x) { // Just on the interval limit
+      return data.begin()->second;
+    }
+    else {
+      // Constant extrapolation (!!)
+      return data.begin()->second;
+    }
+  }
+  if (xf_iterator == data.end()) {
+      // Constant extrapolation (!!)
+      return data.rbegin()->second;
+  }
+  
+  
+  // Ok, we have x_min < x < x_max
+
+  pair<double,double> xf2 = *xf_iterator;
+  pair<double,double> xf1 = *(--xf_iterator);
+  // we now have: xf2.first > x 
+  
+  // Linear interpolation if derivative data is not available:
+  if (ddata.size() != data.size()) {
+    double finterp =  xf1.second + 
+      (xf2.second - xf1.second) / (xf2.first - xf1.first) 
+      * (x - xf1.first);
+    return finterp;
+  }
+  else { // Do Cubic Hermite spline
+    double t = (x - xf1.first)/(xf2.first - xf1.first); // t \in [0,1]
+    double h = xf2.first - xf1.first;
+    double finterp 
+      = xf1.second       * H00(t) 
+      + ddata[xf1.first] * H10(t) * h
+      + xf2.second       * H01(t) 
+      + ddata[xf2.first] * H11(t) * h ;
+    return finterp;
+  }
+  
+}
+
+
+// double 
+// MonotCubicInterpolator::
+// evaluate(double x, double& errorestimate_output) {
+//    cout << "Error: errorestimate not implemented" << endl;
+//    throw("error estimate not implemented");
+//    return x;
+// }
+
+vector<double> 
+MonotCubicInterpolator::
+get_xVector() const 
+{
+  map<double,double>::const_iterator xf_iterator = data.begin();
+  
+  vector<double> outputvector;
+  outputvector.reserve(data.size());
+  for (xf_iterator = data.begin(); xf_iterator != data.end(); ++xf_iterator) {
+    outputvector.push_back(xf_iterator->first);
+  }
+  return outputvector;
+}
+
+
+vector<double>
+MonotCubicInterpolator::
+get_fVector() const 
+{
+
+  map<double,double>::const_iterator xf_iterator = data.begin();
+  
+  vector<double> outputvector;
+  outputvector.reserve(data.size());
+  for (xf_iterator = data.begin(); xf_iterator != data.end(); ++xf_iterator) {
+    outputvector.push_back(xf_iterator->second);
+  }
+  return outputvector;
+}
+
+
+
+string
+MonotCubicInterpolator::
+toString() const 
+{
+  const int precision = 20;
+  std::string dataString;
+  std::stringstream dataStringStream;
+  for (map<double,double>::const_iterator it = data.begin();
+       it != data.end(); ++it) {
+      dataStringStream << setprecision(precision) << it->first;
+    dataStringStream << '\t';
+    dataStringStream << setprecision(precision) << it->second;
+    dataStringStream << '\n';
+  }
+  dataStringStream << "Derivative values:" << endl;
+  for (map<double,double>::const_iterator it = ddata.begin();
+       it != ddata.end(); ++it) {
+    dataStringStream << setprecision(precision) << it->first;
+    dataStringStream << '\t';
+    dataStringStream << setprecision(precision) << it->second;
+    dataStringStream << '\n';
+  }
+  
+  return dataStringStream.str();
+  
+}
+
+
+pair<double,double>
+MonotCubicInterpolator::
+getMissingX() const throw(const char*)
+{
+  if( data.size() < 2) {
+    throw("MonotCubicInterpolator::getMissingX() only one datapoint."); 
+  }
+  
+  // Search for biggest difference value in function-datavalues:
+  
+  map<double,double>::const_iterator  maxfDiffPair_iterator = data.begin();;
+  double maxfDiffValue = 0;
+  
+  map<double,double>::const_iterator xf_iterator;
+  map<double,double>::const_iterator xf_next_iterator;
+  
+  for (xf_iterator = data.begin(), xf_next_iterator = ++(data.begin());
+       xf_next_iterator != data.end(); 
+       ++xf_iterator, ++xf_next_iterator) {
+    double absfDiff = fabs((double)(*xf_next_iterator).second 
+                           - (double)(*xf_iterator).second);
+    if (absfDiff > maxfDiffValue) {
+      maxfDiffPair_iterator = xf_iterator;
+      maxfDiffValue = absfDiff;
+    }
+  }
+  
+  double newXvalue = ((*maxfDiffPair_iterator).first + ((*(++maxfDiffPair_iterator)).first))/2;
+  return make_pair(newXvalue, maxfDiffValue);
+  
+}
+
+
+
+pair<double,double>
+MonotCubicInterpolator::
+getMaximumF() const throw(const char*) {
+  if (data.size() <= 1) {
+    throw ("MonotCubicInterpolator::getMaximumF() empty data.") ;
+  }
+  if (strictlyIncreasing)
+    return *data.rbegin();
+  else if (strictlyDecreasing)
+    return *data.begin();
+  else {
+    pair<double,double> maxf = *data.rbegin() ;
+    map<double,double>::const_iterator xf_iterator;
+    for (xf_iterator = data.begin() ; xf_iterator != data.end(); ++xf_iterator) {
+      if (xf_iterator->second > maxf.second) {
+        maxf = *xf_iterator ;
+      } ;
+    }
+    return maxf ;
+  }
+}
+
+
+pair<double,double>
+MonotCubicInterpolator::
+getMinimumF() const throw(const char*) {
+  if (data.size() <= 1) {
+    throw ("MonotCubicInterpolator::getMinimumF() empty data.") ;
+  }
+  if (strictlyIncreasing)
+    return *data.begin();
+  else if (strictlyDecreasing) {
+    return *data.rbegin();
+  }
+  else {
+    pair<double,double> minf = *data.rbegin() ;
+    map<double,double>::const_iterator xf_iterator;
+    for (xf_iterator = data.begin() ; xf_iterator != data.end(); ++xf_iterator) {
+      if (xf_iterator->second < minf.second) {
+        minf = *xf_iterator ;
+      } ;
+    }
+    return minf ;
+  }
+}
+
+
+void 
+MonotCubicInterpolator::
+computeInternalFunctionData() const { 
+  
+  /* The contents of this function is meaningless if there is only one datapoint */
+  if (data.size() <= 1) {
+    return;
+  }
+  
+  /* We do not check the caching flag if we are instructed
+     to do this computation */
+  
+  /* We compute monotoneness and directions by assuming
+     monotoneness, and setting to false if the function is not for
+     some value */
+  
+  map<double,double>::const_iterator xf_iterator;
+  map<double,double>::const_iterator xf_next_iterator;
+  
+  
+  strictlyMonotone = true; // We assume this is true, and will set to false if not
+  monotone = true;
+  strictlyDecreasing = true;
+  decreasing = true;
+  strictlyIncreasing = true;
+  increasing = true;
+  
+  // Increasing or decreasing??
+  xf_iterator = data.begin();
+  xf_next_iterator = ++(data.begin());
+  /* Cater for non-strictness, search for direction for monotoneness */
+  while (xf_next_iterator != data.end() &&
+         xf_iterator->second == xf_next_iterator->second) {
+    /* Ok, equal values, this is not strict. */
+    strictlyMonotone = false;
+    strictlyIncreasing = false;
+    strictlyDecreasing = false;
+    
+    ++xf_iterator;
+    ++xf_next_iterator;
+  }
+  
+  
+  if (xf_next_iterator != data.end()) {
+
+    if ( xf_iterator->second > xf_next_iterator->second) {
+      // Ok, decreasing, check monotoneness:
+      strictlyDecreasing = true;// if strictlyMonotone == false, this one should not be trusted anyway
+      decreasing = true;
+      strictlyIncreasing = false;
+      increasing = false;
+      while(++xf_iterator, ++xf_next_iterator != data.end()) {
+        if ((*xf_iterator).second <  (*xf_next_iterator).second) {
+          monotone = false;
+          strictlyMonotone = false;
+          strictlyDecreasing = false; // meaningless now
+          break; // out of while loop
+        }
+        if ((*xf_iterator).second <= (*xf_next_iterator).second) {
+          strictlyMonotone = false;
+          strictlyDecreasing = false; // meaningless now
+        }
+      }
+    }
+    else if (xf_iterator->second < xf_next_iterator->second) {
+      // Ok, assume increasing, check monotoneness:
+      strictlyDecreasing = false;
+      strictlyIncreasing = true;
+      decreasing = false;
+      increasing = true;
+      while(++xf_iterator, ++xf_next_iterator != data.end()) {
+        if ((*xf_iterator).second >  (*xf_next_iterator).second) {
+          monotone = false;
+          strictlyMonotone = false;
+          strictlyIncreasing = false; // meaningless now
+          break; // out of while loop
+        }
+        if ((*xf_iterator).second >= (*xf_next_iterator).second) {
+          strictlyMonotone = false;
+          strictlyIncreasing = false; // meaningless now
+        }
+      }
+    }
+    else { 
+      // first two values must be equal if we get 
+      // here, but that should have been taken care
+      // of by the while loop above.
+      throw("Programming logic error.") ;
+    }
+
+  } 
+  computeSimpleDerivatives();
+  
+
+  // If our input data is monotone, we can do monotone cubic
+  // interpolation, so adjust the derivatives if so.
+  //
+  // If input data is not monotone, we should not touch 
+  // the derivatives, as this code should reduce to a
+  // standard cubic interpolation algorithm.
+  if (monotone) {
+    adjustDerivativesForMonotoneness();
+  }
+  
+  strictlyMonotoneCached = true;
+  monotoneCached = true;
+}
+
+//       Checks if the function curve is flat (zero derivative) at the
+//       endpoints, chop off endpoint data points if that is the case.
+//
+//       The notion of "flat" is determined by the input parameter "epsilon"
+//       Values whose difference are less than epsilon are regarded as equal.
+// 
+//       This is implemented to be able to obtain a strictly monotone
+//       curve from a data set that is strictly monotone except at the
+//       endpoints.
+//
+//       Example:
+//         The data points
+//            (1,3), (2,3), (3,4), (4,5), (5,5), (6,5)
+//         will become
+//            (2,3), (3,4), (4,5)
+// 
+//       Assumes at least 3 datapoints. If less than three, this function is a noop.
+void
+MonotCubicInterpolator::
+chopFlatEndpoints(const double epsilon) {
+
+    if (getSize() < 3) {
+        return;
+    }
+    
+    map<double,double>::iterator xf_iterator;
+    map<double,double>::iterator xf_next_iterator;
+
+    // Clear flags:
+    strictlyMonotoneCached = false;
+    monotoneCached = false;
+    
+    // Chop left end:
+    xf_iterator = data.begin();
+    xf_next_iterator = ++(data.begin());
+    // Erase data points that are similar to its right value from the left end.
+    while ((xf_next_iterator != data.end()) && 
+           (fabs(xf_iterator->second - xf_next_iterator->second) < epsilon )) {
+        xf_next_iterator++;
+        data.erase(xf_iterator);
+	xf_iterator++;
+    }
+
+    xf_iterator = data.end(); 
+    xf_iterator--;   // (data.end() points beyond the last element)
+    xf_next_iterator = xf_iterator;
+    xf_next_iterator--;
+    // Erase data points that are similar to its left value from the right end.
+    while ((xf_next_iterator != data.begin()) &&
+           (fabs(xf_iterator->second - xf_next_iterator->second) < epsilon )) {
+        xf_next_iterator--;
+        data.erase(xf_iterator);
+	xf_iterator--;
+    }
+    
+    // Finished chopping, so recompute function data:
+    computeInternalFunctionData();
+}
+
+
+//   
+//       If function is monotone, but not strictly monotone,
+//       this function will remove datapoints from intervals 
+//       with zero derivative so that the curves become
+//       strictly monotone.
+//
+//       Example
+//         The data points
+//           (1,2), (2,3), (3,4), (4,4), (5,5), (6,6)
+//         will become
+//           (1,2), (2,3), (3,4), (5,5), (6,6)
+//       
+//       Assumes at least two datapoints, if one or zero datapoint, this is a noop.
+//   
+//
+void
+MonotCubicInterpolator::
+shrinkFlatAreas(const double epsilon) {
+    
+    if (getSize() < 2) {
+        return;
+    }
+    
+    map<double,double>::iterator xf_iterator;
+    map<double,double>::iterator xf_next_iterator;
+    
+    
+    // Nothing to do if we already are strictly monotone
+    if (isStrictlyMonotone()) {
+        return;
+    }
+    
+    // Refuse to change a curve that is not monotone.
+    if (!isMonotone()) {
+        return;
+    }
+    
+    // Clear flags, they are not to be trusted after we modify the
+    // data
+    strictlyMonotoneCached = false;
+    monotoneCached = false;
+    
+    // Iterate through data values, if two data pairs
+    // have equal values, delete one of the data pair.
+    // Do not trust the source code on which data point is being
+    // removed (x-values of equal y-points might be averaged in the future)
+    xf_iterator = data.begin();
+    xf_next_iterator = ++(data.begin());
+    
+    while (xf_next_iterator != data.end()) {
+        //cout << xf_iterator->first << "," << xf_iterator->second << " " << xf_next_iterator->first << "," << xf_next_iterator->second << "\n";
+        if (fabs(xf_iterator->second - xf_next_iterator->second) < epsilon ) {
+            //cout << "erasing data pair" << xf_next_iterator->first << " " << xf_next_iterator->second << "\n";
+            map <double,double>::iterator xf_tobedeleted_iterator = xf_next_iterator;
+            xf_next_iterator++;
+            data.erase(xf_tobedeleted_iterator);
+        }
+        else {
+            xf_iterator++;
+            xf_next_iterator++;
+        }
+    }
+    
+}
+
+
+void 
+MonotCubicInterpolator::
+computeSimpleDerivatives() const {
+       
+  ddata.clear();
+  
+  // Do endpoints first:
+  map<double,double>::const_iterator xf_prev_iterator;
+  map<double,double>::const_iterator xf_iterator;
+  map<double,double>::const_iterator xf_next_iterator;
+  double diff;
+  
+  // Leftmost interval:
+  xf_iterator = data.begin();
+  xf_next_iterator = ++(data.begin());
+  diff = 
+    (xf_next_iterator->second - xf_iterator->second) /
+    (xf_next_iterator->first  - xf_iterator->first);
+  ddata[xf_iterator->first] = diff ;
+  
+  // Rightmost interval:
+  xf_iterator = --(--(data.end()));
+  xf_next_iterator = --(data.end());
+  diff = 
+    (xf_next_iterator->second - xf_iterator->second) /
+    (xf_next_iterator->first  - xf_iterator->first);
+  ddata[xf_next_iterator->first] = diff ;
+  
+  // If we have more than two intervals, loop over internal points:
+  if (data.size() > 2) {
+    
+    map<double,double>::const_iterator intpoint;
+    for (intpoint = ++data.begin(); intpoint != --data.end(); ++intpoint) {
+      /*
+        diff = (f2 - f1)/(x2-x1)/w + (f3-f1)/(x3-x2)/2 
+        
+        average of the forward and backward difference.
+        Weights are equal, should we weigh with h_i?
+      */
+      
+      map<double,double>::const_iterator lastpoint = intpoint; --lastpoint;
+      map<double,double>::const_iterator nextpoint = intpoint; ++nextpoint;
+      
+      diff = (nextpoint->second - intpoint->second)/
+        (2*(nextpoint->first - intpoint->first)) 
+        +
+        (intpoint->second - lastpoint->second) /
+        (2*(intpoint->first - lastpoint->first));
+
+      ddata[intpoint->first] = diff ;      
+    }
+  }
+}
+
+
+
+void 
+MonotCubicInterpolator::
+adjustDerivativesForMonotoneness() const {
+  map<double,double>::const_iterator point, dpoint;
+  
+  /* Loop over all intervals, ie. loop over all points and look
+     at the interval to the right of the point */
+  for (point = data.begin(), dpoint = ddata.begin();
+       point != --data.end();
+       ++point, ++dpoint) {
+    map<double,double>::const_iterator nextpoint, nextdpoint;
+    nextpoint = point; ++nextpoint;
+    nextdpoint = dpoint; ++nextdpoint;
+    
+    double delta = 
+      (nextpoint->second - point->second) /
+      (nextpoint->first  - point->first);
+    if (fabs(delta) < 1e-14) {
+      ddata[point->first] = 0.0;
+      ddata[nextpoint->first] = 0.0;
+    } else {
+      double alpha = ddata[point->first] / delta;
+      double beta = ddata[nextpoint->first] / delta;
+      
+      if (! isMonotoneCoeff(alpha, beta)) {
+        double tau = 3/sqrt(alpha*alpha + beta*beta);
+        
+        ddata[point->first]     = tau*alpha*delta;
+        ddata[nextpoint->first] = tau*beta*delta;
+      }
+    }
+    
+    
+  }
+  
+  
+}
+
+
+
+
+void 
+MonotCubicInterpolator::
+scaleData(double factor) {
+  map<double,double>::iterator it , itd  ;
+  if (data.size() == ddata.size()) {
+    for (it = data.begin() , itd = ddata.begin() ; it != data.end() ; ++it , ++itd) {
+      it->second  *= factor ;
+      itd->second *= factor ;
+    } ;
+  } else {
+    for (it = data.begin() ; it != data.end() ; ++it ) {
+      it->second  *= factor ;
+    } 
+  } 
+}

dune/common/MonotCubicInterpolator.hpp

@@ -0,0 +1,578 @@
+/* -*-C++-*- */
+
+#ifndef _MONOTCUBICINTERPOLATOR_H
+#define _MONOTCUBICINTERPOLATOR_H
+
+#include <vector>
+#include <map>
+#include <string>
+
+/*
+  MonotCubicInterpolator
+  Copyright (C) 2006 Statoil ASA
+  
+  This program 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 2
+  of the License, or (at your option) any later version.
+  
+  This program 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 this program; if not, write to the Free Software
+  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+*/
+
+/**
+   Class to represent a one-dimensional function f with single-valued
+   argument x. The function is represented by a table of function
+   values. Interpolation between table values is cubic and monotonicity 
+   preserving if input values are monotonous.
+
+   Outside x_min and x_max, the class will extrapolate using the
+   constant f(x_min) or f(x_max).
+   
+   Extra functionality:
+    - Can return (x_1+x_2)/2 where x_1 and x_2 are such that 
+      abs(f(x_1) - f(x_2)) is maximized. This is used to determine where 
+      one should calculate a new value for increased accuracy in the
+      current function
+
+   Monotonicity preserving cubic interpolation algorithm is taken 
+   from Fritsch and Carlson, "Monotone piecewise cubic interpolation", 
+   SIAM J. Numer. Anal. 17, 238--246, no. 2, 
+
+   $Id$
+
+   Algorithm also described here:
+   http://en.wikipedia.org/wiki/Monotone_cubic_interpolation
+
+
+   @author Håvard Berland <havb (at) statoil.com>, December 2006
+   @brief Represents one dimensional function f with single valued argument x that can be interpolated using monotone cubic interpolation
+
+*/
+
+class MonotCubicInterpolator {
+ public:
+   
+   /**
+      @param datafilename A datafile with the x values and the corresponding f(x) values 
+
+      Accepts a filename as input and parses this file for
+      two-column floating point data, interpreting the data as
+      representing function values x and f(x).
+
+      Ignores all lines not conforming to \<whitespace\>\<float\>\<whitespace\>\<float\>\<whatever\>\<newline\>
+   */
+  MonotCubicInterpolator(const std::string & datafilename) throw (const char*) 
+  {
+    if (!read(datafilename)) { 
+      throw("Unable to constuct MonotCubicInterpolator from file.") ; 
+    } ;
+  } ;
+
+
+   /**
+      @param datafilename A datafile with the x values and the corresponding f(x) values 
+
+      Accepts a filename as input and parses this file for
+      two-column floating point data, interpreting the data as
+      representing function values x and f(x).
+
+      Ignores all lines not conforming to \<whitespace\>\<float\>\<whitespace\>\<float\>\<whatever\>\<newline\>
+      
+      All commas in the file will be treated as spaces when parsing.
+
+   */
+ 
+ MonotCubicInterpolator(const char* datafilename) throw (const char*) 
+  {
+    if (!read(std::string(datafilename))) { 
+      throw("Unable to constuct MonotCubicInterpolator from file.") ; 
+    } ;
+  } ;
+
+
+   /** 
+       @param datafilename data file
+       @param XColumn x values
+       @param fColumn f values
+
+       Accepts a filename as input, and parses the chosen columns in
+       that file.
+   */
+   MonotCubicInterpolator(const char* datafilename, int xColumn, int fColumn) throw (const char*) 
+  {
+    if (!read(std::string(datafilename),xColumn,fColumn)) { 
+      throw("Unable to constuct MonotCubicInterpolator from file.") ; 
+    } ;
+  } ;
+   
+   /** 
+       @param datafilename data file
+       @param XColumn x values
+       @param fColumn f values
+
+       Accepts a filename as input, and parses the chosen columns in
+       that file.
+   */
+   MonotCubicInterpolator(const std::string & datafilename, int xColumn, int fColumn) throw (const char*) 
+  {
+    if (!read(datafilename,xColumn,fColumn)) { 
+      throw("Unable to constuct MonotCubicInterpolator from file.") ; 
+    } ;
+  } ;
+   
+   /**
+      @param x vector of x values
+      @param f vector of corresponding f values
+
+      Accepts two equal-length vectors as input for constructing
+      the interpolation object.  First vector is the x-values, the
+      second vector is the function values
+   */
+     MonotCubicInterpolator(const std::vector<double> & x , 
+                            const std::vector<double> & f);
+   
+   /**
+      No input, an empty function object is created.
+      
+      This object must be treated with care until
+      populated.
+   */
+   MonotCubicInterpolator() {;} ;
+
+
+
+   /**
+      @param datafilename A datafile with the x values and the corresponding f(x) values 
+
+      Accepts a filename as input and parses this file for
+      two-column floating point data, interpreting the data as
+      representing function values x and f(x).
+
+      returns true on success
+
+      All commas in file will be treated as spaces when parsing
+
+      Ignores all lines not conforming to \<whitespace\>\<float\>\<whitespace\>\<float\>\<whatever\>\<newline\>
+   */
+  bool read(const std::string & datafilename) {
+    return read(datafilename,1,2) ;
+  } ;
+   
+   /** 
+       @param datafilename data file
+       @param XColumn x values
+       @param fColumn f values
+
+       Accepts a filename as input, and parses the chosen columns in
+       that file.
+   */
+  bool read(const std::string &  datafilename, int xColumn, int fColumn) ; 
+
+
+
+   /**
+      @param x x value
+
+      Returns f(x) for given x (input). Interpolates (monotone cubic 
+      or linearly) if necessary.
+      
+      Extrapolates using the constants f(x_min) or f(x_max) if
+      input x is outside (x_min, x_max)
+
+      @return f(x) for a given x
+   */
+  double operator () (double x) const { return evaluate(x) ; } ;
+   
+   /**
+      @param x x value
+
+      Returns f(x) for given x (input). Interpolates (monotone cubic 
+      or linearly) if necessary.
+      
+      Extrapolates using the constants f(x_min) or f(x_max) if
+      input x is outside (x_min, x_max)
+
+      @return f(x) for a given x
+   */
+   double evaluate(double x) const throw(const char*);
+   
+   /**
+      @param x x value
+      @param errorestimate_output 
+
+      Returns f(x) and an error estimate for given x (input).
+      
+      Interpolates (linearly) if necessary.
+      
+      Throws an exception if extrapolation would be necessary for
+      evaluation. We do not want to do extrapolation (yet).
+      
+      The error estimate for x1 < x < x2 is
+      (x2 - x1)^2/8 * f''(x) where f''(x) is evaluated using
+      the stencil (1 -2  1) using either (x0, x1, x2) or (x1, x2, x3);
+      
+      Throws an exception if the table contains only two x-values.
+      
+      NOT IMPLEMENTED YET!
+   */
+   double evaluate(double x, double & errorestimate_output ) const ; 
+   
+   /**
+      Minimum x-value, returns both x and f in a pair.
+      
+      @return minimum x value
+      @return f(minimum x value)
+   */
+   std::pair<double,double> getMinimumX() const {
+       // Easy since the data is sorted on x:
+       return *data.begin();
+   }
+   
+   /**
+      Maximum x-value, returns both x and f in a pair.
+      
+      @return maximum x value
+      @return f(maximum x value)
+   */
+   std::pair<double,double> getMaximumX() const {
+       // Easy since the data is sorted on x:
+       return *data.rbegin();
+   }
+       
+   /**
+      Maximum f-value, returns both x and f in a pair.
+      
+      @return x value corresponding to maximum f value
+      @return maximum f value
+   */
+   std::pair<double,double> getMaximumF() const throw(const char*) ;
+
+   /**
+      Minimum f-value, returns both x and f in a pair
+      
+      @return x value corresponding to minimal f value
+      @return minimum f value
+   */
+   std::pair<double,double> getMinimumF() const throw(const char*)  ;
+
+   
+   /**
+      Provide a copy of the x-data as a vector
+      
+      Unspecified order, but corresponds to get_fVector.
+      
+      @return x values as a vector
+   */
+   std::vector<double> get_xVector() const ;
+   
+   /**
+      Provide a copy of tghe function data as a vector
+      
+      Unspecified order, but corresponds to get_xVector
+
+      @return f values as a vector
+      
+   */
+   std::vector<double> get_fVector() const ;
+   
+   /**
+      @param factor Scaling constant
+
+      Scale all the function value data by a constant
+   */
+   void scaleData(double factor);
+   
+   /**
+      Determines if the current function-value-data is strictly
+      monotone. This is a utility function for outsiders if they want
+      to invert the data for example.
+
+      @return True if f(x) is strictly monotone, else False
+   */
+   bool isStrictlyMonotone() {
+       
+       /* Use cached value if it can be trusted */
+       if (strictlyMonotoneCached) {
+           return strictlyMonotone;
+       }
+       else {
+           computeInternalFunctionData();
+           return strictlyMonotone;
+       }
+   }
+
+   /**
+      Determines if the current function-value-data is monotone.
+      
+      @return True if f(x) is monotone, else False
+   */
+   bool isMonotone() const {
+       if (monotoneCached) {
+           return monotone;
+       }
+       else {
+           computeInternalFunctionData();
+           return monotone;
+       }
+   }
+   /**
+      Determines if the current function-value-data is strictly
+      increasing. This is a utility function for outsiders if they want
+      to invert the data for example.
+
+      @return True if f(x) is strictly increasing, else False
+   */
+   bool isStrictlyIncreasing() {
+       
+       /* Use cached value if it can be trusted */
+       if (strictlyMonotoneCached) {
+           return (strictlyMonotone && strictlyIncreasing);
+       }
+       else {
+           computeInternalFunctionData();
+           return (strictlyMonotone && strictlyIncreasing);
+       }
+   }
+
+   /**
+      Determines if the current function-value-data is monotone and increasing.
+      
+      @return True if f(x) is monotone and increasing, else False
+   */
+   bool isMonotoneIncreasing() const {
+       if (monotoneCached) {
+           return (monotone && increasing);
+       }
+       else {
+           computeInternalFunctionData();
+           return (monotone && increasing);
+       }
+   }
+   /**
+      Determines if the current function-value-data is strictly
+      decreasing. This is a utility function for outsiders if they want
+      to invert the data for example.
+
+      @return True if f(x) is strictly decreasing, else False
+   */
+   bool isStrictlyDecreasing() {
+       
+       /* Use cached value if it can be trusted */
+       if (strictlyMonotoneCached) {
+           return (strictlyMonotone && strictlyDecreasing);
+       }
+       else {
+           computeInternalFunctionData();
+           return (strictlyMonotone && strictlyDecreasing);
+       }
+   }
+
+   /**
+      Determines if the current function-value-data is monotone and decreasing
+      
+      @return True if f(x) is monotone and decreasing, else False
+   */
+   bool isMonotoneDecreasing() const {
+       if (monotoneCached) {
+           return (monotone && decreasing);
+       }
+       else {
+           computeInternalFunctionData();
+           return (monotone && decreasing);
+       }
+   }
+    
+
+   
+   /**
+      @param newx New x point
+      @param newf New f(x) point
+
+      Adds a new datapoint to the function.
+
+      This causes all the derivatives at all points of the functions
+      to be recomputed and then adjusted for monotone cubic
+      interpolation. If this function ever enters a critical part of
+      any code, the locality of the algorithm for monotone adjustment
+      must be exploited.
+
+   */
+   void addPair(double newx, double newf) throw(const char*);
+   
+   /**
+      Returns an x-value that is believed to yield the best
+      improvement in global accuracy for the interpolation if
+      computed.
+      
+      Searches for the largest jump in f-values, and returns a x
+      value being the average of the two x-values representing the
+      f-value-jump.
+
+      @return New x value beleived to yield the best improvement in global accuracy
+      @return Maximal difference
+   */
+   std::pair<double,double> getMissingX() const throw(const char*) ;
+   
+   /**
+      Constructs a string containing the data in a table 
+
+      @return a string containing the data in a table
+   */
+   std::string toString() const;
+   
+   /**
+     @return Number of datapoint pairs in this object
+   */
+   int getSize() const {
+       return data.size();
+   }
+
+    /**
+       Checks if the function curve is flat at the endpoints, chop off
+       endpoint data points if that is the case.
+      
+       The notion of "flat" is determined by the input parameter "epsilon"
+       Values whose difference are less than epsilon are regarded as equal.
+
+       This is implemented to be able to obtain a strictly monotone
+       curve from a data set that is strictly monotone except at the
+       endpoints.
+
+       Example:
+         The data points
+            (1,3), (2,3), (3,4), (4,5), (5,5), (6,5)
+         will become
+            (2,3), (3,4), (4,5)
+
+       Assumes at least 3 datapoints. If less than three, this function is a noop.
+    */
+    void chopFlatEndpoints(const double);
+
+    /**
+       Wrapper function for chopFlatEndpoints(const double)
+       providing a default epsilon parameter 
+    */
+    void chopFlatEndpoints() {
+        chopFlatEndpoints(1e-14);
+    }
+
+    /**
+       If function is monotone, but not strictly monotone,
+       this function will remove datapoints from intervals 
+       with zero derivative so that the curve become
+       strictly monotone.
+
+       Example
+         The data points
+           (1,2), (2,3), (3,4), (4,4), (5,5), (6,6)
+         will become
+           (1,2), (2,3), (3,4), (5,5), (6,6)
+       
+       Assumes at least two datapoints, if one or zero datapoint, this is a noop.
+    */
+    void shrinkFlatAreas(const double);
+
+    /**
+       Wrapper function for shrinkFlatAreas(const double)
+       providing a default epsilon parameter 
+    */
+    void shrinkFlatAreas() {
+        shrinkFlatAreas(1e-14);
+    };
+
+
+
+private:
+   
+   // Data structure to store x- and f-values
+   std::map<double, double> data;
+   
+   // Data structure to store x- and d-values
+   mutable std::map<double, double> ddata;  
+   
+   
+   // Storage containers for precomputed interpolation data
+   //   std::vector<double> dvalues; // derivatives in Hermite interpolation.
+   
+   // Flag to determine whether the boolean strictlyMonotone can be
+   // trusted.
+   mutable bool strictlyMonotoneCached;
+   mutable bool monotoneCached; /* only monotone, not stricly montone */
+
+   mutable bool strictlyMonotone;
+   mutable bool monotone;
+
+   // if strictlyMonotone is true (and can be trusted), the two next are meaningful
+   mutable bool strictlyDecreasing;
+   mutable bool strictlyIncreasing;
+   mutable bool decreasing;
+   mutable bool increasing;
+   
+
+   /* Hermite basis functions, t \in [0,1] ,
+      notation from:
+      http://en.wikipedia.org/w/index.php?title=Cubic_Hermite_spline&oldid=84495502
+   */
+   
+   double H00(double t) const {
+       return 2*t*t*t - 3*t*t + 1;
+   }
+   double H10(double t) const {
+       return t*t*t - 2*t*t + t;
+   }
+   double H01(double t) const {
+       return -2*t*t*t + 3*t*t;
+   }
+   double H11(double t) const {
+       return t*t*t - t*t;
+   }
+   
+   
+   void computeInternalFunctionData() const ;
+   
+   /** 
+       Computes initial derivative values using centered (second order) difference
+       for internal datapoints, and one-sided derivative for endpoints
+       
+       The internal datastructure map<double,double> ddata is populated by this method.
+   */
+   
+   void computeSimpleDerivatives() const ;
+
+   
+   /**
+      Adjusts the derivative values (ddata) so that we can guarantee that
+      the resulting piecewise Hermite polymial is monotone. This is
+      done according to the algorithm of Fritsch and Carlsson 1980, 
+      see Section 4, especially the two last lines.
+   */
+  void adjustDerivativesForMonotoneness() const ;
+
+   /** 
+       Checks if the coefficient alpha and beta is in 
+       the region that guarantees monotoneness of the 
+       derivative values they represent 
+       
+       See Fritsch and Carlson 1980, Lemma 2, 
+       alternatively Step 5 in Wikipedia's article 
+       on Monotone cubic interpolation.
+   */
+   bool isMonotoneCoeff(double alpha, double beta) const {
+       if ((alpha*alpha + beta*beta) <= 9) {
+         return true;
+       } else {
+         return false;
+       }
+   }
+   
+}; 
+
+
+#endif

dune/common/RootFinders.hpp

@@ -0,0 +1,178 @@
+//===========================================================================
+//
+// File: RootFinders.hpp
+//
+// Created: Thu May  6 19:59:42 2010
+//
+// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
+//            Jostein R Natvig    <jostein.r.natvig@sintef.no>
+//
+// $Date$
+//
+// $Revision$
+//
+//===========================================================================
+
+/*
+  Copyright 2010 SINTEF ICT, Applied Mathematics.
+  Copyright 2010 Statoil ASA.
+
+  This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+  OpenRS is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OpenRS is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPENRS_ROOTFINDERS_HEADER
+#define OPENRS_ROOTFINDERS_HEADER
+
+
+
+//#include <algorithm>
+//#include <limits>
+//#include <cmath>
+//#include <cstdlib>
+#include <dune/common/ErrorMacros.hpp>
+
+namespace Dune
+{
+
+	inline double regulaFalsiStep(const double a,
+                                      const double b,
+                                      const double fa,
+                                      const double fb)
+	{
+	    ASSERT(fa*fb < 0.0);
+	    return (b*fa - a*fb)/(fa - fb);
+	}
+
+
+	/// Implements a modified regula falsi method as described in
+	/// "Improved algorithms of Illinois-type for the numerical
+	/// solution of nonlinear equations"
+	/// by J. A. Ford.
+	/// Current variant is the 'Pegasus' method.
+	template <class Functor>
+	inline double modifiedRegulaFalsi(const Functor& f,
+                                          const double a,
+                                          const double b,
+                                          const int max_iter,
+                                          const double tolerance,
+                                          int& iterations_used)
+	{
+	    using namespace std;
+	    const double macheps = numeric_limits<double>::epsilon();
+	    const double eps = tolerance + macheps*max(max(fabs(a), fabs(b)), 1.0);
+
+	    double x0 = a;
+	    double x1 = b;
+	    double f0 = f(x0);
+	    const double epsF = tolerance + macheps*max(fabs(f0), 1.0);
+	    if (fabs(f0) < epsF) {
+		return x0;
+	    }
+	    double f1 = f(x1);
+	    if (fabs(f1) < epsF) {
+		return x1;
+	    }
+	    if (f0*f1 > 0.0) {
+		THROW("Error in parameters, zero not bracketed: [a, b] = ["
+		      << a << ", " << b << "]    fa = " << f0 << "   fb = " << f1);
+	    }
+	    iterations_used = 0;
+	    // In every iteraton, x1 is the last point computed,
+	    // and x0 is the last point computed that makes it a bracket.
+	    while (fabs(x1 - x0) >= 0.95*eps) {
+		double xnew = regulaFalsiStep(x0, x1, f0, f1);
+		double fnew = f(xnew);
+// 		cout << "xnew = " << xnew << "    fnew = " << fnew << endl;
+		++iterations_used;
+		if (iterations_used > max_iter) {
+		    THROW("Maximum number of iterations exceeded.\n"
+			  << "Current interval is [" << min(x0, x1) << ", "
+			  << max(x0, x1) << "]");
+		}
+		if (fabs(fnew) < epsF) {
+		    return xnew;
+		}
+		// Now we must check which point we must replace.
+		if ((fnew > 0.0) == (f0 > 0.0)) {
+		    // We must replace x0.
+		    x0 = x1;
+		    f0 = f1;
+		} else {
+		    // We must replace x1, this is the case where
+		    // the modification to regula falsi kicks in,
+		    // by modifying f0.
+		    // 1. The classic Illinois method
+//  		    const double gamma = 0.5;
+		    // @afr: The next two methods do not work??!!?
+		    // 2. The method called 'Method 3' in the paper.
+// 		    const double phi0 = f1/f0;
+// 		    const double phi1 = fnew/f1;
+// 		    const double gamma = 1.0 - phi1/(1.0 - phi0);
+		    // 3. The method called 'Method 4' in the paper.
+// 		    const double phi0 = f1/f0;
+// 		    const double phi1 = fnew/f1;
+// 		    const double gamma = 1.0 - phi0 - phi1;
+// 		    cout << "phi0 = " << phi0 <<" phi1 = " << phi1 <<
+// 		    " gamma = " << gamma << endl;
+		    // 4. The 'Pegasus' method
+ 		    const double gamma = f1/(f1 + fnew);
+		    f0 *= gamma;
+		}
+		x1 = xnew;
+		f1 = fnew;
+	    }
+	    return 0.5*(x0 + x1);
+	}
+
+
+	template <class Functor>
+	inline void bracketZero(const Functor& f,
+                                const double x0,
+                                const double dx,
+                                double& a,
+                                double& b)
+        {
+            const int max_iters = 100;
+            double f0 = f(x0);
+            double cur_dx = dx;
+            int i = 0;
+            for (; i < max_iters; ++i) {
+                double x = x0 + cur_dx;
+                double f_new = f(x);
+                if (f0*f_new <= 0.0) {
+                    break;
+                }
+                cur_dx = -2.0*cur_dx;
+            }
+            if (i == max_iters) {
+                THROW("Could not bracket zero in " << max_iters << "iterations.");
+            }
+            if (cur_dx < 0.0) {
+                a = x0 + cur_dx;
+                b = i < 2 ? x0 : x0 + 0.25*cur_dx;
+            } else {
+                a = i < 2 ? x0 : x0 + 0.25*cur_dx;
+                b = x0 + cur_dx;
+            }
+        }
+
+
+} // namespace Dune
+
+
+
+
+#endif // OPENRS_ROOTFINDERS_HEADER

dune/common/SparseTable.hpp

@@ -0,0 +1,239 @@
+//===========================================================================
+//
+// File: SparseTable.hpp
+//
+// Created: Fri Apr 24 09:50:27 2009
+//
+// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
+//
+// $Date$
+//
+// $Revision$
+//
+//===========================================================================
+
+/*
+  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+  Copyright 2009, 2010 Statoil ASA.
+
+  This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+  OpenRS is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OpenRS is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPENRS_SPARSETABLE_HEADER
+#define OPENRS_SPARSETABLE_HEADER
+
+#include <vector>
+#include <numeric>
+#include <algorithm>
+#include <boost/range/iterator_range.hpp>
+#include "ErrorMacros.hpp"
+
+#include <ostream>
+
+namespace Dune
+{
+
+    /// A SparseTable stores a table with rows of varying size
+    /// as efficiently as possible.
+    /// It is supposed to behave similarly to a vector of vectors.
+    /// Its behaviour is similar to compressed row sparse matrices.
+    template <typename T>
+    class SparseTable
+    {
+    public:
+        /// Default constructor. Yields an empty SparseTable.
+        SparseTable()
+        {
+        }
+
+        /// A constructor taking all the data for the table and row sizes.
+        /// \param data_beg The start of the table data.
+        /// \param data_end One-beyond-end of the table data.
+        /// \param rowsize_beg The start of the row length data.
+        /// \param rowsize_end One beyond the end of the row length data.
+        template <typename DataIter, typename IntegerIter>
+        SparseTable(DataIter data_beg, DataIter data_end,
+                    IntegerIter rowsize_beg, IntegerIter rowsize_end)
+            : data_(data_beg, data_end)
+        {
+	    setRowStartsFromSizes(rowsize_beg, rowsize_end);
+        }
+
+
+        /// Sets the table to contain the given data, organized into
+	/// rows as indicated by the given row sizes.
+        /// \param data_beg The start of the table data.
+        /// \param data_end One-beyond-end of the table data.
+        /// \param rowsize_beg The start of the row length data.
+        /// \param rowsize_end One beyond the end of the row length data.
+        template <typename DataIter, typename IntegerIter>
+        void assign(DataIter data_beg, DataIter data_end,
+                    IntegerIter rowsize_beg, IntegerIter rowsize_end)
+        {
+	    data_.assign(data_beg, data_end);
+	    setRowStartsFromSizes(rowsize_beg, rowsize_end);
+        }
+
+
+        /// Request storage for table of given size.
+        /// \param rowsize_beg Start of row size data.
+        /// \param rowsize_end One beyond end of row size data.
+        template <typename IntegerIter>
+        void allocate(IntegerIter rowsize_beg, IntegerIter rowsize_end)
+        {
+            typedef typename std::vector<T>::size_type sz_t;
+
+            sz_t ndata = std::accumulate(rowsize_beg, rowsize_end, sz_t(0));
+            data_.resize(ndata);
+            setRowStartsFromSizes(rowsize_beg, rowsize_end);
+        }
+
+
+        /// Appends a row to the table.
+        template <typename DataIter>
+        void appendRow(DataIter row_beg, DataIter row_end)
+        {
+            data_.insert(data_.end(), row_beg, row_end);
+            if (row_start_.empty()) {
+                row_start_.reserve(2);
+                row_start_.push_back(0);
+            }
+            row_start_.push_back(data_.size());
+        }
+
+        /// True if the table contains no rows.
+        bool empty() const
+        {
+            return row_start_.empty();
+        }
+
+        /// Returns the number of rows in the table.
+        int size() const
+        {
+            return empty() ? 0 : row_start_.size() - 1;
+        }
+
+        /// Allocate storage for table of expected size
+        void reserve(int exptd_nrows, int exptd_ndata)
+        {
+            row_start_.reserve(exptd_nrows + 1);
+            data_.reserve(exptd_ndata);
+        }
+
+        /// Swap contents for other SparseTable<T>
+        void swap(SparseTable<T>& other)
+        {
+            row_start_.swap(other.row_start_);
+            data_.swap(other.data_);
+        }
+
+        /// Returns the number of data elements.
+        int dataSize() const
+        {
+            return data_.size();
+        }
+
+        /// Returns the size of a table row.
+        int rowSize(int row) const
+        {
+            ASSERT(row >= 0 && row < size());
+            return row_start_[row + 1] - row_start_[row];
+        }
+
+        /// Makes the table empty().
+        void clear()
+        {
+            data_.clear();
+            row_start_.clear();
+        }
+
+        /// Defining the row type, returned by operator[].
+        typedef boost::iterator_range<const T*> row_type;
+        typedef boost::iterator_range<T*>       mutable_row_type;
+
+        /// Returns a row of the table.
+        row_type operator[](int row) const
+        {
+            ASSERT(row >= 0 && row < size());
+            const T* start_ptr = data_.empty() ? 0 : &data_[0];
+            return row_type(start_ptr + row_start_[row], start_ptr + row_start_[row + 1]);
+        }
+
+        /// Returns a mutable row of the table.
+        mutable_row_type operator[](int row)
+        {
+            ASSERT(row >= 0 && row < size());
+            T* start_ptr = data_.empty() ? 0 : &data_[0];
+            return mutable_row_type(start_ptr + row_start_[row], start_ptr + row_start_[row + 1]);
+        }
+
+        /// Equality.
+        bool operator==(const SparseTable& other) const
+        {
+            return data_ == other.data_ && row_start_ == other.row_start_;
+        }
+
+        template<class charT, class traits>
+        void print(std::basic_ostream<charT, traits>& os) const
+        {
+            os << "Number of rows: " << size() << '\n';
+
+            os << "Row starts = [";
+            std::copy(row_start_.begin(), row_start_.end(),
+                      std::ostream_iterator<int>(os, " "));
+            os << "\b]\n";
+
+            os << "Data values = [";
+            std::copy(data_.begin(), data_.end(),
+                      std::ostream_iterator<T>(os, " "));
+            os << "\b]\n";
+        }
+
+    private:
+        std::vector<T> data_;
+        // Like in the compressed row sparse matrix format,
+        // row_start_.size() is equal to the number of rows + 1.
+        std::vector<int> row_start_;
+
+	template <class IntegerIter>
+	void setRowStartsFromSizes(IntegerIter rowsize_beg, IntegerIter rowsize_end)
+	{
+            // Since we do not store the row sizes, but cumulative row sizes,
+            // we have to create the cumulative ones.
+            int num_rows = rowsize_end - rowsize_beg;
+            if (num_rows < 1) {
+                THROW("Must have at least one row. Got " << num_rows << " rows.");
+            }
+#ifndef NDEBUG
+            if (*std::min_element(rowsize_beg, rowsize_end) < 0) {
+                THROW("All row sizes must be at least 0.");
+            }
+#endif
+            row_start_.resize(num_rows + 1);
+            row_start_[0] = 0;
+            std::partial_sum(rowsize_beg, rowsize_end, row_start_.begin() + 1);
+            // Check that data_ and row_start_ match.
+            if (int(data_.size()) != row_start_.back()) {
+                THROW("End of row start indices different from data size.");
+            }
+
+	}
+    };
+
+} // namespace Dune
+
+
+#endif // OPENRS_SPARSETABLE_HEADER

dune/common/SparseVector.hpp

@@ -0,0 +1,195 @@
+//===========================================================================
+//
+// File: SparseVector.hpp
+//
+// Created: Mon Jun 29 15:28:59 2009
+//
+// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
+//            Bård Skaflestad     <bard.skaflestad@sintef.no>
+//
+// $Date$
+//
+// $Revision$
+//
+//===========================================================================
+
+/*
+  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+  Copyright 2009, 2010 Statoil ASA.
+
+  This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+  OpenRS is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OpenRS is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+#ifndef OPENRS_SPARSEVECTOR_HEADER
+#define OPENRS_SPARSEVECTOR_HEADER
+
+#include <vector>
+#include <numeric>
+#include <algorithm>
+#include <boost/range/iterator_range.hpp>
+#include "ErrorMacros.hpp"
+
+namespace Dune
+{
+
+    /// A SparseVector stores a vector with possibly many empty elements
+    /// as efficiently as possible.
+    /// It is supposed to behave similarly to a standard vector, but since
+    /// direct indexing is a O(log n) operation instead of O(1), we do not
+    /// supply it as operator[].
+    template <typename T>
+    class SparseVector
+    {
+    public:
+	/// Default constructor. Yields an empty SparseVector.
+	SparseVector()
+	    : size_(0), default_elem_()
+	{
+	}
+
+	/// Constructs a SparseVector with a given size, but no nonzero
+	/// elements.
+	explicit SparseVector(int size)
+	    : size_(size), default_elem_()
+	{
+	}
+
+	/// A constructor taking all the element data for the vector and their indices.
+	/// \param data_beg The start of the element data.
+	/// \param data_end One-beyond-end of the element data.
+	/// \param rowsize_beg The start of the index data.
+	/// \param rowsize_end One beyond the end of the index data.
+	template <typename DataIter, typename IntegerIter>
+	SparseVector(int size,
+		     DataIter data_beg, DataIter data_end,
+		     IntegerIter index_beg, IntegerIter index_end)
+	    : size_(size), data_(data_beg, data_end), indices_(index_beg, index_end),
+	      default_elem_()
+	{
+#ifndef NDEBUG
+	    ASSERT(size >= 0);
+	    ASSERT(indices_.size() == data_.size());
+	    int last_index = -1;
+	    int num_ind = indices_.size();
+	    for (int i = 0; i < num_ind; ++i) {
+		int index = indices_[i];
+		if (index <= last_index || index >= size) {
+		    THROW("Error in SparseVector construction, index is nonincreasing or out of range.");
+		}
+		last_index = index;
+	    }
+#endif
+	}
+
+
+	/// Appends an element to the vector. Note that this function does not
+	/// increase the size() of the vector, it just adds another nonzero element.
+	/// Elements must be added in index order.
+	void addElement(const T& element, int index)
+	{
+	    ASSERT(indices_.empty() || index > indices_.back());
+	    ASSERT(index < size_);
+	    data_.push_back(element);
+	    indices_.push_back(index);
+	}
+
+	/// \return true if the vector has size 0.
+	bool empty() const
+	{
+	    return size_ == 0;
+	}
+
+	/// Returns the size of the vector.
+	/// Recall that most or all of the vector may be default/zero.
+	int size() const
+	{
+	    return size_;
+	}
+
+	/// Returns the number of nonzero data elements.
+	int nonzeroSize() const
+	{
+	    return data_.size();
+	}
+
+	/// Makes the vector empty().
+	void clear()
+	{
+	    data_.clear();
+	    indices_.clear();
+	    size_ = 0;
+	}
+
+	/// Equality.
+	bool operator==(const SparseVector& other) const
+	{
+	    return size_ == other.size_ && data_ == other.data_ && indices_ == other.indices_;
+	}
+
+	/// O(log n) element access.
+	/// \param index the proper vector index
+	/// \return the element with the given index, or the default element if no element in
+	/// the vector has the given index.
+	const T& element(int index) const
+	{
+	    ASSERT(index >= 0);
+	    ASSERT(index < size_);
+	    std::vector<int>::const_iterator lb = std::lower_bound(indices_.begin(), indices_.end(), index);
+	    if (lb != indices_.end() && *lb == index) {
+		return data_[lb - indices_.begin()];
+	    } else {
+		return default_elem_;
+	    }
+	}
+
+	/// O(1) element access.
+	/// \param nzindex an index counting only nonzero elements.
+	/// \return the nzindex'th nonzero element.
+	const T& nonzeroElement(int nzindex) const
+	{
+	    ASSERT(nzindex >= 0);
+	    ASSERT(nzindex < nonzeroSize());
+	    return data_[nzindex];
+	}
+
+	/// O(1) index access.
+	/// \param nzindex an index counting only nonzero elements.
+	/// \return the index of the nzindex'th nonzero element.
+	int nonzeroIndex(int nzindex) const
+	{
+	    ASSERT(nzindex >= 0);
+	    ASSERT(nzindex < nonzeroSize());
+	    return indices_[nzindex];
+	}
+
+    private:
+	// The vectors data_ and indices_ are always the same size.
+	// The indices are supposed to be stored in increasing order,
+	// to be unique, and to be in [0, size_ - 1].
+	// default_elem_ is returned when a default element is requested.
+	int size_;
+	std::vector<T> data_;
+	std::vector<int> indices_;
+	T default_elem_;
+    };
+
+} // namespace Dune
+
+
+
+
+#endif // OPENRS_SPARSEVECTOR_HEADER

dune/common/StopWatch.cpp

@@ -0,0 +1,106 @@
+//===========================================================================
+//
+// File: StopWatch.cpp
+//
+// Created: Thu Jul  2 23:04:51 2009
+//
+// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
+//
+// $Date$
+//
+// $Revision$
+//
+//===========================================================================
+
+/*
+  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+  Copyright 2009, 2010 Statoil ASA.
+
+  This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+  OpenRS is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OpenRS is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#if HAVE_CONFIG_H
+#include "config.h"
+#endif
+#include <boost/date_time/posix_time/posix_time.hpp>
+#include "StopWatch.hpp"
+#include "ErrorMacros.hpp"
+
+namespace Dune
+{
+
+    namespace time
+    {
+
+	StopWatch::StopWatch()
+	    : state_(UnStarted)
+	{
+	}
+
+
+	void StopWatch::start()
+	{
+	    start_time_ = boost::posix_time::microsec_clock::local_time();
+	    last_time_ = start_time_;
+	    state_ = Running;
+	}
+
+	void StopWatch::stop()
+	{
+	    if (state_ != Running) {
+		THROW("Called stop() on a StopWatch that was not running.");
+	    }
+	    stop_time_ = boost::posix_time::microsec_clock::local_time();
+	    state_ = Stopped;
+	}
+
+	double StopWatch::secsSinceLast()
+	{
+	    boost::posix_time::ptime run_time;
+	    if (state_ == Running) {
+		run_time = boost::posix_time::microsec_clock::local_time();
+	    } else if (state_ == Stopped) {
+		run_time = stop_time_;
+	    } else {
+		ASSERT(state_ == UnStarted);
+		THROW("Called secsSinceLast() on a StopWatch that had not been started.");
+	    }
+	    boost::posix_time::time_duration dur = run_time - last_time_;
+	    last_time_ = run_time;
+	    return double(dur.total_microseconds())/1000000.0;
+	}
+
+	double StopWatch::secsSinceStart()
+	{
+	    boost::posix_time::ptime run_time;
+	    if (state_ == Running) {
+		run_time = boost::posix_time::microsec_clock::local_time();
+	    } else if (state_ == Stopped) {
+		run_time = stop_time_;
+	    } else {
+		ASSERT(state_ == UnStarted);
+		THROW("Called secsSinceStart() on a StopWatch that had not been started.");
+	    }
+	    boost::posix_time::time_duration dur = run_time - start_time_;
+	    return double(dur.total_microseconds())/1000000.0;
+	}
+
+    } // namespace time
+
+} // namespace Dune
+
+
+

dune/common/StopWatch.hpp

@@ -0,0 +1,81 @@
+//===========================================================================
+//
+// File: StopWatch.hpp
+//
+// Created: Thu Jul  2 23:04:17 2009
+//
+// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
+//
+// $Date$
+//
+// $Revision$
+//
+//===========================================================================
+
+/*
+  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+  Copyright 2009, 2010 Statoil ASA.
+
+  This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+  OpenRS is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OpenRS is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPENRS_STOPWATCH_HEADER
+#define OPENRS_STOPWATCH_HEADER
+
+#include <boost/date_time/posix_time/posix_time.hpp>
+
+namespace Dune
+{
+
+    namespace time
+    {
+
+	class StopWatch
+	{
+	public:
+	    /// Default constructor. Before the StopWatch is start()-ed,
+	    /// it is an error to call anything other than start().
+	    StopWatch();
+
+	    /// Starts the StopWatch. It is always legal to call
+	    /// start(), even if not stop()-ped.
+	    void start();
+	    /// Stops the StopWatch. The watch no longer runs, until
+	    /// restarted by a call to start().
+	    void stop();
+
+	    /// \ret the number of running seconds that have passed
+	    /// since last call to start(), secsSinceLast() or
+	    /// secsSinceStart()
+	    double secsSinceLast();
+	    /// \ret the number of running seconds that have passed
+	    /// since last call to start().
+	    double secsSinceStart();
+
+	private:
+	    enum StopWatchState { UnStarted, Running, Stopped };
+
+	    StopWatchState state_;
+	    boost::posix_time::ptime start_time_;
+	    boost::posix_time::ptime last_time_;
+	    boost::posix_time::ptime stop_time_;
+	};
+
+    } // namespace time
+
+} // namespace Dune
+
+#endif // OPENRS_STOPWATCH_HEADER

dune/common/Units.hpp

@@ -0,0 +1,182 @@
+//===========================================================================
+//
+// File: Units.hpp
+//
+// Created: Thu Jul  2 09:19:08 2009
+//
+// Author(s): Halvor M Nilsen <hnil@sintef.no>
+//
+// $Date$
+//
+// $Revision$
+//
+//===========================================================================
+
+/*
+  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+  Copyright 2009, 2010 Statoil ASA.
+
+  This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+  OpenRS is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OpenRS is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPENRS_UNITS_HEADER
+#define OPENRS_UNITS_HEADER
+
+namespace Dune
+{
+    namespace prefix
+    {
+        const double micro = 1.0e-6;
+        const double milli = 1.0e-3;
+        const double centi = 1.0e-2;
+        const double deci  = 1.0e-1;
+        const double kilo  = 1.0e3;
+        const double mega  = 1.0e6;
+        const double giga  = 1.0e9;
+    } // namespace prefix
+
+    namespace unit
+    {
+        // Common powers
+        inline double square(double v) { return v * v;     }
+        inline double cubic (double v) { return v * v * v; }
+
+        // --------------------------------------------------------------
+        // Basic (fundamental) units and conversions
+        // --------------------------------------------------------------
+
+        // Length:
+        const double meter =  1;
+        const double inch  =  2.54 * prefix::centi*meter;
+        const double feet  = 12    * inch;
+
+        // Time:
+        const double second =   1;
+        const double minute =  60 * second;
+        const double hour   =  60 * minute;
+        const double day    =  24 * hour;
+        const double year   = 365 * day;
+
+        // Volume
+        const double stb    = 0.158987294928 * cubic(meter);
+
+
+        // Mass:
+        const double kilogram = 1;
+
+        // http://en.wikipedia.org/wiki/Pound_(mass)#Avoirdupois_pound
+        const double pound    = 0.45359237 * kilogram;
+
+        // --------------------------------------------------------------
+        // Standardised constants
+        // --------------------------------------------------------------
+
+        const double gravity = 9.80665 * meter/square(second);
+
+        // --------------------------------------------------------------
+        // Derived units and conversions
+        // --------------------------------------------------------------
+
+        // Force:
+        const double Newton = kilogram*meter / square(second); // == 1
+        const double lbf    = pound * gravity; // Pound-force
+
+        // Pressure:
+        const double Pascal = Newton / square(meter); // == 1
+        const double barsa  = 100000 * Pascal;
+        const double atm    = 101325 * Pascal;
+        const double psia   = lbf / square(inch);
+
+        // Viscosity:
+        const double Pas   = Pascal * second; // == 1
+        const double Poise = prefix::deci*Pas;
+
+        // Permeability:
+        //
+        // A porous medium with a permeability of 1 darcy permits a
+        // flow (flux) of 1 cm³/s of a fluid with viscosity 1 cP (1
+        // mPa·s) under a pressure gradient of 1 atm/cm acting across
+        // an area of 1 cm².
+        //
+        namespace perm_details {
+            const double p_grad   = atm / (prefix::centi*meter);
+            const double area     = square(prefix::centi*meter);
+            const double flux     = cubic (prefix::centi*meter) / second;
+            const double velocity = flux / area;
+            const double visc     = prefix::centi*Poise;
+            const double darcy    = (velocity * visc) / p_grad;
+            //                    == 1e-7 [m^2] / 101325
+            //                    == 9.869232667160130e-13 [m^2]
+        }
+        const double darcy = perm_details::darcy;
+
+        // Unit conversion support.
+        //
+        // Note: Under the penalty of treason will you be
+        //
+        //    using namespace Dune::unit::convert;
+        //
+        // I mean it!
+        //
+        namespace convert {
+            // Convert from external units of measurements to equivalent
+            // internal units of measurements.  Note: The internal units
+            // of measurements are *ALWAYS*, and exclusively, SI.
+            //
+            // Example: Convert a double kx, containing a permeability
+            // value in units of milli-darcy (mD) to the equivalent
+            // value in SI units (m^2).
+            //
+            //    using namespace Dune::unit;
+            //    using namespace Dune::prefix;
+            //    convert::from(kx, milli*darcy);
+            //
+            inline double from(const double q, const double unit)
+            {
+                return q * unit;
+            }
+
+            // Convert from internal units of measurements to equivalent
+            // external units of measurements.  Note: The internal units
+            // of measurements are *ALWAYS*, and exclusively, SI.
+            //
+            // Example: Convert a std::vector<double> p, containing
+            // pressure values in the SI unit Pascal (i.e., unit::Pascal)
+            // to the equivalent values in Psi (unit::psia).
+            //
+            //    using namespace Dune::unit;
+            //    std::transform(p.begin(), p.end(), p.begin(),
+            //                   boost::bind(convert::to, _1, psia));
+            //
+            inline double to(const double q, const double unit)
+            {
+                return q / unit;
+            }
+        } // namespace convert
+    } // namespace unit
+
+    namespace units {
+        //     const double MILLIDARCY = 1.0;//9.86923e-16;
+        //     const double VISCOSITY_UNIT = 1.0;//1e-3;
+        //     const double DAYS2SECONDS = 1.0;//86400;
+        const double MILLIDARCY = 9.86923e-16;
+        const double VISCOSITY_UNIT = 1e-3;
+        const double DAYS2SECONDS = 86400;
+        const double FEET = 0.30479999798832;
+        const double WELL_INDEX_UNIT = VISCOSITY_UNIT/(DAYS2SECONDS*1e5);
+    } // namespace units
+} // namespace Dune
+#endif // OPENRS_UNITS_HEADER

dune/common/linInt.hpp

@@ -0,0 +1,115 @@
+//===========================================================================
+//
+// File: linInt.hpp
+//
+// Created: Tue Feb 16 14:44:10 2010
+//
+// Author(s): Bjørn Spjelkavik    <bsp@sintef.no>
+//            Atgeirr F Rasmussen <atgeirr@sintef.no>
+//            Bård Skaflestad     <bard.skaflestad@sintef.no>
+//
+// $Date$
+//
+// $Revision$
+//
+//===========================================================================
+
+/*
+  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+  Copyright 2009, 2010 Statoil ASA.
+
+  This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+  OpenRS is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OpenRS is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPENRS_LININT_HEADER
+#define OPENRS_LININT_HEADER
+
+
+namespace Dune
+{
+
+    inline int tableIndex(const std::vector<double>& table, double x)
+    {
+	// Returns an index in an ordered table such that x is between
+	// table[j] and table[j+1]. If x is out of range, first or last
+	// interval is returned; Binary search.
+	int n = table.size() - 1;
+	if (n < 2) {
+	    return 0;
+	}
+	int jl = 0;
+	int ju = n;
+	bool ascend = (table[n] > table[0]);
+	while (ju - jl > 1) {
+	    int jm = (ju + jl)/2;   // Compute a midpoint
+	    if ( (x >= table[jm]) == ascend) {
+		jl = jm;     // Replace lower limit
+	    } else {
+		ju = jm;     // Replace upper limit
+	    }
+	}
+	return jl;
+    }
+
+    inline double linearInterpolation(const std::vector<double>& xv,
+                                      const std::vector<double>& yv, double x)
+    {
+	// Returns end point if x is outside xv
+	std::vector<double>::const_iterator lb = lower_bound(xv.begin(), xv.end(), x);
+	int ix2 = lb - xv.begin();
+	if (ix2 == 0) {
+	    return yv[0];
+	} else if (ix2 == int(xv.size())) {
+	    return yv[ix2-1];
+	}
+	int ix1 = ix2 - 1;
+	return  (yv[ix2] - yv[ix1])/(xv[ix2] - xv[ix1])*(x - xv[ix1]) + yv[ix1];
+    }
+
+    inline double linearInterpolDerivative(const std::vector<double>& xv,
+                                           const std::vector<double>& yv, double x)
+    {
+	int ix1 = tableIndex(xv, x);
+	int ix2 = ix1 + 1;
+	return  (yv[ix2] - yv[ix1])/(xv[ix2] - xv[ix1]);
+    }
+
+    inline double linearInterpolationExtrap(const std::vector<double>& xv,
+                                            const std::vector<double>& yv, double x)
+    {
+	// Extrapolates if x is outside xv
+	int ix1 = tableIndex(xv, x);
+	int ix2 = ix1 + 1;
+	return  (yv[ix2] - yv[ix1])/(xv[ix2] - xv[ix1])*(x - xv[ix1]) + yv[ix1];
+    }
+    
+    inline double linearInterpolationExtrap(const std::vector<double>& xv,
+                                            const std::vector<double>& yv,
+                                            double x, int& ix1)
+    {
+	// Extrapolates if x is outside xv
+	ix1 = tableIndex(xv, x);
+	int ix2 = ix1 + 1;
+	return (yv[ix2] - yv[ix1])/(xv[ix2] - xv[ix1])*(x - xv[ix1]) + yv[ix1];
+    }
+
+} // namespace Dune
+
+
+
+
+
+#endif // OPENRS_LININT_HEADER

dune/common/param/Makefile.am

@@ -0,0 +1,21 @@
+# $Date$
+# $Revision$
+
+SUBDIRS = . test
+
+paramdir = $(includedir)/dune/common/param
+
+param_HEADERS = Parameter.hpp ParameterGroup.hpp ParameterGroup_impl.hpp \
+                ParameterMapItem.hpp ParameterRequirement.hpp \
+                ParameterStrings.hpp ParameterTools.hpp ParameterXML.hpp
+
+noinst_LTLIBRARIES = libparam.la
+
+libparam_la_SOURCES = Parameter.cpp ParameterGroup.cpp ParameterTools.cpp \
+                      ParameterXML.cpp
+
+libparam_la_CPPFLAGS = $(XML_CPPFLAGS) $(BOOST_CPPFLAGS)
+libparam_la_LDFLAGS  = $(BOOST_LDFLAGS)
+libparam_la_LIBADD   = $(XML_LIBS) $(BOOST_FILESYSTEM_LIB) $(BOOST_SYSTEM_LIB)
+
+include $(top_srcdir)/am/global-rules

dune/common/param/Parameter.cpp

@@ -0,0 +1,74 @@
+//===========================================================================
+//
+// File: Parameter.cpp
+//
+// Created: Tue Jun  2 19:18:25 2009
+//
+// Author(s): Bård Skaflestad     <bard.skaflestad@sintef.no>
+//            Atgeirr F Rasmussen <atgeirr@sintef.no>
+//
+// $Date$
+//
+// $Revision$
+//
+//===========================================================================
+
+/*
+Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+Copyright 2009, 2010 Statoil ASA.
+
+This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+OpenRS is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+OpenRS is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#if HAVE_CONFIG_H
+#include "config.h"
+#endif
+#include <string>
+#include <dune/common/param/Parameter.hpp>
+
+namespace Dune {
+    namespace parameter {
+	std::string
+        correct_parameter_tag(const ParameterMapItem& item)
+        {
+	    std::string tag = item.getTag();
+	    if (tag != ID_xmltag__param) {
+		std::string error = "The XML tag was '" +
+                                    tag + "' but should be '" +
+                                    ID_xmltag__param + "'.\n";
+		return error;
+	    } else {
+		return "";
+	    }
+	}
+
+	std::string
+        correct_type(const Parameter& parameter,
+                     const std::string& param_type)
+        {
+	    std::string type = parameter.getType();
+	    if ( (type != param_type) &&
+                 (type != ID_param_type__cmdline) ) {
+		std::string error = "The data was of type '" + type +
+                                    "' but should be of type '" +
+                                    param_type + "'.\n";
+		return error;
+	    } else {
+		return "";
+	    }
+	}
+    } // namespace parameter
+} // namespace Dune

dune/common/param/Parameter.hpp

@@ -0,0 +1,214 @@
+//===========================================================================
+//
+// File: Parameter.hpp
+//
+// Created: Tue Jun  2 16:00:21 2009
+//
+// Author(s): Bård Skaflestad     <bard.skaflestad@sintef.no>
+//            Atgeirr F Rasmussen <atgeirr@sintef.no>
+//
+// $Date$
+//
+// $Revision$
+//
+//===========================================================================
+
+/*
+Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+Copyright 2009, 2010 Statoil ASA.
+
+This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+OpenRS is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+OpenRS is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPENRS_PARAMETER_HEADER
+#define OPENRS_PARAMETER_HEADER
+
+#include <string>
+#include <sstream>
+
+#include <dune/common/param/ParameterMapItem.hpp>
+#include <dune/common/param/ParameterStrings.hpp>
+
+namespace Dune {
+    /// See ParameterGroup.hpp for how to use the parameter system
+    namespace parameter {
+	/// @brief
+	/// @todo Doc me!
+	class Parameter : public ParameterMapItem {
+	public:
+	    /// @brief
+	    /// @todo Doc me! 
+	    virtual ~Parameter() {}
+	    /// @brief
+	    /// @todo Doc me!
+	    /// @return
+	    virtual std::string getTag() const {return ID_xmltag__param;}
+	    /// @brief
+	    /// @todo Doc me!
+	    /// @param
+	    Parameter(const std::string& value, const std::string& type)
+                : value_(value), type_(type) {}
+	    /// @brief
+	    /// @todo Doc me!
+	    /// @return
+	    std::string getValue() const {return value_;}
+	    /// @brief
+	    /// @todo Doc me!
+	    /// @return
+	    std::string getType() const {return type_;}
+	private:
+	    std::string value_;
+	    std::string type_;
+	};
+
+	/// @brief
+	/// @todo Doc me!
+	/// @param
+	/// @return
+	std::string correct_parameter_tag(const ParameterMapItem& item);
+	std::string correct_type(const Parameter& parameter,
+                                 const std::string& type);
+
+	/// @brief
+	/// @todo Doc me!
+	/// @tparam
+	/// @param
+	/// @return
+	template<>
+	struct ParameterMapItemTrait<int> {
+	    static int convert(const ParameterMapItem& item,
+                               std::string& conversion_error)
+            {
+		conversion_error = correct_parameter_tag(item);
+		if (conversion_error != "") {
+		    return 0;
+		}
+		const Parameter& parameter = dynamic_cast<const Parameter&>(item);
+		conversion_error = correct_type(parameter, ID_param_type__int);
+		if (conversion_error != "") {
+		    return 0;
+		}
+		std::stringstream stream;
+		stream << parameter.getValue();
+		int value;
+		stream >> value;
+		if (stream.fail()) {
+		    conversion_error = "Conversion to '" +
+                                       ID_param_type__int +
+                                       "' failed. Data was '" +
+                                       parameter.getValue() + "'.\n";
+		    return 0;
+		}
+		return value;
+	    }
+	    static std::string type() {return ID_param_type__int;}
+	};
+
+	/// @brief
+	/// @todo Doc me!
+	/// @tparam
+	/// @param
+	/// @return
+	template<>
+	struct ParameterMapItemTrait<double> {
+	    static double convert(const ParameterMapItem& item,
+                                  std::string& conversion_error)
+            {
+		conversion_error = correct_parameter_tag(item);
+		if (conversion_error != "") {
+		    return 0.0;
+		}
+		const Parameter& parameter = dynamic_cast<const Parameter&>(item);
+		conversion_error = correct_type(parameter, ID_param_type__float);
+		if (conversion_error != "") {
+		    return 0.0;
+		}
+		std::stringstream stream;
+		stream << parameter.getValue();
+		double value;
+		stream >> value;
+		if (stream.fail()) {
+		    conversion_error = "Conversion to '" +
+                                       ID_param_type__float +
+                                       "' failed. Data was '" +
+                                       parameter.getValue() + "'.\n";
+		    return 0.0;
+		}
+		return value;
+	    }
+	    static std::string type() {return ID_param_type__float;}
+	};
+
+	/// @brief
+	/// @todo Doc me!
+	/// @tparam
+	/// @param
+	/// @return
+	template<>
+	struct ParameterMapItemTrait<bool> {
+	    static bool convert(const ParameterMapItem& item,
+                                std::string& conversion_error)
+            {
+		conversion_error = correct_parameter_tag(item);
+		if (conversion_error != "") {
+		    return false;
+		}
+		const Parameter& parameter = dynamic_cast<const Parameter&>(item);
+		conversion_error = correct_type(parameter, ID_param_type__bool);
+		if (conversion_error != "") {
+		    return false;
+		}
+		if (parameter.getValue() == ID_true) {
+		    return true;
+		} else if (parameter.getValue() == ID_false) {
+		    return false;
+		} else {
+		    conversion_error = "Conversion failed. Data was '" +
+                                       parameter.getValue() +
+                                       "', but should be one of '" +
+                                       ID_true + "' or '" + ID_false + "'.\n";
+		    return false;
+		}
+	    }
+	    static std::string type() {return ID_param_type__bool;}
+	};
+
+	/// @brief
+	/// @todo Doc me!
+	/// @tparam
+	/// @param
+	/// @return
+	template<>
+	struct ParameterMapItemTrait<std::string> {
+	    static std::string convert(const ParameterMapItem& item,
+                                       std::string& conversion_error)
+            {
+		conversion_error = correct_parameter_tag(item);
+		if (conversion_error != "") {
+		    return "";
+		}
+		const Parameter& parameter = dynamic_cast<const Parameter&>(item);
+		conversion_error = correct_type(parameter, ID_param_type__string);
+		if (conversion_error != "") {
+		    return "";
+		}
+		return parameter.getValue();
+	    }
+	    static std::string type() {return ID_param_type__string;}
+	};
+    } // namespace parameter
+} // namespace Dune
+#endif  // OPENRS_PARAMETER_HPP

dune/common/param/ParameterGroup.cpp

@@ -0,0 +1,334 @@
+//===========================================================================
+//
+// File: ParameterGroup.cpp
+//
+// Created: Tue Jun  2 19:13:17 2009
+//
+// Author(s): Bård Skaflestad     <bard.skaflestad@sintef.no>
+//            Atgeirr F Rasmussen <atgeirr@sintef.no>
+//
+// $Date$
+//
+// $Revision$
+//
+//===========================================================================
+
+/*
+Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+Copyright 2009, 2010 Statoil ASA.
+
+This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+OpenRS is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+OpenRS is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#if HAVE_CONFIG_H
+#include "config.h"
+#endif
+#include <dune/common/param/ParameterGroup.hpp>
+
+#include <cassert>
+#include <cstdlib>
+#include <fstream>
+#include <iostream>
+
+#include <stdio.h>
+#include <string.h>
+
+//#include "TermColors.hpp" // from utils
+
+#include <dune/common/param/Parameter.hpp>
+#include <dune/common/param/ParameterStrings.hpp>
+#include <dune/common/param/ParameterTools.hpp>
+#include <dune/common/param/ParameterXML.hpp>
+
+namespace Dune {
+    namespace parameter {
+
+	ParameterGroup::ParameterGroup()
+	: path_(ID_path_root), parent_(0), output_is_enabled_(true)
+	{
+	}
+
+	ParameterGroup::~ParameterGroup() {
+	    if (output_is_enabled_) {
+		//TermColors::Normal();
+	    }
+	}
+
+	std::string ParameterGroup::getTag() const {
+	    return ID_xmltag__param_grp;
+	}
+
+
+	bool ParameterGroup::has(const std::string& name) const
+        {
+	    std::pair<std::string, std::string> name_path = split(name);
+	    map_type::const_iterator it = map_.find(name_path.first);
+	    if (it == map_.end()) {
+		return false;
+	    } else if (name_path.second == "") {
+		return true;
+	    } else {
+		ParameterGroup& pg = dynamic_cast<ParameterGroup&>(*(*it).second);
+		return pg.has(name_path.second);
+	    }
+	}
+
+	std::string ParameterGroup::path() const {
+	    return path_;
+	}
+
+	ParameterGroup::ParameterGroup(const std::string& path,
+				       const ParameterGroup* parent)
+	: path_(path), parent_(parent), output_is_enabled_(true)
+	{
+	}
+
+	ParameterGroup
+	ParameterGroup::getGroup(const std::string& name) const
+        {
+	    return get<ParameterGroup>(name);
+	}
+
+	void ParameterGroup::readXML(const std::string& xml_filename)
+        {
+	    fill_xml(*this, xml_filename);
+	}
+
+	namespace {
+	    inline std::istream&
+	    samcode_readline(std::istream& is, std::string& parameter) {
+		return std::getline(is, parameter);
+	    }
+	} // anonymous namespace
+
+	void ParameterGroup::readParam(const std::string& param_filename) {
+	    std::ifstream is(param_filename.c_str());
+	    if (!is) {
+		std::cerr << "Error: Failed to open parameter file '" << param_filename << "'.\n";
+		throw std::exception();
+	    }
+	    std::string parameter;
+	    int lineno = 0;
+	    while (samcode_readline(is, parameter)) {
+		++lineno;
+		int fpos = parameter.find(ID_delimiter_assignment);
+		if (fpos == int(std::string::npos)) {
+		    std::cerr << "WARNING: No '" << ID_delimiter_assignment << "' found on line " << lineno << ".\n";
+		}
+		int pos = fpos + ID_delimiter_assignment.size();
+		int spos = parameter.find(ID_delimiter_assignment, pos);
+		if (spos == int(std::string::npos)) {
+		    std::string name = parameter.substr(0, fpos);
+		    std::string value = parameter.substr(pos, spos);
+		    this->insertParameter(name, value);
+		} else {
+		    std::cerr << "WARNING: To many '" << ID_delimiter_assignment << "' found on line " << lineno << ".\n";
+		}
+	    }
+	    #ifdef MATLAB_MEX_FILE
+	    fclose(is);
+	    #endif
+	}
+
+	void ParameterGroup::writeParam(const std::string& param_filename) const {
+	    std::ofstream os(param_filename.c_str());
+	    if (!os) {
+		std::cerr << "Error: Failed to open parameter file '"
+                          << param_filename << "'.\n";
+		throw std::exception();
+	    }
+	    this->writeParamToStream(os);
+	}
+
+	void ParameterGroup::writeParamToStream(std::ostream& os) const
+        {
+	    if (map_.empty()) {
+		os << this->path() << "/" << "dummy"
+		   << ID_delimiter_assignment << "0\n";
+	    }
+	    for (map_type::const_iterator it = map_.begin(); it != map_.end(); ++it) {
+		if ( (*it).second->getTag() == ID_xmltag__param_grp ) {
+		    ParameterGroup& pg = dynamic_cast<ParameterGroup&>(*(*it).second);
+		    pg.writeParamToStream(os);
+		} else if ( (*it).second->getTag() == ID_xmltag__param) {
+		    os << this->path() << "/" << (*it).first << ID_delimiter_assignment
+		       << dynamic_cast<Parameter&>(*(*it).second).getValue() << "\n";
+		} else {
+		    os << this->path() << "/" << (*it).first << ID_delimiter_assignment
+		       << "<" << (*it).second->getTag() << ">" << "\n";
+		}
+	    }
+	}
+
+	void ParameterGroup::insert(const std::string& name,
+				    const std::tr1::shared_ptr<ParameterMapItem>& data)
+        {
+	    std::pair<std::string, std::string> name_path = split(name);
+	    map_type::const_iterator it = map_.find(name_path.first);
+	    assert(name_path.second == "");
+	    if (it == map_.end()) {
+		map_[name] = data;
+	    } else {
+		if ( (map_[name]->getTag() == data->getTag())  &&
+                     (data->getTag() == ID_xmltag__param_grp) ) {
+		    ParameterGroup& alpha = dynamic_cast<ParameterGroup&>(*(*it).second);
+		    ParameterGroup& beta  = dynamic_cast<ParameterGroup&>(*data);
+		    for (map_type::const_iterator
+                             item = beta.map_.begin(); item != beta.map_.end(); ++item) {
+			alpha.insert((*item).first, (*item).second);
+		    }
+		} else {
+		    std::cout << "WARNING : The '"
+			      << map_[name]->getTag()
+			      << "' element '"
+			      << name
+			      << "' already exist in group '"
+			      << this->path()
+			      << "'. The element will be replaced by a '"
+			      << data->getTag()
+			      << "' element.\n";
+		    map_[name] = data;
+		}
+	    }
+	}
+
+	void ParameterGroup::insertParameter(const std::string& name,
+                                             const std::string& value)
+        {
+	    std::pair<std::string, std::string> name_path = split(name);
+	    while (name_path.first == "") {
+		name_path = split(name_path.second);
+	    }
+	    map_type::const_iterator it = map_.find(name_path.first);
+	    if (it == map_.end()) {
+		if (name_path.second == "") {
+		    std::tr1::shared_ptr<ParameterMapItem> data(new Parameter(value, ID_param_type__cmdline));
+		    map_[name_path.first] = data;
+		} else {
+		    std::tr1::shared_ptr<ParameterMapItem> data(new ParameterGroup(this->path() + ID_delimiter_path + name_path.first,
+										   this));
+		    ParameterGroup& child = dynamic_cast<ParameterGroup&>(*data);
+		    child.insertParameter(name_path.second, value);
+		    map_[name_path.first] = data;
+		}
+	    } else if (name_path.second == "") {
+		if ((*it).second->getTag() != ID_xmltag__param) {
+		    std::cout << "WARNING : The "
+			      << map_[name_path.first]->getTag()
+			      << " element '"
+			      << name
+			      << "' already exist in group '"
+			      << this->path()
+			      << "'. It will be replaced by a "
+			      << ID_xmltag__param
+			      << " element.\n";
+		}
+		std::tr1::shared_ptr<ParameterMapItem> data(new Parameter(value, ID_param_type__cmdline));
+		map_[name_path.first] = data;
+	    } else {
+		ParameterGroup& pg = dynamic_cast<ParameterGroup&>(*(*it).second);
+		pg.insertParameter(name_path.second, value);
+	    }
+	}
+
+#if 0
+	void ParameterGroup::display() const {
+	    std::cout << "Begin: " << this->path() << "\n";
+	    for (map_type::const_iterator it = map_.begin(); it != map_.end(); ++it) {
+		if ( (*it).second->getTag() == ID_xmltag__param_grp ) {
+		    ParameterGroup& pg = dynamic_cast<ParameterGroup&>(*(*it).second);
+		    pg.display();
+		} else if ( (*it).second->getTag() == ID_xmltag__param) {
+		    std::cout << (*it).first
+			      << " ("
+			      << (*it).second->getTag()
+			      << ") has value "
+			      << dynamic_cast<Parameter&>(*(*it).second).getValue()
+			      << "\n";
+		} else {
+		    std::cout << (*it).first
+			      << " ("
+			      << (*it).second->getTag()
+			      << ")\n";
+		}
+	    }
+	    std::cout << "End: " << this->path() << "\n";
+	}
+#endif
+
+        bool ParameterGroup::anyUnused() const
+        {
+            if (!this->used()) {
+                return true;
+            }
+	    for (map_type::const_iterator it = map_.begin(); it != map_.end(); ++it) {
+		if (it->second->getTag() == ID_xmltag__param_grp) {
+		    ParameterGroup& pg = dynamic_cast<ParameterGroup&>(*(it->second));
+		    if (pg.anyUnused()) {
+                        return true;
+                    }
+		} else if (it->second->getTag() == ID_xmltag__param) {
+		    if (!it->second->used()) {
+			return true;
+		    }
+		}
+	    }
+            return false;
+        }
+
+	void ParameterGroup::displayUsage(bool used_params) const
+	{
+	    if (this->used() == used_params) {
+		std::cout << this->path() << '\n';
+	    }
+	    for (map_type::const_iterator it = map_.begin(); it != map_.end(); ++it) {
+		if (it->second->getTag() == ID_xmltag__param_grp) {
+		    ParameterGroup& pg = dynamic_cast<ParameterGroup&>(*(it->second));
+		    pg.displayUsage(used_params);
+		} else if (it->second->getTag() == ID_xmltag__param) {
+		    if (it->second->used() == used_params) {
+			std::cout << path() << '/' << it->first << '\n';
+		    }
+		}
+	    }
+	    std::cout << std::flush;
+	}
+
+	void ParameterGroup::disableOutput() {
+	    this->recursiveSetIsOutputEnabled(false);
+	}
+
+	void ParameterGroup::enableOutput() {
+	    this->recursiveSetIsOutputEnabled(true);
+	}
+
+	bool ParameterGroup::isOutputEnabled() const {
+	    return output_is_enabled_;
+	}
+
+	void ParameterGroup::recursiveSetIsOutputEnabled(bool output_is_enabled)
+        {
+	    output_is_enabled_ = output_is_enabled;
+	    for (map_type::const_iterator it = map_.begin(); it != map_.end(); ++it) {
+		if (it->second->getTag() == ID_xmltag__param_grp) {
+		    ParameterGroup& pg = dynamic_cast<ParameterGroup&>(*(it->second));
+		    pg.recursiveSetIsOutputEnabled(output_is_enabled);
+		}
+	    }
+	}
+
+    } // namespace parameter
+} // namespace Dune

dune/common/param/ParameterGroup.hpp

@@ -0,0 +1,292 @@
+//===========================================================================
+//
+// File: ParameterGroup.hpp
+//
+// Created: Tue Jun  2 19:11:11 2009
+//
+// Author(s): Bård Skaflestad     <bard.skaflestad@sintef.no>
+//            Atgeirr F Rasmussen <atgeirr@sintef.no>
+//
+// $Date$
+//
+// $Revision$
+//
+//===========================================================================
+
+/*
+Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+Copyright 2009, 2010 Statoil ASA.
+
+This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+OpenRS is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+OpenRS is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPENRS_PARAMETERGROUP_HEADER
+#define OPENRS_PARAMETERGROUP_HEADER
+
+#include <exception>
+#include <map>
+#include <string>
+#include <vector>
+#include <tr1/memory>
+
+#include <dune/common/param/ParameterMapItem.hpp>
+#include <dune/common/param/ParameterRequirement.hpp>
+
+namespace Dune {
+    namespace parameter {
+	/// ParameterGroup is a class that is used to provide run-time paramters.
+	/// The standard use of the class is to call create it with the
+	/// (int argc, char** argv) constructor (where the arguments are those
+	/// given by main). This parses the command line, where each token
+	/// either
+	/// A) specifies a parameter (by a "param=value" token).
+	/// B) specifies a xml file to be read (by a "filename.xml" token).
+	/// C) specifies a param file to be read (by a "filename.param" token).
+	/// After the tokens are parsem they are stored in a tree structure
+	/// in the ParameterGroup object; it is worth mentioning that parameters
+	/// are inherited in this tree structure. Thus, if ``grid\_prefix'' is
+	/// given a value in the root node, this value will be visible in all
+	/// parts of the tree (unless the parameter is overwritten in a subtree).
+	/// Applications using this ParameterGroup objects will usually write out
+	/// a message for each node in the three that is used by the application;
+	/// this is one way to determine valid parameters.
+	///
+	/// Parameters specified as "param=value" on the command line
+	///
+	/// To specify a parameter on the command line, you must know where in the tree the
+	/// parameter resides. The syntax for specifying parameters on the command line given
+	/// an application called ``simulator'' is
+	///         simulator param1=value grp/param2=value
+	//  for parameter ``param1'' lying at the root and ``param2'' in the group
+	/// ``grp''. If the same parameters are specified multiple times on the command
+	/// line, only the last will be used. Thus an application named ``simulator'' run with
+	/// the following command
+	///         simulator param=value1 param=value2
+	/// will use value2 as the value of ``param''.
+	///
+	/// XML parameters
+	///
+	/// In the introduction to this section it was mentioned that the parameters for
+	/// the simulator are organized in a tree structure. This is mirrored in the XML
+	/// file by the use of groups; a group is a collection of parameters and subgroups
+	/// enclosed by a \fixed{<ParameterGroup>}-\fixed{</ParameterGroup>} pair. The only
+	/// attribute of relevance in a \fixed{ParameterGroup} tag is that of \fixed{name},
+	/// which is used to navigate in the tree.
+	/// The actual parameters are specified by the \fixed{<Parameter />} tag. Each
+	/// parameter has three attributes: \fixed{name}, \fixed{type}, and \fixed{value}.
+	/// Both name and value should be evident. Type is one of \fixed{bool} (for things
+	/// that are either true or false), \fixed{int} (for integers), \fixed{double} (for
+	/// floating point numbers with double precision), \fixed{string} (for text
+	/// strings), or \fixed{file} (for files and directories relative to the location
+	/// of the XML file).
+	///
+	/// param files
+	///
+	/// A param file consists of multiple lienes, where each line consists of "param=value".
+	/// This syntax is identical to the one for paramters specified on the command line.
+	///
+	///
+	/// If one combines both XML files and parameters, one should note that that if a parameter
+	/// is specified on the command line is also found in a XML file, the parameter
+	/// specified on the command line is the one used.  Thus, if ``parameters.xml''
+	/// specifies that ``stepsize'' is 2.71828 while, the command used to run the
+	/// application ``simulator'' is
+	///         simulator stepsize=3.14159 parameters.xml
+	/// the simulator will run with ``stepsize'' equal to 3.14159.
+	///
+	class ParameterGroup : public ParameterMapItem {
+	public:
+	    struct NotFoundException          : public std::exception {};
+	    struct WrongTypeException         : public std::exception {};
+	    struct ConversionFailedException  : public std::exception {};
+
+	    template<typename T>
+	    struct RequirementFailedException : public std::exception {};
+
+	    ParameterGroup();
+	    ParameterGroup(const std::string& path, const ParameterGroup* parent);
+
+	    // From ParameterMapItem
+	    virtual ~ParameterGroup();
+	    virtual std::string getTag() const;
+
+	    /// \brief A constructor typically used to initialize a
+	    ///        ParameterGroup from command-line arguments.
+	    ///
+	    /// It is required that argv[0] is the program name, while if
+	    /// 0 < i < argc, then argv[i] is either
+	    /// the name of an xml file or parametername=value.
+	    ///
+	    /// \param argc is the number of command-line arguments,
+	    ///        including the name of the executable.
+	    /// \param argv is an array of char*, each of which is a
+	    /// command line argument.
+            template <typename StringArray>
+	    ParameterGroup(int argc, StringArray argv);
+
+	    /// \brief This method checks if there is something with name
+	    ///        \p name in the parameter gropup.
+	    ///
+	    /// \param name is the name of the parameter.
+	    /// \return true if \p name is the name of something in the parameter
+	    ///         group, false otherwise.
+	    bool has(const std::string& name) const;
+
+	    /// \brief This method is used to read a parameter from the
+	    ///        parameter group.
+	    ///
+	    /// NOTE:  If the reading of the parameter fails, then this method
+	    ///        throws an appropriate exception.
+	    ///
+	    /// \param name is the name of the parameter in question.
+	    /// \return The value associated with then name in this parameter
+	    ///         group.
+	    template<typename T>
+	    T get(const std::string& name) const;
+
+	    template<typename T, class Requirement>
+	    T get(const std::string& name, const Requirement&) const;
+
+	    /// \brief This method is used to read a parameter from the
+	    ///        parameter group.
+	    ///
+	    /// NOTE:  If the reading of the parameter fails, then either
+	    ///        a) this method returns \p default_value if there
+	    ///           was no parameter with name \p name in
+	    ///           the parameter group
+	    ///        or
+	    ///        b) this method throws an appropriate exception.
+	    ///
+	    /// \param name is the name of the parameter in question.
+	    /// \param default_value the defualt value of the parameter in
+	    ///        question.
+	    /// \return The value associated with then name in this parameter
+	    ///         group.
+	    template<typename T>
+	    T getDefault(const std::string& name,
+			 const T& default_value) const;
+
+	    template<typename T, class Requirement>
+	    T getDefault(const std::string& name,
+			 const T& default_value,
+			 const Requirement& r) const;
+
+	    /// \brief This method returns the parameter group given by name,
+	    ///        i.e. it is an alias of get<ParameterGroup>().
+	    ///
+	    /// \param name is the name of the parameter group sought.
+	    /// \return the parameter group sought.
+	    ParameterGroup getGroup(const std::string& name) const;
+
+	    /// \brief Disables the output from get, getDefault and getGroup.
+	    ///        By default, such output is enabled.
+	    void disableOutput();
+
+	    /// \brief Enables the output from get, getDefault and getGroup.
+	    ///        By default, such output is enabled.
+	    void enableOutput();
+
+	    /// \brief Returs true if and only if output from get, getDefault
+	    ///        and getGroup is enabled.
+	    ///
+	    /// \return true if and only if output from get, getDefault and
+	    ///         getGroup is enabled.
+	    bool isOutputEnabled() const;
+
+
+	    /// \brief Reads the contents of the  xml file specified by
+	    ///        xml_filename into this ParameterGroup.
+	    ///
+	    /// \param xml_filename is the name of a xml file.
+	    void readXML(const std::string& xml_filename);
+
+	    /// \brief Reads the contents of the param file specified by
+	    ///        param_filename into this ParameterGroup.
+	    ///
+	    /// NOTE: A param file contains lines on the form 'a/b/c=d'.
+	    //        The '/' separates ParameterGroups and Parameters
+	    ///       (e.g. c is a Parameter in the ParameterGroup b,
+	    ///       and b is a ParameterGroup in the ParameterGroup a)
+	    ///       while '=' separates the name from the value (e.g. the
+	    ///       value of the parameter c above will be d).
+	    /// NOTE: A param file does not store any type information about
+	    ///       its values.
+	    ///
+	    /// \param param_filename is the name of a param file.
+	    void readParam(const std::string& param_filename);
+
+	    /// \brief Writes this ParameterGroup into a param file
+	    ///        specified by param_filename.
+	    ///
+	    /// \param param_filename is the name of a param file.
+	    void writeParam(const std::string& param_filename) const;
+
+	    /// \brief Writes this ParameterGroup to a stream.
+	    ///
+	    /// \param stream is the stream to write to.
+	    void writeParamToStream(std::ostream& stream) const;
+
+
+	    /// vki param interface - deprecated
+	    template<typename T>
+	    void get(const char* name, T& value, const T& default_value) const {
+		value = this->getDefault<T>(name, default_value);
+	    }
+
+	    /// vki param interface - deprecated
+	    template<typename T>
+	    void get(const char* name, T& value) const {
+		value = this->get<T>(name);
+	    }
+
+            /// \brief Return true if any parameters are unused.
+            bool anyUnused() const;
+
+	    /// \brief Shows which parameters which are used or unused.
+	    void displayUsage(bool used_params = false) const;
+
+            /// \brief Returns the path of the parameter group.
+	    std::string path() const;
+
+            /// Insert a new item into the group.
+	    void insert(const std::string& name,
+                        const std::tr1::shared_ptr<ParameterMapItem>& data);
+
+            /// Insert a new parameter item into the group.
+	    void insertParameter(const std::string& name, const std::string& value);
+
+	private:
+	    typedef std::tr1::shared_ptr<ParameterMapItem> data_type;
+	    typedef std::pair<std::string, data_type> pair_type;
+	    typedef std::map<std::string, data_type> map_type;
+
+	    std::string path_;
+	    map_type map_;
+	    const ParameterGroup* parent_;
+	    bool output_is_enabled_;
+
+	    template<typename T, class Requirement>
+	    T translate(const pair_type& data, const Requirement& chk) const;
+            template <typename StringArray>
+	    void parseCommandLineArguments(int argc, StringArray argv);
+	    void recursiveSetIsOutputEnabled(bool output_is_enabled);
+	};
+    } // namespace parameter
+} // namespace Dune
+
+#include <dune/common/param/ParameterGroup_impl.hpp>
+
+#endif // OPENRS_PARAMETERGROUP_HEADER

dune/common/param/ParameterGroup_impl.hpp

@@ -0,0 +1,330 @@
+//===========================================================================
+//
+// File: ParameterGroup_impl.hpp
+//
+// Created: Tue Jun  2 19:06:46 2009
+//
+// Author(s): Bård Skaflestad     <bard.skaflestad@sintef.no>
+//            Atgeirr F Rasmussen <atgeirr@sintef.no>
+//
+// $Date$
+//
+// $Revision$
+//
+//===========================================================================
+
+/*
+Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+Copyright 2009, 2010 Statoil ASA.
+
+This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+OpenRS is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+OpenRS is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPENRS_PARAMETERGROUP_IMPL_HEADER
+#define OPENRS_PARAMETERGROUP_IMPL_HEADER
+
+#include <iostream>
+#include <string>
+
+//#include "TermColors.hpp" // from utils
+
+#include <dune/common/param/ParameterStrings.hpp>
+#include <dune/common/param/ParameterTools.hpp>
+#include <dune/common/param/Parameter.hpp>
+
+namespace Dune {
+    namespace parameter {
+
+	template<>
+	struct ParameterMapItemTrait<ParameterGroup> {
+	    static ParameterGroup
+            convert(const ParameterMapItem& item,
+                    std::string& conversion_error)
+            {
+		std::string tag = item.getTag();
+		if (tag != ID_xmltag__param_grp) {
+		    conversion_error = "The XML tag was '" + tag +
+                                       "' but should be '" +
+                                       ID_xmltag__param_grp + "'.\n";
+		    return ParameterGroup("", 0);
+		}
+		conversion_error = "";
+		const ParameterGroup& pg = dynamic_cast<const ParameterGroup&>(item);
+		return pg;
+	    }
+	    static std::string type() {return "ParameterGroup";}
+	};
+
+	namespace {
+	    template <typename T>
+	    inline std::string
+	    to_string(const T& val)
+	    {
+		std::ostringstream os;
+		os << val;
+		return os.str();
+	    }
+
+	    inline std::string
+	    to_string(const bool b) {
+		if (b) {
+		    return ID_true;
+		} else {
+		    return ID_false;
+		}
+	    }
+
+	    inline std::string
+	    to_string(const ParameterGroup&)
+	    {
+		return std::string("<parameter group>");
+	    }
+
+	    std::pair<std::string, std::string>
+            filename_split(const std::string& filename)
+            {
+		int fpos = filename.rfind('.');
+		std::string name = filename.substr(0, fpos);
+		std::string type = filename.substr(fpos+1);
+		return std::make_pair(name, type);
+	    }
+
+	}
+
+
+
+        template <typename StringArray>
+	ParameterGroup::ParameterGroup(int argc, StringArray argv)
+            : path_(ID_path_root), parent_(0), output_is_enabled_(true)
+	{
+	    if (argc < 2) {
+		std::cerr << "Usage: " << argv[0] << " "
+			  << "[paramfilename1.{xml,param}] "
+			  << "[paramfilename2.{xml,param}] "
+			  << "[overridden_arg1=value1] "
+			  << "[overridden_arg2=value2] "
+			  << "[...]" << std::endl;
+		exit(EXIT_FAILURE);
+	    }
+	    this->parseCommandLineArguments(argc, argv);
+	}
+
+        template <typename StringArray>
+	void ParameterGroup::parseCommandLineArguments(int argc, StringArray argv)
+        {
+	    std::vector<std::string> files;
+	    std::vector<std::pair<std::string, std::string> > assignments;
+	    for (int i = 1; i < argc; ++i) {
+		std::string arg(argv[i]);
+		int fpos = arg.find(ID_delimiter_assignment);
+		if (fpos == int(std::string::npos)) {
+		    std::string filename = arg.substr(0, fpos);
+		    files.push_back(filename);
+		    continue;
+		}
+		int pos = fpos + ID_delimiter_assignment.size();
+		int spos = arg.find(ID_delimiter_assignment, pos);
+		if (spos == int(std::string::npos)) {
+		    std::string name = arg.substr(0, fpos);
+		    std::string value = arg.substr(pos, spos);
+		    assignments.push_back(std::make_pair(name, value));
+		    continue;
+		}
+		std::cout << "WARNING: Too many assignements  (' "
+                          << ID_delimiter_assignment
+                          << "') detected in argument " << i << ".\n";
+	    }
+	    for (int i = 0; i < int(files.size()); ++i) {
+		std::pair<std::string, std::string> file_type = filename_split(files[i]);
+		if (file_type.second == "xml") {
+		    this->readXML(files[i]);
+		} else if (file_type.second == "param") {
+		    this->readParam(files[i]);
+		} else {
+		    std::cout << "WARNING: Ignoring file '"
+                              << files[i] << "' with unknown extension.\n"
+			      << "Valid filename extensions are 'xml' and 'param'.\n";
+		}
+	    }
+	    for (int i = 0; i < int(assignments.size()); ++i) {
+		this->insertParameter(assignments[i].first, assignments[i].second);
+	    }
+	}
+
+
+	template<typename T>
+	inline T ParameterGroup::get(const std::string& name) const
+        {
+	    return this->get<T>(name, ParameterRequirementNone());
+	}
+
+	template<typename T, class Requirement>
+	inline T ParameterGroup::get(const std::string& name,
+				     const Requirement& r) const
+        {
+	    setUsed();
+	    std::pair<std::string, std::string> name_path = split(name);
+	    map_type::const_iterator it = map_.find(name_path.first);
+	    if (it == map_.end()) {
+		if (parent_ != 0) {
+		    // If we have a parent, ask it instead.
+		    if (output_is_enabled_) {
+			//TermColors::Red();
+			std::cout << name;
+			//TermColors::Normal();
+			std::cout << " not found at "
+                                  << (path() + ID_delimiter_path)
+                                  << ", asking parent." << std::endl;
+		    }
+		    return parent_->get<T>(name, r);
+		} else {
+		    // We are at the top, name has not been found.
+		    std::cerr << "ERROR: The group '"
+			      << this->path()
+			      << "' does not contain an element named '"
+			      << name
+			      << "'.\n";
+		    throw NotFoundException();
+		}
+	    }
+	    if (name_path.second == "") {
+		T val = this->translate<T>(*it, r);
+		it->second->setUsed();
+		if (output_is_enabled_) {
+		    //TermColors::Green();
+		    std::cout << name;
+		    //TermColors::Normal();
+		    std::cout << " found at " << (path() + ID_delimiter_path)
+			      << ", value is " << to_string(val) << std::endl;
+		}
+		return val;
+	    } else {
+		ParameterGroup& pg = dynamic_cast<ParameterGroup&>(*(*it).second);
+		pg.setUsed();
+		return pg.get<T>(name_path.second, r);
+	    }
+	}
+
+	template<typename T>
+	inline T ParameterGroup::getDefault(const std::string& name,
+					    const T& default_value) const
+        {
+	    return this->getDefault<T>(name, default_value, ParameterRequirementNone());
+	}
+
+	template<typename T, class Requirement>
+	inline T ParameterGroup::getDefault(const std::string& name,
+					    const T& default_value,
+					    const Requirement& r) const
+        {
+	    setUsed();
+	    std::pair<std::string, std::string> name_path = split(name);
+	    map_type::const_iterator it = map_.find(name_path.first);
+	    if (it == map_.end()) {
+		if (parent_ != 0) {
+		    // If we have a parent, ask it instead.
+		    if (output_is_enabled_) {
+			//TermColors::Red();
+			std::cout << name;
+			//TermColors::Normal();
+			std::cout <<  " not found at " << (path() + ID_delimiter_path)
+				  << ", asking parent." << std::endl;
+		    }
+		    return parent_->getDefault<T>(name, default_value, r);
+		} else {
+		    // We check the requirement for the default value
+		    std::string requirement_result = r(default_value);
+		    if (requirement_result != "") {
+			std::cerr << "ERROR: The default value for the "
+				  << " element named '"
+				  << name
+				  << "' in the group '"
+				  << this->path()
+				  << "' failed to meet a requirenemt.\n";
+			std::cerr << "The requirement enforcer returned the following message:\n"
+				  << requirement_result
+				  << "\n";
+			throw RequirementFailedException<Requirement>();
+		    }
+		}
+		if (output_is_enabled_) {
+		    //TermColors::Blue();
+		    std::cout << name;
+		    //TermColors::Normal();
+		    std::cout << " not found. Using default value '"
+                              << to_string(default_value) << "'." << std::endl;
+		}
+		return default_value;
+	    }
+	    if (name_path.second == "") {
+		T val = this->translate<T>(*it, r);
+		it->second->setUsed();
+		if (output_is_enabled_) {
+		    //TermColors::Green();
+		    std::cout << name;
+		    //TermColors::Normal();
+		    std::cout << " found at " << (path() + ID_delimiter_path)
+			      << ", value is '" << to_string(val) << "'." << std::endl;
+		}
+		return val;
+	    } else {
+		ParameterGroup& pg = dynamic_cast<ParameterGroup&>(*(*it).second);
+		pg.setUsed();
+		return pg.getDefault<T>(name_path.second, default_value, r);
+	    }
+	}
+
+	template<typename T, class Requirement>
+	inline T ParameterGroup::translate(const pair_type& named_data,
+					   const Requirement& chk) const
+        {
+	    const std::string& name = named_data.first;
+	    const data_type data = named_data.second;
+	    std::string conversion_error;
+	    T value = ParameterMapItemTrait<T>::convert(*data, conversion_error);
+	    if (conversion_error != "") {
+		std::cerr << "ERROR: Failed to convert the element named '"
+			  << name
+			  << "' in the group '"
+			  << this->path()
+			  << "' to the type '"
+			  << ParameterMapItemTrait<T>::type()
+			  << "'.\n";
+		std::cerr << "The conversion routine returned the following message:\n"
+			  << conversion_error
+			  << "\n";
+		throw WrongTypeException();
+	    }
+	    std::string requirement_result = chk(value);
+	    if (requirement_result != "") {
+		std::cerr << "ERROR: The element named '"
+			  << name
+			  << "' in the group '"
+			  << this->path()
+			  << "' of type '"
+			  << ParameterMapItemTrait<T>::type()
+			  << "' failed to meet a requirenemt.\n";
+		std::cerr << "The requirement enforcer returned the following message:\n"
+			  << requirement_result
+			  << "\n";
+		throw RequirementFailedException<Requirement>();
+	    }
+	    return value;
+	}
+    } // namespace parameter
+} // namespace Dune
+
+#endif // OPENRS_PARAMETERGROUP_IMPL_HEADER

dune/common/param/ParameterMapItem.hpp

@@ -0,0 +1,73 @@
+//===========================================================================
+//
+// File: ParameterMapItem.hpp
+//
+// Created: Tue Jun  2 19:05:54 2009
+//
+// Author(s): Bård Skaflestad     <bard.skaflestad@sintef.no>
+//            Atgeirr F Rasmussen <atgeirr@sintef.no>
+//
+// $Date$
+//
+// $Revision$
+//
+//===========================================================================
+
+/*
+Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+Copyright 2009, 2010 Statoil ASA.
+
+This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+OpenRS is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+OpenRS is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPENRS_PARAMETERMAPITEM_HEADER
+#define OPENRS_PARAMETERMAPITEM_HEADER
+
+#include <string>
+
+namespace Dune {
+    namespace parameter {
+	/// The parameter handlig system is structured as a tree,
+	/// where each node inhertis from ParameterMapItem.
+	///
+	/// The abstract virtual function getTag() is  used to determine
+	/// which derived class the node actually is.
+	struct ParameterMapItem {
+	    /// Default constructor
+	    ParameterMapItem() : used_(false) {}
+
+	    /// Destructor
+	    virtual ~ParameterMapItem() {}
+
+	    /// \brief This function returns a string describing
+	    ///        the ParameterMapItem.
+	    virtual std::string getTag() const = 0;
+	    void setUsed() const { used_ = true; }
+	    bool used() const { return used_; }
+	private:
+	    mutable bool used_;
+	};
+
+	template<typename T>
+	struct ParameterMapItemTrait {
+	    static T convert(const ParameterMapItem&,
+                             std::string& conversion_error);
+	    static std::string type();
+	};
+    } // namespace parameter
+} // namespace Dune
+
+#endif // OPENRS_PARAMETERMAPITEM_HEADER

dune/common/param/ParameterRequirement.hpp

@@ -0,0 +1,243 @@
+//===========================================================================
+//
+// File: ParameterRequirement.hpp
+//
+// Created: Tue Jun  2 19:05:02 2009
+//
+// Author(s): Bård Skaflestad     <bard.skaflestad@sintef.no>
+//            Atgeirr F Rasmussen <atgeirr@sintef.no>
+//
+// $Date$
+//
+// $Revision$
+//
+//===========================================================================
+
+/*
+Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+Copyright 2009, 2010 Statoil ASA.
+
+This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+OpenRS is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+OpenRS is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPENRS_PARAMETERREQUIREMENT_HEADER
+#define OPENRS_PARAMETERREQUIREMENT_HEADER
+
+#include <algorithm>
+#include <cassert>
+#include <sstream>
+#include <string>
+#include <vector>
+
+namespace Dune {
+    namespace parameter {
+	/// @brief
+	/// @todo Doc me!
+	/// @tparam
+	/// @param
+	/// @return
+	struct ParameterRequirementNone {
+	    template<typename T>
+	    std::string operator()(const T&) const {
+		return "";
+	    }
+	};
+
+	/// @brief
+	/// @todo Doc me!
+	/// @param
+	/// @return
+	struct ParameterRequirementProbability {
+	    std::string operator()(double x) const {
+		if ( (x < 0.0) || (x > 1.0) ) {
+		    std::ostringstream stream;
+		    stream << "The value '" << x
+                           << "' is not in the interval [0, 1].";
+		    return stream.str();
+		} else {
+		    return "";
+		}
+	    }
+	};
+
+	/// @brief
+	/// @todo Doc me!
+	/// @tparam
+	/// @param
+	/// @return
+	struct ParameterRequirementPositive {
+	    template<typename T>
+	    std::string operator()(const T& x) const {
+		if (x > 0) {
+		    return "";
+		} else {
+		    std::ostringstream stream;
+		    stream << "The value '" << x << "' is not positive.";
+		    return stream.str();
+		}
+	    }
+	};
+
+	/// @brief
+	/// @todo Doc me!
+	/// @tparam
+	/// @param
+	/// @return
+	struct ParameterRequirementNegative {
+	    template<typename T>
+	    std::string operator()(const T& x) const {
+		if (x < 0) {
+		    return "";
+		} else {
+		    std::ostringstream stream;
+		    stream << "The value '" << x << "' is not negative.";
+		    return stream.str();
+		}
+	    }
+	};
+
+	/// @brief
+	/// @todo Doc me!
+	/// @tparam
+	/// @param
+	/// @return
+	struct ParameterRequirementNonPositive {
+	    template<typename T>
+	    std::string operator()(const T& x) const {
+		if (x > 0) {
+		    std::ostringstream stream;
+		    stream << "The value '" << x << "' is positive.";
+		    return stream.str();
+		} else {
+		    return "";
+		}
+	    }
+	};
+
+	/// @brief
+	/// @todo Doc me!
+	/// @tparam
+	/// @param
+	/// @return
+	struct ParameterRequirementNonNegative {
+	    template<typename T>
+	    std::string operator()(const T& x) const {
+		if (x < 0) {
+		    std::ostringstream stream;
+		    stream << "The value '" << x << "' is negative.";
+		    return stream.str();
+		} else {
+		    return "";
+		}
+	    }
+	};
+
+	/// @brief
+	/// @todo Doc me!
+	/// @tparam
+	/// @param
+	/// @return
+	struct ParameterRequirementNonZero {
+	    template<typename T>
+	    std::string operator()(const T& x) const {
+		if (x != 0) {
+		    return "";
+		} else {
+		    return "The value was zero.";
+		}
+	    }
+	};
+
+	/// @brief
+	/// @todo Doc me!
+	/// @param
+	/// @return
+	struct ParameterRequirementNonEmpty {
+	    std::string operator()(const std::string& x) const {
+		if (x != "") {
+		    return "The string was empty.";
+		} else {
+		    return "";
+		}
+	    }
+	};
+
+	/// @brief
+	/// @todo Doc me!
+	/// @tparam
+	/// @param
+	/// @return
+	template<class Requirement1, class Requirement2>
+	struct ParameterRequirementAnd {
+	    ParameterRequirementAnd(const Requirement1& r1, const Requirement2& r2) :
+	    r1_(r1), r2_(r2) { }
+
+	    template<typename T>
+	    std::string operator()(const T& t) const {
+		std::string e1 = r1_(t);
+		std::string e2 = r2_(t);
+		if (e1 == "") {
+		    return e2;
+		} else if (e2 == "") {
+		    return e1;
+		} else {
+		    return e1 + " AND " + e2;
+		}
+	    }
+	private:
+	    const Requirement1 r1_;
+	    const Requirement2 r2_;
+	};
+
+	/// @brief
+	/// @todo Doc me!
+	/// @param
+	struct ParameterRequirementMemberOf {
+	    ParameterRequirementMemberOf(const std::vector<std::string>& elements)
+	    : elements_(elements) {
+		assert(elements_.size() > 0);
+	    }
+
+	/// @brief
+	/// @todo Doc me!
+	/// @param
+	/// @return
+	    std::string operator()(const std::string& x) const {
+		if (std::find(elements_.begin(), elements_.end(), x) == elements_.end()) {
+		    if (elements_.size() == 1) {
+			return "The string '" + x + "' is not '" + elements_[0] + "'.";
+		    }
+		    std::ostringstream stream;
+		    stream << "The string '" << x << "' is not among '";
+		    for (int i = 0; i < int(elements_.size()) - 2; ++i) {
+			stream << elements_[i] << "', '";
+		    }
+		    stream << elements_[elements_.size() - 2]
+			   << "' and '"
+			   << elements_[elements_.size() - 1]
+			   << "'.";
+		    return stream.str();
+		} else {
+		    return "";
+		}
+	    }
+	private:
+	    const std::vector<std::string> elements_;
+	};
+    } // namespace parameter
+} // namespace Dune
+
+#endif // OPENRS_PARAMETERREQUIREMENT_HEADER

dune/common/param/ParameterStrings.hpp

@@ -0,0 +1,70 @@
+//===========================================================================
+//
+// File: ParameterStrings.hpp
+//
+// Created: Tue Jun  2 19:04:15 2009
+//
+// Author(s): Bård Skaflestad     <bard.skaflestad@sintef.no>
+//            Atgeirr F Rasmussen <atgeirr@sintef.no>
+//
+// $Date$
+//
+// $Revision$
+//
+//===========================================================================
+
+/*
+Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+Copyright 2009, 2010 Statoil ASA.
+
+This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+OpenRS is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+OpenRS is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPENRS_PARAMETERSTRINGS_HEADER
+#define OPENRS_PARAMETERSTRINGS_HEADER
+
+#include <string>
+
+namespace Dune {
+    namespace parameter {
+	const std::string ID_true                   = "true";
+	const std::string ID_false                  = "false";
+
+	const std::string ID_xmltag__param_grp      = "ParameterGroup";
+	const std::string ID_xmltag__param          = "Parameter";
+	const std::string ID_xmltag__file_param_grp = "ParameterGroupFromFile";
+	const std::string ID_xmltag__file_params    = "MergeWithFile";
+
+	const std::string ID_xmlatt__value          = "value";
+	const std::string ID_xmlatt__name           = "name";
+	const std::string ID_xmlatt__type           = "type";
+
+	const std::string ID_param_type__bool       = "bool";
+	const std::string ID_param_type__int        = "int";
+	const std::string ID_param_type__float      = "double";
+	const std::string ID_param_type__string     = "string";
+	const std::string ID_param_type__file       = "file";
+	const std::string ID_param_type__cmdline    = "cmdline";
+
+	//
+
+	const std::string ID_path_root            = "";
+	const std::string ID_delimiter_path       = "/";
+	const std::string ID_delimiter_assignment = "=";
+    } // namespace parameter
+} // namespace Dune
+
+#endif // OPENRS_PARAMETERSTRINGS_HEADER

dune/common/param/ParameterTools.cpp

@@ -0,0 +1,55 @@
+//===========================================================================
+//
+// File: ParameterTools.cpp
+//
+// Created: Tue Jun  2 19:03:09 2009
+//
+// Author(s): Bård Skaflestad     <bard.skaflestad@sintef.no>
+//            Atgeirr F Rasmussen <atgeirr@sintef.no>
+//
+// $Date$
+//
+// $Revision$
+//
+//===========================================================================
+
+/*
+Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+Copyright 2009, 2010 Statoil ASA.
+
+This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+OpenRS is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+OpenRS is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#if HAVE_CONFIG_H
+#include "config.h"
+#endif
+#include <dune/common/param/ParameterTools.hpp>
+#include <dune/common/param/ParameterStrings.hpp>
+
+namespace Dune {
+    namespace parameter {
+	std::pair<std::string, std::string> split(const std::string& name)
+        {
+	    int pos = name.find(ID_delimiter_path);
+	    if (pos == int(std::string::npos)) {
+		return std::make_pair(name, "");
+	    } else {
+		return std::make_pair(name.substr(0, pos),
+                                      name.substr(pos + ID_delimiter_path.size()));
+	    }
+	}
+    } // namespace parameter
+} // namespace Dune

dune/common/param/ParameterTools.hpp

@@ -0,0 +1,48 @@
+//===========================================================================
+//
+// File: ParameterTools.hpp
+//
+// Created: Tue Jun  2 19:02:19 2009
+//
+// Author(s): Bård Skaflestad     <bard.skaflestad@sintef.no>
+//            Atgeirr F Rasmussen <atgeirr@sintef.no>
+//
+// $Date$
+//
+// $Revision$
+//
+//===========================================================================
+
+/*
+Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+Copyright 2009, 2010 Statoil ASA.
+
+This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+OpenRS is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+OpenRS is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPENRS_PARAMETERTOOLS_HEADER
+#define OPENRS_PARAMETERTOOLS_HEADER
+
+#include <string>
+#include <utility>
+
+namespace Dune {
+    namespace parameter {
+	std::pair<std::string, std::string> split(const std::string& name);
+    } // namespace parameter
+} // namespace Dune
+
+#endif // OPENRS_PARAMETERTOOLS_HEADER

dune/common/param/ParameterXML.cpp

@@ -0,0 +1,140 @@
+//===========================================================================
+//
+// File: ParameterXML.cpp
+//
+// Created: Tue Jun  2 18:57:25 2009
+//
+// Author(s): Bård Skaflestad     <bard.skaflestad@sintef.no>
+//            Atgeirr F Rasmussen <atgeirr@sintef.no>
+//
+// $Date$
+//
+// $Revision$
+//
+//===========================================================================
+
+/*
+Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+Copyright 2009, 2010 Statoil ASA.
+
+This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+OpenRS is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+OpenRS is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#if HAVE_CONFIG_H
+#include "config.h"
+#endif
+#include <dune/common/param/ParameterXML.hpp>
+#include <dune/common/param/Parameter.hpp>
+#include <dune/common/param/ParameterStrings.hpp>
+
+#include <exception>
+#include <iostream>
+#include <string>
+#include <boost/filesystem.hpp>
+
+#include <libxml/parser.h> // libxml2
+#include <libxml/tree.h>   // libxml2
+
+
+namespace Dune {
+    namespace parameter {
+
+	namespace libxml2 {
+	    void read_xml(ParameterGroup& pg, const std::string filename);
+	    void fill_tree(ParameterGroup& pg, xmlNodePtr root,
+                           const std::string& xml_file_path);
+	}
+
+	void fill_xml(ParameterGroup& pg, const std::string filename) {
+	    libxml2::read_xml(pg, filename);
+	}
+
+	namespace libxml2 {
+	    std::string getProperty(const std::string& property,
+                                    const xmlNodePtr& node_ptr)
+            {
+		const char* prop_value_ptr = (const char*) xmlGetProp(node_ptr, (const xmlChar*)property.c_str());
+		std::string property_value(prop_value_ptr ? prop_value_ptr : "");
+		xmlFree((xmlChar*) prop_value_ptr);
+		return property_value;
+	    }
+
+	    void read_xml(ParameterGroup& pg, const std::string filename)
+            {
+		// Macro to check that the libxml version in use is compatible
+		// with the version the software has been compiled against
+		LIBXML_TEST_VERSION; // defined by libxml2
+
+		xmlDocPtr doc = xmlParseFile(filename.c_str());
+		if (doc == NULL) {
+		    xmlFreeDoc(doc);
+		    xmlCleanupParser();
+		    std::cerr << "ERROR: Failed to open XML file '" << filename << "'\n";
+		    throw std::exception();
+		}
+		xmlNodePtr root = xmlDocGetRootElement(doc);
+		std::string xml_file_path = boost::filesystem::path(filename).branch_path().string();
+		fill_tree(pg, root, xml_file_path);
+		xmlFreeDoc(doc);
+		xmlCleanupParser();
+	    }
+
+	    void fill_tree(ParameterGroup& pg, xmlNodePtr root,
+                           const std::string& xml_file_path)
+            {
+		//std::cout << "GROUP '" << value << "' BEGIN\n";
+		for (xmlNodePtr cur = root->children; cur; cur = cur->next) {
+		    if (cur->type == XML_ELEMENT_NODE) {
+			const char* tag_ptr = (const char*) cur->name;
+			std::string tag_name(tag_ptr ? tag_ptr : "");
+			if (tag_name == ID_xmltag__file_params) {
+			    std::string relative_filename = getProperty(ID_xmlatt__value, cur);
+			    std::string filename = (boost::filesystem::path(xml_file_path) / relative_filename).string();
+			    fill_xml(pg, filename);
+			    continue;
+			}
+			std::string name = getProperty(ID_xmlatt__name, cur);
+			std::tr1::shared_ptr<ParameterMapItem> data;
+			if (tag_name == ID_xmltag__param) {
+			    std::string value = getProperty(ID_xmlatt__value, cur);
+			    std::string type = getProperty(ID_xmlatt__type, cur);
+			    if (type == ID_param_type__file) {
+				value = (boost::filesystem::path(xml_file_path) / value).string();
+				type  = ID_param_type__string;
+			    }
+			    data.reset(new Parameter(value, type));
+			} else if (tag_name == ID_xmltag__param_grp) {
+			    std::string child_path  = pg.path() + ID_delimiter_path + name;
+			    data.reset(new ParameterGroup(child_path, &pg));
+			    fill_tree(dynamic_cast<ParameterGroup&>(*data), cur, xml_file_path);
+			} else if (tag_name == ID_xmltag__file_param_grp) {
+			    std::string child_path  = pg.path() + ID_delimiter_path + name;
+			    data.reset(new ParameterGroup(child_path, &pg));
+			    std::string relative_filename = getProperty(ID_xmlatt__value, cur);
+			    std::string filename = (boost::filesystem::path(xml_file_path) / relative_filename).string();
+			    fill_xml(dynamic_cast<ParameterGroup&>(*data), filename);
+			} else {
+			    std::cerr << "ERROR: '" << tag_name << "' is an unknown xml tag.\n";
+			    throw std::exception();
+			}
+			pg.insert(name, data);
+		    }
+		}
+		//std::cout << "GROUP '" << id << "' END\n";
+	    }
+	}
+    } // namespace parameter
+} // namespace Dune

dune/common/param/ParameterXML.hpp

@@ -0,0 +1,54 @@
+//===========================================================================
+//
+// File: ParameterXML.hpp
+//
+// Created: Tue Jun  2 18:55:59 2009
+//
+// Author(s): Bård Skaflestad     <bard.skaflestad@sintef.no>
+//            Atgeirr F Rasmussen <atgeirr@sintef.no>
+//
+// $Date$
+//
+// $Revision$
+//
+//===========================================================================
+
+/*
+Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+Copyright 2009, 2010 Statoil ASA.
+
+This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+OpenRS is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+OpenRS is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPENRS_PARAMETERXML_HEADER
+#define OPENRS_PARAMETERXML_HEADER
+
+#include <dune/common/param/ParameterGroup.hpp>
+
+namespace Dune {
+    namespace parameter {
+	/// \brief This function fills a ParameterGroup with a  XML file.
+	///
+	/// NOTE: The function throws if there is an error during reading
+	/// of the XML file.
+	///
+	/// \retval pg is the ParameterGroup to be filled.
+	/// \param filename is the name of an XML file.
+	void fill_xml(ParameterGroup& pg, const std::string filename);
+    } // namespace parameter
+} // namespace Dune
+
+#endif // OPENRS_PARAMETERXML_HEADER

dune/common/param/test/Makefile.am

@@ -0,0 +1,14 @@
+# $Date: 2009-12-08 15:21:36 +0100 (Tue, 08 Dec 2009) $
+# $Revision: 699 $
+
+check_PROGRAMS = param_test
+
+param_test_SOURCES = param_test.cpp
+param_test_CXXFLAGS = $(DUNEMPICPPFLAGS) $(BOOST_CPPFLAGS)
+param_test_LDADD = $(DUNE_LDFLAGS) $(DUNEMPILDFLAGS) $(DUNEMPILIBS) $(BOOST_LDFLAGS) \
+                         $(DUNE_LIBS) $(BOOST_UNIT_TEST_FRAMEWORK_LIB) ../libparam.la
+
+TESTS = $(check_PROGRAMS)
+
+include $(top_srcdir)/am/global-rules
+

dune/common/param/test/param_test.cpp

@@ -0,0 +1,66 @@
+//===========================================================================
+//
+// File: param_test.cpp
+//
+// Created: Sun Dec 13 20:08:36 2009
+//
+// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
+//            Bård Skaflestad     <bard.skaflestad@sintef.no>
+//
+// $Date$
+//
+// $Revision$
+//
+//===========================================================================
+
+/*
+  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+  Copyright 2009, 2010 Statoil ASA.
+
+  This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+  OpenRS is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OpenRS is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+#define BOOST_TEST_DYN_LINK
+#define NVERBOSE // to suppress our messages when throwing
+
+#define BOOST_TEST_MODULE ParameterTest
+#include <boost/test/unit_test.hpp>
+
+#include "../ParameterGroup.hpp"
+#include <cstddef>
+
+using namespace Dune;
+
+BOOST_AUTO_TEST_CASE(commandline_syntax_init)
+{
+    typedef const char* cp;
+    cp argv[] = { "program_command",
+                  "topitem=somestring",
+                  "/slashtopitem=anotherstring",
+                  "/group/item=1",
+                  "/group/anotheritem=2",
+                  "/group/subgroup/item=3",
+                  "/group/subgroup/anotheritem=4",
+                  "/group/item=overridingstring" };
+    const std::size_t argc = sizeof(argv)/sizeof(argv[0]);
+    parameter::ParameterGroup p(argc, argv);
+    BOOST_CHECK(p.get<std::string>("topitem") == "somestring");
+
+    // Tests that only run in debug mode.
+#ifndef NDEBUG
+#endif
+}

dune/common/test/Makefile.am

@@ -0,0 +1,29 @@
+# $Date$
+# $Revision$
+
+check_PROGRAMS = sparsetable_test sparsevector_test monotcubicinterpolator_test
+noinst_PROGRAMS = unit_test
+
+sparsetable_test_SOURCES = sparsetable_test.cpp
+sparsetable_test_CXXFLAGS = $(DUNEMPICPPFLAGS) $(BOOST_CPPFLAGS)
+sparsetable_test_LDADD = $(DUNE_LDFLAGS) $(DUNEMPILDFLAGS) $(BOOST_LDFLAGS) \
+                         $(DUNE_LIBS) $(DUNEMPILIBS) $(BOOST_UNIT_TEST_FRAMEWORK_LIB)
+
+sparsevector_test_SOURCES = sparsevector_test.cpp
+sparsevector_test_CXXFLAGS = $(DUNEMPICPPFLAGS) $(BOOST_CPPFLAGS)
+sparsevector_test_LDADD = $(DUNE_LDFLAGS) $(DUNEMPILDFLAGS) $(BOOST_LDFLAGS) \
+                          $(DUNE_LIBS) $(DUNEMPILIBS) $(BOOST_UNIT_TEST_FRAMEWORK_LIB)
+
+monotcubicinterpolator_test_SOURCES = monotcubicinterpolator_test.cpp
+monotcubicinterpolator_test_CXXFLAGS = $(DUNEMPICPPFLAGS) $(BOOST_CPPFLAGS)
+monotcubicinterpolator_test_LDADD = $(DUNE_LDFLAGS) $(DUNEMPILDFLAGS) \
+                                    $(DUNE_LIBS) $(DUNEMPILIBS) ../libcommon.la
+
+unit_test_SOURCES = unit_test.cpp
+unit_test_CXXFLAGS = $(DUNEMPICPPFLAGS) $(BOOST_CPPFLAGS)
+unit_test_LDADD = $(DUNE_LDFLAGS) $(DUNEMPILDFLAGS) $(DUNE_LIBS) $(DUNEMPILIBS)
+
+TESTS = $(check_PROGRAMS)
+
+include $(top_srcdir)/am/global-rules
+

dune/common/test/monotcubicinterpolator_test.cpp

@@ -0,0 +1,53 @@
+//===========================================================================
+//
+// File: monotcubicinterpolator_test.cpp
+//
+// Created: Tue Dec  8 12:25:30 2009
+//
+// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
+//            Bård Skaflestad     <bard.skaflestad@sintef.no>
+//
+// $Date$
+//
+// $Revision$
+//
+//===========================================================================
+
+/*
+  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+  Copyright 2009, 2010 Statoil ASA.
+
+  This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+  OpenRS is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OpenRS is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <dune/common/MonotCubicInterpolator.hpp>
+
+int main()
+{
+    const int num_v = 3;
+    double xv[num_v] = {0.0, 1.0, 2.0};
+    double fv[num_v] = {10.0, 21.0, 2.0};
+    std::vector<double> x(xv, xv + num_v);
+    std::vector<double> f(fv, fv + num_v);
+    MonotCubicInterpolator interp(x, f);
+    interp.evaluate(-1.0);
+    interp.evaluate(0.0);
+    interp.evaluate(0.0001);
+    interp.evaluate(0.5);
+    interp.evaluate(1.0);
+    interp.evaluate(2.0);
+    interp.evaluate(4.0);
+}

dune/common/test/sparsetable_test.cpp

@@ -0,0 +1,128 @@
+//===========================================================================
+//
+// File: test_sparsetable.cpp
+//
+// Created: Thu May 28 10:01:46 2009
+//
+// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
+//            Bård Skaflestad     <bard.skaflestad@sintef.no>
+//
+// $Date$
+//
+// $Revision$
+//
+//===========================================================================
+
+/*
+  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+  Copyright 2009, 2010 Statoil ASA.
+
+  This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+  OpenRS is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OpenRS is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#define BOOST_TEST_DYN_LINK
+#define NVERBOSE // to suppress our messages when throwing
+
+#define BOOST_TEST_MODULE SparseTableTest
+#include <boost/test/unit_test.hpp>
+
+#include "../SparseTable.hpp"
+
+using namespace Dune;
+
+BOOST_AUTO_TEST_CASE(construction_and_queries)
+{
+    const SparseTable<int> st1;
+    BOOST_CHECK(st1.empty());
+    BOOST_CHECK_EQUAL(st1.size(), 0);
+    BOOST_CHECK_EQUAL(st1.dataSize(), 0);
+
+    // This should be getting us a table like this:
+    // ----------------
+    // 0
+    // <empty row>
+    // 1 2
+    // 3 4 5 6
+    // 7 8 9
+    // ----------------
+    const int num_elem = 10;
+    const int elem[num_elem] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 };
+    const int num_rows = 5;
+    const int rowsizes[num_rows] = { 1, 0, 2, 4, 3 };
+    const SparseTable<int> st2(elem, elem + num_elem, rowsizes, rowsizes + num_rows);
+    BOOST_CHECK(!st2.empty());
+    BOOST_CHECK_EQUAL(st2.size(), num_rows);
+    BOOST_CHECK_EQUAL(st2.dataSize(), num_elem);
+    BOOST_CHECK_EQUAL(st2[0][0], 0);
+    BOOST_CHECK_EQUAL(st2.rowSize(0), 1);
+    BOOST_CHECK(st2[1].empty());
+    BOOST_CHECK_EQUAL(st2.rowSize(1), 0);
+    BOOST_CHECK_EQUAL(st2[3][1], 4);
+    BOOST_CHECK_EQUAL(st2[4][2], 9);
+    BOOST_CHECK(st2[4].size() == rowsizes[4]);
+    const SparseTable<int> st2_again(elem, elem + num_elem, rowsizes, rowsizes + num_rows);
+    BOOST_CHECK(st2 == st2_again);
+    SparseTable<int> st2_byassign;
+    st2_byassign.assign(elem, elem + num_elem, rowsizes, rowsizes + num_rows);
+    BOOST_CHECK(st2 == st2_byassign);
+    const int last_row_size = rowsizes[num_rows - 1];
+    SparseTable<int> st2_append(elem, elem + num_elem - last_row_size, rowsizes, rowsizes + num_rows - 1);
+    BOOST_CHECK_EQUAL(st2_append.dataSize(), num_elem - last_row_size);
+    st2_append.appendRow(elem + num_elem - last_row_size, elem + num_elem);
+    BOOST_CHECK(st2 == st2_append);
+    SparseTable<int> st2_append2;
+    st2_append2.appendRow(elem, elem + 1);
+    st2_append2.appendRow(elem + 1, elem + 1);
+    st2_append2.appendRow(elem + 1, elem + 3);
+    st2_append2.appendRow(elem + 3, elem + 7);
+    st2_append2.appendRow(elem + 7, elem + 10);
+    BOOST_CHECK(st2 == st2_append2);
+    st2_append2.clear();
+    SparseTable<int> st_empty;
+    BOOST_CHECK(st2_append2 == st_empty);
+
+    SparseTable<int> st2_allocate;
+    st2_allocate.allocate(rowsizes, rowsizes + num_rows);
+    BOOST_CHECK_EQUAL(st2_allocate.size(), num_rows);
+    BOOST_CHECK_EQUAL(st2_allocate.dataSize(), num_elem);
+    int s = 0;
+    for (int i = 0; i < num_rows; ++i) {
+        SparseTable<int>::mutable_row_type row = st2_allocate[i];
+        for (int j = 0; j < rowsizes[i]; ++j, ++s)
+            row[j] = elem[s];
+    }
+    BOOST_CHECK(st2 == st2_allocate);
+    
+    // One element too few.
+    BOOST_CHECK_THROW(const SparseTable<int> st3(elem, elem + num_elem - 1, rowsizes, rowsizes + num_rows), std::exception);
+
+    // A few elements too many.
+    BOOST_CHECK_THROW(const SparseTable<int> st4(elem, elem + num_elem, rowsizes, rowsizes + num_rows - 1), std::exception);
+
+    // Need at least one row.
+    BOOST_CHECK_THROW(const SparseTable<int> st5(elem, elem + num_elem, rowsizes, rowsizes), std::exception);
+
+    // Tests that only run in debug mode.
+#ifndef NDEBUG
+    // Do not ask for wrong row numbers.
+    BOOST_CHECK_THROW(st1.rowSize(0), std::exception);
+    BOOST_CHECK_THROW(st2.rowSize(-1), std::exception);
+    BOOST_CHECK_THROW(st2.rowSize(st2.size()), std::exception);
+    // No negative row sizes.
+    const int err_rs[num_rows] = { 1, 0, -1, 7, 3 };
+    BOOST_CHECK_THROW(const SparseTable<int> st6(elem, elem + num_elem, err_rs, err_rs + num_rows), std::exception);
+#endif
+}

dune/common/test/sparsevector_test.cpp

@@ -0,0 +1,106 @@
+//===========================================================================
+//
+// File: sparsevector_test.cpp
+//
+// Created: Mon Jun 29 21:00:53 2009
+//
+// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
+//            Bård Skaflestad     <bard.skaflestad@sintef.no>
+//
+// $Date$
+//
+// $Revision$
+//
+//===========================================================================
+
+/*
+  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+  Copyright 2009, 2010 Statoil ASA.
+
+  This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+  OpenRS is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OpenRS is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+#define BOOST_TEST_DYN_LINK
+#define NVERBOSE // to suppress our messages when throwing
+
+#define BOOST_TEST_MODULE SparseVectorTest
+#include <boost/test/unit_test.hpp>
+
+#include "../SparseVector.hpp"
+
+using namespace Dune;
+
+BOOST_AUTO_TEST_CASE(construction_and_queries)
+{
+    const SparseVector<int> sv1;
+    BOOST_CHECK(sv1.empty());
+    BOOST_CHECK_EQUAL(sv1.size(), 0);
+    BOOST_CHECK_EQUAL(sv1.nonzeroSize(), 0);
+
+    const int size = 100;
+    const int num_elem = 9;
+    const int elem[num_elem] = { 9, 8, 7, 6, 5, 4, 3, 2, 1 };
+    const int indices[num_elem] = { 1, 2, 3, 5, 8, 13, 21, 34, 55 };
+    const SparseVector<int> sv2(size, elem, elem + num_elem, indices, indices + num_elem);
+    BOOST_CHECK(!sv2.empty());
+    BOOST_CHECK_EQUAL(sv2.size(), size);
+    BOOST_CHECK_EQUAL(sv2.element(0), 0);
+    BOOST_CHECK_EQUAL(sv2.element(1), 9);
+    BOOST_CHECK_EQUAL(sv2.element(2), 8);
+    BOOST_CHECK_EQUAL(sv2.element(3), 7);
+    BOOST_CHECK_EQUAL(sv2.element(4), 0);
+    BOOST_CHECK_EQUAL(sv2.element(5), 6);
+    BOOST_CHECK_EQUAL(sv2.element(55), 1);
+    BOOST_CHECK_EQUAL(sv2.element(99), 0);
+    BOOST_CHECK_EQUAL(sv2.nonzeroSize(), num_elem);
+    for (int i = 0; i < num_elem; ++i) {
+	BOOST_CHECK_EQUAL(sv2.nonzeroElement(i), elem[i]);
+    }
+    const SparseVector<int> sv2_again(size, elem, elem + num_elem, indices, indices + num_elem);
+    BOOST_CHECK(sv2 == sv2_again);
+    SparseVector<int> sv2_append(size, elem, elem + num_elem - 1, indices, indices + num_elem - 1);
+    BOOST_CHECK_EQUAL(sv2_append.nonzeroSize(), num_elem - 1);
+    sv2_append.addElement(elem[num_elem - 1], indices[num_elem - 1]);
+    BOOST_CHECK(sv2 == sv2_append);
+    SparseVector<int> sv2_append2(size);
+    for (int i = 0; i < num_elem; ++i) {
+	sv2_append2.addElement(elem[i], indices[i]);
+    }
+    BOOST_CHECK(sv2 == sv2_append2);
+    sv2_append2.clear();
+    SparseVector<int> sv_empty;
+    BOOST_CHECK(sv2_append2 == sv_empty);
+
+    // Tests that only run in debug mode.
+#ifndef NDEBUG
+    // One element too few.
+    BOOST_CHECK_THROW(const SparseVector<int> sv3(size, elem, elem + num_elem - 1, indices, indices + num_elem), std::exception);
+    // One element too many.
+    BOOST_CHECK_THROW(const SparseVector<int> sv4(size, elem, elem + num_elem, indices, indices + num_elem - 1), std::exception);
+    // Indices out of range.
+    BOOST_CHECK_THROW(const SparseVector<int> sv5(4, elem, elem + num_elem, indices, indices + num_elem), std::exception);
+    // Indices not strictly increasing. Cheating by using the elements as indices.
+    BOOST_CHECK_THROW(const SparseVector<int> sv5(size, elem, elem + num_elem, elem, elem + num_elem), std::exception);
+
+    // Do not ask for out of range indices.
+    BOOST_CHECK_THROW(sv1.element(0), std::exception);
+    BOOST_CHECK_THROW(sv2.element(-1), std::exception);
+    BOOST_CHECK_THROW(sv2.element(sv2.size()), std::exception);
+    BOOST_CHECK_THROW(sv2.nonzeroElement(sv2.nonzeroSize()), std::exception);
+#endif
+}
+

dune/common/test/unit_test.cpp

@@ -0,0 +1,78 @@
+//===========================================================================
+//
+// File: unit_test.cpp
+//
+// Created: Fri Jul 17 11:02:33 2009
+//
+// Author(s): Bård Skaflestad     <bard.skaflestad@sintef.no>
+//            Atgeirr F Rasmussen <atgeirr@sintef.no>
+//
+// $Date$
+//
+// $Revision$
+//
+//===========================================================================
+
+/*
+  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
+  Copyright 2009, 2010 Statoil ASA.
+
+  This file is part of The Open Reservoir Simulator Project (OpenRS).
+
+  OpenRS is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OpenRS is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <algorithm>
+#include <iostream>
+#include <iterator>
+#include <vector>
+
+#include <boost/bind.hpp>
+
+#include "../Units.hpp"
+
+using namespace Dune::prefix;
+using namespace Dune::unit;
+
+int main()
+{
+    std::cout << "m  = " << meter << '\n';
+    std::cout << "kg = " << kilogram << '\n';
+    std::cout << "s  = " << second << '\n';
+
+    std::cout << "mD = " << milli*darcy << '\n';
+    std::cout << "MD = " << mega*darcy << '\n';
+
+    std::cout << "MD = " << convert::to(mega*darcy, milli*darcy) << " mD\n";
+
+    std::cout << "1 bar = " << convert::to(convert::from(1.0, barsa), psia) << " psi\n";
+
+    std::cout << "1 atm = " << convert::to(1*atm, barsa) << " bar\n";
+
+    std::vector<double> flux(10, 10000*cubic(meter)/year);
+
+    std::cout << "flux = [";
+    std::copy(flux.begin(), flux.end(),
+              std::ostream_iterator<double>(std::cout, " "));
+    std::cout << "\b] (m^3/s)\n";
+
+    std::transform(flux.begin(), flux.end(), flux.begin(),
+                   boost::bind(convert::to, _1, cubic(meter)/year));
+    std::cout << "     = [";
+    std::copy(flux.begin(), flux.end(),
+              std::ostream_iterator<double>(std::cout, " "));
+    std::cout << "\b] (m^3/year)\n";
+
+    return 0;
+}