Merge branch 'main' of https://github.com/fabricix/MPM-Geomechanics

.gitignore

@@ -22,9 +22,8 @@ doxygen/output.txt
 *.sublime-workspace
 *.code-workspace
 *.sublime-project
-simulation_data.csv
 MPM-Geomechanics
 /make/VS/MPM-Geomechanics.vcxproj.user
-simulation_data.csv
-particle_stress_data.json
-time-energy.csv
+time-energy.csv
+simulation-data.csv
+particle_stress_data.json

inc/Input.h

@@ -238,6 +238,7 @@ using Eigen::Vector3d;
 /// | --------| ----------- | --------- |
 /// | time | define the simulation time | double |
 /// | time_step | define the time step | double |
+/// | critical_time_step_multiplier | define the time step by a fraction of the critical time step | double |
 ///
 
 /// ## Simulation time definition example
@@ -246,6 +247,18 @@ using Eigen::Vector3d;
 /// "time":2,
 /// ```
 
+/// ## Time step definition example
+///
+/// ```
+/// "time_step":1e-4,
+/// ```
+///
+/// The time step of the simulation can be defined by the critical time multiplier:
+/// ```
+/// 	"critical_time_step_multiplier":0.25,
+/// ```
+/// So, in this case the time step will be 0.25 of the critical time step.
+
 /// ## Two phase simulations
 ///	| Keyword | Description | Data type |
 /// | --------| ----------- | --------- |
@@ -458,6 +471,24 @@ using Eigen::Vector3d;
 /// }
 /// ```
 
