STYLE: use clamp/clamp_range instead of min(max(..., upper), lower)

applications/solvers/combustion/XiFoam/bEqn.H

@@ -215,7 +215,7 @@ if (ign.ignited())
           + (
                 scalar(1)
               + (2*XiShapeCoef)
-               *(scalar(0.5) - min(max(b, scalar(0)), scalar(1)))
+               *(scalar(0.5) - clamp(b, zero_one{}))
             )*(XiEqStar - scalar(1.001))
         );
 

applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/constant/constant.C

@@ -99,17 +99,6 @@ Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::mDotAlphal() const
 Foam::Pair<Foam::tmp<Foam::volScalarField>>
 Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::mDot() const
 {
-
-    volScalarField limitedAlpha1
-    (
-        min(max(mixture_.alpha1(), scalar(0)), scalar(1))
-    );
-
-    volScalarField limitedAlpha2
-    (
-        min(max(mixture_.alpha2(), scalar(0)), scalar(1))
-    );
-
     const volScalarField& T = mesh_.lookupObject<volScalarField>("T");
 
     const twoPhaseMixtureEThermo& thermo =
@@ -124,11 +113,15 @@ Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::mDot() const
 
     volScalarField mDotE
     (
-        "mDotE", coeffE_*mixture_.rho1()*limitedAlpha1*max(T - TSat, T0)
+        "mDotE",
+        coeffE_*mixture_.rho1()*clamp(mixture_.alpha1(), zero_one{})
+      * max(T - TSat, T0)
     );
     volScalarField mDotC
     (
-        "mDotC", coeffC_*mixture_.rho2()*limitedAlpha2*max(TSat - T, T0)
+        "mDotC",
+        coeffC_*mixture_.rho2()*clamp(mixture_.alpha2(), zero_one{})
+      * max(TSat - T, T0)
     );
 
     if (mesh_.time().outputTime())
@@ -148,16 +141,6 @@ Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::mDot() const
 Foam::Pair<Foam::tmp<Foam::volScalarField>>
 Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::mDotDeltaT() const
 {
-    volScalarField limitedAlpha1
-    (
-        min(max(mixture_.alpha1(), scalar(0)), scalar(1))
-    );
-
-    volScalarField limitedAlpha2
-    (
-        min(max(mixture_.alpha2(), scalar(0)), scalar(1))
-    );
-
     const volScalarField& T = mesh_.lookupObject<volScalarField>("T");
 
     const twoPhaseMixtureEThermo& thermo =
@@ -170,8 +153,14 @@ Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::mDotDeltaT() const
 
     return Pair<tmp<volScalarField>>
     (
-        coeffC_*mixture_.rho2()*limitedAlpha2*pos(TSat - T),
+        (
-        coeffE_*mixture_.rho1()*limitedAlpha1*pos(T - TSat)
+            coeffC_*mixture_.rho2()*clamp(mixture_.alpha2(), zero_one{})
+          * pos(TSat - T)
+        ),
+        (
+            coeffE_*mixture_.rho1()*clamp(mixture_.alpha1(), zero_one{})
+          * pos(T - TSat)
+        )
     );
 }
 
@@ -201,25 +190,17 @@ Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::TSource() const
 
     const dimensionedScalar& TSat = thermo.TSat();
 
-    dimensionedScalar L = mixture_.Hf2() - mixture_.Hf1();
+    const dimensionedScalar L = mixture_.Hf2() - mixture_.Hf1();
-
-    volScalarField limitedAlpha1
-    (
-        min(max(mixture_.alpha1(), scalar(0)), scalar(1))
-    );
-
-    volScalarField limitedAlpha2
-    (
-        min(max(mixture_.alpha2(), scalar(0)), scalar(1))
-    );
 
     const volScalarField Vcoeff
     (
-        coeffE_*mixture_.rho1()*limitedAlpha1*L*pos(T - TSat)
+        coeffE_*mixture_.rho1()*clamp(mixture_.alpha1(), zero_one{})
+      * L*pos(T - TSat)
     );
     const volScalarField Ccoeff
     (
-        coeffC_*mixture_.rho2()*limitedAlpha2*L*pos(TSat - T)
+        coeffC_*mixture_.rho2()*clamp(mixture_.alpha2(), zero_one{})
+      * L*pos(TSat - T)
     );
 
     TSource =

applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/interfaceHeatResistance/interfaceHeatResistance.C

@@ -167,20 +167,10 @@ Foam::Pair<Foam::tmp<Foam::volScalarField>>
 Foam::temperaturePhaseChangeTwoPhaseMixtures::interfaceHeatResistance::
 mDotAlphal() const
 {
-    volScalarField limitedAlpha1
-    (
-        min(max(mixture_.alpha1(), scalar(0)), scalar(1))
-    );
-
-    volScalarField limitedAlpha2
-    (
-        min(max(mixture_.alpha2(), scalar(0)), scalar(1))
-    );
-
     return Pair<tmp<volScalarField>>
     (
-        (mDotc_/(limitedAlpha2 + SMALL)),
+        (mDotc_/clamp(mixture_.alpha2(), scalarMinMax(SMALL, 1))),
-       -(mDote_/(limitedAlpha1 + SMALL))
+       -(mDote_/clamp(mixture_.alpha1(), scalarMinMax(SMALL, 1)))
     );
 }
 

applications/solvers/multiphase/interCondensatingEvaporatingFoam/temperaturePhaseChangeTwoPhaseMixtures/twoPhaseMixtureEThermo/twoPhaseMixtureEThermo.C

@@ -155,7 +155,7 @@ Foam::tmp<Foam::volScalarField> Foam::twoPhaseMixtureEThermo::Cp() const
 {
     const volScalarField limitedAlpha1
     (
-        min(max(alpha1_, scalar(0)), scalar(1))
+        clamp(alpha1_, zero_one{})
     );
 
     return tmp<volScalarField>
@@ -176,13 +176,11 @@ Foam::tmp<Foam::scalarField> Foam::twoPhaseMixtureEThermo::Cp
     const label patchi
 ) const
 {
-    const volScalarField limitedAlpha1
+    const scalarField alpha1p
     (
-        min(max(alpha1_, scalar(0)), scalar(1))
+        clamp(alpha1_.boundaryField()[patchi], zero_one{})
     );
 
-    const scalarField& alpha1p = limitedAlpha1.boundaryField()[patchi];
-
     return
     (
         alpha1p*Cp1().value() + (scalar(1) - alpha1p)*Cp2().value()
@@ -194,7 +192,7 @@ Foam::tmp<Foam::volScalarField> Foam::twoPhaseMixtureEThermo::rho() const
 {
     const volScalarField limitedAlpha1
     (
-        min(max(alpha1_, scalar(0)), scalar(1))
+        clamp(alpha1_, zero_one{})
     );
 
     return tmp<volScalarField>
@@ -214,13 +212,11 @@ Foam::tmp<Foam::scalarField> Foam::twoPhaseMixtureEThermo::rho
     const label patchi
 ) const
 {
-    const volScalarField limitedAlpha1
+    const scalarField alpha1p
     (
-        min(max(alpha1_, scalar(0)), scalar(1))
+        clamp(alpha1_.boundaryField()[patchi], zero_one{})
     );
 
-    const scalarField& alpha1p = limitedAlpha1.boundaryField()[patchi];
-
     return
     (
         alpha1p*rho1().value() + (scalar(1) - alpha1p)*rho2().value()
@@ -232,7 +228,7 @@ Foam::tmp<Foam::volScalarField> Foam::twoPhaseMixtureEThermo::Cv() const
 {
     const volScalarField limitedAlpha1
     (
-        min(max(alpha1_, scalar(0)), scalar(1))
+        clamp(alpha1_, zero_one{})
     );
 
     return tmp<volScalarField>
@@ -253,13 +249,11 @@ Foam::tmp<Foam::scalarField> Foam::twoPhaseMixtureEThermo::Cv
     const label patchi
 ) const
 {
-    const volScalarField limitedAlpha1
+    const scalarField alpha1p
     (
-        min(max(alpha1_, scalar(0)), scalar(1))
+        clamp(alpha1_.boundaryField()[patchi], zero_one{})
     );
 
-    const scalarField& alpha1p = limitedAlpha1.boundaryField()[patchi];
-
     return
     (
         alpha1p*Cv1().value() + (scalar(1) - alpha1p)*Cv2().value()
@@ -339,7 +333,7 @@ Foam::tmp<Foam::volScalarField> Foam::twoPhaseMixtureEThermo::kappa() const
 {
     const volScalarField limitedAlpha1
     (
-        min(max(alpha1_, scalar(0)), scalar(1))
+        clamp(alpha1_, zero_one{})
     );
 
     return tmp<volScalarField>
@@ -358,13 +352,11 @@ Foam::tmp<Foam::scalarField> Foam::twoPhaseMixtureEThermo::kappa
     const label patchi
 ) const
 {
-    const volScalarField limitedAlpha1
+    const scalarField alpha1p
     (
-        min(max(alpha1_, scalar(0)), scalar(1))
+        clamp(alpha1_.boundaryField()[patchi], zero_one{})
     );
 
-    const scalarField& alpha1p = limitedAlpha1.boundaryField()[patchi];
-
     return (alpha1p*kappa1().value() + (1 - alpha1p)*kappa2().value());
 }
 
@@ -402,13 +394,11 @@ Foam::tmp<Foam::scalarField> Foam::twoPhaseMixtureEThermo::kappaEff
     const label patchi
 ) const
 {
-    const volScalarField limitedAlpha1
+    const scalarField alpha1p
     (
-        min(max(alpha1_, scalar(0)), scalar(1))
+        clamp(alpha1_.boundaryField()[patchi], zero_one{})
     );
 
-    const scalarField& alpha1p = limitedAlpha1.boundaryField()[patchi];
-
     return
         (alpha1p*kappa1().value() + (1 - alpha1p)*kappa2().value()) + kappat;
 
@@ -435,13 +425,11 @@ Foam::tmp<Foam::scalarField> Foam::twoPhaseMixtureEThermo::alphaEff
     const label patchi
 ) const
 {
-    const volScalarField limitedAlpha1
+    const scalarField alpha1p
     (
-        min(max(alpha1_, scalar(0)), scalar(1))
+        clamp(alpha1_.boundaryField()[patchi], zero_one{})
     );
 
-    const scalarField& alpha1p = limitedAlpha1.boundaryField()[patchi];
-
     const scalarField rho
     (
         alpha1p*rho1().value() + (1.0 - alpha1p)*rho2().value()

applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Kunz/Kunz.C

@@ -70,7 +70,7 @@ Foam::Pair<Foam::tmp<Foam::volScalarField>>
 Foam::phaseChangeTwoPhaseMixtures::Kunz::mDotAlphal() const
 {
     const volScalarField& p = alpha1_.db().lookupObject<volScalarField>("p");
-    volScalarField limitedAlpha1(min(max(alpha1_, scalar(0)), scalar(1)));
+    volScalarField limitedAlpha1(clamp(alpha1_, zero_one{}));
 
     return Pair<tmp<volScalarField>>
     (
@@ -85,7 +85,7 @@ Foam::Pair<Foam::tmp<Foam::volScalarField>>
 Foam::phaseChangeTwoPhaseMixtures::Kunz::mDotP() const
 {
     const volScalarField& p = alpha1_.db().lookupObject<volScalarField>("p");
-    volScalarField limitedAlpha1(min(max(alpha1_, scalar(0)), scalar(1)));
+    volScalarField limitedAlpha1(clamp(alpha1_, zero_one{}));
 
     return Pair<tmp<volScalarField>>
     (

applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/Merkle/Merkle.C

@@ -82,7 +82,7 @@ Foam::Pair<Foam::tmp<Foam::volScalarField>>
 Foam::phaseChangeTwoPhaseMixtures::Merkle::mDotP() const
 {
     const volScalarField& p = alpha1_.db().lookupObject<volScalarField>("p");
-    volScalarField limitedAlpha1(min(max(alpha1_, scalar(0)), scalar(1)));
+    volScalarField limitedAlpha1(clamp(alpha1_, zero_one{}));
 
     return Pair<tmp<volScalarField>>
     (

applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.C

@@ -99,7 +99,7 @@ Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::pCoeff
     const volScalarField& p
 ) const
 {
-    volScalarField limitedAlpha1(min(max(alpha1_, scalar(0)), scalar(1)));
+    volScalarField limitedAlpha1(clamp(alpha1_, zero_one{}));
     volScalarField rho
     (
         limitedAlpha1*rho1() + (scalar(1) - limitedAlpha1)*rho2()
@@ -117,7 +117,7 @@ Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::mDotAlphal() const
     const volScalarField& p = alpha1_.db().lookupObject<volScalarField>("p");
     volScalarField pCoeff(this->pCoeff(p));
 
-    volScalarField limitedAlpha1(min(max(alpha1_, scalar(0)), scalar(1)));
+    volScalarField limitedAlpha1(clamp(alpha1_, zero_one{}));
 
     return Pair<tmp<volScalarField>>
     (
@@ -134,7 +134,7 @@ Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::mDotP() const
     const volScalarField& p = alpha1_.db().lookupObject<volScalarField>("p");
     volScalarField pCoeff(this->pCoeff(p));
 
-    volScalarField limitedAlpha1(min(max(alpha1_, scalar(0)), scalar(1)));
+    volScalarField limitedAlpha1(clamp(alpha1_, zero_one{}));
     volScalarField apCoeff(limitedAlpha1*pCoeff);
 
     return Pair<tmp<volScalarField>>

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/vectorTools/vectorTools.H

@@ -114,7 +114,7 @@ namespace vectorTools
         scalar cosPhi = (a & b)/(mag(a)*mag(b) + tolerance);
 
         // Enforce bounding between -1 and 1
-        return min(max(cosPhi, -1), 1);
+        return clamp(cosPhi, -1, 1);
     }
 
     //- Calculate angle between a and b in radians
@@ -129,7 +129,7 @@ namespace vectorTools
         scalar cosPhi = (a & b)/(mag(a)*mag(b) + tolerance);
 
         // Enforce bounding between -1 and 1
-        return acos( min(max(cosPhi, -1), 1) );
+        return acos(clamp(cosPhi, -1, 1));
     }
 
     //- Calculate angle between a and b in degrees

src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.C

@@ -101,7 +101,7 @@ void Foam::adaptiveSolver::solve
     if (err > pow(maxScale_/safeScale_, -1.0/alphaInc_))
     {
         dxTry =
-            min(max(safeScale_*pow(err, -alphaInc_), minScale_), maxScale_)*dx;
+            clamp(safeScale_*pow(err, -alphaInc_), minScale_, maxScale_)*dx;
     }
     else
     {

src/OpenFOAM/meshes/meshShapes/face/face.C

@@ -76,7 +76,7 @@ Foam::label Foam::face::mostConcaveAngle
 
         // NOTE: is -ve angle since left edge pointing in other direction
         scalar edgeCos = (leftEdge & rightEdge);
-        scalar edgeAngle = acos(max(-1.0, min(1.0, edgeCos)));
+        scalar edgeAngle = acos(clamp(edgeCos, -1, 1));
 
         scalar angle;
 
@@ -205,7 +205,7 @@ Foam::label Foam::face::split
             splitEdge.normalise();
 
             const scalar splitCos = splitEdge & rightEdge;
-            const scalar splitAngle = acos(max(-1.0, min(1.0, splitCos)));
+            const scalar splitAngle = acos(clamp(splitCos, -1, 1));
             const scalar angleDiff = fabs(splitAngle - bisectAngle);
 
             if (angleDiff < minDiff)

src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshCheck.C

@@ -101,7 +101,7 @@ bool Foam::polyMesh::checkFaceOrthogonality
                         << " between cells " << own[facei]
                         << " and " << nei[facei]
                         << ": Angle = "
-                        << radToDeg(::acos(min(1.0, max(-1.0, ortho[facei]))))
+                        << radToDeg(::acos(clamp(ortho[facei], -1, 1)))
                         << " deg." << endl;
                 }
 
@@ -130,9 +130,9 @@ bool Foam::polyMesh::checkFaceOrthogonality
             if (debug || report)
             {
                 Info<< "    Mesh non-orthogonality Max: "
-                    << radToDeg(::acos(min(1.0, max(-1.0, minDDotS))))
+                    << radToDeg(::acos(clamp(minDDotS, -1, 1)))
                     << " average: "
-                    << radToDeg(::acos(min(1.0, max(-1.0, sumDDotS/nSummed))))
+                    << radToDeg(::acos(clamp(sumDDotS/nSummed, -1, 1)))
                     << endl;
             }
         }

src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H

@@ -334,7 +334,7 @@ inline Foam::scalar Foam::triangle<Point, PointRef>::circumRadius() const
     }
 
     const scalar a = (d1 + d2)*(d2 + d3)*(d3 + d1) / denom;
-    return 0.5*Foam::sqrt(min(GREAT, max(0, a)));
+    return 0.5*Foam::sqrt(clamp(a, 0, GREAT));
 }
 
 

src/OpenFOAM/primitives/functions/Function1/InputValueMapper/InputValueMapperI.H

@@ -47,7 +47,7 @@ inline Type Foam::Function1Types::InputValueMapper<Type>::value
         }
         case mappingMode::MINMAX:
         {
-            scalar tlim = min(max(t, min_), max_);
+            scalar tlim = clamp(t, min_, max_);
 
             return value_->value(tlim);
         }
@@ -86,8 +86,8 @@ Type Foam::Function1Types::InputValueMapper<Type>::integrate
         }
         case mappingMode::MINMAX:
         {
-            scalar xlim0 = min(max(x1, min_), max_);
+            scalar xlim0 = clamp(x1, min_, max_);
-            scalar xlim1 = min(max(x2, min_), max_);
+            scalar xlim1 = clamp(x2, min_, max_);
 
             Type intValue = value_->integrate(xlim0, xlim1);
 

src/OpenFOAM/primitives/triad/triad.C

@@ -395,7 +395,7 @@ Foam::scalar Foam::diff(const triad& A, const triad& B)
             (tmpA[dir] & tmpB[dir])
            /(Foam::mag(tmpA[dir])*Foam::mag(tmpA[dir]) + SMALL);
 
-        cosPhi = min(max(cosPhi, -1), 1);
+        cosPhi = clamp(cosPhi, -1, 1);
 
         sumDifference += mag(cosPhi - 1);
     }

src/TurbulenceModels/turbulenceModels/RAS/kOmegaSSTLM/kOmegaSSTLM.C

@@ -70,7 +70,7 @@ tmp<volScalarField::Internal> kOmegaSSTLM<BasicTurbulenceModel>::epsilonByk
 ) const
 {
     return
-        min(max(gammaIntEff_, scalar(0.1)), scalar(1))
+        clamp(gammaIntEff_, scalarMinMax(0.1, 1))
        *kOmegaSST<BasicTurbulenceModel>::epsilonByk(F1, gradU);
 }
 

src/TurbulenceModels/turbulenceModels/RAS/realizableKE/realizableKE.C

@@ -63,7 +63,7 @@ tmp<volScalarField> realizableKE<BasicTurbulenceModel>::rCmu
 
     volScalarField phis
     (
-        (1.0/3.0)*acos(min(max(sqrt(6.0)*W, -scalar(1)), scalar(1)))
+        (1.0/3.0)*acos(clamp(sqrt(6.0)*W, scalarMinMax(-1, 1)))
     );
     volScalarField As(sqrt(6.0)*cos(phis));
     volScalarField Us(sqrt(S2/2.0 + magSqr(skew(gradU))));

src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C

@@ -529,11 +529,7 @@ void Foam::edgeInterpolation::makeCorrectionVectors() const
     if (owner.size() > 0)
     {
         scalarField sinAlpha(deltaCoeffs*mag(CorrVecs.internalField()));
-
+        sinAlpha.clamp_range(-1, 1);
-        forAll(sinAlpha, edgeI)
-        {
-            sinAlpha[edgeI] = max(-1, min(sinAlpha[edgeI], 1));
-        }
 
         NonOrthogCoeff = max(Foam::asin(sinAlpha)*180.0/M_PI);
     }

src/finiteArea/interpolation/edgeInterpolation/schemes/Gamma/Gamma.H

@@ -110,7 +110,7 @@ public:
         gradcf = stabilise(gradcf, SMALL);
 
         scalar phict = 1 - 0.5*gradf/gradcf;
-        scalar limiter = min(max(phict/k_, 0), 1);
+        scalar limiter = clamp(phict/k_, 0, 1);
 
         return lerp(udWeight, cdWeight, limiter);
     }

src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/Gamma/Gamma.H

@@ -97,7 +97,7 @@ public:
             faceFlux, phiP, phiN, gradcP, gradcN, d
         );
 
-        return min(max(phict/k_, 0), 1);
+        return clamp(phict/k_, 0, 1);
     }
 };
 

src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/SFCD/SFCD.H

@@ -82,7 +82,7 @@ public:
             faceFlux, phiP, phiN, gradcP, gradcN, d
         );
 
-        scalar limitPhict = min(max(phict, 0), 0.5);
+        scalar limitPhict = clamp(phict, 0, 0.5);
         return limitPhict/(1 - limitPhict);
     }
 };

src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C

@@ -670,7 +670,7 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::processValues
             // Uniformity index
             const scalar ui = 1 - numer/(2*mag(mean*areaTotal) + ROOTVSMALL);
 
-            return min(max(ui, 0), 1);
+            return clamp(ui, 0, 1);
         }
 
         default:
