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 Richard Berger (JKU Linz)
39 Christoph Kloss (DCS Computing GmbH)
40 Alexander Podlozhnyuk (DCS Computing GmbH)
41 Josef Kerbl (DCS Computing GmbH)
42
43 Copyright 2014-2015 JKU Linz
44 Copyright 2015- DCS Computing GmbH
45 ------------------------------------------------------------------------- */
46
47 // include last to ensure correct macros
48 #include "domain_definitions.h"
49
50 static const double SMALL_REGION_NEIGHBOR_LIST = 1.0e-6;
51 static const double BIG_REGION_NEIGHBOR_LIST = 1.0e20;
52
53 /**
54 * @brief Default constructor which will create an empty neighbor list
55 */
56 template<bool INTERPOLATE>
RegionNeighborList(LAMMPS * lmp)57 RegionNeighborList<INTERPOLATE>::RegionNeighborList(LAMMPS *lmp) :
58 IRegionNeighborList(),
59 Pointers(lmp),
60 ncount(0),
61 bbox_set(false)
62 {
63 }
64
65 /**
66 * @brief Determine if the given particle overlaps with any particle in this neighbor list
67 * @param x position of particle to check
68 * @param radius radius of particle to check
69 * @return true if particle has an overlap with a particle in this neighbor list, false otherwise
70 */
71 template<bool INTERPOLATE>
hasOverlap(double * x,double radius)72 bool RegionNeighborList<INTERPOLATE>::hasOverlap(double * x, double radius) const
73 {
74 int ibin = coord2bin(x);
75
76 for(std::vector<int>::const_iterator it = stencil.begin(); it != stencil.end(); ++it)
77 {
78 const int offset = *it;
79 if((ibin+offset < 0) || ((size_t)(ibin+offset) >= bins.size()))
80 {
81
82 error->one(FLERR,"assertion failed");
83 }
84 const std::vector<Particle<INTERPOLATE> > & plist = bins[ibin+offset].particles;
85
86 for(typename std::vector<Particle<INTERPOLATE> >::const_iterator pit = plist.begin(); pit != plist.end(); ++pit)
87 {
88 const Particle<INTERPOLATE> & p = *pit;
89 double del[3];
90 vectorSubtract3D(x, p.x, del);
91 const double rsq = vectorMag3DSquared(del);
92 const double radsum = radius + p.radius;
93 if (rsq <= radsum*radsum) return true;
94 }
95 }
96
97 return false;
98 }
99
100 #ifdef SUPERQUADRIC_ACTIVE_FLAG
101 template<bool INTERPOLATE>
hasOverlap_superquadric(double * x,double radius,double * quaternion,double * shape,double * blockiness)102 bool RegionNeighborList<INTERPOLATE>::hasOverlap_superquadric(double * x, double radius, double *quaternion, double *shape, double *blockiness) const
103 {
104 int ibin = coord2bin(x);
105
106 for(std::vector<int>::const_iterator it = stencil.begin(); it != stencil.end(); ++it) {
107 const int offset = *it;
108 if((ibin+offset < 0) || ((size_t)(ibin+offset) >= bins.size()))
109 {
110
111 error->one(FLERR,"assertion failed");
112 }
113 const std::vector<Particle<INTERPOLATE> > & plist = bins[ibin+offset].particles;
114
115 Superquadric particle1(x, quaternion, shape, blockiness);
116 for(typename std::vector<Particle<INTERPOLATE> >::const_iterator pit = plist.begin(); pit != plist.end(); ++pit) {
117 const Particle<INTERPOLATE> & p = *pit;
118 double del[3];
119 vectorSubtract3D(x, p.x, del);
120 const double rsq = vectorMag3DSquared(del);
121 const double radsum = radius + p.radius;
122 if(check_obb_flag) {
123 double x_copy[3], quaternion_copy[4], shape_copy[3], blockiness_copy[2];
124 vectorCopy3D(p.x, x_copy);
125 vectorCopy4D(p.quaternion, quaternion_copy);
126 vectorCopy3D(p.shape, shape_copy);
127 vectorCopy2D(p.blockiness, blockiness_copy);
128 Superquadric particle2(x_copy, quaternion_copy, shape_copy, blockiness_copy);
129 double contact_point[3];
130
131 if (rsq <= radsum*radsum and (MathExtraLiggghtsNonspherical::capsules_intersect(&particle1, &particle2, contact_point) and
132 MathExtraLiggghtsNonspherical::obb_intersect(&particle1, &particle2))) {
133 return true;
134 }
135 } else
136 if (rsq <= radsum*radsum) return true;
137 }
138 }
139
140 return false;
141 }
142
143 #endif
144 /**
145 * @brief Determine if the given particle overlaps with any particle in this neighbor list
146 * @param x position of particle to check
147 * @param radius radius of particle to check
148 * @param overlap_list list of overlaps, to be populated by this function
149 * @return true if particle has an overlap with a particle in this neighbor list, false otherwise
150 */
151
152 template<bool INTERPOLATE>
hasOverlapWith(double * x,double radius,std::vector<int> & overlap_list)153 bool RegionNeighborList<INTERPOLATE>::hasOverlapWith(double * x, double radius, std::vector<int> &overlap_list) const
154 {
155 int ibin = coord2bin(x);
156
157 bool overlap = false;
158
159 for(std::vector<int>::const_iterator it = stencil.begin(); it != stencil.end(); ++it)
160 {
161 const int offset = *it;
162 if((ibin+offset < 0) || ((size_t)(ibin+offset) >= bins.size()))
163 {
164
165 error->one(FLERR,"assertion failed");
166 }
167 const std::vector<Particle<INTERPOLATE> > & plist = bins[ibin+offset].particles;
168
169 for(typename std::vector<Particle<INTERPOLATE> >::const_iterator pit = plist.begin(); pit != plist.end(); ++pit)
170 {
171 const Particle<INTERPOLATE> & p = *pit;
172 double del[3];
173 vectorSubtract3D(x, p.x, del);
174 const double rsq = vectorMag3DSquared(del);
175 const double radsum = radius + p.radius;
176 if (rsq <= radsum*radsum)
177 {
178 overlap_list.push_back(p.index);
179 overlap = true;
180 }
181 }
182 }
183
184 return overlap;
185 }
186
187 /**
188 * @brief Insert a new particle into neighbor list
189 * @param x position in 3D
190 * @param radius particle radius
191 */
192 template<bool INTERPOLATE>
insert(double * x,double radius,int index)193 void RegionNeighborList<INTERPOLATE>::insert(double * x, double radius,int index)
194 {
195 int quadrant;
196 double wx,wy,wz;
197 int ibin = coord2bin(x,quadrant,wx,wy,wz);
198 if((ibin < 0) || ((size_t)(ibin) >= bins.size()))
199 {
200
201 error->one(FLERR,"assertion failed");
202 }
203
204 bins[ibin].particles.push_back(Particle<INTERPOLATE>(index,x, radius,ibin,quadrant,wx,wy,wz));
205
206 ++ncount;
207 }
208
209 #ifdef SUPERQUADRIC_ACTIVE_FLAG
210 template<bool INTERPOLATE>
insert_superquadric(double * x,double radius,double * quaternion,double * shape,double * blockiness,int index)211 void RegionNeighborList<INTERPOLATE>::insert_superquadric(double * x, double radius, double *quaternion, double *shape, double *blockiness, int index) {
212 int quadrant;
213 double wx,wy,wz;
214 int ibin = coord2bin(x,quadrant,wx,wy,wz);
215 if((ibin < 0) || ((size_t)(ibin) >= bins.size()))
216 {
217
218 error->one(FLERR,"assertion failed");
219 }
220
221 bins[ibin].particles.push_back(Particle<INTERPOLATE>(index,x,radius,quaternion, shape, blockiness, ibin,quadrant,wx,wy,wz));
222 ++ncount;
223 }
224 #endif
225
226 /**
227 * @brief Reset neighbor list and brings it into initial state
228 */
229 template<bool INTERPOLATE>
reset()230 void RegionNeighborList<INTERPOLATE>::reset()
231 {
232 // reset all the settings from setBoundaryBox
233 bins.clear();
234 stencil.clear();
235 bbox_set = false;
236 nbinx = nbiny = nbinz = 0;
237 binsizex = binsizey = binsizez = 0;
238 bininvx = bininvy = bininvz = 0;
239 mbinxlo = mbinylo = mbinzlo = 0;
240 mbinx = mbiny = mbinz = 0;
241
242 // reset particle counter
243 ncount = 0;
244 }
245
246 /**
247 * @brief Return size of buffer for one bin
248 *
249 * \sa pushBinToBuffer
250 */
251 template<bool INTERPOLATE>
getSizeOne()252 int RegionNeighborList<INTERPOLATE>::getSizeOne() const
253 {
254 return 4;
255 }
256
257 /**
258 * @brief Push bin information to buffer
259 *
260 * \sa getSizeOne
261 */
262 template<bool INTERPOLATE>
pushBinToBuffer(int i,double * buf)263 int RegionNeighborList<INTERPOLATE>::pushBinToBuffer(int i, double *buf) const
264 {
265 int m = 0;
266 buf[m++] = this->bins[i].center[0];
267 buf[m++] = this->bins[i].center[1];
268 buf[m++] = this->bins[i].center[2];
269 buf[m++] = this->bins[i].id;
270 return m;
271 }
272
273 /**
274 * @brief Clear bins (remove all inserted particles)
275 */
276 template<bool INTERPOLATE>
clear()277 void RegionNeighborList<INTERPOLATE>::clear()
278 {
279 for(typename std::vector<Bin<INTERPOLATE> >::iterator it = bins.begin(); it != bins.end(); ++it)
280 {
281 (*it).particles.clear();
282 }
283 ncount = 0;
284 }
285
286 /**
287 * @brief Returns the number of particles inserted into the neighbor list
288 * @return number of particles in neighbor list
289 */
290
291 template<bool INTERPOLATE>
count()292 size_t RegionNeighborList<INTERPOLATE>::count() const
293 {
294 return ncount;
295 }
296
297 /**
298 * @brief Set or Update the region bounding box
299 *
300 * This will update internal data structures to ensure they can handle the new
301 * region, which is defined by its bounding box.
302 *
303 * @param bb bounding box of region
304 * @param maxrad largest particle radius
305 * @param extend enable extending the box for ghost particles
306 * @param failsafe limit max. number of bins
307 * @return true if bounding box was set successfully, false bounding box could not
308 * be set and neighbor list is not usable
309 */
310
311 template<bool INTERPOLATE>
setBoundingBox_calc_interpolation_stencil(Bin<INTERPOLATE> & it,int ibin,int ix,int iy,int iz)312 inline void RegionNeighborList<INTERPOLATE>::setBoundingBox_calc_interpolation_stencil(Bin<INTERPOLATE> &it,int ibin,int ix,int iy, int iz) const
313 {
314
315 }
316
317 template<>
setBoundingBox_calc_interpolation_stencil(Bin<true> & it,int ibin,int ix,int iy,int iz)318 inline void RegionNeighborList<true>::setBoundingBox_calc_interpolation_stencil(Bin<true> &it,int ibin,int ix,int iy, int iz) const
319 {
320 Stencil stencilTmp;
321
322 it.stencils.clear();
323
324 int stencil_shift_up_quadrant = 0;
325 int stencil_shift_down_quadrant = 0;
326 int ii,jj,kk;
327 for(int iquad = 0; iquad < 8; iquad++)
328 {
329
330 stencil_shift_up_quadrant = 0;
331 stencil_shift_down_quadrant = 0;
332 const int qx = (iquad & 1) >> 0; // 0 or 1
333 const int qy = (iquad & 2) >> 1; // 0 or 1
334 const int qz = (iquad & 4) >> 2; // 0 or 1
335
336 if(0 == ix && 0 == qx)
337 {
338 stencil_shift_up_quadrant |= 1;
339 ii = 1;
340 }
341 else if(ix == (mbinx-1) && 1 == qx)
342 {
343 stencil_shift_down_quadrant |= 1;
344 ii = -1;
345 }
346 else
347 ii = 0;
348
349 if(0 == iy && 0 == qy)
350 {
351 stencil_shift_up_quadrant |= 2;
352 jj = 1;
353 }
354 else if(iy == (mbiny-1) && 1 == qy)
355 {
356 stencil_shift_down_quadrant |= 2;
357 jj = -1;
358 }
359 else
360 jj = 0;
361
362 if(0 == iz && 0 == qz)
363 {
364 stencil_shift_up_quadrant |= 4;
365 kk = 1;
366 }
367 else if(iz == (mbinz-1) && 1 == qz)
368 {
369 stencil_shift_down_quadrant |= 4;
370 kk = -1;
371 }
372 else
373 kk = 0;
374
375 // generate stencil which will look at 8 relevant bins - self + 7 quadrant neighs
376 // self can be located anywhere within these 8 bins
377 stencilTmp.clear();
378 for (int k = -1+qz; k <= qz; k++)
379 {
380 for (int j = -1+qy; j <= qy; j++)
381 {
382 for (int i = -1+qx; i <= qx; i++)
383 {
384 int isubbin = ((iz+k+kk)*mbiny*mbinx + (iy+j+jj)*mbinx + (ix+i+ii));
385 if( isubbin < 0 || isubbin >= mbins())
386 stencilTmp.push_back(-1);
387 else
388 stencilTmp.push_back(isubbin);
389
390 }
391 }
392 }
393
394 it.stencils.push_back(stencilTmp);
395
396 it.stencil_shift_up.push_back(stencil_shift_up_quadrant);
397 it.stencil_shift_down.push_back(stencil_shift_down_quadrant);
398
399 }
400 }
401
402 template<bool INTERPOLATE>
setBoundingBox(BoundingBox & bb,double maxrad,bool extend,bool failsafe)403 bool RegionNeighborList<INTERPOLATE>::setBoundingBox(BoundingBox & bb, double maxrad, bool extend, bool failsafe)
404 {
405 double extent[3];
406 bb.getExtent(extent);
407
408 if(extent[0] <= 0.0 || extent[1] <= 0.0 || extent[2] <= 0.0) {
409 // empty or invalid region
410 bins.clear();
411 stencil.clear();
412 return false;
413 }
414
415 bb.getBoxBounds(bboxlo, bboxhi);
416
417 // direct comparison of doubles to detect INF
418 if ( bboxlo[0] == -BIG || bboxlo[1] == -BIG || bboxlo[2] == -BIG ||
419 bboxhi[0] == BIG || bboxhi[1] == BIG || bboxhi[2] == BIG )
420 {
421 error->one(FLERR,"'INF' not allowed for definiton of region used for RegionNeighborList.\n"
422 "You may want to use 'EDGE' instead.");
423 }
424
425 // testing code
426 double binsize_optimal = 4*maxrad;
427 double binsizeinv = 1.0/binsize_optimal;
428
429 // test for too many global bins in any dimension due to huge global domain or small maxrad
430 const int max_bins = 8000000;
431 // we may repeat this calculation in case of failsafe
432 bool repeat = false;
433 bigint bbin = -1; // final number of bins
434 double bsubboxlo[3], bsubboxhi[3]; // box bounds
435 do
436 {
437 // create actual bins
438 nbinx = static_cast<int>((extent[0]+binsize_optimal*1e-10)*binsizeinv);
439 nbiny = static_cast<int>((extent[1]+binsize_optimal*1e-10)*binsizeinv);
440 nbinz = static_cast<int>((extent[2]+binsize_optimal*1e-10)*binsizeinv);
441
442 if (nbinx == 0) nbinx = 1;
443 if (nbiny == 0) nbiny = 1;
444 if (nbinz == 0) nbinz = 1;
445
446 binsizex = extent[0]/nbinx;
447 binsizey = extent[1]/nbiny;
448 binsizez = extent[2]/nbinz;
449
450 bininvx = 1.0 / binsizex;
451 bininvy = 1.0 / binsizey;
452 bininvz = 1.0 / binsizez;
453
454 // mbinlo/hi = lowest and highest global bins my ghost atoms could be in
455 // coord = lowest and highest values of coords for my ghost atoms
456 // static_cast(-1.5) = -1, so subract additional -1
457 // add in SMALL for round-off safety
458 bb.getBoxBounds(bsubboxlo, bsubboxhi);
459
460 // list is extended for ghost atoms or
461 // the list covers just exactly the region
462 if (extend)
463 {
464 double coord = bsubboxlo[0] - SMALL_REGION_NEIGHBOR_LIST*extent[0];
465 mbinxlo = static_cast<int> ((coord-bboxlo[0])*bininvx);
466 if (coord < bboxlo[0]) mbinxlo = mbinxlo - 1;
467 coord = bsubboxhi[0] + SMALL_REGION_NEIGHBOR_LIST*extent[0];
468 int mbinxhi = static_cast<int> ((coord-bboxlo[0])*bininvx);
469
470 coord = bsubboxlo[1] - SMALL_REGION_NEIGHBOR_LIST*extent[1];
471 mbinylo = static_cast<int> ((coord-bboxlo[1])*bininvy);
472 if (coord < bboxlo[1]) mbinylo = mbinylo - 1;
473 coord = bsubboxhi[1] + SMALL_REGION_NEIGHBOR_LIST*extent[1];
474 int mbinyhi = static_cast<int> ((coord-bboxlo[1])*bininvy);
475
476 coord = bsubboxlo[2] - SMALL_REGION_NEIGHBOR_LIST*extent[2];
477 mbinzlo = static_cast<int> ((coord-bboxlo[2])*bininvz);
478 if (coord < bboxlo[2]) mbinzlo = mbinzlo - 1;
479 coord = bsubboxhi[2] + SMALL_REGION_NEIGHBOR_LIST*extent[2];
480 int mbinzhi = static_cast<int> ((coord-bboxlo[2])*bininvz);
481
482 // extend bins by 1 to insure stencil extent is included
483
484 mbinxlo = mbinxlo - 1;
485 mbinxhi = mbinxhi + 1;
486 mbinylo = mbinylo - 1;
487 mbinyhi = mbinyhi + 1;
488 mbinzlo = mbinzlo - 1;
489 mbinzhi = mbinzhi + 1;
490
491 mbinx = mbinxhi - mbinxlo + 1;
492 mbiny = mbinyhi - mbinylo + 1;
493 mbinz = mbinzhi - mbinzlo + 1;
494
495 }
496 else
497 {
498 mbinxlo = mbinylo = mbinzlo = 0;
499 mbinx = nbinx;
500 mbiny = nbiny;
501 mbinz = nbinz;
502
503 }
504
505 // final check of number of bins
506 bbin = ((bigint) mbinx) * ((bigint) mbiny) * ((bigint) mbinz);
507
508 if (bbin > max_bins) { // too many bins
509 // check for failsafe mode
510 // check for repeat - ensure that we run the loop only once!
511 if(failsafe && !repeat)
512 {
513 binsizeinv = 1./ (vectorMax3D(extent) / 100.);
514 repeat = true;
515 } else {
516 printf("ERROR: Too many neighbor bins\n");
517 return false;
518 }
519 }
520 else
521 repeat = false;
522 } while (repeat);
523
524 // allocate bins
525 bins.resize(bbin);
526
527 // set cell center and stencils for each bin
528 int cId = 0;
529 for(typename std::vector<Bin<INTERPOLATE> >::iterator it = bins.begin(); it != bins.end(); ++it)
530 {
531 const int ibin = it - bins.begin();
532 int itmp = ibin;
533 const int iz = ibin/(mbinx*mbiny);
534 itmp -= iz*mbinx*mbiny;
535 const int iy = itmp/mbinx;
536 itmp -= iy*mbinx;
537 const int ix = itmp;
538 it->id = cId++;
539 it->center[0] = bsubboxlo[0] + binsizex*((double)(ix+mbinxlo)+0.5);
540 it->center[1] = bsubboxlo[1] + binsizey*((double)(iy+mbinylo)+0.5);
541 it->center[2] = bsubboxlo[2] + binsizez*((double)(iz+mbinzlo)+0.5);
542
543 setBoundingBox_calc_interpolation_stencil(*it,ibin,ix,iy,iz);
544
545 //#ifdef LIGGGHTS_DEBUG
546 // printf("Bin (%d) with indizes: %d %d %d\n",ibin,ix,iy,iz);
547 // printf("Center of bin (%d): %g %g %g\n",ibin,it->center[0],it->center[1],it->center[2]);
548 //#endif
549 }
550
551 // generate stencil which will look at all bins 27 bins
552 for (int k = -1; k <= 1; k++)
553 for (int j = -1; j <= 1; j++)
554 for (int i = -1; i <= 1; i++)
555 stencil.push_back(k*mbiny*mbinx + j*mbinx + i);
556
557 bbox_set = true;
558
559 return true;
560 }
561
562 /**
563 * @brief Compute the closest distance between central bin (0,0,0) and bin (i,j,k)
564 * @param i bin coordinate along x-axis
565 * @param j bin coordinate along y-axis
566 * @param k bin coordinate along z-axis
567 * @return closest distance between central bin (0,0,0) and bin (i,j,k)
568 */
569 template<bool INTERPOLATE>
bin_distance(int i,int j,int k)570 double RegionNeighborList<INTERPOLATE>::bin_distance(int i, int j, int k)
571 {
572 double delx,dely,delz;
573
574 if (i > 0) delx = (i-1)*binsizex;
575 else if (i == 0) delx = 0.0;
576 else delx = (i+1)*binsizex;
577
578 if (j > 0) dely = (j-1)*binsizey;
579 else if (j == 0) dely = 0.0;
580 else dely = (j+1)*binsizey;
581
582 if (k > 0) delz = (k-1)*binsizez;
583 else if (k == 0) delz = 0.0;
584 else delz = (k+1)*binsizez;
585
586 return (delx*delx + dely*dely + delz*delz);
587 }
588
589 /**
590 * @brief Calc local bin index (m) of point x
591 * @param x point in 3D
592 * @return bin index of the given point x
593 */
594
595 template<bool INTERPOLATE>
coord2bin_calc_interpolation_weights(double * x,int ibin,int ix,int iy,int iz,int & quadrant,double & wx,double & wy,double & wz)596 inline void RegionNeighborList<INTERPOLATE>::coord2bin_calc_interpolation_weights(double *x,int ibin,int ix,int iy, int iz,int &quadrant,double &wx,double &wy,double &wz) const
597 {
598 quadrant = 0;
599 wx = wy = wz = 0.;
600 }
601
602 template<>
coord2bin_calc_interpolation_weights(double * x,int ibin,int ix,int iy,int iz,int & quadrant,double & wx,double & wy,double & wz)603 inline void RegionNeighborList<true>::coord2bin_calc_interpolation_weights(double *x,int ibin,int ix,int iy, int iz,int &quadrant,double &wx,double &wy,double &wz) const
604 {
605 double dx = (((x[0]-bboxlo[0])*bininvx) - static_cast<int> ((x[0]-bboxlo[0])*bininvx));//*bininvx;
606 double dy = (((x[1]-bboxlo[1])*bininvy) - static_cast<int> ((x[1]-bboxlo[1])*bininvy));//*bininvy;
607 double dz = (((x[2]-bboxlo[2])*bininvz) - static_cast<int> ((x[2]-bboxlo[2])*bininvz));//*bininvz;
608 quadrant = (dx >= 0.5) * 1 + (dy >= 0.5) * 2 + (dz >= 0.5) * 4;
609
610 if(dx >= 0.5)
611 wx = dx-0.5;
612 else
613 wx = 0.5+dx;
614
615 if(dy >= 0.5)
616 wy = dy-0.5;
617 else
618 wy = 0.5+dy;
619
620 if(dz >= 0.5)
621 wz = dz-0.5;
622 else
623 wz = 0.5+dz;
624
625 int stencil_shift_up = bins[ibin].stencil_shift_up[quadrant];
626 int stencil_shift_down = bins[ibin].stencil_shift_down[quadrant];
627 if( 0 == stencil_shift_up && 0 == stencil_shift_down)
628 return;
629
630 if(stencil_shift_down & 1)
631 wx = 1.;
632 else if(stencil_shift_up & 1)
633 wx = 0.;
634
635 if(stencil_shift_down & 2)
636 wy = 1.;
637 else if(stencil_shift_up & 2)
638 wy = 0.;
639
640 if(stencil_shift_down & 4)
641 wz = 1.;
642 else if(stencil_shift_up & 4)
643 wz = 0.;
644
645 }
646
647 template<bool INTERPOLATE>
coord2bin(double * x,int & quadrant,double & wx,double & wy,double & wz)648 int RegionNeighborList<INTERPOLATE>::coord2bin(double *x,int &quadrant,double &wx,double &wy,double &wz) const
649 {
650 int ix,iy,iz;
651
652 /*if(x[0] < bboxlo[0] || x[0] > bboxhi[0] || x[1] < bboxlo[1] || x[1] > bboxhi[1] || x[2] < bboxlo[2] || x[2] > bboxhi[2])
653 return -1;*/
654
655 if (x[0] >= bboxhi[0])
656 ix = static_cast<int> ((x[0]-bboxhi[0])*bininvx) + nbinx;
657 else if (x[0] >= bboxlo[0]) {
658 ix = static_cast<int> ((x[0]-bboxlo[0])*bininvx);
659 ix = std::min(ix,nbinx-1);
660 } else
661 ix = static_cast<int> ((x[0]-bboxlo[0])*bininvx) - 1;
662
663 if (x[1] >= bboxhi[1])
664 iy = static_cast<int> ((x[1]-bboxhi[1])*bininvy) + nbiny;
665 else if (x[1] >= bboxlo[1]) {
666 iy = static_cast<int> ((x[1]-bboxlo[1])*bininvy);
667 iy = std::min(iy,nbiny-1);
668 } else
669 iy = static_cast<int> ((x[1]-bboxlo[1])*bininvy) - 1;
670
671 if (x[2] >= bboxhi[2])
672 iz = static_cast<int> ((x[2]-bboxhi[2])*bininvz) + nbinz;
673 else if (x[2] >= bboxlo[2]) {
674 iz = static_cast<int> ((x[2]-bboxlo[2])*bininvz);
675 iz = std::min(iz,nbinz-1);
676 } else
677 iz = static_cast<int> ((x[2]-bboxlo[2])*bininvz) - 1;
678
679 int ibin = (iz-mbinzlo)*mbiny*mbinx + (iy-mbinylo)*mbinx + (ix-mbinxlo);
680
681 if(ibin < 0 || ibin >= mbins())
682 return -1; // ibin = -1
683 coord2bin_calc_interpolation_weights(x,ibin,ix,iy,iz,quadrant,wx,wy,wz);
684
685 return ibin;
686 }
687
688 /**
689 * @brief Calc global bin index (n) of point x
690 * @param x point in 3D
691 * @return bin index of the given point x
692 */
693
694 template<bool INTERPOLATE>
isInBoundingBox(double * pos)695 bool RegionNeighborList<INTERPOLATE>::isInBoundingBox(double *pos) const
696 {
697 if(!boundingBoxSet())
698 error->one(FLERR,"internal error: call before bbox is set");
699 if
700 (
701 pos[0] >= bboxlo[0] && pos[0] <= bboxhi[0] &&
702 pos[1] >= bboxlo[1] && pos[1] <= bboxhi[1] &&
703 pos[2] >= bboxlo[2] && pos[2] <= bboxhi[2]
704 ) return true;
705 return false;
706 }
707
708 /**
709 * @brief Set bounding box from region
710 * @param region Region used for creating the neighbor list
711 * @param maxrad largest particle radius
712 * @param extend enable extending the box for ghost particles
713 * @param failsafe limit max. number of bins
714 * @return BoundingBox the bounding box of the region
715 */
716
717 template<bool INTERPOLATE>
setBoundingBoxRegion(const Region & region,double maxrad,bool extend,bool failsafe)718 BoundingBox RegionNeighborList<INTERPOLATE>::setBoundingBoxRegion(const Region ®ion, double maxrad, bool extend, bool failsafe)
719 {
720 // use region to determine and set bounding box
721 BoundingBox bb(region.extent_xlo, region.extent_xhi, region.extent_ylo, region.extent_yhi, region.extent_zlo, region.extent_zhi);
722 if (!setBoundingBox(bb, maxrad, extend, failsafe))
723 error->one(FLERR,"internal error: could not set bounding box");
724
725 return bb;
726 }
727