math.hh 16 KB
Newer Older
1
2
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
3
/*****************************************************************************
4
 *   See the file COPYING for full copying permissions.                      *
5
 *                                                                           *
6
 *   This program is free software: you can redistribute it and/or modify    *
7
 *   it under the terms of the GNU General Public License as published by    *
8
9
 *   the Free Software Foundation, either version 2 of the License, or       *
 *   (at your option) any later version.                                     *
10
 *                                                                           *
11
12
 *   This program is distributed in the hope that it will be useful,         *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
13
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
14
15
16
17
 *   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/>.   *
18
19
20
21
22
23
24
25
 *****************************************************************************/
/*!
 * \file
 * \brief Define some often used mathematical functions
 */
#ifndef DUMUX_MATH_HH
#define DUMUX_MATH_HH

Andreas Lauser's avatar
Andreas Lauser committed
26
27
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
28

29
#include <cmath>
30
#include <algorithm>
31
32
33
34

namespace Dumux
{
/*!
35
 * \ingroup Math
36
 * \brief Calculate the harmonic mean of two scalar values.
37
38
39
 *
 * \param x The first input value
 * \param y The second input value
40
41
42
43
44
45
46
47
48
49
 */
template <class Scalar>
Scalar harmonicMean(Scalar x, Scalar y)
{
    if (x*y <= 0)
        return 0;
    return (2*x*y)/(x + y);
}

/*!
50
 * \ingroup Math
51
 * \brief Calculate the geometric mean of two scalar values.
52
53
54
 *
 * \param x The first input value
 * \param y The second input value
55
56
57
58
59
60
61
62
63
64
 */
template <class Scalar>
Scalar geometricMean(Scalar x, Scalar y)
{
    if (x*y <= 0)
        return 0;
    return std::sqrt(x*y)*((x < 0)?-1:1);
}

/*!
65
 * \ingroup Math
66
67
68
69
 * \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.
70
71
72
73
 *
 * \param K The matrix for storing the result
 * \param Ki The first input matrix
 * \param Kj The second input matrix
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
 */
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];
        }
    }
}

93
/*!
94
 * \ingroup Math
95
96
 * \brief Invert a linear polynomial analytically
 *
97
98
99
100
101
102
103
104
105
 * 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 coefficiont for the linear term
 * \param b The coefficiont for the constant term
106
107
 */
template <class Scalar, class SolContainer>
108
int invertLinearPolynomial(SolContainer &sol,
109
110
111
112
113
114
115
116
117
118
119
                           Scalar a,
                           Scalar b)
{
    if (a == 0.0)
        return 0;

    sol[0] = -b/a;
    return 1;
}

/*!
120
 * \ingroup Math
121
122
 * \brief Invert a quadratic polynomial analytically
 *
123
124
125
 * The polynomial is defined as
 * \f[ p(x) = a\; x^2 + + b\;x + c \f]
 *
126
 * This method returns the number of solutions which are in the real
127
128
129
130
131
132
133
 * 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 coefficiont for the quadratic term
 * \param b The coefficiont for the linear term
 * \param c The coefficiont for the constant term
134
135
 */
template <class Scalar, class SolContainer>
136
int invertQuadraticPolynomial(SolContainer &sol,
137
138
139
140
141
142
143
                              Scalar a,
                              Scalar b,
                              Scalar c)
{
    // check for a line
    if (a == 0.0)
        return invertLinearPolynomial(sol, b, c);
144

145
146
147
148
149
150
151
152
    // 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);
153
154
155
156

    // sort the result
    if (sol[0] > sol[1])
        std::swap(sol[0], sol[1]);
157
158
159
    return 2; // two real roots
}

Thomas Fetzer's avatar
Thomas Fetzer committed
160
//! \cond false
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
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

187
/*!
188
 * \ingroup Math
189
190
 * \brief Invert a cubic polynomial analytically
 *
191
192
193
194
195
196
197
198
199
200
201
202
 * 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 coefficiont for the cubic term
 * \param b The coefficiont for the quadratic term
 * \param c The coefficiont for the linear term
 * \param d The coefficiont for the constant term
203
204
 */