@@ -756,7 +756,7 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::processValues
             // Uniformity index
             const scalar ui = 1 - numer/(2*mag(mean*areaTotal) + ROOTVSMALL);
 
-            return vector(min(max(ui, 0), 1), 0, 0);
+            return vector(clamp(ui, 0, 1), 0, 0);
         }
 
         default:

src/functionObjects/field/stabilityBlendingFactor/stabilityBlendingFactor.C

@@ -417,7 +417,7 @@ bool Foam::functionObjects::stabilityBlendingFactor::init(bool first)
             max
             (
                 indicator,
-                min(max(scalar(0), (Co - Co1_)/(Co2_ - Co1_)), scalar(1))
+                clamp((Co - Co1_)/(Co2_ - Co1_), zero_one{})
             );
 
         if (first)

src/functionObjects/field/stabilityBlendingFactor/stabilityBlendingFactor.H

@@ -152,7 +152,7 @@ Description
     For option 6, the following relation is used:
 
     \f[
-        fCoWeight = min(max((Co - Co1)/(Co2 - Co1), 0), 1)
+        fCoWeight = clamp((Co - Co1)/(Co2 - Co1), 0, 1)
     \f]
 
     where

src/fvOptions/sources/derived/solidificationMeltingSource/solidificationMeltingSource.C