+///
+/// ## Save/Load States
+/// These keywords allow to load and save model states, for example this can be used to load the initial stresses resulting of a previous analysis
+///	| Keyword | Description | Data type |
+/// | --------| ----------- | --------- |
+/// | load_state | activate load model state | bool |
+/// | save_state | activate save model state | bool |
+///
+/// Example
+/// To save the state of the model 
+/// ```
+/// "save_state":true,
+/// ```
+/// And for loading the salved state
+/// ```
+/// "load_state":true,
+/// ```
+
 namespace Input {
 
 	/// \brief Read the input file
@@ -471,6 +502,14 @@ namespace Input {
 	/// \brief Return the file name
 	/// \return File name
 	string getFileName();
+
+	/// \brief Return load state activated
+	/// \return true for load state activated 
+	bool getLoadState();
+
+	/// \brief Return save state activated
+	/// \return true for load state activated
+	bool getSaveState();
 
 	/// \brief Return the simulation time
 	/// \return Simulation time

inc/MPM.h

@@ -94,12 +94,21 @@ class MPM {
 	///
 	void end();
 
+	/// \brief Load model state
+	///
+	void loadState();
+
+	/// \brief Save model state
+	///
+	void saveState();
 private:
 
 	Mesh mesh; //!< grid mesh
 
 	vector<Body*> bodies; //!< bodies discretized by material points
 
+	vector<Particle*> particles; //!< material points list
+
 	vector<Material*> materials; //!< bodies discretized by material points
 
 	Solver* solver; //!< operation for solve the equations in time
@@ -156,6 +165,7 @@ class MPM {
 	///
 	void setNumberPhasesInSimulation();
 
+
 	/// \brief Configure number of threads
 	///
 	void setThreads();

inc/Model.h

@@ -44,6 +44,22 @@ namespace ModelSetup {
 	/// \return True if is Linux system
 	bool getLinuxSystem();
 
+	/// \brief Return true is load state is activated
+	/// \return True is load state is activated
+	bool getLoadState();
+
+	/// \brief Set load state activated
+	/// \param[in] load_state_activated Flag to load an state
+	void setLoadState(bool st);
+
+	/// \brief Return true is save state is activated
+	/// \return True is save state is activated
+	bool getSaveState();
+
+	/// \brief Set save state activated
+	/// \param[in] save_state_activated Flag to save an state
+	void setSaveState(bool st);
+
 	/// \brief Return the number of nodes that a particle contributes
 	/// \return The number of nodes that a particle contributes
 	unsigned getContributionNodesNum();

inc/States.h

@@ -1,27 +1,21 @@
-/*
- * States.h
- *
- *  Created on: 13 de abr de 2021
- *      Author: Fabricio Fernandez <[email protected]>
- */
+
+#include <string>
+#include "Particle/Particle.h"
 
 #ifndef STATES_H_
 #define STATES_H_
 
-/// \class States
-/// \brief Represents model states that
-/// can be loaded before a simulation.
-class States {
+namespace States {
 
-public:
-
-	/// \brief Default constructor
-	///
-	States();
+	/// \brief Save the particles stress into a json file
+	/// \param[in] filename File name
+	/// \param[in] particles List containing pointers to particles
+	void saveParticleStress(const std::string& filename, const std::vector<Particle*>& particles);
 
-	/// \brief Default destructor
-	///
-	virtual ~States();
-};
+	/// \brief Load the particles stress from a json file
+	/// \param[in] filename File name
+	/// \param[in] particles List containing pointers to particles
+	void loadParticleStress(const std::string& filename, std::vector<Particle*>& particles);
+}
 
-#endif /* STATES_H_ */
+#endif /* STATES_H_ */

src/Input.cpp

@@ -65,6 +65,10 @@ inline const json& Input::getJson( ) { return inputFile; }
 
 inline string Input::getFileName( ) { return inputFileName; }
 
+bool Input::getLoadState() { return Input::get_boolean(Input::getJson(), "load_state", false); };
+
+bool Input::getSaveState() { return Input::get_boolean(Input::getJson(), "save_state", false); };
+
 void Input::readInputFile(string filename) {
 
 	try{
@@ -92,7 +96,7 @@ void Input::readInputFile(string filename) {
 double Input::getSimulationTime(){
 
 	try{
-		double time = get_number(inputFile,"time",0);
+		double time = get_number(inputFile,"time",0.0);
 
 		if(time>0) 
 		{ 
@@ -1310,4 +1314,4 @@ vector<Loads::PressureBoundaryForceBox> Input::getPressureBoundaryForceBox() {
 		Warning::printMessage("Error in pressure force box definition");
 		throw;
 	}
-}
+}

src/MPM.cpp

@@ -179,6 +179,10 @@ void MPM::createBodies() {
 
 void MPM::setupParticles() {
 
+	// clear particle list
+	particles.clear();
+
+	// set body id and shape function
 	for (size_t i = 0; i < bodies.size(); ++i){
 
 		vector<Particle*>* particles = bodies.at(i)->getParticles();
@@ -189,8 +193,8 @@ void MPM::setupParticles() {
 			particles->at(j)->setBodyId(bodies.at(i)->getId());
 
 			// set shape function
-			switch(ModelSetup::getInterpolationFunction()){
-
+			switch(ModelSetup::getInterpolationFunction())
+			{
 				case ModelSetup::LINEAR:
 					particles->at(j)->setShape(new ShapeLinear);
 					break;
@@ -203,6 +207,9 @@ void MPM::setupParticles() {
 				Warning::printMessage("Bad definition of shape function in particle");
 				throw (0);
 			}
+
+			// register particle is material point list
+			this->particles.push_back(particles->at(j));
 		}
 	}
 }
@@ -264,6 +271,26 @@ void MPM::setNumberPhasesInSimulation() {
 	ModelSetup::setTwoPhaseActive(Input::getNumberPhases()>1?true:false);
 }
 
+void MPM::loadState()
+{
+	ModelSetup::setLoadState(Input::getLoadState());
+
+	if ( ModelSetup::getLoadState() ) 
+	{
+		States::loadParticleStress("particle_stress_data.json", particles);
+	}
+}
+
+void MPM::saveState()
+{
+	ModelSetup::setSaveState(Input::getSaveState());
+
+	if ( ModelSetup::getSaveState() )
+	{
+		States::saveParticleStress("particle_stress_data.json", particles);
+	}
+}
+
 void MPM::createModel() {
 
 	try{
@@ -300,6 +327,9 @@ void MPM::createModel() {
 		// configures the particles in the model
 		setupParticles();
 
+		// load previous state
+		loadState();
+
 		// configures the loads
 		setupLoads();
 
@@ -332,10 +362,7 @@ void MPM::solve(){
 	else{
 		ModelSetup::setInitialSimulationTime(std::chrono::system_clock::now());
 		solver->Solve();
-		Output::printElapsedTime();
 	}
 }
 
-void MPM::end(){
-
-}
+void MPM::end() { Output::printElapsedTime(); }

src/Main.cpp

@@ -32,6 +32,9 @@ int main(int argc, char **argv)
 	// solve the problem in time
 	mpm.solve();
 
+	// save state
+	mpm.saveState();
+
 	// finish the program
 	mpm.end();
 

src/Model.cpp

@@ -24,6 +24,10 @@ namespace ModelSetup {
 
 	// contact
 	bool contactActive=false; //!< is contact active
+
+	// states
+	bool loadState = false; //!< if state need to be loaded
+	bool saveState = false; //!< if state need to be salved
 
 	// gravity active
 	bool gravityActive=false; //!< is gravity active
@@ -49,15 +53,15 @@ namespace ModelSetup {
 	// time
 	double dt=0.0; //!< time step
 	double time=0.0; //!< simulation time
-	double dt_critical_multiplier=0.25; //!< dafault critical time step fraction
+	double dt_critical_multiplier = 0.25; //!< dafault critical time step fraction
 
 	// initial time simulation
     std::chrono::system_clock::time_point initialSimulationTime;
 
 	// default damping
 	DampingType damping=DampingType::UNDAMPED; //!< damping type
 
-	// default lodal damping value
+	// default local damping value
 	double localDamping=0.0; //!< local damping value
 
 	// input file
@@ -73,6 +77,12 @@ namespace ModelSetup {
 	/// Function members
 	///
 
+    // states 
+	bool getLoadState() { return loadState; }
+	void setLoadState(bool st) { loadState = st; }
+	bool getSaveState() { return saveState; }
+	void setSaveState(bool st) { saveState = st; }
+
 	// contribution nodes
 	unsigned getContributionNodesNum() { return contributionNodes; }
 
@@ -82,8 +92,7 @@ namespace ModelSetup {
 
 	// results
 	unsigned getResultNum() { return resultNumber; }
-	void
-	setResultNum(unsigned d) { resultNumber=d; }
+	void setResultNum(unsigned d) { resultNumber=d; }
 	unsigned getResultSteps() { return static_cast<unsigned int>(floor(time/dt)/resultNumber); }
 
 	// time step
@@ -115,6 +124,7 @@ namespace ModelSetup {
 	bool getTwoPhaseActive() { return twoPhaseCalculationActive; }
 	void setTwoPhaseActive(bool d) { twoPhaseCalculationActive=d; }
 
+    // simulation time
     void setInitialSimulationTime(std::chrono::system_clock::time_point initialTime) { ModelSetup::initialSimulationTime = initialTime; }
     std::chrono::system_clock::time_point getInitialSimulationTime() { return ModelSetup::initialSimulationTime;  }
 
@@ -143,7 +153,8 @@ namespace ModelSetup {
 	// interpolation functions
 	ModelSetup::InterpolationFunctionType getInterpolationFunction() { return interpolationType; }
 	void setInterpolationFunction(ModelSetup::InterpolationFunctionType d) { interpolationType=d; }
-
+
+	// openMP threads
 	void setNumThreads(unsigned nThreads){
 
 		#ifdef _OPENMP

src/Output.cpp

@@ -601,7 +601,7 @@ namespace Output{
     double ienergy = DynamicRelaxation::computeKineticEnergy(bodies);
 
     // Open the simulation CSV file in append mode
-    std::ofstream csv_file("simulation_data.csv", std::ios::app);
+    std::ofstream csv_file("time-energy.csv", std::ios::app);
     if (!csv_file.is_open()) {
         std::cerr << "Error opening the CSV file" << std::endl;
     }
@@ -615,7 +615,7 @@ namespace Output{
 
 	void writeResultInStep(int loopCounter, int resultSteps,vector<Body*>* bodies, double iTime)
 	{
-		if (iTime == 0) { printModelInfo(bodies, iTime); initializeCSVFile("simulation_data.csv");}
+		if (iTime == 0) { printModelInfo(bodies, iTime); initializeCSVFile("time-energy.csv");}
 
 		if (loopCounter%resultSteps==0)
 		{