GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/common/math.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 134 149 89.9%
Functions: 20 25 80.0%
Branches: 142 195 72.8%
Line Branch Exec Source
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 //
7 /*!
8 * \file
9 * \ingroup Core
10 * \brief Define some often used mathematical functions
11 */
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
26 namespace Dumux {
27
28 /*!
29 * \ingroup Core
30 * \brief Calculate the (weighted) arithmetic mean of two scalar values.
31 *
32 * \param x The first input value
33 * \param y The second input value
34 * \param wx The first weight
35 * \param wy The second weight
36 */
37 template <class Scalar>
38 constexpr 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
4/4
✓ Branch 0 taken 126549478 times.
✓ Branch 1 taken 608176 times.
✓ Branch 2 taken 111714016 times.
✓ Branch 3 taken 566378 times.
 238882052 if (x*y <= 0)
43 return 0;
44 238263494 return (x * wx + y * wy)/(wx + wy);
45 }
46
47 /*!
48 * \ingroup Core
49 * \brief Calculate the (weighted) harmonic mean of two scalar values.
50 *
51 * \param x The first input value
52 * \param y The second input value
53 * \param wx The first weight
54 * \param wy The second weight
55 */
56 template <class Scalar>
57 constexpr 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
13/18
✓ Branch 0 taken 1037293536 times.
✓ Branch 1 taken 2020558 times.
✓ Branch 2 taken 686178482 times.
✓ Branch 3 taken 26604 times.
✓ Branch 4 taken 185804488 times.
✓ Branch 5 taken 2020608 times.
✓ Branch 6 taken 278519324 times.
✓ Branch 7 taken 20000 times.
✓ Branch 8 taken 25489377 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 8598102 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 2424240 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 2435216 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 2424240 times.
✗ Branch 17 not taken.
 2233254775 if (x*y <= 0)
