Merge pull request #9 from fabricix/seismic-acceleration

Seismic acceleration

.gitignore

@@ -26,4 +26,5 @@ MPM-Geomechanics
 /make/VS/MPM-Geomechanics.vcxproj.user
 time-energy.csv
 simulation-data.csv
-particle_stress_data.json
+particle_stress_data.json
+base_acceleration.csv

doxygen/Doxyfile

@@ -51,7 +51,7 @@ PROJECT_BRIEF          = "This project provides a platform for developing the Ma
 # pixels and the maximum width should not exceed 200 pixels. Doxygen will copy
 # the logo to the output directory.
 
-PROJECT_LOGO           = logo.png
+PROJECT_LOGO           = 
 
 # The OUTPUT_DIRECTORY tag is used to specify the (relative or absolute) path
 # into which the generated documentation will be written. If a relative path is

inc/Boundary.h

@@ -45,7 +45,7 @@ class Boundary {
 
 	/// \enum BoundaryType
 	/// \brief Determines the type of restrictions to be imposed to the mesh.
-	enum BoundaryType{ FREE, FIXED, SLIDING };
+	enum BoundaryType{ FREE, FIXED, SLIDING, EARTHQUAKE };
 
 	/// \enum BoundaryPlane
 	/// \brief Planes at the mesh boundary

inc/Input.h

@@ -614,6 +614,8 @@ namespace Input {
 	/// \brief Return pore pressure force in particles inside a box
 	/// \return pressure_force  Pore pressure force
 	vector<Loads::PressureBoundaryForceBox> getPressureBoundaryForceBox();
+
+	Loads::SeismicData readSeismicData(const std::string& filename, bool hasHeader);
 };
 
 #endif /* INPUT_H_ */

inc/Interpolation.h

@@ -8,6 +8,9 @@
 #ifndef INTERPOLATION_H_
 #define INTERPOLATION_H_
 
+#include "Eigen/Core"
+using Eigen::Vector3d;
+
 #include<vector>
 using std::vector;
 
@@ -144,6 +147,9 @@ namespace Interpolation {
 	/// \param[in] bodies A list o Body pointers
 	/// \param[in] time_step Time step
 	void particleDeformationGradient(Mesh* mesh, vector<Body*>* bodies, double time_step);
+
+	Eigen::Vector3d interpolateVector(const std::vector<double>& times, const std::vector<Eigen::Vector3d>& values, double itime);
 };
 
-#endif /* INTERPOLATION_H_ */
+#endif /* INTERPOLATION_H_ */
+

inc/Loads.h

@@ -13,6 +13,9 @@
 #include <vector>
 using std::vector;
 
+#include "Eigen/Core"
+using Eigen::Vector3d;
+
 /// \namespace Loads
 /// \brief Operations to manage loads in the MPM model.
 namespace Loads {
@@ -77,7 +80,16 @@ namespace Loads {
         }
     };
 
-	/// \brief Configures the gravity load in particles
+    // TODO doxygen
+    struct SeismicData {
+        std::vector<double> time;
+        std::vector<Eigen::Vector3d> acceleration;
+        std::vector<Eigen::Vector3d> velocity;
+    };
+
+    Loads::SeismicData& getSeismicData();
+
+    /// \brief Configures the gravity load in particles
 	/// \param[in] bodies A list containing Body pointers 
 	void setGravity(vector<Body*>& bodies);
 
@@ -115,6 +127,9 @@ namespace Loads {
     /// \brief Set initial velocity in bodies
     /// \param[in] bodies A list containing Body pointers
     void setInitialVelocity(vector<Body*>& bodies);
+
+    void setSeismicData();
+    SeismicData& Loads::getSeismicData();
 };
 
 #endif /* LOADS_H_ */

inc/MPM.h

@@ -157,7 +157,7 @@ class MPM {
 	///
 	void setupResults();
 
-	/// \brief Configure dampig forces
+	/// \brief Configure damping forces
 	/// 
 	void setupDamping();
 
@@ -169,6 +169,9 @@ class MPM {
 	/// \brief Configure number of threads
 	///
 	void setThreads();
+
+	bool getSeismicAnalysis();
+	void setSeismicAnalysis(bool);
 };
 
 #endif /* MPM_H_ */

inc/Model.h

@@ -201,6 +201,14 @@ namespace ModelSetup {
 	/// \brief Configure the number of threads in the model
 	/// \param[in] nThreads Number of threads
 	void setNumThreads(unsigned nThreads);
+
+	/// \brief Return true if seismic analysis is active
+	/// \param[out] if_seismic_analysis_active 
+	bool getSeismicAnalysis();
+
+	/// \brief Configure if seismic analysis is active
+	/// \param[in] seismic_analysis_active 
+	void setSeismicAnalysis(bool seismic_analysis_active);
 
 	/// \brief Return if two-phase calculation is active
 	/// \return True if two-phase calculation is active
@@ -217,6 +225,16 @@ namespace ModelSetup {
 	/// \brief Get initial simulation time
 	/// \param[out] initialTime
 	std::chrono::system_clock::time_point getInitialSimulationTime( );
+
+	std::string getSeismicFileName();
+
+	/// \brief Get current simulation time
+	/// \param[out] currentTime
+	double getCurrentTime();
+
+    /// \brief Set current simulation time
+	/// \param[in] currentTime 
+	void setCurrentTime(double a);
 };
 
 #endif /* MODEL_H_ */

inc/Update.h

@@ -10,6 +10,7 @@
 
 #include "Mesh/Mesh.h"
 #include "Body/Body.h"
+#include <Loads.h>
 
 /// \namespace Update
 /// \brief Represents operations to update values in nodes and particles.
@@ -101,7 +102,7 @@ namespace Update {
 	/// \param[in] bodies List of Body pointers
 	/// \param[in] time_step Time step
 	void particlePosition(Mesh* mesh, vector<Body*>* bodies, double time_step);
-	
+
 	/// \brief Apply essential boundary condition in 
 	/// terms of force
 	/// \param[in] mesh Mesh reference
@@ -152,6 +153,7 @@ namespace Update {
 	/// \param[in] direction Direction to apply de boundary condition 
 	/// \f$x=0\f$, \f$y=1\f$ , \f$z=2\f$ 
 	void setPlaneMomentumFluid(const Boundary::planeBoundary* boundary, vector<Node*>* nodes, unsigned direction);
+
 };
 
 #endif /* UPDATE_H_ */

src/Input.cpp

@@ -31,6 +31,14 @@ using std::ifstream;
 using std::string;
 using std::to_string;
 
+#include <vector>
+#include <string>
+#include <sstream>
+#include <iostream>
+#include <algorithm>
+using namespace Loads;
+
+
 namespace Input {
 
 	json inputFile; //!< data structure containing all the model information
@@ -939,6 +947,11 @@ static void setRestriction(size_t index,vector<Boundary::BoundaryType>& restrict
 	{
 		restrictions.at(index)=Boundary::BoundaryType::SLIDING;
 	}
+	else if (resPlane=="earthquake")
+	{
+		restrictions.at(index)=Boundary::BoundaryType::EARTHQUAKE;
+		ModelSetup::setSeismicAnalysis(true);
+	}
 	else
 	{
 		throw(0);
@@ -1168,7 +1181,7 @@ vector<Loads::PressureBox> Input::getPrescribedPressureBox() {
 	}
 	catch(...)
 	{
-		Warning::printMessage("Error in prescibed pressure box definition");
+		Warning::printMessage("Error in prescribed pressure box definition");
 		throw;
 	}
 }
@@ -1265,40 +1278,40 @@ vector<Loads::PressureMaterial> Input::getInitialPressureMaterial() {
 
 vector<Loads::PressureBoundaryForceBox> Input::getPressureBoundaryForceBox() {
 
-	try{
+	try {
 		vector<Loads::PressureBoundaryForceBox> pressure_box_list;
 
 		string key = "pressure_force_box";
 
-		if (!inputFile[key].is_null()){
-		
+		if (!inputFile[key].is_null()) {
+
 			json::iterator it;
-			for(it = inputFile[key].begin(); it!=inputFile[key].end();it++) {
+			for (it = inputFile[key].begin(); it != inputFile[key].end(); it++) {
 
 				Loads::PressureBoundaryForceBox ipressure_box;
 
-				Vector3d p1(0,0,0), p2(0,0,0), force(0,0,0);
-				
-				if ((*it)["point_p1"].is_array()) { 
-					p1.x() = ((*it)["point_p1"])[0]; 
-					p1.y() = ((*it)["point_p1"])[1]; 
+				Vector3d p1(0, 0, 0), p2(0, 0, 0), force(0, 0, 0);
+
+				if ((*it)["point_p1"].is_array()) {
+					p1.x() = ((*it)["point_p1"])[0];
+					p1.y() = ((*it)["point_p1"])[1];
 					p1.z() = ((*it)["point_p1"])[2];
 
 					ipressure_box.pointP1 = p1;
 				}
 
-				if ((*it)["point_p2"].is_array()) { 
-					p2.x() = ((*it)["point_p2"])[0]; 
-					p2.y() = ((*it)["point_p2"])[1]; 
-					p2.z() = ((*it)["point_p2"])[2]; 
+				if ((*it)["point_p2"].is_array()) {
+					p2.x() = ((*it)["point_p2"])[0];
+					p2.y() = ((*it)["point_p2"])[1];
+					p2.z() = ((*it)["point_p2"])[2];
 
 					ipressure_box.pointP2 = p2;
 				}
 
-				if ((*it)["pressure_force"].is_array()) { 
-					force.x() = ((*it)["pressure_force"])[0]; 
-					force.y() = ((*it)["pressure_force"])[1]; 
-					force.z() = ((*it)["pressure_force"])[2]; 
+				if ((*it)["pressure_force"].is_array()) {
+					force.x() = ((*it)["pressure_force"])[0];
+					force.y() = ((*it)["pressure_force"])[1];
+					force.z() = ((*it)["pressure_force"])[2];
 
 					ipressure_box.pressureForce = force;
 				}
@@ -1309,9 +1322,73 @@ vector<Loads::PressureBoundaryForceBox> Input::getPressureBoundaryForceBox() {
 
 		return pressure_box_list;
 	}
-	catch(...)
+	catch (...)
 	{
 		Warning::printMessage("Error in pressure force box definition");
 		throw;
 	}
-}
+};
+
+SeismicData Input::readSeismicData(const std::string& filename, bool hasHeader = false) {
+
+	SeismicData data;
+
+	std::ifstream file(filename);
+
+	if (!file.is_open()) {
+		std::cerr << "Error during opening seismic data: " << filename << std::endl;
+		return data;
+	}
+
+	std::string line;
+
+	// ignore headers if we have ones
+	if (hasHeader && std::getline(file, line)) {
+		// headers manipulations if needed 
+	}
+
+	while (std::getline(file, line)) {
+
+		std::stringstream ss(line);
+		std::string item;
+		double t;
+		Eigen::Vector3d acc(0.0, 0.0, 0.0);
+		Eigen::Vector3d vel(0.0, 0.0, 0.0);
+
+		// read time
+		if (!std::getline(ss, item, ',')) continue;
+		t = std::stod(item);
+
+		// ax
+		if (!std::getline(ss, item, ',')) continue;
+		acc.x() = std::stod(item);
+
+		// ay
+		if (!std::getline(ss, item, ',')) continue;
+		acc.y() = std::stod(item);
+
+		// az
+		if (!std::getline(ss, item, ',')) continue;
+		acc.z() = std::stod(item);
+
+		// vx
+		if (!std::getline(ss, item, ',')) continue;
+		vel.x() = std::stod(item);
+
+		// vy
+		if (!std::getline(ss, item, ',')) continue;
+		vel.y() = std::stod(item);
+
+		// vz
+		if (!std::getline(ss, item, ',')) continue;
+		vel.z() = std::stod(item);
+
+		data.time.push_back(t);
+		data.acceleration.push_back(acc);
+		data.velocity.push_back(vel);
+	}
+
+	file.close();
+
+	return data;
+}

src/Interpolation.cpp

@@ -688,3 +688,26 @@ void Interpolation::particleDeformationGradient(Mesh* mesh, vector<Body*>* bodie
 		}
 	}
 }
+
+// Funci�n para interpolar valores de tipo Eigen::Vector3d en el tiempo itime
+Eigen::Vector3d Interpolation::interpolateVector(const std::vector<double>& times, const std::vector<Eigen::Vector3d>& values, double itime) {
+
+	if (itime <= times.front()) return values.front();
+	if (itime >= times.back()) return values.back();
+
+	// Encontrar el indice del tiempo inmediatamente superior a itime
+	auto upper = std::upper_bound(times.begin(), times.end(), itime);
+	size_t idx = std::distance(times.begin(), upper) - 1;
+
+	// Valores adyacentes para la interpolaci�n
+	double t0 = times[idx];
+	double t1 = times[idx + 1];
+	Eigen::Vector3d v0 = values[idx];
+	Eigen::Vector3d v1 = values[idx + 1];
+
+	// Interpolacion lineal para cada componente
+	double factor = (itime - t0) / (t1 - t0);
+	Eigen::Vector3d interpolatedValue = v0 + factor * (v1 - v0);
+
+	return interpolatedValue;
+}