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 
524 BaseParticle* Domain::findParticleInList(unsigned int identification, std::vector<BaseParticle*> particleList)
525 {
526  for (BaseParticle* particle : particleList)
527  {
528  if (particle->getId() == identification)
529  {
530  return particle;
531  }
532  }
533  return nullptr;
534 }
535 
549 std::vector<int> Domain::findNearbyBoundaries(BaseParticle* particle, Mdouble offset)
550 {
551  std::vector<int> boundaryIndex(3);
552  Mdouble interactionDistance = domainHandler_->getInteractionDistance();
553 
554  //for x,y,z directions, find if the particle is close to the left or right wall or not at all
555  for (int d = 0; d < 3; d++)
556  {
557  //Check if the particle is close to Lx
558  if ((particle->getPosition().getComponent(d)) < (domainMin_[d] + interactionDistance + offset))
559  {
560  //Update switch
561  boundaryIndex[d] = -1;
562  }
563  //Check if particle is close to Rx
564  else if ((particle->getPosition().getComponent(d)) >= (domainMax_[d] - interactionDistance - offset))
565  {
566  //Update switch
567  boundaryIndex[d] = 1;
568  }
569  }
570  return boundaryIndex;
571 }
572 
580 {
581  std::vector<unsigned> numberOfDomains = domainHandler_->getNumberOfDomains();
582  std::vector<int> localMeshIndex(3);
583  int localIndex;
584 
585  //disble own list
586  activeBoundaryList_[13] = false;
587 
588  //Special case when numberOfDomains[d] = 1 -> disable all d direction boundaries
589  for (unsigned d = 0; d < 3; d++) //Loop over all dimensions
590  {
591  //If there is only one domain in the d-direction, take action
592  if (numberOfDomains[d] == 1)
593  {
594  //Loop over all locaIndices
595  for (int i = -1; i < 2; i++)
596  {
597  localMeshIndex[Direction::XAXIS] = i;
598  for (int j = -1; j < 2; j++)
599  {
600  localMeshIndex[Direction::YAXIS] = j;
601  for (int k = -1; k < 2; k++)
602  {
603  localMeshIndex[Direction::ZAXIS] = k;
604  //Disable all boundaries that are not in the "middle" of the mesh in the d-direction
605  // i.e. a localMeshIndex[d] = -1 is a neighbour in the d-direction, however there is only one
606  // element in that direction so it can't be a real neighbour.
607  if (localMeshIndex[d] != 0)
608  {
609  //Disable the boundary
610  localIndex = getLocalIndex(localMeshIndex);
611  activeBoundaryList_[localIndex] = false;
612  }
613  }
614  }
615  }
616  }
617  }
618 
619  //Compute the boundaryIndex for all other cases where numberOfDomains is larger than 1
620  std::vector<int> boundaryIndex(3);
621  for (int d = 0; d < 3; d++)
622  {
623  //Check if the domain is on the simulationDomainMin boundary
624  if (globalMeshIndex_[d] == 0)
625  {
626  //Update boundary index
627  boundaryIndex[d] = -1;
628  }
629  //Check if the domain is on the simulationDomainMax boundary
630  else if (globalMeshIndex_[d] == (numberOfDomains[d] - 1))
631  {
632  //Update boundary index
633  boundaryIndex[d] = 1;
634  }
635  }
636 
637  //Disable the boundaries that are on the edge of the simulation domain
638  for (int d = 0; d < 3; d++)
639  {
640  //If boundary is on the edge of the simulation domain
641  if (boundaryIndex[d] != 0)
642  {
643  //Loop over all local indices
644  for (int i = -1; i < 2; i++)
645  {
646  for (int j = -1; j < 2; j++)
647  {
648  for (int k = -1; k < 2; k++)
649  {
650  localMeshIndex = {i, j, k};
651  //Set the localMeshIndex to the boundaryIndex in the d-direction
652  //Since we are only interested in the domains on this boundary
653  localMeshIndex[d] = boundaryIndex[d];
654  //Disable the boundary
655  localIndex = getLocalIndex(localMeshIndex);
656  activeBoundaryList_[localIndex] = false;
657  }
658  }
659  }
660  }
661  }
662 
663  //Create a list of the active boundaries
664  localIndex = 0;
665  for (bool active : activeBoundaryList_)
666  {
667  if (active)
668  {
669  boundaryList_.push_back(localIndex);
670  }
671  localIndex++;
672  }
673 }
674 
683 void Domain::addParticlesToLists(BaseParticle* particle, std::vector<std::vector<BaseParticle*> >& list)
684 {
685  std::vector<int> boundaryIndex = findNearbyBoundaries(particle);
686 
687  //Compute and set complexity of the particle
688  unsigned int complexity = boundaryIndex[0] + 3 * boundaryIndex[1] + 9 * boundaryIndex[2] + 13;
689  unsigned int list_complexity = 0;
690  for (int d = 0; d < 3; d++) //Loop over all directions
691  {
692  list_complexity += std::abs(boundaryIndex[d]);
693  }
694  particle->setCommunicationComplexity(complexity);
695  //particle->setCommunicationComplexity(list_complexity);
696 
697  //Based on the complexity of the particle, add it to the approriate list
698  switch (list_complexity)
699  {
700  //The particle is not close at all
701  case 0:
702  break;
703  //The particle is close to one side
704  // 1 side contribution
705  case 1 :
706  //Add the side contribution
707  list[getLocalIndex(boundaryIndex)].push_back(particle);
708  break;
709  //The particle is close to two neighbouring directions
710  //2 side and 1 rib contrubution
711  case 2 :
712  {
713  //Add the two side contributions
714  for (int d = 0; d < 3; d++)
715  {
716  std::vector<int> localMeshIndex = {0, 0, 0};
717  localMeshIndex[d] = boundaryIndex[d];
718  //Avoid adding the particle in the wrong direction by excluding localMeshIndex[d] = 0
719  if (localMeshIndex[d] != 0)
720  {
721  list[getLocalIndex(localMeshIndex)].push_back(particle);
722  }
723  }
724  }
725  //Add the rib contribution
726  list[getLocalIndex(boundaryIndex)].push_back(particle);
727  break;
728 
729  //The particle is close to three neighbouring directions
730  // 3 side, 3 rib and 1 corner contribution
731  case 3 :
732  {
733  //Add the three side contributions
734  for (int d = 0; d < 3; d++)
735  {
736  std::vector<int> localMeshIndex = {0, 0, 0};
737  //Reset index vector
738  localMeshIndex[d] = boundaryIndex[d];
739  list[getLocalIndex(localMeshIndex)].push_back(particle);
740  }
741 
742  //Add the three rib contributions
743  for (int d = 0; d < 3; d++)
744  {
745  std::vector<int> localMeshIndex = boundaryIndex;
746  //All rib boundary indices are given by the boundaryIndex and setting one of the components to zero
747  localMeshIndex[d] = 0;
748  list[getLocalIndex(localMeshIndex)].push_back(particle);
749  }
750  }
751 
752  //Add the corner contribution
753  list[getLocalIndex(boundaryIndex)].push_back(particle);
754  break;
755 
756  default :
757  logger(INFO, "boundaryIndex : %,%,% | list_complexity: %", boundaryIndex[0], boundaryIndex[1], boundaryIndex[2],
758  list_complexity);
759  logger(ERROR, "Particle is in contact with the wrong number of boundaries");
760  break;
761  }
762 }
763 
768 void Domain::findNewMPIParticles(const ParticleHandler& particleHandler)
769 {
770  //For all particles inside the given domain, loop over all the particles to
771  //see if we have to add the particles to the new particle list
772  for (BaseParticle* particle : particleHandler)
773  {
774  findNewMPIParticle(particle);
775  }
776 }
777 
783 {
784  //If the particle is an MPI particle
785  if (!particle->isInMPIDomain() && !particle->isPeriodicGhostParticle())
786  {
788  }
789 }
790 
791 bool Domain::isInNewBoundaryParticleList(BaseParticle* objectI,int localIndex) const
792 {
793  for (BaseParticle *q : newBoundaryParticleList_[localIndex]) {
794  if (q == objectI)
795  return true;
796  }
797  return false;
798 }
799 
807  {
808  //For all active boundaries
809  for (int localIndex : boundaryList_)
810  {
811  //Check all newly added particles for interactions with already known particles
812  for (BaseParticle* newBoundaryParticle : newBoundaryParticleList_[localIndex])
813  {
814  //Loop over all interactions of the new particle
815  std::vector<BaseInteraction*> interactions = newBoundaryParticle->getInteractions();
816  for (BaseInteraction* interaction : interactions)
817  {
818  //Find out what the new particle is interacting with
819  BaseParticle* particleP = dynamic_cast<BaseParticle*>(interaction->getP());
820  BaseParticle* objectI = dynamic_cast<BaseParticle*>(interaction->getI());
821 
822  //If the P in the interaction structure is the new particle, find I
823  if (newBoundaryParticle == particleP)
824  {
825  //Check if the new particle is interacting with a wall
826  if (!objectI)
827  {
828  //Check if the interaction is still valid
829  BaseWall* wall = dynamic_cast<BaseWall*>(interaction->getI());
830  Mdouble distance;
831  Vec3D normal;
832  wall->getDistanceAndNormal(*particleP, distance, normal);
833  if (distance <= particleP->getMaxInteractionRadius())
834  {
835  //Add the interaction to the list if there are still in contact (hence the if statement)
836  newInteractionList_[localIndex].push_back(interaction);
837  }
838  }
839  else //is I a particle
840  {
841  //check if particle is in mpi domain OR (TODO) if it is in the newBoundaryParticleList_
842  if (objectI->isInMPIDomain() || isInNewBoundaryParticleList(objectI, localIndex) ) {
843  //Is the particle still interacting after the position update?
844  Vec3D branchVector = particleP->getPosition() - objectI->getPosition();
845  //Get the square of the distance between particle i and particle j
846  Mdouble distanceSquared = Vec3D::getLengthSquared(branchVector);
847  Mdouble sumOfInteractionRadii =
848  objectI->getSumOfInteractionRadii(particleP);
849  if (distanceSquared < (sumOfInteractionRadii * sumOfInteractionRadii))
850  {
851  //Add the interaction to the list
852  newInteractionList_[localIndex].push_back(interaction);
853  }
854  }
855  }
856  }
857  else //newBoundaryParticle is I in the interaction, P can only be a particle
858  {
859  //it is "I" in the interaction structure, check if its in the mpi domain
860  if (particleP->isInMPIDomain() || isInNewBoundaryParticleList(particleP, localIndex))
861  {
862  //Is the particle still interacting after the position update?
863  Vec3D branchVector = particleP->getPosition() - objectI->getPosition();
864  //Get the square of the distance between particle i and particle j
865  Mdouble distanceSquared = Vec3D::getLengthSquared(branchVector);
866  Mdouble sumOfInteractionRadii =
867  objectI->getSumOfInteractionRadii(particleP);
868  if (distanceSquared < (sumOfInteractionRadii * sumOfInteractionRadii))
869  {
870  //Add the interaction to the list
871  newInteractionList_[localIndex].push_back(interaction);
872  }
873  }
874  }
875  }
876  }
877  }
878 }
879 
887 {
888  //For all particles
889  for (BaseParticle* particle : newBoundaryParticleList_[localIndex])
890  {
891  //Add the data to the transmission list
892  boundaryParticleDataSend_[localIndex].push_back(copyDataFromParticleToMPIParticle(particle));
893  }
894 }
895 
901 void Domain::collectInteractionData(int localIndex)
902 {
903  //Copy interactions to the data array
904  unsigned int indexInteraction = 0;
905  for (BaseInteraction* interaction : newInteractionList_[localIndex])
906  {
907  interaction->getMPIInteraction(interactionDataSend_[localIndex], indexInteraction);
908  indexInteraction++;
909  }
910 }
911 
918 void Domain::processReceivedBoundaryParticleData(const unsigned index, std::vector<BaseParticle*>& newParticles)
919 {
920 
921  ParticleHandler& particleHandler = getHandler()->getDPMBase()->particleHandler;
922  //for all MPIPparticles that have been received on the other side
923  for (int i = 0; i < numberOfParticlesReceive_[index]; i++)
924  {
925  //copy data from data
927  copyDataFromMPIParticleToParticle(&(boundaryParticleDataReceive_[index][i]), p0, &particleHandler);
928 
929  //Flag that the particle is a ghost particle of a real particle in the neighbour domain
930  p0->setMPIParticle(true);
931  //Flag that the particle is list as a particle in the mpi domain
932  p0->setInMPIDomain(true);
933  //add particle to handler
934  particleHandler.addGhostObject(p0);
935  //Set previous position as current position
936  BaseParticle* pGhost = particleHandler.getLastObject();
937  pGhost->setPreviousPosition(particleHandler.getLastObject()->getPosition());
938  //Add to the newParticle list
939  newParticles.push_back(pGhost);
940  logger(VERBOSE, "Adding mpi particle % at: %", particleHandler.getLastObject()->getId(),
941  particleHandler.getLastObject()->getPosition());
942 
943  //add pointer to the particle to the list and set the status to idle
944  boundaryParticleListNeighbour_[index].push_back(particleHandler.getLastObject());
945  }
946 }
947 
954 void Domain::processSentBoundaryParticles(const unsigned index)
955 {
956  for (BaseParticle* particle : newBoundaryParticleList_[index])
957  {
958  boundaryParticleList_[index].push_back(particle);
959  particle->setInMPIDomain(true);
960  }
961 }
962 
976 void Domain::processReceivedInteractionData(const unsigned localIndex, std::vector<BaseParticle*>& newParticles)
977 {
979  for (unsigned int l = 0; l < numNewInteractionsReceive_[localIndex]; l++)
980  {
981  unsigned int identificationP;
982  unsigned int identificationI;
983  bool isWallInteraction;
984  unsigned timeStamp;
985 
986  //Get the general information required to setup a new interaction
988  l, identificationP, identificationI,
989  isWallInteraction, timeStamp);
990 
991  logger(VERBOSE, "interaction details: %, %, idP %, idI %, wall %, time %", interactionDataReceive_[localIndex],
992  l, identificationP, identificationI,
993  isWallInteraction, timeStamp);
994 
995  //Find the particle in the newParticle list
996  BaseParticle* pGhost = nullptr;
997  int idOther;
998  for (BaseParticle* particle : newParticles)
999  {
1000  if (particle->getId() == identificationP)
1001  {
1002  pGhost = particle;
1003  idOther = identificationI;
1004  break;
1005  }
1006 
1007  if (particle->getId() == identificationI)
1008  {
1009  pGhost = particle;
1010  idOther = identificationP;
1011  break;
1012  }
1013  }
1014  if (pGhost == nullptr)
1015  {
1016  logger(WARN, "In Domain::processReceivedInteractionData: pGhost (id %) is nullptr, the interaction data is not copied. Two particles possibly moved into domain simultaneously.", identificationP);
1017  continue;
1018  }
1019 
1020  //If it is a wall interaction, do stuff
1021  if (isWallInteraction)
1022  {
1023  BaseInteractable* I = getHandler()->getDPMBase()->wallHandler.getObjectById(identificationI);
1024  //Create interactions
1025  BaseInteraction* j = I->getInteractionWith(pGhost, timeStamp, &iH);
1026  if (j!= nullptr) j->setMPIInteraction(interactionDataReceive_[localIndex], l, false);
1027 
1028  }
1029  else
1030  {
1031  //Obtain potential interaction particles
1032  std::vector<BaseParticle*> interactingParticleList;
1033  getHandler()->getDPMBase()->hGridGetInteractingParticleList(pGhost, interactingParticleList);
1034 
1035  //Find the other interacting particle
1036  BaseParticle* otherParticle = nullptr;
1037  for (BaseParticle* p2 : interactingParticleList)
1038  {
1039  if (p2->getId() == idOther)
1040  {
1041  otherParticle = p2;
1042  break;
1043  }
1044  }
1045  if (otherParticle == nullptr) {
1046  //search for the other particle in the newParticles list
1047  for (BaseParticle *particle : newParticles) {
1048  if (particle->getId() == idOther) {
1049  otherParticle = particle;
1050  break;
1051  }
1052  }
1053  if (otherParticle == nullptr) {
1054  logger(WARN,
1055  "In Domain::processReceivedInteractionData: otherParticle (id %) is nullptr, the interaction data with pGhost (id %) is not copied. Two particles possibly moved into domain simultaneously. nt = %",
1056  identificationI, identificationP, getHandler()->getDPMBase()->getNumberOfTimeSteps());
1057  continue;
1058  }
1059  }
1060  //Add the interaction
1061  BaseInteraction* j = pGhost->getInteractionWith(otherParticle, timeStamp, &iH);
1062  if (j!= nullptr) j->setMPIInteraction(interactionDataReceive_[localIndex], l, false);
1063  }
1064  }
1065 }
1066 
1068  std::stringstream s;
1069  std::stringstream m;
1070  for (auto p : getHandler()->getDPMBase()->particleHandler) {
1071  if (p->isMPIParticle())
1072  m << std::setw(4) << p->getId();
1073  else
1074  s << std::setw(4) << p->getId();
1075  }
1076  logger(INFO,"Particles %, MPI %",s.str(),m.str());
1077 }
1078 
1079 
1090  void Domain::sendAndReceiveCount(MercuryMPITag tag, unsigned& countReceive, unsigned& countSend,
1091  unsigned localIndexNeighbour)
1092  {
1093  int globalIndexNeighbour = localIndexToGlobalIndexTable_[localIndexNeighbour];
1094  int processor = localIndexToProcessorList_[localIndexNeighbour];
1095 
1096  //Create communication tags
1097  int tagReceive = globalIndexNeighbour * MAX_PROC + globalIndex_ * 10 + tag;
1098  int tagSend = globalIndex_ * MAX_PROC + globalIndexNeighbour * 10 + tag;
1099 
1100  logger.assert(tagSend > 0, "Send tag is wrong. Tag: %", tagSend);
1101  logger.assert(tagReceive > 0, "Receive tag is wrong. Tag: %", tagReceive);
1102  logger.assert(processor >= 0, "Processor is wrong. processor: %", processor);
1103 
1104  //Communicate the requests
1105  MPIContainer::Instance().receive(countReceive, processor, tagReceive);
1106  MPIContainer::Instance().send(countSend, processor, tagSend);
1107  }
1108 
1118  {
1119  //Find new particles that have entered the communication zone
1120  findNewMPIParticles(getHandler()->getDPMBase()->particleHandler);
1121 
1122  //Find interactions between new particles and other particles in the communication zone
1124 
1125  //Compute number of particles to send
1126  //For all boundaries
1127  for (int localIndex : boundaryList_)
1128  {
1129  numberOfParticlesSend_[localIndex] = newBoundaryParticleList_[localIndex].size();
1130  numNewInteractionsSend_[localIndex] = newInteractionList_[localIndex].size();
1131  }
1132 
1133  //For all active boundaries
1134  for (int localIndex : boundaryList_)
1135  {
1136  //Send and receive the number of new boundary particles
1138  (numberOfParticlesSend_[localIndex]), localIndex);
1139 
1140  //Send and receive the new interactions
1142  (numNewInteractionsSend_[localIndex]), localIndex);
1143  }
1144  }
1145 
1155  {
1156  //Assign the particle to the correct lists if the particle belongs to this domain
1157  if (containsParticle(particle))
1158  {
1159  findNewMPIParticle(particle);
1160  }
1161 
1162  //Note: When inserting a single particle, it has no interactions, so no interactions to be copied either
1163 
1164  //Compute number of particles to send
1165  //For all boundaries
1166  for (int localIndex : boundaryList_)
1167  {
1168  numberOfParticlesSend_[localIndex] = newBoundaryParticleList_[localIndex].size();
1169  numNewInteractionsSend_[localIndex] = newInteractionList_[localIndex].size();
1170  }
1171 
1172 
1173  //For all active boundaries
1174  for (int localIndex : boundaryList_)
1175  {
1176  //Send and recieve the number of new boundary particles
1178  (numberOfParticlesSend_[localIndex]), localIndex);
1179 
1180  //Send and receive the new interactions
1182  (numNewInteractionsSend_[localIndex]), localIndex);
1183  }
1184  }
1185 
1193  {
1194  //For all active boundaries
1195  for (int localIndex : boundaryList_)
1196  {
1197  //make sure enough memory is reserved to receive the data
1198  boundaryParticleDataReceive_[localIndex].resize(numberOfParticlesReceive_[localIndex]);
1199 
1200  //Collect the particle data
1201  collectBoundaryParticleData(localIndex);
1202 
1203  //Send the data
1205  boundaryParticleDataReceive_[localIndex].data(), numberOfParticlesReceive_[localIndex],
1206  boundaryParticleDataSend_[localIndex].data(), numberOfParticlesSend_[localIndex],
1207  localIndex);
1208 
1209 
1210  //create arrays for sending and receiving interaction data
1212  numNewInteractionsSend_[localIndex]);
1214  numNewInteractionsReceive_[localIndex]);
1215 
1216  //Collect the interaction data
1217  collectInteractionData(localIndex);
1218 
1219  //Send the data
1221  interactionDataReceive_[localIndex], numNewInteractionsReceive_[localIndex],
1222  interactionDataSend_[localIndex], numNewInteractionsSend_[localIndex], localIndex);
1223  }
1224  }
1225 
1234  {
1235  //For all boundaries
1236  for (int localIndex : boundaryList_)
1237  {
1238  //copy the received data into the particleHandler and neighbour particle list
1239  std::vector<BaseParticle*> newParticles;
1240  processReceivedBoundaryParticleData(localIndex, newParticles);
1241  //All particles in the current domain that have been send to the other domains need to be flagged as communicating particles
1242  processSentBoundaryParticles(localIndex);
1243  //copy the received interaction data into the standard dpm structure
1244  processReceivedInteractionData(localIndex, newParticles);
1245 
1246  //Delete all send/receive data
1247  boundaryParticleDataSend_[localIndex].clear();
1248  boundaryParticleDataReceive_[localIndex].clear();
1251  interactionDataReceive_[localIndex]);
1252 
1253  //Reset all counters
1254  numberOfParticlesSend_[localIndex] = 0;
1255  numberOfParticlesReceive_[localIndex] = 0;
1256  numNewInteractionsSend_[localIndex] = 0;
1257  numNewInteractionsReceive_[localIndex] = 0;
1258 
1259  //Reset all lists
1260  newBoundaryParticleList_[localIndex].clear();
1261  newInteractionList_[localIndex].clear();
1262  }
1263  }
1264 
1279  void Domain::updateParticles(std::set<BaseParticle*>& ghostParticlesToBeDeleted)
1280  {
1281  int complexityNew;
1282  std::vector<int> boundaryIndex;
1283 
1284  //For all active boundaries
1285  for (int localIndex : boundaryList_)
1286  {
1287  //Step 1A: Remove the particles from the boundaryParticleList_ which require re-assignment
1288  //or have just moved away from the region
1289  for (int p = 0; p < boundaryParticleList_[localIndex].size(); p++)
1290  {
1291  BaseParticle* particle = boundaryParticleList_[localIndex][p];
1292  //Check if the particle is still in the domain, but not in the communication zone
1293  if (containsParticle(particle))
1294  {
1295  //check if the complexity has changed
1296  boundaryIndex = findNearbyBoundaries(particle);
1297  complexityNew = boundaryIndex[0] + 3 * boundaryIndex[1] + 9 * boundaryIndex[2] + 13;
1298  if (particle->getCommunicationComplexity() != complexityNew)
1299  {
1300  logger(VERBOSE, "time: % | global index: % in list % | particle % | CURRENT DOMAIN - CHANGES "
1301  "COMPLEXITY", domainHandler_->getDPMBase()->getTime(), globalIndex_,
1302  localIndexToGlobalIndexTable_[localIndex], particle->getId());
1303  //Flag the particle that it no longer participates in the communication layer
1304  //so it can be reintroduced in the transmission step
1305  particle->setMPIParticle(false);
1306  particle->setInMPIDomain(false);
1307  particle->setCommunicationComplexity(0);
1308  boundaryParticleList_[localIndex][p] = nullptr;
1309  }
1310  }
1311  else
1312  {
1313  logger(VERBOSE, "time: % | global index: % in list % | particle % | CURRENT DOMAIN - TO NEIGHBOURING "
1314  "DOMAIN", domainHandler_->getDPMBase()->getTime(), globalIndex_,
1315  localIndexToGlobalIndexTable_[localIndex], particle->getId());
1316 
1317  ghostParticlesToBeDeleted.insert(particle);
1318  boundaryParticleList_[localIndex][p] = nullptr;
1319  }
1320  }
1321 
1322  //Step 1B: Remove the particles from the boundaryParticleListNeightbour_ which require re-assignment
1323  //or have just moved away from the region
1324  for (int p = 0; p < boundaryParticleListNeighbour_[localIndex].size(); p++)
1325  {
1326  BaseParticle* particle = boundaryParticleListNeighbour_[localIndex][p];
1327  //Check if the particle moved out of the neighbour domain
1329  {
1330  //The particle has moved to this domain
1331  if (containsParticle(particle))
1332  {
1333  logger(VERBOSE,
1334  "time: % | global index: % in list % | particle % | NEIGHBOURING DOMAIN - TO CURRENT DOMAIN",
1336  localIndexToGlobalIndexTable_[localIndex], particle->getId());
1337 
1338  //Flag the particle as not yet communicated
1339  particle->setMPIParticle(false);
1340  particle->setInMPIDomain(false);
1341  particle->setCommunicationComplexity(0);
1342  boundaryParticleListNeighbour_[localIndex][p] = nullptr;
1343  }
1344  //The particle has moved to a different domain
1345  else
1346  {
1347  logger(VERBOSE,
1348  "time: % | global index: % in list % | particle % | NEIGBOURING DOMAIN - TO OTHER DOMAIN",
1350  localIndexToGlobalIndexTable_[localIndex], particle->getId());
1351  //Cruelly destroy the particle without any mercy.
1352  ghostParticlesToBeDeleted.insert(particle);
1353  boundaryParticleListNeighbour_[localIndex][p] = nullptr;
1354  }
1355  }
1356  else
1357  {
1358  //check if the complexity has changed
1359  boundaryIndex = domainHandler_->getObject(
1360  localIndexToGlobalIndexTable_[localIndex])->findNearbyBoundaries(particle);
1361  complexityNew = boundaryIndex[0] + 3 * boundaryIndex[1] + 9 * boundaryIndex[2] + 13;
1362  if (particle->getCommunicationComplexity() != complexityNew)
1363  {
1364  logger(VERBOSE,
1365  "time: % | global index: % in list % | particle % | NEIGHBOURING DOMAIN - CHANGES COMPLEXITY",
1367  localIndexToGlobalIndexTable_[localIndex], particle->getId());
1368  //Cruelly destroy the particle without any mercy.
1369  ghostParticlesToBeDeleted.insert(particle);
1370  boundaryParticleListNeighbour_[localIndex][p] = nullptr;
1371  }
1372  }
1373  }
1374  }
1375  }
1376 
1383  void Domain::updateParticlePosition(int localIndex)
1384  {
1385  //process the updated information
1386  unsigned int index = 0;
1387  for (BaseParticle* particle : boundaryParticleListNeighbour_[localIndex])
1388  {
1389  logger.assert(particle->getId() == updatePositionDataReceive_[localIndex][index].id,
1390  "MPI particle lists are not in sync");
1391 
1392  //set position
1393  particle->setPreviousPosition(particle->getPosition());
1394  particle->setPosition(updatePositionDataReceive_[localIndex][index].position);
1395  particle->setOrientation(updatePositionDataReceive_[localIndex][index].orientation);
1396  if (std::is_base_of<MPILiquidFilmParticle,MPIParticle>())
1397  static_cast<LiquidFilmParticle*>(particle)->setLiquidVolume(updatePositionDataReceive_[localIndex][index].liquidVolume);
1398 
1399  //Update hGrid
1400  Vec3D displacement = particle->getPreviousPosition() - particle->getPosition();
1401  getHandler()->getDPMBase()->hGridUpdateMove(particle, displacement.getLengthSquared());
1402 
1403  index++;
1404  }
1405 }
1406 
1413 {
1414  //process the updated information
1415  unsigned int index = 0;
1416  for (BaseParticle* particle: boundaryParticleListNeighbour_[localIndex])
1417  {
1418  //set velocity
1419  particle->setVelocity(updateVelocityDataReceive_[localIndex][index].velocity);
1420  particle->setAngularVelocity(updateVelocityDataReceive_[localIndex][index].angularVelocity);
1421  index++;
1422  }
1423 }
1424 
1432 {
1433  //For all active boundaries
1434  for (int localIndex : boundaryList_)
1435  {
1436  numberOfParticlesSend_[localIndex] = boundaryParticleList_[localIndex].size();
1437  numberOfParticlesReceive_[localIndex] = boundaryParticleListNeighbour_[localIndex].size();
1438 
1439  //Increase capacity for the receiving data files
1440  updatePositionDataReceive_[localIndex].resize(numberOfParticlesReceive_[localIndex]);
1441  updateVelocityDataReceive_[localIndex].resize(numberOfParticlesReceive_[localIndex]);
1442 
1443 
1444  //Collect data
1445  for (BaseParticle* particle : boundaryParticleList_[localIndex])
1446  {
1447  updatePositionDataSend_[localIndex].push_back(copyPositionFrom(particle));
1448  updateVelocityDataSend_[localIndex].push_back(copyVelocityFrom(particle));
1449  }
1450 
1451  //Send and receive the data
1453  updatePositionDataReceive_[localIndex].data(), numberOfParticlesReceive_[localIndex],
1454  updatePositionDataSend_[localIndex].data(), numberOfParticlesSend_[localIndex],
1455  localIndex);
1457  updateVelocityDataReceive_[localIndex].data(), numberOfParticlesReceive_[localIndex],
1458  updateVelocityDataSend_[localIndex].data(), numberOfParticlesSend_[localIndex],
1459  localIndex);
1460  }
1461 }
1462 
1472 void Domain::finalisePositionAndVelocityUpdate(std::set<BaseParticle*>& ghostParticlesToBeDeleted)
1473 {
1474 
1475  //For all active boundaries
1476  for (int localIndex : boundaryList_)
1477  {
1478  updateParticlePosition(localIndex);
1479  updateParticleVelocity(localIndex);
1480  }
1481 
1482  //Based on the new position, update the particle lists.
1483  //Remove particles that left the communication zone, they will be re-communicated in a later step
1484  updateParticles(ghostParticlesToBeDeleted);
1485 
1486  //For all active boundaries clear the data lists
1487  for (int localIndex : boundaryList_)
1488  {
1489  updatePositionDataSend_[localIndex].clear();
1490  updateVelocityDataSend_[localIndex].clear();
1491  updatePositionDataReceive_[localIndex].clear();
1492  updateVelocityDataReceive_[localIndex].clear();
1493  }
1494 
1495 }
1496 
1501 {
1502  //For all active boundaries
1503  for (int localIndex : boundaryList_)
1504  {
1505  numberOfParticlesSend_[localIndex] = boundaryParticleList_[localIndex].size();
1506  numberOfParticlesReceive_[localIndex] = boundaryParticleListNeighbour_[localIndex].size();
1507 
1508  //Resize the vector to the correct size
1509  updateVelocityDataReceive_[localIndex].resize(numberOfParticlesReceive_[localIndex]);
1510 
1511  //Collect data
1512  for (BaseParticle* particle : boundaryParticleList_[localIndex])
1513  {
1514  updateVelocityDataSend_[localIndex].push_back(copyVelocityFrom(particle));
1515  }
1516 
1517  //Send the data
1519  updateVelocityDataReceive_[localIndex].data(), numberOfParticlesReceive_[localIndex],
1520  updateVelocityDataSend_[localIndex].data(), numberOfParticlesSend_[localIndex],
1521  localIndex);
1522  }
1523 }
1524 
1529 {
1530  //For all active boundaries
1531  for (int localIndex : boundaryList_)
1532  {
1533  updateParticleVelocity(localIndex);
1534  }
1535 
1536  //For all active boundaries clear the data lists
1537  for (int localIndex : boundaryList_)
1538  {
1539  updateVelocityDataSend_[localIndex].clear();
1540  updateVelocityDataReceive_[localIndex].clear();
1541  }
1542 }
1543 
1544 
1550 {
1551  return domainHandler_;
1552 }
1553 
1560 {
1561  //Step 1: For every MPIDomain boundary, create a list of particles that have to be transmitted
1562  //queue send and receive instructions for the number of particles
1565 
1566  //Step 2: queue send and receive of data
1569 
1570  //Step 3: Add the received particles to the particleHandler of the current domain
1572 }
1573 
1582 {
1583  //Step1: check if the particle has to be sent to other processors
1586 
1587  //Step2: queue send and receive data. Note for an inserted particle, no interactions should be required
1590 
1591  //Step3: Add the received particles to the particleHandler of the current domain
1594 }
1595 
1604 void Domain::updateStatus(std::set<BaseParticle*>& ghostParticlesToBeDeleted)
1605 {
1606  //Collect new positions and velocities and send them to the other domains
1609 
1610  //Receive the new positions and velocities from other domains
1611  //and update the mpi flagged particles accordingly. removes and switched particles in the lists
1612  finalisePositionAndVelocityUpdate(ghostParticlesToBeDeleted);
1614 }
1615 
1620 {
1621  //collect new velocity data and send
1624 
1625  //process the received data
1628 }
1629 
1636 {
1637  unsigned int count = 0;
1638  for (auto& index : boundaryParticleListNeighbour_)
1639  {
1640  count += index.size();
1641  }
1642  return count;
1643 }
1644 
1651 {
1652  unsigned int count = 0;
1653  for (auto& index : boundaryParticleListNeighbour_)
1654  {
1655  for (auto& p : index)
1656  {
1657  if (!p->isPeriodicGhostParticle())
1658  {
1659  count++;
1660  }
1661  }
1662  }
1663  return count;
1664 }
1665 
1669 void Domain::flushParticles(std::set<BaseParticle*>& toBeFlushedList)
1670 {
1671  //For all active boundaries
1672  for (int localIndex : boundaryList_)
1673  {
1674  flushParticlesFromList(boundaryParticleList_[localIndex], toBeFlushedList);
1675  flushParticlesFromList(boundaryParticleListNeighbour_[localIndex], toBeFlushedList);
1676  }
1677 }
1678 
1679 void Domain::flushParticlesFromList(std::vector<BaseParticle*>& list, std::set<BaseParticle*>& toBeFlushedList)
1680 {
1681  //Firstly: turn all particles that need to be flushed into nullptrs
1682  for (auto& p : list)
1683  {
1684  if (p != nullptr)
1685  {
1686  BaseParticle* particle1 = p;
1687  for (BaseParticle* particle2 : toBeFlushedList)
1688  {
1689  //If the particle was found in the list, make a nullptr
1690  if (particle1 == particle2)
1691  {
1692  logger(VERBOSE, "Removing particle from mpi domain at: %", particle1->getPosition());
1693  p = nullptr;
1694  }
1695  }
1696  }
1697  }
1698 }
1699 
1705 {
1706  return middle_;
1707 }
1708 
1714 {
1715  //For all active boundaries
1716  for (int i : boundaryList_)
1717  {
1720  }
1721 }
1722 
1728 void Domain::cleanCommunicationList(std::vector<BaseParticle*>& list)
1729 {
1730  for (int i = 0; i < list.size(); i++)
1731  {
1732  if (list[i] == nullptr)
1733  {
1734  list[i] = list.back();
1735  list.pop_back();
1736  i--;
1737  }
1738  }
1739 }
void prepareBoundaryDataTransmission()
Prepares the MPI transmission of particle and interaction data from particles in particleHandler.
Definition: Domain.cc:1117
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:562
std::vector< unsigned > numberOfParticlesReceive_
Counter that keeps track of the number of particles that are being received by this domain...
Definition: Domain.h:517
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:481
void updateParticleVelocity(int localIndex)
Updates the velocity of particles which are flagged as MPIParticles.
Definition: Domain.cc:1412
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
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 updateParticles(std::set< BaseParticle * > &ghostParticlesToBeDeleted)
This step updates all communication lists and particles in the communication zone.
Definition: Domain.cc:1279
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:1581
int getGlobalIndex()
Gets the global index of the domain.
Definition: Domain.cc:346
virtual void hGridGetInteractingParticleList(BaseParticle *obj, std::vector< BaseParticle * > &list)
Creates a list of neighbour particles obtained from the hgrid.
Definition: DPMBase.h:936
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:901
int getRank()
Gets the rank associated with the assigned processorID.
Definition: Domain.cc:328
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
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:791
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:522
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:549
Vec3D getMin() const
Definition: DPMBase.h:623
Vec3D middle_
Middle of the closed domain.
Definition: Domain.h:453
std::vector< double > domainMin_
Minimum domain bounds in the x,y and z direction.
Definition: Domain.h:443
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:470
void debugInformation()
Definition: Domain.cc:1067
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
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:1669
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
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:1559
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:524
Mdouble getSumOfInteractionRadii(const BaseParticle *particle) const
returns the sum of the radii plus the interactionDistance
Definition: BaseParticle.h:376
std::vector< int > localIndexToProcessorList_
look-up table to get the processor of the domain given a local domain index
Definition: Domain.h:476
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:806
void performBoundaryDataTransmission()
Collects data to be transmitted and then performs the transmission of the data.
Definition: Domain.cc:1192
std::vector< std::vector< BaseInteraction * > > newInteractionList_
Array that queues interactions that need to be transmitted.
Definition: Domain.h:507
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:1528
std::vector< std::vector< MPIParticleVelocity > > updateVelocityDataReceive_
Container that keeps a list of MPIParticleVelocities that are being received by this domain...
Definition: Domain.h:557
void updateStatus(std::set< BaseParticle * > &ghostParticlesToBeDeleted)
Updates particles that are not in the current domain and communicates newly added particles...
Definition: Domain.cc:1604
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:464
std::vector< int > boundaryList_
A list of indices of all the active boundaries.
Definition: Domain.h:487
unsigned int getNumberOfTrueMPIParticles()
Obtains the number of particles in the particleHandler that are MPIParticles, but NOT periodic partic...
Definition: Domain.cc:1650
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:567
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:1383
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:542
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:976
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:1713
void preparePositionAndVelocityUpdate()
Function that sends particle position and velocity data for ghost particles to other processors...
Definition: Domain.cc:1431
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:1472
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:768
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:497
ParticleHandler particleHandler
An object of the class ParticleHandler, contains the pointers to all the particles created...
Definition: DPMBase.h:1329
void cleanCommunicationList(std::vector< BaseParticle * > &list)
Removes nullptr's from a given particle list.
Definition: Domain.cc:1728
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:552
DomainHandler * getHandler() const
Gets the domainHandler.
Definition: Domain.cc:1549
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:359
void sendAndReceiveCount(MercuryMPITag tag, unsigned &countReceive, unsigned &countSend, unsigned localIndexNeighbour)
A symmetric communication between two domains exchanging a send/recieve count.
Definition: Domain.cc:1090
std::vector< unsigned > numberOfParticlesSend_
Counter that keeps track of the number of particles that are being send to other domains.
Definition: Domain.h:512
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:918
DomainHandler * domainHandler_
Pointer to the domain's DomainHandler container.
Definition: Domain.h:438
void collectBoundaryParticleData(int localIndex)
collects the data of a particle that has to be communicated to other processors
Definition: Domain.cc:886
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:1233
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:782
InteractionHandler interactionHandler
An object of the class InteractionHandler.
Definition: DPMBase.h:1359
Vec3D getMax() const
Definition: DPMBase.h:629
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:1339
std::vector< std::vector< BaseParticle * > > newBoundaryParticleList_
Array that queues particles that need to be transmitted.
Definition: Domain.h:502
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: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
Vec3D getMiddle() const
Gives the middle of the domain.
Definition: Domain.cc:1704
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:579
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:1619
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:547
int globalIndex_
Global index of the domain in the mesh.
Definition: Domain.h:459
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:572
std::vector< double > domainMax_
Maximum domain bounds in the x,y and z direction.
Definition: Domain.h:448
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:537
virtual void setMPIInteraction(void *interactionDataArray, unsigned int index, bool resetPointers)
virtual void hGridUpdateMove(BaseParticle *, Mdouble)
Definition: DPMBase.cc:1893
void prepareVelocityUpdate()
Function that sends particle velocity data for ghost particles.
Definition: Domain.cc:1500
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:1679
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:492
void processSentBoundaryParticles(unsigned index)
Bookkeep the newly send particles.
Definition: Domain.cc:954
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:527
Mdouble Z
Definition: Vector.h:65
Mdouble getTime() const
Returns the current simulation time.
Definition: DPMBase.cc:797
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:683
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:1635
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:532