62 return 0;
63 2229167005 return (wx + wy) * x * y / (wy * x + wx * y);
64 }
65
66 /*!
67 * \ingroup Core
68 * \brief Calculate the geometric mean of two scalar values.
69 *
70 * \param x The first input value
71 * \param y The second input value
72 * \note as std::sqrt is not constexpr this function is not constexpr
73 */
74 template <class Scalar>
75 Scalar 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
85 /*!
86 * \ingroup Core
87 * \brief Calculate the harmonic mean of a fixed-size matrix.
88 *
89 * This is done by calculating the harmonic mean for each entry
90 * individually. The result is stored the first argument.
91 *
92 * \param K The matrix for storing the result
93 * \param Ki The first input matrix
94 * \param Kj The second input matrix
95 */
96 template <class Scalar, int m, int n>
97 void 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
114 /*!
115 * \ingroup Core
116 * \brief Invert a linear polynomial analytically
117 *
118 * The polynomial is defined as
119 * \f[ p(x) = a\; x + b \f]
120 *
121 * This method Returns the number of solutions which are in the real
122 * numbers, i.e. 1 except if the slope of the line is 0.
123 *
124 * \param sol Container into which the solutions are written
125 * \param a The coefficiont for the linear term
126 * \param b The coefficiont for the constant term
127 */
128 template <class Scalar, class SolContainer>
129 int invertLinearPolynomial(SolContainer &sol,
130 Scalar a,
131 Scalar b)
132 {
133 600000 if (a == 0.0)
134 return 0;
135
136 600000 sol[0] = -b/a;
137 return 1;
138 }
139
140 /*!
141 * \ingroup Core
142 * \brief Invert a quadratic polynomial analytically
143 *
144 * The polynomial is defined as
145 * \f[ p(x) = a\; x^2 + + b\;x + c \f]
146 *
147 * This method returns the number of solutions which are in the real
148 * numbers. The "sol" argument contains the real roots of the parabola
149 * in order with the smallest root first.
150 *
151 * \param sol Container into which the solutions are written
152 * \param a The coefficiont for the quadratic term
153 * \param b The coefficiont for the linear term
154 * \param c The coefficiont for the constant term
155 */
156 template <class Scalar, class SolContainer>
157 1200000 int invertQuadraticPolynomial(SolContainer &sol,
158 Scalar a,
159 Scalar b,
160 Scalar c)
161 {
162 // check for a line
163
2/2
✓ Branch 0 taken 600000 times.
✓ Branch 1 taken 600000 times.
 1200000 if (a == 0.0)
164
1/2
✓ Branch 0 taken 600000 times.
✗ Branch 1 not taken.
 600000 return invertLinearPolynomial(sol, b, c);
165
166 // discriminant
167 600000 Scalar Delta = b*b - 4*a*c;
168
1/2
✓ Branch 0 taken 600000 times.
✗ Branch 1 not taken.
 600000 if (Delta < 0)
169 return 0; // no real roots
170
171 using std::sqrt;
172 600000 Delta = sqrt(Delta);
173 600000 sol[0] = (- b + Delta)/(2*a);
174 600000 sol[1] = (- b - Delta)/(2*a);
175
176 // sort the result
177
2/2
✓ Branch 0 taken 400000 times.
✓ Branch 1 taken 200000 times.
 600000 if (sol[0] > sol[1])
178 {
179 using std::swap;
180 400000 swap(sol[0], sol[1]);
181 }
182 return 2; // two real roots
183 }
184
185 //! \cond false
186 template <class Scalar, class SolContainer>
187 3801593 void invertCubicPolynomialPostProcess_(SolContainer &sol,
188 int numSol,
189 Scalar a, Scalar b, Scalar c, Scalar d,
190 std::size_t numIterations = 1)
191 {
192 28525529 const auto eval = [&](auto x){ return d + x*(c + x*(b + x*a)); };
193 20213642 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
2/2
✓ Branch 0 taken 8311887 times.
✓ Branch 1 taken 3801593 times.
 12113480 for (int i = 0; i < numSol; ++i)
198 {
199 8311887 Scalar x = sol[i];
200 16623774 Scalar fCurrent = eval(x);
201 8311887 const Scalar fOld = fCurrent;
202
2/2
✓ Branch 0 taken 16412049 times.
✓ Branch 1 taken 8311887 times.
 24723936 for (int j = 0; j < numIterations; ++j)
203 {
204 32824098 const Scalar fPrime = evalDeriv(x);
205
1/2
✓ Branch 0 taken 16412049 times.
✗ Branch 1 not taken.
 16412049 if (fPrime == 0.0)
206 break;
207 16412049 x -= fCurrent/fPrime;
208 16412049 fCurrent = eval(x);
209 }
210
211 using std::abs;
212
6/6
✓ Branch 0 taken 5207745 times.
✓ Branch 1 taken 3104142 times.
✓ Branch 2 taken 5207745 times.
✓ Branch 3 taken 3104142 times.
✓ Branch 4 taken 5207745 times.
✓ Branch 5 taken 3104142 times.
 24935661 if (abs(fCurrent) < abs(fOld))
213 5207745 sol[i] = x;
214 }
215 3801593 }
216 //! \endcond
217
218 /*!
219 * \ingroup Core
220 * \brief Invert a cubic polynomial analytically
221 *
222 * The polynomial is defined as
223 * \f[ p(x) = a\; x^3 + + b\;x^3 + c\;x + d \f]
224 *
225 * This method returns the number of solutions which are in the real
226 * numbers. The "sol" argument contains the real roots of the cubic
227 * polynomial in order with the smallest root first.
228 *
229 * \note The closer the roots are to each other the less
230 * precise the inversion becomes. Increase number of post-processing iterations for improved results.
231 *
232 * \param sol Container into which the solutions are written
233 * \param a The coefficient for the cubic term
234 * \param b The coefficient for the quadratic term
235 * \param c The coefficient for the linear term
236 * \param d The coefficient for the constant term
237 * \param numPostProcessIterations The number of iterations to increase precision of the analytical result
238 */
239 template <class Scalar, class SolContainer>
240 3801593 int 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
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3801593 times.
 3801593 if (a == 0)