@@ -147,7 +147,7 @@ void Foam::fv::solidificationMeltingSource::update(const volScalarField& Cp)
         scalar Cpc = Cp[celli];
         scalar alpha1New = alpha1_[celli] + relax_*Cpc*(Tc - Tmelt_)/L_;
 
-        alpha1_[celli] = max(0, min(alpha1New, 1));
+        alpha1_[celli] = clamp(alpha1New, 0, 1);
         deltaT_[i] = Tc - Tmelt_;
     }
 

src/lagrangian/distributionModels/multiNormal/multiNormal.C

@@ -168,7 +168,7 @@ Foam::scalar Foam::distributionModels::multiNormal::sample
     // Note: numerical approximation of the inverse function yields slight
     //       inaccuracies
 
-    return min(max(x, minValue_), maxValue_);
+    return clamp(x, minValue_, maxValue_);
 }
 
 

src/lagrangian/distributionModels/normal/normal.C

@@ -108,7 +108,7 @@ Foam::scalar Foam::distributionModels::normal::sample() const
     // Note: numerical approximation of the inverse function yields slight
     //       inaccuracies
 
-    return min(max(x, minValue_), maxValue_);
+    return clamp(x, minValue_, maxValue_);
 }
 
 

src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C

@@ -300,7 +300,7 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
 
     // Calculate the new temperature and the enthalpy transfer terms
     scalar Tnew = T_ + deltaT;
