version 3.8
math.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3//
4// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5// SPDX-License-Identifier: GPL-3.0-or-later
6//
12#ifndef DUMUX_MATH_HH
13#define DUMUX_MATH_HH
14
15#include <algorithm>
16#include <numeric>
17#include <cmath>
18#include <utility>
19
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>
25
26namespace Dumux {
27
37template <class Scalar>
38constexpr Scalar arithmeticMean(Scalar x, Scalar y, Scalar wx = 1.0, Scalar wy = 1.0) noexcept
39{
40 static_assert(Dune::IsNumber<Scalar>::value, "The arguments x, y, wx, and wy have to be numbers!");
41
42 if (x*y <= 0)
43 return 0;
44 return (x * wx + y * wy)/(wx + wy);
45}
46
56template <class Scalar>
57constexpr Scalar harmonicMean(Scalar x, Scalar y, Scalar wx = 1.0, Scalar wy = 1.0) noexcept
58{
59 static_assert(Dune::IsNumber<Scalar>::value, "The arguments x, y, wx, and wy have to be numbers!");
60
61 if (x*y <= 0)
62 return 0;
63 return (wx + wy) * x * y / (wy * x + wx * y);
64}
65
74template <class Scalar>
75Scalar geometricMean(Scalar x, Scalar y) noexcept
76{
77 static_assert(Dune::IsNumber<Scalar>::value, "The arguments x and y have to be numbers!");
78
79 if (x*y <= 0)
80 return 0;
81 using std::sqrt;
82 return sqrt(x*y)*sign(x);
83}
84
96template <class Scalar, int m, int n>
97void harmonicMeanMatrix(Dune::FieldMatrix<Scalar, m, n> &K,
98 const Dune::FieldMatrix<Scalar, m, n> &Ki,
99 const Dune::FieldMatrix<Scalar, m, n> &Kj)
100{
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])) {
104 K[rowIdx][colIdx] =
105 harmonicMean(Ki[rowIdx][colIdx],
106 Kj[rowIdx][colIdx]);
107 }
108 else
109 K[rowIdx][colIdx] = Ki[rowIdx][colIdx];
110 }
111 }
112}
113
128template <class Scalar, class SolContainer>
129int invertLinearPolynomial(SolContainer &sol,
130 Scalar a,
131 Scalar b)
132{
133 if (a == 0.0)
134 return 0;
135
136 sol[0] = -b/a;
137 return 1;
138}
139
156template <class Scalar, class SolContainer>
157int invertQuadraticPolynomial(SolContainer &sol,
158 Scalar a,
159 Scalar b,
160 Scalar c)
161{
162 // check for a line
163 if (a == 0.0)
164 return invertLinearPolynomial(sol, b, c);
165
166 // discriminant
167 Scalar Delta = b*b - 4*a*c;
168 if (Delta < 0)
169 return 0; // no real roots
170
171 using std::sqrt;
172 Delta = sqrt(Delta);
173 sol[0] = (- b + Delta)/(2*a);
174 sol[1] = (- b - Delta)/(2*a);
175
176 // sort the result
177 if (sol[0] > sol[1])
178 {
179 using std::swap;
180 swap(sol[0], sol[1]);
181 }
182 return 2; // two real roots
183}
184
186template <class Scalar, class SolContainer>
187void invertCubicPolynomialPostProcess_(SolContainer &sol,
188 int numSol,
189 Scalar a, Scalar b, Scalar c, Scalar d,
190 std::size_t numIterations = 1)
191{
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); };
194
195 // do numIterations Newton iterations on the analytic solution
196 // and update result if the precision is increased
197 for (int i = 0; i < numSol; ++i)
198 {
199 Scalar x = sol[i];
200 Scalar fCurrent = eval(x);
201 const Scalar fOld = fCurrent;
202 for (int j = 0; j < numIterations; ++j)
203 {
204 const Scalar fPrime = evalDeriv(x);
205 if (fPrime == 0.0)
206 break;
207 x -= fCurrent/fPrime;
208 fCurrent = eval(x);
209 }
210
211 using std::abs;
212 if (abs(fCurrent) < abs(fOld))
213 sol[i] = x;
214 }
215}
217
239template <class Scalar, class SolContainer>
240int invertCubicPolynomial(SolContainer *sol,
241 Scalar a, Scalar b, Scalar c, Scalar d,
242 std::size_t numPostProcessIterations = 1)
243{
244 // reduces to a quadratic polynomial
245 if (a == 0)
246 return invertQuadraticPolynomial(sol, b, c, d);
247
248 // normalize the polynomial
249 b /= a;
250 c /= a;
251 d /= a;
252 a = 1;
253
254 // get rid of the quadratic term by substituting x = t - b/3
255 Scalar p = c - b*b/3;
256 Scalar q = d + (2*b*b*b - 9*b*c)/27;
257
258 // now we are at the form t^3 + p*t + q = 0. First we handle some
259 // special cases to avoid divisions by zero later...
260 if (p == 0.0 && q == 0.0) {
261 // t^3 = 0, i.e. triple root at t = 0
262 sol[0] = sol[1] = sol[2] = 0.0 - b/3;
263 return 3;
264 }
265 else if (p == 0.0 && q != 0.0) {
266 // t^3 + q = 0,
267 //
268 // i. e. single real root at t=curt(q)
269 using std::cbrt;
270 Scalar t = cbrt(q);
271 sol[0] = t - b/3;
272
273 return 1;
274 }
275 else if (p != 0.0 && q == 0.0) {
276 // t^3 + p*t = 0 = t*(t^2 + p),
277 //
278 // i. e. roots at t = 0, t^2 + p = 0
279 if (p > 0) {
280 sol[0] = 0.0 - b/3;
281 return 1; // only a single real root at t=0
282 }
283
284 // two additional real roots at t = sqrt(-p) and t = -sqrt(-p)
285 using std::sqrt;
286 sol[0] = -sqrt(-p) - b/3;
287 sol[1] = 0.0 - b/3;
288 sol[2] = sqrt(-p) - b/3;
289
290 return 3;
291 }
292
293 // At this point
294 //
295 // t^3 + p*t + q = 0
296 //
297 // with p != 0 and q != 0 holds. Introducing the variables u and v
298 // with the properties
299 //
300 // u + v = t and 3*u*v + p = 0
301 //
302 // leads to
303 //
304 // u^3 + v^3 + q = 0 .
305 //
306 // multiplying both sides with u^3 and taking advantage of the
307 // fact that u*v = -p/3 leads to
308 //
309 // u^6 + q*u^3 - p^3/27 = 0
310 //
311 // Now, substituting u^3 = w yields
312 //
313 // w^2 + q*w - p^3/27 = 0
314 //
315 // This is a quadratic equation with the solutions
316 //
317 // w = -q/2 + sqrt(q^2/4 + p^3/27) and
318 // w = -q/2 - sqrt(q^2/4 + p^3/27)
319 //
320 // Since w is equivalent to u^3 it is sufficient to only look at
321 // one of the two cases. Then, there are still 2 cases: positive
322 // and negative discriminant.
323 Scalar wDisc = q*q/4 + p*p*p/27;
324 if (wDisc >= 0) { // the positive discriminant case:
325 // calculate the cube root of - q/2 + sqrt(q^2/4 + p^3/27)
326 using std::cbrt;
327 using std::sqrt;
328
329 // Choose the root that is safe against the loss of precision
330 // that can cause wDisc to be equal to q*q/4 despite p != 0.
331 // We do not want u to be zero in that case. Mathematically,
332 // we are happy with either root.
333 const Scalar u = [&]{
334 return q > 0 ? cbrt(-0.5*q - sqrt(wDisc)) : cbrt(-0.5*q + sqrt(wDisc));
335 }();
336 // at this point, u != 0 since p^3 = 0 is necessary in order
337 // for u = 0 to hold, so
338 sol[0] = u - p/(3*u) - b/3;
339 // does not produce a division by zero. the remaining two
340 // roots of u are rotated by +- 2/3*pi in the complex plane
341 // and thus not considered here
342 invertCubicPolynomialPostProcess_(sol, 1, a, b, c, d, numPostProcessIterations);
343 return 1;
344 }
345 else { // the negative discriminant case:
346 Scalar uCubedRe = - q/2;
347 using std::sqrt;
348 Scalar uCubedIm = sqrt(-wDisc);
349 // calculate the cube root of - q/2 + sqrt(q^2/4 + p^3/27)
350 using std::cbrt;
351 Scalar uAbs = cbrt(sqrt(uCubedRe*uCubedRe + uCubedIm*uCubedIm));
352 using std::atan2;
353 Scalar phi = atan2(uCubedIm, uCubedRe)/3;
354
355 // calculate the length and the angle of the primitive root
356
357 // with the definitions from above it follows that
358 //
359 // x = u - p/(3*u) - b/3
360 //
361 // where x and u are complex numbers. Rewritten in polar form
362 // this is equivalent to
363 //
364 // x = |u|*e^(i*phi) - p*e^(-i*phi)/(3*|u|) - b/3 .
365 //
366 // Factoring out the e^ terms and subtracting the additional
367 // terms, yields
368 //
369 // x = (e^(i*phi) + e^(-i*phi))*(|u| - p/(3*|u|)) - y - b/3
370 //
371 // with
372 //
373 // y = - |u|*e^(-i*phi) + p*e^(i*phi)/(3*|u|) .
374 //
375 // The crucial observation is the fact that y is the conjugate
376 // of - x + b/3. This means that after taking advantage of the
377 // relation
378 //
379 // e^(i*phi) + e^(-i*phi) = 2*cos(phi)
380 //
381 // the equation
382 //
383 // x = 2*cos(phi)*(|u| - p / (3*|u|)) - conj(x) - 2*b/3
384 //
385 // holds. Since |u|, p, b and cos(phi) are real numbers, it
386 // follows that Im(x) = - Im(x) and thus Im(x) = 0. This
387 // implies
388 //
389 // Re(x) = x = cos(phi)*(|u| - p / (3*|u|)) - b/3 .
390 //
391 // Considering the fact that u is a cubic root, we have three
392 // values for phi which differ by 2/3*pi. This allows to
393 // calculate the three real roots of the polynomial:
394 for (int i = 0; i < 3; ++i) {
395 using std::cos;
396 sol[i] = cos(phi)*(uAbs - p/(3*uAbs)) - b/3;
397 phi += 2*M_PI/3;
398 }
399
400 // post process the obtained solution to increase numerical
401 // precision
402 invertCubicPolynomialPostProcess_(sol, 3, a, b, c, d, numPostProcessIterations);
403
404 // sort the result
405 using std::sort;
406 sort(sol, sol + 3);
407
408 return 3;
409 }
410
411 // NOT REACHABLE!
412 return 0;
413}
414
427template <class Scalar, int dim>
428bool isLarger(const Dune::FieldVector<Scalar, dim> &pos,
429 const Dune::FieldVector<Scalar, dim> &smallerVec)
430{
431 for (int i=0; i < dim; i++)
432 {
433 if (pos[i]<= smallerVec[i])
434 {
435 return false;
436 }
437 }
438 return true;
439}
440
453template <class Scalar, int dim>
454bool isSmaller(const Dune::FieldVector<Scalar, dim> &pos,
455 const Dune::FieldVector<Scalar, dim> &largerVec)
456{
457 for (int i=0; i < dim; i++)
458 {
459 if (pos[i]>= largerVec[i])
460 {
461 return false;
462 }
463 }
464 return true;
465}
466
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)
487{
488 if (isLarger(pos, smallerVec) && isSmaller(pos, largerVec))
489 {
490 return true;
491 }
492 else
493 return false;
494}
495
497namespace InterpolationPolicy { struct Linear; }
498
505template <class Policy = InterpolationPolicy::Linear, class Scalar, class ... Parameter>
506Scalar interpolate(Scalar ip, Parameter&& ... params)
507{ return Policy::interpolate(ip, std::forward<Parameter>(params) ...); }
508
513namespace InterpolationPolicy {
514
519struct Linear
520{
526 template<class Scalar>
527 static constexpr Scalar interpolate(Scalar ip, const std::array<Scalar, 2>& params)
528 {
529 return params[0]*(1.0 - ip) + params[1]*ip;
530 }
531};
532
538{
546 template<class Scalar, class RandomAccessContainer0, class RandomAccessContainer1>
547 static constexpr Scalar interpolate(Scalar ip, const RandomAccessContainer0& range, const RandomAccessContainer1& values)
548 {
549 // check bounds
550 if (ip > range.back()) return values.back();
551 if (ip < range[0]) return values[0];
552
553 // if we are within bounds find the index of the lower bound
554 const auto lookUpIndex = std::distance(range.begin(), std::lower_bound(range.begin(), range.end(), ip));
555 if (lookUpIndex == 0)
556 return values[0];
557
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]}});
560 }
561
562 template<class Scalar, class RandomAccessContainer>
563 static constexpr Scalar interpolate(Scalar ip, const std::pair<RandomAccessContainer, RandomAccessContainer>& table)
564 {
565 const auto& [range, values] = table;
566 return interpolate(ip, range, values);
567 }
568};
569
570} // end namespace InterpolationPolicy
571
581template <class Scalar>
582std::vector<Scalar> linspace(const Scalar begin, const Scalar end,
583 std::size_t samples,
584 bool endPoint = true)
585{
586 using std::max;
587 samples = max(std::size_t{2}, samples); // only makes sense for 2 or more 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;
593 return vec;
594}
595
596
609template <class Scalar>
610Scalar antoine(Scalar temperature,
611 Scalar A,
612 Scalar B,
613 Scalar C)
614{
615 const Scalar ln10 = 2.3025850929940459;
616 using std::exp;
617 return exp(ln10*(A - B/(C + temperature)));
618}
619
628template<class ValueType>
629constexpr int sign(const ValueType& value) noexcept
630{
631 return (ValueType(0) < value) - (value < ValueType(0));
632}
633
641template <class Scalar>
642Dune::FieldVector<Scalar, 3> crossProduct(const Dune::FieldVector<Scalar, 3> &vec1,
643 const Dune::FieldVector<Scalar, 3> &vec2)
644{
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]};
648}
649
657template <class Scalar>
658Scalar crossProduct(const Dune::FieldVector<Scalar, 2> &vec1,
659 const Dune::FieldVector<Scalar, 2> &vec2)
660{ return vec1[0]*vec2[1]-vec1[1]*vec2[0]; }
661
670template <class Scalar>
671Scalar tripleProduct(const Dune::FieldVector<Scalar, 3> &vec1,
672 const Dune::FieldVector<Scalar, 3> &vec2,
673 const Dune::FieldVector<Scalar, 3> &vec3)
674{ return crossProduct<Scalar>(vec1, vec2)*vec3; }
675
682template <class Scalar, int m, int n>
683Dune::FieldMatrix<Scalar, n, m> getTransposed(const Dune::FieldMatrix<Scalar, m, n>& M)
684{
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)
688 T[j][i] = M[i][j];
689
690 return T;
691}
692
699template <class Scalar>
700Dune::DynamicMatrix<Scalar> getTransposed(const Dune::DynamicMatrix<Scalar>& M)
701{
702 std::size_t rows_T = M.M();
703 std::size_t cols_T = M.N();
704
705 Dune::DynamicMatrix<Scalar> M_T(rows_T, cols_T, 0.0);
706
707 for (std::size_t i = 0; i < rows_T; ++i)
708 for (std::size_t j = 0; j < cols_T; ++j)
709 M_T[i][j] = M[j][i];
710
711 return M_T;
712}
713
721template <class Scalar>
722Dune::DynamicMatrix<Scalar> multiplyMatrices(const Dune::DynamicMatrix<Scalar> &M1,
723 const Dune::DynamicMatrix<Scalar> &M2)
724{
725 using size_type = typename Dune::DynamicMatrix<Scalar>::size_type;
726 const size_type rows = M1.N();
727 const size_type cols = M2.M();
728
729 DUNE_ASSERT_BOUNDS(M1.M() == M2.N());
730
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];
736
737 return result;
738}
739
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)
750{
751 using size_type = typename Dune::FieldMatrix<Scalar, rows1, cols2>::size_type;
752
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];
758
759 return result;
760}
761
762
769template <class MatrixType>
770typename Dune::DenseMatrix<MatrixType>::field_type
771trace(const Dune::DenseMatrix<MatrixType>& M)
772{
773 const auto rows = M.N();
774 DUNE_ASSERT_BOUNDS(rows == M.M()); // rows == cols
775
776 using MatType = Dune::DenseMatrix<MatrixType>;
777 typename MatType::field_type trace = 0.0;
778
779 for (typename MatType::size_type i = 0; i < rows; ++i)
780 trace += M[i][i];
781
782 return trace;
783}
784
798template <class MAT, class V>
799typename Dune::DenseVector<V>::derived_type
800mv(const Dune::DenseMatrix<MAT>& M,
801 const Dune::DenseVector<V>& v)
802{
803 typename Dune::DenseVector<V>::derived_type res(v);
804 M.mv(v, res);
805 return res;
806}
807
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)
829{
830 typename Dune::DenseVector<V>::derived_type res(v);
831 res *= m;
832 return res;
833}
834
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)
854{
855 return v1*mv(M, v2);
856}
857
877template <class V1, class FieldScalar, class V2>
878typename std::enable_if_t<Dune::IsNumber<FieldScalar>::value, FieldScalar>
879vtmv(const Dune::DenseVector<V1>& v1,
880 const FieldScalar m,
881 const Dune::DenseVector<V2>& v2)
882{
883 return m*(v1*v2);
884}
885
897template <class Scalar>
898std::array<Scalar, 2> linearRegression(const std::vector<Scalar>& x,
899 const std::vector<Scalar>& y)
900{
901 if (x.size() != y.size())
902 DUNE_THROW(Dune::InvalidStateException, "x and y array must have the same length.");
903
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();
906
907 // calculate temporary variables necessary for slope computation
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); }
911 );
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); }
915 );
916
917 // compute slope and intercept of the regression line
918 const Scalar slope = numerator / denominator;
919 const Scalar intercept = averageY - slope * averageX;
920
921 return {intercept, slope};
922}
923
924} // end namespace Dumux
925
926#endif
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
Definition: adapt.hh:17
interpolate linearly between two given values
Definition: math.hh:520
static constexpr Scalar interpolate(Scalar ip, const std::array< Scalar, 2 > &params)
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