246 1200000 return invertQuadraticPolynomial(sol, b, c, d);
247
248 // normalize the polynomial
249 3801593 b /= a;
250 3801593 c /= a;
251 3801593 d /= a;
252 3801593 a = 1;
253
254 // get rid of the quadratic term by substituting x = t - b/3
255 3801593 Scalar p = c - b*b/3;
256 3801593 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
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 3801593 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 3801593 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
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 3801593 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 3801593 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
2/4
✓ Branch 0 taken 3801593 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 3801593 times.
 3801593 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 3801593 Scalar wDisc = q*q/4 + p*p*p/27;
324
2/2
✓ Branch 0 taken 1546446 times.
✓ Branch 1 taken 2255147 times.
 3801593 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 1546446 const Scalar u = [&]{
334
2/2
✓ Branch 0 taken 1504327 times.
✓ Branch 1 taken 42119 times.
 1546446 return q > 0 ? cbrt(-0.5*q - sqrt(wDisc)) : cbrt(-0.5*q + sqrt(wDisc));
335 1546446 }();
336 // at this point, u != 0 since p^3 = 0 is necessary in order
337 // for u = 0 to hold, so
338 1546446 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 1546446 invertCubicPolynomialPostProcess_(sol, 1, a, b, c, d, numPostProcessIterations);
343 1546446 return 1;
344 }
345 else { // the negative discriminant case:
346 2255147 Scalar uCubedRe = - q/2;
347 using std::sqrt;
348 2255147 Scalar uCubedIm = sqrt(-wDisc);
349 // calculate the cube root of - q/2 + sqrt(q^2/4 + p^3/27)
350 using std::cbrt;
351 2255147 Scalar uAbs = cbrt(sqrt(uCubedRe*uCubedRe + uCubedIm*uCubedIm));
352 using std::atan2;
353 2255147 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
2/2
✓ Branch 0 taken 6765441 times.
✓ Branch 1 taken 2255147 times.
 9020588 for (int i = 0; i < 3; ++i) {
395 using std::cos;
396 6765441 sol[i] = cos(phi)*(uAbs - p/(3*uAbs)) - b/3;
397 6765441 phi += 2*M_PI/3;
398 }
399
400 // post process the obtained solution to increase numerical
401 // precision
402 2255147 invertCubicPolynomialPostProcess_(sol, 3, a, b, c, d, numPostProcessIterations);
403
404 // sort the result
405 using std::sort;
406 4510294 sort(sol, sol + 3);
407
408 2255147 return 3;
409 }
410
411 // NOT REACHABLE!
412 return 0;
413 }
414
415 /*!
416 * \ingroup Core
417 * \brief Comparison of two position vectors
418 *
419 * Compares an current position vector with a reference vector, and returns true
420 * if the position vector is larger.
421 * "Larger" in this case means that all the entries of each spatial dimension are
422 * larger compared to the reference vector.
423 *
424 * \param pos Vector holding the current Position that is to be checked
425 * \param smallerVec Reference vector, holding the minimum values for comparison.
426 */
427 template <class Scalar, int dim>
428 bool 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
441 /*!
442 * \ingroup Core
443 * \brief Comparison of two position vectors
444 *
445 * Compares an current position vector with a reference vector, and returns true
446 * if the position vector is smaller.
447 * "Smaller" in this case means that all the entries of each spatial dimension are
448 * smaller in comparison with the reference vector.
449 *
450 * \param pos Vector holding the current Position that is to be checked
451 * \param largerVec Reference vector, holding the maximum values for comparison.
452 */
453 template <class Scalar, int dim>
454 bool 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
467 /*!
468 * \ingroup Core
469 * \brief Comparison of three position vectors
470 *
471 * Compares an current position vector with two reference vector, and returns true
472 * if the position vector lies in between them.
473 * "Between" in this case means that all the entries of each spatial dimension are
474 * smaller in comparison with the larger reference vector as well as larger compared
475 * to the smaller reference.
476 * This is comfortable to cheack weather the current position is located inside or
477 * outside of a lense with different properties.
478 *
479 * \param pos Vector holding the current Position that is to be checked
480 * \param smallerVec Reference vector, holding the minimum values for comparison.
481 * \param largerVec Reference vector, holding the maximum values for comparison.
482 */
483 template <class Scalar, int dim>
484 bool 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
496 //! forward declaration of the linear interpolation policy (default)
497 namespace InterpolationPolicy { struct Linear; }
498
499 /*!
500 * \ingroup Core
501 * \brief a generic function to interpolate given a set of parameters and an interpolation point
502 * \param params the parameters used for interpolation (depends on the policy used)
503 * \param ip the interpolation point
504 */
505 template <class Policy = InterpolationPolicy::Linear, class Scalar, class ... Parameter>
506 Scalar interpolate(Scalar ip, Parameter&& ... params)
507
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
 99345116 { return Policy::interpolate(ip, std::forward<Parameter>(params) ...); }
