FreeFOAM The Cross-Platform CFD Toolkit
Random_.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 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
13  the Free Software Foundation, either version 3 of the License, or
14  (at your 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, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include <OpenFOAM/Random.H>
27 
28 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
29 
30 namespace Foam
31 {
32 
33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
34 
35 #if INT_MAX != 2147483647
36 # error "INT_MAX != 2147483647"
37 # error "The random number generator may not work!"
38 #endif
39 
40 #ifdef USE_RANDOM
41 # include <climits>
42 #else
43 # include <cstdlib>
44 #endif
45 
46 
47 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
48 
49 // construct given seed
50 Random::Random(const label& seed)
51 {
52  if (seed > 1)
53  {
54  Seed = seed;
55  }
56  else
57  {
58  Seed = 1;
59  }
60 
61 # ifdef USE_RANDOM
62  srandom((unsigned int)Seed);
63 # else
64  srand48(Seed);
65 # endif
66 
67 }
68 
69 
71 {
72 # ifdef USE_RANDOM
73  if (random() > INT_MAX/2)
74 # else
75  if (lrand48() > INT_MAX/2)
76 # endif
77  {
78  return 1;
79  }
80  else
81  {
82  return 0;
83  }
84 }
85 
86 
88 {
89 # ifdef USE_RANDOM
90  return (scalar)random()/INT_MAX;
91 # else
92  return drand48();
93 # endif
94 }
95 
96 
98 {
99  vector rndVec;
100  for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
101  {
102  rndVec.component(cmpt) = scalar01();
103  }
104 
105  return rndVec;
106 }
107 
108 
110 {
111  sphericalTensor rndTen;
112  rndTen.ii() = scalar01();
113 
114  return rndTen;
115 }
116 
117 
119 {
120  symmTensor rndTen;
121  for (direction cmpt=0; cmpt<symmTensor::nComponents; cmpt++)
122  {
123  rndTen.component(cmpt) = scalar01();
124  }
125 
126  return rndTen;
127 }
128 
129 
131 {
132  tensor rndTen;
133  for (direction cmpt=0; cmpt<tensor::nComponents; cmpt++)
134  {
135  rndTen.component(cmpt) = scalar01();
136  }
137 
138  return rndTen;
139 }
140 
141 
142 label Random::integer(const label lower, const label upper)
143 {
144 # ifdef USE_RANDOM
145  return lower + (random() % (upper+1-lower));
146 # else
147  return lower + (lrand48() % (upper+1-lower));
148 # endif
149 }
150 
151 
152 vector Random::position(const vector& start, const vector& end)
153 {
154  vector rndVec(start);
155 
156  for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
157  {
158  rndVec.component(cmpt) +=
159  scalar01()*(end.component(cmpt) - start.component(cmpt));
160  }
161 
162  return rndVec;
163 }
164 
165 
166 void Random::randomise(scalar& s)
167 {
168  s = scalar01();
169 }
170 
171 
173 {
174  v = vector01();
175 }
176 
177 
179 {
180  st = sphericalTensor01();
181 }
182 
183 
185 {
186  st = symmTensor01();
187 }
188 
189 
191 {
192  t = tensor01();
193 }
194 
195 
196 // return a normal Gaussian randon number
197 // with zero mean and unity variance N(0, 1)
198 
200 {
201  static int iset = 0;
202  static scalar gset;
203  scalar fac, rsq, v1, v2;
204 
205  if (iset == 0)
206  {
207  do
208  {
209  v1 = 2.0*scalar01() - 1.0;
210  v2 = 2.0*scalar01() - 1.0;
211  rsq = v1*v1 + v2*v2;
212  } while (rsq >= 1.0 || rsq == 0.0);
213 
214  fac = sqrt(-2.0 * log(rsq)/rsq);
215  gset = v1*fac;
216  iset = 1;
217 
218  return v2*fac;
219  }
220  else
221  {
222  iset = 0;
223 
224  return gset;
225  }
226 }
227 
228 
229 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
230 
231 } // End namespace Foam
232 
233 // ************************ vim: set sw=4 sts=4 et: ************************ //