From 7485f0a2792799e759c0e73bc52e9a78a74ec312 Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Fri, 17 Apr 2020 20:08:32 +0200
Subject: [PATCH] [fluidmatrixinteraction] Use std::clamp for improved
 readability

---
 .../fluidmatrixinteractions/2p/brookscorey.hh | 30 +++++++----------
 .../2p/linearmaterial.hh                      | 11 +++----
 .../2p/vangenuchten.hh                        | 32 ++++++++-----------
 .../3p/parkervangen3p.hh                      |  5 ++-
 .../3p/regularizedparkervangen3p.hh           | 19 +++++------
 .../mp/mplinearmaterial.hh                    |  5 ++-
 .../porositydeformation.hh                    |  5 ++-
 7 files changed, 43 insertions(+), 64 deletions(-)

diff --git a/dumux/material/fluidmatrixinteractions/2p/brookscorey.hh b/dumux/material/fluidmatrixinteractions/2p/brookscorey.hh
index fc04aac05c..cc9827550a 100644
--- a/dumux/material/fluidmatrixinteractions/2p/brookscorey.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/brookscorey.hh
@@ -72,10 +72,9 @@ public:
     static Scalar pc(const Params &params, Scalar swe)
     {
         using std::pow;
-        using std::min;
-        using std::max;
+        using std::clamp;
 
-        swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
+        swe = clamp(swe, 0.0, 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
 
         return params.pe()*pow(swe, -1.0/params.lambda());
     }
@@ -136,10 +135,9 @@ public:
     static Scalar dpc_dswe(const Params &params, Scalar swe)
     {
         using std::pow;
-        using std::min;
-        using std::max;
+        using std::clamp;
 
-        swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
+        swe = clamp(swe, 0.0, 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
 
         return - params.pe()/params.lambda() * pow(swe, -1/params.lambda() - 1);
     }
@@ -184,10 +182,9 @@ public:
     static Scalar krw(const Params &params, Scalar swe)
     {
         using std::pow;
-        using std::min;
-        using std::max;
+        using std::clamp;
 
-        swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
+        swe = clamp(swe, 0.0, 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
 
         return pow(swe, 2.0/params.lambda() + 3);
     }
@@ -210,10 +207,9 @@ public:
     static Scalar dkrw_dswe(const Params &params, Scalar swe)
     {
         using std::pow;
-        using std::min;
-        using std::max;
+        using std::clamp;
 
-        swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
+        swe = clamp(swe, 0.0, 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
 
         return (2.0/params.lambda() + 3)*pow(swe, 2.0/params.lambda() + 2);
     }
@@ -235,10 +231,9 @@ public:
     static Scalar krn(const Params &params, Scalar swe)
     {
         using std::pow;
-        using std::min;
-        using std::max;
+        using std::clamp;
 
-        swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
+        swe = clamp(swe, 0.0, 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
 
         const Scalar exponent = 2.0/params.lambda() + 1;
         const Scalar tmp = 1.0 - swe;
@@ -264,10 +259,9 @@ public:
     static Scalar dkrn_dswe(const Params &params, Scalar swe)
     {
         using std::pow;
-        using std::min;
-        using std::max;
+        using std::clamp;
 
-        swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
+        swe = clamp(swe, 0.0, 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
 
         const auto lambdaInv = 1.0/params.lambda();
         const auto swePow = pow(swe, 2*lambdaInv);
diff --git a/dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh b/dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh
index bb3fd0a91b..e7299fbd19 100644
--- a/dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh
@@ -144,9 +144,8 @@ public:
      */
     static Scalar krw(const Params &params, Scalar swe)
     {
-        using std::max;
-        using std::min;
-        return max(min(swe,1.0),0.0);
+        using std::clamp;
+        return clamp(swe, 0.0, 1.0);
     }
 
     /*!
@@ -160,10 +159,8 @@ public:
      */
     static Scalar krn(const Params &params, Scalar swe)
     {
-        Scalar sne = 1 - swe;
-        using std::max;
-        using std::min;
-        return max(min(sne,1.0),0.0);
+        using std::clamp;
+        return clamp(1.0-swe, 0.0, 1.0);
     }
 };
 } // end namespace Dumux
diff --git a/dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh b/dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh
index 5f12527fbe..36cddb1501 100644
--- a/dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh
@@ -72,10 +72,9 @@ public:
     static Scalar pc(const Params &params, Scalar swe)
     {
         using std::pow;
-        using std::min;
-        using std::max;
+        using std::clamp;
 
-        swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
+        swe = clamp(swe, 0.0, 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
 
         const Scalar pc = pow(pow(swe, -1.0/params.vgm()) - 1, 1.0/params.vgn())/params.vgAlpha();
         return pc;
@@ -142,10 +141,9 @@ public:
     static Scalar dpc_dswe(const Params &params, Scalar swe)
     {
         using std::pow;
-        using std::min;
-        using std::max;
+        using std::clamp;
 
-        swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
+        swe = clamp(swe, 0.0, 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
 
         const Scalar powSwe = pow(swe, -1/params.vgm());
         return - 1.0/params.vgAlpha() * pow(powSwe - 1, 1.0/params.vgn() - 1)/params.vgn()
@@ -193,10 +191,9 @@ public:
     static Scalar krw(const Params &params, Scalar swe)
     {
         using std::pow;
-        using std::min;
-        using std::max;
+        using std::clamp;
 
-        swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
+        swe = clamp(swe, 0.0, 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
 
         const Scalar r = 1.0 - pow(1.0 - pow(swe, 1.0/params.vgm()), params.vgm());
         return pow(swe, params.vgl())*r*r;
@@ -219,10 +216,9 @@ public:
     static Scalar dkrw_dswe(const Params &params, Scalar swe)
     {
         using std::pow;
-        using std::min;
-        using std::max;
+        using std::clamp;
 
-        swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
+        swe = clamp(swe, 0.0, 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
 
         const Scalar x = 1.0 - pow(swe, 1.0/params.vgm());
         const Scalar xToM = pow(x, params.vgm());
@@ -242,15 +238,14 @@ public:
      *
      * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
      *       by clamping the input.
-     * \note See e.g. Dury, Fischer, Schulin (1999) for application of Mualem model to non-wetting rel. perm. 
+     * \note See e.g. Dury, Fischer, Schulin (1999) for application of Mualem model to non-wetting rel. perm.
      */
     static Scalar krn(const Params &params, Scalar swe)
     {
         using std::pow;
-        using std::min;
-        using std::max;
+        using std::clamp;
 
-        swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
+        swe = clamp(swe, 0.0, 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
 
         return pow(1 - swe, params.vgl()) * pow(1 - pow(swe, 1.0/params.vgm()), 2*params.vgm());
     }
@@ -273,10 +268,9 @@ public:
     static Scalar dkrn_dswe(const Params &params, Scalar swe)
     {
         using std::pow;
-        using std::min;
-        using std::max;
+        using std::clamp;
 
-        swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
+        swe = clamp(swe, 0.0, 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
 
         const auto sne = 1.0 - swe;
         const auto x = 1.0 - pow(swe, 1.0/params.vgm());
diff --git a/dumux/material/fluidmatrixinteractions/3p/parkervangen3p.hh b/dumux/material/fluidmatrixinteractions/3p/parkervangen3p.hh
index a516dea681..436269fb09 100644
--- a/dumux/material/fluidmatrixinteractions/3p/parkervangen3p.hh
+++ b/dumux/material/fluidmatrixinteractions/3p/parkervangen3p.hh
@@ -235,13 +235,12 @@ public:
         krn -= pow(1 - pow(ste, 1/params.vgm()), params.vgm());
         krn *= krn;
 
-        using std::max;
-        using std::min;
+        using std::clamp;
         using std::sqrt;
         if (params.krRegardsSnr())
         {
             // regard Snr in the permeability of the n-phase, see Helmig1997
-            Scalar resIncluded = max(min((sn - params.snr()/ (1-params.swr())), 1.0), 0.0);
+            const Scalar resIncluded = clamp(sn - params.snr()/ (1-params.swr()), 0.0, 1.0);
             krn *= sqrt(resIncluded );
         }
         else
diff --git a/dumux/material/fluidmatrixinteractions/3p/regularizedparkervangen3p.hh b/dumux/material/fluidmatrixinteractions/3p/regularizedparkervangen3p.hh
index c82d4fcb1d..dfdef7a824 100644
--- a/dumux/material/fluidmatrixinteractions/3p/regularizedparkervangen3p.hh
+++ b/dumux/material/fluidmatrixinteractions/3p/regularizedparkervangen3p.hh
@@ -25,6 +25,8 @@
 #ifndef REGULARIZED_PARKERVANGEN_3P_HH
 #define REGULARIZED_PARKERVANGEN_3P_HH
 
+#include <algorithm>
+
 #include "parkervangen3p.hh"
 #include "regularizedparkervangen3pparams.hh"
 
@@ -93,13 +95,9 @@ public:
     static Scalar pcgw(const Params &params, Scalar swe)
     {
         //if specified, a constant value is used for regularization
-        if(params.constRegularization())
-        {
-            if(swe < 0.0)
-                swe = 0.0;
-            if(swe > 1.0)
-                swe = 1.0;
-        }
+        using std::clamp;
+        if (params.constRegularization())
+            swe = clamp(swe, 0.0, 1.0);
 
         // retrieve the low and the high threshold saturations for the
         // unregularized capillary pressure curve from the parameters
@@ -333,10 +331,9 @@ public:
      */
     static Scalar krn(const Params &params, Scalar swe, Scalar sn, Scalar ste)
     {
-        using std::max;
-        using std::min;
-        swe = max(min(swe, 1.0), 0.0);
-        ste = max(min(ste, 1.0), 0.0);
+        using std::clamp;
+        swe = clamp(swe, 0.0, 1.0);
+        ste = clamp(ste, 0.0, 1.0);
 
         if(ste - swe <= 0.0)
             return 0.0;
diff --git a/dumux/material/fluidmatrixinteractions/mp/mplinearmaterial.hh b/dumux/material/fluidmatrixinteractions/mp/mplinearmaterial.hh
index 2925710ba2..e5bb06f603 100644
--- a/dumux/material/fluidmatrixinteractions/mp/mplinearmaterial.hh
+++ b/dumux/material/fluidmatrixinteractions/mp/mplinearmaterial.hh
@@ -90,10 +90,9 @@ public:
                                        const FluidState &state,
                                        int wPhaseIdx = 0)
     {
-        using std::max;
-        using std::min;
+        using std::clamp;
         for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
-            values[phaseIdx] = max(min(state.saturation(phaseIdx),1.0),0.0);
+            values[phaseIdx] = clamp(state.saturation(phaseIdx), 0.0, 1.0);
     }
 };
 } // end namespace Dumux
diff --git a/dumux/material/fluidmatrixinteractions/porositydeformation.hh b/dumux/material/fluidmatrixinteractions/porositydeformation.hh
index 2389db780a..e9156a07e7 100644
--- a/dumux/material/fluidmatrixinteractions/porositydeformation.hh
+++ b/dumux/material/fluidmatrixinteractions/porositydeformation.hh
@@ -76,9 +76,8 @@ public:
         for (int dir = 0; dir < FVGridGeom::GridView::dimension; ++dir)
             divU += gradU[dir][dir];
 
-        using std::min;
-        using std::max;
-        return min( maxPoro, max(minPoro, (refPoro+divU)/(1.0+divU)) );
+        using std::clamp;
+        return clamp((refPoro+divU)/(1.0+divU), minPoro, maxPoro);
     }
 
     /*!
-- 
GitLab