1 /* ---------------------------------------------------------------------- 2 This is the 3 4 ██╗ ██╗ ██████╗ ██████╗ ██████╗ ██╗ ██╗████████╗███████╗ 5 ██║ ██║██╔════╝ ██╔════╝ ██╔════╝ ██║ ██║╚══██╔══╝██╔════╝ 6 ██║ ██║██║ ███╗██║ ███╗██║ ███╗███████║ ██║ ███████╗ 7 ██║ ██║██║ ██║██║ ██║██║ ██║██╔══██║ ██║ ╚════██║ 8 ███████╗██║╚██████╔╝╚██████╔╝╚██████╔╝██║ ██║ ██║ ███████║ 9 ╚══════╝╚═╝ ╚═════╝ ╚═════╝ ╚═════╝ ╚═╝ ╚═╝ ╚═╝ ╚══════╝® 10 11 DEM simulation engine, released by 12 DCS Computing Gmbh, Linz, Austria 13 http://www.dcs-computing.com, office@dcs-computing.com 14 15 LIGGGHTS® is part of CFDEM®project: 16 http://www.liggghts.com | http://www.cfdem.com 17 18 Core developer and main author: 19 Christoph Kloss, christoph.kloss@dcs-computing.com 20 21 LIGGGHTS® is open-source, distributed under the terms of the GNU Public 22 License, version 2 or later. It is distributed in the hope that it will 23 be useful, but WITHOUT ANY WARRANTY; without even the implied warranty 24 of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. You should have 25 received a copy of the GNU General Public License along with LIGGGHTS®. 26 If not, see http://www.gnu.org/licenses . See also top-level README 27 and LICENSE files. 28 29 LIGGGHTS® and CFDEM® are registered trade marks of DCS Computing GmbH, 30 the producer of the LIGGGHTS® software and the CFDEM®coupling software 31 See http://www.cfdem.com/terms-trademark-policy for details. 32 33 ------------------------------------------------------------------------- 34 Contributing author and copyright for this file: 35 (if not contributing author is listed, this file has been contributed 36 by the core developer) 37 38 Copyright 2012- DCS Computing GmbH, Linz 39 Copyright 2009-2012 JKU Linz 40 ------------------------------------------------------------------------- */ 41 42 #ifndef LMP_PROBABILITY_DISTRIBUTION_H 43 #define LMP_PROBABILITY_DISTRIBUTION_H 44 45 #include <cmath> 46 #include <stdio.h> 47 #include <string.h> 48 #include "random_park.h" 49 #include "error.h" 50 #include "pointers.h" 51 52 enum{RANDOM_CONSTANT,RANDOM_UNIFORM,RANDOM_GAUSSIAN,RANDOM_LOGNORMAL}; 53 54 namespace LMP_PROBABILITY_NS { 55 class PDF 56 { 57 public: PDF(LAMMPS_NS::Error * error)58 PDF(LAMMPS_NS::Error *error) 59 { 60 mu_ = sigma_ = min_ = max_ = 0.; 61 h1_ = h2_ = 0.; 62 this->error = error; 63 } ~PDF()64 ~PDF(){} 65 66 int rand_style_; 67 68 double mu_,sigma_; 69 double min_,max_; 70 71 // helper 72 double h1_,h2_; 73 74 LAMMPS_NS::Error *error; 75 rand_style()76 inline int rand_style() 77 { return rand_style_; } 78 set_min_max(double min,double max)79 inline void set_min_max(double min,double max) 80 { 81 min_ = min; 82 max_ = max; 83 } 84 set_params(double)85 template<int RAND_STYLE> void set_params(double) 86 { error->all(FLERR,"Faulty usage of Probability::set_params"); } 87 set_params(double,double)88 template<int RAND_STYLE> void set_params(double,double) 89 { error->all(FLERR,"Faulty usage of Probability::set_params"); } 90 }; 91 pdf_max(PDF * pdf)92 inline double pdf_max(PDF *pdf) 93 { 94 return pdf->max_; 95 } 96 pdf_min(PDF * pdf)97 inline double pdf_min(PDF *pdf) 98 { 99 return pdf->min_; 100 } 101 expectancy_value(PDF * pdf)102 template <int RAND_STYLE> inline double expectancy_value(PDF *pdf) 103 { 104 pdf->error->all(FLERR,"Faulty usage of Probability::expectancy"); 105 return 0.; 106 } 107 cubic_expectancy_value(PDF * pdf)108 template <int RAND_STYLE> inline double cubic_expectancy_value(PDF *pdf) 109 { 110 pdf->error->all(FLERR,"Faulty usage of Probability::volume_expectancy"); 111 return 0.; 112 } 113 rand_value(PDF * pdf,LAMMPS_NS::RanPark * rp)114 template <int RAND_STYLE> inline double rand_value(PDF *pdf,LAMMPS_NS::RanPark *rp) 115 { 116 pdf->error->all(FLERR,"Faulty usage of Probability::rand"); 117 return 0.; 118 } 119 120 //------------------------------------------------------------------------------ 121 // CONSTANT 122 //------------------------------------------------------------------------------ 123 124 template<> inline void PDF::set_params<RANDOM_CONSTANT>(double val) 125 { 126 rand_style_ = RANDOM_CONSTANT; 127 mu_ = val; 128 set_min_max(mu_,mu_); 129 } 130 131 template<> inline double cubic_expectancy_value<RANDOM_CONSTANT>(PDF *pdf) 132 { 133 134 return pdf->mu_*pdf->mu_*pdf->mu_; 135 } 136 137 template<> inline double expectancy_value<RANDOM_CONSTANT>(PDF *pdf) 138 { 139 140 return pdf->mu_; 141 } 142 143 template<> inline double rand_value<RANDOM_CONSTANT>(PDF *pdf,LAMMPS_NS::RanPark *rp) 144 { 145 return pdf->mu_; 146 } 147 148 //------------------------------------------------------------------------------ 149 // UNIFORM 150 //------------------------------------------------------------------------------ 151 152 template<> inline void PDF::set_params<RANDOM_UNIFORM>(double min, double max) 153 { 154 rand_style_ = RANDOM_UNIFORM; 155 set_min_max(min,max); 156 h1_ = 2./(1./(min_*min_)-1./(max_*max_)); 157 h2_ = h1_/(2.*min_*min_); 158 } 159 160 template<> inline double cubic_expectancy_value<RANDOM_UNIFORM>(PDF *pdf) 161 { 162 return 0.25*(pdf->max_*pdf->max_*pdf->max_+ 163 pdf->max_*pdf->max_*pdf->min_+ 164 pdf->max_*pdf->min_*pdf->min_+ 165 pdf->min_*pdf->min_*pdf->min_); 166 } 167 168 template<> inline double expectancy_value<RANDOM_UNIFORM>(PDF *pdf) 169 { 170 return sqrt(pdf->h1_/(2.*(pdf->h2_-0.5))); 171 } 172 173 template<> inline double rand_value<RANDOM_UNIFORM>(PDF *pdf,LAMMPS_NS::RanPark *rp) 174 { 175 double rn = rp->uniform(); 176 return sqrt(pdf->h1_/(2.*(pdf->h2_-rn))); 177 } 178 179 //------------------------------------------------------------------------------ 180 // GAUSSIAN 181 //------------------------------------------------------------------------------ 182 183 template<> inline void PDF::set_params<RANDOM_GAUSSIAN>(double mu, double sigma) 184 { 185 rand_style_ = RANDOM_GAUSSIAN; 186 mu_ = mu; 187 sigma_ = sigma; 188 189 // set min-max to +- 3 sigma (99.73% of all values) 190 set_min_max(mu_-3.*sigma_, mu_+3.*sigma_); 191 192 if(min_ < 0.) 193 error->all(FLERR,"Probablity distribution: mu-3*sigma < 0, please increase mu or decrease sigma"); 194 } 195 196 template<> inline double cubic_expectancy_value<RANDOM_GAUSSIAN>(PDF *pdf) 197 { 198 return pdf->mu_*(pdf->mu_*pdf->mu_+3*pdf->sigma_*pdf->sigma_); 199 } 200 201 template<> inline double expectancy_value<RANDOM_GAUSSIAN>(PDF *pdf) 202 { 203 return pdf->mu_; 204 } 205 206 template<> inline double rand_value<RANDOM_GAUSSIAN>(PDF *pdf,LAMMPS_NS::RanPark *rp) 207 { 208 double value; 209 do 210 { 211 value = pdf->mu_ + rp->gaussian() * pdf->sigma_; 212 } while (value < pdf->min_ || value > pdf->max_); 213 return value; 214 } 215 216 //------------------------------------------------------------------------------ 217 // LOGNORMAL 218 //------------------------------------------------------------------------------ 219 220 template<> inline void PDF::set_params<RANDOM_LOGNORMAL>(double mu, double sigma) 221 { 222 error->all(FLERR,"lognormal distribution currently deactivated"); 223 224 rand_style_ = RANDOM_LOGNORMAL; 225 mu_ = mu; 226 sigma_ = sigma; 227 228 // also here, take +- 3 sigma as min/max 229 // change in expectancy considered negligable 230 double min = exp(mu_ - 3. * sigma_); 231 double max = exp(mu_ + 3. * sigma_); 232 set_min_max(min, max); 233 if(min_ < 0.) 234 error->all(FLERR,"Probablity distribution: exp(mu-3*sigma) < 0, please increase mu or decrease sigma"); 235 } 236 237 template<> inline double cubic_expectancy_value<RANDOM_LOGNORMAL>(PDF *pdf) 238 { 239 return exp(3.*pdf->mu_+4.5*pdf->sigma_*pdf->sigma_); 240 } 241 242 template<> inline double expectancy_value<RANDOM_LOGNORMAL>(PDF *pdf) 243 { 244 return exp(pdf->mu_ + 0.5 * pdf->sigma_ * pdf->sigma_); 245 } 246 247 template<> inline double rand_value<RANDOM_LOGNORMAL>(PDF *pdf,LAMMPS_NS::RanPark *rp) 248 { 249 double value; 250 do 251 { 252 value = exp(pdf->mu_ + rp->gaussian() * pdf->sigma_); 253 } while (value < pdf->min_ || value > pdf->max_); 254 return value; 255 } 256 257 //------------------------------------------------------------------------------ 258 // MASTER FUNCTIONS 259 //------------------------------------------------------------------------------ 260 expectancy(PDF * pdf)261 inline double expectancy(PDF *pdf) 262 { 263 if(pdf->rand_style_ == RANDOM_CONSTANT) return expectancy_value<RANDOM_CONSTANT>(pdf); 264 else if(pdf->rand_style_ == RANDOM_UNIFORM) return expectancy_value<RANDOM_UNIFORM>(pdf); 265 else if(pdf->rand_style_ == RANDOM_GAUSSIAN) return expectancy_value<RANDOM_GAUSSIAN>(pdf); 266 else if(pdf->rand_style_ == RANDOM_LOGNORMAL) return expectancy_value<RANDOM_LOGNORMAL>(pdf); 267 else pdf->error->all(FLERR,"Faulty implemantation in Probability::expectancy"); 268 return 0.; 269 } 270 cubic_expectancy(PDF * pdf)271 inline double cubic_expectancy(PDF *pdf) 272 { 273 if(pdf->rand_style_ == RANDOM_CONSTANT) return cubic_expectancy_value<RANDOM_CONSTANT>(pdf); 274 else if(pdf->rand_style_ == RANDOM_UNIFORM) return cubic_expectancy_value<RANDOM_UNIFORM>(pdf); 275 else if(pdf->rand_style_ == RANDOM_GAUSSIAN) return cubic_expectancy_value<RANDOM_GAUSSIAN>(pdf); 276 else if(pdf->rand_style_ == RANDOM_LOGNORMAL) return cubic_expectancy_value<RANDOM_LOGNORMAL>(pdf); 277 else pdf->error->all(FLERR,"Faulty implemantation in Probability::expectancy"); 278 return 0.; 279 } 280 rand(PDF * pdf,LAMMPS_NS::RanPark * rp)281 inline double rand(PDF *pdf,LAMMPS_NS::RanPark *rp) 282 { 283 if(pdf->rand_style_ == RANDOM_CONSTANT) return rand_value<RANDOM_CONSTANT>(pdf,rp); 284 else if(pdf->rand_style_ == RANDOM_UNIFORM) return rand_value<RANDOM_UNIFORM>(pdf,rp); 285 else if(pdf->rand_style_ == RANDOM_GAUSSIAN) return rand_value<RANDOM_GAUSSIAN>(pdf,rp); 286 else if(pdf->rand_style_ == RANDOM_LOGNORMAL) return rand_value<RANDOM_LOGNORMAL>(pdf,rp); 287 else pdf->error->all(FLERR,"Faulty implemantation in Probability::rand"); 288 return 0.; 289 } 290 291 }; 292 293 #endif 294