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 "
134  "ParticleHandler::addObject does not have a species yet."
135  "Please make sure that you have "
136  "set the species somewhere in the driver code.", P->getId());
137  }
138 
139  //Puts the particle in the Particle list
141  if (getDPMBase() != nullptr)
142  {
143  //This places the particle in this grid
145  //This computes where the particle currently is in the grid
147  }
148  //set the particleHandler pointer
149  P->setHandler(this);
150  //compute mass of the particle
151  P->getSpecies()->computeMass(P);
152  //Check if this particle has new extrema
153  checkExtrema(P);
154 }
155 
156 
171 {
172  if (P->getSpecies() == nullptr)
173  {
174  logger(WARN, "WARNING: The particle with ID % that is added in "
175  "ParticleHandler::addObject does not have a species yet."
176  "Please make sure that you have "
177  "set the species somewhere in the driver code.", P->getId());
178  }
179 #ifdef MERCURY_USE_MPI
180  bool insertParticle;
181  //Check if the particle P should be added to the current domain
182  if (NUMBER_OF_PROCESSORS == 1)
183  {
184  insertParticle = true;
185  }
186  else
187  {
188  insertParticle = getDPMBase()->mpiInsertParticleCheck(P);
189  }
190 
191  //Add the particle if it is in the mpi domain or if the domain is not yet defined
192  if(insertParticle)
193  {
194 #endif
195  //Puts the particle in the Particle list
197  if (getDPMBase() != nullptr)
198  {
199  //This places the particle in this grid
201  //This computes where the particle currently is in the grid
203  }
204  //set the particleHandler pointer
205  P->setHandler(this);
206  //compute mass of the particle
207  P->getSpecies()->computeMass(P);
208  //Check if this particle has new extrema
209  checkExtrema(P);
210  if (!P->isSphericalParticle())
211  {
212  getDPMBase()->setRotation(true);
213  }
214 #ifdef MERCURY_USE_MPI
215  P->setPeriodicComplexity(std::vector<int>(0));
216  }
217  else
218  {
219  logger.assert(!P->isMPIParticle(),"Can't add mpi particle as it does not exist");
220  logger.assert(!P->isPeriodicGhostParticle(),"Can't add mpi particle as it does not exist");
221  //Somehwere a really new particle has been added, so to keep the ID's globally unique, we also update
222  //the unique value on all other processors
224  }
225 
226  //Add the particle to the ghost particle lists
228  //Check if this new particle requires an update in the mpi grid (interactionDistance).
230 
231  //Delete the particle that was supposed to be added (but was not)
232  if(!insertParticle)
233  {
234  delete P;
235  }
236 #endif
237 }
238 
247 void ParticleHandler::addObject(int fromProcessor, BaseParticle* p)
248 {
249 #ifdef MERCURY_USE_MPI
250  MPIContainer& communicator = MPIContainer::Instance();
251 
252  //The processor that contains the particle that needs to be copied needs to identify the target, and communicate this
253  MPIParticle pInfo;
254  if (communicator.getProcessorID() == fromProcessor)
255  {
257  }
258 
259  //Broadcast from processor i
260  communicator.broadcast(&pInfo,MercuryMPIType::PARTICLE,fromProcessor);
261  copyDataFromMPIParticleToParticle(&pInfo, p, this);
262 
263  //All processors now know have the same information and we can proceed with the collective add
264  addObject(p);
265 #else
266  logger(WARN, "Function addObject(int fromProcessor, BaseParticle* p) should not be used in serial code");
267 #endif
268 }
269 
280 void ParticleHandler::addGhostObject(int fromProcessor, int toProcessor, BaseParticle* p)
281 {
282 #ifdef MERCURY_USE_MPI
283  MPIContainer& communicator = MPIContainer::Instance();
284  if (fromProcessor == toProcessor)
285  {
286  if (communicator.getProcessorID() == fromProcessor)
287  {
288  addGhostObject(p);
289  }
290 
291  return;
292  }
293 
294  //Create communication information object
296  MPIParticle pInfo;
297  int tag;
298  if (communicator.getProcessorID() == fromProcessor)
299  {
301  tag = fromProcessor*MAX_PROC*10 + toProcessor*10 + MercuryMPITag::PARTICLE_DATA;
302  communicator.send(&pInfo, MercuryMPIType::PARTICLE, 1, toProcessor, tag);
303  }
304 
305  if (communicator.getProcessorID() == toProcessor)
306  {
307  tag = fromProcessor*MAX_PROC*10 + toProcessor*10 + MercuryMPITag::PARTICLE_DATA;
308  communicator.receive(&pInfo, MercuryMPIType::PARTICLE, 1, fromProcessor, tag);
309  }
310 
311  //Sync the communication
312  communicator.sync();
313 
314  //Add the data to the new particle
315  if (communicator.getProcessorID() == toProcessor)
316  {
317  copyDataFromMPIParticleToParticle(&pInfo, p, this);
318  }
319 
320 
321  //Only toProcessor adds the particle, quietly without any check with other processors
322  //IMPORTANT NOTE: When adding a periodic particle in parallel, this is performed just before
323  //finding new particles. For that reason we dont have to add a ghostparticle to the mpi communication lists
324  if (communicator.getProcessorID() == toProcessor)
325  {
326  addGhostObject(p);
327  }
328 #else
329  logger(WARN,
330  "Function addGhostObject(int fromProcessor, int toProcessor, BaseParticle* p) should not be used in serial code");
331 #endif
332 }
333 
342 {
343 #ifdef MERCURY_USE_MPI
344  if (P->getSpecies() == nullptr)
345  {
346  logger(WARN, "[ParticleHandler::adGhostObject(BaseParticle*)] "
347  "WARNING: The particle with ID % that is added in "
348  "ParticleHandler::addObject does not have a species yet."
349  "Please make sure that you have "
350  "set the species somewhere in the driver code.", P->getId());
351  }
352 
353  MPIContainer& communicator = MPIContainer::Instance();
354  //Puts the particle in the Particle list
356  if (getDPMBase() != nullptr)
357  {
358  //This places the particle in this grid
360  //This computes where the particle currently is in the grid
362  }
363  //set the particleHandler pointer
364  P->setHandler(this);
365  //compute mass of the particle
366  P->getSpecies()->computeMass(P) ;
367  //Check if this particle has new extrema
368  checkExtrema(P);
369 #else
370  logger(INFO,
371  "Function ParticleHandler::mpiAddObject(BaseParticle* P) should only be called when compiling with parallel build on");
372 #endif
373 }
374 
384 void ParticleHandler::removeObject(const unsigned int index)
385 {
386 #ifdef MERCURY_USE_MPI
387  MPIContainer& communicator = MPIContainer::Instance();
388  if (communicator.getNumberOfProcessors() > 1 )
389  {
390  logger(WARN, "[ParticleHandler::removeObject(const unsigned int)] Using the function removeObject in parallel could lead to out-of-sync communication, Instead use deletion boundaries");
391  }
392 #endif
393 #ifdef CONTACT_LIST_HGRID
394  getDPMBase()->getPossibleContactList().remove_ParticlePosibleContacts(getObject(id));
395 #endif
398 }
399 
410 void ParticleHandler::removeGhostObject(const unsigned int index)
411 {
412 #ifdef MERCURY_USE_MPI
413 #ifdef CONTACT_LIST_HGRID
414  getDPMBase()->getPossibleContactList().remove_ParticlePosibleContacts(getObject(id));
415 #endif
418 #else
419  logger(ERROR, "This function should only be used interally for mpi routines");
420 #endif
421 }
422 
428 {
429 #ifdef CONTACT_LIST_HGRID
430  getDPMBase()->getPossibleContactList().remove_ParticlePosibleContacts(getLastObject());
431 #endif
434 }
435 
437 {
438  if (getSize() == 0)
439  {
440  logger(DEBUG, "No particles, so cannot compute the smallest particle.");
441 
442  smallestParticle_ = nullptr;
443  return;
444  }
445  Mdouble min = std::numeric_limits<Mdouble>::max();
446  smallestParticle_ = nullptr;
447  for (BaseParticle* const particle : objects_)
448  {
449  if (!(particle->isMPIParticle() || particle->isPeriodicGhostParticle()))
450  {
451  if (particle->getMaxInteractionRadius() < min)
452  {
453  min = particle->getMaxInteractionRadius();
454  smallestParticle_ = particle;
455  }
456  }
457  }
458 }
459 
461 {
462  if (getSize() == 0)
463  {
464  logger(DEBUG, "No particles, so cannot compute the largest particle.");
465  largestParticle_ = nullptr;
466  return;
467  }
468  Mdouble max = -std::numeric_limits<Mdouble>::max();
469  largestParticle_ = nullptr;
470  for (BaseParticle* const particle : objects_)
471  {
472  if (!(particle->isMPIParticle() || particle->isPeriodicGhostParticle()))
473  {
474  if (particle->getMaxInteractionRadius() > max)
475  {
476  max = particle->getMaxInteractionRadius();
477  largestParticle_ = particle;
478  }
479  }
480  }
481 }
482 
488 {
489  return smallestParticle_;
490 }
491 
497 {
498 #ifdef MERCURY_USE_MPI
499  logger(ERROR,"getSmallestParticle should not be used in parallel; use getSmallestInteractionRadius or ParticleSpecies::getSmallestParticleMass instead");
500  return nullptr;
501 #else
502  return getSmallestParticleLocal();
503 #endif
504 }
505 
511 {
512  return largestParticle_;
513 }
514 
519 {
520 #ifdef MERCURY_USE_MPIO
521  logger(ERROR,"getLargestParticle() should not be used in parallel; use getLargestInteractionRadius instead");
522  return nullptr;
523 #else
524  return getLargestParticleLocal();
525 #endif
526 }
527 
529 {
530  Mdouble ene = 0;
531  for (auto p : *this)
532  {
533  if (!(p->isMPIParticle() || p->isPeriodicGhostParticle()))
534  {
535  ene += p->getKineticEnergy();
536  }
537  }
538  return ene;
539 }
540 
542 {
543 #ifdef MERCURY_USE_MPI
544  Mdouble kineticEnergyLocal = getKineticEnergyLocal();
545  Mdouble kineticEnergyGlobal = 0.0;
546 
547  //sum up over all domains
548  MPIContainer& communicator = MPIContainer::Instance();
549  communicator.allReduce(kineticEnergyLocal, kineticEnergyGlobal, MPI_SUM);
550 
551  return kineticEnergyGlobal;
552 #else
553  return getKineticEnergyLocal();
554 #endif
555 }
556 
558 {
559  Mdouble ene = 0;
560  for (auto p : *this)
561  {
562  if (!(p->isMPIParticle() || p->isPeriodicGhostParticle()))
563  {
564  ene += p->getRotationalEnergy();
565  }
566  }
567  return ene;
568 }
569 
571 {
572 #ifdef MERCURY_USE_MPI
573  Mdouble rotationalEnergyLocal = getRotationalEnergyLocal();
574  Mdouble rotationalEnergyGlobal = 0.0;
575 
576  //sum up over all domains
577  MPIContainer& communicator = MPIContainer::Instance();
578  communicator.allReduce(rotationalEnergyLocal, rotationalEnergyGlobal, MPI_SUM);
579 
580  return rotationalEnergyGlobal;
581 #else
582  return getRotationalEnergyLocal();
583 #endif
584 }
585 
587 {
588  Mdouble m = 0;
589  for (auto p : *this)
590  if (!(p->isFixed() || p->isMPIParticle() || p->isPeriodicGhostParticle()))
591  m += p->getMass();
592  return m;
593 }
594 
596 {
597 #ifdef MERCURY_USE_MPI
598  Mdouble massLocal = getMassLocal();
599  Mdouble massGlobal = 0.0;
600 
601  //Sum up over all domains
602  MPIContainer& communicator = MPIContainer::Instance();
603  communicator.allReduce(massLocal, massGlobal, MPI_SUM);
604 
605  return massGlobal;
606 #else
607  return getMassLocal();
608 #endif
609 }
610 
612 {
613  Vec3D com = {0, 0, 0};
614  for (auto p : *this)
615  if (!(p->isFixed() || p->isMPIParticle() || p->isPeriodicGhostParticle()))
616  com += p->getMass() * p->getPosition();
617  return com;
618 }
619 
621 {
622 #ifdef MERCURY_USE_MPI
623  Vec3D massTimesPositionLocal = getMassTimesPositionLocal();
624  Vec3D massTimesPositionGlobal = {0.0, 0.0, 0.0};
625 
626  // Sum up over all domains
627  MPIContainer& communicator = MPIContainer::Instance();
628  communicator.allReduce(massTimesPositionLocal.X, massTimesPositionGlobal.X, MPI_SUM);
629  communicator.allReduce(massTimesPositionLocal.Y, massTimesPositionGlobal.Y, MPI_SUM);
630  communicator.allReduce(massTimesPositionLocal.Z, massTimesPositionGlobal.Z, MPI_SUM);
631 
632  return massTimesPositionGlobal;
633 #else
634  return getMassTimesPositionLocal();
635 #endif
636 }
637 
639 {
640  Mdouble m = getMass();
641  if (m == 0)
642  {
643  Vec3D nanvec = {constants::NaN, constants::NaN, constants::NaN};
644  return nanvec;
645  }
646  else
647  return getMassTimesPosition() / m;
648 }
649 
651 {
652  Vec3D momentum = {0, 0, 0};
653  for (auto p : *this)
654  if (!(p->isFixed() || p->isMPIParticle() || p->isPeriodicGhostParticle()))
655  momentum += p->getMomentum();
656  return getMPISum(momentum);
657 }
658 
660 {
661  Vec3D momentum = {0, 0, 0};
662  for (auto p : *this)
663  if (!(p->isFixed() || p->isMPIParticle() || p->isPeriodicGhostParticle()))
664  momentum += p->getAngularMomentum();
665  return getMPISum(momentum);
666 }
667 
672 {
673  if (getSize() == 0)
674  {
675  logger(WARN, "No particles to set getFastestParticle()");
676  return nullptr;
677  }
678  BaseParticle* p = nullptr;
679  Mdouble maxSpeed = -std::numeric_limits<Mdouble>::max();
680  for (BaseParticle* const pLoop : objects_)
681  {
682  if (!(pLoop->isMPIParticle() || pLoop->isPeriodicGhostParticle()))
683  {
684  if ((pLoop->getVelocity().getLength()) > maxSpeed)
685  {
686  maxSpeed = pLoop->getVelocity().getLength();
687  p = pLoop;
688  }
689  }
690  }
691  return p;
692 }
693 
695 {
696 #ifdef MERCURY_USE_MPI
697  logger(ERROR,"This function should not be used in parallel");
698 #endif
699  return getFastestParticleLocal();
700 }
701 
702 /*
703  * \returns The smallest interaction radius on a local domain
704  */
706 {
707  if (!(getSmallestParticleLocal() == nullptr))
708  {
710  }
711  else
712  {
713  return 0.0;
714  }
715 }
716 
717 /*
718  * \returns The smallest interaction radius
719  */
721 {
722 #ifdef MERCURY_USE_MPI
723  //Compute the local value
724  Mdouble smallestInteractionRadiusLocal = getSmallestInteractionRadiusLocal();
725  Mdouble smallestInteractionRadiusGlobal = 0.0;
726 
727  //Obtain the global value
728  MPIContainer& communicator = MPIContainer::Instance();
729  communicator.allReduce(smallestInteractionRadiusLocal, smallestInteractionRadiusGlobal, MPI_MIN);
730 
731  return smallestInteractionRadiusGlobal;
732 
733 #else
735 #endif
736 }
737 
738 /*
739  * \returns the largest interaction radius on a local domain
740  */
742 {
743  if (!(getLargestParticleLocal() == nullptr))
744  {
746  }
747  else
748  {
749  return 0.0;
750  }
751 }
752 
753 /*
754  * \returns the largest interaction radius
755  */
757 {
758 #ifdef MERCURY_USE_MPI
759  Mdouble largestInteractionRadiusLocal = getLargestInteractionRadiusLocal();
760  Mdouble largestInteractionRadiusGlobal = 0.0;
761 
762  //Obtain the global value
763  MPIContainer& communicator = MPIContainer::Instance();
764  communicator.allReduce(largestInteractionRadiusLocal, largestInteractionRadiusGlobal, MPI_MAX);
765 
766  return largestInteractionRadiusGlobal;
767 #else
769 #endif
770 }
771 
776 {
777  Mdouble sumRadius = 0;
778  for (BaseParticle* const p : objects_)
779  {
780  if (!(p->isMPIParticle() || p->isPeriodicGhostParticle()))
781  {
782  sumRadius += p->getRadius();
783  }
784  }
785  return sumRadius;
786 }
787 
789 {
790 #ifdef MERCURY_USE_MPI
791  Mdouble sumRadiusLocal = getSumRadiusLocal();
792  unsigned numberOfRealParticlesLocal = getNumberOfRealObjectsLocal();
793 
794  Mdouble sumRadiusGlobal = 0.0;
795  unsigned numberOfRealParticlesGlobal = 0;
796 
797  //Sum up over all domains
798  MPIContainer& communicator = MPIContainer::Instance();
799  communicator.allReduce(sumRadiusLocal, sumRadiusGlobal, MPI_SUM);
800  communicator.allReduce(numberOfRealParticlesLocal, numberOfRealParticlesGlobal, MPI_SUM);
801 
802  return sumRadiusGlobal/numberOfRealParticlesGlobal;
803 #else
804  return getSumRadiusLocal() / getSize();
805 #endif
806 }
807 
814 {
815 #ifdef MERCURY_USE_MPI
816  logger(ERROR,"getLowestPositionComponentParticle() not implemented yet in parallel");
817 #endif
818  if (getSize() == 0)
819  {
820  logger(WARN, "No getLowestPositionComponentParticle(const int i) since there are no particles.");
821  return nullptr;
822  }
823  BaseParticle* p = nullptr;
824  Mdouble min = std::numeric_limits<Mdouble>::max();
825  for (BaseParticle* const pLoop : objects_)
826  {
827  if (!(pLoop->isMPIParticle() || pLoop->isPeriodicGhostParticle()))
828  {
829  if (pLoop->getPosition().getComponent(i) < min)
830  {
831  min = pLoop->getPosition().getComponent(i);
832  p = pLoop;
833  }
834  }
835  }
836  return p;
837 }
838 
840 {
841 #ifdef MERCURY_USE_MPI
842  logger(ERROR,"This function should not be used in parallel");
843 #endif
845 }
846 
853 {
854  if (getSize() == 0)
855  {
856  logger(WARN, "No getHighestPositionComponentParticle(const int i) since there are no particles.");
857  return nullptr;
858  }
859  BaseParticle* p = nullptr;
860  Mdouble max = -std::numeric_limits<Mdouble>::max();
861  for (BaseParticle* const pLoop : objects_)
862  {
863  if (!(pLoop->isMPIParticle() || pLoop->isPeriodicGhostParticle()))
864  {
865  if (pLoop->getPosition().getComponent(i) > max)
866  {
867  max = pLoop->getPosition().getComponent(i);
868  p = pLoop;
869  }
870  }
871  }
872 
873  return p;
874 }
875 
877 {
878 #ifdef MERCURY_USE_MPI
879  logger(ERROR,"This function should not be used in parallel");
880 #endif
882 }
883 
890 {
891  if (getSize() == 0)
892  {
893  logger(WARN, "No getLowestVelocityComponentParticle(const int i) since there are no particles");
894  return nullptr;
895  }
896  BaseParticle* p = nullptr;
897  Mdouble min = std::numeric_limits<Mdouble>::max();
898  for (BaseParticle* const pLoop : objects_)
899  {
900  if (!(pLoop->isMPIParticle() || pLoop->isPeriodicGhostParticle()))
901  {
902  if (pLoop->getVelocity().getComponent(i) < min)
903  {
904  min = pLoop->getVelocity().getComponent(i);
905  p = pLoop;
906  }
907  }
908  }
909  return p;
910 }
911 
913 {
914 #ifdef MERCURY_USE_MPI
915  logger(ERROR,"This function should not be used in parallel");
916 #endif
918 }
919 
926 {
927  if (!getSize())
928  {
929  logger(WARN, "No getHighestVelocityComponentParticle(const int i) since there are no particles");
930  return nullptr;
931  }
932  BaseParticle* p = nullptr;
933  Mdouble max = -std::numeric_limits<Mdouble>::max();
934  for (BaseParticle* const pLoop : objects_)
935  {
936  if (!(pLoop->isMPIParticle() || pLoop->isPeriodicGhostParticle()))
937  {
938  if (pLoop->getVelocity().getComponent(i) > max)
939  {
940  max = pLoop->getVelocity().getComponent(i);
941  p = pLoop;
942  }
943  }
944  }
945  return p;
946 }
947 
949 {
950 #ifdef MERCURY_USE_MPI
951  logger(ERROR,"This function should not be used in parallel");
952 #endif
954 }
955 
962 {
963  smallestParticle_ = nullptr;
964  largestParticle_ = nullptr;
966 }
967 
974 {
975  //Define the attribute function
976  std::function<Mdouble(BaseParticle*)> particleAttribute = [](BaseParticle* p) { return p->getPosition().X; };
977 
978  //Obtain the MAX attribute
979  Mdouble positionXGlobal = getParticleAttribute<Mdouble>(particleAttribute, AttributeType::MAX);
980 
981  return positionXGlobal;
982 }
983 
988 {
989  return NFixedParticles_;
990 }
991 
996 {
998 }
999 
1004 {
1005  if (type == "BaseParticle") {
1006  //for backwards compatibility
1007  return new SphericalParticle;
1008  }
1009  else if (type == "SphericalParticle")
1010  {
1011  return new SphericalParticle;
1012  }
1013  else if (type == "LiquidFilmParticle")
1014  {
1015  return new LiquidFilmParticle;
1016  }
1017  else if (type == "SuperQuadricParticle")
1018  {
1019  return new SuperQuadricParticle;
1020  }
1021  else if (type == "ThermalParticle")
1022  {
1023  return new ThermalParticle;
1024  }
1025  else
1026  {
1027  logger(ERROR, "Particle type % not understood in restart file. Particle will not be read.", type);
1028  return nullptr;
1029  }
1030 }
1031 
1036 {
1037  std::string type;
1038  BaseParticle* particle = nullptr;
1039  if (getStorageCapacity() > 0)
1040  {
1041  is >> type;
1042  logger(DEBUG, "ParticleHandler::readAndCreateObject(is): type %", type);
1043  particle = createObject(type);
1044  particle->setHandler(this);
1045  particle->read(is);
1046  }
1047  else //for insertion boundaries
1048  {
1049  is >> type;
1050  logger(DEBUG, "ParticleHandler::readAndCreateObject(is): type %", type);
1051  particle = createObject(type);
1052  particle->setHandler(this);
1053  particle->read(is);
1054  }
1055  particle->setSpecies(getDPMBase()->speciesHandler.getObject(particle->getIndSpecies()));
1056  return particle;
1057 }
1058 
1063 {
1065  addExistingObject(o);
1066 }
1067 
1069 // * \param[in] type The first value of the position.
1070 // * \param[in] is The input stream from which the information is read.
1071 // * \details The old objects did not have their type in the beginning of the line.
1072 // * Instead, the first string of the file was the position in x-direction.
1073 // * Since we already read the first string of the file, we need to give
1074 // * it to this function and convert it to the position in x-direction.
1075 // * The rest of the stream is then read in the usual way.
1076 // */
1077 //void ParticleHandler::readAndCreateOldObject(std::istream& is, const std::string& type)
1078 //{
1079 // //read in next line
1080 // std::stringstream line;
1081 // helpers::getLineFromStringStream(is, line);
1082 // logger(VERBOSE, line.str());
1083 // //std::cout << line.str() << std::endl;
1084 //
1085 // BaseParticle particle;
1086 //
1087 // //Declare all properties of the particle
1088 // unsigned int indSpecies;
1089 // Mdouble radius, inverseMass, inverseInertia;
1090 // Vec3D position, velocity, euler, angularVelocity;
1091 //
1092 // //Read all values
1093 // position.X = atof(type.c_str());
1094 //
1095 // line >> position.Y >> position.Z >> velocity >> radius >> euler >> angularVelocity >> inverseMass >> inverseInertia >> indSpecies;
1096 //
1097 // //Put the values in the particle
1098 // particle.setSpecies(getDPMBase()->speciesHandler.getObject(indSpecies));
1099 // particle.setPosition(position);
1100 // particle.setVelocity(velocity);
1101 // particle.setRadius(radius);
1102 // Quaternion q; q.setEuler(euler);
1103 // particle.setOrientation(q);
1104 // particle.setAngularVelocity(angularVelocity);
1105 // if (inverseMass == 0.0)
1106 // particle.fixParticle();
1107 // else
1108 // {
1109 // particle.setInverseInertia(MatrixSymmetric3D(1,0,0,1,0,1)*inverseInertia);
1110 // }
1111 //
1112 // //Put the particle in the handler
1113 // copyAndAddObject(particle);
1114 //}
1115 
1116 void ParticleHandler::write(std::ostream& os) const
1117 {
1118 #ifdef MERCURY_USE_MPI
1119  os << "Particles " << getNumberOfRealObjectsLocal() << std::endl;
1120  for (BaseParticle* it : *this)
1121  {
1122  if (!it->isPeriodicGhostParticle() && !it->isMPIParticle())
1123  {
1124  os << (*it) << '\n';
1125  }
1126  }
1127 #else
1128  os << "Particles " << getSize() << '\n';
1129  for (BaseParticle* it : *this)
1130  {
1131  os << (*it) << '\n';
1132  }
1133 #endif
1134  //os.flush();
1135 }
1136 
1142 {
1143  if (P == largestParticle_)
1144  {
1145  //if the properties of the largest particle changes
1147  }
1149  {
1150  largestParticle_ = P;
1151  }
1152 
1153  if (P == smallestParticle_)
1154  {
1155  //if the properties of the smallest particle changes
1157  }
1159  {
1160  smallestParticle_ = P;
1161  }
1162 }
1163 
1168 {
1169  if (P == largestParticle_)
1170  {
1172  }
1173  if (P == smallestParticle_)
1174  {
1176  }
1177 }
1178 
1183 void ParticleHandler::computeAllMasses(unsigned int indSpecies)
1184 {
1185  for (BaseParticle* particle : objects_)
1186  {
1187  if (particle->getIndSpecies() == indSpecies)
1188  {
1189  particle->getSpecies()->computeMass(particle);
1190  }
1191  }
1192 }
1193 
1195 {
1196  for (BaseParticle* particle : objects_)
1197  {
1198  particle->getSpecies()->computeMass(particle);
1199  }
1200 }
1201 
1203 {
1204  NFixedParticles_++;
1205 }
1206 
1208 {
1209  NFixedParticles_--;
1210 }
1211 
1215 std::string ParticleHandler::getName() const
1216 {
1217  return "ParticleHandler";
1218 }
1219 
1221 {
1222  Mdouble volume = 0;
1223  for (auto p : *this)
1224  {
1225  if (!(p->isFixed() || p->isPeriodicGhostParticle() || p->isMPIParticle()))
1226  volume += p->getVolume();
1227  }
1228  return volume;
1229 }
1230 
1232 {
1233 #ifdef MERCURY_USE_MPI
1234  Mdouble volumeLocal = getVolumeLocal();
1235  Mdouble volumeGlobal = 0.0;
1236 
1237  // sum up over all domains
1238  MPIContainer& communicator = MPIContainer::Instance();
1239  communicator.allReduce(volumeLocal, volumeGlobal, MPI_SUM);
1240 
1241  return volumeGlobal;
1242 #else
1243  return getVolumeLocal();
1244 #endif
1245 }
1246 
1252 {
1253 #ifdef MERCURY_USE_MPI
1254  const MPIContainer& communicator = MPIContainer::Instance();
1255  if (communicator.getNumberOfProcessors() > 1)
1256  {
1257  unsigned int numberOfFakeParticles = 0;
1258  numberOfFakeParticles += getDPMBase()->domainHandler.getCurrentDomain()->getNumberOfTrueMPIParticles();
1260  logger.assert(numberOfFakeParticles <= getSize(), "More fake particles than getSize()");
1261  return (getSize() - numberOfFakeParticles);
1262  }
1263  else
1264  {
1265  return getSize();
1266  }
1267 #else
1268  return getSize();
1269 #endif
1270 }
1271 
1273 {
1274 #ifdef MERCURY_USE_MPI
1275  MPIContainer& communicator = MPIContainer::Instance();
1276  unsigned int numberOfRealParticles = getNumberOfRealObjectsLocal();
1277 
1278  // \todo MX: use an allreduce here
1279  //Combine the total number of Particles into one number on processor 0
1280  communicator.reduce(numberOfRealParticles, MPI::SUM);
1281 
1282  //Broadcast new number to all the processorsi
1283  communicator.broadcast(numberOfRealParticles);
1284  return numberOfRealParticles;
1285 #else
1286  return getSize();
1287 #endif
1288 }
1289 
1290 /*
1291  * \returns the size of the ParticleHandler. In parallel code it throws an error to avoid confusion with real and fake particles
1292  */
1294 {
1295 #ifdef MERCURY_USE_MPI
1296  MPIContainer& communicator = MPIContainer::Instance();
1297  if (communicator.getNumberOfProcessors() > 1)
1298  {
1299  logger(WARN,"When compiling with MPI please do not use getNumberOfObjects(). Instead use: getNumberOfRealObjectsLocal(), getNumberOfRealObjects() or getSize()");
1300  }
1301 #endif
1302  return getSize();
1303 }
1304 
1311 {
1312  unsigned int numberOfFixedParticles = 0;
1313  for (BaseParticle* particle : *this)
1314  {
1315  if (particle->isFixed())
1316  {
1317  numberOfFixedParticles++;
1318  }
1319  }
1320  return numberOfFixedParticles;
1321 }
1322 
1327 {
1328 #ifdef MERCURY_USE_MPI
1329  unsigned int numberOfFixedParticlesLocal = getNumberOfFixedObjectsLocal();
1330  unsigned int numberOfFixedParticles = 0;
1331 
1332  MPIContainer& communicator = MPIContainer::Instance();
1333  communicator.allReduce(numberOfFixedParticlesLocal, numberOfFixedParticles, MPI_SUM);
1334  return numberOfFixedParticles;
1335 #else
1337 #endif
1338 }
1339 
1341 {
1342  for (auto i: *this)
1343  {
1344  i->actionsAfterTimeStep();
1345  }
1346 }
1347 
1349  double liquidVolume = 0;
1350  for (auto i : objects_) {
1351  auto j = dynamic_cast<LiquidFilmParticle*>(i);
1352  if (j and !j->isMPIParticle()) liquidVolume += j->getLiquidVolume();
1353  }
1354  return getMPISum(liquidVolume);
1355 };
Mdouble getVolume() const
This class contains all information and functions required for communication between processors...
Definition: MpiContainer.h:125
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...
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:1664
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:1657
static MPIContainer & Instance()
fetch the instance to be used for communication
Definition: MpiContainer.h:130
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 basic particle.
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:1692
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
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:213
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.
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:1354
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:1797
unsigned int getNumberOfTrueMPIParticles()
Obtains the number of particles in the particleHandler that are MPIParticles, but NOT periodic partic...
Definition: Domain.cc:1650
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:1349
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:429
Vec3D getMassTimesPositionLocal() const
virtual void hGridRemoveParticle(BaseParticle *obj UNUSED)
Definition: DPMBase.cc:1671
Mdouble getLiquidVolume() const
void sync()
Process all pending asynchronous communication requests before continuing.
Definition: MpiContainer.h:148
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:167
~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:1771
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:538
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