src/foam/primitives/random/Random.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "Random.H"
  27
  28 #ifdef mingw
  29 #   include "rand48.h"
  30 #endif
  31
  32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  33
  34 namespace Foam
  35 {
  36
  37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  38
  39 #if INT_MAX    != 2147483647
  40 #    error "INT_MAX    != 2147483647"
  41 #    error "The random number generator may not work!"
  42 #endif
  43
  44 #ifdef USE_RANDOM
  45 #   include <climits>
  46 #else
  47 #   include <cstdlib>
  48 #endif
  49
  50
  51 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  52
  53 // construct given seed
  54 Random::Random(const label& seed)
  55 {
  56     if (seed > 1)
  57     {
  58         Seed = seed;
  59     }
  60     else
  61     {
  62         Seed = 1;
  63     }
  64
  65 #   ifdef USE_RANDOM
  66         srandom((unsigned int)Seed);
  67 #   else
  68         srand48(Seed);
  69 #   endif
  70
  71 }
  72
  73
  74 int Random::bit()
  75 {
  76 #   ifdef USE_RANDOM
  77     if (random() > INT_MAX/2)
  78 #   else
  79     if (lrand48() > INT_MAX/2)
  80 #   endif
  81     {
  82         return 1;
  83     }
  84     else
  85     {
  86         return 0;
  87     }
  88 }
  89
  90
  91 scalar Random::scalar01()
  92 {
  93 #   ifdef USE_RANDOM
  94         return (scalar)random()/INT_MAX;
  95 #   else
  96         return drand48();
  97 #   endif
  98 }
  99
 100
 101 vector Random::vector01()
 102 {
 103     vector rndVec;
 104     for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
 105     {
 106         rndVec.component(cmpt) = scalar01();
 107     }
 108
 109     return rndVec;
 110 }
 111
 112
 113 sphericalTensor Random::sphericalTensor01()
 114 {
 115     sphericalTensor rndTen;
 116     rndTen.ii() = scalar01();
 117
 118     return rndTen;
 119 }
 120
 121
 122 symmTensor Random::symmTensor01()
 123 {
 124     symmTensor rndTen;
 125     for (direction cmpt=0; cmpt<symmTensor::nComponents; cmpt++)
 126     {
 127         rndTen.component(cmpt) = scalar01();
 128     }
 129
 130     return rndTen;
 131 }
 132
 133
 134 symmTensor4thOrder Random::symmTensor4thOrder01()
 135 {
 136   symmTensor4thOrder rndTen;
 137   for (direction cmpt=0; cmpt<symmTensor4thOrder::nComponents; cmpt++)
 138     {
 139       rndTen.component(cmpt) = scalar01();
 140     }
 141
 142   return rndTen;
 143 }
 144
 145
 146 diagTensor Random::diagTensor01()
 147 {
 148   diagTensor rndTen;
 149   for (direction cmpt=0; cmpt<diagTensor::nComponents; cmpt++)
 150     {
 151       rndTen.component(cmpt) = scalar01();
 152     }
 153
 154   return rndTen;
 155 }
 156
 157
 158 tensor Random::tensor01()
 159 {
 160     tensor rndTen;
 161     for (direction cmpt=0; cmpt<tensor::nComponents; cmpt++)
 162     {
 163         rndTen.component(cmpt) = scalar01();
 164     }
 165
 166     return rndTen;
 167 }
 168
 169
 170 label Random::integer(const label lower, const label upper)
 171 {
 172 #   ifdef USE_RANDOM
 173         return lower + (random() % (upper+1-lower));
 174 #   else
 175         return lower + (lrand48() % (upper+1-lower));
 176 #   endif
 177 }
 178
 179
 180 vector Random::position(const vector& start, const vector& end)
 181 {
 182     vector rndVec(start);
 183
 184     for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
 185     {
 186         rndVec.component(cmpt) +=
 187             scalar01()*(end.component(cmpt) - start.component(cmpt));
 188     }
 189
 190     return rndVec;
 191 }
 192
 193
 194 void Random::randomise(scalar& s)
 195 {
 196      s = scalar01();
 197 }
 198
 199
 200 void Random::randomise(vector& v)
 201 {
 202     v = vector01();
 203 }
 204
 205
 206 void Random::randomise(sphericalTensor& st)
 207 {
 208     st = sphericalTensor01();
 209 }
 210
 211
 212 void Random::randomise(symmTensor& st)
 213 {
 214     st = symmTensor01();
 215 }
 216
 217
 218 void Random::randomise(symmTensor4thOrder& st)
 219 {
 220     st = symmTensor4thOrder01();
 221 }
 222
 223
 224 void Random::randomise(diagTensor& dt)
 225 {
 226     dt = diagTensor01();
 227 }
 228
 229
 230 void Random::randomise(tensor& t)
 231 {
 232     t = tensor01();
 233 }
 234
 235
 236 // return a normal Gaussian randon number
 237 // with zero mean and unity variance N(0, 1)
 238
 239 scalar Random::GaussNormal()
 240 {
 241     static int iset = 0;
 242     static scalar gset;
 243     scalar fac, rsq, v1, v2;
 244
 245     if (iset == 0)
 246     {
 247         do
 248         {
 249             v1 = 2.0*scalar01() - 1.0;
 250             v2 = 2.0*scalar01() - 1.0;
 251             rsq = v1*v1 + v2*v2;
 252         } while (rsq >= 1.0 || rsq == 0.0);
 253
 254         fac = sqrt(-2.0 * log(rsq)/rsq);
 255         gset = v1*fac;
 256         iset = 1;
 257
 258         return v2*fac;
 259     }
 260     else
 261     {
 262         iset = 0;
 263
 264         return gset;
 265     }
 266 }
 267
 268
 269 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 270
 271 } // End namespace Foam
 272
 273 // ************************************************************************* //