src/OpenFOAM/primitives/random/Random.C

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