MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Mercury3D.cc
Go to the documentation of this file.
1 //Copyright (c) 2013-2020, The MercuryDPM Developers Team. All rights reserved.
2 //For the list of developers, see <http://www.MercuryDPM.org/Team>.
3 //
4 //Redistribution and use in source and binary forms, with or without
5 //modification, are permitted provided that the following conditions are met:
6 // * Redistributions of source code must retain the above copyright
7 // notice, this list of conditions and the following disclaimer.
8 // * Redistributions in binary form must reproduce the above copyright
9 // notice, this list of conditions and the following disclaimer in the
10 // documentation and/or other materials provided with the distribution.
11 // * Neither the name MercuryDPM nor the
12 // names of its contributors may be used to endorse or promote products
13 // derived from this software without specific prior written permission.
14 //
15 //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
16 //ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
17 //WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18 //DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
19 //DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
20 //(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
21 //LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
22 //ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23 //(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24 //SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 
26 #include "Mercury3D.h"
27 #include "Particles/BaseParticle.h"
28 
30 {
31  constructor();
32  logger(DEBUG, "Mercury3D::Mercury3D() finished");
33 }
34 
41  : DPMBase(other), MercuryBase(other)
42 {
43  logger(DEBUG, "Mercury3D::Mercury3D(Mercury3D& other) copy constructor finished.");
44 }
45 
53  : DPMBase(other), MercuryBase()
54 {
55  constructor();
56  logger(DEBUG, "Mercury3D::Mercury3D(DPMBase& other) constructor finished");
57 }
58 
60 {
63 }
64 
73 void Mercury3D::hGridFindContactsWithinTargetCell(int x, int y, int z, unsigned int l)
74 {
75  HGrid* const hgrid = getHGrid();
76  const unsigned int bucket = hgrid->computeHashBucketIndex(x, y, z, l);
77 
79  //Check if this function is already applied to this bucket
80  bool bucketIsChecked;
81  #pragma omp critical
82  {
83  bucketIsChecked = hgrid->getBucketIsChecked(bucket);
84  hgrid->setBucketIsChecked(bucket);
85  }
86  if (bucketIsChecked) return;
87 
88  BaseParticle* p1 = hgrid->getFirstBaseParticleInBucket(bucket);
89  while (p1 != nullptr)
90  {
91  BaseParticle* p2 = p1->getHGridNextObject();
92  while (p2 != nullptr)
93  {
96  //Check if the BaseParticle* p1 and BaseParticle* p2 are really in the same cell (i.e. no hashing error has occurred)
97  if (p1->getHGridCell() == (p2->getHGridCell()))
98  {
99  computeInternalForce(p1, p2);
100  }
101  p2 = p2->getHGridNextObject();
102  }
103  p1 = p1->getHGridNextObject();
104  }
105 }
106 
118 void Mercury3D::hGridFindContactsWithTargetCell(int x, int y, int z, unsigned int l, BaseParticle* const obj)
119 {
120  //Check if the object is not in the same cell as being checked, CheckCell_current should handle these cases.
121  //TW a speedcheck revealed that this check costs a 10% performance decrease; it's only a safety check, so I made it an assert.
122  logger.assert(!obj->getHGridCell().equals(x, y, z, l),
123  "hGridFindContactsWithTargetCell should not be called if object is in the same cell");
124 
125  HGrid* const hgrid = getHGrid();
126 
127  // Calculate the bucket
128  const unsigned int bucket = hgrid->computeHashBucketIndex(x, y, z, l);
129 
130  // Loop through all objects in the bucket to find nearby objects
131  for (BaseParticle* p = hgrid->getFirstBaseParticleInBucket(bucket); p != nullptr; p = p->getHGridNextObject())
132  {
133  //This is the most-expensive function in most codes (the two checks in this function slows down granular jet by 15%). It is neccesary, for example: If two cells hash to the same bucket and a particle in one of these cells check for collisions with the other cell. Then due to the hashing collision it also gets all particles in it's own cell and thus generating false collisions.
134  //Check if the BaseParticle *p really is in the target cell (i.e. no hashing error has occurred)
135  //TW speedcheck revealed that this pre-check is cheaper than allowing computeInternalForces to sort out mismatches; even if a large number of hash cells (10*Np) is used.
136  if (p->getHGridCell().equals(x, y, z, l))
137  {
138  if (Vec3D::getDistanceSquared(p->getPosition(),obj->getPosition()) < mathsFunc::square(p->getMaxInteractionRadius()+obj->getMaxInteractionRadius()))
139  computeInternalForce(obj, p);
140  }
141  }
142 }
143 
144 //TODO document
145 //TODO obj kan weg?
146 void Mercury3D::hGridFindParticlesWithTargetCell(int x, int y, int z, unsigned int l, BaseParticle* obj,
147  std::vector<BaseParticle*>& list)
148 {
149  HGrid* const hgrid = getHGrid();
150 
151  // Calculate the bucket
152  const unsigned int bucket = hgrid->computeHashBucketIndex(x, y, z, l);
153 
154  // Loop through all objects in the bucket to find nearby objects
155  BaseParticle* p = hgrid->getFirstBaseParticleInBucket(bucket);
156  while (p != nullptr)
157  {
158  if (p->getHGridCell().equals(x, y, z, l))
159  {
160  list.push_back(p);
161  }
162  p = p->getHGridNextObject();
163  }
164 }
165 
166 void Mercury3D::hGridGetInteractingParticleList(BaseParticle* obj, std::vector<BaseParticle*>& list)
167 {
168  HGrid* hgrid = getHGrid();
169 
171  if (hGridNeedsRebuilding())
172  {
173  hGridRebuild();
174  hgrid = getHGrid();
175  }
176  logger(DEBUG, "hgrid %, object %", hgrid, obj);
177  int occupiedLevelsMask = hgrid->getOccupiedLevelsMask() >> obj->getHGridLevel();
178  for (unsigned int level = 0; level < hgrid->getNumberOfLevels(); level++)
179  {
180  // If no objects in rest of grid, stop now
181  if (occupiedLevelsMask == 0)
182  {
183  break;
184  }
185 
186  // If no objects at this level, go on to the next level
187  if ((occupiedLevelsMask & 1) == 0)
188  {
189  continue;
190  }
191 
192  const Mdouble inv_size = hgrid->getInvCellSize(level);
193  const int xs = static_cast<int>(std::floor(
194  (obj->getPosition().X - obj->getMaxInteractionRadius()) * inv_size - 0.5));
195  const int xe = static_cast<int>(std::floor(
196  (obj->getPosition().X + obj->getMaxInteractionRadius()) * inv_size + 0.5));
197  const int ys = static_cast<int>(std::floor(
198  (obj->getPosition().Y - obj->getMaxInteractionRadius()) * inv_size - 0.5));
199  const int ye = static_cast<int>(std::floor(
200  (obj->getPosition().Y + obj->getMaxInteractionRadius()) * inv_size + 0.5));
201  const int zs = static_cast<int>(std::floor(
202  (obj->getPosition().Z - obj->getMaxInteractionRadius()) * inv_size - 0.5));
203  const int ze = static_cast<int>(std::floor(
204  (obj->getPosition().Z + obj->getMaxInteractionRadius()) * inv_size + 0.5));
205  for (int x = xs; x <= xe; ++x)
206  {
207  for (int y = ys; y <= ye; ++y)
208  {
209  for (int z = zs; z <= ze; ++z)
210  {
211  hGridFindParticlesWithTargetCell(x, y, z, level, obj, list);
212  }
213  }
214  }
215  }
216 }
217 
225 {
226  HGrid* const hgrid = getHGrid();
227  const unsigned int startLevel = obj->getHGridLevel();
228 
229  if (getHGridMethod() == TOPDOWN)
230  {
231  int occupiedLevelsMask = hgrid->getOccupiedLevelsMask();
232  for (unsigned int level = 0; level <= startLevel && occupiedLevelsMask != 0; occupiedLevelsMask >>= 1, level++)
233  {
234  // If no objects at this level, go on to the next level
235  if ((occupiedLevelsMask & 1) == 0)
236  {
237  continue;
238  }
239 
240  if (level == startLevel)
241  {
242  const int x = obj->getHGridX();
243  const int y = obj->getHGridY();
244  const int z = obj->getHGridZ();
245 
246  hGridFindContactsWithinTargetCell(x, y, z, level);
247  hGridFindContactsWithTargetCell(x + 1, y - 1, z, level, obj);
248  hGridFindContactsWithTargetCell(x + 1, y, z, level, obj);
249  hGridFindContactsWithTargetCell(x + 1, y + 1, z, level, obj);
250  hGridFindContactsWithTargetCell(x + 1, y - 1, z + 1, level, obj);
251  hGridFindContactsWithTargetCell(x + 1, y, z + 1, level, obj);
252  hGridFindContactsWithTargetCell(x + 1, y + 1, z + 1, level, obj);
253  hGridFindContactsWithTargetCell(x + 1, y - 1, z - 1, level, obj);
254  hGridFindContactsWithTargetCell(x + 1, y, z - 1, level, obj);
255  hGridFindContactsWithTargetCell(x + 1, y + 1, z - 1, level, obj);
256  hGridFindContactsWithTargetCell(x, y + 1, z, level, obj);
257  hGridFindContactsWithTargetCell(x, y, z - 1, level, obj);
258  hGridFindContactsWithTargetCell(x, y + 1, z - 1, level, obj);
259  hGridFindContactsWithTargetCell(x, y + 1, z + 1, level, obj);
260  }
261  else
262  {
263  const Mdouble inv_size = getHGrid()->getInvCellSize(level);
264  const int xs = static_cast<int>(std::floor(
265  (obj->getPosition().X - obj->getMaxInteractionRadius()) * inv_size - 0.5));
266  const int xe = static_cast<int>(std::floor(
267  (obj->getPosition().X + obj->getMaxInteractionRadius()) * inv_size + 0.5));
268  const int ys = static_cast<int>(std::floor(
269  (obj->getPosition().Y - obj->getMaxInteractionRadius()) * inv_size - 0.5));
270  const int ye = static_cast<int>(std::floor(
271  (obj->getPosition().Y + obj->getMaxInteractionRadius()) * inv_size + 0.5));
272  const int zs = static_cast<int>(std::floor(
273  (obj->getPosition().Z - obj->getMaxInteractionRadius()) * inv_size - 0.5));
274  const int ze = static_cast<int>(std::floor(
275  (obj->getPosition().Z + obj->getMaxInteractionRadius()) * inv_size + 0.5));
276  for (int x = xs; x <= xe; ++x)
277  {
278  for (int y = ys; y <= ye; ++y)
279  {
280  for (int z = zs; z <= ze; ++z)
281  {
282  hGridFindContactsWithTargetCell(x, y, z, level, obj);
283  }
284  }
285  }
286  }
287  }
288  }
289  else
290  {
291  int occupiedLevelsMask = hgrid->getOccupiedLevelsMask() >> obj->getHGridLevel();
292  for (unsigned int level = startLevel; level < hgrid->getNumberOfLevels(); occupiedLevelsMask >>= 1, level++)
293  {
294  // If no objects in rest of grid, stop now
295  if (occupiedLevelsMask == 0)
296  {
297  break;
298  }
299 
300  // If no objects at this level, go on to the next level
301  if ((occupiedLevelsMask & 1) == 0)
302  {
303  continue;
304  }
305 
306  if (level == startLevel)
307  {
308  const int x = obj->getHGridX();
309  const int y = obj->getHGridY();
310  const int z = obj->getHGridZ();
311 
312  hGridFindContactsWithinTargetCell(x, y, z, level);
313  hGridFindContactsWithTargetCell(x + 1, y - 1, z, level, obj);
314  hGridFindContactsWithTargetCell(x + 1, y, z, level, obj);
315  hGridFindContactsWithTargetCell(x + 1, y + 1, z, level, obj);
316  hGridFindContactsWithTargetCell(x + 1, y - 1, z + 1, level, obj);
317  hGridFindContactsWithTargetCell(x + 1, y, z + 1, level, obj);
318  hGridFindContactsWithTargetCell(x + 1, y + 1, z + 1, level, obj);
319  hGridFindContactsWithTargetCell(x + 1, y - 1, z - 1, level, obj);
320  hGridFindContactsWithTargetCell(x + 1, y, z - 1, level, obj);
321  hGridFindContactsWithTargetCell(x + 1, y + 1, z - 1, level, obj);
322  hGridFindContactsWithTargetCell(x, y + 1, z, level, obj);
323  hGridFindContactsWithTargetCell(x, y, z - 1, level, obj);
324  hGridFindContactsWithTargetCell(x, y + 1, z - 1, level, obj);
325  hGridFindContactsWithTargetCell(x, y + 1, z + 1, level, obj);
326  }
327  else
328  {
329  const Mdouble inv_size = hgrid->getInvCellSize(level);
330  const int xs = static_cast<int>(std::floor(
331  (obj->getPosition().X - obj->getMaxInteractionRadius()) * inv_size - 0.5));
332  const int xe = static_cast<int>(std::floor(
333  (obj->getPosition().X + obj->getMaxInteractionRadius()) * inv_size + 0.5));
334  const int ys = static_cast<int>(std::floor(
335  (obj->getPosition().Y - obj->getMaxInteractionRadius()) * inv_size - 0.5));
336  const int ye = static_cast<int>(std::floor(
337  (obj->getPosition().Y + obj->getMaxInteractionRadius()) * inv_size + 0.5));
338  const int zs = static_cast<int>(std::floor(
339  (obj->getPosition().Z - obj->getMaxInteractionRadius()) * inv_size - 0.5));
340  const int ze = static_cast<int>(std::floor(
341  (obj->getPosition().Z + obj->getMaxInteractionRadius()) * inv_size + 0.5));
342  for (int x = xs; x <= xe; ++x)
343  {
344  for (int y = ys; y <= ye; ++y)
345  {
346  for (int z = zs; z <= ze; ++z)
347  {
348  hGridFindContactsWithTargetCell(x, y, z, level, obj);
349  }
350  }
351  }
352  }
353  }
354  }
355 }
356 
362 {
363  HGrid* const hGrid = getHGrid();
364  if (hGrid)
365  {
366  const unsigned int l = obj->getHGridLevel();
367  const Mdouble inv_size = hGrid->getInvCellSize(l);
368 
369  int x = static_cast<int>(std::floor(obj->getPosition().X * inv_size));
370  int y = static_cast<int>(std::floor(obj->getPosition().Y * inv_size));
371  int z = static_cast<int>(std::floor(obj->getPosition().Z * inv_size));
372 
373 #ifdef CONTACT_LIST_HGRID
374  if(obj->getHGridX() != x || obj->getHGridY() != y || obj->getHGridZ() != z)
375  {
376  int bucket = hGrid->computeHashBucketIndex(x, y, z, l);
377 
378  //First the object has to be removed
379  hGridRemoveParticle(obj);
380 
381  //Also remove all contact associated with it
382  getPossibleContactList().remove_ParticlePosibleContacts(obj);
383 
384  //And now reinserted
386  obj->setHGridPrevObject(nullptr);
387  if(hGrid->getFirstBaseParticleInBucket(bucket))
388  {
390  }
391  hGrid->setFirstBaseParticleInBucket(bucket,obj);
392 
393  obj->setHGridX(x);
394  obj->setHGridY(y);
395  obj->setHGridZ(z);
396  InsertObjAgainstGrid(obj);
397  }
398 #else
399  const unsigned int bucket = hGrid->computeHashBucketIndex(x, y, z, l);
400 
401  // this needs to be defined as #pragma omp critical if MercuryBase::hGridActionsBeforeTimeStep is parallelised; however, parallelising it make the code slower, not faster.
402  {
404  obj->setHGridPrevObject(nullptr);
405  if (hGrid->getFirstBaseParticleInBucket(bucket)) {
407  }
408  hGrid->setFirstBaseParticleInBucket(bucket, obj);
409  }
410 
411  obj->setHGridX(x);
412  obj->setHGridY(y);
413  obj->setHGridZ(z);
414 #endif
415  }
416 }
417 
423 {
424  HGrid* const hGrid = getHGrid();
425  if (hGrid)
426  {
427  const unsigned int bucket = hGrid->computeHashBucketIndex(obj->getHGridCell());
428  if (obj->getHGridPrevObject())
429  {
431  }
432  else
433  {
434  if (hGrid->getFirstBaseParticleInBucket(bucket) == obj)
435  {
436  hGrid->setFirstBaseParticleInBucket(bucket, obj->getHGridNextObject());
437  }
438  }
439 
440  if (obj->getHGridNextObject())
441  {
443  }
444  }
445 }
446 
456 bool Mercury3D::hGridHasContactsInTargetCell(int x, int y, int z, unsigned int l, const BaseParticle* obj) const
457 {
458  // Loop through all objects in the bucket to find nearby objects
459  const unsigned int bucket = getHGrid()->computeHashBucketIndex(x, y, z, l);
460 
461  const BaseParticle* p = getHGrid()->getFirstBaseParticleInBucket(bucket);
462  while (p != nullptr)
463  {
464  if (p->getHGridCell().equals(x, y, z, l))
465  {
466  if (areInContact(obj, p))
467  {
468  return true;
469  }
470  }
471  //std::cout << "HERE!" << std::endl;
472  p = p->getHGridNextObject();
473  }
474  return false;
475 }
476 
486 {
487  if (getHGrid() == nullptr || getHGrid()->getNeedsRebuilding())
488  {
489  logger(INFO, "HGrid needs rebuilding for \"bool Mercury3D::hGridHasParticleContacts(BaseParticle *obj)\"");
490  hGridRebuild();
491  }
492 
493  int occupiedLevelsMask = getHGrid()->getOccupiedLevelsMask();
494 
495  for (unsigned int level = 0; level < getHGrid()->getNumberOfLevels(); occupiedLevelsMask >>= 1, level++)
496  {
497  // If no objects in rest of grid, stop now
498  if (occupiedLevelsMask == 0)
499  {
500  logger(VERBOSE, "Level % and higher levels are empty", level);
501  break;
502  }
503 
504  // If no objects at this level, go on to the next level
505  if ((occupiedLevelsMask & 1) == 0)
506  {
507  logger(VERBOSE, "Level % is empty", level);
508  continue;
509  }
510 
511  const Mdouble inv_size = getHGrid()->getInvCellSize(level);
512  const int xs = static_cast<int>(std::floor(
513  (obj->getPosition().X - obj->getMaxInteractionRadius()) * inv_size - 0.5));
514  const int xe = static_cast<int>(std::floor(
515  (obj->getPosition().X + obj->getMaxInteractionRadius()) * inv_size + 0.5));
516  const int ys = static_cast<int>(std::floor(
517  (obj->getPosition().Y - obj->getMaxInteractionRadius()) * inv_size - 0.5));
518  const int ye = static_cast<int>(std::floor(
519  (obj->getPosition().Y + obj->getMaxInteractionRadius()) * inv_size + 0.5));
520  const int zs = static_cast<int>(std::floor(
521  (obj->getPosition().Z - obj->getMaxInteractionRadius()) * inv_size - 0.5));
522  const int ze = static_cast<int>(std::floor(
523  (obj->getPosition().Z + obj->getMaxInteractionRadius()) * inv_size + 0.5));
524 
525  logger(VERBOSE, "Level = % grid cells [%,%] x [%,%] x [%,%]", level, xs, xe, ys, ye, zs, ze);
526  for (int x = xs; x <= xe; ++x)
527  {
528  for (int y = ys; y <= ye; ++y)
529  {
530  for (int z = zs; z <= ze; ++z)
531  {
532  if (hGridHasContactsInTargetCell(x, y, z, level, obj))
533  {
534  return true;
535  }
536  }
537  }
538  }
539  } //end for level
540 
541  return false;
542 }
543 
544 std::vector<BaseParticle*> Mercury3D::hGridFindParticleContacts(const BaseParticle* obj)
545 {
546  if (getHGrid() == nullptr || getHGrid()->getNeedsRebuilding())
547  {
548  logger(INFO, "HGrid needs rebuilding for \"bool Mercury3D::hGridHasParticleContacts(BaseParticle *obj)\"");
549  hGridRebuild();
550  }
551 
552  int occupiedLevelsMask = getHGrid()->getOccupiedLevelsMask();
553 
554  std::vector<BaseParticle*> particlesInContact;
555 
556  for (unsigned int level = 0; level < getHGrid()->getNumberOfLevels(); occupiedLevelsMask >>= 1, level++)
557  {
558  // If no objects in rest of grid, stop now
559  if (occupiedLevelsMask == 0)
560  {
561  logger(VERBOSE, "Level % and higher levels are empty", level);
562  break;
563  }
564 
565  // If no objects at this level, go on to the next level
566  if ((occupiedLevelsMask & 1) == 0)
567  {
568  logger(VERBOSE, "Level % is empty", level);
569  continue;
570  }
571 
572  const Mdouble inv_size = getHGrid()->getInvCellSize(level);
573  const int xs = static_cast<int>(std::floor(
574  (obj->getPosition().X - obj->getMaxInteractionRadius()) * inv_size - 0.5));
575  const int xe = static_cast<int>(std::floor(
576  (obj->getPosition().X + obj->getMaxInteractionRadius()) * inv_size + 0.5));
577  const int ys = static_cast<int>(std::floor(
578  (obj->getPosition().Y - obj->getMaxInteractionRadius()) * inv_size - 0.5));
579  const int ye = static_cast<int>(std::floor(
580  (obj->getPosition().Y + obj->getMaxInteractionRadius()) * inv_size + 0.5));
581  const int zs = static_cast<int>(std::floor(
582  (obj->getPosition().Z - obj->getMaxInteractionRadius()) * inv_size - 0.5));
583  const int ze = static_cast<int>(std::floor(
584  (obj->getPosition().Z + obj->getMaxInteractionRadius()) * inv_size + 0.5));
585 
586  logger(VERBOSE, "Level = % grid cells [%,%] x [%,%] x [%,%]", level, xs, xe, ys, ye, zs, ze);
587  for (int x = xs; x <= xe; ++x)
588  {
589  for (int y = ys; y <= ye; ++y)
590  {
591  for (int z = zs; z <= ze; ++z)
592  {
593  // Loop through all objects in the bucket to find nearby objects
594  const unsigned int bucket = getHGrid()->computeHashBucketIndex(x, y, z, level);
596  while (p != nullptr)
597  {
598  if (p->getHGridCell().equals(x, y, z, level))
599  {
600  if (areInContact(obj, p))
601  {
602  particlesInContact.push_back(p);
603  }
604  }
605  p = p->getHGridNextObject();
606  }
607  }
608  }
609  }
610  } //end for level
611 
612  return particlesInContact;
613 }
614 
616 {
617 
618  // if wall is not local, use the non-hGrid version for finding wall contacts
619  Vec3D min, max;
620  if (w->isLocal(min, max)==false)
621  {
622  return DPMBase::computeWallForces(w);
623  }
624 
625  //compute forces for all particles that are neither fixed or ghosts
626  if (getHGrid() == nullptr || getHGrid()->getNeedsRebuilding())
627  {
628  logger(INFO, "HGrid needs rebuilding for \"bool Mercury3D::hGridHasParticleContacts(BaseParticle *obj)\"");
629  hGridRebuild();
630  }
631 
632  HGrid* const hGrid = getHGrid();
633 
634  int occupiedLevelsMask = hGrid->getOccupiedLevelsMask();
635 
636  for (unsigned int level = 0; level < hGrid->getNumberOfLevels(); occupiedLevelsMask >>= 1, level++)
637  {
638  // If no objects in rest of grid, stop now
639  if (occupiedLevelsMask == 0)
640  {
641  logger(VERBOSE, "Level % and higher levels are empty", level);
642  break;
643  }
644 
645  // If no objects at this level, go on to the next level
646  if ((occupiedLevelsMask & 1) == 0)
647  {
648  logger(VERBOSE, "Level % is empty", level);
649  continue;
650  }
651 
652  const Mdouble inv_size = hGrid->getInvCellSize(level);
653  const int xs = static_cast<int>(std::floor(min.X * inv_size - 0.5));
654  const int xe = static_cast<int>(std::floor(max.X * inv_size + 0.5));
655  const int ys = static_cast<int>(std::floor(min.Y * inv_size - 0.5));
656  const int ye = static_cast<int>(std::floor(max.Y * inv_size + 0.5));
657  const int zs = static_cast<int>(std::floor(min.Z * inv_size - 0.5));
658  const int ze = static_cast<int>(std::floor(max.Z * inv_size + 0.5));
659  //logger(INFO, "Level % grid cells [%,%] x [%,%] x [%,%]", level, xs, xe, ys, ye, zs, ze);
660 
661  for (int x = xs; x <= xe; ++x)
662  {
663  for (int y = ys; y <= ye; ++y)
664  {
665  for (int z = zs; z <= ze; ++z)
666  {
667  // Loop through all objects in the bucket to find nearby objects
668  const unsigned int bucket = hGrid->computeHashBucketIndex(x, y, z, level);
669  BaseParticle* p = hGrid->getFirstBaseParticleInBucket(bucket);
670  while (p != nullptr)
671  {
672  if (!p->isFixed() && p->getPeriodicFromParticle() == nullptr &&
673  p->getHGridCell().equals(x, y, z, level))
674  {
675  //logger(INFO, "t % p % Size % level % cells % % %", getNumberOfTimeSteps(), p->getIndex(), hGrid->getCellSize(level), level, x,y,z);
677  //w->computeForces(p);
678  }
679  p = p->getHGridNextObject();
680  }
681  }
682  }
683  }
684  } //end for level
685 }
686 
687 #ifdef CONTACT_LIST_HGRID
688 
697 void Mercury3D::InsertCell(int x, int y, int z, unsigned int l, BaseParticle *obj)
698 {
699  // Loop through all objects in the bucket to find nearby objects
700  unsigned int bucket = getHGrid()->computeHashBucketIndex(x,y,z,l);
702 
703  while (p!=nullptr)
704  {
705  if ((p->getHGridX() == x) && (p->getHGridY() == y) && (p->getHGridZ() == z) && (p->getHGridLevel() == l) && (obj!=p))
706  {
707  getPossibleContactList().add_PossibleContact(obj,p);
708  }
709  p = p->getHGridNextObject();
710  }
711 }
712 
720 void Mercury3D::InsertObjAgainstGrid(BaseParticle *obj)
721 {
722  Mdouble inv_size;
723  int occupiedLevelsMask_ = getHGrid()->getOccupiedLevelsMask();
724 
725  inv_size=getHGrid()->getInvCellSize(obj->getHGridLevel());
726 
727  double ownXMin = (obj->getHGridX() - 0.5) * getHGrid()->getCellSize(obj->getHGridLevel());
728  double ownXMax = (obj->getHGridX() + 1.5) * getHGrid()->getCellSize(obj->getHGridLevel());
729  double ownYMin = (obj->getHGridY() - 0.5) * getHGrid()->getCellSize(obj->getHGridLevel());
730  double ownYMax = (obj->getHGridY() + 1.5) * getHGrid()->getCellSize(obj->getHGridLevel());
731  double ownZMin = (obj->getHGridZ() - 0.5) * getHGrid()->getCellSize(obj->getHGridLevel());
732  double ownZMax = (obj->getHGridZ() + 1.5) * getHGrid()->getCellSize(obj->getHGridLevel());
733 
734  for (int level = 0; level < getHGrid()->getNumberOfLevels(); occupiedLevelsMask_ >>= 1, level++)
735  {
736  // If no objects in rest of grid, stop now
737  if (occupiedLevelsMask_ == 0)
738  {
739  break;
740  }
741 
742  // If no objects at this level, go on to the next level
743  if ((occupiedLevelsMask_ & 1) == 0)
744  {
745  continue;
746  }
747 
748  // Treat level as a third dimension coordinate
749  inv_size = getHGrid()->getInvCellSize(level);
750 
751  int xs, xe, ys, ye, zs, ze;
752  xs=static_cast<int>(std::floor(ownXMin * inv_size - 0.5));
753  xe=static_cast<int>(std::floor(ownXMax * inv_size + 0.5));
754  ys=static_cast<int>(std::floor(ownYMin * inv_size - 0.5));
755  ye=static_cast<int>(std::floor(ownYMax * inv_size + 0.5));
756  zs=static_cast<int>(std::floor(ownZMin * inv_size - 0.5));
757  ze=static_cast<int>(std::floor(ownZMax * inv_size + 0.5));
758 
759  for(int x=xs; x<=xe; ++x)
760  {
761  for(int y=ys; y<=ye; ++y)
762  {
763  for(int z=zs; z<=ze; ++z)
764  {
765  InsertCell(x, y, z, level, obj);
766  }
767  }
768  }
769  } //end for level
770 }
771 #endif
const HGridCell & getHGridCell() const
Definition: BaseParticle.h:645
bool equals(int x, int y, int z, unsigned int level) const
Checks if the given (x,y,z,level) is the same as the ones in this cell.
Definition: HGridCell.h:39
void computeWallForces(BaseWall *w) override
Compute contacts with a wall.
Definition: Mercury3D.cc:615
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
The DPMBase header includes quite a few header files, defining all the handlers, which are essential...
Definition: DPMBase.h:72
Mdouble X
the vector components
Definition: Vector.h:65
bool hGridNeedsRebuilding()
Gets if the HGrid needs rebuilding before anything else happens.
Definition: MercuryBase.cc:490
void hGridRemoveParticle(BaseParticle *obj) override
Removes a BaseParticle from the HGrid.
Definition: Mercury3D.cc:422
static bool areInContact(const BaseParticle *pI, const BaseParticle *pJ)
Checks if two particle are in contact or is there any positive overlap.
Definition: DPMBase.cc:1621
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
double Mdouble
Definition: GeneralDefine.h:34
Mdouble getMaxInteractionRadius() const
Returns the particle's interaction radius, which might be different from radius_ (e.g., when dealing with wet particles)
Definition: BaseParticle.h:359
void hGridFindContactsWithinTargetCell(int x, int y, int z, unsigned int l)
Finds contacts between particles in the target cell.
Definition: Mercury3D.cc:73
void constructor()
Function that sets the SystemDimension and ParticleDimension to 3.
Definition: Mercury3D.cc:59
bool hGridHasContactsInTargetCell(int x, int y, int z, unsigned int l, const BaseParticle *obj) const
Tests if the BaseParticle has contacts with other Particles in the target cell.
Definition: Mercury3D.cc:456
virtual bool isLocal(Vec3D &min, Vec3D &max) const
Definition: BaseWall.h:137
void setHGridY(const int y)
Sets the particle's HGrid cell Y-coordinate.
Definition: BaseParticle.h:454
void computeInternalForces(BaseParticle *obj) override
Finds contacts with the BaseParticle; avoids multiple checks.
Definition: Mercury3D.cc:224
void setHGridNextObject(BaseParticle *p)
Sets the pointer to the next object in the particle's HGrid cell & level.
Definition: BaseParticle.h:477
void setParticleDimensions(unsigned int particleDimensions)
Sets the particle dimensionality.
Definition: DPMBase.cc:1408
HGridMethod getHGridMethod() const
Gets whether the HGrid in this MercuryBase is BOTTOMUP or TOPDOWN.
Definition: MercuryBase.h:204
BaseParticle * getPeriodicFromParticle() const
Returns the 'original' particle this one's a periodic copy of.
Definition: BaseParticle.h:338
virtual void computeWallForces(BaseWall *w)
Definition: DPMBase.cc:5172
int getHGridY() const
Returns particle's HGrid cell Y-coordinate.
Definition: BaseParticle.h:272
void setSystemDimensions(unsigned int newDim)
Sets the system dimensionality.
Definition: DPMBase.cc:1377
void setHGridX(const int x)
Sets the particle's HGrid cell X-coordinate.
Definition: BaseParticle.h:446
void hGridUpdateParticle(BaseParticle *obj) override
Updates the cell (not the level) of a BaseParticle.
Definition: Mercury3D.cc:361
std::vector< BaseParticle * > hGridFindParticleContacts(const BaseParticle *obj) override
Returns all particles that have a contact with a given particle.
Definition: Mercury3D.cc:544
int getOccupiedLevelsMask() const
Gets the integer that represents which levels are occupied.
Definition: HGrid.h:213
void hGridGetInteractingParticleList(BaseParticle *obj, std::vector< BaseParticle * > &list) override
Obtains all neighbour particles of a given object, obtained from the hgrid.
Definition: Mercury3D.cc:166
This is the base class for both Mercury2D and Mercury3D. Note the actually abstract grid is defined i...
Definition: MercuryBase.h:125
Mercury3D()
This is the default constructor. All it does is set sensible defaults.
Definition: Mercury3D.cc:29
In the HGrid class, here all information about the HGrid is stored.
Definition: HGrid.h:42
int getHGridZ() const
Returns particle's HGrid cell Z-coordinate.
Definition: BaseParticle.h:279
This adds on the hierarchical grid code for 3D problems.
Definition: Mercury3D.h:36
void setFirstBaseParticleInBucket(unsigned int i, BaseParticle *p)
Sets the first particle in bucket i to be the given BaseParticle.
Definition: HGrid.h:117
void setHGridZ(const int z)
Sets the particle's HGrid cell Z-coordinate.
Definition: BaseParticle.h:462
void setBucketIsChecked(unsigned int i)
Sets that the bucket with the given index is checked to true.
Definition: HGrid.h:124
bool hGridHasParticleContacts(const BaseParticle *obj) override
Tests if a BaseParticle has any contacts in the HGrid.
Definition: Mercury3D.cc:485
Basic class for walls.
Definition: BaseWall.h:47
bool isFixed() const override
Is fixed Particle function. It returns whether a Particle is fixed or not, by checking its inverse Ma...
Definition: BaseParticle.h:93
static Mdouble getDistanceSquared(const Vec3D &a, const Vec3D &b)
Calculates the squared distance between two Vec3D: .
Definition: Vector.h:295
Mdouble Y
Definition: Vector.h:65
void hGridRebuild()
This sets up the parameters required for the contact model.
Definition: MercuryBase.cc:203
virtual void computeInternalForce(BaseParticle *, BaseParticle *)
Computes the forces between two particles (internal in the sense that the sum over all these forces i...
Definition: DPMBase.cc:2994
double getCellSize(unsigned int i) const
Gets the size of the cells at the given level.
Definition: HGrid.h:147
BaseParticle * getHGridPrevObject() const
Returns pointer to previous object in particle's HGrid level & cell.
Definition: BaseParticle.h:250
void hGridFindContactsWithTargetCell(int x, int y, int z, unsigned int l, BaseParticle *obj)
Finds contacts between the BaseParticle and the target cell.
Definition: Mercury3D.cc:118
void setHGridPrevObject(BaseParticle *p)
Sets the pointer to the previous object in the particle's HGrid cell & level.
Definition: BaseParticle.h:485
HGrid * getHGrid()
Gets the HGrid used by this problem.
Definition: MercuryBase.h:311
void computeForcesDueToWalls(BaseParticle *, BaseWall *)
Computes the forces on the particles due to the walls (normals are outward normals) ...
Definition: DPMBase.cc:3069
unsigned long getNumberOfLevels() const
Gets the number of levels of this HGrid.
Definition: HGrid.h:206
const BaseParticle * getFirstBaseParticleInBucket(unsigned int i) const
Gets the first BaseParticle in the given bucket, const version.
Definition: HGrid.h:162
void hGridFindParticlesWithTargetCell(int x, int y, int z, unsigned int l, BaseParticle *obj, std::vector< BaseParticle * > &list)
Finds particles within target cell and stores them in a list.
Definition: Mercury3D.cc:146
unsigned int computeHashBucketIndex(int x, int y, int z, unsigned int l) const
Computes hash bucket index in range [0, NUM_BUCKETS-1] for a 3D domain.
Definition: HGrid.h:76
Definition: Vector.h:49
BaseParticle * getHGridNextObject() const
Returns pointer to next object in particle's HGrid level & cell.
Definition: BaseParticle.h:242
T square(const T val)
squares a number
Definition: ExtendedMath.h:104
unsigned int getHGridLevel() const
Returns particle's HGrid level.
Definition: BaseParticle.h:234
Mdouble Z
Definition: Vector.h:65
bool getBucketIsChecked(unsigned int i) const
Gets whether or not the bucket with index i is checked.
Definition: HGrid.h:132
double getInvCellSize(unsigned int i) const
Gets 1/cellSize for the cells on level i.
Definition: HGrid.h:178
int getHGridX() const
Returns particle's HGrid cell X-coordinate.
Definition: BaseParticle.h:265