MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
CGHandler.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.
27 #include "CG/CGHandler.h"
28 #include "DPMBase.h"
29 #include "Walls/InfiniteWall.h"
30 
35  : BaseHandler(ch)
36 {
37  setDPMBase(ch.getDPMBase());
38  //copyContentsFromOtherHandler(ch);
39 }
40 
48 {
49  if (this != &rhs)
50  {
51  clear();
52  setDPMBase(rhs.getDPMBase());
54  }
55  return *this;
56 #ifdef DEBUG_CONSTRUCTOR
57  std::cerr << "CGHandler::operator =(const CGHandler&) finished" << std::endl;
58 #endif
59 }
60 
65 {
66  //Puts the cg object in the list
68  //set the CGHandler pointer
69  cg->setHandler(this);
70 }
71 
75 std::string CGHandler::getName() const
76 {
77  return "CGHandler";
78 }
79 
80 void CGHandler::readAndAddObject(std::istream& is UNUSED)
81 {
82 
83 }
84 
85 void CGHandler::write(std::ostream& os UNUSED) const
86 {
87 
88 }
89 
91 {
92  //initialise all CG objects in the handler
93  for (BaseCG* it : *this)
94  it->initialise();
95 };
96 
98 {
99  //evaluate all CG objects in the handler
100  for (BaseCG* it : *this)
101  {
102  //if we are below timeMax and the next time step should be written
103  if (it->getTimeMin() <= getDPMBase()->getTime() && it->getTimeMax() > getDPMBase()->getTime()
104  && it->statFile.saveCurrentTimeStep(getDPMBase()->getNumberOfTimeSteps()))
105  {
106  //logger(INFO,"evaluate %, nt=%",it->statFile.getName(), getDPMBase()->getNumberOfTimeSteps());
107  it->statFile.setLastSavedTimeStep(getDPMBase()->getNumberOfTimeSteps());
108  it->evaluate();
109  }
110  }
111 };
112 
114 {
115  //enforce that data gets written
116  for (BaseCG* it : *this)
117  {
119  }
120  //evaluate all CG objects in the handler
121  //evaluate();
122  //finish all CG objects in the handler
123  for (BaseCG* it : *this)
124  it->finish();
125 }
126 
127 void CGHandler::restart(std::string name)
128 {
129  DPMBase* dpm = getDPMBase();
130 
131  // determines if the variable name contains the name of a restart file, or a problem name
132  // (in which case the restart file is assumed to be name.restart)
133  std::string::size_type endName = name.find(".restart");
134  if (endName == std::string::npos)
135  {
136  dpm->setName(name);
137  }
138  else
139  {
140  dpm->setName(name.substr(0, endName));
141  dpm->restartFile.setName(name);
142  }
143 
144  // read restart file
145  cgLogger(INFO, "Reading %", dpm->restartFile.getFullName());
146 
148  {
149  if (dpm->speciesHandler.getNumberOfObjects() == 0)
150  {
152  }
153  std::ifstream data(dpm->dataFile.getName());
154  if (data.is_open())
155  {
156  Mdouble N = 0, t = 0;
157  Vec3D min, max;
158  data >> N >> t >> min >> max;
159  data.close();
160  dpm->setDomain(min, max);
161  logger(INFO, "Reading domain size from %: min %, max % ", dpm->dataFile.getName(), min, max);
162  }
163  else
164  {
165  logger(ERROR, "Data file could not be opened");
166  }
167  //adding 10 walls by default
168 // while (dpm->wallHandler.getNumberOfObjects()<10) {
169 // dpm->wallHandler.copyAndAddObject(InfiniteWall(dpm->speciesHandler.getLastObject()));
170 // }
171  //what to do if restart file could not be loaded
172  cgLogger(WARN, "Using default dpm setup as % does not exist", dpm->restartFile.getName());
173  }
174  else
175  {
176  cgLogger(INFO, "Successfully restarted from %, t=%, Np=%, Nc=%", dpm->restartFile.getFullName(), dpm->getTime(),
179  }
180 }
181 
182 void CGHandler::restartAndEvaluateRestartFiles(const std::string& name)
183 {
184  restart(name);
186 }
187 
188 void CGHandler::restartAndEvaluateDataFiles(const std::string& name, bool evaluateFStatFiles)
189 {
190  restart(name);
191  evaluateDataFiles(evaluateFStatFiles);
192 }
193 
195 {
196  //recompute contact point (necessary for old restart files)
197  for (BaseInteraction* const c : getDPMBase()->interactionHandler)
198  {
199  if (c->getContactPoint().isNaN())
200  {
201  const Vec3D& p = c->getP()->getPosition();
202  const Vec3D& i = c->getI()->getPosition();
203  const BaseParticle* const PParticle = dynamic_cast<BaseParticle*>(c->getP());
204  const BaseParticle* const IParticle = dynamic_cast<BaseParticle*>(c->getI());
205  if (IParticle != nullptr)
206  {
207  const Vec3D branchVector = p - i;
208  const Mdouble distance = branchVector.getLength();
209  c->setNormal(branchVector / distance);
210  c->setOverlap(PParticle->getRadius() + IParticle->getRadius() - distance);
211  c->setDistance(distance);
212  c->setContactPoint(p - (PParticle->getRadius() - 0.5 * c->getOverlap()) * c->getNormal());
213  }
214  else
215  {
216  const InfiniteWall* const IWall = dynamic_cast<InfiniteWall*>(c->getI());
217  if (IWall != nullptr)
218  {
219  const Mdouble dist = Vec3D::dot(i - p, IWall->getNormal());
220  const Mdouble r = dynamic_cast<BaseParticle*>(c->getP())->getRadius();
221  const Vec3D branch = (0.5 * (dist + r)) * IWall->getNormal();
222  c->setNormal(-IWall->getNormal());
223  c->setContactPoint(p + branch);
224  }
225  else
226  {
227  const Vec3D IP = p - i;
228  const Mdouble distance = Vec3D::getLength(IP);
229  c->setNormal(IP / distance);
230  c->setContactPoint(i);
231  }
232  }
233  //logger(INFO,"c%",c->getContactPoint());
234  static bool firstTime = true;
235  if (firstTime)
236  {
237  cgLogger(WARN,
238  "recomputing contact point, as contact point information is not available\n"
239  "Contact point is placed assuming spherical particles\n"
240  "Complex walls, non-spherical particles, periodic walls can create errors");
241  firstTime = false;
242  }
243  }
244  }
245 }
246 
248 {
249 
250 #ifdef MERCURY_USE_MPI
251  //Make sure that the number of processors is equal to the number of processors used for the run
252  MPIContainer& communicator = MPIContainer::Instance();
253  std::vector<unsigned> numberOfDomains = this->getDPMBase()->getNumberOfDomains();
254  int numberOfRequiredProcessors = numberOfDomains[0]*numberOfDomains[1]*numberOfDomains[2];
255  if(!(numberOfRequiredProcessors == communicator.getNumberOfProcessors()))
256  {
257  if (communicator.getProcessorID() == 0)
258  {
259  logger(ERROR,"Please re-run the program with % cores",numberOfRequiredProcessors);
260  }
261  else
262  {
263  std::exit(-1);
264  }
265  }
266 #endif
267  // reset counters so reading begins with the first data/fstat file
268  DPMBase* const dpm = getDPMBase();
269 
270  //define and check time limits
271  Mdouble timeMin = getTimeMin();
272  Mdouble timeMax = getTimeMax();
273  if (getDPMBase()->getTime() > timeMax)
274  {
275  logger(ERROR, "initial restart file (t=%) is beyond the maximum cg time (tMax=%)", dpm->getTime(), timeMax);
276  }
277  else
278  while (getDPMBase()->getTime() < timeMin)
279  {
280  cgLogger(INFO, "Skipped %, t = %, because time is below tMin = %", dpm->restartFile.getFullName(),
281  dpm->getTime(), timeMin);
282  //the particle and wall handler is cleared here, because BaseSpecies doesn't delete particles belonging to it
283  dpm->particleHandler.clear();
284  dpm->wallHandler.clear();
285  dpm->readRestartFile();
286  }
287 
288  //call initialise to set up mesh, stat files, etc
289  initialise();
290 
291  // evaluate restart files, starting with the one already read (since interactions where not read)
292  while (dpm->readRestartFile() && getDPMBase()->getTime() < timeMax)
293  {
294  cgLogger(INFO, "Read %, t=%, Np=%, Nc=%", dpm->restartFile.getFullName(), dpm->getTime(),
296 
297  //recompute contact point (necessary for old restart files)
299 
300  evaluate();
301 
302  //the particle and wall handler is cleared here, because BaseSpecies doesn't delete particles belonging to it
303  dpm->particleHandler.clear();
304  dpm->wallHandler.clear();
305 
306  //continue if the next restart file can be read and the max time has not been reached
307  }
308  cgLogger(INFO, "Finished reading from %", dpm->dataFile.getFullName());
309 
310  finish();
311  dpm->dataFile.close();
312  return true;
313 }
314 
315 bool CGHandler::evaluateDataFiles(bool evaluateFStatFiles)
316 {
317 
318 #ifdef MERCURY_USE_MPI
319  //Make sure that the number of processors is equal to the number of processors used for the run
320  MPIContainer& communicator = MPIContainer::Instance();
321  std::vector<unsigned> numberOfDomains = this->getDPMBase()->getNumberOfDomains();
322  int numberOfRequiredProcessors = numberOfDomains[0]*numberOfDomains[1]*numberOfDomains[2];
323  if(!(numberOfRequiredProcessors == communicator.getNumberOfProcessors()))
324  {
325  if (communicator.getProcessorID() == 0)
326  {
327  logger(ERROR,"Please re-run the program with % cores",numberOfRequiredProcessors);
328  }
329  else
330  {
331  std::exit(-1);
332  }
333  }
334 #endif
335  // reset counters so reading begins with the first data/fstat file
336  DPMBase* const dpm = getDPMBase();
337  //dpm->interactionHandler.clear();
340 
341  initialise();
342 
343  // return false if no data file can be read
345  if (!dpm->readNextDataFile()) return false;
346 
347  //define and check time limits
348  const Mdouble timeMin = getTimeMin();
349  const Mdouble timeMax = getTimeMax();
350  if (dpm->getTime() > timeMax)
351  {
352  logger(ERROR, "Time stamp of initial data file (%) is beyond the maximum cg time (tMax=%)", dpm->getTime(),
353  timeMax);
354  }
355  else
356  while (dpm->getTime() < timeMin)
357  {
358  if (evaluateFStatFiles) dpm->readNextFStatFile();
359  cgLogger(INFO, "Skipped %, t = %, because time is below tMin = %", dpm->dataFile.getFullName(),
360  dpm->getTime(), timeMin);
361  dpm->readNextDataFile();
362  }
363 
364  // read data file
365  do
366  {
367  if (evaluateFStatFiles) dpm->readNextFStatFile();
368 
369  cgLogger(INFO, "Read %, t=%, Np=%, Nc=%", dpm->dataFile.getFullName(), dpm->getTime(),
371 
372  evaluate();
373  } while (dpm->readNextDataFile() && getDPMBase()->getTime() <= timeMax);
374  cgLogger(INFO, "Finished reading from %", dpm->dataFile.getFullName());
375 
376  finish();
377  dpm->dataFile.close();
378  return true;
379 }
380 
382 {
383  Mdouble time = constants::inf;
384  for (BaseCG* it : *this)
385  {
386  time = std::min(time, it->getTimeMin());
387  }
388  return time;
389 }
390 
392 {
393  Mdouble time = -constants::inf;
394  for (BaseCG* it : *this)
395  {
396  time = std::max(time, it->getTimeMax());
397  }
398  return time;
399 }
This class contains all information and functions required for communication between processors...
Definition: MpiContainer.h:125
Container that stores all CG objects.
Definition: CGHandler.h:64
void finish()
Contains the code executed after the last time step.
Definition: CGHandler.cc:113
virtual void evaluate()=0
Called after a given number of time steps (statFile::saveCount_) to evaluate the CG fields...
Mdouble getTimeMax() const
Returns timeMax_, the upper limit of the temporal domain.
Definition: BaseCG.cc:180
static MPIContainer & Instance()
fetch the instance to be used for communication
Definition: MpiContainer.h:130
The DPMBase header includes quite a few header files, defining all the handlers, which are essential...
Definition: DPMBase.h:72
bool saveCurrentTimeStep(unsigned int ntimeSteps)
determined if this time step has to be written; if so, opens the output file
Definition: File.cc:313
unsigned int getSize() const
Gets the size of the particleHandler (including mpi and periodic particles)
Definition: BaseHandler.h:655
std::string getName() const final
Definition: CGHandler.cc:75
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
double Mdouble
Definition: GeneralDefine.h:34
void readAndAddObject(std::istream &is) final
Reads objects into the CGHandler from an istream (currently not implemented).
Definition: CGHandler.cc:80
void setDPMBase(DPMBase *DPMBase)
Sets the problem that is solved using this handler.
void setLastSavedTimeStep(unsigned int lastSavedTimeStep)
Sets File::nextSavedTimeStep_.
Definition: File.cc:303
Mdouble getTimeMin()
Definition: CGHandler.cc:381
void readNextFStatFile()
Reads the next fstat file.
Definition: DPMBase.cc:2754
void setCounter(unsigned int counter)
Allows the user to set the file counter according to his need. Sets File::counter_.
Definition: File.cc:232
CGHandler()=default
Default constructor, creates an empty CGHandler.
const std::string getFullName() const
Also allows to access the file name, however with additional information which is the file counter...
Definition: File.cc:171
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
Vec3D getNormal() const
Access function for normal.
void write(std::ostream &os) const
Writes objects into the CGHandler to an ostream (currently not implemented).
Definition: CGHandler.cc:85
unsigned initialFileCounter
Definition: CGHandler.h:149
void evaluate()
Contains the code executed at each time step.
Definition: CGHandler.cc:97
std::size_t getNumberOfProcessors() const
Get the total number of processors participating in this simulation.
static Mdouble dot(const Vec3D &a, const Vec3D &b)
Calculates the dot product of two Vec3D: .
Definition: Vector.cc:76
CGHandler & operator=(const CGHandler &rhs)
Assignment operator that copies the pointer to the DPMBase and all objects.
Definition: CGHandler.cc:47
void clear() override
Empties the whole ParticleHandler by removing all BaseParticle.
bool evaluateDataFiles(bool evaluateFStatFiles=true)
does the same as StatisticsVector::statistics_from_fstat_and_data: loads a restart file (if existing)...
Definition: CGHandler.cc:315
void restart(std::string name)
loads restart file, before evaluateDataFiles is run
Definition: CGHandler.cc:127
const Mdouble inf
Definition: GeneralDefine.h:44
Stores information about interactions between two interactable objects; often particles but could be ...
static Mdouble getLength(const Vec3D &a)
Calculates the length of a Vec3D: .
Definition: Vector.cc:331
void initialise()
Contains the code executed before the first time step.
Definition: CGHandler.cc:90
const unsigned NEVER
Definition: File.h:35
void close()
Closes the file by calling fstream_.close()
Definition: File.cc:408
void setDomain(const Vec3D &min, const Vec3D &max)
Sets the minimum coordinates of the problem domain.
Definition: DPMBase.cc:1059
std::size_t getProcessorID()
Reduces a scalar on all processors to one scalar on a target processor.
File dataFile
An instance of class File to handle in- and output into a .data file.
Definition: DPMBase.h:1370
File fStatFile
An instance of class File to handle in- and output into a .fstat file.
Definition: DPMBase.h:1375
bool evaluateRestartFiles()
Definition: CGHandler.cc:247
void computeContactPoints()
Definition: CGHandler.cc:194
ParticleHandler particleHandler
An object of the class ParticleHandler, contains the pointers to all the particles created...
Definition: DPMBase.h:1329
#define UNUSED
Definition: GeneralDefine.h:39
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.
Definition: BaseHandler.h:379
Mdouble getRadius() const
Returns the particle's radius.
Definition: BaseParticle.h:345
Container to store the pointers to all objects that one creates in a simulation.
Definition: BaseHandler.h:50
bool readNextDataFile(unsigned int format=0)
Reads the next data file with default format=0. However, one can modify the format based on whether t...
Definition: DPMBase.cc:2575
Species< LinearViscoelasticNormalSpecies > LinearViscoelasticSpecies
virtual unsigned int getNumberOfObjects() const
Gets the number of real Object in this BaseHandler. (i.e. no mpi or periodic particles) ...
Definition: BaseHandler.h:648
SpeciesHandler speciesHandler
A handler to that stores the species type i.e. LinearViscoelasticSpecies, etc.
Definition: DPMBase.h:1319
InteractionHandler interactionHandler
An object of the class InteractionHandler.
Definition: DPMBase.h:1359
WallHandler wallHandler
An object of the class WallHandler. Contains pointers to all the walls created.
Definition: DPMBase.h:1339
void restartAndEvaluateDataFiles(const std::string &name, bool evaluateFStatFiles=true)
Definition: CGHandler.cc:188
void setName(const std::string &name)
Allows to set the name of all the files (ene, data, fstat, restart, stat)
Definition: DPMBase.cc:412
std::vector< unsigned > getNumberOfDomains()
returns the number of domains
Definition: DPMBase.cc:5025
Mdouble getLength() const
Calculates the length of this Vec3D: .
Definition: Vector.cc:320
File statFile
File class to handle the output into a .stat file.
Definition: BaseCG.h:369
virtual void addObject(T *object)
Adds a new Object to the BaseHandler.
Definition: BaseHandler.h:431
This is a class defining walls.
Definition: InfiniteWall.h:47
File restartFile
An instance of class File to handle in- and output into a .restart file.
Definition: DPMBase.h:1385
void addObject(BaseCG *cg) final
Definition: CGHandler.cc:64
void decreaseCounter()
Definition: File.h:150
virtual void initialise()=0
Called at the beginning of the DPM simulation to initialise the cg evaluation and to open the statFil...
Mdouble getTimeMax()
Definition: CGHandler.cc:391
Logger< CG_LOGLEVEL > cgLogger("MercuryCG")
void copyContentsFromOtherHandler(const BaseHandler< BaseCG > &BH)
Function that copies the contents (vector of pointers, maxObject_, nextId_, DPMBase_) from one handle...
Definition: Vector.h:49
void setName(const std::string &name)
Sets the file name, e.g. "Name.data".
Definition: File.cc:199
DPMBase * getDPMBase()
Gets the problem that is solved using this handler.
Definition: BaseHandler.h:725
bool readRestartFile(ReadOptions opt=ReadOptions::ReadAll)
Reads all the particle data corresponding to a given, existing . restart file (for more details regar...
Definition: DPMBase.cc:2896
Mdouble getTime() const
Returns the current simulation time.
Definition: DPMBase.cc:797
virtual void finish()=0
Called at the end of the DPM simulation to finish the cg evaluation and to close the statFile...
Base class of all CG objects, needed to store the various CG objects in the CGHandler.
Definition: BaseCG.h:56
const std::string & getName() const
Allows to access the file name, e.g., "problem.data".
Definition: File.cc:166
Mdouble getTimeMin() const
Returns timeMin_, the lower limit of the temporal domain.
Definition: BaseCG.cc:175
void setHandler(CGHandler *handler)
Sets handler_, the pointer to the CGHandler.
Definition: BaseCG.cc:72
virtual void clear()
Empties the whole BaseHandler by removing all Objects and setting all other variables to 0...
void restartAndEvaluateRestartFiles(const std::string &name)
Definition: CGHandler.cc:182