Commit 849cb82f authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

Merge branch 'feature/use-ADL-for-std-functions' into 'master'

Feature/use adl for std functions

* `std::swap` and `dumux/common` was cleaned up
* much more needs to be done
* I found a bunch of other clean-up stuff

See merge request !259
parent 56f43044
......@@ -4,7 +4,6 @@ set(modules
DumuxDoxygen.cmake
DumuxMacros.cmake
DumuxTestMacros.cmake
FindGstat.cmake
FindValgrind.cmake)
FindGstat.cmake)
include(GNUInstallDirs)
install(FILES ${modules} DESTINATION ${DUNE_INSTALL_MODULEDIR})
#
# Module that checks whether Valgrind's header memcheck.h is present
#
# Variables used by this module which you may want to set:
# VALGRIND_ROOT Path list to search for memcheck.h
#
# Sets the follwing variable:
#
# Valgrind_FOUND True if Valgrind was found
# VALGRIND_INCLUDE_DIR Path to Valgrind's include dirs.
# look for header files, only at positions given by the user
find_path(VALGRIND_INCLUDE_DIR
NAMES "valgrind/memcheck.h"
PATHS ${VALGRIND_ROOT}
PATH_SUFFIXES "include"
NO_DEFAULT_PATH
)
# look for header files, including default paths
find_path(VALGRIND_INCLUDE_DIR
NAMES "valgrind/memcheck.h"
PATH_SUFFIXES "include"
)
# handle package arguments
include(FindPackageHandleStandardArgs)
find_package_handle_standard_args(
"Valgrind"
DEFAULT_MSG
VALGRIND_INCLUDE_DIR
)
# set HAVE_VALGRIND for config.h
set(HAVE_VALGRIND ${Valgrind_FOUND})
# register all Valgrind related flags
if(Valgrind_FOUND)
dune_register_package_flags(INCLUDE_DIRS "${VALGRIND_INCLUDE_DIR}")
endif()
......@@ -756,11 +756,11 @@ url = {http://www.sciencedirect.com/science/article/pii/S0169772204001160}
edition = {1}
}
% no doi or link available
@BOOK{kays2005,
title = {{Convective heat and mass transfer}},
publisher = {McGraw-Hill Higher Education},
year = {2005},
isbn={9780071238298},
author = {W. M. Kays and M. E. Crawford and B. Weigand},
edition = {4},
}
......
......@@ -24,6 +24,7 @@
#ifndef DUMUX_BOUNDINGBOXTREE_HH
#define DUMUX_BOUNDINGBOXTREE_HH
#include <algorithm>
#include <memory>
#include <type_traits>
#include <dune/geometry/referenceelements.hh>
......@@ -300,10 +301,7 @@ public:
// we know the points are aligned
// if the dot product is positive and the length in range
// the point is in the interval
if (v1.dot(v2) > 0.0 && v2norm < v1norm*(1 + eps_))
return true;
return false;
return (v1.dot(v2) > 0.0 && v2norm < v1norm*(1 + eps_));
}
/*!
......@@ -546,10 +544,7 @@ public:
// we know the points are aligned
// if the dot product is positive and the length in range
// the point is in the interval
if (v1.dot(v2) > 0.0 && v2norm < v1norm*(1 + eps_))
return true;
return false;
return (v1.dot(v2) > 0.0 && v2norm < v1norm*(1 + eps_));
}
/*!
......@@ -695,10 +690,9 @@ public:
// the point is inside if the length is
// small than the interval length and the
// sign of v1 & v2 are the same
if (std::signbit(v1) == std::signbit(v2)
&& std::abs(v1) < std::abs(v2)*(1 + eps_))
return true;
return false;
using std::abs;
return (std::signbit(v1) == std::signbit(v2)
&& abs(v1) < abs(v2)*(1 + eps_));
}
/*!
......@@ -1175,8 +1169,10 @@ private:
corner = geometry.corner(vLocalIdx);
for (std::size_t dimIdx = 0; dimIdx < dimworld; ++dimIdx)
{
xMin[dimIdx] = std::min(xMin[dimIdx], corner[dimIdx]);
xMax[dimIdx] = std::max(xMax[dimIdx], corner[dimIdx]);
using std::max;
using std::min;
xMin[dimIdx] = min(xMin[dimIdx], corner[dimIdx]);
xMax[dimIdx] = max(xMax[dimIdx], corner[dimIdx]);
}
}
}
......
......@@ -27,6 +27,8 @@
#ifndef DIMENSIONLESS_NUMBERS_HH
#define DIMENSIONLESS_NUMBERS_HH
#include <cmath>
#include <dune/common/exceptions.hh>
namespace Dumux
......@@ -148,13 +150,15 @@ static Scalar nusseltNumberForced(const Scalar reynoldsNumber,
* Dittus, F.W and Boelter, L.M.K, Heat Transfer in Automobile Radiators of the Tubular Type,
* Publications in Engineering, Vol. 2, pages 443-461, 1930
*/
return 0.023 * pow(reynoldsNumber, 0.8) * pow(prandtlNumber,0.33);
using std::pow;
return 0.023 * pow(reynoldsNumber, 0.8) * pow(prandtlNumber,0.33);
}
else if (formulation == NusseltFormulation::WakaoKaguei){
/* example: flow through porous medium *single phase*, fit to many different data
* Wakao and Kaguei, Heat and mass Transfer in Packed Beds, Gordon and Breach Science Publishers, page 293
*/
using std::pow;
return 2. + 1.1 * pow(prandtlNumber,(1./3.)) * pow(reynoldsNumber, 0.6);
}
......@@ -163,6 +167,8 @@ static Scalar nusseltNumberForced(const Scalar reynoldsNumber,
* valid for 0.1<Re<10000, 0.6<Pr/Sc<10000, packed beds of perfect spheres.
*
*/
using std::sqrt;
using std::pow;
Scalar numerator = 0.037 * pow(reynoldsNumber,0.8) * prandtlNumber ;
Scalar reToMin01 = pow(reynoldsNumber,-0.1);
Scalar prTo23 = pow(prandtlNumber, (2./3. ) ) ; // MIND THE pts! :-( otherwise the integer exponent version is chosen
......@@ -258,7 +264,9 @@ static Scalar sherwoodNumber(const Scalar reynoldsNumber,
/* example: flow through porous medium *single phase*
* Wakao and Kaguei, Heat and mass Transfer in Packed Beds, Gordon and Breach Science Publishers, page 156
*/
return 2. + 1.1 * pow(schmidtNumber,(1./3.)) * pow(reynoldsNumber, 0.6);
using std::cbrt;
using std::pow;
return 2. + 1.1 * cbrt(schmidtNumber) * pow(reynoldsNumber, 0.6);
}
else {
......
......@@ -24,16 +24,13 @@
#ifndef DUMUX_EIGENVALUES_HH
#define DUMUX_EIGENVALUES_HH
#include <algorithm>
#include <cmath>
namespace Dumux
{
#include "math.hh"
template<class ValueType>
int sign(const ValueType& value)
namespace Dumux
{
return (value < 0 ? -1 : 1);
}
template<int dim, class Matrix>
void identityMatrix(Matrix& matrix)
......@@ -54,11 +51,13 @@ double calcOffDiagonalNorm(Matrix& matrix)
norm += matrix[i][j] * matrix[i][j];
}
return std::sqrt(norm);
using std::sqrt;
return sqrt(norm);
}
//! Function to calculate eigenvalues of n x n matrices
/*
/*!
* \briefFunction to calculate eigenvalues of n x n matrices
*
* \param eigVel Vector for storing the eigenvalues
* \param matrix n x n matrices for which eigenvalues have to be calculated
* \param relativeTolerance tolerance for the relative convergence criterion (default: 0.01)
......@@ -76,10 +75,12 @@ bool calculateEigenValues(EVVectorType &eigVel, MatrixType& matrix, double relat
eigVel[0] = (-b + sqrt(b * b - 4.0 * c)) / 2.0;
eigVel[1] = (-b - sqrt(b * b - 4.0 * c)) / 2.0;
if (std::isnan(eigVel[0]) || std::isinf(eigVel[0]))
using std::isnan;
using std::isinf;
if (isnan(eigVel[0]) || isinf(eigVel[0]))
return false;
if (std::isnan(eigVel[1]) || std::isinf(eigVel[1]))
if (isnan(eigVel[1]) || isinf(eigVel[1]))
return false;
return true;
......@@ -103,9 +104,10 @@ bool calculateEigenValues(EVVectorType &eigVel, MatrixType& matrix, double relat
double theta = (evMatrix[i][i] - evMatrix[j][j])
/ (2 * evMatrix[i][j]);
double t = sign(theta)
/ (std::abs(theta) + std::sqrt(1 + theta * theta));
double c = 1 / std::sqrt(1 + t * t);
using std::abs;
using std::sqrt;
double t = sign(theta) / (abs(theta) + sqrt(1 + theta * theta));
double c = 1 / sqrt(1 + t * t);
double s = c * t;
rotationMatrix[i][i] = c;
......@@ -129,7 +131,9 @@ bool calculateEigenValues(EVVectorType &eigVel, MatrixType& matrix, double relat
for (int i = 0; i < dim; i++)
{
eigVel[i] = evMatrix[i][i];
if (std::isnan(eigVel[i]) || std::isinf(eigVel[i]))
using std::isinf;
using std::isnan;
if (isnan(eigVel[i]) || isinf(eigVel[i]))
return false;
}
......@@ -168,9 +172,10 @@ bool calculateEigenValues(EVVectorType &eigVel, MatrixType& eigVec, MatrixType&
double theta = (evMatrix[i][i] - evMatrix[j][j])
/ (2 * evMatrix[i][j]);
double t = sign(theta)
/ (std::abs(theta) + std::sqrt(1 + theta * theta));
double c = 1 / std::sqrt(1 + t * t);
using std::abs;
using std::sqrt;
double t = sign(theta) / (abs(theta) + sqrt(1 + theta * theta));
double c = 1 / sqrt(1 + t * t);
double s = c * t;
rotationMatrix[i][i] = c;
......@@ -195,7 +200,9 @@ bool calculateEigenValues(EVVectorType &eigVel, MatrixType& eigVec, MatrixType&
for (int i = 0; i < dim; i++)
{
eigVel[i] = evMatrix[i][i];
if (std::isnan(eigVel[i]) || std::isinf(eigVel[i]))
using std::isinf;
using std::isnan;
if (isnan(eigVel[i]) || isinf(eigVel[i]))
return false;
for (int j = 0; j < dim; j++)
{
......
......@@ -23,13 +23,13 @@
#ifndef DUMUX_MATH_HH
#define DUMUX_MATH_HH
#include <algorithm>
#include <cmath>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/dynmatrix.hh>
#include <cmath>
#include <algorithm>
namespace Dumux
{
/*!
......@@ -59,7 +59,8 @@ Scalar geometricMean(Scalar x, Scalar y)
{
if (x*y <= 0)
return 0;
return std::sqrt(x*y)*((x < 0)?-1:1);
using std::sqrt;
return sqrt(x*y)*sign(x);
}
/*!
......@@ -148,13 +149,17 @@ int invertQuadraticPolynomial(SolContainer &sol,
if (Delta < 0)
return 0; // no real roots
Delta = std::sqrt(Delta);
using std::sqrt;
Delta = 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]);
{
using std::swap;
swap(sol[0], sol[1]);
}
return 2; // two real roots
}
......@@ -179,7 +184,8 @@ void invertCubicPolynomialPostProcess_(SolContainer &sol,
x -= fOld/fPrime;
Scalar fNew = d + x*(c + x*(b + x*a));
if (std::abs(fNew) < std::abs(fOld))
using std::abs;
if (abs(fNew) < abs(fOld))
sol[i] = x;
}
}
......@@ -234,9 +240,8 @@ int invertCubicPolynomial(SolContainer *sol,
// 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);
using std::cbrt;
Scalar t = cbrt(q);
sol[0] = t - b/3;
return 1;
......@@ -251,9 +256,10 @@ int invertCubicPolynomial(SolContainer *sol,
}
// two additional real roots at t = sqrt(-p) and t = -sqrt(-p)
sol[0] = -std::sqrt(-p) - b/3;
using std::sqrt;
sol[0] = -sqrt(-p) - b/3;
sol[1] = 0.0 - b/3;
sol[2] = std::sqrt(-p) - b/3;
sol[2] = sqrt(-p) - b/3;
return 3;
}
......@@ -291,9 +297,9 @@ int invertCubicPolynomial(SolContainer *sol,
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);
using std::cbrt;
using std::sqrt;
Scalar u = cbrt(-q/2 + sqrt(wDisc));
// at this point, u != 0 since p^3 = 0 is necessary in order
// for u = 0 to hold, so
......@@ -306,10 +312,13 @@ int invertCubicPolynomial(SolContainer *sol,
}
else { // the negative discriminant case:
Scalar uCubedRe = - q/2;
Scalar uCubedIm = std::sqrt(-wDisc);
using std::sqrt;
Scalar uCubedIm = 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;
using std::cbrt;
Scalar uAbs = cbrt(sqrt(uCubedRe*uCubedRe + uCubedIm*uCubedIm));
using std::atan2;
Scalar phi = atan2(uCubedIm, uCubedRe)/3;
// calculate the length and the angle of the primitive root
......@@ -351,7 +360,8 @@ int invertCubicPolynomial(SolContainer *sol,
// 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;
using std::cos;
sol[i] = cos(phi)*(uAbs - p/(3*uAbs)) - b/3;
phi += 2*M_PI/3;
}
......@@ -360,7 +370,8 @@ int invertCubicPolynomial(SolContainer *sol,
invertCubicPolynomialPostProcess_(sol, 3, a, b, c, d);
// sort the result
std::sort(sol, sol + 3);
using std::sort;
sort(sol, sol + 3);
return 3;
}
......@@ -470,7 +481,21 @@ Scalar antoine(Scalar temperature,
Scalar C)
{
const Scalar ln10 = 2.3025850929940459;
return std::exp(ln10*(A - B/(C + temperature)));
using std::exp;
return exp(ln10*(A - B/(C + temperature)));
}
/*!
* \brief Sign or signum function.
*
* Returns 1 for a positive argument.
* Returns -1 for a negative argument.
* Returns 0 if the argument is zero.
*/
template<class ValueType>
int sign(const ValueType& value)
{
return (ValueType(0) < value) - (value < ValueType(0));
}
/*!
......
......@@ -23,14 +23,15 @@
#ifndef DUMUX_SPLINE_COMMON__HH
#define DUMUX_SPLINE_COMMON__HH
#include <algorithm>
#include <iostream>
#include <cassert>
#include <dune/common/exceptions.hh>
#include <dune/common/float_cmp.hh>
#include "valgrind.hh"
#include "math.hh"
#include "valgrind.hh"
namespace Dumux
{
......@@ -55,19 +56,19 @@ public:
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); };
{ return x_(0); }
/*!
* \brief Return the x value of the rightmost sampling point.
*/
Scalar xMax() const
{ return x_(numSamples_() - 1); };
{ return x_(numSamples_() - 1); }
/*!
* \brief Prints k tuples of the format (x, y, dx/dy, isMonotonic)
......@@ -88,8 +89,10 @@ public:
*/
void printCSV(Scalar xi0, Scalar xi1, int k) const
{
Scalar x0 = std::min(xi0, xi1);
Scalar x1 = std::max(xi0, xi1);
using std::max;
using std::min;
Scalar x0 = min(xi0, xi1);
Scalar x1 = max(xi0, xi1);
const int n = numSamples_() - 1;
for (int i = 0; i <= k; ++i) {
double x = i*(x1 - x0)/k + x0;
......@@ -116,7 +119,7 @@ public:
else {
y = eval(x);
dy_dx = evalDerivative(x);
mono = monotonic(std::max<Scalar>(x_(0), x), std::min<Scalar>(x_(n), x_p1));
mono = monotonic(max<Scalar>(x_(0), x), min<Scalar>(x_(n), x_p1));
}
std::cout << x << " " << y << " " << dy_dx << " " << mono << "\n";
......@@ -190,7 +193,7 @@ public:
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
......@@ -222,7 +225,7 @@ public:
"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
......@@ -240,9 +243,10 @@ public:
// make sure that x0 is smaller than x1
if (x0 > x1)
std::swap(x0, x1);
assert(x0 < x1);
{
using std::swap;
swap(x0, x1);
}
// corner case where the whole spline is a constant
if (moment_(0) == 0 &&
......@@ -287,7 +291,7 @@ public:
protected:
// this is an internal class, so everything is protected!
SplineCommon_()
{ Valgrind::SetUndefined(asImp_()); };
{ Valgrind::SetUndefined(asImp_()); }
/*!
* \brief Set the sampling point vectors.
......@@ -551,7 +555,8 @@ protected:
// not exhibit any extrema.
return (x0*(x0*3*a + 2*b) + c > 0) ? 1 : -1;
}
disc = std::sqrt(disc);
using std::sqrt;
disc = sqrt(disc);
Scalar xE1 = (-2*b + disc)/(6*a);
Scalar xE2 = (-2*b - disc)/(6*a);
......@@ -562,7 +567,7 @@ protected:
// to determine whether we're monotonically increasing
// or decreasing
x0 = x1;
return (x0*(x0*3*a + 2*b) + c > 0) ? 1 : -1;
return sign(x0*(x0*3*a + 2*b) + c);
}
if ((x0 < xE1 && xE1 < x1) ||
(x0 < xE2 && xE2 < x1))
......@@ -573,7 +578,7 @@ protected:
// 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;
return sign(x0*(x0*3*a + 2*b) + c);
}
/*!
......@@ -590,8 +595,9 @@ protected:
b_(segIdx) - b,
c_(segIdx) - c,
d_(segIdx) - d);
x0 = std::max(x_(segIdx), x0);
x1 = std::max(x_(segIdx+1), x1);
using std::max;
x0 = max(x_(segIdx), x0);
x1 = max(x_(segIdx+1), x1);
// filter the intersections outside of the specified intervall
int k = 0;
......
......@@ -79,7 +79,7 @@ public:
yMin_ = yMin;
yMax_ = yMax;
};
}
/*!
* \brief Return the position on the x-axis of the i-th interval.
......@@ -89,7 +89,7 @@ public:
assert(0 <= i && i < m_);
return xMin_ + i*(xMax_ - xMin_)/(m_ - 1);
};
}
/*!
* \brief Return the position on the y-axis of the j-th interval.
......@@ -99,7 +99,7 @@ public:
assert(0 <= j && j < n_);
return yMin_ + j*(yMax_ - yMin_)/(n_ - 1);
};
}
/*!
* \brief Return the interval index of a given position on the x-axis.
......@@ -112,7 +112,7 @@ public:
Scalar xToI(Scalar x) const
{
return (x - xMin_)/(xMax_ - xMin_)*m_;
};
}
/*!
......@@ -126,7 +126,7 @@ public:
Scalar yToJ(Scalar y) const
{
return (y - yMin_)/(yMax_ - yMin_)*n_;
};
}
/*!
......@@ -140,7 +140,7 @@ public:
assert(0 <= j && j < n_);
return samples_[j*m_ + i];
};
}
/*!
* \brief Set the value of the sample point which is at the
......@@ -153,7 +153,7 @@ public:
assert(0 <= j && j < n_);
samples_[j*m_ + i] = value;
};
}