Merge branch 'feature/improve-math-functions' into 'master'

Feature/improve math functions

See merge request !807

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;

test/common/math/CMakeLists.txt

 # 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)

test/common/math/test_vtmv.cc → 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;