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 <mpi.h>
43 #include <string.h>
44 #include <stdio.h>
45 #include <cmath>
46 #include <algorithm>
47 #include "fix_contact_history_mesh.h"
48 #include "atom.h"
49 #include "fix_mesh_surface.h"
50 #include "comm.h"
51 #include "neighbor.h"
52 #include "neigh_list.h"
53 #include "force.h"
54 #include "pair.h"
55 #include "update.h"
56 #include "modify.h"
57 #include "memory.h"
58 #include "error.h"
59
60 #if defined(_OPENMP)
61 #include "omp.h"
62 #endif
63
64 using namespace LAMMPS_NS;
65 using namespace FixConst;
66
67 /* ---------------------------------------------------------------------- */
68
FixContactHistoryMesh(LAMMPS * lmp,int narg,char ** arg)69 FixContactHistoryMesh::FixContactHistoryMesh(LAMMPS *lmp, int narg, char **arg) :
70 FixContactHistory(lmp, narg, arg),
71 ipage1_(0),
72 dpage1_(0),
73 ipage2_(0),
74 dpage2_(0),
75 keeppage_(0),
76 intersectpage_(0),
77 keepflag_(0),
78 intersectflag_(0),
79 mesh_(0),
80 fix_neighlist_mesh_(0),
81 fix_nneighs_(0),
82 build_neighlist_(true),
83 numpages_(0)
84 {
85
86 // parse args
87 Fix *f = modify->find_fix_id(arg[iarg_++]);
88 if(!f || strncmp(f->style,"mesh/surface",12) )
89 error->fix_error(FLERR,this,"wrong ID for fix mesh/surface");
90 mesh_ = (static_cast<FixMeshSurface*>(f))->triMesh();
91 fix_neighlist_mesh_ = (static_cast<FixMeshSurface*>(f))->meshNeighlist();
92
93 swap_ = new double[dnum_];
94
95 // initial allocation of keepflag
96 keepflag_ = (bool **) memory->srealloc(keepflag_,atom->nmax*sizeof(bool *),
97 "contact_history:keepflag");
98 // initial allocation of intersectflag
99 intersectflag_ = (bool **) memory->srealloc(intersectflag_,atom->nmax*sizeof(bool *),
100 "contact_history:intersectflag");
101 }
102
103 /* ---------------------------------------------------------------------- */
104
~FixContactHistoryMesh()105 FixContactHistoryMesh::~FixContactHistoryMesh()
106 {
107 // delete locally stored arrays
108
109 if(ipage1_) delete [] ipage1_;
110 if(dpage1_) delete [] dpage1_;
111 if(ipage2_) delete [] ipage2_;
112 if(dpage2_) delete [] dpage2_;
113
114 if(keeppage_) {
115 for(int i = 0; i < numpages_; i++) {
116 delete keeppage_[i];
117 keeppage_[i] = NULL;
118 }
119 delete [] keeppage_;
120 keeppage_ = NULL;
121 }
122 if(intersectpage_) {
123 for(int i = 0; i < numpages_; i++) {
124 delete intersectpage_[i];
125 intersectpage_[i] = NULL;
126 }
127 delete [] intersectpage_;
128 intersectpage_ = NULL;
129 }
130
131 ipage_ = 0;
132 dpage_ = 0;
133
134 delete [] swap_;
135
136 if(keepflag_) memory->sfree(keepflag_);
137 if(intersectflag_) memory->sfree(intersectflag_);
138 }
139
140 /* ---------------------------------------------------------------------- */
141
setmask()142 int FixContactHistoryMesh::setmask()
143 {
144 int mask = 0;
145 mask |= PRE_FORCE;
146 mask |= MIN_PRE_FORCE;
147 mask |= PRE_NEIGHBOR;
148 mask |= PRE_EXCHANGE;
149 mask |= MIN_PRE_EXCHANGE;
150 return mask;
151 }
152
153 /* ---------------------------------------------------------------------- */
154
init()155 void FixContactHistoryMesh::init()
156 {
157 FixContactHistory::init();
158
159 char *fix_nneighs_name = new char[strlen(mesh_->mesh_id())+1+14];
160 sprintf(fix_nneighs_name,"n_neighs_mesh_%s",mesh_->mesh_id());
161
162 fix_nneighs_ = static_cast<FixPropertyAtom*>(modify->find_fix_property(fix_nneighs_name,"property/atom","scalar",0,0,this->style));
163
164 delete [] fix_nneighs_name;
165 }
166
167 /* ----------------------------------------------------------------------
168 create pages if first time or if neighbor pgsize/oneatom has changed
169 note that latter could cause shear history info to be discarded
170 ------------------------------------------------------------------------- */
171
allocate_pages()172 void FixContactHistoryMesh::allocate_pages()
173 {
174 if ((ipage_ == NULL) || (pgsize_ != neighbor->pgsize) || (oneatom_ != neighbor->oneatom)) {
175
176 bool use_first = ipage_ == ipage2_;
177
178 delete [] ipage1_;
179 delete [] dpage1_;
180 delete [] ipage2_;
181 delete [] dpage2_;
182
183 if(keeppage_) {
184 for(int i = 0; i < numpages_; i++) {
185 delete keeppage_[i];
186 keeppage_[i] = NULL;
187 }
188 delete [] keeppage_;
189 keeppage_ = NULL;
190 }
191 if(intersectpage_) {
192 for(int i = 0; i < numpages_; i++) {
193 delete intersectpage_[i];
194 intersectpage_[i] = NULL;
195 }
196 delete [] intersectpage_;
197 intersectpage_ = NULL;
198 }
199
200 pgsize_ = neighbor->pgsize;
201 oneatom_ = neighbor->oneatom;
202 numpages_ = comm->nthreads;
203 ipage1_ = new MyPage<int>[numpages_];
204 dpage1_ = new MyPage<double>[numpages_];
205 ipage2_ = new MyPage<int>[numpages_];
206 dpage2_ = new MyPage<double>[numpages_];
207 keeppage_ = new MyPage<bool>*[numpages_];
208 intersectpage_ = new MyPage<bool>*[numpages_];
209 for (int i = 0; i < numpages_; i++) {
210 ipage1_[i].init(oneatom_,pgsize_);
211 dpage1_[i].init(oneatom_*std::max(1,dnum_),pgsize_);
212 ipage2_[i].init(oneatom_,pgsize_);
213 dpage2_[i].init(oneatom_*std::max(1,dnum_),pgsize_);
214 }
215
216 #if defined(_OPENMP)
217 #pragma omp parallel
218 {
219 const int tid = omp_get_thread_num();
220 // make sure page is allocated in memory near core
221 keeppage_[tid] = new MyPage<bool>();
222 keeppage_[tid]->init(oneatom_,pgsize_);
223 intersectpage_[tid] = new MyPage<bool>();
224 intersectpage_[tid]->init(oneatom_,pgsize_);
225 }
226 #else
227 for (int i = 0; i < numpages_; i++) {
228 keeppage_[i] = new MyPage<bool>();
229 keeppage_[i]->init(oneatom_,pgsize_);
230 intersectpage_[i] = new MyPage<bool>();
231 intersectpage_[i]->init(oneatom_,pgsize_);
232 }
233 #endif
234
235 if(use_first)
236 {
237 ipage_ = ipage1_;
238 dpage_ = dpage1_;
239 }
240 else
241 {
242 ipage_ = ipage2_;
243 dpage_ = dpage2_;
244 }
245 }
246 }
247
248 /* ---------------------------------------------------------------------- */
249
setup_pre_neighbor()250 void FixContactHistoryMesh::setup_pre_neighbor()
251 {
252
253 }
254
255 /* ---------------------------------------------------------------------- */
256
setup_pre_exchange()257 void FixContactHistoryMesh::setup_pre_exchange()
258 {
259 pre_exchange();
260 }
261 /* ---------------------------------------------------------------------- */
262
min_setup_pre_exchange()263 void FixContactHistoryMesh::min_setup_pre_exchange()
264 {
265 pre_exchange();
266 }
267
268 /* ----------------------------------------------------------------------
269 sort contacthistory so need
270 ------------------------------------------------------------------------- */
271
pre_exchange()272 void FixContactHistoryMesh::pre_exchange()
273 {
274
275 if(!recent_restart)
276 sort_contacts();
277
278 // set maxtouch = max # of partners of any owned atom
279 // bump up comm->maxexchange_fix if necessary
280
281 const int nlocal = atom->nlocal;
282
283 maxtouch_ = 0;
284 for (int i = 0; i < nlocal; i++) maxtouch_ = MAX(maxtouch_,npartner_[i]);
285 comm->maxexchange_fix = MAX(comm->maxexchange_fix,(dnum_+1)*maxtouch_+1);
286 }
287
288 /* ----------------------------------------------------------------------
289 need to execute here since neighlist is refreshed in setup_pre_force(int foo)
290 if this is not done, data will get out-of-sync
291
292 need not care about overwriting hist values since always have pointer
293 to most current data
294 ------------------------------------------------------------------------- */
295
setup_pre_force(int dummy)296 void FixContactHistoryMesh::setup_pre_force(int dummy)
297 {
298 pre_neighbor();
299 pre_force(0);
300 }
301
302 /* ---------------------------------------------------------------------- */
303
pre_neighbor()304 void FixContactHistoryMesh::pre_neighbor()
305 {
306 build_neighlist_ = true;
307
308 }
309
310 /* ----------------------------------------------------------------------
311 swap 2 pages data structures where contacthistory is stored
312 do this here because exchange and sort can rag the data structure
313 ------------------------------------------------------------------------- */
314
pre_force(int dummy)315 void FixContactHistoryMesh::pre_force(int dummy)
316 {
317 if(!build_neighlist_)
318 return;
319 build_neighlist_ = false;
320
321 const int nlocal = atom->nlocal;
322
323 cleanUpContactJumps();
324
325 int nneighs_next;
326 int *partner_prev;
327 double *contacthistory_prev;
328
329 MyPage<int> *ipage_next = (ipage_ == ipage1_) ? ipage2_ : ipage1_;
330 MyPage<double> *dpage_next = (dpage_ == dpage1_) ? dpage2_ : dpage1_;
331
332 ipage_next->reset();
333 dpage_next->reset();
334
335 for (int i = 0; i < nlocal; i++)
336 {
337 nneighs_next = fix_nneighs_->get_vector_atom_int(i);
338
339 partner_prev = partner_[i];
340 contacthistory_prev = contacthistory_[i];
341
342 partner_[i] = ipage_next->get(nneighs_next);
343 if (!partner_[i])
344 error->one(FLERR,"mesh neighbor list overflow, boost neigh_modify one and/or page");
345 vectorInitializeN(partner_[i],nneighs_next,-1);
346
347 contacthistory_[i] = dpage_next->get(nneighs_next*dnum_);
348 if(!contacthistory_[i])
349 error->one(FLERR,"mesh neighbor list overflow, boost neigh_modify one");
350 vectorZeroizeN(contacthistory_[i],nneighs_next*dnum_);
351
352 if(npartner_[i] > nneighs_next)
353 {
354
355 error->one(FLERR,"internal error");
356 }
357
358 for(int ipartner = 0; ipartner < npartner_[i]; ipartner++)
359 {
360 if(partner_prev[ipartner] < 0)
361 error->one(FLERR,"internal error");
362
363 partner_[i][ipartner] = partner_prev[ipartner];
364 vectorCopyN(&(contacthistory_prev[ipartner*dnum_]),&(contacthistory_[i][ipartner*dnum_]),dnum_);
365
366 }
367 }
368
369 ipage_ = ipage_next;
370 dpage_ = dpage_next;
371 }
372
373 /* ---------------------------------------------------------------------- */
374
sort_contacts()375 void FixContactHistoryMesh::sort_contacts()
376 {
377 const int nlocal = atom->nlocal;
378 int nneighs, first_empty, last_filled;
379
380 for(int i = 0; i < nlocal; i++)
381 {
382 nneighs = fix_nneighs_->get_vector_atom_int(i);
383
384 if(0 == nneighs)
385 continue;
386
387 do
388 {
389 first_empty = last_filled = -1;
390
391 for(int j = 0; j < nneighs; j++)
392 {
393 if(-1 == first_empty && -1 == partner_[i][j])
394 first_empty = j;
395 if(partner_[i][j] >= 0)
396 last_filled = j;
397 }
398
399 if(first_empty > -1 && last_filled > -1 && first_empty < last_filled)
400 swap(i,first_empty,last_filled,true);
401 }
402 while(first_empty > -1 && last_filled > -1 && first_empty < last_filled);
403 }
404
405 }
406
407 /* ----------------------------------------------------------------------
408 mark all contacts for deletion, mark all as not intersecting
409 ------------------------------------------------------------------------- */
410
markAllContacts()411 void FixContactHistoryMesh::markAllContacts()
412 {
413 const int nlocal = atom->nlocal;
414
415 keeppage_[0]->reset(true);
416 intersectpage_[0]->reset(false);
417
418 for(int i = 0; i < nlocal; i++)
419 {
420 const int nneighs = fix_nneighs_->get_vector_atom_int(i);
421 keepflag_[i] = keeppage_[0]->get(nneighs);
422 intersectflag_[i] = intersectpage_[0]->get(nneighs);
423 if (!keepflag_[i] || !intersectflag_[i])
424 error->one(FLERR,"mesh contact history overflow, boost neigh_modify one");
425 }
426 }
427
428 /* ---------------------------------------------------------------------- */
429
cleanUpContacts()430 void FixContactHistoryMesh::cleanUpContacts()
431 {
432 cleanUpContacts(0, atom->nlocal);
433 }
434
435 /* ---------------------------------------------------------------------- */
436
cleanUpContacts(int ifrom,int ito)437 void FixContactHistoryMesh::cleanUpContacts(int ifrom, int ito)
438 {
439
440 for(int i = ifrom; i < ito; i++)
441 {
442 const int nneighs = fix_nneighs_->get_vector_atom_int(i);
443
444 for(int j = 0; j < nneighs; j++)
445 {
446 // delete values
447
448 if(!keepflag_[i][j])
449 {
450 if(partner_[i][j] > -1)
451 {
452 npartner_[i]--;
453
454 }
455 partner_[i][j] = -1;
456 intersectflag_[i][j] = false;
457 vectorZeroizeN(&(contacthistory_[i][j*dnum_]),dnum_);
458 }
459 }
460
461 }
462 }
463
464 /* ---------------------------------------------------------------------- */
465
cleanUpContactJumps()466 void FixContactHistoryMesh::cleanUpContactJumps()
467 {
468 const int nlocal = atom->nlocal;
469
470 for(int i = 0; i < nlocal; i++)
471 {
472 int ipartner = 0;
473 while (ipartner < npartner_[i])
474 {
475
476 if(partner_[i][ipartner] < 0)
477 {
478
479 error->one(FLERR,"internal error");
480 }
481
482 int iTri = -1;
483 const int nTri_id = mesh_->map_size(partner_[i][ipartner]);
484 for (int j = 0; j < nTri_id; j++)
485 {
486 iTri = mesh_->map(partner_[i][ipartner], j);
487 if (iTri != -1 && fix_neighlist_mesh_->contactInList(iTri,i))
488 break;
489 else
490 iTri = -1;
491 }
492
493 if(iTri == -1)
494 {
495
496 partner_[i][ipartner] = -1;
497 vectorZeroizeN(&(contacthistory_[i][ipartner*dnum_]),dnum_);
498 swap(i,ipartner,npartner_[i]-1,false);
499 npartner_[i]--;
500
501 }
502 else
503 ipartner++;
504 }
505 }
506 }
507
508 /* ---------------------------------------------------------------------- */
509
reset_history()510 void FixContactHistoryMesh::reset_history()
511 {
512 const int nlocal = atom->nlocal;
513
514 for(int i = 0; i < nlocal; i++)
515 {
516 const int nneighs = fix_nneighs_->get_vector_atom_int(i);
517
518 // zeroize values
519 for(int j = 0; j < nneighs; j++)
520 vectorZeroizeN(&(contacthistory_[i][j*dnum_]),dnum_);
521 }
522 }
523
524 /* ---------------------------------------------------------------------- */
525
min_setup_pre_force(int dummy)526 void FixContactHistoryMesh::min_setup_pre_force(int dummy)
527 {
528 if (*computeflag_) pre_force(0);
529 *computeflag_ = 0;
530 }
531
532 /* ---------------------------------------------------------------------- */
533
min_pre_force(int dummy)534 void FixContactHistoryMesh::min_pre_force(int dummy)
535 {
536 pre_force(0);
537 }
538
539 /* ----------------------------------------------------------------------
540 allocate local atom-based arrays
541 ------------------------------------------------------------------------- */
542
grow_arrays(int nmax)543 void FixContactHistoryMesh::grow_arrays(int nmax)
544 {
545 FixContactHistory::grow_arrays(nmax);
546 keepflag_ = (bool **) memory->srealloc(keepflag_,nmax*sizeof(bool *),
547 "contact_history:keepflag");
548 intersectflag_ = (bool **) memory->srealloc(intersectflag_,nmax*sizeof(bool *),
549 "contact_history:keepflag");
550 }
551
552 /* ----------------------------------------------------------------------
553 copy values within local atom-based arrays
554 ------------------------------------------------------------------------- */
555
copy_arrays(int i,int j,int delflag)556 void FixContactHistoryMesh::copy_arrays(int i, int j, int delflag)
557 {
558 // just copy pointers for partner and shearpartner
559 // b/c can't overwrite chunk allocation inside ipage,dpage
560 // incoming atoms in unpack_exchange just grab new chunks
561 // so are orphaning chunks for migrating atoms
562 // OK, b/c will reset ipage,dpage on next reneighboring
563
564 FixContactHistory::copy_arrays(i,j,delflag);
565 keepflag_[j] = keepflag_[i];
566 intersectflag_[j] = intersectflag_[i];
567 }
568
569 /* ----------------------------------------------------------------------
570 unpack values in local atom-based arrays from exchange with another proc
571 ------------------------------------------------------------------------- */
572
unpack_exchange(int nlocal,double * buf)573 int FixContactHistoryMesh::unpack_exchange(int nlocal, double *buf)
574 {
575 // allocate new chunks from ipage,dpage for incoming values
576
577 int m = 0;
578 const int nneighs = fix_nneighs_->get_vector_atom_int(nlocal);
579
580 npartner_[nlocal] = ubuf(buf[m++]).i;
581 maxtouch_ = MAX(maxtouch_,npartner_[nlocal]);
582 int nalloc = std::max(nneighs,npartner_[nlocal]);
583 partner_[nlocal] = ipage_->get(nalloc);
584 contacthistory_[nlocal] = dpage_->get(dnum_*nalloc);
585
586 if (!partner_[nlocal] || !contacthistory_[nlocal])
587 error->one(FLERR,"mesh contact history overflow, boost neigh_modify one");
588
589 for (int n = 0; n < npartner_[nlocal]; n++) {
590 partner_[nlocal][n] = ubuf(buf[m++]).i;
591
592 for (int d = 0; d < dnum_; d++) {
593 contacthistory_[nlocal][n*dnum_+d] = buf[m++];
594
595 }
596 }
597
598 for (int n = npartner_[nlocal]; n < nneighs; n++) {
599 partner_[nlocal][n] = -1;
600
601 for (int d = 0; d < dnum_; d++) {
602 contacthistory_[nlocal][n*dnum_+d] = 0.;
603 }
604 }
605
606 return m;
607 }
608
609 /* ----------------------------------------------------------------------
610 unpack values from atom->extra array to restart the fix
611 ------------------------------------------------------------------------- */
612
unpack_restart(int nlocal,int nth)613 void FixContactHistoryMesh::unpack_restart(int nlocal, int nth)
614 {
615 // ipage = NULL if being called from granular pair style init()
616
617 if (ipage_ == NULL) allocate_pages();
618
619 // skip to Nth set of extra values
620
621 double **extra = atom->extra;
622
623 int m = 0;
624 for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);
625 m++;
626
627 // allocate new chunks from ipage,dpage for incoming values
628
629 int d;
630
631 npartner_[nlocal] = ubuf(extra[nlocal][m++]).i;
632 maxtouch_ = MAX(maxtouch_,npartner_[nlocal]);
633 partner_[nlocal] = ipage_->get(npartner_[nlocal]);
634 contacthistory_[nlocal] = dpage_->get(npartner_[nlocal]*dnum_);
635
636 if (!partner_[nlocal] || !contacthistory_[nlocal])
637 error->one(FLERR,"mesh contact history overflow, boost neigh_modify one");
638
639 for (int n = 0; n < npartner_[nlocal]; n++) {
640 partner_[nlocal][n] = ubuf(extra[nlocal][m++]).i;
641
642 for (d = 0; d < dnum_; d++) {
643 contacthistory_[nlocal][n*dnum_+d] = extra[nlocal][m++];
644
645 }
646 }
647 }
648
649 /* ----------------------------------------------------------------------
650 pack state of Fix into one write
651 ------------------------------------------------------------------------- */
652
write_restart(FILE * fp)653 void FixContactHistoryMesh::write_restart(FILE *fp)
654 {
655 FixContactHistory::write_restart(fp);
656
657 sort_contacts();
658 }
659
660 /* ----------------------------------------------------------------------
661 memory usage of local atom-based arrays
662 ------------------------------------------------------------------------- */
663
memory_usage()664 double FixContactHistoryMesh::memory_usage()
665 {
666 int nmax = atom->nmax;
667 double bytes = nmax * sizeof(int);
668 bytes += nmax * sizeof(int *);
669 bytes += nmax * sizeof(double *);
670
671 for (int i = 0; i < numpages_; i++) {
672 bytes += ipage1_[i].size();
673 bytes += dpage1_[i].size();
674 bytes += ipage2_[i].size();
675 bytes += dpage2_[i].size();
676 bytes += keeppage_[i]->size();
677 bytes += intersectpage_[i]->size();
678 }
679
680 return bytes;
681 }
682
get_contact(const int i,const int j)683 int FixContactHistoryMesh::get_contact(const int i, const int j)
684 {
685 int count = -1;
686 int offset = 0;
687 for (; offset < nneighs(i); offset++)
688 {
689 // check whether real partner
690 if (partner_[i][offset] != -1)
691 {
692 count++;
693 if (j == count) // found j-th real contact partner
694 break;
695 }
696 }
697
698 if (j != count)
699 error->fix_error(FLERR, this, "could not find suitable partner");
700
701 return offset;
702 }
703