508
509 /*!
510 * \ingroup Core
511 * \brief Interpolation policies
512 */
513 namespace InterpolationPolicy {
514
515 /*!
516 * \ingroup Core
517 * \brief interpolate linearly between two given values
518 */
519 struct Linear
520 {
521 /*!
522 * \brief interpolate linearly between two given values
523 * \param ip the interpolation point in [0,1]
524 * \param params array with the lower and upper bound
525 */
526 template<class Scalar>
527 static constexpr Scalar interpolate(Scalar ip, const std::array<Scalar, 2>& params)
528 {
529 94626512 return params[0]*(1.0 - ip) + params[1]*ip;
530 }
531 };
532
533 /*!
534 * \ingroup Core
535 * \brief interpolate linearly in a piecewise linear function (tabularized function)
536 */
537 struct LinearTable
538 {
539 /*!
540 * \brief interpolate linearly in a piecewise linear function (tabularized function)
541 * \param ip the interpolation point
542 * \param range positions of values
543 * \param values values to interpolate from
544 * \note if the interpolation point is out of bounds this will return the bounds
545 */
546 template<class Scalar, class RandomAccessContainer0, class RandomAccessContainer1>
547 99347470 static constexpr Scalar interpolate(Scalar ip, const RandomAccessContainer0& range, const RandomAccessContainer1& values)
548 {
549 // check bounds
550
4/4
✓ Branch 0 taken 4718223 times.
✓ Branch 1 taken 94629247 times.
✓ Branch 2 taken 4718223 times.
✓ Branch 3 taken 94629247 times.
 198694940 if (ip > range.back()) return values.back();
551
2/2
✓ Branch 0 taken 2734 times.
✓ Branch 1 taken 94626513 times.
 94629247 if (ip < range[0]) return values[0];
552
553 // if we are within bounds find the index of the lower bound
554
4/4
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 94626512 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 94626512 times.
 473132565 const auto lookUpIndex = std::distance(range.begin(), std::lower_bound(range.begin(), range.end(), ip));
555
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 94626512 times.
 94626513 if (lookUpIndex == 0)
556 1 return values[0];
557
558 189253024 const auto ipLinear = (ip - range[lookUpIndex-1])/(range[lookUpIndex] - range[lookUpIndex-1]);
559 378506048 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 4692 const auto& [range, values] = table;
566 2346 return interpolate(ip, range, values);
567 }
568 };
569
570 } // end namespace InterpolationPolicy
571
572 /*!
573 * \ingroup Core
574 * \brief Generates linearly spaced vectors
575 *
576 * \param begin The first value in the vector
577 * \param end The last value in the vector
578 * \param samples The size of the vector
579 * \param endPoint if the range is including the interval's end point or not
580 */
581 template <class Scalar>
582 50 std::vector<Scalar> linspace(const Scalar begin, const Scalar end,
583 std::size_t samples,
584 bool endPoint = true)
585 {
586 using std::max;
587
1/2
✓ Branch 0 taken 50 times.
✗ Branch 1 not taken.
 50 samples = max(std::size_t{2}, samples); // only makes sense for 2 or more samples
588
2/2
✓ Branch 0 taken 49 times.
✓ Branch 1 taken 1 times.
 50 const Scalar divisor = endPoint ? samples-1 : samples;
589 50 const Scalar delta = (end-begin)/divisor;
590 100 std::vector<Scalar> vec(samples);
591
2/2
✓ Branch 0 taken 8048755 times.
✓ Branch 1 taken 50 times.
 8048805 for (std::size_t i = 0; i < samples; ++i)
592 16097510 vec[i] = begin + i*delta;
593 50 return vec;
594 }
595
596
597 /*!
598 * \ingroup Core
599 * \brief Evaluates the Antoine equation used to calculate the vapour
600 * pressure of various liquids.
601 *
602 * See http://en.wikipedia.org/wiki/Antoine_equation
603 *
604 * \param temperature The temperature [K] of the fluid
605 * \param A The first coefficient for the Antoine equation
606 * \param B The first coefficient for the Antoine equation
607 * \param C The first coefficient for the Antoine equation
608 */
609 template <class Scalar>
610 Scalar 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
620 /*!
621 * \ingroup Core
622 * \brief Sign or signum function.
623 *
624 * Returns 1 for a positive argument.
625 * Returns -1 for a negative argument.
626 * Returns 0 if the argument is zero.
627 */
628 template<class ValueType>
629 ✗ constexpr int sign(const ValueType& value) noexcept
630 {
631
13/13
✓ Branch 0 taken 27 times.
✓ Branch 1 taken 159884759 times.
✓ Branch 2 taken 101242334 times.
✓ Branch 3 taken 45830437 times.
✓ Branch 4 taken 124074537 times.
✓ Branch 5 taken 69539 times.
✓ Branch 6 taken 148002987 times.
✓ Branch 7 taken 252988 times.
✓ Branch 8 taken 57861382 times.
✓ Branch 9 taken 1152762 times.
✓ Branch 10 taken 1835 times.
✓ Branch 11 taken 1024472 times.
✓ Branch 12 taken 15 times.
 1177269101 return (ValueType(0) < value) - (value < ValueType(0));
632 }
633
634 /*!
635 * \ingroup Core
636 * \brief Cross product of two vectors in three-dimensional Euclidean space
637 *
638 * \param vec1 The first vector
639 * \param vec2 The second vector
640 */
641 template <class Scalar>
642 39365508 Dune::FieldVector<Scalar, 3> crossProduct(const Dune::FieldVector<Scalar, 3> &vec1,
643 const Dune::FieldVector<Scalar, 3> &vec2)
644 {
645 157462032 return {vec1[1]*vec2[2]-vec1[2]*vec2[1],
646 118096524 vec1[2]*vec2[0]-vec1[0]*vec2[2],
647 236193048 vec1[0]*vec2[1]-vec1[1]*vec2[0]};
648 }
649
650 /*!
651 * \ingroup Core
652 * \brief Cross product of two vectors in two-dimensional Euclidean space retuning scalar
653 *
654 * \param vec1 The first vector
655 * \param vec2 The second vector
656 */
657 template <class Scalar>
658 Scalar crossProduct(const Dune::FieldVector<Scalar, 2> &vec1,
659 const Dune::FieldVector<Scalar, 2> &vec2)
660
20/20
✓ Branch 0 taken 13628205 times.
✓ Branch 1 taken 21035116 times.
✓ Branch 2 taken 13628205 times.
✓ Branch 3 taken 21035116 times.
✓ Branch 4 taken 13628205 times.
✓ Branch 5 taken 21035116 times.
✓ Branch 6 taken 13628205 times.
✓ Branch 7 taken 21035116 times.
✓ Branch 8 taken 13628205 times.
✓ Branch 9 taken 21035116 times.
✓ Branch 10 taken 200 times.
✓ Branch 11 taken 200 times.
✓ Branch 12 taken 200 times.
✓ Branch 13 taken 200 times.
✓ Branch 14 taken 200 times.
✓ Branch 15 taken 200 times.
✓ Branch 16 taken 200 times.
✓ Branch 17 taken 200 times.
✓ Branch 18 taken 200 times.
✓ Branch 19 taken 200 times.
 363630790 { return vec1[0]*vec2[1]-vec1[1]*vec2[0]; }
