41 #include <sys/types.h>
45 #define GetCurrentDir _getcwd
51 #define GetCurrentDir getcwd
56 std::transform(s.begin(), s.end(), s.begin(), ::tolower);
334 ans.
disp = -2.0 * mass *
log(r) / tc;
342 ans.
dispt = 2. * sqrt(1.0 / 7.0 * mass * ans.
kt);
355 Mdouble DispCorrection =
exp(-disp / mass / w) * asin(w / w0);
357 return radius * sqrt(8.0 * k / mass) / DispCorrection;
377 if (numberOfSaves > 0 && timeMax > 0 && timeStep > 0)
379 return static_cast<unsigned int>(ceil(
380 (timeMax + timeStep) / timeStep / static_cast<double>(numberOfSaves - 1)));
385 "[Helpers::getSaveCountFromNumberOfSavesAndTimeMaxAndTimeStep()] numberOfSaves: %, timeMax: %, "
386 "timestep: %\n Arguments need to be positive",
387 numberOfSaves, timeMax, timeStep);
425 std::string line_string;
426 getline(in, line_string);
427 out.str(std::move(line_string));
448 file.open(filename.c_str(), std::ios::out);
451 logger(
WARN,
"Error in writeToFile: file could not be opened");
461 std::stringstream ss;
462 for (
int i=0;
i<argc; ++
i) {
463 ss << argv[
i] <<
' ';
472 logger(
WARN,
"[helpers::gnuplot] is not supported on Cygwin");
475 logger(
WARN,
"[helpers::gnuplot] is not supported on Windows");
477 FILE* pipe = popen(
"gnuplot -persist",
"w");
478 fprintf(pipe,
"%s", command.c_str());
487 file.open(filename.c_str(), std::ios::app);
490 logger(
INFO,
"Error in writeToFile: file could not be opened");
504 struct stat stFileInfo;
510 intStat = stat(strFilename.c_str(), &stFileInfo);
541 file.open(filename.c_str(), mode);
557 return mass0 * mass1 / (mass0 + mass1);
563 file.open(filename.c_str(), std::ios::in);
566 logger(
ERROR,
"Error in readArrayFromFile: file could not be opened");
569 std::vector<double> v;
571 for (
int i = 0;
i < n * m;
i++)
585 file.open(filename.c_str(), std::ios::in);
587 logger(
ERROR,
"Error in readArrayFromFile: file could not be opened");
589 for (
unsigned i = 0;
i < nLines;
i++)
591 if (file.eof())
break;
600 const Mdouble logValue = log10(value);
601 const int factor = std::pow(10, precision - std::ceil(logValue));
607 std::ostringstream stm;
608 stm <<
round(value,precision);
621 class LoadingTest :
public DPMBase
630 : species(species), displacement(displacement), velocity(velocity), radius(radius)
642 setMax({radius, radius, radius + radius});
643 setMin({-radius, -radius, 0});
647 speciesHandler.copyAndAddObject(*species);
653 particleHandler.copyAndAddObject(p);
658 wallHandler.copyAndAddObject(w);
664 logger.assert_debug(p,
"Empty particle handler");
668 if (
getTime() <= displacement / velocity)
679 } test(species, displacement, velocity, radius);
682 writeToFile(test.getName() +
".gnu",
"plot '" + test.getName() +
".fstat' u 7:9 w lp");
683 logger(
INFO,
"finished loading test: run 'gnuplot %.gnu' to view output", test.getName());
696 class LoadingTest :
public DPMBase
700 Mdouble tangentialDisplacement;
707 : species(species), displacement(displacement), tangentialDisplacement(tangentialDisplacement),
708 velocity(velocity), radius(radius)
714 setTimeMax(4.0 * tangentialDisplacement / velocity);
720 setMax({radius, radius, radius + radius});
721 setMin({-radius, -radius, 0});
725 speciesHandler.copyAndAddObject(*species);
731 particleHandler.copyAndAddObject(p);
736 wallHandler.copyAndAddObject(w);
742 logger.assert_debug(p,
"Empty particle handler");
746 bool moveRight =
static_cast<int>(
getTime() / (2.0*tangentialDisplacement / velocity) +0.5)%2==0;
750 p->
setPosition({tangentialDisplacement - velocity *
getTime(), 0, radius - displacement});
755 p->
setPosition({-2*tangentialDisplacement + velocity *
getTime(), 0, radius - displacement});
759 } test(species, displacement, tangentialDisplacement, velocity, radius);
762 writeToFile(test.getName() +
".gnu",
"plot '" + test.getName() +
".fstat' u 8:($10*$14) w lp");
763 logger(
INFO,
"finished tangential loading test: run 'gnuplot %.gnu' to view output", test.getName());
775 class ObjectivenessTest :
public DPMBase
779 Mdouble tangentialDisplacement;
786 : species(species), displacement(displacement), tangentialDisplacement(tangentialDisplacement),
787 velocity(velocity), radius(radius)
805 speciesHandler.copyAndAddObject(*species);
811 particleHandler.copyAndAddObject(p);
813 particleHandler.copyAndAddObject(p);
820 logger.assert_debug(p,
"Empty particle handler");
821 logger.assert_debug(q,
"Empty particle handler");
824 if (
getTime() <= tangentialDisplacement / velocity)
828 p->
setPosition({-tangentialDisplacement + velocity *
getTime(), radius - displacement, 0});
831 q->
setPosition({tangentialDisplacement - velocity *
getTime(), -radius + displacement, 0});
835 Mdouble angle = velocity / (radius - displacement) * (
getTime() - tangentialDisplacement / velocity);
851 } test(species, displacement, tangentialDisplacement, velocity, radius);
854 writeToFile(test.getName() +
".gnu",
"set size ratio -1; plot '" + test.getName() +
".fstat' u 14:15 every 2 w lp");
855 logger(
INFO,
"finished objectiveness test: run 'gnuplot %.gnu' to view output", test.getName());
862 int len = is.tellg();
867 is.seekg(len, std::ios_base::beg);
868 logger(
INFO,
"helpers::compare: Next stream value (%) is not %", dummy, s);
876 for (
unsigned i=0;
i<argc; ++
i) {
877 if (varName == argv[
i]) {
878 logger(
INFO,
"readFromCommandLine: % set to true",varName.substr(1));
883 logger(
INFO,
"readFromCommandLine: % set to default value false",varName.substr(1));
889 Mdouble dx = (Max - Min) / static_cast<Mdouble>(numberOfBins - 1);
891 std::vector<Mdouble> linearVector(numberOfBins);
892 typename std::vector<Mdouble>::iterator x;
893 for (x = linearVector.begin(), val = Min; x != linearVector.end(); ++x, val += dx)
898 linearVector.pop_back();
899 linearVector.push_back(Max);
904 void checkTemplate(T real, T ideal,
double error, std::string whatIsChecked)
907 whatIsChecked +
": % (should be %) ", real, ideal);
908 logger(
INFO, whatIsChecked +
": % (correct)", real);
911 void helpers::check(
double real,
double ideal,
double error, std::string whatIsChecked)
933 char cCurrentPath[FILENAME_MAX];
935 logger(
WARN,
"Get current dir failed: %", cCurrentPath);
937 return std::string(cCurrentPath);
943 static auto start = std::chrono::steady_clock::now();
944 auto end = std::chrono::steady_clock::now();
945 std::chrono::duration<double> diff = end - start;
959 auto pos = is.tellg();
970 std::string helpers::readFromCommandLine<std::string>(
int argc,
char *argv[], std::string varName, std::string value)
972 for (
unsigned i=0;
i<argc-1; ++
i) {
973 if (varName == argv[
i]) {
975 logger(
INFO,
"readFromCommandLine: % set to % ", varName.substr(1), value);
980 logger(
INFO,
"readFromCommandLine: % set to default value % ", varName.substr(1), value);
989 nError = _mkdir(path.c_str());
991 nError = mkdir(path.c_str(),nMode);
1000 return constants::pi*radius*sqrt(density/shearModulus)/(0.1631*poisson+0.8766);
MERCURY_DEPRECATED KAndDispAndKtAndDispt computeDisptFromCollisionTimeAndRestitutionCoefficientAndTangentialRestitutionCoefficientAndEffectiveMass(Mdouble tc, Mdouble r, Mdouble beta, Mdouble mass)
Set disp, k, dispt and kt such that is matches a given collision time tc and a normal and tangential ...
bool writeToFile(std::string filename, std::string filecontent)
Writes a string to a file.
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
The DPMBase header includes quite a few header files, defining all the handlers, which are essential...
std::string lower(std::string s)
A spherical particle is the most simple particle used in MercuryDPM.
void setVelocity(const Vec3D &velocity)
set the velocity of the BaseInteractable.
void setTimeMax(Mdouble newTMax)
Sets a new value for the maximum simulation duration.
Mdouble getEffectiveMass(Mdouble mass0, Mdouble mass1)
Calculates the effective mass of a particle pair, i.e. half the harmonic mean of two particle masses...
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
Definition of different loggers with certain modules. A user can define its own custom logger here...
virtual void setRadius(Mdouble radius)
Sets the particle's radius_ (and adjusts the mass_ accordingly, based on the particle's species) ...
void setMax(const Vec3D &max)
Sets the maximum coordinates of the problem domain.
LL< Log::INFO > INFO
Info log level.
Mdouble getRayleighTime(Mdouble radius, Mdouble shearModulus, Mdouble poisson, Mdouble density)
Mdouble exp(Mdouble Exponent)
Mdouble beta(Mdouble z, Mdouble w)
This is the beta function, returns the approximation based on cmath's implementation of ln(gamma) ...
void check(double real, double ideal, double error, std::string errorMessage)
const std::complex< Mdouble > i
void setParticleDimensions(unsigned int particleDimensions)
Sets the particle dimensionality.
virtual const Vec3D & getAngularVelocity() const
Returns the angular velocity of this interactable.
void setSpecies(const ParticleSpecies *species)
Set disp and k such that is matches a given collision time tc and restitution coefficient r for a col...
void normalAndTangentialLoadingTest(const ParticleSpecies *species, Mdouble displacement, Mdouble tangentialDisplacement, Mdouble velocity, Mdouble radius, std::string name)
Used to test the tangential loading/unloading/reloading properties of a frictional contact law...
void setSystemDimensions(unsigned int newDim)
Sets the system dimensionality.
Mdouble log(Mdouble Power)
LL< Log::ERROR > ERROR
Error log level.
bool isEqual(Mdouble v1, Mdouble v2, Mdouble absError)
Compares the difference of two Mdouble with an absolute error, useful in UnitTests.
const unsigned unsignedMax
file will not be created/read
LL< Log::WARN > WARN
Warning log level.
void getLineFromStringStream(std::istream &in, std::stringstream &out)
Reads a line from one stringstream into another, and prepares the latter for reading in...
all data will be written into/ read from a single file called name_
void setMin(const Vec3D &min)
Sets the minimum coordinates of the problem domain.
void setSaveCount(unsigned int saveCount)
Sets File::saveCount_ for all files (ene, data, fstat, restart, stat)
unsigned int getSaveCountFromNumberOfSavesAndTimeMaxAndTimeStep(unsigned int numberOfSaves, Mdouble timeMax, Mdouble timeStep)
Returns the correct saveCount if the total number of saves, the final time and the time step is known...
void loadingTest(const ParticleSpecies *species, Mdouble displacement, Mdouble velocity, Mdouble radius, std::string name)
Used to test the loading/unloading properties of a contact law.
void gnuplot(std::string command)
Plots to a gnuplot window.
bool createDirectory(std::string)
virtual void actionsBeforeTimeStep()
A virtual function which allows to define operations to be executed before the new time step...
bool isNext(std::istream &is, const std::string name)
reads next value in stream as a string and compares it with name.
std::string to_string(const T &n)
void setPosition(const Vec3D &position)
Sets the position of this BaseInteractable.
Mdouble round(const Mdouble value, unsigned precision)
bool readFromCommandLine(int argc, char *argv[], std::string varName)
Returns true if command line arguments contain varName, false else Usage example: if (readFromCommand...
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...
std::vector< Mdouble > linspace(Mdouble a, Mdouble b, int N)
void setTimeStep(Mdouble newDt)
Sets a new value for the simulation time step.
bool compare(std::istream &is, std::string s)
Checks if the next argument in the input stream is a certain string.
This is a class defining walls.
void writeCommandLineToFile(const std::string filename, const int argc, char *const argv[])
Writes a string to a file.
void setOrientation(const Quaternion &orientation)
Sets the orientation of this BaseInteractable.
bool fileExists(std::string strFilename)
Function to check if a file exists, is used to check if a run has already need done.
MERCURY_DEPRECATED Mdouble getMaximumVelocity(Mdouble k, Mdouble disp, Mdouble radius, Mdouble mass)
Calculates the maximum relative velocity allowed for a normal collision of two particles of radius r ...
Implementation of a 3D matrix.
virtual const Vec3D & getVelocity() const
Returns the velocity of this interactable.
T square(const T val)
squares a number
bool openFile(std::fstream &file, std::string filename, std::fstream::openmode mode)
Provides a simple interface for opening a file.
void checkTemplate(T real, T ideal, double error, std::string whatIsChecked)
const Quaternion & getOrientation() const
Returns the orientation of this BaseInteractable.
void objectivenessTest(const ParticleSpecies *species, Mdouble displacement, Mdouble tangentialDisplacement, Mdouble velocity, Mdouble radius, std::string name)
Used to test the tangential loading/unloading/reloading properties of a frictional contact law...
bool addToFile(std::string filename, std::string filecontent)
Adds a string to an existing file.
void setAngularVelocity(const Vec3D &angularVelocity)
set the angular velocity of the BaseInteractble.
Mdouble getTime() const
Returns the current simulation time.
Implementation of a 3D symmetric matrix.
void more(std::string filename, unsigned nLines=constants::unsignedMax)
Mdouble getTimeMax() const
Returns the maximum simulation duration.
std::vector< double > readArrayFromFile(std::string filename, int &n, int &m)
virtual void setupInitialConditions()
This function allows to set the initial conditions for our problem to be solved, by default particle ...
void setSpecies(const ParticleSpecies *species)
Defines the species of the current wall.
void setFileType(FileType fileType)
Sets File::fileType_ for all files (ene, data, fstat, restart, stat)