[Initial condition] Add functions to initialize velocity and accelerations

include/io/io_mesh.h

@@ -70,6 +70,12 @@ class IOMesh {
   virtual std::vector<std::tuple<mpm::Index, double>>
       read_particles_scalar_properties(const std::string& scalar_file) = 0;
 
+  //! Read particle vector properties
+  //! \param[in] vector_file file name with particle vector properties
+  //! \retval Vector of particles vector properties
+  virtual std::vector<std::tuple<mpm::Index, Eigen::Matrix<double, Tdim, 1>>>
+      read_particles_vector_properties(const std::string& vector_file) = 0;
+
   //! Read pressure constraints file
   //! \param[in] pressure_constraints_files file name with pressure constraints
   virtual std::vector<std::tuple<mpm::Index, double>> read_pressure_constraints(

include/io/io_mesh_ascii.h

@@ -57,6 +57,12 @@ class IOMeshAscii : public IOMesh<Tdim> {
   std::vector<std::tuple<mpm::Index, double>> read_particles_scalar_properties(
       const std::string& scalar_file) override;
 
+  //! Read particle vector properties
+  //! \param[in] vector_file file name with particle vector properties
+  //! \retval Vector of particles vector properties
+  std::vector<std::tuple<mpm::Index, Eigen::Matrix<double, Tdim, 1>>>
+      read_particles_vector_properties(const std::string& vector_file) override;
+
   //! Read pressure constraints file
   //! \param[in] pressure_constraints_files file name with pressure
   //! constraints

include/io/io_mesh_ascii.tcc

@@ -300,6 +300,55 @@ std::vector<std::tuple<mpm::Index, double>>
   return scalar_properties;
 }
 
+//! Return particles vector properties
+template <unsigned Tdim>
+std::vector<std::tuple<mpm::Index, Eigen::Matrix<double, Tdim, 1>>>
+    mpm::IOMeshAscii<Tdim>::read_particles_vector_properties(
+        const std::string& vector_file) {
+
+  // Particles vector properties
+  std::vector<std::tuple<mpm::Index, VectorDim>> vector_properties;
+
+  // input file stream
+  std::fstream file;
+  file.open(vector_file.c_str(), std::ios::in);
+
+  try {
+    if (file.is_open() && file.good()) {
+      // Line
+      std::string line;
+      while (std::getline(file, line)) {
+        boost::algorithm::trim(line);
+        std::istringstream istream(line);
+        // ignore comment lines (# or !) or blank lines
+        if ((line.find('#') == std::string::npos) &&
+            (line.find('!') == std::string::npos) && (line != "")) {
+          // ID
+          mpm::Index id;
+          // Vector
+          Eigen::Matrix<double, Tdim, 1> vector;
+          while (istream.good()) {
+            // Read stream
+            istream >> id;
+            // Read to vector
+            for (unsigned i = 0; i < Tdim; ++i) istream >> vector[i];
+
+            vector_properties.emplace_back(std::make_tuple(id, vector));
+          }
+        }
+      }
+    } else {
+      throw std::runtime_error("File not open or not good!");
+    }
+    file.close();
+  } catch (std::exception& exception) {
+    console_->error("Read particle {} #{}: {}\n", __FILE__, __LINE__,
+                    exception.what());
+    file.close();
+  }
+  return vector_properties;
+}
+
 //! Read pressure constraints file
 template <unsigned Tdim>
 std::vector<std::tuple<mpm::Index, double>>

include/mesh/mesh.h

@@ -317,10 +317,22 @@ class Mesh {
       const std::map<mpm::Index, Eigen::Matrix<double, Tdim, 1>>& euler_angles);
 
   //! Assign particles volumes
-  //! \param[in] particle_volumes Volume at dir on particle
+  //! \param[in] particle_volumes Volume on particle
   bool assign_particles_volumes(
       const std::vector<std::tuple<mpm::Index, double>>& particle_volumes);
 
+  //! Assign particles velocities
+  //! \param[in] particle_velocities Initial particle velocity
+  bool assign_particles_velocities(
+      const std::vector<std::tuple<mpm::Index, Eigen::Matrix<double, Tdim, 1>>>&
+          particle_velocities);
+
+  //! Assign particles accelerations
+  //! \param[in] particle_accelerations Initial particle velocity
+  bool assign_particles_accelerations(
+      const std::vector<std::tuple<mpm::Index, Eigen::Matrix<double, Tdim, 1>>>&
+          particle_accelerations);
+
   //! Create particles tractions
   //! \param[in] mfunction Math function if defined
   //! \param[in] setid Particle set id

include/mesh/mesh.tcc

@@ -1313,6 +1313,78 @@ bool mpm::Mesh<Tdim>::assign_particles_volumes(
   return status;
 }
 
+//! Assign particle velocities
+template <unsigned Tdim>
+bool mpm::Mesh<Tdim>::assign_particles_velocities(
+    const std::vector<std::tuple<mpm::Index, Eigen::Matrix<double, Tdim, 1>>>&
+        particle_velocities) {
+  bool status = true;
+
+  try {
+    if (!particles_.size())
+      throw std::runtime_error(
+          "No particles have been assigned in mesh, cannot assign particles "
+          "velocities");
+
+    if (particles_.size() < particle_velocities.size())
+      throw std::runtime_error(
+          "Number of particles in mesh and initial velocities don't match");
+
+    // Loop over particle velocities
+    for (const auto& particle_vel : particle_velocities) {
+      // Particle id
+      mpm::Index pid = std::get<0>(particle_vel);
+      // Velocity vector
+      VectorDim pvel = std::get<1>(particle_vel);
+
+      if (map_particles_.find(pid) != map_particles_.end()) {
+        map_particles_[pid]->assign_velocity(pvel);
+      }
+    }
+
+  } catch (std::exception& exception) {
+    console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what());
+    status = false;
+  }
+  return status;
+}
+
+//! Assign particle accelerations
+template <unsigned Tdim>
+bool mpm::Mesh<Tdim>::assign_particles_accelerations(
+    const std::vector<std::tuple<mpm::Index, Eigen::Matrix<double, Tdim, 1>>>&
+        particle_accelerations) {
+  bool status = true;
+
+  try {
+    if (!particles_.size())
+      throw std::runtime_error(
+          "No particles have been assigned in mesh, cannot assign particles "
+          "accelerations");
+
+    if (particles_.size() < particle_accelerations.size())
+      throw std::runtime_error(
+          "Number of particles in mesh and initial accelerations don't match");
+
+    // Loop over particle accelerations
+    for (const auto& particle_acc : particle_accelerations) {
+      // Particle id
+      mpm::Index pid = std::get<0>(particle_acc);
+      // Acceleration vector
+      VectorDim pacc = std::get<1>(particle_acc);
+
+      if (map_particles_.find(pid) != map_particles_.end()) {
+        map_particles_[pid]->assign_acceleration(pacc);
+      }
+    }
+
+  } catch (std::exception& exception) {
+    console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what());
+    status = false;
+  }
+  return status;
+}
+
 //! Compute and assign rotation matrix to nodes
 template <unsigned Tdim>
 bool mpm::Mesh<Tdim>::compute_nodal_rotation_matrices(

include/solvers/mpm_base.h

@@ -203,6 +203,20 @@ class MPMBase : public MPM {
       const Json& mesh_prop,
       const std::shared_ptr<mpm::IOMesh<Tdim>>& particle_io);
 
+  //! Particles velocities
+  //! \param[in] mesh_prop Mesh properties
+  //! \param[in] particle_io Particle IO handle
+  void particles_velocities(
+      const Json& mesh_prop,
+      const std::shared_ptr<mpm::IOMesh<Tdim>>& particle_io);
+
+  //! Particles accelerations
+  //! \param[in] mesh_prop Mesh properties
+  //! \param[in] particle_io Particle IO handle
+  void particles_accelerations(
+      const Json& mesh_prop,
+      const std::shared_ptr<mpm::IOMesh<Tdim>>& particle_io);
+
   // Particles pore pressures
   //! \param[in] mesh_prop Mesh properties
   //! \param[in] particle_io Particle IO handle

include/solvers/mpm_base.tcc

@@ -425,6 +425,12 @@ void mpm::MPMBase<Tdim>::initialise_particles() {
   // Read and assign particles stresses
   this->particles_stresses(mesh_props, particle_io);
 
+  // Read and assign particles initial velocity
+  this->particles_velocities(mesh_props, particle_io);
+
+  // Read and assign particles initial acceleration
+  this->particles_accelerations(mesh_props, particle_io);
+
   // Read and assign particles initial pore pressure
   this->particles_pore_pressures(mesh_props, particle_io);
 
@@ -1444,6 +1450,7 @@ void mpm::MPMBase<Tdim>::particles_volumes(
     if (mesh_props.find("particles_volumes") != mesh_props.end()) {
       std::string fparticles_volumes =
           mesh_props["particles_volumes"].template get<std::string>();
+
       if (!io_->file_name(fparticles_volumes).empty()) {
         bool particles_volumes =
             mesh_->assign_particles_volumes(particle_io->read_particles_volumes(
@@ -1454,12 +1461,117 @@ void mpm::MPMBase<Tdim>::particles_volumes(
       }
     } else
       throw std::runtime_error("Particle volumes JSON data not found");
+
   } catch (std::exception& exception) {
     console_->warn("#{}: Particle volumes are undefined; {}", __LINE__,
                    exception.what());
   }
 }
 
+// Read and assign particle velocities
+template <unsigned Tdim>
+void mpm::MPMBase<Tdim>::particles_velocities(
+    const Json& mesh_props,
+    const std::shared_ptr<mpm::IOMesh<Tdim>>& particle_io) {
+  try {
+    // Read particle initial velocities
+    if (mesh_props.find("particles_velocities") != mesh_props.end()) {
+      // Get generator type
+      const std::string type = mesh_props["particles_velocities"]["type"]
+                                   .template get<std::string>();
+      if (type == "file") {
+        std::string fparticle_vel =
+            mesh_props["particles_velocities"]["location"]
+                .template get<std::string>();
+        if (!io_->file_name(fparticle_vel).empty()) {
+
+          // Get particle initial velocities
+          const auto particles_vel =
+              particle_io->read_particles_vector_properties(
+                  io_->file_name(fparticle_vel));
+
+          // Assign particle velocities
+          if (!mesh_->assign_particles_velocities(particles_vel))
+            throw std::runtime_error(
+                "Particles velocities are not properly assigned");
+        }
+      } else if (type == "isotropic") {
+        Eigen::Matrix<double, Tdim, 1> in_vel;
+        in_vel.setZero();
+        if (mesh_props["particles_velocities"]["values"].is_array() &&
+            mesh_props["particles_velocities"]["values"].size() ==
+                in_vel.size()) {
+          for (unsigned i = 0; i < in_vel.size(); ++i) {
+            in_vel[i] = mesh_props["particles_velocities"]["values"].at(i);
+          }
+          mesh_->iterate_over_particles(
+              std::bind(&mpm::ParticleBase<Tdim>::assign_velocity,
+                        std::placeholders::_1, in_vel));
+        } else {
+          throw std::runtime_error("Initial velocity dimension is invalid");
+        }
+      }
+    } else
+      throw std::runtime_error("Particle velocities JSON data not found");
+
+  } catch (std::exception& exception) {
+    console_->warn("#{}: Particle velocities are undefined {} ", __LINE__,
+                   exception.what());
+  }
+}
+
+// Read and assign particle accelerations
+template <unsigned Tdim>
+void mpm::MPMBase<Tdim>::particles_accelerations(
+    const Json& mesh_props,
+    const std::shared_ptr<mpm::IOMesh<Tdim>>& particle_io) {
+  try {
+    // Read particle initial accelerations
+    if (mesh_props.find("particles_accelerations") != mesh_props.end()) {
+      // Get generator type
+      const std::string type = mesh_props["particles_accelerations"]["type"]
+                                   .template get<std::string>();
+      if (type == "file") {
+        std::string fparticle_acc =
+            mesh_props["particles_accelerations"]["location"]
+                .template get<std::string>();
+        if (!io_->file_name(fparticle_acc).empty()) {
+
+          // Get particle initial accelerations
+          const auto particles_acc =
+              particle_io->read_particles_vector_properties(
+                  io_->file_name(fparticle_acc));
+
+          // Assign particle accelerations
+          if (!mesh_->assign_particles_accelerations(particles_acc))
+            throw std::runtime_error(
+                "Particles accelerations are not properly assigned");
+        }
+      } else if (type == "isotropic") {
+        Eigen::Matrix<double, Tdim, 1> in_acc;
+        in_acc.setZero();
+        if (mesh_props["particles_accelerations"]["values"].is_array() &&
+            mesh_props["particles_accelerations"]["values"].size() ==
+                in_acc.size()) {
+          for (unsigned i = 0; i < in_acc.size(); ++i) {
+            in_acc[i] = mesh_props["particles_accelerations"]["values"].at(i);
+          }
+          mesh_->iterate_over_particles(
+              std::bind(&mpm::ParticleBase<Tdim>::assign_acceleration,
+                        std::placeholders::_1, in_acc));
+        } else {
+          throw std::runtime_error("Initial acceleration dimension is invalid");
+        }
+      }
+    } else
+      throw std::runtime_error("Particle accelerations JSON data not found");
+
+  } catch (std::exception& exception) {
+    console_->warn("#{}: Particle accelerations are undefined {} ", __LINE__,
+                   exception.what());
+  }
+}
+
 // Particle velocity constraints
 template <unsigned Tdim>
 void mpm::MPMBase<Tdim>::particle_velocity_constraints() {