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 // this header for signal handler
55 #include <csignal>
56 // this header for memset function
57 #include <cstring>
58 
59 #ifdef MERCURY_USE_MPI
60 #include <mpi.h>
61 #endif
62 
63 //This is OMP related
64 #ifdef MERCURY_USE_OMP
65 #include <omp.h>
66 #endif
67 
79 [[noreturn]] void logWriteAndDie(const std::string& module, std::string message)
80 {
81  std::cerr << "A fatal error has occured"
82  << "\n Module :" << module
83  << "\n Message :" << message << std::endl;
84 
85  std::exit(-1);
86 }
87 
96 std::ostream& operator<<(std::ostream& os, DPMBase& md)
97 {
98  md.write(os);
99  return os;
100 }
101 
121 DPMBase::DPMBase(const DPMBase& other) : wallVTKWriter_(other.wallVTKWriter_),
122  interactionVTKWriter_(other.interactionVTKWriter_),
123  boundaryVTKWriter_(other.boundaryVTKWriter_)
124 {
125  setName(other.getName());
126  runNumber_ = other.runNumber_;
129  gravity_ = other.gravity_;
130 /* xMin_ = other.xMin_;
131  xMax_ = other.xMax_;
132  yMin_ = other.yMin_;
133  yMax_ = other.yMax_;
134  zMin_ = other.zMin_;
135  zMax_ = other.zMax_;*/
136  min_ = other.min_;
137  max_ = other.max_;
139  time_ = other.time_;
140  timeStep_ = other.timeStep_;
142  timeMax_ = other.timeMax_;
143  restartVersion_ = other.restartVersion_; //to read new and old restart data
144  restarted_ = other.restarted_; //to see if it was restarted or not
145  append_ = other.append_;
146  rotation_ = other.rotation_;
147  xBallsColourMode_ = other.xBallsColourMode_; // sets the xballs argument cmode (see xballs.txt)
148  xBallsVectorScale_ = other.xBallsVectorScale_; // sets the xballs argument vscale (see xballs.txt)
149  xBallsScale_ = other.xBallsScale_; // sets the xballs argument scale (see xballs.txt)
150  xBallsAdditionalArguments_ = other.xBallsAdditionalArguments_; // std::string where additional xballs argument can be specified (see xballs.txt)
154 
155 //effectively saying "if there exists a CONTACT_LIST_HGRID, copy it, if not, ignore.
156 #ifdef CONTACT_LIST_HGRID
157  possibleContactList=other.possibleContactList;
158 #endif
159  random = other.random;
160 
165  wallHandler.setDPMBase(this);
168  //Initialise the handlers
171 
172  //setting contents equal to the other handlers!
175  cgHandler = other.cgHandler;
176  //cgHandler = other.cgHandler.copy(); //todo
177  //cgHandler.setDPMBase(this);
178  wallHandler = other.wallHandler;
181  vtkWriter_ = other.vtkWriter_;
186 }
187 
193 DPMBase::DPMBase() : wallVTKWriter_(wallHandler), interactionVTKWriter_(interactionHandler), boundaryVTKWriter_(boundaryHandler)
194 {
195  constructor();
196 }
197 
207 {
208  // sofStop function
209  setSoftStop();
210  //constructor();
211  dataFile.getFstream().precision(10);
212  fStatFile.getFstream().precision(10);
213  eneFile.getFstream().precision(10);
214  restartFile.getFstream().precision(
215  std::numeric_limits<double>::digits10); //highly accurate, so the overlap is accurate
216  statFile.getFstream().precision(10);
217  statFile.getFstream().setf(std::ios::left);
218  interactionFile.getFstream().precision(10);
219  name_ = ""; // needs to be user-specified, otherwise checkSettings throws error
220  //by default, the fileType of all files is ONE_FILE. However, by default we don't want an interaction file since it
221  // is very large.
223 
224  runNumber_ = 0;
225 
226  //Decomposition direction for MPI
227  numberOfDomains_ = {1, 1, 1};
228 
229  //Check if MPI is already initialised
230  initialiseMPI();
231 
232  //This sets the maximum number of particles
237  cgHandler.setDPMBase(this);
239  wallHandler.setDPMBase(this);
245 
246  //set defaults for DPMBase parameters
249  setRestarted(false);
250  setGravity(Vec3D(0, 0, 0));
251 
252  //This is the parameter of the numerical part
253  setTime(0);
254  numberOfTimeSteps_ = 0;
255  setTimeMax(0);
256  timeStep_ = 0; // needs to be user-specified, otherwise checkSettings throws error
257  setSaveCount(20);
258 
259  //This sets the default xballs domain
260  min_ = Vec3D(0, 0, 0);
261  max_ = Vec3D(0, 0, 0); // needs to be user-specified, otherwise checkSettings throws error
262 
263  //sets the default write particles data in VTK format flag to false
264  writeParticlesVTK_ = false;
267  vtkWriter_ = nullptr;
268 
269  setName(""); // needs to be user-specified, otherwise checkSettings throws error
270 
271  //Default mode is energy with no scale of the vectors
272  xBallsColourMode_ = 0;
273  xBallsVectorScale_ = -1;
274  xBallsScale_ = -1;
276  setAppend(false);
277 
278  //The default random seed is 0
280 
281  logger(DEBUG, "DPMBase problem constructor finished");
282 
283  readSpeciesFromDataFile_ = false;
284 
286 
287  //Set number of elements to write to the screen if a user wants to output write information to the terminal
288  nToWrite_ = 4;
289 }
290 
296 {
297  delete vtkWriter_;
298 }
299 
304 {
305  return dataFile;
306 }
307 
312 {
313  return eneFile;
314 }
315 
320 {
321  return fStatFile;
322 }
323 
328 {
329  return restartFile;
330 }
331 
336 {
337  return statFile;
338 }
339 
344 {
345  return interactionFile;
346 }
347 
348 
353 {
354  return dataFile;
355 }
356 
360 const File& DPMBase::getEneFile() const
361 {
362  return eneFile;
363 }
364 
369 {
370  return fStatFile;
371 }
372 
377 {
378  return restartFile;
379 }
380 
385 {
386  return statFile;
387 }
391 const File& DPMBase::getInteractionFile() const
393 {
394  return interactionFile;
395 }
396 
397 const std::string& DPMBase::getName() const
398 {
399  return name_;
400 }
401 
406 void DPMBase::setSaveCount(unsigned int saveCount)
407 {
408  dataFile.setSaveCount(saveCount);
409  fStatFile.setSaveCount(saveCount);
410  restartFile.setSaveCount(saveCount);
411  statFile.setSaveCount(saveCount);
412  eneFile.setSaveCount(saveCount);
413  for (auto cg : cgHandler)
414  cg->statFile.setSaveCount(saveCount);
415 }
416 
420 void DPMBase::setName(const std::string& name)
421 {
422  if (NUMBER_OF_PROCESSORS > 1)
423  {
424  name_ = name; // was before this->name_ = name
431  }
432  else
433  {
434  name_ = name; // was before this->name_ = name
435  dataFile.setName(name_ + ".data");
436  fStatFile.setName(name_ + ".fstat");
437  restartFile.setName(name_ + ".restart");
438  statFile.setName(name_ + ".stat");
439  eneFile.setName(name_ + ".ene");
440  interactionFile.setName(name_ + ".interaction");
441  }
442 }
443 
447 void DPMBase::setName(const char* name)
448 {
449  setName(std::string(name));
450 }
451 
458 {
459  dataFile.setFileType(fileType);
460  fStatFile.setFileType(fileType);
461  restartFile.setFileType(fileType);
462  statFile.setFileType(fileType);
463  eneFile.setFileType(fileType);
464 }
465 
470 {
471  dataFile.setCounter(0);
474  statFile.setCounter(0);
475  eneFile.setCounter(0);
482 }
483 
487 void DPMBase::setOpenMode(std::fstream::openmode openMode)
488 {
489  dataFile.setOpenMode(openMode);
490  fStatFile.setOpenMode(openMode);
491  restartFile.setOpenMode(openMode);
492  statFile.setOpenMode(openMode);
493  eneFile.setOpenMode(openMode);
494  interactionFile.setOpenMode(openMode);
495 }
496 
501 {
502  dataFile.close();
503  fStatFile.close();
504  restartFile.close();
505  statFile.close();
506  eneFile.close();
508 }
509 
516 void DPMBase::setLastSavedTimeStep(unsigned int nextSavedTimeStep)
517 {
518  dataFile.setLastSavedTimeStep(nextSavedTimeStep);
519  fStatFile.setLastSavedTimeStep(nextSavedTimeStep);
520  restartFile.setLastSavedTimeStep(nextSavedTimeStep);
521  //statFile.setLastSavedTimeStep(nextSavedTimeStep); //this one is not force-written
522  eneFile.setLastSavedTimeStep(nextSavedTimeStep);
523 }
524 
537 {
539 
540  if (!getRestarted())
541  {
543  }
544 }
545 
551 {
552  int counter;
553 
554  FILE* counter_file;
555  //checking if there exists already a file named "COUNTER_DONOTDEL" which can be opened for
556  //input and output (returns "true" if no such file exists).
557  if ((counter_file = fopen("COUNTER_DONOTDEL", "r+")) == nullptr)
558  {
559  //if a file does not already exist, checks whether a new file can be created
560  //retutns "true" if a new file CANNOT be created
561  if ((counter_file = fopen("COUNTER_DONOTDEL", "w")) == nullptr)
562  {
563  //If true, outputs an error message and ends the program
564  fprintf(stderr, "\n\n\tERROR :: Counter File NOT found, please re-create\n\n");
565  fclose(counter_file);
566  exit(-1);
567  }
568  //alternatively, if a new file CAN be created...
569  else
570  {
571  //starts the new counter file, writing to it the value "1"
572  fprintf(counter_file, "1");
573  fprintf(stderr, "Counter File created\n");
574  fclose(counter_file);
575  return 1;
576  }
577  }
578  //alternatively, if a counter file DOES already exist...
579  else
580  {
581  //...checks if there exists only 1 value in the file (as would be expected from a COUNTER_DONOTDEL file...
582  if (fscanf(counter_file, "%d", &counter) != 1)
583  {
584  //...and if not, returns an error.
585  fprintf(stderr, "\n\n\tERROR :: Counter File found, but something went wrong with reading it\n\n");
586  fclose(counter_file);
587  exit(-1);
588  }
589  else
590  {
591  //...otherwise, if all has been successful, returns the current value of the file!
592  fclose(counter_file);
593  return counter;
594  }
595  }
596 
597 }
598 
603 void DPMBase::setRunNumber(int runNumber)
604 {
605  runNumber_ = runNumber;
606 }
607 
615 {
616  return runNumber_;
617 }
618 
626 {
627  //opening two filestreams - counter_file and counter_file2
628  std::fstream counter_file, counter_file2;
629  //declares an integer, temp_counter
630  int temp_counter;
631  //attempts to open the COUNTER_DONOTDEL text file
632  counter_file.open("COUNTER_DONOTDEL", std::ios::in);
633  //gives error message if file could not be successfully opened and ends the program
634  if (counter_file.fail())
635  {
636  fprintf(stderr, "\n\n\tERROR :: Counter File NOT found, please re-create\n\n");
637  counter_file.close();
638  exit(0);
639  }
640  // if opened successfully, reads in the counter corresponding to the current run number
641  //and stored it in the "temp_counter" variable
642  counter_file >> temp_counter;
643  counter_file.close();
644  //Increments the temp_counter
645  temp_counter++;
646  //opens an output stream to the COUNTER_DONOTDEL file
647  counter_file2.open("COUNTER_DONOTDEL", std::ios::out);
648  if (counter_file2.fail())
649  {
650  fprintf(stderr, "\n\n\tERROR :: Counter File NOT found, please re-create2\n\n");
651  counter_file2.close();
652  exit(0);
653  }
654  //writes the new valuer of the counter to COUNTER_DONOTDEL
655  counter_file2 << temp_counter;
656 
657  counter_file2.close();
658 }
659 
667 std::vector<int> DPMBase::get1DParametersFromRunNumber(int sizeX) const
668 {
669  // Declare a vector of integers capable of storing 2 values
670  std::vector<int> temp(2);
671 
672  // Declare and initialise for the current simulation run number
673  int counter = getRunNumber();
674 
675  // Give studyNum value 0 if study is incomplete, otherwise value > 0
676  int studyNum = (counter-1)/sizeX;
677  counter = counter - sizeX*studyNum;
678 
679  int i = ((counter - 1) % sizeX) + 1;
680  logger(INFO,"StudyNum: % \t Counter: % \t i: %", studyNum, counter, i);
681  temp[0] = studyNum;
682  temp[1] = i;
683 
684  return temp;
685 }
686 
695 std::vector<int> DPMBase::get2DParametersFromRunNumber(int sizeX, int sizeY) const
696 {
697  //declares a vector of integers capable of storing 3 values,
698  std::vector<int> temp(3);
699  //declares and initialises an integer variable named "counter"
700  //with the current counter number, runNumber_
701  int counter = getRunNumber();
702  //calculates the total size of the study, i.e. the number of points
703  //in the 2D parameter space explored
704  int studySize = sizeX * sizeY;
705  //(counter - 1) / studySize gives a fraction comparing the number of runs conducted so far
706  //to the total size of the study, i.e. the total number of runs that need to be performed.
707  //since studyNum is an integer, will declare zero until an adequate number of runs has been performed,
708  //at which point it will equal 1
709  int studyNum = (counter - 1) / studySize;
710 
711  counter = counter - studySize * studyNum;
712  int i = ((counter - 1) % sizeX) + 1;
713  int j = ((counter - i) / sizeX) + 1;
714  logger(INFO,"StudyNum: % \t Counter: % \t i: % \t j: %", studyNum, counter, i, j);
715 
716  temp[0] = studyNum;
717  temp[1] = i;
718  temp[2] = j;
719 
720  return (temp);
721 }
722 
732 std::vector<int> DPMBase::get3DParametersFromRunNumber(int sizeX, int sizeY, int sizeZ) const
733 {
734  //declares a vector of integers capable of storing 4 values,
735  std::vector<int> temp(4);
736  //declares and initialises an integer variable named "counter"
737  //with the current counter number, runNumber_
738  int counter = getRunNumber();
739  //calculates the total size of the study, i.e. the number of points
740  //in the 3D parameter space explored
741  int studySize = sizeX * sizeY * sizeZ;
742  //(counter - 1) / studySize gives a fraction comparing the number of runs conducted so far
743  //to the total size of the study, i.e. the total number of runs that need to be performed.
744  //since studyNum is an integer, will declare zero until an adequate number of runs has been performed,
745  //at which point it will equal 1
746  int studyNum = (counter - 1) / studySize;
747 
748  counter = counter - studySize * studyNum;
749  int i = ((counter-1) % sizeX) + 1;
750  int j = static_cast<int>(std::floor((counter-1)/sizeX)) % sizeY + 1;
751  int k = static_cast<int>(std::floor((counter-1)/(sizeX*sizeY))) % sizeZ + 1;
752  logger(INFO,"StudyNum: % \t Counter: % \t i: % \t j: % \t k: %", studyNum, counter, i, j, k);
753 
754  temp[0] = studyNum;
755  temp[1] = i;
756  temp[2] = j;
757  temp[3] = k;
758 
759  return (temp);
760 }
761 
772 int DPMBase::launchNewRun(const char* name, bool quick UNUSED)
773 {
774  //defines an (empty) stringstream named "com"
775  std::stringstream com("");
776  //adds the name of the code to run (fed in as an argument)
777  //to the "com" string and appends the string with " &"
778  com << name << " &";
779  //converts the stringstream "com" to a standard string, and then
780  //converts this string to a C string
781  //the string is then fed to the "system" function, which will run the named command
782  return system(com.str().c_str());
783 }
784 
796 void DPMBase::solve(int argc, char* argv[])
797 {
798  readArguments(argc, argv);
799  solve();
800 }
801 
806 {
807  return time_;
808 }
809 
814 {
815  return time_ + timeStep_;
816 }
817 
821 unsigned int DPMBase::getNumberOfTimeSteps() const
822 {
823  return numberOfTimeSteps_;
824 }
825 
834 {
835  Mdouble diff = time_ - time;
836  time_ = time;
837  //this sets the interaction timestamp, so each interaction has the right time
838  for (auto i: interactionHandler)
839  {
840  i->setTimeStamp(i->getTimeStamp() - diff);
841  }
842 }
843 
844 
850 void DPMBase::setNToWrite(int nToWrite)
851 {
852  nToWrite_ = nToWrite;
853 }
854 
862 {
863  return nToWrite_;
864 }
865 
871 {
872  if (newTMax >= 0)
873  {
874  timeMax_ = newTMax;
875  }
876  else
877  {
878  logger(ERROR, "Error in setTimeMax, new timeMax=% is not positive", newTMax);
879  }
880 }
881 
886 {
887  return timeMax_;
888 }
895 #ifdef CONTACT_LIST_HGRID
896 PossibleContactList& DPMBase::getPossibleContactList()
897 {
898  return possibleContactList;
899 }
900 #endif
901 
909 {
910  writeWallsVTK_ = writeWallsVTK;
911 }
912 
919 void DPMBase::setWallsWriteVTK(bool writeVTK)
920 {
922 }
923 
925 {
927 }
934 void DPMBase::setParticlesWriteVTK(bool writeParticlesVTK)
935 {
936  writeParticlesVTK_ = writeParticlesVTK;
937  if (writeParticlesVTK_)
938  {
940  }
941  delete vtkWriter_;
943 }
944 
948 void DPMBase::setSuperquadricParticlesWriteVTK(bool writeParticlesVTK)
949 {
950  writeSuperquadricParticlesVTK_ = writeParticlesVTK;
952  {
953  writeParticlesVTK_ = false;
954  }
955  delete vtkWriter_;
957 }
958 
966 {
967  return writeWallsVTK_;
968 }
969 
977 {
978  return writeParticlesVTK_;
979 }
980 
985 {
987 }
988 
1002 {
1003  if (newXMin <= getXMax())
1004  {
1005  min_.x() = newXMin;
1006  }
1007  else
1008  {
1009  logger(WARN, "Warning in setXMin(%): xMax=%", newXMin, getXMax());
1010  }
1011 }
1012 
1026 {
1027  if (newYMin <= getYMax())
1028  {
1029  min_.y() = newYMin;
1030  }
1031  else
1032  {
1033  logger(WARN, "Warning in setYMin(%): yMax=%", newYMin, getYMax());
1034  }
1035 }
1036 
1050 {
1051 
1052  if (newZMin <= getZMax())
1053  {
1054  min_.z() = newZMin;
1055  }
1056  else
1057  {
1058  logger(WARN, "Warning in setZMin(%): zMax=%", newZMin, getZMax());
1059  }
1060 
1061 }
1062 
1073 void DPMBase::setMax(const Vec3D& newMax)
1074 {
1075  if (min_.x() > newMax.x() ||
1076  min_.y() > newMax.y() ||
1077  min_.z() > newMax.z())
1078  {
1079  logger(WARN, "Warning in setMax: upper bound is smaller"
1080  " than lower bound. (%,%,%) > (%,%,%)",
1081  min_.x(), min_.y(), min_.z(), newMax.x(), newMax.y(), newMax.z());
1082  }
1083  else
1084  {
1085  max_ = newMax;
1086  }
1087 }
1088 
1089 void DPMBase::setDomain(const Vec3D& min, const Vec3D& max)
1090 {
1091 
1092  logger.assert_debug(min.X <= max.X, "lower x-bound (%) is larger than upper x-bound (%)", min.X, max.X);
1093  logger.assert_debug(min.Y <= max.Y, "lower x-bound (%) is larger than upper x-bound (%)", min.Y, max.Y);
1094  logger.assert_debug(min.Z <= max.Z, "lower x-bound (%) is larger than upper x-bound (%)", min.Z, max.Z);
1095  min_ = min;
1096  max_ = max;
1097 }
1098 
1109 void DPMBase::setMin(const Vec3D& newMin)
1110 {
1111  if (max_.x() < newMin.x() ||
1112  max_.y() < newMin.y() ||
1113  max_.z() < newMin.z())
1114  {
1115  logger(WARN, "Warning in setMin: lower bound is larger"
1116  " than upper bound. (%,%,%) < (%,%,%)",
1117  max_.x(), max_.y(), max_.z(), newMin.x(), newMin.y(), newMin.z());
1118  }
1119  else
1120  {
1121  min_ = newMin;
1122  }
1123 }
1124 
1129 void DPMBase::setMin(const Mdouble xMin, const Mdouble yMin, const Mdouble zMin)
1130 {
1131  setMin(Vec3D(xMin, yMin, zMin));
1132 }
1133 
1138 void DPMBase::setMax(const Mdouble xMax, const Mdouble yMax, const Mdouble zMax)
1139 {
1140  setMax(Vec3D(xMax, yMax, zMax));
1141 }
1142 
1157 {
1158 
1159  if (newXMax >= getXMin())
1160  {
1161  max_.x() = newXMax;
1162  }
1163  else
1164  {
1165  logger(WARN, "Warning in setXMax(%): xMax=%", newXMax, getXMin());
1166  }
1167 
1168 }
1169 
1183 {
1184 
1185  if (newYMax >= getYMin())
1186  {
1187  max_.y() = newYMax;
1188  }
1189  else
1190  {
1191  logger(WARN, "Warning in setYMax(%): yMax=%", newYMax, getYMin());
1192  }
1193 
1194 }
1195 
1209 {
1210  if (newZMax >= getZMin())
1211  {
1212  max_.z() = newZMax;
1213  }
1214  else
1215  {
1216  logger(WARN, "Warning in setZMax(%): zMax=%", newZMax, getZMin());
1217  }
1218 }
1219 
1226 {
1227  if (timeStep > 0.0)
1228  {
1229  timeStep_ = timeStep;
1230  }
1231  else
1232  {
1233  logger(ERROR, "Error in setTimeStep: new timeStep % is not positive", timeStep);
1234  }
1235 }
1236 
1242 {
1243  return timeStep_;
1244 }
1245 
1246 
1247 /* Allows user to set the number of omp threads */
1248 void DPMBase::setNumberOfOMPThreads(int numberOfOMPThreads)
1249 {
1250  logger.assert_always(numberOfOMPThreads > 0, "Number of OMP threads must be positive");
1251  numberOfOMPThreads_ = numberOfOMPThreads;
1252 
1253 #ifdef MERCURY_USE_OMP
1254  if(numberOfOMPThreads > omp_get_max_threads()) {
1255  logger(INFO, "Number of omp threads set to the maximum number of threads allowed: %",
1256  omp_get_max_threads());
1257  numberOfOMPThreads_ = numberOfOMPThreads = omp_get_max_threads();
1258  }
1259 #pragma omp parallel num_threads(getNumberOfOMPThreads())
1260  {
1261  if (omp_get_thread_num()==0)
1262  logger(INFO, "Using % of % omp threads; testing thread",
1263  omp_get_num_threads(), omp_get_max_threads(), Flusher::NO_FLUSH);
1264  }
1265 #pragma omp parallel num_threads(getNumberOfOMPThreads())
1266  {
1267  logger(INFO) " %", std::to_string(omp_get_thread_num(), Flusher::NO_FLUSH);
1268  }
1269  logger(INFO,"");
1270 
1271 #else
1272  logger(WARN, "You are setting the number of omp threads to %, but OMP is not turned on", getNumberOfOMPThreads());
1273 #endif
1274 }
1275 
1276 /* Returns the number of omp threads */
1278 {
1279  //logger.assert_debug(numberOfOMPThreads_,"You need to set the number of OMP threads");
1280  return numberOfOMPThreads_;
1281 }
1282 
1292 {
1293  xBallsColourMode_ = newCMode;
1294 }
1295 
1302 {
1303  return xBallsColourMode_;
1304 }
1305 
1311 void DPMBase::setXBallsVectorScale(double newVScale)
1312 {
1313  xBallsVectorScale_ = newVScale;
1314 }
1315 
1322 {
1323  return xBallsVectorScale_;
1324 }
1325 
1338 void DPMBase::setXBallsAdditionalArguments(std::string xBallsAdditionalArguments)
1339 {
1340  xBallsAdditionalArguments_ = xBallsAdditionalArguments;
1341 }
1342 
1347 {
1349 }
1350 
1355 {
1356  xBallsScale_ = newScale;
1357 }
1358 
1364 {
1365  return xBallsScale_;
1366 }
1367 
1374 void DPMBase::setGravity(Vec3D newGravity)
1375 {
1376  gravity_ = newGravity;
1377 }
1378 
1383 {
1384  return gravity_;
1385 }
1386 
1394 void DPMBase::setDimension(unsigned int newDim)
1395 {
1396  setSystemDimensions(newDim);
1397  setParticleDimensions(newDim);
1398 }
1399 
1408 void DPMBase::setSystemDimensions(unsigned int newDim)
1409 {
1410  if (newDim >= 1 && newDim <= 3)
1411  systemDimensions_ = newDim;
1412  else
1413  {
1414  logger(ERROR, "Error in setSystemDimensions; newDim % is not 1, 2 or 3", newDim);
1415  }
1416 }
1417 
1421 unsigned int DPMBase::getSystemDimensions() const
1422 {
1423  return systemDimensions_;
1424 }
1425 
1439 void DPMBase::setParticleDimensions(unsigned int particleDimensions)
1440 {
1441  if (particleDimensions >= 1 && particleDimensions <= 3)
1442  {
1443  particleDimensions_ = particleDimensions;
1445  }
1446  else
1447  {
1448  logger(WARN, "Error in setParticleDimensions; particleDimensions % is not 1, 2 or 3, setting it to "
1449  "systemDimension now (%)", particleDimensions, systemDimensions_);
1451  }
1452 }
1453 
1459 {
1460  return particleDimensions_;
1461 }
1462 
1466 std::string DPMBase::getRestartVersion() const
1467 {
1468  return restartVersion_;
1469 }
1470 
1475 void DPMBase::setRestartVersion(std::string newRV)
1476 {
1477  restartVersion_ = newRV;
1478 }
1479 
1485 {
1486  return restarted_;
1487 }
1488 
1492 void DPMBase::setRestarted(bool newRestartedFlag)
1493 {
1494  restarted_ = newRestartedFlag;
1495  //setAppend(new_);
1496 }
1497 
1502 {
1503  return append_;
1504 }
1505 
1513 void DPMBase::setAppend(bool newAppendFlag)
1514 {
1515  append_ = newAppendFlag;
1516 }
1517 
1522 {
1523  Mdouble elasticEnergy = 0.0;
1524  // JMFT: Note that we do count the elastic energy of fixed particles here.
1525  for (const BaseInteraction* c : interactionHandler)
1526  {
1527  elasticEnergy += c->getElasticEnergy();
1528  }
1529  return elasticEnergy;
1530 }
1531 
1536 {
1537  Mdouble kineticEnergy = 0;
1538  for (const BaseParticle* const p : particleHandler)
1539  {
1540  if (!(p->isFixed()))
1541  {
1542  kineticEnergy += .5 * p->getMass() * p->getVelocity().getLengthSquared();
1543  }
1544  }
1545  return kineticEnergy;
1546 }
1547 
1553 {
1554  Mdouble gravitationalEnergy = 0;
1555  for (const BaseParticle* const p : particleHandler)
1556  {
1557  // Don't consider fixed particles. 'Fixed' particles aren't necessarily
1558  // stationary; it just means their position is prescribed.
1559  if (!(p->isFixed()))
1560  {
1561  gravitationalEnergy += p->getMass() * Vec3D::dot((-getGravity()), p->getPosition());
1562  }
1563  }
1564  return gravitationalEnergy;
1565 }
1566 
1569 {
1570  Mdouble ene_rot = 0;
1571  for (std::vector<BaseParticle*>::const_iterator it = particleHandler.begin(); it != particleHandler.end(); ++it)
1572  {
1573  // See above.
1574  if (!(*it)->isFixed())
1575  {
1576  // ene_rot += .5 * (*it)->getInertia() * (*it)->getAngularVelocity().getLengthSquared();
1577  }
1578  }
1579  return ene_rot;
1580 }
1581 
1584 }
1585 
1590 {
1591  /*
1592  double mass_sum = 0;
1593  for (std::vector<BaseParticle*>::const_iterator it = particleHandler.begin(); it != particleHandler.end(); ++it)
1594  if (!(*it)->isFixed())
1595  mass_sum += (*it)->getMass();
1596  return mass_sum;
1597  */
1598  return particleHandler.getMass();
1599 }
1600 
1606 {
1608 }
1609 
1616 {
1617  return particleHandler.getMomentum();
1618  /*
1619  Vec3D total_momentum = Vec3D(0,0,0);
1620  for (std::vector<BaseParticle*>::const_iterator it = particleHandler.begin(); it != particleHandler.end(); ++it)
1621  if (!(*it)->isFixed())
1622  total_momentum += (*it)->getMass() * (*it)->getVelocity();
1623  return total_momentum;
1624  */
1625 }
1626 
1630 /* JMFT: This information is placed on the last entry on each line of the .data
1631  * file. That space is, in general, reserved for 'additional' information.
1632  */
1634 {
1635 // return p.getSpecies()->getId(); // was getIndex()
1636  return p.getInfo();
1637 }
1638 
1653 {
1654  return (pI != pJ && pI->isInContactWith(pJ));
1655 }
1656 
1661 {
1662 }
1663 
1668 {
1669 }
1670 
1675 {
1676 }
1677 
1682 {
1683 }
1684 
1689 {
1690 }
1691 
1696 {
1697 }
1698 
1703 {
1704 }
1705 
1710 {
1711  return true;
1712 }
1713 
1724 {
1725 #ifdef MERCURY_USE_MPI
1726  //If only one core is used (i.e. domainHandler is empty) then the result is always true
1727  if (domainHandler.getSize() == 0)
1728  {
1729  return true;
1730  }
1731  //Get the current domain
1733 
1734  //Check if the particle is in the current domain
1735  if(domain->containsParticle(P))
1736  {
1737  //When adding a particle inside the domain, this should always be true
1738  return true;
1739  }
1740  else
1741  {
1742  //MPI particles that are inserted in the communication zone should still be inserted
1743  return (P->isMPIParticle());
1744  }
1745 #else
1746  return false;
1747 #endif
1748 }
1749 
1756 {
1757 
1758  bool insideCommunicationZone = false;
1759 #ifdef MERCURY_USE_MPI
1760  MPIContainer& communicator = MPIContainer::Instance();
1761 
1762  //Check for the current domain if the particle is within the communication domain
1763  int val = domainHandler.getCurrentDomain()->isInCommunicationZone(particle);
1764 
1765  //The root gathers all results
1766  int *list = nullptr;
1767  if (PROCESSOR_ID == 0)
1768  {
1769  list = new int [NUMBER_OF_PROCESSORS];
1770  }
1771  communicator.gather(val,list);
1772 
1773  //Compute the global value
1774  //if on any processor the val is true, we have to do the communcation step
1776  int result = 0;
1777  if (PROCESSOR_ID == 0)
1778  {
1779  for (int i = 0; i< NUMBER_OF_PROCESSORS; i++)
1780  {
1781  if (list[i] == 1)
1782  {
1783  result = 1;
1784  break;
1785  }
1786  }
1787  }
1788 
1789  //The root now tells the other processors what the global value for the interaction is
1790  communicator.broadcast(result);
1791 
1792  //Convert the result back to bool
1793  insideCommunicationZone = result;
1794 #endif
1795  return insideCommunicationZone;
1796 }
1797 
1803 {
1804 #ifdef MERCURY_USE_MPI
1805  //mpi particles only exist when there is more than one domain
1806  if (domainHandler.getSize() > 0)
1807  {
1808  //Add the particle to the mpi domain
1810  }
1811 
1812  //If periodic boundaries are present..
1813  if (periodicBoundaryHandler.getSize() > 0)
1814  {
1816  }
1817 #endif
1818 }
1819 
1829 {
1830 #ifdef MERCURY_USE_MPI
1831  if (NUMBER_OF_PROCESSORS == 1) { return; }
1832 
1833  //Check if the interactionRadius of the BaseParticle is larger than given in the domain
1836  {
1837  logger(VERBOSE,"Processor % | Updating mpi grid. Old interactionDistance: %, new interactionDistance %.",
1839 
1840  //Update the interactionDistance in the domain and periodicBoundaryHandler
1843 
1844  //Find new ghost particless
1847  }
1848 #endif
1849 }
1850 
1851 
1856 {
1857 }
1858 
1863 {
1864 }
1865 
1870 {
1871 }
1872 
1877 {
1879 }
1880 
1885 {
1886  //cgHandler.evaluate();
1887 }
1888 
1890 {
1892  {
1893  c->gatherContactStatistics();
1894  }
1895 }
1896 
1900 void DPMBase::gatherContactStatistics(unsigned int index1 UNUSED, int index2 UNUSED, Vec3D Contact UNUSED,
1901  Mdouble delta UNUSED, Mdouble ctheta UNUSED, Mdouble fdotn UNUSED,
1902  Mdouble fdott UNUSED, Vec3D P1_P2_normal_ UNUSED, Vec3D P1_P2_tangential UNUSED)
1903 {
1904 }
1905 
1910 {
1911 }
1912 
1917 {
1918  cgHandler.finish();
1919 }
1920 
1925 {
1926 }
1927 
1932 {
1933 }
1934 
1939 {
1940 }
1941 
1951 void DPMBase::setFixedParticles(unsigned int n)
1952 {
1953  for (unsigned int i = 0; i < std::min(particleHandler.getSize(), n); i++)
1955 }
1956 
1962 {
1963 #ifdef MERCURY_USE_MPI
1964  MPIContainer& communicator = MPIContainer::Instance();
1965  if (communicator.getProcessorID() == 0)
1966  {
1967 #endif
1968  logger(INFO, "t=%3.6, tmax=%3.6", getTime(), getTimeMax());
1969 #ifdef MERCURY_USE_MPI
1970  }
1971 #endif
1972 }
1973 
1982 {
1983  return continueFlag_ != 0;
1984 }
1985 
1990 {
1991 }
1992 
2005 void DPMBase::writeEneHeader(std::ostream& os) const
2006 {
2007  //only write if we don't restart
2008  if (getAppend())
2009  return;
2010 
2012 
2015  long width = os.precision() + 6;
2016  os << std::setw(width)
2017  << "time " << std::setw(width)
2018  << "gravitEnergy " << std::setw(width) //gravitational potential energy
2019  << "traKineticEnergy " << std::setw(width) //translational kinetic energy
2020  << "rotKineticEnergy " << std::setw(width) //rotational kE
2021  << "elasticEnergy " << std::setw(width)
2022  << "centerOfMassX " << std::setw(width)
2023  << "centerOfMassY " << std::setw(width)
2024  << "centerOfMassZ\n";
2025 }
2026 
2034 void DPMBase::writeFstatHeader(std::ostream& os) const
2035 {
2036 
2037  // line #1: time, volume fraction
2038  // line #2: wall box: wx0, wy0, wz0, wx1, wy1, wz1
2039  // line #3: radii-min-max & moments: rad_min, rad_max, r1, r2, r3, r4
2040  os << "#"
2041  << " " << getTime()
2042  << " " << 1 //marker that these are fstat files with contact point instead of center point
2043  << '\n';
2044  os << "#"
2045  << " " << getXMin()
2046  << " " << getYMin()
2047  << " " << getZMin()
2048  << " " << getXMax()
2049  << " " << getYMax()
2050  << " " << getZMax()
2051  << '\n';
2052  os << "#"
2053  << " ";
2054 
2055  if (!(particleHandler.getSmallestParticleLocal() == nullptr))
2056  {
2058  }
2059  else
2060  {
2061  os << std::numeric_limits<double>::quiet_NaN();
2062  }
2063  os << " ";
2064  if (!(particleHandler.getLargestParticleLocal() == nullptr))
2065  {
2067  }
2068  else
2069  {
2070  os << std::numeric_limits<double>::quiet_NaN();
2071  }
2072 
2073  os << " " << 0
2074  << " " << 0
2075  << " " << 0
2076  << " " << 0
2077  << '\n';
2078  //B: write data
2080  {
2081  c->writeToFStat(os, getTime());
2082  }
2083  //os << std::flush;
2084 }
2085 
2095 void DPMBase::writeEneTimeStep(std::ostream& os) const
2096 {
2099  writeEneHeader(os);
2100 
2101  const Mdouble m = particleHandler.getMass();
2103  //Ensure the numbers fit into a constant width column: for this we need the precision given by the operating system,
2104  //plus a few extra characters for characters like a minus and scientific notation.
2105  const static int width = os.precision() + 6;
2106  os << std::setw(width) << getTime()
2107  << " " << std::setw(width) << -Vec3D::dot(getGravity(), com)
2108  << " " << std::setw(width) << particleHandler.getKineticEnergy()
2109  << " " << std::setw(width) << particleHandler.getRotationalEnergy()
2110  << " " << std::setw(width) << getElasticEnergy()
2111  // we need to write x, y and z coordinates separately, otherwise the width of the columns is incorrect
2112  << " " << std::setw(width)
2113  << (m == 0 ? constants::NaN : com.X / m) //set to nan because 0/0 implementation in gcc and clang differs
2114  << " " << std::setw(width) << (m == 0 ? constants::NaN : com.Y / m)
2115  << " " << std::setw(width) << (m == 0 ? constants::NaN : com.Z / m)
2116  << std::endl;
2117 }
2118 
2120 {
2121  static bool writeWall = true;
2122  if (getWallsWriteVTK() == FileType::ONE_FILE && writeWall)
2123  {
2125  writeWall=false;
2129  } // else do nothing
2130 
2132  {
2133  vtkWriter_->writeVTK();
2134  } // else do nothing
2135 
2137  {
2139  }
2140 
2142  {
2144  }
2145 
2146  //only write once
2147  bool writePython = getParticlesWriteVTK() || getWallsWriteVTK() != FileType::NO_FILE ||
2149  if (writePython && getTime() == 0)
2150  {
2152  }
2153 }
2154 
2156 {
2157 #ifdef MERCURY_USE_MPI
2158  if (PROCESSOR_ID == 0)
2159  {
2160  logger(INFO, "Writing python script for paraview visualisation");
2161 #else
2162  logger(INFO, "Writing python script for paraview visualisation");
2163 #endif
2164 
2165  std::string script = "#script to visualise the output of data2pvd of MercuryDPM in paraview.\n"
2166  "#usage: change the path below to your own path, open paraview\n"
2167  "#Tools->Python Shell->Run Script->VisualisationScript.py\n"
2168  "#or run paraview --script=VisualisationScript.py \n"
2169  "\n"
2170  "from paraview.simple import *\n"
2171  "import os\n"
2172  "import glob\n"
2173  "os.chdir('" + helpers::getPath() + "')\n\n";
2174 
2175 #ifdef MERCURY_USE_MPI
2176  for (int i = 0; i < NUMBER_OF_PROCESSORS; i++)
2177  {
2178 #endif
2179  if (getParticlesWriteVTK())
2180  {
2181  script += "#Load data in any order\n";
2182 #ifdef MERCURY_USE_MPI
2183  if (NUMBER_OF_PROCESSORS > 1)
2184  {
2185  script += "Data = glob.glob('./" + getName() + "Processor_" + std::to_string(i) + "_Particle_*.vtu')\n";
2186  }
2187  else
2188  {
2189  script += "Data = glob.glob('./" + getName() + "Particle_*.vtu')\n";
2190  }
2191 #else
2192  script += "Data = glob.glob('./" + getName() + "Particle_*.vtu')\n";
2193 #endif
2194  script += "\n"
2195  "#Find the maximum timestep\n"
2196  "maxTime = 0\n"
2197  "for fileName in Data:\n"
2198  "\ttokens1 = fileName.split('.')\n"
2199  "\ttokens2 = tokens1[1].split('_')\n"
2200  "\tif int(tokens2[-1]) > maxTime:\n"
2201  "\t\tmaxTime = int(tokens2[-1])\n"
2202  "print str(maxTime)\n"
2203  "\n"
2204  "#Create correct order of time steps\n"
2205  "DataSorted = []\n"
2206  "for x in range(0,maxTime+1):\n";
2207 #ifdef MERCURY_USE_MPI
2208  if (NUMBER_OF_PROCESSORS > 1)
2209  {
2210  script += "\tDataSorted.append('./" + getName() + "Processor_" + std::to_string(i) + "_Particle_' + " + "str(x)" + " + '.vtu')\n";
2211  }
2212  else
2213  {
2214  script += "\tDataSorted.append('./" + getName() + "Particle_' + " + "str(x)" + " + '.vtu')\n";
2215  }
2216 #else
2217  script += "\tDataSorted.append('./" + getName() + "Particle_' + " + "str(x)" + " + '.vtu')\n";
2218 #endif
2219 
2220  script += "\n"
2221  "#Load the data and visualise it in paraview\n"
2222  "particles = XMLUnstructuredGridReader(FileName=DataSorted)\n"
2223  "glyphP = Glyph(particles)\n"
2224  "glyphP.GlyphType = 'Sphere'\n"
2225  "glyphP.Scalars = 'Radius'\n"
2226  "glyphP.Vectors = 'None'\n"
2227  "glyphP.ScaleMode = 'scalar'\n"
2228  "glyphP.ScaleFactor = 2\n"
2229  "glyphP.GlyphMode = 'All Points'\n"
2230  "Show(glyphP)\n\n";
2231  }
2233  {
2234  script += "walls = XMLUnstructuredGridReader(FileName=glob.glob('./"
2235  + getName() + "Wall_*.vtu'))\n"
2236  "Show(walls)\n\n";
2237  }
2240  {
2241  script += "interactions = XMLUnstructuredGridReader(FileName=glob.glob('./"
2242  + getName() + "Interaction_*.vtu'))\n"
2243  "glyphI = Glyph(interactions)\n"
2244  "glyphI.GlyphType = 'Sphere'\n"
2245  "glyphI.Scalars = 'Cylinder'\n"
2246  "glyphI.Vectors = 'None'\n"
2247  "glyphI.ScaleMode = 'scalar'\n"
2248  "glyphI.ScaleFactor = 10\n" //5 times too large
2249  "glyphI.GlyphMode = 'All Points'\n"
2250  "Show(glyphI)\n\n";
2251  }
2252 #ifdef MERCURY_USE_MPI
2253  } // end of loop over number of processors
2254 #endif
2255  script += "Render()\n"
2256  "ResetCamera()\n";
2257 
2258  helpers::writeToFile(getName() + ".py", script);
2259 #ifdef MERCURY_USE_MPI
2260  } // end of communicator is root statement
2261 #endif
2262 }
2263 
2264 
2268 void DPMBase::outputXBallsData(std::ostream& os) const
2269 {
2270 
2271 
2272  //Set the correct formation based of dimension if the formation is not specified by the user
2273 
2274  unsigned int format;
2275  switch (getSystemDimensions())
2276  {
2277  case 2:
2278  format = 8;
2279  break;
2280  case 3:
2281  format = 14;
2282  break;
2283  default:
2284  logger(ERROR, "Unknown system dimension");
2285  }
2286 
2287  unsigned int numberOfParticles = particleHandler.getNumberOfRealObjectsLocal();
2288 
2289  // This outputs the location of walls and how many particles there are to file this is required by the xballs plotting
2290  if (format != 14) // dim = 1 or 2
2291  {
2292  os << numberOfParticles
2293  << " " << getTime()
2294  << " " << getXMin()
2295  << " " << getYMin()
2296  << " " << getXMax()
2297  << " " << getYMax()
2298  << " " << std::endl;
2299  }
2300  else
2301  {
2302  //dim==3
2303  os << numberOfParticles
2304  << " " << getTime()
2305  << " " << getXMin()
2306  << " " << getYMin()
2307  << " " << getZMin()
2308  << " " << getXMax()
2309  << " " << getYMax()
2310  << " " << getZMax()
2311  << " " << std::endl;
2312  }
2313 
2314  // This outputs the particle data
2315  for (unsigned int i = 0; i < particleHandler.getSize(); i++)
2316  {
2317 #ifdef MERCURY_USE_MPI
2319  {
2320  outputXBallsDataParticle(i, format, os);
2321  }
2322 #else
2323  outputXBallsDataParticle(i, format, os);
2324 #endif
2325  }
2326 #ifdef DEBUG_OUTPUT
2327  logger(DEBUG, "Have output the properties of the problem to disk ");
2328 #endif
2329 }
2330 
2346 bool DPMBase::readDataFile(std::string fileName, unsigned int format)
2347 {
2348  //default value: dataFile.getFullName()
2349  if (!fileName.compare(""))
2350  fileName = dataFile.getFullName();
2351 
2352  std::string oldFileName = dataFile.getName();
2353  unsigned oldCounter = dataFile.getCounter();
2354  //Updates the name of the data file to the user-input from the argument.
2355  dataFile.setName(fileName);
2356  //opens a filestream of the input type
2357  dataFile.open(std::fstream::in);
2358  //Checks if the file has been successfully opened...
2359  if (!dataFile.getFstream().is_open() || dataFile.getFstream().bad())
2360  {
2361  //...and if not, ends the function and returns "false"
2362  logger(WARN, "Loading data file % failed.", fileName);
2363  return false;
2364  }
2365 
2366  //retrieves and saves the "FileType" of the file
2367  FileType fileTypeData = dataFile.getFileType();
2369  readNextDataFile(format);
2370  dataFile.setFileType(fileTypeData);
2371  dataFile.close();
2372  dataFile.setName(oldFileName);
2373  dataFile.setCounter(oldCounter);
2374  return true;
2375 }
2376 
2382 bool DPMBase::readParAndIniFiles(const std::string fileName)
2383 {
2384  //Opens the par.ini file
2385  std::fstream file;
2386  file.open(fileName, std::fstream::in);
2387  if (!file.is_open() || file.bad())
2388  {
2389  //std::cout << "Loading par.ini file " << filename << " failed" << std::endl;
2390  return false;
2391  }
2392 
2393  Mdouble doubleValue;
2394  int integerValue;
2395 
2396  // inputfile par.ini
2397  // line 1 =============================================================
2398  // Example: 1 1 0
2399  // 1: integer (0|1) switches from non-periodic to periodic
2400  // integer (5|6) does 2D integration only (y-coordinates fixed)
2401  // and switches from non-periodic to periodic
2402  // integer (11) uses a quarter system with circular b.c.
2403  file >> integerValue;
2404  //~ std::cout << "11" << integerValue << std::endl;
2405  if (integerValue == 0)
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  w0.set(Vec3D(0, 0, -1), Vec3D(0, 0, getZMin()));
2419  w0.set(Vec3D(0, 0, 1), Vec3D(0, 0, getZMax()));
2421  }
2422  else if (integerValue == 1)
2423  {
2424  PeriodicBoundary b0;
2425  b0.set(Vec3D(1, 0, 0), getXMin(), getXMax());
2427  b0.set(Vec3D(0, 1, 0), getYMin(), getYMax());
2429  b0.set(Vec3D(0, 0, 1), getZMin(), getZMax());
2431  }
2432  else if (integerValue == 5)
2433  {
2434  InfiniteWall w0;
2436  w0.set(Vec3D(-1, 0, 0), Vec3D(-getXMin(), 0, 0));
2438  w0.set(Vec3D(1, 0, 0), Vec3D(getXMax(), 0, 0));
2440  w0.set(Vec3D(0, -1, 0), Vec3D(0, -getYMin(), 0));
2442  w0.set(Vec3D(0, 1, 0), Vec3D(0, getYMax(), 0));
2444 
2445  }
2446  else if (integerValue == 6)
2447  {
2448  PeriodicBoundary b0;
2449  b0.set(Vec3D(1, 0, 0), getXMin(), getXMax());
2451  b0.set(Vec3D(0, 1, 0), getYMin(), getYMax());
2453  b0.set(Vec3D(0, 0, 1), getZMin(), getZMax());
2455  }
2456  else
2457  {
2458  logger(ERROR, "Error in par.ini: line 1, value 1 is ", integerValue);
2459  }
2460 
2461  // 2: integer (0|1) dont use | use the search pattern for linked cells
2462  file >> integerValue; //ignore
2463 
2464  // 3: real - gravity in z direction: positive points upwards
2465  file >> doubleValue;
2466  setGravity(Vec3D(0.0, 0.0, doubleValue));
2467 
2468  // line 2 =============================================================
2469  // Example: -10000 .5e-2
2470  // 1: time end of simulation - (negative resets start time to zero
2471  // and uses -time as end-time)
2472  file >> doubleValue;
2473  if (doubleValue < 0)
2474  setTime(0);
2475  setTimeMax(fabs(doubleValue));
2476 
2477  // 2: time-step of simulation
2478  file >> doubleValue;
2479  setTimeStep(doubleValue);
2480 
2481  // line 3 =============================================================
2482  // Example: 1e-1 100
2483  file >> doubleValue;
2484  if (doubleValue >= 0)
2485  {
2486  // 1: time-step for output on time-series protocoll file -> "ene"
2487  unsigned int savecount = static_cast<unsigned int>(round(doubleValue / getTimeStep()));
2488  setSaveCount(savecount);
2489 
2490  // 2: time-step for output on film (coordinate) file -> "c3d"
2491  // (fstat-output is coupled to c3d-output time-step)
2492  file >> doubleValue;
2493  savecount = static_cast<unsigned int>(round(doubleValue / getTimeStep()));
2494  dataFile.setSaveCount(savecount);
2495  fStatFile.setSaveCount(savecount);
2496  }
2497  else
2498  {
2499  // or: ---------------------------------------------------------------
2500  // 1: negative number is multiplied to the previous log-output-time
2501  // 2: requires initial log-output time
2502  // 3: negative number is multiplied to the previous film-output-time
2503  // 4: requires initial film-output time
2504  logger(ERROR, "Error in par.ini: line 3, value 1 is ", doubleValue);
2505  }
2506 
2507  // line 4 =============================================================
2508  // Example: 2000 1e5 1e3 1e2
2509  // 1: particle density (mass=4/3*constants::pi*density*rad^3)
2510  file >> doubleValue;
2511 
2512  //clear species handler
2516 
2518  S->setDensity(doubleValue);
2519 
2520  // 2: linear spring constant
2521  file >> doubleValue;
2522  S->setStiffness(doubleValue);
2523 
2524  // 3: linear dashpot constant
2525  file >> doubleValue;
2526  S->setDissipation(doubleValue);
2527 
2528  // 4: background damping dashpot constant
2529  file >> doubleValue;
2530  if (doubleValue != 0.0)
2531  logger(WARN, "Warning in par.ini: ignored background damping %", doubleValue);
2532 
2533  // line 5 =============================================================
2534  // Example: 0 0
2535  // 1: growth rate: d(radius) = xgrow * dt
2536  file >> doubleValue;
2537  if (doubleValue != 0.0)
2538  logger(WARN, "Warning in par.ini: ignored growth rate ", doubleValue);
2539 
2540  // 2: target volume_fraction
2541  file >> doubleValue;
2542  if (doubleValue != 0.0)
2543  logger(WARN, "Warning in par.ini: ignored target volume_fraction %", doubleValue);
2544 
2545  file.close();
2546  //std::cout << "Loaded par.ini file " << filename << std::endl;
2547  return true;
2548 }
2549 
2564 {
2566  {
2567  while (true)// This true corresponds to the if s
2568  {
2569  dataFile.open();
2570  //check if file exists and contains data
2571  int N;
2572  dataFile.getFstream() >> N;
2573  if (dataFile.getFstream().eof() || dataFile.getFstream().peek() == -1)
2574  {
2575  logger(WARN, "file % not found", dataFile.getName());
2576  return false;
2577  }
2578  //check if tmin condition is satisfied
2579  Mdouble t;
2580  dataFile.getFstream() >> t;
2581  if (t > tMin)
2582  {
2583  //set_file_counter(get_file_counter()-1);
2584  return true;
2585  }
2586  if (verbose)
2587  logger(VERBOSE, "Jumping file counter: %", dataFile.getCounter());
2588  }
2589  }
2590  return true;
2591 }
2592 
2600 bool DPMBase::readNextDataFile(unsigned int format)
2601 {
2602  dataFile.open(std::fstream::in);
2603  //logger(INFO,"Reading %",dataFile.getFullName());
2604  //fStatFile.open();
2605  //Set the correct format based of dimension if the format is not explicitly specified by the user
2606  if (format == 0)
2607  {
2608  //checking the dimensionality of the system
2610  switch (getSystemDimensions())
2611  {
2612  case 1:
2613  //if 2D, sets format 8 (data file has 8 columns for 2D)
2614  case 2:
2615  format = 8;
2616  break;
2617  case 3:
2618  //if 3D, sets format 14 (data file has 14 columns for 3D)
2619  format = 14;
2620  break;
2621  }
2622  //end case
2623  }
2624  //end if
2625 
2626  // read in the particle number (as a double)
2627  double doubleN = -1;
2628  dataFile.getFstream() >> doubleN;
2629 
2630  // If N cannot be read, we reached the end of the file
2631  if (doubleN == -1) return false;
2632 
2633  // converting N to an integer; skipping the line if there is a problem (this happens when there is a corrupt data file)
2634  unsigned N = doubleN;
2635  while (doubleN != N) {
2636  std::string dummy;
2637  getline(dataFile.getFstream(),dummy,'\n');
2638  logger(WARN,"Skipping bad line in data file: % %",doubleN, dummy);
2639  dataFile.getFstream() >> doubleN;
2640  N = doubleN;
2641  }
2642 
2643  //store the parameters you want to preserve:
2644  const size_t nHistory = std::min(N,particleHandler.getSize());
2645  std::vector<const ParticleSpecies*> species(nHistory);
2646  std::vector<bool> fix(nHistory);
2647  for (size_t i=0; i<nHistory; ++i) {
2649  species[i] = p->getSpecies();
2650  fix[i] = p->isFixed();
2651  //store from reading the restart file which particles are fixed and which species each particle has
2652  }
2653 
2654  BaseParticle* p;
2655  if (particleHandler.getSize() == 0) {
2656  logger.assert_always(speciesHandler.getSize()>0,"readData: species needs to be set first");
2658  p->unfix();
2659  } else {
2660  p = particleHandler.getObject(0)->copy();
2661  }
2662 
2663  //empty the particle handler
2666 
2667  //now fill it again with the read-in particles
2668  std::stringstream line;
2670  Mdouble radius;
2671  Vec3D position, velocity, angle, angularVelocity;
2672  size_t indSpecies;
2673  //read all other data available for the time step
2674  if (format==7) {
2675  line >> time_ >> min_.x() >> min_.y() >> min_.z() >> max_.x() >> max_.y() >> max_.z();
2676  for (size_t i = 0; i < N; ++i) {
2678  line >> position.X >> position.Z >> position.Y >> velocity.X >> velocity.Z
2679  >> velocity.Y >> radius
2680  >> indSpecies;
2681  p->setPosition(position);
2682  p->setVelocity(velocity);
2683  p->setOrientation({1, 0, 0, 0});
2684  p->setAngularVelocity({0.0, 0.0, 0.0});
2685  p->setRadius(radius);
2687  p->setSpecies(speciesHandler.getObject(indSpecies));
2688  else if (i < nHistory)
2689  p->setSpecies(species[i]);
2690  if (i < nHistory && fix[i])
2691  p->fixParticle();
2693  p->unfix();
2694  }
2695  } else if (format==8) {
2696  line >> time_ >> min_.x() >> min_.y() >> max_.x() >> max_.y();
2697  min_.z() = 0.0; max_.z() = 0.0;//For 2d functions we define Z to be of no interest
2698  for (size_t i = 0; i < N; ++i) {
2700  line >> position.X >> position.Y >> velocity.X >> velocity.Y >> radius >> angle.Z >> angularVelocity.Z >> indSpecies;
2701  Quaternion q;
2702  q.setAngleZ(angle.Z);
2703  p->setPosition(position);
2704  p->setVelocity(velocity);
2705  p->setOrientation(q);
2706  p->setAngularVelocity(-angularVelocity);
2707  p->setRadius(radius);
2709  else if (i < nHistory) p->setSpecies(species[i]);
2710  if (i < nHistory && fix[i]) p->fixParticle();
2712  p->unfix();
2713  } //end for all particles
2714  } else if (format==14) {
2715  //This is a 3D format_
2716  //@TODO: Check bounds or get rid of this function
2717  line >> time_ >> min_.x() >> min_.y() >> min_.z() >> max_.x() >> max_.y() >> max_.z();
2718  for (size_t i = 0; i < N; ++i) {
2720  line >> position >> velocity >> radius >> angle >> angularVelocity >> indSpecies;
2721 #ifdef MERCURY_USE_MPI
2722  //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
2723  bool isMPIParticle = false;
2724  bool isPeriodicGhostParticle = false;
2725  if (NUMBER_OF_PROCESSORS > 1)
2726  {
2728  line >> isMPIParticle >> isPeriodicGhostParticle;
2729  }
2730 #endif
2731  p->setPosition(position);
2732  p->setVelocity(velocity);
2733  p->setOrientationViaEuler(angle);
2734  p->setAngularVelocity(angularVelocity);
2735  p->setRadius(radius);
2737  if (indSpecies<speciesHandler.getSize()) {
2738  p->setSpecies(speciesHandler.getObject(indSpecies));
2739  } else {
2740  logger(WARN, "Read in bad species data; species is not set");
2741  }
2742  } else if (i < nHistory)
2743  p->setSpecies(species[i]);
2744  if (i < nHistory && fix[i])
2745  p->fixParticle();
2746 #ifdef MERCURY_USE_MPI
2748  {
2749  p->setMPIParticle(isMPIParticle);
2750  p->setPeriodicGhostParticle(isPeriodicGhostParticle);
2751  }
2752 #endif
2754  p->unfix();
2755  } //end read into existing particles logger(INFO, "read % particles", particleHandler.getNumberOfObjects());
2756  } else if (format==15) {
2757  line >> time_ >> min_.x() >> min_.y() >> min_.z() >> max_.z() >> max_.y() >> max_.z();
2758  for (size_t i = 0; i < N; ++i) {
2760  line >> position >> velocity >> radius >> angle >> angularVelocity >> indSpecies >> indSpecies;
2761  Quaternion q;
2762  q.setEuler(angle);
2763  p->setPosition(position);
2764  p->setVelocity(velocity);
2765  p->setOrientation(q);
2766  p->setAngularVelocity(angularVelocity);
2767  p->setRadius(radius);
2769  else if (i < nHistory) p->setSpecies(species[i]);
2770  if (i < nHistory && fix[i]) p->fixParticle();
2772  } //end for all particles
2773  } //end if format
2774 
2776  return true;
2777 }
2778 
2780 {
2781  fStatFile.open(std::fstream::in);
2782  std::string line;
2783  std::fstream& in = fStatFile.getFstream();
2784  // read the first three lines
2785  getline(in, line);
2786  getline(in, line);
2787  getline(in, line);
2788  Mdouble time;
2789  unsigned int indexP;
2790  int indexI; //could be negative
2791  unsigned counter = 0;
2793  while ((in.peek() != -1) && (in.peek() != '#'))
2794  {
2795  /* # 1: time
2796  # 2: particle Number i
2797  # 3: contact partner j (particles >= 0, walls < 0)
2798  # 4: x-position \
2799  # 5: y-position > of the contact point (I hope)
2800  # 6: z-position /
2801  # 7: delta = overlap at the contact
2802  # 8: ctheta = length of the tangential spring
2803  # 9: P1_P2_normal force |f^n|
2804  # 10: remaining (tangential) force |f^t|=|f-f^n|
2805  # 11-13: P1_P2_normal unit vector nx, ny, nz
2806  # 14-16: tangential unit vector tx, ty, tz
2807  */
2808  in >> time >> indexP >> indexI;
2809  BaseParticle* P = particleHandler.getObject(indexP);
2810  BaseInteraction* C;
2811  if (indexI >= 0)
2812  {
2813  //read only one of the two fstat lines reported
2814  if (indexI >= indexP)
2815  {
2816  // particle pair contact
2817  BaseParticle* I = particleHandler.getObject(static_cast<const unsigned int>(indexI));
2819  C->setFStatData(in, P, I);
2820  // skip next line
2821  //in.ignore(256, '\n');
2822  }
2823  }
2824  else
2825  {
2826  // wall-particle contact
2827  while (wallHandler.getNumberOfObjects() <= -indexI - 1)
2828  {
2830  logger(WARN, "Added new wall because .fstat file indicates contact with wall % that doesn't exist",
2831  -indexI - 1);
2832  }
2833  BaseWall* I = wallHandler.getObject(static_cast<const unsigned int>(-indexI - 1));
2835  C->setFStatData(in, P, I);
2836  }
2837  counter++;
2838  //skip the rest of the line, to allow additional output and to ignore the second time a particle-particle contact is printed
2839  in.ignore(256, '\n');
2840  }
2841  //logger(INFO,"read % contacts at t = % (N=%)",counter,timeStamp,interactionHandler.getNumberOfObjects());
2842  //interactionHandler.write(std::cout);
2843  //logger(INFO,"normal % %",interactionHandler.getObject(0)->getNormal(),interactionHandler.getObject(0)->getContactPoint());
2844  //logger(INFO,"normal % %",interactionHandler.getLastObject()->getNormal(),interactionHandler.getLastObject()->getContactPoint());
2845 }
2846 
2855 {
2857  {
2858  //logger(DEBUG, "Writing restart file %th time step",getNumberOfTimeSteps());
2860  restartFile.close();
2861  }
2862 }
2863 
2865 {
2867  {
2869  dataFile.close();
2870  }
2871 }
2872 
2874 {
2876  {
2877  //If the file type is "multiple files, writes a header for each individual files. If not, only writes for the first time step
2879  eneFile.close();
2880  }
2881 }
2882 
2884 {
2886  {
2888  //fStatFile.getFstream().ignore(2000,'\t');
2889  fStatFile.close();
2890  }
2891 }
2892 
2900  logger.assert_always(speciesHandler.getSize()>0,"There needs to be at least one species");
2902  SphericalParticle p(s);
2904  Mdouble r = cbrt(getTotalVolume())/N;
2905  b.set(p, 100, getMin(), getMax(), {0, 0, 0}, {0, 0, 0});
2906  b.insertParticles(this);
2907  logger(INFO,"Inserted % particles",particleHandler.getSize());
2908  //setTimeMax(0);
2909  //solve();
2910 }
2911 
2912 
2922 {
2923  //Assuming a filename corresponding to "restartFile" has already been established,
2924  //opens an input filestream to the relevant file
2925  if (restartFile.open(std::fstream::in))
2926  {
2927  //reads the input stream line-by-line
2928  read(restartFile.getFstream(), opt);
2929  logger(INFO, "Loaded restart file %", restartFile.getFullName());
2930  restartFile.close();
2931  //sets the flag such that other functions can tell that the file has been restarted
2932  //e.g. does not run "setUpInitialConditions" or add headers to the .ene files etc.
2933  setRestarted(true);
2934  return true;
2935  }
2936  else /* if the file could not be opened */
2937  {
2938  logger(INFO, "% could not be loaded.", restartFile.getFullName());
2939  return false;
2940  }
2941 }
2942 
2950 int DPMBase::readRestartFile(std::string fileName, ReadOptions opt)
2951 {
2952  //add ".restart" if necessary
2953  if (fileName.find(".restart") == std::string::npos)
2954  {
2955  fileName.append(".restart");
2956  }
2957 
2958  //If padded or numbered files are used, we need to extract the counter and remove it from the filename
2959  //First find the last point in the restart name
2960  unsigned int pos = fileName.find('.');
2961  while (fileName.find('.', pos + 1) != std::string::npos)
2962  {
2963  pos = fileName.find('.', pos + 1);
2964  }
2965  //If the next char after the last . is a digit we are using numbered files
2966  std::string counter;
2967  if (isdigit(fileName[pos + 1]))
2968  {
2969  for (int i = pos + 1; i < fileName.length(); i++)
2970  {
2971  counter.push_back(fileName[i]);
2972  }
2973  //Set counter in restart file
2974  restartFile.setCounter(std::stoi(counter));
2975  logger(INFO, "Counter: %", std::stoi(counter));
2976  }
2977 
2978 #ifdef MERCURY_USE_MPI
2979  //Correct for the processor number
2980  if (NUMBER_OF_PROCESSORS > 1 && !helpers::fileExists(fileName))
2981  {
2982  //Modify file name
2983  const unsigned int length = fileName.length();
2984  if (isdigit(fileName[pos + 1]))
2985  {
2986  for (int i = pos + 1; i < length + 1; i++)
2987  {
2988  fileName.pop_back();
2989  }
2990  }
2991  fileName.append(std::to_string(PROCESSOR_ID));
2992  if (counter.size() > 0)
2993  {
2994  fileName.append(".");
2995  fileName.append(counter);
2996  }
2997  }
2998 #endif
2999 
3000  restartFile.setName(fileName);
3001 
3002  logger(INFO, "Restarting from %", fileName);
3003  return readRestartFile(opt);
3004 }
3005 
3020 {
3021  //Does not compute forces if particles are fixed
3022  //this is necessary because the rough bottom allows overlapping fixed particles
3023  if (P1->isFixed() && P2->isFixed())
3024  {
3025  return;
3026  }
3027 //Ensures that interactions between the "ghost" particles used to implement periodic behaviour
3028  //are not included in calculations
3029  //i.e. ends the function if both particles are "ghosts".
3030  if ((P1->getPeriodicFromParticle() != nullptr) && (P2->getPeriodicFromParticle() != nullptr))
3031  {
3032  return;
3033  }
3034 //if statement below ensures that the PI has the lower id than PJ
3035  BaseParticle* PI, * PJ;
3036  if (P1->getId() > P2->getId())
3037  {
3038  PI = P2;
3039  PJ = P1;
3040  }
3041  else
3042  {
3043  PI = P1;
3044  PJ = P2;
3045  }
3046  //checks if the two particles are interacting
3047  //("getInteractionWith" returns the relevant pointer if PI and PJ are interacting,
3048  //zero if not)
3049  //if statement above ensures that the PI has the lower id than PJ
3052  if (i!= nullptr) {
3053  //calculates the force corresponding to the interaction
3054  i->computeForce();
3055 
3056  //Applies the relevant calculated forces to PI and PJ
3057  PI->addForce(i->getForce());
3058  PJ->addForce(-i->getForce());
3059 
3060  //checks if particle rotation is turned on...
3061  if (getRotation()) {
3062  //...and, if so, performs equivalent calculations for the torque as were
3063  //performed for the force.
3064  PI->addTorque(i->getTorque() - Vec3D::cross(PI->getPosition() - i->getContactPoint(), i->getForce()));
3065  PJ->addTorque(-i->getTorque() + Vec3D::cross(PJ->getPosition() - i->getContactPoint(), i->getForce()));
3066  }
3067  }
3068 }
3069 
3075 {
3076  //Checks that the current particle is not "fixed"
3077  //and hence infinitely massive!
3078  if (!CI->isFixed())
3079  {
3080  // Applying the force due to gravity (F = m.g)
3081  CI->addForce(getGravity() * CI->getMass());
3082  // Still calls this in compute External Forces.
3083  // computeForcesDueToWalls(CI);
3084  }
3085 }
3086 
3095 {
3096  //No need to compute interactions between periodic particle images and walls
3097  if (pI->getPeriodicFromParticle() != nullptr)
3098  return;
3099 
3100  //Checks if the particle is interacting with the current wall
3103  if (i!=nullptr) {
3104  //...calculates the forces between the two objects...
3105  i->computeForce();
3106 
3107  //...and applies them to each of the two objects (wall and particle).
3108  pI->addForce(i->getForce());
3109  w->addForce(-i->getForce());
3110 
3111  //If the rotation flag is on, also applies the relevant torques
3112  //(getRotation() returns a boolean).
3113  if (getRotation()) // getRotation() returns a boolean.
3114  {
3115  pI->addTorque(i->getTorque() - Vec3D::cross(pI->getPosition() - i->getContactPoint(), i->getForce()));
3117  w->addTorque(-i->getTorque() + Vec3D::cross(w->getPosition() - i->getContactPoint(), i->getForce()));
3118  }
3119  }
3120 }
3121 
3134 {
3135  //cycling through all particles, p, in the particleHandler
3136  //for_each(particleHandler.begin(), particleHandler.end(), [this](BaseParticle* p)
3137  //for (BaseParticle* p : particleHandler) {
3138 
3139  #pragma omp parallel for num_threads(getNumberOfOMPThreads()) //schedule(dynamic)
3140  for (int k = 0; k < particleHandler.getSize(); ++k) {
3142 #ifdef MERCURY_USE_MPI
3143  //MPI particles are not integrated, they are purely ghost particles and get their new velocity and position from an MPI update
3144  if (!(p->isMPIParticle() || p->isPeriodicGhostParticle()))
3145  {
3147  }
3148 #else
3149  //using the particle p's internal "integrateBeforeForceComputation" function
3150  //to update the relevant parameters concerning the particle's position and motion
3152 #endif
3153  }
3154  //});
3155  //cycling through all walls, w, in the wallHandler
3156  //for_each(wallHandler.begin(), wallHandler.end(), [this](BaseWall* w)
3157  //for (BaseWall* w : wallHandler) {
3158  #pragma omp parallel for num_threads(getNumberOfOMPThreads()) //schedule(dynamic)
3159  for (int k = 0; k < wallHandler.getSize(); k++) {
3160  BaseWall *w = wallHandler.getObject(k);
3161  //using the wall's internal "integrateBeforeForceComputation" function
3162  //to update the relevant parameters concerning its position and motion
3164  }
3165  //});
3166 }
3167 
3178 {
3179 
3180  //Cycling over all boundaries within the system...
3181  for (BaseBoundary* b : boundaryHandler)
3182  {
3183  //check all boundaries...
3184  b->checkBoundaryAfterParticlesMove(particleHandler);
3185 
3186 
3187 #ifdef MERCURY_USE_MPI
3188  //When ghost particles are deleted by deletion boundaries they need to be removed
3189  //from their communication lists to avoid segfaults
3190  if (NUMBER_OF_PROCESSORS > 1)
3191  {
3192  //Flush deleted particles from mpi communication zones
3193  getCurrentDomain()->flushParticles(b->getParticlesToBeDeleted());
3195  periodicBoundaryHandler.flushParticles(b->getParticlesToBeDeleted());
3197  }
3198 
3199  //Delete particles that were in communication zone
3200  for (auto p_it = b->getParticlesToBeDeleted().begin(); p_it != b->getParticlesToBeDeleted().end(); p_it++)
3201  {
3202  particleHandler.removeGhostObject((*p_it)->getIndex());
3203  }
3204 #endif
3205  }
3206 }
3207 
3220 {
3221  //cycling through all particles, p, in the particleHandler
3222  //for_each(particleHandler.begin(), particleHandler.end(), [this](BaseParticle* p){
3223  #pragma omp parallel for num_threads(getNumberOfOMPThreads()) //schedule(dynamic)
3224  for (int k = 0; k < particleHandler.getSize(); ++k) {
3226 #ifdef MERCURY_USE_MPI
3227  //MPI particles do not require integration - they are updated by the communication step
3228  if (!(p->isMPIParticle() || p->isPeriodicGhostParticle()))
3229  {
3231  }
3232 #else
3233  //using the particle p's internal "integrateAfterForceComputation" function
3234  //to update the relevant parameters concerning the particle's position and motion
3236 #endif
3237  }
3238  //});
3239  //cycling through all walls, w, in the wallHandler
3240  //for_each(wallHandler.begin(), wallHandler.end(), [this](BaseWall* w){
3241  #pragma omp parallel for num_threads(getNumberOfOMPThreads()) //schedule(dynamic)
3242  for (int k = 0; k < wallHandler.getSize(); k++) {
3243  BaseWall *w = wallHandler.getObject(k);
3244  //using the wall's internal "integrateAfterForceComputation" function
3245  //to update the relevant parameters concerning its position and motion
3247  }
3248  //});
3249 }
3250 
3254 //void DPMBase::statisticsFromRestartData(const char *name)
3255 //{
3256 // ///todo{Check this whole function}
3257 // //This function loads all MD data
3258 // readRestartFile();
3259 //
3260 // //This creates the file statistics will be saved to
3261 // std::stringstream ss("");
3262 // ss << name << ".stat";
3263 // statFile.setName(ss.str());
3264 // statFile.setOpenMode(std::fstream::out);
3265 // statFile.open();
3266 //
3267 // // Sets up the initial conditions for the simulation
3268 // // setupInitialConditions();
3269 // // Setup the previous position arrays and mass of each particle.
3270 // computeParticleMasses();
3271 // // Other routines required to jump-start the simulation
3272 // actionsBeforeTimeLoop();
3273 // initialiseStatistics();
3274 // hGridActionsBeforeTimeLoop();
3275 // writeEneHeader(eneFile.getFstream());
3276 //
3277 // while (readDataFile(dataFile.getName().c_str()))
3278 // {
3279 // hGridActionsBeforeTimeLoop();
3280 // actionsBeforeTimeStep();
3281 // checkAndDuplicatePeriodicParticles();
3282 // hGridActionsBeforeTimeStep();
3287 // computeAllForces();
3288 // removeDuplicatePeriodicParticles();
3289 // actionsAfterTimeStep();
3290 // writeEneTimeStep(eneFile.getFstream());
3291 // std::cout << std::setprecision(6) << getTime() << std::endl;
3292 // }
3293 //
3294 // dataFile.close();
3295 // statFile.close();
3296 //}
3303 {
3304  //Resetting all forces on both particles and walls to zero
3305  #pragma omp parallel num_threads(getNumberOfOMPThreads())
3306  {
3307  #pragma omp for
3308  for (int k = 0; k < particleHandler.getSize(); ++k) {
3310  }
3311  #pragma omp for
3312  for (int k = 0; k < wallHandler.getSize(); k++) {
3314  }
3315  }
3316  logger(DEBUG,"All forces set to zero");
3317 
3318  // for omp simulations, reset the newObjects_ variable (used for reduction)
3320 
3321  // compute all internal and external forces; for omp simulations, this can be done in parallel
3322  #pragma omp parallel num_threads(getNumberOfOMPThreads())
3323  {
3324  //logger(INFO, "Number of omp threads = %", getNumberOfOMPThreads());
3326  #pragma omp for schedule(dynamic)
3327  for (int k = 0; k < particleHandler.getSize(); ++k) {
3329  //computing both internal forces (e.g. due to collisions)
3330  //and external forces (e.g. gravity)
3331  //(compute internal forces compares the current particle p
3332  //with all others in the handler!)
3334  // body forces
3336  }
3337 
3338  // wall-forces
3339  #pragma omp for schedule(dynamic)
3340  for (int k = 0; k < wallHandler.getSize(); k++) {
3341  BaseWall *w = wallHandler.getObject(k);
3342  computeWallForces(w);
3343  }
3344 
3345  }
3346 
3347 #ifdef CONTACT_LIST_HGRID
3348  PossibleContact* Curr=possibleContactList.getFirstPossibleContact();
3349  while(Curr)
3350  {
3351  computeInternalForces(Curr->getP1(),Curr->getP2());
3352  Curr=Curr->getNext();
3353  }
3354 #endif
3355 
3356  // Check wall forces
3357  #pragma omp for schedule(dynamic)
3358  for (int k = 0; k < wallHandler.getNumberOfObjects(); k++) {
3359  BaseWall *w = wallHandler.getObject(k);
3361  }
3362 
3364 
3365  // for omp simulations, sum up all forces and add all newObjects_ (needed since both are using reduction)
3366  #ifdef MERCURY_USE_OMP
3367  if (getNumberOfOMPThreads()>1) {
3369  }
3370  //Resetting all forces on both particles and walls to zero
3371  #pragma omp parallel num_threads(getNumberOfOMPThreads())
3372  {
3373  #pragma omp for
3374  for (int k = 0; k < particleHandler.getSize(); k++) {
3376  }
3377  #pragma omp for
3378  for (int k = 0; k < wallHandler.getSize(); k++) {
3380  } //end reset forces loop
3381  }
3382  #endif
3383 
3384  //end outer loop over contacts.
3385 }
3386 
3394 {
3395  for (auto it = particleHandler.begin(); (*it) != i; ++it)
3396  {
3397  computeInternalForce(i, *it);
3398  }
3399 }
3400 
3401 
3409 void DPMBase::write(std::ostream& os, bool writeAllParticles) const
3410 {
3411  os << "MercuryDPM " << getVersion();
3412  //which outputs basic information regarding the various files (.data, .fstat etc. etc.)
3413  //only writes the run number if it is different from 0
3414  if (runNumber_ != 0)
3415  os << " runNumber " << runNumber_;
3416  os << " name " << name_;
3417  os << " revision " << getSVNRevision();
3418  os << " repository " << getSVNURL() << '\n';
3419  os << "dataFile " << dataFile << '\n';
3420  os << "fStatFile " << fStatFile << '\n';
3421  os << "eneFile " << eneFile << '\n';
3422  os << "restartFile " << restartFile << '\n';
3423  os << "statFile " << statFile << '\n';
3424  os << "interactionFile " << interactionFile << '\n';
3425  //Outputs the "domain" corresponding to the system for
3426  //use with XBalls, as well as other information regarding the system as a whole
3427  os << "xMin " << getXMin()
3428  << " xMax " << getXMax()
3429  << " yMin " << getYMin()
3430  << " yMax " << getYMax()
3431  << " zMin " << getZMin()
3432  << " zMax " << getZMax() << '\n'
3433  << "timeStep " << getTimeStep()
3434  << " time " << getTime()
3435  << " ntimeSteps " << numberOfTimeSteps_
3436  << " timeMax " << getTimeMax() << '\n'
3437  << "systemDimensions " << getSystemDimensions()
3438  << " particleDimensions " << getParticleDimensions()
3439  << " gravity " << getGravity();
3440  os << " writeVTK " << writeParticlesVTK_
3441  << " " << writeWallsVTK_
3442  << " " << interactionHandler.getWriteVTK()
3443  << " " << (vtkWriter_?vtkWriter_->getFileCounter():0)
3444  << " " << wallVTKWriter_.getFileCounter()
3446  << " " << boundaryVTKWriter_.getFileCounter();
3447  os << " random ";
3448  random.write(os);
3449 #ifdef MERCURY_USE_OMP
3450  //Write number of OMP threads
3451  if(getNumberOfOMPThreads() > 1) {
3452  os << " numberOfOMPThreads " << getNumberOfOMPThreads();
3453  }
3454 #endif
3455 #ifdef MERCURY_USE_MPI
3456  //Check if we are dealing with multiple cores
3457  if (NUMBER_OF_PROCESSORS > 1 )
3458  {
3459  os << " numberOfProcessors " << NUMBER_OF_PROCESSORS
3461  }
3462 #endif
3463  //only write xBallsArguments if they are nonzero
3465  os << " xBallsArguments " << getXBallsAdditionalArguments();
3466  os << '\n';
3467  //writes all species (including mixed species) to an output stream
3468 
3469  speciesHandler.write(os);
3470 
3471  //outputs the number of walls in the system
3472  os << "Walls " << wallHandler.getNumberOfObjects() << std::endl;
3473  if (writeAllParticles || wallHandler.getSize() < 9) {
3474  for (BaseWall* w : wallHandler)
3475  os << (*w) << std::endl;
3476  } else {
3477  for (int i=0; i<2; ++i)
3478  os << *wallHandler.getObject(i) << std::endl;
3479  os << "...\n";
3480  }
3481 
3482  //outputs the number of boundaries in the system
3483  os << "Boundaries " << boundaryHandler.getNumberOfObjects() << std::endl;
3484  if (writeAllParticles || boundaryHandler.getSize() < 9)
3485  {
3486  for (BaseBoundary* b: boundaryHandler)
3487  os << (*b) << std::endl;
3488  }
3489  else
3490  {
3491  for (int i = 0; i < 2; ++i)
3492  os << *boundaryHandler.getObject(i) << std::endl;
3493  os << "...\n";
3494  }
3495 
3496  if (writeAllParticles || particleHandler.getSize() < getNToWrite())
3497  {
3498  //if the "writeAllParticles" bool == true, or there are fewer than 4 particles
3499  //calls the particleHandler version of the "write" function and also
3500  //outputs to file all relevant particle information for all particles in the system
3501  particleHandler.write(os);
3502  }
3503  else
3504  {
3505  //otherwise, only prints out limited information
3506  os << "Particles " << particleHandler.getSize() << '\n';
3507  for (unsigned int i = 0; i < getNToWrite(); i++)
3508  os << *(particleHandler.getObject(i)) << '\n';
3509  os << "..." << '\n';
3510  }
3511  // Similarly, print out interaction details (all of them, or up to nToWrite of them)
3512  if (writeAllParticles || interactionHandler.getNumberOfObjects() < getNToWrite())
3513  {
3515  }
3516  else
3517  {
3518  os << "Interactions " << interactionHandler.getNumberOfObjects() << '\n';
3519  for (unsigned int i = 0; i < getNToWrite(); i++)
3520  os << *(interactionHandler.getObject(i)) << '\n';
3521  os << "..." << '\n';
3522  }
3523 }
3524 
3530 void DPMBase::read(std::istream& is, ReadOptions opt)
3531 {
3532 #ifdef MERCURY_USE_MPI
3533  int previousNumberOfProcessors;
3534 #endif
3535  //Declares...
3536  std::string dummy;
3537  //...and reads in a dummy variable from the start of the stream "is"
3538  is >> dummy;
3539  //compare the string read in to the phrase "restart_version" to see if the stream corresponds
3540  //to a restart file (all restart files begin with this phrase)
3541  //if both strings match, strcmp(dummy.c_str(), "restart_version") returns 0 (here read as "false")
3542  if (dummy != "restart_version" && dummy != "MercuryDPM")
3543  {
3544  //If the strings do not match, if statement is fulfilled and the error logged
3545  //Note: only very old files did not have a restart_version
3546  logger(FATAL, "Error in DPMBase::read(is): this is not a valid restart file");
3547  }
3548  else
3549  {
3550  //reads in the restart version (earlier versions of Mercury possess different file formats!)
3551  is >> restartVersion_;
3552  //checking which version the current data file corresponds to, and reads the data in
3553  //accordingly
3554  if (restartVersion_ == "1.0" || restartVersion_ == "0.14")
3555  {
3556  //reads in and saves the relevant values from the data file to the current instance of DPMBase
3557  std::stringstream line;
3558 
3559  // Store path (if restart file is nonlocal)
3560  auto slash = restartFile.getName().rfind('/');
3561  std::string path;
3562  if (slash != std::string::npos)
3563  {
3564  path = restartFile.getName().substr(0, slash + 1);
3565  }
3566  if (!path.empty())
3567  {
3568  logger(INFO, "Adding path information (%) to file names", path);
3569  }
3570 
3571  //line 1
3573  //discards the whitespace (ws) at the start of the stream
3574  line >> std::ws;
3575  //uses the "peek" function to access the stream's first
3576  //non-whitespace character, and check if it is an "r"
3577  if (line.peek() == 'r')
3578  //if so, reads in the current run number
3579  line >> dummy >> runNumber_;
3580  //In either case, then uses the "Files" version of the read function
3581  //to read in the rest of the relevant information.
3582  line >> dummy >> name_;
3583  setName(name_);
3584 
3585  //Read line 2-7 (definition of i/o files)
3587  line >> dummy >> dataFile;
3589  line >> dummy >> fStatFile;
3591  line >> dummy >> eneFile;
3593  line >> dummy >> restartFile;
3595  line >> dummy >> statFile;
3596 
3597  // Add the file path from the restart file to the file names
3598  dataFile.setName(path + dataFile.getName());
3599  fStatFile.setName(path + fStatFile.getName());
3600  eneFile.setName(path + eneFile.getName());
3601  restartFile.setName(path + restartFile.getName());
3602  statFile.setName(path + statFile.getName());
3603 
3604  // Get current position
3605  //check if the next line starts with 'interactionFile'; otherwise, skip interaction
3606  if (helpers::compare(is, "interactionFile"))
3607  {
3609  line >> interactionFile;
3610  interactionFile.setName(path + interactionFile.getName());
3611  }
3612 
3614  line >> dummy >> min_.x()
3615  >> dummy >> max_.x()
3616  >> dummy >> min_.y()
3617  >> dummy >> max_.y()
3618  >> dummy >> min_.z()
3619  >> dummy >> max_.z();
3620 
3622  line >> dummy >> timeStep_
3623  >> dummy >> time_
3624  >> dummy >> numberOfTimeSteps_
3625  >> dummy >> timeMax_;
3626 
3628  line >> dummy >> systemDimensions_
3629  >> dummy >> particleDimensions_
3630  >> dummy >> gravity_;
3631 
3632  line >> dummy;
3633  if (!dummy.compare("writeVTK"))
3634  {
3635  FileType writeInteractionsVTK = FileType::NO_FILE;
3636  unsigned particlesCounter, wallCounter, interactionCounter;
3637  bool writeBoundaryVTK;
3638  line >> writeParticlesVTK_ >> writeWallsVTK_ >> writeInteractionsVTK >> writeBoundaryVTK >> particlesCounter >> wallCounter >> interactionCounter;
3639  line.clear();//because the number of arguments in writeVTK has changed
3640  line >> dummy;
3643  interactionHandler.setWriteVTK(writeInteractionsVTK);
3644  boundaryHandler.setWriteVTK(writeBoundaryVTK);
3645  vtkWriter_->setFileCounter(particlesCounter);
3646  wallVTKWriter_.setFileCounter(particlesCounter);
3647  interactionVTKWriter_.setFileCounter(particlesCounter);
3648  boundaryVTKWriter_.setFileCounter(particlesCounter);
3649  }
3650  if (!dummy.compare("random"))
3651  {
3652  random.read(line);
3653  line >> dummy;
3654  }
3655 
3656 #ifdef MERCURY_USE_OMP
3657  //Read the number of OMP threads
3658  if (!dummy.compare("numberOfOMPThreads")) {
3659  int numberOfOMPThreads;
3660  line >> numberOfOMPThreads;
3661  setNumberOfOMPThreads(numberOfOMPThreads);
3662  //logger(INFO," Check the number of OMP threads = % ", getNumberOfOMPThreads());
3663  }
3664 #endif
3665 #ifdef MERCURY_USE_MPI
3666  if (!dummy.compare("numberOfProcessors"))
3667  {
3668  line >> previousNumberOfProcessors
3669  >> dummy >> numberOfDomains_[Direction::XAXIS]
3672  }
3673  else
3674  {
3675  logger(INFO,"Reading a serial restart file");
3676  //numberOfDomains_ = {1,1,1};
3677  }
3678 #endif
3679  if (!dummy.compare("xBallsArguments")) {
3681  setXBallsAdditionalArguments(line.str());
3682  }
3683 
3684  speciesHandler.read(is);
3685 
3686 #ifdef MERCURY_USE_MPI
3687  //Initialise MPI structures and perform domain decomposition
3688  decompose();
3689 #endif
3690 
3691  //reading in the various relevant handlers
3692  unsigned int N;
3693  is >> dummy >> N;
3694  if (dummy.compare("Walls"))
3695  logger(ERROR, "DPMBase::read(is): Error during restart: 'Walls' argument could not be found.");
3696  wallHandler.clear();
3699  for (unsigned int i = 0; i < N; i++)
3700  {
3703  }
3704 
3705  is >> dummy >> N;
3708  if (dummy.compare("Boundaries"))
3709  logger(ERROR, "DPMBase::read(is): Error during restart: 'Boundaries' argument could not be found.");
3711  for (unsigned int i = 0; i < N; i++)
3712  {
3715  }
3716 
3718 
3719  is >> dummy >> N;
3720  is.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
3721  //display a message if a large amount o fparticles is read
3722  if (N>2.5e5) logger(INFO, "Reading % particles (may take a while)",N);
3723  logger.assert_always(dummy.compare("Particles")==0, "DPMBase::read(is): Error during restart: 'Particles' argument could not be found. %",dummy);
3726  for (unsigned int i = 0; i < N; i++)
3727  {
3728  //ParticleHandler::readAndCreateObject reads line-by-line
3730  //skip the remaining data in line
3731  is.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
3733  //particleHandler.getLastObject()->computeMass();
3734  }
3735 #ifdef MERCURY_USE_MPI
3736  //Interaction distances of the domainHandler and periodicBoundaryHandler need to be set
3737  Mdouble interactionRadius = particleHandler.getLargestInteractionRadius();
3738  domainHandler.setInteractionDistance(2.0*interactionRadius);
3739  periodicBoundaryHandler.setInteractionDistance(2.0*interactionRadius);
3740 
3741  if (NUMBER_OF_PROCESSORS > 1)
3742  {
3743  //Create ghost particles
3746  }
3747 #endif
3748  //Add interactions to particles and ghost particles
3749  if (opt==ReadOptions::ReadNoInteractions) return;
3751  }
3752  //reading in for older versions of the Mercury restart file.
3753  else if (!restartVersion_.compare("3"))
3754  {
3755  logger(INFO, "DPMBase::read(is): restarting from an old restart file (restart_version %).",
3756  restartVersion_);
3757  readOld(is);
3758  }
3759  //returning an error if there is no restart file to read in due to the use of outdated files.
3760  else
3761  {
3762  //only very old files did not have a restart_version
3763  logger(FATAL,
3764  "Error in DPMBase::read(is): restart_version % cannot be read; use an older version of Mercury to upgrade the file",
3765  restartVersion_);
3766  }
3767  }
3768 }
3769 
3773 void DPMBase::readOld(std::istream& is)
3774 {
3775  std::string dummy;
3776  is >> dummy >> dummy;
3777  setName(dummy);
3778 
3779  unsigned int saveCountData, saveCountEne, saveCountStat, saveCountFStat;
3780  unsigned int fileTypeFstat, fileTypeData, fileTypeEne, fileTypeRestart;
3781  is >> dummy >> min_.x()
3782  >> dummy >> max_.x()
3783  >> dummy >> min_.y()
3784  >> dummy >> max_.y()
3785  >> dummy >> min_.z()
3786  >> dummy >> max_.z()
3787  >> dummy >> timeStep_
3788  >> dummy >> time_
3789  >> dummy >> timeMax_
3790  >> dummy >> saveCountData
3791  >> dummy >> saveCountEne
3792  >> dummy >> saveCountStat
3793  >> dummy >> saveCountFStat
3794  >> dummy >> systemDimensions_
3795  >> dummy >> gravity_
3796  >> dummy >> fileTypeFstat
3797  >> dummy >> fileTypeData
3798  >> dummy >> fileTypeEne;
3799  dataFile.setSaveCount(saveCountData);
3800  eneFile.setSaveCount(saveCountEne);
3801  statFile.setSaveCount(saveCountStat);
3802  fStatFile.setSaveCount(saveCountFStat);
3803 
3804  fStatFile.setFileType(static_cast<FileType>(fileTypeFstat));
3805  dataFile.setFileType(static_cast<FileType>(fileTypeData));
3806  eneFile.setFileType(static_cast<FileType>(fileTypeEne));
3807 
3808  //this is optional to allow restart files with and without restartFile.getFileType()
3809  is >> dummy;
3810  if (!strcmp(dummy.c_str(), "options_restart"))
3811  {
3812  is >> fileTypeRestart;
3813  restartFile.setFileType(static_cast<FileType>(fileTypeRestart));
3814  }
3815 
3816  speciesHandler.read(is);
3817  wallHandler.read(is);
3818  boundaryHandler.read(is);
3819  particleHandler.read(is);
3820  setRestarted(true);
3821 }
3822 
3831 {
3832  //check if name is set
3833  logger.assert_always(getName() != "",
3834  "File name not set: use setName()");
3835  //check if time step is set
3836  logger.assert_always(getTimeStep() != 0,
3837  "Time step undefined: use setTimeStep()");
3838  //check if domain is set
3839  logger.assert_always(getXMax() > getXMin(),
3840  "Domain size not set: use setXMin() and setXMax()");
3841  logger.assert_always(getYMax() > getYMin(),
3842  "Domain size not set: use setYMin() and setYMax()");
3843  logger.assert_always(systemDimensions_ == 3 ? (getZMax() > getZMin()) : (getZMax() >= getZMin()),
3844  "Domain size not set: use setZMin() and setZMax()", systemDimensions_);
3845 
3846  //check for species parameters
3847  logger.assert_always(speciesHandler.getNumberOfObjects() > 0,
3848  "No species defined: use speciesHandler.copyAndAddObject()");
3849  for (BaseParticle* p : particleHandler)
3850  {
3851  logger.assert_always(p->getSpecies() != nullptr, "particle % has no species", p->getId());
3852  }
3853  for (BaseWall* w : wallHandler)
3854  {
3855  logger.assert_always(w->getSpecies() != nullptr, "% with index % has no species", w->getName(), w->getId());
3856  }
3857 }
3858 
3860 {
3862  writeOutputFiles();
3863 }
3864 
3881 {
3882  //writing fstat data if .saveCurrentTimeStep(numberOfTimeSteps_) is true
3884  {
3885  writeFStatFile();
3886  }
3887 
3888  //writing .ene data if .saveCurrentTimeStep(numberOfTimeSteps_) is true
3890  {
3891  writeEneFile();
3892  }
3893  //writing .data data if .saveCurrentTimeStep(numberOfTimeSteps_) is true
3895  {
3897  if (getRestarted() || dataFile.getCounter() == 0)
3899  writeDataFile();
3900  } else {
3902  }
3903  printTime();
3904  writeVTKFiles();
3905  }
3906  cgHandler.evaluate();
3907 
3908 
3909  //write restart file last, otherwise the output counters are wrong
3911  {
3912  writeRestartFile();
3913  }
3914 }
3915 
3920 {
3921 #ifdef MERCURY_USE_MPI
3922 
3923  //If running in parallel build, but just running with one core - no domain decomposition required
3924  int numberOfRequiredProcessors = numberOfDomains_[Direction::XAXIS]*
3927  if (NUMBER_OF_PROCESSORS != numberOfRequiredProcessors)
3928  {
3929  logger(ERROR,"The domain decompositions expects % processors, but only % are requested.\n"
3930  "Either run your process using \"mpirun -np % [executable]\", "
3931  "or change the domain decomposition to e.g. setNumberOfDomains({%,1,1}).", numberOfRequiredProcessors, NUMBER_OF_PROCESSORS, numberOfRequiredProcessors, NUMBER_OF_PROCESSORS);
3932  }
3933 
3934  if (NUMBER_OF_PROCESSORS == 1) {return;}
3935 
3936  //Check if the simulation domain has been set
3937  logger.assert_always(getXMax() - getXMin() > 0,"Please set your simulation domain (setXMax(),setXmin()) before calling solve()");
3938  logger.assert_always(getYMax() - getYMin() > 0,"Please set your simulation domain (setYMax(),setYmin()) before calling solve()");
3939  logger.assert_always(getZMax() - getZMin() > 0,"Please set your simulation domain (setZMax(),setZmin()) before calling solve()");
3940 
3941  //Grab simulation domains
3942  std::vector<Mdouble> simulationMin{getXMin(), getYMin(), getZMin()};
3943  std::vector<Mdouble> simulationMax{getXMax(), getYMax(), getZMax()};
3944 
3945  //Check if the user input decomposition is correct
3946  logger.assert_always(numberOfDomains_[Direction::XAXIS] > 0,"Number of domain in x-direction incorrect: %",numberOfDomains_[Direction::XAXIS]);
3947  logger.assert_always(numberOfDomains_[Direction::YAXIS] > 0,"Number of domain in y-direction incorrect: %",numberOfDomains_[Direction::YAXIS]);
3948  logger.assert_always(numberOfDomains_[Direction::ZAXIS] > 0,"Number of domain in z-direction incorrect: %",numberOfDomains_[Direction::ZAXIS]);
3949 
3950  //Open domain decomposition, closed is not implemented
3951  bool open = true;
3952 
3953  //Check if the number of domains is equal to the number of processors
3954  logger.assert_always(numberOfDomains_[Direction::XAXIS]*numberOfDomains_[Direction::YAXIS]*numberOfDomains_[Direction::ZAXIS] == NUMBER_OF_PROCESSORS,
3955  "Number of Processors is not equal to number of domains. Processors %, domains, %",
3957  numberOfDomains_[Direction::XAXIS]*numberOfDomains_[Direction::YAXIS]*numberOfDomains_[Direction::ZAXIS]);
3958 
3959  //Create all processor domains
3960 
3961  domainHandler.setDPMBase(this);
3962  domainHandler.createMesh(simulationMin, simulationMax, numberOfDomains_, open);
3963  logger(VERBOSE,"Number of domains: % | Number of processors: %",domainHandler.getNumberOfObjects(), NUMBER_OF_PROCESSORS);
3964  //logger.assert_always(domainHandler.getNumberOfObjects() == numberOfProcessors, "Invalid decomposition: Number of domains and processors are different");
3965 
3966  //Tell the current processor to which domain it belongs
3967  for (Domain* domain : domainHandler)
3968  {
3969  if (domain->getRank() == PROCESSOR_ID)
3970  {
3971  logger(VERBOSE,"processor: %, domain index: %",PROCESSOR_ID, domain->getIndex());
3972  domainHandler.setCurrentDomainIndex(domain->getIndex());
3973  }
3974  }
3975 
3976  //Define the mpi transfer types, which requires a definition of the species already
3978  logger.assert_always(speciesHandler.getNumberOfObjects() > 0, "Please create a particle species before calling solve()");
3980 
3981  //Make sure all processors are done with decomposition before proceeding
3982  logger(VERBOSE,"processor %: #real particles: %, #total particles: %", PROCESSOR_ID, particleHandler.getNumberOfRealObjects(), particleHandler.getSize());
3984 #endif
3985 }
3986 
3987 
4003 {
4004  logger(DEBUG, "Entered solve");
4005 #ifdef CONTACT_LIST_HGRID
4006  logger(INFO,"Using CONTACT_LIST_HGRID");
4007 #endif
4008 
4013  if (!getRestarted())
4014  {
4015  // If the simulation is "new" (i.e. not restarted):
4016  // - set time, nTimeSteps to zero
4017  // - reset the file counter etc.
4018  // - decompose the domain based on XMin, XMax, ....
4019  // - run user-defined setupInitialConditions
4020  numberOfTimeSteps_ = 0;
4021  setTime(0.0);
4022  resetFileCounter();
4023  decompose();
4024  //\todo tw there was a function combining the next two lines, why is it back to the old version?
4025  //setLastSavedTimeStep(NEVER); //reset the counter
4026  //this is to ensure that the interaction time stamps agree with the resetting of the time value
4027  for (auto& i : interactionHandler)
4028  i->setTimeStamp(0);
4030  logger(DEBUG, "Have created the particles initial conditions");
4031  }
4032  else
4033  {
4034  // If the simulation is "restarted" (i.e. not restarted):
4035 
4036  // - run wall-defined actionsOnRestart
4037  for (auto w: wallHandler)
4038  {
4039  w->actionsOnRestart();
4040  }
4041 
4042  // - run user-defined actionsOnRestart
4043  actionsOnRestart();
4044  }
4045 
4046  // Check that the code has been correctly set up,
4047  // i.e. system dimensions, particles and time steps are sensibly implemented
4048  checkSettings();
4049 
4050  // If the simulation is "new" and the runNumber is used, append the run number to the problem name
4051  if (getRunNumber() > 0 && !getRestarted())
4052  {
4053  std::stringstream name;
4054  name << getName() << "." << getRunNumber();
4055  setName(name.str());
4056  }
4057 
4058  //If append is true, files are appended, not overwritten
4059  if (getAppend())
4060  {
4061  setOpenMode(std::fstream::out | std::fstream::app);
4062  //Restart files should always be overwritten.
4063  restartFile.setOpenMode(std::fstream::out);
4064  }
4065  else
4066  {
4067  setOpenMode(std::fstream::out);
4068  }
4069 
4070  //sets the hgrid, writes headers to the .stat output file
4072 
4073  if (getInteractionFile().getFileType() == FileType::ONE_FILE)
4074  {
4075  logger(WARN, "Warning: interaction file will take up a lot of disk space!");
4077  }
4078 
4079  // Sets the mass of all particle.
4083 
4084  // Other initialisations
4085  //max_radius = getLargestParticle()->getRadius();
4086 
4090 
4091  // Performs a first force computation
4093 
4094 #ifdef MERCURY_USE_MPI
4095  if (NUMBER_OF_PROCESSORS > 1)
4096  {
4097  //Find new mpi particles
4099  //Periodic particles in parallel
4101  }
4102 #endif
4103 
4105  computeAllForces();
4108  logger(DEBUG, "Have computed the initial values for the forces ");
4109 
4110  // Can be used to measure simulation time
4111  clock_.tic();
4112 
4113  // This is the main loop over advancing time
4114  while (getTime() < getTimeMax() && continueSolve())
4115  {
4117  }
4118 
4119  // Can be used to measure simulation time
4120  clock_.toc();
4121 
4122  //force writing of the last time step
4124 
4125  //end loop over interaction count
4127 
4128  //To make sure getTime gets the correct time for outputting statistics
4129  finishStatistics();
4130 
4131  closeFiles();
4132 }
4133 
4139 {
4140  logger(DEBUG, "starting computeOneTimeStep()");
4141 
4142  logger(DEBUG, "about to call writeOutputFiles()");
4143  writeOutputFiles(); //everything is written at the beginning of the time step!
4144 
4145  logger(DEBUG, "about to call hGridActionsBeforeIntegration()");
4147 
4148  //Computes the half-time step velocity and full time step position and updates the particles accordingly
4149  logger(DEBUG, "about to call integrateBeforeForceComputation()");
4150 
4152  //New positions require the MPI and parallel periodic boundaries to do things
4153  logger(DEBUG, "about to call performGhostParticleUpdate()");
4155  // Some walls need to be aware of the new positions
4157 
4161 
4162  logger(DEBUG, "about to call checkInteractionWithBoundaries()");
4163  checkInteractionWithBoundaries(); // INSERTION boundaries handled
4164 
4165  logger(DEBUG, "about to call hGridActionsAfterIntegration()");
4167 
4168  // Compute forces
4170  // INSERTION/DELETION boundary flag change
4171  for (BaseBoundary* b : boundaryHandler)
4172  {
4173  b->checkBoundaryBeforeTimeStep(this);
4174  }
4175 
4176  logger(DEBUG, "about to call actionsBeforeTimeStep()");
4178 
4179  logger(DEBUG, "about to call checkAndDuplicatePeriodicParticles()");
4181 
4182  logger(DEBUG, "about to call hGridActionsBeforeTimeStep()");
4184 
4185  //Creates and updates interactions and computes forces based on these
4186  logger(DEBUG, "about to call computeAllForces()");
4187  computeAllForces();
4188 
4189  logger(DEBUG, "about to call removeDuplicatePeriodicParticles()");
4191 
4192  logger(DEBUG, "about to call actionsAfterTimeStep()");
4194 
4195  //Computes new velocities and updates the particles accordingly
4196  logger(DEBUG, "about to call integrateAfterForceComputation()");
4198 
4199  //erase interactions that have not been used during the last time step
4200  //logger(DEBUG, "about to call interactionHandler.eraseOldInteractions(getNumberOfTimeSteps())");
4202  logger(DEBUG, "about to call interactionHandler.actionsAfterTimeStep()");
4205 
4206  time_ += timeStep_;
4208 
4209  logger(DEBUG, "finished computeOneTimeStep()");
4210 }
4211 
4224 bool DPMBase::readArguments(int argc, char* argv[])
4225 {
4226  bool isRead = true;
4227  // Cycles over every second element. (Most flags will contain both name labels and actual data.
4228  // Those that don't will have to do i--; some examples in readNextArgument.)
4229  for (int i = 1; i < argc; i += 2)
4230  {
4231  logger(INFO, "interpreting input argument %\n", argv[i], Flusher::NO_FLUSH);
4232  for (int j = i + 1; j < argc; j++)
4233  {
4234  //looks for the next string that starts with a minus sign
4235  //i.e. the next flag, as each flag may take 0, 1 , 2, 3... arguments
4236  //and we need to make sure all are read in!
4237  if (argv[j][0] == '-')
4238  break;
4239  logger(INFO, " %\n", argv[j], Flusher::NO_FLUSH);
4240  }
4241  logger(INFO, "");
4242  //if "isRead"is true and "readNextArgument" is also true...
4243  //(i.e. checking if any argument is false)
4244  isRead &= readNextArgument(i, argc, argv);
4245 
4246  // If the read was unsuccessful, raise an error and quit. (JMFT: Used to just be a warning.)
4247  if (!isRead)
4248  {
4249  logger(ERROR, "Warning: not all arguments read correctly!");
4250  }
4251  }
4252  return isRead;
4253 }
4254 
4256 {
4257  //logger(INFO,"ID %",PROCESSOR_ID);
4258  //if (PROCESSOR_ID!=0) return;
4259 
4260  logger(INFO,"Removing old files named %.*",getName());
4261  std::ostringstream filename;
4262 
4263  // add processor id to file extension for mpi jobs
4264  std::string p = (NUMBER_OF_PROCESSORS > 1)?std::to_string(PROCESSOR_ID):"";
4265  // all the file extensions that should be deleted
4266  std::vector<std::string> ext{".restart"+p, ".stat"+p, ".fstat"+p, ".data"+p, ".ene"+p, ".xballs"};
4267  for (const auto& j : ext)
4268  {
4269  // remove files with given extension for FileType::ONE_FILE
4270  filename.str("");
4271  filename << getName() << j;
4272  if (!remove(filename.str().c_str()))
4273  {
4274  logger(INFO," File % successfully deleted",filename.str());
4275  }
4276  // remove files with given extension for FileType::MULTIPLE_FILES
4277  unsigned k = 0;
4278  filename.str("");
4279  filename << getName() << j << '.' << k;
4280  while (!remove(filename.str().c_str()))
4281  {
4282  if (k<3) logger(INFO," File % successfully deleted",filename.str());
4283  filename.clear();
4284  filename << getName() << j << '.' << ++k;
4285  }
4286  // remove files with given extension for FileType::MULTIPLE_FILES_PADDED
4287  k = 0;
4288  filename.str("");
4289  filename << getName() << j << '.' << to_string_padded(k);
4290  while (!remove(filename.str().c_str()))
4291  {
4292  if (k<3) logger(INFO," File % successfully deleted",filename.str());
4293  filename.clear();
4294  filename << getName() << j << '.' << to_string_padded(++k);
4295  }
4296  }
4297  // remove vtk files
4298  // add processor id to file extension for mpi jobs
4299  std::string q = (NUMBER_OF_PROCESSORS > 1)?("Processor_"+std::to_string(PROCESSOR_ID)+"_"):"";
4300  // all the file extensions that should be deleted
4301  ext = {"Wall_", q+"Particle_", q+"Interaction_"};
4302  for (const auto& j : ext)
4303  {
4304  // remove files with given extension for FileType::ONE_FILE
4305  filename.str("");
4306  filename << getName() << j << ".vtu";
4307  if (!remove(filename.str().c_str()))
4308  {
4309  logger(INFO," File % successfully deleted",filename.str());
4310  }
4311  // remove files with given extension for FileType::MULTIPLE_FILES
4312  unsigned k = 0;
4313  filename.str("");
4314  filename << getName() << j << k << ".vtu";
4315  while (!remove(filename.str().c_str()))
4316  {
4317  if (k<3) logger(INFO," File % successfully deleted",filename.str());
4318  filename.str("");
4319  filename << getName() << j << ++k << ".vtu";
4320  }
4321  //std::cout << " File " << filename.str() << " not found" << std::endl;
4322  }
4323 }
4324 
4349 bool DPMBase::readNextArgument(int& i, int argc, char* argv[])
4350 {
4351  // The argument argv[i] identifies the label of the flag, and subsequent arguments (usually 1)
4352  // contain the content.
4353  //
4354  // For example...
4355  // Checks if the "-name" flag has been passed
4356  // The strcmp returns 0 if "argv[i]" is "-name" (i.e. !strcmp(argv[i], "-name") --> 1)
4357  // In this case, the setName function is run with the relevant input (i.e. the value which
4358  // immediately follows the "-name" flag
4359  if (!strcmp(argv[i], "-name"))
4360  {
4361  setName(argv[i + 1]);
4362  }
4363  // The above process is repeated for all viable flags.
4364  else if (!strcmp(argv[i], "-xmin"))
4365  {
4366  setXMin(atof(argv[i + 1]));
4367  }
4368  else if (!strcmp(argv[i], "-ymin"))
4369  {
4370  setYMin(atof(argv[i + 1]));
4371  }
4372  else if (!strcmp(argv[i], "-zmin"))
4373  {
4374  setZMin(atof(argv[i + 1]));
4375  }
4376  else if (!strcmp(argv[i], "-xmax"))
4377  {
4378  setXMax(atof(argv[i + 1]));
4379  }
4380  else if (!strcmp(argv[i], "-ymax"))
4381  {
4382  setYMax(atof(argv[i + 1]));
4383  }
4384  else if (!strcmp(argv[i], "-zmax"))
4385  {
4386  setZMax(atof(argv[i + 1]));
4387  //} else if (!strcmp(argv[i],"-svn")) {
4388  // std::cout << "svn version " << SVN_VERSION << std::endl;
4389  // i--;
4390  }
4391  else if (!strcmp(argv[i], "-dt"))
4392  {
4393  Mdouble old = getTimeStep();
4394  setTimeStep(atof(argv[i + 1]));
4395  logger(INFO, " reset dt from % to %", old, getTimeStep());
4396  }
4397 // else if (!strcmp(argv[i], "-Hertz"))
4398 // {
4399 // speciesHandler.getObject(0)->setForceType(ForceType::HERTZ);
4400 // i--;
4401 // }
4402  else if (!strcmp(argv[i], "-tmax"))
4403  {
4404  Mdouble old = getTimeMax();
4405  setTimeMax(atof(argv[i + 1]));
4406  logger(INFO, " reset timeMax from % to %", old, getTimeMax());
4407  }
4408  else if (!strcmp(argv[i], "-saveCount"))
4409  {
4410  Mdouble old = dataFile.getSaveCount();
4411  setSaveCount(static_cast<unsigned int>(atoi(argv[i + 1])));
4412  logger(INFO, " reset saveCount from & to %", old, dataFile.getSaveCount());
4413  }
4414  else if (!strcmp(argv[i], "-saveCountData"))
4415  {
4416  dataFile.setSaveCount(static_cast<unsigned int>(atoi(argv[i + 1])));
4417  }
4418  else if (!strcmp(argv[i], "-saveCountFStat"))
4419  {
4420  fStatFile.setSaveCount(static_cast<unsigned int>(atoi(argv[i + 1])));
4421  }
4422  else if (!strcmp(argv[i], "-saveCountStat"))
4423  {
4424  statFile.setSaveCount(static_cast<unsigned int>(atoi(argv[i + 1])));
4425  }
4426  else if (!strcmp(argv[i], "-saveCountEne"))
4427  {
4428  eneFile.setSaveCount(static_cast<unsigned int>(atoi(argv[i + 1])));
4429  }
4430  else if (!strcmp(argv[i], "-saveCountRestart"))
4431  {
4432  restartFile.setSaveCount(static_cast<unsigned int>(atoi(argv[i + 1])));
4433  }
4434  else if (!strcmp(argv[i], "-dim"))
4435  {
4436  setSystemDimensions(static_cast<unsigned int>(atoi(argv[i + 1])));
4437  }
4438  else if (!strcmp(argv[i], "-gravity"))
4439  {
4441  setGravity(Vec3D(atof(argv[i + 1]), atof(argv[i + 2]), atof(argv[i + 3])));
4442  i += 2;
4443  }
4444  else if (!strcmp(argv[i], "-fileType"))
4445  { //uses int input
4446  setFileType(static_cast<FileType>(atoi(argv[i + 1])));
4447  }
4448  else if (!strcmp(argv[i], "-fileTypeFStat"))
4449  { //uses int input
4450  fStatFile.setFileType(static_cast<FileType>(atoi(argv[i + 1])));
4451  }
4452  else if (!strcmp(argv[i], "-fileTypeRestart"))
4453  {
4454  restartFile.setFileType(static_cast<FileType>(atoi(argv[i + 1])));
4455  }
4456  else if (!strcmp(argv[i], "-fileTypeData"))
4457  {
4458  dataFile.setFileType(static_cast<FileType>(atoi(argv[i + 1])));
4459  }
4460  else if (!strcmp(argv[i], "-fileTypeStat"))
4461  {
4462  statFile.setFileType(static_cast<FileType>(atoi(argv[i + 1])));
4463  }
4464  else if (!strcmp(argv[i], "-fileTypeEne"))
4465  {
4466  eneFile.setFileType(static_cast<FileType>(atoi(argv[i + 1])));
4467  }
4468  else if (!strcmp(argv[i], "-auto_number"))
4469  {
4470  autoNumber();
4471  i--;
4472  }
4473 // else if (!strcmp(argv[i], "-number_of_saves"))
4474 // {
4475 // set_number_of_saves_all(atof(argv[i + 1]));
4476 // }
4477  else if (!strcmp(argv[i], "-restart") || !strcmp(argv[i], "-r"))
4478  {
4482  std::string filename;
4483 
4484  //use default filename if no argument is given
4485  if (i + 1 >= argc || argv[i + 1][0] == '-')
4486  {
4487  i--;
4488  filename = getName();
4489  logger(INFO, "%", getName());
4490  }
4491  else
4492  {
4493  filename = argv[i + 1];
4494  }
4495 
4496  //add ".restart" if necessary
4497  if (filename.find(".restart") == std::string::npos)
4498  {
4499  filename = filename + ".restart";
4500  }
4501 
4502  readRestartFile(filename);
4503  }
4504  else if (!strcmp(argv[i], "-clean") || !strcmp(argv[i], "-c"))
4505  {
4506  logger(INFO, "Remove old %.* files", getName());
4507  removeOldFiles();
4508  i--;
4509  }
4510  else if (!strcmp(argv[i], "-data"))
4511  {
4512  std::string filename = argv[i + 1];
4513  readDataFile(filename);
4514  }
4515  else if (!strcmp(argv[i], "-readSpeciesFromDataFile"))
4516  {
4517  readSpeciesFromDataFile_ = true;
4518  i--;
4519  logger(INFO, "Last column of data file will be interpreted as species index");
4520  }
4521 // else if (!strcmp(argv[i], "-k"))
4522 // {
4523 // Mdouble old = getSpecies(0)->getStiffness();
4524 // getSpecies(0)->setStiffness(atof(argv[i + 1]));
4525 // std::cout << " reset k from " << old << " to " << getSpecies(0)->getStiffness() << std::endl;
4526 // }
4527 // else if (!strcmp(argv[i], "-dissipation") || !strcmp(argv[i], "-disp"))
4528 // {
4529 // Mdouble old = getSpecies(0)->getDissipation();
4530 // getSpecies(0)->setDissipation(atof(argv[i + 1]));
4531 // std::cout << " reset getDissipation() from " << old << " to " << getSpecies(0)->getDissipation() << std::endl;
4532 // }
4533 // else if (!strcmp(argv[i], "-kt"))
4534 // {
4535 // Mdouble old = getSpecies(0)->getSlidingStiffness();
4536 // getSpecies(0)->setSlidingStiffness(atof(argv[i + 1]));
4537 // std::cout << " reset kt from " << old << " to " << getSpecies(0)->getSlidingStiffness() << std::endl;
4538 // }
4539 // else if (!strcmp(argv[i], "-dispt"))
4540 // {
4541 // Mdouble old = getSpecies(0)->getSlidingDissipation();
4542 // getSpecies(0)->setSlidingDissipation(atof(argv[i + 1]));
4543 // std::cout << " reset dispt from " << old << " to " << getSpecies(0)->getSlidingDissipation() << std::endl;
4544 // }
4545 // else if (!strcmp(argv[i], "-krolling"))
4546 // {
4547 // Mdouble old = getSpecies(0)->getRollingStiffness();
4548 // getSpecies(0)->setRollingStiffness(atof(argv[i + 1]));
4549 // std::cout << " reset krolling from " << old << " to " << getSpecies(0)->getRollingStiffness() << std::endl;
4550 // }
4551 // else if (!strcmp(argv[i], "-disprolling"))
4552 // {
4553 // Mdouble old = getSpecies(0)->getRollingDissipation();
4554 // getSpecies(0)->setRollingDissipation(atof(argv[i + 1]));
4555 // std::cout << " reset disprolling from " << old << " to " << getSpecies(0)->getRollingDissipation() << std::endl;
4556 // }
4557 // else if (!strcmp(argv[i], "-mu"))
4558 // {
4559 // Mdouble old = getSpecies(0)->getSlidingFrictionCoefficient();
4560 // getSpecies(0)->setSlidingFrictionCoefficient(atof(argv[i + 1]));
4561 // std::cout << " reset mu from " << old << " to " << getSpecies(0)->getSlidingFrictionCoefficient() << std::endl;
4562 // }
4563 // else if (!strcmp(argv[i], "-murolling"))
4564 // {
4565 // Mdouble old = getSpecies(0)->getRollingFrictionCoefficient();
4566 // getSpecies(0)->setRollingFrictionCoefficient(atof(argv[i + 1]));
4567 // std::cout << " reset murolling from " << old << " to " << getSpecies(0)->getRollingFrictionCoefficient() << std::endl;
4568 // }
4569  else if (!strcmp(argv[i], "-randomise") || !strcmp(argv[i], "-randomize"))
4570  {
4571  random.randomise();
4572  i--;
4573  }
4574 // else if (!strcmp(argv[i], "-k0"))
4575 // {
4576 // Mdouble old = speciesHandler.getObject(0)->getAdhesionStiffness();
4577 // speciesHandler.getObject(0)->setAdhesionStiffness(atof(argv[i + 1]));
4578 // std::cout << " reset k0 from " << old << " to " << speciesHandler.getObject(0)->getAdhesionStiffness() << std::endl;
4579 // }
4580 // else if (!strcmp(argv[i], "-f0"))
4581 // {
4582 // Mdouble old = speciesHandler.getObject(0)->getBondForceMax();
4583 // speciesHandler.getObject(0)->setBondForceMax(atof(argv[i + 1]));
4584 // std::cout << " reset f0 from " << old << " to " << speciesHandler.getObject(0)->getBondForceMax() << std::endl;
4585 // }
4586 // else if (!strcmp(argv[i], "-AdhesionForceType"))
4587 // {
4588 // AdhesionForceType old = speciesHandler.getObject(0)->getAdhesionForceType();
4589 // speciesHandler.getObject(0)->setAdhesionForceType(argv[i + 1]);
4590 // std::cout << " reset AdhesionForceType from "
4591 // << static_cast<signed char>(old) << " to "
4592 // << static_cast<signed char>(speciesHandler.getObject(0)->getAdhesionForceType()) << std::endl;
4593 // }
4594  else if (!strcmp(argv[i], "-append"))
4595  {
4596  setAppend(true);
4597  i--;
4598  }
4599  else if (!strcmp(argv[i], "-fixedParticles"))
4600  {
4601  setFixedParticles(static_cast<unsigned int>(atoi(argv[i + 1])));
4602  }
4603 // else if (!strcmp(argv[i], "-rho"))
4604 // {
4605 // Mdouble old = speciesHandler.getObject(0)->getDensity();
4606 // speciesHandler.getObject(0)->setDensity(atof(argv[i + 1]));
4607 // std::cout << " reset rho from " << old << " to " << speciesHandler.getObject(0)->getDensity() << std::endl;
4608 // }
4609 // else if (!strcmp(argv[i], "-dim_particle"))
4610 // {
4611 // setParticleDimensions(atoi(argv[i + 1]));
4612 // }
4613  else if (!strcmp(argv[i], "-counter"))
4614  {
4615  setRunNumber(atoi(argv[i + 1]));
4616  }
4617  else
4618  {
4619  //returns false if the flag passed does not match any of the currently used flags.
4620  return false;
4621  }
4622  //returns true if argv is found
4623  return true;
4624 }
4625 
4639 {
4640 #ifdef MERCURY_USE_MPI
4641  if (NUMBER_OF_PROCESSORS == 1)
4642  {
4644  }
4645 
4646  int localInteraction = checkParticleForInteractionLocal(p);
4647  //The root gathers all values and computes the global value
4648  int *interactionList = nullptr;
4649  if (PROCESSOR_ID == 0)
4650  {
4651  interactionList = new int [NUMBER_OF_PROCESSORS];
4652  }
4653 
4654  //Gather all local values
4655  MPIContainer::Instance().gather(localInteraction,interactionList);
4656 
4657  //Compute the global value
4658  int globalInteraction = 1;
4659  if (PROCESSOR_ID == 0)
4660  {
4661  for (int i = 0; i < NUMBER_OF_PROCESSORS; i++)
4662  {
4663  if (interactionList[i] == 0)
4664  {
4665  globalInteraction = 0;
4666  break;
4667  }
4668  }
4669  }
4670  //The root now tells the other processors what the global value for the interaction is
4671  MPIContainer::Instance().broadcast(globalInteraction);
4672 
4673  //Convert the result back to bool
4674  bool interaction = globalInteraction;
4675 #else
4676  bool interaction = checkParticleForInteractionLocalPeriodic(p);
4677 #endif
4678  return interaction;
4679 }
4680 
4688 {
4689  //A vector of ghost particles of the particle that is to be inserted (empty if no periodic boundaries are present)
4690  std::vector<Vec3D> pPeriodic;
4691  for (BaseBoundary* b : boundaryHandler)
4692  {
4693  PeriodicBoundary* pb = dynamic_cast<PeriodicBoundary*>(b);
4694  if (pb && particleHandler.getNumberOfObjects() > 0 )
4695  {
4697  for (int i = pPeriodic.size() - 1; i >= 0; --i)
4698  {
4699  if (pb->getDistance(pPeriodic[i]) < maxDistance)
4700  {
4701  pPeriodic.push_back(pPeriodic[i]);
4702  pb->shiftPosition(pPeriodic.back());
4703  }
4704  }
4705  if (pb->getDistance(p) < maxDistance)
4706  {
4707  pPeriodic.push_back(p.getPosition());
4708  pb->shiftPosition(pPeriodic.back());
4709  }
4710  }
4711  }
4712  //check the particle AND the ghost particles for intersection problems
4713  bool insertable = checkParticleForInteractionLocal(p);
4714  if (!pPeriodic.empty()) {
4715  BaseParticle* q = p.copy();
4716  for (const Vec3D& pos : pPeriodic) {
4717  q->setPosition(pos);
4718  insertable &= checkParticleForInteractionLocal(*q);
4719  }
4720  delete q;
4721  }
4722  return insertable;
4723 }
4724 
4738 {
4739  Mdouble distance;
4740  Vec3D normal;
4741 
4742  //Check if it has no collision with walls
4743  for (BaseWall* w : wallHandler)
4744  {
4745  //returns false if the function getDistanceAndNormal returns true,
4746  //i.e. if there exists an interaction between wall and particle
4747  //\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.)
4748  if (w->getDistanceAndNormal(p, distance, normal))
4749  {
4750  //std::cout<<"failure: Collision with wall: "<<**it<<std::endl;
4751  return false;
4752  }
4753  else
4754  {
4755  //std::cout<<"No collision with wall: "<<**it<<std::endl;
4756  }
4757  }
4758 
4759  //Check if it has no collision with other particles
4760  for (BaseParticle* q : particleHandler)
4761  {
4762  //returns false if the particle separation is less than the relevant sum of interaction radii
4763  //(i.e. another particle is in contact with P)
4764  if (Vec3D::getDistanceSquared(q->getPosition(), p.getPosition())
4766  {
4767  //std::cout<<"failure: Collision with particle "<<**it<<std::endl;
4768  return false;
4769  }
4770  else
4771  {
4772  //std::cout<<"No collision with particle "<<**it<<std::endl;
4773  }
4774  }
4775  return true;
4777 }
4778 
4787 void DPMBase::importParticlesAs(ParticleHandler& particleH, InteractionHandler& interactionH, const ParticleSpecies* species )
4788 {
4789  size_t nParticlesPreviouslyIn = particleHandler.getSize();
4790  int l = 0;
4791  for (auto k = particleH.begin(); k != particleH.end(); ++k) {
4792  auto p = particleHandler.copyAndAddObject( *k );
4793  p->setSpecies(species);
4794  l++;
4795  }
4796 
4797  for (std::vector<BaseInteraction*>::const_iterator i = interactionH.begin(); i != interactionH.end(); ++i) {
4798  if ( (*i)->getP()->getInvMass() != 0.0 && (*i)->getI()->getInvMass() != 0.0 ) {
4800  j->importP(particleHandler.getObject(nParticlesPreviouslyIn + j->getP()->getIndex()));
4801  j->importI(particleHandler.getObject(nParticlesPreviouslyIn + j->getI()->getIndex()));
4802  j->setTimeStamp(getNumberOfTimeSteps());
4803  }
4804  }
4805 }
4806 
4822 {
4823 #ifdef MERCURY_USE_MPI
4824  if (NUMBER_OF_PROCESSORS == 1)
4826  {
4827 #endif
4828  //Looping from the *back* of the particle handler, where all the periodic or "ghost" particles are stored.
4829  //The loop ends when there exist no more periodic particles (i.e. "getPeriodicFromParticle()" returns a null pointer)
4830  for (unsigned int i = particleHandler.getSize();
4831  i >= 1 && particleHandler.getObject(i - 1)->getPeriodicFromParticle() != nullptr; i--)
4832  {
4833  while (!particleHandler.getObject(i - 1)->getInteractions().empty())
4834  {
4836  particleHandler.getObject(i - 1)->getInteractions().front()->getIndex());
4837  }
4839  }
4840 
4841  // OMP parallelism
4842  /*#pragma omp parallel for num_threads(getNumberOfOMPThreads()) //schedule(dynamic)
4843  for (unsigned int i = particleHandler.getSize(); i >= 1 ; i--)
4844  {
4845  if (particleHandler.getObject(i - 1)->getPeriodicFromParticle() != nullptr) {
4846  while (!particleHandler.getObject(i - 1)->getInteractions().empty()) {
4847  interactionHandler.removeObjectKeepingPeriodics(
4848  particleHandler.getObject(i - 1)->getInteractions().front()->getIndex());
4849  }
4850  particleHandler.removeObject(i - 1);
4851  }
4852  }*/
4853 
4854 #ifdef MERCURY_USE_MPI
4855  }
4856 #endif
4857 }
4858 
4865 {
4866  //Looping over all boundaries in the boundaryHandler
4867  for (BaseBoundary* boundary : boundaryHandler)
4868  {
4869  //Calls the createPeriodicParticles() function which checks if a particle is adequately
4870  //close to a periodic particle that a periodic (ghost) particle should be created and,
4871  //if so, creates one and adds it to the system (hence the necessity to keep "N" variable).
4872  //
4873  // (The loop is over all boundaries, but if a boundary is not a PeriodicBoundary, then
4874  // this does nothing.)
4875  boundary->createPeriodicParticles(particleHandler);
4876  }
4877 
4878  // OMP parallelism
4879  /*#pragma omp parallel for num_threads(getNumberOfOMPThreads()) //schedule(dynamic)
4880  for (int k = 0; k < boundaryHandler.getNumberOfObjects(); k++)
4881  {
4882  //Calls the createPeriodicParticles() function which checks if a particle is adequately
4883  //close to a periodic particle that a periodic (ghost) particle should be created and,
4884  //if so, creates one and adds it to the system (hence the necessity to keep "N" variable).
4885  //
4886  // (The loop is over all boundaries, but if a boundary is not a PeriodicBoundary, then
4887  // this does nothing.)
4888 
4889  BaseBoundary* boundary = boundaryHandler.getObject(k);
4890  #pragma omp critical
4891  boundary->createPeriodicParticles(particleHandler);
4892  }*/
4893 }
4894 
4899 {
4900 #ifdef MERCURY_USE_MPI
4901  //MPIContainer& communicator = MPIContainer::Instance();
4902  if (NUMBER_OF_PROCESSORS == 1) {return;}
4903 
4904  //Update the postion and velocity data of ghosts and perform some bookkeeping
4905  std::set<BaseParticle*> particlesToBeDeleted;
4906  domainHandler.updateStatus(particlesToBeDeleted);
4907  periodicBoundaryHandler.updateStatus(particlesToBeDeleted);
4908 
4909  //Delete particles
4910  deleteGhostParticles(particlesToBeDeleted);
4911 
4912  //Add new particles
4915 #endif
4916 }
4917 
4921 void DPMBase::deleteGhostParticles(std::set<BaseParticle*>& particlesToBeDeleted)
4922 {
4923  //Flush mixed particles from lists (MP particles are located in both structures)
4924  if (periodicBoundaryHandler.getSize() > 0)
4925  {
4926  //Flush particles from boundaries
4927  domainHandler.getCurrentDomain()->flushParticles(particlesToBeDeleted);
4928  periodicBoundaryHandler.flushParticles(particlesToBeDeleted);
4929  }
4930 
4931  //Clean communication lists
4934 
4935  //Delete the particles
4936  for (auto particle_it : particlesToBeDeleted)
4937  {
4938  particleHandler.removeGhostObject(particle_it->getIndex());
4939  }
4940 }
4941 
4942 
4943 //This function takes a base particle and then copies all the particle information from the root particle to the other particles
4944 // on neighbouring domains.
4945 void DPMBase::synchroniseParticle(BaseParticle* p, unsigned fromProcessor)
4946 {
4947 #ifdef MERCURY_USE_MPI
4948  MPIContainer& communicator = MPIContainer::Instance();
4949 
4950  //The processor that contains the particle that needs to be copied needs to identify the target, and communicate this
4951  MPIParticle pInfo;
4952  if (communicator.getProcessorID() == fromProcessor)
4953  {
4955  }
4956 
4957  //Broadcast from processor i
4958  communicator.broadcast(&pInfo,MercuryMPIType::PARTICLE,fromProcessor);
4960 #endif
4961 }
4962 
4964 {
4965 #ifdef MERCURY_USE_MPI
4966  if (NUMBER_OF_PROCESSORS == 1) {return;}
4967  //TODO If required, I can implement this for periodic particles, first discuss with Thomas if it is actually requiredf
4968  //periodicDomainHandler.updateVelocity()
4969  //domainHandler.updateVelocity();
4970 #endif
4971 }
4972 
4976 void DPMBase::signalHandler(int signal) {
4977  switch (signal) {
4978  case SIGINT:
4979  logger(INFO, "SIGINT has been captured!\nMercuryDPM is writing the files, then it will stop!");
4980  // continue Flag must be set to false here!
4981  continueFlag_ = false;
4982  //exit(SIGTERM); // this will interrupt the simulation
4983  return;
4984  case SIGTERM:
4985  logger(INFO, "\nSIGTERM has been captured!\nMercuryDPM is writing the files, then it will stop!");
4986  // continue Flag must be set to false here!
4987  continueFlag_ = false;
4988  return;
4989  case SIGKILL:
4990  logger(INFO, "\nSIGKILL has been captured!\nMercuryDPM is writing the files, then it will stop!");
4991  continueFlag_ = false;
4992  return;
4993  default:
4994  logger(INFO, "No Signal to Capture!");
4995  }
4996 }
4997 
4998 
5003  logger(INFO,"Initiated soft stop");
5004  struct sigaction act{};
5005  memset(&act, 0, sizeof(act));
5006 
5007  act.sa_handler = &signalHandler;
5008 
5009  sigaction(SIGINT, &act, nullptr);
5010  sigaction(SIGTERM, &act, nullptr);
5011  sigaction(SIGKILL, &act, nullptr);
5012 }
5013 
5019 {
5020  logger(INFO, "Interactions currently in the handler:\n", Flusher::NO_FLUSH);
5021  //looping over all individual objects in the interactionHandler
5023  {
5024  p->write(std::cout);
5025  logger(INFO, "\nInteraction % % between % and %",
5026  p->getName(), p->getId(), p->getP()->getId(), p->getI()->getId());
5027  }
5028 }
5029 
5039 {
5040  return getTime() <= time && getTime() + getTimeStep() > time;
5041 }
5042 
5048 void DPMBase::setNumberOfDomains(std::vector<unsigned> numberOfDomains)
5049 {
5050 #ifdef MERCURY_USE_MPI
5051  numberOfDomains_ = numberOfDomains;
5052  logger(INFO, "Split domain into a %x%x% grid",numberOfDomains[0],numberOfDomains[1],numberOfDomains[2]);
5053 #else
5054  logger(WARN, "Setting number of domains, but code is not compiled with MPI on");
5055 #endif
5056 }
5057 
5059  //one-d problems
5060  if (domainSplit == DomainSplit::X) {
5062  return;
5063  } else if (domainSplit == DomainSplit::Y) {
5065  return;
5066  } else if (domainSplit == DomainSplit::Z) {
5068  return;
5069  }
5070  //two-d problems
5071  // split into axb grid with a the largest integer that divides NUMBER_OF_PROCESSORS and is smaller than sqrt(NUMBER_OF_PROCESSORS)
5072  unsigned a;
5073  for (unsigned n = floor(sqrt(NUMBER_OF_PROCESSORS));n>0; n--) {
5074  if (NUMBER_OF_PROCESSORS % n == 0) {
5075  a = n;
5076  break;
5077  }
5078  }
5079  if (domainSplit == DomainSplit::XY) {
5081  return;
5082  } else if (domainSplit == DomainSplit::XZ) {
5084  return;
5085  } else if (domainSplit == DomainSplit::YZ) {
5087  return;
5088  }
5089  //three-d problems
5090  // split into axbxc grid with
5091  // - a the largest integer that divides NUMBER_OF_PROCESSORS and is smaller than cbrt(NUMBER_OF_PROCESSORS)
5092  // - b the largest integer that divides NUMBER_OF_PROCESSORS/a and is smaller than sqrt(NUMBER_OF_PROCESSORS/a)
5093  unsigned b;
5094  for (unsigned n = floor(cbrt(NUMBER_OF_PROCESSORS));n>0; n--) {
5095  if (NUMBER_OF_PROCESSORS % n == 0) {
5096  a = n;
5097  break;
5098  }
5099  }
5100  for (unsigned n = floor(sqrt(NUMBER_OF_PROCESSORS/a));n>0; n--) {
5101  if (NUMBER_OF_PROCESSORS % (n*a) == 0) {
5102  b = n;
5103  break;
5104  }
5105  }
5107 }
5108 
5109 
5115 std::vector<unsigned> DPMBase::getNumberOfDomains()
5116 {
5117  return numberOfDomains_;
5118 }
5119 
5125 {
5127 }
5128 
5130 {
5131  return vtkWriter_;
5132 }
5133 
5141 {
5142  Vec3D meanVelocity = getTotalMomentum() / getTotalMass();
5143 
5144  //correct the mean velocity to zero
5145  for (auto& p : particleHandler)
5146  {
5147  p->addVelocity(-meanVelocity);
5148  }
5149 }
5150 
5158 {
5159  Vec3D V_mean;
5160  Mdouble Ek_mean_goal = 0, Ek_fluct_factor = 0;
5161  RNG rng;
5162 
5163  //assign random velocity to each particle
5164  for (auto& p : particleHandler)
5165  {
5166  p->setVelocity(Vec3D(rng.getRandomNumber(-1, 1), rng.getRandomNumber(-1, 1), rng.getRandomNumber(-1, 1)));
5167  }
5168 
5169  //calculate the mean velocity in the system now
5170  Ek_mean_goal = 0.5 * getTotalMass() * V_mean_goal.getLengthSquared();
5171  V_mean = getTotalMomentum() / getTotalMass();
5172  //check if the user input mean kinetic energy is larger than the total kinetic energy input, then return error
5173  logger.assert_always(0.5 * getTotalMass() * V_mean_goal.getLengthSquared() < Ek_goal,
5174  "Too large mean velocity input, Kinetic energy from mean velocity part is larger than the "
5175  "total kinetic energy you want to set");
5176 
5177  //correct the mean velocity to zero
5178  for (auto& p : particleHandler)
5179  {
5180  p->addVelocity(-V_mean);
5181  }
5182 
5183  //set the new fluctuating velocity based on the goal fluctuating kinetic energy
5184  Ek_fluct_factor = std::sqrt((Ek_goal - Ek_mean_goal) / getKineticEnergy());
5185  for (auto& p : particleHandler)
5186  {
5187  p->setVelocity(Ek_fluct_factor * p->getVelocity());
5188  }
5189 
5190  //correct the mean velocity finally to the user set values
5191  V_mean = getTotalMomentum() / getTotalMass();
5192  for (auto& p : particleHandler)
5193  {
5194  p->addVelocity(V_mean_goal - V_mean);
5195  }
5196 
5197  //check the final mean velocity and kinetic energy
5198  logger(INFO, "In DPMBase::setMeanVelocityAndKineticEnergy,\nV_mean_final %\n Ek_final %",
5200 }
5201 
5206 {
5207  return (getXMax() - getXMin()) * (getYMax() - getYMin()) * (getZMax() - getZMin());
5208 }
5209 
5215 {
5216  Matrix3D F; //set the kinetic energy tensor, this is in terms of Sum(m*v^2)
5217  Vec3D J; //set the momentum tensor
5218 
5219  //calculate stress for kinetic part
5220  for (const auto& p : particleHandler)
5221  {
5222  F += Matrix3D::dyadic(p->getVelocity(), p->getVelocity()) * p->getMass();
5223  J += p->getVelocity() * p->getMass();
5224  }
5225 
5226  Matrix3D stressKinetic = F - Matrix3D::dyadic(J, J) / getTotalMass();
5227  stressKinetic /= getTotalVolume();
5228  return stressKinetic;
5229 }
5230 
5237 {
5238  //stress components calculation variables
5239  Matrix3D stressStatic;
5240 
5241  //calculate the static stress tensor based on all the interactions
5242  for (const auto i : interactionHandler)
5243  {
5244  stressStatic += Matrix3D::dyadic(i->getForce(), i->getNormal()) * i->getDistance();
5245  }
5246 
5247  stressStatic /= getTotalVolume();
5248  return stressStatic;
5249 }
5250 
5257 {
5258  return getKineticStress() + getStaticStress();
5259 }
5260 
5261 
5263 {
5264  //compute forces for all particles that are neither fixed or ghosts
5265  for (auto p : particleHandler)
5266  {
5267  if (!p->isFixed() && p->getPeriodicFromParticle() == nullptr)
5268  {
5269  //w->computeForces(p);
5271  }
5272  }
5273 }
5274 
5282 void DPMBase::setLogarithmicSaveCount(const Mdouble logarithmicSaveCountBase)
5283 {
5284  dataFile.setlogarithmicSaveCount(logarithmicSaveCountBase);
5285  fStatFile.setlogarithmicSaveCount(logarithmicSaveCountBase);
5286  restartFile.setlogarithmicSaveCount(logarithmicSaveCountBase);
5287  statFile.setlogarithmicSaveCount(logarithmicSaveCountBase);
5288  eneFile.setlogarithmicSaveCount(logarithmicSaveCountBase);
5289 }
5290 
5297 void DPMBase::handleParticleRemoval(unsigned int id)
5298 {
5299  for (auto w: wallHandler)
5300  {
5301  w->handleParticleRemoval(id);
5302  }
5303 }
5304 
5313 {
5314  for (auto w: wallHandler)
5315  {
5316  w->handleParticleAddition(id, p);
5317  }
5318 }
5319 
5320 // initialise static member variables
5321 volatile sig_atomic_t DPMBase::continueFlag_ = true;
5322 
virtual void computeExternalForces(BaseParticle *)
Computes the external forces, such as gravity, acting on particles.
Definition: DPMBase.cc:3074
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:129
Mdouble timeMax_
Stores the duration of the simulation.
Definition: DPMBase.h:1286
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:772
void setFileCounter(unsigned fileCounter)
Definition: BaseVTKWriter.h:61
virtual 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:833
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 handleParticleRemoval(unsigned int id)
Handles the removal of particles from the particleHandler.
Definition: DPMBase.cc:5297
virtual void actionsBeforeTimeLoop()
A virtual function. Allows one to carry out any operations before the start of the time loop...
Definition: DPMBase.cc:1660
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:1311
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:3409
FileType getWallsWriteVTK() const
Returns whether walls are written in a VTK file.
Definition: DPMBase.cc:965
virtual void hGridUpdateParticle(BaseParticle *obj UNUSED)
Definition: DPMBase.cc:1695
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:1156
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:445
void setWallsWriteVTK(FileType writeWallsVTK)
Sets whether walls are written into a VTK file.
Definition: DPMBase.cc:908
unsigned getFileCounter() const
Definition: BaseVTKWriter.h:56
void write(std::ostream &os) const
void solve()
The work horse of the code.
Definition: DPMBase.cc:4002
virtual void hGridInsertParticle(BaseParticle *obj UNUSED)
Definition: DPMBase.cc:1688
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:134
The DPMBase header includes quite a few header files, defining all the handlers, which are essential...
Definition: DPMBase.h:75
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:516
Mdouble getRotationalEnergy() const
JMFT Returns the global rotational energy stored in the system.
Definition: DPMBase.cc:1568
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:5124
void set(BaseParticle *particleToCopy, unsigned int maxFailed, Vec3D posMin, Vec3D posMax, Vec3D velMin, Vec3D velMax)
Sets the properties of the InsertionBoundary for mutliple different particle types.
bool getWriteVTK() const
virtual void integrateAfterForceComputation()
Update particles' and walls' positions and velocities after force computation.
Definition: DPMBase.cc:3219
void initialiseMercuryMPITypes(const SpeciesHandler &speciesHandler)
Creates the MPI types required for communication of Mercury data through the MPI interface.
Definition: MpiContainer.cc:74
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:4864
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:312
A spherical particle is the most simple particle used in MercuryDPM.
void addParticle(BaseParticle *particle)
Initialises a single particle which is added during the simulation.
Definition: Domain.cc:1610
std::string name_
the name of the problem, used, e.g., for the output files
Definition: DPMBase.h:1364
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:1652
int getXBallsColourMode() const
Get the xballs colour mode (CMode).
Definition: DPMBase.cc:1301
Matrix3D getStaticStress() const
Calculate the static stress tensor in the system averaged over the whole volume.
Definition: DPMBase.cc:5236
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:2119
void setTimeMax(Mdouble newTMax)
Sets a new value for the maximum simulation duration.
Definition: DPMBase.cc:870
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:206
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:281
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:1421
bool readArguments(int argc, char *argv[])
Can interpret main function input arguments that are passed by the driver codes.
Definition: DPMBase.cc:4224
FileType getWriteVTK() const
FileType getFileType() const
Gets the file type e.g. NOFILE, ONEFILE and MULTIPLE FILES. File::fileType_.
Definition: File.cc:207
bool mpiInsertParticleCheck(BaseParticle *P)
Function that checks if the mpi particle should really be inserted by the current domain...
Definition: DPMBase.cc:1723
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:343
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
Definition of different loggers with certain modules. A user can define its own custom logger here...
Vec3D max_
Definition: DPMBase.h:1266
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:1250
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:931
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:362
void splitDomain(DomainSplit domainSplit)
Definition: DPMBase.cc:5058
void setMax(const Vec3D &max)
Sets the maximum coordinates of the problem domain.
Definition: DPMBase.cc:1073
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:1025
virtual void processStatistics(bool)
Definition: DPMBase.cc:1909
Mdouble getZMax() const
If the length of the problem domain in z-direction is ZMax - ZMin, then getZMax() returns ZMax...
Definition: DPMBase.h:631
Mdouble getXMin() const
If the length of the problem domain in x-direction is XMax - XMin, then getXMin() returns XMin...
Definition: DPMBase.h:600
Vec3D getCentreOfMass() const
JMFT: Return the centre of mass of the system, excluding fixed particles.
Definition: DPMBase.cc:1605
void setLastSavedTimeStep(unsigned int lastSavedTimeStep)
Sets File::nextSavedTimeStep_.
Definition: File.cc:302
void readNextFStatFile()
Reads the next fstat file.
Definition: DPMBase.cc:2779
void setCounter(unsigned int counter)
Allows the user to set the file counter according to his need. Sets File::counter_.
Definition: File.cc:231
bool readSpeciesFromDataFile_
Determines if the last column of the data file is interpreted as the info parameter during restart...
Definition: DPMBase.h:1374
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:3919
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:1458
Vec3D getMin() const
Definition: DPMBase.h:637
void setParticlesWriteVTK(bool writeParticlesVTK)
Sets whether particles are written in a VTK file.
Definition: DPMBase.cc:934
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:170
virtual void initialiseStatistics()
Definition: DPMBase.cc:1876
void setLogarithmicSaveCount(Mdouble logarithmicSaveCountBase)
Sets File::logarithmicSaveCount_ for all files (ene, data, fstat, restart, stat)
Definition: DPMBase.cc:5282
const std::complex< Mdouble > i
Definition: ExtendedMath.h:51
void flushParticles(std::set< BaseParticle * > &toBeDeletedList)
Particles that are going to be deleted from the simulation are flushed out of the communication bound...
Definition: Domain.cc:1698
WallVTKWriter wallVTKWriter_
Definition: DPMBase.h:1327
Mdouble getYMin() const
If the length of the problem domain in y-direction is YMax - YMin, then getYMin() returns YMin...
Definition: DPMBase.h:613
void setXBallsColourMode(int newCMode)
Set the xballs output mode.
Definition: DPMBase.cc:1291
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:1462
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:3773
Mdouble getSumOfInteractionRadii(const BaseParticle *particle) const
returns the sum of the radii plus the interactionDistance
Definition: BaseParticle.h:379
const std::string & getName() const
Returns the name of the file. Does not allow to change it though.
Definition: DPMBase.cc:397
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:1329
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:1359
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:1439
virtual ~DPMBase()
virtual destructor
Definition: DPMBase.cc:295
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:1869
void writePythonFileForVTKVisualisation() const
Definition: DPMBase.cc:2155
Mdouble getInteractionDistance()
Gets the interaction distance of the domain handler.
virtual void handleParticleAddition(unsigned int id, BaseParticle *p)
Definition: DPMBase.cc:5312
BaseParticle * getPeriodicFromParticle() const
Returns the 'original' particle this one's a periodic copy of.
Definition: BaseParticle.h:341
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:1420
void setZMax(Mdouble newZMax)
Sets the value of ZMax, the upper bound of the problem domain in the z-direction. ...
Definition: DPMBase.cc:1208
virtual void computeWallForces(BaseWall *w)
Definition: DPMBase.cc:5262
static void signalHandler(int signal)
signal handler function.
Definition: DPMBase.cc:4976
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:1394
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:1260
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:311
void tic()
This is like a start button of a stopwatch. Assigns the variable start with the current number of clo...
Definition: MercuryTime.h:55
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)
Fill a certain domain with particles.
void synchroniseParticle(BaseParticle *, unsigned fromProcessor=0)
Definition: DPMBase.cc:4945
virtual BaseParticle * copy() const =0
Particle copy method. It calls to copy constructor of this Particle, useful for polymorphism.
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:1513
int getNumberOfOMPThreads() const
Definition: DPMBase.cc:1277
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:3880
Species< LinearViscoelasticNormalSpecies, SlidingFrictionSpecies > LinearViscoelasticSlidingFrictionSpecies
void setSystemDimensions(unsigned int newDim)
Sets the system dimensionality.
Definition: DPMBase.cc:1408
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 setSoftStop()
function for setting sigaction constructor.
Definition: DPMBase.cc:5002
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:5157
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:1374
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:2268
virtual bool getHGridUpdateEachTimeStep() const
Definition: DPMBase.cc:1709
bool rotation_
A flag to turn on/off particle rotation. true will enable particle rotation. false will disable parti...
Definition: DPMBase.h:1311
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:976
#define MERCURY_DEPRECATED
Definition: GeneralDefine.h:37
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:1828
void addNewParticles()
unsigned int numberOfTimeSteps_
Stores the number of time steps.
Definition: DPMBase.h:1276
bool openWrite(unsigned)
First sets openmode to write (and append in some cases), then calls open().
Definition: File.cc:381
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:1298
void initialise()
Initialises the communication list vectors as they can not be determined on compile time...
void actionsAfterParticleGhostUpdate()
Calls the method actionsAfterParticleGhostUpdate of every wall in the handler.
Definition: WallHandler.cc:479
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:607
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:5256
unsigned int getSaveCount() const
Gets File::saveCount_.
Definition: File.cc:255
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:407
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:2563
virtual void writeEneTimeStep(std::ostream &os) const
Write the global kinetic, potential energy, etc. in the system.
Definition: DPMBase.cc:2095
void forceWriteOutputFiles()
Writes output files immediately, even if the current time step was not meant to be written...
Definition: DPMBase.cc:3859
void setYMax(Mdouble newYMax)
Sets the value of YMax, the upper bound of the problem domain in the y-direction. ...
Definition: DPMBase.cc:1182
unsigned int systemDimensions_
The dimensions of the simulation i.e. 2D or 3D.
Definition: DPMBase.h:1245
void readAndAddObject(std::istream &is) final
Reads BaseBoundary into the BoundaryHandler from restart data.
void writeFStatFile()
Definition: DPMBase.cc:2883
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:428
void setDomain(const Vec3D &min, const Vec3D &max)
Sets the minimum coordinates of the problem domain.
Definition: DPMBase.cc:1089
This is a class that generates random numbers i.e. named the Random Number Generator (RNG)...
Definition: RNG.h:52
MERCURY_DEPRECATED[[noreturn]] void logWriteAndDie(const std::string &module, std::string message)
todo strcmp relies on this, should be changed to more modern version
Definition: DPMBase.cc:79
Mdouble getElasticEnergy() const
Returns the global elastic energy within the system.
Definition: DPMBase.cc:1521
Matrix3D getKineticStress() const
Calculate the kinetic stress tensor in the system averaged over the whole volume. ...
Definition: DPMBase.cc:5214
Mdouble toc()
This is like a stop button of a stopwatch. Assigns the variable finish to the current value of ticks ...
Definition: MercuryTime.h:66
Mdouble getGravitationalEnergy() const
Returns the global gravitational potential energy stored in the system.
Definition: DPMBase.cc:1552
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:625
unsigned int getCounter() const
In case of multiple files, File::getCounter() returns the the number (FILE::Counter_) of the next fil...
Definition: File.cc:223
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:1436
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:2382
MERCURY_DEPRECATED File & getDataFile()
The non const version. Allows one to edit the File::dataFile.
Definition: DPMBase.cc:303
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:3530
void setlogarithmicSaveCount(const Mdouble logarithmicSaveCountBase)
the function to set the user input base of logarithmic saving count
Definition: File.cc:283
virtual void hGridActionsBeforeIntegration()
This function has to be called before integrateBeforeForceComputation.
Definition: DPMBase.cc:1931
std::fstream & getFstream()
Allows to access the member variable File::fstream_.
Definition: File.cc:153
File fStatFile
An instance of class File to handle in- and output into a .fstat file.
Definition: DPMBase.h:1441
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:1501
virtual bool readNextArgument(int &i, int argc, char *argv[])
Interprets the i^th command-line argument.
Definition: DPMBase.cc:4349
static volatile sig_atomic_t continueFlag_
Stores whether code should be stopped.
Definition: DPMBase.h:1237
virtual bool checkParticleForInteraction(const BaseParticle &P)
Checks whether a particle P has any interaction with walls or other particles.
Definition: DPMBase.cc:4638
void cleanCommunicationLists()
Removes nullptrs from boundaryParticleList_ and boundaryParticleListNeighbour_.
Definition: Domain.cc:1742
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:4963
bool getRestarted() const
Returns the flag denoting if the simulation was restarted or not.
Definition: DPMBase.cc:1484
Mdouble getMass() const
Returns the particle's mass.
Definition: BaseParticle.h:322
bool writeSuperquadricParticlesVTK_
Definition: DPMBase.h:1323
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:423
CGHandler cgHandler
Object of the class cgHandler.
Definition: DPMBase.h:1431
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:821
int xBallsColourMode_
XBalls is a package to view the particle data. As an alternative MercuryDPM also supports ParaView...
Definition: DPMBase.h:1339
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:1410
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:1415
ParticleVtkWriter * vtkWriter_
Definition: DPMBase.h:1325
BoundaryVTKWriter boundaryVTKWriter_
Definition: DPMBase.h:1331
void fillDomainWithParticles(unsigned N=50)
Inserts particles in the whole domain.
Definition: DPMBase.cc:2899
virtual Mdouble getInfo(const BaseParticle &P) const
A virtual function that returns some user-specified information about a particle. ...
Definition: DPMBase.cc:1633
Mdouble xBallsScale_
sets the xballs argument scale (see XBalls/xballs.txt)
Definition: DPMBase.h:1349
void writeEneFile()
Definition: DPMBase.cc:2873
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:1001
ParticleHandler particleHandler
An object of the class ParticleHandler, contains the pointers to all the particles created...
Definition: DPMBase.h:1395
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:317
int nToWrite_
number of elements to write to a screen
Definition: DPMBase.h:1379
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:1535
Mdouble & x()
RW reference to X.
Definition: Vector.h:344
Mdouble time_
Stores the current simulation time.
Definition: DPMBase.h:1271
void setRestartVersion(std::string newRV)
Sets restart_version.
Definition: DPMBase.cc:1475
void setMin(const Vec3D &min)
Sets the minimum coordinates of the problem domain.
Definition: DPMBase.cc:1109
void setSaveCount(unsigned int saveCount)
Sets File::saveCount_ for all files (ene, data, fstat, restart, stat)
Definition: DPMBase.cc:406
void closeFiles()
Closes all files (ene, data, fstat, restart, stat) that were opened to read or write.
Definition: DPMBase.cc:500
Mdouble getTotalVolume() const
Get the total volume of the cuboid system.
Definition: DPMBase.cc:5205
Vec3D getGravity() const
Returns the gravitational acceleration.
Definition: DPMBase.cc:1382
void setOpenMode(std::fstream::openmode openMode)
Sets File::openMode_ for all files (ene, data, fstat, restart, stat)
Definition: DPMBase.cc:487
int getNToWrite() const
get the number of elements to write to the
Definition: DPMBase.cc:861
void checkSettings()
Checks if the essentials are set properly to go ahead with solving the problem.
Definition: DPMBase.cc:3830
#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:1667
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:398
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:3302
Mdouble getRadius() const
Returns the particle's radius.
Definition: BaseParticle.h:348
bool mpiIsInCommunicationZone(BaseParticle *particle)
Checks if the position of the particle is in an mpi communication zone or not.
Definition: DPMBase.cc:1755
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:732
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:2600
void setInteractionsWriteVTK(bool)
Sets whether interactions are written into a VTK file.
Definition: DPMBase.cc:924
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:1049
virtual void actionsAfterSolve()
A virtual function which allows to define operations to be executed after the solve().
Definition: DPMBase.cc:1862
void setSaveCount(unsigned int saveCount)
Sets File::saveCount_.
Definition: File.cc:273
void clear() override
Empties the whole BaseHandler by removing all Objects and setting all other variables to 0...
int numberOfOMPThreads_
Definition: DPMBase.h:1240
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:1961
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:1385
InteractionHandler interactionHandler
An object of the class InteractionHandler.
Definition: DPMBase.h:1425
FileType writeWallsVTK_
A flag to turn on/off the vtk writer for walls.
Definition: DPMBase.h:1316
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:1248
const int getSVNRevision()
virtual void checkInteractions(InteractionHandler *interactionHandler, unsigned int timeStamp)
Check if all interactions are valid.
Definition: BaseWall.h:219
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:5018
Vec3D getMax() const
Definition: DPMBase.h:643
virtual void computeAdditionalForces()
A virtual function which allows to define operations to be executed prior to the OMP force collect...
Definition: DPMBase.h:1068
void setXBallsAdditionalArguments(std::string newXBArgs)
Set the additional arguments for xballs.
Definition: DPMBase.cc:1338
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:1951
Mdouble getYMax() const
If the length of the problem domain in y-direction is YMax - YMin, then getYMax() returns XMax...
Definition: DPMBase.h:619
Mdouble Y
Definition: Vector.h:65
Mdouble xBallsVectorScale_
sets the xballs argument vscale (see XBalls/xballs.txt)
Definition: DPMBase.h:1344
virtual void actionsBeforeTimeStep()
A virtual function which allows to define operations to be executed before the new time step...
Definition: DPMBase.cc:1855
ParticleVtkWriter * getVtkWriter() const
Definition: DPMBase.cc:5129
BaseParticle * getSmallestParticleLocal() const
Gets a pointer to the smallest BaseParticle (by interactionRadius) in this ParticleHandler of the loc...
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:3019
Mdouble timeStep_
Stores the simulation time step.
Definition: DPMBase.h:1281
void setFileType(FileType fileType)
Sets the type of file needed to write into or read from. File::fileType_.
Definition: File.cc:215
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:1405
void setWriteVTK(FileType f)
RNG random
This is a random generator, often used for setting up the initial conditions etc...
Definition: DPMBase.h:1390
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:441
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:3393
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:1702
void setPosition(const Vec3D &position)
Sets the position of this BaseInteractable.
Mdouble round(const Mdouble value, unsigned precision)
Definition: Helpers.cc:598
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:5140
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:667
bool writeParticlesVTK_
A flag to turn on/off the vtk writer for particles.
Definition: DPMBase.h:1321
void sync()
Process all pending asynchronous communication requests before continuing.
Definition: MpiContainer.h:152
void setName(const std::string &name)
Allows to set the name of all the files (ene, data, fstat, restart, stat)
Definition: DPMBase.cc:420
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:1674
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:2864
virtual void writeEneHeader(std::ostream &os) const
Writes a header with a certain format for ENE file.
Definition: DPMBase.cc:2005
Domain * getCurrentDomain()
Gets the domain assigned to the processor.
void autoNumber()
The autoNumber() function calls three functions: setRunNumber(), readRunNumberFromFile() and incremen...
Definition: DPMBase.cc:536
Vec3D getTotalMomentum() const
JMFT: Return the total momentum of the system, excluding fixed particles.
Definition: DPMBase.cc:1615
std::vector< unsigned > getNumberOfDomains()
returns the number of domains
Definition: DPMBase.cc:5115
bool open()
Checks if the file stream fstream_ has any issues while opening. Alongside, it also increments the ne...
Definition: File.cc:347
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:4255
int getRunNumber() const
This returns the current value of the counter (runNumber_)
Definition: DPMBase.cc:614
Mdouble getZMin() const
If the length of the problem domain in z-direction is ZMax - ZMin, then getZMin() returns ZMin...
Definition: DPMBase.h:625
void setSuperquadricParticlesWriteVTK(bool writeSuperquadricParticlesVTK)
Definition: DPMBase.cc:948
void insertGhostParticle(BaseParticle *P)
This function inserts a particle in the mpi communication boundaries.
Definition: DPMBase.cc:1802
void gatherContactStatistics()
Definition: DPMBase.cc:1889
MERCURY_DEPRECATED File & getFStatFile()
The non const version. Allows to edit the File::fStatFile.
Definition: DPMBase.cc:319
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:1884
File statFile
An instance of class File to handle in- and output into a .stat file.
Definition: DPMBase.h:1456
void boundaryActionsBeforeTimeLoop()
Definition: