Merge remote-tracking branch 'origin/master' into develop

src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.C

@@ -615,8 +615,8 @@ void Foam::isoAdvection::limitFluxes()
     const scalar aTol = 1.0e-12;          // Note: tolerances
     const scalar maxAlphaMinus1 = 1;      // max(alphaNew - 1);
     const scalar minAlpha = -1;           // min(alphaNew);
-    const label nUndershoots = 20;        // sum(neg(alphaNew + aTol));
-    const label nOvershoots = 20;         // sum(pos(alphaNew - 1 - aTol));
+    const label nUndershoots = 20;        // sum(neg0(alphaNew + aTol));
+    const label nOvershoots = 20;         // sum(pos0(alphaNew - 1 - aTol));
     cellIsBounded_ = false;
 
     // Loop number of bounding steps
@@ -682,8 +682,8 @@ void Foam::isoAdvection::limitFluxes()
             scalarField alphaNew(alpha1In_ - fvc::surfaceIntegrate(dVf_)());
             label maxAlphaMinus1 = max(alphaNew - 1);
             scalar minAlpha = min(alphaNew);
-            label nUndershoots = sum(neg(alphaNew + aTol));
-            label nOvershoots = sum(pos(alphaNew - 1 - aTol));
+            label nUndershoots = sum(neg0(alphaNew + aTol));
+            label nOvershoots = sum(pos0(alphaNew - 1 - aTol));
             Info<< "After bounding number " << n + 1 << " of time "
                 << mesh_.time().value() << ":" << endl;
             Info<< "nOvershoots = " << nOvershoots << " with max(alphaNew-1) = "
@@ -792,7 +792,7 @@ void Foam::isoAdvection::boundFromAbove
                         fluidToPassOn*mag(phi[fi]*dt)/dVftot;
 
                     nFacesToPassFluidThrough +=
-                        pos(dVfmax[fi] - fluidToPassThroughFace);
+                        pos0(dVfmax[fi] - fluidToPassThroughFace);
 
                     fluidToPassThroughFace =
                         min(fluidToPassThroughFace, dVfmax[fi]);
@@ -1014,9 +1014,9 @@ void Foam::isoAdvection::applyBruteForceBounding()
     {
         alpha1_ =
             alpha1_
-           *pos(alpha1_ - snapAlphaTol)
-           *neg(alpha1_ - (1.0 - snapAlphaTol))
-          + pos(alpha1_ - (1.0 - snapAlphaTol));
+           *pos0(alpha1_ - snapAlphaTol)
+           *neg0(alpha1_ - (1.0 - snapAlphaTol))
+          + pos0(alpha1_ - (1.0 - snapAlphaTol));
 
         alpha1Changed = true;
     }

src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutFace/isoCutFace.C

@@ -362,7 +362,7 @@ Foam::scalar Foam::isoCutFace::timeIntegratedArea
     {
         // If all face cuttings were in the past and cell is filling up (Un0>0)
         // then face must be full during whole time interval
-        tIntArea = magSf*dt*pos(Un0);
+        tIntArea = magSf*dt*pos0(Un0);
         return tIntArea;
     }
 
@@ -373,7 +373,7 @@ Foam::scalar Foam::isoCutFace::timeIntegratedArea
         // If all cuttings are in the future but non of them within [0,dt] then
         // if cell is filling up (Un0 > 0) face must be empty during whole time
         // interval
-        tIntArea = magSf*dt*(1 - pos(Un0));
+        tIntArea = magSf*dt*(1 - pos0(Un0));
         return tIntArea;
     }
 
@@ -398,7 +398,7 @@ Foam::scalar Foam::isoCutFace::timeIntegratedArea
         // If Un0 > 0 cell is filling up and it must initially be empty.
         // If Un0 < 0 cell must initially be full(y immersed in fluid A).
         time = firstTime;
-        initialArea = magSf*(1.0 - pos(Un0));
+        initialArea = magSf*(1.0 - pos0(Un0));
         tIntArea = initialArea*time;
         cutPoints(fPts, pTimes, time, FIIL);
     }
@@ -470,7 +470,7 @@ Foam::scalar Foam::isoCutFace::timeIntegratedArea
     {
         // FIIL will leave the face at lastTime and face will be fully in fluid
         // A or fluid B in the time interval from lastTime to dt.
-        tIntArea += magSf*(dt - lastTime)*pos(Un0);
+        tIntArea += magSf*(dt - lastTime)*pos0(Un0);
     }
 
     return tIntArea;

src/lagrangian/distributionModels/binned/binned.C

@@ -47,10 +47,6 @@ const char* Foam::distributionModels::binned::header =
 void Foam::distributionModels::binned::initialise()
 {
     const label nSample(xy_.size());
-    forAll(xy_, bini)
-    {
-        xy_[bini][1] /= scalar(nSample);
-    }
 
     // Convert values to integral values
     for (label bini = 1; bini < nSample; ++bini)
@@ -58,6 +54,13 @@ void Foam::distributionModels::binned::initialise()
         xy_[bini][1] += xy_[bini - 1][1];
     }
 
+    // Normalise
+    scalar sumProb = xy_.last()[1];
+    forAll(xy_, bini)
+    {
+        xy_[bini][1] /= sumProb;
+    }
+
     // Calculate the mean value
     label bini = 0;
     forAll(xy_, i)
@@ -188,16 +191,7 @@ Foam::scalar Foam::distributionModels::binned::sample() const
     {
         if (xy_[i][1] > y)
         {
-            scalar d1 = y - xy_[i][1];
-            scalar d2 = xy_[i+1][1] - y;
-            if (d1 < d2)
-            {
-                return xy_[i][0];
-            }
-            else
-            {
-                return xy_[i+1][0];
-            }
+            return xy_[i][0];
         }
     }