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 This file is from LAMMPS, but has been modified. Copyright for
36 modification:
37
38 Copyright 2012- DCS Computing GmbH, Linz
39 Copyright 2009-2012 JKU Linz
40 Copyright 2016- CFDEMresearch GmbH, Linz
41
42 Copyright of original file:
43 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
44 http://lammps.sandia.gov, Sandia National Laboratories
45 Steve Plimpton, sjplimp@sandia.gov
46
47 Copyright (2003) Sandia Corporation. Under the terms of Contract
48 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
49 certain rights in this software. This software is distributed under
50 the GNU General Public License.
51 ------------------------------------------------------------------------- */
52
53 #include <cmath>
54 #include <stdio.h>
55 #include <stdlib.h>
56 #include <string.h>
57 #include "atom.h"
58 #include "atom_vec.h"
59 #include "domain.h"
60 #include "modify.h"
61 #include "force.h"
62 #include "update.h"
63 #include "modify.h"
64 #include "fix.h"
65 #include "fix_contact_history.h"
66 #include "comm.h"
67 #include "neighbor.h"
68 #include "neigh_list.h"
69 #include "neigh_request.h"
70 #include "memory.h"
71 #include "error.h"
72 #include "properties.h"
73 #include "fix_rigid.h"
74 #include "fix_particledistribution_discrete.h"
75 #include "fix_property_global.h"
76 #include "fix_property_atom.h"
77 #include "fix_contact_property_atom.h"
78 #include "compute_pair_gran_local.h"
79 #include "pair_gran.h"
80
81 using namespace LAMMPS_NS;
82
83 /* ---------------------------------------------------------------------- */
84
PairGran(LAMMPS * lmp)85 PairGran::PairGran(LAMMPS *lmp) : Pair(lmp)
86 {
87 single_enable = 0;
88 no_virial_fdotr_compute = 1;
89
90 suffix = NULL;
91 neighprev = 0;
92
93 history = 0;
94 dnum_pairgran = 0;
95 dnum_all = 0;
96 fix_history = NULL;
97
98 nfix = 0;
99 fix_dnum = NULL;
100 dnum_index = NULL;
101
102 mass_rigid = NULL;
103
104 nmax = 0;
105
106 cpl_enable = 1;
107 cpl_ = NULL;
108
109 energytrack_enable = 0;
110 fppaCPEn = fppaCDEn = fppaCPEt = fppaCDEVt = fppaCDEFt = fppaCTFW = fppaDEH = NULL;
111 CPEn = CDEn = CPEt = CDEVt = CDEFt = CTFW = DEH = NULL;
112
113 computeflag_ = 0;
114
115 needs_neighlist = true;
116
117 fix_contact_forces_ = NULL;
118 store_contact_forces_ = false;
119 store_contact_forces_every_ = 1;
120 fix_contact_forces_stress_ = NULL;
121 store_contact_forces_stress_ = false;
122 fix_store_multicontact_data_ = NULL;
123 store_multicontact_data_ = false;
124
125 fix_relax_ = NULL;
126
127 if(modify->n_fixes_style("multisphere/advanced"))
128 do_store_contact_forces();
129 }
130
131 /* ---------------------------------------------------------------------- */
132
~PairGran()133 PairGran::~PairGran()
134 {
135 if (fix_history) modify->delete_fix("contacthistory");
136 if (fix_contact_forces_) modify->delete_fix(fix_contact_forces_->id);
137 if (fix_contact_forces_stress_) modify->delete_fix(fix_contact_forces_stress_->id);
138 if (fix_store_multicontact_data_) modify->delete_fix(fix_store_multicontact_data_->id);
139
140 if (suffix) delete[] suffix;
141
142 if (allocated) {
143 memory->destroy(setflag);
144 memory->destroy(cutsq);
145 delete [] onerad_dynamic;
146 delete [] onerad_frozen;
147 delete [] maxrad_dynamic;
148 delete [] maxrad_frozen;
149 }
150
151 // tell cpl that pair gran is deleted
152 if(cpl_) cpl_->reference_deleted();
153
154 //unregister energy terms as property/atom
155 if (fppaCPEn) modify->delete_fix("CPEn");
156 if (fppaCDEn) modify->delete_fix("CDEn");
157 if (fppaCPEt) modify->delete_fix("CPEt");
158 if (fppaCDEVt) modify->delete_fix("CDEVt");
159 if (fppaCDEFt) modify->delete_fix("CDEFt");
160 if (fppaCTFW) modify->delete_fix("CTFW");
161 if (fppaDEH) modify->delete_fix("DEH");
162
163 if(fix_dnum) delete []fix_dnum;
164 if(dnum_index) delete []dnum_index;
165 }
166
167 /* ---------------------------------------------------------------------- */
168
updatePtrs()169 void PairGran::updatePtrs()
170 {
171 if(fppaCPEn) CPEn = fppaCPEn->vector_atom;
172 if(fppaCDEn) CDEn = fppaCDEn->vector_atom;
173 if(fppaCPEt) CPEt = fppaCPEt->vector_atom;
174 if(fppaCDEVt) CDEVt = fppaCDEVt->vector_atom;
175 if(fppaCDEFt) CDEFt = fppaCDEFt->vector_atom;
176 if(fppaCTFW) CTFW = fppaCTFW->vector_atom;
177 if(fppaDEH) DEH = fppaDEH->vector_atom;
178 }
179
180 /* ----------------------------------------------------------------------
181 allocate all arrays
182 ------------------------------------------------------------------------- */
183
allocate()184 void PairGran::allocate()
185 {
186 allocated = 1;
187 int n = atom->ntypes; //ensured elsewhere that this is high enough
188
189 memory->create(setflag,n+1,n+1,"pair:setflag");
190 for (int i = 1; i <= n; i++)
191 for (int j = i; j <= n; j++)
192 setflag[i][j] = 0;
193
194 memory->create(cutsq,n+1,n+1,"pair:cutsq");
195
196 onerad_dynamic = new double[n+1];
197 onerad_frozen = new double[n+1];
198 maxrad_dynamic = new double[n+1];
199 maxrad_frozen = new double[n+1];
200 }
201
202 /* ----------------------------------------------------------------------
203 set coeffs for one or more type pairs
204 ------------------------------------------------------------------------- */
205
coeff(int narg,char ** arg)206 void PairGran::coeff(int narg, char **arg)
207 {
208 if (narg > 2) error->all(FLERR,"Incorrect args for pair coefficients");
209 if (!allocated) allocate();
210
211 int ilo,ihi,jlo,jhi;
212 force->bounds(arg[0],atom->ntypes,ilo,ihi);
213 force->bounds(arg[1],atom->ntypes,jlo,jhi);
214
215 int count = 0;
216 for (int i = ilo; i <= ihi; i++) {
217 for (int j = MAX(jlo,i); j <= jhi; j++) {
218 setflag[i][j] = 1;
219 count++;
220 }
221 }
222
223 if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
224 }
225
226 /* ----------------------------------------------------------------------
227 init specific to this pair style
228 ------------------------------------------------------------------------- */
229
init_style()230 void PairGran::init_style()
231 {
232 int i;
233
234 // error and warning checks
235
236 if(0 == comm->me && 0 != neighbor->delay)
237 error->warning(FLERR,"It is heavily recommended to use 'neigh_modify delay 0' with granular pair styles");
238
239 if(strcmp(update->unit_style,"metal") ==0 || strcmp(update->unit_style,"real") == 0)
240 error->all(FLERR,"Cannot use a non-consistent unit system with pair gran. Please use si,cgs or lj.");
241
242 if (!atom->sphere_flag && !atom->superquadric_flag && !atom->shapetype_flag)
243 error->all(FLERR,"Pair granular requires atom style sphere, superquadric or convexhull");
244 if (comm->ghost_velocity == 0)
245 error->all(FLERR,"Pair granular requires ghost atoms store velocity");
246
247 // check if a fix rigid is registered - important for damp
248
249 fix_rigid = static_cast<FixRigid*>(modify->find_fix_style_strict("rigid",0));
250
251 if(modify->n_fixes_style("rigid") > 1)
252 error->warning(FLERR,"Pair gran does currently not support more than one fix rigid. This may result in under-damping.");
253
254 dt = update->dt;
255
256 // if shear history is stored:
257 // check if newton flag is valid
258 // if first init, create Fix needed for storing shear history
259
260 if (history && force->newton_pair == 1)
261 error->all(FLERR,"Pair granular with shear history requires newton pair off");
262
263 dnum_all = dnum_pairgran;
264
265 int dnum_extra = 0;
266 nfix = modify->nfix;
267
268 if(fix_dnum) delete []fix_dnum;
269 if(dnum_index) delete []dnum_index;
270 fix_dnum = new Fix*[nfix];
271 dnum_index = new int[nfix];
272
273 for(int ifix = 0; ifix < nfix; ifix++)
274 {
275 fix_dnum[ifix] = modify->fix[ifix];
276 dnum_index[ifix] = dnum_pairgran + dnum_extra;
277
278 dnum_extra += fix_dnum[ifix]->n_history_extra();
279 }
280
281 if(dnum_extra && !history)
282 error->all(FLERR,"Cannot use extra dnum with non-history pair style");
283 else
284 dnum_all += dnum_extra;
285
286 // init contact history
287 if(history)
288 {
289 if (!fix_history)
290 {
291
292 fix_history = static_cast<FixContactHistory*>(modify->find_fix_style_strict("contacthistory",0));
293 if(fix_history)
294 error->all(FLERR,"Can not use more than one pair style that uses contact history");
295
296 char **fixarg = new char*[4+2*dnum_all];
297 fixarg[3] = new char[3];
298
299 fixarg[0] = (char *) "contacthistory";
300 fixarg[1] = (char *) "all";
301 fixarg[2] = (char *) "contacthistory";
302 sprintf(fixarg[3],"%d",dnum_all);
303 history_args(&fixarg[4]);
304
305 // add history args from fixes
306 for(int ifix = 0; ifix < nfix; ifix++)
307 {
308 if(fix_dnum[ifix]->n_history_extra())
309 if(!fix_dnum[ifix]->history_args(&fixarg[4+2*dnum_index[ifix]]))
310 error->all(FLERR,"Missing history_args() implementation for fix that requests extra history");
311 }
312
313 modify->add_fix(4+2*dnum_all,fixarg);
314
315 if(modify->n_fixes_style_strict("contacthistory") != 1) error->all(FLERR,"Pair granular with shear history requires exactly one fix of style contacthistory");
316 fix_history = static_cast<FixContactHistory*>(modify->find_fix_style_strict("contacthistory",0));
317
318 delete [] fixarg[3];
319 delete [] fixarg;
320 }
321 }
322
323 // register per-particle properties for energy tracking
324 if(energytrack_enable)
325 {
326 if(comm->nprocs > 1) error->all(FLERR,"check communication settings");
327 char **fixarg = new char*[9];
328 if (fppaCPEn == NULL) {
329 fixarg[0]=(char *) "CPEn";
330 fixarg[1]=(char *) "all";
331 fixarg[2]=(char *) "property/atom";
332 fixarg[3]=(char *) "CPEn";
333 fixarg[4]=(char *) "scalar";
334 fixarg[5]=(char *) "yes";
335 fixarg[6]=(char *) "no";
336 fixarg[7]=(char *) "no";
337 fixarg[8]=(char *) "0.0";
338 modify->add_fix(9,fixarg);
339 fppaCPEn=static_cast<FixPropertyAtom*>(modify->find_fix_property("CPEn","property/atom","scalar",0,0,force->pair_style));
340 }
341 if (fppaCDEn == NULL) {
342 fixarg[0]=(char *) "CDEn";
343 fixarg[1]=(char *) "all";
344 fixarg[2]=(char *) "property/atom";
345 fixarg[3]=(char *) "CDEn";
346 fixarg[4]=(char *) "scalar";
347 fixarg[5]=(char *) "yes";
348 fixarg[6]=(char *) "no";
349 fixarg[7]=(char *) "no";
350 fixarg[8]=(char *) "0.0";
351 modify->add_fix(9,fixarg);
352 fppaCDEn=static_cast<FixPropertyAtom*>(modify->find_fix_property("CDEn","property/atom","scalar",0,0,force->pair_style));
353 }
354 if (fppaCPEt == NULL) {
355 fixarg[0]=(char *) "CPEt";
356 fixarg[1]=(char *) "all";
357 fixarg[2]=(char *) "property/atom";
358 fixarg[3]=(char *) "CPEt";
359 fixarg[4]=(char *) "scalar";
360 fixarg[5]=(char *) "yes";
361 fixarg[6]=(char *) "no";
362 fixarg[7]=(char *) "no";
363 fixarg[8]=(char *) "0.0";
364 modify->add_fix(9,fixarg);
365 fppaCPEt=static_cast<FixPropertyAtom*>(modify->find_fix_property("CPEt","property/atom","scalar",0,0,force->pair_style));
366 }
367 if (fppaCDEVt == NULL) {
368 fixarg[0]=(char *) "CDEVt";
369 fixarg[1]=(char *) "all";
370 fixarg[2]=(char *) "property/atom";
371 fixarg[3]=(char *) "CDEVt";
372 fixarg[4]=(char *) "scalar";
373 fixarg[5]=(char *) "yes";
374 fixarg[6]=(char *) "no";
375 fixarg[7]=(char *) "no";
376 fixarg[8]=(char *) "0.0";
377 modify->add_fix(9,fixarg);
378 fppaCDEVt=static_cast<FixPropertyAtom*>(modify->find_fix_property("CDEVt","property/atom","scalar",0,0,force->pair_style));
379 }
380 if (fppaCDEFt == NULL) {
381 fixarg[0]=(char *) "CDEFt";
382 fixarg[1]=(char *) "all";
383 fixarg[2]=(char *) "property/atom";
384 fixarg[3]=(char *) "CDEFt";
385 fixarg[4]=(char *) "scalar";
386 fixarg[5]=(char *) "yes";
387 fixarg[6]=(char *) "no";
388 fixarg[7]=(char *) "no";
389 fixarg[8]=(char *) "0.0";
390 modify->add_fix(9,fixarg);
391 fppaCDEFt=static_cast<FixPropertyAtom*>(modify->find_fix_property("CDEFt","property/atom","scalar",0,0,force->pair_style));
392 }
393 if (fppaCTFW == NULL) {
394 fixarg[0]=(char *) "CTFW";
395 fixarg[1]=(char *) "all";
396 fixarg[2]=(char *) "property/atom";
397 fixarg[3]=(char *) "CTFW";
398 fixarg[4]=(char *) "scalar";
399 fixarg[5]=(char *) "yes";
400 fixarg[6]=(char *) "no";
401 fixarg[7]=(char *) "no";
402 fixarg[8]=(char *) "0.0";
403 modify->add_fix(9,fixarg);
404 fppaCTFW=static_cast<FixPropertyAtom*>(modify->find_fix_property("CTFW","property/atom","scalar",0,0,force->pair_style));
405 }
406 if (fppaDEH == NULL) {
407 fixarg[0]=(char *) "DEH";
408 fixarg[1]=(char *) "all";
409 fixarg[2]=(char *) "property/atom";
410 fixarg[3]=(char *) "DEH";
411 fixarg[4]=(char *) "scalar";
412 fixarg[5]=(char *) "yes";
413 fixarg[6]=(char *) "no";
414 fixarg[7]=(char *) "no";
415 fixarg[8]=(char *) "0.0";
416 modify->add_fix(9,fixarg);
417 fppaDEH=static_cast<FixPropertyAtom*>(modify->find_fix_property("DEH","property/atom","scalar",0,0,force->pair_style));
418 }
419 delete []fixarg;
420
421 if(force->newton_pair == 1) error->all(FLERR,"Have to implement rev comm of energy terms");
422 }
423
424 // register per-particle properties for energy tracking
425 if(store_contact_forces_ && fix_contact_forces_ == NULL)
426 {
427 char **fixarg = new char*[19];
428 fixarg[0]=(char *) "contactforces";
429 fixarg[1]=(char *) "all";
430 fixarg[2]=(char *) "contactproperty/atom";
431 fixarg[3]=(char *) "contactforces";
432 fixarg[4]=(char *) "6";
433 fixarg[5]=(char *) "fx";
434 fixarg[6]=(char *) "0";
435 fixarg[7]=(char *) "fy";
436 fixarg[8]=(char *) "0";
437 fixarg[9]=(char *) "fz";
438 fixarg[10]=(char *) "0";
439 fixarg[11]=(char *) "tx";
440 fixarg[12]=(char *) "0";
441 fixarg[13]=(char *) "ty";
442 fixarg[14]=(char *) "0";
443 fixarg[15]=(char *) "tz";
444 fixarg[16]=(char *) "0";
445 fixarg[17]=(char *) "reset";
446 fixarg[18]=(char *) "yes";
447 modify->add_fix(19,fixarg);
448 fix_contact_forces_ = static_cast<FixContactPropertyAtom*>(modify->find_fix_id("contactforces"));
449 delete []fixarg;
450 }
451
452 // register per-particle properties for energy tracking
453 if(store_contact_forces_stress_ && fix_contact_forces_stress_ == NULL)
454 {
455 char **fixarg = new char*[15];
456 fixarg[0]=(char *) "contactforces_stress_";
457 fixarg[1]=(char *) "all";
458 fixarg[2]=(char *) "contactproperty/atom";
459 fixarg[3]=(char *) "contactforces_stress_";
460 fixarg[4]=(char *) "4";
461 fixarg[5]=(char *) "fx";
462 fixarg[6]=(char *) "0";
463 fixarg[7]=(char *) "fy";
464 fixarg[8]=(char *) "0";
465 fixarg[9]=(char *) "fz";
466 fixarg[10]=(char *) "0";
467 fixarg[11]=(char *) "j";
468 fixarg[12]=(char *) "0";
469 fixarg[13]=(char *) "reset";
470 fixarg[14]=(char *) "yes";
471 modify->add_fix(15,fixarg);
472 fix_contact_forces_stress_ = static_cast<FixContactPropertyAtom*>(modify->find_fix_id("contactforces_stress_"));
473 delete []fixarg;
474 }
475
476 if(store_multicontact_data_ && fix_store_multicontact_data_ == NULL)
477 {
478 // create a new per contact property which will contain the data for the computation according to Brodu et. al. 2016
479 // surfPosIJ will contain the position of the contact surface ij, realtive to position i
480 // normalForce will contain the normal component of the contact force
481 char **fixarg = new char*[15];
482 fixarg[0]=(char *) "multicontactData_";
483 fixarg[1]=(char *) "all";
484 fixarg[2]=(char *) "contactproperty/atom";
485 fixarg[3]=(char *) "multicontactData_";
486 fixarg[4]=(char *) "4";
487 fixarg[5]=(char *) "surfPosIJ_x";
488 fixarg[6]=(char *) "0";
489 fixarg[7]=(char *) "surfPosIJ_y";
490 fixarg[8]=(char *) "0";
491 fixarg[9]=(char *) "surfPosIJ_z";
492 fixarg[10]=(char *) "0";
493 fixarg[11]=(char *) "normalForce";
494 fixarg[12]=(char *) "0";
495 fixarg[13]=(char *) "reset";
496 fixarg[14]=(char *) "no";
497 modify->add_fix(15,fixarg);
498 fix_store_multicontact_data_ = static_cast<FixContactPropertyAtom*>(modify->find_fix_id("multicontactData_"));
499 delete []fixarg;
500 }
501
502 fix_sum_normal_force_ =
503 static_cast<FixPropertyAtom*>
504 (
505 modify->find_fix_property("sum_normal_force_","property/atom","scalar",0,0, "pair/gran", false)
506 );
507
508 // need a gran neigh list and optionally a granular history neigh list
509
510 if(needs_neighlist)
511 {
512 int irequest = neighbor->request(this);
513 neighbor->requests[irequest]->pairgran_hashcode = hashcode();
514 neighbor->requests[irequest]->half = 0;
515 neighbor->requests[irequest]->gran = 1;
516
517 if (history) {
518 irequest = neighbor->request(this);
519 neighbor->requests[irequest]->pairgran_hashcode = hashcode();
520 neighbor->requests[irequest]->id = 1;
521 neighbor->requests[irequest]->half = 0;
522 neighbor->requests[irequest]->granhistory = 1;
523 neighbor->requests[irequest]->dnum = dnum_all;
524 }
525 }
526
527 // check for freeze Fix and set freeze_group_bit_
528
529 for (i = 0; i < modify->nfix; i++)
530 if (strcmp(modify->fix[i]->style,"freeze") == 0) break;
531 if (i < modify->nfix) freeze_group_bit_ = modify->fix[i]->groupbit;
532 else freeze_group_bit_ = 0;
533
534 // set maxrad_dynamic and maxrad_frozen for each type
535 for (i = 1; i <= atom->ntypes; i++)
536 onerad_dynamic[i] = onerad_frozen[i] = 0.0;
537
538 // include future particles as dynamic
539
540 for (i = 0; i < modify->nfix; i++)
541 {
542
543 if(!modify->fix[i]->use_rad_for_cut_neigh_and_ghost())
544 continue;
545
546 for(int j=1;j<=atom->ntypes;j++)
547 {
548 int pour_type = 0;
549 double pour_maxrad = 0.0;
550 pour_type = j;
551 pour_maxrad = modify->fix[i]->max_rad(pour_type);
552 onerad_dynamic[pour_type] = MAX(onerad_dynamic[pour_type],pour_maxrad);
553 }
554 }
555
556 // further dynamic and frozen
557
558 double *radius = atom->radius;
559 int *mask = atom->mask;
560 int *type = atom->type;
561 int nlocal = atom->nlocal;
562
563 for (i = 0; i < nlocal; i++)
564 if (mask[i] & freeze_group_bit_)
565 onerad_frozen[type[i]] = MAX(onerad_frozen[type[i]],radius[i]);
566 else
567 onerad_dynamic[type[i]] = MAX(onerad_dynamic[type[i]],radius[i]);
568
569 MPI_Allreduce(&onerad_dynamic[1],&maxrad_dynamic[1],atom->ntypes,MPI_DOUBLE,MPI_MAX,world);
570 MPI_Allreduce(&onerad_frozen[1],&maxrad_frozen[1],atom->ntypes,MPI_DOUBLE,MPI_MAX,world);
571
572 init_granular();
573 }
574
575 /* ----------------------------------------------------------------------
576 register and unregister callback to compute
577 ------------------------------------------------------------------------- */
578
register_compute_pair_local(ComputePairGranLocal * ptr,int & dnum_compute)579 void PairGran::register_compute_pair_local(ComputePairGranLocal *ptr,int &dnum_compute)
580 {
581 if(cpl_ != NULL) error->all(FLERR,"Pair gran allows only one compute of type pair/local");
582 cpl_ = ptr;
583 dnum_compute = dnum_pairgran; //history values
584 }
585
unregister_compute_pair_local(ComputePairGranLocal * ptr)586 void PairGran::unregister_compute_pair_local(ComputePairGranLocal *ptr)
587 {
588 if(cpl_ != ptr) error->all(FLERR,"Illegal situation in PairGran::unregister_compute_pair_local");
589 cpl_ = NULL;
590 }
591
592 /* ----------------------------------------------------------------------
593 return index for extra dnum
594 ------------------------------------------------------------------------- */
595
fix_extra_dnum_index(class Fix * fix)596 int PairGran::fix_extra_dnum_index(class Fix *fix)
597 {
598 for(int ifix = 0; ifix < nfix; ifix++)
599 {
600 if(fix == fix_dnum[ifix])
601 return dnum_index[ifix];
602 }
603 error->all(FLERR,"Internal error - illegal fix_extra_dnum_index() call");
604 return 0;
605 }
606
607 /* ----------------------------------------------------------------------
608 neighbor callback to inform pair style of neighbor list to use
609 optional granular history list
610 ------------------------------------------------------------------------- */
611
init_list(int id,NeighList * ptr)612 void PairGran::init_list(int id, NeighList *ptr)
613 {
614
615 if (id == 0) list = ptr;
616 else if (id == 1) listgranhistory = ptr;
617 }
618
619 /* ----------------------------------------------------------------------
620 init for one type pair i,j and corresponding j,i
621 ------------------------------------------------------------------------- */
622
init_one(int i,int j)623 double PairGran::init_one(int i, int j)
624 {
625 if (!allocated) allocate();
626
627 // cutoff = sum of max I,J radii for
628 // dynamic/dynamic & dynamic/frozen interactions, but not frozen/frozen
629
630 double cutoff = maxrad_dynamic[i]+maxrad_dynamic[j];
631 cutoff = MAX(cutoff,maxrad_frozen[i]+maxrad_dynamic[j]);
632 cutoff = MAX(cutoff,maxrad_dynamic[i]+maxrad_frozen[j]);
633
634 return cutoff * neighbor->contactDistanceFactor;
635 }
636
637 /* ----------------------------------------------------------------------
638 compute as called via force
639 ------------------------------------------------------------------------- */
640
compute(int eflag,int vflag)641 void PairGran::compute(int eflag, int vflag)
642 {
643 if(forceoff()) return;
644
645 // update rigid body info for owned & ghost atoms if using FixRigid masses
646 // body[i] = which body atom I is in, -1 if none
647 // mass_body = mass of each rigid body
648
649 if (fix_rigid && neighbor->ago == 0) {
650 int tmp;
651 int *body = (int *) fix_rigid->extract("body",tmp);
652 double *mass_body = (double *) fix_rigid->extract("masstotal",tmp);
653 if (atom->nmax > nmax) {
654 memory->destroy(mass_rigid);
655 nmax = atom->nmax;
656 memory->create(mass_rigid,nmax,"pair:mass_rigid");
657 }
658 int nlocal = atom->nlocal;
659 for (int i = 0; i < nlocal; i++)
660 if (body[i] >= 0) mass_rigid[i] = mass_body[body[i]];
661 else mass_rigid[i] = 0.0;
662 comm->forward_comm_pair(this);
663 }
664
665 computeflag_ = 1;
666 shearupdate_ = 1;
667 if (update->setupflag) shearupdate_ = 0;
668
669 compute_force(eflag,vflag,0);
670 }
671
672 /* ----------------------------------------------------------------------
673 compute as called via compute pair gran local
674 ------------------------------------------------------------------------- */
675
compute_pgl(int eflag,int vflag)676 void PairGran::compute_pgl(int eflag, int vflag)
677 {
678
679 // update rigid body info for owned & ghost atoms if using FixRigid masses
680 // body[i] = which body atom I is in, -1 if none
681 // mass_body = mass of each rigid body
682
683 if (fix_rigid && neighbor->ago == 0) {
684 int tmp;
685 int *body = (int *) fix_rigid->extract("body",tmp);
686 double *mass_body = (double *) fix_rigid->extract("masstotal",tmp);
687 if (atom->nmax > nmax) {
688 memory->destroy(mass_rigid);
689 nmax = atom->nmax;
690 memory->create(mass_rigid,nmax,"pair:mass_rigid");
691 }
692 int nlocal = atom->nlocal;
693 for (int i = 0; i < nlocal; i++)
694 if (body[i] >= 0) mass_rigid[i] = mass_body[body[i]];
695 else mass_rigid[i] = 0.0;
696 comm->forward_comm_pair(this);
697 }
698
699 bool reset_computeflag = (computeflag_ == 1) ? true : false;
700
701 computeflag_ = 0;
702 shearupdate_ = 0;
703
704 compute_force(eflag,vflag,1);
705
706 if(reset_computeflag)
707 computeflag_ = 1;
708 }
709
710 /* ---------------------------------------------------------------------- */
711
pack_comm(int n,int * list,double * buf,int pbc_flag,int * pbc)712 int PairGran::pack_comm(int n, int *list,double *buf, int pbc_flag, int *pbc)
713 {
714 int i,j,m;
715
716 m = 0;
717 for (i = 0; i < n; i++) {
718 j = list[i];
719 buf[m++] = mass_rigid[j];
720 }
721 return 1;
722 }
723
724 /* ---------------------------------------------------------------------- */
725
unpack_comm(int n,int first,double * buf)726 void PairGran::unpack_comm(int n, int first, double *buf)
727 {
728 int i,m,last;
729
730 m = 0;
731 last = first + n;
732 for (i = first; i < last; i++)
733 mass_rigid[i] = static_cast<int> (buf[m++]);
734 }
735
736 /* ----------------------------------------------------------------------
737 proc 0 writes to restart file
738 ------------------------------------------------------------------------- */
739
write_restart(FILE * fp)740 void PairGran::write_restart(FILE *fp)
741 {
742 write_restart_settings(fp);
743
744 int i,j;
745 for (i = 1; i <= atom->ntypes; i++)
746 for (j = i; j <= atom->ntypes; j++)
747 fwrite(&setflag[i][j],sizeof(int),1,fp);
748 }
749
750 /* ----------------------------------------------------------------------
751 proc 0 reads from restart file, bcasts
752 ------------------------------------------------------------------------- */
753
read_restart(FILE * fp,const int major,const int minor)754 void PairGran::read_restart(FILE *fp, const int major, const int minor)
755 {
756 read_restart_settings(fp, major, minor);
757 allocate();
758
759 int i,j;
760 int me = comm->me;
761 for (i = 1; i <= atom->ntypes; i++)
762 for (j = i; j <= atom->ntypes; j++) {
763 if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
764 MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
765 }
766 }
767
768 /* ---------------------------------------------------------------------- */
769
reset_dt()770 void PairGran::reset_dt()
771 {
772 dt = update->dt;
773 }
774
775 /* ---------------------------------------------------------------------- */
776
extract(const char * str,int & dim)777 void *PairGran::extract(const char *str, int &dim)
778 {
779 dim = 0;
780 if (strcmp(str,"computeflag") == 0) return (void *) &computeflag_;
781 return NULL;
782 }
783
784 /* ----------------------------------------------------------------------
785 memory usage of local atom-based arrays
786 ------------------------------------------------------------------------- */
787
memory_usage()788 double PairGran::memory_usage()
789 {
790 double bytes = nmax * sizeof(double);
791 return bytes;
792 }
793
forceoff()794 bool PairGran::forceoff() {
795 return false;
796 }
797
cpl_add_pair(LCM::SurfacesIntersectData & sidata,LCM::ForceData & i_forces)798 void PairGran::cpl_add_pair(LCM::SurfacesIntersectData & sidata, LCM::ForceData & i_forces)
799 {
800 const double fx = i_forces.delta_F[0];
801 const double fy = i_forces.delta_F[1];
802 const double fz = i_forces.delta_F[2];
803 const double tor1 = i_forces.delta_torque[0];
804 const double tor2 = i_forces.delta_torque[1];
805 const double tor3 = i_forces.delta_torque[2];
806 const double * const contact_point =
807 #ifdef NONSPHERICAL_ACTIVE_FLAG
808 atom->shapetype_flag ? sidata.contact_point : NULL;
809 #else
810 NULL;
811 #endif
812 cpl_->add_pair(sidata.i, sidata.j, fx,fy,fz,tor1,tor2,tor3,sidata.contact_history, contact_point);
813 }
814
cpl_pair_finalize()815 void PairGran::cpl_pair_finalize()
816 {
817 cpl_->pair_finalize();
818 }
819