36namespace distributionModels
46void Foam::distributionModels::general::initialise()
48 static scalar eps = ROOTVSMALL;
57 for (label i = 1; i < nEntries_; ++i)
59 scalar
k = (xy_[i][1] - xy_[i-1][1])/(xy_[i][0] - xy_[i-1][0] + eps);
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;
67 integral_[i] = xy_[i][1];
72 integral_[i] =
area + integral_[i-1];
74 y1 =
sqr(xy_[i][0])*(1.0/3.0*
k*xy_[i][0] + 0.5*d);
75 y0 =
sqr(xy_[i-1][0])*(1.0/3.0*
k*xy_[i-1][0] + 0.5*d);
76 meanValue_ +=
y1 -
y0;
81 const scalar sumArea = integral_.
last();
83 for (label i = 0; i < nEntries_; ++i)
85 xy_[i][1] /= sumArea + eps;
86 integral_[i] /= sumArea + eps;
89 meanValue_ /= sumArea + eps;
90 meanValue_ = cumulative_ ? (
maxValue_ - meanValue_) : meanValue_;
103 xy_(distributionModelDict_.
lookup(
"distribution")),
104 nEntries_(xy_.size()),
106 integral_(nEntries_),
107 cumulative_(distributionModelDict_.getOrDefault(
"cumulative", false))
115 if (cumulative_ && xy_[0][1] != 0)
118 <<
type() <<
"distribution: "
119 <<
"Elements in the second column for cumulative "
120 <<
"distribution functions must start from zero." <<
nl
121 <<
"First element = " << xy_[0][1]
125 for (label i = 0; i < nEntries_; ++i)
127 if (i > 0 && xy_[i][0] <= xy_[i-1][0])
130 <<
type() <<
"distribution: "
131 <<
"Elements in the first column must "
132 <<
"be specified in an ascending order." <<
nl
133 <<
"Please see the row i = " << i <<
nl
134 <<
"xy[i-1] = " << xy_[i-1] <<
nl
135 <<
"xy[i] = " << xy_[i]
139 if (xy_[i][0] < 0 || xy_[i][1] < 0)
142 <<
type() <<
"distribution: "
143 <<
"Input pairs cannot contain any negative element." <<
nl
144 <<
"Please see the row i = " << i <<
nl
145 <<
"xy[i] = " << xy_[i]
149 if (cumulative_ && i > 0 && xy_[i][1] < xy_[i-1][1])
152 <<
type() <<
"distribution: "
153 <<
"Elements in the second column for cumulative "
154 <<
"distribution functions must be non-decreasing." <<
nl
155 <<
"Please see the row i = " << i <<
nl
156 <<
"xy[i-1] = " << xy_[i-1] <<
nl
157 <<
"xy[i] = " << xy_[i]
169 const scalar binWidth,
190 nEntries_ = bin1 - bin0;
195 <<
"Data cannot be binned - zero bins generated" <<
nl
196 <<
" Bin width : " << binWidth <<
nl
197 <<
" Sample data : " << sampleData
203 xy_.setSize(nEntries_);
206 for (label bini = 0; bini < nEntries_; ++bini)
208 xy_[bini][0] = (bin0 + bini)*binWidth;
215 label bini = floor(sampleData[i]/binWidth) - bin0;
227 nEntries_(
p.nEntries_),
228 meanValue_(
p.meanValue_),
229 integral_(
p.integral_),
230 cumulative_(
p.cumulative_)
238 const scalar u = rndGen_.sample01<scalar>();
242 while (integral_[
n] <= u)
247 const scalar
k = (xy_[
n][1] - xy_[
n-1][1])/(xy_[
n][0] - xy_[
n-1][0]);
248 const scalar d = xy_[
n-1][1] -
k*xy_[
n-1][0];
256 u + xy_[
n-1][0]*(0.5*
k*xy_[
n-1][0] + d) - integral_[
n-1];
261 const scalar
p = 2.0*d/
k;
262 const scalar q = -2.0*
alpha/
k;
263 const scalar sqrtEr =
sqrt(0.25*
p*
p - q);
265 const scalar x1 = -0.5*
p + sqrtEr;
266 const scalar x2 = -0.5*
p - sqrtEr;
267 if ((x1 >= xy_[
n-1][0]) && (x1 <= xy_[
n][0]))
322 xy_.setSize(
x.size());
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
virtual bool check(const char *operation) const
Check IOstream status for given operation.
An Istream is an abstract base class for all input systems (streams, files, token lists etc)....
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
void setSize(const label n)
Alias for resize()
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
T & last()
Return the last element of the list.
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
A library of runtime-selectable doubly-truncated probability distribution models. Returns random samp...
virtual void check() const
Check that the distribution model is valid.
scalar minValue_
Minimum of the distribution.
scalar maxValue_
Maximum of the distribution.
Particle-size distribution model wherein random samples are drawn from a given arbitrary probability ...
virtual scalar meanValue() const
Return the arithmetic mean of the distribution data.
virtual tmp< scalarField > y() const
Probabilities.
virtual tmp< scalarField > x() const
Bin boundaries.
virtual void writeData(Ostream &os) const
Write data to stream.
virtual scalar sample() const
Sample the distribution.
virtual void readDict(const dictionary &dict)
Read data from dictionary.
virtual void readData(Istream &os)
Read data from stream.
Foam::dictionary writeDict() const
Write to dictionary.
Lookup type of boundary radiation properties.
A class for managing temporary objects.
A class for handling words, derived from Foam::string.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
OBJstream os(runTime.globalPath()/outputName)
const word dictName("faMeshDefinition")
#define WarningInFunction
Report a warning using Foam::Warning.
const wordList area
Standard area field types (scalar, vector, tensor, etc)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedScalar y0(const dimensionedScalar &ds)
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensionedScalar y1(const dimensionedScalar &ds)
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Istream & operator>>(Istream &, directionInfo &)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
errorManipArg< error, int > exit(error &err, const int errNo=1)
constexpr char nl
The newline '\n' character (0x0a)
#define forAll(list, i)
Loop across all elements in list.