diff --git a/dumux/common/math.hh b/dumux/common/math.hh index 905d5aab727e7e9ad821e103f188bc6b9f7c3a8c..0c2dd6e89d8257472a0d8bb8f8d9ad54aac517cf 100644 --- a/dumux/common/math.hh +++ b/dumux/common/math.hh @@ -42,8 +42,10 @@ namespace Dumux * \param y The second input value */ template <class Scalar> -Scalar harmonicMean(Scalar x, Scalar y) +constexpr Scalar harmonicMean(Scalar x, Scalar y) noexcept { + static_assert(Dune::IsNumber<Scalar>::value, "The arguments x and y have to be numbers!"); + if (x*y <= 0) return 0; return (2*x*y)/(x + y); @@ -55,10 +57,13 @@ Scalar harmonicMean(Scalar x, Scalar y) * * \param x The first input value * \param y The second input value + * \note as std::sqrt is not constexpr this function is not constexpr */ template <class Scalar> -Scalar geometricMean(Scalar x, Scalar y) +Scalar geometricMean(Scalar x, Scalar y) noexcept { + static_assert(Dune::IsNumber<Scalar>::value, "The arguments x and y have to be numbers!"); + if (x*y <= 0) return 0; using std::sqrt; @@ -496,7 +501,7 @@ Scalar antoine(Scalar temperature, * Returns 0 if the argument is zero. */ template<class ValueType> -int sign(const ValueType& value) +constexpr int sign(const ValueType& value) noexcept { return (ValueType(0) < value) - (value < ValueType(0)); } @@ -632,20 +637,21 @@ Dune::FieldMatrix<Scalar, rows1, cols2> multiplyMatrices(const Dune::FieldMatrix /*! * \ingroup Common - * \brief Trace of dynamic matrix + * \brief Trace of a dense matrix * - * \param M The dynamic matrix + * \param M The dense matrix */ -template <class Scalar> -Scalar trace(const Dune::DynamicMatrix<Scalar>& M) +template <class MatrixType> +typename Dune::DenseMatrix<MatrixType>::field_type +trace(const Dune::DenseMatrix<MatrixType>& M) { - std::size_t rows_T = M.M(); - - DUNE_ASSERT_BOUNDS(rows_T == M.N()); + const auto rows = M.N(); + DUNE_ASSERT_BOUNDS(rows == M.M()); // rows == cols - Scalar trace = 0.0; + using MatType = Dune::DenseMatrix<MatrixType>; + typename MatType::field_type trace = 0.0; - for (std::size_t i = 0; i < rows_T; ++i) + for (typename MatType::size_type i = 0; i < rows; ++i) trace += M[i][i]; return trace; diff --git a/test/common/math/CMakeLists.txt b/test/common/math/CMakeLists.txt index 08a8c62e5572013f0a57c3d09b17687828a37277..f5836ca31593b4e2dbd52e693fa2e3ca4e3237f7 100644 --- a/test/common/math/CMakeLists.txt +++ b/test/common/math/CMakeLists.txt @@ -1,7 +1,7 @@ # build the test for the property system -dune_add_test(SOURCES test_vtmv.cc) +dune_add_test(SOURCES test_math.cc) #install sources install(FILES -test_vtmv.cc -DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/test/common/math) \ No newline at end of file +test_math.cc +DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/test/common/math) diff --git a/test/common/math/test_vtmv.cc b/test/common/math/test_math.cc similarity index 71% rename from test/common/math/test_vtmv.cc rename to test/common/math/test_math.cc index d59d6be3071a097bdd35888e54975843d9a9a151..7c62d3e0c3a94dcd2545a5ad96065a5ae5db9e15 100644 --- a/test/common/math/test_vtmv.cc +++ b/test/common/math/test_math.cc @@ -19,8 +19,11 @@ /*! * \file * - * \brief This file tests the vtmv function with - * vtmv(Vector, FieldScalar, Vector) and vtmv(Vector, Matrix, Vector). + * \brief This file tests several math functions: + * the vtmv function with vtmv(Vector, FieldScalar, Vector) and vtmv(Vector, Matrix, Vector). + * the trace function + * the harmonicMean function + * \todo test more math functions! * * We declare some vectors and matrices and test the combinations of them * against a previously calculated result. @@ -55,6 +58,11 @@ int main() try K[2][2] = k; Dune::DynamicMatrix<double> K_dyn(K); + + ////////////////////////////////////////////////////////////////// + ///// Dumux::vtmv //////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////// + //! Set reference result. Should be -15 for all combinations const double reference = -15; @@ -90,6 +98,37 @@ int main() try //! Test vtmv function with DynamicVector v1_dyn and DynamicMatrix K_dyn and DynamicVector v2_dyn if (!Dune::FloatCmp::eq(reference, Dumux::vtmv(v1_dyn, K_dyn, v2_dyn), 1e-6)) DUNE_THROW(Dune::Exception, "vtmv-result does not match reference"); + + + ////////////////////////////////////////////////////////////////// + ///// Dumux::trace /////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////// + + const auto t1 = Dumux::trace(K); + const auto t2 = Dumux::trace(K_dyn); + if (!Dune::FloatCmp::eq(t1, t2, 1e-30)) + DUNE_THROW(Dune::Exception, "Traces do not match!"); + + + ////////////////////////////////////////////////////////////////// + ///// Dumux::harmonicMean //////////////////////////////////////// + ////////////////////////////////////////////////////////////////// + + constexpr auto mean = Dumux::harmonicMean(4.0, 5.0); + static_assert( (mean - 4.0*5.0*2.0/(4.0+5.0)) < 1e-30 && (mean - 4.0*5.0*2.0/(4.0+5.0)) > -1e-30 , "Wrong harmonic mean!"); + + ////////////////////////////////////////////////////////////////// + ///// Dumux::sign //////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////// + + static_assert(Dumux::sign(0.0) == 0, "Wrong sign!"); + static_assert(Dumux::sign(-0.0) == 0, "Wrong sign!"); + static_assert(Dumux::sign(0) == 0, "Wrong sign!"); + static_assert(Dumux::sign(-0) == 0, "Wrong sign!"); + static_assert(Dumux::sign(1) == 1, "Wrong sign!"); + static_assert(Dumux::sign(2.0) == 1, "Wrong sign!"); + static_assert(Dumux::sign(-2) == -1, "Wrong sign!"); + static_assert(Dumux::sign(-3.5) == -1, "Wrong sign!"); } catch (Dune::Exception& e) { std::cerr << e << std::endl;