OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #798 from andlaus/remove_unused_code

Browse Source 
Remove some unused code
This commit is contained in:
Andreas Lauser 2015-04-28 17:45:09 +02:00
commit b18d609396
11 changed files with 1 additions and 4720 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
CMakeLists_files.cmake
View File

@ -160,8 +160,6 @@ list (APPEND TEST_SOURCE_FILES
tests/test_EclipseWriter.cpp
tests/test_EclipseWriteRFTHandler.cpp
tests/test_compressedpropertyaccess.cpp
tests/test_spline.cpp
tests/test_propertysystem.cpp
tests/test_dgbasis.cpp
tests/test_cartgrid.cpp
tests/test_ug.cpp
@ -416,7 +414,6 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/transport/reorder/reordersequence.h
opm/core/transport/reorder/tarjan.h
opm/core/utility/Average.hpp
opm/core/utility/ClassName.hpp
opm/core/utility/CompressedPropertyAccess.hpp
opm/core/utility/DataMap.hpp
opm/core/utility/ErrorMacros.hpp
@ -427,17 +424,13 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/utility/MonotCubicInterpolator.hpp
opm/core/utility/NonuniformTableLinear.hpp
opm/core/utility/NullStream.hpp
opm/core/utility/PolynomialUtils.hpp
opm/core/utility/RegionMapping.hpp
opm/core/utility/RootFinders.hpp
opm/core/utility/SparseTable.hpp
opm/core/utility/SparseVector.hpp
opm/core/utility/Spline.hpp
opm/core/utility/StopWatch.hpp
opm/core/utility/TridiagonalMatrix.hpp
opm/core/utility/UniformTableLinear.hpp
opm/core/utility/Units.hpp
opm/core/utility/Unused.hpp
opm/core/utility/VelocityInterpolation.hpp
opm/core/utility/WachspressCoord.hpp
opm/core/utility/buildUniformMonotoneTable.hpp
@ -458,7 +451,6 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/utility/parameters/tinyxml/tinyxml.h
opm/core/utility/platform_dependent/disable_warnings.h
opm/core/utility/platform_dependent/reenable_warnings.h
opm/core/utility/PropertySystem.hpp
opm/core/utility/share_obj.hpp
opm/core/utility/thresholdPressures.hpp
opm/core/wells/InjectionSpecification.hpp

2
cmake/Modules/ewoms-prereqs.cmake
View File

@ -22,9 +22,7 @@ set (ewoms_DEPS
"dune-geometry REQUIRED"
"dune-grid REQUIRED"
"dune-istl REQUIRED"
"opm-core REQUIRED"
"opm-material REQUIRED"
"opm-parser"
"dune-alugrid"
"dune-cornerpoint"
# valgrind client requests

2
cmake/Modules/opm-material-prereqs.cmake
View File

@ -15,7 +15,7 @@ set (opm-material_DEPS
# compile with C++0x/11 support if available
"CXX11Features REQUIRED"
# prerequisite OPM modules
"opm-core REQUIRED"
"opm-parser"
# DUNE dependency
"dune-common REQUIRED"
)

84
opm/core/utility/ClassName.hpp
View File

@ -1,84 +0,0 @@
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set ts=4 sw=2 et sts=2:
/*****************************************************************************
* Copyright (C) 2012-2013 by Andreas Lauser *
* Copyright (C) 2010 by Oliver Sander *
* Copyright (C) 2011 by Martin Nolte *
* *
* 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, version 2. *
* *
* This file is based on works developed by the DUNE project, see *
* <http://www.dune-project.org>. The following license exception *
* applies to it: *
* *
* As a special exception, you may use the DUNE library without *
* restriction. Specifically, if other files instantiate templates or *
* use macros or inline functions from one or more of the DUNE source *
* files, or you compile one or more of the DUNE source files and link *
* them with other files to produce an executable, this does not by *
* itself cause the resulting executable to be covered by the GNU *
* General Public License. This exception does not however invalidate *
* any other reasons why the executable file might be covered by the *
* GNU General Public License. *
* *
* 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 <http://www.gnu.org/licenses/>. *
*****************************************************************************/
#ifndef OPM_CLASSNAME_HPP
#define OPM_CLASSNAME_HPP
/** \file
* \brief A free function to provide the demangled class name
* of a given object or type as a string
*/
#include <cstdlib>
#include <string>
#include <typeinfo>
#if HAVE_CXA_DEMANGLE
#include <cxxabi.h>
#endif
namespace Opm {
/** \brief Provide the demangled class name of a given object as a string */
template <class T>
std::string className()
{
std::string className = typeid( T ).name();
#if HAVE_CXA_DEMANGLE
int status;
char *demangled = abi::__cxa_demangle( className.c_str(), 0, 0, &status );
if( demangled )
{
className = demangled;
std::free( demangled );
}
#endif // HAVE_CXA_DEMANGLE
return className;
}
#if HAVE_QUAD
// specialize for quad precision floating point values to avoid
// needing a type_info structure
template <>
std::string className<__float128>()
{ return "quad"; }
#endif
/** \brief Provide the demangled class name of a given object as a string */
template <class T>
std::string className(const T &)
{
return className<T>();
}
} // namespace Opm
#endif // OPM_CLASSNAME_HPP

310
opm/core/utility/PolynomialUtils.hpp
View File

@ -1,310 +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-2013 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 <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \brief Define some often used mathematical functions
*/
#ifndef OPM_MATH_HH
#define OPM_MATH_HH
#include <cmath>
#include <algorithm>
namespace Opm
{
/*!
* \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 <class Scalar, class SolContainer>
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 <class Scalar, class SolContainer>
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 <class Scalar, class SolContainer>
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 <class Scalar, class SolContainer>
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;
}
}
#endif

1134
opm/core/utility/PropertySystem.hpp
View File

File diff suppressed because it is too large Load Diff

1762
opm/core/utility/Spline.hpp
View File

File diff suppressed because it is too large Load Diff

849
opm/core/utility/TridiagonalMatrix.hpp
View File

@ -1,849 +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-2013 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 *
* MERCHANTBILITY 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 <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \copydetails Opm::TridiagonalMatrix
*/
#ifndef OPM_TRIDIAGONAL_MATRIX_HH
#define OPM_TRIDIAGONAL_MATRIX_HH
#include <iostream>
#include <vector>
#include <algorithm>
#include <cmath>
#include <assert.h>
namespace Opm {
/*!
* \brief Provides a tridiagonal matrix that also supports non-zero
* entries in the upper right and lower left
*
* The entries in the lower left and upper right are supported to make
* implementing periodic systems easy.
*
* The API of this class is designed to be close to the one used by
* the DUNE matrix classes.
*/
template <class Scalar>
class TridiagonalMatrix
{
struct TridiagRow_ {
TridiagRow_(TridiagonalMatrix &m, size_t rowIdx)
: matrix_(m)
, rowIdx_(rowIdx)
{};
Scalar &operator[](size_t colIdx)
{ return matrix_.at(rowIdx_, colIdx); }
Scalar operator[](size_t colIdx) const
{ return matrix_.at(rowIdx_, colIdx); }
/*!
* \brief Prefix increment operator
*/
TridiagRow_ &operator++()
{ ++ rowIdx_; return *this; }
/*!
* \brief Prefix decrement operator
*/
TridiagRow_ &operator--()
{ -- rowIdx_; return *this; }
/*!
* \brief Comparision operator
*/
bool operator==(const TridiagRow_ &other) const
{ return other.rowIdx_ == rowIdx_ && &other.matrix_ == &matrix_; }
/*!
* \brief Comparision operator
*/
bool operator!=(const TridiagRow_ &other) const
{ return !operator==(other); }
/*!
* \brief Dereference operator
*/
TridiagRow_ &operator*()
{ return *this; }
/*!
* \brief Return the row index of the this row.
*
* 0 is the first row.
*/
size_t index() const
{ return rowIdx_; }
private:
TridiagonalMatrix &matrix_;
mutable size_t rowIdx_;
};
public:
typedef Scalar FieldType;
typedef TridiagRow_ RowType;
typedef size_t SizeType;
typedef TridiagRow_ iterator;
typedef TridiagRow_ const_iterator;
explicit TridiagonalMatrix(int numRows = 0)
{
resize(numRows);
};
TridiagonalMatrix(int numRows, Scalar value)
{
resize(numRows);
this->operator=(value);
};
/*!
* \brief Copy constructor.
*/
TridiagonalMatrix(const TridiagonalMatrix &source)
{ *this = source; };
/*!
* \brief Return the number of rows/columns of the matrix.
*/
size_t size() const
{ return diag_[0].size(); }
/*!
* \brief Return the number of rows of the matrix.
*/
size_t rows() const
{ return size(); }
/*!
* \brief Return the number of columns of the matrix.
*/
size_t cols() const
{ return size(); }
/*!
* \brief Change the number of rows of the matrix.
*/
void resize(size_t n)
{
if (n == size())
return;
for (int diagIdx = 0; diagIdx < 3; ++ diagIdx)
diag_[diagIdx].resize(n);
}
/*!
* \brief Access an entry.
*/
Scalar &at(size_t rowIdx, size_t colIdx)
{
size_t n = size();
// special cases
if (n > 2) {
if (rowIdx == 0 && colIdx == n - 1)
return diag_[2][0];
if (rowIdx == n - 1 && colIdx == 0)
return diag_[0][n - 1];
}
int diagIdx = 1 + colIdx - rowIdx;
// make sure that the requested column is in range
assert(0 <= diagIdx && diagIdx < 3);
return diag_[diagIdx][colIdx];
}
/*!
* \brief Access an entry.
*/
Scalar at(size_t rowIdx, size_t colIdx) const
{
int n = size();
// special cases
if (rowIdx == 0 && colIdx == n - 1)
return diag_[2][0];
if (rowIdx == n - 1 && colIdx == 0)
return diag_[0][n - 1];
int diagIdx = 1 + colIdx - rowIdx;
// make sure that the requested column is in range
assert(0 <= diagIdx && diagIdx < 3);
return diag_[diagIdx][colIdx];
}
/*!
* \brief Assignment operator from another tridiagonal matrix.
*/
TridiagonalMatrix &operator=(const TridiagonalMatrix &source)
{
for (int diagIdx = 0; diagIdx < 3; ++ diagIdx)
diag_[diagIdx] = source.diag_[diagIdx];
return *this;
}
/*!
* \brief Assignment operator from a Scalar.
*/
TridiagonalMatrix &operator=(Scalar value)
{
for (int diagIdx = 0; diagIdx < 3; ++ diagIdx)
diag_[diagIdx].assign(size(), value);
return *this;
}
/*!
* \begin Iterator for the first row
*/
iterator begin()
{ return TridiagRow_(*this, 0); }
/*!
* \begin Const iterator for the first row
*/
const_iterator begin() const
{ return TridiagRow_(const_cast<TridiagonalMatrix&>(*this), 0); }
/*!
* \begin Const iterator for the next-to-last row
*/
const_iterator end() const
{ return TridiagRow_(const_cast<TridiagonalMatrix&>(*this), size()); }
/*!
* \brief Row access operator.
*/
TridiagRow_ operator[](size_t rowIdx)
{ return TridiagRow_(*this, rowIdx); }
/*!
* \brief Row access operator.
*/
const TridiagRow_ operator[](size_t rowIdx) const
{ return TridiagRow_(*this, rowIdx); }
/*!
* \brief Multiplication with a Scalar
*/
TridiagonalMatrix &operator*=(Scalar alpha)
{
int n = size();
for (int diagIdx = 0; diagIdx < 3; ++ diagIdx) {
for (int i = 0; i < n; ++i) {
diag_[diagIdx][i] *= alpha;
}
}
return *this;
}
/*!
* \brief Division by a Scalar
*/
TridiagonalMatrix &operator/=(Scalar alpha)
{
alpha = 1.0/alpha;
int n = size();
for (int diagIdx = 0; diagIdx < 3; ++ diagIdx) {
for (int i = 0; i < n; ++i) {
diag_[diagIdx][i] *= alpha;
}
}
return *this;
}
/*!
* \brief Subtraction operator
*/
TridiagonalMatrix &operator-=(const TridiagonalMatrix &other)
{ return axpy(-1.0, other); }
/*!
* \brief Addition operator
*/
TridiagonalMatrix &operator+=(const TridiagonalMatrix &other)
{ return axpy(1.0, other); }
/*!
* \brief Multiply and add the matrix entries of another
* tridiagonal matrix.
*
* This means that
* \code
* A.axpy(alpha, B)
* \endcode
* is equivalent to
* \code
* A += alpha*C
* \endcode
*/
TridiagonalMatrix &axpy(Scalar alpha, const TridiagonalMatrix &other)
{
assert(size() == other.size());
int n = size();
for (int diagIdx = 0; diagIdx < 3; ++ diagIdx)
for (int i = 0; i < n; ++ i)
diag_[diagIdx][i] += alpha * other[diagIdx][i];
return *this;
}
/*!
* \brief Matrix-vector product
*
* This means that
* \code
* A.mv(x, y)
* \endcode
* is equivalent to
* \code
* y = A*x
* \endcode
*/
template<class Vector>
void mv(const Vector &source, Vector &dest) const
{
assert(source.size() == size());
assert(dest.size() == size());
assert(size() > 1);
// deal with rows 1 .. n-2
int n = size();
for (int i = 1; i < n - 1; ++ i) {
dest[i] =
diag_[0][i - 1]*source[i-1] +
diag_[1][i]*source[i] +
diag_[2][i + 1]*source[i + 1];
}
// rows 0 and n-1
dest[0] =
diag_[1][0]*source[0] +
diag_[2][1]*source[1] +
diag_[2][0]*source[n - 1];
dest[n - 1] =
diag_[0][n-1]*source[0] +
diag_[0][n-2]*source[n-2] +
diag_[1][n-1]*source[n-1];
}
/*!
* \brief Additive matrix-vector product
*
* This means that
* \code
* A.umv(x, y)
* \endcode
* is equivalent to
* \code
* y += A*x
* \endcode
*/
template<class Vector>
void umv(const Vector &source, Vector &dest) const
{
assert(source.size() == size());
assert(dest.size() == size());
assert(size() > 1);
// deal with rows 1 .. n-2
int n = size();
for (int i = 1; i < n - 1; ++ i) {
dest[i] +=
diag_[0][i - 1]*source[i-1] +
diag_[1][i]*source[i] +
diag_[2][i + 1]*source[i + 1];
}
// rows 0 and n-1
dest[0] +=
diag_[1][0]*source[0] +
diag_[2][1]*source[1] +
diag_[2][0]*source[n - 1];
dest[n - 1] +=
diag_[0][n-1]*source[0] +
diag_[0][n-2]*source[n-2] +
diag_[1][n-1]*source[n-1];
}
/*!
* \brief Subtractive matrix-vector product
*
* This means that
* \code
* A.mmv(x, y)
* \endcode
* is equivalent to
* \code
* y -= A*x
* \endcode
*/
template<class Vector>
void mmv(const Vector &source, Vector &dest) const
{
assert(source.size() == size());
assert(dest.size() == size());
assert(size() > 1);
// deal with rows 1 .. n-2
int n = size();
for (int i = 1; i < n - 1; ++ i) {
dest[i] -=
diag_[0][i - 1]*source[i-1] +
diag_[1][i]*source[i] +
diag_[2][i + 1]*source[i + 1];
}
// rows 0 and n-1
dest[0] -=
diag_[1][0]*source[0] +
diag_[2][1]*source[1] +
diag_[2][0]*source[n - 1];
dest[n - 1] -=
diag_[0][n-1]*source[0] +
diag_[0][n-2]*source[n-2] +
diag_[1][n-1]*source[n-1];
}
/*!
* \brief Scaled additive matrix-vector product
*
* This means that
* \code
* A.usmv(x, y)
* \endcode
* is equivalent to
* \code
* y += alpha*(A*x)
* \endcode
*/
template<class Vector>
void usmv(Scalar alpha, const Vector &source, Vector &dest) const
{
assert(source.size() == size());
assert(dest.size() == size());
assert(size() > 1);
// deal with rows 1 .. n-2
int n = size();
for (int i = 1; i < n - 1; ++ i) {
dest[i] +=
alpha*(
diag_[0][i - 1]*source[i-1] +
diag_[1][i]*source[i] +
diag_[2][i + 1]*source[i + 1]);
}
// rows 0 and n-1
dest[0] +=
alpha*(
diag_[1][0]*source[0] +
diag_[2][1]*source[1] +
diag_[2][0]*source[n - 1]);
dest[n - 1] +=
alpha*(
diag_[0][n-1]*source[0] +
diag_[0][n-2]*source[n-2] +
diag_[1][n-1]*source[n-1]);
}
/*!
* \brief Transposed matrix-vector product
*
* This means that
* \code
* A.mtv(x, y)
* \endcode
* is equivalent to
* \code
* y = A^T*x
* \endcode
*/
template<class Vector>
void mtv(const Vector &source, Vector &dest) const
{
assert(source.size() == size());
assert(dest.size() == size());
assert(size() > 1);
// deal with rows 1 .. n-2
int n = size();
for (int i = 1; i < n - 1; ++ i) {
dest[i] =
diag_[2][i + 1]*source[i-1] +
diag_[1][i]*source[i] +
diag_[0][i - 1]*source[i + 1];
}
// rows 0 and n-1
dest[0] =
diag_[1][0]*source[0] +
diag_[0][1]*source[1] +
diag_[0][n-1]*source[n - 1];
dest[n - 1] =
diag_[2][0]*source[0] +
diag_[2][n-1]*source[n-2] +
diag_[1][n-1]*source[n-1];
}
/*!
* \brief Transposed additive matrix-vector product
*
* This means that
* \code
* A.umtv(x, y)
* \endcode
* is equivalent to
* \code
* y += A^T*x
* \endcode
*/
template<class Vector>
void umtv(const Vector &source, Vector &dest) const
{
assert(source.size() == size());
assert(dest.size() == size());
assert(size() > 1);
// deal with rows 1 .. n-2
int n = size();
for (int i = 1; i < n - 1; ++ i) {
dest[i] +=
diag_[2][i + 1]*source[i-1] +
diag_[1][i]*source[i] +
diag_[0][i - 1]*source[i + 1];
}
// rows 0 and n-1
dest[0] +=
diag_[1][0]*source[0] +
diag_[0][1]*source[1] +
diag_[0][n-1]*source[n - 1];
dest[n - 1] +=
diag_[2][0]*source[0] +
diag_[2][n-1]*source[n-2] +
diag_[1][n-1]*source[n-1];
}
/*!
* \brief Transposed subtractive matrix-vector product
*
* This means that
* \code
* A.mmtv(x, y)
* \endcode
* is equivalent to
* \code
* y -= A^T*x
* \endcode
*/
template<class Vector>
void mmtv (const Vector &source, Vector &dest) const
{
assert(source.size() == size());
assert(dest.size() == size());
assert(size() > 1);
// deal with rows 1 .. n-2
int n = size();
for (int i = 1; i < n - 1; ++ i) {
dest[i] -=
diag_[2][i + 1]*source[i-1] +
diag_[1][i]*source[i] +
diag_[0][i - 1]*source[i + 1];
}
// rows 0 and n-1
dest[0] -=
diag_[1][0]*source[0] +
diag_[0][1]*source[1] +
diag_[0][n-1]*source[n - 1];
dest[n - 1] -=
diag_[2][0]*source[0] +
diag_[2][n-1]*source[n-2] +
diag_[1][n-1]*source[n-1];
}
/*!
* \brief Transposed scaled additive matrix-vector product
*
* This means that
* \code
* A.umtv(alpha, x, y)
* \endcode
* is equivalent to
* \code
* y += alpha*A^T*x
* \endcode
*/
template<class Vector>
void usmtv(Scalar alpha, const Vector &source, Vector &dest) const
{
assert(source.size() == size());
assert(dest.size() == size());
assert(size() > 1);
// deal with rows 1 .. n-2
int n = size();
for (int i = 1; i < n - 1; ++ i) {
dest[i] +=
alpha*(
diag_[2][i + 1]*source[i-1] +
diag_[1][i]*source[i] +
diag_[0][i - 1]*source[i + 1]);
}
// rows 0 and n-1
dest[0] +=
alpha*(
diag_[1][0]*source[0] +
diag_[0][1]*source[1] +
diag_[0][n-1]*source[n - 1]);
dest[n - 1] +=
alpha*(
diag_[2][0]*source[0] +
diag_[2][n-1]*source[n-2] +
diag_[1][n-1]*source[n-1]);
}
/*!
* \brief Calculate the frobenius norm
*
* i.e., the square root of the sum of all squared entries. This
* corresponds to the euclidean norm for vectors.
*/
Scalar frobeniusNorm() const
{ return std::sqrt(frobeniusNormSquared()); }
/*!
* \brief Calculate the squared frobenius norm
*
* i.e., the sum of all squared entries.
*/
Scalar frobeniusNormSquared() const
{
Scalar result = 0;
int n = size();
for (int i = 0; i < n; ++ i)
for (int diagIdx = 0; diagIdx < 3; ++ diagIdx)
result += diag_[diagIdx][i];
return result;
}
/*!
* \brief Calculate the infinity norm
*
* i.e., the maximum of the sum of the absolute values of all rows.
*/
Scalar infinityNorm() const
{
Scalar result=0;
// deal with rows 1 .. n-2
int n = size();
for (int i = 1; i < n - 1; ++ i) {
result = std::max(result,
std::abs(diag_[0][i - 1]) +
std::abs(diag_[1][i]) +
std::abs(diag_[2][i + 1]));
}
// rows 0 and n-1
result = std::max(result,
std::abs(diag_[1][0]) +
std::abs(diag_[2][1]) +
std::abs(diag_[2][0]));
result = std::max(result,
std::abs(diag_[0][n-1]) +
std::abs(diag_[0][n-2]) +
std::abs(diag_[1][n-2]));
return result;
}
/*!
* \brief Calculate the solution for a linear system of equations
*
* i.e., calculate x, so that it solves Ax = b, where A is a
* tridiagonal matrix.
*/
template <class XVector, class BVector>
void solve(XVector &x, const BVector &b) const
{
if (size() > 2 && diag_[2][0] != 0)
solveWithUpperRight_(x, b);
else
solveWithoutUpperRight_(x, b);
}
/*!
* \brief Print the matrix to a given output stream.
*/
void print(std::ostream &os = std::cout) const
{
int n = size();
// row 0
os << at(0, 0) << "\t"
<< at(0, 1) << "\t";
if (n > 3)
os << "\t";
if (n > 2)
os << at(0, n-1);
os << "\n";
// row 1 .. n - 2
for (int rowIdx = 1; rowIdx < n-1; ++rowIdx) {
if (rowIdx > 1)
os << "\t";
if (rowIdx == n - 2)
os << "\t";
os << at(rowIdx, rowIdx - 1) << "\t"
<< at(rowIdx, rowIdx) << "\t"
<< at(rowIdx, rowIdx + 1) << "\n";
}
// row n - 1
if (n > 2)
os << at(n-1, 0) << "\t";
if (n > 3)
os << "\t";
if (n > 4)
os << "\t";
os << at(n-1, n-2) << "\t"
<< at(n-1, n-1) << "\n";
}
private:
template <class XVector, class BVector>
void solveWithUpperRight_(XVector &x, const BVector &b) const
{
size_t n = size();
std::vector<Scalar> lowerDiag(diag_[0]), mainDiag(diag_[1]), upperDiag(diag_[2]), lastColumn(n);
std::vector<Scalar> bStar(n);
std::copy(b.begin(), b.end(), bStar.begin());
lastColumn[0] = upperDiag[0];
// forward elimination
for (size_t i = 1; i < n; ++i) {
Scalar alpha = lowerDiag[i - 1]/mainDiag[i - 1];
lowerDiag[i - 1] -= alpha * mainDiag[i - 1];
mainDiag[i] -= alpha * upperDiag[i];
bStar[i] -= alpha * bStar[i - 1];
};
// deal with the last row if the entry on the lower left is not zero
if (lowerDiag[n - 1] != 0.0 && size() > 2) {
Scalar lastRow = lowerDiag[n - 1];
for (size_t i = 0; i < n - 1; ++i) {
Scalar alpha = lastRow/mainDiag[i];
lastRow = - alpha*upperDiag[i + 1];
bStar[n - 1] -= alpha * bStar[i];
}
mainDiag[n-1] += lastRow;
}
// backward elimination
x[n - 1] = bStar[n - 1]/mainDiag[n-1];
for (int i = n - 2; i >= 0; --i) {
x[i] = (bStar[i] - x[i + 1]*upperDiag[i+1] - x[n-1]*lastColumn[i])/mainDiag[i];
}
}
template <class XVector, class BVector>
void solveWithoutUpperRight_(XVector &x, const BVector &b) const
{
size_t n = size();
std::vector<Scalar> lowerDiag(diag_[0]), mainDiag(diag_[1]), upperDiag(diag_[2]);
std::vector<Scalar> bStar(n);
std::copy(b.begin(), b.end(), bStar.begin());
// forward elimination
for (size_t i = 1; i < n; ++i) {
Scalar alpha = lowerDiag[i - 1]/mainDiag[i - 1];
lowerDiag[i - 1] -= alpha * mainDiag[i - 1];
mainDiag[i] -= alpha * upperDiag[i];
bStar[i] -= alpha * bStar[i - 1];
};
// deal with the last row if the entry on the lower left is not zero
if (lowerDiag[n - 1] != 0.0 && size() > 2) {
Scalar lastRow = lowerDiag[n - 1];
for (size_t i = 0; i < n - 1; ++i) {
Scalar alpha = lastRow/mainDiag[i];
lastRow = - alpha*upperDiag[i + 1];
bStar[n - 1] -= alpha * bStar[i];
}
mainDiag[n-1] += lastRow;
}
// backward elimination
x[n - 1] = bStar[n - 1]/mainDiag[n-1];
for (int i = n - 2; i >= 0; --i) {
x[i] = (bStar[i] - x[i + 1]*upperDiag[i+1])/mainDiag[i];
}
}
mutable std::vector<Scalar> diag_[3];
};
} // namespace Opm
template <class Scalar>
std::ostream &operator<<(std::ostream &os, const Opm::TridiagonalMatrix<Scalar> &mat)
{
mat.print(os);
return os;
}
#endif

28
opm/core/utility/Unused.hpp
View File

@ -1,28 +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-2013 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 <http://www.gnu.org/licenses/>. *
*****************************************************************************/
#ifndef OPM_CORE_UNUSED_HH
#define OPM_CORE_UNUSED_HH
#ifndef HAS_ATTRIBUTE_UNUSED
#define OPM_UNUSED
#else
#define OPM_UNUSED __attribute__((unused))
#endif
#endif

196
tests/test_propertysystem.cpp
View File

@ -1,196 +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 <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
*
* \brief This file tests the properties system.
*
* We define a few type tags and property tags, then we attach values
* to (TypeTag, PropertyTag) tuples and finally we use them in the
* main function and print some diagnostic messages.
*/
#include "config.h"
#include <opm/core/utility/PropertySystem.hpp>
#include <iostream>
namespace Opm {
namespace Properties {
///////////////////
// Define some hierarchy of type tags:
//
// Vehicle -- CompactCar -- Sedan -_
// \ \.
// \ +- Pickup ---_
// \ / \.
// +-- Truck ---------------^ \.
// \ \.
// +- Tank ----------------------------------+- HummerH1
///////////////////
NEW_TYPE_TAG(Vehicle);
NEW_TYPE_TAG(CompactCar, INHERITS_FROM(Vehicle));
NEW_TYPE_TAG(Truck, INHERITS_FROM(Vehicle));
NEW_TYPE_TAG(Tank, INHERITS_FROM(Vehicle));
NEW_TYPE_TAG(Sedan, INHERITS_FROM(CompactCar));
NEW_TYPE_TAG(Pickup, INHERITS_FROM(Sedan, Truck));
NEW_TYPE_TAG(HummerH1, INHERITS_FROM(Sedan, Pickup, Tank));
///////////////////
// Define the property tags:
// TopSpeed, NumSeats, CanonCaliber, GasUsage, AutomaticTransmission, Payload
///////////////////
NEW_PROP_TAG(TopSpeed); // [km/h]
NEW_PROP_TAG(NumSeats); // []
NEW_PROP_TAG(CanonCaliber); // [mm]
NEW_PROP_TAG(GasUsage); // [l/100km]
NEW_PROP_TAG(AutomaticTransmission); // true/false
NEW_PROP_TAG(Payload); // [t]
///////////////////
// Make the AutomaticTransmission default to false
SET_BOOL_PROP(Vehicle, AutomaticTransmission, false);
///////////////////
// Define some values for the properties on the type tags:
//
// (CompactCar, TopSpeed) = GasUsage*35
// (CompactCar, NumSeats) = 5
// (CompactCar, GasUsage) = 4
//
// (Truck, TopSpeed) = 100
// (Truck, NumSeats) = 2
// (Truck, GasUsage) = 12
// (Truck, Payload) = 35
//
// (Tank, TopSpeed) = 60
// (Tank, GasUsage) = 65
// (Tank, CanonCaliber) = 120
//
// (Sedan, GasUsage) = 7
// (Sedan, AutomaticTransmission) = true
//
// (Pickup, TopSpeed) = 120
// (Pickup, Payload) = 5
//
// (HummmerH1, TopSpeed) = (Pickup, TopSpeed)
///////////////////
SET_INT_PROP(CompactCar, TopSpeed, GET_PROP_VALUE(TypeTag, GasUsage) * 30);
SET_INT_PROP(CompactCar, NumSeats, 5);
SET_INT_PROP(CompactCar, GasUsage, 4);
SET_INT_PROP(Truck, TopSpeed, 100);
SET_INT_PROP(Truck, NumSeats, 2);
SET_INT_PROP(Truck, GasUsage, 12);
SET_INT_PROP(Truck, Payload, 35);
SET_INT_PROP(Tank, TopSpeed, 60);
SET_INT_PROP(Tank, GasUsage, 65);
SET_INT_PROP(Tank, CanonCaliber, 120);
SET_INT_PROP(Sedan, GasUsage, 7);
SET_BOOL_PROP(Sedan, AutomaticTransmission, true);
SET_INT_PROP(Pickup, TopSpeed, 120);
SET_INT_PROP(Pickup, Payload, 5);
SET_INT_PROP(HummerH1, TopSpeed, GET_PROP_VALUE(TTAG(Pickup), TopSpeed));
///////////////////
// Unmount the canon from the Hummer
UNSET_PROP(HummerH1, CanonCaliber);
} // namespace Properties
} // namespace Opm
int main()
{
// print all properties for all type tags
std::cout << "---------------------------------------\n";
std::cout << "-- Property values\n";
std::cout << "---------------------------------------\n";
std::cout << "---------- Values for CompactCar ----------\n";
std::cout << "(CompactCar, TopSpeed) = " << GET_PROP_VALUE(TTAG(CompactCar), TopSpeed) << "\n";
std::cout << "(CompactCar, NumSeats) = " << GET_PROP_VALUE(TTAG(CompactCar), NumSeats) << "\n";
std::cout << "(CompactCar, GasUsage) = " << GET_PROP_VALUE(TTAG(CompactCar), GasUsage) << "\n";
std::cout << "(CompactCar, AutomaticTransmission) = " << GET_PROP_VALUE(TTAG(CompactCar), AutomaticTransmission) << "\n";
std::cout << "---------- Values for Truck ----------\n";
std::cout << "(Truck, TopSpeed) = " << GET_PROP_VALUE(TTAG(Truck), TopSpeed) << "\n";
std::cout << "(Truck, NumSeats) = " << GET_PROP_VALUE(TTAG(Truck), NumSeats) << "\n";
std::cout << "(Truck, GasUsage) = " << GET_PROP_VALUE(TTAG(Truck), GasUsage) << "\n";
std::cout << "(Truck, Payload) = " << GET_PROP_VALUE(TTAG(Truck), Payload) << "\n";
std::cout << "(Truck, AutomaticTransmission) = " << GET_PROP_VALUE(TTAG(Truck), AutomaticTransmission) << "\n";
std::cout << "---------- Values for Tank ----------\n";
std::cout << "(Tank, TopSpeed) = " << GET_PROP_VALUE(TTAG(Tank), TopSpeed) << "\n";
std::cout << "(Tank, GasUsage) = " << GET_PROP_VALUE(TTAG(Tank), GasUsage) << "\n";
std::cout << "(Tank, AutomaticTransmission) = " << GET_PROP_VALUE(TTAG(Tank), AutomaticTransmission) << "\n";
std::cout << "(Tank, CanonCaliber) = " << GET_PROP_VALUE(TTAG(Tank), CanonCaliber) << "\n";
std::cout << "---------- Values for Sedan ----------\n";
std::cout << "(Sedan, TopSpeed) = " << GET_PROP_VALUE(TTAG(Sedan), TopSpeed) << "\n";
std::cout << "(Sedan, NumSeats) = " << GET_PROP_VALUE(TTAG(Sedan), NumSeats) << "\n";
std::cout << "(Sedan, GasUsage) = " << GET_PROP_VALUE(TTAG(Sedan), GasUsage) << "\n";
std::cout << "(Sedan, AutomaticTransmission) = " << GET_PROP_VALUE(TTAG(Sedan), AutomaticTransmission) << "\n";
std::cout << "---------- Values for Pickup ----------\n";
std::cout << "(Pickup, TopSpeed) = " << GET_PROP_VALUE(TTAG(Pickup), TopSpeed) << "\n";
std::cout << "(Pickup, NumSeats) = " << GET_PROP_VALUE(TTAG(Pickup), NumSeats) << "\n";
std::cout << "(Pickup, GasUsage) = " << GET_PROP_VALUE(TTAG(Pickup), GasUsage) << "\n";
std::cout << "(Pickup, Payload) = " << GET_PROP_VALUE(TTAG(Pickup), Payload) << "\n";
std::cout << "(Pickup, AutomaticTransmission) = " << GET_PROP_VALUE(TTAG(Pickup), AutomaticTransmission) << "\n";
std::cout << "---------- Values for HummerH1 ----------\n";
std::cout << "(HummerH1, TopSpeed) = " << GET_PROP_VALUE(TTAG(HummerH1), TopSpeed) << "\n";
std::cout << "(HummerH1, NumSeats) = " << GET_PROP_VALUE(TTAG(HummerH1), NumSeats) << "\n";
std::cout << "(HummerH1, GasUsage) = " << GET_PROP_VALUE(TTAG(HummerH1), GasUsage) << "\n";
std::cout << "(HummerH1, Payload) = " << GET_PROP_VALUE(TTAG(HummerH1), Payload) << "\n";
std::cout << "(HummerH1, AutomaticTransmission) = " << GET_PROP_VALUE(TTAG(HummerH1), AutomaticTransmission) << "\n";
// CanonCaliber is explcitly unset for the Hummer -> this would not compile:
// std::cout << "(HummerH1, CanonCaliber) = " << GET_PROP_VALUE(TTAG(HummerH1), CanonCaliber) << "\n";
std::cout << "\n";
std::cout << "---------------------------------------\n";
std::cout << "-- Diagnostic messages\n";
std::cout << "---------------------------------------\n";
std::cout << "---- All properties for Sedan ---\n";
Opm::Properties::printValues<TTAG(Sedan)>();
std::cout << "---- Message for (HummerH1, CanonCaliber) ---\n"
<< PROP_DIAGNOSTIC(TTAG(HummerH1), CanonCaliber);
std::cout << "---- Message for (HummerH1, GasUsage) ---\n"
<< PROP_DIAGNOSTIC(TTAG(HummerH1), GasUsage);
std::cout << "---- Message for (HummerH1, AutomaticTransmission) ---\n"
<< PROP_DIAGNOSTIC(TTAG(HummerH1), AutomaticTransmission);
return 0;
}

346
tests/test_spline.cpp
View File

@ -1,346 +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 <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
*
* \brief This is a program to test the polynomial spline interpolation.
*
* It just prints some function to stdout. You can look at the result
* using the following commands:
*
----------- snip -----------
./test_spline > 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 "Monotonical"
----------- snap -----------
*/
#include "config.h"
#include <array>
#include <opm/core/utility/Spline.hpp>
#define GCC_VERSION (__GNUC__ * 10000 \
+ __GNUC_MINOR__ * 100 \
+ __GNUC_PATCHLEVEL__)
template <class Spline>
void testCommon(const Spline &sp,
const double *x,
const double *y)
{
static double eps = 1e-10;
static double epsFD = 1e-7;
int n = sp.numSamples();
for (int i = 0; i < n; ++i) {
// sure that we hit all sampling points
double y0 = (i>0)?sp.eval(x[i]-eps):y[0];
double y1 = sp.eval(x[i]);
double y2 = (i<n-1)?sp.eval(x[i]+eps):y[n-1];
if (std::abs(y0 - y[i]) > 100*eps || std::abs(y2 - y[i]) > 100*eps)
OPM_THROW(std::runtime_error,
"Spline seems to be discontinuous at sampling point " << i << "!");
if (std::abs(y1 - y[i]) > eps)
OPM_THROW(std::runtime_error,
"Spline does not capture sampling point " << i << "!");
// make sure the derivative is continuous (assuming that the
// second derivative is smaller than 1000)
double d1 = sp.evalDerivative(x[i]);
double d0 = (i>0)?sp.evalDerivative(x[i]-eps):d1;
double d2 = (i<n-1)?sp.evalDerivative(x[i]+eps):d1;
if (std::abs(d1 - d0) > 1000*eps || std::abs(d2 - d0) > 1000*eps)
OPM_THROW(std::runtime_error,
"Spline seems to exhibit a discontinuous derivative at sampling point " << i << "!");
}
// make sure the derivatives are consistent with the curve
int np = 3*n;
for (int i = 0; i < np; ++i) {
double xval = sp.xMin() + (sp.xMax() - sp.xMin())*i/np;
// first derivative
double y1 = sp.eval(xval+epsFD);
double y0 = sp.eval(xval);
double mFD = (y1 - y0)/epsFD;
double m = sp.evalDerivative(xval);
if (std::abs( mFD - m ) > 1000*epsFD)
OPM_THROW(std::runtime_error,
"Derivative of spline seems to be inconsistent with cuve"
" (" << mFD << " - " << m << " = " << mFD - m << ")!");
// second derivative
y1 = sp.evalDerivative(xval+epsFD);
y0 = sp.evalDerivative(xval);
mFD = (y1 - y0)/epsFD;
m = sp.evalSecondDerivative(xval);
if (std::abs( mFD - m ) > 1000*epsFD)
OPM_THROW(std::runtime_error,
"Second derivative of spline seems to be inconsistent with cuve"
" (" << mFD << " - " << m << " = " << mFD - m << ")!");
// Third derivative
y1 = sp.evalSecondDerivative(xval+epsFD);
y0 = sp.evalSecondDerivative(xval);
mFD = (y1 - y0)/epsFD;
m = sp.evalThirdDerivative(xval);
if (std::abs( mFD - m ) > 1000*epsFD)
OPM_THROW(std::runtime_error,
"Third derivative of spline seems to be inconsistent with cuve"
" (" << mFD << " - " << m << " = " << mFD - m << ")!");
}
}
template <class Spline>
void testFull(const Spline &sp,
const double *x,
const double *y,
double m0,
double m1)
{
// test the common properties of splines
testCommon(sp, x, y);
static double eps = 1e-5;
int n = sp.numSamples();
// make sure the derivative at both end points is correct
double d0 = sp.evalDerivative(x[0]);
double d1 = sp.evalDerivative(x[n-1]);
if (std::abs(d0 - m0) > eps)
OPM_THROW(std::runtime_error,
"Invalid derivative at beginning of interval: is "
<< d0 << " ought to be " << m0);
if (std::abs(d1 - m1) > eps)
OPM_THROW(std::runtime_error,
"Invalid derivative at end of interval: is "
<< d1 << " ought to be " << m1);
}
template <class Spline>
void testNatural(const Spline &sp,
const double *x,
const double *y)
{
// test the common properties of splines
testCommon(sp, x, y);
static double eps = 1e-5;
int n = sp.numSamples();
// make sure the second derivatives at both end points are 0
double d0 = sp.evalDerivative(x[0]);
double d1 = sp.evalDerivative(x[0] + eps);
double d2 = sp.evalDerivative(x[n-1] - eps);
double d3 = sp.evalDerivative(x[n-1]);
if (std::abs(d1 - d0)/eps > 1000*eps)
OPM_THROW(std::runtime_error,
"Invalid second derivative at beginning of interval: is "
<< (d1 - d0)/eps << " ought to be 0");
if (std::abs(d3 - d2)/eps > 1000*eps)
OPM_THROW(std::runtime_error,
"Invalid second derivative at end of interval: is "
<< (d3 - d2)/eps << " ought to be 0");
}
template <class Spline>
void testMonotonic(const Spline &sp,
const double *x,
const double *y)
{
// test the common properties of splines
testCommon(sp, x, y);
int n = sp.numSamples();
for (int i = 0; i < n - 1; ++ i) {
// make sure that the spline is monotonic for each interval
// between sampling points
if (!sp.monotonic(x[i], x[i + 1]))
OPM_THROW(std::runtime_error,
"Spline says it is not monotonic in interval "
<< i << " where it should be");
// test the intersection methods
double d = (y[i] + y[i+1])/2;
double interX = sp.intersectInterval(x[i], x[i+1],
/*a=*/0, /*b=*/0, /*c=*/0, d);
double interY = sp.eval(interX);
if (std::abs(interY - d) > 1e-5)
OPM_THROW(std::runtime_error,
"Spline::intersectInterval() seems to be broken: "
<< sp.eval(interX) << " - " << d << " = " << sp.eval(interX) - d << "!");
}
// make sure the spline says to be monotonic on the (extrapolated)
// left and right sides
if (!sp.monotonic(x[0] - 1.0, (x[0] + x[1])/2, /*extrapolate=*/true))
OPM_THROW(std::runtime_error,
"Spline says it is not monotonic on left side where it should be");
if (!sp.monotonic((x[n - 2]+ x[n - 1])/2, x[n-1] + 1.0, /*extrapolate=*/true))
OPM_THROW(std::runtime_error,
"Spline says it is not monotonic on right side where it should be");
for (int i = 0; i < n - 2; ++ i) {
// make sure that the spline says that it is non-monotonic for
// if extrema are within the queried interval
if (sp.monotonic((x[i] + x[i + 1])/2, (x[i + 1] + x[i + 2])/2))
OPM_THROW(std::runtime_error,
"Spline says it is monotonic in interval "
<< i << " where it should not be");
}
}
void testAll()
{
double x[] = { 0, 5, 7.5, 8.75, 9.375 };
double y[] = { 10, 0, 10, 0, 10 };
double m0 = 10;
double m1 = -10;
double points[][2] =
{
{x[0], y[0]},
{x[1], y[1]},
{x[2], y[2]},
{x[3], y[3]},
{x[4], y[4]},
};
#if GCC_VERSION >= 40500
std::initializer_list<const std::pair<double, double> > pointsInitList =
{
{x[0], y[0]},
{x[1], y[1]},
{x[2], y[2]},
{x[3], y[3]},
{x[4], y[4]},
};
#endif
std::vector<double> xVec;
std::vector<double> yVec;
std::vector<double*> pointVec;
for (int i = 0; i < 5; ++i) {
xVec.push_back(x[i]);
yVec.push_back(y[i]);
pointVec.push_back(points[i]);
}
/////////
// test spline with two sampling points
/////////
// full spline
{ Opm::Spline<double> sp(x[0], x[1], y[0], y[1], m0, m1); sp.set(x[0],x[1],y[0],y[1],m0, m1); testFull(sp, x, y, m0, m1); };
{ Opm::Spline<double> sp(2, x, y, m0, m1); sp.setXYArrays(2, x, y, m0, m1); testFull(sp, x, y, m0, m1); };
{ Opm::Spline<double> sp(2, points, m0, m1); sp.setArrayOfPoints(2, points, m0, m1); testFull(sp, x, y, m0, m1); };
/////////
// test variable length splines
/////////
// full spline
{ Opm::Spline<double> sp(5, x, y, m0, m1); sp.setXYArrays(5,x,y,m0, m1); testFull(sp, x, y, m0, m1); };
{ Opm::Spline<double> sp(xVec, yVec, m0, m1); sp.setXYContainers(xVec,yVec,m0, m1); testFull(sp, x, y, m0, m1); };
{ Opm::Spline<double> sp; sp.setArrayOfPoints(5,points,m0, m1); testFull(sp, x, y, m0, m1); };
{ Opm::Spline<double> sp; sp.setContainerOfPoints(pointVec,m0, m1); testFull(sp, x, y, m0, m1); };
#if GCC_VERSION >= 40500
{ Opm::Spline<double> sp; sp.setContainerOfTuples(pointsInitList,m0, m1); testFull(sp, x, y, m0, m1); };
#endif
// natural spline
{ Opm::Spline<double> sp(5, x, y); sp.setXYArrays(5,x,y); testNatural(sp, x, y); };
{ Opm::Spline<double> sp(xVec, yVec); sp.setXYContainers(xVec,yVec); testNatural(sp, x, y); };
{ Opm::Spline<double> sp; sp.setArrayOfPoints(5,points); testNatural(sp, x, y); };
{ Opm::Spline<double> sp; sp.setContainerOfPoints(pointVec); testNatural(sp, x, y); };
#if GCC_VERSION >= 40500
{ Opm::Spline<double> sp; sp.setContainerOfTuples(pointsInitList); testNatural(sp, x, y); };
#endif
}
void plot()
{
const int numSamples = 5;
const int n = numSamples - 1;
typedef std::array<double, numSamples> FV;
double x_[] = { 0, 5, 7.5, 8.75, 10 };
double y_[] = { 10, 0, 10, 0, 10 };
double m1 = 10;
double m2 = -10;
FV &xs = *reinterpret_cast<FV*>(x_);
FV &ys = *reinterpret_cast<FV*>(y_);
Opm::Spline<double> spFull(xs, ys, m1, m2);
Opm::Spline<double> spNatural(xs, ys);
Opm::Spline<double> spPeriodic(xs, ys, /*type=*/Opm::Spline<double>::Periodic);
Opm::Spline<double> spMonotonic(xs, ys, /*type=*/Opm::Spline<double>::Monotonic);
testMonotonic(spMonotonic, x_, y_);
spFull.printCSV(x_[0] - 1.00001,
x_[n] + 1.00001,
1000);
std::cout << "\n";
spNatural.printCSV(x_[0] - 1.00001,
x_[n] + 1.00001,
1000);
std::cout << "\n";
spPeriodic.printCSV(x_[0] - 1.00001,
x_[n] + 1.00001,
1000);
std::cout << "\n";
spMonotonic.printCSV(x_[0] - 1.00001,
x_[n] + 1.00001,
1000);
std::cout << "\n";
}
int main()
{
try {
testAll();
plot();
}
catch (const std::exception &e) {
std::cout << "Caught OPM exception: " << e.what() << "\n";
}
return 0;
}