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 Christoph Kloss (DCS Computing GmbH, Linz)
36 Christoph Kloss (JKU Linz)
37 Richard Berger (JKU Linz)
38
39 Copyright 2012- DCS Computing GmbH, Linz
40 Copyright 2009-2015 JKU Linz
41 ------------------------------------------------------------------------- */
42 #include <cmath>
43 #include <algorithm>
44 #include <stdlib.h>
45 #include <string.h>
46 #include "fix_insert_stream.h"
47 #include "fix_mesh_surface.h"
48 #include "atom.h"
49 #include "atom_vec.h"
50 #include "force.h"
51 #include "update.h"
52 #include "comm.h"
53 #include "modify.h"
54 #include "vector_liggghts.h"
55 #include "domain.h"
56 #include "random_park.h"
57 #include "memory.h"
58 #include "error.h"
59 #include "fix_property_atom.h"
60 #include "fix_property_atom_tracer_stream.h"
61 #include "fix_particledistribution_discrete.h"
62 #include "fix_multisphere.h"
63 #include "fix_template_sphere.h"
64 #include "particleToInsert.h"
65 #include "tri_mesh_planar.h"
66
67 enum{FACE_NONE,FACE_MESH,FACE_CIRCLE};
68
69 using namespace LAMMPS_NS;
70 using namespace FixConst;
71
72 #define FIX_INSERT_NTRY_SUBBOX 500
73 #define FIX_INSERT_STREAM_TINY 1e-14
74
75 /* ---------------------------------------------------------------------- */
76
FixInsertStream(LAMMPS * lmp,int narg,char ** arg)77 FixInsertStream::FixInsertStream(LAMMPS *lmp, int narg, char **arg) :
78 FixInsert(lmp, narg, arg),
79 recalc_release_ms(false),
80 dt_ratio(0.),
81 save_template_(false),
82 fix_template_(NULL)
83 {
84 // set defaults first, then parse args
85 init_defaults();
86
87 bool hasargs = true;
88 while(iarg < narg && hasargs)
89 {
90 hasargs = false;
91
92 if (strcmp(arg[iarg],"insertion_face") == 0)
93 {
94
95 if (iarg+2 > narg) error->fix_error(FLERR,this,"not enough arguments");
96 int f_i = modify->find_fix(arg[iarg+1]);
97 if (f_i == -1) error->fix_error(FLERR,this,"Could not find fix mesh/surface id you provided");
98 if (strncmp(modify->fix[f_i]->style,"mesh",4))
99 error->fix_error(FLERR,this,"The fix belonging to the id you provided is not of type mesh");
100 ins_face = (static_cast<FixMeshSurface*>(modify->fix[f_i]))->triMesh();
101 ins_face->useAsInsertionMesh(false);
102 face_style = FACE_MESH;
103 iarg += 2;
104 hasargs = true;
105 } else if (strcmp(arg[iarg],"extrude_length") == 0) {
106 if (iarg+2 > narg) error->fix_error(FLERR,this,"not enough arguments");
107 extrude_length = atof(arg[iarg+1]);
108 if(extrude_length < 0. ) error->fix_error(FLERR,this,"invalid extrude_length");
109 iarg += 2;
110 hasargs = true;
111 } else if (strcmp(arg[iarg],"duration") == 0 || strcmp(arg[iarg],"duration_time") == 0) {
112 if (iarg+2 > narg) error->fix_error(FLERR,this,"not enough arguments");
113 if(strcmp(arg[iarg],"duration_time") == 0)
114 duration = static_cast<int>(atof(arg[iarg+1])/update->dt);
115 else
116 duration = atoi(arg[iarg+1]);
117 if(duration < 1 ) error->fix_error(FLERR,this,"'duration' can not be < 1");
118 iarg += 2;
119 hasargs = true;
120 } else if (strcmp(arg[iarg],"parallel") == 0) {
121 if (iarg+2 > narg) error->fix_error(FLERR,this,"not enough arguments");
122 if(strcmp("yes",arg[iarg+1]) == 0)
123 parallel = true;
124 else if(strcmp("no",arg[iarg+1]) == 0)
125 parallel = false;
126 else error->fix_error(FLERR,this,"expecting 'yes' or 'no' for 'parallel'");
127 iarg += 2;
128 hasargs = true;
129 } else if (strcmp(arg[iarg],"ntry_mc") == 0) {
130 if (iarg+2 > narg) error->fix_error(FLERR,this,"not enough arguments");
131 ntry_mc = atoi(arg[iarg+1]);
132 if(ntry_mc < 1000) error->fix_error(FLERR,this,"ntry_mc must be > 1000");
133 iarg += 2;
134 hasargs = true;
135 }
136 else if (strcmp(arg[iarg], "save_template") == 0)
137 {
138 if (iarg+2 > narg)
139 error->fix_error(FLERR,this,"not enough arguments");
140
141 if(strcmp("yes",arg[iarg+1]) == 0)
142 save_template_ = true;
143 else if(strcmp("no",arg[iarg+1]) == 0)
144 save_template_ = false;
145 else
146 error->fix_error(FLERR,this,"expecting 'yes' or 'no' for 'save_template'");
147 iarg += 2;
148 hasargs = true;
149 }
150 else if (0 == strcmp(style,"insert/stream"))
151 error->fix_error(FLERR,this,"unknown keyword or wrong keyword order");
152 }
153
154 fix_release = NULL;
155 i_am_integrator = false;
156
157 tracer = NULL;
158 ntracer = 0;
159
160 ins_fraction = 0.;
161 do_ins_fraction_calc = true;
162
163 nevery = 1;
164 }
165
166 /* ---------------------------------------------------------------------- */
167
~FixInsertStream()168 FixInsertStream::~FixInsertStream()
169 {
170 if(tracer) delete []tracer;
171 }
172
173 /* ---------------------------------------------------------------------- */
174
post_create()175 void FixInsertStream::post_create()
176 {
177 FixInsert::post_create();
178
179 // only register property if I am the first fix/insert/stream in the simulation
180
181 if(modify->n_fixes_style(style) == 1)
182 {
183 const char* fixarg[22];
184 fixarg[0]="release_fix_insert_stream";
185 fixarg[1]="all";
186 fixarg[2]="property/atom";
187 fixarg[3]="release_fix_insert_stream";
188 fixarg[4]="vector";
189 fixarg[5]="yes";
190 fixarg[6]="yes";
191 fixarg[7]="no";
192 fixarg[8]="0.";
193 fixarg[9]="0.";
194 fixarg[10]="0.";
195 fixarg[11]="0.";
196 fixarg[12]="0.";
197 fixarg[13]="0.";
198 fixarg[14]="0.";
199 fixarg[15]="0.";
200 fixarg[16]="0.";
201 fixarg[17]="0.";
202 fixarg[18]="0.";
203 fixarg[19]="0.";
204 fixarg[20]="0.";
205 fixarg[21]="0.";
206 modify->add_fix_property_atom(22,const_cast<char**>(fixarg),style);
207
208 fix_release = static_cast<FixPropertyAtom*>(modify->find_fix_property("release_fix_insert_stream","property/atom","vector",14,0,style));
209 if(!fix_release) error->fix_error(FLERR,this,"Internal error in fix insert/stream");
210
211 if( modify->fix_restart_in_progress())
212 recalc_release_restart();
213 }
214
215 if (save_template_)
216 {
217 fix_template_ = static_cast<FixPropertyAtom*>(modify->find_fix_property("insertion_template_", "property/atom", "scalar", 1, 0, style, false));
218 if (!fix_template_)
219 {
220 const char *fixarg[9];
221 fixarg[0] = "insertion_template_";
222 fixarg[1] = "all";
223 fixarg[2] = "property/atom";
224 fixarg[3] = "insertion_template_";
225 fixarg[4] = "scalar";
226 fixarg[5] = "yes"; // restart
227 fixarg[6] = "yes"; // ghost
228 fixarg[7] = "no"; // reverse
229 fixarg[8] = "-1.0";
230 fix_template_ = modify->add_fix_property_atom(9, const_cast<char**>(fixarg), style);
231 }
232 fix_distribution->save_templates(fix_template_);
233 }
234 }
235
236 /* ---------------------------------------------------------------------- */
237
pre_delete(bool unfixflag)238 void FixInsertStream::pre_delete(bool unfixflag)
239 {
240 // delete if I am the last fix of this style to be deleted
241 if(unfixflag && modify->n_fixes_style(style) == 1)
242 modify->delete_fix("release_fix_insert_stream");
243 }
244
245 /* ---------------------------------------------------------------------- */
246
init_defaults()247 void FixInsertStream::init_defaults()
248 {
249 face_style = FACE_NONE;
250 extrude_length = 0.;
251
252 extrude_length_min = extrude_length_max = 0.;
253
254 duration = 0;
255
256 parallel = false;
257
258 ntry_mc = 100000;
259
260 vel_normal_to_face = false;
261 }
262
263 /* ---------------------------------------------------------------------- */
264
register_tracer_callback(FixPropertyAtomTracerStream * tr)265 void FixInsertStream::register_tracer_callback(FixPropertyAtomTracerStream* tr)
266 {
267 // just return if I already have this callback
268 for(int i = 0; i < ntracer; i++)
269 if(tracer[i] == tr) return;
270
271 FixPropertyAtomTracerStream** tracer_new = new FixPropertyAtomTracerStream*[ntracer+1];
272
273 for(int i = 0; i < ntracer; i++)
274 tracer_new[i] = tracer[i];
275
276 tracer_new[ntracer] = tr;
277 ntracer++;
278 delete []tracer;
279 tracer = tracer_new;
280 }
281
282 /* ----------------------------------------------------------------------
283 calculate ninsert, insert_every, ninsert_per, massinsert, flowrates etc
284 also perform error checks
285 ------------------------------------------------------------------------- */
286
calc_insertion_properties()287 void FixInsertStream::calc_insertion_properties()
288 {
289 double dt,dot,extrude_vec[3],t1[3],t2[3];
290
291 // error check on insertion face
292 if(face_style == FACE_NONE)
293 error->fix_error(FLERR,this,"must define an insertion face");
294
295 // check properties of insertion face
296 if(face_style == FACE_MESH)
297 {
298 // check if face planar
299 if(!ins_face->isPlanar())
300 error->fix_error(FLERR,this,"command requires a planar face for insertion");
301
302 if(all_in_flag)
303 {
304 if(!dynamic_cast<TriMeshPlanar*>(ins_face))
305 error->fix_error(FLERR,this,"using all_in yes requires you to use a fix mesh/surface/planar");
306 }
307
308 // get normal vector of face 0
309 ins_face->surfaceNorm(0,normalvec);
310
311 // flip normal vector so dot product with v_insert is > 0
312 dot = vectorDot3D(v_insert,normalvec);
313 if(dot < 0) vectorScalarMult3D(normalvec,-1.);
314
315 // calc v normal
316 dot = vectorDot3D(v_insert,normalvec);
317 vectorCopy3D(normalvec,v_normal);
318 vectorScalarMult3D(v_normal,dot);
319
320 double diff[3];
321 vectorSubtract3D(v_insert,v_normal,diff);
322
323 if(vectorMag3DSquared(diff) < 1e-6)
324 vel_normal_to_face = true;
325 else
326 vel_normal_to_face = false;
327
328 // error check on v normal
329 if(vectorMag3D(v_normal) < 1.e-3)
330 error->fix_error(FLERR,this,"insertion velocity projected on face normal is < 1e-3");
331
332 // get reference point on face
333 ins_face->node(0,0,p_ref);
334 }
335 else error->fix_error(FLERR,this,"FixInsertStream::calc_insertion_properties(): Implementation missing");
336
337 // error check on insertion velocity
338 if(vectorMag3D(v_insert) < 1e-5)
339 error->fix_error(FLERR,this,"insertion velocity too low");
340
341 // further error-checks
342 if(insert_every == -1 && extrude_length == 0.)
343 error->fix_error(FLERR,this,"must define either 'insert_every' or 'extrude_length'");
344 if(insert_every > -1 && extrude_length > 0.)
345 error->fix_error(FLERR,this,"must not provide both 'insert_every' and 'extrude_length'");
346 if(extrude_length > 0. && duration > 0)
347 error->fix_error(FLERR,this,"must not provide both 'extrude_length' and 'duration'");
348
349 dt = update->dt;
350
351 // if extrude_length given, calculate insert_every
352 if(insert_every == -1)
353 {
354 // no duration allowed here (checked before)
355
356 if(extrude_length < 3.*max_r_bound() && (all_in_flag || check_ol_flag))
357 error->fix_error(FLERR,this,"'extrude_length' is too small");
358 // add TINY for resolving round-off
359 insert_every = static_cast<int>((extrude_length+FIX_INSERT_STREAM_TINY)/(dt*vectorMag3D(v_normal)));
360
361 if(insert_every == 0)
362 error->fix_error(FLERR,this,"insertion velocity too high or extrude_length too low");
363 }
364 // if insert_every given, calculate extrude_length
365 // take into account duration can be != insert_every
366 else
367 {
368 if(insert_every < 1) error->fix_error(FLERR,this,"'insert_every' must be > 0");
369
370 // duration = insert_every by default (if already > 0, defined directly)
371 if(duration == 0) duration = insert_every;
372 else if (duration > insert_every) error->fix_error(FLERR,this,"'duration' > 'insert_every' not allowed");
373
374 extrude_length = static_cast<double>(duration) * dt * vectorMag3D(v_normal);
375
376 if(extrude_length < 2.*max_r_bound())
377 error->fix_error(FLERR,this,"'insert_every' or 'vel' is too small, or (bounding) radius of inserted particles too large");
378 }
379
380 // ninsert - if ninsert not defined directly, calculate it
381 if(ninsert == 0 && ninsert_exists)
382 {
383 if(massinsert/fix_distribution->mass_expect() > 2.e9)
384 error->fix_error(FLERR,this,"you are attempting to insert more than 2e9 particles. Reduce the mass to be inserted or increase the particle diameter");
385
386 if(massinsert > 0.) ninsert = static_cast<int>((massinsert+FIX_INSERT_STREAM_TINY) / fix_distribution->mass_expect());
387 else error->fix_error(FLERR,this,"must define either 'nparticles' or 'mass'");
388 }
389
390 // flow rate
391 if(nflowrate == 0.)
392 {
393 if(massflowrate == 0.) error->fix_error(FLERR,this,"must define either 'massrate' or 'particlerate'");
394 nflowrate = massflowrate / fix_distribution->mass_expect();
395 }
396 else massflowrate = nflowrate * fix_distribution->mass_expect();
397
398 // ninsert_per and massinsert
399 ninsert_per = nflowrate*(static_cast<double>(insert_every)*dt);
400 if(ninsert_exists) massinsert = static_cast<double>(ninsert) * fix_distribution->mass_expect();
401
402 // calculate bounding box of extruded face
403 if(face_style == FACE_MESH)
404 {
405 // get bounding box for face
406 ins_face->getGlobalBoundingBox().getBoxBounds(ins_vol_xmin,ins_vol_xmax);
407
408 // get bounding box for extruded face - store in t1,t2
409 vectorScalarMult3D(normalvec,-extrude_length,extrude_vec);
410 vectorAdd3D(ins_vol_xmin,extrude_vec,t1);
411 vectorAdd3D(ins_vol_xmax,extrude_vec,t2);
412
413 // take min and max
414 vectorComponentMin3D(ins_vol_xmin,t1,ins_vol_xmin);
415 vectorComponentMax3D(ins_vol_xmax,t2,ins_vol_xmax);
416
417 }
418 else error->fix_error(FLERR,this,"Missing implementation in calc_insertion_properties()");
419
420 extrude_length_min = 0.;
421 extrude_length_max = extrude_length;
422 }
423
424 /* ---------------------------------------------------------------------- */
425
setmask()426 int FixInsertStream::setmask()
427 {
428 int mask = FixInsert::setmask();
429 mask |= END_OF_STEP;
430 return mask;
431 }
432
433 /* ---------------------------------------------------------------------- */
434
init()435 void FixInsertStream::init()
436 {
437
438 FixInsert::init();
439
440 if(fix_multisphere && v_randomSetting != RANDOM_CONSTANT)
441 error->fix_error(FLERR,this,"Currently only fix insert/stream with multisphere particles only supports constant velocity");
442
443 fix_release = static_cast<FixPropertyAtom*>(modify->find_fix_property("release_fix_insert_stream","property/atom","vector",5,0,style));
444 if(!fix_release) error->fix_error(FLERR,this,"Internal error if fix insert/stream");
445 fix_release->set_internal();
446
447 i_am_integrator = modify->i_am_first_of_style(this);
448
449 // error check on insertion face
450 if(face_style == FACE_NONE)
451 error->fix_error(FLERR,this,"must define an insertion face");
452
453 if(ins_face->isMoving() || ins_face->isScaling())
454 error->fix_error(FLERR,this,"cannot translate, rotate, scale mesh which is used for particle insertion");
455
456 if(recalc_release_ms)
457 {
458 recalc_release_ms = false;
459
460 if(fix_multisphere && dt_ratio > 0.)
461 fix_multisphere->data().recalc_n_steps(dt_ratio);
462
463 }
464 }
465
466 /* ---------------------------------------------------------------------- */
467
setup_pre_exchange()468 void FixInsertStream::setup_pre_exchange()
469 {
470
471 }
472
473 /* ---------------------------------------------------------------------- */
474
insertion_fraction()475 double FixInsertStream::insertion_fraction()
476 {
477
478 // have to re-calculate insertion fraction for my subbox
479 // in case subdomains of simulation box are changing
480
481 if(domain->box_change || do_ins_fraction_calc || ins_face->isMoving())
482 calc_ins_fraction();
483
484 return ins_fraction;
485 }
486
487 /* ----------------------------------------------------------------------
488 calculate insertion fraction for my subbox
489 has to be called at initialization and before every insertion in case
490 box is changing
491 ------------------------------------------------------------------------- */
492
calc_ins_fraction()493 void FixInsertStream::calc_ins_fraction()
494 {
495
496 do_ins_fraction_calc = false;
497
498 double pos[3], boxedgevec[3], dot;
499 int n_in_local = 0, n_test = ntry_mc;
500
501 for(int i = 0; i < n_test; i++)
502 {
503 generate_random_global(pos);
504
505 if(domain->is_in_subdomain(pos))
506 n_in_local++;
507 }
508
509 ins_fraction = static_cast<double>(n_in_local)/static_cast<double>(n_test);
510
511 // also calculate min and max extrusion
512 // this can speed up insertion if extrusion volume extends across multiple procs
513
514 if(parallel)
515 {
516 extrude_length_min = extrude_length;
517 extrude_length_max = 0.;
518
519 for(int ix = 0; ix < 2; ix++)
520 for(int iy = 0; iy < 2; iy++)
521 for(int iz = 0; iz < 2; iz++)
522 {
523 vectorConstruct3D
524 (
525 boxedgevec,
526 (ix == 0 ? domain->sublo[0] : domain->subhi[0]) - p_ref[0],
527 (iy == 0 ? domain->sublo[1] : domain->subhi[1]) - p_ref[1],
528 (iz == 0 ? domain->sublo[2] : domain->subhi[2]) - p_ref[2]
529 );
530
531 dot = -vectorDot3D(boxedgevec,normalvec);
532
533 if(dot > 0. && dot < extrude_length)
534 {
535 extrude_length_max = std::max(extrude_length_max,dot);
536 extrude_length_min = std::min(extrude_length_min,dot);
537 }
538 else if(dot < 0.)
539 extrude_length_min = 0.;
540 else if(dot >= extrude_length)
541 extrude_length_max = extrude_length;
542 }
543 if(extrude_length_min == extrude_length)
544 extrude_length_min = 0.;
545 if(extrude_length_max == 0.)
546 extrude_length_max = extrude_length;
547 }
548
549 double ins_fraction_all;
550 MPI_Sum_Scalar(ins_fraction,ins_fraction_all,world);
551 if(ins_fraction_all < 0.9 || ins_fraction_all > 1.1)
552 error->fix_error(FLERR,this,"insertion volume could not be distributed properly in parallel. "
553 "Bad decomposition or insertion face extrusion is too small or outside domain");
554 }
555
556 /* ---------------------------------------------------------------------- */
557
pre_insert()558 bool FixInsertStream::pre_insert()
559 {
560 if((!domain->is_in_domain(ins_vol_xmin) || !domain->is_in_domain(ins_vol_xmax)) && comm->me == 0)
561 error->warning(FLERR,"Fix insert/stream: Extruded insertion face extends outside domain, may not insert all particles correctly");
562
563 return true;
564 }
565
566 /* ---------------------------------------------------------------------- */
567
is_nearby(int i)568 inline int FixInsertStream::is_nearby(int i)
569 {
570 double pos_rel[3], pos_projected[3], t[3];
571 double **x = atom->x;
572
573 vectorSubtract3D(x[i],p_ref,pos_rel);
574 double dist_normal = vectorDot3D(pos_rel,normalvec);
575
576 // on wrong side of extrusion
577 if(dist_normal > maxrad) return 0;
578
579 // on right side of extrusion, but too far away
580 // 3*maxrad b/c extrude_length+rad is max extrusion for overlapcheck yes
581 if(dist_normal < -(extrude_length + 3.*maxrad)) return 0;
582
583 // on right side of extrusion, within extrude_length
584 // check if projection is on face or not
585
586 vectorScalarMult3D(normalvec,dist_normal,t);
587 vectorAdd3D(x[i],t,pos_projected);
588
589 //TODO also should check if projected point is NEAR surface
590
591 return ins_face->isOnSurface(pos_projected);
592 }
593
594 /* ---------------------------------------------------------------------- */
595
getBoundingBox()596 BoundingBox FixInsertStream::getBoundingBox() {
597 BoundingBox bb = ins_face->getGlobalBoundingBox();
598
599 const double cut = 3*maxrad;
600 const double delta = -(extrude_length + 2*cut);
601 bb.extrude(delta, normalvec);
602 bb.shrinkToSubbox(domain->sublo, domain->subhi);
603
604 const double extend = 3*maxrad /*cut*/ + 2.*fix_distribution->max_r_bound();
605 bb.extendByDelta(extend);
606
607 return bb;
608 }
609
610 /* ----------------------------------------------------------------------
611 generate random positions on insertion face
612 extrude by random length in negative face normal direction
613 currently only implemented for all_in_flag = 0
614 since would be tedious to check/implement otherwise
615 ------------------------------------------------------------------------- */
616
generate_random(double * pos,double rad)617 inline void FixInsertStream::generate_random(double *pos, double rad)
618 {
619 double r, ext[3];
620
621 // generate random position on the mesh
622
623 if(all_in_flag)
624 ins_face->generateRandomOwnedGhostWithin(pos,rad);
625
626 else
627 ins_face->generateRandomOwnedGhost(pos);
628
629 // extrude the position
630
631 if(check_ol_flag)
632 r = -1.*(random->uniform()*(extrude_length_max ) + rad + extrude_length_min);
633 else
634 r = -1.*(random->uniform()*(extrude_length_max - 2.*rad) + rad + extrude_length_min);
635
636 vectorScalarMult3D(normalvec,r,ext);
637 vectorAdd3D(pos,ext,pos);
638
639 }
640
641 /* ----------------------------------------------------------------------
642 generate random positions on shallow copy insertion face
643 extrude by random length in negative face normal direction
644 currently only implemented for all_in_flag = 0
645 since would be tedious to check/implement otherwise
646 ------------------------------------------------------------------------- */
647
generate_random_global(double * pos)648 inline void FixInsertStream::generate_random_global(double *pos)
649 {
650 double r, ext[3];
651
652 // generate random position on the mesh
653
654 ins_face->generateRandomOwnedGhost(pos);
655
656 // extrude the position
657 r = -1.*(random->uniform()*extrude_length);
658 vectorScalarMult3D(normalvec,r,ext);
659 vectorAdd3D(pos,ext,pos);
660
661 }
662
663 /* ----------------------------------------------------------------------
664 generate random positions within extruded face
665 perform overlap check via xnear if requested
666 returns # bodies and # spheres that could actually be inserted
667 ------------------------------------------------------------------------- */
668
x_v_omega(int ninsert_this_local,int & ninserted_this_local,int & ninserted_spheres_this_local,double & mass_inserted_this_local)669 void FixInsertStream::x_v_omega(int ninsert_this_local,int &ninserted_this_local, int &ninserted_spheres_this_local, double &mass_inserted_this_local)
670 {
671 ninserted_this_local = ninserted_spheres_this_local = 0;
672 mass_inserted_this_local = 0.;
673
674 int nins;
675 double pos[3];
676 ParticleToInsert *pti;
677
678 double omega_tmp[] = {0.,0.,0.};
679
680 int ntry = 0;
681 int maxtry = ninsert_this_local * maxattempt;
682
683 // no overlap check
684 // insert with v_normal, no omega
685 if(!check_ol_flag)
686 {
687 for(int itotal = 0; itotal < ninsert_this_local; itotal++)
688 {
689 pti = fix_distribution->pti_list[ninserted_this_local];
690 double rad_to_insert = pti->r_bound_ins;
691
692 do
693 {
694 generate_random(pos,rad_to_insert);
695 ntry++;
696 }
697
698 while(ntry < maxtry && (!domain->is_in_subdomain(pos)));
699
700 if(ntry < maxtry)
701 {
702 // randomize quat here
703
704 if(quat_random_)
705 MathExtraLiggghts::random_unit_quat(random,quat_insert);
706
707 nins = pti->set_x_v_omega(pos,v_normal,omega_tmp,quat_insert);
708
709 ninserted_spheres_this_local += nins;
710 mass_inserted_this_local += pti->mass_ins;
711 ninserted_this_local++;
712 }
713 }
714 }
715 // overlap check
716 // account for maxattempt
717 // pti checks against xnear and adds self contributions
718 else
719 {
720 while(ntry < maxtry && ninserted_this_local < ninsert_this_local)
721 {
722 pti = fix_distribution->pti_list[ninserted_this_local];
723 double rad_to_insert = pti->r_bound_ins;
724
725 nins = 0;
726 while(nins == 0 && ntry < maxtry)
727 {
728 do
729 {
730 generate_random(pos,rad_to_insert);
731 ntry++;
732
733 }
734 while(ntry < maxtry && ((!domain->is_in_subdomain(pos)) || (domain->dist_subbox_borders(pos) < rad_to_insert)));
735
736 if(ntry < maxtry)
737 {
738 // randomize quat here
739 if(quat_random_)
740 MathExtraLiggghts::random_unit_quat(random,quat_insert);
741
742 nins = pti->check_near_set_x_v_omega(pos,v_normal,omega_tmp,quat_insert,neighList);
743 }
744 }
745
746 if(nins > 0)
747 {
748 ninserted_spheres_this_local += nins;
749 mass_inserted_this_local += pti->mass_ins;
750 ninserted_this_local++;
751 }
752 }
753 }
754
755 }
756
757 /* ---------------------------------------------------------------------- */
758
finalize_insertion(int ninserted_spheres_this_local)759 void FixInsertStream::finalize_insertion(int ninserted_spheres_this_local)
760 {
761 // nins particles have been inserted on this proc, set initial position, insertion step and release step according to pos
762
763 int n_steps = -1;
764 int step = update->ntimestep;
765 int ilo = atom->nlocal - ninserted_spheres_this_local;
766 int ihi = atom->nlocal;
767
768 double pos_rel[3], dist_normal;
769 double **x = atom->x;
770 double dt = update->dt;
771
772 double **release_data = fix_release->array_atom;
773
774 Multisphere *multisphere = NULL;
775 if(fix_multisphere) multisphere = &fix_multisphere->data();
776
777 for(int i = ilo; i < ihi; i++)
778 {
779
780 if(multisphere)
781 n_steps = fix_multisphere->calc_n_steps(i,p_ref,normalvec,v_normal);
782 if(!multisphere || n_steps == -1)
783 {
784 vectorSubtract3D(p_ref,x[i],pos_rel);
785 dist_normal = vectorDot3D(pos_rel,normalvec);
786 n_steps = static_cast<int>((dist_normal+FIX_INSERT_STREAM_TINY)/(vectorMag3D(v_normal)*dt));
787 }
788
789 // first 3 values is original position to integrate
790 vectorCopy3D(x[i],release_data[i]);
791
792 // 4th value is insertion step
793 release_data[i][3] = static_cast<double>(step);
794
795 // 5th value is step to release
796 release_data[i][4] = static_cast<double>(step + n_steps);
797
798 // 6-8th value is integration velocity
799 vectorCopy3D(v_normal,&release_data[i][5]);
800
801 // set initial conditions
802 // randomize vel, omega here
803 double v_toInsert[3],omega_toInsert[3];
804
805 vectorCopy3D(v_insert,v_toInsert);
806 vectorCopy3D(omega_insert,omega_toInsert);
807
808 // could randomize vel, omega here
809 generate_random_velocity(v_toInsert);
810
811 // 9-11th value is velocity, 12-14 is omega
812 vectorCopy3D(v_toInsert,&release_data[i][8]);
813 vectorCopy3D(omega_toInsert,&release_data[i][11]);
814
815 }
816
817 if(ntracer)
818 for(int i = 0; i < ntracer; i++)
819 tracer[i]->mark_tracers(ilo,ihi);
820 }
821
822 /* ---------------------------------------------------------------------- */
823
end_of_step()824 void FixInsertStream::end_of_step()
825 {
826 int r_step, i_step;
827
828 int step = update->ntimestep;
829 int nlocal = atom->nlocal;
830 double **release_data = fix_release->array_atom;
831 double time_elapsed, dist_elapsed[3], v_integrate[3], *v_toInsert, *omega_toInsert;
832 double dt = update->dt;
833
834 double **x = atom->x;
835 double **v = atom->v;
836 double **f = atom->f;
837 double **omega = atom->omega;
838 double **torque = atom->torque;
839 int *mask = atom->mask;
840
841 // only one fix handles the integration
842 if(!i_am_integrator) return;
843
844 for(int i = 0; i < nlocal; i++)
845 {
846 if (mask[i] & groupbit)
847 {
848 if(MathExtraLiggghts::compDouble(release_data[i][3],0.,1.e-13))
849 continue;
850
851 i_step = static_cast<int>(release_data[i][3]+FIX_INSERT_STREAM_TINY);
852 r_step = static_cast<int>(release_data[i][4]+FIX_INSERT_STREAM_TINY);
853 vectorCopy3D(&release_data[i][5],v_integrate);
854
855 if(step > r_step) continue;
856 else if (r_step == step)
857 {
858
859 if(fix_multisphere && fix_multisphere->belongs_to(i) >= 0)
860 {
861
862 v_toInsert = &release_data[i][8];
863 omega_toInsert = &release_data[i][11];
864 fix_multisphere->release(i,v_toInsert,omega_toInsert);
865 continue;
866 }
867
868 // integrate with constant vel and set v,omega
869
870 time_elapsed = (step - i_step) * dt;
871
872 // particle moves with v_integrate
873 vectorScalarMult3D(v_integrate,time_elapsed,dist_elapsed);
874 double *x_ins = release_data[i];
875
876 // set x,v,omega
877 // zero out force, torque
878
879 vectorAdd3D(x_ins,dist_elapsed,x[i]);
880
881 vectorZeroize3D(f[i]);
882 vectorZeroize3D(torque[i]);
883
884 v_toInsert = &release_data[i][8];
885 omega_toInsert = &release_data[i][11];
886
887 vectorCopy3D(v_toInsert,v[i]);
888 vectorCopy3D(omega_toInsert,omega[i]);
889
890 }
891 // step < r_step, only true for inserted particles
892 // b/c r_step is 0 for all other particles
893 // integrate with constant vel
894 else
895 {
896
897 time_elapsed = (step - i_step) * dt;
898
899 // particle moves with v_integrate
900 vectorScalarMult3D(v_integrate,time_elapsed,dist_elapsed);
901 double *x_ins = release_data[i];
902
903 // set x,v,omega
904 vectorAdd3D(x_ins,dist_elapsed,x[i]);
905
906 vectorCopy3D(v_integrate,v[i]);
907 vectorZeroize3D(omega[i]);
908
909 // zero out force, torque
910 vectorZeroize3D(f[i]);
911 vectorZeroize3D(torque[i]);
912 }
913 }
914 }
915
916 }
917
918 /* ---------------------------------------------------------------------- */
919
reset_timestep(bigint newstep,bigint oldstep)920 void FixInsertStream::reset_timestep(bigint newstep,bigint oldstep)
921 {
922 FixInsert::reset_timestep(newstep,oldstep);
923
924 reset_releasedata(newstep,oldstep);
925 }
926
927 /* ---------------------------------------------------------------------- */
928
reset_releasedata(bigint newstep,bigint oldstep)929 void FixInsertStream::reset_releasedata(bigint newstep,bigint oldstep)
930 {
931
932 int nlocal = atom->nlocal;
933 double **x = atom->x;
934 double **release_data = fix_release->array_atom;
935
936 for(int i = 0; i < nlocal; i++)
937 {
938
939 // first 3 values is original position to integrate
940 vectorCopy3D(x[i],release_data[i]);
941
942 // 4th value is insertion step
943 release_data[i][3] -= static_cast<double>(oldstep-newstep);
944
945 // 5th value is step to release
946 release_data[i][4] -= static_cast<double>(oldstep-newstep);
947
948 // 6-8th value is integration velocity
949 vectorCopy3D(v_normal,&release_data[i][5]);
950 }
951 }
952
953 /* ---------------------------------------------------------------------- */
954
recalc_release_restart()955 void FixInsertStream::recalc_release_restart()
956 {
957
958 const int nlocal = atom->nlocal;
959 double **x = atom->x;
960 double **release_data = fix_release->array_atom;
961 const double dt = update->dt;
962 double change_ratio = -1.;
963
964 for(int i = 0; i < nlocal; i++)
965 {
966 if(release_data[i][4] > update->ntimestep)
967 {
968 double dx[3];
969 vectorSubtract3D(x[i],release_data[i],dx);
970 const double dt_old = (vectorMag3D(dx)) / (vectorMag3D(&release_data[i][5]) * (static_cast<double>(update->ntimestep) - release_data[i][3]));
971
972 change_ratio = dt_old / dt;
973
974 release_data[i][4] = static_cast<double>(update->ntimestep) + static_cast<double>(static_cast<int>(change_ratio*(release_data[i][4] - static_cast<double>(update->ntimestep))));
975
976 }
977 }
978
979 recalc_release_ms = true;
980 dt_ratio = change_ratio;
981 }
982