20#include <dune/common/typetraits.hh>
21#include <dune/common/fvector.hh>
22#include <dune/common/fmatrix.hh>
23#include <dune/common/dynmatrix.hh>
24#include <dune/common/float_cmp.hh>
37template <
class Scalar>
38constexpr Scalar
arithmeticMean(Scalar x, Scalar y, Scalar wx = 1.0, Scalar wy = 1.0) noexcept
40 static_assert(Dune::IsNumber<Scalar>::value,
"The arguments x, y, wx, and wy have to be numbers!");
44 return (x * wx + y * wy)/(wx + wy);
56template <
class Scalar>
57constexpr Scalar
harmonicMean(Scalar x, Scalar y, Scalar wx = 1.0, Scalar wy = 1.0) noexcept
59 static_assert(Dune::IsNumber<Scalar>::value,
"The arguments x, y, wx, and wy have to be numbers!");
63 return (wx + wy) * x * y / (wy * x + wx * y);
74template <
class Scalar>
77 static_assert(Dune::IsNumber<Scalar>::value,
"The arguments x and y have to be numbers!");
82 return sqrt(x*y)*
sign(x);
96template <
class Scalar,
int m,
int n>
98 const Dune::FieldMatrix<Scalar, m, n> &Ki,
99 const Dune::FieldMatrix<Scalar, m, n> &Kj)
101 for (
int rowIdx=0; rowIdx < m; rowIdx++){
102 for (
int colIdx=0; colIdx< n; colIdx++){
103 if (Dune::FloatCmp::ne<Scalar>(Ki[rowIdx][colIdx], Kj[rowIdx][colIdx])) {
109 K[rowIdx][colIdx] = Ki[rowIdx][colIdx];
128template <
class Scalar,
class SolContainer>
156template <
class Scalar,
class SolContainer>
167 Scalar Delta = b*b - 4*a*c;
173 sol[0] = (- b + Delta)/(2*a);
174 sol[1] = (- b - Delta)/(2*a);
180 swap(sol[0], sol[1]);
186template <
class Scalar,
class SolContainer>
187void invertCubicPolynomialPostProcess_(SolContainer &sol,
189 Scalar a, Scalar b, Scalar c, Scalar d,
190 std::size_t numIterations = 1)
192 const auto eval = [&](
auto x){
return d + x*(c + x*(b + x*a)); };
193 const auto evalDeriv = [&](
auto x){
return c + x*(2*b + x*3*a); };
197 for (
int i = 0; i < numSol; ++i)
200 Scalar fCurrent = eval(x);
201 const Scalar fOld = fCurrent;
202 for (
int j = 0; j < numIterations; ++j)
204 const Scalar fPrime = evalDeriv(x);
207 x -= fCurrent/fPrime;
212 if (abs(fCurrent) < abs(fOld))
239template <
class Scalar,
class SolContainer>
241 Scalar a, Scalar b, Scalar c, Scalar d,
242 std::size_t numPostProcessIterations = 1)
255 Scalar p = c - b*b/3;
256 Scalar q = d + (2*b*b*b - 9*b*c)/27;
260 if (p == 0.0 && q == 0.0) {
262 sol[0] = sol[1] = sol[2] = 0.0 - b/3;
265 else if (p == 0.0 && q != 0.0) {
275 else if (p != 0.0 && q == 0.0) {
286 sol[0] = -sqrt(-p) - b/3;
288 sol[2] = sqrt(-p) - b/3;
323 Scalar wDisc = q*q/4 + p*p*p/27;
333 const Scalar u = [&]{
334 return q > 0 ? cbrt(-0.5*q - sqrt(wDisc)) : cbrt(-0.5*q + sqrt(wDisc));
338 sol[0] = u - p/(3*u) - b/3;
342 invertCubicPolynomialPostProcess_(sol, 1, a, b, c, d, numPostProcessIterations);
346 Scalar uCubedRe = - q/2;
348 Scalar uCubedIm = sqrt(-wDisc);
351 Scalar uAbs = cbrt(sqrt(uCubedRe*uCubedRe + uCubedIm*uCubedIm));
353 Scalar phi = atan2(uCubedIm, uCubedRe)/3;
394 for (
int i = 0; i < 3; ++i) {
396 sol[i] = cos(phi)*(uAbs - p/(3*uAbs)) - b/3;
402 invertCubicPolynomialPostProcess_(sol, 3, a, b, c, d, numPostProcessIterations);
427template <
class Scalar,
int dim>
428bool isLarger(
const Dune::FieldVector<Scalar, dim> &pos,
429 const Dune::FieldVector<Scalar, dim> &smallerVec)
431 for (
int i=0; i < dim; i++)
433 if (pos[i]<= smallerVec[i])
453template <
class Scalar,
int dim>
454bool isSmaller(
const Dune::FieldVector<Scalar, dim> &pos,
455 const Dune::FieldVector<Scalar, dim> &largerVec)
457 for (
int i=0; i < dim; i++)
459 if (pos[i]>= largerVec[i])
483template <
class Scalar,
int dim>
484bool isBetween(
const Dune::FieldVector<Scalar, dim> &pos,
485 const Dune::FieldVector<Scalar, dim> &smallerVec,
486 const Dune::FieldVector<Scalar, dim> &largerVec)
497namespace InterpolationPolicy {
struct Linear; }
513namespace InterpolationPolicy {
526 template<
class Scalar>
527 static constexpr Scalar
interpolate(Scalar ip,
const std::array<Scalar, 2>& params)
529 return params[0]*(1.0 - ip) + params[1]*ip;
546 template<
class Scalar,
class RandomAccessContainer0,
class RandomAccessContainer1>
547 static constexpr Scalar
interpolate(Scalar ip,
const RandomAccessContainer0& range,
const RandomAccessContainer1& values)
550 if (ip > range.back())
return values.back();
551 if (ip < range[0])
return values[0];
554 const auto lookUpIndex =
std::distance(range.begin(), std::lower_bound(range.begin(), range.end(), ip));
555 if (lookUpIndex == 0)
558 const auto ipLinear = (ip - range[lookUpIndex-1])/(range[lookUpIndex] - range[lookUpIndex-1]);
559 return Dumux::interpolate<Linear>(ipLinear, std::array<Scalar, 2>{{values[lookUpIndex-1], values[lookUpIndex]}});
562 template<
class Scalar,
class RandomAccessContainer>
563 static constexpr Scalar
interpolate(Scalar ip,
const std::pair<RandomAccessContainer, RandomAccessContainer>& table)
565 const auto& [range, values] = table;
581template <
class Scalar>
582std::vector<Scalar>
linspace(
const Scalar begin,
const Scalar end,
584 bool endPoint =
true)
587 samples = max(std::size_t{2}, samples);
588 const Scalar divisor = endPoint ? samples-1 : samples;
589 const Scalar delta = (end-begin)/divisor;
590 std::vector<Scalar> vec(samples);
591 for (std::size_t i = 0; i < samples; ++i)
592 vec[i] = begin + i*delta;
609template <
class Scalar>
615 const Scalar ln10 = 2.3025850929940459;
628template<
class ValueType>
629constexpr int sign(
const ValueType& value)
noexcept
631 return (ValueType(0) < value) - (value < ValueType(0));
641template <
class Scalar>
642Dune::FieldVector<Scalar, 3>
crossProduct(
const Dune::FieldVector<Scalar, 3> &vec1,
643 const Dune::FieldVector<Scalar, 3> &vec2)
645 return {vec1[1]*vec2[2]-vec1[2]*vec2[1],
646 vec1[2]*vec2[0]-vec1[0]*vec2[2],
647 vec1[0]*vec2[1]-vec1[1]*vec2[0]};
657template <
class Scalar>
659 const Dune::FieldVector<Scalar, 2> &vec2)
660{
return vec1[0]*vec2[1]-vec1[1]*vec2[0]; }
670template <
class Scalar>
672 const Dune::FieldVector<Scalar, 3> &vec2,
673 const Dune::FieldVector<Scalar, 3> &vec3)
674{
return crossProduct<Scalar>(vec1, vec2)*vec3; }
682template <
class Scalar,
int m,
int n>
683Dune::FieldMatrix<Scalar, n, m>
getTransposed(
const Dune::FieldMatrix<Scalar, m, n>& M)
685 Dune::FieldMatrix<Scalar, n, m> T;
686 for (std::size_t i = 0; i < m; ++i)
687 for (std::size_t j = 0; j < n; ++j)
699template <
class Scalar>
700Dune::DynamicMatrix<Scalar>
getTransposed(
const Dune::DynamicMatrix<Scalar>& M)
702 std::size_t rows_T = M.M();
703 std::size_t cols_T = M.N();
705 Dune::DynamicMatrix<Scalar> M_T(rows_T, cols_T, 0.0);
707 for (std::size_t i = 0; i < rows_T; ++i)
708 for (std::size_t j = 0; j < cols_T; ++j)
721template <
class Scalar>
723 const Dune::DynamicMatrix<Scalar> &M2)
725 using size_type =
typename Dune::DynamicMatrix<Scalar>::size_type;
726 const size_type rows = M1.N();
727 const size_type cols = M2.M();
729 DUNE_ASSERT_BOUNDS(M1.M() == M2.N());
731 Dune::DynamicMatrix<Scalar> result(rows, cols, 0.0);
732 for (size_type i = 0; i < rows; i++)
733 for (size_type j = 0; j < cols; j++)
734 for (size_type k = 0; k < M1.M(); k++)
735 result[i][j] += M1[i][k]*M2[k][j];
747template <
class Scalar,
int rows1,
int cols1,
int cols2>
748Dune::FieldMatrix<Scalar, rows1, cols2>
multiplyMatrices(
const Dune::FieldMatrix<Scalar, rows1, cols1> &M1,
749 const Dune::FieldMatrix<Scalar, cols1, cols2> &M2)
751 using size_type =
typename Dune::FieldMatrix<Scalar, rows1, cols2>::size_type;
753 Dune::FieldMatrix<Scalar, rows1, cols2> result(0.0);
754 for (size_type i = 0; i < rows1; i++)
755 for (size_type j = 0; j < cols2; j++)
756 for (size_type k = 0; k < cols1; k++)
757 result[i][j] += M1[i][k]*M2[k][j];
769template <
class MatrixType>
770typename Dune::DenseMatrix<MatrixType>::field_type
771trace(
const Dune::DenseMatrix<MatrixType>& M)
773 const auto rows = M.N();
774 DUNE_ASSERT_BOUNDS(rows == M.M());
776 using MatType = Dune::DenseMatrix<MatrixType>;
777 typename MatType::field_type
trace = 0.0;
779 for (
typename MatType::size_type i = 0; i < rows; ++i)
798template <
class MAT,
class V>
799typename Dune::DenseVector<V>::derived_type
800mv(
const Dune::DenseMatrix<MAT>& M,
801 const Dune::DenseVector<V>& v)
803 typename Dune::DenseVector<V>::derived_type res(v);
825template <
class FieldScalar,
class V>
826typename std::enable_if_t<Dune::IsNumber<FieldScalar>::value,
827 typename Dune::DenseVector<V>::derived_type>
828mv(
const FieldScalar m,
const Dune::DenseVector<V>& v)
830 typename Dune::DenseVector<V>::derived_type res(v);
849template <
class V1,
class MAT,
class V2>
850typename Dune::DenseMatrix<MAT>::value_type
851vtmv(
const Dune::DenseVector<V1>& v1,
852 const Dune::DenseMatrix<MAT>& M,
853 const Dune::DenseVector<V2>& v2)
877template <
class V1,
class FieldScalar,
class V2>
878typename std::enable_if_t<Dune::IsNumber<FieldScalar>::value, FieldScalar>
879vtmv(
const Dune::DenseVector<V1>& v1,
881 const Dune::DenseVector<V2>& v2)
897template <
class Scalar>
899 const std::vector<Scalar>& y)
901 if (x.size() != y.size())
902 DUNE_THROW(Dune::InvalidStateException,
"x and y array must have the same length.");
904 const Scalar averageX = std::accumulate(x.begin(), x.end(), 0.0)/x.size();
905 const Scalar averageY = std::accumulate(y.begin(), y.end(), 0.0)/y.size();
908 const Scalar numerator = std::inner_product(
909 x.begin(), x.end(), y.begin(), 0.0, std::plus<Scalar>(),
910 [&](
auto xx,
auto yy) { return (xx - averageX) * (yy - averageY); }
912 const Scalar denominator = std::inner_product(
913 x.begin(), x.end(), x.begin(), 0.0, std::plus<Scalar>(),
914 [&](
auto xx,
auto yy) { return (xx - averageX) * (yy - averageX); }
918 const Scalar slope = numerator / denominator;
919 const Scalar intercept = averageY - slope * averageX;
921 return {intercept, slope};
constexpr Scalar harmonicMean(Scalar x, Scalar y, Scalar wx=1.0, Scalar wy=1.0) noexcept
Calculate the (weighted) harmonic mean of two scalar values.
Definition: math.hh:57
int invertCubicPolynomial(SolContainer *sol, Scalar a, Scalar b, Scalar c, Scalar d, std::size_t numPostProcessIterations=1)
Invert a cubic polynomial analytically.
Definition: math.hh:240
Dune::DenseMatrix< MAT >::value_type vtmv(const Dune::DenseVector< V1 > &v1, const Dune::DenseMatrix< MAT > &M, const Dune::DenseVector< V2 > &v2)
Evaluates the scalar product of a vector v2, projected by a matrix M, with a vector v1.
Definition: math.hh:851
Dune::DynamicMatrix< Scalar > multiplyMatrices(const Dune::DynamicMatrix< Scalar > &M1, const Dune::DynamicMatrix< Scalar > &M2)
Multiply two dynamic matrices.
Definition: math.hh:722
Scalar geometricMean(Scalar x, Scalar y) noexcept
Calculate the geometric mean of two scalar values.
Definition: math.hh:75
int invertLinearPolynomial(SolContainer &sol, Scalar a, Scalar b)
Invert a linear polynomial analytically.
Definition: math.hh:129
bool isBetween(const Dune::FieldVector< Scalar, dim > &pos, const Dune::FieldVector< Scalar, dim > &smallerVec, const Dune::FieldVector< Scalar, dim > &largerVec)
Comparison of three position vectors.
Definition: math.hh:484
void harmonicMeanMatrix(Dune::FieldMatrix< Scalar, m, n > &K, const Dune::FieldMatrix< Scalar, m, n > &Ki, const Dune::FieldMatrix< Scalar, m, n > &Kj)
Calculate the harmonic mean of a fixed-size matrix.
Definition: math.hh:97
Dune::FieldVector< Scalar, 3 > crossProduct(const Dune::FieldVector< Scalar, 3 > &vec1, const Dune::FieldVector< Scalar, 3 > &vec2)
Cross product of two vectors in three-dimensional Euclidean space.
Definition: math.hh:642
Scalar interpolate(Scalar ip, Parameter &&... params)
a generic function to interpolate given a set of parameters and an interpolation point
Definition: math.hh:506
bool isSmaller(const Dune::FieldVector< Scalar, dim > &pos, const Dune::FieldVector< Scalar, dim > &largerVec)
Comparison of two position vectors.
Definition: math.hh:454
Dune::DenseVector< V >::derived_type mv(const Dune::DenseMatrix< MAT > &M, const Dune::DenseVector< V > &v)
Returns the result of the projection of a vector v with a Matrix M.
Definition: math.hh:800
Dune::DenseMatrix< MatrixType >::field_type trace(const Dune::DenseMatrix< MatrixType > &M)
Trace of a dense matrix.
Definition: math.hh:771
Dune::FieldMatrix< Scalar, n, m > getTransposed(const Dune::FieldMatrix< Scalar, m, n > &M)
Transpose a FieldMatrix.
Definition: math.hh:683
bool isLarger(const Dune::FieldVector< Scalar, dim > &pos, const Dune::FieldVector< Scalar, dim > &smallerVec)
Comparison of two position vectors.
Definition: math.hh:428
int invertQuadraticPolynomial(SolContainer &sol, Scalar a, Scalar b, Scalar c)
Invert a quadratic polynomial analytically.
Definition: math.hh:157
std::array< Scalar, 2 > linearRegression(const std::vector< Scalar > &x, const std::vector< Scalar > &y)
Returns the [intercept, slope] of the regression line fitted to a set of (x, y) data points.
Definition: math.hh:898
std::vector< Scalar > linspace(const Scalar begin, const Scalar end, std::size_t samples, bool endPoint=true)
Generates linearly spaced vectors.
Definition: math.hh:582
Scalar tripleProduct(const Dune::FieldVector< Scalar, 3 > &vec1, const Dune::FieldVector< Scalar, 3 > &vec2, const Dune::FieldVector< Scalar, 3 > &vec3)
Triple product of three vectors in three-dimensional Euclidean space retuning scalar.
Definition: math.hh:671
constexpr int sign(const ValueType &value) noexcept
Sign or signum function.
Definition: math.hh:629
Scalar antoine(Scalar temperature, Scalar A, Scalar B, Scalar C)
Evaluates the Antoine equation used to calculate the vapour pressure of various liquids.
Definition: math.hh:610
constexpr Scalar arithmeticMean(Scalar x, Scalar y, Scalar wx=1.0, Scalar wy=1.0) noexcept
Calculate the (weighted) arithmetic mean of two scalar values.
Definition: math.hh:38
static ctype distance(const Dune::FieldVector< ctype, dimWorld > &a, const Dune::FieldVector< ctype, dimWorld > &b)
Compute the shortest distance between two points.
Definition: distance.hh:282
std::string temperature() noexcept
I/O name of temperature for equilibrium models.
Definition: name.hh:39
interpolate linearly between two given values
Definition: math.hh:520
static constexpr Scalar interpolate(Scalar ip, const std::array< Scalar, 2 > ¶ms)
interpolate linearly between two given values
Definition: math.hh:527
interpolate linearly in a piecewise linear function (tabularized function)
Definition: math.hh:538
static constexpr Scalar interpolate(Scalar ip, const std::pair< RandomAccessContainer, RandomAccessContainer > &table)
Definition: math.hh:563
static constexpr Scalar interpolate(Scalar ip, const RandomAccessContainer0 &range, const RandomAccessContainer1 &values)
interpolate linearly in a piecewise linear function (tabularized function)
Definition: math.hh:547