[math] Add smoothMin/Max to common/math.hh

dumux/common/math.hh

@@ -111,6 +111,35 @@ void harmonicMeanMatrix(Dune::FieldMatrix<Scalar, m, n> &K,
     }
 }
 
+/*!
+ * \ingroup Core
+ * \brief A smoothed minimum function (using cubic interpolation)
+ * \note The returned value is guaranteed to be smaller or equal to the minimum of the two input values.
+ * \note The smoothing kicks in if the difference between the two input values is smaller than the smoothing parameter.
+ * \param a The first input value
+ * \param b The second input value
+ * \param k The smoothing parameter (k > 0)
+ */
+template<class Scalar>
+Scalar smoothMin(const Scalar a, const Scalar b, const Scalar k)
+{
+    using std::max; using std::min; using std::abs;
+    const auto h = max(k-abs(a-b), 0.0 )/k;
+    return min(a, b) - h*h*h*k*(1.0/6.0);
+}
+
+/*!
+ * \ingroup Core
+ * \brief A smoothed maximum function (using cubic interpolation)
+ *
+ * \param a The first input value
+ * \param b The second input value
+ * \param k The smoothing parameter (k > 0)
+ */
+template<class Scalar>
+Scalar smoothMax(const Scalar a, const Scalar b, const Scalar k)
+{ return -smoothMin(-a, -b, k); }
+
 /*!
  * \ingroup Core
  * \brief Invert a linear polynomial analytically

test/common/math/test_math.cc

@@ -377,5 +377,25 @@ int main()
     if (!Dune::FloatCmp::eq(roots[0], -3.863448718244389e-02, 1e-14))
         DUNE_THROW(Dune::Exception, "Root of cubic equation does not match reference");
 
+    // Test smoothMin and smoothMax
+    {
+        const double a = 1.0;
+        const double b = 2.0;
+        const double k = 0.5;
+        if (!Dune::FloatCmp::eq(Dumux::smoothMin(a, b, k), a))
+            DUNE_THROW(Dune::Exception, "smoothMin does not work as expected");
+        if (!Dune::FloatCmp::eq(Dumux::smoothMax(a, b, k), b))
+            DUNE_THROW(Dune::Exception, "smoothMax does not work as expected");
+    }
+    {
+        const double a = 1.0;
+        const double b = 2.0;
+        const double k = 2.0;
+        if (Dumux::smoothMin(a, b, k) >= a)
+            DUNE_THROW(Dune::Exception, "smoothMin does not work as expected");
+        if (Dumux::smoothMax(a, b, k) <= b)
+            DUNE_THROW(Dune::Exception, "smoothMax does not work as expected");
+    }
+
     return 0;
 }