make the tests compile using OPM's shiny new CMake based build system
This commit is contained in:
Andreas Lauser
parent
c61b0bd63d
commit
1fa4ad0c5d
162
opm/common/dynamictabulated2dfunction.hh
Normal file
162
opm/common/dynamictabulated2dfunction.hh
Normal file
@ -0,0 +1,162 @@
|
||||
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
||||
// vi: set et ts=4 sw=4 sts=4:
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 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
|
||||
*
|
||||
* \copydoc Opm::DynamicTabulated2DFunction
|
||||
*/
|
||||
#ifndef OPM_DYNAMIC_TABULATED_2D_FUNCTION_HH
|
||||
#define OPM_DYNAMIC_TABULATED_2D_FUNCTION_HH
|
||||
|
||||
#include <opm/common/exceptions.hh>
|
||||
#include <opm/common/tabulated2dfunction.hh>
|
||||
|
||||
#include <vector>
|
||||
|
||||
#include <assert.h>
|
||||
|
||||
namespace Opm {
|
||||
|
||||
/*!
|
||||
* \copydoc Opm::Tabulated2DFunction
|
||||
*
|
||||
* This class can be used when the sampling points are calculated at
|
||||
* run time.
|
||||
*/
|
||||
template <class Scalar>
|
||||
class DynamicTabulated2DFunction
|
||||
: public Tabulated2DFunction<Scalar, DynamicTabulated2DFunction<Scalar> >
|
||||
{
|
||||
public:
|
||||
DynamicTabulated2DFunction()
|
||||
{ }
|
||||
|
||||
/*!
|
||||
* \brief Constructor where the tabulation parameters are already
|
||||
* provided.
|
||||
*/
|
||||
DynamicTabulated2DFunction(Scalar xMin, Scalar xMax, int m,
|
||||
Scalar yMin, Scalar yMax, int n)
|
||||
{
|
||||
resize(xMin, xMax, m, yMin, yMax, n);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Resize the tabulation to a new range.
|
||||
*/
|
||||
void resize(Scalar xMin, Scalar xMax, int m,
|
||||
Scalar yMin, Scalar yMax, int n)
|
||||
{
|
||||
samples_.resize(m*n);
|
||||
|
||||
m_ = m;
|
||||
n_ = n;
|
||||
|
||||
xMin_ = xMin;
|
||||
xMax_ = xMax;
|
||||
|
||||
yMin_ = yMin;
|
||||
yMax_ = yMax;
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns the minimum of the X coordinate of the sampling points.
|
||||
*/
|
||||
Scalar xMin() const
|
||||
{ return xMin_; }
|
||||
|
||||
/*!
|
||||
* \brief Returns the maximum of the X coordinate of the sampling points.
|
||||
*/
|
||||
Scalar xMax() const
|
||||
{ return xMax_; }
|
||||
|
||||
/*!
|
||||
* \brief Returns the minimum of the Y coordinate of the sampling points.
|
||||
*/
|
||||
Scalar yMin() const
|
||||
{ return yMin_; }
|
||||
|
||||
/*!
|
||||
* \brief Returns the maximum of the Y coordinate of the sampling points.
|
||||
*/
|
||||
Scalar yMax() const
|
||||
{ return yMax_; }
|
||||
|
||||
/*!
|
||||
* \brief Returns the number of sampling points in X direction.
|
||||
*/
|
||||
int numX() const
|
||||
{ return m_; }
|
||||
|
||||
/*!
|
||||
* \brief Returns the number of sampling points in Y direction.
|
||||
*/
|
||||
int numY() const
|
||||
{ return n_; }
|
||||
|
||||
/*!
|
||||
* \brief Get the value of the sample point which is at the
|
||||
* intersection of the \f$i\f$-th interval of the x-Axis
|
||||
* and the \f$j\f$-th of the y-Axis.
|
||||
*/
|
||||
Scalar getSamplePoint(int i, int j) const
|
||||
{
|
||||
assert(0 <= i && i < m_);
|
||||
assert(0 <= j && j < n_);
|
||||
|
||||
return samples_[j*m_ + i];
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Set the value of the sample point which is at the
|
||||
* intersection of the \f$i\f$-th interval of the x-Axis
|
||||
* and the \f$j\f$-th of the y-Axis.
|
||||
*/
|
||||
void setSamplePoint(int i, int j, Scalar value)
|
||||
{
|
||||
assert(0 <= i && i < m_);
|
||||
assert(0 <= j && j < n_);
|
||||
|
||||
samples_[j*m_ + i] = value;
|
||||
}
|
||||
|
||||
private:
|
||||
// the vector which contains the values of the sample points
|
||||
// f(x_i, y_j). don't use this directly, use getSamplePoint(i,j)
|
||||
// instead!
|
||||
std::vector<Scalar> samples_;
|
||||
|
||||
// the number of sample points in x direction
|
||||
int m_;
|
||||
|
||||
// the number of sample points in y direction
|
||||
int n_;
|
||||
|
||||
// the range of the tabulation on the x axis
|
||||
Scalar xMin_;
|
||||
Scalar xMax_;
|
||||
|
||||
// the range of the tabulation on the y axis
|
||||
Scalar yMin_;
|
||||
Scalar yMax_;
|
||||
};
|
||||
}
|
||||
|
||||
#endif
|
||||
77
opm/common/exceptions.hh
Normal file
77
opm/common/exceptions.hh
Normal file
@ -0,0 +1,77 @@
|
||||
// -*- 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 <http://www.gnu.org/licenses/>. *
|
||||
*****************************************************************************/
|
||||
/*!
|
||||
* \file
|
||||
* \brief Some exceptions thrown by opm-material.
|
||||
*/
|
||||
#ifndef OPM_EXCEPTIONS_HH
|
||||
#define OPM_EXCEPTIONS_HH
|
||||
|
||||
#include <dune/common/exceptions.hh>
|
||||
|
||||
#include <string>
|
||||
|
||||
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
|
||||
412
opm/common/fixedlengthspline_.hh
Normal file
412
opm/common/fixedlengthspline_.hh
Normal file
@ -0,0 +1,412 @@
|
||||
// -*- 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
|
||||
* \copydoc Opm::FixedLengthSpline_
|
||||
*/
|
||||
#ifndef OPM_FIXED_LENGTH_SPLINE_HH
|
||||
#define OPM_FIXED_LENGTH_SPLINE_HH
|
||||
|
||||
#include <dune/common/fvector.hh>
|
||||
#include <dune/common/fmatrix.hh>
|
||||
#include <dune/istl/btdmatrix.hh>
|
||||
|
||||
#include "splinecommon_.hh"
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
/*!
|
||||
* \brief The common code for all 3rd order polynomial splines with
|
||||
* more than two sampling points.
|
||||
*/
|
||||
template<class ScalarT, int nSamples>
|
||||
class FixedLengthSpline_
|
||||
: public SplineCommon_<ScalarT,
|
||||
FixedLengthSpline_<ScalarT, nSamples> >
|
||||
{
|
||||
friend class SplineCommon_<ScalarT, FixedLengthSpline_<ScalarT, nSamples> >;
|
||||
|
||||
typedef ScalarT Scalar;
|
||||
typedef Dune::FieldVector<Scalar, nSamples> Vector;
|
||||
typedef Dune::BlockVector<Dune::FieldVector<Scalar, 1> > BlockVector;
|
||||
typedef Dune::BTDMatrix<Dune::FieldMatrix<Scalar, 1, 1> > 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 <class ScalarArrayX, class ScalarArrayY>
|
||||
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 <class ScalarContainerX, class ScalarContainerY>
|
||||
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 <class PointArray>
|
||||
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 <class XYContainer>
|
||||
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 <class XYContainer>
|
||||
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 <class ScalarArrayX, class ScalarArrayY>
|
||||
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 <class ScalarContainerX, class ScalarContainerY>
|
||||
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 <class PointArray>
|
||||
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 <class XYContainer>
|
||||
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 <class XYContainer>
|
||||
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
|
||||
392
opm/common/math.hh
Normal file
392
opm/common/math.hh
Normal file
@ -0,0 +1,392 @@
|
||||
// -*- 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 Define some often used mathematical functions
|
||||
*/
|
||||
#ifndef OPM_MATH_HH
|
||||
#define OPM_MATH_HH
|
||||
|
||||
#include <dune/common/fvector.hh>
|
||||
#include <dune/common/fmatrix.hh>
|
||||
|
||||
#include <cmath>
|
||||
#include <algorithm>
|
||||
|
||||
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 <class Scalar>
|
||||
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 <class Scalar>
|
||||
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 <class Scalar, int m, int n>
|
||||
void harmonicMeanMatrix(Dune::FieldMatrix<Scalar, m, n> &K,
|
||||
const Dune::FieldMatrix<Scalar, m, n> &Ki,
|
||||
const Dune::FieldMatrix<Scalar, m, n> &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 <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;
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Cross product of two vectors in three-dimensional Euclidean space
|
||||
*
|
||||
* \param vec1 The first vector
|
||||
* \param vec2 The second vector
|
||||
*/
|
||||
template <class Scalar>
|
||||
Dune::FieldVector<Scalar, 3> crossProduct(const Dune::FieldVector<Scalar, 3> &vec1,
|
||||
const Dune::FieldVector<Scalar, 3> &vec2)
|
||||
{
|
||||
Dune::FieldVector<Scalar, 3> 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
|
||||
532
opm/common/spline.hh
Normal file
532
opm/common/spline.hh
Normal file
@ -0,0 +1,532 @@
|
||||
// -*- 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
|
||||
* \copydoc Opm::Spline
|
||||
*/
|
||||
#ifndef OPM_SPLINE_HH
|
||||
#define OPM_SPLINE_HH
|
||||
|
||||
#include "fixedlengthspline_.hh"
|
||||
#include "variablelengthspline_.hh"
|
||||
#include "splinecommon_.hh"
|
||||
|
||||
#include <dune/common/fvector.hh>
|
||||
#include <dune/common/fmatrix.hh>
|
||||
|
||||
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 Scalar, int numSamples = 2>
|
||||
class Spline : public FixedLengthSpline_<Scalar, numSamples>
|
||||
{
|
||||
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 <class ScalarArray>
|
||||
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 <class PointArray>
|
||||
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 <class ScalarArray>
|
||||
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 <class PointArray>
|
||||
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 Scalar>
|
||||
class Spline<Scalar, /*numSamples=*/-1> : public VariableLengthSpline_<Scalar>
|
||||
{
|
||||
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 <class ScalarArrayX, class ScalarArrayY>
|
||||
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 <class PointArray>
|
||||
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 <class ScalarContainer>
|
||||
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 <class PointContainer>
|
||||
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 <class ScalarArray>
|
||||
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 <class PointArray>
|
||||
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 <class ScalarContainerX, class ScalarContainerY>
|
||||
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 <class PointContainer>
|
||||
Spline(const PointContainer &points,
|
||||
Scalar m0,
|
||||
Scalar m1)
|
||||
{ this->setContainerOfPoints(points, m0, m1); }
|
||||
};
|
||||
|
||||
/*!
|
||||
* \brief Do not allow splines with zero sampling points
|
||||
*/
|
||||
template<class Scalar>
|
||||
class Spline<Scalar, /*numSamples=*/0>
|
||||
// Splines with zero sampling points do not make sense!
|
||||
{ private: Spline() { } };
|
||||
|
||||
/*!
|
||||
* \brief Do not allow splines with one sampling point
|
||||
*/
|
||||
template<class Scalar>
|
||||
class Spline<Scalar, /*numSamples=*/1>
|
||||
// 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 Scalar>
|
||||
class Spline<Scalar, 2> : public SplineCommon_<Scalar, Spline<Scalar, 2> >
|
||||
{
|
||||
friend class SplineCommon_<Scalar, Spline<Scalar, 2> >;
|
||||
typedef Dune::FieldVector<Scalar, 2> Vector;
|
||||
typedef Dune::FieldMatrix<Scalar, 2, 2> 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 <class ScalarArrayX, class ScalarArrayY>
|
||||
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 <class PointArray>
|
||||
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 <class ScalarContainer>
|
||||
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 <class ScalarContainerX, class ScalarContainerY>
|
||||
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 <class PointArray>
|
||||
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 <class PointContainer>
|
||||
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 <class TupleContainer>
|
||||
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
|
||||
685
opm/common/splinecommon_.hh
Normal file
685
opm/common/splinecommon_.hh
Normal file
@ -0,0 +1,685 @@
|
||||
// -*- 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
|
||||
* \copydoc Opm::SplineCommon_
|
||||
*/
|
||||
#ifndef OPM_SPLINE_COMMON__HH
|
||||
#define OPM_SPLINE_COMMON__HH
|
||||
|
||||
#include "valgrind.hh"
|
||||
#include "math.hh"
|
||||
|
||||
#include <dune/common/exceptions.hh>
|
||||
|
||||
#include <tuple>
|
||||
#include <algorithm>
|
||||
#include <iostream>
|
||||
#include <cassert>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
/*!
|
||||
* \brief The common code for all 3rd order polynomial splines.
|
||||
*/
|
||||
template<class ScalarT, class ImplementationT>
|
||||
class SplineCommon_
|
||||
{
|
||||
typedef ScalarT Scalar;
|
||||
typedef ImplementationT Implementation;
|
||||
|
||||
Implementation &asImp_()
|
||||
{ return *static_cast<Implementation*>(this); }
|
||||
const Implementation &asImp_() const
|
||||
{ return *static_cast<const Implementation*>(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<Scalar>(x_(0), x), std::min<Scalar>(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"); //<<a<<"x^3 + "<<b<"x^2 + "<<c<"x + "<<d);
|
||||
}
|
||||
}
|
||||
|
||||
if (nSol != 1)
|
||||
DUNE_THROW(Dune::MathError,
|
||||
"Spline has no intersection"); //<<a<"x^3 + " <<b<"x^2 + "<<c<"x + "<<d<<"!");
|
||||
|
||||
return tmpSol[0];
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns 1 if the spline is monotonically increasing, -1
|
||||
* if the spline is mononously decreasing and 0 if the
|
||||
* spline is not monotonous in the interval (x0, x1).
|
||||
*
|
||||
* In the corner case where the whole spline is flat, it returns
|
||||
* 2.
|
||||
*/
|
||||
int monotonic(Scalar x0, Scalar x1) const
|
||||
{
|
||||
assert(applies(x0));
|
||||
assert(applies(x1));
|
||||
assert(x0 != x1);
|
||||
|
||||
// make sure that x0 is smaller than x1
|
||||
if (x0 > 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 <class DestVector, class SourceVector>
|
||||
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 <class DestVector, class ListIterator>
|
||||
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 <class DestVector, class ListIterator>
|
||||
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 <class Vector, class Matrix>
|
||||
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 <class Vector, class Matrix>
|
||||
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 <class Vector, class Matrix>
|
||||
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
|
||||
108
opm/common/statictabulated2dfunction.hh
Normal file
108
opm/common/statictabulated2dfunction.hh
Normal file
@ -0,0 +1,108 @@
|
||||
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
||||
// vi: set et ts=4 sw=4 sts=4:
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 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
|
||||
*
|
||||
* \copydoc Opm::StaticTabulated2DFunction
|
||||
*/
|
||||
#ifndef OPM_STATIC_TABULATED_2D_FUNCTION_HH
|
||||
#define OPM_STATIC_TABULATED_2D_FUNCTION_HH
|
||||
|
||||
#include <opm/common/exceptions.hh>
|
||||
#include <opm/common/tabulated2dfunction.hh>
|
||||
|
||||
#include <assert.h>
|
||||
|
||||
namespace Opm {
|
||||
/*!
|
||||
* \copydoc Opm::Tabulated2DFunction
|
||||
*
|
||||
* This class can be used when the sampling points are calculated at
|
||||
* compile time.
|
||||
*/
|
||||
template <class Traits>
|
||||
class StaticTabulated2DFunction
|
||||
: public Tabulated2DFunction<typename Traits::Scalar, StaticTabulated2DFunction<Traits> >
|
||||
{
|
||||
typedef typename Traits::Scalar Scalar;
|
||||
|
||||
public:
|
||||
StaticTabulated2DFunction()
|
||||
{ }
|
||||
|
||||
/*!
|
||||
* \brief Returns the minimum of the X coordinate of the sampling points.
|
||||
*/
|
||||
Scalar xMin() const
|
||||
{ return Traits::xMin; }
|
||||
|
||||
/*!
|
||||
* \brief Returns the maximum of the X coordinate of the sampling points.
|
||||
*/
|
||||
Scalar xMax() const
|
||||
{ return Traits::xMax; }
|
||||
|
||||
/*!
|
||||
* \brief Returns the minimum of the Y coordinate of the sampling points.
|
||||
*/
|
||||
Scalar yMin() const
|
||||
{ return Traits::yMin; }
|
||||
|
||||
/*!
|
||||
* \brief Returns the maximum of the Y coordinate of the sampling points.
|
||||
*/
|
||||
Scalar yMax() const
|
||||
{ return Traits::yMax; }
|
||||
|
||||
/*!
|
||||
* \brief Returns the number of sampling points in X direction.
|
||||
*/
|
||||
int numX() const
|
||||
{ return Traits::numX; }
|
||||
|
||||
/*!
|
||||
* \brief Returns the number of sampling points in Y direction.
|
||||
*/
|
||||
int numY() const
|
||||
{ return Traits::numY; }
|
||||
|
||||
/*!
|
||||
* \brief Get the value of the sample point which is at the
|
||||
* intersection of the \f$i\f$-th interval of the x-Axis
|
||||
* and the \f$j\f$-th of the y-Axis.
|
||||
*/
|
||||
Scalar getSamplePoint(int i, int j) const
|
||||
{
|
||||
#if !defined NDEBUG
|
||||
if (i < 0 || i >= Traits::numX ||
|
||||
j < 0 || j >= Traits::numY) {
|
||||
DUNE_THROW(NumericalProblem,
|
||||
"Attempt to access element ("
|
||||
<< i << ", " << j
|
||||
<< ") on a " << Traits::name << " table of size ("
|
||||
<< Traits::numX << ", " << Traits::numY
|
||||
<< ")\n");
|
||||
};
|
||||
#endif
|
||||
return Traits::vals[i][j];
|
||||
}
|
||||
};
|
||||
}
|
||||
|
||||
#endif
|
||||
143
opm/common/tabulated2dfunction.hh
Normal file
143
opm/common/tabulated2dfunction.hh
Normal file
@ -0,0 +1,143 @@
|
||||
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
||||
// vi: set et ts=4 sw=4 sts=4:
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2011-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
|
||||
*
|
||||
* \copydoc Opm::Tabulated2DFunction
|
||||
*/
|
||||
#ifndef OPM_TABULATED_2D_FUNCTION_HH
|
||||
#define OPM_TABULATED_2D_FUNCTION_HH
|
||||
|
||||
#include <opm/common/exceptions.hh>
|
||||
|
||||
#include <assert.h>
|
||||
|
||||
namespace Opm {
|
||||
/*!
|
||||
* \brief A generic class that represents tabulated 2 dimensional functions
|
||||
*
|
||||
* This class can be used to tabulate a two dimensional function
|
||||
* \f$f(x, y)\f$ over the range \f$[x_{min}, x_{max}] \times [y_{min},
|
||||
* y_{max}]\f$. For this, the ranges of the \f$x\f$ and \f$y\f$ axes are
|
||||
* divided into \f$m\f$ and \f$n\f$ sub-intervals and the values of
|
||||
* \f$f(x_i, y_j)\f$ need to be provided. Here, \f$x_i\f$ and
|
||||
* \f$y_j\f$ are the largest positions of the \f$i\f$-th and
|
||||
* \f$j\f$-th intervall. Between these sampling points this tabulation
|
||||
* class uses linear interpolation.
|
||||
*
|
||||
* If the class is queried for a value outside of the tabulated range,
|
||||
* a \c Opm::NumericalProblem exception is thrown.
|
||||
*/
|
||||
template <class Scalar, class Implementation>
|
||||
class Tabulated2DFunction
|
||||
{
|
||||
public:
|
||||
Tabulated2DFunction()
|
||||
{ }
|
||||
|
||||
/*!
|
||||
* \brief Return the position on the x-axis of the i-th interval.
|
||||
*/
|
||||
Scalar iToX(int i) const
|
||||
{
|
||||
assert(0 <= i && i < asImp_().numX());
|
||||
|
||||
return asImp_().xMin() + i*(asImp_().xMax() - asImp_().xMin())/(asImp_().numX() - 1);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Return the position on the y-axis of the j-th interval.
|
||||
*/
|
||||
Scalar jToY(int j) const
|
||||
{
|
||||
assert(0 <= j && j < asImp_().numY());
|
||||
|
||||
return asImp_().yMin() + j*(asImp_().yMax() - asImp_().yMin())/(asImp_().numY() - 1);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Return the interval index of a given position on the x-axis.
|
||||
*
|
||||
* This method returns a *floating point* number. The integer part
|
||||
* should be interpreted as interval, the decimal places are the
|
||||
* position of the x value between the i-th and the (i+1)-th
|
||||
* sample point.
|
||||
*/
|
||||
Scalar xToI(Scalar x) const
|
||||
{ return (x - asImp_().xMin())/(asImp_().xMax() - asImp_().xMin())*(asImp_().numX() - 1); }
|
||||
|
||||
/*!
|
||||
* \brief Return the interval index of a given position on the y-axis.
|
||||
*
|
||||
* This method returns a *floating point* number. The integer part
|
||||
* should be interpreted as interval, the decimal places are the
|
||||
* position of the y value between the j-th and the (j+1)-th
|
||||
* sample point.
|
||||
*/
|
||||
Scalar yToJ(Scalar y) const
|
||||
{ return (y - asImp_().yMin())/(asImp_().yMax() - asImp_().yMin())*(asImp_().numY() - 1); }
|
||||
|
||||
/*!
|
||||
* \brief Returns true iff a coordinate lies in the tabulated range
|
||||
*/
|
||||
bool applies(Scalar x, Scalar y) const
|
||||
{ return asImp_().xMin() <= x && x <= asImp_().xMax() && asImp_().yMin() <= y && y <= asImp_().yMax(); }
|
||||
|
||||
/*!
|
||||
* \brief Evaluate the function at a given (x,y) position.
|
||||
*
|
||||
* If this method is called for a value outside of the tabulated
|
||||
* range, a \c Opm::NumericalProblem exception is thrown.
|
||||
*/
|
||||
Scalar eval(Scalar x, Scalar y) const
|
||||
{
|
||||
#ifndef NDEBUG
|
||||
if (!applies(x,y))
|
||||
{
|
||||
DUNE_THROW(NumericalProblem,
|
||||
"Attempt to get tabulated value for ("
|
||||
<< x << ", " << y
|
||||
<< ") on a table of extend "
|
||||
<< asImp_().xMin() << " to " << asImp_().xMax() << " times "
|
||||
<< asImp_().yMin() << " to " << asImp_().yMax());
|
||||
};
|
||||
#endif
|
||||
|
||||
Scalar alpha = xToI(x);
|
||||
Scalar beta = yToJ(y);
|
||||
|
||||
int i = std::max(0, std::min(asImp_().numX(), static_cast<int>(alpha)));
|
||||
int j = std::max(0, std::min(asImp_().numY(), static_cast<int>(beta)));
|
||||
|
||||
alpha -= i;
|
||||
beta -= j;
|
||||
|
||||
// bi-linear interpolation
|
||||
Scalar s1 = asImp_().getSamplePoint(i, j)*(1.0 - alpha) + asImp_().getSamplePoint(i + 1, j)*alpha;
|
||||
Scalar s2 = asImp_().getSamplePoint(i, j + 1)*(1.0 - alpha) + asImp_().getSamplePoint(i + 1, j + 1)*alpha;
|
||||
return s1*(1.0 - beta) + s2*beta;
|
||||
}
|
||||
|
||||
private:
|
||||
const Implementation &asImp_() const
|
||||
{ return *static_cast<const Implementation*>(this); }
|
||||
};
|
||||
}
|
||||
|
||||
#endif
|
||||
285
opm/common/valgrind.hh
Normal file
285
opm/common/valgrind.hh
Normal file
@ -0,0 +1,285 @@
|
||||
// -*- 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 Some templates to wrap the valgrind macros
|
||||
*/
|
||||
#ifndef OPM_VALGRIND_HH
|
||||
#define OPM_VALGRIND_HH
|
||||
|
||||
#if ! HAVE_VALGRIND && ! defined(DOXYGEN)
|
||||
namespace Valgrind
|
||||
{
|
||||
bool boolBlubb(bool value) { return value; }
|
||||
void voidBlubb() { }
|
||||
|
||||
#define SetUndefined(t) voidBlubb()
|
||||
#define SetDefined(t) voidBlubb()
|
||||
#define CheckDefined(t) boolBlubb(true)
|
||||
#define SetNoAccess(t) voidBlubb()
|
||||
#define IsRunning() boolBlubb(false)
|
||||
}
|
||||
|
||||
#else
|
||||
|
||||
#include <valgrind/memcheck.h>
|
||||
|
||||
namespace Valgrind
|
||||
{
|
||||
/*!
|
||||
* \ingroup Valgrind
|
||||
* \brief Returns whether the program is running under Valgrind or not.
|
||||
*/
|
||||
inline bool IsRunning()
|
||||
{
|
||||
#if !defined NDEBUG && HAVE_VALGRIND
|
||||
return RUNNING_ON_VALGRIND;
|
||||
#else
|
||||
return false;
|
||||
#endif
|
||||
}
|
||||
|
||||
|
||||
/*!
|
||||
* \ingroup Valgrind
|
||||
* \brief Make valgrind complain if any of the memory occupied by an object
|
||||
* is undefined.
|
||||
*
|
||||
* Please note that this does not check whether the destinations of an
|
||||
* object's pointers or references are defined. Also, for performance
|
||||
* reasons the compiler might insert "padding bytes" between within
|
||||
* the objects which leads to false positives.
|
||||
*
|
||||
* Example:
|
||||
*
|
||||
* \code
|
||||
* int i;
|
||||
* Opm::Valgrind::CheckDefined(i); // Valgrind complains!
|
||||
* \endcode
|
||||
*
|
||||
* \tparam T The type of the object which ought to be checked
|
||||
*
|
||||
* \param value the object which valgrind should check
|
||||
*
|
||||
* \return true iff there are no undefined bytes in the memory
|
||||
* occupied by the object.
|
||||
*/
|
||||
template <class T>
|
||||
inline bool CheckDefined(const T &value)
|
||||
{
|
||||
#if !defined NDEBUG && HAVE_VALGRIND
|
||||
unsigned int tmp = VALGRIND_CHECK_MEM_IS_DEFINED(&value, sizeof(T));
|
||||
return tmp == 0;
|
||||
#else
|
||||
return true;
|
||||
#endif
|
||||
}
|
||||
/*!
|
||||
* \ingroup Valgrind
|
||||
*
|
||||
* * \brief Make valgrind complain if any of the the memory occupied
|
||||
* by a C-style array objects is undefined.
|
||||
*
|
||||
* Please note that this does not check whether the destinations of an
|
||||
* object's pointers or references are defined. Also, for performance
|
||||
* reasons the compiler might insert "padding bytes" between within
|
||||
* the objects which leads to false positives.
|
||||
*
|
||||
* Example:
|
||||
*
|
||||
* \code
|
||||
* int i[2];
|
||||
* Opm::Valgrind::CheckDefined(i, 2); // Valgrind complains!
|
||||
* \endcode
|
||||
*
|
||||
* \tparam T The type of the object which ought to be checked
|
||||
*
|
||||
* \param value Pointer to the first object of the array.
|
||||
* \param size The size of the array in number of objects
|
||||
*
|
||||
* \return true iff there are no undefined bytes in the memory
|
||||
* occupied by the array.
|
||||
*/
|
||||
template <class T>
|
||||
inline bool CheckDefined(const T *value, int size)
|
||||
{
|
||||
#if !defined NDEBUG && HAVE_VALGRIND
|
||||
unsigned int tmp = VALGRIND_CHECK_MEM_IS_DEFINED(value, size*sizeof(T));
|
||||
return tmp == 0;
|
||||
#else
|
||||
return true;
|
||||
#endif
|
||||
}
|
||||
|
||||
/*!
|
||||
* \ingroup Valgrind
|
||||
* \brief Make the memory on which an object resides undefined in
|
||||
* valgrind runs.
|
||||
*
|
||||
* Example:
|
||||
*
|
||||
* \code
|
||||
* int i = 0;
|
||||
* Opm::Valgrind::SetUndefined(i);
|
||||
* Opm::Valgrind::CheckDefined(i); // Valgrind complains!
|
||||
* \endcode
|
||||
*
|
||||
* \tparam T The type of the object which ought to be set to undefined
|
||||
*
|
||||
* \param value The object which's memory valgrind should be told is undefined
|
||||
*/
|
||||
template <class T>
|
||||
inline void SetUndefined(const T &value)
|
||||
{
|
||||
#if !defined NDEBUG && HAVE_VALGRIND
|
||||
auto DUNE_UNUSED result = VALGRIND_MAKE_MEM_UNDEFINED(&value, sizeof(T));
|
||||
#endif
|
||||
}
|
||||
|
||||
/*!
|
||||
* \ingroup Valgrind
|
||||
* \brief Make the memory on which an array of object resides
|
||||
* undefined in valgrind runs.
|
||||
*
|
||||
* Example:
|
||||
*
|
||||
* \code
|
||||
* int i[3] = {0, 1, 3};
|
||||
* Opm::Valgrind::SetUndefined(&i[1], 2);
|
||||
* Opm::Valgrind::CheckDefined(i, 3); // Valgrind complains!
|
||||
* \endcode
|
||||
*
|
||||
* \tparam T The type of the object which ought to be set to undefined
|
||||
*
|
||||
* \param value Pointer to the first object of the array.
|
||||
* \param size The size of the array in number of objects
|
||||
*/
|
||||
template <class T>
|
||||
inline void SetUndefined(const T *value, int size)
|
||||
{
|
||||
#if !defined NDEBUG && HAVE_VALGRIND
|
||||
auto DUNE_UNUSED result = VALGRIND_MAKE_MEM_UNDEFINED(value, size*sizeof(T));
|
||||
#endif
|
||||
}
|
||||
|
||||
/*!
|
||||
* \ingroup Valgrind
|
||||
* \brief Make the memory on which an object resides defined.
|
||||
*
|
||||
* Example:
|
||||
*
|
||||
* \code
|
||||
* int i;
|
||||
* Opm::Valgrind::SetDefined(i);
|
||||
* Opm::Valgrind::CheckDefined(i); // Valgrind does not complain!
|
||||
* \endcode
|
||||
*
|
||||
* \tparam T The type of the object which valgrind should consider as defined
|
||||
*
|
||||
* \param value The object which's memory valgrind should consider as defined
|
||||
*/
|
||||
template <class T>
|
||||
inline void SetDefined(const T &value)
|
||||
{
|
||||
#if !defined NDEBUG && HAVE_VALGRIND
|
||||
auto DUNE_UNUSED result = VALGRIND_MAKE_MEM_DEFINED(&value, sizeof(T));
|
||||
#endif
|
||||
}
|
||||
|
||||
/*!
|
||||
* \ingroup Valgrind
|
||||
* \brief Make the memory on which a C-style array of objects resides
|
||||
* defined.
|
||||
*
|
||||
* Example:
|
||||
*
|
||||
* \code
|
||||
* int i[3];
|
||||
* Opm::Valgrind::SetDefined(i, 3);
|
||||
* Opm::Valgrind::CheckDefined(i, 3); // Valgrind does not complain!
|
||||
* \endcode
|
||||
*
|
||||
* \tparam T The type of the object which valgrind should consider as defined
|
||||
*
|
||||
* \param value Pointer to the first object of the array.
|
||||
* \param n The size of the array in number of objects
|
||||
*/
|
||||
template <class T>
|
||||
inline void SetDefined(const T *value, int n)
|
||||
{
|
||||
#if !defined NDEBUG && HAVE_VALGRIND
|
||||
auto DUNE_UNUSED result = VALGRIND_MAKE_MEM_DEFINED(value, n*sizeof(T));
|
||||
#endif
|
||||
}
|
||||
|
||||
/*!
|
||||
* \ingroup Valgrind
|
||||
* \brief Make valgrind complain if an object's memory is accessed.
|
||||
*
|
||||
* Example:
|
||||
*
|
||||
* \code
|
||||
* int i = 1;
|
||||
* Opm::Valgrind::SetNoAccess(i);
|
||||
* int j = i; // Valgrind complains!
|
||||
* \endcode
|
||||
*
|
||||
* \tparam T The type of the object which valgrind should complain if accessed
|
||||
*
|
||||
* \param value The object which's memory valgrind should complain if accessed
|
||||
*/
|
||||
template <class T>
|
||||
inline void SetNoAccess(const T &value)
|
||||
{
|
||||
#if !defined NDEBUG && HAVE_VALGRIND
|
||||
auto DUNE_UNUSED result = VALGRIND_MAKE_MEM_NOACCESS(&value, sizeof(T));
|
||||
#endif
|
||||
}
|
||||
|
||||
/*!
|
||||
* \ingroup Valgrind
|
||||
* \brief Make valgrind complain if the memory of a C-style array of
|
||||
* objects is accessed.
|
||||
*
|
||||
* Example:
|
||||
*
|
||||
* \code
|
||||
* int i[3] = {0, 1, 2};
|
||||
* Opm::Valgrind::SetNoAccess(i, 2);
|
||||
* int j = i[1]; // Valgrind complains!
|
||||
* \endcode
|
||||
*
|
||||
* \param value Pointer to the first object of the array.
|
||||
* \param size The size of the array in number of objects
|
||||
*
|
||||
* \param value The object which's memory valgrind should complain if accessed
|
||||
*/
|
||||
template <class T>
|
||||
inline void SetNoAccess(const T *value, int size)
|
||||
{
|
||||
#if !defined NDEBUG && HAVE_VALGRIND
|
||||
auto DUNE_UNUSED result = VALGRIND_MAKE_MEM_NOACCESS(value, size*sizeof(T));
|
||||
#endif
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
|
||||
#endif
|
||||
479
opm/common/variablelengthspline_.hh
Normal file
479
opm/common/variablelengthspline_.hh
Normal file
@ -0,0 +1,479 @@
|
||||
// -*- 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
|
||||
* \copydoc Opm::VariableLengthSpline_
|
||||
*/
|
||||
#ifndef OPM_VARIABLE_LENGTH_SPLINE_HH
|
||||
#define OPM_VARIABLE_LENGTH_SPLINE_HH
|
||||
|
||||
#include <dune/common/fvector.hh>
|
||||
#include <dune/common/fmatrix.hh>
|
||||
#include <dune/istl/bvector.hh>
|
||||
#include <dune/istl/btdmatrix.hh>
|
||||
|
||||
#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 ScalarT>
|
||||
class VariableLengthSpline_
|
||||
: public SplineCommon_<ScalarT,
|
||||
VariableLengthSpline_<ScalarT> >
|
||||
{
|
||||
friend class SplineCommon_<ScalarT, VariableLengthSpline_<ScalarT> >;
|
||||
|
||||
typedef ScalarT Scalar;
|
||||
typedef Dune::BlockVector<Dune::FieldVector<Scalar, 1> > Vector;
|
||||
typedef Dune::BTDMatrix<Dune::FieldMatrix<Scalar, 1, 1> > 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 <class ScalarArrayX, class ScalarArrayY>
|
||||
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 <class ScalarContainerX, class ScalarContainerY>
|
||||
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 <class PointArray>
|
||||
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 <class XYContainer>
|
||||
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 <class XYContainer>
|
||||
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 <class ScalarArrayX, class ScalarArrayY>
|
||||
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 <class ScalarContainerX, class ScalarContainerY>
|
||||
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 <class PointArray>
|
||||
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 <class XYContainer>
|
||||
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 <class XYContainer>
|
||||
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
|
||||
@ -21,8 +21,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::BinaryCoeff::Air_Mesitylene
|
||||
*/
|
||||
#ifndef EWOMS_BINARY_COEFF_AIR_MESITYLENE_HH
|
||||
#define EWOMS_BINARY_COEFF_AIR_MESITYLENE_HH
|
||||
#ifndef OPM_BINARY_COEFF_AIR_MESITYLENE_HH
|
||||
#define OPM_BINARY_COEFF_AIR_MESITYLENE_HH
|
||||
|
||||
#include <opm/material/components/air.hh>
|
||||
#include <opm/material/components/mesitylene.hh>
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::BinaryCoeff::Air_Xylene
|
||||
*/
|
||||
#ifndef EWOMS_BINARY_COEFF_AIR_XYLENE_HH
|
||||
#define EWOMS_BINARY_COEFF_AIR_XYLENE_HH
|
||||
#ifndef OPM_BINARY_COEFF_AIR_XYLENE_HH
|
||||
#define OPM_BINARY_COEFF_AIR_XYLENE_HH
|
||||
|
||||
#include <opm/material/components/air.hh>
|
||||
#include <opm/material/components/xylene.hh>
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
*
|
||||
* \copydoc Opm::BinaryCoeff::Brine_CO2
|
||||
*/
|
||||
#ifndef EWOMS_BINARY_COEFF_BRINE_CO2_HH
|
||||
#define EWOMS_BINARY_COEFF_BRINE_CO2_HH
|
||||
#ifndef OPM_BINARY_COEFF_BRINE_CO2_HH
|
||||
#define OPM_BINARY_COEFF_BRINE_CO2_HH
|
||||
|
||||
#include <opm/material/components/brine.hh>
|
||||
#include <opm/material/components/h2o.hh>
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::BinaryCoeff::fullerMethod
|
||||
*/
|
||||
#ifndef EWOMS_FULLERMETHOD_HH
|
||||
#define EWOMS_FULLERMETHOD_HH
|
||||
#ifndef OPM_FULLERMETHOD_HH
|
||||
#define OPM_FULLERMETHOD_HH
|
||||
|
||||
#include <opm/common/math.hh>
|
||||
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::BinaryCoeff::H2O_Air
|
||||
*/
|
||||
#ifndef EWOMS_BINARY_COEFF_H2O_AIR_HH
|
||||
#define EWOMS_BINARY_COEFF_H2O_AIR_HH
|
||||
#ifndef OPM_BINARY_COEFF_H2O_AIR_HH
|
||||
#define OPM_BINARY_COEFF_H2O_AIR_HH
|
||||
|
||||
#include <cmath>
|
||||
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
*
|
||||
* \copydoc Opm::BinaryCoeff::H2O_CO2
|
||||
*/
|
||||
#ifndef EWOMS_BINARY_COEFF_H2O_CO2_HH
|
||||
#define EWOMS_BINARY_COEFF_H2O_CO2_HH
|
||||
#ifndef OPM_BINARY_COEFF_H2O_CO2_HH
|
||||
#define OPM_BINARY_COEFF_H2O_CO2_HH
|
||||
|
||||
#include <opm/material/binarycoefficients/henryiapws.hh>
|
||||
#include <opm/material/binarycoefficients/fullermethod.hh>
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::BinaryCoeff::H2O_Mesitylene
|
||||
*/
|
||||
#ifndef EWOMS_BINARY_COEFF_H2O_MESITYLENE_HH
|
||||
#define EWOMS_BINARY_COEFF_H2O_MESITYLENE_HH
|
||||
#ifndef OPM_BINARY_COEFF_H2O_MESITYLENE_HH
|
||||
#define OPM_BINARY_COEFF_H2O_MESITYLENE_HH
|
||||
|
||||
#include <opm/material/components/h2o.hh>
|
||||
#include <opm/material/components/mesitylene.hh>
|
||||
@ -51,9 +51,9 @@ public:
|
||||
// after Sanders
|
||||
Scalar sanderH = 1.7e-1; // [M/atm]
|
||||
//conversion to our Henry definition
|
||||
Scalar ewomsH = sanderH / 101.325; // has now [(mol/m^3)/Pa]
|
||||
ewomsH *= 18.02e-6; // multiplied by molar volume of reference phase = water
|
||||
return 1.0/ewomsH; // [Pa]
|
||||
Scalar opmH = sanderH / 101.325; // has now [(mol/m^3)/Pa]
|
||||
opmH *= 18.02e-6; // multiplied by molar volume of reference phase = water
|
||||
return 1.0/opmH; // [Pa]
|
||||
}
|
||||
|
||||
/*!
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::BinaryCoeff::H2O_N2
|
||||
*/
|
||||
#ifndef EWOMS_BINARY_COEFF_H2O_N2_HH
|
||||
#define EWOMS_BINARY_COEFF_H2O_N2_HH
|
||||
#ifndef OPM_BINARY_COEFF_H2O_N2_HH
|
||||
#define OPM_BINARY_COEFF_H2O_N2_HH
|
||||
|
||||
#include "henryiapws.hh"
|
||||
#include "fullermethod.hh"
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::BinaryCoeff::H2O_Xylene
|
||||
*/
|
||||
#ifndef EWOMS_BINARY_COEFF_H2O_XYLENE_HH
|
||||
#define EWOMS_BINARY_COEFF_H2O_XYLENE_HH
|
||||
#ifndef OPM_BINARY_COEFF_H2O_XYLENE_HH
|
||||
#define OPM_BINARY_COEFF_H2O_XYLENE_HH
|
||||
|
||||
#include <opm/material/components/h2o.hh>
|
||||
#include <opm/material/components/xylene.hh>
|
||||
@ -52,9 +52,9 @@ public:
|
||||
// after Sanders
|
||||
Scalar sanderH = 1.5e-1; //[M/atm]
|
||||
//conversion to our Henry definition
|
||||
Scalar ewomsH = sanderH / 101.325; // has now [(mol/m^3)/Pa]
|
||||
ewomsH *= 18.02e-6; //multiplied by molar volume of reference phase = water
|
||||
return 1.0/ewomsH; // [Pa]
|
||||
Scalar opmH = sanderH / 101.325; // has now [(mol/m^3)/Pa]
|
||||
opmH *= 18.02e-6; //multiplied by molar volume of reference phase = water
|
||||
return 1.0/opmH; // [Pa]
|
||||
}
|
||||
|
||||
/*!
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
*
|
||||
* \brief The IAPWS formulation of Henry coefficients in water.
|
||||
*/
|
||||
#ifndef EWOMS_HENRY_IAPWS_HH
|
||||
#define EWOMS_HENRY_IAPWS_HH
|
||||
#ifndef OPM_HENRY_IAPWS_HH
|
||||
#define OPM_HENRY_IAPWS_HH
|
||||
|
||||
#include <opm/material/components/h2o.hh>
|
||||
|
||||
@ -81,4 +81,4 @@ inline Scalar henryIAPWS(Scalar E,
|
||||
}
|
||||
}
|
||||
|
||||
#endif // EWOMS_HENRY_IAPWS_HH
|
||||
#endif // OPM_HENRY_IAPWS_HH
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::Air
|
||||
*/
|
||||
#ifndef EWOMS_AIR_HH
|
||||
#define EWOMS_AIR_HH
|
||||
#ifndef OPM_AIR_HH
|
||||
#define OPM_AIR_HH
|
||||
|
||||
#include <dune/common/exceptions.hh>
|
||||
#include <opm/material/components/component.hh>
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
*
|
||||
* \copydoc Opm::Brine
|
||||
*/
|
||||
#ifndef EWOMS_BRINE_HH
|
||||
#define EWOMS_BRINE_HH
|
||||
#ifndef OPM_BRINE_HH
|
||||
#define OPM_BRINE_HH
|
||||
|
||||
#include <opm/material/components/component.hh>
|
||||
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
*
|
||||
* \copydoc Opm::CO2
|
||||
*/
|
||||
#ifndef EWOMS_CO2_HH
|
||||
#define EWOMS_CO2_HH
|
||||
#ifndef OPM_CO2_HH
|
||||
#define OPM_CO2_HH
|
||||
|
||||
#include <opm/common/exceptions.hh>
|
||||
#include <opm/material/components/component.hh>
|
||||
|
||||
@ -22,8 +22,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::Component
|
||||
*/
|
||||
#ifndef EWOMS_COMPONENT_HH
|
||||
#define EWOMS_COMPONENT_HH
|
||||
#ifndef OPM_COMPONENT_HH
|
||||
#define OPM_COMPONENT_HH
|
||||
|
||||
#include <dune/common/stdstreams.hh>
|
||||
|
||||
|
||||
@ -22,8 +22,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::DNAPL
|
||||
*/
|
||||
#ifndef EWOMS_DNAPL_HH
|
||||
#define EWOMS_DNAPL_HH
|
||||
#ifndef OPM_DNAPL_HH
|
||||
#define OPM_DNAPL_HH
|
||||
|
||||
#include <opm/material/idealgas.hh>
|
||||
#include "component.hh"
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::H2O
|
||||
*/
|
||||
#ifndef EWOMS_H2O_HH
|
||||
#define EWOMS_H2O_HH
|
||||
#ifndef OPM_H2O_HH
|
||||
#define OPM_H2O_HH
|
||||
|
||||
#include <cmath>
|
||||
#include <cassert>
|
||||
|
||||
@ -22,8 +22,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::IAPWS::Common
|
||||
*/
|
||||
#ifndef EWOMS_IAPWS_COMMON_HH
|
||||
#define EWOMS_IAPWS_COMMON_HH
|
||||
#ifndef OPM_IAPWS_COMMON_HH
|
||||
#define OPM_IAPWS_COMMON_HH
|
||||
|
||||
#include <opm/material/constants.hh>
|
||||
|
||||
|
||||
@ -22,8 +22,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::IAPWS::Region1
|
||||
*/
|
||||
#ifndef EWOMS_IAPWS_REGION1_HH
|
||||
#define EWOMS_IAPWS_REGION1_HH
|
||||
#ifndef OPM_IAPWS_REGION1_HH
|
||||
#define OPM_IAPWS_REGION1_HH
|
||||
|
||||
#include <cmath>
|
||||
|
||||
|
||||
@ -22,8 +22,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::IAPWS::Region2
|
||||
*/
|
||||
#ifndef EWOMS_IAPWS_REGION2_HH
|
||||
#define EWOMS_IAPWS_REGION2_HH
|
||||
#ifndef OPM_IAPWS_REGION2_HH
|
||||
#define OPM_IAPWS_REGION2_HH
|
||||
|
||||
#include <cmath>
|
||||
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::IAPWS::Region4
|
||||
*/
|
||||
#ifndef EWOMS_IAPWS_REGION4_HH
|
||||
#define EWOMS_IAPWS_REGION4_HH
|
||||
#ifndef OPM_IAPWS_REGION4_HH
|
||||
#define OPM_IAPWS_REGION4_HH
|
||||
|
||||
#include <cmath>
|
||||
|
||||
|
||||
@ -22,8 +22,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::LNAPL
|
||||
*/
|
||||
#ifndef EWOMS_LNAPL_HH
|
||||
#define EWOMS_LNAPL_HH
|
||||
#ifndef OPM_LNAPL_HH
|
||||
#define OPM_LNAPL_HH
|
||||
|
||||
#include "component.hh"
|
||||
|
||||
|
||||
@ -22,8 +22,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::Mesitylene
|
||||
*/
|
||||
#ifndef EWOMS_MESITYLENE_HH
|
||||
#define EWOMS_MESITYLENE_HH
|
||||
#ifndef OPM_MESITYLENE_HH
|
||||
#define OPM_MESITYLENE_HH
|
||||
|
||||
#include <opm/material/idealgas.hh>
|
||||
#include <opm/material/components/component.hh>
|
||||
|
||||
@ -22,8 +22,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::N2
|
||||
*/
|
||||
#ifndef EWOMS_N2_HH
|
||||
#define EWOMS_N2_HH
|
||||
#ifndef OPM_N2_HH
|
||||
#define OPM_N2_HH
|
||||
|
||||
#include <opm/material/idealgas.hh>
|
||||
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::NullComponent
|
||||
*/
|
||||
#ifndef EWOMS_NULL_COMPONENT_HH
|
||||
#define EWOMS_NULL_COMPONENT_HH
|
||||
#ifndef OPM_NULL_COMPONENT_HH
|
||||
#define OPM_NULL_COMPONENT_HH
|
||||
|
||||
#include "component.hh"
|
||||
|
||||
|
||||
@ -23,8 +23,8 @@
|
||||
*
|
||||
* \copydoc Opm::SimpleCO2
|
||||
*/
|
||||
#ifndef EWOMS_SIMPLE_CO2_HH
|
||||
#define EWOMS_SIMPLE_CO2_HH
|
||||
#ifndef OPM_SIMPLE_CO2_HH
|
||||
#define OPM_SIMPLE_CO2_HH
|
||||
|
||||
#include <opm/material/idealgas.hh>
|
||||
#include <opm/material/components/component.hh>
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::SimpleH2O
|
||||
*/
|
||||
#ifndef EWOMS_SIMPLE_H2O_HH
|
||||
#define EWOMS_SIMPLE_H2O_HH
|
||||
#ifndef OPM_SIMPLE_H2O_HH
|
||||
#define OPM_SIMPLE_H2O_HH
|
||||
|
||||
#include <opm/material/idealgas.hh>
|
||||
|
||||
|
||||
@ -23,8 +23,8 @@
|
||||
*
|
||||
* \copydoc Opm::TabulatedComponent
|
||||
*/
|
||||
#ifndef EWOMS_TABULATED_COMPONENT_HH
|
||||
#define EWOMS_TABULATED_COMPONENT_HH
|
||||
#ifndef OPM_TABULATED_COMPONENT_HH
|
||||
#define OPM_TABULATED_COMPONENT_HH
|
||||
|
||||
#include <cmath>
|
||||
#include <limits>
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::Unit
|
||||
*/
|
||||
#ifndef EWOMS_UNIT_HH
|
||||
#define EWOMS_UNIT_HH
|
||||
#ifndef OPM_UNIT_HH
|
||||
#define OPM_UNIT_HH
|
||||
|
||||
#include "component.hh"
|
||||
|
||||
|
||||
@ -22,8 +22,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::Xylene
|
||||
*/
|
||||
#ifndef EWOMS_XYLENE_HH
|
||||
#define EWOMS_XYLENE_HH
|
||||
#ifndef OPM_XYLENE_HH
|
||||
#define OPM_XYLENE_HH
|
||||
|
||||
#include <cmath>
|
||||
#include <opm/material/idealgas.hh>
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::Constants
|
||||
*/
|
||||
#ifndef EWOMS_CONSTANTS_HH
|
||||
#define EWOMS_CONSTANTS_HH
|
||||
#ifndef OPM_CONSTANTS_HH
|
||||
#define OPM_CONSTANTS_HH
|
||||
|
||||
#include <cmath>
|
||||
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::CompositionFromFugacities
|
||||
*/
|
||||
#ifndef EWOMS_COMPOSITION_FROM_FUGACITIES_HH
|
||||
#define EWOMS_COMPOSITION_FROM_FUGACITIES_HH
|
||||
#ifndef OPM_COMPOSITION_FROM_FUGACITIES_HH
|
||||
#define OPM_COMPOSITION_FROM_FUGACITIES_HH
|
||||
|
||||
#include <dune/common/fvector.hh>
|
||||
#include <dune/common/fmatrix.hh>
|
||||
|
||||
169
opm/material/constraintsolvers/computefromreferencephase.hh
Normal file
169
opm/material/constraintsolvers/computefromreferencephase.hh
Normal file
@ -0,0 +1,169 @@
|
||||
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
||||
// vi: set et ts=4 sw=4 sts=4:
|
||||
/*****************************************************************************
|
||||
* Copyright (C) 2011-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
|
||||
* \copydoc Opm::ComputeFromReferencePhase
|
||||
*/
|
||||
#ifndef OPM_COMPUTE_FROM_REFERENCE_PHASE_HH
|
||||
#define OPM_COMPUTE_FROM_REFERENCE_PHASE_HH
|
||||
|
||||
#include <opm/material/constraintsolvers/compositionfromfugacities.hh>
|
||||
|
||||
#include <opm/common/exceptions.hh>
|
||||
#include <opm/common/valgrind.hh>
|
||||
|
||||
#include <dune/common/fvector.hh>
|
||||
|
||||
namespace Opm {
|
||||
|
||||
/*!
|
||||
* \brief Computes all quantities of a generic fluid state if a
|
||||
* reference phase has been specified.
|
||||
*
|
||||
* This makes it is possible to specify just one phase and let the
|
||||
* remaining ones be calculated by the constraint solver. This
|
||||
* constraint solver assumes thermodynamic equilibrium. It assumes the
|
||||
* following quantities to be set:
|
||||
*
|
||||
* - composition (mole+mass fractions) of the *reference* phase
|
||||
* - temperature of the *reference* phase
|
||||
* - saturations of *all* phases
|
||||
* - pressures of *all* phases
|
||||
*
|
||||
* after calling the solve() method the following quantities are
|
||||
* calculated in addition:
|
||||
*
|
||||
* - temperature of *all* phases
|
||||
* - density, molar density, molar volume of *all* phases
|
||||
* - composition in mole and mass fractions and molaries of *all* phases
|
||||
* - mean molar masses of *all* phases
|
||||
* - fugacity coefficients of *all* components in *all* phases
|
||||
* - if the setViscosity parameter is true, also dynamic viscosities of *all* phases
|
||||
* - if the setEnthalpy parameter is true, also specific enthalpies and internal energies of *all* phases
|
||||
*/
|
||||
template <class Scalar, class FluidSystem>
|
||||
class ComputeFromReferencePhase
|
||||
{
|
||||
enum { numPhases = FluidSystem::numPhases };
|
||||
enum { numComponents = FluidSystem::numComponents };
|
||||
typedef Opm::CompositionFromFugacities<Scalar, FluidSystem> CompositionFromFugacities;
|
||||
typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
|
||||
|
||||
public:
|
||||
/*!
|
||||
* \brief Computes all quantities of a generic fluid state if a
|
||||
* reference phase has been specified.
|
||||
*
|
||||
* This makes it is possible to specify just one phase and let the
|
||||
* remaining ones be calculated by the constraint solver. This
|
||||
* constraint solver assumes thermodynamic equilibrium. It assumes the
|
||||
* following quantities to be set:
|
||||
*
|
||||
* - composition (mole+mass fractions) of the *reference* phase
|
||||
* - temperature of the *all* phases
|
||||
* - saturations of *all* phases
|
||||
* - pressures of *all* phases
|
||||
*
|
||||
* after calling the solve() method the following quantities are
|
||||
* calculated in addition:
|
||||
*
|
||||
* - temperature of *all* phases
|
||||
* - density, molar density, molar volume of *all* phases
|
||||
* - composition in mole and mass fractions and molaries of *all* phases
|
||||
* - mean molar masses of *all* phases
|
||||
* - fugacity coefficients of *all* components in *all* phases
|
||||
* - if the setViscosity parameter is true, also dynamic viscosities of *all* phases
|
||||
* - if the setEnthalpy parameter is true, also specific enthalpies and internal energies of *all* phases
|
||||
*
|
||||
* \param fluidState Thermodynamic state of the fluids
|
||||
* \param paramCache Container for cache parameters
|
||||
* \param refPhaseIdx The phase index of the reference phase
|
||||
* \param setViscosity Specify whether the dynamic viscosity of
|
||||
* each phase should also be set.
|
||||
* \param setEnthalpy Specify whether the specific
|
||||
* enthalpy/internal energy of each phase
|
||||
* should also be set.
|
||||
*/
|
||||
template <class FluidState, class ParameterCache>
|
||||
static void solve(FluidState &fluidState,
|
||||
ParameterCache ¶mCache,
|
||||
int refPhaseIdx,
|
||||
bool setViscosity,
|
||||
bool setEnthalpy)
|
||||
{
|
||||
|
||||
// compute the density and enthalpy of the
|
||||
// reference phase
|
||||
paramCache.updatePhase(fluidState, refPhaseIdx);
|
||||
fluidState.setDensity(refPhaseIdx,
|
||||
FluidSystem::density(fluidState,
|
||||
paramCache,
|
||||
refPhaseIdx));
|
||||
|
||||
if (setEnthalpy)
|
||||
fluidState.setEnthalpy(refPhaseIdx,
|
||||
FluidSystem::enthalpy(fluidState,
|
||||
paramCache,
|
||||
refPhaseIdx));
|
||||
|
||||
if (setViscosity)
|
||||
fluidState.setViscosity(refPhaseIdx,
|
||||
FluidSystem::viscosity(fluidState,
|
||||
paramCache,
|
||||
refPhaseIdx));
|
||||
|
||||
// compute the fugacities of all components in the reference phase
|
||||
for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
|
||||
fluidState.setFugacityCoefficient(refPhaseIdx,
|
||||
compIdx,
|
||||
FluidSystem::fugacityCoefficient(fluidState,
|
||||
paramCache,
|
||||
refPhaseIdx,
|
||||
compIdx));
|
||||
}
|
||||
|
||||
// compute all quantities for the non-reference phases
|
||||
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||
if (phaseIdx == refPhaseIdx)
|
||||
continue; // reference phase is already calculated
|
||||
|
||||
ComponentVector fugVec;
|
||||
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
|
||||
fugVec[compIdx] = fluidState.fugacity(refPhaseIdx, compIdx);
|
||||
|
||||
CompositionFromFugacities::solve(fluidState, paramCache, phaseIdx, fugVec);
|
||||
|
||||
if (setViscosity)
|
||||
fluidState.setViscosity(phaseIdx,
|
||||
FluidSystem::viscosity(fluidState,
|
||||
paramCache,
|
||||
phaseIdx));
|
||||
|
||||
if (setEnthalpy)
|
||||
fluidState.setEnthalpy(phaseIdx,
|
||||
FluidSystem::enthalpy(fluidState,
|
||||
paramCache,
|
||||
phaseIdx));
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
} // end namespace Opm
|
||||
|
||||
#endif
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::ImmiscibleFlash
|
||||
*/
|
||||
#ifndef EWOMS_IMMISCIBLE_FLASH_HH
|
||||
#define EWOMS_IMMISCIBLE_FLASH_HH
|
||||
#ifndef OPM_IMMISCIBLE_FLASH_HH
|
||||
#define OPM_IMMISCIBLE_FLASH_HH
|
||||
|
||||
#include <opm/common/exceptions.hh>
|
||||
#include <opm/common/valgrind.hh>
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::MiscibleMultiPhaseComposition
|
||||
*/
|
||||
#ifndef EWOMS_MISCIBLE_MULTIPHASE_COMPOSITION_HH
|
||||
#define EWOMS_MISCIBLE_MULTIPHASE_COMPOSITION_HH
|
||||
#ifndef OPM_MISCIBLE_MULTIPHASE_COMPOSITION_HH
|
||||
#define OPM_MISCIBLE_MULTIPHASE_COMPOSITION_HH
|
||||
|
||||
#include <dune/common/fvector.hh>
|
||||
#include <dune/common/fmatrix.hh>
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::NcpFlash
|
||||
*/
|
||||
#ifndef EWOMS_NCP_FLASH_HH
|
||||
#define EWOMS_NCP_FLASH_HH
|
||||
#ifndef OPM_NCP_FLASH_HH
|
||||
#define OPM_NCP_FLASH_HH
|
||||
|
||||
#include <dune/common/fvector.hh>
|
||||
#include <dune/common/fmatrix.hh>
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
*
|
||||
* \copydoc Opm::PengRobinson
|
||||
*/
|
||||
#ifndef EWOMS_PENG_ROBINSON_HH
|
||||
#define EWOMS_PENG_ROBINSON_HH
|
||||
#ifndef OPM_PENG_ROBINSON_HH
|
||||
#define OPM_PENG_ROBINSON_HH
|
||||
|
||||
#include <opm/material/fluidstates/temperatureoverlayfluidstate.hh>
|
||||
#include <opm/material/idealgas.hh>
|
||||
@ -30,6 +30,8 @@
|
||||
#include <opm/common/math.hh>
|
||||
#include <opm/common/dynamictabulated2dfunction.hh>
|
||||
|
||||
#include <dune/common/unused.hh>
|
||||
|
||||
#include <csignal>
|
||||
|
||||
namespace Opm {
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::PengRobinsonMixture
|
||||
*/
|
||||
#ifndef EWOMS_PENG_ROBINSON_MIXTURE_HH
|
||||
#define EWOMS_PENG_ROBINSON_MIXTURE_HH
|
||||
#ifndef OPM_PENG_ROBINSON_MIXTURE_HH
|
||||
#define OPM_PENG_ROBINSON_MIXTURE_HH
|
||||
|
||||
#include "pengrobinson.hh"
|
||||
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::PengRobinsonParams
|
||||
*/
|
||||
#ifndef EWOMS_PENG_ROBINSON_PARAMS_HH
|
||||
#define EWOMS_PENG_ROBINSON_PARAMS_HH
|
||||
#ifndef OPM_PENG_ROBINSON_PARAMS_HH
|
||||
#define OPM_PENG_ROBINSON_PARAMS_HH
|
||||
|
||||
#include <opm/common/valgrind.hh>
|
||||
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::PengRobinsonParamsMixture
|
||||
*/
|
||||
#ifndef EWOMS_PENG_ROBINSON_PARAMS_MIXTURE_HH
|
||||
#define EWOMS_PENG_ROBINSON_PARAMS_MIXTURE_HH
|
||||
#ifndef OPM_PENG_ROBINSON_PARAMS_MIXTURE_HH
|
||||
#define OPM_PENG_ROBINSON_PARAMS_MIXTURE_HH
|
||||
|
||||
#include <algorithm>
|
||||
#include <opm/material/constants.hh>
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::BrooksCorey
|
||||
*/
|
||||
#ifndef EWOMS_BROOKS_COREY_HH
|
||||
#define EWOMS_BROOKS_COREY_HH
|
||||
#ifndef OPM_BROOKS_COREY_HH
|
||||
#define OPM_BROOKS_COREY_HH
|
||||
|
||||
#include "brookscoreyparams.hh"
|
||||
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::BrooksCoreyParams
|
||||
*/
|
||||
#ifndef EWOMS_BROOKS_COREY_PARAMS_HH
|
||||
#define EWOMS_BROOKS_COREY_PARAMS_HH
|
||||
#ifndef OPM_BROOKS_COREY_PARAMS_HH
|
||||
#define OPM_BROOKS_COREY_PARAMS_HH
|
||||
|
||||
#include <opm/common/valgrind.hh>
|
||||
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::EffToAbsLaw
|
||||
*/
|
||||
#ifndef EWOMS_EFF_TO_ABS_LAW_HH
|
||||
#define EWOMS_EFF_TO_ABS_LAW_HH
|
||||
#ifndef OPM_EFF_TO_ABS_LAW_HH
|
||||
#define OPM_EFF_TO_ABS_LAW_HH
|
||||
|
||||
#include "efftoabslawparams.hh"
|
||||
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::EffToAbsLawParams
|
||||
*/
|
||||
#ifndef EWOMS_EFF_TO_ABS_LAW_PARAMS_HH
|
||||
#define EWOMS_EFF_TO_ABS_LAW_PARAMS_HH
|
||||
#ifndef OPM_EFF_TO_ABS_LAW_PARAMS_HH
|
||||
#define OPM_EFF_TO_ABS_LAW_PARAMS_HH
|
||||
|
||||
namespace Opm {
|
||||
/*!
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::LinearMaterial
|
||||
*/
|
||||
#ifndef EWOMS_LINEAR_MATERIAL_HH
|
||||
#define EWOMS_LINEAR_MATERIAL_HH
|
||||
#ifndef OPM_LINEAR_MATERIAL_HH
|
||||
#define OPM_LINEAR_MATERIAL_HH
|
||||
|
||||
#include "linearmaterialparams.hh"
|
||||
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::ParkerLenhard
|
||||
*/
|
||||
#ifndef EWOMS_PARKER_LENHARD_HH
|
||||
#define EWOMS_PARKER_LENHARD_HH
|
||||
#ifndef OPM_PARKER_LENHARD_HH
|
||||
#define OPM_PARKER_LENHARD_HH
|
||||
|
||||
#include "parkerlenhardparams.hh"
|
||||
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::ParkerLenhardParams
|
||||
*/
|
||||
#ifndef EWOMS_PARKER_LENHARD_PARAMS_HH
|
||||
#define EWOMS_PARKER_LENHARD_PARAMS_HH
|
||||
#ifndef OPM_PARKER_LENHARD_PARAMS_HH
|
||||
#define OPM_PARKER_LENHARD_PARAMS_HH
|
||||
|
||||
#include <opm/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh>
|
||||
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::RegularizedBrooksCoreyParams
|
||||
*/
|
||||
#ifndef EWOMS_REGULARIZED_BROOKS_COREY_PARAMS_HH
|
||||
#define EWOMS_REGULARIZED_BROOKS_COREY_PARAMS_HH
|
||||
#ifndef OPM_REGULARIZED_BROOKS_COREY_PARAMS_HH
|
||||
#define OPM_REGULARIZED_BROOKS_COREY_PARAMS_HH
|
||||
|
||||
#include "brookscoreyparams.hh"
|
||||
|
||||
|
||||
@ -22,8 +22,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::RegularizedLinearMaterial
|
||||
*/
|
||||
#ifndef EWOMS_REGULARIZED_LINEAR_MATERIAL_HH
|
||||
#define EWOMS_REGULARIZED_LINEAR_MATERIAL_HH
|
||||
#ifndef OPM_REGULARIZED_LINEAR_MATERIAL_HH
|
||||
#define OPM_REGULARIZED_LINEAR_MATERIAL_HH
|
||||
|
||||
#include "linearmaterial.hh"
|
||||
#include "regularizedlinearmaterialparams.hh"
|
||||
|
||||
@ -22,8 +22,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::RegularizedVanGenuchten
|
||||
*/
|
||||
#ifndef EWOMS_REGULARIZED_VAN_GENUCHTEN_HH
|
||||
#define EWOMS_REGULARIZED_VAN_GENUCHTEN_HH
|
||||
#ifndef OPM_REGULARIZED_VAN_GENUCHTEN_HH
|
||||
#define OPM_REGULARIZED_VAN_GENUCHTEN_HH
|
||||
|
||||
#include "vangenuchten.hh"
|
||||
#include "regularizedvangenuchtenparams.hh"
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::RegularizedVanGenuchtenParams
|
||||
*/
|
||||
#ifndef EWOMS_REGULARIZED_VAN_GENUCHTEN_PARAMS_HH
|
||||
#define EWOMS_REGULARIZED_VAN_GENUCHTEN_PARAMS_HH
|
||||
#ifndef OPM_REGULARIZED_VAN_GENUCHTEN_PARAMS_HH
|
||||
#define OPM_REGULARIZED_VAN_GENUCHTEN_PARAMS_HH
|
||||
|
||||
#include "vangenuchtenparams.hh"
|
||||
|
||||
|
||||
@ -22,8 +22,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::VanGenuchten
|
||||
*/
|
||||
#ifndef EWOMS_VAN_GENUCHTEN_HH
|
||||
#define EWOMS_VAN_GENUCHTEN_HH
|
||||
#ifndef OPM_VAN_GENUCHTEN_HH
|
||||
#define OPM_VAN_GENUCHTEN_HH
|
||||
|
||||
#include "vangenuchtenparams.hh"
|
||||
|
||||
|
||||
@ -22,8 +22,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::ThreePParkerVanGenuchten
|
||||
*/
|
||||
#ifndef EWOMS_3P_PARKER_VAN_GENUCHTEN_HH
|
||||
#define EWOMS_3P_PARKER_VAN_GENUCHTEN_HH
|
||||
#ifndef OPM_3P_PARKER_VAN_GENUCHTEN_HH
|
||||
#define OPM_3P_PARKER_VAN_GENUCHTEN_HH
|
||||
|
||||
#include "3pparkervangenuchtenparams.hh"
|
||||
|
||||
@ -61,8 +61,7 @@ public:
|
||||
Sw = wetting phase saturation, or,
|
||||
sum of wetting phase saturations
|
||||
alpha : VanGenuchten-alpha
|
||||
this function is just copied from MUFTE/pml/constrel3p3cni.c
|
||||
that is why variable names do not yet fulfill eWoms rules, TODO Change */
|
||||
this function is copied from MUFTE/pml/constrel3p3cni.c */
|
||||
|
||||
Scalar r,Se,x,vg_m;
|
||||
Scalar pc,pc_prime,Se_regu;
|
||||
@ -105,8 +104,7 @@ public:
|
||||
Sw = wetting phase saturation, or,
|
||||
sum of wetting phase saturations
|
||||
alpha : VanGenuchten-alpha
|
||||
this function is just copied from MUFTE/pml/constrel3p3cni.c
|
||||
that is why variable names do not yet fulfill eWoms rules, TODO Change */
|
||||
this function is just copied from MUFTE/pml/constrel3p3cni.c */
|
||||
|
||||
Scalar r,Se,x,vg_m;
|
||||
Scalar pc,pc_prime,Se_regu;
|
||||
@ -148,8 +146,7 @@ public:
|
||||
/*
|
||||
St = sum of wetting (liquid) phase saturations
|
||||
alpha : VanGenuchten-alpha
|
||||
this function is just copied from MUFTE/pml/constrel3p3cni.c
|
||||
that is why variable names do not yet fulfill eWoms rules, TODO Change */
|
||||
this function is just copied from MUFTE/pml/constrel3p3cni.c */
|
||||
|
||||
Scalar r,Se,x,vg_m;
|
||||
Scalar pc,pc_prime,Se_regu;
|
||||
|
||||
@ -22,8 +22,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::ParkerVanGen3PParams
|
||||
*/
|
||||
#ifndef EWOMS_3P_PARKER_VAN_GENUCHTEN_PARAMS_HH
|
||||
#define EWOMS_3P_PARKER_VAN_GENUCHTEN_PARAMS_HH
|
||||
#ifndef OPM_3P_PARKER_VAN_GENUCHTEN_PARAMS_HH
|
||||
#define OPM_3P_PARKER_VAN_GENUCHTEN_PARAMS_HH
|
||||
|
||||
#include <dune/common/fvector.hh>
|
||||
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::TwoPAdapter
|
||||
*/
|
||||
#ifndef EWOMS_MP_2P_ADAPTER_HH
|
||||
#define EWOMS_MP_2P_ADAPTER_HH
|
||||
#ifndef OPM_MP_2P_ADAPTER_HH
|
||||
#define OPM_MP_2P_ADAPTER_HH
|
||||
|
||||
#include <algorithm>
|
||||
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::ThreePAdapter
|
||||
*/
|
||||
#ifndef EWOMS_MP_3P_ADAPTER_HH
|
||||
#define EWOMS_MP_3P_ADAPTER_HH
|
||||
#ifndef OPM_MP_3P_ADAPTER_HH
|
||||
#define OPM_MP_3P_ADAPTER_HH
|
||||
|
||||
#include <dune/common/exceptions.hh>
|
||||
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::MpLinearMaterial
|
||||
*/
|
||||
#ifndef EWOMS_MP_LINEAR_MATERIAL_HH
|
||||
#define EWOMS_MP_LINEAR_MATERIAL_HH
|
||||
#ifndef OPM_MP_LINEAR_MATERIAL_HH
|
||||
#define OPM_MP_LINEAR_MATERIAL_HH
|
||||
|
||||
#include "mplinearmaterialparams.hh"
|
||||
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::NullMaterialLaw
|
||||
*/
|
||||
#ifndef EWOMS_NULL_MATERIAL_LAW_HH
|
||||
#define EWOMS_NULL_MATERIAL_LAW_HH
|
||||
#ifndef OPM_NULL_MATERIAL_LAW_HH
|
||||
#define OPM_NULL_MATERIAL_LAW_HH
|
||||
|
||||
#include "nullmateriallawparams.hh"
|
||||
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::NullMaterialLawParams
|
||||
*/
|
||||
#ifndef EWOMS_NULL_MATERIAL_LAW_PARAMS_HH
|
||||
#define EWOMS_NULL_MATERIAL_LAW_PARAMS_HH
|
||||
#ifndef OPM_NULL_MATERIAL_LAW_PARAMS_HH
|
||||
#define OPM_NULL_MATERIAL_LAW_PARAMS_HH
|
||||
|
||||
namespace Opm {
|
||||
/*!
|
||||
|
||||
@ -23,8 +23,8 @@
|
||||
* multi-phase, multi-component fluid system assuming
|
||||
* thermodynamic equilibrium.
|
||||
*/
|
||||
#ifndef EWOMS_COMPOSITIONAL_FLUID_STATE_HH
|
||||
#define EWOMS_COMPOSITIONAL_FLUID_STATE_HH
|
||||
#ifndef OPM_COMPOSITIONAL_FLUID_STATE_HH
|
||||
#define OPM_COMPOSITIONAL_FLUID_STATE_HH
|
||||
|
||||
#include "modularfluidstate.hh"
|
||||
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
*
|
||||
* \brief Modules for the ModularFluidState which represent composition.
|
||||
*/
|
||||
#ifndef EWOMS_FLUID_STATE_COMPOSITION_MODULES_HH
|
||||
#define EWOMS_FLUID_STATE_COMPOSITION_MODULES_HH
|
||||
#ifndef OPM_FLUID_STATE_COMPOSITION_MODULES_HH
|
||||
#define OPM_FLUID_STATE_COMPOSITION_MODULES_HH
|
||||
|
||||
#include <opm/common/valgrind.hh>
|
||||
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
*
|
||||
* \brief Modules for the ModularFluidState which represent density.
|
||||
*/
|
||||
#ifndef EWOMS_FLUID_STATE_DENSITY_MODULES_HH
|
||||
#define EWOMS_FLUID_STATE_DENSITY_MODULES_HH
|
||||
#ifndef OPM_FLUID_STATE_DENSITY_MODULES_HH
|
||||
#define OPM_FLUID_STATE_DENSITY_MODULES_HH
|
||||
|
||||
#include <opm/common/valgrind.hh>
|
||||
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
*
|
||||
* \brief Modules for the ModularFluidState which represent enthalpy.
|
||||
*/
|
||||
#ifndef EWOMS_FLUID_STATE_ENTHALPY_MODULES_HH
|
||||
#define EWOMS_FLUID_STATE_ENTHALPY_MODULES_HH
|
||||
#ifndef OPM_FLUID_STATE_ENTHALPY_MODULES_HH
|
||||
#define OPM_FLUID_STATE_ENTHALPY_MODULES_HH
|
||||
|
||||
#include <opm/common/valgrind.hh>
|
||||
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
*
|
||||
* \brief Modules for the ModularFluidState which represent fugacity/chemical potential.
|
||||
*/
|
||||
#ifndef EWOMS_FLUID_STATE_FUGACITY_MODULES_HH
|
||||
#define EWOMS_FLUID_STATE_FUGACITY_MODULES_HH
|
||||
#ifndef OPM_FLUID_STATE_FUGACITY_MODULES_HH
|
||||
#define OPM_FLUID_STATE_FUGACITY_MODULES_HH
|
||||
|
||||
#include <opm/common/valgrind.hh>
|
||||
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
*
|
||||
* \brief Modules for the ModularFluidState which represent pressure.
|
||||
*/
|
||||
#ifndef EWOMS_FLUID_STATE_PRESSURE_MODULES_HH
|
||||
#define EWOMS_FLUID_STATE_PRESSURE_MODULES_HH
|
||||
#ifndef OPM_FLUID_STATE_PRESSURE_MODULES_HH
|
||||
#define OPM_FLUID_STATE_PRESSURE_MODULES_HH
|
||||
|
||||
#include <opm/common/valgrind.hh>
|
||||
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
*
|
||||
* \brief Modules for the ModularFluidState which represent saturation.
|
||||
*/
|
||||
#ifndef EWOMS_FLUID_STATE_SATURATION_MODULES_HH
|
||||
#define EWOMS_FLUID_STATE_SATURATION_MODULES_HH
|
||||
#ifndef OPM_FLUID_STATE_SATURATION_MODULES_HH
|
||||
#define OPM_FLUID_STATE_SATURATION_MODULES_HH
|
||||
|
||||
#include <opm/common/valgrind.hh>
|
||||
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
*
|
||||
* \brief Modules for the ModularFluidState which represent temperature.
|
||||
*/
|
||||
#ifndef EWOMS_FLUID_STATE_TEMPERATURE_MODULES_HH
|
||||
#define EWOMS_FLUID_STATE_TEMPERATURE_MODULES_HH
|
||||
#ifndef OPM_FLUID_STATE_TEMPERATURE_MODULES_HH
|
||||
#define OPM_FLUID_STATE_TEMPERATURE_MODULES_HH
|
||||
|
||||
#include <opm/common/valgrind.hh>
|
||||
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
*
|
||||
* \brief Modules for the ModularFluidState which represent viscosity.
|
||||
*/
|
||||
#ifndef EWOMS_FLUID_STATE_VISCOSITY_MODULES_HH
|
||||
#define EWOMS_FLUID_STATE_VISCOSITY_MODULES_HH
|
||||
#ifndef OPM_FLUID_STATE_VISCOSITY_MODULES_HH
|
||||
#define OPM_FLUID_STATE_VISCOSITY_MODULES_HH
|
||||
|
||||
#include <opm/common/valgrind.hh>
|
||||
|
||||
|
||||
@ -23,8 +23,8 @@
|
||||
* multi-phase, multi-component fluid system assuming
|
||||
* thermodynamic equilibrium.
|
||||
*/
|
||||
#ifndef EWOMS_IMMISCIBLE_FLUID_STATE_HH
|
||||
#define EWOMS_IMMISCIBLE_FLUID_STATE_HH
|
||||
#ifndef OPM_IMMISCIBLE_FLUID_STATE_HH
|
||||
#define OPM_IMMISCIBLE_FLUID_STATE_HH
|
||||
|
||||
#include "modularfluidstate.hh"
|
||||
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::ModularFluidState
|
||||
*/
|
||||
#ifndef EWOMS_MODULAR_FLUID_STATE_HH
|
||||
#define EWOMS_MODULAR_FLUID_STATE_HH
|
||||
#ifndef OPM_MODULAR_FLUID_STATE_HH
|
||||
#define OPM_MODULAR_FLUID_STATE_HH
|
||||
|
||||
#include "fluidstatepressuremodules.hh"
|
||||
#include "fluidstatetemperaturemodules.hh"
|
||||
|
||||
@ -23,8 +23,8 @@
|
||||
* multi-phase, multi-component fluid system _not_ assuming
|
||||
* thermodynamic equilibrium.
|
||||
*/
|
||||
#ifndef EWOMS_NON_EQUILIBRIUM_FLUID_STATE_HH
|
||||
#define EWOMS_NON_EQUILIBRIUM_FLUID_STATE_HH
|
||||
#ifndef OPM_NON_EQUILIBRIUM_FLUID_STATE_HH
|
||||
#define OPM_NON_EQUILIBRIUM_FLUID_STATE_HH
|
||||
|
||||
#include "modularfluidstate.hh"
|
||||
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::PressureOverlayFluidState
|
||||
*/
|
||||
#ifndef EWOMS_PRESSURE_OVERLAY_FLUID_STATE_HH
|
||||
#define EWOMS_PRESSURE_OVERLAY_FLUID_STATE_HH
|
||||
#ifndef OPM_PRESSURE_OVERLAY_FLUID_STATE_HH
|
||||
#define OPM_PRESSURE_OVERLAY_FLUID_STATE_HH
|
||||
|
||||
#include <opm/common/valgrind.hh>
|
||||
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::SaturationOverlayFluidState
|
||||
*/
|
||||
#ifndef EWOMS_SATURATION_OVERLAY_FLUID_STATE_HH
|
||||
#define EWOMS_SATURATION_OVERLAY_FLUID_STATE_HH
|
||||
#ifndef OPM_SATURATION_OVERLAY_FLUID_STATE_HH
|
||||
#define OPM_SATURATION_OVERLAY_FLUID_STATE_HH
|
||||
|
||||
#include <opm/common/valgrind.hh>
|
||||
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::TemperatureOverlayFluidState
|
||||
*/
|
||||
#ifndef EWOMS_TEMPERATURE_OVERLAY_FLUID_STATE_HH
|
||||
#define EWOMS_TEMPERATURE_OVERLAY_FLUID_STATE_HH
|
||||
#ifndef OPM_TEMPERATURE_OVERLAY_FLUID_STATE_HH
|
||||
#define OPM_TEMPERATURE_OVERLAY_FLUID_STATE_HH
|
||||
|
||||
#include <opm/common/valgrind.hh>
|
||||
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::FluidSystems::OneP
|
||||
*/
|
||||
#ifndef EWOMS_1P_FLUIDSYSTEM_HH
|
||||
#define EWOMS_1P_FLUIDSYSTEM_HH
|
||||
#ifndef OPM_1P_FLUIDSYSTEM_HH
|
||||
#define OPM_1P_FLUIDSYSTEM_HH
|
||||
|
||||
#include "basefluidsystem.hh"
|
||||
#include "nullparametercache.hh"
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::FluidSystems::TwoPImmiscible
|
||||
*/
|
||||
#ifndef EWOMS_2P_IMMISCIBLE_FLUID_SYSTEM_HH
|
||||
#define EWOMS_2P_IMMISCIBLE_FLUID_SYSTEM_HH
|
||||
#ifndef OPM_2P_IMMISCIBLE_FLUID_SYSTEM_HH
|
||||
#define OPM_2P_IMMISCIBLE_FLUID_SYSTEM_HH
|
||||
|
||||
#include <limits>
|
||||
#include <cassert>
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::BaseFluidSystem
|
||||
*/
|
||||
#ifndef EWOMS_BASE_FLUID_SYSTEM_HH
|
||||
#define EWOMS_BASE_FLUID_SYSTEM_HH
|
||||
#ifndef OPM_BASE_FLUID_SYSTEM_HH
|
||||
#define OPM_BASE_FLUID_SYSTEM_HH
|
||||
|
||||
#include "nullparametercache.hh"
|
||||
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::FluidSystems::BlackOil
|
||||
*/
|
||||
#ifndef EWOMS_BLACK_OIL_FLUID_SYSTEM_HH
|
||||
#define EWOMS_BLACK_OIL_FLUID_SYSTEM_HH
|
||||
#ifndef OPM_BLACK_OIL_FLUID_SYSTEM_HH
|
||||
#define OPM_BLACK_OIL_FLUID_SYSTEM_HH
|
||||
|
||||
#include <opm/common/exceptions.hh>
|
||||
#include <opm/common/spline.hh>
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
*
|
||||
* \copydoc Opm::FluidSystems::BrineCO2
|
||||
*/
|
||||
#ifndef EWOMS_BRINE_CO2_SYSTEM_HH
|
||||
#define EWOMS_BRINE_CO2_SYSTEM_HH
|
||||
#ifndef OPM_BRINE_CO2_SYSTEM_HH
|
||||
#define OPM_BRINE_CO2_SYSTEM_HH
|
||||
|
||||
#include "basefluidsystem.hh"
|
||||
#include "nullparametercache.hh"
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::GasPhase
|
||||
*/
|
||||
#ifndef EWOMS_GAS_PHASE_HH
|
||||
#define EWOMS_GAS_PHASE_HH
|
||||
#ifndef OPM_GAS_PHASE_HH
|
||||
#define OPM_GAS_PHASE_HH
|
||||
|
||||
namespace Opm {
|
||||
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::FluidSystems::H2OAir
|
||||
*/
|
||||
#ifndef EWOMS_H2O_AIR_SYSTEM_HH
|
||||
#define EWOMS_H2O_AIR_SYSTEM_HH
|
||||
#ifndef OPM_H2O_AIR_SYSTEM_HH
|
||||
#define OPM_H2O_AIR_SYSTEM_HH
|
||||
|
||||
#include "basefluidsystem.hh"
|
||||
#include "nullparametercache.hh"
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::FluidSystems::H2OAirMesitylene
|
||||
*/
|
||||
#ifndef EWOMS_H2O_AIR_MESITYLENE_FLUID_SYSTEM_HH
|
||||
#define EWOMS_H2O_AIR_MESITYLENE_FLUID_SYSTEM_HH
|
||||
#ifndef OPM_H2O_AIR_MESITYLENE_FLUID_SYSTEM_HH
|
||||
#define OPM_H2O_AIR_MESITYLENE_FLUID_SYSTEM_HH
|
||||
|
||||
#include "basefluidsystem.hh"
|
||||
#include "nullparametercache.hh"
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::FluidSystems::H2OAirXylene
|
||||
*/
|
||||
#ifndef EWOMS_H2O_AIR_XYLENE_FLUID_SYSTEM_HH
|
||||
#define EWOMS_H2O_AIR_XYLENE_FLUID_SYSTEM_HH
|
||||
#ifndef OPM_H2O_AIR_XYLENE_FLUID_SYSTEM_HH
|
||||
#define OPM_H2O_AIR_XYLENE_FLUID_SYSTEM_HH
|
||||
|
||||
#include <opm/material/idealgas.hh>
|
||||
#include <opm/material/components/air.hh>
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::FluidSystems::H2ON2
|
||||
*/
|
||||
#ifndef EWOMS_H2O_N2_FLUID_SYSTEM_HH
|
||||
#define EWOMS_H2O_N2_FLUID_SYSTEM_HH
|
||||
#ifndef OPM_H2O_N2_FLUID_SYSTEM_HH
|
||||
#define OPM_H2O_N2_FLUID_SYSTEM_HH
|
||||
|
||||
#include "basefluidsystem.hh"
|
||||
#include "nullparametercache.hh"
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::FluidSystems::H2ON2LiquidPhase
|
||||
*/
|
||||
#ifndef EWOMS_H2O_N2_LIQUIDPHASE_FLUID_SYSTEM_HH
|
||||
#define EWOMS_H2O_N2_LIQUIDPHASE_FLUID_SYSTEM_HH
|
||||
#ifndef OPM_H2O_N2_LIQUIDPHASE_FLUID_SYSTEM_HH
|
||||
#define OPM_H2O_N2_LIQUIDPHASE_FLUID_SYSTEM_HH
|
||||
|
||||
#include "basefluidsystem.hh"
|
||||
#include "nullparametercache.hh"
|
||||
|
||||
@ -21,8 +21,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::LiquidPhase
|
||||
*/
|
||||
#ifndef EWOMS_LIQUID_PHASE_HH
|
||||
#define EWOMS_LIQUID_PHASE_HH
|
||||
#ifndef OPM_LIQUID_PHASE_HH
|
||||
#define OPM_LIQUID_PHASE_HH
|
||||
|
||||
namespace Opm {
|
||||
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::NullParameterCache
|
||||
*/
|
||||
#ifndef EWOMS_NULL_PARAMETER_CACHE_HH
|
||||
#define EWOMS_NULL_PARAMETER_CACHE_HH
|
||||
#ifndef OPM_NULL_PARAMETER_CACHE_HH
|
||||
#define OPM_NULL_PARAMETER_CACHE_HH
|
||||
|
||||
#include "parametercachebase.hh"
|
||||
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::ParameterCacheBase
|
||||
*/
|
||||
#ifndef EWOMS_PARAMETER_CACHE_BASE_HH
|
||||
#define EWOMS_PARAMETER_CACHE_BASE_HH
|
||||
#ifndef OPM_PARAMETER_CACHE_BASE_HH
|
||||
#define OPM_PARAMETER_CACHE_BASE_HH
|
||||
|
||||
namespace Opm {
|
||||
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::FluidSystems::Spe5
|
||||
*/
|
||||
#ifndef EWOMS_SPE5_FLUID_SYSTEM_HH
|
||||
#define EWOMS_SPE5_FLUID_SYSTEM_HH
|
||||
#ifndef OPM_SPE5_FLUID_SYSTEM_HH
|
||||
#define OPM_SPE5_FLUID_SYSTEM_HH
|
||||
|
||||
#include "basefluidsystem.hh"
|
||||
#include "spe5parametercache.hh"
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::Spe5ParameterCache
|
||||
*/
|
||||
#ifndef EWOMS_SPE5_PARAMETER_CACHE_HH
|
||||
#define EWOMS_SPE5_PARAMETER_CACHE_HH
|
||||
#ifndef OPM_SPE5_PARAMETER_CACHE_HH
|
||||
#define OPM_SPE5_PARAMETER_CACHE_HH
|
||||
|
||||
#include <cassert>
|
||||
|
||||
|
||||
@ -20,8 +20,8 @@
|
||||
* \file
|
||||
* \copydoc Opm::DummyHeatConductionLaw
|
||||
*/
|
||||
#ifndef EWOMS_DUMMY_HEATCONDUCTION_LAW_HH
|
||||
#define EWOMS_DUMMY_HEATCONDUCTION_LAW_HH
|
||||
#ifndef OPM_DUMMY_HEATCONDUCTION_LAW_HH
|
||||
#define OPM_DUMMY_HEATCONDUCTION_LAW_HH
|
||||
|
||||
#include <dune/common/exceptions.hh>
|
||||
|
||||
|
||||
Some files were not shown because too many files have changed in this diff Show More
Loading…
Reference in New Issue
Block a user