FreeFOAM The Cross-Platform CFD Toolkit
general.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 "general.H"
28 
29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
30 
31 namespace Foam
32 {
33  namespace pdfs
34  {
35  defineTypeNameAndDebug(general, 0);
36  addToRunTimeSelectionTable(pdf, general, dictionary);
37  }
38 }
39 
40 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
41 
43 :
44  pdf(typeName, dict, rndGen),
45  xy_(pdfDict_.lookup("distribution")),
46  nEntries_(xy_.size()),
47  minValue_(xy_[0][0]),
48  maxValue_(xy_[nEntries_-1][0]),
49  integral_(nEntries_)
50 {
51  check();
52 
53  // normalize the cumulative pdf
54 
55  integral_[0] = 0.0;
56  for (label i=1; i<nEntries_; i++)
57  {
58 
59  scalar k = (xy_[i][1] - xy_[i-1][1])/(xy_[i][0] - xy_[i-1][0]);
60  scalar d = xy_[i-1][1] - k*xy_[i-1][0];
61  scalar y1 = xy_[i][0]*(0.5*k*xy_[i][0] + d);
62  scalar y0 = xy_[i-1][0]*(0.5*k*xy_[i-1][0] + d);
63  scalar area = y1 - y0;
64 
65  integral_[i] = area + integral_[i-1];
66  }
67 
68  for (label i=0; i<nEntries_; i++)
69  {
70  xy_[i][1] /= integral_[nEntries_-1];
71  integral_[i] /= integral_[nEntries_-1];
72  }
73 
74 }
75 
76 
77 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
78 
80 {}
81 
82 
83 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
84 
85 Foam::scalar Foam::pdfs::general::sample() const
86 {
87  scalar y = rndGen_.scalar01();
88 
89  // find the interval where y is in the table
90  label n=1;
91  while (integral_[n] <= y)
92  {
93  n++;
94  }
95 
96  scalar k = (xy_[n][1] - xy_[n-1][1])/(xy_[n][0] - xy_[n-1][0]);
97  scalar d = xy_[n-1][1] - k*xy_[n-1][0];
98 
99  scalar alpha = y + xy_[n-1][0]*(0.5*k*xy_[n-1][0] + d) - integral_[n-1];
100  scalar x = 0.0;
101 
102  // if k is small it is a linear equation, otherwise it is of second order
103  if (mag(k) > SMALL)
104  {
105  scalar p = 2.0*d/k;
106  scalar q = -2.0*alpha/k;
107  scalar sqrtEr = sqrt(0.25*p*p - q);
108 
109  scalar x1 = -0.5*p + sqrtEr;
110  scalar x2 = -0.5*p - sqrtEr;
111  if ((x1 >= xy_[n-1][0]) && (x1 <= xy_[n][0]))
112  {
113  x = x1;
114  }
115  else
116  {
117  x = x2;
118  }
119  }
120  else
121  {
122  x = alpha/d;
123  }
124 
125  return x;
126 }
127 
128 
129 Foam::scalar Foam::pdfs::general::minValue() const
130 {
131  return minValue_;
132 }
133 
134 
135 Foam::scalar Foam::pdfs::general::maxValue() const
136 {
137  return maxValue_;
138 }
139 
140 
141 // ************************ vim: set sw=4 sts=4 et: ************************ //