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
36 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
37 http://lammps.sandia.gov, Sandia National Laboratories
38 Steve Plimpton, sjplimp@sandia.gov
39
40 Copyright (2003) Sandia Corporation. Under the terms of Contract
41 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
42 certain rights in this software. This software is distributed under
43 the GNU General Public License.
44 ------------------------------------------------------------------------- */
45
46 #include "neighbor.h"
47 #include "neigh_list.h"
48 #include "atom.h"
49
50 using namespace LAMMPS_NS;
51
52 /* ----------------------------------------------------------------------
53 routines to create a stencil = list of bin offsets
54 stencil = bins whose closest corner to central bin is within cutoff
55 sx,sy,sz = bin bounds = furthest the stencil could possibly extend
56 3d creates xyz stencil, 2d creates xy stencil
57 for half list with newton off:
58 stencil is all surrounding bins including self
59 regardless of triclinic
60 for half list with newton on:
61 stencil is bins to the "upper right" of central bin
62 stencil does not include self
63 for half list with triclinic:
64 stencil is all bins in z-plane of self and above, but not below
65 in 2d is all bins in y-plane of self and above, but not below
66 stencil includes self
67 for full list:
68 stencil is all surrounding bins including self
69 regardless of newton on/off or triclinic
70 for multi:
71 create one stencil for each atom type
72 stencil is same bin bounds as newton on/off, triclinic, half/full
73 cutoff is not cutneighmaxsq, but max cutoff for that atom type
74 ------------------------------------------------------------------------- */
75
76 /* ---------------------------------------------------------------------- */
77
stencil_half_bin_2d_no_newton(NeighList * list,int sx,int sy,int sz)78 void Neighbor::stencil_half_bin_2d_no_newton(NeighList *list,
79 int sx, int sy, int sz)
80 {
81 int i,j;
82 int *stencil = list->stencil;
83 int nstencil = 0;
84
85 for (j = -sy; j <= sy; j++)
86 for (i = -sx; i <= sx; i++)
87 if (bin_distance(i,j,0) < cutneighmaxsq)
88 stencil[nstencil++] = j*mbinx + i;
89 list->nstencil = nstencil;
90 }
91
92 /* ---------------------------------------------------------------------- */
93
stencil_half_ghost_bin_2d_no_newton(NeighList * list,int sx,int sy,int sz)94 void Neighbor::stencil_half_ghost_bin_2d_no_newton(NeighList *list,
95 int sx, int sy, int sz)
96 {
97
98 int *stencil = list->stencil;
99 int **stencilxyz = list->stencilxyz;
100 int nstencil = 0;
101
102 for (int j = -sy; j <= sy; j++)
103 for (int i = -sx; i <= sx; i++)
104 if (bin_distance(i,j,0) < cutneighmaxsq) {
105 stencilxyz[nstencil][0] = i;
106 stencilxyz[nstencil][1] = j;
107 stencilxyz[nstencil][2] = 0;
108 stencil[nstencil++] = j*mbinx + i;
109 }
110
111 list->nstencil = nstencil;
112 }
113
114 /* ---------------------------------------------------------------------- */
115
stencil_half_bin_3d_no_newton(NeighList * list,int sx,int sy,int sz)116 void Neighbor::stencil_half_bin_3d_no_newton(NeighList *list,
117 int sx, int sy, int sz)
118 {
119 int i,j,k;
120 int *stencil = list->stencil;
121 int nstencil = 0;
122
123 for (k = -sz; k <= sz; k++)
124 for (j = -sy; j <= sy; j++)
125 for (i = -sx; i <= sx; i++)
126 if (bin_distance(i,j,k) < cutneighmaxsq)
127 stencil[nstencil++] = k*mbiny*mbinx + j*mbinx + i;
128
129 list->nstencil = nstencil;
130 }
131
132 /* ---------------------------------------------------------------------- */
133
stencil_half_ghost_bin_3d_no_newton(NeighList * list,int sx,int sy,int sz)134 void Neighbor::stencil_half_ghost_bin_3d_no_newton(NeighList *list,
135 int sx, int sy, int sz)
136 {
137 int i,j,k;
138 int *stencil = list->stencil;
139 int **stencilxyz = list->stencilxyz;
140 int nstencil = 0;
141
142 for (k = -sz; k <= sz; k++)
143 for (j = -sy; j <= sy; j++)
144 for (i = -sx; i <= sx; i++)
145 if (bin_distance(i,j,k) < cutneighmaxsq) {
146 stencilxyz[nstencil][0] = i;
147 stencilxyz[nstencil][1] = j;
148 stencilxyz[nstencil][2] = k;
149 stencil[nstencil++] = k*mbiny*mbinx + j*mbinx + i;
150 }
151
152 list->nstencil = nstencil;
153 }
154
155 /* ---------------------------------------------------------------------- */
156
stencil_half_bin_2d_newton(NeighList * list,int sx,int sy,int sz)157 void Neighbor::stencil_half_bin_2d_newton(NeighList *list,
158 int sx, int sy, int sz)
159 {
160 int i,j;
161 int *stencil = list->stencil;
162 int nstencil = 0;
163
164 for (j = 0; j <= sy; j++)
165 for (i = -sx; i <= sx; i++)
166 if (j > 0 || (j == 0 && i > 0))
167 if (bin_distance(i,j,0) < cutneighmaxsq)
168 stencil[nstencil++] = j*mbinx + i;
169
170 list->nstencil = nstencil;
171 }
172
173 /* ---------------------------------------------------------------------- */
174
stencil_half_bin_3d_newton(NeighList * list,int sx,int sy,int sz)175 void Neighbor::stencil_half_bin_3d_newton(NeighList *list,
176 int sx, int sy, int sz)
177 {
178 int i,j,k;
179 int *stencil = list->stencil;
180 int nstencil = 0;
181
182 for (k = 0; k <= sz; k++)
183 for (j = -sy; j <= sy; j++)
184 for (i = -sx; i <= sx; i++)
185 if (k > 0 || j > 0 || (j == 0 && i > 0))
186 if (bin_distance(i,j,k) < cutneighmaxsq)
187 stencil[nstencil++] = k*mbiny*mbinx + j*mbinx + i;
188
189 list->nstencil = nstencil;
190 }
191
192 /* ---------------------------------------------------------------------- */
193
stencil_half_bin_2d_newton_tri(NeighList * list,int sx,int sy,int sz)194 void Neighbor::stencil_half_bin_2d_newton_tri(NeighList *list,
195 int sx, int sy, int sz)
196 {
197 int i,j;
198 int *stencil = list->stencil;
199 int nstencil = 0;
200
201 for (j = 0; j <= sy; j++)
202 for (i = -sx; i <= sx; i++)
203 if (bin_distance(i,j,0) < cutneighmaxsq)
204 stencil[nstencil++] = j*mbinx + i;
205
206 list->nstencil = nstencil;
207 }
208
209 /* ---------------------------------------------------------------------- */
210
stencil_half_bin_3d_newton_tri(NeighList * list,int sx,int sy,int sz)211 void Neighbor::stencil_half_bin_3d_newton_tri(NeighList *list,
212 int sx, int sy, int sz)
213 {
214 int i,j,k;
215 int *stencil = list->stencil;
216 int nstencil = 0;
217
218 for (k = 0; k <= sz; k++)
219 for (j = -sy; j <= sy; j++)
220 for (i = -sx; i <= sx; i++)
221 if (bin_distance(i,j,k) < cutneighmaxsq)
222 stencil[nstencil++] = k*mbiny*mbinx + j*mbinx + i;
223
224 list->nstencil = nstencil;
225 }
226
227 /* ---------------------------------------------------------------------- */
228
stencil_half_multi_2d_no_newton(NeighList * list,int sx,int sy,int sz)229 void Neighbor::stencil_half_multi_2d_no_newton(NeighList *list,
230 int sx, int sy, int sz)
231 {
232 int i,j,n;
233 double rsq,typesq;
234 int *s;
235 double *distsq;
236
237 int *nstencil_multi = list->nstencil_multi;
238 int **stencil_multi = list->stencil_multi;
239 double **distsq_multi = list->distsq_multi;
240
241 int ntypes = atom->ntypes;
242 for (int itype = 1; itype <= ntypes; itype++) {
243 typesq = cuttypesq[itype];
244 s = stencil_multi[itype];
245 distsq = distsq_multi[itype];
246 n = 0;
247 for (j = -sy; j <= sy; j++)
248 for (i = -sx; i <= sx; i++) {
249 rsq = bin_distance(i,j,0);
250 if (rsq < typesq) {
251 distsq[n] = rsq;
252 s[n++] = j*mbinx + i;
253 }
254 }
255 nstencil_multi[itype] = n;
256 }
257 }
258
259 /* ---------------------------------------------------------------------- */
260
stencil_half_multi_3d_no_newton(NeighList * list,int sx,int sy,int sz)261 void Neighbor::stencil_half_multi_3d_no_newton(NeighList *list,
262 int sx, int sy, int sz)
263 {
264 int i,j,k,n;
265 double rsq,typesq;
266 int *s;
267 double *distsq;
268
269 int *nstencil_multi = list->nstencil_multi;
270 int **stencil_multi = list->stencil_multi;
271 double **distsq_multi = list->distsq_multi;
272
273 int ntypes = atom->ntypes;
274 for (int itype = 1; itype <= ntypes; itype++) {
275 typesq = cuttypesq[itype];
276 s = stencil_multi[itype];
277 distsq = distsq_multi[itype];
278 n = 0;
279 for (k = -sz; k <= sz; k++)
280 for (j = -sy; j <= sy; j++)
281 for (i = -sx; i <= sx; i++) {
282 rsq = bin_distance(i,j,k);
283 if (rsq < typesq) {
284 distsq[n] = rsq;
285 s[n++] = k*mbiny*mbinx + j*mbinx + i;
286 }
287 }
288 nstencil_multi[itype] = n;
289 }
290 }
291
292 /* ---------------------------------------------------------------------- */
293
stencil_half_multi_2d_newton(NeighList * list,int sx,int sy,int sz)294 void Neighbor::stencil_half_multi_2d_newton(NeighList *list,
295 int sx, int sy, int sz)
296 {
297 int i,j,n;
298 double rsq,typesq;
299 int *s;
300 double *distsq;
301
302 int *nstencil_multi = list->nstencil_multi;
303 int **stencil_multi = list->stencil_multi;
304 double **distsq_multi = list->distsq_multi;
305
306 int ntypes = atom->ntypes;
307 for (int itype = 1; itype <= ntypes; itype++) {
308 typesq = cuttypesq[itype];
309 s = stencil_multi[itype];
310 distsq = distsq_multi[itype];
311 n = 0;
312 for (j = 0; j <= sy; j++)
313 for (i = -sx; i <= sx; i++)
314 if (j > 0 || (j == 0 && i > 0)) {
315 rsq = bin_distance(i,j,0);
316 if (rsq < typesq) {
317 distsq[n] = rsq;
318 s[n++] = j*mbinx + i;
319 }
320 }
321 nstencil_multi[itype] = n;
322 }
323 }
324
325 /* ---------------------------------------------------------------------- */
326
stencil_half_multi_3d_newton(NeighList * list,int sx,int sy,int sz)327 void Neighbor::stencil_half_multi_3d_newton(NeighList *list,
328 int sx, int sy, int sz)
329 {
330 int i,j,k,n;
331 double rsq,typesq;
332 int *s;
333 double *distsq;
334
335 int *nstencil_multi = list->nstencil_multi;
336 int **stencil_multi = list->stencil_multi;
337 double **distsq_multi = list->distsq_multi;
338
339 int ntypes = atom->ntypes;
340 for (int itype = 1; itype <= ntypes; itype++) {
341 typesq = cuttypesq[itype];
342 s = stencil_multi[itype];
343 distsq = distsq_multi[itype];
344 n = 0;
345 for (k = 0; k <= sz; k++)
346 for (j = -sy; j <= sy; j++)
347 for (i = -sx; i <= sx; i++)
348 if (k > 0 || j > 0 || (j == 0 && i > 0)) {
349 rsq = bin_distance(i,j,k);
350 if (rsq < typesq) {
351 distsq[n] = rsq;
352 s[n++] = k*mbiny*mbinx + j*mbinx + i;
353 }
354 }
355 nstencil_multi[itype] = n;
356 }
357 }
358
359 /* ---------------------------------------------------------------------- */
360
stencil_half_multi_2d_newton_tri(NeighList * list,int sx,int sy,int sz)361 void Neighbor::stencil_half_multi_2d_newton_tri(NeighList *list,
362 int sx, int sy, int sz)
363 {
364 int i,j,n;
365 double rsq,typesq;
366 int *s;
367 double *distsq;
368
369 int *nstencil_multi = list->nstencil_multi;
370 int **stencil_multi = list->stencil_multi;
371 double **distsq_multi = list->distsq_multi;
372
373 int ntypes = atom->ntypes;
374 for (int itype = 1; itype <= ntypes; itype++) {
375 typesq = cuttypesq[itype];
376 s = stencil_multi[itype];
377 distsq = distsq_multi[itype];
378 n = 0;
379 for (j = 0; j <= sy; j++)
380 for (i = -sx; i <= sx; i++) {
381 rsq = bin_distance(i,j,0);
382 if (rsq < typesq) {
383 distsq[n] = rsq;
384 s[n++] = j*mbinx + i;
385 }
386 }
387 nstencil_multi[itype] = n;
388 }
389 }
390
391 /* ---------------------------------------------------------------------- */
392
stencil_half_multi_3d_newton_tri(NeighList * list,int sx,int sy,int sz)393 void Neighbor::stencil_half_multi_3d_newton_tri(NeighList *list,
394 int sx, int sy, int sz)
395 {
396 int i,j,k,n;
397 double rsq,typesq;
398 int *s;
399 double *distsq;
400
401 int *nstencil_multi = list->nstencil_multi;
402 int **stencil_multi = list->stencil_multi;
403 double **distsq_multi = list->distsq_multi;
404
405 int ntypes = atom->ntypes;
406 for (int itype = 1; itype <= ntypes; itype++) {
407 typesq = cuttypesq[itype];
408 s = stencil_multi[itype];
409 distsq = distsq_multi[itype];
410 n = 0;
411 for (k = 0; k <= sz; k++)
412 for (j = -sy; j <= sy; j++)
413 for (i = -sx; i <= sx; i++) {
414 rsq = bin_distance(i,j,k);
415 if (rsq < typesq) {
416 distsq[n] = rsq;
417 s[n++] = k*mbiny*mbinx + j*mbinx + i;
418 }
419 }
420 nstencil_multi[itype] = n;
421 }
422 }
423
424 /* ---------------------------------------------------------------------- */
425
stencil_full_bin_2d(NeighList * list,int sx,int sy,int sz)426 void Neighbor::stencil_full_bin_2d(NeighList *list,
427 int sx, int sy, int sz)
428 {
429 int i,j;
430 int *stencil = list->stencil;
431 int nstencil = 0;
432
433 for (j = -sy; j <= sy; j++)
434 for (i = -sx; i <= sx; i++)
435 if (bin_distance(i,j,0) < cutneighmaxsq)
436 stencil[nstencil++] = j*mbinx + i;
437
438 list->nstencil = nstencil;
439 }
440
441 /* ---------------------------------------------------------------------- */
442
stencil_full_ghost_bin_2d(NeighList * list,int sx,int sy,int sz)443 void Neighbor::stencil_full_ghost_bin_2d(NeighList *list,
444 int sx, int sy, int sz)
445 {
446 int i,j;
447 int *stencil = list->stencil;
448 int **stencilxyz = list->stencilxyz;
449 int nstencil = 0;
450
451 for (j = -sy; j <= sy; j++)
452 for (i = -sx; i <= sx; i++)
453 if (bin_distance(i,j,0) < cutneighmaxsq) {
454 stencilxyz[nstencil][0] = i;
455 stencilxyz[nstencil][1] = j;
456 stencilxyz[nstencil][2] = 0;
457 stencil[nstencil++] = j*mbinx + i;
458 }
459
460 list->nstencil = nstencil;
461 }
462
463 /* ---------------------------------------------------------------------- */
464
stencil_full_bin_3d(NeighList * list,int sx,int sy,int sz)465 void Neighbor::stencil_full_bin_3d(NeighList *list,
466 int sx, int sy, int sz)
467 {
468 int i,j,k;
469 int *stencil = list->stencil;
470 int nstencil = 0;
471
472 for (k = -sz; k <= sz; k++)
473 for (j = -sy; j <= sy; j++)
474 for (i = -sx; i <= sx; i++)
475 if (bin_distance(i,j,k) < cutneighmaxsq)
476 stencil[nstencil++] = k*mbiny*mbinx + j*mbinx + i;
477
478 list->nstencil = nstencil;
479 }
480
481 /* ---------------------------------------------------------------------- */
482
stencil_full_ghost_bin_3d(NeighList * list,int sx,int sy,int sz)483 void Neighbor::stencil_full_ghost_bin_3d(NeighList *list,
484 int sx, int sy, int sz)
485 {
486 int i,j,k;
487 int *stencil = list->stencil;
488 int **stencilxyz = list->stencilxyz;
489 int nstencil = 0;
490
491 for (k = -sz; k <= sz; k++)
492 for (j = -sy; j <= sy; j++)
493 for (i = -sx; i <= sx; i++)
494 if (bin_distance(i,j,k) < cutneighmaxsq) {
495 stencilxyz[nstencil][0] = i;
496 stencilxyz[nstencil][1] = j;
497 stencilxyz[nstencil][2] = k;
498 stencil[nstencil++] = k*mbiny*mbinx + j*mbinx + i;
499 }
500
501 list->nstencil = nstencil;
502 }
503
504 /* ---------------------------------------------------------------------- */
505
stencil_full_multi_2d(NeighList * list,int sx,int sy,int sz)506 void Neighbor::stencil_full_multi_2d(NeighList *list,
507 int sx, int sy, int sz)
508 {
509 int i,j,n;
510 double rsq,typesq;
511 int *s;
512 double *distsq;
513
514 int *nstencil_multi = list->nstencil_multi;
515 int **stencil_multi = list->stencil_multi;
516 double **distsq_multi = list->distsq_multi;
517
518 int ntypes = atom->ntypes;
519 for (int itype = 1; itype <= ntypes; itype++) {
520 typesq = cuttypesq[itype];
521 s = stencil_multi[itype];
522 distsq = distsq_multi[itype];
523 n = 0;
524 for (j = -sy; j <= sy; j++)
525 for (i = -sx; i <= sx; i++) {
526 rsq = bin_distance(i,j,0);
527 if (rsq < typesq) {
528 distsq[n] = rsq;
529 s[n++] = j*mbinx + i;
530 }
531 }
532 nstencil_multi[itype] = n;
533 }
534 }
535
536 /* ---------------------------------------------------------------------- */
537
stencil_full_multi_3d(NeighList * list,int sx,int sy,int sz)538 void Neighbor::stencil_full_multi_3d(NeighList *list,
539 int sx, int sy, int sz)
540 {
541 int i,j,k,n;
542 double rsq,typesq;
543 int *s;
544 double *distsq;
545
546 int *nstencil_multi = list->nstencil_multi;
547 int **stencil_multi = list->stencil_multi;
548 double **distsq_multi = list->distsq_multi;
549
550 int ntypes = atom->ntypes;
551 for (int itype = 1; itype <= ntypes; itype++) {
552 typesq = cuttypesq[itype];
553 s = stencil_multi[itype];
554 distsq = distsq_multi[itype];
555 n = 0;
556 for (k = -sz; k <= sz; k++)
557 for (j = -sy; j <= sy; j++)
558 for (i = -sx; i <= sx; i++) {
559 rsq = bin_distance(i,j,k);
560 if (rsq < typesq) {
561 distsq[n] = rsq;
562 s[n++] = k*mbiny*mbinx + j*mbinx + i;
563 }
564 }
565 nstencil_multi[itype] = n;
566 }
567 }
568