-    Tnew = min(max(Tnew, cloud.constProps().TMin()), cloud.constProps().TMax());
+    Tnew = clamp(Tnew, cloud.constProps().TMin(), cloud.constProps().TMax());
 
     dhsTrans -= m*Cp_*deltaTcp;
 

src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/DistortedSphereDrag/DistortedSphereDragForce.C

@@ -84,7 +84,7 @@ Foam::forceSuSp Foam::DistortedSphereDragForce<CloudType>::calcCoupled
 ) const
 {
     // Limit the drop distortion to y=0 (sphere) and y=1 (disk)
-    const scalar y = min(max(p.y(), scalar(0)), scalar(1));
+    const scalar y = clamp(p.y(), 0, 1);
 
     // (LMR:Eq. 10)
     return

src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C

@@ -361,9 +361,7 @@ Foam::scalar Foam::SprayParcel<ParcelType>::chi
         }
     }
 
-    chi = min(1.0, max(chi, 0.0));
+    return clamp(chi, 0, 1);
-
-    return chi;
 }
 
 

src/lagrangian/spray/submodels/BreakupModel/ETAB/ETAB.C

@@ -177,8 +177,7 @@ bool Foam::ETAB<CloudType>::update
                     Kbr =k2_*omega*sqrtWe;
                 }
 
