diff --git a/.gitignore b/.gitignore index 575d7ba..acc035f 100644 --- a/.gitignore +++ b/.gitignore @@ -26,4 +26,5 @@ MPM-Geomechanics /make/VS/MPM-Geomechanics.vcxproj.user time-energy.csv simulation-data.csv -particle_stress_data.json \ No newline at end of file +particle_stress_data.json +base_acceleration.csv \ No newline at end of file diff --git a/doxygen/Doxyfile b/doxygen/Doxyfile index 8a0b00f..b675d3e 100644 --- a/doxygen/Doxyfile +++ b/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 diff --git a/doxygen/MPM.png b/doxygen/MPM.png new file mode 100644 index 0000000..3e2b370 Binary files /dev/null and b/doxygen/MPM.png differ diff --git a/inc/Boundary.h b/inc/Boundary.h index ba7426f..5d221a9 100644 --- a/inc/Boundary.h +++ b/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 diff --git a/inc/Input.h b/inc/Input.h index 6ce6a88..c40e3d7 100644 --- a/inc/Input.h +++ b/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 getPressureBoundaryForceBox(); + + Loads::SeismicData readSeismicData(const std::string& filename, bool hasHeader); }; #endif /* INPUT_H_ */ diff --git a/inc/Interpolation.h b/inc/Interpolation.h index 8a78d0a..b825561 100644 --- a/inc/Interpolation.h +++ b/inc/Interpolation.h @@ -8,6 +8,9 @@ #ifndef INTERPOLATION_H_ #define INTERPOLATION_H_ +#include "Eigen/Core" +using Eigen::Vector3d; + #include 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* bodies, double time_step); + + Eigen::Vector3d interpolateVector(const std::vector& times, const std::vector& values, double itime); }; -#endif /* INTERPOLATION_H_ */ \ No newline at end of file +#endif /* INTERPOLATION_H_ */ + diff --git a/inc/Loads.h b/inc/Loads.h index 4a5f9ae..537a540 100644 --- a/inc/Loads.h +++ b/inc/Loads.h @@ -13,6 +13,9 @@ #include 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 time; + std::vector acceleration; + std::vector velocity; + }; + + Loads::SeismicData& getSeismicData(); + + /// \brief Configures the gravity load in particles /// \param[in] bodies A list containing Body pointers void setGravity(vector& bodies); @@ -115,6 +127,9 @@ namespace Loads { /// \brief Set initial velocity in bodies /// \param[in] bodies A list containing Body pointers void setInitialVelocity(vector& bodies); + + void setSeismicData(); + SeismicData& Loads::getSeismicData(); }; #endif /* LOADS_H_ */ diff --git a/inc/MPM.h b/inc/MPM.h index 33283f5..620403f 100644 --- a/inc/MPM.h +++ b/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_ */ \ No newline at end of file diff --git a/inc/Model.h b/inc/Model.h index c1dca70..81290ff 100644 --- a/inc/Model.h +++ b/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_ */ \ No newline at end of file diff --git a/inc/Update.h b/inc/Update.h index 3b13c04..67486fa 100644 --- a/inc/Update.h +++ b/inc/Update.h @@ -10,6 +10,7 @@ #include "Mesh/Mesh.h" #include "Body/Body.h" +#include /// \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* 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* nodes, unsigned direction); + }; #endif /* UPDATE_H_ */ diff --git a/src/Input.cpp b/src/Input.cpp index 3b30cae..9692a01 100644 --- a/src/Input.cpp +++ b/src/Input.cpp @@ -31,6 +31,14 @@ using std::ifstream; using std::string; using std::to_string; +#include +#include +#include +#include +#include +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& 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 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 Input::getInitialPressureMaterial() { vector Input::getPressureBoundaryForceBox() { - try{ + try { vector 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 Input::getPressureBoundaryForceBox() { return pressure_box_list; } - catch(...) + catch (...) { Warning::printMessage("Error in pressure force box definition"); throw; } -} \ No newline at end of file +}; + +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; +} diff --git a/src/Interpolation.cpp b/src/Interpolation.cpp index c73ef4d..6856004 100644 --- a/src/Interpolation.cpp +++ b/src/Interpolation.cpp @@ -688,3 +688,26 @@ void Interpolation::particleDeformationGradient(Mesh* mesh, vector* bodie } } } + +// Funci�n para interpolar valores de tipo Eigen::Vector3d en el tiempo itime +Eigen::Vector3d Interpolation::interpolateVector(const std::vector& times, const std::vector& 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; +} diff --git a/src/Loads.cpp b/src/Loads.cpp index 244943d..ae59e38 100644 --- a/src/Loads.cpp +++ b/src/Loads.cpp @@ -8,13 +8,19 @@ #include "Loads.h" #include "Model.h" #include "Geometry.h" +#include "Input.h" namespace Loads { // list of particles with prescribed pore pressure vector prescribedPorePressureParticlesList; + + // seismic data + SeismicData seismicRecord; } +Loads::SeismicData& Loads::getSeismicData(){ return seismicRecord ;} + void Loads::setGravity(vector& bodies) { // return if gravity is not activated @@ -230,4 +236,12 @@ void Loads::setInitialVelocity(vector& bodies) { particle->setVelocity(body->getInitialVelocity()); } } +} + +void Loads::setSeismicData() +{ + if(!ModelSetup::getSeismicAnalysis()) return; + + // set seismic data + seismicRecord = Input::readSeismicData(ModelSetup::getSeismicFileName(),false); } \ No newline at end of file diff --git a/src/MPM.cpp b/src/MPM.cpp index ecef567..ace133a 100644 --- a/src/MPM.cpp +++ b/src/MPM.cpp @@ -143,7 +143,6 @@ void MPM::setupMesh() { mesh.setBoundaryRestrictions(Input::getMeshBoundaryConditions()); if (ModelSetup::getTwoPhaseActive()){ - // configure the mesh boundary conditions of fluid mesh.setBoundaryRestrictionsFluid(Input::getMeshBoundaryConditionsFluid()); } @@ -236,6 +235,8 @@ void MPM::setupLoads() { // set initial velocity Loads::setInitialVelocity(bodies); + // set seismic acceleration + Loads::setSeismicData(); } void MPM::setupDamping() { diff --git a/src/Model.cpp b/src/Model.cpp index 176c680..0ebb653 100644 --- a/src/Model.cpp +++ b/src/Model.cpp @@ -51,9 +51,10 @@ namespace ModelSetup { unsigned resultNumber=10; //!< number of results to write // time - double dt=0.0; //!< time step - double time=0.0; //!< simulation 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 currentTime = 0.0; //!< current simulation time // initial time simulation std::chrono::system_clock::time_point initialSimulationTime; @@ -73,6 +74,9 @@ namespace ModelSetup { // default interpolation functions InterpolationFunctionType interpolationType = InterpolationFunctionType::GIMP; //!< interpolation function type + // seismic analysis + bool seismicAnalysisActive = false; //!< is seismic analysis + /// /// Function members /// @@ -120,13 +124,15 @@ namespace ModelSetup { setGravityActive((gravity.x()!=0.0||gravity.y()!=0.0||gravity.z()!=0.0)?true:false); } - // coupled analisys + // coupled analysis 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; } + double getCurrentTime() { return currentTime; } + void setCurrentTime(double a) { currentTime = a; } // axisymmetric analisys bool getAxisymetricActive() { return axisymetricActive; } @@ -161,4 +167,10 @@ namespace ModelSetup { omp_set_num_threads((nThreads>0&&nThreads<=(unsigned)omp_get_num_procs())?nThreads:omp_get_num_procs()); #endif } + + // Seismic analysis + string seismic_file_name = "base_acceleration.csv"; + string getSeismicFileName() {return seismic_file_name;} + bool getSeismicAnalysis() {return seismicAnalysisActive;} + void setSeismicAnalysis(bool a) {seismicAnalysisActive = a;} } diff --git a/src/Output.cpp b/src/Output.cpp index 9f94262..b1e4f42 100644 --- a/src/Output.cpp +++ b/src/Output.cpp @@ -243,6 +243,16 @@ namespace Output{ partFile<<"\n"; } + if (isFieldRequired("velocity")) { + + // particle velocity + partFile << "\n"; + for (int i = 0; i < nPoints; ++i) { + partFile << scientific << particles->at(i)->getVelocity() << "\n"; + } + partFile << "\n"; + } + if (isFieldRequired("stress")){ // particle stress xx diff --git a/src/SolverExplicitTwoPhaseUSL.cpp b/src/SolverExplicitTwoPhaseUSL.cpp index 462cdb5..7e81a8b 100644 --- a/src/SolverExplicitTwoPhaseUSL.cpp +++ b/src/SolverExplicitTwoPhaseUSL.cpp @@ -163,7 +163,7 @@ void SolverExplicitTwoPhaseUSL::Solve( ){ DynamicRelaxation::setStaticSolution(bodies,loopCounter); // advance in time - iTime+=dt; + ModelSetup::setCurrentTime(iTime += dt); } // write the Eulerian mesh diff --git a/src/SolverExplicitUSL.cpp b/src/SolverExplicitUSL.cpp index 2e3ffae..2385ba4 100644 --- a/src/SolverExplicitUSL.cpp +++ b/src/SolverExplicitUSL.cpp @@ -19,37 +19,39 @@ using std::cout; #include "Output.h" #include "DynamicRelaxation.h" -SolverExplicitUSL::SolverExplicitUSL():Solver() { } +SolverExplicitUSL::SolverExplicitUSL() : Solver() {} -void SolverExplicitUSL::Solve( ){ +void SolverExplicitUSL::Solve() +{ // initialize simulation variables - double time=ModelSetup::getTime(); - double dt=ModelSetup::getTimeStep(); - int resultSteps=ModelSetup::getResultSteps(); - double iTime=0.0; - int loopCounter=0; + double time = ModelSetup::getTime(); + double dt = ModelSetup::getTimeStep(); + int resultSteps = ModelSetup::getResultSteps(); + double iTime = 0.0; + int loopCounter = 0; // solve in time - while(iTime* bodies, double dt) { particles->at(i)->updatePressure(dt); } } - // set up prescribed pore pressure Loads::updatePrescribedPorePressure(bodies); } @@ -295,6 +295,8 @@ void Update::particlePosition(Mesh* mesh, vector* bodies, double dt) { void Update::setPlaneMomentum(const Boundary::planeBoundary* plane, vector* nodes, unsigned dir) { + Eigen::Vector3d interpolatedVelocity = ModelSetup::getSeismicAnalysis() ? Interpolation::interpolateVector(Loads::getSeismicData().time, Loads::getSeismicData().velocity, ModelSetup::getCurrentTime()) : Vector3d{0,0,0} ; + // for each boundary node #pragma omp parallel for shared(plane, nodes, dir) for (int i = 0; i < plane->nodes.size(); ++i){ @@ -305,47 +307,59 @@ void Update::setPlaneMomentum(const Boundary::planeBoundary* plane, vectorgetActive()) { // witch type of restriction - switch(plane->restriction) { - + switch(plane->restriction) + { // free condition case Boundary::BoundaryType::FREE: + { break; - + } // fixed condition case Boundary::BoundaryType::FIXED: - + { // set all momentum components as zero nodeI->setMomentum(Vector3d::Zero()); break; - + } // sliding restriction case Boundary::BoundaryType::SLIDING: - + { // get current boundary nodal momentum Vector3d momentum = nodeI->getMomentum(); - - // witch direction of the normal vector - switch(dir) { + // witch direction of the normal vector + switch (dir) + { // normal pointed to x - case Update::Direction::X : - momentum.x()=0.0; + case Update::Direction::X: + { + momentum.x() = 0.0; break; - + } // normal pointed to y - case Update::Direction::Y : - momentum.y()=0.0; + case Update::Direction::Y: + { + momentum.y() = 0.0; break; - + } // normal pointed to z - case Update::Direction::Z : - momentum.z()=0.0; + case Update::Direction::Z: + { + momentum.z() = 0.0; break; + } } // set the boundary nodal momentum nodeI->setMomentum(momentum); break; + } + // earthquake boundary condition + case Boundary::BoundaryType::EARTHQUAKE: + { + nodeI->setMomentum(nodeI->getMass() * interpolatedVelocity); + break; + } } } } @@ -378,7 +392,7 @@ void Update::setPlaneMomentumFluid(const Boundary::planeBoundary* plane, vector< // perpendicular restriction case Boundary::BoundaryType::SLIDING: - + { // get current boundary nodal momentum of fluid Vector3d momentum = *(nodeI->getMomentumFluid()); @@ -404,6 +418,7 @@ void Update::setPlaneMomentumFluid(const Boundary::planeBoundary* plane, vector< // set the boundary nodal momentum nodeI->setMomentumFluid(momentum); break; + } } } } @@ -414,22 +429,14 @@ void Update::boundaryConditionsMomentumFluid(Mesh* mesh) { // get nodes vector* nodes = mesh->getNodes(); - // plane X0 is the plane passing for the origin and points to -x axes - setPlaneMomentumFluid(mesh->getBoundary()->getPlaneX0(), nodes, Update::Direction::X); + // set p=0 in fixed nodes - // plane Y0 is the plane passing for the origin and points to -y axes + setPlaneMomentumFluid(mesh->getBoundary()->getPlaneX0(), nodes, Update::Direction::X); setPlaneMomentumFluid(mesh->getBoundary()->getPlaneY0(), nodes, Update::Direction::Y); - - // plane Z0 is the plane passing for the origin and points to -z axes setPlaneMomentumFluid(mesh->getBoundary()->getPlaneZ0(), nodes, Update::Direction::Z); - // plane Xn is the plane passing for the maximum x coordinate and points to x axes setPlaneMomentumFluid(mesh->getBoundary()->getPlaneXn(), nodes, Update::Direction::X); - - // plane Yn is the plane passing for the maximum y coordinate and points to y axes setPlaneMomentumFluid(mesh->getBoundary()->getPlaneYn(), nodes, Update::Direction::Y); - - // plane Zn is the plane passing for the maximum z coordinate and points to z axes setPlaneMomentumFluid(mesh->getBoundary()->getPlaneZn(), nodes, Update::Direction::Z); } @@ -438,27 +445,22 @@ void Update::boundaryConditionsMomentum(Mesh* mesh) { // get nodes vector* nodes = mesh->getNodes(); - // plane X0 is the plane passing for the origin and points to -x axes + // set p = 0 in fixed direction + setPlaneMomentum(mesh->getBoundary()->getPlaneX0(), nodes, Update::Direction::X); - - // plane Y0 is the plane passing for the origin and points to -y axes setPlaneMomentum(mesh->getBoundary()->getPlaneY0(), nodes, Update::Direction::Y); - - // plane Z0 is the plane passing for the origin and points to -z axes setPlaneMomentum(mesh->getBoundary()->getPlaneZ0(), nodes, Update::Direction::Z); - - // plane Xn is the plane passing for the maximum x coordinate and points to x axes - setPlaneMomentum(mesh->getBoundary()->getPlaneXn(), nodes, Update::Direction::X); - // plane Yn is the plane passing for the maximum y coordinate and points to y axes + setPlaneMomentum(mesh->getBoundary()->getPlaneXn(), nodes, Update::Direction::X); setPlaneMomentum(mesh->getBoundary()->getPlaneYn(), nodes, Update::Direction::Y); - - // plane Zn is the plane passing for the maximum z coordinate and points to z axes setPlaneMomentum(mesh->getBoundary()->getPlaneZn(), nodes, Update::Direction::Z); -} + +} void Update::setPlaneForce(const Boundary::planeBoundary* plane, vector* nodes, unsigned dir) { + Eigen::Vector3d interpolatedAcceleration = ModelSetup::getSeismicAnalysis() ? Interpolation::interpolateVector(Loads::getSeismicData().time, Loads::getSeismicData().acceleration, ModelSetup::getCurrentTime()) : Vector3d{ 0,0,0 }; + // get boundary nodes #pragma omp parallel for shared(plane, nodes, dir) for (int i = 0; i < plane->nodes.size(); ++i) { @@ -469,74 +471,79 @@ void Update::setPlaneForce(const Boundary::planeBoundary* plane, vector* if (nodeI->getActive()) { // witch type of restriction - switch(plane->restriction) { - + switch(plane->restriction) + { // free condition case Boundary::BoundaryType::FREE: + { break; - + } // fixed condition case Boundary::BoundaryType::FIXED: - + { // set all force component as zero nodeI->setTotalForce(Vector3d::Zero()); break; - + } // perpendicular restriction case Boundary::BoundaryType::SLIDING: - + { // get current boundary nodal force Vector3d force = nodeI->getTotalForce(); // witch direction of the normal vector - switch(dir) { - + switch(dir) + { // normal pointed to x case Update::Direction::X : + { force.x()=0.0; break; - + } // normal pointed to y case Update::Direction::Y : + { force.y()=0.0; break; - + } // normal pointed to z case Update::Direction::Z : + { force.z()=0.0; break; + } } // set boundary nodal force nodeI->setTotalForce(force); break; + } + // earthquake boundary condition + case Boundary::BoundaryType::EARTHQUAKE: + { + nodeI->setTotalForce(nodeI->getMass() * interpolatedAcceleration); + break; + } } } } } void Update::boundaryConditionsForce(Mesh* mesh) { - + // get nodes vector* nodes = mesh->getNodes(); - // plane X0 is the plane passing for the origin and points to -x axes - setPlaneForce(mesh->getBoundary()->getPlaneX0(), nodes, Update::Direction::X); + // set f = 0 in fixed direction - // plane Y0 is the plane passing for the origin and points to -y axes + setPlaneForce(mesh->getBoundary()->getPlaneX0(), nodes, Update::Direction::X); setPlaneForce(mesh->getBoundary()->getPlaneY0(), nodes, Update::Direction::Y); - - // plane Z0 is the plane passing for the origin and points to -z axes setPlaneForce(mesh->getBoundary()->getPlaneZ0(), nodes, Update::Direction::Z); - - // plane Xn is the plane passing for the maximum x coordinate and points to x axes - setPlaneForce(mesh->getBoundary()->getPlaneXn(), nodes, Update::Direction::X); - // plane Yn is the plane passing for the maximum x coordinate and points to y axes + setPlaneForce(mesh->getBoundary()->getPlaneXn(), nodes, Update::Direction::X); setPlaneForce(mesh->getBoundary()->getPlaneYn(), nodes, Update::Direction::Y); - - // plane Zn is the plane passing for the maximum x coordinate and points to z axes setPlaneForce(mesh->getBoundary()->getPlaneZn(), nodes, Update::Direction::Z); + } void Update::setPlaneForceFluid(const Boundary::planeBoundary* plane, vector* nodes, unsigned dir) { @@ -601,23 +608,14 @@ void Update::boundaryConditionsForceFluid(Mesh* mesh) { // get nodes vector* nodes = mesh->getNodes(); - - // plane X0 is the plane passing for the origin and points to -x axes + + // set f=0 in fixed nodes setPlaneForceFluid(mesh->getBoundary()->getPlaneX0(), nodes, Update::Direction::X); - - // plane Y0 is the plane passing for the origin and points to -y axes setPlaneForceFluid(mesh->getBoundary()->getPlaneY0(), nodes, Update::Direction::Y); - - // plane Z0 is the plane passing for the origin and points to -z axes setPlaneForceFluid(mesh->getBoundary()->getPlaneZ0(), nodes, Update::Direction::Z); - // plane Xn is the plane passing for the maximum x coordinate and points to x axes setPlaneForceFluid(mesh->getBoundary()->getPlaneXn(), nodes, Update::Direction::X); - - // plane Yn is the plane passing for the maximum x coordinate and points to y axes setPlaneForceFluid(mesh->getBoundary()->getPlaneYn(), nodes, Update::Direction::Y); - - // plane Zn is the plane passing for the maximum x coordinate and points to z axes setPlaneForceFluid(mesh->getBoundary()->getPlaneZn(), nodes, Update::Direction::Z); } @@ -640,4 +638,4 @@ void Update::contributionNodes(Mesh* mesh, vector* bodies) { particles->at(i)->updateContributionNodes(mesh); } } -} +} \ No newline at end of file diff --git a/tests/seismic-load/seismic-column.json b/tests/seismic-load/seismic-column.json new file mode 100644 index 0000000..5593a81 --- /dev/null +++ b/tests/seismic-load/seismic-column.json @@ -0,0 +1,69 @@ +{ + "stress_scheme_update":"USL", + + "shape_function":"GIMP", + + "time":10, + + "time_step":0.0005, + + "gravity":[0.0,0.0,0.0], + + "n_threads":1, + + "damping": + { + "type":"local", + "value":0.0 + }, + + "results": + { + "print":100, + "fields":["id","displacement","velocity","material","active","body"] + }, + + "n_phases":1, + + "mesh": + { + "cells_dimension":[0.1,0.1,0.1], + "cells_number":[10,10,15], + "origin":[0,0,0], + + "boundary_conditions": + { + "plane_X0":"sliding", + "plane_Y0":"sliding", + "plane_Z0":"earthquake", + "plane_Xn":"sliding", + "plane_Yn":"sliding", + "plane_Zn":"sliding" + } + }, + + "material": + { + "elastic_1": + { + "type":"elastic", + "id":1, + "young":10e6, + "density":2500, + "poisson":0.25 + } + }, + + "body": + { + "columns_1": + { + "type":"cuboid", + "id":1, + "point_p1":[0.2,0.2,0], + "point_p2":[0.5,0.5,1.0], + "material_id":1 + } + + } + } \ No newline at end of file