32 logger(
DEBUG,
"BaseCluster::BaseCluster() finished");
39 logger(
DEBUG,
"BaseCluster::BaseClusterr() finished");
77 if (cTOTS < 45 && cTOTS > 0)
78 logger(
WARN,
"collisionTimeOverTimeStep = % is too small: consider setting it greater or equal than 50.", cTOTS);
80 logger(
ERROR,
"collisionTimeOverTimeStep = % must be grater than zero.", cTOTS);
99 logger(
ERROR,
"radiusParticle must be greater than zero. radiusParticle = %", rP);
119 logger(
ERROR,
"sizeDispersityParticle must be greater or equal than 1. sizeDispersityParticle = %", sDP);
136 logger(
ERROR,
"nParticles must be grater than zero. nParticles = %", nP);
152 "relativeClusterRadius is smaller than 0. relativeClusterRadius = %",
185 logger(
WARN,
"idCluster = % is less than zero.", iC);
201 if (vDM < 0 || vDM > 1)
202 logger(
ERROR,
"velocityDampingModulus must be grater than zero and less than 1. velocityDampingModulus = %", vDM);
219 logger(
ERROR,
"nInternalStructurePoints_ must be grater than zero. nInternalStructurePoints_ = %", gL);
236 logger(
ERROR,
"energyRatioTolerance must be grater than zero. energyRatioTolerance = %", eRT);
430 logger(
WARN,
"ParticleHandler was not empty");
442 logger(
ERROR,
"Please set exactly two values among radiusParticle, radiusCluster and numberOfParticles."
443 " radiusParticle = %, radiusCluster = %, numberOfParticles = %.",
450 logger(
VERBOSE,
"clusterRadius is small compared to the radius of a single particle:"
451 " consider setting clusterRadius >= 2 * radiusParticle."
516 #ifdef MERCURY_USE_MPI
548 #ifdef MERCURY_USE_MPI
587 std::ostringstream name;
702 #ifdef MERCURY_USE_MPI
745 logger(
WARN,
"Dissipation process not completed: final energyRatioTollerance_ = %."
746 "Try to increase energyRatioTollerance_ or decreasing velocityDampingModulus.",
793 writeAllParticles =
true;
798 "stage " <<
stage_ <<
" " <<
842 std::stringstream line;
912 std::ostringstream cdatName;
913 cdatName <<
getName() <<
".cdat";
914 cdatFile_.open(cdatName.str(), std::ios::app);
919 std::ostringstream overlName;
920 overlName <<
getName() <<
".overl";
921 overlFile_.open(overlName.str(), std::ios::app);
955 case 3: printTime <<
"Dissipating energy: ";
958 default: printTime <<
"Final values: ";
1020 std::ostringstream printSpecies;
1026 printSpecies <<
"loadingStiffness: " << std::scientific <<
particleSpecies_ -> getLoadingStiffness()
1028 <<
"unloadingStiffnessMax: " <<
particleSpecies_ -> getUnloadingStiffnessMax()
1046 Mdouble initialSolidFraction = 0.1;
1048 std::ostringstream printDomainLimits;
1049 printDomainLimits <<
"Cubic size " <<
boxSize_ << std::endl;
1065 std::ostringstream printTimeStep;
1066 printTimeStep <<
"timeStep: " << std::setprecision(4) <<
getTimeStep() << std::endl
1067 <<
"cT/tS, at least: " << std::fixed << std::setprecision(1)
1079 int nParticlesInserted = 0;
1090 nParticlesInserted++;
1094 logger(
ERROR,
"Cannot insert all particles, try to decrase the value of initialSolidFraction in "
1095 "BaseCluster::setDomainLimits();\n"
1096 "Inserted %/% particles.", nParticlesInserted,
nParticles_);
1099 logger(
VERBOSE,
"PARTICLE INSERTION TERMINATED SUCCESSFULLY\n");
1138 std::ostringstream cdatName;
1139 cdatName <<
getName() <<
".cdat";
1141 cdatFile_.open(cdatName.str(), std::ios::out);
1143 cdatFile_ <<
"CLUSTER DATA AND INFORMATION" << std::endl << std::endl;
1169 cdatFile_ <<
"progress" << std::setw(16) <<
"ElastEne" << std::setw(23) <<
"Ekin/ElastEne" << std::setw(18) <<
"coord_number" << std::setw(14) <<
"meanRadius"
1170 << std::setw(19) <<
"solidFraction" << std::setw(16) <<
"forceModulus" << std::setw(22) <<
"AveFOnOverl" << std::setw(15) <<
"dMin" << std::setw(17) <<
"dMean"
1171 << std::setw(14) <<
"dMax" << std::setw(24) <<
"Mass Centre" << std::endl;
1181 std::ostringstream overlName;
1182 overlName <<
getName() <<
".overl";
1184 overlFile_.open(overlName.str(), std::ios::out);
1185 overlFile_ <<
"Overlap Vs Normal Force" << std::endl;
1200 int insertionFailCounter = 0;
1202 Vec3D particlePosition;
1211 while (insertionFailCounter < 1000)
1229 insertionFailCounter++;
1258 Mdouble solidVolumeInsideRadius = 0;
1265 Vec3D distanceFromCenterOfMass;
1267 Mdouble furthestParticleDistance = 0;
1283 centerOfMass_ += ((*p)->getVolume()) * ((*p)->getPosition());
1293 distanceFromCenterOfMass = (*p)->getPosition() -
centerOfMass_;
1295 if (distanceFromCenterOfMass.
getLength() > furthestParticleDistance)
1296 furthestParticleDistance = distanceFromCenterOfMass.
getLength();
1302 distanceFromCenterOfMass = (*p)->getPosition() -
centerOfMass_;
1330 solidVolumeInsideRadius += (*p) -> getVolume();
1343 relativeOverlap = ((*i) -> getOverlap())/(
particleHandler.
getObject((*i) -> getP() -> getIndex()) -> getRadius()) +
1345 relativeOverlap /= 2;
1392 std::setprecision(2)<<
1398 std::setprecision(2) <<
1405 std::setprecision(2) <<
1409 std::setprecision(3) <<
1414 std::setprecision(5) <<
1422 std::setprecision(2) <<
1462 forceOnOverlap = ((*i) -> getForce()).getLength();
1463 relativeOverlap = ((*i) -> getOverlap())/(
particleHandler.
getObject((*i) -> getP() -> getIndex()) -> getRadius()) +
1465 relativeOverlap = relativeOverlap / 2;
1466 overlFile_ << std::setprecision(2) << std::scientific << std::setw(18) << forceOnOverlap
1467 << std::defaultfloat << std::fixed << std::setprecision(4) << std::setw(9) << relativeOverlap;
1543 std::vector<int> temporaryRowVector;
1547 temporaryRowVector.push_back(0);
1554 adjacencyMatrix_[(*i) -> getP() -> getIndex()][(*i) -> getI() -> getIndex()] = 1;
1555 adjacencyMatrix_[(*i) -> getI() -> getIndex()][(*i) -> getP() -> getIndex()] = 1;
1564 std::ostringstream amatName;
1565 amatName <<
getName() <<
".amat";
1567 amatFile_.open(amatName.str(), std::ios::out);
1568 amatFile_ <<
"ADJACENCY MATRIX" << std::endl << std::endl;
1609 Mdouble fictitiousGridPointRadiusRatio = 1.0e-5;
1614 Mdouble nPointsInsideComponentsForMCTest = 0;
1618 for (
int i = 0;
i < nMonteCarloSamplingPoints; ++
i)
1625 mcPoint.
X = rad *
sin(phi) *
cos(theta);
1626 mcPoint.
Y = rad *
sin(phi) *
sin(theta);
1627 mcPoint.
Z = rad *
cos(phi);
1634 nPointsInsideComponentsForMCTest++;
1635 intStructFile_ << std::scientific << std::setprecision(5) << std::setw(12) << mcPoint.
X
1636 << std::setw(13) << mcPoint.
Y << std::setw(13) << mcPoint.
Z << std::setw(6) << 0 << std::endl;
1647 Mdouble accordance = (theoVal - diff)/theoVal;
1654 Mdouble accordanceAvOverl = (theoValAvOverl - diffAvOverl)/theoValAvOverl;
1656 intStructFile_ <<
"n_points_inside_boundary: " << std::scientific << nMonteCarloSamplingPoints << std::endl;
1657 intStructFile_ <<
"n_points_inside_components: " << nPointsInsideComponentsForMCTest << std::endl;
1659 <<
", accordance with theoretical values: " << 100*accordance <<
"%." << std::endl
1660 <<
"Accordance with average overlap: " << 100*accordanceAvOverl <<
"%." << std::endl
1661 <<
"It is very important to notice that this formula is accurate only if sliding friction" << std::endl
1662 <<
"is set to 0.5 and relative tangential stiffness is set to 0.3 while creating the cluster." << std::endl
1663 <<
"Different values do not guarantee accuracy." << std::endl << std::endl;
1670 std::ostringstream printResults;
1671 printResults <<
"n_points_inside_boundary: " << std::scientific << nMonteCarloSamplingPoints << std::endl;
1672 printResults <<
"n_points_inside_components: " << nPointsInsideComponentsForMCTest << std::endl;
1674 <<
", accordance with theoretical values: " << 100*accordance <<
"%." << std::endl
1675 <<
"Accordance with average overlap: " << 100*accordanceAvOverl <<
"%." << std::endl
1676 <<
"It is very important to notice that this formula is accurate only if sliding friction" << std::endl
1677 <<
"is set to 0.5 and relative tangential stiffness is set to 0.3 while creating the cluster." << std::endl
1678 <<
"Different values do not guarantee accuracy." << std::endl << std::endl;
1690 std::ostringstream gnuplotname;
1691 gnuplotname <<
getName() <<
".gnuplot";
1696 std::string titleLine=R
"(set title "Overlap Vs Force")";
1697 std::string xLabel=R
"(set xlabel "Overlap")";
1698 std::string yLabel=R"(set ylabel "Force")";
1706 gnuplotFile_ <<
"\"" <<
getName() <<
".overl" <<
"\"" <<
" using " << 2*
i+4 <<
":" << 2*
i+3 <<
" title \"\" with lines lt 1 dashtype 2, ";
1721 std::ostringstream intStructName;
1722 intStructName <<
getName() <<
".struct";
std::vector< std::vector< int > > adjacencyMatrix_
Vec3D getPosition() const
This returns the value of position_, which is the position in which the cluster will be inserted...
void setSpecies()
Sets species of particles.
BaseParticle * getLargestParticle() const
Returns the pointer of the largest particle in the particle handler. When mercury is running in paral...
void computeInternalStructure()
This computes the internal structure of the cluster.
void setCollisionTimeOverTimeStep(Mdouble cTOTS)
This sets the collisionTimeOverTimeStep number (which is the ratio between collision time and time st...
Mdouble minRelativeOverlap_
Mdouble getCollisionTimeOverTimeStep() const
This returns the value of the ratio between collision time and time step.
Mdouble X
the vector components
void makeDataAnalysis()
This functions computes some important cluster information needed by the program. ...
bool isEneOutputOn() const
This returns the bool variable that defines whether the cluster ene output is written or not...
void writeToCdatFile()
This writes on the cluster data output file.
A spherical particle is the most simple particle used in MercuryDPM.
void setVelocity(const Vec3D &velocity)
set the velocity of the BaseInteractable.
LinearPlasticViscoelasticFrictionSpecies * particleSpecies_
unsigned int getSize() const
Gets the size of the particleHandler (including mpi and periodic particles)
void setTimeMax(Mdouble newTMax)
Sets a new value for the maximum simulation duration.
void setNumberOfInternalStructurePoints(int gL)
This sets the value of the number of particles used to compute the internal structure.
LinearPlasticViscoelasticFrictionSpecies * getParticleSpecies() const
This returns the species of the particle.
void actionsAfterSolve() override
Overrides DPMBase actionsAfterSolve(): in this cluster data file and cluster overlap file are closed ...
void createAdjacencyMatrix()
This calculates the adjacency matrix of the cluster.
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 doOverlOutput(bool iOOO)
This sets the bool variable that defines whether the cluster overlap output will be written or not...
Mdouble velocityDampingInterval_
void doIntStrucOutput(bool iISOO)
This sets the bool variable that defines whether the cluster internal structure output will be writte...
void doAmatOutput(bool iAOO)
This sets the bool variable that defines whether the cluster adjacency matrix output will be written ...
unsigned int getClusterId() const
This returns the value of the cluster ID.
void doEneOutput(bool isEneOutputOn)
This sets the bool variable that defines whether the cluster ene output will be written or not...
void insertParticles()
Inserts particles inside the domain.
void setParticlesWriteVTK(bool writeParticlesVTK)
Sets whether particles are written in a VTK file.
void writeToOverlFile()
This writes on the cluster overlap output file.
void decreaseForce()
This linearly decreases values of forceModulus (stage = 2).
const std::complex< Mdouble > i
Mdouble meanRelativeOverlap_
void setRadii()
Sets all radii according to particleRadius and sizeDispersityParticle.
const std::vector< T * >::const_iterator end() const
Gets the end of the const_iterator over all BaseBoundary in this BaseHandler.
void setSizeDispersityParticle(Mdouble sDP)
This sets the value of particles' dispersity in size.
const std::string & getName() const
Returns the name of the file. Does not allow to change it though.
bool isVtkOutputOn() const
This returns the bool variable that defines whether the cluster vtk output is written or not...
void setZero()
Sets all elements to zero.
void dampForce()
This damps values of forceModulus (stage = 3).
bool setNumberOfParticles_
Mdouble maxRelativeOverlap_
bool isIntStrucOutputOn() const
This returns the bool variable that defines whether the cluster internal structure output is written ...
int getNumberOfInternalStructurePoints() const
This returns the value of the number of particles used to compute internal structure.
void setSpecies(const ParticleSpecies *species)
Mdouble dissipationDuration_
Mdouble getRandomNumber()
This is a random generating routine can be used for initial positions.
Mdouble getMassFromRadius(Mdouble radius) const
bool particleInsertionSuccessful(int n)
This function tries to insert the n-th particle (returns true if it manage to do that). It is inside insertParticles().
Mdouble getMeanClusterRadius()
this returns meanClusterRadius (radius of an ideal perfectly spherical cluster, there's no setter)...
bool isAmatOutputOn() const
This returns the bool variable that defines whether the cluster adjacency matrix output is written or...
std::ofstream gnuplotFile_
void increaseForce()
This linearly increases the value of forceModulus (stage = 1).
void clear() override
Empties the whole ParticleHandler by removing all BaseParticle.
void doCdatOutput(bool iCOO)
This sets the bool variable that defines whether the cluster data output will be written or not...
void makeAmatFile()
This creates the adjacency matrix file.
bool isCdatOutputOn() const
This returns the bool variable that defines whether the cluster data output (which is NOT the mercury...
static Mdouble getLength(const Vec3D &a)
Calculates the length of a Vec3D: .
void read(std::istream &is, ReadOptions opt=ReadOptions::ReadAll) override
Overrides DPMBase read(): in this all variables needed by the program for restarting are read...
const std::vector< T * >::const_iterator begin() const
Gets the begin of the const_iterator over all Object in this BaseHandler.
Mdouble getSizeDispersityParticle() const
This returns the value of particles' dispersity in size.
void setVelocity(Vec3D v)
This sets the value of velocity after creation.
bool checkParticleForInteraction(const BaseParticle &P) final
Checks if given BaseParticle has an interaction with a BaseWall or other BaseParticle.
Mdouble meanCoordinationNumber_
void setDomain(const Vec3D &min, const Vec3D &max)
Sets the minimum coordinates of the problem domain.
Mdouble getElasticEnergy() const
Returns the global elastic energy within the system.
void setDomainLimits()
Sets domain limits.
File dataFile
An instance of class File to handle in- and output into a .data file.
file will not be created/read
void makeCdatFile()
Creates the cluster data output file.
Mdouble maximumForceModulus_
void write(std::ostream &os, bool writeAllParticles=true) const override
Writes all data into a restart file.
Mdouble fileOutputTimeInterval_
File fStatFile
An instance of class File to handle in- and output into a .fstat file.
void write(std::ostream &os, bool writeAllParticles) const override
Overrides DPMBase write(): in this all variables needed by the program for restarting are written...
std::ofstream intStructFile_
Mdouble getFinalMassFraction()
This gets the final value obtained for the mass fraction;.
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...
void actionsOnRestart() override
Overrides DPMBase actionsOnRestart(): in this all variables needed by the program for restarting are ...
void setVelocityDampingModulus(Mdouble vDM)
This sets the value of the velocity damping modulus.
void makeOverlFile()
Creates the cluster overlap output file.
ParticleHandler particleHandler
An object of the class ParticleHandler, contains the pointers to all the particles created...
all data will be written into/ read from a single file called name_
void makeGnuplotFile()
This creates the gnuplot file needed for printing force vs overlaps values.
T * getObject(const unsigned int id)
Gets a pointer to the Object at the specified index in the BaseHandler.
void doVtkOutput(bool iVOO)
This sets the bool variable that defines whether the cluster vtk output will be written or not...
Mdouble getKineticEnergy() const
Returns the global kinetic energy stored in the system.
void printTime() const override
Overrides DPMBase printTime(): this way variables of interest are shown.
void doFStatOutput(bool isfStatOutputOn)
This sets the bool variable that defines whether the cluster fStat output will be written or not...
void setSaveCount(unsigned int saveCount)
Sets File::saveCount_ for all files (ene, data, fstat, restart, stat)
Mdouble collisionTimeOverTimeStep_
Mdouble energyRatioTolerance_
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.
Mdouble radiusForSolidFraction_
void writeAmatFile()
This writes on the adjacency matrix file.
Mdouble forceDampingModulus_
Mdouble getRadius() const
Returns the particle's radius.
Mdouble getAverageOverlap()
this returns the average overlap.
void clear() override
Empties the whole BaseHandler by removing all Objects and setting all other variables to 0...
Mdouble sizeDispersityParticle_
SpeciesHandler speciesHandler
A handler to that stores the species type i.e. LinearViscoelasticSpecies, etc.
InteractionHandler interactionHandler
An object of the class InteractionHandler.
Mdouble totalParticleVolume_
void setXBallsAdditionalArguments(std::string newXBArgs)
Set the additional arguments for xballs.
~BaseCluster() final
Default destructor.
void setFileType(FileType fileType)
Sets the type of file needed to write into or read from. File::fileType_.
RNG random
This is a random generator, often used for setting up the initial conditions etc...
Vec3D getVelocity()
This gets the value of velocity after creation.
Mdouble solidFractionIntStruct_
void setGroupId(unsigned groupId)
void dampVelocities()
This damps values of each particle velocity (stage = 1, stage = 2, stage = 3).
void setPosition(const Vec3D &position)
Sets the position of this BaseInteractable.
Mdouble round(const Mdouble value, unsigned precision)
void setNumberOfParticles(int nP)
This sets the value of the number of particles in the cluster.
void setName(const std::string &name)
Allows to set the name of all the files (ene, data, fstat, restart, stat)
void applyCentralForce()
This applies force on each particle.
Mdouble meanClusterRadius_
bool isRestartOutputOn() const
This returns the bool variable that defines whether the cluster restart output is written or not...
Mdouble velocityDampingModulus_
Mdouble getEnergyRatioTolerance() const
This returns the value of the value of the energy ratio threshold under which the process can be cons...
Mdouble forceTuningDuration_
void makeIntenalStructureFile()
This creates the file needed for writing down datas from computeInternalStructure().
void setParticleSpecies(LinearPlasticViscoelasticFrictionSpecies *particleSpecies)
This sets the species of the particle.
void setTimeStep(Mdouble newDt)
Sets a new value for the simulation time step.
Mdouble getRadiusParticle() const
This returns the value of particles' radius if there's no dispersity in size. In case of dispersity !...
File eneFile
An instance of class File to handle in- and output into a .ene file.
File restartFile
An instance of class File to handle in- and output into a .restart file.
int nInternalStructurePoints_
void setEnergyRatioTolerance(Mdouble eRT)
This sets the value of the value of the energy ratio threshold under which the process can be conside...
Contains material and contact force properties.
void doRestartOutput(bool isRestartOutputOn)
This sets the bool variable that defines whether the cluster restart output will be written or not...
void setRadiusParticle(Mdouble rP)
This sets the value of particles' radius if there's no dispersity in size.
void actionsAfterTimeStep() override
Overrides DPMBase actionsAfterTimeStep(): in this compression and decompression are computed...
void calculateTimeStep()
Calculates the time step.
Mdouble getTimeStep() const
Returns the simulation time step.
Mdouble getVelocityDampingModulus() const
This returns the value of the velocity damping modulus.
bool isFStatOutputOn() const
This returns the bool variable that defines whether the cluster fStat output is written or not...
void read(std::istream &is, ReadOptions opt=ReadOptions::ReadAll) override
Reads the MercuryBase from an input stream, for example a restart file.
bool readRestartFile(ReadOptions opt=ReadOptions::ReadAll)
Reads all the particle data corresponding to a given, existing . restart file (for more details regar...
void setupInitialConditions() override
Overrides DPMBase setupInitialConditions(): in this initial conditions for the problem are set...
Mdouble getTime() const
Returns the current simulation time.
Mdouble forceTuningInterval_
bool isOverlOutputOn() const
This returns the bool variable that defines whether the cluster overlap output is written or not...
void setPosition(Vec3D p)
This sets the value of position_, which is the position in which the cluster will be inserted...
void setClusterId(unsigned int iC)
This sets the value of the cluster ID.
BaseCluster()
Default constructor.
int getNumberOfParticles() const
This returns the value of the number of particles in the cluster.
void setRadiusCluster(Mdouble rCR)
This sets the desired value of the cluster radius (there is no getter of this value, but there is a getter of the actual mean cluster radius obtained, getMeanClusterRadius)
std::vector< Mdouble > radii_