diff --git a/dumux/common/math.hh b/dumux/common/math.hh index 314bc2d2f2b15509ff9ee9c94321396c0bbf4c5c..623dfb1dabbeeabf75987ab380d0efd81337e112 100644 --- a/dumux/common/math.hh +++ b/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 diff --git a/test/common/math/test_math.cc b/test/common/math/test_math.cc index 0f06b0a035352f110e93c09fdc7d4c3334e45dd7..e7404d55215d5334d5d8d0d380307eebab1de504 100644 --- a/test/common/math/test_math.cc +++ b/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; }