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 #include <cmath>
43 #include <stdio.h>
44 #include <stdlib.h>
45 #include <string.h>
46 #include "fix_template_sphere.h"
47 #include "fix_template_multiplespheres.h"
48 #include "fix_particledistribution_discrete.h"
49 #include "atom.h"
50 #include "atom_vec.h"
51 #include "domain.h"
52 #include "region.h"
53 #include "update.h"
54 #include "modify.h"
55 #include "output.h"
56 #include "memory.h"
57 #include "error.h"
58 #include "random_park.h"
59 #include "particleToInsert.h"
60 #include "comm.h"
61
62 using namespace LAMMPS_NS;
63 using namespace FixConst;
64
65 #define LMP_DEBUGMODE_SPHERE false
66
67 /* ---------------------------------------------------------------------- */
68
FixParticledistributionDiscrete(LAMMPS * lmp,int narg,char ** arg)69 FixParticledistributionDiscrete::FixParticledistributionDiscrete(LAMMPS *lmp, int narg, char **arg) :
70 Fix(lmp, narg, arg),
71 fix_template_(NULL)
72 {
73 restart_global = 1;
74
75 mass_based = true;
76
77 if(strstr(arg[2],"numberbased"))
78 mass_based = false;
79
80 // random number generator, same for all procs
81
82 if (narg < 7)
83 error->fix_error(FLERR,this,"not enough arguments");
84 random = new RanPark(lmp, arg[3], true);
85 seed = random->getSeed();
86 ntemplates = atoi(arg[4]);
87 if(ntemplates < 1)
88 error->fix_error(FLERR,this,"illegal number of templates");
89
90 templates = new FixTemplateSphere*[ntemplates];
91 distweight = new double[ntemplates];
92 cumweight = new double[ntemplates];
93 parttogen = new int[ntemplates];
94 distorder = new int[ntemplates];
95
96 iarg = 5;
97
98 int itemp=0;
99
100 if(narg != iarg+2*ntemplates)
101 error->fix_error(FLERR,this,"# of templates does not match # of arguments");
102
103 // parse further args
104 do {
105 if(itemp == ntemplates) break;
106 if(narg < iarg+1)
107 error->fix_error(FLERR,this,"not enough arguments");
108 int ifix = modify->find_fix(arg[iarg]);
109
110 if(ifix < 0)
111 error->fix_error(FLERR,this,"invalid ID for fix particletemplate provided");
112
113 if(strncmp(modify->fix[ifix]->style,"particletemplate/",16))
114 error->fix_error(FLERR,this,"fix is not of type particletemplate");
115
116 templates[itemp] = static_cast<FixTemplateSphere*>(modify->fix[ifix]);
117 distweight[itemp] = atof(arg[iarg+1]);
118 if (distweight[itemp] < 0) error->fix_error(FLERR,this,"invalid weight");
119 itemp++;
120 iarg += 2;
121 } while (iarg < narg);
122
123 // check for double use of template which is not allowed
124 for(int i = 0; i < ntemplates; i++)
125 for(int j = 0; j < i; j++)
126 if(templates[i] == templates[j])
127 error->fix_error(FLERR,this,"cannot use the same template twice");
128
129 // normalize distribution
130 double weightsum = 0;
131 for(int i = 0; i < ntemplates; i++)
132 weightsum += distweight[i];
133
134 if(comm->me == 0 && fabs(weightsum-1.) > 0.00001)
135 error->warning(FLERR,"particledistribution/discrete: sum of distribution weights != 1, normalizing distribution");
136
137 for(int i = 0; i < ntemplates; i++)
138 distweight[i]/=weightsum;
139
140 if(mass_based && comm->me == 0 && screen)
141 {
142 fprintf(screen,"Fix particledistribution/discrete (id %s): distribution based on mass%%:\n",this->id);
143 for(int i = 0; i < ntemplates; i++)
144 fprintf(screen," %s: d=%e (max. bounding sphere) mass%%=%f%%\n",templates[i]->id,2.*templates[i]->max_r_bound(),100.*distweight[i]);
145 }
146
147 // convert distribution from mass% to number%
148 // do not do if already number-based
149 if(mass_based)
150 {
151 for(int i=0;i<ntemplates; i++)
152 distweight[i]=distweight[i]/templates[i]->massexpect();
153
154 weightsum=0;
155 for(int i=0;i<ntemplates; i++)
156 weightsum+=distweight[i];
157
158 for(int i=0;i<ntemplates; i++)
159 distweight[i]/=weightsum;
160 }
161
162 if(comm->me == 0 && screen)
163 {
164 fprintf(screen,"Fix particledistribution/discrete (id %s): distribution based on number%%:\n",this->id);
165 for(int i = 0; i < ntemplates; i++)
166 fprintf(screen," %s: d=%e (max. bounding sphere) number%%=%f%%\n",templates[i]->id,2.*templates[i]->max_r_bound(),100.*distweight[i]);
167 }
168
169 cumweight[0] = distweight[0];
170 for(int i = 1; i < ntemplates; i++)
171 cumweight[i] = distweight[i]+cumweight[i-1];
172
173 volexpect=0.;massexpect=0.;
174
175 for(int i = 0; i < ntemplates; i++)
176 {
177 volexpect += templates[i]->volexpect() * distweight[i];
178 massexpect += templates[i]->massexpect() * distweight[i];
179 }
180
181 //get min/maxtype
182 maxtype = 0;
183 mintype = 10000;
184 for(int i = 0; i < ntemplates; i++)
185 {
186 if(templates[i]->maxtype() > maxtype)
187 maxtype = templates[i]->maxtype();
188 if(templates[i]->mintype() < mintype)
189 mintype = templates[i]->mintype();
190 }
191
192 // check which template has the most spheres
193 maxnspheres = 0;
194 for(int i = 0; i < ntemplates;i++)
195 if(templates[i]->number_spheres() > maxnspheres)
196 maxnspheres=templates[i]->number_spheres();
197
198 // sort the distributions by insertion volume (in descending order)
199 // use bubble sort
200 for(int i = 0; i < ntemplates; i++)
201 distorder[i]=i;
202
203 bool swaped;
204 int n = ntemplates;
205 do
206 {
207 swaped = false;
208 for(int i = 0; i < ntemplates-1; i++)
209 {
210 if(templates[distorder[i]]->volexpect() < templates[distorder[i+1]]->volexpect())
211 {
212 //swap
213 int tmp = distorder[i];
214 distorder[i] = distorder[i+1];
215 distorder[i+1] = tmp;
216 swaped = true;
217 }
218 }
219 n--;
220 } while(swaped && n > 0);
221
222 pti = templates[distorder[0]]->pti;
223
224 pti_list = NULL;
225 n_pti = n_pti_max = 0;
226
227 //calc max radius and bounding sphere radius
228
229 maxrad = maxrbound = 0.;
230 minrad = 1000.;
231
232 for(int i = 0; i < ntemplates;i++)
233 if(templates[i]->max_r_bound() > maxrbound)
234 maxrbound = templates[i]->max_r_bound();
235
236 for(int i = 0; i < ntemplates;i++)
237 if(templates[i]->max_rad() > maxrad)
238 maxrad = templates[i]->max_rad();
239
240 for(int i = 0; i < ntemplates;i++)
241 if(templates[i]->min_rad() < minrad)
242 minrad = templates[i]->min_rad();
243
244 }
245
246 /* ---------------------------------------------------------------------- */
247
~FixParticledistributionDiscrete()248 FixParticledistributionDiscrete::~FixParticledistributionDiscrete()
249 {
250 delete []templates;
251 delete []distweight;
252 delete []cumweight;
253 delete []parttogen;
254 delete []distorder;
255 if(pti_list) delete []pti_list;
256 delete random;
257 }
258
259 /* ----------------------------------------------------------------------*/
260
setmask()261 int FixParticledistributionDiscrete::setmask()
262 {
263 int mask = 0;
264 return mask;
265 }
266
267 /* ----------------------------------------------------------------------
268 prepares the fix for a series of randomize_single() commands
269 typically called once per insertion step
270 ------------------------------------------------------------------------- */
271
random_init_single(int ntotal)272 int FixParticledistributionDiscrete::random_init_single(int ntotal)
273 {
274 ninsert = ntotal;
275 ninserted = 0;
276
277 for(int i = 0; i < ntemplates; i++)
278 parttogen[i] = static_cast<int>(static_cast<double>(ninsert) * distweight[i] + random->uniform());
279
280 ninsert = 0;
281 for(int i = 0; i < ntemplates; i++)
282 ninsert += parttogen[i];
283 return ninsert;
284 }
285
286 /* ----------------------------------------------------------------------
287 request one template to generate one pti
288 ------------------------------------------------------------------------- */
289
randomize_single()290 Region* FixParticledistributionDiscrete::randomize_single()
291 {
292 if(ntemplates == 1){
293 templates[0]->randomize_single();
294
295 return templates[0]->region();
296 }
297
298 //choose a template from the discrete distribution, beginning from large to small particles
299 int chosen = 0;
300 int chosendist = distorder[chosen];
301 int ntoinsert = parttogen[chosendist];
302 while(ninserted >= ntoinsert && chosen < ntemplates-1)
303 {
304 chosen++;
305 chosendist = distorder[chosen];
306 ntoinsert += parttogen[chosendist];
307 }
308
309 templates[chosendist]->randomize_single();
310
311 pti = templates[chosendist]->pti;
312
313 ninserted++;
314
315 return templates[chosendist]->region();
316
317 }
318
319 /* ----------------------------------------------------------------------
320 prepares the fix for a series of randomize_list() command
321 also prepares templates
322 - deletes their old lists if present and allocates new lists
323 typically only called once before first insertion step
324
325 allocates for max # particles
326
327 can be called by multiple fix insert commands, so check first if max #
328 particles to be inserted is exceeded and only re-allocate in this case
329 ------------------------------------------------------------------------- */
330
random_init_list(int ntotal)331 void FixParticledistributionDiscrete::random_init_list(int ntotal)
332 {
333 int parttogen_max_i, n_pti_max_requested;
334 int nprocs = comm->nprocs;
335
336 ntotal += 2 * ntemplates;
337
338 // number of requested pti
339 n_pti_max_requested = 0;
340
341 for(int i = 0; i < ntemplates; i++)
342 {
343 parttogen_max_i = static_cast<int>(static_cast<double>(ntotal) * distweight[i] + static_cast<double>(1.01)*(ntemplates+nprocs));
344 n_pti_max_requested += parttogen_max_i;
345
346 // re-allocated if need more ptis in this template than allocated so far
347 if(parttogen_max_i > templates[i]->n_pti_max)
348 {
349 templates[i]->delete_ptilist();
350 templates[i]->init_ptilist(parttogen_max_i);
351 }
352 }
353
354 // re-allocate if need more total ptis in distribution than allocated so far
355 if(n_pti_max_requested > n_pti_max)
356 {
357 n_pti_max = n_pti_max_requested;
358 if(pti_list) delete []pti_list;
359 pti_list = new ParticleToInsert*[n_pti_max];
360
361 }
362
363 }
364
direct_init_list(const int * const parttogen,FixPropertyAtom * const fix_release)365 void FixParticledistributionDiscrete::direct_init_list(const int * const parttogen, FixPropertyAtom * const fix_release)
366 {
367 int n_pti_max_requested = 0;
368 for (int i = 0; i < ntemplates; i++)
369 {
370 if (parttogen[i] > templates[i]->n_pti_max)
371 {
372 templates[i]->delete_ptilist();
373 templates[i]->init_ptilist(parttogen[i], true, fix_release);
374 }
375 n_pti_max_requested += parttogen[i];
376 }
377
378 // re-allocate if need more total ptis in distribution than allocated so far
379 if(n_pti_max_requested > n_pti_max)
380 {
381 n_pti_max = n_pti_max_requested;
382 if(pti_list) delete []pti_list;
383 pti_list = new ParticleToInsert*[n_pti_max];
384
385 }
386 }
387
388 /* ----------------------------------------------------------------------
389 tell all templates to generate their pti_list, wire their pti_list to
390 the list in this fix. returns number of particles to be inserted.
391 typically called once per insertion step
392
393 for exact_number = 1, truncate distribution so to exactly meet
394 requested # particles
395 for exact_number = 0, use random gen to fulfill distribution
396 ------------------------------------------------------------------------- */
397
randomize_list(int ntotal,int insert_groupbit,int exact_number)398 int FixParticledistributionDiscrete::randomize_list(int ntotal,int insert_groupbit,int exact_number)
399 {
400 if(ntotal > n_pti_max)
401 {
402
403 error->one(FLERR,"Faulty implementation: FixParticledistributionDiscrete::randomize_list() called for more particles than defined in random_init_list()");
404 }
405
406 ninsert = ntotal;
407 ninserted = 0;
408
409 // use random generator so long-time average of insertion will represent distribution correctly
410 if(exact_number == 0)
411 {
412
413 for(int i = 0; i < ntemplates; i++)
414 {
415 parttogen[i] = static_cast<int>(static_cast<double>(ninsert) * distweight[i] + random->uniform());
416
417 }
418 }
419 // truncate distribution so # particles to insert is met exactly
420 else
421 {
422 int ninsert_truncated = 0, j;
423 double *remainder = new double[ntemplates], rsum, r;
424
425 // distribute particles and calculate remainder
426 for(int i = 0; i < ntemplates; i++)
427 {
428 parttogen[i] = static_cast<int>(static_cast<double>(ninsert) * distweight[i]);
429 ninsert_truncated += parttogen[i];
430 remainder[i] = static_cast<double>(ninsert) * distweight[i] - static_cast<double>(parttogen[i]);
431
432 }
433
434 int ninsert_gap = ninsert - ninsert_truncated;
435
436 // distribute remaining ninsert_gap particles
437 for(int i = 0; i < ninsert_gap; i++)
438 {
439 r = random->uniform() * static_cast<double>(ninsert_gap);
440 j = 0;
441 rsum = remainder[0];
442
443 while(rsum < r && j < (ntemplates-1) )
444 {
445 j++;
446 rsum += remainder[j];
447 }
448
449 parttogen[j]++;
450 }
451
452 delete []remainder;
453 }
454
455 // count total particle number to be inserted, let templates generate a pti_list
456 ninsert = 0;
457 for(int i = 0; i < ntemplates; i++)
458 {
459 ninsert += parttogen[i];
460 templates[i]->randomize_ptilist(parttogen[i],groupbit | insert_groupbit,distorder[i]);
461 }
462
463 // wire lists, make sure in correct order (large to small particles)
464
465 n_pti = 0;
466 for(int i = 0; i < ntemplates; i++)
467 {
468 int chosendist = distorder[i];
469 for (int j = 0; j < parttogen[chosendist]; j++)
470 {
471 pti_list[n_pti + j] = templates[chosendist]->pti_list[j];
472 }
473 n_pti += parttogen[chosendist];
474 }
475
476 if(n_pti != ninsert)
477 error->one(FLERR,"Internal error in FixParticledistributionDiscrete::randomize_list");
478
479 ninserted = ninsert;
480 return ninsert;
481 }
482
483 /* ---------------------------------------------------------------------- */
484
direct_set_ptlist(const int itemplate,const int i,const void * const data,const int distribution_groupbit)485 void FixParticledistributionDiscrete::direct_set_ptlist(const int itemplate, const int i, const void * const data, const int distribution_groupbit)
486 {
487 templates[itemplate]->direct_set_ptlist(i, data, distribution_groupbit | groupbit, distorder[itemplate]);
488 }
489
490 /* ---------------------------------------------------------------------- */
491
update_ptlist_pointer(const int * ext_parttogen)492 int FixParticledistributionDiscrete::update_ptlist_pointer(const int * ext_parttogen)
493 {
494 n_pti = 0;
495 ninsert = 0;
496 for (int i = 0; i < ntemplates; i++)
497 {
498 ninsert += ext_parttogen[i];
499 int chosendist = distorder[i];
500 parttogen[chosendist] = ext_parttogen[chosendist];
501 for (int j = 0; j < parttogen[chosendist]; j++)
502 {
503 pti_list[n_pti + j] = templates[chosendist]->pti_list[j];
504 }
505 n_pti += parttogen[chosendist];
506 }
507
508 if(n_pti != ninsert)
509 error->one(FLERR,"Internal error in FixParticledistributionDiscrete::update_ptlist_ptr");
510
511 ninserted = ninsert;
512 return ninsert;
513 }
514
515 /* ----------------------------------------------------------------------
516 preparations before insertion
517 ------------------------------------------------------------------------- */
518
pre_insert(int n,class FixPropertyAtom * fp,double val)519 void FixParticledistributionDiscrete::pre_insert(int n,class FixPropertyAtom *fp,double val)
520 {
521 // allow fixes to e.g. update some pointers before set_arrays is called
522 // set_arrays called in ParticleToInsert::insert()
523
524 int nfix = modify->nfix;
525 Fix **fix = modify->fix;
526
527 for (int j = 0; j < nfix; j++)
528 if (fix[j]->create_attribute) fix[j]->pre_set_arrays();
529
530 // set fix property as desired by fix insert
531 // loop to n, not n_pti
532 if(fp)
533 {
534
535 for(int i = 0; i < ntemplates; i++)
536 {
537 if( dynamic_cast<FixTemplateMultiplespheres*>(templates[i]) &&
538 dynamic_cast<FixTemplateMultiplespheres*>(templates[i])->is_bonded()
539 )
540 error->one(FLERR,"'bonded' and setting values for a fix property upon insertion can not be used together");
541 }
542
543 for(int i = 0; i < n; i++)
544 {
545 if (!pti_list[i]->fix_property)
546 {
547 pti_list[i]->fix_property = new FixPropertyAtom*[1];
548 if (pti_list[i]->fix_property_value)
549 error->one(FLERR, "Internal error (fix property pti list)");
550 pti_list[i]->fix_property_value = new double*[1];
551 pti_list[i]->fix_property_value[0] = new double[1];
552 if (pti_list[i]->fix_property_nentry)
553 error->one(FLERR, "Internal error (fix property pti list)");
554 pti_list[i]->fix_property_nentry = new int[1];
555 }
556 pti_list[i]->fix_property[0] = fp;
557 pti_list[i]->fix_property_value[0][0] = val;
558 pti_list[i]->n_fix_property = 1;
559 pti_list[i]->fix_property_nentry[0] = 1;
560 }
561 }
562
563 for(int i = 0; i < n; i++)
564 pti_list[i]->setFixTemplate(fix_template_);
565 }
566
567 /* ----------------------------------------------------------------------
568 set particle properties - only pti needs to know which properties to set
569 loop to n, not n_pti, since not all particles may have been inserted
570 ------------------------------------------------------------------------- */
571
insert(int n)572 int FixParticledistributionDiscrete::insert(int n)
573 {
574 int ninserted_spheres_local = 0;
575 for(int i = 0; i < n; i++)
576 {
577
578 ninserted_spheres_local += pti_list[i]->insert();
579 }
580 return ninserted_spheres_local;
581 }
582
583 /* ----------------------------------------------------------------------
584 wrap up insertion
585 ------------------------------------------------------------------------- */
586
finalize_insertion()587 void FixParticledistributionDiscrete::finalize_insertion()
588 {
589 for(int i = 0; i < ntemplates; i++)
590 templates[i]->finalize_insertion();
591 }
592
593 /* ----------------------------------------------------------------------*/
594
vol_expect()595 double FixParticledistributionDiscrete::vol_expect()
596 {
597 return volexpect;
598 }
599
600 /* ----------------------------------------------------------------------*/
601
mass_expect()602 double FixParticledistributionDiscrete::mass_expect()
603 {
604 return massexpect;
605 }
606
607 /* ----------------------------------------------------------------------*/
608
min_rad(int type)609 double FixParticledistributionDiscrete::min_rad(int type)
610 {
611 //get minrad
612 double minrad_type = 1000.;
613 for(int i = 0; i < ntemplates;i++)
614 {
615 if(
616 (type >= templates[i]->mintype() && type <= templates[i]->maxtype()) &&
617 (templates[i]->min_rad() < minrad_type)
618 )
619 minrad_type = templates[i]->min_rad();
620 }
621
622 return minrad_type;
623 }
624
625 /* ----------------------------------------------------------------------*/
626
max_rad(int type)627 double FixParticledistributionDiscrete::max_rad(int type)
628 {
629 //get maxrad
630 double maxrad_type = 0.;
631 for(int i = 0; i < ntemplates;i++)
632 {
633
634 if(!templates[i]->use_rad_for_cut_neigh_and_ghost())
635 continue;
636
637 if(
638 (type >= templates[i]->mintype() && type <= templates[i]->maxtype()) &&
639 (templates[i]->max_rad() > maxrad_type)
640 )
641 maxrad_type = templates[i]->max_rad();
642 }
643
644 return maxrad_type;
645 }
646
647 /* ----------------------------------------------------------------------*/
648
max_type()649 int FixParticledistributionDiscrete::max_type()
650 {
651 return maxtype;
652 }
653
654 /* ----------------------------------------------------------------------*/
655
min_type()656 int FixParticledistributionDiscrete::min_type()
657 {
658 return mintype;
659 }
660
661 /* ----------------------------------------------------------------------*/
662
max_nspheres()663 int FixParticledistributionDiscrete::max_nspheres()
664 {
665 return maxnspheres;
666 }
667
668 /* ----------------------------------------------------------------------
669 pack entire state of Fix into one write
670 ------------------------------------------------------------------------- */
671
write_restart(FILE * fp)672 void FixParticledistributionDiscrete::write_restart(FILE *fp)
673 {
674 int n = 0;
675 double list[1];
676 list[n++] = static_cast<int>(random->state());
677
678 if (comm->me == 0) {
679 int size = n * sizeof(double);
680 fwrite(&size,sizeof(int),1,fp);
681 fwrite(list,sizeof(double),n,fp);
682 }
683 }
684
685 /* ----------------------------------------------------------------------
686 use state info from restart file to restart the Fix
687 ------------------------------------------------------------------------- */
688
restart(char * buf)689 void FixParticledistributionDiscrete::restart(char *buf)
690 {
691 int n = 0;
692 double *list = (double *) buf;
693
694 seed = static_cast<int> (list[n++]) + comm->me;
695
696 random->reset(seed);
697 }
698
699 /* ----------------------------------------------------------------------
700 generate a hash for unique identification
701 ------------------------------------------------------------------------- */
generate_hash()702 unsigned int FixParticledistributionDiscrete::generate_hash()
703 {
704 unsigned int hash = 0;
705 unsigned int start = seed*420001; // it's magic
706 add_hash_value(ntemplates, start, hash);
707 for (int i=0; i<ntemplates; i++)
708 {
709 add_hash_value(distweight[i], start, hash);
710 add_hash_value(distorder[i], start, hash);
711 add_hash_value((int)templates[i]->generate_hash(), start, hash);
712 }
713 add_hash_value(maxtype, start, hash);
714 add_hash_value(mintype, start, hash);
715 add_hash_value(volexpect, start, hash);
716 add_hash_value(massexpect, start, hash);
717 add_hash_value(maxnspheres, start, hash);
718 return hash;
719 }
720
add_hash_value(const int value,unsigned int & start,unsigned int & hash)721 void FixParticledistributionDiscrete::add_hash_value(const int value, unsigned int &start, unsigned int &hash)
722 {
723 if (value >= 0)
724 hash = hash*start + (unsigned int)value;
725 else
726 hash = hash*start + (unsigned int)(-value) + INT_MAX;
727 start = start*seed;
728 }
729
add_hash_value(double value,unsigned int & start,unsigned int & hash)730 void FixParticledistributionDiscrete::add_hash_value(double value, unsigned int &start, unsigned int &hash)
731 {
732 int ivalue = 0;
733 if (value < 0)
734 value *= -1.0;
735 if (value > 1e-50)
736 {
737 while (value > 1e6)
738 value *= 1e-6;
739
740 while (value < 1e0)
741 value *= 1e6;
742 int high = (int) value;
743 double remainder = (value - (double)high)*1e6;
744 int low = (int) remainder;
745 ivalue = high + low;
746 }
747 add_hash_value(ivalue, start, hash);
748 }
749