-                scalar rWetmp = 1.0/Wetmp;
+                scalar cosdtbu = clamp((1.0 - 1.0/Wetmp), -1, 1);
-                scalar cosdtbu = max(-1.0, min(1.0, 1.0 - rWetmp));
                 scalar dtbu = romega*acos(cosdtbu);
                 scalar decay = exp(-Kbr*dtbu);
 

src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatch.C

@@ -326,8 +326,8 @@ void Foam::cyclicACMIPolyPatch::resetAMI(const UList<point>& points) const
     // Note:
     // - assumes that the non-overlap patches are decomposed using the same
     //   decomposition as the coupled patches (per side)
-    srcMask_ = min(scalar(1), max(scalar(0), AMI.srcWeightsSum()));
+    srcMask_ = clamp(AMI.srcWeightsSum(), zero_one{});
-    tgtMask_ = min(scalar(1), max(scalar(0), AMI.tgtWeightsSum()));
+    tgtMask_ = clamp(AMI.tgtWeightsSum(), zero_one{});
 
     if (debug)
     {

src/meshTools/searchableSurfaces/searchableCone/searchableCone.C

@@ -113,11 +113,11 @@ void Foam::searchableCone::findNearestAndNormal
     v.normalise();
 
     // Nearest and normal on disk at point1
-    point disk1Point(point1_ + min(max(magV, innerRadius1_), radius1_)*v);
+    point disk1Point(point1_ + clamp(magV, innerRadius1_, radius1_)*v);
     vector disk1Normal(-unitDir_);
 
     // Nearest and normal on disk at point2
-    point disk2Point(point2_ + min(max(magV, innerRadius2_), radius2_)*v);
+    point disk2Point(point2_ + clamp(magV, innerRadius2_, radius2_)*v);
     vector disk2Normal(unitDir_);
 
     // Nearest and normal on cone. Initialise to far-away point so if not

src/meshTools/triSurface/faceTriangulation/faceTriangulation.C

@@ -83,7 +83,7 @@ void Foam::faceTriangulation::calcHalfAngle
 )
 {
     // truncate cos to +-1 to prevent negative numbers
-    scalar cos = max(-1, min(1, e0 & e1));
+    scalar cos = clamp((e0 & e1), -1, 1);
 
     scalar sin = (e0 ^ e1) & normal;
 

src/phaseSystemModels/multiphaseInter/phasesSystem/MassTransferPhaseSystem/MassTransferPhaseSystem.C

@@ -600,11 +600,8 @@ void Foam::MassTransferPhaseSystem<BasePhaseSystem>::alphaTransfer
 
         if (includeDivU)
         {
-            SuPhase1 +=
+            SuPhase1 += fvc::div(phi)*clamp(alpha1, zero_one{});
-                fvc::div(phi)*min(max(alpha1, scalar(0)), scalar(1));
+            SuPhase2 += fvc::div(phi)*clamp(alpha2, zero_one{});
-
-            SuPhase2 +=
-                fvc::div(phi)*min(max(alpha2, scalar(0)), scalar(1));
         }
 
         // NOTE: dmdtNet is distributed in terms =

src/phaseSystemModels/multiphaseInter/phasesSystem/massTransferModels/Lee/Lee.C

@@ -57,7 +57,7 @@ Foam::meltingEvaporationModels::Lee<Thermo, OtherThermo>::Kexp
     {
         const volScalarField from
         (
-            min(max(this->pair().from(), scalar(0)), scalar(1))
+            clamp(this->pair().from(), zero_one{})
         );
 
         const volScalarField coeff
@@ -95,9 +95,9 @@ Foam::meltingEvaporationModels::Lee<Thermo, OtherThermo>::KSp
 {
     if (this->modelVariable_ == variable)
     {
-        volScalarField from
+        const volScalarField from
         (
-            min(max(this->pair().from(), scalar(0)), scalar(1))
+            clamp(this->pair().from(), zero_one{})
         );
 
         const volScalarField coeff
@@ -123,7 +123,7 @@ Foam::meltingEvaporationModels::Lee<Thermo, OtherThermo>::KSp
     }
     else
     {
-        return tmp<volScalarField> ();
+        return nullptr;
     }
 }
 
@@ -138,9 +138,9 @@ Foam::meltingEvaporationModels::Lee<Thermo, OtherThermo>::KSu
 {
     if (this->modelVariable_ == variable)
     {
-        volScalarField from
+        const volScalarField from
         (
-            min(max(this->pair().from(), scalar(0)), scalar(1))
+            clamp(this->pair().from(), zero_one{})
         );
 
         const volScalarField coeff
@@ -165,7 +165,7 @@ Foam::meltingEvaporationModels::Lee<Thermo, OtherThermo>::KSu
     }
     else
     {
-        return tmp<volScalarField> ();
+        return nullptr;
     }
 }
 

src/phaseSystemModels/multiphaseInter/phasesSystem/multiphaseInterSystem/multiphaseInterSystem.C

@@ -1277,12 +1277,12 @@ Foam::tmp<Foam::volVectorField> Foam::multiphaseInterSystem::nVolHatfv
 {
     const volScalarField alpha1m
     (
-        min(max(alpha1, scalar(0)), scalar(1))
+        clamp(alpha1, zero_one{})
     );
 
     const volScalarField alpha2m
     (
-        min(max(alpha2, scalar(0)), scalar(1))
+        clamp(alpha2, zero_one{})
     );
 
     const volVectorField gradAlphaf
@@ -1311,12 +1311,12 @@ Foam::tmp<Foam::surfaceVectorField> Foam::multiphaseInterSystem::nHatfv
 
     const volScalarField alpha1b
     (
-        min(max(alpha1, scalar(0)), scalar(1))
+        clamp(alpha1, zero_one{})
     );
 
     const volScalarField alpha2b
     (
-        min(max(alpha2, scalar(0)), scalar(1))
+        clamp(alpha2, zero_one{})
     );
 
     surfaceVectorField gradAlphaf

src/phaseSystemModels/multiphaseInter/phasesSystem/multiphaseSystem/multiphaseSystem.C

@@ -286,8 +286,7 @@ void Foam::multiphaseInter::multiphaseSystem::solveAlphas()
 
             // Add alpha*div(U)
             //const volScalarField& alpha = phase;
-            //Su_[phase.name()] +=
+            //Su_[phase.name()] += fvc::div(phi)*clamp(alpha, zero_one{});
-            //    fvc::div(phi)*min(max(alpha, scalar(0)), scalar(1));
         }
 
         // Fill Su and Sp

src/phaseSystemModels/reactingEuler/multiphaseSystem/interfacialModels/virtualMassModels/Lamb/Lamb.C

@@ -69,7 +69,7 @@ Foam::virtualMassModels::Lamb::~Lamb()
 
 Foam::tmp<Foam::volScalarField> Foam::virtualMassModels::Lamb::Cvm() const
 {
-    volScalarField E(min(max(pair_.E(), SMALL), 1 - SMALL));
+    volScalarField E(clamp(pair_.E(), scalarMinMax(SMALL, 1 - SMALL)));
     volScalarField rtOmEsq(sqrt(1 - sqr(E)));
 
     return

src/phaseSystemModels/reactingEuler/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/radialModel/SinclairJackson/SinclairJacksonRadial.C

@@ -90,7 +90,11 @@ Foam::kineticTheoryModels::radialModels::SinclairJackson::g0prime
 {
     volScalarField aByaMax
     (
-        cbrt(min(max(alpha, scalar(1e-3)), alphaMinFriction)/alphaMax)
+        cbrt
+        (
+            clamp(alpha, scalarMinMax(1e-3, alphaMinFriction.value()))
+          / alphaMax
+        )
     );
 
     return (1.0/(3*alphaMax))/sqr(aByaMax - sqr(aByaMax));

src/phaseSystemModels/twoPhaseEuler/phaseCompressibleTurbulenceModels/kineticTheoryModels/radialModel/SinclairJackson/SinclairJacksonRadial.C

@@ -90,7 +90,11 @@ Foam::kineticTheoryModels::radialModels::SinclairJackson::g0prime
 {
     volScalarField aByaMax
     (
-        cbrt(min(max(alpha, scalar(1e-3)), alphaMinFriction)/alphaMax)
+        cbrt
+        (
+            clamp(alpha, scalarMinMax(1e-3, alphaMinFriction.value()))
+          / alphaMax
+        )
     );
 
     return (1.0/(3*alphaMax))/sqr(aByaMax - sqr(aByaMax));

src/phaseSystemModels/twoPhaseEuler/twoPhaseSystem/interfacialModels/virtualMassModels/Lamb/Lamb.C

@@ -69,7 +69,7 @@ Foam::virtualMassModels::Lamb::~Lamb()
 
 Foam::tmp<Foam::volScalarField> Foam::virtualMassModels::Lamb::Cvm() const
 {
-    volScalarField E(min(max(pair_.E(), SMALL), 1 - SMALL));
+    volScalarField E(clamp(pair_.E(), scalarMinMax(SMALL, 1 - SMALL)));
     volScalarField rtOmEsq(sqrt(1 - sqr(E)));
 
     return

src/sixDoFRigidBodyMotion/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C

@@ -233,7 +233,7 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
         g_ = g.value();
     }
 
-    // scalar ramp = min(max((t.value() - 5)/10, 0), 1);
+    // scalar ramp = clamp((t.value() - 5)/10, 0, 1);
     scalar ramp = 1.0;
 
     motion_.update

src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.C

@@ -219,7 +219,7 @@ void Foam::sixDoFRigidBodyMotionSolver::solve()
         coeffDict().readIfPresent("g", g);
     }
 
-    // const scalar ramp = min(max((this->db().time().value() - 5)/10, 0), 1);
+    // const scalar ramp = clamp((this->db().time().value() - 5)/10, 0, 1);
     const scalar ramp = 1.0;
 
     if (test_)

src/thermoTools/derivedFvPatchFields/outletMappedUniformInletHeatAddition/outletMappedUniformInletHeatAdditionFvPatchField.C

@@ -170,7 +170,7 @@ void Foam::outletMappedUniformInletHeatAdditionFvPatchField::updateCoeffs()
 
         scalar totalPhiCp = gSum(outletPatchPhi)*gAverage(Cpf);
 
-        operator==(min(max(averageOutletField + Q_/totalPhiCp, TMin_), TMax_));
+        operator==(clamp(averageOutletField + Q_/totalPhiCp, TMin_, TMax_));
     }
     else
     {

src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H

@@ -123,12 +123,10 @@ inline Foam::scalar Foam::janafThermo<EquationOfState>::limit
             << Tlow_ << " -> " << Thigh_ << ";  T = " << T
             << endl;
 
-        return min(max(T, Tlow_), Thigh_);
+        return clamp(T, Tlow_, Thigh_);
-    }
-    else
-    {
-        return T;
     }
+
+    return T;
 }
 
 

src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvection.C

@@ -631,7 +631,7 @@ void Foam::isoAdvection::applyBruteForceBounding()
 
     if (dict_.getOrDefault("clip", true))
     {
-        alpha1_ = min(scalar(1), max(scalar(0), alpha1_));
+        alpha1_.clamp_range(zero_one{});
         alpha1Changed = true;
     }
 

src/transportModels/incompressible/incompressibleTwoPhaseMixture/incompressibleTwoPhaseMixture.C

@@ -48,7 +48,7 @@ void Foam::incompressibleTwoPhaseMixture::calcNu()
     const volScalarField limitedAlpha1
     (
         "limitedAlpha1",
-        min(max(alpha1_, scalar(0)), scalar(1))
+        clamp(alpha1_, zero_one{})
     );
 
     // Average kinematic viscosity calculated from dynamic viscosity
@@ -128,7 +128,7 @@ Foam::incompressibleTwoPhaseMixture::mu() const
 {
     const volScalarField limitedAlpha1
     (
-        min(max(alpha1_, scalar(0)), scalar(1))
+        clamp(alpha1_, zero_one{})
     );
 
     return tmp<volScalarField>::New
@@ -153,7 +153,7 @@ Foam::incompressibleTwoPhaseMixture::muf() const
 {
     const surfaceScalarField alpha1f
     (
-        min(max(fvc::interpolate(alpha1_), scalar(0)), scalar(1))
+        clamp(fvc::interpolate(alpha1_), zero_one{})
     );
 
     return tmp<surfaceScalarField>::New
@@ -170,7 +170,7 @@ Foam::incompressibleTwoPhaseMixture::nuf() const
 {
     const surfaceScalarField alpha1f
     (
-        min(max(fvc::interpolate(alpha1_), scalar(0)), scalar(1))
+        clamp(fvc::interpolate(alpha1_), zero_one{})
     );
 
     return tmp<surfaceScalarField>::New

src/transportModels/interfaceProperties/interfaceCompression/interfaceCompression.H

@@ -71,10 +71,11 @@ public:
         //return min(max(4*min(phiP*(1 - phiP), phiN*(1 - phiN)), 0), 1);
 
         // Quartic compression scheme
-        return
+        return clamp
-            min(max(
+        (
             1 - max(sqr(1 - 4*phiP*(1 - phiP)), sqr(1 - 4*phiN*(1 - phiN))),
-            0), 1);
+            0, 1
+        );
     }
 };
 

src/waveModels/derivedPointPatchFields/waveMaker/waveMakerPointPatchVectorField.C

@@ -89,7 +89,7 @@ Foam::scalar Foam::waveMakerPointPatchVectorField::timeCoeff
     const scalar t
 ) const
 {
-    return max(0, min(t/rampTime_, 1));
+    return clamp(t/rampTime_, 0, 1);
 }
 
 

src/waveModels/waveGenerationModels/base/irregularWaveModel/irregularWaveModel.C

@@ -46,7 +46,7 @@ Foam::scalar Foam::waveModels::irregularWaveModel::timeCoeff
     const scalar t
 ) const
 {
-    return max(0, min(t/rampTime_, 1));
+    return clamp(t/rampTime_, 0, 1);
 }