MercuryDPM  Alpha
 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-2014, The MercuryDPM Developers Team. All rights reserved.
2 //For the list of developers, see <http://www.MercuryDPM.org/Team>.
3 //
4 //Redistribution and use in source and binary forms, with or without
5 //modification, are permitted provided that the following conditions are met:
6 // * Redistributions of source code must retain the above copyright
7 // notice, this list of conditions and the following disclaimer.
8 // * Redistributions in binary form must reproduce the above copyright
9 // notice, this list of conditions and the following disclaimer in the
10 // documentation and/or other materials provided with the distribution.
11 // * Neither the name MercuryDPM nor the
12 // names of its contributors may be used to endorse or promote products
13 // derived from this software without specific prior written permission.
14 //
15 //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
16 //ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
17 //WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18 //DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
19 //DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
20 //(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
21 //LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
22 //ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23 //(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24 //SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 
26 
27 #include "DPMBase.h"
28 #include <iostream>
29 #include <iomanip>
30 #include <algorithm>
31 #include <cmath>
32 #include <fstream>
33 #include <cstdlib>
34 #include <limits>
35 #include <string>
36 #include <sstream>
37 #include <cstdio>
39 #include <cstring>
40 //This is only used to change the file permission of the xball script create, at some point this code may be moved from this file to a different file.
41 #include <sys/types.h>
42 #include <sys/stat.h>
45 
46 #include "Species/Species.h"
49 #include "Particles/BaseParticle.h"
50 
51 //This is part of this class and just separates out the stuff to do with xballs.
52 #include "CMakeDefinitions.h"
53 #include "DPMBaseXBalls.icc"
54 #include "Logger.h"
55 #include "Particles/BaseParticle.h"
56 #include "Walls/BaseWall.h"
57 #include "Walls/InfiniteWall.h"
59 
65 void logWriteAndDie(std::string module, std::string message)
66 {
67  std::cerr << "A fatal error has occured"
68  << "\n Module :" << module
69  << "\n Message :" << message << std::endl;
70 
71  std::exit(-1);
72 }
73 
78 std::ostream& operator<<(std::ostream& os, DPMBase& md)
79 {
80  md.write(os);
81  return os;
82 }
83 
88  : FilesAndRunNumber(other)
89 {
90  constructor();
91 }
92 
97 {
101  gravity_ = other.gravity_;
102  xMin_ = other.xMin_;
103  xMax_ = other.xMax_;
104  yMin_ = other.yMin_;
105  yMax_ = other.yMax_;
106  zMin_ = other.zMin_;
107  zMax_ = other.zMax_;
108  time_ = other.time_;
109  timeStep_ = other.timeStep_;
110  ntimeSteps_ = other.ntimeSteps_;
111  timeMax_ = other.timeMax_;
112  restartVersion_ = other.restartVersion_; //to read new and old restart data
113  restarted_ = other.restarted_; //to see if it was restarted or not
114  append_ = other.append_;
115  rotation_ = other.rotation_;
116  xBallsColourMode_ = other.xBallsColourMode_; // sets the xballs argument cmode (see xballs.txt)
117  xBallsVectorScale_ = other.xBallsVectorScale_; // sets the xballs argument vscale (see xballs.txt)
118  xBallsScale_ = other.xBallsScale_; // sets the xballs argument scale (see xballs.txt)
119  xBallsAdditionalArguments_ = other.xBallsAdditionalArguments_; // std::string where additional xballs argument can be specified (see xballs.txt)
122 #ifdef CONTACT_LIST_HGRID
123  possibleContactList=other.possibleContactList;
124 #endif
125  random = other.random;
126 
131  wallHandler.setDPMBase(this);
132 
135  wallHandler = other.wallHandler;
138 }
139 
144 {
145  constructor();
146 }
147 
152 {
153 
154 }
155 
160 void DPMBase::solve(int argc, char* argv[])
161 {
162  readArguments(argc, argv);
163  solve();
164 }
165 
170 {
171  return time_;
172 }
173 
177 unsigned int DPMBase::getNtimeSteps() const
178 {
179  return ntimeSteps_;
180 }
181 
186 {
187  Mdouble diff = time_ - time;
188  time_ = time;
189  //this sets the interaction timestamp, so each interaction has the right time
190  for (auto i : interactionHandler)
191  {
192  i->setTimeStamp(i->getTimeStamp() - diff);
193  }
194 }
195 
200 {
201  if (newTMax >= 0)
202  {
203  timeMax_ = newTMax;
204  }
205  else
206  {
207  logger(ERROR, "Error in setTimeMax, newTMax=%", newTMax);
208  }
209 }
210 
215 {
216  return timeMax_;
217 }
221 #ifdef CONTACT_LIST_HGRID
222 PossibleContactList& DPMBase::getPossibleContactList()
223 {
224  return possibleContactList;
225 }
226 #endif
227 
231 void DPMBase::setRotation(bool newRotFlag)
232 {
233  rotation_ = newRotFlag;
234 }
235 
240 {
241  return rotation_;
242 }
243 
248 {
249  writeWallsVTK_ = writeWallsVTK;
250 }
251 
255 void DPMBase::setParticlesWriteVTK(bool writeParticlesVTK)
256 {
257  writeParticlesVTK_ = writeParticlesVTK;
258 }
259 
264 {
265  return writeWallsVTK_;
266 }
267 
272 {
273  return writeParticlesVTK_;
274 }
275 
280 {
281  return xMin_;
282 }
283 
288 {
289  return xMax_;
290 }
291 
296 {
297  return yMin_;
298 }
299 
304 {
305  return yMax_;
306 }
307 
312 {
313  return zMin_;
314 }
315 
320 {
321  return zMax_;
322 }
323 
328 {
329  return addL_;
330 }
331 
336 {
337  return addS_;
338 }
339 
343 void DPMBase::setAddSmall(int addS)
344 {
345  addS_ = addS;
346 }
347 
351 void DPMBase::setAddLarge(int addL)
352 {
353  addL_ = addL;
354 }
355 
360 {
361  if (xMin > getXMax())
362  {
363  logger(WARN, "Warning in setXMin(%): value cannot be larger than xMax=%", xMin, getXMax());
364  }
365  xMin_ = xMin;
366 }
367 
372 {
373  if (yMin > getYMax())
374  {
375  logger(WARN, "Warning in setYMin(%): value cannot be larger than yMax=%", yMin, getYMax());
376  }
377  yMin_ = yMin;
378 }
379 
384 {
385  if (zMin > getZMax())
386  {
387  logger(WARN, "Warning in setZMin(%): value cannot be larger than zMax=%", zMin, getZMax());
388  }
389  zMin_ = zMin;
390 }
391 
396 {
397  setXMax(max.X);
398  setYMax(max.Y);
399  setZMax(max.Z);
400 }
401 
406 {
407  setXMin(min.X);
408  setYMin(min.Y);
409  setZMin(min.Z);
410 }
411 
416 {
417  if (xMax < getXMin())
418  {
419  logger(WARN, "Warning in setXMax(%): value cannot be smaller than xMin=%", xMax, getXMin());
420  }
421  xMax_ = xMax;
422 }
423 
428 {
429  if (yMax < getYMin())
430  {
431  logger(WARN, "Warning in setYMax(%): value cannot be smaller than yMin=%", yMax, getYMin());
432  }
433  yMax_ = yMax;
434 }
435 
440 {
441  if (zMax < getZMin())
442  {
443  logger(WARN, "Warning in setZMax(%): value cannot be smaller than zMin=%", zMax, getZMin());
444  }
445  zMax_ = zMax;
446 }
447 
452 {
453  if (timeStep >= 0.0)
454  {
455  timeStep_ = timeStep;
456  } else
457  {
458  logger(ERROR, "Error in setTimeStep, proposed time step %", timeStep);
459  }
460 }
461 
466 {
467  return timeStep_;
468 }
469 
474 {
475  xBallsColourMode_ = newCMode;
476 }
477 
482 {
483  return xBallsColourMode_;
484 }
485 
489 void DPMBase::setXBallsVectorScale(double newVScale)
490 {
491  xBallsVectorScale_ = newVScale;
492 }
493 
498 {
499  return xBallsVectorScale_;
500 }
501 
505 void DPMBase::setXBallsAdditionalArguments(std::string newXBArgs)
506 {
507  xBallsAdditionalArguments_ = newXBArgs.c_str();
508 }
509 
514 {
516 }
517 
522 {
523  xBallsScale_ = newScale;
524 }
525 
530 {
531  return xBallsScale_;
532 }
533 
537 void DPMBase::setGravity(Vec3D newGravity)
538 {
539  gravity_ = newGravity;
540 }
541 
546 {
547  return gravity_;
548 }
549 
553 void DPMBase::setDimension(unsigned int newDim)
554 {
555  setSystemDimensions(newDim);
556  setParticleDimensions(newDim);
557 }
558 
562 void DPMBase::setSystemDimensions(unsigned int newDim)
563 {
564  if (newDim >= 1 && newDim <= 3)
565  systemDimensions_ = newDim;
566  else
567  {
568  logger(ERROR,"Error in setSystemDimensions, input argument is %", newDim);
569  }
570 }
571 
575 unsigned int DPMBase::getSystemDimensions() const
576 {
577  return systemDimensions_;
578 }
579 
583 void DPMBase::setParticleDimensions(unsigned int particleDimensions)
584 {
585  if (particleDimensions >= 1 && particleDimensions <= 3)
586  {
587  particleDimensions_ = particleDimensions;
589  } else
590  {
591  logger(WARN, "Error in setParticleDimensions, input argument is %", particleDimensions);
592  }
593 }
594 
599 unsigned int DPMBase::getParticleDimensions() const
600 {
601  return particleDimensions_;
602 }
603 
608 std::string DPMBase::getRestartVersion() const
609 {
610  return restartVersion_;
611 }
612 
617 void DPMBase::setRestartVersion(std::string newRV)
618 {
619  restartVersion_ = newRV;
620 }
621 
627 {
628  return restarted_;
629 }
630 
634 void DPMBase::setRestarted(bool newRestartedFlag)
635 {
636  restarted_ = newRestartedFlag;
637  //setAppend(new_);
638 }
639 
643 bool DPMBase::getAppend() const
644 {
645  return append_;
646 }
647 
651 void DPMBase::setAppend(bool newAppendFlag)
652 {
653  append_ = newAppendFlag;
654 }
655 
660 {
661  Mdouble elasticEnergy = 0.0;
662  for (const BaseInteraction* c : interactionHandler)
663  {
664  elasticEnergy += c->getElasticEnergy();
665  }
666  return elasticEnergy;
667 }
668 
673 {
674  Mdouble kineticEnergy = 0;
675  for (const BaseParticle* const p : particleHandler)
676  {
677  if (!(p->isFixed()))
678  {
679  kineticEnergy += .5 * p->getMass() * p->getVelocity().getLengthSquared();
680  }
681  }
682  return kineticEnergy;
683 }
684 
685 //finding the global **potential** energy
687 {
688  Mdouble gravitationalEnergy = 0;
689  for (const BaseParticle* const p : particleHandler)
690  {
691  if (!(p->isFixed()))
692  {
693  gravitationalEnergy += p->getMass() * Vec3D::dot((-getGravity()), p->getPosition());
694  }
695  }
696  return gravitationalEnergy;
697 
698 }
699 
703 double DPMBase::getInfo(const BaseParticle& p) const
704 {
705  return p.getSpecies()->getId(); // was getIndex()
706 }
707 
713 bool DPMBase::areInContact(const BaseParticle* pI, const BaseParticle* pJ) const
714 {
715  return (pI != pJ && Vec3D::getDistanceSquared(pI->getPosition(), pJ->getPosition()) <
717 }
718 
723 {
724 }
725 
730 {
731 }
732 
737 {
738 }
739 
744 {
745 }
746 
751 {
752 }
753 
758 {
759 }
760 
765 {
766 }
767 
772 {
773  return true;
774 }
775 
780 {
781 }
782 
787 {
788 }
789 
794 {
795 }
796 
801 {
802 }
803 
808 {
809 }
810 
812 {
814  {
815  c->gatherContactStatistics();
816  }
817 }
818 
822 void DPMBase::gatherContactStatistics(unsigned int index1 UNUSED, int index2 UNUSED, Vec3D Contact UNUSED,
823  Mdouble delta UNUSED, Mdouble ctheta UNUSED, Mdouble fdotn UNUSED,
824  Mdouble fdott UNUSED, Vec3D P1_P2_normal_ UNUSED, Vec3D P1_P2_tangential UNUSED)
825 {
826 }
827 
832 {
833 }
834 
839 {
840 }
841 
846 {
847 }
848 
853 {
854 }
855 
860 {
861 }
862 
867 {
868  for (auto it = particleHandler.begin(); (*it) != i; ++it)
869  {
870  computeInternalForces(i, *it);
871  }
872 }
873 
877 void DPMBase::setFixedParticles(unsigned int n)
878 {
879  for (unsigned int i = 0; i < std::min(particleHandler.getNumberOfObjects(), n); i++)
881 }
882 
886 void DPMBase::printTime() const
887 {
888  std::cout << "t=" << std::setprecision(3) << std::left << std::setw(6) << getTime()
889  << ", tmax=" << std::setprecision(3) << std::left << std::setw(6) << getTimeMax()
890  << std::endl;
891  std::cout.flush();
892 }
893 
898 {
899  return true;
900 }
901 
906 {
907  //These are the particle parameters like dissipation etc..
908  setRestarted(false);
909 
910  //This sets the maximum number of particles
915  wallHandler.setDPMBase(this);
917 
921 
922  // Set gravitational acceleration
923  gravity_ = Vec3D(0.0, -9.8, 0.0);
924 
925  //This is the parameter of the numerical part
926  time_ = 0.0;
927  timeStep_ = 0.0; // if dt is not user-specified, this is set in actionsBeforeTimeLoop()
928  ntimeSteps_ = 0;
929  setSaveCount(20);
930  timeMax_ = 1.0;
931 
932  //This sets the default xballs domain
933  xMin_ = 0.0;
934  xMax_ = 0.01;
935  yMin_ = 0.0;
936  yMax_ = 0.01;
937  zMin_ = 0.0;
938  zMax_ = 0.0;
939 
940  //sets the default write particles data in VTK format flag to false
941  writeParticlesVTK_ = false;
943 
945 
946  setName("");
947 
948  //Default mode is energy with no scale of the vectors
949  xBallsColourMode_ = 0;
950  xBallsVectorScale_ = -1;
951  xBallsScale_ = -1;
953  setAppend(false);
954 //TODO
955  //The default random seed is 0
957 
958 #ifdef DEBUG_OUTPUT
959  std::cerr << "MD problem constructor finished " << std::endl;
960 #endif
961 }
962 
967 {
968 }
969 
973 void DPMBase::writeEneHeader(std::ostream& os) const
974 {
975  //only write if we don't restart
976  if (getAppend())
977  return;
978 
980  long width = os.precision() + 6;
981  os << std::setw(width) << "t" << " " << std::setw(width)
982  << "ene_gra" << " " << std::setw(width)
983  << "ene_kin" << " " << std::setw(width)
984  << "ene_rot" << " " << std::setw(width)
985  << "ene_ela" << " " << std::setw(width)
986  << "X_COM" << " " << std::setw(width)
987  << "Y_COM" << " " << std::setw(width)
988  << "Z_COM" << std::endl;
989 }
990 
994 void DPMBase::writeFstatHeader(std::ostream& os) const
995 {
996 
997  // line #1: time, volume fraction
998  // line #2: wall box: wx0, wy0, wz0, wx1, wy1, wz1
999  // line #3: radii-min-max & moments: rad_min, rad_max, r1, r2, r3, r4
1000  os << "#"
1001  << " " << getTime()
1002  << " " << 0
1003  << std::endl;
1004  os << "#"
1005  << " " << getXMin()
1006  << " " << getYMin()
1007  << " " << getZMin()
1008  << " " << getXMax()
1009  << " " << getYMax()
1010  << " " << getZMax()
1011  << std::endl;
1012  os << "#"
1013  << " ";
1015  {
1017  } else
1018  {
1019  os << std::numeric_limits<double>::quiet_NaN();
1020  }
1021  os << " ";
1023  {
1025  } else
1026  {
1027  os << std::numeric_limits<double>::quiet_NaN();
1028  }
1029  os << " " << 0
1030  << " " << 0
1031  << " " << 0
1032  << " " << 0
1033  << std::endl;
1034  //B: write data
1036  {
1037  c->writeToFStat(os);
1038  }
1039 }
1040 
1044 void DPMBase::writeEneTimestep(std::ostream& os) const
1045 {
1046  Mdouble ene_kin = 0, ene_elastic = 0, ene_rot = 0, ene_gra = 0, mass_sum = 0, x_masslength = 0, y_masslength = 0, z_masslength = 0;
1047 
1048  for (BaseParticle* p : particleHandler)
1049  if (!p->isFixed())
1050  {
1051  ene_kin += .5 * p->getMass() * p->getVelocity().getLengthSquared();
1052  ene_rot += .5 * p->getInertia() * p->getAngularVelocity().getLengthSquared();
1053  ene_gra -= p->getMass() * Vec3D::dot(getGravity(), p->getPosition());
1054  mass_sum += p->getMass();
1055  x_masslength += p->getMass() * p->getPosition().X;
1056  y_masslength += p->getMass() * p->getPosition().Y;
1057  z_masslength += p->getMass() * p->getPosition().Z;
1058  } //end for loop over Particles
1059 
1060  ene_elastic = getElasticEnergy();
1061 
1063  long width = os.precision() + 6;
1064  os << std::setw(width) << getTime()
1065  << " " << std::setw(width) << ene_gra
1066  << " " << std::setw(width) << ene_kin
1067  << " " << std::setw(width) << ene_rot
1068  << " " << std::setw(width) << ene_elastic
1069  << " " << std::setw(width)
1070  << (mass_sum != 0.0 ? x_masslength / mass_sum : std::numeric_limits<double>::quiet_NaN())
1071  << " " << std::setw(width)
1072  << (mass_sum != 0.0 ? y_masslength / mass_sum : std::numeric_limits<double>::quiet_NaN())
1073  << " " << std::setw(width)
1074  << (mass_sum != 0.0 ? z_masslength / mass_sum : std::numeric_limits<double>::quiet_NaN()) << std::endl;
1075 
1076  //sliding = sticking = 0;
1077 }
1078 
1079 
1080 void DPMBase::writeVTK() const {
1081  if ((getWallsWriteVTK() == FileType::ONE_FILE && getTime() == 0.0)
1084  {
1086  }
1087  if (getParticlesWriteVTK() == true)
1088  {
1090  }
1092 
1094  {
1095  logger(INFO,"Writing python script for paraview visualisation");
1096  helpers::writeToFile(getName()+".paraview.py",
1097  "#script to visualise the output of data2pvd of MercuryDPM in paraview.\n"
1098  "#usage: change the path below to your own path, open paraview\n"
1099  "#Tools->Python Shell->Run Script->VisualisationScript.py\n"
1100  "#or run paraview --script=VisualisationScript.py \n"
1101  "\n"
1102  "try: paraview.simple\n"
1103  "except: from paraview.simple import *\n"
1104  "import glob\n"
1105  //"import os\n"
1106  //"os.chdir(\"/Users/weinhartt/Mercury/Alpha/cmake-build-debug/Drivers/USER/Thomas/Mainz\")\n"
1107  "\n"
1108  "particles = XMLUnstructuredGridReader(FileName=glob.glob('./"
1109  + getName() + "Particle_*.vtu'))\n"
1110  "\n"
1111  "walls = XMLUnstructuredGridReader(FileName=glob.glob('./"
1112  + getName() + "Wall_*.vtu'))\n"
1113  "\n"
1114  "interactions = XMLUnstructuredGridReader(FileName=glob.glob('./"
1115  + getName() + "Interaction_*.vtu'))\n"
1116  "\n"
1117  "#Make the glyph, see visualisation guide on the docs\n"
1118  "glyphP = Glyph(particles)\n"
1119  "glyphP.GlyphType = 'Sphere'\n"
1120  "glyphP.Scalars = 'Radius'\n"
1121  "glyphP.Vectors = 'None'\n"
1122  "glyphP.ScaleMode = 'scalar'\n"
1123  "glyphP.ScaleFactor = 2\n"
1124  "glyphP.GlyphMode = 'All Points'\n"
1125  "\n"
1126  "#Make the glyph, see visualisation guide on the docs\n"
1127  "glyphI = Glyph(interactions)\n"
1128  "glyphI.GlyphType = 'Sphere'\n"
1129  "glyphI.Scalars = 'ContactRadius'\n"
1130  "glyphI.Vectors = 'None'\n"
1131  "glyphI.ScaleMode = 'scalar'\n"
1132  "glyphI.ScaleFactor = 10\n" //5 times too large
1133  "glyphI.GlyphMode = 'All Points'\n"
1134  "\n"
1135  "#show on screen\n"
1136  "Show(glyphI)\n"
1137  "Show(glyphP)\n"
1138  "Show(walls)\n"
1139  "Render()\n"
1140  "ResetCamera()");
1141  }
1142 }
1143 
1147 void DPMBase::outputXBallsData(std::ostream& os) const
1148 {
1149  //Set the correct formation based of dimension if the formation is not specified by the user
1150  unsigned int format;
1151  switch (getSystemDimensions())
1152  {
1153  case 2:
1154  format = 8;
1155  break;
1156  case 3:
1157  format = 14;
1158  break;
1159  default:
1160  std::cerr << "Unknown systemdimension" << std::endl;
1161  exit(-1);
1162  }
1163 
1164  // This outputs the location of walls and how many particles there are to file this is required by the xballs plotting
1165  if (format != 14) // dim = 1 or 2
1166  {
1168  << " " << getTime()
1169  << " " << getXMin()
1170  << " " << getYMin()
1171  << " " << getXMax()
1172  << " " << getYMax()
1173  << " " << std::endl;
1174  } else
1175  {
1176  //dim==3
1178  << " " << getTime()
1179  << " " << getXMin()
1180  << " " << getYMin()
1181  << " " << getZMin()
1182  << " " << getXMax()
1183  << " " << getYMax()
1184  << " " << getZMax()
1185  << " " << std::endl;
1186  }
1187  // This outputs the particle data
1188  for (unsigned int i = 0; i < particleHandler.getNumberOfObjects(); i++)
1189  {
1190  outputXBallsDataParticle(i, format, os);
1191  }
1192 #ifdef DEBUG_OUTPUT
1193  std::cerr << "Have output the properties of the problem to disk " << std::endl;
1194 #endif
1195 }
1196 
1202 bool DPMBase::readDataFile(const std::string fileName, unsigned int format)
1203 {
1204  std::string oldFileName = dataFile.getName();
1205  //This opens the file the data will be recalled from
1206  dataFile.setName(fileName);
1207  dataFile.open(std::fstream::in);
1208  if (!dataFile.getFstream().is_open() || dataFile.getFstream().bad())
1209  {
1210  std::cout << "Loading data file " << fileName << " failed" << std::endl;
1211  return false;
1212  }
1213 
1214  //dataFile.getFileType() is ignored
1215  FileType fileTypeData = dataFile.getFileType();
1217  readNextDataFile(format);
1218  dataFile.setFileType(fileTypeData);
1219 
1220  dataFile.close();
1221  dataFile.setName(oldFileName);
1222 
1223  //std::cout << "Loaded data file " << filename << " (t=" << getTime() << ")" << std::endl;
1224  return true;
1225 }
1226 
1231 bool DPMBase::readParAndIniFiles(const std::string fileName)
1232 {
1233  //Opens the par.ini file
1234  std::fstream file;
1235  file.open(fileName, std::fstream::in);
1236  if (!file.is_open() || file.bad())
1237  {
1238  //std::cout << "Loading par.ini file " << filename << " failed" << std::endl;
1239  return false;
1240  }
1241 
1242  Mdouble doubleValue;
1243  int integerValue;
1244 
1245  // inputfile par.ini
1246  // line 1 =============================================================
1247  // Example: 1 1 0
1248  // 1: integer (0|1) switches from non-periodic to periodic
1249  // integer (5|6) does 2D integration only (y-coordinates fixed)
1250  // and switches from non-periodic to periodic
1251  // integer (11) uses a quarter system with circular b.c.
1252  file >> integerValue;
1253  //~ std::cout << "11" << integerValue << std::endl;
1254  if (integerValue == 0)
1255  {
1256  InfiniteWall w0;
1258  w0.set(Vec3D(-1, 0, 0), Vec3D(getXMin(), 0, 0));
1260  w0.set(Vec3D(1, 0, 0), Vec3D(getXMax(), 0, 0));
1262  w0.set(Vec3D(0, -1, 0), Vec3D(0, getYMin(), 0));
1264  w0.set(Vec3D(0, 1, 0), Vec3D(0, getYMax(), 0));
1266  w0.set(Vec3D(0, 0, -1), Vec3D(0, 0, getZMin()));
1268  w0.set(Vec3D(0, 0, 1), Vec3D(0, 0, getZMax()));
1270  } else if (integerValue == 1)
1271  {
1272  PeriodicBoundary b0;
1273  b0.set(Vec3D(1, 0, 0), getXMin(), getXMax());
1275  b0.set(Vec3D(0, 1, 0), getYMin(), getYMax());
1277  b0.set(Vec3D(0, 0, 1), getZMin(), getZMax());
1279  } else if (integerValue == 5)
1280  {
1281  InfiniteWall w0;
1283  w0.set(Vec3D(-1, 0, 0), Vec3D(-getXMin(), 0, 0));
1285  w0.set(Vec3D(1, 0, 0), Vec3D(getXMax(), 0, 0));
1287  w0.set(Vec3D(0, -1, 0), Vec3D(0, -getYMin(), 0));
1289  w0.set(Vec3D(0, 1, 0), Vec3D(0, getYMax(), 0));
1291 
1292  } else if (integerValue == 6)
1293  {
1294  PeriodicBoundary b0;
1295  b0.set(Vec3D(1, 0, 0), getXMin(), getXMax());
1297  b0.set(Vec3D(0, 1, 0), getYMin(), getYMax());
1299  b0.set(Vec3D(0, 0, 1), getZMin(), getZMax());
1301  } else
1302  {
1303  std::cerr << "Error in par.ini: line 1, value 1 is " << integerValue << std::endl;
1304  exit(-1);
1305  }
1306 
1307  // 2: integer (0|1) dont use | use the search pattern for linked cells
1308  file >> integerValue; //ignore
1309 
1310  // 3: real - gravity in z direction: positive points upwards
1311  file >> doubleValue;
1312  setGravity(Vec3D(0.0, 0.0, doubleValue));
1313 
1314  // line 2 =============================================================
1315  // Example: -10000 .5e-2
1316  // 1: time end of simulation - (negative resets start time to zero
1317  // and uses -time as end-time)
1318  file >> doubleValue;
1319  if (doubleValue < 0)
1320  setTime(0);
1321  setTimeMax(fabs(doubleValue));
1322 
1323  // 2: time-step of simulation
1324  file >> doubleValue;
1325  setTimeStep(doubleValue);
1326 
1327  // line 3 =============================================================
1328  // Example: 1e-1 100
1329  file >> doubleValue;
1330  if (doubleValue >= 0)
1331  {
1332  // 1: time-step for output on time-series protocoll file -> "ene"
1333  unsigned int savecount = static_cast<unsigned int>(round(doubleValue / getTimeStep()));
1334  setSaveCount(savecount);
1335 
1336  // 2: time-step for output on film (coordinate) file -> "c3d"
1337  // (fstat-output is coupled to c3d-output time-step)
1338  file >> doubleValue;
1339  savecount = static_cast<unsigned int>(round(doubleValue / getTimeStep()));
1340  dataFile.setSaveCount(savecount);
1341  fStatFile.setSaveCount(savecount);
1342  } else
1343  {
1344  // or: ---------------------------------------------------------------
1345  // 1: negative number is multiplied to the previous log-output-time
1346  // 2: requires initial log-output time
1347  // 3: negative number is multiplied to the previous film-output-time
1348  // 4: requires initial film-output time
1349  std::cerr << "Error in par.ini: line 3, value 1 is " << doubleValue << std::endl;
1350  exit(-1);
1351  }
1352 
1353  // line 4 =============================================================
1354  // Example: 2000 1e5 1e3 1e2
1355  // 1: particle density (mass=4/3*constants::pi*density*rad^3)
1356  file >> doubleValue;
1357 
1358  //clear species handler
1362 
1364  S->setDensity(doubleValue);
1365 
1366  // 2: linear spring constant
1367  file >> doubleValue;
1368  S->setStiffness(doubleValue);
1369 
1370  // 3: linear dashpot constant
1371  file >> doubleValue;
1372  S->setDissipation(doubleValue);
1373 
1374  // 4: background damping dashpot constant
1375  file >> doubleValue;
1376  if (doubleValue != 0.0)
1377  std::cerr << "Warning in par.ini: ignored background damping " << doubleValue << std::endl;
1378 
1379  // line 5 =============================================================
1380  // Example: 0 0
1381  // 1: growth rate: d(radius) = xgrow * dt
1382  file >> doubleValue;
1383  if (doubleValue != 0.0)
1384  std::cerr << "Warning in par.ini: ignored growth rate " << doubleValue << std::endl;
1385 
1386  // 2: target volume_fraction
1387  file >> doubleValue;
1388  if (doubleValue != 0.0)
1389  std::cerr << "Warning in par.ini: ignored target volume_fraction " << doubleValue << std::endl;
1390 
1391  file.close();
1392  //std::cout << "Loaded par.ini file " << filename << std::endl;
1393  return true;
1394 }
1395 
1402 {
1404  {
1405  while (true)// This true corresponds to the if s
1406  {
1408  //check if file exists and contains data
1409  int N;
1410  dataFile.getFstream() >> N;
1411  if (dataFile.getFstream().eof() || dataFile.getFstream().peek() == -1)
1412  {
1413  std::cout << "file " << dataFile.getName() << " not found" << std::endl;
1414  return false;
1415  }
1416  //check if tmin condition is satisfied
1417  Mdouble t;
1418  dataFile.getFstream() >> t;
1419  if (t > tMin)
1420  {
1421  //set_file_counter(get_file_counter()-1);
1422  return true;
1423  }
1424  if (verbose)
1425  std::cout << "Jumping file counter: " << dataFile.getCounter() << std::endl;
1426  }
1427  }
1428  return true;
1429 }
1430 
1434 bool DPMBase::readNextDataFile(unsigned int format)
1435 {
1436  dataFile.openNextFile(std::fstream::in);
1437  //fStatFile.openNextFile();
1438  //Set the correct formation based of dimension if the formation is not specified by the user
1439  if (format == 0)
1440  {
1441  switch (getSystemDimensions())
1442  {
1443  case 1:
1444  case 2:
1445  format = 8;
1446  break;
1447  case 3:
1448  format = 14;
1449  break;
1450  }
1451  //end case
1452 
1453  }
1454  //end if
1455 
1456  unsigned int N = 0;
1457  Mdouble dummy;
1458  dataFile.getFstream() >> N;
1459  //read first parameter and check if we reached the end of the file
1460  if (N == 0)
1461  {
1462  //std::cout << "Reached the end of the data file" << std::endl;
1463  return false;
1464  }
1465 
1466  //create particles if N is not fixed
1468  {
1470  {
1471  BaseParticle p;
1474  std::cout << "Warning in readNextDataFile: particleHandler is empty, "
1475  "so newly created particles will be of BaseParticle type." << std::endl;
1476  }
1478  {
1479  while (particleHandler.getNumberOfObjects() != N)
1480  {
1482  //std::cout << "S" << particleHandler.getNumberOfObjects() << std::endl;
1483  }
1484  } else
1485  {
1486  while (particleHandler.getNumberOfObjects() != N)
1487  {
1490  //std::cout << particleHandler.getNumberOfObjects() << N << std::endl;
1491  }
1492  }
1493  }
1494  //for (auto p: particleHandler)
1495  // std::cout << *p << std::endl;
1496 #ifdef DEBUG_OUTPUT
1497  std::cout << "Number of particles read from file "<<N << std::endl;
1498 #endif
1499 
1500  //read all other data available for the time step
1501  switch (format)
1502  {
1503  //This is the format_ that Stefan's and Vitaly's code saves in - note the axis rotation_
1504  case 7:
1505  {
1506  std::stringstream line;
1508  line >> dummy;
1510  line >> dummy;
1511  setTime(dummy);
1513  line >> dummy;
1514  setXMin(dummy);
1516  line >> dummy;
1517  setYMin(dummy);
1519  line >> dummy;
1520  setZMin(dummy);
1522  line >> dummy;
1523  setXMax(dummy);
1525  line >> dummy;
1526  setYMax(dummy);
1528  line >> dummy;
1529  setZMax(dummy);
1530  //std::cout << " time " << t << std::endl;
1531  Mdouble radius;
1532  Vec3D position, velocity;
1533  for (BaseParticle* p : particleHandler)
1534  {
1536  line >> position.X >> position.Z >> position.Y >> velocity.X >> velocity.Z >> velocity.Y >> radius
1537  >> dummy;
1538  p->setPosition(position);
1539  p->setVelocity(velocity);
1540  p->setOrientation(Vec3D(0.0, 0.0, 0.0));
1541  p->setAngularVelocity(Vec3D(0.0, 0.0, 0.0));
1542  p->setRadius(radius);
1543  }
1544  break;
1545  }
1546  //This is a 2D format_
1547  case 8:
1548  {
1549  std::stringstream line;
1551  line >> time_ >> xMin_ >> yMin_ >> xMax_ >> yMax_;
1552  zMin_ = 0.0;
1553  zMax_ = 0.0;
1554  Mdouble radius;
1555  Vec3D position, velocity, angle, angularVelocity;
1556  for (BaseParticle* p : particleHandler)
1557  {
1559  line >> position.X >> position.Y >> velocity.X >> velocity.Y >> radius >> angle.Z >> angularVelocity.Z
1560  >> dummy;
1561  p->setPosition(position);
1562  p->setVelocity(velocity);
1563  p->setOrientation(-angle);
1564  p->setAngularVelocity(-angularVelocity);
1565  p->setRadius(radius);
1566  } //end for all particles
1567  break;
1568  }
1569  //This is a 3D format_
1570  case 14:
1571  {
1572  std::stringstream line;
1574  line >> time_ >> xMin_ >> yMin_ >> zMin_ >> xMax_ >> yMax_ >> zMax_;
1575  Mdouble radius;
1576  Vec3D position, velocity, angle, angularVelocity;
1577  for (BaseParticle* p : particleHandler)
1578  {
1580  line >> position >> velocity >> radius >> angle >> angularVelocity;
1581  p->setPosition(position);
1582  p->setVelocity(velocity);
1583  p->setOrientation(angle);
1584  p->setAngularVelocity(angularVelocity);
1585  p->setRadius(radius);
1586  } //end read into existing particles
1587  logger(INFO, "read % particles", particleHandler.getNumberOfObjects());
1588  break;
1589  }
1590  //This is a 3D format_
1591  case 15:
1592  {
1593  std::stringstream line;
1595  line >> time_ >> xMin_ >> yMin_ >> zMin_ >> xMax_ >> yMax_ >> zMax_;
1596  Mdouble radius;
1597  Vec3D position, velocity, angle, angularVelocity;
1598  for (BaseParticle* p : particleHandler)
1599  {
1601  line >> position >> velocity >> radius >> angle >> angularVelocity >> dummy >> dummy;
1602  p->setPosition(position);
1603  p->setVelocity(velocity);
1604  p->setOrientation(angle);
1605  p->setAngularVelocity(angularVelocity);
1606  p->setRadius(radius);
1607  } //end for all particles
1608  break;
1609  }
1610  } //end switch statement
1612  return true;
1613 }
1614 
1619 {
1622 }
1623 
1629 bool DPMBase::readRestartFile(bool restarted)
1630 {
1631  if (restartFile.open(std::fstream::in))
1632  {
1634  restartFile.close();
1635  setRestarted(restarted);
1636  return true;
1637  } else {
1638  //std::cerr << restartFile.getName() << " could not be loaded" << std::endl;
1639  return false;
1640  }
1641 }
1642 
1647 int DPMBase::readRestartFile(std::string fileName)
1648 {
1649  restartFile.setName(fileName);
1650  return readRestartFile();
1651 }
1652 
1658 {
1659  //this is because the rough bottom allows overlapping fixed particles
1660  if (P1->isFixed() && P2->isFixed())
1661  {
1662  return;
1663  }
1664 
1666  {
1667  return;
1668  }
1669 
1670  BaseParticle* PI, * PJ;
1671  if (P1->getId() > P2->getId())
1672  {
1673  PI = P2;
1674  PJ = P1;
1675  } else
1676  {
1677  PI = P1;
1678  PJ = P2;
1679  }
1680 
1681  //if statement above ensures that the PI has the lower id than PJ
1682  std::vector<BaseInteraction*> interactions = PJ->getInteractionWith(PI, getTime() + getTimeStep(),
1684  for (auto i : interactions)
1685  {
1686  i->computeForce();
1687 
1688  PI->addForce(i->getForce());
1689  PJ->addForce(-i->getForce());
1690 
1691  if (getRotation())
1692  {
1693  PI->addTorque(i->getTorque() - Vec3D::cross(PI->getPosition() - i->getContactPoint(), i->getForce()));
1694  PJ->addTorque(-i->getTorque() + Vec3D::cross(PJ->getPosition() - i->getContactPoint(), i->getForce()));
1695  }
1696  }
1697 }
1698 
1704 {
1705  if (!CI->isFixed())
1706  {
1707  // Gravitational force
1708  CI->addForce(getGravity() * CI->getMass());
1709  // Still calls this in compute External Forces.
1711  }
1712 }
1713 
1718 {
1719  //No need to compute interactions between periodic particle images and walls
1720  if (pI->getPeriodicFromParticle() != nullptr)
1721  return;
1722 
1723  for (BaseWall* w : wallHandler)
1724  {
1725  std::vector<BaseInteraction*> interactions = w->getInteractionWith(pI, getTime() + getTimeStep(),
1727 
1728  for (auto i : interactions)
1729  {
1730  i->computeForce();
1731 
1732  pI->addForce(i->getForce());
1733  w->addForce(-i->getForce());
1734 
1735  if (getRotation()) // getRotation() returns a boolean.
1736  {
1737  pI->addTorque(i->getTorque() - Vec3D::cross(pI->getPosition() - i->getContactPoint(), i->getForce()));
1739  w->addTorque(-i->getTorque() + Vec3D::cross(w->getPosition() - i->getContactPoint(), i->getForce()));
1740  }
1741  }
1742  }
1743 }
1744 
1749 {
1750  for_each(particleHandler.begin(), particleHandler.end(), [this](BaseParticle* p)
1751  {
1752  p->integrateBeforeForceComputation(getTime(), getTimeStep());
1753  });
1754  for_each(wallHandler.begin(), wallHandler.end(), [this](BaseWall* w)
1755  {
1756  w->integrateBeforeForceComputation(getTime(), getTimeStep());
1757  });
1758 }
1759 
1764 {
1765  for (BaseBoundary* b : boundaryHandler)
1766  {
1767  for (auto p = particleHandler.begin(); p != particleHandler.end(); ++p)
1768  {
1769  if (b->checkBoundaryAfterParticleMoved(*p, particleHandler))
1770  //Returns true if the particle is deleted
1771  --p;
1772  }
1773  }
1774 }
1775 
1780 {
1781  for_each(particleHandler.begin(), particleHandler.end(), [this](BaseParticle* p)
1782  {
1783  p->integrateAfterForceComputation(getTime(), getTimeStep());
1784  });
1785  for_each(wallHandler.begin(), wallHandler.end(), [this](BaseWall* w)
1786  {
1787  w->integrateAfterForceComputation(getTime(), getTimeStep());
1788  });
1789 }
1790 
1794 //void DPMBase::statisticsFromRestartData(const char *name)
1795 //{
1796 // ///todo{Check this whole function}
1797 // //This function loads all MD data
1798 // readRestartFile();
1799 //
1800 // //This creates the file statistics will be saved to
1801 // std::stringstream ss("");
1802 // ss << name << ".stat";
1803 // statFile.setName(ss.str());
1804 // statFile.setOpenMode(std::fstream::out);
1805 // statFile.open();
1806 //
1807 // // Sets up the initial conditions for the simulation
1808 // // setupInitialConditions();
1809 // // Setup the previous position arrays and mass of each particle.
1810 // computeParticleMasses();
1811 // // Other routines required to jump-start the simulation
1812 // actionsBeforeTimeLoop();
1813 // initialiseStatistics();
1814 // hGridActionsBeforeTimeLoop();
1815 // writeEneHeader(eneFile.getFstream());
1816 //
1817 // while (readDataFile(dataFile.getName().c_str()))
1818 // {
1819 // hGridActionsBeforeTimeLoop();
1820 // actionsBeforeTimeStep();
1821 // checkAndDuplicatePeriodicParticles();
1822 // hGridActionsBeforeTimeStep();
1827 // computeAllForces();
1828 // removeDuplicatePeriodicParticles();
1829 // actionsAfterTimeStep();
1830 // writeEneTimestep(eneFile.getFstream());
1831 // std::cout << std::setprecision(6) << getTime() << std::endl;
1832 // }
1833 //
1834 // dataFile.close();
1835 // statFile.close();
1836 //}
1837 
1842 {
1844  for (BaseParticle* const p : particleHandler)
1845  {
1846  p->setForce(Vec3D(0.0, 0.0, 0.0));
1847  p->setTorque(Vec3D(0.0, 0.0, 0.0));
1848  }
1849  for (BaseWall* const w : wallHandler)
1850  {
1851  w->setForce(Vec3D(0.0, 0.0, 0.0));
1852  w->setTorque(Vec3D(0.0, 0.0, 0.0));
1853  } //end reset forces loop
1854 
1855 #ifdef DEBUG_OUTPUT
1856  std::cerr << "Have all forces set to zero " << std::endl;
1857 #endif
1858 
1860 
1861  for (BaseParticle* p : particleHandler)
1862  {
1863 
1866  //end inner loop over contacts.
1867 
1869 
1870  }
1871 
1872 #ifdef CONTACT_LIST_HGRID
1873  //possibleContactList.write(std::cout);
1874  PossibleContact* Curr=possibleContactList.getFirstPossibleContact();
1875  while(Curr)
1876  {
1877  computeInternalForces(Curr->getP1(),Curr->getP2());
1878  Curr=Curr->getNext();
1879  }
1880 #endif
1881  //end outer loop over contacts.
1882 
1883 }
1884 
1889 {
1890  broadPhase(i);
1891 }
1892 
1897 void DPMBase::write(std::ostream& os, bool writeAllParticles) const
1898 {
1899  os << "restart_version " << "1.0";
1901  os << "xMin " << getXMin()
1902  << " xMax " << getXMax()
1903  << " yMin " << getYMin()
1904  << " yMax " << getYMax()
1905  << " zMin " << getZMin()
1906  << " zMax " << getZMax() << std::endl
1907  << "timeStep " << getTimeStep()
1908  << " time " << getTime()
1909  << " ntimeSteps " << ntimeSteps_
1910  << " timeMax " << getTimeMax() << std::endl
1911  << "systemDimensions " << getSystemDimensions()
1912  << " particleDimensions " << getParticleDimensions()
1913  << " gravity " << getGravity();
1914  os << " writeVTK " << writeParticlesVTK_ << " " << writeWallsVTK_ << " " << interactionHandler.getWriteVTK();
1915  os << " random "; random.write(os);
1916  //only write xBallsArguments if they are nonzero
1918  os << " xBallsArguments " << getXBallsAdditionalArguments();
1919  os << std::endl;
1920 
1921  speciesHandler.write(os);
1922  os << "Walls " << wallHandler.getNumberOfObjects() << std::endl;
1923  for (BaseWall* w : wallHandler)
1924  os << (*w) << std::endl;
1925  os << "Boundaries " << boundaryHandler.getNumberOfObjects() << std::endl;
1926  for (BaseBoundary* b : boundaryHandler)
1927  os << (*b) << std::endl;
1928  if (writeAllParticles || particleHandler.getNumberOfObjects() < 4)
1929  {
1930  particleHandler.write(os);
1931  } else
1932  {
1933  os << "Particles " << particleHandler.getNumberOfObjects() << std::endl;
1934  for (unsigned int i = 0; i < 2; i++)
1935  os << *particleHandler.getObject(i) << std::endl;
1936  os << "..." << std::endl;
1937  }
1938  if (writeAllParticles || interactionHandler.getNumberOfObjects() < 4)
1939  {
1941  } else
1942  {
1943  os << "Interactions " << interactionHandler.getNumberOfObjects() << std::endl;
1944  for (unsigned int i = 0; i < 2; i++)
1945  os << *interactionHandler.getObject(i) << std::endl;
1946  os << "..." << std::endl;
1947  }
1950 }
1951 
1955 void DPMBase::read(std::istream& is)
1956 {
1957  std::string dummy;
1958  is >> dummy;
1959  if (strcmp(dummy.c_str(), "restart_version"))
1960  {
1961  //only very old files did not have a restart_version
1962  logger.log(Log::FATAL, "Error in DPMBase::read(is): this is not a valid restart file");
1963  } else
1964  {
1965  is >> restartVersion_;
1966  if (!restartVersion_.compare("1.0"))
1967  {
1969  is >> dummy >> xMin_
1970  >> dummy >> xMax_
1971  >> dummy >> yMin_
1972  >> dummy >> yMax_
1973  >> dummy >> zMin_
1974  >> dummy >> zMax_
1975  >> dummy >> timeStep_
1976  >> dummy >> time_
1977  >> dummy >> ntimeSteps_
1978  >> dummy >> timeMax_
1979  >> dummy >> systemDimensions_
1980  >> dummy >> particleDimensions_
1981  >> dummy >> gravity_;
1982 
1983  std::stringstream line(std::stringstream::in | std::stringstream::out);
1985  line >> dummy;
1986  if (!dummy.compare("writeVTK")) {
1987  FileType writeInteractionsVTK = FileType::NO_FILE;
1988  line >> writeParticlesVTK_ >> writeWallsVTK_ >> writeInteractionsVTK >> dummy;
1989  interactionHandler.setWriteVTK(writeInteractionsVTK);
1990  }
1991  if (!dummy.compare("random"))
1992  {
1994  random.read(line);
1995  line >> dummy;
1997  setXBallsAdditionalArguments(line.str());
1998  }
1999 
2000  speciesHandler.read(is);
2001 
2002  unsigned int N;
2003  //std::stringstream line(std::stringstream::in | std::stringstream::out);
2004  is >> dummy >> N;
2005  if (dummy.compare("Walls"))
2006  logger(ERROR, "DPMBase::read(is): Error during restart: 'Walls' argument could not be found.");
2007  wallHandler.clear();
2010 
2011  for (unsigned int i = 0; i < N; i++)
2012  {
2014  wallHandler.readObject(line);
2015  }
2016 
2017  is >> dummy >> N;
2020  if (dummy.compare("Boundaries"))
2021  logger(ERROR, "DPMBase::read(is): Error during restart: 'Boundaries' argument could not be found.");
2023  for (unsigned int i = 0; i < N; i++)
2024  {
2027  }
2028 
2029  is >> dummy >> N;
2030  if (dummy.compare("Particles"))
2031  logger(ERROR, "DPMBase::read(is): Error during restart: 'Boundaries' argument could not be found.");
2035  for (unsigned int i = 0; i < N; i++)
2036  {
2040  //particleHandler.getLastObject()->computeMass();
2041  }
2042 
2044 
2045  } else if (!restartVersion_.compare("3"))
2046  {
2047  logger.log(Log::INFO, "DPMBase::read(is): restarting from an old restart file (restart_version %).",
2048  restartVersion_);
2049  readOld(is);
2050  } else
2051  {
2052  //only very old files did not have a restart_version
2053  logger.log(Log::FATAL,
2054  "Error in DPMBase::read(is): restart_version % cannot be read; use an older version of Mercury to upgrade the file",
2055  restartVersion_);
2056  }
2057  }
2058 }
2059 
2063 void DPMBase::readOld(std::istream& is)
2064 {
2065  std::string dummy;
2066  is >> dummy >> dummy;
2067  setName(dummy);
2068 
2069  unsigned int saveCountData, saveCountEne, saveCountStat, saveCountFStat;
2070  unsigned int fileTypeFstat, fileTypeData, fileTypeEne, fileTypeRestart;
2071  is >> dummy >> xMin_
2072  >> dummy >> xMax_
2073  >> dummy >> yMin_
2074  >> dummy >> yMax_
2075  >> dummy >> zMin_
2076  >> dummy >> zMax_
2077  >> dummy >> timeStep_
2078  >> dummy >> time_
2079  >> dummy >> timeMax_
2080  >> dummy >> saveCountData
2081  >> dummy >> saveCountEne
2082  >> dummy >> saveCountStat
2083  >> dummy >> saveCountFStat
2084  >> dummy >> systemDimensions_
2085  >> dummy >> gravity_
2086  >> dummy >> fileTypeFstat
2087  >> dummy >> fileTypeData
2088  >> dummy >> fileTypeEne;
2089  dataFile.setSaveCount(saveCountData);
2090  eneFile.setSaveCount(saveCountEne);
2091  statFile.setSaveCount(saveCountStat);
2092  fStatFile.setSaveCount(saveCountFStat);
2093 
2094  fStatFile.setFileType(static_cast<FileType>(fileTypeFstat));
2095  dataFile.setFileType(static_cast<FileType>(fileTypeData));
2096  eneFile.setFileType(static_cast<FileType>(fileTypeEne));
2097 
2098  //this is optional to allow restart files with and without restartFile.getFileType()
2099  is >> dummy;
2100  if (!strcmp(dummy.c_str(), "options_restart"))
2101  {
2102  is >> fileTypeRestart;
2103  restartFile.setFileType(static_cast<FileType>(fileTypeRestart));
2104  }
2105 
2106  speciesHandler.read(is);
2107  wallHandler.read(is);
2108  boundaryHandler.read(is);
2109  particleHandler.read(is);
2110 }
2111 
2116 {
2118  {
2119  std::cerr << "Error in solve(): at least one species has to be defined" << std::endl;
2120  exit(-1);
2121  }
2122  if (getParticleDimensions() == 0)
2123  {
2124  std::cerr << "Error in solve(): particleDimensions not specified" << std::endl;
2125  exit(-1);
2126  }
2127  if (getSystemDimensions() == 0)
2128  {
2129  std::cerr << "Error in solve(): systemDimensions not specified" << std::endl;
2130  exit(-1);
2131  }
2132  if (getTimeStep() == 0.0)
2133  {
2134  std::cerr << "Error in solve(): timeStep not specified" << std::endl;
2135  exit(-1);
2136  }
2137 }
2138 
2140 {
2143 
2145  {
2150  }
2151 
2153  {
2154  printTime();
2158  writeVTK();
2159  }
2160 
2161 // if (statFile.saveCurrentTimestep(ntimeSteps_))
2162 // {
2163 // outputStatistics();
2164 // gatherContactStatistics();
2165 // processStatistics(true);
2166 // }
2167 
2168  //write restart file last, otherwise the output cunters are wrong
2170  {
2171  writeRestartFile();
2172  restartFile.close(); //overwrite old restart file if FileType::ONE_FILE
2173  }
2174 }
2175 
2189 {
2190  logger(DEBUG, "Entered solve");
2191 #ifdef CONTACT_LIST_HGRID
2192  logger(INFO,"Using CONTACT_LIST_HGRID");
2193 #endif
2194 
2199  if (!getRestarted())
2200  {
2201  ntimeSteps_ = 0;
2202  resetFileCounter();
2203  setTime(0.0);
2204  setNextSavedTimeStep(0); //reset the counter
2205  //this is to ensure that the interaction time stamps agree with the resetting of the time value
2206  for (auto& i : interactionHandler)
2207  i->setTimeStamp(0);
2209  logger(DEBUG, "Have created the particles initial conditions");
2210  }
2211  else
2212  {
2213  actionsOnRestart();
2214  }
2215 
2216  checkSettings();
2217 
2218  if (getRunNumber() > 0 && !getRestarted())
2219  {
2220  std::stringstream name;
2221  name << getName() << "." << getRunNumber();
2222  setName(name.str());
2223  }
2224 
2225  //append_ determines if files have to be appended or not
2226  if (getAppend())
2227  {
2228  setOpenMode(std::fstream::out | std::fstream::app);
2229  restartFile.setOpenMode(std::fstream::out); //Restart files should always be overwritten.
2230  }
2231  else
2232  {
2233  setOpenMode(std::fstream::out);
2234  }
2235 
2237 
2238  // Setup the mass of each particle.
2240 
2241  // Other initialisations
2242  //max_radius = getLargestParticle()->getRadius();
2245 
2246  // do a first force computation
2249  time_ -= timeStep_; //to enforce that old time stamps are created
2250  ntimeSteps_--;
2251  computeAllForces();
2252  ntimeSteps_++;
2253  time_ += timeStep_;
2256  logger(DEBUG, "Have computed the initial values for the forces ");
2257 
2258  // This is the main loop over advancing time
2259  while (getTime() < getTimeMax() && continueSolve())
2260  {
2261  writeOutputFiles(); //everything is written at the beginning of the timestep!
2262  // Loop over all particles doing the time integration step
2265  checkInteractionWithBoundaries(); // INSERTION boundaries handled
2267 
2268  // Compute forces
2269 
2271  // INSERTION/DELETION boundary flag change
2272  for (BaseBoundary* b : boundaryHandler)
2273  {
2274  b->checkBoundaryBeforeTimeStep(this);
2275  }
2276 
2278 
2280 
2282 
2283  computeAllForces();
2284 
2286 
2288 
2289  // Loop over all particles doing the time integration step
2292 
2293  checkInteractionWithBoundaries(); // DELETION boundaries handled
2295 
2296  //erase interactions that have not been used during the last timestep
2299 
2300  time_ += timeStep_;
2301  ntimeSteps_++;
2302  }
2303  //force writing of the last time step
2305  writeOutputFiles();
2306 
2307  //end loop over interaction count
2309 
2310  //To make sure getTime gets the correct time for outputting statistics
2311  finishStatistics();
2312 
2313  closeFiles();
2314 }
2315 
2320 bool DPMBase::readArguments(int argc, char* argv[])
2321 {
2322  bool isRead = true;
2323  for (int i = 1; i < argc; i += 2)
2324  {
2325  std::cout << "interpreting input argument " << argv[i];
2326  for (int j = i + 1; j < argc; j++)
2327  {
2328  if (argv[j][0] == '-')
2329  break;
2330  std::cout << " " << argv[j];
2331  }
2332  std::cout << std::endl;
2333  isRead &= readNextArgument(i, argc, argv);
2334  if (!isRead)
2335  {
2336  std::cerr << "Warning: not all arguments read correctly!" << std::endl;
2337  exit(-1);
2338  }
2339  }
2340  return isRead;
2341 }
2342 
2349 bool DPMBase::readNextArgument(int& i, int argc, char* argv[])
2350 {
2352  if (!strcmp(argv[i], "-name"))
2353  {
2354  setName(argv[i + 1]);
2355  } else if (!strcmp(argv[i], "-xmin"))
2356  {
2357  setXMin(atof(argv[i + 1]));
2358  } else if (!strcmp(argv[i], "-ymin"))
2359  {
2360  setYMin(atof(argv[i + 1]));
2361  } else if (!strcmp(argv[i], "-zmin"))
2362  {
2363  setZMin(atof(argv[i + 1]));
2364  } else if (!strcmp(argv[i], "-xmax"))
2365  {
2366  setXMax(atof(argv[i + 1]));
2367  } else if (!strcmp(argv[i], "-ymax"))
2368  {
2369  setYMax(atof(argv[i + 1]));
2370  } else if (!strcmp(argv[i], "-zmax"))
2371  {
2372  setZMax(atof(argv[i + 1]));
2373  //} else if (!strcmp(argv[i],"-svn")) {
2374  // std::cout << "svn version " << SVN_VERSION << std::endl;
2375  // i--;
2376  } else if (!strcmp(argv[i], "-addS"))
2377  {
2378  setAddSmall(atof(argv[i + 1]));
2379  } else if (!strcmp(argv[i], "-addL"))
2380  {
2381  setAddLarge(atof(argv[i + 1]));
2382  } else if (!strcmp(argv[i], "-dt"))
2383  {
2384  Mdouble old = getTimeStep();
2385  setTimeStep(atof(argv[i + 1]));
2386  std::cout << " reset dt from " << old << " to " << getTimeStep() << std::endl;
2387  }
2388 // else if (!strcmp(argv[i], "-Hertz"))
2389 // {
2390 // speciesHandler.getObject(0)->setForceType(ForceType::HERTZ);
2391 // i--;
2392 // }
2393  else if (!strcmp(argv[i], "-tmax"))
2394  {
2395  Mdouble old = getTimeMax();
2396  setTimeMax(atof(argv[i + 1]));
2397  std::cout << " reset timeMax from " << old << " to " << getTimeMax() << std::endl;
2398  } else if (!strcmp(argv[i], "-saveCount"))
2399  {
2400  Mdouble old = dataFile.getSaveCount();
2401  setSaveCount(static_cast<unsigned int>(atoi(argv[i + 1])));
2403  std::cout << " reset saveCount from " << old << " to " << dataFile.getSaveCount() << std::endl;
2404  } else if (!strcmp(argv[i], "-saveCountData"))
2405  {
2406  dataFile.setSaveCount(static_cast<unsigned int>(atoi(argv[i + 1])));
2407  } else if (!strcmp(argv[i], "-saveCountFStat"))
2408  {
2409  fStatFile.setSaveCount(static_cast<unsigned int>(atoi(argv[i + 1])));
2410  } else if (!strcmp(argv[i], "-saveCountStat"))
2411  {
2412  statFile.setSaveCount(static_cast<unsigned int>(atoi(argv[i + 1])));
2413  } else if (!strcmp(argv[i], "-saveCountEne"))
2414  {
2415  eneFile.setSaveCount(static_cast<unsigned int>(atoi(argv[i + 1])));
2416  } else if (!strcmp(argv[i], "-saveCountRestart"))
2417  {
2418  restartFile.setSaveCount(static_cast<unsigned int>(atoi(argv[i + 1])));
2419  } else if (!strcmp(argv[i], "-dim"))
2420  {
2421  setSystemDimensions(static_cast<unsigned int>(atoi(argv[i + 1])));
2422  } else if (!strcmp(argv[i], "-gravity"))
2423  {
2425  setGravity(Vec3D(atof(argv[i + 1]), atof(argv[i + 2]), atof(argv[i + 3])));
2426  i += 2;
2427  } else if (!strcmp(argv[i], "-fileType"))
2428  { //uses int input
2429  setFileType(static_cast<FileType>(atoi(argv[i + 1])));
2430  } else if (!strcmp(argv[i], "-fileTypeFStat"))
2431  { //uses int input
2432  fStatFile.setFileType(static_cast<FileType>(atoi(argv[i + 1])));
2433  } else if (!strcmp(argv[i], "-fileTypeRestart"))
2434  {
2435  restartFile.setFileType(static_cast<FileType>(atoi(argv[i + 1])));
2436  } else if (!strcmp(argv[i], "-fileTypeData"))
2437  {
2438  dataFile.setFileType(static_cast<FileType>(atoi(argv[i + 1])));
2439  } else if (!strcmp(argv[i], "-fileTypeStat"))
2440  {
2441  statFile.setFileType(static_cast<FileType>(atoi(argv[i + 1])));
2442  } else if (!strcmp(argv[i], "-fileTypeEne"))
2443  {
2444  eneFile.setFileType(static_cast<FileType>(atoi(argv[i + 1])));
2445  } else if (!strcmp(argv[i], "-auto_number"))
2446  {
2447  autoNumber();
2448  i--;
2449  }
2450 // else if (!strcmp(argv[i], "-number_of_saves"))
2451 // {
2452 // set_number_of_saves_all(atof(argv[i + 1]));
2453 // }
2454  else if (!strcmp(argv[i], "-restart") || !strcmp(argv[i], "-r"))
2455  {
2459  std::string filename;
2460 
2461  //use default filename if no argument is given
2462  if (i + 1 >= argc || argv[i + 1][0] == '-')
2463  {
2464  i--;
2465  filename = getName();
2466  std::cout << getName() << std::endl;
2467  } else
2468  filename = argv[i + 1];
2469 
2470  //add ".restart" if necessary
2471  if (filename.find(".restart") == std::string::npos)
2472  {
2473  filename = filename + ".restart";
2474  }
2475 
2476  std::cout << "restart from " << filename << std::endl;
2477  readRestartFile(filename);
2478  } else if (!strcmp(argv[i], "-clean") || !strcmp(argv[i], "-c"))
2479  {
2480  std::cout << "Remove old " << getName() << ".* files" << std::endl;
2482 // std::string filename;
2483  std::ostringstream filename;
2484  std::vector<std::string> ext{".restart", ".stat", ".fstat", ".data", ".ene", ".xballs"};
2485  for (unsigned int j = 0; j < ext.size(); j++)
2486  {
2487  // remove files with given extension for FileType::ONE_FILE
2488  filename.clear();
2489  filename << getName() << ext[j];
2490  if (!remove(filename.str().c_str()))
2491  {
2492  std::cout << " File " << filename.str() << " successfully deleted" << std::endl;
2493  }
2494  // remove files with given extension for FileType::MULTIPLE_FILES
2495  unsigned int k = 0;
2496  filename.clear();
2497  filename << getName() << ext[j] << '.' << k;
2498  while (!remove(filename.str().c_str()))
2499  {
2500  std::cout << " File " << filename.str() << " successfully deleted" << std::endl;
2501  filename.clear();
2502  filename << getName() << ext[j] << '.' << ++k;
2503  }
2504  // remove files with given extension for FileType::MULTIPLE_FILES_PADDED
2505  k = 0;
2506  filename.clear();
2507  filename << getName() << ext[j] << '.' << to_string_padded(k);
2508  while (!remove(filename.str().c_str()))
2509  {
2510  std::cout << " File " << filename.str() << " successfully deleted" << std::endl;
2511  filename.clear();
2512  filename << getName() << ext[j] << '.' << to_string_padded(++k);
2513  }
2514  }
2515  i--;
2516  } else if (!strcmp(argv[i], "-data"))
2517  {
2518  std::string filename = argv[i + 1];
2519  readDataFile(filename.c_str());
2520  }
2521 // else if (!strcmp(argv[i], "-k"))
2522 // {
2523 // Mdouble old = getSpecies(0)->getStiffness();
2524 // getSpecies(0)->setStiffness(atof(argv[i + 1]));
2525 // std::cout << " reset k from " << old << " to " << getSpecies(0)->getStiffness() << std::endl;
2526 // }
2527 // else if (!strcmp(argv[i], "-dissipation") || !strcmp(argv[i], "-disp"))
2528 // {
2529 // Mdouble old = getSpecies(0)->getDissipation();
2530 // getSpecies(0)->setDissipation(atof(argv[i + 1]));
2531 // std::cout << " reset getDissipation() from " << old << " to " << getSpecies(0)->getDissipation() << std::endl;
2532 // }
2533 // else if (!strcmp(argv[i], "-kt"))
2534 // {
2535 // Mdouble old = getSpecies(0)->getSlidingStiffness();
2536 // getSpecies(0)->setSlidingStiffness(atof(argv[i + 1]));
2537 // std::cout << " reset kt from " << old << " to " << getSpecies(0)->getSlidingStiffness() << std::endl;
2538 // }
2539 // else if (!strcmp(argv[i], "-dispt"))
2540 // {
2541 // Mdouble old = getSpecies(0)->getSlidingDissipation();
2542 // getSpecies(0)->setSlidingDissipation(atof(argv[i + 1]));
2543 // std::cout << " reset dispt from " << old << " to " << getSpecies(0)->getSlidingDissipation() << std::endl;
2544 // }
2545 // else if (!strcmp(argv[i], "-krolling"))
2546 // {
2547 // Mdouble old = getSpecies(0)->getRollingStiffness();
2548 // getSpecies(0)->setRollingStiffness(atof(argv[i + 1]));
2549 // std::cout << " reset krolling from " << old << " to " << getSpecies(0)->getRollingStiffness() << std::endl;
2550 // }
2551 // else if (!strcmp(argv[i], "-disprolling"))
2552 // {
2553 // Mdouble old = getSpecies(0)->getRollingDissipation();
2554 // getSpecies(0)->setRollingDissipation(atof(argv[i + 1]));
2555 // std::cout << " reset disprolling from " << old << " to " << getSpecies(0)->getRollingDissipation() << std::endl;
2556 // }
2557 // else if (!strcmp(argv[i], "-mu"))
2558 // {
2559 // Mdouble old = getSpecies(0)->getSlidingFrictionCoefficient();
2560 // getSpecies(0)->setSlidingFrictionCoefficient(atof(argv[i + 1]));
2561 // std::cout << " reset mu from " << old << " to " << getSpecies(0)->getSlidingFrictionCoefficient() << std::endl;
2562 // }
2563 // else if (!strcmp(argv[i], "-murolling"))
2564 // {
2565 // Mdouble old = getSpecies(0)->getRollingFrictionCoefficient();
2566 // getSpecies(0)->setRollingFrictionCoefficient(atof(argv[i + 1]));
2567 // std::cout << " reset murolling from " << old << " to " << getSpecies(0)->getRollingFrictionCoefficient() << std::endl;
2568 // }
2569  else if (!strcmp(argv[i], "-randomise") || !strcmp(argv[i], "-randomize"))
2570  {
2571  random.randomise();
2572  i--;
2573  }
2574 // else if (!strcmp(argv[i], "-k0"))
2575 // {
2576 // Mdouble old = speciesHandler.getObject(0)->getAdhesionStiffness();
2577 // speciesHandler.getObject(0)->setAdhesionStiffness(atof(argv[i + 1]));
2578 // std::cout << " reset k0 from " << old << " to " << speciesHandler.getObject(0)->getAdhesionStiffness() << std::endl;
2579 // }
2580 // else if (!strcmp(argv[i], "-f0"))
2581 // {
2582 // Mdouble old = speciesHandler.getObject(0)->getBondForceMax();
2583 // speciesHandler.getObject(0)->setBondForceMax(atof(argv[i + 1]));
2584 // std::cout << " reset f0 from " << old << " to " << speciesHandler.getObject(0)->getBondForceMax() << std::endl;
2585 // }
2586 // else if (!strcmp(argv[i], "-AdhesionForceType"))
2587 // {
2588 // AdhesionForceType old = speciesHandler.getObject(0)->getAdhesionForceType();
2589 // speciesHandler.getObject(0)->setAdhesionForceType(argv[i + 1]);
2590 // std::cout << " reset AdhesionForceType from "
2591 // << static_cast<signed char>(old) << " to "
2592 // << static_cast<signed char>(speciesHandler.getObject(0)->getAdhesionForceType()) << std::endl;
2593 // }
2594  else if (!strcmp(argv[i], "-append"))
2595  {
2596  setAppend(true);
2597  i--;
2598  } else if (!strcmp(argv[i], "-fixedParticles"))
2599  {
2600  setFixedParticles(static_cast<unsigned int>(atoi(argv[i + 1])));
2601  }
2602 // else if (!strcmp(argv[i], "-rho"))
2603 // {
2604 // Mdouble old = speciesHandler.getObject(0)->getDensity();
2605 // speciesHandler.getObject(0)->setDensity(atof(argv[i + 1]));
2606 // std::cout << " reset rho from " << old << " to " << speciesHandler.getObject(0)->getDensity() << std::endl;
2607 // }
2608 // else if (!strcmp(argv[i], "-dim_particle"))
2609 // {
2610 // setParticleDimensions(atoi(argv[i + 1]));
2611 // }
2612  else if (!strcmp(argv[i], "-counter"))
2613  {
2614  setRunNumber(atoi(argv[i + 1]));
2615  } else
2616  return false;
2617  return true; //returns true if argv is found
2618 }
2619 
2633 {
2634 #ifdef MERCURY_USE_MPI
2635  int localInteraction = checkParticleForInteractionLocal(p);
2636 
2637  MPIContainer& communicator = MPIContainer::Instance();
2638  int numberOfProcessors = communicator.getNumberOfProcessors();
2639 
2640  //The root gathers all values and computes the global value
2641  int *interactionList = nullptr;
2642  if (communicator.getProcessorID() == 0)
2643  {
2644  interactionList = new int [numberOfProcessors];
2645  }
2646 
2647  //Gather all local values
2648  communicator.gather(localInteraction,interactionList);
2649 
2650  //Compute the global value
2651  int globalInteraction = 1;
2652  if (communicator.getProcessorID() == 0)
2653  {
2654  for (int i = 0; i<numberOfProcessors; i++)
2655  {
2656  if (interactionList[i] == 0)
2657  {
2658  globalInteraction = 0;
2659  break;
2660  }
2661  }
2662  }
2663  //The root now tells the other processors what the global value for the interaction is
2664  communicator.broadcast(globalInteraction);
2665 
2666  //Convert the result back to bool
2667  bool interaction = globalInteraction;
2668 #else
2669  bool interaction = checkParticleForInteractionLocalPeriodic(p);
2670 #endif
2671  return interaction;
2672 }
2673 
2681  //A vector of ghost particles of the particle that is to be inserted (empty if no periodic boundaries are present)
2682  std::vector<BaseParticle> pPeriodic;
2683  for (BaseBoundary* b : boundaryHandler) {
2684  PeriodicBoundary* pb = dynamic_cast<PeriodicBoundary*>(b);
2685  if (pb) {
2686  for (int i=pPeriodic.size()-1; i>=0; --i) {
2687  if (particleHandler.getNumberOfObjects()>0 && pb->getDistance(pPeriodic[i]) < pPeriodic[i].getInteractionRadius() + particleHandler.getLargestParticle()->getInteractionRadius()) {
2688  pPeriodic.push_back(pPeriodic[i]);
2689  pb->shiftPosition(&pPeriodic.back());
2690  }
2691  }
2693  pPeriodic.push_back(p);
2694  pb->shiftPosition(&pPeriodic.back());
2695  }
2696  }
2697  }
2698  //check the particle AND the ghost particles for intersection problems
2699  bool insertable = checkParticleForInteractionLocal(p);
2700  for (BaseParticle& pP : pPeriodic) {
2701  insertable &= checkParticleForInteractionLocal(pP);
2702  }
2703  return insertable;
2704 }
2705 
2719 {
2720  Mdouble distance;
2721  Vec3D normal;
2722 
2723  //Check if it has no collision with walls
2724  for (BaseWall* w : wallHandler)
2725  {
2726  //returns false if the function getDistanceAndNormal returns true,
2727  //i.e. if there exists an interaction between wall and particle
2728  //\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.)
2729  if (w->getDistanceAndNormal(p, distance, normal))
2730  {
2731  //std::cout<<"failure: Collision with wall: "<<**it<<std::endl;
2732  return false;
2733  }
2734  else
2735  {
2736  //std::cout<<"No collision with wall: "<<**it<<std::endl;
2737  }
2738  }
2739 
2740  //Check if it has no collision with other particles
2741  for (BaseParticle* q : particleHandler)
2742  {
2743  //returns false if the particle separation is less than the relevant sum of interaction radii
2744  //(i.e. another particle is in contact with P)
2745  if (Vec3D::getDistanceSquared(q->getPosition(), p.getPosition())
2746  < mathsFunc::square(q->getInteractionRadius() + p.getInteractionRadius()))
2747  {
2748  //std::cout<<"failure: Collision with particle "<<**it<<std::endl;
2749  return false;
2750  }
2751  else
2752  {
2753  //std::cout<<"No collision with particle "<<**it<<std::endl;
2754  }
2755  }
2756  return true;
2758 }
2759 
2766 {
2767  for (unsigned int i = particleHandler.getNumberOfObjects();
2768  i >= 1 && particleHandler.getObject(i - 1)->getPeriodicFromParticle() != nullptr; i--)
2769  {
2770  while (particleHandler.getObject(i - 1)->getInteractions().size() > 0)
2771  {
2773  particleHandler.getObject(i - 1)->getInteractions().front()->getIndex());
2774  }
2776  }
2777 }
2778 
2783 {
2784  for (BaseBoundary* b : boundaryHandler) // For all pointers of type BaseBoundary in Boundary handler
2785  {
2786  unsigned int N = particleHandler.getNumberOfObjects(); //Required because the number of particles increases
2787  for (unsigned int i = 0; i < N; i++)
2788  {
2789  b->createPeriodicParticles(particleHandler.getObject(i), particleHandler);
2790  }
2791  }
2792 }
2793 
2799 {
2800  std::cout << "Interactions currently in the handler:" << std::endl;
2802  {
2803  p->write(std::cout);
2804  std::cout << std::endl;
2805  std::cout << "Interaction " << p->getName() << " " << p->getId() << " between " << p->getP()->getId() << " and "
2806  << p->getI()->getId() << std::endl;
2807  }
2808 }
2809 
2814 {
2815  return getTime() <= time && getTime() + getTimeStep() > time;
2816 }
2817 
2819 {
2820  return (getXMax() - getXMin()) * (getYMax() - getYMin()) * (getZMax() - getZMin());
2821 }
2822 
void writeVTK() const
Writes all particles into a vtk file format (unstructured grid), consisting of particle positions...
Mdouble timeMax_
Stores the duration of the simulation.
Definition: DPMBase.h:914
void setTime(Mdouble time)
Access function for the time.
Definition: DPMBase.cc:185
BaseParticle * getLargestParticle() const
Gets a pointer to the largest BaseParticle (by interactionRadius) in this ParticleHandler.
virtual void actionsBeforeTimeLoop()
A virtual function. Allows one to carry out any operations before the start of the time loop...
Definition: DPMBase.cc:722
unsigned int getId() const
Returns the unique identifier of any particular object.
Definition: BaseObject.cc:116
void setXBallsVectorScale(double newVScale)
Set the scale of vectors in xballs.
Definition: DPMBase.cc:489
virtual void write(std::ostream &os, bool writeAllParticles=true) const
Loads all MD data and plots statistics for all timesteps in the .data file.
Definition: DPMBase.cc:1897
FileType getWallsWriteVTK() const
Returns a flag indicating if particle rotation is enabled or disabled.
Definition: DPMBase.cc:263
virtual void hGridUpdateParticle(BaseParticle *obj UNUSED)
no implementation but can be overidden in its derived classes.
Definition: DPMBase.cc:757
void set(Vec3D normal, Mdouble distanceLeft, Mdouble distanceRight)
Defines a periodic wall.
void setXMax(Mdouble newXMax)
If the length of the problem domain in x-direction is XMax - XMin, this method sets XMax...
Definition: DPMBase.cc:415
bool writeToFile(std::string filename, std::string filecontent)
Writes a string to a file.
Definition: Helpers.cc:418
void setWallsWriteVTK(FileType writeWallsVTK)
Allows to set the flag for enabling or disabling particle rotation in the simulations.
Definition: DPMBase.cc:247
void write(std::ostream &os) const
void solve()
The work horse of the code.
Definition: DPMBase.cc:2188
Mdouble yMax_
If the length of the problem domain in y-direction is YMax - XMin, the above variable stores YMax...
Definition: DPMBase.h:866
virtual void hGridInsertParticle(BaseParticle *obj UNUSED)
no implementation but can be overidden in its derived classes.
Definition: DPMBase.cc:750
int addL_
Stores the number of large particles that are to be added on restart.
Definition: DPMBase.h:886
The DPMBase header includes quite a few header files, defining all the handlers, which are essential...
Definition: DPMBase.h:65
Mdouble X
the vector components
Definition: Vector.h:52
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, ..
void writeVTK() const
Writes all particles into a vtk file format (unstructured grid), consisting of particle positions...
virtual void integrateAfterForceComputation()
Integration is done after force computations. We apply the Velocity verlet scheme. See http://en.wikipedia.org/wiki/Verlet_integration#Velocity_Verlet.
Definition: DPMBase.cc:1779
void checkAndDuplicatePeriodicParticles()
In case of periodic boundaries, the below methods checks and adds particles when necessary into the p...
Definition: DPMBase.cc:2782
void read(std::istream &is)
Definition: RNG.cc:52
int getXBallsColourMode() const
Get the xball colour mode (CMode)
Definition: DPMBase.cc:481
void setNextSavedTimeStep(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: Files.cc:263
void setTimeMax(Mdouble newTMax)
Allows the upper time limit to be changed.
Definition: DPMBase.cc:199
bool areInContact(const BaseParticle *pI, const BaseParticle *pJ) const
Checks if two particle are in contact or is there any positive overlap.
Definition: DPMBase.cc:713
void constructor()
A function which initialises the member variables to default values, so that the problem can be solve...
Definition: DPMBase.cc:905
void write(std::ostream &os) const
Writes the InteractionHandler to an output stream, for example a restart file.
unsigned int getSystemDimensions() const
Returns the dimension of the simulation. Note there is also a particle dimension. ...
Definition: DPMBase.cc:575
bool readArguments(int argc, char *argv[])
Can interpret main function input arguments that are passed by the driver codes.
Definition: DPMBase.cc:2320
FileType getWriteVTK() const
FileType getFileType() const
Gets the file type e.g. NOFILE, ONEFILE and MULTIPLE FILES. File::fileType_.
Definition: File.cc:203
BaseParticle * getP1()
Gets a pointer to the first BaseParticle in this PossibleContact.
int getAddLarge() const
returns the number of large particles that are to be added on restart.
Definition: DPMBase.cc:327
void setAddLarge(int new_addL)
sets the number of large particles that are to be added on restart.
Definition: DPMBase.cc:351
virtual void removeObject(unsigned const int index)
Removes a BaseParticle from the ParticleHandler.
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
unsigned int particleDimensions_
determines if 2D or 3D particle volume is used for mass calculations
Definition: DPMBase.h:837
bool readRestartFile(bool restarted=true)
Reads all the particle data corresponding to the current saved time step. Which is what the restart f...
Definition: DPMBase.cc:1629
void addForce(const Vec3D &addForce)
Adds an amount to the force on this BaseInteractable.
Mdouble yMin_
If the length of the problem domain in y-direction is YMax - YMin, the above variable stores YMin...
Definition: DPMBase.h:860
void setDPMBase(DPMBase *DPMBase)
Sets the problem that is solved using this handler.
Definition: BaseHandler.h:536
void setYMin(Mdouble newYMin)
If the length of the problem domain in y-direction is YMax - YMin, this method sets YMin...
Definition: DPMBase.cc:371
virtual void computeInternalForces(BaseParticle *i)
Computes the forces between particles (internal in the sense that the sum over all these forces is ze...
Definition: DPMBase.cc:1888
Mdouble getZMax() const
If the length of the problem domain in z-direction is ZMax - ZMin, then getZMax() returns ZMax...
Definition: DPMBase.cc:319
Mdouble getXMin() const
If the length of the problem domain in x-direction is XMax - XMin, then getXMin() returns XMin...
Definition: DPMBase.cc:279
const ParticleSpecies * getSpecies() const
Returns a pointer to the species of this BaseInteractable.
void setStorageCapacity(const unsigned int N)
Sets the storage capacity of this BaseHandler.
Definition: BaseHandler.h:501
void setAddSmall(int new_addS)
sets the number of small particles that are to be added on restart.
Definition: DPMBase.cc:343
unsigned int getParticleDimensions() const
Returns the particle dimensions.
Definition: DPMBase.cc:599
bool readDataFile(const std::string fileName, unsigned int format=0)
This allows particle data to be reloaded from data files.
Definition: DPMBase.cc:1202
void setParticlesWriteVTK(bool writeParticlesVTK)
Allows to set the flag for enabling or disabling particle rotation in the simulations.
Definition: DPMBase.cc:255
void computeAllMasses(unsigned int indSpecies)
Computes the mass for all BaseParticle of the given species in this ParticleHandler.
virtual void initialiseStatistics()
no implementation but can be overidden in its derived classes.
Definition: DPMBase.cc:800
File restartFile
An instance of class File to handle in- and output into a .restart file.
Definition: Files.h:219
Mdouble elasticEnergy_
used in force calculations
Definition: DPMBase.h:921
Mdouble getYMin() const
If the length of the problem domain in y-direction is YMax - YMin, then getYMin() returns YMin...
Definition: DPMBase.cc:295
void setXBallsColourMode(int newCMode)
Set the xball output mode.
Definition: DPMBase.cc:473
const std::vector< T * >::const_iterator end() const
Gets the end of the const_iterator over all BaseBoundary in this BaseHandler.
Definition: BaseHandler.h:522
double Mdouble
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:2063
void closeFiles()
Closes all files (ene, data, fstat, restart, stat) that were opened to read or write.
Definition: Files.cc:252
void setParticleDimensions(unsigned int particleDimensions)
Allows the dimension of the particle (f.e. for mass) to be changed. e.g. discs or spheres...
Definition: DPMBase.cc:583
virtual ~DPMBase()
virtual destructor
Definition: DPMBase.cc:151
virtual void actionsAfterTimeStep()
A virtual function which allows to define operations to be executed after time step.
Definition: DPMBase.cc:793
void readObject(std::istream &is)
Reads BaseParticle into the ParticleHandler from restart data.
void unfix()
Unfix Particle function, which required a reference to the Species vector. It unfixes a Particle by c...
void setZMax(Mdouble newZMax)
If the length of the problem domain in z-direction is XMax - XMin, this method sets ZMax...
Definition: DPMBase.cc:439
void write(std::ostream &os) const
Accepts an output stream read function, which accepts an input stream std::ostream.
bool openNextFile()
This function should be called before a data corresponding to the new time step is written or read...
Definition: File.cc:231
void setSpecies(const ParticleSpecies *species)
void setDimension(unsigned int newDim)
Sets the system and particle dimension.
Definition: DPMBase.cc:553
static Mdouble dot(const Vec3D &a, const Vec3D &b)
Calculates the dot product of two Vec3D: .
Definition: Vector.cc:167
BaseParticle * getP2()
Gets a pointer to the second BaseParticle in this PossibleContact.
void setAppend(bool newAppendFlag)
Allows to set the append option.
Definition: DPMBase.cc:651
virtual void writeOutputFiles()
Writes the simulation data onto all the files i.e. .data, .ene, .fstat ...
Definition: DPMBase.cc:2139
void setSystemDimensions(unsigned int newDim)
Allows for the dimension of the simulation to be changed.
Definition: DPMBase.cc:562
void writeVTK() const
Definition: DPMBase.cc:1080
Species< LinearViscoelasticNormalSpecies, SlidingFrictionSpecies > LinearViscoelasticSlidingFrictionSpecies
std::function< void(std::string, std::string)> onFatal
Definition: Logger.h:185
void setFileType(FileType fileType)
Sets File::fileType_ for all files (ene, data, fstat, restart, stat)
Definition: Files.cc:170
void read(std::istream &is)
Reads all objects from restart data.
Definition: BaseHandler.h:407
virtual void write(std::ostream &os) const
Write all the species and mixed species to an output stream.
void setGravity(Vec3D newGravity)
Allows to modify the gravity vector.
Definition: DPMBase.cc:537
T square(T val)
squares a number
Definition: ExtendedMath.h:91
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:1147
void logWriteAndDie(std::string module, std::string message)
todo strcmp relies on this, should be changed to more modern version
Definition: DPMBase.cc:65
virtual bool getHGridUpdateEachTimeStep() const
Definition: DPMBase.cc:771
bool rotation_
A flag to turn on/off particle rotation.
Definition: DPMBase.h:944
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
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:91
void eraseOldInteractions(Mdouble lastTimeStep)
erases interactions which have an old timestamp.
void setRunNumber(int runNumber)
This sets the counter/Run number, overriding the defaults.
bool getParticlesWriteVTK() const
Returns a flag indicating if particle rotation is enabled or disabled.
Definition: DPMBase.cc:271
Mdouble zMax_
If the length of the problem domain in z-direction is ZMax - ZMin, the above variable stores ZMax...
Definition: DPMBase.h:878
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:35
bool restarted_
A bool to check if the simulation was restarted or not, ie. if setupInitialConditionsShould be run an...
Definition: DPMBase.h:933
const std::vector< T * >::const_iterator begin() const
Gets the begin of the const_iterator over all Object in this BaseHandler.
Definition: BaseHandler.h:508
Defines a pair of periodic walls. Inherits from BaseBoundary.
BaseParticle * getPeriodicFromParticle() const
Returns the 'original' particle this one's a periodic copy of.
Mdouble getXMax() const
If the length of the problem domain in x-direction is XMax - XMin, then getXMax() returns XMax...
Definition: DPMBase.cc:287
void clear()
Empties the whole ParticleHandler by removing all BaseParticle.
void setRotation(bool newRotFlag)
Allows to set the flag for enabling or disabling particle rotation in the simulations.
Definition: DPMBase.cc:231
void readObject(std::istream &is)
Reads BaseWall into the WallHandler from restart data.
Definition: WallHandler.cc:109
unsigned int getSaveCount() const
Gets File::saveCount_.
Definition: File.cc:275
PossibleContact * getNext()
Gets the next PossibleContact in the general linked list of PossibleContact.
void setOpenMode(std::fstream::openmode openMode)
Sets File::openMode_ for all files (ene, data, fstat, restart, stat)
Definition: Files.cc:201
virtual void processStatistics(bool usethese UNUSED)
no implementation but can be overidden in its derived classes.
Definition: DPMBase.cc:831
void close()
Closes the file by calling fstream_.close()
Definition: File.cc:362
bool findNextExistingDataFile(Mdouble tMin, bool verbose=true)
Useful when fileType is chosen as Multiple Files or Multiple files with padded.
Definition: DPMBase.cc:1401
void setYMax(Mdouble newYMax)
If the length of the problem domain in y-direction is YMax - YMin, this method sets YMax...
Definition: DPMBase.cc:427
unsigned int systemDimensions_
The dimensions of the simulation i.e. 2D or 3D.
Definition: DPMBase.h:832
void readObject(std::istream &is)
Reads BaseBoundary into the BoundaryHandler from restart data.
File eneFile
An instance of class File to handle in- and output into a .ene file.
Definition: Files.h:214
Mdouble zMin_
If the length of the problem domain in z-direction is ZMax - ZMin, the above variable stores ZMin...
Definition: DPMBase.h:872
Mdouble getElasticEnergy() const
Returns the global elastic energy within the system.
Definition: DPMBase.cc:659
Mdouble xMax_
If the length of the problem domain in x-direction is XMax - XMin, the above variable stores XMax...
Definition: DPMBase.h:854
Mdouble getGravitationalEnergy() const
Returns the global gravitational potential energy stored in the system.
Definition: DPMBase.cc:686
unsigned int getCounter() const
In case of multiple files, File::getCounter() returns the the number (FILE::Counter_) of the next fil...
Definition: File.cc:217
void removeObjectKeepingPeriodics(unsigned const int id)
Removes interactions of periodic particles when the periodic particles get deleted (see DPMBase::remo...
file will not be created/read
unsigned int ntimeSteps_
Stores the number of time steps.
Definition: DPMBase.h:904
virtual void writeEneTimestep(std::ostream &os) const
This function enables one to write the global energy available in the system after each time step...
Definition: DPMBase.cc:1044
std::fstream & getFstream()
Allows to access the member variable File::fstream_.
Definition: File.cc:151
virtual void hGridActionsBeforeIntegration()
no implementation but can be overidden in its derived classes.
Definition: DPMBase.cc:852
BaseParticle * getSmallestParticle() const
Gets a pointer to the smallest BaseParticle (by interactionRadius) in this ParticleHandler.
void write(std::ostream &os) const
Definition: RNG.cc:71
bool getAppend() const
Returns the flag denoting if the append option is on or off.
Definition: DPMBase.cc:643
virtual bool readNextArgument(int &i, int argc, char *argv[])
Interprets the i^th command-line argument.
Definition: DPMBase.cc:2349
virtual bool checkParticleForInteraction(const BaseParticle &P)
Checks if the particle having any interaction with walls or other particles.
Definition: DPMBase.cc:2632
void setName(const std::string &name)
Allows to set the name of all the files (ene, data, fstat, restart, stat)
Definition: Files.cc:150
Mdouble getDistance(const BaseParticle &p) const
Returns the distance of the wall to the particle.
void resetFileCounter()
Resets the file counter for each file i.e. for ene, data, fstat, restart, stat)
Definition: Files.cc:190
bool getRestarted() const
Returns the flag denoting if the simulation was restarted or not.
Definition: DPMBase.cc:626
Mdouble getMass() const
Returns the particle's mass_.
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:396
virtual void addObject(ParticleSpecies *const S)
Adds a new ParticleSpecies to the SpeciesHandler.
int xBallsColourMode_
XBalls is a package to view the particle data. As an alternative MercuryDPM also supports Paraview...
Definition: DPMBase.h:962
BoundaryHandler boundaryHandler
An object of the class BoundaryHandler which concerns insertion and deletion of particles into or fro...
Definition: DPMBase.h:1011
Mdouble xBallsScale_
sets the xballs argument scale (see XBalls/xballs.txt)
Definition: DPMBase.h:972
void setXMin(Mdouble newXMin)
If the length of the problem domain in x-direction is XMax - XMin, this method sets XMin...
Definition: DPMBase.cc:359
ParticleHandler particleHandler
An object of the class ParticleHandler, contains the pointers to all the particles created...
Definition: DPMBase.h:1001
int addS_
Stores the number of small particles that are to be added on restart.
Definition: DPMBase.h:892
void removeLastObject()
Removes the last BaseParticle from the ParticleHandler.
all data will be written into/ read from a single file called name_
Mdouble domainSize() const
Definition: DPMBase.cc:2818
T * getObject(const unsigned int id)
Gets a pointer to the Object at the specified index in the BaseHandler.
Definition: BaseHandler.h:451
unsigned int getNtimeSteps() const
Returns the current counter of time steps.
Definition: DPMBase.cc:177
void shiftPosition(BaseParticle *p) const
shifts the particle
Mdouble getKineticEnergy() const
Returns the global kinetic energy stored in the system.
Definition: DPMBase.cc:672
virtual void outputXBallsDataParticle(const unsigned int i, const unsigned int format, std::ostream &os) const
This function writes out the particle locations into an output stream in a format the XBalls program ...
Mdouble time_
Stores the current simulation time.
Definition: DPMBase.h:899
void setRestartVersion(std::string newRV)
Sets restart_version.
Definition: DPMBase.cc:617
Vec3D getGravity() const
Returns the gravity vector.
Definition: DPMBase.cc:545
void checkSettings()
Checks if the essentials are set properly to go ahead with solving the problem.
Definition: DPMBase.cc:2115
#define UNUSED
Definition: GeneralDefine.h:39
virtual void hGridActionsBeforeTimeLoop()
A virtual function that allows one to carry out hGrid operations before the start of the time loop...
Definition: DPMBase.cc:729
File fStatFile
An instance of class File to handle in- and output into a .fstat file.
Definition: Files.h:209
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:295
void read(std::istream &is)
Accepts an input stream std::istream.
bool readParAndIniFiles(const std::string fileName)
Allows the user to read par.ini files (useful to read MDCLR files)
Definition: DPMBase.cc:1231
void setSaveCount(unsigned int saveCount)
Sets File::saveCount_ for all files (ene, data, fstat, restart, stat)
Definition: Files.cc:139
static Vec3D cross(const Vec3D &a, const Vec3D &b)
Calculates the cross product of two Vec3D: .
Definition: Vector.cc:255
virtual void computeAllForces()
Computes all the forces acting on the particles by using the setTorque and setForce methods...
Definition: DPMBase.cc:1841
File dataFile
An instance of class File to handle in- and output into a .data file.
Definition: Files.h:204
Mdouble getRadius() const
Returns the particle's radius_.
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:1434
Basic class for walls.
Definition: BaseWall.h:44
void setZMin(Mdouble newZMin)
If the length of the problem domain in z-direction is ZMax - ZMin, this method sets ZMin...
Definition: DPMBase.cc:383
virtual void actionsAfterSolve()
A virtual function which allows to define operations to be executed after the solve().
Definition: DPMBase.cc:786
void setSaveCount(unsigned int saveCount)
Sets File::saveCount_.
Definition: File.cc:283
virtual void printTime() const
Displays the current simulation time onto your screen for example.
Definition: DPMBase.cc:886
unsigned int getNumberOfObjects() const
Gets the number of Object in this BaseHandler.
Definition: BaseHandler.h:487
SpeciesHandler speciesHandler
A handler to that stores the species type i.e. elastic, linear visco-elastic... et cetera...
Definition: DPMBase.h:991
InteractionHandler interactionHandler
An object of the class InteractionHandler.
Definition: DPMBase.h:1016
FileType writeWallsVTK_
A flag to turn on/off the vtk writer for walls.
Definition: DPMBase.h:949
virtual void writeXBallsScript() const
This writes a script which can be used to load the xballs problem to display the data just generated...
void outputInteractionDetails() const
Displays the interaction details corresponding to the pointer objects in the interaction handler...
Definition: DPMBase.cc:2798
void setXBallsAdditionalArguments(std::string newXBArgs)
Set the additional arguments for xballs.
Definition: DPMBase.cc:505
void setFixedParticles(unsigned int n)
Definition: DPMBase.cc:877
Mdouble getYMax() const
If the length of the problem domain in y-direction is YMax - YMin, then getYMax() returns XMax...
Definition: DPMBase.cc:303
Mdouble Y
Definition: Vector.h:52
Mdouble xBallsVectorScale_
sets the xballs argument vscale (see XBalls/xballs.txt)
Definition: DPMBase.h:967
Mdouble xMin_
If the length of the problem domain in x-direction is XMax - XMin, the above variable stores XMin...
Definition: DPMBase.h:848
virtual void actionsBeforeTimeStep()
A virtual function which allows to define operations to be executed before the new time step...
Definition: DPMBase.cc:779
LoggerOutput * loggerOutput
Declaration of the output functions. If the output needs to be redirected, please swap the loggerOutp...
Definition: Logger.cc:143
Mdouble timeStep_
Stores the simulation time step.
Definition: DPMBase.h:909
void setFileType(FileType fileType)
Sets the type of file needed to write into or read from. File::fileType_.
Definition: File.cc:210
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 ...
Definition: RNG.cc:46
WallHandler wallHandler
An object of the class WallHandler. Contains pointers to all the walls created.
Definition: DPMBase.h:1006
void setWriteVTK(FileType f)
RNG random
This is a random generator, often used for setting up the initial conditions etc...
Definition: DPMBase.h:996
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:61
virtual void hGridRemoveParticle(BaseParticle *obj UNUSED)
no implementation but can be overidden in its derived classes.
Definition: DPMBase.cc:764
virtual double getInfo(const BaseParticle &P) const
A virtual method that allows the user to overrride and set what is written into the info column in th...
Definition: DPMBase.cc:703
bool writeParticlesVTK_
A flag to turn on/off the vtk writer for particles.
Definition: DPMBase.h:954
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:736
void autoNumber()
The autoNumber() function is the trigger. It calls three functions. setRunNumber(), readRunNumberFromFile() and incrementRunNumberInFile().
Manages the linked list of PossibleContact.
virtual void writeEneHeader(std::ostream &os) const
Writes a header with a certain format for ENE file.
Definition: DPMBase.cc:973
bool open()
Checks if the file stream fstream_ has any issues while opening. Alongside, it also increments the ne...
Definition: File.cc:320
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...
Definition: InfiniteWall.cc:79
Mdouble getZMin() const
If the length of the problem domain in z-direction is ZMax - ZMin, then getZMin() returns ZMin...
Definition: DPMBase.cc:311
int getAddSmall() const
returns the number of small particles that are to be added on restart.
Definition: DPMBase.cc:335
bool saveCurrentTimestep(unsigned int ntimeSteps)
Definition: File.cc:309
void gatherContactStatistics()
Definition: DPMBase.cc:811
virtual void outputStatistics()
no implementation but can be overidden in its derived classes.
Definition: DPMBase.cc:807
void setTimeStep(Mdouble newDt)
Allows the time step dt to be changed.
Definition: DPMBase.cc:451
bool compare(std::istream &is, std::string s)
Definition: Helpers.cc:754
void removeDuplicatePeriodicParticles()
Removes periodic duplicate Particles.
Definition: DPMBase.cc:2765
bool checkParticleForInteractionLocalPeriodic(const BaseParticle &P)
Definition: DPMBase.cc:2680
void setMin(Vec3D min)
Sets the values xMin, yMin, zMin of the problem domain, which is [xMin,xMax]x[yMin,yMax]x[zMin,zMax].
Definition: DPMBase.cc:405
virtual void integrateBeforeForceComputation()
This is were the integration is done, at the moment it is velocity Verlet integration and is done bef...
Definition: DPMBase.cc:1748
virtual void writeRestartFile()
Stores all the particle data for current save time step. Calls the write function.
Definition: DPMBase.cc:1618
std::string xBallsAdditionalArguments_
A string of additional arguments for xballs can be specified (see XBalls/xballs.txt). e.g. "-solidf -v0".
Definition: DPMBase.h:977
double getXBallsVectorScale() const
Returns the scale of vectors used in xballs.
Definition: DPMBase.cc:497
This is a class defining walls.
Definition: InfiniteWall.h:47
virtual bool continueSolve() const
Definition: DPMBase.cc:897
each time-step will be written into/read from separate files numbered consecutively: name_...
It is publicly inherited from class Files. It defines an awesome feature that is ideal when doing a p...
static Mdouble getDistanceSquared(const Vec3D &a, const Vec3D &b)
Calculates the squared distance between two Vec3D: .
Definition: Vector.cc:280
void addTorque(const Vec3D &addTorque)
Adds an amount to the torque on this BaseInteractable.
std::ostream & operator<<(std::ostream &os, DPMBase &md)
Definition: DPMBase.cc:78
virtual void hGridUpdateMove(BaseParticle *, Mdouble)
no implementation but can be overidden in its derived classes.
Definition: DPMBase.cc:845
virtual void computeForcesDueToWalls(BaseParticle *PI)
Computes the forces on the particles due to the walls (normals are outward normals) ...
Definition: DPMBase.cc:1717
Class that describes a possible contact between two BaseParticle.
std::string restartVersion_
Previous versions of MercuryDPM had a different restart file format, the below member variable allows...
Definition: DPMBase.h:928
const std::list< BaseInteraction * > & getInteractions() const
Returns a reference to the list of interactions in this BaseInteractable.
T * getLastObject()
Gets a pointer to the last Object in this BaseHandler.
Definition: BaseHandler.h:473
void setRestarted(bool newRestartedFlag)
Allows to set the flag stating if the simulation is to be restarted or not.
Definition: DPMBase.cc:634
Implementation of a 3D vector (by Vitaliy).
Definition: Vector.h:45
void setName(const std::string &name)
Sets the file name, e.g. "Name.data".
Definition: File.cc:196
std::string getXBallsAdditionalArguments() const
Returns the additional arguments for xballs.
Definition: DPMBase.cc:513
virtual void hGridActionsAfterIntegration()
no implementation but can be overidden in its derived classes.
Definition: DPMBase.cc:859
Mdouble getTimeStep() const
Allows the time step dt to be accessed.
Definition: DPMBase.cc:465
File statFile
An instance of class File to handle in- and output into a .stat file.
Definition: Files.h:224
double getXBallsScale() const
Returns the scale of the view in xballs.
Definition: DPMBase.cc:529
bool isFixed() const
Is fixed Particle function. It returns whether a Particle is fixed or not, by checking its inverse Ma...
virtual void hGridActionsBeforeTimeStep()
A virtual function that allows one to set or execute hGrid parameters or operations before every simu...
Definition: DPMBase.cc:743
void setXBallsScale(Mdouble newScale)
Sets the scale of the view (either normal, zoom in or zoom out) to display in xballs. The default is fit to screen.
Definition: DPMBase.cc:521
Mdouble Z
Definition: Vector.h:52
Mdouble getInteractionRadius() const
Returns the particle's interaction radius, which might be different from radius_ (e.g., when dealing with wet particles)
bool getRotation() const
Returns a flag indicating if particle rotation is enabled or disabled.
Definition: DPMBase.cc:239
virtual void finishStatistics()
no implementation but can be overidden in its derived classes.
Definition: DPMBase.cc:838
virtual void broadPhase(BaseParticle *i)
By broad one means to screen and determine an approximate number of particles that a given particle c...
Definition: DPMBase.cc:866
Mdouble getTime() const
Access function for the time.
Definition: DPMBase.cc:169
const std::string & getName() const
Returns the name of the file. Does not allow to change it though.
Definition: Files.cc:132
virtual void checkInteractionWithBoundaries()
There are a range of boundaries one could implement depending on ones' problem. This methods checks f...
Definition: DPMBase.cc:1763
virtual void writeFstatHeader(std::ostream &os) const
Writes a header with a certain format for FStat file.
Definition: DPMBase.cc:994
virtual void read(std::istream &is)
Reads all data from a restart file, e.g. domain data and particle data.
Definition: DPMBase.cc:1955
Mdouble getTimeMax() const
Allows the user to access the total simulation time during the simulation. Cannot change it though...
Definition: DPMBase.cc:214
int getRunNumber() const
This returns the current value of the counter (runNumber_)
const std::string & getName() const
Allows to access the file name, e.g., "problem.data".
Definition: File.cc:163
DPMBase()
Constructor that calls the "void constructor()".
Definition: DPMBase.cc:143
void writeVTK() const
Writes all walls into a vtk format, consisting of points (edges) and cells (faces).
Definition: WallHandler.cc:239
void setOpenMode(std::fstream::openmode openMode)
Allows the user to Sets File::openMode_.
Definition: File.cc:268
Vec3D gravity_
Gravity vector.
Definition: DPMBase.h:842
void clear()
Empties the whole BaseHandler by removing all Objects and setting all other variables to 0...
Definition: BaseHandler.h:392
virtual bool checkParticleForInteractionLocal(const BaseParticle &P)
Checks if a particle P has any interaction with walls or other particles in the local domain...
Definition: DPMBase.cc:2718
virtual void setupInitialConditions()
This function allows to set the initial conditions for our problem to be solved, by default particle ...
Definition: DPMBase.cc:966
void setSpecies(const ParticleSpecies *species)
Define the species of this wall.
Definition: BaseWall.cc:113
void fixParticle()
Fix Particle function. It fixes a Particle by setting its inverse mass and inertia and velocities to ...
bool append_
A flag to determine if the file has to be appended or not. See DPMBase::Solve() for example...
Definition: DPMBase.h:939
void setMax(Vec3D max)
Sets the values xMax, yMax, zMax of the problem domain, which is [xMin,xMax]x[yMin,yMax]x[zMin,zMax].
Definition: DPMBase.cc:395
virtual void computeExternalForces(BaseParticle *PI)
Computes the external forces acting on particles (e.g. gravitational)
Definition: DPMBase.cc:1703
std::vector< BaseInteraction * > getInteractionWith(BaseParticle *P, Mdouble timeStamp, InteractionHandler *interactionHandler)
Checks if particle is in interaction with given particle P, and if so, returns pointer to the associa...
std::string getRestartVersion() const
This is to take into account for different Mercury versions. Returns the version of the restart file...
Definition: DPMBase.cc:608
bool isTimeEqualTo(Mdouble time) const
Checks if the input variable "time" is the current time in the simulation.
Definition: DPMBase.cc:2813