template <class Scalar, class SolContainer>
205
int invertCubicPolynomial(SolContainer *sol,
206
207
                          Scalar a,
                          Scalar b,
208
                          Scalar c,
209
210
211
                          Scalar d)
{
    // reduces to a quadratic polynomial
212
    if (a == 0)
213
        return invertQuadraticPolynomial(sol, b, c, d);
214

215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
    // 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;

    // now we are at the form t^3 + p*t + q = 0. First we handle some
    // special cases to avoid divisions by zero later...
    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) {
233
        // t^3 + q = 0,
234
235
236
237
238
239
        //
        // 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;
240

241
242
243
        return 1;
    }
    else if (p != 0.0 && q == 0.0) {
244
        // t^3 + p*t = 0 = t*(t^2 + p),
245
        //
246
        // i. e. roots at t = 0, t^2 + p = 0
247
248
        if (p > 0) {
            sol[0] = 0.0 - b/3;
249
            return 1; // only a single real root at t=0
250
        }
251

252
        // two additional real roots at t = sqrt(-p) and t = -sqrt(-p)
253
        sol[0] = -std::sqrt(-p) - b/3;
254
255
        sol[1] = 0.0 - b/3;
        sol[2] = std::sqrt(-p) - b/3;
256

257
258
259
260
261
262
263
264
        return 3;
    }

    // At this point
    //
    // t^3 + p*t + q = 0
    //
    // with p != 0 and q != 0 holds. Introducing the variables u and v
265
    // with the properties
266
    //
267
    //   u + v = t       and       3*u*v + p = 0
268
    //
269
    // leads to
270
271
272
273
274
275
276
277
278
279
    //
    // 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
    //
280
    // w^2 + q*w - p^3/27 = 0
281
    //
282
    // This is a quadratic equation with the solutions
283
    //
284
285
    // w = -q/2 + sqrt(q^2/4 + p^3/27) and
    // w = -q/2 - sqrt(q^2/4 + p^3/27)
286
    //
287
288
289
    // 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.
290
    Scalar wDisc = q*q/4 + p*p*p/27;
291
292
    if (wDisc >= 0) { // the positive discriminant case:
        // calculate the cube root of - q/2 + sqrt(q^2/4 + p^3/27)
293
294
295
296
297
        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
298
        // for u = 0 to hold, so
299
        sol[0] = u - p/(3*u) - b/3;
300
301
302
        // 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
303
        invertCubicPolynomialPostProcess_(sol, 1, a, b, c, d);
304
305
        return 1;
    }
306
307
    else { // the negative discriminant case:
        Scalar uCubedRe = - q/2;
308
        Scalar uCubedIm = std::sqrt(-wDisc);
309
310
        // 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);
311
        Scalar phi = std::atan2(uCubedIm, uCubedRe)/3;
312
313

        // calculate the length and the angle of the primitive root
314
315

        // with the definitions from above it follows that
316
317
318
        //
        // x = u - p/(3*u) - b/3
        //
319
        // where x and u are complex numbers. Rewritten in polar form
320
321
        // this is equivalent to
        //
322
        // x = |u|*e^(i*phi) - p*e^(-i*phi)/(3*|u|) - b/3 .
323
        //
324
325
        // Factoring out the e^ terms and subtracting the additional
        // terms, yields
326
        //
327
        // x = (e^(i*phi) + e^(-i*phi))*(|u| - p/(3*|u|)) - y - b/3
328
329
        //
        // with
330
        //
331
        // y = - |u|*e^(-i*phi) + p*e^(i*phi)/(3*|u|) .
332
        //
333
334
335
        // The crucial observation is the fact that y is the conjugate
        // of - x + b/3. This means that after taking advantage of the
        // relation
336
        //
337
        // e^(i*phi) + e^(-i*phi) = 2*cos(phi)
338
        //
339
        // the equation
340
        //
341
        // x = 2*cos(phi)*(|u| - p / (3*|u|)) - conj(x) - 2*b/3
342
        //
343
344
345
        // 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
346
        //
347
        // Re(x) = x = cos(phi)*(|u| - p / (3*|u|)) - b/3 .
348
349
350
351
352
        //
        // 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) {
353
            sol[i] = std::cos(phi)*(uAbs - p/(3*uAbs)) - b/3;
354
355
            phi += 2*M_PI/3;
        }
356

357
358
359
360
        // post process the obtained solution to increase numerical
        // precision
        invertCubicPolynomialPostProcess_(sol, 3, a, b, c, d);

361
362
        // sort the result
        std::sort(sol, sol + 3);
363

364
365
366
367
368
369
        return 3;
    }

    // NOT REACHABLE!
    return 0;
}
370

371
/*!
372
 * \ingroup Math
373
374
375
376
377
378
379
 * \brief Comparison of two position vectors
 *
 * Compares an current position vector with a reference vector, and returns true
 * if the position vector is larger.
 * "Larger" in this case means that all the entries of each spacial dimension are
 * larger compared to the reference vector.
 *
380
381
 * \param pos Vektor holding the current Position that is to be checked
 * \param smallerVec Reference vector, holding the minimum values for comparison.
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
 */
template <class Scalar, int dim>
bool isLarger(const Dune::FieldVector<Scalar, dim> &pos,
              const Dune::FieldVector<Scalar, dim> &smallerVec)
{
    for (int i=0; i < dim; i++)
    {
        if (pos[i]<= smallerVec[i])
        {
            return false;
        }
    }
    return true;
}