661
662 /*!
663 * \ingroup Core
664 * \brief Triple product of three vectors in three-dimensional Euclidean space retuning scalar
665 *
666 * \param vec1 The first vector
667 * \param vec2 The second vector
668 * \param vec3 The third vector
669 */
670 template <class Scalar>
671 Scalar tripleProduct(const Dune::FieldVector<Scalar, 3> &vec1,
672 const Dune::FieldVector<Scalar, 3> &vec2,
673 const Dune::FieldVector<Scalar, 3> &vec3)
674 10272607 { return crossProduct<Scalar>(vec1, vec2)*vec3; }
675
676 /*!
677 * \ingroup Core
678 * \brief Transpose a FieldMatrix
679 *
680 * \param M The matrix to be transposed
681 */
682 template <class Scalar, int m, int n>
683 Dune::FieldMatrix<Scalar, n, m> getTransposed(const Dune::FieldMatrix<Scalar, m, n>& M)
684 {
685 2185722 Dune::FieldMatrix<Scalar, n, m> T;
686
6/6
✓ Branch 0 taken 1753604 times.
✓ Branch 1 taken 876802 times.
✓ Branch 2 taken 38758 times.
✓ Branch 3 taken 12920 times.
✓ Branch 4 taken 2592000 times.
✓ Branch 5 taken 1296000 times.
 6570084 for (std::size_t i = 0; i < m; ++i)
687
6/6
✓ Branch 0 taken 3507208 times.
✓ Branch 1 taken 1753604 times.
✓ Branch 2 taken 116270 times.
✓ Branch 3 taken 38758 times.
✓ Branch 4 taken 5184000 times.
✓ Branch 5 taken 2592000 times.
 13191840 for (std::size_t j = 0; j < n; ++j)
688 44037390 T[j][i] = M[i][j];
689
690 return T;
691 }
692
693 /*!
694 * \ingroup Core
695 * \brief Transpose a DynamicMatrix
696 *
697 * \param M The matrix to be transposed
698 */
699 template <class Scalar>
700 Dune::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
714 /*!
715 * \ingroup Core
716 * \brief Multiply two dynamic matrices
717 *
718 * \param M1 The first dynamic matrix
719 * \param M2 The second dynamic matrix (to be multiplied to M1 from the right side)
720 */
721 template <class Scalar>
722 52461048 Dune::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 52461048 const size_type rows = M1.N();
727 52461048 const size_type cols = M2.M();
728
729 DUNE_ASSERT_BOUNDS(M1.M() == M2.N());
730
731 52461048 Dune::DynamicMatrix<Scalar> result(rows, cols, 0.0);
732
2/2
✓ Branch 1 taken 207548397 times.
✓ Branch 2 taken 52461048 times.
 260009445 for (size_type i = 0; i < rows; i++)
