MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ParticleHandler.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 
27 #include <limits>
28 #include <Math/Helpers.h>
31 #include "ParticleHandler.h"
32 #include "DPMBase.h"
33 #include "SpeciesHandler.h"
36 #include "Particles/BaseParticle.h"
38 #include "MpiContainer.h"
39 #include "MpiDataClass.h"
40 #include "BoundaryHandler.h"
41 #include "DomainHandler.h"
42 #include "Domain.h"
43 
49 {
50  largestParticle_ = nullptr;
51  smallestParticle_ = nullptr;
52  logger(DEBUG, "ParticleHandler::ParticleHandler() finished");
53 }
54 
63 {
64  clear();
65  setDPMBase(PH.getDPMBase());
66  largestParticle_ = nullptr;
67  smallestParticle_ = nullptr;
69  if (!objects_.empty())
70  {
73  }
74  logger(DEBUG, "ParticleHandler::ParticleHandler(const ParticleHandler &PH) finished");
75 }
76 
84 {
85  if (this != &rhs)
86  {
87  clear();
88  largestParticle_ = nullptr;
89  smallestParticle_ = nullptr;
91  if (!objects_.empty())
92  {
95  }
96  }
97  logger(DEBUG, "ParticleHandler::operator = (const ParticleHandler& rhs) finished");
98  return *this;
99 }
100 
106 {
107  //First reset the pointers, such that they are not checked twice when removing particles
108  largestParticle_ = nullptr;
109  smallestParticle_ = nullptr;
110  logger(DEBUG, "ParticleHandler::~ParticleHandler() finished");
111 }
112 
130 {
131  if (P->getSpecies() == nullptr)
132  {
133  logger(WARN, "WARNING: The particle with ID % that is added in ParticleHandler::addObject does not have a "
134  "species yet. Please make sure that you have "
135  "set the species somewhere in the driver code.", P->getId());
136  }
137 
138  //Puts the particle in the Particle list
140  if (getDPMBase() != nullptr)
141  {
142  //This places the particle in this grid
144  //This computes where the particle currently is in the grid
146  }
147  //set the particleHandler pointer
148  P->setHandler(this);
149  //compute mass of the particle
150  P->getSpecies()->computeMass(P);
151  //Check if this particle has new extrema
152  checkExtrema(P);
153 }
154 
155 
170 {
171  if (P->getSpecies() == nullptr)
172  {
173  logger(WARN, "WARNING: The particle with ID % that is added in "
174  "ParticleHandler::addObject does not have a species yet. "
175  "Please make sure that you have "
176  "set the species somewhere in the driver code.", P->getId());
177  }
178 #ifdef MERCURY_USE_MPI
179  bool insertParticle;
180  //Check if the particle P should be added to the current domain
181  if (NUMBER_OF_PROCESSORS == 1)
182  {
183  insertParticle = true;
184  }
185  else
186  {
187  insertParticle = getDPMBase()->mpiInsertParticleCheck(P);
188  }
189 
190  //Add the particle if it is in the mpi domain or if the domain is not yet defined
191  if(insertParticle)
192  {
193 #endif
194  //Puts the particle in the Particle list
196  if (getDPMBase() != nullptr)
197  {
198  //This places the particle in this grid
200  //This computes where the particle currently is in the grid
202 
203  // Broadcasts the existance of a new particle
205  }
206  //set the particleHandler pointer
207  P->setHandler(this);
208  //compute mass of the particle
209  P->getSpecies()->computeMass(P);
210  //Check if this particle has new extrema
211  checkExtrema(P);
212  if (!P->isSphericalParticle())
213  {
214  getDPMBase()->setRotation(true);
215  }
216 #ifdef MERCURY_USE_MPI
217  P->setPeriodicComplexity(std::vector<int>(0));
218  }
219  else
220  {
221  logger.assert_debug(!P->isMPIParticle(),"Can't add mpi particle as it does not exist");
222  logger.assert_debug(!P->isPeriodicGhostParticle(),"Can't add mpi particle as it does not exist");
223  //Somehwere a really new particle has been added, so to keep the ID's globally unique, we also update
224  //the unique value on all other processors
226  }
227 
228  //Add the particle to the ghost particle lists
230  //Check if this new particle requires an update in the mpi grid (interactionDistance).
232 
233  //Delete the particle that was supposed to be added (but was not)
234  if(!insertParticle)
235  {
236  delete P;
237  }
238 #endif
239 }
240 
249 void ParticleHandler::addObject(int fromProcessor, BaseParticle* p)
250 {
251 #ifdef MERCURY_USE_MPI
252  MPIContainer& communicator = MPIContainer::Instance();
253 
254  //The processor that contains the particle that needs to be copied needs to identify the target, and communicate this
255  MPIParticle pInfo;
256  if (communicator.getProcessorID() == fromProcessor)
257  {
259  }
260 
261  //Broadcast from processor i
262  communicator.broadcast(&pInfo,MercuryMPIType::PARTICLE,fromProcessor);
263  copyDataFromMPIParticleToParticle(&pInfo, p, this);
264 
265  //All processors now know have the same information and we can proceed with the collective add
266  addObject(p);
267 #else
268  logger(WARN, "Function addObject(int fromProcessor, BaseParticle* p) should not be used in serial code");
269 #endif
270 }
271 
282 void ParticleHandler::addGhostObject(int fromProcessor, int toProcessor, BaseParticle* p)
283 {
284 #ifdef MERCURY_USE_MPI
285  MPIContainer& communicator = MPIContainer::Instance();
286  if (fromProcessor == toProcessor)
287  {
288  if (communicator.getProcessorID() == fromProcessor)
289  {
290  addGhostObject(p);
291  }
292 
293  return;
294  }
295 
296  //Create communication information object
298  MPIParticle pInfo;
299  int tag;
300  if (communicator.getProcessorID() == fromProcessor)
301  {
303  tag = fromProcessor*MAX_PROC*10 + toProcessor*10 + MercuryMPITag::PARTICLE_DATA;
304  communicator.send(&pInfo, MercuryMPIType::PARTICLE, 1, toProcessor, tag);
305  }
306 
307  if (communicator.getProcessorID() == toProcessor)
308  {
309  tag = fromProcessor*MAX_PROC*10 + toProcessor*10 + MercuryMPITag::PARTICLE_DATA;
310  communicator.receive(&pInfo, MercuryMPIType::PARTICLE, 1, fromProcessor, tag);
311  }
312 
313  //Sync the communication
314  communicator.sync();
315 
316  //Add the data to the new particle
317  if (communicator.getProcessorID() == toProcessor)
318  {
319  copyDataFromMPIParticleToParticle(&pInfo, p, this);
320  }
321 
322 
323  //Only toProcessor adds the particle, quietly without any check with other processors
324  //IMPORTANT NOTE: When adding a periodic particle in parallel, this is performed just before
325  //finding new particles. For that reason we dont have to add a ghostparticle to the mpi communication lists
326  if (communicator.getProcessorID() == toProcessor)
327  {
328  addGhostObject(p);
329  }
330 #else
331  logger(WARN,
332  "Function addGhostObject(int fromProcessor, int toProcessor, BaseParticle* p) should not be used in serial code");
333 #endif
334 }
335 
344 {
345 #ifdef MERCURY_USE_MPI
346  if (P->getSpecies() == nullptr)
347  {
348  logger(WARN, "[ParticleHandler::adGhostObject(BaseParticle*)] "
349  "WARNING: The particle with ID % that is added in "
350  "ParticleHandler::addObject does not have a species yet."
351  " Please make sure that you have "
352  "set the species somewhere in the driver code.", P->getId());
353  }
354 
355  MPIContainer& communicator = MPIContainer::Instance();
356  //Puts the particle in the Particle list
358  if (getDPMBase() != nullptr)
359  {
360  //This places the particle in this grid
362  //This computes where the particle currently is in the grid
364 
365  // broadcast particle addition
367  }
368  //set the particleHandler pointer
369  P->setHandler(this);
370  //compute mass of the particle
371  P->getSpecies()->computeMass(P) ;
372  //Check if this particle has new extrema
373  checkExtrema(P);
374 #else
375  logger(INFO,
376  "Function ParticleHandler::mpiAddObject(BaseParticle* P) should only be called when compiling with parallel build on");
377 #endif
378 }
379 
389 void ParticleHandler::removeObject(const unsigned int index)
390 {
391 #ifdef MERCURY_USE_MPI
392  MPIContainer& communicator = MPIContainer::Instance();
393  if (communicator.getNumberOfProcessors() > 1 )
394  {
395  logger(WARN, "[ParticleHandler::removeObject(const unsigned int)] Using the function removeObject in parallel could lead to out-of-sync communication, Instead use deletion boundaries");
396  }
397 #endif
398 #ifdef CONTACT_LIST_HGRID
399  getDPMBase()->getPossibleContactList().remove_ParticlePosibleContacts(getObject(id));
400 #endif
402  // broadcast the particle removal
403  getDPMBase()->handleParticleRemoval(getObject(index)->getId());
405 }
406 
417 void ParticleHandler::removeGhostObject(const unsigned int index)
418 {
419 #ifdef MERCURY_USE_MPI
420 #ifdef CONTACT_LIST_HGRID
421  getDPMBase()->getPossibleContactList().remove_ParticlePosibleContacts(getObject(id));
422 #endif
424  // broadcast the particle removal
425  getDPMBase()->handleParticleRemoval(getObject(index)->getId());
427 #else
428  logger(ERROR, "This function should only be used interally for mpi routines");
429 #endif
430 }
431 
437 {
438 #ifdef CONTACT_LIST_HGRID
439  getDPMBase()->getPossibleContactList().remove_ParticlePosibleContacts(getLastObject());
440 #endif
443 }
444 
446 {
447  if (getSize() == 0)
448  {
449  logger(DEBUG, "No particles, so cannot compute the smallest particle.");
450 
451  smallestParticle_ = nullptr;
452  return;
453  }
454  Mdouble min = std::numeric_limits<Mdouble>::max();
455  smallestParticle_ = nullptr;
456  for (BaseParticle* const particle : objects_)
457  {
458  //if (!(particle->isMPIParticle() || particle->isPeriodicGhostParticle()))
459  {
460  if (particle->getMaxInteractionRadius() < min)
461  {
462  min = particle->getMaxInteractionRadius();
463  smallestParticle_ = particle;
464  }
465  }
466  }
467 }
468 
470 {
471  if (getSize() == 0)
472  {
473  logger(DEBUG, "No particles, so cannot compute the largest particle.");
474  largestParticle_ = nullptr;
475  return;
476  }
477  Mdouble max = -std::numeric_limits<Mdouble>::max();
478  largestParticle_ = nullptr;
479  for (BaseParticle* const particle : objects_)
480  {
481  //if (!(particle->isMPIParticle() || particle->isPeriodicGhostParticle()))
482  {
483  if (particle->getMaxInteractionRadius() > max)
484  {
485  max = particle->getMaxInteractionRadius();
486  largestParticle_ = particle;
487  }
488  }
489  }
490 }
491 
497 {
498  return smallestParticle_;
499 }
500 
506 {
507 #ifdef MERCURY_USE_MPI
508  logger(WARN,"getSmallestParticle should not be used in parallel; use getSmallestInteractionRadius or "
509  "ParticleSpecies::getSmallestParticleMass instead");
510  return nullptr;
511 #else
512  return getSmallestParticleLocal();
513 #endif
514 }
515 
521 {
522  return largestParticle_;
523 }
524 
529 {
530 #ifdef MERCURY_USE_MPIO
531  logger(WARN,"getLargestParticle() should not be used in parallel; use getLargestInteractionRadius instead");
532  return nullptr;
533 #else
534  return getLargestParticleLocal();
535 #endif
536 }
537 
539 {
540  Mdouble ene = 0;
541  for (auto p : *this)
542  {
543  if (!(p->isMPIParticle() || p->isPeriodicGhostParticle()))
544  {
545  ene += p->getKineticEnergy();
546  }
547  }
548  return ene;
549 }
550 
552 {
553 #ifdef MERCURY_USE_MPI
554  Mdouble kineticEnergyLocal = getKineticEnergyLocal();
555  Mdouble kineticEnergyGlobal = 0.0;
556 
557  //sum up over all domains
558  MPIContainer& communicator = MPIContainer::Instance();
559  communicator.allReduce(kineticEnergyLocal, kineticEnergyGlobal, MPI_SUM);
560 
561  return kineticEnergyGlobal;
562 #else
563  return getKineticEnergyLocal();
564 #endif
565 }
566 
568 {
569  Mdouble ene = 0;
570  for (auto p : *this)
571  {
572  if (!(p->isMPIParticle() || p->isPeriodicGhostParticle()))
573  {
574  ene += p->getRotationalEnergy();
575  }
576  }
577  return ene;
578 }
579 
581 {
582 #ifdef MERCURY_USE_MPI
583  Mdouble rotationalEnergyLocal = getRotationalEnergyLocal();
584  Mdouble rotationalEnergyGlobal = 0.0;
585 
586  //sum up over all domains
587  MPIContainer& communicator = MPIContainer::Instance();
588  communicator.allReduce(rotationalEnergyLocal, rotationalEnergyGlobal, MPI_SUM);
589 
590  return rotationalEnergyGlobal;
591 #else
592  return getRotationalEnergyLocal();
593 #endif
594 }
595 
597 {
598  Mdouble m = 0;
599  for (auto p : *this)
600  if (!(p->isFixed() || p->isMPIParticle() || p->isPeriodicGhostParticle()))
601  m += p->getMass();
602  return m;
603 }
604 
606 {
607 #ifdef MERCURY_USE_MPI
608  Mdouble massLocal = getMassLocal();
609  Mdouble massGlobal = 0.0;
610 
611  //Sum up over all domains
612  MPIContainer& communicator = MPIContainer::Instance();
613  communicator.allReduce(massLocal, massGlobal, MPI_SUM);
614 
615  return massGlobal;
616 #else
617  return getMassLocal();
618 #endif
619 }
620 
622 {
623  Vec3D com = {0, 0, 0};
624  for (auto p : *this)
625  if (!(p->isFixed() || p->isMPIParticle() || p->isPeriodicGhostParticle()))
626  com += p->getMass() * p->getPosition();
627  return com;
628 }
629 
631 {
632 #ifdef MERCURY_USE_MPI
633  Vec3D massTimesPositionLocal = getMassTimesPositionLocal();
634  Vec3D massTimesPositionGlobal = {0.0, 0.0, 0.0};
635 
636  // Sum up over all domains
637  MPIContainer& communicator = MPIContainer::Instance();
638  communicator.allReduce(massTimesPositionLocal.X, massTimesPositionGlobal.X, MPI_SUM);
639  communicator.allReduce(massTimesPositionLocal.Y, massTimesPositionGlobal.Y, MPI_SUM);
640  communicator.allReduce(massTimesPositionLocal.Z, massTimesPositionGlobal.Z, MPI_SUM);
641 
642  return massTimesPositionGlobal;
643 #else
644  return getMassTimesPositionLocal();
645 #endif
646 }
647 
649 {
650  Mdouble m = getMass();
651  if (m == 0)
652  {
653  Vec3D nanvec = {constants::NaN, constants::NaN, constants::NaN};
654  return nanvec;
655  }
656  else
657  return getMassTimesPosition() / m;
658 }
659 
661 {
662  Vec3D momentum = {0, 0, 0};
663  for (auto p : *this)
664  if (!(p->isFixed() || p->isMPIParticle() || p->isPeriodicGhostParticle()))
665  momentum += p->getMomentum();
666  return getMPISum(momentum);
667 }
668 
670 {
671  Vec3D momentum = {0, 0, 0};
672  for (auto p : *this)
673  if (!(p->isFixed() || p->isMPIParticle() || p->isPeriodicGhostParticle()))
674  momentum += p->getAngularMomentum();
675  return getMPISum(momentum);
676 }
677 
682 {
683  if (getSize() == 0)
684  {
685  logger(WARN, "No particles to set getFastestParticle()");
686  return nullptr;
687  }
688  BaseParticle* p = nullptr;
689  Mdouble maxSpeed = -std::numeric_limits<Mdouble>::max();
690  for (BaseParticle* const pLoop : objects_)
691  {
692  if (!(pLoop->isMPIParticle() || pLoop->isPeriodicGhostParticle()))
693  {
694  if ((pLoop->getVelocity().getLength()) > maxSpeed)
695  {
696  maxSpeed = pLoop->getVelocity().getLength();
697  p = pLoop;
698  }
699  }
700  }
701  return p;
702 }
703 
705 {
706 #ifdef MERCURY_USE_MPI
707  logger(ERROR,"This function should not be used in parallel");
708 #endif
709  return getFastestParticleLocal();
710 }
711 
712 /*
713  * \returns The smallest interaction radius on a local domain
714  */
716 {
717  if (!(getSmallestParticleLocal() == nullptr))
718  {
720  }
721  else
722  {
723  return 0.0;
724  }
725 }
726 
727 /*
728  * \returns The smallest interaction radius
729  */
731 {
732 #ifdef MERCURY_USE_MPI
733  //Compute the local value
734  Mdouble smallestInteractionRadiusLocal = getSmallestInteractionRadiusLocal();
735  Mdouble smallestInteractionRadiusGlobal = 0.0;
736 
737  //Obtain the global value
738  MPIContainer& communicator = MPIContainer::Instance();
739  communicator.allReduce(smallestInteractionRadiusLocal, smallestInteractionRadiusGlobal, MPI_MIN);
740 
741  return smallestInteractionRadiusGlobal;
742 
743 #else
745 #endif
746 }
747 
748 /*
749  * \returns the largest interaction radius on a local domain
750  */
752 {
753  if (!(getLargestParticleLocal() == nullptr))
754  {
756  }
757  else
758  {
759  return 0.0;
760  }
761 }
762 
763 /*
764  * \returns the largest interaction radius
765  */
767 {
768 #ifdef MERCURY_USE_MPI
769  Mdouble largestInteractionRadiusLocal = getLargestInteractionRadiusLocal();
770  Mdouble largestInteractionRadiusGlobal = 0.0;
771 
772  //Obtain the global value
773  MPIContainer& communicator = MPIContainer::Instance();
774  communicator.allReduce(largestInteractionRadiusLocal, largestInteractionRadiusGlobal, MPI_MAX);
775 
776  return largestInteractionRadiusGlobal;
777 #else
779 #endif
780 }
781 
786 {
787  Mdouble sumRadius = 0;
788  for (BaseParticle* const p : objects_)
789  {
790  if (!(p->isMPIParticle() || p->isPeriodicGhostParticle()))
791  {
792  sumRadius += p->getRadius();
793  }
794  }
795  return sumRadius;
796 }
797 
799 {
800 #ifdef MERCURY_USE_MPI
801  Mdouble sumRadiusLocal = getSumRadiusLocal();
802  unsigned numberOfRealParticlesLocal = getNumberOfRealObjectsLocal();
803 
804  Mdouble sumRadiusGlobal = 0.0;
805  unsigned numberOfRealParticlesGlobal = 0;
806 
807  //Sum up over all domains
808  MPIContainer& communicator = MPIContainer::Instance();
809  communicator.allReduce(sumRadiusLocal, sumRadiusGlobal, MPI_SUM);
810  communicator.allReduce(numberOfRealParticlesLocal, numberOfRealParticlesGlobal, MPI_SUM);
811 
812  return sumRadiusGlobal/numberOfRealParticlesGlobal;
813 #else
814  return getSumRadiusLocal() / getSize();
815 #endif
816 }
817 
824 {
825 #ifdef MERCURY_USE_MPI
826  logger(ERROR,"getLowestPositionComponentParticle() not implemented yet in parallel");
827 #endif
828  if (getSize() == 0)
829  {
830  logger(WARN, "No getLowestPositionComponentParticle(const int i) since there are no particles.");
831  return nullptr;
832  }
833  BaseParticle* p = nullptr;
834  Mdouble min = std::numeric_limits<Mdouble>::max();
835  for (BaseParticle* const pLoop : objects_)
836  {
837  if (!(pLoop->isMPIParticle() || pLoop->isPeriodicGhostParticle()))
838  {
839  if (pLoop->getPosition().getComponent(i) < min)
840  {
841  min = pLoop->getPosition().getComponent(i);
842  p = pLoop;
843  }
844  }
845  }
846  return p;
847 }
848 
850 {
851 #ifdef MERCURY_USE_MPI
852  logger(ERROR,"This function should not be used in parallel");
853 #endif
855 }
856 
863 {
864  if (getSize() == 0)
865  {
866  logger(WARN, "No getHighestPositionComponentParticle(const int i) since there are no particles.");
867  return nullptr;
868  }
869  BaseParticle* p = nullptr;
870  Mdouble max = -std::numeric_limits<Mdouble>::max();
871  for (BaseParticle* const pLoop : objects_)
872  {
873  if (!(pLoop->isMPIParticle() || pLoop->isPeriodicGhostParticle()))
874  {
875  if (pLoop->getPosition().getComponent(i) > max)
876  {
877  max = pLoop->getPosition().getComponent(i);
878  p = pLoop;
879  }
880  }
881  }
882 
883  return p;
884 }
885 
887 {
888 #ifdef MERCURY_USE_MPI
889  logger(ERROR,"This function should not be used in parallel");
890 #endif
892 }
893 
900 {
901  if (getSize() == 0)
902  {
903  logger(WARN, "No getLowestVelocityComponentParticle(const int i) since there are no particles");
904  return nullptr;
905  }
906  BaseParticle* p = nullptr;
907  Mdouble min = std::numeric_limits<Mdouble>::max();
908  for (BaseParticle* const pLoop : objects_)
909  {
910  if (!(pLoop->isMPIParticle() || pLoop->isPeriodicGhostParticle()))
911  {
912  if (pLoop->getVelocity().getComponent(i) < min)
913  {
914  min = pLoop->getVelocity().getComponent(i);
915  p = pLoop;
916  }
917  }
918  }
919  return p;
920 }
921 
923 {
924 #ifdef MERCURY_USE_MPI
925  logger(ERROR,"This function should not be used in parallel");
926 #endif
928 }
929 
936 {
937  if (!getSize())
938  {
939  logger(WARN, "No getHighestVelocityComponentParticle(const int i) since there are no particles");
940  return nullptr;
941  }
942  BaseParticle* p = nullptr;
943  Mdouble max = -std::numeric_limits<Mdouble>::max();
944  for (BaseParticle* const pLoop : objects_)
945  {
946  if (!(pLoop->isMPIParticle() || pLoop->isPeriodicGhostParticle()))
947  {
948  if (pLoop->getVelocity().getComponent(i) > max)
949  {
950  max = pLoop->getVelocity().getComponent(i);
951  p = pLoop;
952  }
953  }
954  }
955  return p;
956 }
957 
959 {
960 #ifdef MERCURY_USE_MPI
961  logger(ERROR,"This function should not be used in parallel");
962 #endif
964 }
965 
972 {
973  smallestParticle_ = nullptr;
974  largestParticle_ = nullptr;
976 }
977 
984 {
985  //Define the attribute function
986  std::function<Mdouble(BaseParticle*)> particleAttribute = [](BaseParticle* p) { return p->getPosition().X; };
987 
988  //Obtain the MAX attribute
989  Mdouble positionXGlobal = getParticleAttribute<Mdouble>(particleAttribute, AttributeType::MAX);
990 
991  return positionXGlobal;
992 }
993 
998 {
999  return NFixedParticles_;
1000 }
1001 
1006 {
1008 }
1009 
1014 {
1015  if (type == "BaseParticle") {
1016  //for backwards compatibility
1017  return new SphericalParticle;
1018  }
1019  else if (type == "SphericalParticle")
1020  {
1021  return new SphericalParticle;
1022  }
1023  else if (type == "LiquidFilmParticle")
1024  {
1025  return new LiquidFilmParticle;
1026  }
1027  else if (type == "SuperQuadricParticle")
1028  {
1029  return new SuperQuadricParticle;
1030  }
1031  else if (type == "ThermalParticle")
1032  {
1033  return new ThermalParticle;
1034  }
1035  else
1036  {
1037  logger(WARN, "Particle type % not understood in restart file. Particle will not be read.", type);
1038  return nullptr;
1039  }
1040 }
1041 
1046 {
1047  std::string type;
1048  BaseParticle* particle = nullptr;
1049  if (getStorageCapacity() > 0)
1050  {
1051  is >> type;
1052  logger(DEBUG, "ParticleHandler::readAndCreateObject(is): type %", type);
1053  particle = createObject(type);
1054  particle->setHandler(this);
1055  particle->read(is);
1056  }
1057  else //for insertion boundaries
1058  {
1059  is >> type;
1060  logger(DEBUG, "ParticleHandler::readAndCreateObject(is): type %", type);
1061  particle = createObject(type);
1062  particle->setHandler(this);
1063  particle->read(is);
1064  }
1065  particle->setSpecies(getDPMBase()->speciesHandler.getObject(particle->getIndSpecies()));
1066  return particle;
1067 }
1068 
1073 {
1075  addExistingObject(o);
1076 }
1077 
1079 // * \param[in] type The first value of the position.
1080 // * \param[in] is The input stream from which the information is read.
1081 // * \details The old objects did not have their type in the beginning of the line.
1082 // * Instead, the first string of the file was the position in x-direction.
1083 // * Since we already read the first string of the file, we need to give
1084 // * it to this function and convert it to the position in x-direction.
1085 // * The rest of the stream is then read in the usual way.
1086 // */
1087 //void ParticleHandler::readAndCreateOldObject(std::istream& is, const std::string& type)
1088 //{
1089 // //read in next line
1090 // std::stringstream line;
1091 // helpers::getLineFromStringStream(is, line);
1092 // logger(VERBOSE, line.str());
1093 // //std::cout << line.str() << std::endl;
1094 //
1095 // BaseParticle particle;
1096 //
1097 // //Declare all properties of the particle
1098 // unsigned int indSpecies;
1099 // Mdouble radius, inverseMass, inverseInertia;
1100 // Vec3D position, velocity, euler, angularVelocity;
1101 //
1102 // //Read all values
1103 // position.X = atof(type.c_str());
1104 //
1105 // line >> position.Y >> position.Z >> velocity >> radius >> euler >> angularVelocity >> inverseMass >> inverseInertia >> indSpecies;
1106 //
1107 // //Put the values in the particle
1108 // particle.setSpecies(getDPMBase()->speciesHandler.getObject(indSpecies));
1109 // particle.setPosition(position);
1110 // particle.setVelocity(velocity);
1111 // particle.setRadius(radius);
1112 // Quaternion q; q.setEuler(euler);
1113 // particle.setOrientation(q);
1114 // particle.setAngularVelocity(angularVelocity);
1115 // if (inverseMass == 0.0)
1116 // particle.fixParticle();
1117 // else
1118 // {
1119 // particle.setInverseInertia(MatrixSymmetric3D(1,0,0,1,0,1)*inverseInertia);
1120 // }
1121 //
1122 // //Put the particle in the handler
1123 // copyAndAddObject(particle);
1124 //}
1125 
1126 void ParticleHandler::write(std::ostream& os) const
1127 {
1128 #ifdef MERCURY_USE_MPI
1129  os << "Particles " << getNumberOfRealObjectsLocal() << std::endl;
1130  for (BaseParticle* it : *this)
1131  {
1132  if (!it->isPeriodicGhostParticle() && !it->isMPIParticle())
1133  {
1134  os << (*it) << '\n';
1135  }
1136  }
1137 #else
1138  os << "Particles " << getSize() << '\n';
1139  for (BaseParticle* it : *this)
1140  {
1141  os << (*it) << '\n';
1142  }
1143 #endif
1144  //os.flush();
1145 }
1146 
1152 {
1153  if (P == largestParticle_)
1154  {
1155  //if the properties of the largest particle changes
1157  }
1159  {
1160  largestParticle_ = P;
1161  }
1162 
1163  if (P == smallestParticle_)
1164  {
1165  //if the properties of the smallest particle changes
1167  }
1169  {
1170  smallestParticle_ = P;
1171  }
1172 }
1173 
1178 {
1179  if (P == largestParticle_)
1180  {
1182  }
1183  if (P == smallestParticle_)
1184  {
1186  }
1187 }
1188 
1193 void ParticleHandler::computeAllMasses(unsigned int indSpecies)
1194 {
1195  for (BaseParticle* particle : objects_)
1196  {
1197  if (particle->getIndSpecies() == indSpecies)
1198  {
1199  particle->getSpecies()->computeMass(particle);
1200  }
1201  }
1202 }
1203 
1205 {
1206  for (BaseParticle* particle : objects_)
1207  {
1208  particle->getSpecies()->computeMass(particle);
1209  }
1210 }
1211 
1213 {
1214  NFixedParticles_++;
1215 }
1216 
1218 {
1219  NFixedParticles_--;
1220 }
1221 
1225 std::string ParticleHandler::getName() const
1226 {
1227  return "ParticleHandler";
1228 }
1229 
1231 {
1232  Mdouble volume = 0;
1233  for (auto p : *this)
1234  {
1235  if (!(p->isFixed() || p->isPeriodicGhostParticle() || p->isMPIParticle()))
1236  volume += p->getVolume();
1237  }
1238  return volume;
1239 }
1240 
1242 {
1243 #ifdef MERCURY_USE_MPI
1244  Mdouble volumeLocal = getVolumeLocal();
1245  Mdouble volumeGlobal = 0.0;
1246 
1247  // sum up over all domains
1248  MPIContainer& communicator = MPIContainer::Instance();
1249  communicator.allReduce(volumeLocal, volumeGlobal, MPI_SUM);
1250 
1251  return volumeGlobal;
1252 #else
1253  return getVolumeLocal();
1254 #endif
1255 }
1256 
1262 {
1263 #ifdef MERCURY_USE_MPI
1264  const MPIContainer& communicator = MPIContainer::Instance();
1265  if (communicator.getNumberOfProcessors() > 1)
1266  {
1267  unsigned int numberOfFakeParticles = 0;
1268  numberOfFakeParticles += getDPMBase()->domainHandler.getCurrentDomain()->getNumberOfTrueMPIParticles();
1270  logger.assert_debug(numberOfFakeParticles <= getSize(), "More fake particles than getSize()");
1271  return (getSize() - numberOfFakeParticles);
1272  }
1273  else
1274  {
1275  return getSize();
1276  }
1277 #else
1278  return getSize();
1279 #endif
1280 }
1281 
1283 {
1284 #ifdef MERCURY_USE_MPI
1285  MPIContainer& communicator = MPIContainer::Instance();
1286  unsigned int numberOfRealParticles = getNumberOfRealObjectsLocal();
1287 
1288  // \todo MX: use an allreduce here
1289  //Combine the total number of Particles into one number on processor 0
1290  communicator.reduce(numberOfRealParticles, MPI_SUM);
1291 
1292  //Broadcast new number to all the processorsi
1293  communicator.broadcast(numberOfRealParticles);
1294  return numberOfRealParticles;
1295 #else
1296  return getSize();
1297 #endif
1298 }
1299 
1300 /*
1301  * \returns the size of the ParticleHandler. In parallel code it throws an error to avoid confusion with real and fake particles
1302  */
1304 {
1305 //#ifdef MERCURY_USE_MPI
1306 // MPIContainer& communicator = MPIContainer::Instance();
1307 // if (communicator.getNumberOfProcessors() > 1)
1308 // {
1309 // logger(WARN,"When compiling with MPI please do not use getNumberOfObjects(). Instead use: getNumberOfRealObjectsLocal(), getNumberOfRealObjects() or getSize()");
1310 // }
1311 //#endif
1312  return getSize();
1313 }
1314 
1321 {
1322  unsigned int numberOfFixedParticles = 0;
1323  for (BaseParticle* particle : *this)
1324  {
1325  if (particle->isFixed())
1326  {
1327  numberOfFixedParticles++;
1328  }
1329  }
1330  return numberOfFixedParticles;
1331 }
1332 
1337 {
1338 #ifdef MERCURY_USE_MPI
1339  unsigned int numberOfFixedParticlesLocal = getNumberOfFixedObjectsLocal();
1340  unsigned int numberOfFixedParticles = 0;
1341 
1342  MPIContainer& communicator = MPIContainer::Instance();
1343  communicator.allReduce(numberOfFixedParticlesLocal, numberOfFixedParticles, MPI_SUM);
1344  return numberOfFixedParticles;
1345 #else
1347 #endif
1348 }
1349 
1351 {
1352  for (auto i: *this)
1353  {
1354  i->actionsAfterTimeStep();
1355  }
1356 }
1357 
1359  double liquidVolume = 0;
1360  for (auto i : objects_) {
1361  auto j = dynamic_cast<LiquidFilmParticle*>(i);
1362  if (j and !j->isMPIParticle()) liquidVolume += j->getLiquidVolume();
1363  }
1364  return getMPISum(liquidVolume);
1365 };
Mdouble getVolume() const
This class contains all information and functions required for communication between processors...
Definition: MpiContainer.h:129
void read(std::istream &is) override
Particle read function, which accepts an std::istream as input.
void copyDataFromParticleToMPIParticle(BaseParticle *p)
BaseParticle * getLargestParticle() const
Returns the pointer of the largest particle in the particle handler. When mercury is running in paral...
virtual void handleParticleRemoval(unsigned int id)
Handles the removal of particles from the particleHandler.
Definition: DPMBase.cc:5249
unsigned int getId() const
Returns the unique identifier of any particular object.
Definition: BaseObject.h:125
virtual void hGridUpdateParticle(BaseParticle *obj UNUSED)
Definition: DPMBase.cc:1688
void checkExtremaOnDelete(BaseParticle *P)
Checks if the extrema of this ParticleHandler needs updating when a particle is deleted.
void write(std::ostream &os) const
virtual void hGridInsertParticle(BaseParticle *obj UNUSED)
Definition: DPMBase.cc:1681
static MPIContainer & Instance()
fetch the instance to be used for communication
Definition: MpiContainer.h:134
Mdouble X
the vector components
Definition: Vector.h:65
void removeObject(unsigned int index) override
Removes a BaseParticle from the ParticleHandler.
const Mdouble NaN
Definition: GeneralDefine.h:43
A spherical particle is the most simple particle used in MercuryDPM.
void increaseId()
Definition: BaseHandler.h:251
void addObject(BaseParticle *P) override
Adds a BaseParticle to the ParticleHandler.
unsigned int getSize() const
Gets the size of the particleHandler (including mpi and periodic particles)
virtual bool isSphericalParticle() const
Definition: BaseParticle.h:642
bool mpiInsertParticleCheck(BaseParticle *P)
Function that checks if the mpi particle should really be inserted by the current domain...
Definition: DPMBase.cc:1716
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
Definition of different loggers with certain modules. A user can define its own custom logger here...
double Mdouble
Definition: GeneralDefine.h:34
const ParticleSpecies * getSpecies() const
Returns a pointer to the species of this BaseInteractable.
unsigned int getNumberOfObjects() const override
Returns the number of objects in the container. In parallel code this practice is forbidden to avoid ...
void checkExtrema(BaseParticle *P)
Checks if the extrema of this ParticleHandler needs updating.
std::string getName() const override
Returns the name of the handler, namely the string "ParticleHandler".
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 setDPMBase(DPMBase *DPMBase)
Sets the problem that is solved using this handler.
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
std::enable_if< std::is_scalar< T >::value, void >::type receive(T &t, int from, int tag)
asynchronously receive a scalar from some other processor.
Definition: MpiContainer.h:221
void setHandler(ParticleHandler *handler)
Sets the pointer to the particle's ParticleHandler.
void addGhostObject(int fromPrcessor, int toProcessor, BaseParticle *p)
Adds a ghost particle located at fromProcessor to toProcessor.
virtual void addGhostObject(T *O)
Adds a new Object to the BaseHandler. called by the to avoid increasing the id.
Definition: BaseHandler.h:444
BaseParticle * getHighestVelocityComponentParticleLocal(int i) const
Gets a pointer to the particle with the highest velocity in direction i in this ParticleHandler.
Mdouble getLargestInteractionRadiusLocal() const
Returns the largest interaction radius of the current domain.
Mdouble getSmallestInteractionRadius() const
Returns the smallest interaction radius.
std::size_t getNumberOfProcessors() const
Get the total number of processors participating in this simulation.
Mdouble getSumRadiusLocal() const
void computeMass(BaseParticle *p) const
Compute Particle mass function, which required a reference to the Species vector. It computes the Par...
void computeAllMasses()
Computes the mass for all BaseParticle in this ParticleHandler.
virtual void handleParticleAddition(unsigned int id, BaseParticle *p)
Definition: DPMBase.cc:5264
BaseParticle * getHighestPositionComponentParticle(int i) const
Gets a pointer to the particle with the highest coordinates in direction i in this ParticleHandler...
DomainHandler domainHandler
An object of the class DomainHandler which deals with parallel code.
Definition: DPMBase.h:1401
void setSpecies(const ParticleSpecies *species)
Mdouble getMass() const
unsigned int getNumberOfFixedParticles() const
Gets the number of particles that are fixed.
void clear() override
Empties the whole ParticleHandler by removing all BaseParticle.
BaseParticle * smallestParticle_
A pointer to the smallest BaseParticle (by interactionRadius) in this ParticleHandler.
void removeLastObject()
Removes the last Object from the BaseHandler.
Definition: BaseHandler.h:511
void removedFixedParticle()
Decrement of the number of fixed particles.
void updateGhostGrid(BaseParticle *P)
Checks if the Domain/periodic interaction distance needs to be updated and updates it accordingly...
Definition: DPMBase.cc:1821
unsigned int getNumberOfTrueMPIParticles()
Obtains the number of particles in the particleHandler that are MPIParticles, but NOT periodic partic...
Definition: Domain.cc:1679
unsigned int getNumberOfRealObjects() const
Returns the number of real objects (on all processors)
void copyDataFromMPIParticleToParticle(MPIParticle *bP, BaseParticle *p, ParticleHandler *particleHandler)
Copies data from an MPIParticle class to a BaseParticle and sets the particleHandler and species...
std::size_t getProcessorID()
Reduces a scalar on all processors to one scalar on a target processor.
virtual void addExistingObject(T *O)
Adds an existing object to the BaseHandler without changing the id of the object. ...
Definition: BaseHandler.h:417
BaseParticle * getHighestVelocityComponentParticle(int i) const
Gets a pointer to the particle with the highest velocity in direction i in this ParticleHandler.
#define MAX_PROC
Definition: GeneralDefine.h:51
BaseParticle * getSmallestParticle() const
Gets a pointer to the smallest BaseParticle (by interactionRadius) in this ParticleHandler of the loc...
Vec3D getMassTimesPosition() const
Mdouble getSmallestInteractionRadiusLocal() const
Returns the smallest interaction radius of the current domain.
PeriodicBoundaryHandler periodicBoundaryHandler
Internal handler that deals with periodic boundaries, especially in a parallel build.
Definition: DPMBase.h:1396
unsigned int getNumberOfPeriodicGhostParticles()
Returns the number of particles that are flagged is periodicGhostParticle.
BaseParticle * getFastestParticleLocal() const
Gets a pointer to the fastest BaseParticle in this ParticleHandler.
virtual void removeObject(unsigned const int index)
Removes an Object from the BaseHandler.
Definition: BaseHandler.h:472
void removeLastObject()
Removes the last BaseParticle from the ParticleHandler.
void addedFixedParticle()
Increment of the number of fixed particles.
double getLiquidFilmVolume() const
BaseParticle * getObject(const unsigned int id)
Gets a pointer to the Object at the specified index in the BaseHandler.
Data class to send a particle over MPI.
Definition: MpiDataClass.h:81
void readAndAddObject(std::istream &is) override
Create a new particle, based on the information from old-style restart data.
ParticleHandler()
Default constructor, it creates an empty ParticleHandler.
BaseParticle * largestParticle_
A pointer to the largest BaseParticle (by interactionRadius) in this ParticleHandler.
BaseParticle * getLowestVelocityComponentParticleLocal(int i) const
Gets a pointer to the particle with the lowest velocity in direction i in this ParticleHandler.
Vec3D getMPISum(Vec3D &val)
Sums the values over all processors using MPI_reduce.
Mdouble getKineticEnergy() const
std::vector< BaseParticle * > objects_
The actual list of Object pointers.
Definition: BaseHandler.h:302
BaseParticle * getHighestPositionComponentParticleLocal(int i) const
Gets a pointer to the particle with the highest coordinates in direction i in this ParticleHandler...
Container to store the pointers to all objects that one creates in a simulation.
Definition: BaseHandler.h:50
unsigned int getNumberOfRealObjectsLocal() const
Returns the number of real objects on a local domain. MPI particles and periodic particles are neglec...
void removeGhostObject(unsigned int index)
Removes a BaseParticle from the ParticleHandler without a global check, this is only to be done for m...
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
unsigned int getNumberOfFixedObjects() const
Computes the number of fixed particles in the whole simulation.
unsigned int getNumberOfFixedObjectsLocal() const
Computes the number of Fixed particles on a local domain.
Mdouble getHighestPositionX() const
Function returns the highest position in the x-direction.
unsigned int NFixedParticles_
Number of fixed particles.
Vec3D getAngularMomentum() const
bool isPeriodicGhostParticle() const
Indicates if this particle is a ghost in the periodic boundary.
Container to store all BaseParticle.
#define NUMBER_OF_PROCESSORS
For the MPI communication routines this quantity is often required. defining this macro makes the cod...
Definition: GeneralDefine.h:62
Mdouble Y
Definition: Vector.h:65
BaseParticle * getSmallestParticleLocal() const
Gets a pointer to the smallest BaseParticle (by interactionRadius) in this ParticleHandler of the loc...
std::enable_if< std::is_scalar< T >::value, void >::type broadcast(T &t, int fromProcessor=0)
Broadcasts a scalar from the root to all other processors.
Definition: MpiContainer.h:441
Vec3D getMassTimesPositionLocal() const
virtual void hGridRemoveParticle(BaseParticle *obj UNUSED)
Definition: DPMBase.cc:1695
Mdouble getLiquidVolume() const
void sync()
Process all pending asynchronous communication requests before continuing.
Definition: MpiContainer.h:152
std::enable_if< std::is_scalar< T >::value, void >::type send(T &t, int to, int tag)
Asynchronously send a scalar to some other processor.
Definition: MpiContainer.h:171
~ParticleHandler() override
Destructor, it destructs the ParticleHandler and all BaseParticle it contains.
Mdouble getKineticEnergyLocal() const
void computeSmallestParticle()
Computes the smallest particle (by interaction radius) and sets it in smallestParticle_.
Domain * getCurrentDomain()
Gets the domain assigned to the processor.
BaseParticle * getLowestPositionComponentParticle(int i) const
Gets a pointer to the particle with the lowest coordinates in direction i in this ParticleHandler...
void setPeriodicComplexity(std::vector< int > complexity)
Set the periodic communication complexity of the particle.
BaseParticle * getLowestVelocityComponentParticle(int i) const
Gets a pointer to the particle with the lowest velocity in direction i in this ParticleHandler.
void insertGhostParticle(BaseParticle *P)
This function inserts a particle in the mpi communication boundaries.
Definition: DPMBase.cc:1795
bool isMPIParticle() const
Indicates if this particle is a ghost in the MPI domain.
void computeLargestParticle()
Computes the largest particle (by interaction radius) and sets it in largestParticle_.
Mdouble getMeanRadius() const
virtual void addObject(T *object)
Adds a new Object to the BaseHandler.
Definition: BaseHandler.h:431
Mdouble getRotationalEnergyLocal() const
BaseParticle * getLowestPositionComponentParticleLocal(int i) const
Gets a pointer to the particle with the lowest coordinates in direction i in this ParticleHandler...
void setRotation(bool rotation)
Sets whether particle rotation is enabled or disabled.
Definition: DPMBase.h:550
Mdouble getMassLocal() const
Vec3D getMomentum() const
unsigned int getIndSpecies() const
Returns the index of the species associated with the interactable object.
BaseParticle * getLastObject()
Gets a pointer to the last Object in this BaseHandler.
static BaseParticle * createObject(const std::string &type)
Reads BaseParticle into the ParticleHandler from restart data.
unsigned int getStorageCapacity() const
Gets the storage capacity of this BaseHandler.
void copyContentsFromOtherHandler(const BaseHandler< BaseParticle > &BH)
Function that copies the contents (vector of pointers, maxObject_, nextId_, DPMBase_) from one handle...
Definition: Vector.h:49
Mdouble getVolumeLocal() const
DPMBase * getDPMBase()
Gets the problem that is solved using this handler.
Definition: BaseHandler.h:725
ParticleHandler & operator=(const ParticleHandler &rhs)
Assignment operator.
unsigned int getNumberOfUnfixedParticles() const
Gets the number of particles that are not fixed.
BaseParticle * getLargestParticleLocal() const
Gets a pointer to the largest BaseParticle (by interactionRadius) in the ParticleHandler of the local...
Mdouble getRotationalEnergy() const
Vec3D getCentreOfMass() const
Mdouble Z
Definition: Vector.h:65
void addExistingObject(BaseParticle *P) override
Adds a BaseParticle to the ParticleHandler.
Mdouble getLargestInteractionRadius() const
Returns the largest interaction radius.
virtual void clear()
Empties the whole BaseHandler by removing all Objects and setting all other variables to 0...
Definition: BaseHandler.h:528
BaseParticle * readAndCreateObject(std::istream &is)
Create a new particle, based on the information provided in a restart file.
BaseParticle * getFastestParticle() const