/*!
398
 * \ingroup Math
399
400
401
402
403
404
405
 * \brief Comparison of two position vectors
 *
 * Compares an current position vector with a reference vector, and returns true
 * if the position vector is smaller.
 * "Smaller" in this case means that all the entries of each spacial dimension are
 * smaller in comparison with the reference vector.
 *
406
407
 * \param pos Vektor holding the current Position that is to be checked
 * \param largerVec Reference vector, holding the maximum values for comparison.
408
409
410
 */
template <class Scalar, int dim>
bool isSmaller(const Dune::FieldVector<Scalar, dim> &pos,
Andreas Lauser's avatar
Andreas Lauser committed
411
               const Dune::FieldVector<Scalar, dim> &largerVec)
412
413
414
415
416
417
418
419
420
421
422
423
{
    for (int i=0; i < dim; i++)
    {
        if (pos[i]>= largerVec[i])
        {
            return false;
        }
    }
    return true;
}

/*!
424
 * \ingroup Math
425
426
427
428
429
430
431
432
433
434
 * \brief Comparison of three position vectors
 *
 * Compares an current position vector with two reference vector, and returns true
 * if the position vector lies in between them.
 * "Between" in this case means that all the entries of each spacial dimension are
 * smaller in comparison with the larger reference vector as well as larger campared
 * to the smaller reference.
 * This is comfortable to cheack weather the current position is located inside or
 * outside of a lense with different properties.
 *
435
436
437
 * \param pos Vektor holding the current Position that is to be checked
 * \param smallerVec Reference vector, holding the minimum values for comparison.
 * \param largerVec Reference vector, holding the maximum values for comparison.
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
 */
template <class Scalar, int dim>
bool isBetween(const Dune::FieldVector<Scalar, dim> &pos,
              const Dune::FieldVector<Scalar, dim> &smallerVec,
              const Dune::FieldVector<Scalar, dim> &largerVec)
{
   if (isLarger(pos, smallerVec) && isSmaller(pos, largerVec))
       {
           return true;
       }
   else
       return false;
}


/*!
454
 * \ingroup Math
455
456
457
458
 * \brief Evaluates the Antoine equation used to calculate the vapour
 *        pressure of various liquids.
 *
 * See http://en.wikipedia.org/wiki/Antoine_equation
459
460
461
462
463
 *
 * \param temperature The temperature [K] of the fluid
 * \param A The first coefficient for the Antoine equation
 * \param B The first coefficient for the Antoine equation
 * \param C The first coefficient for the Antoine equation
464
465
466
467
468
469
470
471
472
473
474
 */
template <class Scalar>
Scalar antoine(Scalar temperature,
               Scalar A,
               Scalar B,
               Scalar C)
{
    const Scalar ln10 = 2.3025850929940459;
    return std::exp(ln10*(A - B/(C + temperature)));
}

475
/*!
476
 * \brief Cross product of two vectors in three-dimensional Euclidean space
477
478
479
480
 *
 * \param vec1 The first vector
 * \param vec2 The second vector
 */
481
482
483
template <class Scalar>
Dune::FieldVector<Scalar, 3> crossProduct(const Dune::FieldVector<Scalar, 3> &vec1,
                                          const Dune::FieldVector<Scalar, 3> &vec2)
484
{
485
486
487
    return {vec1[1]*vec2[2]-vec1[2]*vec2[1],
            vec1[2]*vec2[0]-vec1[0]*vec2[2],
            vec1[0]*vec2[1]-vec1[1]*vec2[0]};
488
489
}

490
491
492
493
494
495
496
497
498
499
/*!
 * \brief Cross product of two vectors in two-dimensional Euclidean space retuning scalar
 *
 * \param vec1 The first vector
 * \param vec2 The second vector
 */
template <class Scalar>
Scalar crossProduct(const Dune::FieldVector<Scalar, 2> &vec1,
                    const Dune::FieldVector<Scalar, 2> &vec2)
{   return vec1[0]*vec2[1]-vec1[1]*vec2[0]; }
500

Dennis Gläser's avatar
Dennis Gläser committed
501
502
503
504
505
506
507
508
509
510
511
512
/*!
 * \brief Triple product of three vectors in three-dimensional Euclidean space retuning scalar
 *
 * \param vec1 The first vector
 * \param vec2 The second vector
 * \param vec3 The third vector
 */
template <class Scalar>
Scalar tripleProduct(const Dune::FieldVector<Scalar, 3> &vec1,
                    const Dune::FieldVector<Scalar, 3> &vec2,
                    const Dune::FieldVector<Scalar, 3> &vec3)
{   return crossProduct<Scalar>(vec1, vec2)*vec3; }
513
} // end namespace Dumux
514
515

#endif