src/lagrangian/distributionModels/general/general.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   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"
  27 #include "addToRunTimeSelectionTable.H"
  28
  29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  30
  31 namespace Foam
  32 {
  33     namespace distributionModels
  34     {
  35         defineTypeNameAndDebug(general, 0);
  36         addToRunTimeSelectionTable(distributionModel, general, dictionary);
  37     }
  38 }
  39
  40 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  41
  42 Foam::distributionModels::general::general
  43 (
  44     const dictionary& dict,
  45     cachedRandom& rndGen
  46 )
  47 :
  48     distributionModel(typeName, dict, rndGen),
  49     xy_(distributionModelDict_.lookup("distribution")),
  50     nEntries_(xy_.size()),
  51     minValue_(xy_[0][0]),
  52     maxValue_(xy_[nEntries_-1][0]),
  53     integral_(nEntries_)
  54 {
  55     check();
  56
  57     // normalize the cumulative distributionModel
  58
  59     integral_[0] = 0.0;
  60     for (label i=1; i<nEntries_; i++)
  61     {
  62
  63         scalar k = (xy_[i][1] - xy_[i-1][1])/(xy_[i][0] - xy_[i-1][0]);
  64         scalar d = xy_[i-1][1] - k*xy_[i-1][0];
  65         scalar y1 = xy_[i][0]*(0.5*k*xy_[i][0] + d);
  66         scalar y0 = xy_[i-1][0]*(0.5*k*xy_[i-1][0] + d);
  67         scalar area = y1 - y0;
  68
  69         integral_[i] = area + integral_[i-1];
  70     }
  71
  72     for (label i=0; i<nEntries_; i++)
  73     {
  74         xy_[i][1] /= integral_[nEntries_-1];
  75         integral_[i] /= integral_[nEntries_-1];
  76     }
  77
  78 }
  79
  80
  81 Foam::distributionModels::general::general(const general& p)
  82 :
  83     distributionModel(p),
  84     xy_(p.xy_),
  85     nEntries_(p.nEntries_),
  86     minValue_(p.minValue_),
  87     maxValue_(p.maxValue_),
  88     integral_(p.integral_)
  89 {}
  90
  91
  92 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
  93
  94 Foam::distributionModels::general::~general()
  95 {}
  96
  97
  98 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  99
 100 Foam::scalar Foam::distributionModels::general::sample() const
 101 {
 102     scalar y = rndGen_.sample01<scalar>();
 103
 104     // find the interval where y is in the table
 105     label n=1;
 106     while (integral_[n] <= y)
 107     {
 108         n++;
 109     }
 110
 111     scalar k = (xy_[n][1] - xy_[n-1][1])/(xy_[n][0] - xy_[n-1][0]);
 112     scalar d = xy_[n-1][1] - k*xy_[n-1][0];
 113
 114     scalar alpha = y + xy_[n-1][0]*(0.5*k*xy_[n-1][0] + d) - integral_[n-1];
 115     scalar x = 0.0;
 116
 117     // if k is small it is a linear equation, otherwise it is of second order
 118     if (mag(k) > SMALL)
 119     {
 120         scalar p = 2.0*d/k;
 121         scalar q = -2.0*alpha/k;
 122         scalar sqrtEr = sqrt(0.25*p*p - q);
 123
 124         scalar x1 = -0.5*p + sqrtEr;
 125         scalar x2 = -0.5*p - sqrtEr;
 126         if ((x1 >= xy_[n-1][0]) && (x1 <= xy_[n][0]))
 127         {
 128             x = x1;
 129         }
 130         else
 131         {
 132             x = x2;
 133         }
 134     }
 135     else
 136     {
 137         x = alpha/d;
 138     }
 139
 140     return x;
 141 }
 142
 143
 144 Foam::scalar Foam::distributionModels::general::minValue() const
 145 {
 146     return minValue_;
 147 }
 148
 149
 150 Foam::scalar Foam::distributionModels::general::maxValue() const
 151 {
 152     return maxValue_;
 153 }
 154
 155
 156 // ************************************************************************* //