733
2/2
✓ Branch 0 taken 822834623 times.
✓ Branch 1 taken 207548397 times.
 1030383020 for (size_type j = 0; j < cols; j++)
734
2/2
✓ Branch 1 taken 3395029373 times.
✓ Branch 2 taken 822834623 times.
 4217863996 for (size_type k = 0; k < M1.M(); k++)
735 23765205611 result[i][j] += M1[i][k]*M2[k][j];
736
737 52461048 return result;
738 }
739
740 /*!
741 * \ingroup Core
742 * \brief Multiply two field matrices
743 *
744 * \param M1 The first field matrix
745 * \param M2 The second field matrix (to be multiplied to M1 from the right side)
746 */
747 template <class Scalar, int rows1, int cols1, int cols2>
748 706929 Dune::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 706929 Dune::FieldMatrix<Scalar, rows1, cols2> result(0.0);
754
2/2
✓ Branch 0 taken 2827716 times.
✓ Branch 1 taken 706929 times.
 3534645 for (size_type i = 0; i < rows1; i++)
755
2/2
✓ Branch 0 taken 11310864 times.
✓ Branch 1 taken 2827716 times.
 14138580 for (size_type j = 0; j < cols2; j++)
756
2/2
✓ Branch 0 taken 45243456 times.
✓ Branch 1 taken 11310864 times.
 56554320 for (size_type k = 0; k < cols1; k++)
