MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Domain.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 "Domain.h"
27 #include "DPMBase.h"
28 #include "Logger.h"
29 #include "MpiContainer.h"
30 #include "MpiDataClass.h"
32 #include "Walls/BaseWall.h"
33 #include "InteractionHandler.h"
34 #include "DomainHandler.h"
35 #include "Math/Vector.h"
36 #include <limits>
37 #include <utility>
38 #include <vector>
39 #include <set>
41 
49 {
50  constructor();
51 #ifdef DEBUG_CONSTRUCTOR
52  std::cout<<"Domain::Domain() finished"<<std::endl;
53 #endif
54 }
55 
62 Domain::Domain(std::vector<unsigned> globalMeshIndex) : globalMeshIndex_(std::move(globalMeshIndex))
63 {
64  constructor();
65 #ifdef DEBUG_CONSTRUCTOR
66  std::cout<<"Domain::Domain() finished"<<std::endl;
67 #endif
68 }
69 
77  : BaseObject(b)
78 {
79  rank_ = b.rank_;
84  middle_ = b.middle_;
85 
86  //A cube has 3^3=27 neighbours
87  unsigned long numberOfNeighbours = 27;
88 
89  //Create all lists
90  localIndexToGlobalIndexTable_ = std::vector<int>(numberOfNeighbours);
91  localIndexToProcessorList_ = std::vector<int>(numberOfNeighbours);
92  boundaryParticleList_ = std::vector<std::vector<BaseParticle*> >(numberOfNeighbours, std::vector<BaseParticle*>(0));
93  boundaryParticleListNeighbour_ = std::vector<std::vector<BaseParticle*> >(numberOfNeighbours,
94  std::vector<BaseParticle*>(0));
95  newBoundaryParticleList_ = std::vector<std::vector<BaseParticle*> >(numberOfNeighbours,
96  std::vector<BaseParticle*>(0));
97  newInteractionList_ = std::vector<std::vector<BaseInteraction*> >(numberOfNeighbours,
98  std::vector<BaseInteraction*>(0));
99  numberOfParticlesSend_ = std::vector<unsigned>(numberOfNeighbours);
100  numberOfParticlesReceive_ = std::vector<unsigned>(numberOfNeighbours);
101  numNewInteractionsSend_ = std::vector<unsigned>(numberOfNeighbours);
102  numNewInteractionsReceive_ = std::vector<unsigned>(numberOfNeighbours);
103  boundaryParticleDataSend_ = std::vector<std::vector<MPIParticle> >(numberOfNeighbours);
104  boundaryParticleDataReceive_ = std::vector<std::vector<MPIParticle> >(numberOfNeighbours);
105  updatePositionDataSend_ = std::vector<std::vector<MPIParticlePosition> >(numberOfNeighbours);
106  updatePositionDataReceive_ = std::vector<std::vector<MPIParticlePosition> >(numberOfNeighbours);
107  updateVelocityDataSend_ = std::vector<std::vector<MPIParticleVelocity> >(numberOfNeighbours);
108  updateVelocityDataReceive_ = std::vector<std::vector<MPIParticleVelocity> >(numberOfNeighbours);
109  interactionDataSend_ = std::vector<void*>(numberOfNeighbours);
110  interactionDataReceive_ = std::vector<void*>(numberOfNeighbours);
111  activeBoundaryList_ = std::vector<bool>(numberOfNeighbours, true);
112  boundaryList_ = std::vector<int>(0);
113 #ifdef DEBUG_CONSTRUCTOR
114  std::cout<<"Domain::Domain(const Domain &b) finished"<<std::endl;
115 #endif
116 }
117 
122 {
123 #ifdef DEBUG_DESTRUCTOR
124  std::cout << "Domain::~Domain() finished"<<std::endl;
125 #endif
126 }
127 
134 {
136  domainHandler_ = nullptr;
137  domainMin_ = {-constants::inf, -constants::inf, -constants::inf};
138  domainMax_ = {constants::inf, constants::inf, constants::inf};
139 
140  //A cube has 3^3=27 neighbours
141  int numberOfNeighbours = 27;
142 
143  //Create all lists
144  localIndexToGlobalIndexTable_ = std::vector<int>(numberOfNeighbours);
145  localIndexToProcessorList_ = std::vector<int>(numberOfNeighbours);
146  boundaryParticleList_ = std::vector<std::vector<BaseParticle*> >(numberOfNeighbours, std::vector<BaseParticle*>(0));
147  boundaryParticleListNeighbour_ = std::vector<std::vector<BaseParticle*> >(numberOfNeighbours,
148  std::vector<BaseParticle*>(0));
149  newBoundaryParticleList_ = std::vector<std::vector<BaseParticle*> >(numberOfNeighbours,
150  std::vector<BaseParticle*>(0));
151  newInteractionList_ = std::vector<std::vector<BaseInteraction*> >(numberOfNeighbours,
152  std::vector<BaseInteraction*>(0));
153  numberOfParticlesSend_ = std::vector<unsigned>(numberOfNeighbours);
154  numberOfParticlesReceive_ = std::vector<unsigned>(numberOfNeighbours);
155  numNewInteractionsSend_ = std::vector<unsigned>(numberOfNeighbours);
156  numNewInteractionsReceive_ = std::vector<unsigned>(numberOfNeighbours);
157  boundaryParticleDataSend_ = std::vector<std::vector<MPIParticle> >(numberOfNeighbours);
158  boundaryParticleDataReceive_ = std::vector<std::vector<MPIParticle> >(numberOfNeighbours);
159  updatePositionDataSend_ = std::vector<std::vector<MPIParticlePosition> >(numberOfNeighbours);
160  updatePositionDataReceive_ = std::vector<std::vector<MPIParticlePosition> >(numberOfNeighbours);
161  updateVelocityDataSend_ = std::vector<std::vector<MPIParticleVelocity> >(numberOfNeighbours);
162  updateVelocityDataReceive_ = std::vector<std::vector<MPIParticleVelocity> >(numberOfNeighbours);
163  interactionDataSend_ = std::vector<void*>(numberOfNeighbours);
164  interactionDataReceive_ = std::vector<void*>(numberOfNeighbours);
165  activeBoundaryList_ = std::vector<bool>(numberOfNeighbours, true);
166  boundaryList_ = std::vector<int>(0);
167 }
168 
176 {
177  return new Domain(*this);
178 }
179 
186 void Domain::read(std::istream& is)
187 {
188  logger(WARN, "[Domain::read] should not be called");
189  //BaseObject::read(is);
190 }
191 
197 void Domain::write(std::ostream& os) const
198 {
199  logger(WARN, "[Domain::write] should not be called");
200  //BaseObject::write(os);
201 }
202 
207 std::string Domain::getName() const
208 {
209  return "Domain";
210 }
211 
219 void Domain::setRange(Direction direction, Mdouble min, Mdouble max)
220 {
221  if (min > max)
222  {
223  logger(ERROR, "[MercuryMPI ERROR]: min is larger than max. (%,%)", min, max);
224  }
225 
226  double maxClosed;
227  double minClosed;
228  if (min == -constants::inf)
229  {
230  minClosed = getHandler()->getDPMBase()->getMin().getComponent(direction);
231  }
232 
233  if (max == constants::inf)
234  {
235  maxClosed = getHandler()->getDPMBase()->getMax().getComponent(direction);
236  }
237 
238 
239  switch (direction)
240  {
241  case Direction::XAXIS :
242  domainMin_[0] = min;
243  domainMax_[0] = max;
244  middle_.X = minClosed + (maxClosed - minClosed) / 2.0;
245  break;
246  case Direction::YAXIS :
247  domainMin_[1] = min;
248  domainMax_[1] = max;
249  middle_.Y = minClosed + (maxClosed - minClosed) / 2.0;
250  break;
251  case Direction::ZAXIS :
252  domainMin_[2] = min;
253  domainMax_[2] = max;
254  middle_.Z = minClosed + (maxClosed - minClosed) / 2.0;
255  break;
256  default :
257  logger(ERROR, "Direction is not a valid direction. (%)", direction);
258  break;
259  }
260 }
261 
268 void Domain::setBounds(std::vector<double> domainMin, std::vector<double> domainMax, bool computeMiddle)
269 {
270  domainMin_ = domainMin;
271  domainMax_ = domainMax;
272 
273  //Compute the middle of the closed domain
274  if (computeMiddle)
275  {
276  double minClosed;
277  double maxClosed;
278  for (int i = 0; i < 3; i++)
279  {
280  minClosed = domainMin_[i];
281  maxClosed = domainMax_[i];
282  if (domainMin_[i] == -constants::inf)
283  {
284  minClosed = getHandler()->getDPMBase()->getMin().getComponent(i);
285  }
286 
287  if (domainMax_[i] == constants::inf)
288  {
289  maxClosed = getHandler()->getDPMBase()->getMax().getComponent(i);
290  }
291 
292  middle_.setComponent(i, minClosed + (maxClosed - minClosed) / 2.0);
293  }
294  }
295 }
296 
301 void Domain::setRank(int rank)
302 {
303  rank_ = rank;
304 }
305 
310 std::vector<double> Domain::getDomainMin()
311 {
312  return domainMin_;
313 }
314 
319 std::vector<double> Domain::getDomainMax()
320 {
321  return domainMax_;
322 }
323 
329 {
330  return rank_;
331 }
332 
337 void Domain::setHandler(DomainHandler* domainHandler)
338 {
339  domainHandler_ = domainHandler;
340 }
341 
347 {
348  return globalIndex_;
349 }
350 
355 std::vector<unsigned> Domain::getGlobalMeshIndex()
356 {
357  return globalMeshIndex_;
358 }
359 
365 void Domain::setGlobalMeshIndex(std::vector<unsigned> globalMeshIndex)
366 {
367  globalMeshIndex_ = globalMeshIndex;
368 }
369 
378 void Domain::disableBoundary(unsigned localIndex)
379 {
380  activeBoundaryList_[localIndex] = false;
381 }
382 
387 std::vector<bool> Domain::getActiveBoundaryList()
388 {
389  return activeBoundaryList_;
390 }
391 
401 {
402  for (unsigned i = 0; i < 3; i++)
403  {
404  if (!(((domainMin_[i] + offset) < particle->getPosition().getComponent(i))
405  &&
406  ((domainMax_[i] - offset) >= particle->getPosition().getComponent(i)))
407  )
408  {
409  return false;
410  }
411  }
412  return true;
413 }
414 
421 {
423 }
424 
431 {
433 }
434 
442 {
443  //First check if the particle is actually in the domain
444  //Secondly check if the particle is in the inner domain
445  return containsParticle(particle) && !(isInInnerDomain(particle));
446 }
447 
452 {
453  int localIndex;
454  int globalIndex;
455  //Get the global decomposition vector, (nx,ny,nz) and compute the mesh multipliers: (1,nx,nx*ny)
456  std::vector<unsigned> numberOfDomains = domainHandler_->getNumberOfDomains();
457  std::vector<unsigned> globalMeshMultiplier = {1,
458  numberOfDomains[Direction::XAXIS],
459  numberOfDomains[Direction::XAXIS] *
460  numberOfDomains[Direction::YAXIS]};
461 
462  //Compute the globalIndex of this domain
464  globalMeshMultiplier[Direction::YAXIS] * globalMeshIndex_[Direction::YAXIS] +
465  globalMeshMultiplier[Direction::ZAXIS] * globalMeshIndex_[Direction::ZAXIS];
466 
467  //Create the lookup table from localIndex to globalIndex of the neihbours
468  for (int i = -1; i < 2; i++)
469  {
470  for (int j = -1; j < 2; j++)
471  {
472  for (int k = -1; k < 2; k++)
473  {
474  //Get the local index
475  localIndex = i + 3 * j + 9 * k + 13;
476  //Compute the global index
477  globalIndex = globalMeshMultiplier[Direction::XAXIS] * (globalMeshIndex_[Direction::XAXIS] + i) +
478  globalMeshMultiplier[Direction::YAXIS] * (globalMeshIndex_[Direction::YAXIS] + j) +
479  globalMeshMultiplier[Direction::ZAXIS] * (globalMeshIndex_[Direction::ZAXIS] + k);
480  //Fill in the look-up table
481  localIndexToGlobalIndexTable_[localIndex] = globalIndex;
482  }
483  }
484  }
485  //Create the lookup table from localIndex to processorRank
488 }
489 
498 int Domain::getLocalIndex(const int i, const int j, const int k)
499 {
500  return i + 3 * j + 9 * k + 13;
501 }
502 
509 int Domain::getLocalIndex(const std::vector<int> localMeshIndex)
510 {
511  return localMeshIndex[Direction::XAXIS] + 3 * localMeshIndex[Direction::YAXIS] +
512  9 * localMeshIndex[Direction::ZAXIS] + 13;
513 }
514 
515 std::vector<int> Domain::getLocalIndexInverse(int localIndex) {
516  localIndex -= 13;
517  return {localIndex%3, (localIndex/3)%3, (localIndex/9)%3};
518 }
519 
529 BaseParticle* Domain::findParticleInList(unsigned int identification, std::vector<BaseParticle*> particleList)
530 {
531  for (BaseParticle* particle : particleList)
532  {
533  if (particle->getId() == identification)
534  {
535  return particle;
536  }
537  }
538  return nullptr;
539 }
540 
554 std::vector<int> Domain::findNearbyBoundaries(BaseParticle* particle, Mdouble offset)
555 {
556  std::vector<int> boundaryIndex(3);
557  Mdouble interactionDistance = domainHandler_->getInteractionDistance();
558 
559  //for x,y,z directions, find if the particle is close to the left or right wall or not at all
560  for (int d = 0; d < 3; d++)
561  {
562  //Check if the particle is close to Lx
563  if ((particle->getPosition().getComponent(d)) < (domainMin_[d] + interactionDistance + offset))
564  {
565  //Update switch
566  boundaryIndex[d] = -1;
567  }
568  //Check if particle is close to Rx
569  else if ((particle->getPosition().getComponent(d)) >= (domainMax_[d] - interactionDistance - offset))
570  {
571  //Update switch
572  boundaryIndex[d] = 1;
573  }
574  }
575  return boundaryIndex;
576 }
577 
578 bool Domain::inBoundary(BaseParticle* particle, int localIndex_) {
579  std::vector<int> boundaryIndex = findNearbyBoundaries(particle);
580  std::vector<int> localIndex = getLocalIndexInverse(localIndex_);
581  // if the particle is in the boundary zone into which we want to transmit
582  for (int i = 0; i < 3; ++i)
583  {
584  if (localIndex[i]!=0 && boundaryIndex[i]!=localIndex[i]) {
585  return false;
586  }
587  }
588  return true;
589 }
590 
598 {
599  std::vector<unsigned> numberOfDomains = domainHandler_->getNumberOfDomains();
600  std::vector<int> localMeshIndex(3);
601  int localIndex;
602 
603  //disble own list
604  activeBoundaryList_[13] = false;
605 
606  //Special case when numberOfDomains[d] = 1 -> disable all d direction boundaries
607  for (unsigned d = 0; d < 3; d++) //Loop over all dimensions
608  {
609  //If there is only one domain in the d-direction, take action
610  if (numberOfDomains[d] == 1)
611  {
612  //Loop over all locaIndices
613  for (int i = -1; i < 2; i++)
614  {
615  localMeshIndex[Direction::XAXIS] = i;
616  for (int j = -1; j < 2; j++)
617  {
618  localMeshIndex[Direction::YAXIS] = j;
619  for (int k = -1; k < 2; k++)
620  {
621  localMeshIndex[Direction::ZAXIS] = k;
622  //Disable all boundaries that are not in the "middle" of the mesh in the d-direction
623  // i.e. a localMeshIndex[d] = -1 is a neighbour in the d-direction, however there is only one
624  // element in that direction so it can't be a real neighbour.
625  if (localMeshIndex[d] != 0)
626  {
627  //Disable the boundary
628  localIndex = getLocalIndex(localMeshIndex);
629  activeBoundaryList_[localIndex] = false;
630  }
631  }
632  }
633  }
634  }
635  }
636 
637  //Compute the boundaryIndex for all other cases where numberOfDomains is larger than 1
638  std::vector<int> boundaryIndex(3);
639  for (int d = 0; d < 3; d++)
640  {
641  //Check if the domain is on the simulationDomainMin boundary
642  if (globalMeshIndex_[d] == 0)
643  {
644  //Update boundary index
645  boundaryIndex[d] = -1;
646  }
647  //Check if the domain is on the simulationDomainMax boundary
648  else if (globalMeshIndex_[d] == (numberOfDomains[d] - 1))
649  {
650  //Update boundary index
651  boundaryIndex[d] = 1;
652  }
653  }
654 
655  //Disable the boundaries that are on the edge of the simulation domain
656  for (int d = 0; d < 3; d++)
657  {
658  //If boundary is on the edge of the simulation domain
659  if (boundaryIndex[d] != 0)
660  {
661  //Loop over all local indices
662  for (int i = -1; i < 2; i++)
663  {
664  for (int j = -1; j < 2; j++)
665  {
666  for (int k = -1; k < 2; k++)
667  {
668  localMeshIndex = {i, j, k};
669  //Set the localMeshIndex to the boundaryIndex in the d-direction
670  //Since we are only interested in the domains on this boundary
671  localMeshIndex[d] = boundaryIndex[d];
672  //Disable the boundary
673  localIndex = getLocalIndex(localMeshIndex);
674  activeBoundaryList_[localIndex] = false;
675  }
676  }
677  }
678  }
679  }
680 
681  //Create a list of the active boundaries
682  localIndex = 0;
683  for (bool active : activeBoundaryList_)
684  {
685  if (active)
686  {
687  boundaryList_.push_back(localIndex);
688  }
689  localIndex++;
690  }
691 }
692 
701 void Domain::addParticlesToLists(BaseParticle* particle, std::vector<std::vector<BaseParticle*> >& list)
702 {
703  std::vector<int> boundaryIndex = findNearbyBoundaries(particle);
704 
705  //Compute and set complexity of the particle
706  unsigned int complexity = boundaryIndex[0] + 3 * boundaryIndex[1] + 9 * boundaryIndex[2] + 13;
707  unsigned int list_complexity = 0;
708  for (int d = 0; d < 3; d++) //Loop over all directions
709  {
710  list_complexity += std::abs(boundaryIndex[d]);
711  }
712  particle->setCommunicationComplexity(complexity);
713  //particle->setCommunicationComplexity(list_complexity);
714 
715  //Based on the complexity of the particle, add it to the approriate list
716  switch (list_complexity)
717  {
718  //The particle is not close at all
719  case 0:
720  break;
721  //The particle is close to one side
722  // 1 side contribution
723  case 1 :
724  //Add the side contribution
725  list[getLocalIndex(boundaryIndex)].push_back(particle);
726  break;
727  //The particle is close to two neighbouring directions
728  //2 side and 1 rib contrubution
729  case 2 :
730  {
731  //Add the two side contributions
732  for (int d = 0; d < 3; d++)
733  {
734  std::vector<int> localMeshIndex = {0, 0, 0};
735  localMeshIndex[d] = boundaryIndex[d];
736  //Avoid adding the particle in the wrong direction by excluding localMeshIndex[d] = 0
737  if (localMeshIndex[d] != 0)
738  {
739  list[getLocalIndex(localMeshIndex)].push_back(particle);
740  }
741  }
742  }
743  //Add the rib contribution
744  list[getLocalIndex(boundaryIndex)].push_back(particle);
745  break;
746 
747  //The particle is close to three neighbouring directions
748  // 3 side, 3 rib and 1 corner contribution
749  case 3 :
750  {
751  //Add the three side contributions
752  for (int d = 0; d < 3; d++)
753  {
754  std::vector<int> localMeshIndex = {0, 0, 0};
755  //Reset index vector
756  localMeshIndex[d] = boundaryIndex[d];
757  list[getLocalIndex(localMeshIndex)].push_back(particle);
758  }
759 
760  //Add the three rib contributions
761  for (int d = 0; d < 3; d++)
762  {
763  std::vector<int> localMeshIndex = boundaryIndex;
764  //All rib boundary indices are given by the boundaryIndex and setting one of the components to zero
765  localMeshIndex[d] = 0;
766  list[getLocalIndex(localMeshIndex)].push_back(particle);
767  }
768  }
769 
770  //Add the corner contribution
771  list[getLocalIndex(boundaryIndex)].push_back(particle);
772  break;
773 
774  default :
775  logger(INFO, "boundaryIndex : %,%,% | list_complexity: %", boundaryIndex[0], boundaryIndex[1], boundaryIndex[2],
776  list_complexity);
777  logger(ERROR, "Particle is in contact with the wrong number of boundaries");
778  break;
779  }
780 }
781 
786 void Domain::findNewMPIParticles(const ParticleHandler& particleHandler)
787 {
788  //For all particles inside the given domain, loop over all the particles to
789  //see if we have to add the particles to the new particle list
790  for (BaseParticle* particle : particleHandler)
791  {
792  findNewMPIParticle(particle);
793  }
794 }
795 
801 {
802  //If the particle is an MPI particle
803  if (!particle->isInMPIDomain() && !particle->isPeriodicGhostParticle())
804  {
806  }
807 }
808 
809 bool Domain::isInNewBoundaryParticleList(BaseParticle* objectI,int localIndex) const
810 {
811  for (BaseParticle *q : newBoundaryParticleList_[localIndex]) {
812  if (q == objectI)
813  return true;
814  }
815  return false;
816 }
817 
825  {
826  //For all active boundaries
827  for (int localIndex : boundaryList_)
828  {
829  //Check all newly added particles for interactions with already known particles
830  for (BaseParticle* newBoundaryParticle : newBoundaryParticleList_[localIndex])
831  {
832  //Loop over all interactions of the new particle
833  std::vector<BaseInteraction*> interactions = newBoundaryParticle->getInteractions();
834  for (BaseInteraction* interaction : interactions)
835  {
836  //Find out what the new particle is interacting with
837  BaseParticle* particleP = dynamic_cast<BaseParticle*>(interaction->getP());
838  BaseParticle* objectI = dynamic_cast<BaseParticle*>(interaction->getI());
839 
840  //If the P in the interaction structure is the new particle, find I
841  if (newBoundaryParticle == particleP)
842  {
843  //Check if the new particle is interacting with a wall
844  if (!objectI)
845  {
846  //Check if the interaction is still valid
847  BaseWall* wall = dynamic_cast<BaseWall*>(interaction->getI());
848  Mdouble distance;
849  Vec3D normal;
850  wall->getDistanceAndNormal(*particleP, distance, normal);
851  if (distance <= particleP->getMaxInteractionRadius())
852  {
853  //Add the interaction to the list if there are still in contact (hence the if statement)
854  newInteractionList_[localIndex].push_back(interaction);
855  }
856  }
857  else //is I a particle
858  {
859  //check if particle is in mpi domain OR (TODO) if it is in the newBoundaryParticleList_
860  if (objectI->isInMPIDomain() || isInNewBoundaryParticleList(objectI, localIndex))
861  {
862  // iff the interaction partner of the new particle is in the communication zone, transmit the interaction
863  if (inBoundary(objectI,localIndex))
864  {
865  //Is the particle still interacting after the position update?
866  Vec3D branchVector = particleP->getPosition() - objectI->getPosition();
867  //Get the square of the distance between particle i and particle j
868  Mdouble distanceSquared = Vec3D::getLengthSquared(branchVector);
869  Mdouble sumOfInteractionRadii =
870  objectI->getSumOfInteractionRadii(particleP);
871  if (distanceSquared < (sumOfInteractionRadii * sumOfInteractionRadii))
872  {
873  //Add the interaction to the list
874  newInteractionList_[localIndex].push_back(interaction);
875  }
876  }
877  }
878  }
879  }
880  else //newBoundaryParticle is I in the interaction, P can only be a particle
881  {
882  //it is "I" in the interaction structure, check if its in the mpi domain
883  if (particleP->isInMPIDomain() || isInNewBoundaryParticleList(particleP, localIndex))
884  {
885  //Is the particle still interacting after the position update?
886  Vec3D branchVector = particleP->getPosition() - objectI->getPosition();
887  //Get the square of the distance between particle i and particle j
888  Mdouble distanceSquared = Vec3D::getLengthSquared(branchVector);
889  Mdouble sumOfInteractionRadii =
890  objectI->getSumOfInteractionRadii(particleP);
891  if (distanceSquared < (sumOfInteractionRadii * sumOfInteractionRadii))
892  {
893  //Add the interaction to the list
894  newInteractionList_[localIndex].push_back(interaction);
895  }
896  }
897  }
898  }
899  }
900  }
901 }
902 
910 {
911  //For all particles
912  for (BaseParticle* particle : newBoundaryParticleList_[localIndex])
913  {
914  //Add the data to the transmission list
915  boundaryParticleDataSend_[localIndex].push_back(copyDataFromParticleToMPIParticle(particle));
916  }
917 }
918 
924 void Domain::collectInteractionData(int localIndex)
925 {
926  //Copy interactions to the data array
927  unsigned int indexInteraction = 0;
928  for (BaseInteraction* interaction : newInteractionList_[localIndex])
929  {
930  interaction->getMPIInteraction(interactionDataSend_[localIndex], indexInteraction);
931  indexInteraction++;
932  }
933 }
934 
941 void Domain::processReceivedBoundaryParticleData(const unsigned index, std::vector<BaseParticle*>& newParticles)
942 {
943 
944  ParticleHandler& particleHandler = getHandler()->getDPMBase()->particleHandler;
945  //for all MPIPparticles that have been received on the other side
946  for (int i = 0; i < numberOfParticlesReceive_[index]; i++)
947  {
948  //copy data from data
950  copyDataFromMPIParticleToParticle(&(boundaryParticleDataReceive_[index][i]), p0, &particleHandler);
951 
952  //Flag that the particle is a ghost particle of a real particle in the neighbour domain
953  p0->setMPIParticle(true);
954  //Flag that the particle is list as a particle in the mpi domain
955  p0->setInMPIDomain(true);
956  //add particle to handler
957  particleHandler.addGhostObject(p0);
958  //Set previous position as current position
959  BaseParticle* pGhost = particleHandler.getLastObject();
960  pGhost->setPreviousPosition(particleHandler.getLastObject()->getPosition());
961  //Add to the newParticle list
962  newParticles.push_back(pGhost);
963  logger(VERBOSE, "Adding mpi particle % at: %", particleHandler.getLastObject()->getId(),
964  particleHandler.getLastObject()->getPosition());
965 
966  //add pointer to the particle to the list and set the status to idle
967  boundaryParticleListNeighbour_[index].push_back(particleHandler.getLastObject());
968  }
969 }
970 
977 void Domain::processSentBoundaryParticles(const unsigned index)
978 {
979  for (BaseParticle* particle : newBoundaryParticleList_[index])
980  {
981  boundaryParticleList_[index].push_back(particle);
982  particle->setInMPIDomain(true);
983  }
984 }
985 
999 void Domain::processReceivedInteractionData(const unsigned localIndex, std::vector<BaseParticle*>& newParticles)
1000 {
1002  for (unsigned int l = 0; l < numNewInteractionsReceive_[localIndex]; l++)
1003  {
1004  unsigned int identificationP;
1005  unsigned int identificationI;
1006  bool isWallInteraction;
1007  unsigned timeStamp;
1008 
1009  //Get the general information required to setup a new interaction
1011  l, identificationP, identificationI,
1012  isWallInteraction, timeStamp);
1013 
1014  logger(VERBOSE, "interaction details: %, %, idP %, idI %, wall %, time %", interactionDataReceive_[localIndex],
1015  l, identificationP, identificationI,
1016  isWallInteraction, timeStamp);
1017 
1018  //Find the particle in the newParticle list
1019  BaseParticle* pGhost = nullptr;
1020  int idOther;
1021  for (BaseParticle* particle : newParticles)
1022  {
1023  if (particle->getId() == identificationP)
1024  {
1025  pGhost = particle;
1026  idOther = identificationI;
1027  break;
1028  }
1029 
1030  if (particle->getId() == identificationI)
1031  {
1032  pGhost = particle;
1033  idOther = identificationP;
1034  break;
1035  }
1036  }
1037  if (pGhost == nullptr)
1038  {
1039  logger(WARN, "In Domain::processReceivedInteractionData: pGhost (id %) is nullptr, the interaction data is not copied. Two particles possibly moved into domain simultaneously.", identificationP);
1040  continue;
1041  }
1042 
1043  //If it is a wall interaction, do stuff
1044  if (isWallInteraction)
1045  {
1046  BaseInteractable* I = getHandler()->getDPMBase()->wallHandler.getObjectById(identificationI);
1047  //Create interactions
1048  BaseInteraction* j = I->getInteractionWith(pGhost, timeStamp, &iH);
1049  if (j!= nullptr) j->setMPIInteraction(interactionDataReceive_[localIndex], l, false);
1050 
1051  }
1052  else
1053  {
1054  //Obtain potential interaction particles
1055  std::vector<BaseParticle*> interactingParticleList;
1056  getHandler()->getDPMBase()->hGridGetInteractingParticleList(pGhost, interactingParticleList);
1057 
1058  //Find the other interacting particle
1059  BaseParticle* otherParticle = nullptr;
1060  for (BaseParticle* p2 : interactingParticleList)
1061  {
1062  if (p2->getId() == idOther)
1063  {
1064  otherParticle = p2;
1065  break;
1066  }
1067  }
1068  if (otherParticle == nullptr) {
1069  //search for the other particle in the newParticles list
1070  for (BaseParticle *particle : newParticles) {
1071  if (particle->getId() == idOther) {
1072  otherParticle = particle;
1073  break;
1074  }
1075  }
1076  if (otherParticle == nullptr) {
1077  logger(WARN,
1078  "In Domain::processReceivedInteractionData: otherParticle (id %) is nullptr, the interaction data with pGhost (id %) is not copied. Two particles possibly moved into domain simultaneously. nt = %",
1079  identificationI, identificationP, getHandler()->getDPMBase()->getNumberOfTimeSteps());
1080  continue;
1081  }
1082  }
1083  //Add the interaction
1084  // flipped order of other and ghost, so it is in the same order as received
1085  BaseInteraction* j;
1086  if (otherParticle->getId() < pGhost->getId()) {
1087  j = pGhost->getInteractionWith(otherParticle, timeStamp, &iH);
1088  } else {
1089  j = otherParticle->getInteractionWith(pGhost, timeStamp, &iH);
1090  }
1091  if (j != nullptr) j->setMPIInteraction(interactionDataReceive_[localIndex], l, false);
1092  }
1093  }
1094 }
1095 
1097  std::stringstream s;
1098  std::stringstream m;
1099  for (auto p : getHandler()->getDPMBase()->particleHandler) {
1100  if (p->isMPIParticle())
1101  m << std::setw(4) << p->getId();
1102  else
1103  s << std::setw(4) << p->getId();
1104  }
1105  logger(INFO,"Particles %, MPI %",s.str(),m.str());
1106 }
1107 
1108 
1119  void Domain::sendAndReceiveCount(MercuryMPITag tag, unsigned& countReceive, unsigned& countSend,
1120  unsigned localIndexNeighbour)
1121  {
1122  int globalIndexNeighbour = localIndexToGlobalIndexTable_[localIndexNeighbour];
1123  int processor = localIndexToProcessorList_[localIndexNeighbour];
1124 
1125  //Create communication tags
1126  int tagReceive = globalIndexNeighbour * MAX_PROC + globalIndex_ * 10 + tag;
1127  int tagSend = globalIndex_ * MAX_PROC + globalIndexNeighbour * 10 + tag;
1128 
1129  logger.assert_debug(tagSend > 0, "Send tag is wrong. Tag: %", tagSend);
1130  logger.assert_debug(tagReceive > 0, "Receive tag is wrong. Tag: %", tagReceive);
1131  logger.assert_debug(processor >= 0, "Processor is wrong. processor: %", processor);
1132 
1133  //Communicate the requests
1134  MPIContainer::Instance().receive(countReceive, processor, tagReceive);
1135  MPIContainer::Instance().send(countSend, processor, tagSend);
1136  }
1137 
1147  {
1148  //Find new particles that have entered the communication zone
1149  findNewMPIParticles(getHandler()->getDPMBase()->particleHandler);
1150 
1151  //Find interactions between new particles and other particles in the communication zone
1153 
1154  //Compute number of particles to send
1155  //For all boundaries
1156  for (int localIndex : boundaryList_)
1157  {
1158  numberOfParticlesSend_[localIndex] = newBoundaryParticleList_[localIndex].size();
1159  numNewInteractionsSend_[localIndex] = newInteractionList_[localIndex].size();
1160  }
1161 
1162  //For all active boundaries
1163  for (int localIndex : boundaryList_)
1164  {
1165  //Send and receive the number of new boundary particles
1167  (numberOfParticlesSend_[localIndex]), localIndex);
1168 
1169  //Send and receive the new interactions
1171  (numNewInteractionsSend_[localIndex]), localIndex);
1172  }
1173  }
1174 
1184  {
1185  //Assign the particle to the correct lists if the particle belongs to this domain
1186  if (containsParticle(particle))
1187  {
1188  findNewMPIParticle(particle);
1189  }
1190 
1191  //Note: When inserting a single particle, it has no interactions, so no interactions to be copied either
1192 
1193  //Compute number of particles to send
1194  //For all boundaries
1195  for (int localIndex : boundaryList_)
1196  {
1197  numberOfParticlesSend_[localIndex] = newBoundaryParticleList_[localIndex].size();
1198  numNewInteractionsSend_[localIndex] = newInteractionList_[localIndex].size();
1199  }
1200 
1201 
1202  //For all active boundaries
1203  for (int localIndex : boundaryList_)
1204  {
1205  //Send and recieve the number of new boundary particles
1207  (numberOfParticlesSend_[localIndex]), localIndex);
1208 
1209  //Send and receive the new interactions
1211  (numNewInteractionsSend_[localIndex]), localIndex);
1212  }
1213  }
1214 
1222  {
1223  //For all active boundaries
1224  for (int localIndex : boundaryList_)
1225  {
1226  //make sure enough memory is reserved to receive the data
1227  boundaryParticleDataReceive_[localIndex].resize(numberOfParticlesReceive_[localIndex]);
1228 
1229  //Collect the particle data
1230  collectBoundaryParticleData(localIndex);
1231 
1232  //Send the data
1234  boundaryParticleDataReceive_[localIndex].data(), numberOfParticlesReceive_[localIndex],
1235  boundaryParticleDataSend_[localIndex].data(), numberOfParticlesSend_[localIndex],
1236  localIndex);
1237 
1238 
1239  //create arrays for sending and receiving interaction data
1241  numNewInteractionsSend_[localIndex]);
1243  numNewInteractionsReceive_[localIndex]);
1244 
1245  //Collect the interaction data
1246  collectInteractionData(localIndex);
1247 
1248  //Send the data
1250  interactionDataReceive_[localIndex], numNewInteractionsReceive_[localIndex],
1251  interactionDataSend_[localIndex], numNewInteractionsSend_[localIndex], localIndex);
1252  }
1253  }
1254 
1263  {
1264  //For all boundaries
1265  for (int localIndex : boundaryList_)
1266  {
1267  //copy the received data into the particleHandler and neighbour particle list
1268  std::vector<BaseParticle*> newParticles;
1269  processReceivedBoundaryParticleData(localIndex, newParticles);
1270  //All particles in the current domain that have been send to the other domains need to be flagged as communicating particles
1271  processSentBoundaryParticles(localIndex);
1272  //copy the received interaction data into the standard dpm structure
1273  processReceivedInteractionData(localIndex, newParticles);
1274 
1275  //Delete all send/receive data
1276  boundaryParticleDataSend_[localIndex].clear();
1277  boundaryParticleDataReceive_[localIndex].clear();
1280  interactionDataReceive_[localIndex]);
1281 
1282  //Reset all counters
1283  numberOfParticlesSend_[localIndex] = 0;
1284  numberOfParticlesReceive_[localIndex] = 0;
1285  numNewInteractionsSend_[localIndex] = 0;
1286  numNewInteractionsReceive_[localIndex] = 0;
1287 
1288  //Reset all lists
1289  newBoundaryParticleList_[localIndex].clear();
1290  newInteractionList_[localIndex].clear();
1291  }
1292  }
1293 
1308  void Domain::updateParticles(std::set<BaseParticle*>& ghostParticlesToBeDeleted)
1309  {
1310  int complexityNew;
1311  std::vector<int> boundaryIndex;
1312 
1313  //For all active boundaries
1314  for (int localIndex : boundaryList_)
1315  {
1316  //Step 1A: Remove the particles from the boundaryParticleList_ which require re-assignment
1317  //or have just moved away from the region
1318  for (int p = 0; p < boundaryParticleList_[localIndex].size(); p++)
1319  {
1320  BaseParticle* particle = boundaryParticleList_[localIndex][p];
1321  //Check if the particle is still in the domain, but not in the communication zone
1322  if (containsParticle(particle))
1323  {
1324  //check if the complexity has changed
1325  boundaryIndex = findNearbyBoundaries(particle);
1326  complexityNew = boundaryIndex[0] + 3 * boundaryIndex[1] + 9 * boundaryIndex[2] + 13;
1327  if (particle->getCommunicationComplexity() != complexityNew)
1328  {
1329  logger(VERBOSE, "time: % | global index: % in list % | particle % | CURRENT DOMAIN - CHANGES "
1330  "COMPLEXITY", domainHandler_->getDPMBase()->getTime(), globalIndex_,
1331  localIndexToGlobalIndexTable_[localIndex], particle->getId());
1332  //Flag the particle that it no longer participates in the communication layer
1333  //so it can be reintroduced in the transmission step
1334  particle->setMPIParticle(false);
1335  particle->setInMPIDomain(false);
1336  particle->setCommunicationComplexity(0);
1337  boundaryParticleList_[localIndex][p] = nullptr;
1338  }
1339  }
1340  else
1341  {
1342  logger(VERBOSE, "time: % | global index: % in list % | particle % | CURRENT DOMAIN - TO NEIGHBOURING "
1343  "DOMAIN", domainHandler_->getDPMBase()->getTime(), globalIndex_,
1344  localIndexToGlobalIndexTable_[localIndex], particle->getId());
1345 
1346  ghostParticlesToBeDeleted.insert(particle);
1347  boundaryParticleList_[localIndex][p] = nullptr;
1348  }
1349  }
1350 
1351  //Step 1B: Remove the particles from the boundaryParticleListNeightbour_ which require re-assignment
1352  //or have just moved away from the region
1353  for (int p = 0; p < boundaryParticleListNeighbour_[localIndex].size(); p++)
1354  {
1355  BaseParticle* particle = boundaryParticleListNeighbour_[localIndex][p];
1356  //Check if the particle moved out of the neighbour domain
1358  {
1359  //The particle has moved to this domain
1360  if (containsParticle(particle))
1361  {
1362  logger(VERBOSE,
1363  "time: % | global index: % in list % | particle % | NEIGHBOURING DOMAIN - TO CURRENT DOMAIN",
1365  localIndexToGlobalIndexTable_[localIndex], particle->getId());
1366 
1367  //Flag the particle as not yet communicated
1368  particle->setMPIParticle(false);
1369  particle->setInMPIDomain(false);
1370  particle->setCommunicationComplexity(0);
1371  boundaryParticleListNeighbour_[localIndex][p] = nullptr;
1372  }
1373  //The particle has moved to a different domain
1374  else
1375  {
1376  logger(VERBOSE,
1377  "time: % | global index: % in list % | particle % | NEIGBOURING DOMAIN - TO OTHER DOMAIN",
1379  localIndexToGlobalIndexTable_[localIndex], particle->getId());
1380  //Cruelly destroy the particle without any mercy.
1381  ghostParticlesToBeDeleted.insert(particle);
1382  boundaryParticleListNeighbour_[localIndex][p] = nullptr;
1383  }
1384  }
1385  else
1386  {
1387  //check if the complexity has changed
1388  boundaryIndex = domainHandler_->getObject(
1389  localIndexToGlobalIndexTable_[localIndex])->findNearbyBoundaries(particle);
1390  complexityNew = boundaryIndex[0] + 3 * boundaryIndex[1] + 9 * boundaryIndex[2] + 13;
1391  if (particle->getCommunicationComplexity() != complexityNew)
1392  {
1393  logger(VERBOSE,
1394  "time: % | global index: % in list % | particle % | NEIGHBOURING DOMAIN - CHANGES COMPLEXITY",
1396  localIndexToGlobalIndexTable_[localIndex], particle->getId());
1397  //Cruelly destroy the particle without any mercy.
1398  ghostParticlesToBeDeleted.insert(particle);
1399  boundaryParticleListNeighbour_[localIndex][p] = nullptr;
1400  }
1401  }
1402  }
1403  }
1404  }
1405 
1412  void Domain::updateParticlePosition(int localIndex)
1413  {
1414  //process the updated information
1415  unsigned int index = 0;
1416  for (BaseParticle* particle : boundaryParticleListNeighbour_[localIndex])
1417  {
1418  logger.assert_debug(particle->getId() == updatePositionDataReceive_[localIndex][index].id,
1419  "MPI particle lists are not in sync");
1420 
1421  //set position
1422  particle->setPreviousPosition(particle->getPosition());
1423  particle->setPosition(updatePositionDataReceive_[localIndex][index].position);
1424  particle->setOrientation(updatePositionDataReceive_[localIndex][index].orientation);
1425  if (std::is_base_of<MPILiquidFilmParticle,MPIParticle>())
1426  static_cast<LiquidFilmParticle*>(particle)->setLiquidVolume(updatePositionDataReceive_[localIndex][index].liquidVolume);
1427 
1428  //Update hGrid
1429  Vec3D displacement = particle->getPreviousPosition() - particle->getPosition();
1430  getHandler()->getDPMBase()->hGridUpdateMove(particle, displacement.getLengthSquared());
1431 
1432  index++;
1433  }
1434 }
1435 
1442 {
1443  //process the updated information
1444  unsigned int index = 0;
1445  for (BaseParticle* particle: boundaryParticleListNeighbour_[localIndex])
1446  {
1447  //set velocity
1448  particle->setVelocity(updateVelocityDataReceive_[localIndex][index].velocity);
1449  particle->setAngularVelocity(updateVelocityDataReceive_[localIndex][index].angularVelocity);
1450  index++;
1451  }
1452 }
1453 
1461 {
1462  //For all active boundaries
1463  for (int localIndex : boundaryList_)
1464  {
1465  numberOfParticlesSend_[localIndex] = boundaryParticleList_[localIndex].size();
1466  numberOfParticlesReceive_[localIndex] = boundaryParticleListNeighbour_[localIndex].size();
1467 
1468  //Increase capacity for the receiving data files
1469  updatePositionDataReceive_[localIndex].resize(numberOfParticlesReceive_[localIndex]);
1470  updateVelocityDataReceive_[localIndex].resize(numberOfParticlesReceive_[localIndex]);
1471 
1472 
1473  //Collect data
1474  for (BaseParticle* particle : boundaryParticleList_[localIndex])
1475  {
1476  updatePositionDataSend_[localIndex].push_back(copyPositionFrom(particle));
1477  updateVelocityDataSend_[localIndex].push_back(copyVelocityFrom(particle));
1478  }
1479 
1480  //Send and receive the data
1482  updatePositionDataReceive_[localIndex].data(), numberOfParticlesReceive_[localIndex],
1483  updatePositionDataSend_[localIndex].data(), numberOfParticlesSend_[localIndex],
1484  localIndex);
1486  updateVelocityDataReceive_[localIndex].data(), numberOfParticlesReceive_[localIndex],
1487  updateVelocityDataSend_[localIndex].data(), numberOfParticlesSend_[localIndex],
1488  localIndex);
1489  }
1490 }
1491 
1501 void Domain::finalisePositionAndVelocityUpdate(std::set<BaseParticle*>& ghostParticlesToBeDeleted)
1502 {
1503 
1504  //For all active boundaries
1505  for (int localIndex : boundaryList_)
1506  {
1507  updateParticlePosition(localIndex);
1508  updateParticleVelocity(localIndex);
1509  }
1510 
1511  //Based on the new position, update the particle lists.
1512  //Remove particles that left the communication zone, they will be re-communicated in a later step
1513  updateParticles(ghostParticlesToBeDeleted);
1514 
1515  //For all active boundaries clear the data lists
1516  for (int localIndex : boundaryList_)
1517  {
1518  updatePositionDataSend_[localIndex].clear();
1519  updateVelocityDataSend_[localIndex].clear();
1520  updatePositionDataReceive_[localIndex].clear();
1521  updateVelocityDataReceive_[localIndex].clear();
1522  }
1523 
1524 }
1525 
1530 {
1531  //For all active boundaries
1532  for (int localIndex : boundaryList_)
1533  {
1534  numberOfParticlesSend_[localIndex] = boundaryParticleList_[localIndex].size();
1535  numberOfParticlesReceive_[localIndex] = boundaryParticleListNeighbour_[localIndex].size();
1536 
1537  //Resize the vector to the correct size
1538  updateVelocityDataReceive_[localIndex].resize(numberOfParticlesReceive_[localIndex]);
1539 
1540  //Collect data
1541  for (BaseParticle* particle : boundaryParticleList_[localIndex])
1542  {
1543  updateVelocityDataSend_[localIndex].push_back(copyVelocityFrom(particle));
1544  }
1545 
1546  //Send the data
1548  updateVelocityDataReceive_[localIndex].data(), numberOfParticlesReceive_[localIndex],
1549  updateVelocityDataSend_[localIndex].data(), numberOfParticlesSend_[localIndex],
1550  localIndex);
1551  }
1552 }
1553 
1558 {
1559  //For all active boundaries
1560  for (int localIndex : boundaryList_)
1561  {
1562  updateParticleVelocity(localIndex);
1563  }
1564 
1565  //For all active boundaries clear the data lists
1566  for (int localIndex : boundaryList_)
1567  {
1568  updateVelocityDataSend_[localIndex].clear();
1569  updateVelocityDataReceive_[localIndex].clear();
1570  }
1571 }
1572 
1573 
1579 {
1580  return domainHandler_;
1581 }
1582 
1589 {
1590  //Step 1: For every MPIDomain boundary, create a list of particles that have to be transmitted
1591  //queue send and receive instructions for the number of particles
1594 
1595  //Step 2: queue send and receive of data
1598 
1599  //Step 3: Add the received particles to the particleHandler of the current domain
1601 }
1602 
1611 {
1612  //Step1: check if the particle has to be sent to other processors
1615 
1616  //Step2: queue send and receive data. Note for an inserted particle, no interactions should be required
1619 
1620  //Step3: Add the received particles to the particleHandler of the current domain
1623 }
1624 
1633 void Domain::updateStatus(std::set<BaseParticle*>& ghostParticlesToBeDeleted)
1634 {
1635  //Collect new positions and velocities and send them to the other domains
1638 
1639  //Receive the new positions and velocities from other domains
1640  //and update the mpi flagged particles accordingly. removes and switched particles in the lists
1641  finalisePositionAndVelocityUpdate(ghostParticlesToBeDeleted);
1643 }
1644 
1649 {
1650  //collect new velocity data and send
1653 
1654  //process the received data
1657 }
1658 
1665 {
1666  unsigned int count = 0;
1667  for (auto& index : boundaryParticleListNeighbour_)
1668  {
1669  count += index.size();
1670  }
1671  return count;
1672 }
1673 
1680 {
1681  unsigned int count = 0;
1682  for (auto& index : boundaryParticleListNeighbour_)
1683  {
1684  for (auto& p : index)
1685  {
1686  if (!p->isPeriodicGhostParticle())
1687  {
1688  count++;
1689  }
1690  }
1691  }
1692  return count;
1693 }
1694 
1698 void Domain::flushParticles(std::set<BaseParticle*>& toBeFlushedList)
1699 {
1700  //For all active boundaries
1701  for (int localIndex : boundaryList_)
1702  {
1703  flushParticlesFromList(boundaryParticleList_[localIndex], toBeFlushedList);
1704  flushParticlesFromList(boundaryParticleListNeighbour_[localIndex], toBeFlushedList);
1705  }
1706 }
1707 
1708 void Domain::flushParticlesFromList(std::vector<BaseParticle*>& list, std::set<BaseParticle*>& toBeFlushedList)
1709 {
1710  //Firstly: turn all particles that need to be flushed into nullptrs
1711  for (auto& p : list)
1712  {
1713  if (p != nullptr)
1714  {
1715  BaseParticle* particle1 = p;
1716  for (BaseParticle* particle2 : toBeFlushedList)
1717  {
1718  //If the particle was found in the list, make a nullptr
1719  if (particle1 == particle2)
1720  {
1721  logger(VERBOSE, "Removing particle from mpi domain at: %", particle1->getPosition());
1722  p = nullptr;
1723  }
1724  }
1725  }
1726  }
1727 }
1728 
1734 {
1735  return middle_;
1736 }
1737 
1743 {
1744  //For all active boundaries
1745  for (int i : boundaryList_)
1746  {
1749  }
1750 }
1751 
1757 void Domain::cleanCommunicationList(std::vector<BaseParticle*>& list)
1758 {
1759  for (int i = 0; i < list.size(); i++)
1760  {
1761  if (list[i] == nullptr)
1762  {
1763  list[i] = list.back();
1764  list.pop_back();
1765  i--;
1766  }
1767  }
1768 }
void prepareBoundaryDataTransmission()
Prepares the MPI transmission of particle and interaction data from particles in particleHandler.
Definition: Domain.cc:1146
unsigned int getId() const
Returns the unique identifier of any particular object.
Definition: BaseObject.h:125
std::vector< void * > interactionDataSend_
Container that keeps a void array of all the interaction data that are being send to other domains...
Definition: Domain.h:572
std::vector< unsigned > numberOfParticlesReceive_
Counter that keeps track of the number of particles that are being received by this domain...
Definition: Domain.h:527
bool isInInnerDomain(BaseParticle *particle)
Check if the particle is in the current domain but not in the communication zone. ...
Definition: Domain.cc:430
void disableBoundary(unsigned localIndex)
Disables a boundary of the domain with a neighbouring domain.
Definition: Domain.cc:378
std::vector< bool > activeBoundaryList_
A list of flags corresponding to an inactive or active boundary.
Definition: Domain.h:491
void updateParticleVelocity(int localIndex)
Updates the velocity of particles which are flagged as MPIParticles.
Definition: Domain.cc:1441
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
static MPIContainer & Instance()
fetch the instance to be used for communication
Definition: MpiContainer.h:134
std::vector< int > getLocalIndexInverse(int localMeshIndex)
Definition: Domain.cc:515
Mdouble X
the vector components
Definition: Vector.h:65
void updateParticles(std::set< BaseParticle * > &ghostParticlesToBeDeleted)
This step updates all communication lists and particles in the communication zone.
Definition: Domain.cc:1308
void setBounds(std::vector< double > domainLeft, std::vector< double > domainRight, bool computeMiddle)
Sets the domain bounds.
Definition: Domain.cc:268
void addParticle(BaseParticle *particle)
Initialises a single particle which is added during the simulation.
Definition: Domain.cc:1610
int getGlobalIndex()
Gets the global index of the domain.
Definition: Domain.cc:346
bool inBoundary(BaseParticle *particle, int localIndex)
Definition: Domain.cc:578
virtual void hGridGetInteractingParticleList(BaseParticle *obj, std::vector< BaseParticle * > &list)
Creates a list of neighbour particles obtained from the hgrid.
Definition: DPMBase.h:961
void setInMPIDomain(bool flag)
Flags the status of the particle if wether it is in the communication zone or not.
void collectInteractionData(int localIndex)
Collects the data of an interaction that has to be communicated to other processors.
Definition: Domain.cc:924
int getRank()
Gets the rank associated with the assigned processorID.
Definition: Domain.cc:328
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
void constructor()
contructor of a domain
Definition: Domain.cc:133
bool isInNewBoundaryParticleList(BaseParticle *object, int localIndex) const
Definition: Domain.cc:809
void * createMPIInteractionDataArray(unsigned int numberOfInteractions) const
creates an empty MPIInteractionDataArray
BaseInteraction * getInteractionWith(BaseParticle *P, unsigned timeStamp, InteractionHandler *interactionHandler) override
Checks if particle is in interaction with given particle P, and if so, returns vector of pointer to t...
T * getObjectById(const unsigned int id)
Gets a pointer to the Object at the specified index in the BaseHandler.
Definition: BaseHandler.h:565
std::vector< unsigned > numNewInteractionsSend_
Counter that keeps track of the number of interactions that are being send to other domains...
Definition: Domain.h:532
std::vector< int > findNearbyBoundaries(BaseParticle *particle, Mdouble offset=0)
This function finds if a given particle is close to a given boundary.
Definition: Domain.cc:554
Vec3D getMin() const
Definition: DPMBase.h:637
Vec3D middle_
Middle of the closed domain.
Definition: Domain.h:463
std::vector< double > domainMin_
Minimum domain bounds in the x,y and z direction.
Definition: Domain.h:453
std::vector< unsigned > getGlobalMeshIndex()
Gets the global mesh index of the domain.
Definition: Domain.cc:355
It is an abstract base class due to the purely virtual functions declared below. Even if the function...
Definition: BaseObject.h:50
std::vector< int > localIndexToGlobalIndexTable_
look-up table to get the global index given a local domain index
Definition: Domain.h:480
void debugInformation()
Definition: Domain.cc:1096
const std::complex< Mdouble > i
Definition: ExtendedMath.h:51
void flushParticles(std::set< BaseParticle * > &toBeDeletedList)
Particles that are going to be deleted from the simulation are flushed out of the communication bound...
Definition: Domain.cc:1698
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
virtual BaseInteraction * getInteractionWith(BaseParticle *P, unsigned timeStamp, InteractionHandler *interactionHandler)=0
Returns the interaction between this object and a given BaseParticle.
void addNewParticles()
Initialises the MPIParticles by communicating newly found particles.
Definition: Domain.cc:1588
void addGhostObject(int fromPrcessor, int toProcessor, BaseParticle *p)
Adds a ghost particle located at fromProcessor to toProcessor.
void setCommunicationComplexity(unsigned complexity)
Set the communication complexity of the particle.
BaseParticle * findParticleInList(unsigned int identification, std::vector< BaseParticle * > particleList)
Searches for a particle with a specific id in a list of particles.
Definition: Domain.cc:529
Mdouble getSumOfInteractionRadii(const BaseParticle *particle) const
returns the sum of the radii plus the interactionDistance
Definition: BaseParticle.h:379
std::vector< int > localIndexToProcessorList_
look-up table to get the processor of the domain given a local domain index
Definition: Domain.h:486
bool containsParticle(BaseParticle *particle, Mdouble offset=0.0)
Check to see if a given particle is within the current domain.
Definition: Domain.cc:400
Domain()
Default Domain constructor.
Definition: Domain.cc:48
void findNewMPIInteractions()
Finds interactions that have to be send over to another domain.
Definition: Domain.cc:824
void performBoundaryDataTransmission()
Collects data to be transmitted and then performs the transmission of the data.
Definition: Domain.cc:1221
std::vector< std::vector< BaseInteraction * > > newInteractionList_
Array that queues interactions that need to be transmitted.
Definition: Domain.h:517
Mdouble getInteractionDistance()
Gets the interaction distance of the domain handler.
void deleteMPIInteractionDataArray(void *dataArray)
deletes an MPIInteractionDataArray
bool isInGreaterDomain(BaseParticle *particle)
Check to see if a given particle is in the current domain or in neighbouring communication zones...
Definition: Domain.cc:420
void finaliseVelocityUpdate()
Processes particle velocity data for ghost particles.
Definition: Domain.cc:1557
std::vector< std::vector< MPIParticleVelocity > > updateVelocityDataReceive_
Container that keeps a list of MPIParticleVelocities that are being received by this domain...
Definition: Domain.h:567
void updateStatus(std::set< BaseParticle * > &ghostParticlesToBeDeleted)
Updates particles that are not in the current domain and communicates newly added particles...
Definition: Domain.cc:1633
Mdouble getLengthSquared() const
Calculates the squared length of this Vec3D: .
Definition: Vector.cc:184
const Mdouble inf
Definition: GeneralDefine.h:44
Stores information about interactions between two interactable objects; often particles but could be ...
std::vector< unsigned > globalMeshIndex_
Vector containing the global mesh indices i,j,k.
Definition: Domain.h:474
std::vector< int > boundaryList_
A list of indices of all the active boundaries.
Definition: Domain.h:497
unsigned int getNumberOfTrueMPIParticles()
Obtains the number of particles in the particleHandler that are MPIParticles, but NOT periodic partic...
Definition: Domain.cc:1679
void setComponent(int index, double val)
Sets the requested component of this Vec3D to the requested value.
Definition: Vector.cc:217
void copyDataFromMPIParticleToParticle(MPIParticle *bP, BaseParticle *p, ParticleHandler *particleHandler)
Copies data from an MPIParticle class to a BaseParticle and sets the particleHandler and species...
std::vector< void * > interactionDataReceive_
Container that keeps a void array of all the interaction data that is being received by this domain...
Definition: Domain.h:577
void setRange(Direction direction, Mdouble min, Mdouble max)
Sets the domain range in a given direction.
Definition: Domain.cc:219
void getInteractionDetails(void *interactionData, unsigned int index, unsigned int &identificationP, unsigned int &identificationI, bool &isWallInteraction, unsigned &timeStamp)
reads the basic interaction details from an MPIInteractionDataArray
void updateParticlePosition(int localIndex)
Updates the position of particles which are flagged as MPIParticles.
Definition: Domain.cc:1412
void write(std::ostream &os) const override
This function does nothing.
Definition: Domain.cc:197
std::vector< std::vector< MPIParticlePosition > > updatePositionDataSend_
Container that keeps a list of MPIParticlePositions that are being send to other domains.
Definition: Domain.h:552
const int intMax
Definition: GeneralDefine.h:45
#define MAX_PROC
Definition: GeneralDefine.h:51
void processReceivedInteractionData(unsigned index, std::vector< BaseParticle * > &newParticles)
Processes the received interactions from newly added mpi particles.
Definition: Domain.cc:999
static BaseParticle * newParticle()
MPIParticle copyDataFromParticleToMPIParticle(BaseParticle *p)
Copies data from a SuperQuadricParticle to an MPIParticle class and returns this. ...
void cleanCommunicationLists()
Removes nullptrs from boundaryParticleList_ and boundaryParticleListNeighbour_.
Definition: Domain.cc:1742
void preparePositionAndVelocityUpdate()
Function that sends particle position and velocity data for ghost particles to other processors...
Definition: Domain.cc:1460
bool isInMPIDomain()
Indicates if the particle is in the communication zone of the mpi domain.
virtual bool getDistanceAndNormal(const BaseParticle &P, Mdouble &distance, Vec3D &normal_return) const =0
Pure virtual function that computes the distance of a BaseParticle to this wall and returns the norma...
void finalisePositionAndVelocityUpdate(std::set< BaseParticle * > &ghostParticlesToBeDeleted)
processes position and velocity data for ghost particles
Definition: Domain.cc:1501
void findNewMPIParticles(const ParticleHandler &particleHandler)
Function that finds new particles in the particle handler that should be added to the communication l...
Definition: Domain.cc:786
void setPreviousPosition(const Vec3D &pos)
Sets the particle's position in the previous time step.
MercuryMPITag
An enum that facilitates the creation of unique communication tags in the parallel code...
Definition: MpiContainer.h:76
std::vector< double > getDomainMin()
Gets the minimum domain bounds.
Definition: Domain.cc:310
Mdouble getComponent(int index) const
Returns the requested component of this Vec3D.
Definition: Vector.cc:194
std::vector< std::vector< BaseParticle * > > boundaryParticleListNeighbour_
a list of ghost particles on the current domain, which are real on the neighbour domain ...
Definition: Domain.h:507
ParticleHandler particleHandler
An object of the class ParticleHandler, contains the pointers to all the particles created...
Definition: DPMBase.h:1395
void cleanCommunicationList(std::vector< BaseParticle * > &list)
Removes nullptr's from a given particle list.
Definition: Domain.cc:1757
Container to store Interaction objects.
T * getObject(const unsigned int id)
Gets a pointer to the Object at the specified index in the BaseHandler.
Definition: BaseHandler.h:613
std::vector< std::vector< MPIParticleVelocity > > updateVelocityDataSend_
Container that keeps a list of MPIParticleVelocities that are being send to other domains...
Definition: Domain.h:562
DomainHandler * getHandler() const
Gets the domainHandler.
Definition: Domain.cc:1578
void sendAndReceiveMPIData(MercuryMPITag tag, MercuryMPIType type, T *receiveData, unsigned receiveCount, T *sendData, unsigned sendCount, unsigned localIndexNeighbour)
Function that sends transmissionData/positionData/velocityData to other processors.
Definition: Domain.h:369
void sendAndReceiveCount(MercuryMPITag tag, unsigned &countReceive, unsigned &countSend, unsigned localIndexNeighbour)
A symmetric communication between two domains exchanging a send/recieve count.
Definition: Domain.cc:1119
std::vector< unsigned > numberOfParticlesSend_
Counter that keeps track of the number of particles that are being send to other domains.
Definition: Domain.h:522
void processReceivedBoundaryParticleData(unsigned index, std::vector< BaseParticle * > &newParticles)
Function that copies the mpi data format of a base particle to a real particle and adds it to the par...
Definition: Domain.cc:941
DomainHandler * domainHandler_
Pointer to the domain's DomainHandler container.
Definition: Domain.h:448
void collectBoundaryParticleData(int localIndex)
collects the data of a particle that has to be communicated to other processors
Definition: Domain.cc:909
void read(std::istream &is) override
This function does nothing.
Definition: Domain.cc:186
Basic class for walls.
Definition: BaseWall.h:47
std::vector< double > getDomainMax()
Gets the maximum domain bounds.
Definition: Domain.cc:319
void finaliseBoundaryDataTransmission()
This function processes the transmitted data.
Definition: Domain.cc:1262
void setRank(int rank)
Sets the rank associated with the assigned processorID.
Definition: Domain.cc:301
void findNewMPIParticle(BaseParticle *particle)
Function that check if a given particle should be added to the communication lists.
Definition: Domain.cc:800
InteractionHandler interactionHandler
An object of the class InteractionHandler.
Definition: DPMBase.h:1425
Vec3D getMax() const
Definition: DPMBase.h:643
bool isPeriodicGhostParticle() const
Indicates if this particle is a ghost in the periodic boundary.
Container to store all BaseParticle.
Mdouble Y
Definition: Vector.h:65
MPIParticlePosition copyPositionFrom(BaseParticle *particle)
Copies the position from a particle to an MPIParticlePosition class.
~Domain() override
Destructor, destroys the domain.
Definition: Domain.cc:121
WallHandler wallHandler
An object of the class WallHandler. Contains pointers to all the walls created.
Definition: DPMBase.h:1405
std::vector< std::vector< BaseParticle * > > newBoundaryParticleList_
Array that queues particles that need to be transmitted.
Definition: Domain.h:512
The simulation can be subdivided into Domain's used in parallel code.
Definition: Domain.h:63
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
Vec3D getMiddle() const
Gives the middle of the domain.
Definition: Domain.cc:1733
void createLookUpTable()
Create a look up table between local index system to global index system.
Definition: Domain.cc:451
bool isInCommunicationZone(BaseParticle *particle)
Check if the particle is in the communication zone of the current domain.
Definition: Domain.cc:441
unsigned getCommunicationComplexity()
Obtains the communication complexity of the particle.
void setMPIParticle(bool flag)
Flags the mpi particle status.
void disableBoundaries()
disables all domain boundaries that have no neighbour
Definition: Domain.cc:597
Container to store all Domain.
Definition: DomainHandler.h:46
static Mdouble getLengthSquared(const Vec3D &a)
Calculates the squared length of a Vec3D: .
Definition: Vector.h:316
void updateVelocity()
Updates MPI particle velocity at the half-time step.
Definition: Domain.cc:1648
std::string getName() const override
Returns the name of the object.
Definition: Domain.cc:207
Direction
An enum that indicates the direction in Cartesian coordinates.
Definition: GeneralDefine.h:76
std::vector< std::vector< MPIParticlePosition > > updatePositionDataReceive_
Container that keeps a list of MPIParticlePositions that are being received by this domain...
Definition: Domain.h:557
int globalIndex_
Global index of the domain in the mesh.
Definition: Domain.h:469
std::vector< unsigned > getNumberOfDomains()
Gets the number of domains in the domain handler.
int rank_
Rank of the domain which identifies to which processor it belongs.
Definition: Domain.h:582
std::vector< double > domainMax_
Maximum domain bounds in the x,y and z direction.
Definition: Domain.h:458
Defines the basic properties that a interactable object can have.
std::vector< std::vector< MPIParticle > > boundaryParticleDataReceive_
Container that keeps a list of MPIParticles that are being received by this domain.
Definition: Domain.h:547
virtual void setMPIInteraction(void *interactionDataArray, unsigned int index, bool resetPointers)
virtual void hGridUpdateMove(BaseParticle *, Mdouble)
Definition: DPMBase.cc:1924
void prepareVelocityUpdate()
Function that sends particle velocity data for ghost particles.
Definition: Domain.cc:1529
virtual Domain * copy() const
Function that creates a copy of this current domain, using the copy constructor.
Definition: Domain.cc:175
void flushParticlesFromList(std::vector< BaseParticle * > &list, std::set< BaseParticle * > &toBeDeletedList)
Particles that are going to be deleted from the simulation are flushed out of a give communcation bou...
Definition: Domain.cc:1708
T * getLastObject()
Gets a pointer to the last Object in this BaseHandler.
Definition: BaseHandler.h:634
std::vector< std::vector< BaseParticle * > > boundaryParticleList_
A list of boundary particles in the communication zone that are ghost particles on other domains...
Definition: Domain.h:502
void processSentBoundaryParticles(unsigned index)
Bookkeep the newly send particles.
Definition: Domain.cc:977
Definition: Vector.h:49
int getLocalIndex(int i, int j, int k)
return the local index of a domain given local mesh indices i,j and k
Definition: Domain.cc:498
MPIParticleVelocity copyVelocityFrom(BaseParticle *particle)
Copies the velocity from a particle to an MPIParticleVelocity class.
DPMBase * getDPMBase()
Gets the problem that is solved using this handler.
Definition: BaseHandler.h:725
std::vector< unsigned > numNewInteractionsReceive_
Counter that keeps track of the number of interactions that are being received by this domain...
Definition: Domain.h:537
Mdouble Z
Definition: Vector.h:65
Mdouble getTime() const
Returns the current simulation time.
Definition: DPMBase.cc:805
void setHandler(DomainHandler *handler)
Sets the domainHandler.
Definition: Domain.cc:337
void addParticlesToLists(BaseParticle *particle, std::vector< std::vector< BaseParticle * > > &list)
Function that adds the particles to the approriate boundary list.
Definition: Domain.cc:701
std::vector< bool > getActiveBoundaryList()
Returns a list of boundaries that are active in mpi communication.
Definition: Domain.cc:387
unsigned int getNumberOfMPIParticles()
Obtains the number of particles in the particleHandler that are MPIParticles.
Definition: Domain.cc:1664
void setGlobalMeshIndex(std::vector< unsigned > globalMeshIndex)
Sets the global mesh index of theh domain.
Definition: Domain.cc:365
std::vector< std::vector< MPIParticle > > boundaryParticleDataSend_
Container that keeps a list of MPIParticles that are being send to other domains. ...
Definition: Domain.h:542