MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
DPMBase.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 #include "DPMBase.h"
26 #include <iostream>
27 #include <iomanip>
28 #include <algorithm>
29 #include <fstream>
30 #include <cstdlib>
31 #include <limits>
32 #include <string>
33 #include <cstdio>
35 #include <cstring>
40 #include "CMakeDefinitions.h"
41 #include "DPMBaseXBalls.icc" //This is part of this class and just separates out the stuff to do with xballs.
42 #include "Logger.h"
44 #include "Walls/BaseWall.h"
45 #include "Walls/InfiniteWall.h"
49 
50 //This is MPI related
51 #include "MpiContainer.h"
52 #include "MpiDataClass.h"
53 #include "Domain.h"
54 #ifdef MERCURY_USE_MPI
55 #include <mpi.h>
56 #endif
57 
58 //This is OMP related
59 #ifdef MERCURY_USE_OMP
60 #include <omp.h>
61 #endif
62 
73 [[noreturn]] void logWriteAndDie(const std::string& module, std::string message)
74 {
75  std::cerr << "A fatal error has occured"
76  << "\n Module :" << module
77  << "\n Message :" << message << std::endl;
78 
79  std::exit(-1);
80 }
81 
90 std::ostream& operator<<(std::ostream& os, DPMBase& md)
91 {
92  md.write(os);
93  return os;
94 }
95 
115 DPMBase::DPMBase(const DPMBase& other) : wallVTKWriter_(other.wallVTKWriter_),
116  interactionVTKWriter_(other.interactionVTKWriter_),
117  boundaryVTKWriter_(other.boundaryVTKWriter_)
118 {
119  setName(other.getName());
120  runNumber_ = other.runNumber_;
123  gravity_ = other.gravity_;
124 /* xMin_ = other.xMin_;
125  xMax_ = other.xMax_;
126  yMin_ = other.yMin_;
127  yMax_ = other.yMax_;
128  zMin_ = other.zMin_;
129  zMax_ = other.zMax_;*/
130  min_ = other.min_;
131  max_ = other.max_;
133  time_ = other.time_;
134  timeStep_ = other.timeStep_;
136  timeMax_ = other.timeMax_;
137  restartVersion_ = other.restartVersion_; //to read new and old restart data
138  restarted_ = other.restarted_; //to see if it was restarted or not
139  append_ = other.append_;
140  rotation_ = other.rotation_;
141  xBallsColourMode_ = other.xBallsColourMode_; // sets the xballs argument cmode (see xballs.txt)
142  xBallsVectorScale_ = other.xBallsVectorScale_; // sets the xballs argument vscale (see xballs.txt)
143  xBallsScale_ = other.xBallsScale_; // sets the xballs argument scale (see xballs.txt)
144  xBallsAdditionalArguments_ = other.xBallsAdditionalArguments_; // std::string where additional xballs argument can be specified (see xballs.txt)
148 
149 //effectively saying "if there exists a CONTACT_LIST_HGRID, copy it, if not, ignore.
150 #ifdef CONTACT_LIST_HGRID
151  possibleContactList=other.possibleContactList;
152 #endif
153  random = other.random;
154 
159  wallHandler.setDPMBase(this);
162  //Initialise the handlers
165 
166  //setting contents equal to the other handlers!
169  cgHandler = other.cgHandler;
170  //cgHandler = other.cgHandler.copy(); //todo
171  //cgHandler.setDPMBase(this);
172  wallHandler = other.wallHandler;
175  vtkWriter_ = other.vtkWriter_;
180 }
181 
187 DPMBase::DPMBase() : wallVTKWriter_(wallHandler), interactionVTKWriter_(interactionHandler), boundaryVTKWriter_(boundaryHandler)
188 {
189  constructor();
190 }
191 
201 {
202  //constructor();
203  dataFile.getFstream().precision(10);
204  fStatFile.getFstream().precision(10);
205  eneFile.getFstream().precision(10);
206  restartFile.getFstream().precision(
207  std::numeric_limits<double>::digits10); //highly accurate, so the overlap is accurate
208  statFile.getFstream().precision(10);
209  statFile.getFstream().setf(std::ios::left);
210  interactionFile.getFstream().precision(10);
211  name_ = ""; // needs to be user-specified, otherwise checkSettings throws error
212  //by default, the fileType of all files is ONE_FILE. However, by default we don't want an interaction file since it
213  // is very large.
215 
216  runNumber_ = 0;
217 
218  //Decomposition direction for MPI
219  numberOfDomains_ = {1, 1, 1};
220 
221  //Check if MPI is already initialised
222  initialiseMPI();
223 
224  //This sets the maximum number of particles
229  cgHandler.setDPMBase(this);
231  wallHandler.setDPMBase(this);
237 
238  //set defaults for DPMBase parameters
241  setRestarted(false);
242  setGravity(Vec3D(0, 0, 0));
243 
244  //This is the parameter of the numerical part
245  setTime(0);
246  numberOfTimeSteps_ = 0;
247  setTimeMax(0);
248  timeStep_ = 0; // needs to be user-specified, otherwise checkSettings throws error
249  setSaveCount(20);
250 
251  //This sets the default xballs domain
252  min_ = Vec3D(0, 0, 0);
253  max_ = Vec3D(0, 0, 0); // needs to be user-specified, otherwise checkSettings throws error
254 
255  //sets the default write particles data in VTK format flag to false
256  writeParticlesVTK_ = false;
259  vtkWriter_ = nullptr;
260 
261  //defines logger behaviour
263 
264  setName(""); // needs to be user-specified, otherwise checkSettings throws error
265 
266  //Default mode is energy with no scale of the vectors
267  xBallsColourMode_ = 0;
268  xBallsVectorScale_ = -1;
269  xBallsScale_ = -1;
271  setAppend(false);
272 
273  //The default random seed is 0
275 
276  logger(DEBUG, "DPMBase problem constructor finished");
277 
278  readSpeciesFromDataFile_ = false;
279 
281 }
282 
288 {
289  delete vtkWriter_;
290 }
291 
296 {
297  return dataFile;
298 }
299 
304 {
305  return eneFile;
306 }
307 
312 {
313  return fStatFile;
314 }
315 
320 {
321  return restartFile;
322 }
323 
328 {
329  return statFile;
330 }
331 
336 {
337  return interactionFile;
338 }
339 
340 
345 {
346  return dataFile;
347 }
348 
352 const File& DPMBase::getEneFile() const
353 {
354  return eneFile;
355 }
356 
361 {
362  return fStatFile;
363 }
364 
369 {
370  return restartFile;
371 }
372 
377 {
378  return statFile;
379 }
383 const File& DPMBase::getInteractionFile() const
385 {
386  return interactionFile;
387 }
388 
389 const std::string& DPMBase::getName() const
390 {
391  return name_;
392 }
393 
398 void DPMBase::setSaveCount(unsigned int saveCount)
399 {
400  dataFile.setSaveCount(saveCount);
401  fStatFile.setSaveCount(saveCount);
402  restartFile.setSaveCount(saveCount);
403  statFile.setSaveCount(saveCount);
404  eneFile.setSaveCount(saveCount);
405  for (auto cg : cgHandler)
406  cg->statFile.setSaveCount(saveCount);
407 }
408 
412 void DPMBase::setName(const std::string& name)
413 {
414  if (NUMBER_OF_PROCESSORS > 1)
415  {
416  name_ = name; // was before this->name_ = name
423  }
424  else
425  {
426  name_ = name; // was before this->name_ = name
427  dataFile.setName(name_ + ".data");
428  fStatFile.setName(name_ + ".fstat");
429  restartFile.setName(name_ + ".restart");
430  statFile.setName(name_ + ".stat");
431  eneFile.setName(name_ + ".ene");
432  interactionFile.setName(name_ + ".interaction");
433  }
434 }
435 
439 void DPMBase::setName(const char* name)
440 {
441  setName(std::string(name));
442 }
443 
450 {
451  dataFile.setFileType(fileType);
452  fStatFile.setFileType(fileType);
453  restartFile.setFileType(fileType);
454  statFile.setFileType(fileType);
455  eneFile.setFileType(fileType);
456 }
457 
462 {
463  dataFile.setCounter(0);
466  statFile.setCounter(0);
467  eneFile.setCounter(0);
474 }
475 
479 void DPMBase::setOpenMode(std::fstream::openmode openMode)
480 {
481  dataFile.setOpenMode(openMode);
482  fStatFile.setOpenMode(openMode);
483  restartFile.setOpenMode(openMode);
484  statFile.setOpenMode(openMode);
485  eneFile.setOpenMode(openMode);
486  interactionFile.setOpenMode(openMode);
487 }
488 
493 {
494  dataFile.close();
495  fStatFile.close();
496  restartFile.close();
497  statFile.close();
498  eneFile.close();
500 }
501 
508 void DPMBase::setLastSavedTimeStep(unsigned int nextSavedTimeStep)
509 {
510  dataFile.setLastSavedTimeStep(nextSavedTimeStep);
511  fStatFile.setLastSavedTimeStep(nextSavedTimeStep);
512  restartFile.setLastSavedTimeStep(nextSavedTimeStep);
513  //statFile.setLastSavedTimeStep(nextSavedTimeStep); //this one is not force-written
514  eneFile.setLastSavedTimeStep(nextSavedTimeStep);
515 }
516 
529 {
531 
532  if (!getRestarted())
533  {
535  }
536 }
537 
543 {
544  int counter;
545 
546  FILE* counter_file;
547  //checking if there exists already a file named "COUNTER_DONOTDEL" which can be opened for
548  //input and output (returns "true" if no such file exists).
549  if ((counter_file = fopen("COUNTER_DONOTDEL", "r+")) == nullptr)
550  {
551  //if a file does not already exist, checks whether a new file can be created
552  //retutns "true" if a new file CANNOT be created
553  if ((counter_file = fopen("COUNTER_DONOTDEL", "w")) == nullptr)
554  {
555  //If true, outputs an error message and ends the program
556  fprintf(stderr, "\n\n\tERROR :: Counter File NOT found, please re-create\n\n");
557  fclose(counter_file);
558  exit(-1);
559  }
560  //alternatively, if a new file CAN be created...
561  else
562  {
563  //starts the new counter file, writing to it the value "1"
564  fprintf(counter_file, "1");
565  fprintf(stderr, "Counter File created\n");
566  fclose(counter_file);
567  return 1;
568  }
569  }
570  //alternatively, if a counter file DOES already exist...
571  else
572  {
573  //...checks if there exists only 1 value in the file (as would be expected from a COUNTER_DONOTDEL file...
574  if (fscanf(counter_file, "%d", &counter) != 1)
575  {
576  //...and if not, returns an error.
577  fprintf(stderr, "\n\n\tERROR :: Counter File found, but something went wrong with reading it\n\n");
578  fclose(counter_file);
579  exit(-1);
580  }
581  else
582  {
583  //...otherwise, if all has been successful, returns the current value of the file!
584  fclose(counter_file);
585  return counter;
586  }
587  }
588 
589 }
590 
595 void DPMBase::setRunNumber(int runNumber)
596 {
597  runNumber_ = runNumber;
598 }
599 
607 {
608  return runNumber_;
609 }
610 
618 {
619  //opening two filestreams - counter_file and counter_file2
620  std::fstream counter_file, counter_file2;
621  //declares an integer, temp_counter
622  int temp_counter;
623  //attempts to open the COUNTER_DONOTDEL text file
624  counter_file.open("COUNTER_DONOTDEL", std::ios::in);
625  //gives error message if file could not be successfully opened and ends the program
626  if (counter_file.fail())
627  {
628  fprintf(stderr, "\n\n\tERROR :: Counter File NOT found, please re-create\n\n");
629  counter_file.close();
630  exit(0);
631  }
632  // if opened successfully, reads in the counter corresponding to the current run number
633  //and stored it in the "temp_counter" variable
634  counter_file >> temp_counter;
635  counter_file.close();
636  //Increments the temp_counter
637  temp_counter++;
638  //opens an output stream to the COUNTER_DONOTDEL file
639  counter_file2.open("COUNTER_DONOTDEL", std::ios::out);
640  if (counter_file2.fail())
641  {
642  fprintf(stderr, "\n\n\tERROR :: Counter File NOT found, please re-create2\n\n");
643  counter_file2.close();
644  exit(0);
645  }
646  //writes the new valuer of the counter to COUNTER_DONOTDEL
647  counter_file2 << temp_counter;
648 
649  counter_file2.close();
650 }
651 
659 std::vector<int> DPMBase::get1DParametersFromRunNumber(int sizeX) const
660 {
661  // Declare a vector of integers capable of storing 2 values
662  std::vector<int> temp(2);
663 
664  // Declare and initialise for the current simulation run number
665  int counter = getRunNumber();
666 
667  // Give studyNum value 0 if study is incomplete, otherwise value > 0
668  int studyNum = (counter-1)/sizeX;
669  counter = counter - sizeX*studyNum;
670 
671  int i = ((counter - 1) % sizeX) + 1;
672  logger(INFO,"StudyNum: % \t Counter: % \t i: %", studyNum, counter, i);
673  temp[0] = studyNum;
674  temp[1] = i;
675 
676  return temp;
677 }
678 
687 std::vector<int> DPMBase::get2DParametersFromRunNumber(int sizeX, int sizeY) const
688 {
689  //declares a vector of integers capable of storing 3 values,
690  std::vector<int> temp(3);
691  //declares and initialises an integer variable named "counter"
692  //with the current counter number, runNumber_
693  int counter = getRunNumber();
694  //calculates the total size of the study, i.e. the number of points
695  //in the 2D parameter space explored
696  int studySize = sizeX * sizeY;
697  //(counter - 1) / studySize gives a fraction comparing the number of runs conducted so far
698  //to the total size of the study, i.e. the total number of runs that need to be performed.
699  //since studyNum is an integer, will declare zero until an adequate number of runs has been performed,
700  //at which point it will equal 1
701  int studyNum = (counter - 1) / studySize;
702 
703  counter = counter - studySize * studyNum;
704  int i = ((counter - 1) % sizeX) + 1;
705  int j = ((counter - i) / sizeX) + 1;
706  logger(INFO,"StudyNum: % \t Counter: % \t i: % \t j: %", studyNum, counter, i, j);
707 
708  temp[0] = studyNum;
709  temp[1] = i;
710  temp[2] = j;
711 
712  return (temp);
713 }
714 
724 std::vector<int> DPMBase::get3DParametersFromRunNumber(int sizeX, int sizeY, int sizeZ) const
725 {
726  //declares a vector of integers capable of storing 4 values,
727  std::vector<int> temp(4);
728  //declares and initialises an integer variable named "counter"
729  //with the current counter number, runNumber_
730  int counter = getRunNumber();
731  //calculates the total size of the study, i.e. the number of points
732  //in the 3D parameter space explored
733  int studySize = sizeX * sizeY * sizeZ;
734  //(counter - 1) / studySize gives a fraction comparing the number of runs conducted so far
735  //to the total size of the study, i.e. the total number of runs that need to be performed.
736  //since studyNum is an integer, will declare zero until an adequate number of runs has been performed,
737  //at which point it will equal 1
738  int studyNum = (counter - 1) / studySize;
739 
740  counter = counter - studySize * studyNum;
741  int i = ((counter-1) % sizeX) + 1;
742  int j = static_cast<int>(std::floor((counter-1)/sizeX)) % sizeY + 1;
743  int k = static_cast<int>(std::floor((counter-1)/(sizeX*sizeY))) % sizeZ + 1;
744  logger(INFO,"StudyNum: % \t Counter: % \t i: % \t j: % \t k: %", studyNum, counter, i, j, k);
745 
746  temp[0] = studyNum;
747  temp[1] = i;
748  temp[2] = j;
749  temp[3] = k;
750 
751  return (temp);
752 }
753 
764 int DPMBase::launchNewRun(const char* name, bool quick UNUSED)
765 {
766  //defines an (empty) stringstream named "com"
767  std::stringstream com("");
768  //adds the name of the code to run (fed in as an argument)
769  //to the "com" string and appends the string with " &"
770  com << name << " &";
771  //converts the stringstream "com" to a standard string, and then
772  //converts this string to a C string
773  //the string is then fed to the "system" function, which will run the named command
774  return system(com.str().c_str());
775 }
776 
788 void DPMBase::solve(int argc, char* argv[])
789 {
790  readArguments(argc, argv);
791  solve();
792 }
793 
798 {
799  return time_;
800 }
801 
806 {
807  return time_ + timeStep_;
808 }
809 
813 unsigned int DPMBase::getNumberOfTimeSteps() const
814 {
815  return numberOfTimeSteps_;
816 }
817 
826 {
827  Mdouble diff = time_ - time;
828  time_ = time;
829  //this sets the interaction timestamp, so each interaction has the right time
830  for (auto i : interactionHandler)
831  {
832  i->setTimeStamp(i->getTimeStamp() - diff);
833  }
834 }
835 
841 {
842  if (newTMax >= 0)
843  {
844  timeMax_ = newTMax;
845  }
846  else
847  {
848  logger(ERROR, "Error in setTimeMax, new timeMax=% is not positive", newTMax);
849  }
850 }
851 
856 {
857  return timeMax_;
858 }
865 #ifdef CONTACT_LIST_HGRID
866 PossibleContactList& DPMBase::getPossibleContactList()
867 {
868  return possibleContactList;
869 }
870 #endif
871 
879 {
880  writeWallsVTK_ = writeWallsVTK;
881 }
882 
889 void DPMBase::setWallsWriteVTK(bool writeVTK)
890 {
892 }
893 
895 {
897 }
904 void DPMBase::setParticlesWriteVTK(bool writeParticlesVTK)
905 {
906  writeParticlesVTK_ = writeParticlesVTK;
907  if (writeParticlesVTK_)
908  {
910  }
911  delete vtkWriter_;
913 }
914 
918 void DPMBase::setSuperquadricParticlesWriteVTK(bool writeParticlesVTK)
919 {
920  writeSuperquadricParticlesVTK_ = writeParticlesVTK;
922  {
923  writeParticlesVTK_ = false;
924  }
925  delete vtkWriter_;
927 }
928 
936 {
937  return writeWallsVTK_;
938 }
939 
947 {
948  return writeParticlesVTK_;
949 }
950 
955 {
957 }
958 
972 {
973  if (newXMin <= getXMax())
974  {
975  min_.x() = newXMin;
976  }
977  else
978  {
979  logger(WARN, "Warning in setXMin(%): xMax=%", newXMin, getXMax());
980  }
981 }
982 
996 {
997  if (newYMin <= getYMax())
998  {
999  min_.y() = newYMin;
1000  }
1001  else
1002  {
1003  logger(WARN, "Warning in setYMin(%): yMax=%", newYMin, getYMax());
1004  }
1005 }
1006 
1020 {
1021 
1022  if (newZMin <= getZMax())
1023  {
1024  min_.z() = newZMin;
1025  }
1026  else
1027  {
1028  logger(WARN, "Warning in setZMin(%): zMax=%", newZMin, getZMax());
1029  }
1030 
1031 }
1032 
1043 void DPMBase::setMax(const Vec3D& newMax)
1044 {
1045  if (min_.x() > newMax.x() ||
1046  min_.y() > newMax.y() ||
1047  min_.z() > newMax.z())
1048  {
1049  logger(WARN, "Warning in setMax: upper bound is smaller"
1050  " than lower bound. (%,%,%) > (%,%,%)",
1051  min_.x(), min_.y(), min_.z(), newMax.x(), newMax.y(), newMax.z());
1052  }
1053  else
1054  {
1055  max_ = newMax;
1056  }
1057 }
1058 
1059 void DPMBase::setDomain(const Vec3D& min, const Vec3D& max)
1060 {
1061 
1062  logger.assert(min.X <= max.X, "lower x-bound (%) is larger than upper x-bound (%)", min.X, max.X);
1063  logger.assert(min.Y <= max.Y, "lower x-bound (%) is larger than upper x-bound (%)", min.Y, max.Y);
1064  logger.assert(min.Z <= max.Z, "lower x-bound (%) is larger than upper x-bound (%)", min.Z, max.Z);
1065  min_ = min;
1066  max_ = max;
1067 }
1068 
1079 void DPMBase::setMin(const Vec3D& newMin)
1080 {
1081  if (max_.x() < newMin.x() ||
1082  max_.y() < newMin.y() ||
1083  max_.z() < newMin.z())
1084  {
1085  logger(WARN, "Warning in setMin: lower bound is larger"
1086  " than upper bound. (%,%,%) < (%,%,%)",
1087  max_.x(), max_.y(), max_.z(), newMin.x(), newMin.y(), newMin.z());
1088  }
1089  else
1090  {
1091  min_ = newMin;
1092  }
1093 }
1094 
1099 void DPMBase::setMin(const Mdouble xMin, const Mdouble yMin, const Mdouble zMin)
1100 {
1101  setMin(Vec3D(xMin, yMin, zMin));
1102 }
1103 
1108 void DPMBase::setMax(const Mdouble xMax, const Mdouble yMax, const Mdouble zMax)
1109 {
1110  setMax(Vec3D(xMax, yMax, zMax));
1111 }
1112 
1127 {
1128 
1129  if (newXMax >= getXMin())
1130  {
1131  max_.x() = newXMax;
1132  }
1133  else
1134  {
1135  logger(WARN, "Warning in setXMax(%): xMax=%", newXMax, getXMin());
1136  }
1137 
1138 }
1139 
1153 {
1154 
1155  if (newYMax >= getYMin())
1156  {
1157  max_.y() = newYMax;
1158  }
1159  else
1160  {
1161  logger(WARN, "Warning in setYMax(%): yMax=%", newYMax, getYMin());
1162  }
1163 
1164 }
1165 
1179 {
1180  if (newZMax >= getZMin())
1181  {
1182  max_.z() = newZMax;
1183  }
1184  else
1185  {
1186  logger(WARN, "Warning in setZMax(%): zMax=%", newZMax, getZMin());
1187  }
1188 }
1189 
1196 {
1197  if (timeStep > 0.0)
1198  {
1199  timeStep_ = timeStep;
1200  }
1201  else
1202  {
1203  logger(ERROR, "Error in setTimeStep: new timeStep % is not positive", timeStep);
1204  }
1205 }
1206 
1212 {
1213  return timeStep_;
1214 }
1215 
1216 
1217 /* Allows user to set the number of omp threads */
1218 void DPMBase::setNumberOfOMPThreads(int numberOfOMPThreads)
1219 {
1220  logger.assert_always(numberOfOMPThreads>0, "Number of OMP threads must be positive");
1221  numberOfOMPThreads_ = numberOfOMPThreads;
1222 
1223  #ifdef MERCURY_USE_OMP
1224  if(numberOfOMPThreads > omp_get_max_threads()) {
1225  logger(INFO, "Number of omp threads set to the maximum number of threads allowed: %",
1226  omp_get_max_threads());
1227  numberOfOMPThreads_ = numberOfOMPThreads = omp_get_max_threads();
1228  }
1229  #pragma omp parallel num_threads(getNumberOfOMPThreads())
1230  {
1231  if (omp_get_thread_num()==0)
1232  std::cout << "Using " << omp_get_num_threads() << " of " << omp_get_max_threads() << " omp threads; testing thread";
1233  }
1234  #pragma omp parallel num_threads(getNumberOfOMPThreads())
1235  {
1236  std::cout << ' ' + std::to_string(omp_get_thread_num());
1237  }
1238  std::cout << '\n';
1239 
1240  #else
1241  logger(WARN, "You are setting the number of omp threads to %, but OMP is not turned on", getNumberOfOMPThreads());
1242  #endif
1243 }
1244 
1245 /* Returns the number of omp threads */
1247 {
1248  //logger.assert(numberOfOMPThreads_,"You need to set the number of OMP threads");
1249  return numberOfOMPThreads_;
1250 }
1251 
1261 {
1262  xBallsColourMode_ = newCMode;
1263 }
1264 
1271 {
1272  return xBallsColourMode_;
1273 }
1274 
1280 void DPMBase::setXBallsVectorScale(double newVScale)
1281 {
1282  xBallsVectorScale_ = newVScale;
1283 }
1284 
1291 {
1292  return xBallsVectorScale_;
1293 }
1294 
1307 void DPMBase::setXBallsAdditionalArguments(std::string xBallsAdditionalArguments)
1308 {
1309  xBallsAdditionalArguments_ = xBallsAdditionalArguments;
1310 }
1311 
1316 {
1318 }
1319 
1324 {
1325  xBallsScale_ = newScale;
1326 }
1327 
1333 {
1334  return xBallsScale_;
1335 }
1336 
1343 void DPMBase::setGravity(Vec3D newGravity)
1344 {
1345  gravity_ = newGravity;
1346 }
1347 
1352 {
1353  return gravity_;
1354 }
1355 
1363 void DPMBase::setDimension(unsigned int newDim)
1364 {
1365  setSystemDimensions(newDim);
1366  setParticleDimensions(newDim);
1367 }
1368 
1377 void DPMBase::setSystemDimensions(unsigned int newDim)
1378 {
1379  if (newDim >= 1 && newDim <= 3)
1380  systemDimensions_ = newDim;
1381  else
1382  {
1383  logger(ERROR, "Error in setSystemDimensions; newDim % is not 1, 2 or 3", newDim);
1384  }
1385 }
1386 
1390 unsigned int DPMBase::getSystemDimensions() const
1391 {
1392  return systemDimensions_;
1393 }
1394 
1408 void DPMBase::setParticleDimensions(unsigned int particleDimensions)
1409 {
1410  if (particleDimensions >= 1 && particleDimensions <= 3)
1411  {
1412  particleDimensions_ = particleDimensions;
1414  }
1415  else
1416  {
1417  logger(WARN, "Error in setParticleDimensions; particleDimensions % is not 1, 2 or 3, setting it to "
1418  "systemDimension now (%)", particleDimensions, systemDimensions_);
1420  }
1421 }
1422 
1428 {
1429  return particleDimensions_;
1430 }
1431 
1435 std::string DPMBase::getRestartVersion() const
1436 {
1437  return restartVersion_;
1438 }
1439 
1444 void DPMBase::setRestartVersion(std::string newRV)
1445 {
1446  restartVersion_ = newRV;
1447 }
1448 
1454 {
1455  return restarted_;
1456 }
1457 
1461 void DPMBase::setRestarted(bool newRestartedFlag)
1462 {
1463  restarted_ = newRestartedFlag;
1464  //setAppend(new_);
1465 }
1466 
1471 {
1472  return append_;
1473 }
1474 
1482 void DPMBase::setAppend(bool newAppendFlag)
1483 {
1484  append_ = newAppendFlag;
1485 }
1486 
1491 {
1492  Mdouble elasticEnergy = 0.0;
1493  // JMFT: Note that we do count the elastic energy of fixed particles here.
1494  for (const BaseInteraction* c : interactionHandler)
1495  {
1496  elasticEnergy += c->getElasticEnergy();
1497  }
1498  return elasticEnergy;
1499 }
1500 
1505 {
1506  Mdouble kineticEnergy = 0;
1507  for (const BaseParticle* const p : particleHandler)
1508  {
1509  if (!(p->isFixed()))
1510  {
1511  kineticEnergy += .5 * p->getMass() * p->getVelocity().getLengthSquared();
1512  }
1513  }
1514  return kineticEnergy;
1515 }
1516 
1522 {
1523  Mdouble gravitationalEnergy = 0;
1524  for (const BaseParticle* const p : particleHandler)
1525  {
1526  // Don't consider fixed particles. 'Fixed' particles aren't necessarily
1527  // stationary; it just means their position is prescribed.
1528  if (!(p->isFixed()))
1529  {
1530  gravitationalEnergy += p->getMass() * Vec3D::dot((-getGravity()), p->getPosition());
1531  }
1532  }
1533  return gravitationalEnergy;
1534 }
1535 
1538 {
1539  Mdouble ene_rot = 0;
1540  for (std::vector<BaseParticle*>::const_iterator it = particleHandler.begin(); it != particleHandler.end(); ++it)
1541  {
1542  // See above.
1543  if (!(*it)->isFixed())
1544  {
1545  // ene_rot += .5 * (*it)->getInertia() * (*it)->getAngularVelocity().getLengthSquared();
1546  }
1547  }
1548  return ene_rot;
1549 }
1550 
1553 }
1554 
1559 {
1560  /*
1561  double mass_sum = 0;
1562  for (std::vector<BaseParticle*>::const_iterator it = particleHandler.begin(); it != particleHandler.end(); ++it)
1563  if (!(*it)->isFixed())
1564  mass_sum += (*it)->getMass();
1565  return mass_sum;
1566  */
1567  return particleHandler.getMass();
1568 }
1569 
1575 {
1577 }
1578 
1585 {
1586  return particleHandler.getMomentum();
1587  /*
1588  Vec3D total_momentum = Vec3D(0,0,0);
1589  for (std::vector<BaseParticle*>::const_iterator it = particleHandler.begin(); it != particleHandler.end(); ++it)
1590  if (!(*it)->isFixed())
1591  total_momentum += (*it)->getMass() * (*it)->getVelocity();
1592  return total_momentum;
1593  */
1594 }
1595 
1599 /* JMFT: This information is placed on the last entry on each line of the .data
1600  * file. That space is, in general, reserved for 'additional' information.
1601  */
1603 {
1604 // return p.getSpecies()->getId(); // was getIndex()
1605  return p.getInfo();
1606 }
1607 
1622 {
1623  return (pI != pJ && pI->isInContactWith(pJ));
1624 }
1625 
1630 {
1631 }
1632 
1637 {
1638 }
1639 
1644 {
1645 }
1646 
1651 {
1652 }
1653 
1658 {
1659 }
1660 
1665 {
1666 }
1667 
1672 {
1673 }
1674 
1679 {
1680  return true;
1681 }
1682 
1693 {
1694 #ifdef MERCURY_USE_MPI
1695  //If only one core is used (i.e. domainHandler is empty) then the result is always true
1696  if (domainHandler.getSize() == 0)
1697  {
1698  return true;
1699  }
1700  //Get the current domain
1702 
1703  //Check if the particle is in the current domain
1704  if(domain->containsParticle(P))
1705  {
1706  //When adding a particle inside the domain, this should always be true
1707  return true;
1708  }
1709  else
1710  {
1711  //MPI particles that are inserted in the communication zone should still be inserted
1712  return (P->isMPIParticle());
1713  }
1714 #else
1715  return false;
1716 #endif
1717 }
1718 
1725 {
1726 
1727  bool insideCommunicationZone = false;
1728 #ifdef MERCURY_USE_MPI
1729  MPIContainer& communicator = MPIContainer::Instance();
1730 
1731  //Check for the current domain if the particle is within the communication domain
1732  int val = domainHandler.getCurrentDomain()->isInCommunicationZone(particle);
1733 
1734  //The root gathers all results
1735  int *list = nullptr;
1736  if (PROCESSOR_ID == 0)
1737  {
1738  list = new int [NUMBER_OF_PROCESSORS];
1739  }
1740  communicator.gather(val,list);
1741 
1742  //Compute the global value
1743  //if on any processor the val is true, we have to do the communcation step
1745  int result = 0;
1746  if (PROCESSOR_ID == 0)
1747  {
1748  for (int i = 0; i< NUMBER_OF_PROCESSORS; i++)
1749  {
1750  if (list[i] == 1)
1751  {
1752  result = 1;
1753  break;
1754  }
1755  }
1756  }
1757 
1758  //The root now tells the other processors what the global value for the interaction is
1759  communicator.broadcast(result);
1760 
1761  //Convert the result back to bool
1762  insideCommunicationZone = result;
1763 #endif
1764  return insideCommunicationZone;
1765 }
1766 
1772 {
1773 #ifdef MERCURY_USE_MPI
1774  //mpi particles only exist when there is more than one domain
1775  if (domainHandler.getSize() > 0)
1776  {
1777  //Add the particle to the mpi domain
1779  }
1780 
1781  //If periodic boundaries are present..
1782  if (periodicBoundaryHandler.getSize() > 0)
1783  {
1785  }
1786 #endif
1787 }
1788 
1798 {
1799 #ifdef MERCURY_USE_MPI
1800  if (NUMBER_OF_PROCESSORS == 1) { return; }
1801 
1802  //Check if the interactionRadius of the BaseParticle is larger than given in the domain
1805  {
1806  logger(VERBOSE,"Processor % | Updating mpi grid. Old interactionDistance: %, new interactionDistance %.",
1808 
1809  //Update the interactionDistance in the domain and periodicBoundaryHandler
1812 
1813  //Find new ghost particless
1816  }
1817 #endif
1818 }
1819 
1820 
1825 {
1826 }
1827 
1832 {
1833 }
1834 
1839 {
1840 }
1841 
1846 {
1848 }
1849 
1854 {
1855  //cgHandler.evaluate();
1856 }
1857 
1859 {
1861  {
1862  c->gatherContactStatistics();
1863  }
1864 }
1865 
1869 void DPMBase::gatherContactStatistics(unsigned int index1 UNUSED, int index2 UNUSED, Vec3D Contact UNUSED,
1870  Mdouble delta UNUSED, Mdouble ctheta UNUSED, Mdouble fdotn UNUSED,
1871  Mdouble fdott UNUSED, Vec3D P1_P2_normal_ UNUSED, Vec3D P1_P2_tangential UNUSED)
1872 {
1873 }
1874 
1879 {
1880 }
1881 
1886 {
1887  cgHandler.finish();
1888 }
1889 
1894 {
1895 }
1896 
1901 {
1902 }
1903 
1908 {
1909 }
1910 
1920 void DPMBase::setFixedParticles(unsigned int n)
1921 {
1922  for (unsigned int i = 0; i < std::min(particleHandler.getSize(), n); i++)
1924 }
1925 
1931 {
1932 #ifdef MERCURY_USE_MPI
1933  MPIContainer& communicator = MPIContainer::Instance();
1934  if (communicator.getProcessorID() == 0)
1935  {
1936 #endif
1937  std::cout << "t=" << std::setprecision(3) << std::left << std::setw(6) << getTime()
1938  << ", tmax=" << std::setprecision(3) << std::left << std::setw(6) << getTimeMax()
1939  << std::endl;
1940  std::cout.flush();
1941 #ifdef MERCURY_USE_MPI
1942  }
1943 #endif
1944 }
1945 
1954 {
1955  return true;
1956 }
1957 
1962 {
1963 }
1964 
1977 void DPMBase::writeEneHeader(std::ostream& os) const
1978 {
1979  //only write if we don't restart
1980  if (getAppend())
1981  return;
1982 
1984 
1987  long width = os.precision() + 6;
1988  os << std::setw(width)
1989  << "time " << std::setw(width)
1990  << "gravitEnergy " << std::setw(width) //gravitational potential energy
1991  << "traKineticEnergy " << std::setw(width) //translational kinetic energy
1992  << "rotKineticEnergy " << std::setw(width) //rotational kE
1993  << "elasticEnergy " << std::setw(width)
1994  << "centerOfMassX " << std::setw(width)
1995  << "centerOfMassY " << std::setw(width)
1996  << "centerOfMassZ\n";
1997 }
1998 
2006 void DPMBase::writeFstatHeader(std::ostream& os) const
2007 {
2008 
2009  // line #1: time, volume fraction
2010  // line #2: wall box: wx0, wy0, wz0, wx1, wy1, wz1
2011  // line #3: radii-min-max & moments: rad_min, rad_max, r1, r2, r3, r4
2012  os << "#"
2013  << " " << getTime()
2014  << " " << 1 //marker that these are fstat files with contact point instead of center point
2015  << '\n';
2016  os << "#"
2017  << " " << getXMin()
2018  << " " << getYMin()
2019  << " " << getZMin()
2020  << " " << getXMax()
2021  << " " << getYMax()
2022  << " " << getZMax()
2023  << '\n';
2024  os << "#"
2025  << " ";
2026 
2027  if (!(particleHandler.getSmallestParticleLocal() == nullptr))
2028  {
2030  }
2031  else
2032  {
2033  os << std::numeric_limits<double>::quiet_NaN();
2034  }
2035  os << " ";
2036  if (!(particleHandler.getLargestParticleLocal() == nullptr))
2037  {
2039  }
2040  else
2041  {
2042  os << std::numeric_limits<double>::quiet_NaN();
2043  }
2044 
2045  os << " " << 0
2046  << " " << 0
2047  << " " << 0
2048  << " " << 0
2049  << '\n';
2050  //B: write data
2052  {
2053  c->writeToFStat(os, getTime());
2054  }
2055  //os << std::flush;
2056 }
2057 
2067 void DPMBase::writeEneTimeStep(std::ostream& os) const
2068 {
2071  writeEneHeader(os);
2072 
2073  const Mdouble m = particleHandler.getMass();
2075  //Ensure the numbers fit into a constant width column: for this we need the precision given by the operating system,
2076  //plus a few extra characters for characters like a minus and scientific notation.
2077  const static int width = os.precision() + 6;
2078  os << std::setw(width) << getTime()
2079  << " " << std::setw(width) << -Vec3D::dot(getGravity(), com)
2080  << " " << std::setw(width) << particleHandler.getKineticEnergy()
2081  << " " << std::setw(width) << particleHandler.getRotationalEnergy()
2082  << " " << std::setw(width) << getElasticEnergy()
2083  // we need to write x, y and z coordinates separately, otherwise the width of the columns is incorrect
2084  << " " << std::setw(width)
2085  << (m == 0 ? constants::NaN : com.X / m) //set to nan because 0/0 implementation in gcc and clang differs
2086  << " " << std::setw(width) << (m == 0 ? constants::NaN : com.Y / m)
2087  << " " << std::setw(width) << (m == 0 ? constants::NaN : com.Z / m)
2088  << std::endl;
2089 }
2090 
2092 {
2093  static bool writeWall = true;
2094  if (getWallsWriteVTK() == FileType::ONE_FILE && writeWall)
2095  {
2097  writeWall=false;
2101  } // else do nothing
2102 
2104  {
2105  vtkWriter_->writeVTK();
2106  } // else do nothing
2107 
2109  {
2111  }
2112 
2114  {
2116  }
2117 
2118  //only write once
2119  bool writePython = getParticlesWriteVTK() || getWallsWriteVTK() != FileType::NO_FILE ||
2121  if (writePython && getTime() == 0)
2122  {
2124  }
2125 }
2126 
2128 {
2129 #ifdef MERCURY_USE_MPI
2130  if (PROCESSOR_ID == 0)
2131  {
2132  logger(INFO, "Writing python script for paraview visualisation");
2133 #else
2134  logger(INFO, "Writing python script for paraview visualisation");
2135 #endif
2136 
2137  std::string script = "#script to visualise the output of data2pvd of MercuryDPM in paraview.\n"
2138  "#usage: change the path below to your own path, open paraview\n"
2139  "#Tools->Python Shell->Run Script->VisualisationScript.py\n"
2140  "#or run paraview --script=VisualisationScript.py \n"
2141  "\n"
2142  "from paraview.simple import *\n"
2143  "import os\n"
2144  "import glob\n"
2145  "os.chdir('" + helpers::getPath() + "')\n\n";
2146 
2147 #ifdef MERCURY_USE_MPI
2148  for (int i = 0; i < NUMBER_OF_PROCESSORS; i++)
2149  {
2150 #endif
2151  if (getParticlesWriteVTK())
2152  {
2153  script += "#Load data in any order\n";
2154 #ifdef MERCURY_USE_MPI
2155  if (NUMBER_OF_PROCESSORS > 1)
2156  {
2157  script += "Data = glob.glob('./" + getName() + "Processor_" + std::to_string(i) + "_Particle_*.vtu')\n";
2158  }
2159  else
2160  {
2161  script += "Data = glob.glob('./" + getName() + "Particle_*.vtu')\n";
2162  }
2163 #else
2164  script += "Data = glob.glob('./" + getName() + "Particle_*.vtu')\n";
2165 #endif
2166  script += "\n"
2167  "#Find the maximum timestep\n"
2168  "maxTime = 0\n"
2169  "for fileName in Data:\n"
2170  "\ttokens1 = fileName.split('.')\n"
2171  "\ttokens2 = tokens1[1].split('_')\n"
2172  "\tif int(tokens2[-1]) > maxTime:\n"
2173  "\t\tmaxTime = int(tokens2[-1])\n"
2174  "print str(maxTime)\n"
2175  "\n"
2176  "#Create correct order of time steps\n"
2177  "DataSorted = []\n"
2178  "for x in range(0,maxTime+1):\n";
2179 #ifdef MERCURY_USE_MPI
2180  if (NUMBER_OF_PROCESSORS > 1)
2181  {
2182  script += "\tDataSorted.append('./" + getName() + "Processor_" + std::to_string(i) + "_Particle_' + " + "str(x)" + " + '.vtu')\n";
2183  }
2184  else
2185  {
2186  script += "\tDataSorted.append('./" + getName() + "Particle_' + " + "str(x)" + " + '.vtu')\n";
2187  }
2188 #else
2189  script += "\tDataSorted.append('./" + getName() + "Particle_' + " + "str(x)" + " + '.vtu')\n";
2190 #endif
2191 
2192  script += "\n"
2193  "#Load the data and visualise it in paraview\n"
2194  "particles = XMLUnstructuredGridReader(FileName=DataSorted)\n"
2195  "glyphP = Glyph(particles)\n"
2196  "glyphP.GlyphType = 'Sphere'\n"
2197  "glyphP.Scalars = 'Radius'\n"
2198  "glyphP.Vectors = 'None'\n"
2199  "glyphP.ScaleMode = 'scalar'\n"
2200  "glyphP.ScaleFactor = 2\n"
2201  "glyphP.GlyphMode = 'All Points'\n"
2202  "Show(glyphP)\n\n";
2203  }
2205  {
2206  script += "walls = XMLUnstructuredGridReader(FileName=glob.glob('./"
2207  + getName() + "Wall_*.vtu'))\n"
2208  "Show(walls)\n\n";
2209  }
2212  {
2213  script += "interactions = XMLUnstructuredGridReader(FileName=glob.glob('./"
2214  + getName() + "Interaction_*.vtu'))\n"
2215  "glyphI = Glyph(interactions)\n"
2216  "glyphI.GlyphType = 'Sphere'\n"
2217  "glyphI.Scalars = 'Cylinder'\n"
2218  "glyphI.Vectors = 'None'\n"
2219  "glyphI.ScaleMode = 'scalar'\n"
2220  "glyphI.ScaleFactor = 10\n" //5 times too large
2221  "glyphI.GlyphMode = 'All Points'\n"
2222  "Show(glyphI)\n\n";
2223  }
2224 #ifdef MERCURY_USE_MPI
2225  } // end of loop over number of processors
2226 #endif
2227  script += "Render()\n"
2228  "ResetCamera()\n";
2229 
2230  helpers::writeToFile(getName() + ".py", script);
2231 #ifdef MERCURY_USE_MPI
2232  } // end of communicator is root statement
2233 #endif
2234 }
2235 
2236 
2240 void DPMBase::outputXBallsData(std::ostream& os) const
2241 {
2242 
2243 
2244  //Set the correct formation based of dimension if the formation is not specified by the user
2245 
2246  unsigned int format;
2247  switch (getSystemDimensions())
2248  {
2249  case 2:
2250  format = 8;
2251  break;
2252  case 3:
2253  format = 14;
2254  break;
2255  default:
2256  std::cerr << "Unknown system dimension" << std::endl;
2257  exit(-1);
2258  }
2259 
2260  unsigned int numberOfParticles = particleHandler.getNumberOfRealObjectsLocal();
2261 
2262  // This outputs the location of walls and how many particles there are to file this is required by the xballs plotting
2263  if (format != 14) // dim = 1 or 2
2264  {
2265  os << numberOfParticles
2266  << " " << getTime()
2267  << " " << getXMin()
2268  << " " << getYMin()
2269  << " " << getXMax()
2270  << " " << getYMax()
2271  << " " << std::endl;
2272  }
2273  else
2274  {
2275  //dim==3
2276  os << numberOfParticles
2277  << " " << getTime()
2278  << " " << getXMin()
2279  << " " << getYMin()
2280  << " " << getZMin()
2281  << " " << getXMax()
2282  << " " << getYMax()
2283  << " " << getZMax()
2284  << " " << std::endl;
2285  }
2286 
2287  // This outputs the particle data
2288  for (unsigned int i = 0; i < particleHandler.getSize(); i++)
2289  {
2290 #ifdef MERCURY_USE_MPI
2292  {
2293  outputXBallsDataParticle(i, format, os);
2294  }
2295 #else
2296  outputXBallsDataParticle(i, format, os);
2297 #endif
2298  }
2299 #ifdef DEBUG_OUTPUT
2300  std::cerr << "Have output the properties of the problem to disk " << std::endl;
2301 #endif
2302 }
2303 
2319 bool DPMBase::readDataFile(std::string fileName, unsigned int format)
2320 {
2321  //default value: dataFile.getFullName()
2322  if (!fileName.compare(""))
2323  fileName = dataFile.getFullName();
2324 
2325  std::string oldFileName = dataFile.getName();
2326  unsigned oldCounter = dataFile.getCounter();
2327  //Updates the name of the data file to the user-input from the argument.
2328  dataFile.setName(fileName);
2329  //opens a filestream of the input type
2330  dataFile.open(std::fstream::in);
2331  //Checks if the file has been successfully opened...
2332  if (!dataFile.getFstream().is_open() || dataFile.getFstream().bad())
2333  {
2334  //...and if not, ends the function and returns "false"
2335  logger(WARN, "Loading data file % failed.", fileName);
2336  return false;
2337  }
2338 
2339  //retrieves and saves the "FileType" of the file
2340  FileType fileTypeData = dataFile.getFileType();
2342  readNextDataFile(format);
2343  dataFile.setFileType(fileTypeData);
2344  dataFile.close();
2345  dataFile.setName(oldFileName);
2346  dataFile.setCounter(oldCounter);
2347  return true;
2348 }
2349 
2355 bool DPMBase::readParAndIniFiles(const std::string fileName)
2356 {
2357  //Opens the par.ini file
2358  std::fstream file;
2359  file.open(fileName, std::fstream::in);
2360  if (!file.is_open() || file.bad())
2361  {
2362  //std::cout << "Loading par.ini file " << filename << " failed" << std::endl;
2363  return false;
2364  }
2365 
2366  Mdouble doubleValue;
2367  int integerValue;
2368 
2369  // inputfile par.ini
2370  // line 1 =============================================================
2371  // Example: 1 1 0
2372  // 1: integer (0|1) switches from non-periodic to periodic
2373  // integer (5|6) does 2D integration only (y-coordinates fixed)
2374  // and switches from non-periodic to periodic
2375  // integer (11) uses a quarter system with circular b.c.
2376  file >> integerValue;
2377  //~ std::cout << "11" << integerValue << std::endl;
2378  if (integerValue == 0)
2379  {
2380  InfiniteWall w0;
2382  w0.set(Vec3D(-1, 0, 0), Vec3D(getXMin(), 0, 0));
2384  w0.set(Vec3D(1, 0, 0), Vec3D(getXMax(), 0, 0));
2386  w0.set(Vec3D(0, -1, 0), Vec3D(0, getYMin(), 0));
2388  w0.set(Vec3D(0, 1, 0), Vec3D(0, getYMax(), 0));
2390  w0.set(Vec3D(0, 0, -1), Vec3D(0, 0, getZMin()));
2392  w0.set(Vec3D(0, 0, 1), Vec3D(0, 0, getZMax()));
2394  }
2395  else if (integerValue == 1)
2396  {
2397  PeriodicBoundary b0;
2398  b0.set(Vec3D(1, 0, 0), getXMin(), getXMax());
2400  b0.set(Vec3D(0, 1, 0), getYMin(), getYMax());
2402  b0.set(Vec3D(0, 0, 1), getZMin(), getZMax());
2404  }
2405  else if (integerValue == 5)
2406  {
2407  InfiniteWall w0;
2409  w0.set(Vec3D(-1, 0, 0), Vec3D(-getXMin(), 0, 0));
2411  w0.set(Vec3D(1, 0, 0), Vec3D(getXMax(), 0, 0));
2413  w0.set(Vec3D(0, -1, 0), Vec3D(0, -getYMin(), 0));
2415  w0.set(Vec3D(0, 1, 0), Vec3D(0, getYMax(), 0));
2417 
2418  }
2419  else if (integerValue == 6)
2420  {
2421  PeriodicBoundary b0;
2422  b0.set(Vec3D(1, 0, 0), getXMin(), getXMax());
2424  b0.set(Vec3D(0, 1, 0), getYMin(), getYMax());
2426  b0.set(Vec3D(0, 0, 1), getZMin(), getZMax());
2428  }
2429  else
2430  {
2431  std::cerr << "Error in par.ini: line 1, value 1 is " << integerValue << std::endl;
2432  exit(-1);
2433  }
2434 
2435  // 2: integer (0|1) dont use | use the search pattern for linked cells
2436  file >> integerValue; //ignore
2437 
2438  // 3: real - gravity in z direction: positive points upwards
2439  file >> doubleValue;
2440  setGravity(Vec3D(0.0, 0.0, doubleValue));
2441 
2442  // line 2 =============================================================
2443  // Example: -10000 .5e-2
2444  // 1: time end of simulation - (negative resets start time to zero
2445  // and uses -time as end-time)
2446  file >> doubleValue;
2447  if (doubleValue < 0)
2448  setTime(0);
2449  setTimeMax(fabs(doubleValue));
2450 
2451  // 2: time-step of simulation
2452  file >> doubleValue;
2453  setTimeStep(doubleValue);
2454 
2455  // line 3 =============================================================
2456  // Example: 1e-1 100
2457  file >> doubleValue;
2458  if (doubleValue >= 0)
2459  {
2460  // 1: time-step for output on time-series protocoll file -> "ene"
2461  unsigned int savecount = static_cast<unsigned int>(round(doubleValue / getTimeStep()));
2462  setSaveCount(savecount);
2463 
2464  // 2: time-step for output on film (coordinate) file -> "c3d"
2465  // (fstat-output is coupled to c3d-output time-step)
2466  file >> doubleValue;
2467  savecount = static_cast<unsigned int>(round(doubleValue / getTimeStep()));
2468  dataFile.setSaveCount(savecount);
2469  fStatFile.setSaveCount(savecount);
2470  }
2471  else
2472  {
2473  // or: ---------------------------------------------------------------
2474  // 1: negative number is multiplied to the previous log-output-time
2475  // 2: requires initial log-output time
2476  // 3: negative number is multiplied to the previous film-output-time
2477  // 4: requires initial film-output time
2478  std::cerr << "Error in par.ini: line 3, value 1 is " << doubleValue << std::endl;
2479  exit(-1);
2480  }
2481 
2482  // line 4 =============================================================
2483  // Example: 2000 1e5 1e3 1e2
2484  // 1: particle density (mass=4/3*constants::pi*density*rad^3)
2485  file >> doubleValue;
2486 
2487  //clear species handler
2491 
2493  S->setDensity(doubleValue);
2494 
2495  // 2: linear spring constant
2496  file >> doubleValue;
2497  S->setStiffness(doubleValue);
2498 
2499  // 3: linear dashpot constant
2500  file >> doubleValue;
2501  S->setDissipation(doubleValue);
2502 
2503  // 4: background damping dashpot constant
2504  file >> doubleValue;
2505  if (doubleValue != 0.0)
2506  std::cerr << "Warning in par.ini: ignored background damping " << doubleValue << std::endl;
2507 
2508  // line 5 =============================================================
2509  // Example: 0 0
2510  // 1: growth rate: d(radius) = xgrow * dt
2511  file >> doubleValue;
2512  if (doubleValue != 0.0)
2513  std::cerr << "Warning in par.ini: ignored growth rate " << doubleValue << std::endl;
2514 
2515  // 2: target volume_fraction
2516  file >> doubleValue;
2517  if (doubleValue != 0.0)
2518  std::cerr << "Warning in par.ini: ignored target volume_fraction " << doubleValue << std::endl;
2519 
2520  file.close();
2521  //std::cout << "Loaded par.ini file " << filename << std::endl;
2522  return true;
2523 }
2524 
2539 {
2541  {
2542  while (true)// This true corresponds to the if s
2543  {
2544  dataFile.open();
2545  //check if file exists and contains data
2546  int N;
2547  dataFile.getFstream() >> N;
2548  if (dataFile.getFstream().eof() || dataFile.getFstream().peek() == -1)
2549  {
2550  std::cout << "file " << dataFile.getName() << " not found" << std::endl;
2551  return false;
2552  }
2553  //check if tmin condition is satisfied
2554  Mdouble t;
2555  dataFile.getFstream() >> t;
2556  if (t > tMin)
2557  {
2558  //set_file_counter(get_file_counter()-1);
2559  return true;
2560  }
2561  if (verbose)
2562  std::cout << "Jumping file counter: " << dataFile.getCounter() << std::endl;
2563  }
2564  }
2565  return true;
2566 }
2567 
2575 bool DPMBase::readNextDataFile(unsigned int format)
2576 {
2577  dataFile.open(std::fstream::in);
2578  //logger(INFO,"Reading %",dataFile.getFullName());
2579  //fStatFile.open();
2580  //Set the correct format based of dimension if the format is not explicitly specified by the user
2581  if (format == 0)
2582  {
2583  //checking the dimensionality of the system
2585  switch (getSystemDimensions())
2586  {
2587  case 1:
2588  //if 2D, sets format 8 (data file has 8 columns for 2D)
2589  case 2:
2590  format = 8;
2591  break;
2592  case 3:
2593  //if 3D, sets format 14 (data file has 14 columns for 3D)
2594  format = 14;
2595  break;
2596  }
2597  //end case
2598  }
2599  //end if
2600 
2601  // read in the particle number (as a double)
2602  double doubleN = -1;
2603  dataFile.getFstream() >> doubleN;
2604 
2605  // If N cannot be read, we reached the end of the file
2606  if (doubleN == -1) return false;
2607 
2608  // converting N to an integer; skipping the line if there is a problem (this happens when there is a corrupt data file)
2609  unsigned N = doubleN;
2610  while (doubleN != N) {
2611  std::string dummy;
2612  getline(dataFile.getFstream(),dummy,'\n');
2613  logger(WARN,"Skipping bad line in data file: % %",doubleN, dummy);
2614  dataFile.getFstream() >> doubleN;
2615  N = doubleN;
2616  }
2617 
2618  //store the parameters you want to preserve:
2619  const size_t nHistory = std::min(N,particleHandler.getSize());
2620  std::vector<const ParticleSpecies*> species(nHistory);
2621  std::vector<bool> fix(nHistory);
2622  for (size_t i=0; i<nHistory; ++i) {
2624  species[i] = p->getSpecies();
2625  fix[i] = p->isFixed();
2626  //store from reading the restart file which particles are fixed and which species each particle has
2627  }
2628 
2629  BaseParticle* p;
2630  if (particleHandler.getSize() == 0) {
2631  logger.assert_always(speciesHandler.getSize()>0,"readData: species needs to be set first");
2633  p->unfix();
2634  } else {
2635  p = particleHandler.getObject(0)->copy();
2636  }
2637 
2638  //empty the particle handler
2641 
2642  //now fill it again with the read-in particles
2643  std::stringstream line;
2645  Mdouble radius;
2646  Vec3D position, velocity, angle, angularVelocity;
2647  size_t indSpecies;
2648  //read all other data available for the time step
2649  if (format==7) {
2650  line >> time_ >> min_.x() >> min_.y() >> min_.z() >> max_.x() >> max_.y() >> max_.z();
2651  for (size_t i = 0; i < N; ++i) {
2653  line >> position.X >> position.Z >> position.Y >> velocity.X >> velocity.Z
2654  >> velocity.Y >> radius
2655  >> indSpecies;
2656  p->setPosition(position);
2657  p->setVelocity(velocity);
2658  p->setOrientation({1, 0, 0, 0});
2659  p->setAngularVelocity({0.0, 0.0, 0.0});
2660  p->setRadius(radius);
2662  p->setSpecies(speciesHandler.getObject(indSpecies));
2663  else if (i < nHistory)
2664  p->setSpecies(species[i]);
2665  if (i < nHistory && fix[i])
2666  p->fixParticle();
2668  p->unfix();
2669  }
2670  } else if (format==8) {
2671  line >> time_ >> min_.x() >> min_.y() >> max_.x() >> max_.y();
2672  min_.z() = 0.0; max_.z() = 0.0;//For 2d functions we define Z to be of no interest
2673  for (size_t i = 0; i < N; ++i) {
2675  line >> position.X >> position.Y >> velocity.X >> velocity.Y >> radius >> angle.Z >> angularVelocity.Z >> indSpecies;
2676  Quaternion q;
2677  q.setAngleZ(angle.Z);
2678  p->setPosition(position);
2679  p->setVelocity(velocity);
2680  p->setOrientation(q);
2681  p->setAngularVelocity(-angularVelocity);
2682  p->setRadius(radius);
2684  else if (i < nHistory) p->setSpecies(species[i]);
2685  if (i < nHistory && fix[i]) p->fixParticle();
2687  p->unfix();
2688  } //end for all particles
2689  } else if (format==14) {
2690  //This is a 3D format_
2691  //@TODO: Check bounds or get rid of this function
2692  line >> time_ >> min_.x() >> min_.y() >> min_.z() >> max_.x() >> max_.y() >> max_.z();
2693  for (size_t i = 0; i < N; ++i) {
2695  line >> position >> velocity >> radius >> angle >> angularVelocity >> indSpecies;
2696 #ifdef MERCURY_USE_MPI
2697  //This is required for the CG tool. When reading the data file it is neseccary to know if a particle is an MPIParticle or not
2698  bool isMPIParticle = false;
2699  bool isPeriodicGhostParticle = false;
2700  if (NUMBER_OF_PROCESSORS > 1)
2701  {
2703  line >> isMPIParticle >> isPeriodicGhostParticle;
2704  }
2705 #endif
2706  p->setPosition(position);
2707  p->setVelocity(velocity);
2708  p->setOrientationViaEuler(angle);
2709  p->setAngularVelocity(angularVelocity);
2710  p->setRadius(radius);
2712  if (indSpecies<speciesHandler.getSize()) {
2713  p->setSpecies(speciesHandler.getObject(indSpecies));
2714  } else {
2715  logger(WARN, "Read in bad species data; species is not set");
2716  }
2717  } else if (i < nHistory)
2718  p->setSpecies(species[i]);
2719  if (i < nHistory && fix[i])
2720  p->fixParticle();
2721 #ifdef MERCURY_USE_MPI
2723  {
2724  p->setMPIParticle(isMPIParticle);
2725  p->setPeriodicGhostParticle(isPeriodicGhostParticle);
2726  }
2727 #endif
2729  p->unfix();
2730  } //end read into existing particles logger(INFO, "read % particles", particleHandler.getNumberOfObjects());
2731  } else if (format==15) {
2732  line >> time_ >> min_.x() >> min_.y() >> min_.z() >> max_.z() >> max_.y() >> max_.z();
2733  for (size_t i = 0; i < N; ++i) {
2735  line >> position >> velocity >> radius >> angle >> angularVelocity >> indSpecies >> indSpecies;
2736  Quaternion q;
2737  q.setEuler(angle);
2738  p->setPosition(position);
2739  p->setVelocity(velocity);
2740  p->setOrientation(q);
2741  p->setAngularVelocity(angularVelocity);
2742  p->setRadius(radius);
2744  else if (i < nHistory) p->setSpecies(species[i]);
2745  if (i < nHistory && fix[i]) p->fixParticle();
2747  } //end for all particles
2748  } //end if format
2749 
2751  return true;
2752 }
2753 
2755 {
2756  fStatFile.open(std::fstream::in);
2757  std::string line;
2758  std::fstream& in = fStatFile.getFstream();
2759  // read the first three lines
2760  getline(in, line);
2761  getline(in, line);
2762  getline(in, line);
2763  Mdouble time;
2764  unsigned int indexP;
2765  int indexI; //could be negative
2766  unsigned counter = 0;
2768  while ((in.peek() != -1) && (in.peek() != '#'))
2769  {
2770  /* # 1: time
2771  # 2: particle Number i
2772  # 3: contact partner j (particles >= 0, walls < 0)
2773  # 4: x-position \
2774  # 5: y-position > of the contact point (I hope)
2775  # 6: z-position /
2776  # 7: delta = overlap at the contact
2777  # 8: ctheta = length of the tangential spring
2778  # 9: P1_P2_normal force |f^n|
2779  # 10: remaining (tangential) force |f^t|=|f-f^n|
2780  # 11-13: P1_P2_normal unit vector nx, ny, nz
2781  # 14-16: tangential unit vector tx, ty, tz
2782  */
2783  in >> time >> indexP >> indexI;
2784  BaseParticle* P = particleHandler.getObject(indexP);
2785  BaseInteraction* C;
2786  if (indexI >= 0)
2787  {
2788  //read only one of the two fstat lines reported
2789  if (indexI >= indexP)
2790  {
2791  // particle pair contact
2792  BaseParticle* I = particleHandler.getObject(static_cast<const unsigned int>(indexI));
2794  C->setFStatData(in, P, I);
2795  // skip next line
2796  //in.ignore(256, '\n');
2797  }
2798  }
2799  else
2800  {
2801  // wall-particle contact
2802  while (wallHandler.getNumberOfObjects() <= -indexI - 1)
2803  {
2805  logger(WARN, "Added new wall because .fstat file indicates contact with wall % that doesn't exist",
2806  -indexI - 1);
2807  }
2808  BaseWall* I = wallHandler.getObject(static_cast<const unsigned int>(-indexI - 1));
2810  C->setFStatData(in, P, I);
2811  }
2812  counter++;
2813  //skip the rest of the line, to allow additional output and to ignore the second time a particle-particle contact is printed
2814  in.ignore(256, '\n');
2815  }
2816  //logger(INFO,"read % contacts at t = % (N=%)",counter,timeStamp,interactionHandler.getNumberOfObjects());
2817  //interactionHandler.write(std::cout);
2818  //logger(INFO,"normal % %",interactionHandler.getObject(0)->getNormal(),interactionHandler.getObject(0)->getContactPoint());
2819  //logger(INFO,"normal % %",interactionHandler.getLastObject()->getNormal(),interactionHandler.getLastObject()->getContactPoint());
2820 }
2821 
2830 {
2832  {
2833  //logger(DEBUG, "Writing restart file %th time step",getNumberOfTimeSteps());
2835  restartFile.close();
2836  }
2837 }
2838 
2840 {
2842  {
2844  dataFile.close();
2845  }
2846 }
2847 
2849 {
2851  {
2852  //If the file type is "multiple files, writes a header for each individual files. If not, only writes for the first time step
2854  eneFile.close();
2855  }
2856 }
2857 
2859 {
2861  {
2863  //fStatFile.getFstream().ignore(2000,'\t');
2864  fStatFile.close();
2865  }
2866 }
2867 
2875  logger.assert_always(speciesHandler.getSize()>0,"There needs to be at least one species");
2877  SphericalParticle p(s);
2879  Mdouble r = cbrt(getTotalVolume())/N;
2880  b.set(p,100,getMin(),getMax(),{0,0,0},{0,0,0},r,r);
2881  b.insertParticles(this);
2882  logger(INFO,"Inserted % particles",particleHandler.getSize());
2883  //setTimeMax(0);
2884  //solve();
2885 }
2886 
2887 
2897 {
2898  //Assuming a filename corresponding to "restartFile" has already been established,
2899  //opens an input filestream to the relevant file
2900  if (restartFile.open(std::fstream::in))
2901  {
2902  //reads the input stream line-by-line
2903  read(restartFile.getFstream(), opt);
2904  logger(INFO, "Loaded restart file %", restartFile.getFullName());
2905  restartFile.close();
2906  //sets the flag such that other functions can tell that the file has been restarted
2907  //e.g. does not run "setUpInitialConditions" or add headers to the .ene files etc.
2908  setRestarted(true);
2909  return true;
2910  }
2911  else /* if the file could not be opened */
2912  {
2913  logger(INFO, "% could not be loaded.", restartFile.getFullName());
2914  return false;
2915  }
2916 }
2917 
2925 int DPMBase::readRestartFile(std::string fileName, ReadOptions opt)
2926 {
2927  //add ".restart" if necessary
2928  if (fileName.find(".restart") == std::string::npos)
2929  {
2930  fileName.append(".restart");
2931  }
2932 
2933  //If padded or numbered files are used, we need to extract the counter and remove it from the filename
2934  //First find the last point in the restart name
2935  unsigned int pos = fileName.find('.');
2936  while (fileName.find('.', pos + 1) != std::string::npos)
2937  {
2938  pos = fileName.find('.', pos + 1);
2939  }
2940  //If the next char after the last . is a digit we are using numbered files
2941  std::string counter;
2942  if (isdigit(fileName[pos + 1]))
2943  {
2944  for (int i = pos + 1; i < fileName.length(); i++)
2945  {
2946  counter.push_back(fileName[i]);
2947  }
2948  //Set counter in restart file
2949  restartFile.setCounter(std::stoi(counter));
2950  logger(INFO, "Counter: %", std::stoi(counter));
2951  }
2952 
2953 #ifdef MERCURY_USE_MPI
2954  //Correct for the processor number
2955  if (NUMBER_OF_PROCESSORS > 1 && !helpers::fileExists(fileName))
2956  {
2957  //Modify file name
2958  const unsigned int length = fileName.length();
2959  if (isdigit(fileName[pos + 1]))
2960  {
2961  for (int i = pos + 1; i < length + 1; i++)
2962  {
2963  fileName.pop_back();
2964  }
2965  }
2966  fileName.append(std::to_string(PROCESSOR_ID));
2967  if (counter.size() > 0)
2968  {
2969  fileName.append(".");
2970  fileName.append(counter);
2971  }
2972  }
2973 #endif
2974 
2975  restartFile.setName(fileName);
2976 
2977  logger(INFO, "Restarting from %", fileName);
2978  return readRestartFile(opt);
2979 }
2980 
2995 {
2996  //Does not compute forces if particles are fixed
2997  //this is necessary because the rough bottom allows overlapping fixed particles
2998  if (P1->isFixed() && P2->isFixed())
2999  {
3000  return;
3001  }
3002 //Ensures that interactions between the "ghost" particles used to implement periodic behaviour
3003  //are not included in calculations
3004  //i.e. ends the function if both particles are "ghosts".
3005  if ((P1->getPeriodicFromParticle() != nullptr) && (P2->getPeriodicFromParticle() != nullptr))
3006  {
3007  return;
3008  }
3009 //if statement below ensures that the PI has the lower id than PJ
3010  BaseParticle* PI, * PJ;
3011  if (P1->getId() > P2->getId())
3012  {
3013  PI = P2;
3014  PJ = P1;
3015  }
3016  else
3017  {
3018  PI = P1;
3019  PJ = P2;
3020  }
3021  //checks if the two particles are interacting
3022  //("getInteractionWith" returns the relevant pointer if PI and PJ are interacting,
3023  //zero if not)
3024  //if statement above ensures that the PI has the lower id than PJ
3027  if (i!= nullptr) {
3028  //calculates the force corresponding to the interaction
3029  i->computeForce();
3030 
3031  //Applies the relevant calculated forces to PI and PJ
3032  PI->addForce(i->getForce());
3033  PJ->addForce(-i->getForce());
3034 
3035  //checks if particle rotation is turned on...
3036  if (getRotation()) {
3037  //...and, if so, performs equivalent calculations for the torque as were
3038  //performed for the force.
3039  PI->addTorque(i->getTorque() - Vec3D::cross(PI->getPosition() - i->getContactPoint(), i->getForce()));
3040  PJ->addTorque(-i->getTorque() + Vec3D::cross(PJ->getPosition() - i->getContactPoint(), i->getForce()));
3041  }
3042  }
3043 }
3044 
3050 {
3051  //Checks that the current particle is not "fixed"
3052  //and hence infinitely massive!
3053  if (!CI->isFixed())
3054  {
3055  // Applying the force due to gravity (F = m.g)
3056  CI->addForce(getGravity() * CI->getMass());
3057  // Still calls this in compute External Forces.
3058  // computeForcesDueToWalls(CI);
3059  }
3060 }
3061 
3070 {
3071  //No need to compute interactions between periodic particle images and walls
3072  if (pI->getPeriodicFromParticle() != nullptr)
3073  return;
3074 
3075  //Checks if the particle is interacting with the current wall
3078  if (i!=nullptr) {
3079  //...calculates the forces between the two objects...
3080  i->computeForce();
3081 
3082  //...and applies them to each of the two objects (wall and particle).
3083  pI->addForce(i->getForce());
3084  w->addForce(-i->getForce());
3085 
3086  //If the rotation flag is on, also applies the relevant torques
3087  //(getRotation() returns a boolean).
3088  if (getRotation()) // getRotation() returns a boolean.
3089  {
3090  pI->addTorque(i->getTorque() - Vec3D::cross(pI->getPosition() - i->getContactPoint(), i->getForce()));
3092  w->addTorque(-i->getTorque() + Vec3D::cross(w->getPosition() - i->getContactPoint(), i->getForce()));
3093  }
3094  }
3095 }
3096 
3109 {
3110  //cycling through all particles, p, in the particleHandler
3111  //for_each(particleHandler.begin(), particleHandler.end(), [this](BaseParticle* p)
3112  //for (BaseParticle* p : particleHandler) {
3113 
3114  #pragma omp parallel for num_threads(getNumberOfOMPThreads()) //schedule(dynamic)
3115  for (int k = 0; k < particleHandler.getNumberOfObjects(); ++k) {
3117 #ifdef MERCURY_USE_MPI
3118  //MPI particles are not integrated, they are purely ghost particles and get their new velocity and position from an MPI update
3119  if (!(p->isMPIParticle() || p->isPeriodicGhostParticle()))
3120  {
3122  }
3123 #else
3124  //using the particle p's internal "integrateBeforeForceComputation" function
3125  //to update the relevant parameters concerning the particle's position and motion
3127 #endif
3128  }
3129  //});
3130  //cycling through all walls, w, in the wallHandler
3131  //for_each(wallHandler.begin(), wallHandler.end(), [this](BaseWall* w)
3132  //for (BaseWall* w : wallHandler) {
3133  #pragma omp parallel for num_threads(getNumberOfOMPThreads()) //schedule(dynamic)
3134  for (int k = 0; k < wallHandler.getNumberOfObjects(); k++) {
3135  BaseWall *w = wallHandler.getObject(k);
3136  //using the wall's internal "integrateBeforeForceComputation" function
3137  //to update the relevant parameters concerning its position and motion
3139  }
3140  //});
3141 }
3142 
3153 {
3154 
3155  //Cycling over all boundaries within the system...
3156  for (BaseBoundary* b : boundaryHandler)
3157  {
3158  //check all boundaries...
3159  b->checkBoundaryAfterParticlesMove(particleHandler);
3160 
3161 
3162 #ifdef MERCURY_USE_MPI
3163  //When ghost particles are deleted by deletion boundaries they need to be removed
3164  //from their communication lists to avoid segfaults
3165  if (NUMBER_OF_PROCESSORS > 1)
3166  {
3167  //Flush deleted particles from mpi communication zones
3168  getCurrentDomain()->flushParticles(b->getParticlesToBeDeleted());
3170  periodicBoundaryHandler.flushParticles(b->getParticlesToBeDeleted());
3172  }
3173 
3174  //Delete particles that were in communication zone
3175  for (auto p_it = b->getParticlesToBeDeleted().begin(); p_it != b->getParticlesToBeDeleted().end(); p_it++)
3176  {
3177  particleHandler.removeGhostObject((*p_it)->getIndex());
3178  }
3179 #endif
3180  }
3181 }
3182 
3195 {
3196  //cycling through all particles, p, in the particleHandler
3197  //for_each(particleHandler.begin(), particleHandler.end(), [this](BaseParticle* p){
3198  #pragma omp parallel for num_threads(getNumberOfOMPThreads()) //schedule(dynamic)
3199  for (int k = 0; k < particleHandler.getNumberOfObjects(); ++k) {
3201 #ifdef MERCURY_USE_MPI
3202  //MPI particles do not require integration - they are updated by the communication step
3203  if (!(p->isMPIParticle() || p->isPeriodicGhostParticle()))
3204  {
3206  }
3207 #else
3208  //using the particle p's internal "integrateAfterForceComputation" function
3209  //to update the relevant parameters concerning the particle's position and motion
3211 #endif
3212  }
3213  //});
3214  //cycling through all walls, w, in the wallHandler
3215  //for_each(wallHandler.begin(), wallHandler.end(), [this](BaseWall* w){
3216  #pragma omp parallel for num_threads(getNumberOfOMPThreads()) //schedule(dynamic)
3217  for (int k = 0; k < wallHandler.getNumberOfObjects(); k++) {
3218  BaseWall *w = wallHandler.getObject(k);
3219  //using the wall's internal "integrateAfterForceComputation" function
3220  //to update the relevant parameters concerning its position and motion
3222  }
3223  //});
3224 }
3225 
3229 //void DPMBase::statisticsFromRestartData(const char *name)
3230 //{
3231 // ///todo{Check this whole function}
3232 // //This function loads all MD data
3233 // readRestartFile();
3234 //
3235 // //This creates the file statistics will be saved to
3236 // std::stringstream ss("");
3237 // ss << name << ".stat";
3238 // statFile.setName(ss.str());
3239 // statFile.setOpenMode(std::fstream::out);
3240 // statFile.open();
3241 //
3242 // // Sets up the initial conditions for the simulation
3243 // // setupInitialConditions();
3244 // // Setup the previous position arrays and mass of each particle.
3245 // computeParticleMasses();
3246 // // Other routines required to jump-start the simulation
3247 // actionsBeforeTimeLoop();
3248 // initialiseStatistics();
3249 // hGridActionsBeforeTimeLoop();
3250 // writeEneHeader(eneFile.getFstream());
3251 //
3252 // while (readDataFile(dataFile.getName().c_str()))
3253 // {
3254 // hGridActionsBeforeTimeLoop();
3255 // actionsBeforeTimeStep();
3256 // checkAndDuplicatePeriodicParticles();
3257 // hGridActionsBeforeTimeStep();
3262 // computeAllForces();
3263 // removeDuplicatePeriodicParticles();
3264 // actionsAfterTimeStep();
3265 // writeEneTimeStep(eneFile.getFstream());
3266 // std::cout << std::setprecision(6) << getTime() << std::endl;
3267 // }
3268 //
3269 // dataFile.close();
3270 // statFile.close();
3271 //}
3278 {
3279  //Resetting all forces on both particles and walls to zero
3280  #pragma omp parallel num_threads(getNumberOfOMPThreads())
3281  {
3282  #pragma omp for
3283  for (int k = 0; k < particleHandler.getNumberOfObjects(); ++k) {
3285  }
3286  #pragma omp for
3287  for (int k = 0; k < wallHandler.getNumberOfObjects(); k++) {
3289  }
3290  }
3291  logger(DEBUG,"All forces set to zero");
3292 
3293  // for omp simulations, reset the newObjects_ variable (used for reduction)
3295 
3296  // compute all internal and external forces; for omp simulations, this can be done in parallel
3297  #pragma omp parallel num_threads(getNumberOfOMPThreads())
3298  {
3299  //logger(INFO, "Number of omp threads = %", getNumberOfOMPThreads());
3301  #pragma omp for schedule(dynamic)
3302  for (int k = 0; k < particleHandler.getNumberOfObjects(); ++k) {
3304  //computing both internal forces (e.g. due to collisions)
3305  //and external forces (e.g. gravity)
3306  //(compute internal forces compares the current particle p
3307  //with all others in the handler!)
3309  // body forces
3311  }
3312 
3313  // wall-forces
3314  #pragma omp for schedule(dynamic)
3315  for (int k = 0; k < wallHandler.getNumberOfObjects(); k++) {
3316  BaseWall *w = wallHandler.getObject(k);
3317  computeWallForces(w);
3318  }
3319 
3320  }
3321 
3322 #ifdef CONTACT_LIST_HGRID
3323  PossibleContact* Curr=possibleContactList.getFirstPossibleContact();
3324  while(Curr)
3325  {
3326  computeInternalForces(Curr->getP1(),Curr->getP2());
3327  Curr=Curr->getNext();
3328  }
3329 #endif
3330 
3331  // for omp simulations, sum up all forces and add all newObjects_ (needed since both are using reduction)
3332  #ifdef MERCURY_USE_OMP
3333  if (getNumberOfOMPThreads()>1) {
3335  }
3336  //Resetting all forces on both particles and walls to zero
3337  #pragma omp parallel num_threads(getNumberOfOMPThreads())
3338  {
3339  #pragma omp for
3340  for (int k = 0; k < particleHandler.getNumberOfObjects(); k++) {
3342  }
3343  #pragma omp for
3344  for (int k = 0; k < wallHandler.getNumberOfObjects(); k++) {
3346  } //end reset forces loop
3347  }
3348  #endif
3349 
3350  //end outer loop over contacts.
3351 }
3352 
3360 {
3361  for (auto it = particleHandler.begin(); (*it) != i; ++it)
3362  {
3363  computeInternalForce(i, *it);
3364  }
3365 }
3366 
3367 
3375 void DPMBase::write(std::ostream& os, bool writeAllParticles) const
3376 {
3377  os << "MercuryDPM " << getVersion();
3378  //which outputs basic information regarding the various files (.data, .fstat etc. etc.)
3379  //only writes the run number if it is different from 0
3380  if (runNumber_ != 0)
3381  os << " runNumber " << runNumber_;
3382  os << " name " << name_;
3383  os << " revision " << getSVNRevision();
3384  os << " repository " << getSVNURL() << '\n';
3385  os << "dataFile " << dataFile << '\n';
3386  os << "fStatFile " << fStatFile << '\n';
3387  os << "eneFile " << eneFile << '\n';
3388  os << "restartFile " << restartFile << '\n';
3389  os << "statFile " << statFile << '\n';
3390  os << "interactionFile " << interactionFile << '\n';
3391  //Outputs the "domain" corresponding to the system for
3392  //use with XBalls, as well as other information regarding the system as a whole
3393  os << "xMin " << getXMin()
3394  << " xMax " << getXMax()
3395  << " yMin " << getYMin()
3396  << " yMax " << getYMax()
3397  << " zMin " << getZMin()
3398  << " zMax " << getZMax() << '\n'
3399  << "timeStep " << getTimeStep()
3400  << " time " << getTime()
3401  << " ntimeSteps " << numberOfTimeSteps_
3402  << " timeMax " << getTimeMax() << '\n'
3403  << "systemDimensions " << getSystemDimensions()
3404  << " particleDimensions " << getParticleDimensions()
3405  << " gravity " << getGravity();
3406  os << " writeVTK " << writeParticlesVTK_
3407  << " " << writeWallsVTK_
3408  << " " << interactionHandler.getWriteVTK()
3409  << " " << (vtkWriter_?vtkWriter_->getFileCounter():0)
3410  << " " << wallVTKWriter_.getFileCounter()
3412  << " " << boundaryVTKWriter_.getFileCounter();
3413  os << " random ";
3414  random.write(os);
3415 #ifdef MERCURY_USE_OMP
3416  //Write number of OMP threads
3417  if(getNumberOfOMPThreads() > 1) {
3418  os << " numberOfOMPThreads " << getNumberOfOMPThreads();
3419  }
3420 #endif
3421 #ifdef MERCURY_USE_MPI
3422  //Check if we are dealing with multiple cores
3423  if (NUMBER_OF_PROCESSORS > 1 )
3424  {
3425  os << " numberOfProcessors " << NUMBER_OF_PROCESSORS
3427  }
3428 #endif
3429 
3430  //only write xBallsArguments if they are nonzero
3432  os << " xBallsArguments " << getXBallsAdditionalArguments();
3433  os << '\n';
3434  //writes all species (including mixed species) to an output stream
3435 
3436  speciesHandler.write(os);
3437 
3438  //outputs the number of walls in the system
3439  os << "Walls " << wallHandler.getNumberOfObjects() << std::endl;
3440  if (writeAllParticles || wallHandler.getSize() < 9) {
3441  for (BaseWall* w : wallHandler)
3442  os << (*w) << std::endl;
3443  } else {
3444  for (int i=0; i<2; ++i)
3445  os << *wallHandler.getObject(i) << std::endl;
3446  os << "...\n";
3447  }
3448 
3449  //outputs the number of boundaries in the system
3450  os << "Boundaries " << boundaryHandler.getNumberOfObjects() << std::endl;
3451  if (writeAllParticles || boundaryHandler.getSize() < 9) {
3452  for (BaseBoundary* b : boundaryHandler)
3453  os << (*b) << std::endl;
3454  } else {
3455  for (int i=0; i<2; ++i)
3456  os << *boundaryHandler.getObject(i) << std::endl;
3457  os << "...\n";
3458  }
3459 
3460  int nToWrite = 4; // \todo JMFT: Don't hardcode this here, but put it in the argument
3461 
3462  if (writeAllParticles || particleHandler.getSize() < nToWrite)
3463  {
3464  //if the "writeAllParticles" bool == true, or there are fewer than 4 particles
3465  //calls the particleHandler version of the "write" function and also
3466  //outputs to file all relevant particle information for all particles in the system
3467  particleHandler.write(os);
3468  }
3469  else
3470  {
3471  //otherwise, only prints out limited information
3472  os << "Particles " << particleHandler.getSize() << '\n';
3473  for (unsigned int i = 0; i < nToWrite; i++)
3474  os << *(particleHandler.getObject(i)) << '\n';
3475  os << "..." << '\n';
3476  }
3477  // Similarly, print out interaction details (all of them, or up to nToWrite of them)
3478  if (writeAllParticles || interactionHandler.getNumberOfObjects() < nToWrite)
3479  {
3481  }
3482  else
3483  {
3484  os << "Interactions " << interactionHandler.getNumberOfObjects() << '\n';
3485  for (unsigned int i = 0; i < nToWrite; i++)
3486  os << *(interactionHandler.getObject(i)) << '\n';
3487  os << "..." << '\n';
3488  }
3489 }
3490 
3496 void DPMBase::read(std::istream& is, ReadOptions opt)
3497 {
3498 #ifdef MERCURY_USE_MPI
3499  int previousNumberOfProcessors;
3500 #endif
3501  //Declares...
3502  std::string dummy;
3503  //...and reads in a dummy variable from the start of the stream "is"
3504  is >> dummy;
3505  //compare the string read in to the phrase "restart_version" to see if the stream corresponds
3506  //to a restart file (all restart files begin with this phrase)
3507  //if both strings match, strcmp(dummy.c_str(), "restart_version") returns 0 (here read as "false")
3508  if (dummy != "restart_version" && dummy != "MercuryDPM")
3509  {
3510  //If the strings do not match, if statement is fulfilled and the error logged
3511  //Note: only very old files did not have a restart_version
3512  logger(FATAL, "Error in DPMBase::read(is): this is not a valid restart file");
3513  }
3514  else
3515  {
3516  //reads in the restart version (earlier versions of Mercury possess different file formats!)
3517  is >> restartVersion_;
3518  //checking which version the current data file corresponds to, and reads the data in
3519  //accordingly
3520  if (restartVersion_ == "1.0" || restartVersion_ == "0.14")
3521  {
3522  //reads in and saves the relevant values from the data file to the current instance of DPMBase
3523  std::stringstream line;
3524 
3525  // Store path (if restart file is nonlocal)
3526  auto slash = restartFile.getName().rfind('/');
3527  std::string path;
3528  if (slash != std::string::npos)
3529  {
3530  path = restartFile.getName().substr(0, slash + 1);
3531  }
3532  if (!path.empty())
3533  {
3534  logger(INFO, "Adding path information (%) to file names", path);
3535  }
3536 
3537  //line 1
3539  //discards the whitespace (ws) at the start of the stream
3540  line >> std::ws;
3541  //uses the "peek" function to access the stream's first
3542  //non-whitespace character, and check if it is an "r"
3543  if (line.peek() == 'r')
3544  //if so, reads in the current run number
3545  line >> dummy >> runNumber_;
3546  //In either case, then uses the "Files" version of the read function
3547  //to read in the rest of the relevant information.
3548  line >> dummy >> name_;
3549  setName(name_);
3550 
3551  //Read line 2-7 (definition of i/o files)
3553  line >> dummy >> dataFile;
3555  line >> dummy >> fStatFile;
3557  line >> dummy >> eneFile;
3559  line >> dummy >> restartFile;
3561  line >> dummy >> statFile;
3562 
3563  // Add the file path from the restart file to the file names
3564  dataFile.setName(path + dataFile.getName());
3565  fStatFile.setName(path + fStatFile.getName());
3566  eneFile.setName(path + eneFile.getName());
3567  restartFile.setName(path + restartFile.getName());
3568  statFile.setName(path + statFile.getName());
3569 
3570  // Get current position
3571  //check if the next line starts with 'interactionFile'; otherwise, skip interaction
3572  if (helpers::compare(is, "interactionFile"))
3573  {
3575  line >> interactionFile;
3576  interactionFile.setName(path + interactionFile.getName());
3577  }
3578 
3580  line >> dummy >> min_.x()
3581  >> dummy >> max_.x()
3582  >> dummy >> min_.y()
3583  >> dummy >> max_.y()
3584  >> dummy >> min_.z()
3585  >> dummy >> max_.z();
3586 
3588  line >> dummy >> timeStep_
3589  >> dummy >> time_
3590  >> dummy >> numberOfTimeSteps_
3591  >> dummy >> timeMax_;
3592 
3594  line >> dummy >> systemDimensions_
3595  >> dummy >> particleDimensions_
3596  >> dummy >> gravity_;
3597 
3598  line >> dummy;
3599  if (!dummy.compare("writeVTK"))
3600  {
3601  FileType writeInteractionsVTK = FileType::NO_FILE;
3602  unsigned particlesCounter, wallCounter, interactionCounter;
3603  bool writeBoundaryVTK;
3604  line >> writeParticlesVTK_ >> writeWallsVTK_ >> writeInteractionsVTK >> writeBoundaryVTK >> particlesCounter >> wallCounter >> interactionCounter;
3605  line.clear();//because the number of arguments in writeVTK has changed
3606  line >> dummy;
3609  interactionHandler.setWriteVTK(writeInteractionsVTK);
3610  boundaryHandler.setWriteVTK(writeBoundaryVTK);
3611  vtkWriter_->setFileCounter(particlesCounter);
3612  wallVTKWriter_.setFileCounter(particlesCounter);
3613  interactionVTKWriter_.setFileCounter(particlesCounter);
3614  boundaryVTKWriter_.setFileCounter(particlesCounter);
3615  }
3616  if (!dummy.compare("random"))
3617  {
3618  random.read(line);
3619  line >> dummy;
3620  }
3621 
3622 #ifdef MERCURY_USE_OMP
3623  //Read the number of OMP threads
3624  if (!dummy.compare("numberOfOMPThreads")) {
3625  int numberOfOMPThreads;
3626  line >> numberOfOMPThreads;
3627  setNumberOfOMPThreads(numberOfOMPThreads);
3628  //logger(INFO," Check the number of OMP threads = % ", getNumberOfOMPThreads());
3629  }
3630 #endif
3631 #ifdef MERCURY_USE_MPI
3632  if (!dummy.compare("numberOfProcessors"))
3633  {
3634  line >> previousNumberOfProcessors
3635  >> dummy >> numberOfDomains_[Direction::XAXIS]
3638  }
3639  else
3640  {
3641  logger(INFO,"Reading a serial restart file");
3642  //numberOfDomains_ = {1,1,1};
3643  }
3644 #endif
3645  if (!dummy.compare("xBallsArguments")) {
3647  setXBallsAdditionalArguments(line.str());
3648  }
3649 
3650  speciesHandler.read(is);
3651 
3652 #ifdef MERCURY_USE_MPI
3653  //Initialise MPI structures and perform domain decomposition
3654  decompose();
3655 #endif
3656 
3657  //reading in the various relevant handlers
3658  unsigned int N;
3659  is >> dummy >> N;
3660  if (dummy.compare("Walls"))
3661  logger(ERROR, "DPMBase::read(is): Error during restart: 'Walls' argument could not be found.");
3662  wallHandler.clear();
3665  for (unsigned int i = 0; i < N; i++)
3666  {
3669  }
3670 
3671  is >> dummy >> N;
3674  if (dummy.compare("Boundaries"))
3675  logger(ERROR, "DPMBase::read(is): Error during restart: 'Boundaries' argument could not be found.");
3677  for (unsigned int i = 0; i < N; i++)
3678  {
3681  }
3682 
3684 
3685  is >> dummy >> N;
3686  is.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
3687  //display a message if a large amount o fparticles is read
3688  if (N>2.5e5) logger(INFO, "Reading % particles (may take a while)",N);
3689  logger.assert_always(dummy.compare("Particles")==0, "DPMBase::read(is): Error during restart: 'Particles' argument could not be found. %",dummy);
3692  for (unsigned int i = 0; i < N; i++)
3693  {
3694  //ParticleHandler::readAndCreateObject reads line-by-line
3696  //skip the remaining data in line
3697  is.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
3699  //particleHandler.getLastObject()->computeMass();
3700  }
3701 #ifdef MERCURY_USE_MPI
3702  //Interaction distances of the domainHandler and periodicBoundaryHandler need to be set
3703  Mdouble interactionRadius = particleHandler.getLargestInteractionRadius();
3704  domainHandler.setInteractionDistance(2.0*interactionRadius);
3705  periodicBoundaryHandler.setInteractionDistance(2.0*interactionRadius);
3706 
3707  if (NUMBER_OF_PROCESSORS > 1)
3708  {
3709  //Create ghost particles
3712  }
3713 #endif
3714  //Add interactions to particles and ghost particles
3715  if (opt==ReadOptions::ReadNoInteractions) return;
3717  }
3718  //reading in for older versions of the Mercury restart file.
3719  else if (!restartVersion_.compare("3"))
3720  {
3721  logger(INFO, "DPMBase::read(is): restarting from an old restart file (restart_version %).",
3722  restartVersion_);
3723  readOld(is);
3724  }
3725  //returning an error if there is no restart file to read in due to the use of outdated files.
3726  else
3727  {
3728  //only very old files did not have a restart_version
3729  logger(FATAL,
3730  "Error in DPMBase::read(is): restart_version % cannot be read; use an older version of Mercury to upgrade the file",
3731  restartVersion_);
3732  }
3733  }
3734 }
3735 
3739 void DPMBase::readOld(std::istream& is)
3740 {
3741  std::string dummy;
3742  is >> dummy >> dummy;
3743  setName(dummy);
3744 
3745  unsigned int saveCountData, saveCountEne, saveCountStat, saveCountFStat;
3746  unsigned int fileTypeFstat, fileTypeData, fileTypeEne, fileTypeRestart;
3747  is >> dummy >> min_.x()
3748  >> dummy >> max_.x()
3749  >> dummy >> min_.y()
3750  >> dummy >> max_.y()
3751  >> dummy >> min_.z()
3752  >> dummy >> max_.z()
3753  >> dummy >> timeStep_
3754  >> dummy >> time_
3755  >> dummy >> timeMax_
3756  >> dummy >> saveCountData
3757  >> dummy >> saveCountEne
3758  >> dummy >> saveCountStat
3759  >> dummy >> saveCountFStat
3760  >> dummy >> systemDimensions_
3761  >> dummy >> gravity_
3762  >> dummy >> fileTypeFstat
3763  >> dummy >> fileTypeData
3764  >> dummy >> fileTypeEne;
3765  dataFile.setSaveCount(saveCountData);
3766  eneFile.setSaveCount(saveCountEne);
3767  statFile.setSaveCount(saveCountStat);
3768  fStatFile.setSaveCount(saveCountFStat);
3769 
3770  fStatFile.setFileType(static_cast<FileType>(fileTypeFstat));
3771  dataFile.setFileType(static_cast<FileType>(fileTypeData));
3772  eneFile.setFileType(static_cast<FileType>(fileTypeEne));
3773 
3774  //this is optional to allow restart files with and without restartFile.getFileType()
3775  is >> dummy;
3776  if (!strcmp(dummy.c_str(), "options_restart"))
3777  {
3778  is >> fileTypeRestart;
3779  restartFile.setFileType(static_cast<FileType>(fileTypeRestart));
3780  }
3781 
3782  speciesHandler.read(is);
3783  wallHandler.read(is);
3784  boundaryHandler.read(is);
3785  particleHandler.read(is);
3786  setRestarted(true);
3787 }
3788 
3797 {
3798  //check if name is set
3799  logger.assert_always(getName() != "",
3800  "File name not set: use setName()");
3801  //check if time step is set
3802  logger.assert_always(getTimeStep() != 0,
3803  "Time step undefined: use setTimeStep()");
3804  //check if domain is set
3805  logger.assert_always(getXMax() > getXMin(),
3806  "Domain size not set: use setXMin() and setXMax()");
3807  logger.assert_always(getYMax() > getYMin(),
3808  "Domain size not set: use setYMin() and setYMax()");
3809  logger.assert_always(systemDimensions_ == 3 ? (getZMax() > getZMin()) : (getZMax() >= getZMin()),
3810  "Domain size not set: use setZMin() and setZMax()", systemDimensions_);
3811 
3812  //check for species parameters
3813  logger.assert_always(speciesHandler.getNumberOfObjects() > 0,
3814  "No species defined: use speciesHandler.copyAndAddObject()");
3815  for (BaseParticle* p : particleHandler)
3816  {
3817  logger.assert_always(p->getSpecies() != nullptr, "particle % has no species", p->getId());
3818  }
3819  for (BaseWall* w : wallHandler)
3820  {
3821  logger.assert_always(w->getSpecies() != nullptr, "% with index % has no species", w->getName(), w->getId());
3822  }
3823 }
3824 
3826 {
3828  writeOutputFiles();
3829 }
3830 
3847 {
3848  //writing fstat data if .saveCurrentTimeStep(numberOfTimeSteps_) is true
3850  {
3851  writeFStatFile();
3852  }
3853 
3854  //writing .ene data if .saveCurrentTimeStep(numberOfTimeSteps_) is true
3856  {
3857  writeEneFile();
3858  }
3859  //writing .data data if .saveCurrentTimeStep(numberOfTimeSteps_) is true
3861  {
3863  if (getRestarted() || dataFile.getCounter() == 0)
3865  writeDataFile();
3866  } else {
3868  }
3869  printTime();
3870  writeVTKFiles();
3871  }
3872  cgHandler.evaluate();
3873 
3874 
3875  //write restart file last, otherwise the output counters are wrong
3877  {
3878  writeRestartFile();
3879  }
3880 }
3881 
3886 {
3887 #ifdef MERCURY_USE_MPI
3888 
3889  //If running in parallel build, but just running with one core - no domain decomposition required
3890  int numberOfRequiredProcessors = numberOfDomains_[Direction::XAXIS]*
3893  if (NUMBER_OF_PROCESSORS != numberOfRequiredProcessors)
3894  {
3895  logger(ERROR,"The domain decompositions expects % processors, but only % are requested.\n"
3896  "Either run your process using \"mpirun -np % [executable]\", "
3897  "or change the domain decomposition to e.g. setNumberOfDomains({%,1,1}).", numberOfRequiredProcessors, NUMBER_OF_PROCESSORS, numberOfRequiredProcessors, NUMBER_OF_PROCESSORS);
3898  }
3899 
3900  if (NUMBER_OF_PROCESSORS == 1) {return;}
3901 
3902  //Check if the simulation domain has been set
3903  logger.assert_always(getXMax() - getXMin() > 0,"Please set your simulation domain (setXMax(),setXmin()) before calling solve()");
3904  logger.assert_always(getYMax() - getYMin() > 0,"Please set your simulation domain (setYMax(),setYmin()) before calling solve()");
3905  logger.assert_always(getZMax() - getZMin() > 0,"Please set your simulation domain (setZMax(),setZmin()) before calling solve()");
3906 
3907  //Grab simulation domains
3908  std::vector<Mdouble> simulationMin{getXMin(), getYMin(), getZMin()};
3909  std::vector<Mdouble> simulationMax{getXMax(), getYMax(), getZMax()};
3910 
3911  //Check if the user input decomposition is correct
3912  logger.assert_always(numberOfDomains_[Direction::XAXIS] > 0,"Number of domain in x-direction incorrect: %",numberOfDomains_[Direction::XAXIS]);
3913  logger.assert_always(numberOfDomains_[Direction::YAXIS] > 0,"Number of domain in y-direction incorrect: %",numberOfDomains_[Direction::YAXIS]);
3914  logger.assert_always(numberOfDomains_[Direction::ZAXIS] > 0,"Number of domain in z-direction incorrect: %",numberOfDomains_[Direction::ZAXIS]);
3915 
3916  //Open domain decomposition, closed is not implemented
3917  bool open = true;
3918 
3919  //Check if the number of domains is equal to the number of processors
3920  logger.assert_always(numberOfDomains_[Direction::XAXIS]*numberOfDomains_[Direction::YAXIS]*numberOfDomains_[Direction::ZAXIS] == NUMBER_OF_PROCESSORS,
3921  "Number of Processors is not equal to number of domains. Processors %, domains, %",
3923  numberOfDomains_[Direction::XAXIS]*numberOfDomains_[Direction::YAXIS]*numberOfDomains_[Direction::ZAXIS]);
3924 
3925  //Create all processor domains
3926 
3927  domainHandler.setDPMBase(this);
3928  domainHandler.createMesh(simulationMin, simulationMax, numberOfDomains_, open);
3929  logger(VERBOSE,"Number of domains: % | Number of processors: %",domainHandler.getNumberOfObjects(), NUMBER_OF_PROCESSORS);
3930  //logger.assert_always(domainHandler.getNumberOfObjects() == numberOfProcessors, "Invalid decomposition: Number of domains and processors are different");
3931 
3932  //Tell the current processor to which domain it belongs
3933  for (Domain* domain : domainHandler)
3934  {
3935  if (domain->getRank() == PROCESSOR_ID)
3936  {
3937  logger(VERBOSE,"processor: %, domain index: %",PROCESSOR_ID, domain->getIndex());
3938  domainHandler.setCurrentDomainIndex(domain->getIndex());
3939  }
3940  }
3941 
3942  //Define the mpi transfer types, which requires a definition of the species already
3944  logger.assert_always(speciesHandler.getNumberOfObjects() > 0, "Please create a particle species before calling solve()");
3946 
3947  //Make sure all processors are done with decomposition before proceeding
3948  logger(VERBOSE,"processor %: #real particles: %, #total particles: %", PROCESSOR_ID, particleHandler.getNumberOfRealObjects(), particleHandler.getSize());
3950 #endif
3951 }
3952 
3953 
3969 {
3970  logger(DEBUG, "Entered solve");
3971 #ifdef CONTACT_LIST_HGRID
3972  logger(INFO,"Using CONTACT_LIST_HGRID");
3973 #endif
3974 
3979  if (!getRestarted())
3980  {
3981  // If the simulation is "new" (i.e. not restarted):
3982  // - set time, nTimeSteps to zero
3983  // - reset the file counter etc.
3984  // - decompose the domain based on XMin, XMax, ....
3985  // - run user-defined setupInitialConditions
3986  numberOfTimeSteps_ = 0;
3987  setTime(0.0);
3988  resetFileCounter();
3989  decompose();
3990  //\todo tw there was a function combining the next two lines, why is it back to the old version?
3991  //setLastSavedTimeStep(NEVER); //reset the counter
3992  //this is to ensure that the interaction time stamps agree with the resetting of the time value
3993  for (auto& i : interactionHandler)
3994  i->setTimeStamp(0);
3996  logger(DEBUG, "Have created the particles initial conditions");
3997  }
3998  else
3999  {
4000  // If the simulation is "restarted" (i.e. not restarted):
4001  // - run user-defined actionsOnRestart
4002  actionsOnRestart();
4003  }
4004 
4005  // Check that the code has been correctly set up,
4006  // i.e. system dimensions, particles and time steps are sensibly implemented
4007  checkSettings();
4008 
4009  // If the simulation is "new" and the runNumber is used, append the run number to the problem name
4010  if (getRunNumber() > 0 && !getRestarted())
4011  {
4012  std::stringstream name;
4013  name << getName() << "." << getRunNumber();
4014  setName(name.str());
4015  }
4016 
4017  //If append is true, files are appended, not overwritten
4018  if (getAppend())
4019  {
4020  setOpenMode(std::fstream::out | std::fstream::app);
4021  //Restart files should always be overwritten.
4022  restartFile.setOpenMode(std::fstream::out);
4023  }
4024  else
4025  {
4026  setOpenMode(std::fstream::out);
4027  }
4028 
4029  //sets the hgrid, writes headers to the .stat output file
4031 
4032  if (getInteractionFile().getFileType() == FileType::ONE_FILE)
4033  {
4034  logger(WARN, "Warning: interaction file will take up a lot of disk space!");
4036  }
4037 
4038  // Sets the mass of all particle.
4042 
4043  // Other initialisations
4044  //max_radius = getLargestParticle()->getRadius();
4045 
4049 
4050  // Performs a first force computation
4052 
4053 #ifdef MERCURY_USE_MPI
4054  if (NUMBER_OF_PROCESSORS > 1)
4055  {
4056  //Find new mpi particles
4058  //Periodic particles in parallel
4060  }
4061 #endif
4062 
4064  computeAllForces();
4067  logger(DEBUG, "Have computed the initial values for the forces ");
4068 
4069  // This is the main loop over advancing time
4070  while (getTime() < getTimeMax() && continueSolve())
4071  {
4073  }
4074  //force writing of the last time step
4076 
4077  //end loop over interaction count
4079 
4080  //To make sure getTime gets the correct time for outputting statistics
4081  finishStatistics();
4082 
4083  closeFiles();
4084 }
4085 
4091 {
4092  logger(DEBUG, "starting computeOneTimeStep()");
4093 
4094  logger(DEBUG, "about to call writeOutputFiles()");
4095  writeOutputFiles(); //everything is written at the beginning of the time step!
4096 
4097  logger(DEBUG, "about to call hGridActionsBeforeIntegration()");
4099 
4100  //Computes the half-time step velocity and full time step position and updates the particles accordingly
4101  logger(DEBUG, "about to call integrateBeforeForceComputation()");
4102 
4104  //New positions require the MPI and parallel periodic boundaries to do things
4105  logger(DEBUG, "about to call performGhostParticleUpdate()");
4107 
4111 
4112  logger(DEBUG, "about to call checkInteractionWithBoundaries()");
4113  checkInteractionWithBoundaries(); // INSERTION boundaries handled
4114 
4115  logger(DEBUG, "about to call hGridActionsAfterIntegration()");
4117 
4118  // Compute forces
4120  // INSERTION/DELETION boundary flag change
4121  for (BaseBoundary* b : boundaryHandler)
4122  {
4123  b->checkBoundaryBeforeTimeStep(this);
4124  }
4125 
4126  logger(DEBUG, "about to call actionsBeforeTimeStep()");
4128 
4129  logger(DEBUG, "about to call checkAndDuplicatePeriodicParticles()");
4131 
4132  logger(DEBUG, "about to call hGridActionsBeforeTimeStep()");
4134 
4135  //Creates and updates interactions and computes forces based on these
4136  logger(DEBUG, "about to call computeAllForces()");
4137  computeAllForces();
4138 
4139  logger(DEBUG, "about to call removeDuplicatePeriodicParticles()");
4141 
4142  logger(DEBUG, "about to call actionsAfterTimeStep()");
4144 
4145  //Computes new velocities and updates the particles accordingly
4146  logger(DEBUG, "about to call integrateAfterForceComputation()");
4148 
4149  //erase interactions that have not been used during the last time step
4150  //logger(DEBUG, "about to call interactionHandler.eraseOldInteractions(getNumberOfTimeSteps())");
4152  logger(DEBUG, "about to call interactionHandler.actionsAfterTimeStep()");
4155 
4156  time_ += timeStep_;
4158 
4159  logger(DEBUG, "finished computeOneTimeStep()");
4160 }
4161 
4174 bool DPMBase::readArguments(int argc, char* argv[])
4175 {
4176  bool isRead = true;
4177  // Cycles over every second element. (Most flags will contain both name labels and actual data.
4178  // Those that don't will have to do i--; some examples in readNextArgument.)
4179  for (int i = 1; i < argc; i += 2)
4180  {
4181  std::cout << "interpreting input argument " << argv[i];
4182  for (int j = i + 1; j < argc; j++)
4183  {
4184  //looks for the next string that starts with a minus sign
4185  //i.e. the next flag, as each flag may take 0, 1 , 2, 3... arguments
4186  //and we need to make sure all are read in!
4187  if (argv[j][0] == '-')
4188  break;
4189  std::cout << " " << argv[j];
4190  }
4191  std::cout << std::endl;
4192  //if "isRead"is true and "readNextArgument" is also true...
4193  //(i.e. checking if any argument is false)
4194  isRead &= readNextArgument(i, argc, argv);
4195 
4196  // If the read was unsuccessful, raise an error and quit. (JMFT: Used to just be a warning.)
4197  if (!isRead)
4198  {
4199  logger(ERROR, "Warning: not all arguments read correctly!");
4200  }
4201  }
4202  return isRead;
4203 }
4204 
4206 {
4207  //logger(INFO,"ID %",PROCESSOR_ID);
4208  //if (PROCESSOR_ID!=0) return;
4209 
4210  logger(INFO,"Removing old files named %.*",getName());
4211  std::ostringstream filename;
4212 
4213  // add processor id to file extension for mpi jobs
4214  std::string p = (NUMBER_OF_PROCESSORS > 1)?std::to_string(PROCESSOR_ID):"";
4215  // all the file extensions that should be deleted
4216  std::vector<std::string> ext{".restart"+p, ".stat"+p, ".fstat"+p, ".data"+p, ".ene"+p, ".xballs"};
4217  for (const auto& j : ext)
4218  {
4219  // remove files with given extension for FileType::ONE_FILE
4220  filename.str("");
4221  filename << getName() << j;
4222  if (!remove(filename.str().c_str()))
4223  {
4224  logger(INFO," File % successfully deleted",filename.str());
4225  }
4226  // remove files with given extension for FileType::MULTIPLE_FILES
4227  unsigned k = 0;
4228  filename.str("");
4229  filename << getName() << j << '.' << k;
4230  while (!remove(filename.str().c_str()))
4231  {
4232  if (k<3) logger(INFO," File % successfully deleted",filename.str());
4233  filename.clear();
4234  filename << getName() << j << '.' << ++k;
4235  }
4236  // remove files with given extension for FileType::MULTIPLE_FILES_PADDED
4237  k = 0;
4238  filename.str("");
4239  filename << getName() << j << '.' << to_string_padded(k);
4240  while (!remove(filename.str().c_str()))
4241  {
4242  if (k<3) logger(INFO," File % successfully deleted",filename.str());
4243  filename.clear();
4244  filename << getName() << j << '.' << to_string_padded(++k);
4245  }
4246  }
4247  // remove vtk files
4248  // add processor id to file extension for mpi jobs
4249  std::string q = (NUMBER_OF_PROCESSORS > 1)?("Processor_"+std::to_string(PROCESSOR_ID)+"_"):"";
4250  // all the file extensions that should be deleted
4251  ext = {"Wall_", q+"Particle_", q+"Interaction_"};
4252  for (const auto& j : ext)
4253  {
4254  // remove files with given extension for FileType::ONE_FILE
4255  filename.str("");
4256  filename << getName() << j << ".vtu";
4257  if (!remove(filename.str().c_str()))
4258  {
4259  logger(INFO," File % successfully deleted",filename.str());
4260  }
4261  // remove files with given extension for FileType::MULTIPLE_FILES
4262  unsigned k = 0;
4263  filename.str("");
4264  filename << getName() << j << k << ".vtu";
4265  while (!remove(filename.str().c_str()))
4266  {
4267  if (k<3) logger(INFO," File % successfully deleted",filename.str());
4268  filename.str("");
4269  filename << getName() << j << ++k << ".vtu";
4270  }
4271  //std::cout << " File " << filename.str() << " not found" << std::endl;
4272  }
4273 }
4274 
4299 bool DPMBase::readNextArgument(int& i, int argc, char* argv[])
4300 {
4301  // The argument argv[i] identifies the label of the flag, and subsequent arguments (usually 1)
4302  // contain the content.
4303  //
4304  // For example...
4305  // Checks if the "-name" flag has been passed
4306  // The strcmp returns 0 if "argv[i]" is "-name" (i.e. !strcmp(argv[i], "-name") --> 1)
4307  // In this case, the setName function is run with the relevant input (i.e. the value which
4308  // immediately follows the "-name" flag
4309  if (!strcmp(argv[i], "-name"))
4310  {
4311  setName(argv[i + 1]);
4312  }
4313  // The above process is repeated for all viable flags.
4314  else if (!strcmp(argv[i], "-xmin"))
4315  {
4316  setXMin(atof(argv[i + 1]));
4317  }
4318  else if (!strcmp(argv[i], "-ymin"))
4319  {
4320  setYMin(atof(argv[i + 1]));
4321  }
4322  else if (!strcmp(argv[i], "-zmin"))
4323  {
4324  setZMin(atof(argv[i + 1]));
4325  }
4326  else if (!strcmp(argv[i], "-xmax"))
4327  {
4328  setXMax(atof(argv[i + 1]));
4329  }
4330  else if (!strcmp(argv[i], "-ymax"))
4331  {
4332  setYMax(atof(argv[i + 1]));
4333  }
4334  else if (!strcmp(argv[i], "-zmax"))
4335  {
4336  setZMax(atof(argv[i + 1]));
4337  //} else if (!strcmp(argv[i],"-svn")) {
4338  // std::cout << "svn version " << SVN_VERSION << std::endl;
4339  // i--;
4340  }
4341  else if (!strcmp(argv[i], "-dt"))
4342  {
4343  Mdouble old = getTimeStep();
4344  setTimeStep(atof(argv[i + 1]));
4345  std::cout << " reset dt from " << old << " to " << getTimeStep() << std::endl;
4346  }
4347 // else if (!strcmp(argv[i], "-Hertz"))
4348 // {
4349 // speciesHandler.getObject(0)->setForceType(ForceType::HERTZ);
4350 // i--;
4351 // }
4352  else if (!strcmp(argv[i], "-tmax"))
4353  {
4354  Mdouble old = getTimeMax();
4355  setTimeMax(atof(argv[i + 1]));
4356  std::cout << " reset timeMax from " << old << " to " << getTimeMax() << std::endl;
4357  }
4358  else if (!strcmp(argv[i], "-saveCount"))
4359  {
4360  Mdouble old = dataFile.getSaveCount();
4361  setSaveCount(static_cast<unsigned int>(atoi(argv[i + 1])));
4362  std::cout << " reset saveCount from " << old << " to " << dataFile.getSaveCount() << std::endl;
4363  }
4364  else if (!strcmp(argv[i], "-saveCountData"))
4365  {
4366  dataFile.setSaveCount(static_cast<unsigned int>(atoi(argv[i + 1])));
4367  }
4368  else if (!strcmp(argv[i], "-saveCountFStat"))
4369  {
4370  fStatFile.setSaveCount(static_cast<unsigned int>(atoi(argv[i + 1])));
4371  }
4372  else if (!strcmp(argv[i], "-saveCountStat"))
4373  {
4374  statFile.setSaveCount(static_cast<unsigned int>(atoi(argv[i + 1])));
4375  }
4376  else if (!strcmp(argv[i], "-saveCountEne"))
4377  {
4378  eneFile.setSaveCount(static_cast<unsigned int>(atoi(argv[i + 1])));
4379  }
4380  else if (!strcmp(argv[i], "-saveCountRestart"))
4381  {
4382  restartFile.setSaveCount(static_cast<unsigned int>(atoi(argv[i + 1])));
4383  }
4384  else if (!strcmp(argv[i], "-dim"))
4385  {
4386  setSystemDimensions(static_cast<unsigned int>(atoi(argv[i + 1])));
4387  }
4388  else if (!strcmp(argv[i], "-gravity"))
4389  {
4391  setGravity(Vec3D(atof(argv[i + 1]), atof(argv[i + 2]), atof(argv[i + 3])));
4392  i += 2;
4393  }
4394  else if (!strcmp(argv[i], "-fileType"))
4395  { //uses int input
4396  setFileType(static_cast<FileType>(atoi(argv[i + 1])));
4397  }
4398  else if (!strcmp(argv[i], "-fileTypeFStat"))
4399  { //uses int input
4400  fStatFile.setFileType(static_cast<FileType>(atoi(argv[i + 1])));
4401  }
4402  else if (!strcmp(argv[i], "-fileTypeRestart"))
4403  {
4404  restartFile.setFileType(static_cast<FileType>(atoi(argv[i + 1])));
4405  }
4406  else if (!strcmp(argv[i], "-fileTypeData"))
4407  {
4408  dataFile.setFileType(static_cast<FileType>(atoi(argv[i + 1])));
4409  }
4410  else if (!strcmp(argv[i], "-fileTypeStat"))
4411  {
4412  statFile.setFileType(static_cast<FileType>(atoi(argv[i + 1])));
4413  }
4414  else if (!strcmp(argv[i], "-fileTypeEne"))
4415  {
4416  eneFile.setFileType(static_cast<FileType>(atoi(argv[i + 1])));
4417  }
4418  else if (!strcmp(argv[i], "-auto_number"))
4419  {
4420  autoNumber();
4421  i--;
4422  }
4423 // else if (!strcmp(argv[i], "-number_of_saves"))
4424 // {
4425 // set_number_of_saves_all(atof(argv[i + 1]));
4426 // }
4427  else if (!strcmp(argv[i], "-restart") || !strcmp(argv[i], "-r"))
4428  {
4432  std::string filename;
4433 
4434  //use default filename if no argument is given
4435  if (i + 1 >= argc || argv[i + 1][0] == '-')
4436  {
4437  i--;
4438  filename = getName();
4439  std::cout << getName() << std::endl;
4440  }
4441  else
4442  {
4443  filename = argv[i + 1];
4444  }
4445 
4446  //add ".restart" if necessary
4447  if (filename.find(".restart") == std::string::npos)
4448  {
4449  filename = filename + ".restart";
4450  }
4451 
4452  readRestartFile(filename);
4453  }
4454  else if (!strcmp(argv[i], "-clean") || !strcmp(argv[i], "-c"))
4455  {
4456  std::cout << "Remove old " << getName() << ".* files" << std::endl;
4457  removeOldFiles();
4458  i--;
4459  }
4460  else if (!strcmp(argv[i], "-data"))
4461  {
4462  std::string filename = argv[i + 1];
4463  readDataFile(filename);
4464  }
4465  else if (!strcmp(argv[i], "-readSpeciesFromDataFile"))
4466  {
4467  readSpeciesFromDataFile_ = true;
4468  i--;
4469  logger(INFO, "Last column of data file will be interpreted as species index");
4470  }
4471 // else if (!strcmp(argv[i], "-k"))
4472 // {
4473 // Mdouble old = getSpecies(0)->getStiffness();
4474 // getSpecies(0)->setStiffness(atof(argv[i + 1]));
4475 // std::cout << " reset k from " << old << " to " << getSpecies(0)->getStiffness() << std::endl;
4476 // }
4477 // else if (!strcmp(argv[i], "-dissipation") || !strcmp(argv[i], "-disp"))
4478 // {
4479 // Mdouble old = getSpecies(0)->getDissipation();
4480 // getSpecies(0)->setDissipation(atof(argv[i + 1]));
4481 // std::cout << " reset getDissipation() from " << old << " to " << getSpecies(0)->getDissipation() << std::endl;
4482 // }
4483 // else if (!strcmp(argv[i], "-kt"))
4484 // {
4485 // Mdouble old = getSpecies(0)->getSlidingStiffness();
4486 // getSpecies(0)->setSlidingStiffness(atof(argv[i + 1]));
4487 // std::cout << " reset kt from " << old << " to " << getSpecies(0)->getSlidingStiffness() << std::endl;
4488 // }
4489 // else if (!strcmp(argv[i], "-dispt"))
4490 // {
4491 // Mdouble old = getSpecies(0)->getSlidingDissipation();
4492 // getSpecies(0)->setSlidingDissipation(atof(argv[i + 1]));
4493 // std::cout << " reset dispt from " << old << " to " << getSpecies(0)->getSlidingDissipation() << std::endl;
4494 // }
4495 // else if (!strcmp(argv[i], "-krolling"))
4496 // {
4497 // Mdouble old = getSpecies(0)->getRollingStiffness();
4498 // getSpecies(0)->setRollingStiffness(atof(argv[i + 1]));
4499 // std::cout << " reset krolling from " << old << " to " << getSpecies(0)->getRollingStiffness() << std::endl;
4500 // }
4501 // else if (!strcmp(argv[i], "-disprolling"))
4502 // {
4503 // Mdouble old = getSpecies(0)->getRollingDissipation();
4504 // getSpecies(0)->setRollingDissipation(atof(argv[i + 1]));
4505 // std::cout << " reset disprolling from " << old << " to " << getSpecies(0)->getRollingDissipation() << std::endl;
4506 // }
4507 // else if (!strcmp(argv[i], "-mu"))
4508 // {
4509 // Mdouble old = getSpecies(0)->getSlidingFrictionCoefficient();
4510 // getSpecies(0)->setSlidingFrictionCoefficient(atof(argv[i + 1]));
4511 // std::cout << " reset mu from " << old << " to " << getSpecies(0)->getSlidingFrictionCoefficient() << std::endl;
4512 // }
4513 // else if (!strcmp(argv[i], "-murolling"))
4514 // {
4515 // Mdouble old = getSpecies(0)->getRollingFrictionCoefficient();
4516 // getSpecies(0)->setRollingFrictionCoefficient(atof(argv[i + 1]));
4517 // std::cout << " reset murolling from " << old << " to " << getSpecies(0)->getRollingFrictionCoefficient() << std::endl;
4518 // }
4519  else if (!strcmp(argv[i], "-randomise") || !strcmp(argv[i], "-randomize"))
4520  {
4521  random.randomise();
4522  i--;
4523  }
4524 // else if (!strcmp(argv[i], "-k0"))
4525 // {
4526 // Mdouble old = speciesHandler.getObject(0)->getAdhesionStiffness();
4527 // speciesHandler.getObject(0)->setAdhesionStiffness(atof(argv[i + 1]));
4528 // std::cout << " reset k0 from " << old << " to " << speciesHandler.getObject(0)->getAdhesionStiffness() << std::endl;
4529 // }
4530 // else if (!strcmp(argv[i], "-f0"))
4531 // {
4532 // Mdouble old = speciesHandler.getObject(0)->getBondForceMax();
4533 // speciesHandler.getObject(0)->setBondForceMax(atof(argv[i + 1]));
4534 // std::cout << " reset f0 from " << old << " to " << speciesHandler.getObject(0)->getBondForceMax() << std::endl;
4535 // }
4536 // else if (!strcmp(argv[i], "-AdhesionForceType"))
4537 // {
4538 // AdhesionForceType old = speciesHandler.getObject(0)->getAdhesionForceType();
4539 // speciesHandler.getObject(0)->setAdhesionForceType(argv[i + 1]);
4540 // std::cout << " reset AdhesionForceType from "
4541 // << static_cast<signed char>(old) << " to "
4542 // << static_cast<signed char>(speciesHandler.getObject(0)->getAdhesionForceType()) << std::endl;
4543 // }
4544  else if (!strcmp(argv[i], "-append"))
4545  {
4546  setAppend(true);
4547  i--;
4548  }
4549  else if (!strcmp(argv[i], "-fixedParticles"))
4550  {
4551  setFixedParticles(static_cast<unsigned int>(atoi(argv[i + 1])));
4552  }
4553 // else if (!strcmp(argv[i], "-rho"))
4554 // {
4555 // Mdouble old = speciesHandler.getObject(0)->getDensity();
4556 // speciesHandler.getObject(0)->setDensity(atof(argv[i + 1]));
4557 // std::cout << " reset rho from " << old << " to " << speciesHandler.getObject(0)->getDensity() << std::endl;
4558 // }
4559 // else if (!strcmp(argv[i], "-dim_particle"))
4560 // {
4561 // setParticleDimensions(atoi(argv[i + 1]));
4562 // }
4563  else if (!strcmp(argv[i], "-counter"))
4564  {
4565  setRunNumber(atoi(argv[i + 1]));
4566  }
4567  else
4568  {
4569  //returns false if the flag passed does not match any of the currently used flags.
4570  return false;
4571  }
4572  //returns true if argv is found
4573  return true;
4574 }
4575 
4589 {
4590 #ifdef MERCURY_USE_MPI
4591  if (NUMBER_OF_PROCESSORS == 1)
4592  {
4594  }
4595 
4596  int localInteraction = checkParticleForInteractionLocal(p);
4597  //The root gathers all values and computes the global value
4598  int *interactionList = nullptr;
4599  if (PROCESSOR_ID == 0)
4600  {
4601  interactionList = new int [NUMBER_OF_PROCESSORS];
4602  }
4603 
4604  //Gather all local values
4605  MPIContainer::Instance().gather(localInteraction,interactionList);
4606 
4607  //Compute the global value
4608  int globalInteraction = 1;
4609  if (PROCESSOR_ID == 0)
4610  {
4611  for (int i = 0; i < NUMBER_OF_PROCESSORS; i++)
4612  {
4613  if (interactionList[i] == 0)
4614  {
4615  globalInteraction = 0;
4616  break;
4617  }
4618  }
4619  }
4620  //The root now tells the other processors what the global value for the interaction is
4621  MPIContainer::Instance().broadcast(globalInteraction);
4622 
4623  //Convert the result back to bool
4624  bool interaction = globalInteraction;
4625 #else
4626  bool interaction = checkParticleForInteractionLocalPeriodic(p);
4627 #endif
4628  return interaction;
4629 }
4630 
4638 {
4639  //A vector of ghost particles of the particle that is to be inserted (empty if no periodic boundaries are present)
4640  std::vector<Vec3D> pPeriodic;
4641  for (BaseBoundary* b : boundaryHandler)
4642  {
4643  PeriodicBoundary* pb = dynamic_cast<PeriodicBoundary*>(b);
4644  if (pb && particleHandler.getNumberOfObjects() > 0 )
4645  {
4647  for (int i = pPeriodic.size() - 1; i >= 0; --i)
4648  {
4649  if (pb->getDistance(pPeriodic[i]) < maxDistance)
4650  {
4651  pPeriodic.push_back(pPeriodic[i]);
4652  pb->shiftPosition(pPeriodic.back());
4653  }
4654  }
4655  if (pb->getDistance(p) < maxDistance)
4656  {
4657  pPeriodic.push_back(p.getPosition());
4658  pb->shiftPosition(pPeriodic.back());
4659  }
4660  }
4661  }
4662  //check the particle AND the ghost particles for intersection problems
4663  bool insertable = checkParticleForInteractionLocal(p);
4664  if (!pPeriodic.empty()) {
4665  BaseParticle* q = p.copy();
4666  for (const Vec3D& pos : pPeriodic) {
4667  q->setPosition(pos);
4668  insertable &= checkParticleForInteractionLocal(*q);
4669  }
4670  delete q;
4671  }
4672  return insertable;
4673 }
4674 
4688 {
4689  Mdouble distance;
4690  Vec3D normal;
4691 
4692  //Check if it has no collision with walls
4693  for (BaseWall* w : wallHandler)
4694  {
4695  //returns false if the function getDistanceAndNormal returns true,
4696  //i.e. if there exists an interaction between wall and particle
4697  //\todo TW getDistanceAndNormal(p,distance,normal) should ideally be replaced by a inContact(p) function, as it doesn't require distance and normal for anything (and walls now can have multiple contacts, soon particles can have it too.)
4698  if (w->getDistanceAndNormal(p, distance, normal))
4699  {
4700  //std::cout<<"failure: Collision with wall: "<<**it<<std::endl;
4701  return false;
4702  }
4703  else
4704  {
4705  //std::cout<<"No collision with wall: "<<**it<<std::endl;
4706  }
4707  }
4708 
4709  //Check if it has no collision with other particles
4710  for (BaseParticle* q : particleHandler)
4711  {
4712  //returns false if the particle separation is less than the relevant sum of interaction radii
4713  //(i.e. another particle is in contact with P)
4714  if (Vec3D::getDistanceSquared(q->getPosition(), p.getPosition())
4716  {
4717  //std::cout<<"failure: Collision with particle "<<**it<<std::endl;
4718  return false;
4719  }
4720  else
4721  {
4722  //std::cout<<"No collision with particle "<<**it<<std::endl;
4723  }
4724  }
4725  return true;
4727 }
4728 
4737 void DPMBase::importParticlesAs(ParticleHandler& particleH, InteractionHandler& interactionH, const ParticleSpecies* species )
4738 {
4739  size_t nParticlesPreviouslyIn = particleHandler.getSize();
4740  int l = 0;
4741  for (auto k = particleH.begin(); k != particleH.end(); ++k) {
4742  auto p = particleHandler.copyAndAddObject( *k );
4743  p->setSpecies(species);
4744  l++;
4745  }
4746 
4747  for (std::vector<BaseInteraction*>::const_iterator i = interactionH.begin(); i != interactionH.end(); ++i) {
4748  if ( (*i)->getP()->getInvMass() != 0.0 && (*i)->getI()->getInvMass() != 0.0 ) {
4750  j->importP(particleHandler.getObject(nParticlesPreviouslyIn + j->getP()->getIndex()));
4751  j->importI(particleHandler.getObject(nParticlesPreviouslyIn + j->getI()->getIndex()));
4752  j->setTimeStamp(getNumberOfTimeSteps());
4753  }
4754  }
4755 }
4756 
4772 {
4773 #ifdef MERCURY_USE_MPI
4774  if (NUMBER_OF_PROCESSORS == 1)
4776  {
4777 #endif
4778  //Looping from the *back* of the particle handler, where all the periodic or "ghost" particles are stored.
4779  //The loop ends when there exist no more periodic particles (i.e. "getPeriodicFromParticle()" returns a null pointer)
4780  for (unsigned int i = particleHandler.getSize();
4781  i >= 1 && particleHandler.getObject(i - 1)->getPeriodicFromParticle() != nullptr; i--)
4782  {
4783  while (!particleHandler.getObject(i - 1)->getInteractions().empty())
4784  {
4786  particleHandler.getObject(i - 1)->getInteractions().front()->getIndex());
4787  }
4789  }
4790 
4791  // OMP parallelism
4792  /*#pragma omp parallel for num_threads(getNumberOfOMPThreads()) //schedule(dynamic)
4793  for (unsigned int i = particleHandler.getSize(); i >= 1 ; i--)
4794  {
4795  if (particleHandler.getObject(i - 1)->getPeriodicFromParticle() != nullptr) {
4796  while (!particleHandler.getObject(i - 1)->getInteractions().empty()) {
4797  interactionHandler.removeObjectKeepingPeriodics(
4798  particleHandler.getObject(i - 1)->getInteractions().front()->getIndex());
4799  }
4800  particleHandler.removeObject(i - 1);
4801  }
4802  }*/
4803 
4804 #ifdef MERCURY_USE_MPI
4805  }
4806 #endif
4807 }
4808 
4815 {
4816  //Looping over all boundaries in the boundaryHandler
4817  for (BaseBoundary* boundary : boundaryHandler)
4818  {
4819  //Calls the createPeriodicParticles() function which checks if a particle is adequately
4820  //close to a periodic particle that a periodic (ghost) particle should be created and,
4821  //if so, creates one and adds it to the system (hence the necessity to keep "N" variable).
4822  //
4823  // (The loop is over all boundaries, but if a boundary is not a PeriodicBoundary, then
4824  // this does nothing.)
4825  boundary->createPeriodicParticles(particleHandler);
4826  }
4827 
4828  // OMP parallelism
4829  /*#pragma omp parallel for num_threads(getNumberOfOMPThreads()) //schedule(dynamic)
4830  for (int k = 0; k < boundaryHandler.getNumberOfObjects(); k++)
4831  {
4832  //Calls the createPeriodicParticles() function which checks if a particle is adequately
4833  //close to a periodic particle that a periodic (ghost) particle should be created and,
4834  //if so, creates one and adds it to the system (hence the necessity to keep "N" variable).
4835  //
4836  // (The loop is over all boundaries, but if a boundary is not a PeriodicBoundary, then
4837  // this does nothing.)
4838 
4839  BaseBoundary* boundary = boundaryHandler.getObject(k);
4840  #pragma omp critical
4841  boundary->createPeriodicParticles(particleHandler);
4842  }*/
4843 }
4844 
4849 {
4850 #ifdef MERCURY_USE_MPI
4851  //MPIContainer& communicator = MPIContainer::Instance();
4852  if (NUMBER_OF_PROCESSORS == 1) {return;}
4853 
4854  //Update the postion and velocity data of ghosts and perform some bookkeeping
4855  std::set<BaseParticle*> particlesToBeDeleted;
4856  domainHandler.updateStatus(particlesToBeDeleted);
4857  periodicBoundaryHandler.updateStatus(particlesToBeDeleted);
4858 
4859  //Delete particles
4860  deleteGhostParticles(particlesToBeDeleted);
4861 
4862  //Add new particles
4865 #endif
4866 }
4867 
4871 void DPMBase::deleteGhostParticles(std::set<BaseParticle*>& particlesToBeDeleted)
4872 {
4873  //Flush mixed particles from lists (MP particles are located in both structures)
4874  if (periodicBoundaryHandler.getSize() > 0)
4875  {
4876  //Flush particles from boundaries
4877  domainHandler.getCurrentDomain()->flushParticles(particlesToBeDeleted);
4878  periodicBoundaryHandler.flushParticles(particlesToBeDeleted);
4879  }
4880 
4881  //Clean communication lists
4884 
4885  //Delete the particles
4886  for (auto particle_it : particlesToBeDeleted)
4887  {
4888  particleHandler.removeGhostObject(particle_it->getIndex());
4889  }
4890 }
4891 
4892 
4893 //This function takes a base particle and then copies all the particle information from the root particle to the other particles
4894 // on neighbouring domains.
4895 void DPMBase::synchroniseParticle(BaseParticle* p, unsigned fromProcessor)
4896 {
4897 #ifdef MERCURY_USE_MPI
4898  MPIContainer& communicator = MPIContainer::Instance();
4899 
4900  //The processor that contains the particle that needs to be copied needs to identify the target, and communicate this
4901  MPIParticle pInfo;
4902  if (communicator.getProcessorID() == fromProcessor)
4903  {
4905  }
4906 
4907  //Broadcast from processor i
4908  communicator.broadcast(&pInfo,MercuryMPIType::PARTICLE,fromProcessor);
4910 #endif
4911 }
4912 
4914 {
4915 #ifdef MERCURY_USE_MPI
4916  if (NUMBER_OF_PROCESSORS == 1) {return;}
4917  //TODO If required, I can implement this for periodic particles, first discuss with Thomas if it is actually requiredf
4918  //periodicDomainHandler.updateVelocity()
4919  //domainHandler.updateVelocity();
4920 #endif
4921 }
4922 
4928 {
4929  std::cout << "Interactions currently in the handler:" << std::endl;
4930  //looping over all individual objects in the interactionHandler
4932  {
4933  p->write(std::cout);
4934  std::cout << std::endl;
4935  std::cout << "Interaction " << p->getName() << " " << p->getId() << " between " << p->getP()->getId() << " and "
4936  << p->getI()->getId() << std::endl;
4937  }
4938 }
4939 
4949 {
4950  return getTime() <= time && getTime() + getTimeStep() > time;
4951 }
4952 
4958 void DPMBase::setNumberOfDomains(std::vector<unsigned> numberOfDomains)
4959 {
4960 #ifdef MERCURY_USE_MPI
4961  numberOfDomains_ = numberOfDomains;
4962  logger(INFO, "Split domain into a %x%x% grid",numberOfDomains[0],numberOfDomains[1],numberOfDomains[2]);
4963 #else
4964  logger(WARN, "Setting number of domains, but code is not compiled with MPI on");
4965 #endif
4966 }
4967 
4969  //one-d problems
4970  if (domainSplit == DomainSplit::X) {
4972  return;
4973  } else if (domainSplit == DomainSplit::Y) {
4975  return;
4976  } else if (domainSplit == DomainSplit::Z) {
4978  return;
4979  }
4980  //two-d problems
4981  // split into axb grid with a the largest integer that divides NUMBER_OF_PROCESSORS and is smaller than sqrt(NUMBER_OF_PROCESSORS)
4982  unsigned a;
4983  for (unsigned n = floor(sqrt(NUMBER_OF_PROCESSORS));n>0; n--) {
4984  if (NUMBER_OF_PROCESSORS % n == 0) {
4985  a = n;
4986  break;
4987  }
4988  }
4989  if (domainSplit == DomainSplit::XY) {
4991  return;
4992  } else if (domainSplit == DomainSplit::XZ) {
4994  return;
4995  } else if (domainSplit == DomainSplit::YZ) {
4997  return;
4998  }
4999  //three-d problems
5000  // split into axbxc grid with
5001  // - a the largest integer that divides NUMBER_OF_PROCESSORS and is smaller than cbrt(NUMBER_OF_PROCESSORS)
5002  // - b the largest integer that divides NUMBER_OF_PROCESSORS/a and is smaller than sqrt(NUMBER_OF_PROCESSORS/a)
5003  unsigned b;
5004  for (unsigned n = floor(cbrt(NUMBER_OF_PROCESSORS));n>0; n--) {
5005  if (NUMBER_OF_PROCESSORS % n == 0) {
5006  a = n;
5007  break;
5008  }
5009  }
5010  for (unsigned n = floor(sqrt(NUMBER_OF_PROCESSORS/a));n>0; n--) {
5011  if (NUMBER_OF_PROCESSORS % (n*a) == 0) {
5012  b = n;
5013  break;
5014  }
5015  }
5017 }
5018 
5019 
5025 std::vector<unsigned> DPMBase::getNumberOfDomains()
5026 {
5027  return numberOfDomains_;
5028 }
5029 
5035 {
5037 }
5038 
5040 {
5041  return vtkWriter_;
5042 }
5043 
5051 {
5052  Vec3D meanVelocity = getTotalMomentum() / getTotalMass();
5053 
5054  //correct the mean velocity to zero
5055  for (auto& p : particleHandler)
5056  {
5057  p->addVelocity(-meanVelocity);
5058  }
5059 }
5060 
5068 {
5069  Vec3D V_mean;
5070  Mdouble Ek_mean_goal = 0, Ek_fluct_factor = 0;
5071  RNG rng;
5072 
5073  //assign random velocity to each particle
5074  for (auto& p : particleHandler)
5075  {
5076  p->setVelocity(Vec3D(rng.getRandomNumber(-1, 1), rng.getRandomNumber(-1, 1), rng.getRandomNumber(-1, 1)));
5077  }
5078 
5079  //calculate the mean velocity in the system now
5080  Ek_mean_goal = 0.5 * getTotalMass() * V_mean_goal.getLengthSquared();
5081  V_mean = getTotalMomentum() / getTotalMass();
5082  //check if the user input mean kinetic energy is larger than the total kinetic energy input, then return error
5083  logger.assert_always(0.5 * getTotalMass() * V_mean_goal.getLengthSquared() < Ek_goal,
5084  "Too large mean velocity input, Kinetic energy from mean velocity part is larger than the "
5085  "total kinetic energy you want to set");
5086 
5087  //correct the mean velocity to zero
5088  for (auto& p : particleHandler)
5089  {
5090  p->addVelocity(-V_mean);
5091  }
5092 
5093  //set the new fluctuating velocity based on the goal fluctuating kinetic energy
5094  Ek_fluct_factor = std::sqrt((Ek_goal - Ek_mean_goal) / getKineticEnergy());
5095  for (auto& p : particleHandler)
5096  {
5097  p->setVelocity(Ek_fluct_factor * p->getVelocity());
5098  }
5099 
5100  //correct the mean velocity finally to the user set values
5101  V_mean = getTotalMomentum() / getTotalMass();
5102  for (auto& p : particleHandler)
5103  {
5104  p->addVelocity(V_mean_goal - V_mean);
5105  }
5106 
5107  //check the final mean velocity and kinetic energy
5108  logger(INFO, "In DPMBase::setMeanVelocityAndKineticEnergy,\nV_mean_final %\n Ek_final %\n",
5110 }
5111 
5116 {
5117  return (getXMax() - getXMin()) * (getYMax() - getYMin()) * (getZMax() - getZMin());
5118 }
5119 
5125 {
5126  Matrix3D F; //set the kinetic energy tensor, this is in terms of Sum(m*v^2)
5127  Vec3D J; //set the momentum tensor
5128 
5129  //calculate stress for kinetic part
5130  for (const auto& p : particleHandler)
5131  {
5132  F += Matrix3D::dyadic(p->getVelocity(), p->getVelocity()) * p->getMass();
5133  J += p->getVelocity() * p->getMass();
5134  }
5135 
5136  Matrix3D stressKinetic = F - Matrix3D::dyadic(J, J) / getTotalMass();
5137  stressKinetic /= getTotalVolume();
5138  return stressKinetic;
5139 }
5140 
5147 {
5148  //stress components calculation variables
5149  Matrix3D stressStatic;
5150 
5151  //calculate the static stress tensor based on all the interactions
5152  for (const auto i : interactionHandler)
5153  {
5154  stressStatic += Matrix3D::dyadic(i->getForce(), i->getNormal()) * i->getDistance();
5155  }
5156 
5157  stressStatic /= getTotalVolume();
5158  return stressStatic;
5159 }
5160 
5167 {
5168  return getKineticStress() + getStaticStress();
5169 }
5170 
5171 
5173 {
5174  //compute forces for all particles that are neither fixed or ghosts
5175  for (auto p : particleHandler)
5176  {
5177  if (!p->isFixed() && p->getPeriodicFromParticle() == nullptr)
5178  {
5179  //w->computeForces(p);
5181  }
5182  }
5183 }
5184 
5192 void DPMBase::setLogarithmicSaveCount(const Mdouble logarithmicSaveCountBase)
5193 {
5194  dataFile.setlogarithmicSaveCount(logarithmicSaveCountBase);
5195  fStatFile.setlogarithmicSaveCount(logarithmicSaveCountBase);
5196  restartFile.setlogarithmicSaveCount(logarithmicSaveCountBase);
5197  statFile.setlogarithmicSaveCount(logarithmicSaveCountBase);
5198  eneFile.setlogarithmicSaveCount(logarithmicSaveCountBase);
5199 }
5200 
virtual void computeExternalForces(BaseParticle *)
Computes the external forces, such as gravity, acting on particles.
Definition: DPMBase.cc:3049
Implementation of a 3D quaternion (by Vitaliy).
Definition: Quaternion.h:62
This class contains all information and functions required for communication between processors...
Definition: MpiContainer.h:125
Mdouble timeMax_
Stores the duration of the simulation.
Definition: DPMBase.h:1224
int launchNewRun(const char *name, bool quick=false)
This launches a code from within this code. Please pass the name of the code to run.
Definition: DPMBase.cc:764
void setFileCounter(unsigned fileCounter)
Definition: BaseVTKWriter.h:61
void shiftPosition(BaseParticle *p) const override
shifts the particle
void setTime(Mdouble time)
Sets a new value for the current simulation time.
Definition: DPMBase.cc:825
void copyDataFromParticleToMPIParticle(BaseParticle *p)
BaseParticle * getLargestParticle() const
Returns the pointer of the largest particle in the particle handler. When mercury is running in paral...
virtual void actionsBeforeTimeLoop()
A virtual function. Allows one to carry out any operations before the start of the time loop...
Definition: DPMBase.cc:1629
unsigned int getId() const
Returns the unique identifier of any particular object.
Definition: BaseObject.h:125
const Vec3D & getTorque() const
Gets the current torque (vector) between the two interacting objects.
void addNewParticle(BaseParticle *particle)
Adds a new particle to the periodic list.
void setXBallsVectorScale(double newVScale)
Set the scale of vectors in xballs.
Definition: DPMBase.cc:1280
void finish()
Contains the code executed after the last time step.
Definition: CGHandler.cc:113
virtual void write(std::ostream &os, bool writeAllParticles=true) const
Loads all MD data and plots statistics for all time steps in the .data file.
Definition: DPMBase.cc:3375
FileType getWallsWriteVTK() const
Returns whether walls are written in a VTK file.
Definition: DPMBase.cc:935
virtual void hGridUpdateParticle(BaseParticle *obj UNUSED)
Definition: DPMBase.cc:1664
void set(Vec3D normal, Mdouble distanceLeft, Mdouble distanceRight)
Defines a PeriodicBoundary by its normal and positions.
void setXMax(Mdouble newXMax)
Sets the value of XMax, the upper bound of the problem domain in the x-direction. ...
Definition: DPMBase.cc:1126
void eraseOldInteractions(unsigned)
erases interactions which have an old timestamp.
bool writeToFile(std::string filename, std::string filecontent)
Writes a string to a file.
Definition: Helpers.cc:446
void setWallsWriteVTK(FileType writeWallsVTK)
Sets whether walls are written into a VTK file.
Definition: DPMBase.cc:878
unsigned getFileCounter() const
Definition: BaseVTKWriter.h:56
void write(std::ostream &os) const
void solve()
The work horse of the code.
Definition: DPMBase.cc:3968
virtual void hGridInsertParticle(BaseParticle *obj UNUSED)
Definition: DPMBase.cc:1657
void writeVTK() const override
writes a vtk file
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
The DPMBase header includes quite a few header files, defining all the handlers, which are essential...
Definition: DPMBase.h:72
void setLastSavedTimeStep(unsigned int nextSavedTimeStep)
Sets the next time step for all the files (ene, data, fstat, restart, stat) at which the data is to b...
Definition: DPMBase.cc:508
Mdouble getRotationalEnergy() const
JMFT Returns the global rotational energy stored in the system.
Definition: DPMBase.cc:1537
Mdouble X
the vector components
Definition: Vector.h:65
void removeObject(unsigned int index) override
Removes a BaseParticle from the ParticleHandler.
each time-step will be written into/read from separate files numbered consecutively, with numbers padded by zeros to a minimum of four digits: name_.0000, name_.0001, ..
Domain * getCurrentDomain()
Function that returns a pointer to the domain corresponding to the processor.
Definition: DPMBase.cc:5034
bool getWriteVTK() const
virtual void integrateAfterForceComputation()
Update particles' and walls' positions and velocities after force computation.
Definition: DPMBase.cc:3194
void initialiseMercuryMPITypes(const SpeciesHandler &speciesHandler)
Creates the MPI types required for communication of Mercury data through the MPI interface.
Definition: MpiContainer.cc:70
const Mdouble NaN
Definition: GeneralDefine.h:43
void checkAndDuplicatePeriodicParticles()
For simulations using periodic boundaries, checks and adds particles when necessary into the particle...
Definition: DPMBase.cc:4814
void read(std::istream &is)
Definition: RNG.cc:59
bool saveCurrentTimeStep(unsigned int ntimeSteps)
determined if this time step has to be written; if so, opens the output file
Definition: File.cc:313
A basic particle.
void addParticle(BaseParticle *particle)
Initialises a single particle which is added during the simulation.
Definition: Domain.cc:1581
std::string name_
the name of the problem, used, e.g., for the output files
Definition: DPMBase.h:1302
void setVelocity(const Vec3D &velocity)
set the velocity of the BaseInteractable.
static bool areInContact(const BaseParticle *pI, const BaseParticle *pJ)
Checks if two particle are in contact or is there any positive overlap.
Definition: DPMBase.cc:1621
int getXBallsColourMode() const
Get the xballs colour mode (CMode).
Definition: DPMBase.cc:1270
Matrix3D getStaticStress() const
Calculate the static stress tensor in the system averaged over the whole volume.
Definition: DPMBase.cc:5146
unsigned int getSize() const
Gets the size of the particleHandler (including mpi and periodic particles)
Definition: BaseHandler.h:655
void writeVTKFiles() const
Definition: DPMBase.cc:2091
void setTimeMax(Mdouble newTMax)
Sets a new value for the maximum simulation duration.
Definition: DPMBase.cc:840
virtual void writeVTK() const =0
void constructor()
A function which initialises the member variables to default values, so that the problem can be solve...
Definition: DPMBase.cc:200
void readAndAddObject(std::istream &is) final
Create a new wall in the WallHandler, based on the information provided in a restart file...
Definition: WallHandler.cc:276
void write(std::ostream &os) const
Writes the InteractionHandler to an output stream, for example a restart file.
unsigned int getSystemDimensions() const
Returns the system dimensionality.
Definition: DPMBase.cc:1390
bool readArguments(int argc, char *argv[])
Can interpret main function input arguments that are passed by the driver codes.
Definition: DPMBase.cc:4174
FileType getWriteVTK() const
FileType getFileType() const
Gets the file type e.g. NOFILE, ONEFILE and MULTIPLE FILES. File::fileType_.
Definition: File.cc:208
bool mpiInsertParticleCheck(BaseParticle *P)
Function that checks if the mpi particle should really be inserted by the current domain...
Definition: DPMBase.cc:1692
virtual bool isInContactWith(const BaseParticle *P) const
Get whether or not this particle is in contact with the given particle.
BaseParticle * getP1()
Gets a pointer to the first BaseParticle in this PossibleContact.
File & getInteractionFile()
Return a reference to the file InteractionsFile.
Definition: DPMBase.cc:335
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
Vec3D max_
Definition: DPMBase.h:1204
double Mdouble
Definition: GeneralDefine.h:34
const ParticleSpecies * getSpecies() const
Returns a pointer to the species of this BaseInteractable.
const std::string getSVNURL()
unsigned int particleDimensions_
determines if 2D or 3D particle volume is used for mass calculations
Definition: DPMBase.h:1188
virtual void setRadius(Mdouble radius)
Sets the particle's radius_ (and adjusts the mass_ accordingly, based on the particle's species) ...
void addForce(const Vec3D &addForce)
Adds an amount to the force on this BaseInteractable.
unsigned int getNumberOfObjects() const override
Returns the number of objects in the container. In parallel code this practice is forbidden to avoid ...
std::string getPath()
Definition: Helpers.cc:918
Mdouble getMaxInteractionRadius() const
Returns the particle's interaction radius, which might be different from radius_ (e.g., when dealing with wet particles)
Definition: BaseParticle.h:359
void splitDomain(DomainSplit domainSplit)
Definition: DPMBase.cc:4968
void setMax(const Vec3D &max)
Sets the maximum coordinates of the problem domain.
Definition: DPMBase.cc:1043
void setDPMBase(DPMBase *DPMBase)
Sets the problem that is solved using this handler.
Definition: BaseHandler.h:718
void resetForceTorque(int numberOfOMPthreads)
void setYMin(Mdouble newYMin)
Sets the value of YMin, the lower bound of the problem domain in the y-direction. ...
Definition: DPMBase.cc:995
virtual void processStatistics(bool)
Definition: DPMBase.cc:1878
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
Vec3D getCentreOfMass() const
JMFT: Return the centre of mass of the system, excluding fixed particles.
Definition: DPMBase.cc:1574
void setLastSavedTimeStep(unsigned int lastSavedTimeStep)
Sets File::nextSavedTimeStep_.
Definition: File.cc:303
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
bool readSpeciesFromDataFile_
Determines if the last column of the data file is interpreted as the info parameter during restart...
Definition: DPMBase.h:1312
void read(std::istream &is)
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...
virtual void decompose()
Sends particles from processorId to the root processor.
Definition: DPMBase.cc:3885
void setStorageCapacity(const unsigned int N)
Sets the storage capacity of this BaseHandler.
Definition: BaseHandler.h:669
unsigned int getParticleDimensions() const
Returns the particle dimensionality.
Definition: DPMBase.cc:1427
Vec3D getMin() const
Definition: DPMBase.h:623
void setParticlesWriteVTK(bool writeParticlesVTK)
Sets whether particles are written in a VTK file.
Definition: DPMBase.cc:904
void computeAllMasses(unsigned int indSpecies)
Computes the mass for all BaseParticle of the given species in this ParticleHandler.
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
virtual void initialiseStatistics()
Definition: DPMBase.cc:1845
void setLogarithmicSaveCount(Mdouble logarithmicSaveCountBase)
Sets File::logarithmicSaveCount_ for all files (ene, data, fstat, restart, stat)
Definition: DPMBase.cc:5192
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
WallVTKWriter wallVTKWriter_
Definition: DPMBase.h:1265
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 setXBallsColourMode(int newCMode)
Set the xballs output mode.
Definition: DPMBase.cc:1260
const std::vector< T * >::const_iterator end() const
Gets the end of the const_iterator over all BaseBoundary in this BaseHandler.
Definition: BaseHandler.h:704
File interactionFile
File class to handle in- and output into .interactions file. This file hold information about interac...
Definition: DPMBase.h:1396
virtual void readOld(std::istream &is)
Reads all data from a restart file, e.g. domain data and particle data; old version.
Definition: DPMBase.cc:3739
Mdouble getSumOfInteractionRadii(const BaseParticle *particle) const
returns the sum of the radii plus the interactionDistance
Definition: BaseParticle.h:376
const std::string & getName() const
Returns the name of the file. Does not allow to change it though.
Definition: DPMBase.cc:389
bool containsParticle(BaseParticle *particle, Mdouble offset=0.0)
Check to see if a given particle is within the current domain.
Definition: Domain.cc:400
void removeObjectKeepingPeriodics(unsigned int id)
Removes interactions of periodic particles when the periodic particles get deleted (see DPMBase::remo...
void integrateBeforeForceComputation(double time, double timeStep)
This is part of integrate routine for objects with infinite mass.
InteractionVTKWriter interactionVTKWriter_
Definition: DPMBase.h:1267
void evaluate()
Contains the code executed at each time step.
Definition: CGHandler.cc:97
int runNumber_
This stores the run number for saving.
Definition: DPMBase.h:1297
virtual void outputXBallsDataParticle(unsigned int i, unsigned int format, std::ostream &os) const
This function writes out the particle locations into an output stream in a format the XBalls program ...
void setParticleDimensions(unsigned int particleDimensions)
Sets the particle dimensionality.
Definition: DPMBase.cc:1408
virtual ~DPMBase()
virtual destructor
Definition: DPMBase.cc:287
void setOrientationViaEuler(Vec3D eulerAngle)
Sets the orientation of this BaseInteractable by defining the euler angles.
virtual void actionsAfterTimeStep()
A virtual function which allows to define operations to be executed after time step.
Definition: DPMBase.cc:1838
void writePythonFileForVTKVisualisation() const
record when the simulation started
Definition: DPMBase.cc:2127
Mdouble getInteractionDistance()
Gets the interaction distance of the domain handler.
BaseParticle * getPeriodicFromParticle() const
Returns the 'original' particle this one's a periodic copy of.
Definition: BaseParticle.h:338
void unfix()
Unfix Particle function, which required a reference to the Species vector. It unfixes a Particle by c...
DomainHandler domainHandler
An object of the class DomainHandler which deals with parallel code.
Definition: DPMBase.h:1354
void setZMax(Mdouble newZMax)
Sets the value of ZMax, the upper bound of the problem domain in the z-direction. ...
Definition: DPMBase.cc:1178
virtual void computeWallForces(BaseWall *w)
Definition: DPMBase.cc:5172
void flushParticles(std::set< BaseParticle * > &particlesToBeFlushed)
Removes particles from the periodiocParticleList_ and periociGhostList_.
void setSpecies(const ParticleSpecies *species)
void setDimension(unsigned int newDim)
Sets both the system dimensions and the particle dimensionality.
Definition: DPMBase.cc:1363
Mdouble getRandomNumber()
This is a random generating routine can be used for initial positions.
Definition: RNG.cc:143
std::vector< unsigned > numberOfDomains_
Vector containing the number of domains in x-,y- and z-direction, required for parallel computations...
Definition: DPMBase.h:1198
const Vec3D & getContactPoint() const
Gets constant reference to contact point (vector).
MERCURY_DEPRECATED File & getEneFile()
The non const version. Allows to edit the File::eneFile.
Definition: DPMBase.cc:303
Mdouble getMass() const
static Mdouble dot(const Vec3D &a, const Vec3D &b)
Calculates the dot product of two Vec3D: .
Definition: Vector.cc:76
void insertParticles(DPMBase *md)
void synchroniseParticle(BaseParticle *, unsigned fromProcessor=0)
Definition: DPMBase.cc:4895
virtual BaseParticle * copy() const =0
Particle copy method. It calls to copy constructor of this Particle, useful for polymorfism.
BaseParticle * getP2()
Gets a pointer to the second BaseParticle in this PossibleContact.
void setAppend(bool newAppendFlag)
Sets whether the "append" option is on or off.
Definition: DPMBase.cc:1482
int getNumberOfOMPThreads() const
Definition: DPMBase.cc:1246
void setFStatData(std::fstream &fstat, BaseParticle *P, BaseWall *I)
virtual void writeOutputFiles()
Writes simulation data to all the main Mercury files: .data, .ene, .fstat, .xballs and ...
Definition: DPMBase.cc:3846
Species< LinearViscoelasticNormalSpecies, SlidingFrictionSpecies > LinearViscoelasticSlidingFrictionSpecies
void setSystemDimensions(unsigned int newDim)
Sets the system dimensionality.
Definition: DPMBase.cc:1377
BaseInteraction * getInteractionWith(BaseParticle *p, unsigned timeStamp, InteractionHandler *interactionHandler) override
Returns the interaction between this wall and a given particle, nullptr if there is no interaction...
Definition: BaseWall.cc:369
void setMeanVelocityAndKineticEnergy(Vec3D V_mean_goal, Mdouble Ek_goal)
This function will help you set a fixed kinetic energy and mean velocity in your system.
Definition: DPMBase.cc:5067
std::function< void(std::string, std::string)> onFatal
Definition: Logger.h:195
void read(std::istream &is)
Reads all objects from restart data.
Definition: BaseHandler.h:543
void clear() override
Empties the whole ParticleHandler by removing all BaseParticle.
virtual void write(std::ostream &os) const
Write all the species and mixed species to an output stream.
void setGravity(Vec3D newGravity)
Sets a new value for the gravitational acceleration.
Definition: DPMBase.cc:1343
Mdouble getDistance(const BaseParticle &p) const override
Returns the distance of the edge to the particle.
void initialiseMPI()
Inialises the MPI library.
virtual void outputXBallsData(std::ostream &os) const
This function writes the location of the walls and particles in a format the XBalls program can read...
Definition: DPMBase.cc:2240
virtual bool getHGridUpdateEachTimeStep() const
Definition: DPMBase.cc:1678
bool rotation_
A flag to turn on/off particle rotation. true will enable particle rotation. false will disable parti...
Definition: DPMBase.h:1249
Mdouble & z()
RW reference to Z.
Definition: Vector.h:368
Stores information about interactions between two interactable objects; often particles but could be ...
void randomise()
sets the random variables such that they differ for each run
Definition: RNG.cc:98
bool getParticlesWriteVTK() const
Returns whether particles are written in a VTK file.
Definition: DPMBase.cc:946
void initialise()
Contains the code executed before the first time step.
Definition: CGHandler.cc:90
void updateGhostGrid(BaseParticle *P)
Checks if the Domain/periodic interaction distance needs to be updated and updates it accordingly...
Definition: DPMBase.cc:1797
void addNewParticles()
unsigned int numberOfTimeSteps_
Stores the number of time steps.
Definition: DPMBase.h:1214
bool openWrite(unsigned)
First sets openmode to write (and append in some cases), then calls open().
Definition: File.cc:382
FileType
With FileType options, one is able to choose if data is to be read/written from/into no or single or ...
Definition: File.h:40
bool restarted_
A bool to check if the simulation was restarted or not, ie. if setupInitialConditionsShould be run an...
Definition: DPMBase.h:1236
void initialise()
Initialises the communication list vectors as they can not be determined on compile time...
const std::vector< T * >::const_iterator begin() const
Gets the begin of the const_iterator over all Object in this BaseHandler.
Definition: BaseHandler.h:690
Defines a pair of periodic walls. Inherits from BaseBoundary.
Mdouble getXMax() const
If the length of the problem domain in x-direction is XMax - XMin, then getXMax() returns XMax...
Definition: DPMBase.h:593
static Matrix3D dyadic(const Vec3D &a, const Vec3D &b)
Calculates the dyadic product of a two Vec3D: .
Definition: Matrix.cc:323
unsigned int getNumberOfRealObjects() const
Returns the number of real objects (on all processors)
const unsigned NEVER
Definition: File.h:35
Matrix3D getTotalStress() const
Calculate the total stress tensor in the system averaged over the whole volume.
Definition: DPMBase.cc:5166
unsigned int getSaveCount() const
Gets File::saveCount_.
Definition: File.cc:256
PossibleContact * getNext()
Gets the next PossibleContact in the general linked list of PossibleContact.
void writeVTK() const override
writes a vtk file
void close()
Closes the file by calling fstream_.close()
Definition: File.cc:408
void setAngleZ(Mdouble psi)
Converts the rotation angle in the XY plane into a quaternion (for Mercury2D). See Wikipedia for deta...
Definition: Quaternion.cc:492
void copyDataFromMPIParticleToParticle(MPIParticle *bP, BaseParticle *p, ParticleHandler *particleHandler)
Copies data from an MPIParticle class to a BaseParticle and sets the particleHandler and species...
bool findNextExistingDataFile(Mdouble tMin, bool verbose=true)
Finds and opens the next data file, if such a file exists.
Definition: DPMBase.cc:2538
virtual void writeEneTimeStep(std::ostream &os) const
Write the global kinetic, potential energy, etc. in the system.
Definition: DPMBase.cc:2067
void forceWriteOutputFiles()
Writes output files immediately, even if the current time step was not meant to be written...
Definition: DPMBase.cc:3825
void setYMax(Mdouble newYMax)
Sets the value of YMax, the upper bound of the problem domain in the y-direction. ...
Definition: DPMBase.cc:1152
unsigned int systemDimensions_
The dimensions of the simulation i.e. 2D or 3D.
Definition: DPMBase.h:1183
void readAndAddObject(std::istream &is) final
Reads BaseBoundary into the BoundaryHandler from restart data.
void writeFStatFile()
Definition: DPMBase.cc:2858
void gather(T &send_t, T *receive_t)
Gathers a scaler from all processors to a vector of scalars on the root.
Definition: MpiContainer.h:416
void setDomain(const Vec3D &min, const Vec3D &max)
Sets the minimum coordinates of the problem domain.
Definition: DPMBase.cc:1059
This is a class that generates random numbers i.e. named the Random Number Generator (RNG)...
Definition: RNG.h:52
Mdouble getElasticEnergy() const
Returns the global elastic energy within the system.
Definition: DPMBase.cc:1490
Matrix3D getKineticStress() const
Calculate the kinetic stress tensor in the system averaged over the whole volume. ...
Definition: DPMBase.cc:5124
Mdouble getGravitationalEnergy() const
Returns the global gravitational potential energy stored in the system.
Definition: DPMBase.cc:1521
std::size_t getProcessorID()
Reduces a scalar on all processors to one scalar on a target processor.
static void incrementRunNumberInFile()
Increment the run Number (counter value) stored in the file_counter (COUNTER_DONOTDEL) by 1 and store...
Definition: DPMBase.cc:617
unsigned int getCounter() const
In case of multiple files, File::getCounter() returns the the number (FILE::Counter_) of the next fil...
Definition: File.cc:224
void integrateBeforeForceComputation(double time, double timeStep)
First step of Velocity Verlet integration.
File dataFile
An instance of class File to handle in- and output into a .data file.
Definition: DPMBase.h:1370
file will not be created/read
bool readParAndIniFiles(std::string fileName)
Allows the user to read par.ini files (useful to read files produced by the MDCLR simulation code - e...
Definition: DPMBase.cc:2355
MERCURY_DEPRECATED File & getDataFile()
The non const version. Allows one to edit the File::dataFile.
Definition: DPMBase.cc:295
void setWriteVTK(bool writeVTK)
virtual void read(std::istream &is, ReadOptions opt=ReadOptions::ReadAll)
Reads all data from a restart file, e.g. domain data and particle data.
Definition: DPMBase.cc:3496
void setlogarithmicSaveCount(const Mdouble logarithmicSaveCountBase)
the function to set the user input base of logarithmic saving count
Definition: File.cc:284
virtual void hGridActionsBeforeIntegration()
This function has to be called before integrateBeforeForceComputation.
Definition: DPMBase.cc:1900
std::fstream & getFstream()
Allows to access the member variable File::fstream_.
Definition: File.cc:154
void set(BaseParticle *particleToCopy, unsigned int maxFailed, Vec3D posMin, Vec3D posMax, Vec3D velMin, Vec3D velMax, double radMin, double radMax)
Sets the properties of the CubeInsertionBoundary.
File fStatFile
An instance of class File to handle in- and output into a .fstat file.
Definition: DPMBase.h:1375
void write(std::ostream &os) const
Definition: RNG.cc:78
bool getAppend() const
Returns whether the "append" option is on or off.
Definition: DPMBase.cc:1470
virtual bool readNextArgument(int &i, int argc, char *argv[])
Interprets the i^th command-line argument.
Definition: DPMBase.cc:4299
virtual bool checkParticleForInteraction(const BaseParticle &P)
Checks whether a particle P has any interaction with walls or other particles.
Definition: DPMBase.cc:4588
void cleanCommunicationLists()
Removes nullptrs from boundaryParticleList_ and boundaryParticleListNeighbour_.
Definition: Domain.cc:1713
Vec3D getMassTimesPosition() const
It's an insertion boundary which has cuboidal shape (yes, 'CuboidalInsertionBoundary' would have been...
void performGhostVelocityUpdate()
updates the final time-step velocity of the ghost particles
Definition: DPMBase.cc:4913
bool getRestarted() const
Returns the flag denoting if the simulation was restarted or not.
Definition: DPMBase.cc:1453
Mdouble getMass() const
Returns the particle's mass.
Definition: BaseParticle.h:322
bool writeSuperquadricParticlesVTK_
Definition: DPMBase.h:1261
void getLineFromStringStream(std::istream &in, std::stringstream &out)
Reads a line from one stringstream into another, and prepares the latter for reading in...
Definition: Helpers.cc:424
CGHandler cgHandler
Object of the class cgHandler.
Definition: DPMBase.h:1365
unsigned int getNumberOfTimeSteps() const
Returns the current counter of time-steps, i.e. the number of time-steps that the simulation has unde...
Definition: DPMBase.cc:813
int xBallsColourMode_
XBalls is a package to view the particle data. As an alternative MercuryDPM also supports ParaView...
Definition: DPMBase.h:1277
void updateStatus(std::set< BaseParticle * > &ghostParticlesToBeDeleted)
Updates the positions/velocity of ghost particles and accordingly the status of these particles...
void setEuler(const Vec3D &e)
Convert Euler angles to a quaternion. See Wikipedia for details.
Definition: Quaternion.cc:473
BoundaryHandler boundaryHandler
An object of the class BoundaryHandler which concerns insertion and deletion of particles into or fro...
Definition: DPMBase.h:1344
void addNewParticles()
Adds new particles to the periodic particle lists.
PeriodicBoundaryHandler periodicBoundaryHandler
Internal handler that deals with periodic boundaries, especially in a parallel build.
Definition: DPMBase.h:1349
ParticleVtkWriter * vtkWriter_
Definition: DPMBase.h:1263
BoundaryVTKWriter boundaryVTKWriter_
Definition: DPMBase.h:1269
void fillDomainWithParticles(unsigned N=50)
Inserts particles in the whole domain.
Definition: DPMBase.cc:2874
virtual Mdouble getInfo(const BaseParticle &P) const
A virtual function that returns some user-specified information about a particle. ...
Definition: DPMBase.cc:1602
Mdouble xBallsScale_
sets the xballs argument scale (see XBalls/xballs.txt)
Definition: DPMBase.h:1287
void writeEneFile()
Definition: DPMBase.cc:2848
const std::string getVersion()
void setXMin(Mdouble newXMin)
Sets the value of XMin, the lower bound of the problem domain in the x-direction. ...
Definition: DPMBase.cc:971
ParticleHandler particleHandler
An object of the class ParticleHandler, contains the pointers to all the particles created...
Definition: DPMBase.h:1329
all data will be written into/ read from a single file called name_
Container to store Interaction objects.
BaseInteraction * addInteraction(BaseInteractable *P, BaseInteractable *I, unsigned timeStamp)
T * getObject(const unsigned int id)
Gets a pointer to the Object at the specified index in the BaseHandler.
Definition: BaseHandler.h:613
bool saveCurrentTimeStepNoFileTypeCheck(unsigned int ntimeSteps)
Definition: File.cc:318
Data class to send a particle over MPI.
Definition: MpiDataClass.h:81
void readAndAddObject(std::istream &is) override
Create a new particle, based on the information from old-style restart data.
Mdouble getKineticEnergy() const
Returns the global kinetic energy stored in the system.
Definition: DPMBase.cc:1504
Mdouble & x()
RW reference to X.
Definition: Vector.h:344
Mdouble time_
Stores the current simulation time.
Definition: DPMBase.h:1209
void setRestartVersion(std::string newRV)
Sets restart_version.
Definition: DPMBase.cc:1444
void setMin(const Vec3D &min)
Sets the minimum coordinates of the problem domain.
Definition: DPMBase.cc:1079
void setSaveCount(unsigned int saveCount)
Sets File::saveCount_ for all files (ene, data, fstat, restart, stat)
Definition: DPMBase.cc:398
void closeFiles()
Closes all files (ene, data, fstat, restart, stat) that were opened to read or write.
Definition: DPMBase.cc:492
Mdouble getTotalVolume() const
Get the total volume of the cuboid system.
Definition: DPMBase.cc:5115
Vec3D getGravity() const
Returns the gravitational acceleration.
Definition: DPMBase.cc:1351
void setOpenMode(std::fstream::openmode openMode)
Sets File::openMode_ for all files (ene, data, fstat, restart, stat)
Definition: DPMBase.cc:479
void checkSettings()
Checks if the essentials are set properly to go ahead with solving the problem.
Definition: DPMBase.cc:3796
#define UNUSED
Definition: GeneralDefine.h:39
#define PROCESSOR_ID
Definition: GeneralDefine.h:63
virtual void hGridActionsBeforeTimeLoop()
A virtual function that allows one to carry out hGrid operations before the start of the time loop...
Definition: DPMBase.cc:1636
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 getKineticEnergy() const
bool openWriteNoAppend(unsigned)
Definition: File.cc:399
const Vec3D & getForce() const
Gets the current force (vector) between the two interacting objects.
static Vec3D cross(const Vec3D &a, const Vec3D &b)
Calculates the cross product of two Vec3D: .
Definition: Vector.cc:163
virtual void computeAllForces()
Computes all the forces acting on the particles using the BaseInteractable::setForce() and BaseIntera...
Definition: DPMBase.cc:3277
Mdouble getRadius() const
Returns the particle's radius.
Definition: BaseParticle.h:345
bool mpiIsInCommunicationZone(BaseParticle *particle)
Checks if the position of the particle is in an mpi communication zone or not.
Definition: DPMBase.cc:1724
std::vector< int > get3DParametersFromRunNumber(int size_x, int size_y, int size_z) const
This turns a counter into 3 indices, which is a useful feature for performing a 3D parameter study...
Definition: DPMBase.cc:724
unsigned int getNumberOfRealObjectsLocal() const
Returns the number of real objects on a local domain. MPI particles and periodic particles are neglec...
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
void setInteractionsWriteVTK(bool)
Sets whether interactions are written into a VTK file.
Definition: DPMBase.cc:894
void removeGhostObject(unsigned int index)
Removes a BaseParticle from the ParticleHandler without a global check, this is only to be done for m...
Basic class for walls.
Definition: BaseWall.h:47
bool isFixed() const override
Is fixed Particle function. It returns whether a Particle is fixed or not, by checking its inverse Ma...
Definition: BaseParticle.h:93
void setZMin(Mdouble newZMin)
Sets the value of ZMin, the lower bound of the problem domain in the z-direction. ...
Definition: DPMBase.cc:1019
virtual void actionsAfterSolve()
A virtual function which allows to define operations to be executed after the solve().
Definition: DPMBase.cc:1831
void setSaveCount(unsigned int saveCount)
Sets File::saveCount_.
Definition: File.cc:274
void clear() override
Empties the whole BaseHandler by removing all Objects and setting all other variables to 0...
int numberOfOMPThreads_
Definition: DPMBase.h:1178
void writeVTK() const override
writes a vtk file
virtual void printTime() const
Displays the current simulation time and the maximum simulation duration.
Definition: DPMBase.cc:1930
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
FileType writeWallsVTK_
A flag to turn on/off the vtk writer for walls.
Definition: DPMBase.h:1254
virtual void writeXBallsScript() const
This writes a script which can be used to load the xballs problem to display the data just generated...
void setNumberOfOMPThreads(int numberOfOMPThreads)
Definition: DPMBase.cc:1218
const int getSVNRevision()
Mdouble & y()
RW reference to Y.
Definition: Vector.h:356
void outputInteractionDetails() const
Displays the interaction details corresponding to the pointer objects in the interaction handler...
Definition: DPMBase.cc:4927
Vec3D getMax() const
Definition: DPMBase.h:629
void setXBallsAdditionalArguments(std::string newXBArgs)
Set the additional arguments for xballs.
Definition: DPMBase.cc:1307
bool isPeriodicGhostParticle() const
Indicates if this particle is a ghost in the periodic boundary.
Container to store all BaseParticle.
#define NUMBER_OF_PROCESSORS
For the MPI communication routines this quantity is often required. defining this macro makes the cod...
Definition: GeneralDefine.h:62
static Mdouble getDistanceSquared(const Vec3D &a, const Vec3D &b)
Calculates the squared distance between two Vec3D: .
Definition: Vector.h:295
void setFixedParticles(unsigned int n)
Sets a number, n, of particles in the particleHandler as "fixed particles".
Definition: DPMBase.cc:1920
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
Mdouble xBallsVectorScale_
sets the xballs argument vscale (see XBalls/xballs.txt)
Definition: DPMBase.h:1282
virtual void actionsBeforeTimeStep()
A virtual function which allows to define operations to be executed before the new time step...
Definition: DPMBase.cc:1824
ParticleVtkWriter * getVtkWriter() const
Definition: DPMBase.cc:5039
BaseParticle * getSmallestParticleLocal() const
Gets a pointer to the smallest BaseParticle (by interactionRadius) in this ParticleHandler of the loc...
LoggerOutput * loggerOutput
Declaration of the output functions. If the output needs to be redirected, please swap the loggerOutp...
Definition: Logger.cc:168
virtual void computeInternalForce(BaseParticle *, BaseParticle *)
Computes the forces between two particles (internal in the sense that the sum over all these forces i...
Definition: DPMBase.cc:2994
Mdouble timeStep_
Stores the simulation time step.
Definition: DPMBase.h:1219
void setFileType(FileType fileType)
Sets the type of file needed to write into or read from. File::fileType_.
Definition: File.cc:216
void setRandomSeed(unsigned long int new_seed)
This is the seed for the random number generator (note the call to seed_LFG is only required really i...
Definition: RNG.cc:53
WallHandler wallHandler
An object of the class WallHandler. Contains pointers to all the walls created.
Definition: DPMBase.h:1339
void setWriteVTK(FileType f)
RNG random
This is a random generator, often used for setting up the initial conditions etc...
Definition: DPMBase.h:1324
std::enable_if< std::is_scalar< T >::value, void >::type broadcast(T &t, int fromProcessor=0)
Broadcasts a scalar from the root to all other processors.
Definition: MpiContainer.h:429
std::string to_string_padded(unsigned int value)
Pads the number This function tries to pad the number to 4 digits, which is used when you create mult...
Definition: File.cc:44
virtual void computeInternalForces(BaseParticle *)
Computes the internal forces on particle i (internal in the sense that the sum over all these forces ...
Definition: DPMBase.cc:3359
The simulation can be subdivided into Domain's used in parallel code.
Definition: Domain.h:63
std::string to_string(const T &n)
Definition: Helpers.h:227
virtual void hGridRemoveParticle(BaseParticle *obj UNUSED)
Definition: DPMBase.cc:1671
void setPosition(const Vec3D &position)
Sets the position of this BaseInteractable.
Mdouble round(const Mdouble value, unsigned precision)
Definition: Helpers.cc:600
void setMeanVelocity(Vec3D V_mean_goal)
This function will help you set a fixed kinetic energy and mean velocity in your system.
Definition: DPMBase.cc:5050
std::vector< int > get1DParametersFromRunNumber(int size_x) const
This turns a counter into 1 index, which is a useful feature for performing 1D parameter study...
Definition: DPMBase.cc:659
bool writeParticlesVTK_
A flag to turn on/off the vtk writer for particles.
Definition: DPMBase.h:1259
void sync()
Process all pending asynchronous communication requests before continuing.
Definition: MpiContainer.h:148
void setName(const std::string &name)
Allows to set the name of all the files (ene, data, fstat, restart, stat)
Definition: DPMBase.cc:412
virtual void actionsOnRestart()
A virtual function where the users can add extra code which is executed only when the code is restart...
Definition: DPMBase.cc:1643
bool isInCommunicationZone(BaseParticle *particle)
Check if the particle is in the communication zone of the current domain.
Definition: Domain.cc:441
Manages the linked list of PossibleContact.
void writeDataFile()
Definition: DPMBase.cc:2839
virtual void writeEneHeader(std::ostream &os) const
Writes a header with a certain format for ENE file.
Definition: DPMBase.cc:1977
Domain * getCurrentDomain()
Gets the domain assigned to the processor.
void autoNumber()
The autoNumber() function calls three functions: setRunNumber(), readRunNumberFromFile() and incremen...
Definition: DPMBase.cc:528
Vec3D getTotalMomentum() const
JMFT: Return the total momentum of the system, excluding fixed particles.
Definition: DPMBase.cc:1584
std::vector< unsigned > getNumberOfDomains()
returns the number of domains
Definition: DPMBase.cc:5025
bool open()
Checks if the file stream fstream_ has any issues while opening. Alongside, it also increments the ne...
Definition: File.cc:348
void logWriteAndDie(const std::string &module, std::string message)
todo strcmp relies on this, should be changed to more modern version
Definition: DPMBase.cc:73
void set(Vec3D normal, Vec3D point)
Defines a standard wall, given an outward normal vector s.t. normal*x=normal*point for all x of the w...
void removeOldFiles() const
Definition: DPMBase.cc:4205
int getRunNumber() const
This returns the current value of the counter (runNumber_)
Definition: DPMBase.cc:606
Mdouble getZMin() const
If the length of the problem domain in z-direction is ZMax - ZMin, then getZMin() returns ZMin...
Definition: DPMBase.h:611
void setSuperquadricParticlesWriteVTK(bool writeSuperquadricParticlesVTK)
Definition: DPMBase.cc:918
void insertGhostParticle(BaseParticle *P)
This function inserts a particle in the mpi communication boundaries.
Definition: DPMBase.cc:1771
void gatherContactStatistics()
Definition: DPMBase.cc:1858
MERCURY_DEPRECATED File & getFStatFile()
The non const version. Allows to edit the File::fStatFile.
Definition: DPMBase.cc:311
bool isMPIParticle() const
Indicates if this particle is a ghost in the MPI domain.
void setMPIParticle(bool flag)
Flags the mpi particle status.
virtual void outputStatistics()
Definition: DPMBase.cc:1853
File statFile
An instance of class File to handle in- and output into a .stat file.
Definition: DPMBase.h:1390
void boundaryActionsBeforeTimeLoop()
void setTimeStep(Mdouble newDt)
Sets a new value for the simulation time step.
Definition: DPMBase.cc:1195
const std::vector< BaseInteraction * > & getInteractions() const
Returns a list of interactions which belong to this interactable.
static Mdouble getLengthSquared(const Vec3D &a)
Calculates the squared length of a Vec3D: .
Definition: Vector.h:316
void setPeriodicGhostParticle(bool flag)
Flags the status of the particle to be a ghost in periodic boundary or not.
bool compare(std::istream &is, std::string s)
Checks if the next argument in the input stream is a certain string.
Definition: Helpers.cc:860
void removeDuplicatePeriodicParticles()
Removes periodic duplicate Particles.
Definition: DPMBase.cc:4771
bool checkParticleForInteractionLocalPeriodic(const BaseParticle &P)
Definition: DPMBase.cc:4637
bool getSuperquadricParticlesWriteVTK() const
Definition: DPMBase.cc:954
File eneFile
An instance of class File to handle in- and output into a .ene file.
Definition: DPMBase.h:1380
virtual void integrateBeforeForceComputation()
Update particles' and walls' positions and velocities before force computation.
Definition: DPMBase.cc:3108
static int readRunNumberFromFile()
Read the run number or the counter from the counter file (COUNTER_DONOTDEL)
Definition: DPMBase.cc:542
virtual void writeRestartFile()
Stores all the particle data for current save time step to a "restart" file, which is a file simply i...
Definition: DPMBase.cc:2829
std::string xBallsAdditionalArguments_
A string of additional arguments for xballs can be specified (see XBalls/xballs.txt). e.g. "-solidf -v0".
Definition: DPMBase.h:1292
void computeForcesDueToWalls(BaseParticle *, BaseWall *)
Computes the forces on the particles due to the walls (normals are outward normals) ...
Definition: DPMBase.cc:3069
double getXBallsVectorScale() const
Returns the scale of vectors used in xballs.
Definition: DPMBase.cc:1290
Mdouble getTotalMass() const
JMFT: Return the total mass of the system, excluding fixed particles.
Definition: DPMBase.cc:1558
This is a class defining walls.
Definition: InfiniteWall.h:47
virtual bool continueSolve() const
A virtual function for deciding whether to continue the simulation, based on a user-specified criteri...
Definition: DPMBase.cc:1953
void setInteractionDistance(Mdouble interactionDistance)
Sets the interaction distance of the domain handler.
File restartFile
An instance of class File to handle in- and output into a .restart file.
Definition: DPMBase.h:1385
each time-step will be written into/read from separate files numbered consecutively: name_...
Vec3D getMomentum() const
void importParticlesAs(ParticleHandler &particleHandler, InteractionHandler &interactionHandler, const ParticleSpecies *species)
Copies particles, interactions assigning species from a local simulation to a global one...
Definition: DPMBase.cc:4737
Mdouble getNextTime() const
Returns the current simulation time.
Definition: DPMBase.cc:805
void setOrientation(const Quaternion &orientation)
Sets the orientation of this BaseInteractable.
void addTorque(const Vec3D &addTorque)
Adds an amount to the torque on this BaseInteractable.
bool readDataFile(std::string fileName="", unsigned int format=0)
This allows particle data to be reloaded from data files.
Definition: DPMBase.cc:2319
std::ostream & operator<<(std::ostream &os, DPMBase &md)
Definition: DPMBase.cc:90
void updateStatus(std::set< BaseParticle * > &particlesToBeDeleted)
virtual void hGridUpdateMove(BaseParticle *, Mdouble)
Definition: DPMBase.cc:1893
Class that describes a possible contact between two BaseParticle.
bool fileExists(std::string strFilename)
Function to check if a file exists, is used to check if a run has already need done.
Definition: Helpers.cc:503
void integrateAfterForceComputation(double time, double timeStep)
Second step of Velocity Verlet integration.
std::string restartVersion_
Previous versions of MercuryDPM had a different restart file format, the below member variable allows...
Definition: DPMBase.h:1231
void setNumberOfDomains(std::vector< unsigned > direction)
Sets the number of domains in x-,y- and z-direction. Required for parallel computations.
Definition: DPMBase.cc:4958
DomainSplit
Definition: DPMBase.h:912
T * getLastObject()
Gets a pointer to the last Object in this BaseHandler.
Definition: BaseHandler.h:634
virtual void computeForce()
Virtual function that contains the force law between the two objects interacting. ...
void setRestarted(bool newRestartedFlag)
Allows to set the flag stating if the simulation is to be restarted or not.
Definition: DPMBase.cc:1461
void deleteGhostParticles(std::set< BaseParticle * > &particlesToBeDeleted)
Definition: DPMBase.cc:4871
Implementation of a 3D matrix.
Definition: Matrix.h:37
Definition: Vector.h:49
MERCURY_DEPRECATED File & getStatFile()
The non const version. Allows to edit the File::statFile.
Definition: DPMBase.cc:327
T square(const T val)
squares a number
Definition: ExtendedMath.h:104
void setName(const std::string &name)
Sets the file name, e.g. "Name.data".
Definition: File.cc:199
std::string getXBallsAdditionalArguments() const
Returns the additional arguments for xballs.
Definition: DPMBase.cc:1315
virtual void hGridActionsAfterIntegration()
This function has to be called after integrateBeforeForceComputation.
Definition: DPMBase.cc:1907
Mdouble getTimeStep() const
Returns the simulation time step.
Definition: DPMBase.cc:1211
double getXBallsScale() const
Returns the scale of the view in xballs.
Definition: DPMBase.cc:1332
Definition: File.h:80
void createMesh(std::vector< Mdouble > &simulationMin, std::vector< Mdouble > &simulationMax, std::vector< unsigned > &numberOfDomains, bool open)
Creates a Cartesian square mesh in 3D.
BaseParticle * getLargestParticleLocal() const
Gets a pointer to the largest BaseParticle (by interactionRadius) in the ParticleHandler of the local...
Mdouble getRotationalEnergy() const
std::vector< int > get2DParametersFromRunNumber(int size_x, int size_y) const
This turns a counter into 2 indices which is a very useful feature for performing a 2D study...
Definition: DPMBase.cc:687