MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
DomainHandler.cc
Go to the documentation of this file.
1 //Copyright (c) 2013-2020, The MercuryDPM Developers Team. All rights reserved.
2 //For the list of developers, see <http://www.MercuryDPM.org/Team>.
3 //
4 //Redistribution and use in source and binary forms, with or without
5 //modification, are permitted provided that the following conditions are met:
6 // * Redistributions of source code must retain the above copyright
7 // notice, this list of conditions and the following disclaimer.
8 // * Redistributions in binary form must reproduce the above copyright
9 // notice, this list of conditions and the following disclaimer in the
10 // documentation and/or other materials provided with the distribution.
11 // * Neither the name MercuryDPM nor the
12 // names of its contributors may be used to endorse or promote products
13 // derived from this software without specific prior written permission.
14 //
15 //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
16 //ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
17 //WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18 //DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
19 //DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
20 //(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
21 //LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
22 //ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23 //(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24 //SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 
26 
27 #include <iostream>
28 #include <string>
29 #include <cstdlib>
30 #include <Math/Helpers.h>
31 #include "DomainHandler.h"
33 #include "DPMBase.h"
34 #include <limits>
35 
41 {
44  logger(DEBUG, "DomainHandler::DomainHandler() finished");
45 }
46 
53  : BaseHandler<Domain>()
54 {
56  logger(DEBUG, "DomainHandler::DomainHandler(const DomainHandler &DH) finished");
57 }
58 
65 {
66  if (this != &rhs)
67  {
68  clear();
70  }
71  logger(DEBUG, "DomainHandler DomainHandler::operator =(const DomainHandler& rhs)");
72  return *this;
73 }
74 
79 {
80  logger(DEBUG, "DomainHandler::~DomainHandler() finished");
81 }
82 
93 void DomainHandler::createMesh(std::vector<Mdouble>& simulationMin, std::vector<Mdouble>& simulationMax,
94  std::vector<unsigned>& numberOfDomains, bool open)
95 {
96  //Clear the objects in the list, we're gonna make a new mesh, owyeah
97  objects_.clear();
98 
99  //Create the mesh
100  setNumberOfDomains(numberOfDomains);
101 
102  std::vector<unsigned> globalMeshIndex(3);
103  //Recursive function creating the domains
104  int dimCounter = 3;
105  //Compute the mesh size in each direction
106  std::vector<Mdouble> meshSize;
107  for (int d = 0; d < 3; d++)
108  {
109  meshSize.push_back((simulationMax[d] - simulationMin[d]) / numberOfDomains[d]);
110  }
111  //Create the domains
112  createDomains(simulationMin, simulationMax, globalMeshIndex, numberOfDomains, dimCounter, meshSize, open);
113 
114  //Create lookup tables
115  for (Domain* domain : objects_)
116  {
117  domain->setRank(domain->getGlobalIndex());
118  domain->createLookUpTable();
119  }
120 
121  //Disable boundaries that dont require communication
122  for (Domain* domain : objects_)
123  {
124  domain->disableBoundaries();
125  }
126 
127  //Output the result
128  std::string meshType;
129  if (open)
130  {
131  meshType = "open";
132  }
133  else
134  {
135  meshType = "closed";
136  }
137  if (PROCESSOR_ID == 0)
138  {
139  logger(INFO, "A simulation mesh has been created with % number of domains and % boundaries.", this->getSize(),
140  meshType);
141  }
142 }
143 
156 void DomainHandler::createDomains(std::vector<Mdouble> domainMin, std::vector<Mdouble> domainMax,
157  std::vector<unsigned>& globalMeshIndex, std::vector<unsigned>& numberOfDomains,
158  int dimCounter, std::vector<Mdouble>& meshSize, bool open)
159 {
160  //
161  if (dimCounter == 0)
162  {
163  //Create a new domain
164  Domain domain(globalMeshIndex);
165  domain.setHandler(this);
166  const std::vector<double>& domainBoundMin = domainMin;
167  const std::vector<double>& domainBoundMax = domainMax;
168  domain.setBounds(domainBoundMin, domainBoundMax, true);
169  this->copyAndAddObject(domain);
170  }
171  else
172  {
173  dimCounter--;
174  //Compute size of a domain in the dimCounter'th direction.
175  Mdouble boundLeft = domainMin[dimCounter];
176  //Starting with the bound left, create the number of domains in the given dimCounter-direction
177  for (int i = 0; i < numberOfDomains[dimCounter]; i++)
178  {
179  globalMeshIndex[dimCounter] = i;
180  domainMin[dimCounter] = boundLeft + i * meshSize[dimCounter];
181  domainMax[dimCounter] = boundLeft + (i + 1) * meshSize[dimCounter];
182  if ((i == 0) && (open))
183  {
184  domainMin[dimCounter] = -constants::inf;
185  }
186 
187  if ((i == numberOfDomains[dimCounter] - 1) && (open))
188  {
189  domainMax[dimCounter] = constants::inf;
190  }
191 
192  //Start recursively a new createDomains function
193  createDomains(domainMin, domainMax, globalMeshIndex, numberOfDomains, dimCounter, meshSize, open);
194  }
195  }
196 }
197 
204 {
205  //Puts the particle in the Particle list
207  //set the particleHandler pointer
208  D->setHandler(this);
209 }
210 
215 void DomainHandler::readAndAddObject(std::istream& is)
216 {
217 }
218 
223 void DomainHandler::readOldObject(std::istream& is)
224 {
225 }
226 
230 std::string DomainHandler::getName() const
231 {
232  return "DomainHandler";
233 }
234 
240 void DomainHandler::setCurrentDomainIndex(unsigned int index)
241 {
242  currentDomainIndex_ = index;
243 }
244 
251 {
253 }
254 
256 {
258 }
259 
266 {
267  return currentDomainIndex_;
268 }
269 
276 void DomainHandler::setNumberOfDomains(std::vector<unsigned>& numberOfDomains)
277 {
278  numberOfDomains_ = numberOfDomains;
279 }
280 
287 std::vector<unsigned> DomainHandler::getNumberOfDomains()
288 {
289  return numberOfDomains_;
290 }
291 
300 {
301  //Update the interaction distance
302  interactionDistance_ = interactionDistance;
303 
304  //Check if the domainsize is not too small
305  Domain* domain = getCurrentDomain();
306  logger.assert_always((domain->getDomainMax()[0] - domain->getDomainMin()[0]) > 2 * interactionDistance_,
307  "Size of the domain in x-direction is smaller than communication zone. Size: %, communication zone: %",
308  (domain->getDomainMax()[0] - domain->getDomainMin()[0]), 2 * interactionDistance_);
309  logger.assert_always((domain->getDomainMax()[1] - domain->getDomainMin()[1]) > 2 * interactionDistance_,
310  "Size of the domain in y-direction is smaller than communication zone. Size: %, communication zone: %",
311  (domain->getDomainMax()[1] - domain->getDomainMin()[1]), 2 * interactionDistance);
312  logger.assert_always((domain->getDomainMax()[2] - domain->getDomainMin()[2]) > 2 * interactionDistance_,
313  "Size of the domain in z-direction is smaller than communication zone. Size: %, communication zone: %",
314  (domain->getDomainMax()[2] - domain->getDomainMin()[2]), 2 * interactionDistance_);
315 }
316 
325 {
326  return interactionDistance_;
327 }
328 
331 {
332  //Step 1: obtain values i,j,k by looking at the position
333  //TODO this could possibly be stored in the domainHandler to save computational power
334  std::vector<int> decompositionVector(3);
335 
336  int i, j, k;
340 
341  //x-direction
342  if ((particle->getPosition().X - getDPMBase()->getXMin()) <= 0)
343  {
344  i = 0;
345  }
346  else
347  {
348  //Compute the relative domain cell it is in
349  i = floor((particle->getPosition().X - getDPMBase()->getXMin()) / dx);
350 
351  //Make sure we dont go over our number of domain bounds
352  if (i >= getDPMBase()->getNumberOfDomains()[0] - 1)
353  {
354  i = getDPMBase()->getNumberOfDomains()[0] - 1;
355  }
356  if (i < 0)
357  {
358  i = 0;
359  }
360  }
361 
362  //y-direction
363  if ((particle->getPosition().Y - getDPMBase()->getYMin()) <= 0)
364  {
365  j = 0;
366  }
367  else
368  {
369  j = floor((particle->getPosition().Y - getDPMBase()->getYMin()) / dy);
370 
371  //Make sure we dont go over our number of domain bounds
372  if (j >= getDPMBase()->getNumberOfDomains()[1] - 1)
373  {
374  j = getDPMBase()->getNumberOfDomains()[1] - 1;
375  }
376  if (i < 0)
377  {
378  j = 0;
379  }
380  }
381 
382  //z-direction
383  if ((particle->getPosition().Z - getDPMBase()->getZMin()) <= 0)
384  {
385  k = 0;
386  }
387  else
388  {
389  k = floor((particle->getPosition().Z - getDPMBase()->getZMin()) / dz);
390  //Make sure we dont go over our number of domain bounds
391  if (k >= getDPMBase()->getNumberOfDomains()[2] - 1)
392  {
393  k = getDPMBase()->getNumberOfDomains()[2] - 1;
394  }
395  if (k < 0)
396  {
397  k = 0;
398  }
399 
400  }
401 
402  //Step 2: obtain the processor number
403  int globalIndex = i +
404  getDPMBase()->getNumberOfDomains()[0] * j +
406  return globalIndex;
407 }
408 
410 {
411  return globalIndex;//ToProcessorList_[globalIndex];
412 }
413 
415 {
416  return getObject(globalIndex);
417 }
418 
419 void DomainHandler::updateStatus(std::set<BaseParticle*>& particlesToBeDeleted)
420 {
421  getCurrentDomain()->updateStatus(particlesToBeDeleted);
422 }
423 
425 {
427 }
428 
430 {
432 }
433 
435 {
436  //Create a single domain
437  Domain domain;
438  domain.setHandler(this);
439  std::vector<double> domainBoundMin = {-constants::inf, -constants::inf, -constants::inf};
440  std::vector<double> domainBoundMax = {constants::inf, constants::inf, constants::inf};
441  domain.setBounds(domainBoundMin, domainBoundMax, false);
442  this->copyAndAddObject(domain);
444 }
Domain * getParticleDomain(int globalIndex)
DomainHandler & operator=(const DomainHandler &rhs)
Assignment operator.
void setNumberOfDomains(std::vector< unsigned > &numberOfdomains)
Sets the number of domains in the domain handler.
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
Mdouble X
the vector components
Definition: Vector.h:65
void setBounds(std::vector< double > domainLeft, std::vector< double > domainRight, bool computeMiddle)
Sets the domain bounds.
Definition: Domain.cc:268
DomainHandler()
Default constructor, it creates an empty DomainHandler.
unsigned int getSize() const
Gets the size of the particleHandler (including mpi and periodic particles)
~DomainHandler() final
Destructor, it destructs the DomainHandler and all Domain it contains.
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
double Mdouble
Definition: GeneralDefine.h:34
Mdouble getZMax() const
If the length of the problem domain in z-direction is ZMax - ZMin, then getZMax() returns ZMax...
Definition: DPMBase.h:617
Mdouble getXMin() const
If the length of the problem domain in x-direction is XMax - XMin, then getXMin() returns XMin...
Definition: DPMBase.h:586
void readAndAddObject(std::istream &is) final
reads a domain object
void readOldObject(std::istream &is)
reads an old domain object
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
Mdouble getYMin() const
If the length of the problem domain in y-direction is YMax - YMin, then getYMin() returns YMin...
Definition: DPMBase.h:599
void addNewParticles()
Initialises the MPIParticles by communicating newly found particles.
Definition: Domain.cc:1559
unsigned int currentDomainIndex_
Index to the particular domain of this process.
unsigned int getCurrentDomainIndex() const
Gets the domain index assigned to the processor.
Mdouble getInteractionDistance()
Gets the interaction distance of the domain handler.
void createDomains(std::vector< Mdouble > domainMin, std::vector< Mdouble > domainMax, std::vector< unsigned > &globalMeshIndex, std::vector< unsigned > &numberOfDomains, int dimCounter, std::vector< Mdouble > &meshSize, bool open)
Recursive function that creates the domains for a 3D mesh.
std::string getName() const final
returns the name of the class
void updateStatus(std::set< BaseParticle * > &ghostParticlesToBeDeleted)
Updates particles that are not in the current domain and communicates newly added particles...
Definition: Domain.cc:1604
const Mdouble inf
Definition: GeneralDefine.h:44
void addNewParticles()
Mdouble getXMax() const
If the length of the problem domain in x-direction is XMax - XMin, then getXMax() returns XMax...
Definition: DPMBase.h:593
void setCurrentDomainIndex(unsigned int index)
This sets a domain to the processor.
std::vector< double > getDomainMin()
Gets the minimum domain bounds.
Definition: Domain.cc:310
Domain * getObject(const unsigned int id)
Gets a pointer to the Object at the specified index in the BaseHandler.
void addObject(Domain *D) final
Adds a Domain to the DomainHandler.
#define PROCESSOR_ID
Definition: GeneralDefine.h:63
std::enable_if<!std::is_pointer< U >::value, U * >::type copyAndAddObject(const U &object)
Creates a copy of a Object and adds it to the BaseHandler.
std::vector< Domain * > objects_
The actual list of Object pointers.
Definition: BaseHandler.h:302
Container to store the pointers to all objects that one creates in a simulation.
Definition: BaseHandler.h:50
std::vector< double > getDomainMax()
Gets the maximum domain bounds.
Definition: Domain.cc:319
Mdouble getYMax() const
If the length of the problem domain in y-direction is YMax - YMin, then getYMax() returns XMax...
Definition: DPMBase.h:605
Mdouble Y
Definition: Vector.h:65
int getParticleDomainGlobalIndex(BaseParticle *particle)
The simulation can be subdivided into Domain's used in parallel code.
Definition: Domain.h:63
Domain * getCurrentDomain()
Gets the domain assigned to the processor.
std::vector< unsigned > getNumberOfDomains()
returns the number of domains
Definition: DPMBase.cc:5025
int getParticleProcessor(int globalIndex)
Mdouble getZMin() const
If the length of the problem domain in z-direction is ZMax - ZMin, then getZMin() returns ZMin...
Definition: DPMBase.h:611
Container to store all Domain.
Definition: DomainHandler.h:46
virtual void addObject(T *object)
Adds a new Object to the BaseHandler.
Definition: BaseHandler.h:431
void updateVelocity()
Updates MPI particle velocity at the half-time step.
Definition: Domain.cc:1619
void setInteractionDistance(Mdouble interactionDistance)
Sets the interaction distance of the domain handler.
std::vector< unsigned > getNumberOfDomains()
Gets the number of domains in the domain handler.
std::vector< unsigned > numberOfDomains_
vector containing the number of domains in Cartesian direction
Mdouble interactionDistance_
Interaction distance between a domain boundary and the communication zone boundary.
void updateStatus(std::set< BaseParticle * > &particlesToBeDeleted)
void updateVelocity()
void copyContentsFromOtherHandler(const BaseHandler< Domain > &BH)
Function that copies the contents (vector of pointers, maxObject_, nextId_, DPMBase_) from one handle...
DPMBase * getDPMBase()
Gets the problem that is solved using this handler.
void createMesh(std::vector< Mdouble > &simulationMin, std::vector< Mdouble > &simulationMax, std::vector< unsigned > &numberOfDomains, bool open)
Creates a Cartesian square mesh in 3D.
Mdouble Z
Definition: Vector.h:65
void setHandler(DomainHandler *handler)
Sets the domainHandler.
Definition: Domain.cc:337
virtual void clear()
Empties the whole BaseHandler by removing all Objects and setting all other variables to 0...