757 316704192 result[i][j] += M1[i][k]*M2[k][j];
758
759 706929 return result;
760 }
761
762
763 /*!
764 * \ingroup Core
765 * \brief Trace of a dense matrix
766 *
767 * \param M The dense matrix
768 */
769 template <class MatrixType>
770 typename Dune::DenseMatrix<MatrixType>::field_type
771 trace(const Dune::DenseMatrix<MatrixType>& M)
772 {
773 1 const auto rows = M.N();
774 DUNE_ASSERT_BOUNDS(rows == M.M()); // rows == cols
775
776 using MatType = Dune::DenseMatrix<MatrixType>;
777 1 typename MatType::field_type trace = 0.0;
778
779
4/4
✓ Branch 0 taken 1204579 times.
✓ Branch 1 taken 3585729 times.
✓ Branch 2 taken 651 times.
✓ Branch 3 taken 1297 times.
 4792255 for (typename MatType::size_type i = 0; i < rows; ++i)
780 10761090 trace += M[i][i];
781
782 return trace;
783 }
784
785 /*!
786 * \ingroup Core
787 * \brief Returns the result of the projection of
788 * a vector v with a Matrix M.
789 *
790 * Note: We use DenseVector and DenseMatrix here so that
791 * it can be used with the statically and dynamically
792 * allocated Dune Vectors/Matrices. Size mismatch
793 * assertions are done in the respective Dune classes.
794 *
795 * \param M The matrix
796 * \param v The vector
797 */
798 template <class MAT, class V>
799 typename Dune::DenseVector<V>::derived_type
800 88756348 mv(const Dune::DenseMatrix<MAT>& M,
801 const Dune::DenseVector<V>& v)
802 {
803
2/4
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
 97215157 typename Dune::DenseVector<V>::derived_type res(v);
804 97215150 M.mv(v, res);
805 88756348 return res;
806 }
807
808 /*!
809 * \ingroup Core
810 * \brief Returns the result of a vector v multiplied by a scalar m.
811 *
812 * Note: We use DenseVector and DenseMatrix here so that
813 * it can be used with the statically and dynamically
814 * allocated Dune Vectors/Matrices. Size mismatch
815 * assertions are done in the respective Dune classes.
816 *
817 * \note We need the enable_if to make sure that only Scalars
818 * fit here. Matrix types are forced to use the above
819 * mv(DenseMatrix, DenseVector) instead.
820 *
821 * \param m The scale factor
822 * \param v The vector
823 */
824
825 template <class FieldScalar, class V>
826 typename std::enable_if_t<Dune::IsNumber<FieldScalar>::value,
827 typename Dune::DenseVector<V>::derived_type>
828 1764 mv(const FieldScalar m, const Dune::DenseVector<V>& v)
829 {
830 1766 typename Dune::DenseVector<V>::derived_type res(v);
831 1764 res *= m;
832 1764 return res;
833 }
834
835 /*!
836 * \ingroup Core
837 * \brief Evaluates the scalar product of a vector v2, projected by
838 * a matrix M, with a vector v1.
839 *
840 * Note: We use DenseVector and DenseMatrix here so that
841 * it can be used with the statically and dynamically
842 * allocated Dune Vectors/Matrices. Size mismatch
843 * assertions are done in the respective Dune classes.
844 *
845 * \param v1 The first vector
846 * \param M The matrix
847 * \param v2 The second vector
848 */
849 template <class V1, class MAT, class V2>
850 typename Dune::DenseMatrix<MAT>::value_type
851 6 vtmv(const Dune::DenseVector<V1>& v1,
852 const Dune::DenseMatrix<MAT>& M,
853 const Dune::DenseVector<V2>& v2)
854 {
855
5/10
✓ Branch 1 taken 3914614 times.
✓ Branch 2 taken 82758 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 1 times.
 103970622 return v1*mv(M, v2);
856 }
857
858 /*!
859 * \ingroup Core
860 * \brief Evaluates the scalar product of a vector v2, scaled by
861 * a scalar m, with a vector v1.
862 *
863 * Note: We use DenseVector and DenseMatrix here so that
864 * it can be used with the statically and dynamically
865 * allocated Dune Vectors/Matrices. Size mismatch
866 * assertions are done in the respective Dune classes.
867 *
868 * \note We need the enable_if to make sure that only Scalars
869 * fit here. Matrix types are forced to use the above
870 * vtmv(DenseVector, DenseMatrix, DenseVector) instead.
871 *
872 * \param v1 The first vector
873 * \param m The scale factor
874 * \param v2 The second vector
875 */
876
877 template <class V1, class FieldScalar, class V2>
878 typename std::enable_if_t<Dune::IsNumber<FieldScalar>::value, FieldScalar>
879 ✗ vtmv(const Dune::DenseVector<V1>& v1,
880 const FieldScalar m,
881 const Dune::DenseVector<V2>& v2)
882 {
883
8/12
✓ Branch 0 taken 474041684 times.
✓ Branch 1 taken 7891822 times.
✓ Branch 2 taken 278150781 times.
✓ Branch 3 taken 1679216 times.
✓ Branch 4 taken 875901 times.
✓ Branch 5 taken 1083504 times.
✓ Branch 6 taken 770556 times.
✓ Branch 7 taken 49896 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 3023799312 return m*(v1*v2);
884 }
885
886 /*!
887 * \ingroup Core
888 * \brief Returns the [intercept, slope] of the regression line
889 * fitted to a set of (x, y) data points.
890 *
891 * Note: We use least-square regression method to find
892 * the regression line.
893 *
894 * \param x x-values of the data set
895 * \param y y-values of the data set
896 */
897 template <class Scalar>
898 3 std::array<Scalar, 2> linearRegression(const std::vector<Scalar>& x,
899 const std::vector<Scalar>& y)
900 {
901
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 3 times.
 9 if (x.size() != y.size())
902 ✗ DUNE_THROW(Dune::InvalidStateException, "x and y array must have the same length.");
903
904 3 const Scalar averageX = std::accumulate(x.begin(), x.end(), 0.0)/x.size();
905 9 const Scalar averageY = std::accumulate(y.begin(), y.end(), 0.0)/y.size();
906
907 // calculate temporary variables necessary for slope computation
908 12 const Scalar numerator = std::inner_product(
909 x.begin(), x.end(), y.begin(), 0.0, std::plus<Scalar>(),
910 208 [&](auto xx, auto yy) { return (xx - averageX) * (yy - averageY); }
911 );
912 3 const Scalar denominator = std::inner_product(
913 x.begin(), x.end(), x.begin(), 0.0, std::plus<Scalar>(),
914 208 [&](auto xx, auto yy) { return (xx - averageX) * (yy - averageX); }
915 );
916
917 // compute slope and intercept of the regression line
918 3 const Scalar slope = numerator / denominator;
919 3 const Scalar intercept = averageY - slope * averageX;
920
921 3 return {intercept, slope};
922 }
923
924 } // end namespace Dumux
925
926 #endif
927