diff --git a/opm/common/exceptions.hh b/opm/common/exceptions.hh
deleted file mode 100644
index 3004455fb..000000000
--- a/opm/common/exceptions.hh
+++ /dev/null
@@ -1,77 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * Copyright (C) 2008-2012 by Andreas Lauser *
- * *
- * 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, see . *
- *****************************************************************************/
-/*!
- * \file
- * \brief Some exceptions thrown by opm-material.
- */
-#ifndef OPM_EXCEPTIONS_HH
-#define OPM_EXCEPTIONS_HH
-
-#include
-
-#include
-
-namespace Opm {
-/*!
- * \ingroup Exception
- * \brief Exception thrown if a fixable numerical problem occurs.
- *
- * (e.g. time step too big, etc.)
- */
-class NumericalProblem : public Dune::Exception
-{
-public:
- // copy constructor
- NumericalProblem(const NumericalProblem &v)
- : Dune::Exception(v)
- {}
-
- // default constructor
- NumericalProblem()
- {}
-
- // constructor with error message
- NumericalProblem(const std::string &s)
- { this->message(s); }
-};
-
-/*!
- * \ingroup Exception
- * \brief Exception thrown if a run-time parameter is not specified correctly.
- */
-class ParameterException : public Dune::Exception
-{
-public:
- // copy constructor
- ParameterException(const ParameterException &v)
- : Dune::Exception(v)
- {}
-
- // default constructor
- ParameterException()
- {}
-
- // constructor with error message
- ParameterException(const std::string &s)
- { this->message(s); }
-};
-
-} // namespace Opm
-
-#endif
diff --git a/opm/common/fixedlengthspline_.hh b/opm/common/fixedlengthspline_.hh
deleted file mode 100644
index 8e090580c..000000000
--- a/opm/common/fixedlengthspline_.hh
+++ /dev/null
@@ -1,412 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * Copyright (C) 2010-2012 by Andreas Lauser *
- * *
- * 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, see . *
- *****************************************************************************/
-/*!
- * \file
- * \copydoc Opm::FixedLengthSpline_
- */
-#ifndef OPM_FIXED_LENGTH_SPLINE_HH
-#define OPM_FIXED_LENGTH_SPLINE_HH
-
-#include
-#include
-#include
-
-#include "splinecommon_.hh"
-
-namespace Opm
-{
-/*!
- * \brief The common code for all 3rd order polynomial splines with
- * more than two sampling points.
- */
-template
-class FixedLengthSpline_
- : public SplineCommon_ >
-{
- friend class SplineCommon_ >;
-
- typedef ScalarT Scalar;
- typedef Dune::FieldVector Vector;
- typedef Dune::BlockVector > BlockVector;
- typedef Dune::BTDMatrix > BTDMatrix;
-
-protected:
- FixedLengthSpline_()
- : m_(nSamples)
- {}
-
-public:
- /*!
- * \brief Returns the number of sampling points.
- */
- size_t numSamples() const
- { return nSamples; }
-
- ///////////////////////////////////////
- ///////////////////////////////////////
- ///////////////////////////////////////
- // Full splines //
- ///////////////////////////////////////
- ///////////////////////////////////////
- ///////////////////////////////////////
-
- /*!
- * \brief Set the sampling points and the boundary slopes of a
- * full spline using C-style arrays.
- *
- * This method uses separate array-like objects for the values of
- * the X and Y coordinates. In this context 'array-like' means
- * that an access to the members is provided via the []
- * operator. (e.g. C arrays, std::vector, std::array, etc.) Each
- * array must be of size 'numberSamples' at least. Also, the number of
- * sampling points must be larger than 1.
- */
- template
- void setXYArrays(size_t numberSamples,
- const ScalarArrayX &x,
- const ScalarArrayY &y,
- Scalar m0, Scalar m1)
- {
- assert(numberSamples == numSamples());
-
- for (size_t i = 0; i < numberSamples; ++i) {
- xPos_[i] = x[i];
- yPos_[i] = y[i];
- }
- makeFullSpline_(m0, m1);
- }
-
- /*!
- * \brief Set the sampling points and the boundary slopes of a
- * full spline using STL-compatible containers.
- *
- * This method uses separate STL-compatible containers for the
- * values of the sampling points' X and Y
- * coordinates. "STL-compatible" means that the container provides
- * access to iterators using the begin(), end() methods and also
- * provides a size() method. Also, the number of entries in the X
- * and the Y containers must be equal and larger than 1.
- */
- template
- void setXYContainers(const ScalarContainerX &x,
- const ScalarContainerY &y,
- Scalar m0, Scalar m1)
- {
- assert(x.size() == y.size());
- assert(x.size() == numSamples());
-
- std::copy(x.begin(), x.end(), xPos_.begin());
- std::copy(y.begin(), y.end(), yPos_.begin());
- makeFullSpline_(m0, m1);
- }
-
- /*!
- * \brief Set the sampling points and the boundary slopes of a
- * full spline using a C-style array.
- *
- * This method uses a single array of sampling points, which are
- * seen as an array-like object which provides access to the X and
- * Y coordinates. In this context 'array-like' means that an
- * access to the members is provided via the [] operator. (e.g. C
- * arrays, std::vector, std::array, etc.) The array containing
- * the sampling points must be of size 'numberSamples' at least. Also,
- * the number of sampling points must be larger than 1.
- */
- template
- void setArrayOfPoints(size_t numberSamples,
- const PointArray &points,
- Scalar m0,
- Scalar m1)
- {
- // a spline with no or just one sampling points? what an
- // incredible bad idea!
- assert(numberSamples == numSamples());
-
- for (size_t i = 0; i < numberSamples; ++i) {
- xPos_[i] = points[i][0];
- yPos_[i] = points[i][1];
- }
- makeFullSpline_(m0, m1);
- }
-
- /*!
- * \brief Set the sampling points and the boundary slopes of a
- * full spline using a STL-compatible container of
- * array-like objects.
- *
- * This method uses a single STL-compatible container of sampling
- * points, which are assumed to be array-like objects storing the
- * X and Y coordinates. "STL-compatible" means that the container
- * provides access to iterators using the begin(), end() methods
- * and also provides a size() method. Also, the number of entries
- * in the X and the Y containers must be equal and larger than 1.
- */
- template
- void setContainerOfPoints(const XYContainer &points,
- Scalar m0,
- Scalar m1)
- {
- assert(points.size() == numSamples());
-
- typename XYContainer::const_iterator it = points.begin();
- typename XYContainer::const_iterator endIt = points.end();
- for (int i = 0; it != endIt; ++i, ++it) {
- xPos_[i] = (*it)[0];
- yPos_[i] = (*it)[1];
- }
-
- // make a full spline
- makeFullSpline_(m0, m1);
- }
-
- /*!
- * \brief Set the sampling points and the boundary slopes of a
- * full spline using a STL-compatible container of
- * tuple-like objects.
- *
- * This method uses a single STL-compatible container of sampling
- * points, which are assumed to be tuple-like objects storing the
- * X and Y coordinates. "tuple-like" means that the objects
- * provide access to the x values via std::get<0>(obj) and to the
- * y value via std::get<1>(obj) (e.g. std::tuple or
- * std::pair). "STL-compatible" means that the container provides
- * access to iterators using the begin(), end() methods and also
- * provides a size() method. Also, the number of entries in the X
- * and the Y containers must be equal and larger than 1.
- */
- template
- void setContainerOfTuples(const XYContainer &points,
- Scalar m0,
- Scalar m1)
- {
- assert(points.size() == numSamples());
-
- typename XYContainer::const_iterator it = points.begin();
- typename XYContainer::const_iterator endIt = points.end();
- for (int i = 0; it != endIt; ++i, ++it) {
- xPos_[i] = std::get<0>(*it);
- yPos_[i] = std::get<1>(*it);
- }
-
- // make a full spline
- makeFullSpline_(m0, m1);
- }
-
- ///////////////////////////////////////
- ///////////////////////////////////////
- ///////////////////////////////////////
- // Natural splines //
- ///////////////////////////////////////
- ///////////////////////////////////////
- ///////////////////////////////////////
- /*!
- * \brief Set the sampling points natural spline using C-style arrays.
- *
- * This method uses separate array-like objects for the values of
- * the X and Y coordinates. In this context 'array-like' means
- * that an access to the members is provided via the []
- * operator. (e.g. C arrays, std::vector, std::array, etc.) Each
- * array must be of size 'numberSamples' at least. Also, the number of
- * sampling points must be larger than 1.
- */
- template
- void setXYArrays(size_t numberSamples,
- const ScalarArrayX &x,
- const ScalarArrayY &y)
- {
- assert(numberSamples == numSamples());
-
- for (size_t i = 0; i < numberSamples; ++i) {
- xPos_[i] = x[i];
- yPos_[i] = y[i];
- }
-
- makeNaturalSpline_();
- }
-
- /*!
- * \brief Set the sampling points of a natural spline using
- * STL-compatible containers.
- *
- * This method uses separate STL-compatible containers for the
- * values of the sampling points' X and Y
- * coordinates. "STL-compatible" means that the container provides
- * access to iterators using the begin(), end() methods and also
- * provides a size() method. Also, the number of entries in the X
- * and the Y containers must be equal and larger than 1.
- */
- template
- void setXYContainers(const ScalarContainerX &x,
- const ScalarContainerY &y)
- {
- assert(x.size() == y.size());
- assert(x.size() > 1);
-
- setNumSamples_(x.size());
- std::copy(x.begin(), x.end(), xPos_.begin());
- std::copy(y.begin(), y.end(), yPos_.begin());
- makeNaturalSpline_();
- }
-
-
- /*!
- * \brief Set the sampling points of a natural spline using a
- * C-style array.
- *
- * This method uses a single array of sampling points, which are
- * seen as an array-like object which provides access to the X and
- * Y coordinates. In this context 'array-like' means that an
- * access to the members is provided via the [] operator. (e.g. C
- * arrays, std::vector, std::array, etc.) The array containing
- * the sampling points must be of size 'numberSamples' at least. Also,
- * the number of sampling points must be larger than 1.
- */
- template
- void setArrayOfPoints(int numberSamples,
- const PointArray &points)
- {
- assert(numberSamples == numSamples());
-
- for (int i = 0; i < numberSamples; ++i) {
- xPos_[i] = points[i][0];
- yPos_[i] = points[i][1];
- }
- makeNaturalSpline_();
- }
-
- /*!
- * \brief Set the sampling points of a natural spline using a
- * STL-compatible container of array-like objects.
- *
- * This method uses a single STL-compatible container of sampling
- * points, which are assumed to be array-like objects storing the
- * X and Y coordinates. "STL-compatible" means that the container
- * provides access to iterators using the begin(), end() methods
- * and also provides a size() method. Also, the number of entries
- * in the X and the Y containers must be equal and larger than 1.
- */
- template
- void setContainerOfPoints(const XYContainer &points)
- {
- assert(points.size() == numSamples());
-
- typename XYContainer::const_iterator it = points.begin();
- typename XYContainer::const_iterator endIt = points.end();
- for (int i = 0; it != endIt; ++ i, ++it) {
- xPos_[i] = (*it)[0];
- yPos_[i] = (*it)[1];
- }
-
- // make a natural spline
- makeNaturalSpline_();
- }
-
- /*!
- * \brief Set the sampling points of a natural spline using a
- * STL-compatible container of tuple-like objects.
- *
- * This method uses a single STL-compatible container of sampling
- * points, which are assumed to be tuple-like objects storing the
- * X and Y coordinates. "tuple-like" means that the objects
- * provide access to the x values via std::get<0>(obj) and to the
- * y value via std::get<1>(obj) (e.g. std::tuple or
- * std::pair). "STL-compatible" means that the container provides
- * access to iterators using the begin(), end() methods and also
- * provides a size() method. Also, the number of entries in the X
- * and the Y containers must be equal and larger than 1.
- */
- template
- void setContainerOfTuples(const XYContainer &points)
- {
- assert(points.size() == numSamples());
-
- typename XYContainer::const_iterator it = points.begin();
- typename XYContainer::const_iterator endIt = points.end();
- for (int i = 0; it != endIt; ++i, ++it) {
- xPos_[i] = std::get<0>(*it);
- yPos_[i] = std::get<1>(*it);
- }
-
- // make a natural spline
- makeNaturalSpline_();
- }
-
-protected:
- /*!
- * \brief Create a full spline from the already set sampling points.
- *
- * Also creates temporary matrix and right hand side vector.
- */
- void makeFullSpline_(Scalar m0, Scalar m1)
- {
- BTDMatrix M(numSamples());
- BlockVector d(numSamples());
-
- // create linear system of equations
- this->makeFullSystem_(M, d, m0, m1);
-
- // solve for the moments
- M.solve(m_, d);
- }
-
- /*!
- * \brief Create a natural spline from the already set sampling points.
- *
- * Also creates temporary matrix and right hand side vector.
- */
- void makeNaturalSpline_()
- {
- BTDMatrix M(numSamples());
- BlockVector d(numSamples());
-
- // create linear system of equations
- this->makeNaturalSystem_(M, d);
-
- // solve for the moments
- M.solve(m_, d);
- }
-
- /*!
- * \brief Returns the x coordinate of the i-th sampling point.
- */
- Scalar x_(int i) const
- { return xPos_[i]; }
-
- /*!
- * \brief Returns the y coordinate of the i-th sampling point.
- */
- Scalar y_(int i) const
- { return yPos_[i]; }
-
- /*!
- * \brief Returns the moment (i.e. second derivative) of the
- * spline at the i-th sampling point.
- */
- Scalar moment_(int i) const
- { return m_[i]; }
-
- Vector xPos_;
- Vector yPos_;
- BlockVector m_;
-};
-
-}
-
-#endif
diff --git a/opm/common/math.hh b/opm/common/math.hh
deleted file mode 100644
index 55b2376d2..000000000
--- a/opm/common/math.hh
+++ /dev/null
@@ -1,392 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * Copyright (C) 2010-2012 by Andreas Lauser *
- * *
- * 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, see . *
- *****************************************************************************/
-/*!
- * \file
- * \brief Define some often used mathematical functions
- */
-#ifndef OPM_MATH_HH
-#define OPM_MATH_HH
-
-#include
-#include
-
-#include
-#include
-
-namespace Opm
-{
-/*!
- * \ingroup Math
- * \brief Calculate the harmonic mean of two scalar values.
- *
- * \param x The first input value
- * \param y The second input value
- */
-template
-Scalar harmonicMean(Scalar x, Scalar y)
-{
- if (x*y <= 0)
- return 0;
- return (2*x*y)/(x + y);
-}
-
-/*!
- * \ingroup Math
- * \brief Calculate the geometric mean of two scalar values.
- *
- * \param x The first input value
- * \param y The second input value
- */
-template
-Scalar geometricMean(Scalar x, Scalar y)
-{
- if (x*y <= 0)
- return 0;
- return std::sqrt(x*y)*((x < 0)?-1:1);
-}
-
-/*!
- * \ingroup Math
- * \brief Calculate the harmonic mean of a fixed-size matrix.
- *
- * This is done by calculating the harmonic mean for each entry
- * individually. The result is stored the first argument.
- *
- * \param K The matrix for storing the result
- * \param Ki The first input matrix
- * \param Kj The second input matrix
- */
-template
-void harmonicMeanMatrix(Dune::FieldMatrix &K,
- const Dune::FieldMatrix &Ki,
- const Dune::FieldMatrix &Kj)
-{
- for (int rowIdx=0; rowIdx < m; rowIdx++){
- for (int colIdx=0; colIdx< n; colIdx++){
- if (Ki[rowIdx][colIdx] != Kj[rowIdx][colIdx]) {
- K[rowIdx][colIdx] =
- harmonicMean(Ki[rowIdx][colIdx],
- Kj[rowIdx][colIdx]);
- }
- else
- K[rowIdx][colIdx] = Ki[rowIdx][colIdx];
- }
- }
-}
-
-/*!
- * \ingroup Math
- * \brief Invert a linear polynomial analytically
- *
- * The polynomial is defined as
- * \f[ p(x) = a\; x + b \f]
- *
- * This method Returns the number of solutions which are in the real
- * numbers, i.e. 1 except if the slope of the line is 0.
- *
- * \param sol Container into which the solutions are written
- * \param a The coefficient for the linear term
- * \param b The coefficient for the constant term
- */
-template
-int invertLinearPolynomial(SolContainer &sol,
- Scalar a,
- Scalar b)
-{
- if (a == 0.0)
- return 0;
-
- sol[0] = -b/a;
- return 1;
-}
-
-/*!
- * \ingroup Math
- * \brief Invert a quadratic polynomial analytically
- *
- * The polynomial is defined as
- * \f[ p(x) = a\; x^2 + + b\;x + c \f]
- *
- * This method teturns the number of solutions which are in the real
- * numbers. The "sol" argument contains the real roots of the parabola
- * in order with the smallest root first.
- *
- * \param sol Container into which the solutions are written
- * \param a The coefficient for the quadratic term
- * \param b The coefficient for the linear term
- * \param c The coefficient for the constant term
- */
-template
-int invertQuadraticPolynomial(SolContainer &sol,
- Scalar a,
- Scalar b,
- Scalar c)
-{
- // check for a line
- if (a == 0.0)
- return invertLinearPolynomial(sol, b, c);
-
- // discriminant
- Scalar Delta = b*b - 4*a*c;
- if (Delta < 0)
- return 0; // no real roots
-
- Delta = std::sqrt(Delta);
- sol[0] = (- b + Delta)/(2*a);
- sol[1] = (- b - Delta)/(2*a);
-
- // sort the result
- if (sol[0] > sol[1])
- std::swap(sol[0], sol[1]);
- return 2; // two real roots
-}
-
-//! \cond SKIP_THIS
-template
-void invertCubicPolynomialPostProcess_(SolContainer &sol,
- int numSol,
- Scalar a,
- Scalar b,
- Scalar c,
- Scalar d)
-{
- // do one Newton iteration on the analytic solution if the
- // precision is increased
- for (int i = 0; i < numSol; ++i) {
- Scalar x = sol[i];
- Scalar fOld = d + x*(c + x*(b + x*a));
-
- Scalar fPrime = c + x*(2*b + x*3*a);
- if (fPrime == 0.0)
- continue;
- x -= fOld/fPrime;
-
- Scalar fNew = d + x*(c + x*(b + x*a));
- if (std::abs(fNew) < std::abs(fOld))
- sol[i] = x;
- }
-}
-//! \endcond
-
-/*!
- * \ingroup Math
- * \brief Invert a cubic polynomial analytically
- *
- * The polynomial is defined as
- * \f[ p(x) = a\; x^3 + + b\;x^3 + c\;x + d \f]
- *
- * This method teturns the number of solutions which are in the real
- * numbers. The "sol" argument contains the real roots of the cubic
- * polynomial in order with the smallest root first.
- *
- * \param sol Container into which the solutions are written
- * \param a The coefficient for the cubic term
- * \param b The coefficient for the quadratic term
- * \param c The coefficient for the linear term
- * \param d The coefficient for the constant term
- */
-template
-int invertCubicPolynomial(SolContainer *sol,
- Scalar a,
- Scalar b,
- Scalar c,
- Scalar d)
-{
- // reduces to a quadratic polynomial
- if (a == 0)
- return invertQuadraticPolynomial(sol, b, c, d);
-
- // normalize the polynomial
- b /= a;
- c /= a;
- d /= a;
- a = 1;
-
- // get rid of the quadratic term by subsituting x = t - b/3
- Scalar p = c - b*b/3;
- Scalar q = d + (2*b*b*b - 9*b*c)/27;
-
- if (p != 0.0 && q != 0.0) {
- // At this point
- //
- // t^3 + p*t + q = 0
- //
- // with p != 0 and q != 0 holds. Introducing the variables u and v
- // with the properties
- //
- // u + v = t and 3*u*v + p = 0
- //
- // leads to
- //
- // u^3 + v^3 + q = 0 .
- //
- // multiplying both sides with u^3 and taking advantage of the
- // fact that u*v = -p/3 leads to
- //
- // u^6 + q*u^3 - p^3/27 = 0
- //
- // Now, substituting u^3 = w yields
- //
- // w^2 + q*w - p^3/27 = 0
- //
- // This is a quadratic equation with the solutions
- //
- // w = -q/2 + sqrt(q^2/4 + p^3/27) and
- // w = -q/2 - sqrt(q^2/4 + p^3/27)
- //
- // Since w is equivalent to u^3 it is sufficient to only look at
- // one of the two cases. Then, there are still 2 cases: positive
- // and negative discriminant.
- Scalar wDisc = q*q/4 + p*p*p/27;
- if (wDisc >= 0) { // the positive discriminant case:
- // calculate the cube root of - q/2 + sqrt(q^2/4 + p^3/27)
- Scalar u = - q/2 + std::sqrt(wDisc);
- if (u < 0) u = - std::pow(-u, 1.0/3);
- else u = std::pow(u, 1.0/3);
-
- // at this point, u != 0 since p^3 = 0 is necessary in order
- // for u = 0 to hold, so
- sol[0] = u - p/(3*u) - b/3;
- // does not produce a division by zero. the remaining two
- // roots of u are rotated by +- 2/3*pi in the complex plane
- // and thus not considered here
- invertCubicPolynomialPostProcess_(sol, 1, a, b, c, d);
- return 1;
- }
- else { // the negative discriminant case:
- Scalar uCubedRe = - q/2;
- Scalar uCubedIm = std::sqrt(-wDisc);
- // calculate the cube root of - q/2 + sqrt(q^2/4 + p^3/27)
- Scalar uAbs = std::pow(std::sqrt(uCubedRe*uCubedRe + uCubedIm*uCubedIm), 1.0/3);
- Scalar phi = std::atan2(uCubedIm, uCubedRe)/3;
-
- // calculate the length and the angle of the primitive root
-
- // with the definitions from above it follows that
- //
- // x = u - p/(3*u) - b/3
- //
- // where x and u are complex numbers. Rewritten in polar form
- // this is equivalent to
- //
- // x = |u|*e^(i*phi) - p*e^(-i*phi)/(3*|u|) - b/3 .
- //
- // Factoring out the e^ terms and subtracting the additional
- // terms, yields
- //
- // x = (e^(i*phi) + e^(-i*phi))*(|u| - p/(3*|u|)) - y - b/3
- //
- // with
- //
- // y = - |u|*e^(-i*phi) + p*e^(i*phi)/(3*|u|) .
- //
- // The crucial observation is the fact that y is the conjugate
- // of - x + b/3. This means that after taking advantage of the
- // relation
- //
- // e^(i*phi) + e^(-i*phi) = 2*cos(phi)
- //
- // the equation
- //
- // x = 2*cos(phi)*(|u| - p / (3*|u|)) - conj(x) - 2*b/3
- //
- // holds. Since |u|, p, b and cos(phi) are real numbers, it
- // follows that Im(x) = - Im(x) and thus Im(x) = 0. This
- // implies
- //
- // Re(x) = x = cos(phi)*(|u| - p / (3*|u|)) - b/3 .
- //
- // Considering the fact that u is a cubic root, we have three
- // values for phi which differ by 2/3*pi. This allows to
- // calculate the three real roots of the polynomial:
- for (int i = 0; i < 3; ++i) {
- sol[i] = std::cos(phi)*(uAbs - p/(3*uAbs)) - b/3;
- phi += 2*M_PI/3;
- }
-
- // post process the obtained solution to increase numerical
- // precision
- invertCubicPolynomialPostProcess_(sol, 3, a, b, c, d);
-
- // sort the result
- std::sort(sol, sol + 3);
-
- return 3;
- }
- }
- // Handle some (pretty unlikely) special cases required to avoid
- // divisions by zero in the code above...
- else if (p == 0.0 && q == 0.0) {
- // t^3 = 0, i.e. triple root at t = 0
- sol[0] = sol[1] = sol[2] = 0.0 - b/3;
- return 3;
- }
- else if (p == 0.0 && q != 0.0) {
- // t^3 + q = 0,
- //
- // i. e. single real root at t=curt(q)
- Scalar t;
- if (-q > 0) t = std::pow(-q, 1./3);
- else t = - std::pow(q, 1./3);
- sol[0] = t - b/3;
-
- return 1;
- }
-
- assert(p != 0.0 && q == 0.0);
-
- // t^3 + p*t = 0 = t*(t^2 + p),
- //
- // i. e. roots at t = 0, t^2 + p = 0
- if (p > 0) {
- sol[0] = 0.0 - b/3;
- return 1; // only a single real root at t=0
- }
-
- // two additional real roots at t = sqrt(-p) and t = -sqrt(-p)
- sol[0] = -std::sqrt(-p) - b/3;;
- sol[1] = 0.0 - b/3;
- sol[2] = std::sqrt(-p) - b/3;
-
- return 3;
-}
-
-/*!
- * \brief Cross product of two vectors in three-dimensional Euclidean space
- *
- * \param vec1 The first vector
- * \param vec2 The second vector
- */
-template
-Dune::FieldVector crossProduct(const Dune::FieldVector &vec1,
- const Dune::FieldVector &vec2)
-{
- Dune::FieldVector result;
-
- result[0] = vec1[1]*vec2[2] - vec1[2]*vec2[1];
- result[1] = vec1[2]*vec2[0] - vec1[0]*vec2[2];
- result[2] = vec1[0]*vec2[1] - vec1[1]*vec2[0];
-
- return result;
-}
-
-}
-
-
-#endif
diff --git a/opm/common/spline.hh b/opm/common/spline.hh
deleted file mode 100644
index 60f6d25ed..000000000
--- a/opm/common/spline.hh
+++ /dev/null
@@ -1,532 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * Copyright (C) 2010-2012 by Andreas Lauser *
- * *
- * 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, see . *
- *****************************************************************************/
-/*!
- * \file
- * \copydoc Opm::Spline
- */
-#ifndef OPM_SPLINE_HH
-#define OPM_SPLINE_HH
-
-#include "fixedlengthspline_.hh"
-#include "variablelengthspline_.hh"
-#include "splinecommon_.hh"
-
-#include
-#include
-
-namespace Opm
-{
-/*!
- * \ingroup Spline
- * \brief A 3rd order polynomial spline.
-
- * This class implements a spline \f$s(x)\f$ for which, given \f$n\f$ sampling
- * points \f$x_1, \dots, x_n\f$, the following conditions hold
- *\f{align*}{
- s(x_i) & = y_i \quad \forall i \in \{1, \dots, n \}\\
- s'(x_1) & = m_1 \\
- s'(x_n) & = m_n
- \f}
-*
-* for any given boundary slopes \f$m_1\f$ and \f$m_n\f$. Alternatively, natural
-* splines are supported which are defined by
-*\f{align*}{
- s(x_i) & = y_i \quad \forall i \in \{1, \dots, n \} \\
- s''(x_1) & = 0 \\
- s''(x_n) & = 0
-\f}
- */
-template
-class Spline : public FixedLengthSpline_
-{
-public:
- /*!
- * \brief Default constructor for a spline.
- *
- * To specfiy the acutal curve, use one of the set() methods.
- */
- Spline()
- { }
-
- /*!
- * \brief Convenience constructor for a full spline
- *
- * \param x An array containing the \f$x\f$ values of the spline's sampling points
- * \param y An array containing the \f$y\f$ values of the spline's sampling points
- */
- template
- Spline(const ScalarArray &x,
- const ScalarArray &y)
- { this->setXYArrays(numSamples, x, y); }
-
- /*!
- * \brief Convenience constructor for a full spline
- *
- * \param points An array of \f$(x,y)\f$ tuples of the spline's sampling points
- */
- template
- Spline(const PointArray &points)
- { this->setArrayOfPoints(numSamples, points); }
-
- /*!
- * \brief Convenience constructor for a full spline
- *
- * \param x An array containing the \f$x\f$ values of the spline's sampling points
- * \param y An array containing the \f$y\f$ values of the spline's sampling points
- * \param m0 The slope of the spline at \f$x_0\f$
- * \param m1 The slope of the spline at \f$x_n\f$
- */
- template
- Spline(const ScalarArray &x,
- const ScalarArray &y,
- Scalar m0,
- Scalar m1)
- { this->setXYArrays(numSamples, x, y, m0, m1); }
-
- /*!
- * \brief Convenience constructor for a full spline
- *
- * \param points An array of \f$(x,y)\f$ tuples of the spline's sampling points
- * \param m0 The slope of the spline at \f$x_0\f$
- * \param m1 The slope of the spline at \f$x_n\f$
- */
- template
- Spline(const PointArray &points,
- Scalar m0,
- Scalar m1)
- { this->setArrayOfPoints(numSamples, points, m0, m1); }
-};
-
-/*!
- * \brief Specialization of a spline with the number of sampling
- * points only known at run time.
- *
- * This class implements a spline \f$s(x)\f$ for which, given \f$n\f$ sampling
- * points \f$x_1, \dots, x_n\f$, the following conditions hold
- *\f{align*}{
- s(x_i) & = y_i \quad \forall i \in \{1, \dots, n \}\\
- s'(x_1) & = m_1 \\
- s'(x_n) & = m_n
- \f}
- *
- * for any given boundary slopes \f$m_1\f$ and \f$m_n\f$. Alternatively, natural
- * splines are supported which are defined by
- *\f{align*}{
- s(x_i) & = y_i \quad \forall i \in \{1, \dots, n \} \\
- s''(x_1) & = 0 \\
- s''(x_n) & = 0
- \f}
-*/
-template
-class Spline : public VariableLengthSpline_
-{
-public:
- /*!
- * \brief Default constructor for a spline.
- *
- * To specfiy the acutal curve, use one of the set() methods.
- */
- Spline()
- { }
-
- /*!
- * \brief Convenience constructor for a natural spline
- *
- * \param nSamples The number of sampling points (must be > 2)
- * \param x An array containing the \f$x\f$ values of the spline's sampling points
- * \param y An array containing the \f$y\f$ values of the spline's sampling points
- */
- template
- Spline(int nSamples,
- const ScalarArrayX &x,
- const ScalarArrayY &y)
- { this->setXYArrays(nSamples, x, y); }
-
- /*!
- * \brief Convenience constructor for a natural spline
- *
- * \param nSamples The number of sampling points (must be > 2)
- * \param points An array of \f$(x,y)\f$ tuples of the spline's sampling points
- */
- template
- Spline(int nSamples,
- const PointArray &points)
- { this->setArrayOfPoints(nSamples, points); }
-
- /*!
- * \brief Convenience constructor for a natural spline
- *
- * \param x An array containing the \f$x\f$ values of the spline's sampling points (must have a size() method)
- * \param y An array containing the \f$y\f$ values of the spline's sampling points (must have a size() method)
- */
- template
- Spline(const ScalarContainer &x,
- const ScalarContainer &y)
- { this->setXYContainers(x, y); }
-
- /*!
- * \brief Convenience constructor for a natural spline
- *
- * \param points An array of \f$(x,y)\f$ tuples of the spline's sampling points (must have a size() method)
- */
- template
- Spline(const PointContainer &points)
- { this->setContainerOfPoints(points); }
-
- /*!
- * \brief Convenience constructor for a full spline
- *
- * \param nSamples The number of sampling points (must be >= 2)
- * \param x An array containing the \f$x\f$ values of the spline's sampling points
- * \param y An array containing the \f$y\f$ values of the spline's sampling points
- * \param m0 The slope of the spline at \f$x_0\f$
- * \param m1 The slope of the spline at \f$x_n\f$
- */
- template
- Spline(int nSamples,
- const ScalarArray &x,
- const ScalarArray &y,
- Scalar m0,
- Scalar m1)
- { this->setXYArrays(nSamples, x, y, m0, m1); }
-
- /*!
- * \brief Convenience constructor for a full spline
- *
- * \param nSamples The number of sampling points (must be >= 2)
- * \param points An array containing the \f$x\f$ and \f$x\f$ values of the spline's sampling points
- * \param m0 The slope of the spline at \f$x_0\f$
- * \param m1 The slope of the spline at \f$x_n\f$
- */
- template
- Spline(int nSamples,
- const PointArray &points,
- Scalar m0,
- Scalar m1)
- { this->setArrayOfPoints(nSamples, points, m0, m1); }
-
- /*!
- * \brief Convenience constructor for a full spline
- *
- * \param x An array containing the \f$x\f$ values of the spline's sampling points (must have a size() method)
- * \param y An array containing the \f$y\f$ values of the spline's sampling points (must have a size() method)
- * \param m0 The slope of the spline at \f$x_0\f$
- * \param m1 The slope of the spline at \f$x_n\f$
- */
- template
- Spline(const ScalarContainerX &x,
- const ScalarContainerY &y,
- Scalar m0,
- Scalar m1)
- { this->setXYContainers(x, y, m0, m1); }
-
- /*!
- * \brief Convenience constructor for a full spline
- *
- * \param points An array of \f$(x,y)\f$ tuples of the spline's sampling points (must have a size() method)
- * \param m0 The slope of the spline at \f$x_0\f$
- * \param m1 The slope of the spline at \f$x_n\f$
- */
- template
- Spline(const PointContainer &points,
- Scalar m0,
- Scalar m1)
- { this->setContainerOfPoints(points, m0, m1); }
-};
-
-/*!
- * \brief Do not allow splines with zero sampling points
- */
-template
-class Spline
-// Splines with zero sampling points do not make sense!
-{ private: Spline() { } };
-
-/*!
- * \brief Do not allow splines with one sampling point
- */
-template
-class Spline
-// Splines with one sampling point do not make sense!
-{ private: Spline() { } };
-
-/*!
- * \brief Spline for two sampling points.
- *
- * For this type of spline there is no natural spline.
- */
-template
-class Spline : public SplineCommon_ >
-{
- friend class SplineCommon_ >;
- typedef Dune::FieldVector Vector;
- typedef Dune::FieldMatrix Matrix;
-
-public:
- Spline()
- {}
-
- /*!
- * \brief Convenience constructor for a full spline
- *
- * \param x An array containing the \f$x\f$ values of the spline's sampling points
- * \param y An array containing the \f$y\f$ values of the spline's sampling points
- * \param m0 The slope of the spline at \f$x_0\f$
- * \param m1 The slope of the spline at \f$x_n\f$
- */
- template
- Spline(const ScalarArrayX &x,
- const ScalarArrayY &y,
- Scalar m0, Scalar m1)
- { setXYArrays(2, x, y, m0, m1); }
-
- /*!
- * \brief Convenience constructor for a full spline
- *
- * \param points An array of \f$(x,y)\f$ tuples of the spline's sampling points
- * \param m0 The slope of the spline at \f$x_0\f$
- * \param m1 The slope of the spline at \f$x_n\f$
- */
- template
- Spline(const PointArray &points,
- Scalar m0,
- Scalar m1)
- { this->setArrayOfPoints(2, points, m0, m1); }
-
- /*!
- * \brief Convenience constructor for a full spline
- *
- * \param x0 The \f$x\f$ value of the first sampling point
- * \param x1 The \f$x\f$ value of the second sampling point
- * \param y0 The \f$y\f$ value of the first sampling point
- * \param y1 The \f$y\f$ value of the second sampling point
- * \param m0 The slope of the spline at \f$x_0\f$
- * \param m1 The slope of the spline at \f$x_n\f$
- */
- Spline(Scalar x0, Scalar x1,
- Scalar y0, Scalar y1,
- Scalar m0, Scalar m1)
- {
- set(x0, x1,
- y0, y1,
- m0, m1);
- }
-
- /*!
- * \brief Returns the number of sampling points.
- */
- int numSamples() const
- { return 2; }
-
- /*!
- * \brief Set the sampling points and the boundary slopes of the
- * spline.
- *
- * \param x0 The \f$x\f$ value of the first sampling point
- * \param x1 The \f$x\f$ value of the second sampling point
- * \param y0 The \f$y\f$ value of the first sampling point
- * \param y1 The \f$y\f$ value of the second sampling point
- * \param m0 The slope of the spline at \f$x_0\f$
- * \param m1 The slope of the spline at \f$x_1\f$
- */
- void set(Scalar x0, Scalar x1,
- Scalar y0, Scalar y1,
- Scalar m0, Scalar m1)
- {
- Matrix M(numSamples());
- Vector d;
- assignXY_(x0, x1, y0, y1);
- this->makeFullSystem_(M, d, m0, m1);
-
- // solve for the moments
- M.solve(m_, d);
- }
-
- /*!
- * \brief Set the sampling points and the boundary slopes of the
- * spline.
- *
- * \param nSamples The number of sampling points (must be >= 2)
- * \param x An array containing the \f$x\f$ values of the sampling points
- * \param y An array containing the \f$y\f$ values of the sampling points
- * \param m0 The slope of the spline at \f$x_0\f$
- * \param m1 The slope of the spline at \f$x_1\f$
- */
- template
- void setXYArrays(int nSamples,
- const ScalarContainer &x,
- const ScalarContainer &y,
- Scalar m0, Scalar m1)
- {
- assert(nSamples == 2);
- set(x[0], x[1], y[0], y[1], m0, m1);
- }
-
- /*!
- * \brief Set the sampling points and the boundary slopes of the
- * spline.
- *
- * \param x An array containing the \f$x\f$ values of the sampling points
- * \param y An array containing the \f$y\f$ values of the sampling points
- * \param m0 The slope of the spline at \f$x_0\f$
- * \param m1 The slope of the spline at \f$x_1\f$
- */
- template
- void setXYContainers(const ScalarContainerX &x,
- const ScalarContainerY &y,
- Scalar m0, Scalar m1)
- {
- assert(x.size() == y.size());
- assert(x.size() == 2);
-
- Matrix M(numSamples());
- Vector d;
-
- typename ScalarContainerX::const_iterator xIt0 = x.begin();
- typename ScalarContainerX::const_iterator xIt1 = xIt0;
- ++xIt1;
- typename ScalarContainerY::const_iterator yIt0 = y.begin();
- typename ScalarContainerY::const_iterator yIt1 = yIt0;
- ++yIt1;
- set(*xIt0, *xIt1, *yIt0, *yIt1);
- }
-
- /*!
- * \brief Set the sampling points and the boundary slopes of the
- * spline.
- *
- * \param nSamples The number of sampling points (must be >= 2)
- * \param points An array of \f$(x,y)\f$ tuples of the spline's sampling points
- * \param m0 The slope of the spline at \f$x_0\f$
- * \param m1 The slope of the spline at \f$x_1\f$
- */
- template
- void setArrayOfPoints(int nSamples,
- const PointArray &points,
- Scalar m0,
- Scalar m1)
- {
- assert(nSamples == 2);
-
- set(points[0][0],
- points[1][0],
- points[0][1],
- points[1][1],
- m0, m1);
- }
-
- /*!
- * \brief Set the sampling points and the boundary slopes from an
- * STL-like container of points.
- *
- * \param points An array of \f$(x,y)\f$ tuples of the spline's sampling points
- * \param m0 The slope of the spline at \f$x_0\f$
- * \param m1 The slope of the spline at \f$x_1\f$
- */
- template
- void setContainerOfPoints(const PointContainer &points,
- Scalar m0,
- Scalar m1)
- {
- assert(points.size() == 2);
-
- Matrix M;
- Vector d;
- typename PointContainer::const_iterator it0 = points.begin();
- typename PointContainer::const_iterator it1 = it0;
- ++it1;
-
- set((*it0)[0],
- (*it0)[1],
- (*it1)[0],
- (*it1)[1],
- m0, m1);
- }
-
- /*!
- * \brief Set the sampling points and the boundary slopes from an
- * STL-like container of tuples.
- *
- * \param tuples An array of \f$(x,y)\f$ tuples of the spline's sampling points
- * \param m0 The slope of the spline at \f$x_0\f$
- * \param m1 The slope of the spline at \f$x_1\f$
- */
- template
- void setContainerOfTuples(const TupleContainer &tuples,
- Scalar m0,
- Scalar m1)
- {
- assert(tuples.size() == 2);
-
- typename TupleContainer::const_iterator it0 = tuples.begin();
- typename TupleContainer::const_iterator it1 = it0;
- ++it1;
-
- set(std::get<0>(*it0),
- std::get<1>(*it0),
- std::get<0>(*it1),
- std::get<1>(*it1),
- m0, m1);
- }
-
-protected:
- void assignXY_(Scalar x0, Scalar x1,
- Scalar y0, Scalar y1)
- {
- if (x0 > x1) {
- xPos_[0] = x1;
- xPos_[1] = x0;
- yPos_[0] = y1;
- yPos_[1] = y0;
- }
- else {
- xPos_[0] = x0;
- xPos_[1] = x1;
- yPos_[0] = y0;
- yPos_[1] = y1;
- }
- }
-
- /*!
- * \brief Returns the x coordinate of the i-th sampling point.
- */
- Scalar x_(int i) const
- { return xPos_[i]; }
-
- /*!
- * \brief Returns the y coordinate of the i-th sampling point.
- */
- Scalar y_(int i) const
- { return yPos_[i]; }
-
- /*!
- * \brief Returns the moment (i.e. second derivative) of the
- * spline at the i-th sampling point.
- */
- Scalar moment_(int i) const
- { return m_[i]; }
-
- Vector xPos_;
- Vector yPos_;
- Vector m_;
-};
-
-}
-
-#endif
diff --git a/opm/common/splinecommon_.hh b/opm/common/splinecommon_.hh
deleted file mode 100644
index 11f79e15b..000000000
--- a/opm/common/splinecommon_.hh
+++ /dev/null
@@ -1,685 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * Copyright (C) 2010-2012 by Andreas Lauser *
- * *
- * 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, see . *
- *****************************************************************************/
-/*!
- * \file
- * \copydoc Opm::SplineCommon_
- */
-#ifndef OPM_SPLINE_COMMON__HH
-#define OPM_SPLINE_COMMON__HH
-
-#include "valgrind.hh"
-#include "math.hh"
-
-#include
-
-#include
-#include
-#include
-#include
-
-namespace Opm
-{
-/*!
- * \brief The common code for all 3rd order polynomial splines.
- */
-template
-class SplineCommon_
-{
- typedef ScalarT Scalar;
- typedef ImplementationT Implementation;
-
- Implementation &asImp_()
- { return *static_cast(this); }
- const Implementation &asImp_() const
- { return *static_cast(this); }
-
-public:
- /*!
- * \brief Return true iff the given x is in range [x1, xn].
- */
- bool applies(Scalar x) const
- {
- return x_(0) <= x && x <= x_(numSamples_() - 1);
- }
-
- /*!
- * \brief Return the x value of the leftmost sampling point.
- */
- Scalar xMin() const
- { return x_(0); }
-
- /*!
- * \brief Return the x value of the rightmost sampling point.
- */
- Scalar xMax() const
- { return x_(numSamples_() - 1); }
-
- /*!
- * \brief Prints k tuples of the format (x, y, dx/dy, isMonotonic)
- * to stdout.
- *
- * If the spline does not apply for parts of [x0, x1] it is
- * extrapolated using a straight line. The result can be inspected
- * using the following commands:
- *
- ----------- snip -----------
- ./yourProgramm > spline.csv
- gnuplot
-
- gnuplot> plot "spline.csv" using 1:2 w l ti "Curve", \
- "spline.csv" using 1:3 w l ti "Derivative", \
- "spline.csv" using 1:4 w p ti "Monotonic"
- ----------- snap -----------
- */
- void printCSV(Scalar xi0, Scalar xi1, int k) const
- {
- Scalar x0 = std::min(xi0, xi1);
- Scalar x1 = std::max(xi0, xi1);
- const int n = numSamples_() - 1;
- for (int i = 0; i <= k; ++i) {
- double x = i*(x1 - x0)/k + x0;
- double x_p1 = x + (x1 - x0)/k;
- double y;
- double dy_dx;
- double mono = 1;
- if (!applies(x)) {
- if (x < 0) {
- dy_dx = evalDerivative(x_(0));
- y = (x - x_(0))*dy_dx + y_(0);
- mono = (dy_dx>0)?1:-1;
- }
- else if (x > x_(n)) {
- dy_dx = evalDerivative(x_(n));
- y = (x - x_(n))*dy_dx + y_(n);
- mono = (dy_dx>0)?1:-1;
- }
- else {
- DUNE_THROW(Dune::InvalidStateException,
- "The sampling points given to a spline must be sorted by their x value!");
- }
- }
- else {
- y = eval(x);
- dy_dx = evalDerivative(x);
- mono = monotonic(std::max(x_(0), x), std::min(x_(n), x_p1));
- }
-
- std::cout << x << " " << y << " " << dy_dx << " " << mono << "\n";
- }
- }
-
- /*!
- * \brief Evaluate the spline at a given position.
- *
- * \param x The value on the abscissa where the spline ought to be
- * evaluated
- * \param extrapolate If this parameter is set to true, the spline
- * will be extended beyond its range by
- * straight lines, if false calling extrapolate
- * for \f$ x \not [x_{min}, x_{max}]\f$ will
- * cause a failed assertation.
- */
- Scalar eval(Scalar x, bool extrapolate=false) const
- {
- assert(extrapolate || applies(x));
-
- // handle extrapolation
- if (extrapolate) {
- if (x < xMin()) {
- Scalar m = evalDerivative(xMin(), /*segmentIdx=*/0);
- Scalar y0 = y_(0);
- return y0 + m*(x - xMin());
- }
- else if (x > xMax()) {
- Scalar m = evalDerivative(xMax(), /*segmentIdx=*/numSamples_()-2);
- Scalar y0 = y_(numSamples_() - 1);
- return y0 + m*(x - xMax());
- }
- }
-
- return eval_(x, segmentIdx_(x));
- }
-
- /*!
- * \brief Evaluate the spline's derivative at a given position.
- *
- * \param x The value on the abscissa where the spline's
- * derivative ought to be evaluated
- *
- * \param extrapolate If this parameter is set to true, the spline
- * will be extended beyond its range by
- * straight lines, if false calling extrapolate
- * for \f$ x \not [x_{min}, x_{max}]\f$ will
- * cause a failed assertation.
-
- */
- Scalar evalDerivative(Scalar x, bool extrapolate=false) const
- {
- assert(extrapolate || applies(x));
- if (extrapolate) {
- if (x < xMin())
- evalDerivative_(xMin(), 0);
- else if (x > xMax())
- evalDerivative_(xMax(), numSamples_() - 2);
- }
-
- return evalDerivative_(x, segmentIdx_(x));
- }
-
- /*!
- * \brief Find the intersections of the spline with a cubic
- * polynomial in the whole intervall, throws
- * Dune::MathError exception if there is more or less than
- * one solution.
- */
- Scalar intersect(Scalar a, Scalar b, Scalar c, Scalar d) const
- {
- return intersectIntervall(xMin(), xMax(), a, b, c, d);
- }
-
- /*!
- * \brief Find the intersections of the spline with a cubic
- * polynomial in a sub-intervall of the spline, throws
- * Dune::MathError exception if there is more or less than
- * one solution.
- */
- Scalar intersectInterval(Scalar x0, Scalar x1,
- Scalar a, Scalar b, Scalar c, Scalar d) const
- {
- assert(applies(x0) && applies(x1));
-
- Scalar tmpSol[3];
- int nSol = 0;
- int iFirst = segmentIdx_(x0);
- int iLast = segmentIdx_(x1);
- for (int i = iFirst; i <= iLast; ++i)
- {
- nSol += intersectSegment_(tmpSol, i, a, b, c, d, x0, x1);
-
- if (nSol > 1) {
- DUNE_THROW(Dune::MathError,
- "Spline has more than one intersection"); //< x1)
- std::swap(x0, x1);
-
- assert(x0 < x1);
-
- // corner case where the whole spline is a constant
- if (moment_(0) == 0 &&
- moment_(1) == 0 &&
- y_(0) == y_(1))
- {
- // actually the is monotonically increasing as well as
- // monotonously decreasing
- return 2;
- }
-
-
- int i = segmentIdx_(x0);
- if (x_(i) <= x0 && x1 <= x_(i+1)) {
- // the interval to check is completely included in a
- // single segment
- return monotonic_(i, x0, x1);
- }
-
- // make sure that the segments which are completly in the
- // interval [x0, x1] all exhibit the same monotonicity.
- int iEnd = segmentIdx_(x1);
- int r = monotonic_(i, x0, x_(1));
- for (; i < iEnd - 1; ++i)
- if (r != monotonic_(i, x_(i), x_(i + 1)))
- return 0;
-
- // check for the last segment
- if (x_(iEnd) < x1 && r != monotonic_(iEnd, x_(iEnd), x1))
- { return 0; }
-
- return r;
- }
-
- /*!
- * \brief Same as monotonic(x0, x1), but with the entire range of the
- * spline as interval.
- */
- int monotonic() const
- { return monotonic(xMin(), xMax()); }
-
-protected:
- // this is an internal class, so everything is protected!
- SplineCommon_()
- { Valgrind::SetUndefined(asImp_()); }
-
- /*!
- * \brief Set the sampling point vectors.
- *
- * This takes care that the order of the x-values is ascending,
- * although the input must be ordered already!
- */
- template
- void assignSamplingPoints_(DestVector &destX,
- DestVector &destY,
- const SourceVector &srcX,
- const SourceVector &srcY,
- int numSamples)
- {
- assert(numSamples >= 2);
-
- // copy sample points, make sure that the first x value is
- // smaller than the last one
- for (int i = 0; i < numSamples; ++i) {
- int idx = i;
- if (srcX[0] > srcX[numSamples - 1])
- idx = numSamples - i - 1;
- destX[i] = srcX[idx];
- destY[i] = srcY[idx];
- }
- }
-
- template
- void assignFromArrayList_(DestVector &destX,
- DestVector &destY,
- const ListIterator &srcBegin,
- const ListIterator &srcEnd,
- int numSamples)
- {
- assert(numSamples >= 2);
-
- // find out wether the x values are in reverse order
- ListIterator it = srcBegin;
- ++it;
- bool reverse = false;
- if ((*srcBegin)[0] > (*it)[0])
- reverse = true;
- --it;
-
- // loop over all sampling points
- for (int i = 0; it != srcEnd; ++i, ++it) {
- int idx = i;
- if (reverse)
- idx = numSamples - i - 1;
- destX[i] = (*it)[0];
- destY[i] = (*it)[1];
- }
- }
-
- /*!
- * \brief Set the sampling points.
- *
- * Here we assume that the elements of the source vector have an
- * [] operator where v[0] is the x value and v[1] is the y value
- * if the sampling point.
- */
- template
- void assignFromTupleList_(DestVector &destX,
- DestVector &destY,
- ListIterator srcBegin,
- ListIterator srcEnd,
- int numSamples)
- {
- assert(numSamples >= 2);
-
- // copy sample points, make sure that the first x value is
- // smaller than the last one
-
- // find out wether the x values are in reverse order
- ListIterator it = srcBegin;
- ++it;
- bool reverse = false;
- if (std::get<0>(*srcBegin) > std::get<0>(*it))
- reverse = true;
- --it;
-
- // loop over all sampling points
- for (int i = 0; it != srcEnd; ++i, ++it) {
- int idx = i;
- if (reverse)
- idx = numSamples - i - 1;
- destX[i] = std::get<0>(*it);
- destY[i] = std::get<1>(*it);
- }
- }
-
-
- /*!
- * \brief Make the linear system of equations Mx = d which results
- * in the moments of the periodic spline.
- *
- * When solving Mx = d, it should be noted that x[0] is trash and
- * needs to be set to x[n-1]
- */
- template
- void makePeriodicSystem_(Matrix &M, Vector &d)
- {
- // a periodic spline only makes sense if the first and the
- // last sampling point have the same y value
- assert(y_(0) == y_(numSamples_() - 1));
-
- makeNaturalSystem(M, d);
-
- int n = numSamples_() - 1;
-
- Scalar lambda_n = h_(1)/(h_(n) + h_(1));
- Scalar lambda_1 = M[1][2];
- Scalar mu_1 = 1 - lambda_1;
- Scalar mu_n = 1 - lambda_n;
-
- Scalar d_n =
- 6 / (h_(n) + h_(1))
- *
- ( (y_(1) - y_(n))/h_(1) - (y_(n) - y_(n-1))/h_(n));
-
- // insert lambda_n and mu_1
- M[1][n] = mu_1;
- M[n][1] = lambda_n;
- M[n][n-1] = mu_n;
- d[n] = d_n;
- }
-
- /*!
- * \brief Make the linear system of equations Mx = d which results
- * in the moments of the full spline.
- */
- template
- void makeFullSystem_(Matrix &M, Vector &d, Scalar m0, Scalar m1)
- {
- makeNaturalSystem_(M, d);
-
- int n = numSamples_() - 1;
- // first row
- M[0][1] = 1;
- d[0] = 6/h_(1) * ( (y_(1) - y_(0))/h_(1) - m0);
-
- // last row
- M[n][n - 1] = 1;
-
- // right hand side
- d[n] =
- 6/h_(n)
- *
- (m1 - (y_(n) - y_(n - 1))/h_(n));
- }
-
- /*!
- * \brief Make the linear system of equations Mx = d which results
- * in the moments of the natural spline.
- */
- template
- void makeNaturalSystem_(Matrix &M, Vector &d)
- {
- M = 0;
- d = 0;
-
- // See: J. Stoer: "Numerische Mathematik 1", 9th edition,
- // Springer, 2005, p. 111
- const int n = numSamples_() - 1;
-
- // second to next to last rows
- for (int i = 1; i < n; ++i) {
- const Scalar lambda_i = h_(i + 1) / (h_(i) + h_(i+1));
- const Scalar mu_i = 1 - lambda_i;
- const Scalar d_i =
- 6 / (h_(i) + h_(i+1))
- *
- ( (y_(i+1) - y_(i))/h_(i+1) - (y_(i) - y_(i-1))/h_(i));
-
- M[i][i-1] = mu_i;
- M[i][i] = 2;
- M[i][i + 1] = lambda_i;
- d[i] = d_i;
- };
-
- // first row
- M[0][0] = 2;
-
- // last row
- M[n][n] = 2;
- // right hand side
- d[0] = 0.0;
- d[n] = 0.0;
- }
-
- // evaluate the spline at a given the position and given the
- // segment index
- Scalar eval_(Scalar x, int i) const
- {
- // See: J. Stoer: "Numerische Mathematik 1", 9th edition,
- // Springer, 2005, p. 109
- Scalar h_i1 = h_(i + 1);
- Scalar x_i = x - x_(i);
- Scalar x_i1 = x_(i+1) - x;
-
- Scalar A_i =
- (y_(i+1) - y_(i))/h_i1
- -
- h_i1/6*(moment_(i+1) - moment_(i));
- Scalar B_i = y_(i) - moment_(i)* (h_i1*h_i1) / 6;
-
- return
- moment_(i)* x_i1*x_i1*x_i1 / (6 * h_i1)
- +
- moment_(i + 1)* x_i*x_i*x_i / (6 * h_i1)
- +
- A_i*x_i
- +
- B_i;
- }
-
- // evaluate the derivative of a spline given the actual position
- // and the segment index
- Scalar evalDerivative_(Scalar x, int i) const
- {
- // See: J. Stoer: "Numerische Mathematik 1", 9th edition,
- // Springer, 2005, p. 109
- Scalar h_i1 = h_(i + 1);
- Scalar x_i = x - x_(i);
- Scalar x_i1 = x_(i+1) - x;
-
- Scalar A_i =
- (y_(i+1) - y_(i))/h_i1
- -
- h_i1/6*(moment_(i+1) - moment_(i));
-
- return
- -moment_(i) * x_i1*x_i1 / (2 * h_i1)
- +
- moment_(i + 1) * x_i*x_i / (2 * h_i1)
- +
- A_i;
- }
-
-
- // returns the monotonicality of an interval of a spline segment
- //
- // The return value have the following meaning:
- //
- // 1: spline is monotonously increasing in the specified interval
- // 0: spline is not monotonic in the specified interval
- // -1: spline is monotonously decreasing in the specified interval
- int monotonic_(int i, Scalar x0, Scalar x1) const
- {
- // shift the interval so that it is consistent with the
- // definitions by Stoer
- x0 = x0 - x_(i);
- x1 = x1 - x_(i);
-
- Scalar a = a_(i);
- Scalar b = b_(i);
- Scalar c = c_(i);
-
- Scalar disc = 4*b*b - 12*a*c;
- if (disc < 0) {
- // discriminant is smaller than 0, i.e. the segment does
- // not exhibit any extrema.
- return (x0*(x0*3*a + 2*b) + c > 0) ? 1 : -1;
- }
- disc = std::sqrt(disc);
- Scalar xE1 = (-2*b + disc)/(6*a);
- Scalar xE2 = (-2*b - disc)/(6*a);
-
- if (disc == 0) {
- // saddle point -> no extrema
- if (xE1 == x0)
- // make sure that we're not picking the saddle point
- // to determine whether we're monotonically increasing
- // or decreasing
- x0 = x1;
- return (x0*(x0*3*a + 2*b) + c > 0) ? 1 : -1;
- };
- if ((x0 < xE1 && xE1 < x1) ||
- (x0 < xE2 && xE2 < x1))
- {
- // there is an extremum in the range (x0, x1)
- return 0;
- }
- // no extremum in range (x0, x1)
- x0 = (x0 + x1)/2; // pick point in the middle of the interval
- // to avoid extrema on the boundaries
- return (x0*(x0*3*a + 2*b) + c > 0) ? 1 : -1;
- }
-
- /*!
- * \brief Find all the intersections of a segment of the spline
- * with a cubic polynomial within a specified interval.
- */
- int intersectSegment_(Scalar *sol,
- int segIdx,
- Scalar a, Scalar b, Scalar c, Scalar d,
- Scalar x0 = -1e100, Scalar x1 = 1e100) const
- {
- int n = Opm::invertCubicPolynomial(sol,
- a_(segIdx) - a,
- b_(segIdx) - b,
- c_(segIdx) - c,
- d_(segIdx) - d);
- x0 = std::max(x_(segIdx), x0);
- x1 = std::max(x_(segIdx+1), x1);
-
- // filter the intersections outside of the specified intervall
- int k = 0;
- for (int j = 0; j < n; ++j) {
- sol[j] += x_(segIdx); // add the offset of the intervall. For details see Stoer
- if (x0 <= sol[j] && sol[j] <= x1) {
- sol[k] = sol[j];
- ++k;
- }
- }
- return k;
- }
-
- // find the segment index for a given x coordinate
- int segmentIdx_(Scalar x) const
- {
- // bisection
- int iLow = 0;
- int iHigh = numSamples_() - 1;
-
- while (iLow + 1 < iHigh) {
- int i = (iLow + iHigh) / 2;
- if (x_(i) > x)
- iHigh = i;
- else
- iLow = i;
- };
- return iLow;
- }
-
- /*!
- * \brief Returns x[i] - x[i - 1]
- */
- Scalar h_(int i) const
- {
- assert(x_(i) > x_(i-1)); // the sampling points must be given
- // in ascending order
- return x_(i) - x_(i - 1);
- }
-
- /*!
- * \brief Returns the y coordinate of the i-th sampling point.
- */
- Scalar x_(int i) const
- { return asImp_().x_(i); }
-
- /*!
- * \brief Returns the y coordinate of the i-th sampling point.
- */
- Scalar y_(int i) const
- { return asImp_().y_(i); }
-
- /*!
- * \brief Returns the moment (i.e. second derivative) of the
- * spline at the i-th sampling point.
- */
- Scalar moment_(int i) const
- { return asImp_().moment_(i); }
-
- // returns the coefficient in front of the x^0 term. In Stoer this
- // is delta.
- Scalar a_(int i) const
- { return (moment_(i+1) - moment_(i))/(6*h_(i+1)); }
-
- // returns the coefficient in front of the x^2 term In Stoer this
- // is gamma.
- Scalar b_(int i) const
- { return moment_(i)/2; }
-
- // returns the coefficient in front of the x^1 term. In Stoer this
- // is beta.
- Scalar c_(int i) const
- { return (y_(i+1) - y_(i))/h_(i + 1) - h_(i+1)/6*(2*moment_(i) + moment_(i+1)); }
-
- // returns the coefficient in front of the x^0 term. In Stoer this
- // is alpha.
- Scalar d_(int i) const
- { return y_(i); }
-
- /*!
- * \brief Returns the number of sampling points.
- */
- int numSamples_() const
- { return asImp_().numSamples(); }
-};
-
-}
-
-#endif
diff --git a/opm/common/variablelengthspline_.hh b/opm/common/variablelengthspline_.hh
deleted file mode 100644
index 6e743e253..000000000
--- a/opm/common/variablelengthspline_.hh
+++ /dev/null
@@ -1,479 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * Copyright (C) 2010-2012 by Andreas Lauser *
- * *
- * 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, see . *
- *****************************************************************************/
-/*!
- * \file
- * \copydoc Opm::VariableLengthSpline_
- */
-#ifndef OPM_VARIABLE_LENGTH_SPLINE_HH
-#define OPM_VARIABLE_LENGTH_SPLINE_HH
-
-#include
-#include
-#include
-#include
-
-#include "splinecommon_.hh"
-
-namespace Opm
-{
-/*!
- * \brief The common code for all 3rd order polynomial splines with
- * where the number of sampling points only known at run-time.
- */
-template
-class VariableLengthSpline_
- : public SplineCommon_ >
-{
- friend class SplineCommon_ >;
-
- typedef ScalarT Scalar;
- typedef Dune::BlockVector > Vector;
- typedef Dune::BTDMatrix > BTDMatrix;
-
-public:
- /*!
- * \brief Returns the number of sampling points.
- */
- int numSamples() const
- { return xPos_.size(); }
-
-
- ///////////////////////////////////////
- ///////////////////////////////////////
- ///////////////////////////////////////
- // Full splines //
- ///////////////////////////////////////
- ///////////////////////////////////////
- ///////////////////////////////////////
-
- /*!
- * \brief Set the sampling points and the boundary slopes of a
- * full spline using C-style arrays.
- *
- * This method uses separate array-like objects for the values of
- * the X and Y coordinates. In this context 'array-like' means
- * that an access to the members is provided via the []
- * operator. (e.g. C arrays, std::vector, std::array, etc.) Each
- * array must be of size 'nSamples' at least. Also, the number of
- * sampling points must be larger than 1.
- */
- template
- void setXYArrays(int nSamples,
- const ScalarArrayX &x,
- const ScalarArrayY &y,
- Scalar m0, Scalar m1)
- {
- assert(nSamples > 1);
-
- setNumSamples_(nSamples);
- for (int i = 0; i < nSamples; ++i) {
- xPos_[i] = x[i];
- yPos_[i] = y[i];
- }
-
- if (xPos_[0] > xPos_[numSamples() - 1])
- reverseSamplingPoints_();
-
- makeFullSpline_(m0, m1);
- }
-
- /*!
- * \brief Set the sampling points and the boundary slopes of a
- * full spline using STL-compatible containers.
- *
- * This method uses separate STL-compatible containers for the
- * values of the sampling points' X and Y
- * coordinates. "STL-compatible" means that the container provides
- * access to iterators using the begin(), end() methods and also
- * provides a size() method. Also, the number of entries in the X
- * and the Y containers must be equal and larger than 1.
- */
- template
- void setXYContainers(const ScalarContainerX &x,
- const ScalarContainerY &y,
- Scalar m0, Scalar m1)
- {
- assert(x.size() == y.size());
- assert(x.size() > 1);
-
- setNumSamples_(x.size());
-
- std::copy(x.begin(), x.end(), xPos_.begin());
- std::copy(y.begin(), y.end(), yPos_.begin());
-
- if (xPos_[0] > xPos_[numSamples() - 1])
- reverseSamplingPoints_();
-
- makeFullSpline_(m0, m1);
- }
-
- /*!
- * \brief Set the sampling points and the boundary slopes of a
- * full spline using a C-style array.
- *
- * This method uses a single array of sampling points, which are
- * seen as an array-like object which provides access to the X and
- * Y coordinates. In this context 'array-like' means that an
- * access to the members is provided via the [] operator. (e.g. C
- * arrays, std::vector, std::array, etc.) The array containing
- * the sampling points must be of size 'nSamples' at least. Also,
- * the number of sampling points must be larger than 1.
- */
- template
- void setArrayOfPoints(int nSamples,
- const PointArray &points,
- Scalar m0,
- Scalar m1)
- {
- // a spline with no or just one sampling points? what an
- // incredible bad idea!
- assert(nSamples > 1);
-
- setNumSamples_(nSamples);
- for (int i = 0; i < nSamples; ++i) {
- xPos_[i] = points[i][0];
- yPos_[i] = points[i][1];
- }
-
- if (xPos_[0] > xPos_[numSamples() - 1])
- reverseSamplingPoints_();
-
- makeFullSpline_(m0, m1);
- }
-
- /*!
- * \brief Set the sampling points and the boundary slopes of a
- * full spline using a STL-compatible container of
- * array-like objects.
- *
- * This method uses a single STL-compatible container of sampling
- * points, which are assumed to be array-like objects storing the
- * X and Y coordinates. "STL-compatible" means that the container
- * provides access to iterators using the begin(), end() methods
- * and also provides a size() method. Also, the number of entries
- * in the X and the Y containers must be equal and larger than 1.
- */
- template
- void setContainerOfPoints(const XYContainer &points,
- Scalar m0,
- Scalar m1)
- {
- // a spline with no or just one sampling points? what an
- // incredible bad idea!
- assert(points.size() > 1);
-
- setNumSamples_(points.size());
- typename XYContainer::const_iterator it = points.begin();
- typename XYContainer::const_iterator endIt = points.end();
- for (int i = 0; it != endIt; ++i, ++it) {
- xPos_[i] = (*it)[0];
- yPos_[i] = (*it)[1];
- }
-
- if (xPos_[0] > xPos_[numSamples() - 1])
- reverseSamplingPoints_();
-
- // make a full spline
- makeFullSpline_(m0, m1);
- }
-
- /*!
- * \brief Set the sampling points and the boundary slopes of a
- * full spline using a STL-compatible container of
- * tuple-like objects.
- *
- * This method uses a single STL-compatible container of sampling
- * points, which are assumed to be tuple-like objects storing the
- * X and Y coordinates. "tuple-like" means that the objects
- * provide access to the x values via std::get<0>(obj) and to the
- * y value via std::get<1>(obj) (e.g. std::tuple or
- * std::pair). "STL-compatible" means that the container provides
- * access to iterators using the begin(), end() methods and also
- * provides a size() method. Also, the number of entries in the X
- * and the Y containers must be equal and larger than 1.
- */
- template
- void setContainerOfTuples(const XYContainer &points,
- Scalar m0,
- Scalar m1)
- {
- // resize internal arrays
- setNumSamples_(points.size());
- typename XYContainer::const_iterator it = points.begin();
- typename XYContainer::const_iterator endIt = points.end();
- for (int i = 0; it != endIt; ++i, ++it) {
- xPos_[i] = std::get<0>(*it);
- yPos_[i] = std::get<1>(*it);
- }
-
- if (xPos_[0] > xPos_[numSamples() - 1])
- reverseSamplingPoints_();
-
- // make a full spline
- makeFullSpline_(m0, m1);
- }
-
- ///////////////////////////////////////
- ///////////////////////////////////////
- ///////////////////////////////////////
- // Natural splines //
- ///////////////////////////////////////
- ///////////////////////////////////////
- ///////////////////////////////////////
- /*!
- * \brief Set the sampling points natural spline using C-style arrays.
- *
- * This method uses separate array-like objects for the values of
- * the X and Y coordinates. In this context 'array-like' means
- * that an access to the members is provided via the []
- * operator. (e.g. C arrays, std::vector, std::array, etc.) Each
- * array must be of size 'nSamples' at least. Also, the number of
- * sampling points must be larger than 1.
- */
- template
- void setXYArrays(int nSamples,
- const ScalarArrayX &x,
- const ScalarArrayY &y)
- {
- assert(nSamples > 1);
-
- setNumSamples_(nSamples);
- for (int i = 0; i < nSamples; ++i) {
- xPos_[i] = x[i];
- yPos_[i] = y[i];
- }
-
- if (xPos_[0] > xPos_[numSamples() - 1])
- reverseSamplingPoints_();
-
- makeNaturalSpline_();
- }
-
- /*!
- * \brief Set the sampling points of a natural spline using
- * STL-compatible containers.
- *
- * This method uses separate STL-compatible containers for the
- * values of the sampling points' X and Y
- * coordinates. "STL-compatible" means that the container provides
- * access to iterators using the begin(), end() methods and also
- * provides a size() method. Also, the number of entries in the X
- * and the Y containers must be equal and larger than 1.
- */
- template
- void setXYContainers(const ScalarContainerX &x,
- const ScalarContainerY &y)
- {
- assert(x.size() == y.size());
- assert(x.size() > 1);
-
- setNumSamples_(x.size());
- std::copy(x.begin(), x.end(), xPos_.begin());
- std::copy(y.begin(), y.end(), yPos_.begin());
-
- if (xPos_[0] > xPos_[numSamples() - 1])
- reverseSamplingPoints_();
-
- makeNaturalSpline_();
- }
-
- /*!
- * \brief Set the sampling points of a natural spline using a
- * C-style array.
- *
- * This method uses a single array of sampling points, which are
- * seen as an array-like object which provides access to the X and
- * Y coordinates. In this context 'array-like' means that an
- * access to the members is provided via the [] operator. (e.g. C
- * arrays, std::vector, std::array, etc.) The array containing
- * the sampling points must be of size 'nSamples' at least. Also,
- * the number of sampling points must be larger than 1.
- */
- template
- void setArrayOfPoints(int nSamples,
- const PointArray &points)
- {
- // a spline with no or just one sampling points? what an
- // incredible bad idea!
- assert(nSamples > 1);
-
- setNumSamples_(nSamples);
- for (int i = 0; i < nSamples; ++i) {
- xPos_[i] = points[i][0];
- yPos_[i] = points[i][1];
- }
-
- if (xPos_[0] > xPos_[numSamples() - 1])
- reverseSamplingPoints_();
-
- makeNaturalSpline_();
- }
-
- /*!
- * \brief Set the sampling points of a natural spline using a
- * STL-compatible container of array-like objects.
- *
- * This method uses a single STL-compatible container of sampling
- * points, which are assumed to be array-like objects storing the
- * X and Y coordinates. "STL-compatible" means that the container
- * provides access to iterators using the begin(), end() methods
- * and also provides a size() method. Also, the number of entries
- * in the X and the Y containers must be equal and larger than 1.
- */
- template
- void setContainerOfPoints(const XYContainer &points)
- {
- // a spline with no or just one sampling points? what an
- // incredible bad idea!
- assert(points.size() > 1);
-
- setNumSamples_(points.size());
- typename XYContainer::const_iterator it = points.begin();
- typename XYContainer::const_iterator endIt = points.end();
- for (int i = 0; it != endIt; ++ i, ++it) {
- xPos_[i] = (*it)[0];
- yPos_[i] = (*it)[1];
- }
-
- if (xPos_[0] > xPos_[numSamples() - 1])
- reverseSamplingPoints_();
-
- // make a natural spline
- makeNaturalSpline_();
- }
-
- /*!
- * \brief Set the sampling points of a natural spline using a
- * STL-compatible container of tuple-like objects.
- *
- * This method uses a single STL-compatible container of sampling
- * points, which are assumed to be tuple-like objects storing the
- * X and Y coordinates. "tuple-like" means that the objects
- * provide access to the x values via std::get<0>(obj) and to the
- * y value via std::get<1>(obj) (e.g. std::tuple or
- * std::pair). "STL-compatible" means that the container provides
- * access to iterators using the begin(), end() methods and also
- * provides a size() method. Also, the number of entries in the X
- * and the Y containers must be equal and larger than 1.
- */
- template
- void setContainerOfTuples(const XYContainer &points)
- {
- // resize internal arrays
- setNumSamples_(points.size());
- typename XYContainer::const_iterator it = points.begin();
- typename XYContainer::const_iterator endIt = points.end();
- for (int i = 0; it != endIt; ++i, ++it) {
- xPos_[i] = std::get<0>(*it);
- yPos_[i] = std::get<1>(*it);
- }
-
- if (xPos_[0] > xPos_[numSamples() - 1])
- reverseSamplingPoints_();
-
- // make a natural spline
- makeNaturalSpline_();
- }
-
-protected:
- /*!
- * \brief Reverse order of the elements in the arrays which
- * contain the sampling points.
- */
- void reverseSamplingPoints_()
- {
- // reverse the arrays
- int n = numSamples();
- for (int i = 0; i <= (n - 1)/2; ++i) {
- std::swap(xPos_[i], xPos_[n - i - 1]);
- std::swap(yPos_[i], yPos_[n - i - 1]);
- }
- }
-
- /*!
- * \brief Resizes the internal vectors to store the sample points.
- */
- void setNumSamples_(int nSamples)
- {
- xPos_.resize(nSamples);
- yPos_.resize(nSamples);
- m_.resize(nSamples);
- }
-
- /*!
- * \brief Create a natural spline from the already set sampling points.
- *
- * Also creates temporary matrix and right hand side vector.
- */
- void makeFullSpline_(Scalar m0, Scalar m1)
- {
- BTDMatrix M(numSamples());
- Vector d(numSamples());
-
- // create linear system of equations
- this->makeFullSystem_(M, d, m0, m1);
-
- // solve for the moments
- M.solve(m_, d);
- }
-
- /*!
- * \brief Create a natural spline from the already set sampling points.
- *
- * Also creates temporary matrix and right hand side vector.
- */
- void makeNaturalSpline_()
- {
- BTDMatrix M(numSamples());
- Vector d(numSamples());
-
- // create linear system of equations
- this->makeNaturalSystem_(M, d);
-
- // solve for the moments
- M.solve(m_, d);
- }
-
- /*!
- * \brief Returns the x coordinate of the i-th sampling point.
- */
- Scalar x_(int i) const
- { return xPos_[i]; }
-
- /*!
- * \brief Returns the y coordinate of the i-th sampling point.
- */
- Scalar y_(int i) const
- { return yPos_[i]; }
-
- /*!
- * \brief Returns the moment (i.e. second derivative) of the
- * spline at the i-th sampling point.
- */
- Scalar moment_(int i) const
- { return m_[i]; }
-
- Vector xPos_;
- Vector yPos_;
- Vector m_;
-};
-}
-
-#endif
diff --git a/opm/material/binarycoefficients/air_mesitylene.hh b/opm/material/binarycoefficients/air_mesitylene.hh
index db9c166a9..73df039a2 100644
--- a/opm/material/binarycoefficients/air_mesitylene.hh
+++ b/opm/material/binarycoefficients/air_mesitylene.hh
@@ -43,9 +43,7 @@ public:
*/
template
static Scalar henry(Scalar temperature)
- { DUNE_THROW(Dune::NotImplemented,
- "Henry coefficient of air in mesitylene");
- }
+ { OPM_THROW(std::runtime_error, "Not implemented: Henry coefficient of air in mesitylene"); }
/*!
* \brief Binary diffusion coefficent [m^2/s] for air and mesitylene.
@@ -96,8 +94,7 @@ public:
*/
template
static Scalar liquidDiffCoeff(Scalar temperature, Scalar pressure)
- { DUNE_THROW(Dune::NotImplemented,
- "Binary liquid diffusion coefficients of air and mesitylene");
+ { OPM_THROW(std::runtime_error, "Not implemented: Binary liquid diffusion coefficients of air and mesitylene");
}
};
diff --git a/opm/material/binarycoefficients/air_xylene.hh b/opm/material/binarycoefficients/air_xylene.hh
index 82a19c7e4..f1e8ecce6 100644
--- a/opm/material/binarycoefficients/air_xylene.hh
+++ b/opm/material/binarycoefficients/air_xylene.hh
@@ -43,9 +43,7 @@ public:
*/
template
static Scalar henry(Scalar temperature)
- { DUNE_THROW(Dune::NotImplemented,
- "Henry coefficient of air in xylene");
- }
+ { OPM_THROW(std::runtime_error, "Not implemented: Henry coefficient of air in xylene"); }
/*!
* \brief Binary diffusion coefficent [m^2/s] for air and xylene.
@@ -96,8 +94,7 @@ public:
*/
template
static Scalar liquidDiffCoeff(Scalar temperature, Scalar pressure)
- { DUNE_THROW(Dune::NotImplemented,
- "Binary liquid diffusion coefficients of air and xylene");
+ { OPM_THROW(std::runtime_error, "Not implemented: Binary liquid diffusion coefficients of air and xylene");
}
};
diff --git a/opm/material/binarycoefficients/fullermethod.hh b/opm/material/binarycoefficients/fullermethod.hh
index 958fb1b6f..3086c7df0 100644
--- a/opm/material/binarycoefficients/fullermethod.hh
+++ b/opm/material/binarycoefficients/fullermethod.hh
@@ -23,12 +23,13 @@
#ifndef OPM_FULLERMETHOD_HH
#define OPM_FULLERMETHOD_HH
-#include
+#include
+
+#include
+
+namespace Opm {
+namespace BinaryCoeff {
-namespace Opm
-{
-namespace BinaryCoeff
-{
/*!
* \ingroup Binarycoefficients
* \brief Estimate binary diffusion coefficents \f$\mathrm{[m^2/s]}\f$ in gases according to
@@ -53,7 +54,7 @@ inline Scalar fullerMethod(const Scalar *M, // molar masses [g/mol]
const Scalar pressure) // [Pa]
{
// "effective" molar mass in [g/m^3]
- Scalar Mab = harmonicMean(M[0], M[1]);
+ Scalar Mab = Opm::utils::harmonicAverage(M[0], M[1]);
// Fuller's method
Scalar tmp = std::pow(SigmaNu[0], 1./3) + std::pow(SigmaNu[1], 1./3);
diff --git a/opm/material/binarycoefficients/h2o_co2.hh b/opm/material/binarycoefficients/h2o_co2.hh
index fc0bde485..8c8113bd9 100644
--- a/opm/material/binarycoefficients/h2o_co2.hh
+++ b/opm/material/binarycoefficients/h2o_co2.hh
@@ -87,10 +87,7 @@ public:
*/
template
static Scalar liquidDiffCoeff(Scalar temperature, Scalar pressure)
- {
- DUNE_THROW(Dune::NotImplemented,
- "Binary liquid diffusion coefficients of CO2 and CH4");
- }
+ { OPM_THROW(std::runtime_error, "Not implemented: Binary liquid diffusion coefficients of CO2 and CH4"); }
};
}
diff --git a/opm/material/components/air.hh b/opm/material/components/air.hh
index 4705fa0e9..6c6333adb 100644
--- a/opm/material/components/air.hh
+++ b/opm/material/components/air.hh
@@ -24,10 +24,12 @@
#ifndef OPM_AIR_HH
#define OPM_AIR_HH
-#include
#include
#include
+#include
+#include
+
namespace Opm {
/*!
@@ -156,8 +158,8 @@ public:
Scalar r;
if(temperature < 273.15 || temperature > 660.)
{
- DUNE_THROW(Dune::NotImplemented,
- "ConstrelAir: Temperature out of range at ");
+ OPM_THROW(MaterialLawProblem,
+ "ConstrelAir: Temperature out of range at ");
}
r = 1.496*1.E-6*std::pow(temperature,1.5)/(temperature+120.);
return (r);
diff --git a/opm/material/components/co2.hh b/opm/material/components/co2.hh
index 9fb4cbbad..d1a9ea849 100644
--- a/opm/material/components/co2.hh
+++ b/opm/material/components/co2.hh
@@ -24,7 +24,8 @@
#ifndef OPM_CO2_HH
#define OPM_CO2_HH
-#include
+#include
+#include
#include
#include
#include
@@ -163,8 +164,8 @@ public:
#ifndef NDEBUG
if ((temperature < criticalTemperature() or pressure < criticalPressure()) and !warningPrinted)
{
- Dune::dwarn << "The tables used for the CO2 exhibit subcritical as well as critical values: "
- << "Double-check your results if you cross the critical-subcritical line during the simulation!\n";
+ std::cout << "The tables used for the CO2 exhibit subcritical as well as critical values: "
+ << "Double-check your results if you cross the critical-subcritical line during the simulation!\n";
warningPrinted=true;
}
#endif
@@ -192,8 +193,8 @@ public:
#ifndef NDEBUG
if ((temperature < criticalTemperature() or pressure < criticalPressure()) and !warningPrinted)
{
- Dune::dwarn << "The tables used for the CO2 exhibit subcritical as well as critical values: "
- << "Double-check your results if you cross the critical-subcritical line during the simulation!\n";
+ std::cout << "The tables used for the CO2 exhibit subcritical as well as critical values: "
+ << "Double-check your results if you cross the critical-subcritical line during the simulation!\n";
warningPrinted=true;
}
#endif
diff --git a/opm/material/components/component.hh b/opm/material/components/component.hh
index cc64d274b..d1578e44e 100644
--- a/opm/material/components/component.hh
+++ b/opm/material/components/component.hh
@@ -25,10 +25,10 @@
#ifndef OPM_COMPONENT_HH
#define OPM_COMPONENT_HH
-#include
+#include
+#include
-namespace Opm
-{
+namespace Opm {
/*!
* \ingroup Components
@@ -57,61 +57,61 @@ public:
*/
static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp,
Scalar pressMin, Scalar pressMax, unsigned nPress)
- { Dune::dwarn << "No init routine defined - make sure that this is not necessary!" << std::endl; }
+ { }
/*!
* \brief Returns true iff the gas phase is assumed to be compressible
*/
static bool gasIsCompressible()
- { DUNE_THROW(Dune::NotImplemented, "Component::gasIsCompressible()"); }
+ { OPM_THROW(std::runtime_error, "Not implemented: Component::gasIsCompressible()"); }
/*!
* \brief Returns true iff the gas phase is assumed to be ideal
*/
static bool gasIsIdeal()
- { DUNE_THROW(Dune::NotImplemented, "Component::gasIsIdeal()"); }
+ { OPM_THROW(std::runtime_error, "Not implemented: Component::gasIsIdeal()"); }
/*!
* \brief Returns true iff the liquid phase is assumed to be compressible
*/
static bool liquidIsCompressible()
- { DUNE_THROW(Dune::NotImplemented, "Component::liquidIsCompressible()"); }
+ { OPM_THROW(std::runtime_error, "Not implemented: Component::liquidIsCompressible()"); }
/*!
* \brief A human readable name for the component.
*/
static const char *name()
- { DUNE_THROW(Dune::NotImplemented, "Component::name()"); }
+ { OPM_THROW(std::runtime_error, "Not implemented: Component::name()"); }
/*!
* \brief The molar mass in \f$\mathrm{[kg]}\f$ of the component.
*/
static Scalar molarMass()
- { DUNE_THROW(Dune::NotImplemented, "Component::molarMass()"); }
+ { OPM_THROW(std::runtime_error, "Not implemented: Component::molarMass()"); }
/*!
* \brief Returns the critical temperature in \f$\mathrm{[K]}\f$ of the component.
*/
static Scalar criticalTemperature()
- { DUNE_THROW(Dune::NotImplemented, "Component::criticalTemperature()"); }
+ { OPM_THROW(std::runtime_error, "Not implemented: Component::criticalTemperature()"); }
/*!
* \brief Returns the critical pressure in \f$\mathrm{[Pa]}\f$ of the component.
*/
static Scalar criticalPressure()
- { DUNE_THROW(Dune::NotImplemented, "Component::criticalPressure()"); }
+ { OPM_THROW(std::runtime_error, "Not implemented: Component::criticalPressure()"); }
/*!
* \brief Returns the temperature in \f$\mathrm{[K]}\f$ at the component's triple point.
*/
static Scalar tripleTemperature()
- { DUNE_THROW(Dune::NotImplemented, "Component::tripleTemperature()"); }
+ { OPM_THROW(std::runtime_error, "Not implemented: Component::tripleTemperature()"); }
/*!
* \brief Returns the pressure in \f$\mathrm{[Pa]}\f$ at the component's triple point.
*/
static Scalar triplePressure()
- { DUNE_THROW(Dune::NotImplemented, "Component::triplePressure()"); }
+ { OPM_THROW(std::runtime_error, "Not implemented: Component::triplePressure()"); }
/*!
* \brief The vapor pressure in \f$\mathrm{[Pa]}\f$ of the component at a given
@@ -120,7 +120,7 @@ public:
* \param T temperature of the component in \f$\mathrm{[K]}\f$
*/
static Scalar vaporPressure(Scalar T)
- { DUNE_THROW(Dune::NotImplemented, "Component::vaporPressure()"); }
+ { OPM_THROW(std::runtime_error, "Not implemented: Component::vaporPressure()"); }
/*!
* \brief The density in \f$\mathrm{[kg/m^3]}\f$ of the component at a given pressure in \f$\mathrm{[Pa]}\f$ and temperature in \f$\mathrm{[K]}\f$.
@@ -129,7 +129,7 @@ public:
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar gasDensity(Scalar temperature, Scalar pressure)
- { DUNE_THROW(Dune::NotImplemented, "Component::gasDensity()"); }
+ { OPM_THROW(std::runtime_error, "Not implemented: Component::gasDensity()"); }
/*!
* \brief The density \f$\mathrm{[kg/m^3]}\f$ of the liquid component at a given pressure in \f$\mathrm{[Pa]}\f$ and temperature in \f$\mathrm{[K]}\f$.
@@ -138,7 +138,7 @@ public:
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar liquidDensity(Scalar temperature, Scalar pressure)
- { DUNE_THROW(Dune::NotImplemented, "Component::liquidDensity()"); }
+ { OPM_THROW(std::runtime_error, "Not implemented: Component::liquidDensity()"); }
/*!
* \brief Specific enthalpy \f$\mathrm{[J/kg]}\f$ of the pure component in gas.
@@ -147,7 +147,7 @@ public:
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static const Scalar gasEnthalpy(Scalar temperature, Scalar pressure)
- { DUNE_THROW(Dune::NotImplemented, "Component::gasEnthalpy()"); }
+ { OPM_THROW(std::runtime_error, "Not implemented: Component::gasEnthalpy()"); }
/*!
* \brief Specific enthalpy \f$\mathrm{[J/kg]}\f$ of the pure component in liquid.
@@ -156,7 +156,7 @@ public:
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static const Scalar liquidEnthalpy(Scalar temperature, Scalar pressure)
- { DUNE_THROW(Dune::NotImplemented, "Component::liquidEnthalpy()"); }
+ { OPM_THROW(std::runtime_error, "Not implemented: Component::liquidEnthalpy()"); }
/*!
* \brief Specific internal energy \f$\mathrm{[J/kg]}\f$ of the pure component in gas.
@@ -165,7 +165,7 @@ public:
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static const Scalar gasInternalEnergy(Scalar temperature, Scalar pressure)
- { DUNE_THROW(Dune::NotImplemented, "Component::gasInternalEnergy()"); }
+ { OPM_THROW(std::runtime_error, "Not implemented: Component::gasInternalEnergy()"); }
/*!
* \brief Specific internal energy \f$\mathrm{[J/kg]}\f$ of pure the pure component in liquid.
@@ -174,7 +174,7 @@ public:
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static const Scalar liquidInternalEnergy(Scalar temperature, Scalar pressure)
- { DUNE_THROW(Dune::NotImplemented, "Component::liquidInternalEnergy()"); }
+ { OPM_THROW(std::runtime_error, "Not implemented: Component::liquidInternalEnergy()"); }
/*!
* \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of the pure component at a given pressure in \f$\mathrm{[Pa]}\f$ and
@@ -184,7 +184,7 @@ public:
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar gasViscosity(Scalar temperature, Scalar pressure)
- { DUNE_THROW(Dune::NotImplemented, "Component::gasViscosity()"); }
+ { OPM_THROW(std::runtime_error, "Not implemented: Component::gasViscosity()"); }
/*!
* \brief The dynamic liquid viscosity \f$\mathrm{[Pa*s]}\f$ of the pure component.
@@ -193,31 +193,31 @@ public:
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
static Scalar liquidViscosity(Scalar temperature, Scalar pressure)
- { DUNE_THROW(Dune::NotImplemented, "Component::liquidViscosity()"); }
+ { OPM_THROW(std::runtime_error, "Not implemented: Component::liquidViscosity()"); }
/*!
* \brief Thermal conductivity of the component [W/(m^2 K/m)] as a gas.
*/
static Scalar gasThermalConductivity(Scalar temperature, Scalar pressure)
- { DUNE_THROW(Dune::NotImplemented, "Component::gasThermalConductivity()"); }
+ { OPM_THROW(std::runtime_error, "Not implemented: Component::gasThermalConductivity()"); }
/*!
* \brief Thermal conductivity of the component [W/(m^2 K/m)] as a liquid.
*/
static Scalar liquidThermalConductivity(Scalar temperature, Scalar pressure)
- { DUNE_THROW(Dune::NotImplemented, "Component::liquidThermalConductivity()"); }
+ { OPM_THROW(std::runtime_error, "Not implemented: Component::liquidThermalConductivity()"); }
/*!
* \brief Specific isobaric heat capacity of the component [J/kg] as a gas.
*/
static Scalar gasHeatCapacity(Scalar temperature, Scalar pressure)
- { DUNE_THROW(Dune::NotImplemented, "Component::gasHeatCapacity()"); }
+ { OPM_THROW(std::runtime_error, "Not implemented: Component::gasHeatCapacity()"); }
/*!
* \brief Specific isobaric heat capacity of the component [J/kg] as a liquid.
*/
static Scalar liquidHeatCapacity(Scalar temperature, Scalar pressure)
- { DUNE_THROW(Dune::NotImplemented, "Component::liquidHeatCapacity()"); }
+ { OPM_THROW(std::runtime_error, "Not implemented: Component::liquidHeatCapacity()"); }
};
} // end namepace
diff --git a/opm/material/components/h2o.hh b/opm/material/components/h2o.hh
index 71a58c4fe..03b58b426 100644
--- a/opm/material/components/h2o.hh
+++ b/opm/material/components/h2o.hh
@@ -28,8 +28,9 @@
#include
#include
-#include
-#include
+#include
+#include
+#include
#include "component.hh"
@@ -172,7 +173,7 @@ public:
{
if (!Region2::isValid(temperature, pressure))
{
- DUNE_THROW(NumericalProblem,
+ OPM_THROW(NumericalProblem,
"Enthalpy of steam is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
}
@@ -222,7 +223,7 @@ public:
{
if (!Region1::isValid(temperature, pressure))
{
- DUNE_THROW(NumericalProblem,
+ OPM_THROW(NumericalProblem,
"Enthalpy of water is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
}
@@ -263,7 +264,7 @@ public:
{
if (!Region2::isValid(temperature, pressure))
{
- DUNE_THROW(NumericalProblem,
+ OPM_THROW(NumericalProblem,
"Heat capacity of steam is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
}
@@ -299,7 +300,7 @@ public:
{
if (!Region1::isValid(temperature, pressure))
{
- DUNE_THROW(NumericalProblem,
+ OPM_THROW(NumericalProblem,
"heat Capacity of water is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
}
@@ -332,7 +333,7 @@ public:
{
if (!Region1::isValid(temperature, pressure))
{
- DUNE_THROW(NumericalProblem,
+ OPM_THROW(NumericalProblem,
"Internal Energy of water is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
}
@@ -388,7 +389,7 @@ public:
{
if (!Region2::isValid(temperature, pressure))
{
- DUNE_THROW(NumericalProblem,
+ OPM_THROW(NumericalProblem,
"Internal Energy of steam is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
}
@@ -462,7 +463,7 @@ public:
{
if (!Region1::isValid(temperature, pressure))
{
- DUNE_THROW(NumericalProblem,
+ OPM_THROW(NumericalProblem,
"Heat capacity of water is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
}
@@ -495,7 +496,7 @@ public:
{
if (!Region2::isValid(temperature, pressure))
{
- DUNE_THROW(NumericalProblem,
+ OPM_THROW(NumericalProblem,
"Heat capacity of steam is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
}
@@ -541,7 +542,7 @@ public:
{
if (!Region2::isValid(temperature, pressure))
{
- DUNE_THROW(NumericalProblem,
+ OPM_THROW(NumericalProblem,
"Density of steam is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
}
@@ -662,7 +663,7 @@ public:
{
if (!Region1::isValid(temperature, pressure))
{
- DUNE_THROW(NumericalProblem,
+ OPM_THROW(NumericalProblem,
"Density of water is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
}
@@ -761,7 +762,7 @@ public:
{
if (!Region2::isValid(temperature, pressure))
{
- DUNE_THROW(NumericalProblem,
+ OPM_THROW(NumericalProblem,
"Viscosity of steam is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
}
@@ -785,7 +786,7 @@ public:
{
if (!Region1::isValid(temperature, pressure))
{
- DUNE_THROW(NumericalProblem,
+ OPM_THROW(NumericalProblem,
"Viscosity of water is only implemented for temperatures below 623.15K and "
"pressures below 100MPa. (T = " << temperature << ", p=" << pressure);
};
@@ -818,7 +819,7 @@ public:
|| (pressure <= 150e6 && ((523.15 < temperature) || (temperature > 673.15)) )
|| (pressure <= 100e6 && ((673.15 < temperature) || (temperature > 1073.15)) ) )
{
- DUNE_THROW(NumericalProblem,
+ OPM_THROW(NumericalProblem,
"Liquid thermal conductivity of H2O for "
<< "T="< 673.15)) )
|| (pressure <= 100e6 && ((673.15 < temperature) || (temperature > 1073.15)) ) )
{
- DUNE_THROW(NumericalProblem,
+ OPM_THROW(NumericalProblem,
"Gas thermal conductivity of H2O for "
" T="<
#include
-#include
+#include
+#include
namespace Opm
{
@@ -109,7 +110,7 @@ public:
Scalar temperature = iT * (tempMax_ - tempMin_)/(nTemp_ - 1) + tempMin_;
try { vaporPressure_[iT] = RawComponent::vaporPressure(temperature); }
- catch (Dune::NotImplemented) { vaporPressure_[iT] = NaN; }
+ catch (std::exception) { vaporPressure_[iT] = NaN; }
catch (NumericalProblem e) { vaporPressure_[iT] = NaN; };
Scalar pgMax = maxGasPressure_(iT);
@@ -122,23 +123,23 @@ public:
unsigned i = iT + iP*nTemp_;
try { gasEnthalpy_[i] = RawComponent::gasEnthalpy(temperature, pressure); }
- catch (Dune::NotImplemented) { gasEnthalpy_[i] = NaN; }
+ catch (std::exception) { gasEnthalpy_[i] = NaN; }
catch (NumericalProblem) { gasEnthalpy_[i] = NaN; };
try { gasHeatCapacity_[i] = RawComponent::gasHeatCapacity(temperature, pressure); }
- catch (Dune::NotImplemented) { gasHeatCapacity_[i] = NaN; }
+ catch (std::exception) { gasHeatCapacity_[i] = NaN; }
catch (NumericalProblem) { gasHeatCapacity_[i] = NaN; };
try { gasDensity_[i] = RawComponent::gasDensity(temperature, pressure); }
- catch (Dune::NotImplemented) { gasDensity_[i] = NaN; }
+ catch (std::exception) { gasDensity_[i] = NaN; }
catch (NumericalProblem) { gasDensity_[i] = NaN; };
try { gasViscosity_[i] = RawComponent::gasViscosity(temperature, pressure); }
- catch (Dune::NotImplemented) { gasViscosity_[i] = NaN; }
+ catch (std::exception) { gasViscosity_[i] = NaN; }
catch (NumericalProblem) { gasViscosity_[i] = NaN; };
try { gasThermalConductivity_[i] = RawComponent::gasThermalConductivity(temperature, pressure); }
- catch (Dune::NotImplemented) { gasThermalConductivity_[i] = NaN; }
+ catch (std::exception) { gasThermalConductivity_[i] = NaN; }
catch (NumericalProblem) { gasThermalConductivity_[i] = NaN; };
};
@@ -150,23 +151,23 @@ public:
unsigned i = iT + iP*nTemp_;
try { liquidEnthalpy_[i] = RawComponent::liquidEnthalpy(temperature, pressure); }
- catch (Dune::NotImplemented) { liquidEnthalpy_[i] = NaN; }
+ catch (std::exception) { liquidEnthalpy_[i] = NaN; }
catch (NumericalProblem) { liquidEnthalpy_[i] = NaN; };
try { liquidHeatCapacity_[i] = RawComponent::liquidHeatCapacity(temperature, pressure); }
- catch (Dune::NotImplemented) { liquidHeatCapacity_[i] = NaN; }
+ catch (std::exception) { liquidHeatCapacity_[i] = NaN; }
catch (NumericalProblem) { liquidHeatCapacity_[i] = NaN; };
try { liquidDensity_[i] = RawComponent::liquidDensity(temperature, pressure); }
- catch (Dune::NotImplemented) { liquidDensity_[i] = NaN; }
+ catch (std::exception) { liquidDensity_[i] = NaN; }
catch (NumericalProblem) { liquidDensity_[i] = NaN; };
try { liquidViscosity_[i] = RawComponent::liquidViscosity(temperature, pressure); }
- catch (Dune::NotImplemented) { liquidViscosity_[i] = NaN; }
+ catch (std::exception) { liquidViscosity_[i] = NaN; }
catch (NumericalProblem) { liquidViscosity_[i] = NaN; };
try { liquidThermalConductivity_[i] = RawComponent::liquidThermalConductivity(temperature, pressure); }
- catch (Dune::NotImplemented) { liquidThermalConductivity_[i] = NaN; }
+ catch (std::exception) { liquidThermalConductivity_[i] = NaN; }
catch (NumericalProblem) { liquidThermalConductivity_[i] = NaN; };
}
}
@@ -595,13 +596,13 @@ private:
#if 0 && !defined NDEBUG
if(!(0 <= alphaT && alphaT <= 1.0))
- DUNE_THROW(NumericalProblem, "Temperature out of range: "
+ OPM_THROW(NumericalProblem, "Temperature out of range: "
<< "T=" << T << " range: [" << tempMin_ << ", " << tempMax_ << "]");
if(!(0 <= alphaP1 && alphaP1 <= 1.0))
- DUNE_THROW(NumericalProblem, "First liquid pressure out of range: "
+ OPM_THROW(NumericalProblem, "First liquid pressure out of range: "
<< "p=" << p << " range: [" << minLiquidPressure_(tempIdx_(T)) << ", " << maxLiquidPressure_(tempIdx_(T)) << "]");
if(!(0 <= alphaP2 && alphaP2 <= 1.0))
- DUNE_THROW(NumericalProblem, "Second liquid pressure out of range: "
+ OPM_THROW(NumericalProblem, "Second liquid pressure out of range: "
<< "p=" << p << " range: [" << minLiquidPressure_(tempIdx_(T) + 1) << ", " << maxLiquidPressure_(tempIdx_(T) + 1) << "]");
#endif
@@ -633,13 +634,13 @@ private:
#if 0 && !defined NDEBUG
if(!(0 <= alphaT && alphaT <= 1.0))
- DUNE_THROW(NumericalProblem, "Temperature out of range: "
+ OPM_THROW(NumericalProblem, "Temperature out of range: "
<< "T=" << T << " range: [" << tempMin_ << ", " << tempMax_ << "]");
if(!(0 <= alphaP1 && alphaP1 <= 1.0))
- DUNE_THROW(NumericalProblem, "First gas pressure out of range: "
+ OPM_THROW(NumericalProblem, "First gas pressure out of range: "
<< "p=" << p << " range: [" << minGasPressure_(tempIdx_(T)) << ", " << maxGasPressure_(tempIdx_(T)) << "]");
if(!(0 <= alphaP2 && alphaP2 <= 1.0))
- DUNE_THROW(NumericalProblem, "Second gas pressure out of range: "
+ OPM_THROW(NumericalProblem, "Second gas pressure out of range: "
<< "p=" << p << " range: [" << minGasPressure_(tempIdx_(T) + 1) << ", " << maxGasPressure_(tempIdx_(T) + 1) << "]");
#endif
diff --git a/opm/material/components/xylene.hh b/opm/material/components/xylene.hh
index d87c5cf25..42643ba8f 100644
--- a/opm/material/components/xylene.hh
+++ b/opm/material/components/xylene.hh
@@ -73,17 +73,13 @@ public:
* \brief Returns the temperature \f$\mathrm{[K]}\f$ at xylene's triple point.
*/
static const Scalar tripleTemperature()
- {
- DUNE_THROW(Dune::NotImplemented, "tripleTemperature for xylene");
- }
+ { OPM_THROW(std::runtime_error, "Not implemented: tripleTemperature for xylene"); }
/*!
* \brief Returns the pressure \f$\mathrm{[Pa]}\f$ at xylene's triple point.
*/
static const Scalar triplePressure()
- {
- DUNE_THROW(Dune::NotImplemented, "triplePressure for xylene");
- }
+ { OPM_THROW(std::runtime_error, "Not implemented: triplePressure for xylene"); }
/*!
* \brief The saturation vapor pressure in \f$\mathrm{[Pa]}\f$ of pure xylene
diff --git a/opm/material/constraintsolvers/compositionfromfugacities.hh b/opm/material/constraintsolvers/compositionfromfugacities.hh
index 61bac0d8b..b83bb4de6 100644
--- a/opm/material/constraintsolvers/compositionfromfugacities.hh
+++ b/opm/material/constraintsolvers/compositionfromfugacities.hh
@@ -26,8 +26,9 @@
#include
#include
-#include
-#include
+#include
+#include
+#include
#include
@@ -165,12 +166,12 @@ public:
}
}
- DUNE_THROW(NumericalProblem,
- "Calculating the " << FluidSystem::phaseName(phaseIdx)
- << "Phase composition failed. Initial {x} = {"
- << xInit
- << "}, {fug_t} = {" << targetFug << "}, p = " << fluidState.pressure(phaseIdx)
- << ", T = " << fluidState.temperature(phaseIdx));
+ OPM_THROW(Opm::NumericalProblem,
+ "Calculating the " << FluidSystem::phaseName(phaseIdx)
+ << "Phase composition failed. Initial {x} = {"
+ << xInit
+ << "}, {fug_t} = {" << targetFug << "}, p = " << fluidState.pressure(phaseIdx)
+ << ", T = " << fluidState.temperature(phaseIdx));
}
diff --git a/opm/material/constraintsolvers/computefromreferencephase.hh b/opm/material/constraintsolvers/computefromreferencephase.hh
index 8bb312902..97d3b58df 100644
--- a/opm/material/constraintsolvers/computefromreferencephase.hh
+++ b/opm/material/constraintsolvers/computefromreferencephase.hh
@@ -25,8 +25,9 @@
#include
-#include
-#include
+#include
+#include
+#include
#include
diff --git a/opm/material/constraintsolvers/immiscibleflash.hh b/opm/material/constraintsolvers/immiscibleflash.hh
index c2e3b2d3c..7f74f17f3 100644
--- a/opm/material/constraintsolvers/immiscibleflash.hh
+++ b/opm/material/constraintsolvers/immiscibleflash.hh
@@ -23,8 +23,9 @@
#ifndef OPM_IMMISCIBLE_FLASH_HH
#define OPM_IMMISCIBLE_FLASH_HH
-#include
-#include
+#include
+#include
+#include
#include
#include
@@ -214,10 +215,10 @@ public:
std::cout << "\n";
*/
- DUNE_THROW(NumericalProblem,
- "Flash calculation failed."
- " {c_alpha^kappa} = {" << globalMolarities << "}, T = "
- << fluidState.temperature(/*phaseIdx=*/0));
+ OPM_THROW(Opm::NumericalProblem,
+ "Flash calculation failed."
+ " {c_alpha^kappa} = {" << globalMolarities << "}, T = "
+ << fluidState.temperature(/*phaseIdx=*/0));
}
diff --git a/opm/material/constraintsolvers/misciblemultiphasecomposition.hh b/opm/material/constraintsolvers/misciblemultiphasecomposition.hh
index 924256fe9..b44223707 100644
--- a/opm/material/constraintsolvers/misciblemultiphasecomposition.hh
+++ b/opm/material/constraintsolvers/misciblemultiphasecomposition.hh
@@ -26,8 +26,9 @@
#include
#include
-#include
-#include
+#include
+#include
+#include
namespace Opm {
@@ -242,8 +243,8 @@ public:
M.solve(x, b);
}
catch (const Dune::FMatrixError &e) {
- DUNE_THROW(NumericalProblem,
- "Numerical problem in MiscibleMultiPhaseComposition::solve(): " << NumericalProblem(e.what()) << "; M="<
+#include
namespace Opm {
/*!
diff --git a/opm/material/fluidmatrixinteractions/2p/regularizedlinearmaterial.hh b/opm/material/fluidmatrixinteractions/2p/regularizedlinearmaterial.hh
index e37b4b73e..7c31c0cad 100644
--- a/opm/material/fluidmatrixinteractions/2p/regularizedlinearmaterial.hh
+++ b/opm/material/fluidmatrixinteractions/2p/regularizedlinearmaterial.hh
@@ -28,7 +28,7 @@
#include "linearmaterial.hh"
#include "regularizedlinearmaterialparams.hh"
-#include
+#include
namespace Opm {
/*!
diff --git a/opm/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh b/opm/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh
index eb555e7a9..7a92f09d7 100644
--- a/opm/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh
+++ b/opm/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh
@@ -28,7 +28,7 @@
#include "vangenuchten.hh"
#include "regularizedvangenuchtenparams.hh"
-#include
+#include
#include
diff --git a/opm/material/fluidmatrixinteractions/3p/3pparkervangenuchten.hh b/opm/material/fluidmatrixinteractions/3p/3pparkervangenuchten.hh
index 9239d7264..bace198e1 100644
--- a/opm/material/fluidmatrixinteractions/3p/3pparkervangenuchten.hh
+++ b/opm/material/fluidmatrixinteractions/3p/3pparkervangenuchten.hh
@@ -52,7 +52,7 @@ public:
*/
static Scalar pC(const Params ¶ms, Scalar Sw)
{
- DUNE_THROW(Dune::NotImplemented, "Capillary pressures for three phases is not so simple! Use pCGN, pCNW, and pcGW");
+ OPM_THROW(Dune::NotImplemented, "Capillary pressures for three phases is not so simple! Use pCGN, pCNW, and pcGW");
}
static Scalar pCGW(const Params ¶ms, Scalar Sw)
@@ -208,7 +208,7 @@ public:
*/
static Scalar Sw(const Params ¶ms, Scalar pC)
{
- DUNE_THROW(Dune::NotImplemented, "Sw(pc) for three phases not implemented! Do it yourself!");
+ OPM_THROW(Dune::NotImplemented, "Sw(pc) for three phases not implemented! Do it yourself!");
}
/*!
@@ -218,7 +218,7 @@ public:
*/
static Scalar dpC_dSw(const Params ¶ms, Scalar Sw)
{
- DUNE_THROW(Dune::NotImplemented, "dpC/dSw for three phases not implemented! Do it yourself!");
+ OPM_THROW(Dune::NotImplemented, "dpC/dSw for three phases not implemented! Do it yourself!");
}
/*!
@@ -227,7 +227,7 @@ public:
*/
static Scalar dSw_dpC(const Params ¶ms, Scalar pC)
{
- DUNE_THROW(Dune::NotImplemented, "dSw/dpC for three phases not implemented! Do it yourself!");
+ OPM_THROW(Dune::NotImplemented, "dSw/dpC for three phases not implemented! Do it yourself!");
}
/*!
diff --git a/opm/material/fluidmatrixinteractions/3p/3pparkervangenuchtenparams.hh b/opm/material/fluidmatrixinteractions/3p/3pparkervangenuchtenparams.hh
index a414e1d28..6eb0e8a59 100644
--- a/opm/material/fluidmatrixinteractions/3p/3pparkervangenuchtenparams.hh
+++ b/opm/material/fluidmatrixinteractions/3p/3pparkervangenuchtenparams.hh
@@ -27,7 +27,7 @@
#include
-#include
+#include
namespace Opm {
/*!
@@ -128,7 +128,7 @@ public:
return Sgr_;
break;
};
- DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
+ OPM_THROW(std::logic_error, "Invalid phase index " << phaseIdx);
}
/*!
diff --git a/opm/material/fluidmatrixinteractions/mp/3padapter.hh b/opm/material/fluidmatrixinteractions/mp/3padapter.hh
index 877b712a6..4f916cf05 100644
--- a/opm/material/fluidmatrixinteractions/mp/3padapter.hh
+++ b/opm/material/fluidmatrixinteractions/mp/3padapter.hh
@@ -20,12 +20,12 @@
* \file
* \copydoc Opm::ThreePAdapter
*/
-#ifndef OPM_MP_3P_ADAPTER_HH
-#define OPM_MP_3P_ADAPTER_HH
+#ifndef OPM_3P_ADAPTER_HPP
+#define OPM_3P_ADAPTER_HPP
-#include
-
-#include
+#include
+#include
+#include
#include
@@ -80,7 +80,7 @@ public:
static void saturations(ContainerT &values,
const Params ¶ms,
const FluidState &fluidState)
- { DUNE_THROW(Dune::NotImplemented, "Inverse capillary pressure curves"); }
+ { OPM_THROW(std::runtime_error, "Not implemented: Inverse capillary pressure curves"); }
/*!
* \brief The relative permeability of all phases.
diff --git a/opm/material/fluidmatrixinteractions/mp/mplinearmaterial.hh b/opm/material/fluidmatrixinteractions/mp/mplinearmaterial.hh
index 233796873..3834b44f3 100644
--- a/opm/material/fluidmatrixinteractions/mp/mplinearmaterial.hh
+++ b/opm/material/fluidmatrixinteractions/mp/mplinearmaterial.hh
@@ -25,9 +25,10 @@
#include "mplinearmaterialparams.hh"
-#include
+#include
-#include
+#include
+#include
#include
@@ -93,7 +94,7 @@ public:
const Params ¶ms,
const FluidState &state)
{
- DUNE_THROW(Dune::NotImplemented, "MpLinearMaterial::saturations()");
+ OPM_THROW(std::runtime_error, "Not implemented: MpLinearMaterial::saturations()");
}
/*!
diff --git a/opm/material/fluidmatrixinteractions/mp/nullmateriallaw.hh b/opm/material/fluidmatrixinteractions/mp/nullmateriallaw.hh
index e7008cab7..fcbd6e5e9 100644
--- a/opm/material/fluidmatrixinteractions/mp/nullmateriallaw.hh
+++ b/opm/material/fluidmatrixinteractions/mp/nullmateriallaw.hh
@@ -25,7 +25,8 @@
#include "nullmateriallawparams.hh"
-#include
+#include
+#include
#include
@@ -69,7 +70,7 @@ public:
const Params ¶ms,
const FluidState &state)
{
- DUNE_THROW(Dune::NotImplemented, "MpLinearMaterial::saturations()");
+ OPM_THROW(std::runtime_error, "Not implemented: MpLinearMaterial::saturations()");
}
/*!
diff --git a/opm/material/fluidstates/compositionalfluidstate.hh b/opm/material/fluidstates/compositionalfluidstate.hh
index ef5b8a582..cd3ab54e4 100644
--- a/opm/material/fluidstates/compositionalfluidstate.hh
+++ b/opm/material/fluidstates/compositionalfluidstate.hh
@@ -28,7 +28,7 @@
#include "modularfluidstate.hh"
-#include
+#include
#include
namespace Opm {
diff --git a/opm/material/fluidstates/fluidstatecompositionmodules.hh b/opm/material/fluidstates/fluidstatecompositionmodules.hh
index 0618ffd8a..dbd463b3e 100644
--- a/opm/material/fluidstates/fluidstatecompositionmodules.hh
+++ b/opm/material/fluidstates/fluidstatecompositionmodules.hh
@@ -24,9 +24,10 @@
#ifndef OPM_FLUID_STATE_COMPOSITION_MODULES_HH
#define OPM_FLUID_STATE_COMPOSITION_MODULES_HH
-#include
+#include
-#include
+#include
+#include
#include
#include
@@ -238,7 +239,7 @@ protected:
/*!
* \brief Module for the modular fluid state which does not store the
- * compositions but throws Dune::InvalidState instead.
+ * compositions but throws std::logic_error instead.
*/
template
+#include
-#include
+#include
+#include
#include
@@ -105,7 +106,7 @@ protected:
/*!
* \brief Module for the modular fluid state which does not the
- * densitys but throws Dune::InvalidState instead.
+ * densitys but throws std::logic_error instead.
*/
template
+#include
-#include
+#include
+#include
#include
@@ -99,7 +100,7 @@ protected:
/*!
* \brief Module for the modular fluid state which does not store the
- * enthalpies but throws Dune::InvalidState instead.
+ * enthalpies but throws std::logic_error instead.
*/
template
+#include
-#include
+#include
+#include
#include
#include
@@ -173,7 +174,7 @@ protected:
/*!
* \brief Module for the modular fluid state which does not store the
- * fugacitys but throws Dune::InvalidState instead.
+ * fugacitys but throws std::logic_error instead.
*/
template
+#include
-#include
+#include
+#include
#include
@@ -91,7 +92,7 @@ protected:
/*!
* \brief Module for the modular fluid state which does not the
- * pressures but throws Dune::InvalidState instead.
+ * pressures but throws std::logic_error instead.
*/
template
+#include
-#include
+#include
+#include
#include
@@ -90,7 +91,7 @@ protected:
/*!
* \brief Module for the modular fluid state which does not the
- * saturations but throws Dune::InvalidState instead.
+ * saturations but throws std::logic_error instead.
*/
template
+#include
-#include
+#include
+#include
#include
#include
@@ -149,7 +150,7 @@ protected:
/*!
* \brief Module for the modular fluid state which does not the
- * temperatures but throws Dune::InvalidState instead.
+ * temperatures but throws std::logic_error instead.
*/
template
+#include
-#include
+#include
+#include