Feature/multiphase material

include/generators/injection.h

@@ -7,7 +7,7 @@ struct Injection {
   // Particle type
   std::string particle_type;
   // Material id
-  unsigned material_id;
+  std::vector<unsigned> material_ids;
   // Cell id
   int cell_set_id;
   // Start

include/mesh.h

@@ -195,8 +195,8 @@ class Mesh {
   //! \retval status Create particle status
   bool create_particles(const std::string& particle_type,
                         const std::vector<VectorDim>& coordinates,
-                        unsigned material_id, unsigned pset_id,
-                        bool check_duplicates = true);
+                        const std::vector<unsigned>& material_ids,
+                        unsigned pset_id, bool check_duplicates = true);
 
   //! Add a particle to the mesh
   //! \param[in] particle A shared pointer to particle
@@ -341,8 +341,8 @@ class Mesh {
   //! \retval point Material point coordinates
   bool generate_material_points(unsigned nquadratures,
                                 const std::string& particle_type,
-                                unsigned material_id, int cset_id,
-                                unsigned pset_id);
+                                const std::vector<unsigned>& material_ids,
+                                int cset_id, unsigned pset_id);
 
   //! Initialise material models
   //! \param[in] materials Material models

include/mesh.tcc

@@ -465,10 +465,9 @@ mpm::Index mpm::Mesh<Tdim>::nnodes_rank() {
 
 //! Create cells from node lists
 template <unsigned Tdim>
-bool mpm::Mesh<Tdim>::generate_material_points(unsigned nquadratures,
-                                               const std::string& particle_type,
-                                               unsigned material_id,
-                                               int cset_id, unsigned pset_id) {
+bool mpm::Mesh<Tdim>::generate_material_points(
+    unsigned nquadratures, const std::string& particle_type,
+    const std::vector<unsigned>& material_ids, int cset_id, unsigned pset_id) {
   bool status = true;
   try {
     if (cells_.size() > 0) {
@@ -477,7 +476,9 @@ bool mpm::Mesh<Tdim>::generate_material_points(unsigned nquadratures,
       unsigned before_generation = this->nparticles();
       bool checks = false;
       // Get material
-      auto material = materials_.at(material_id);
+      std::vector<std::shared_ptr<mpm::Material<Tdim>>> materials;
+      for (auto m_id : material_ids)
+        materials.emplace_back(materials_.at(m_id));
 
       // If set id is -1, use all cells
       auto cset = (cset_id == -1) ? this->cells_ : cell_sets_.at(cset_id);
@@ -501,7 +502,8 @@ bool mpm::Mesh<Tdim>::generate_material_points(unsigned nquadratures,
           status = this->add_particle(particle, checks);
           if (status) {
             map_particles_[pid]->assign_cell(*citr);
-            map_particles_[pid]->assign_material(material);
+            for (unsigned phase = 0; phase < materials.size(); phase++)
+              map_particles_[pid]->assign_material(materials[phase], phase);
             pids.emplace_back(pid);
           } else
             throw std::runtime_error("Generate particles in mesh failed");
@@ -535,13 +537,15 @@ bool mpm::Mesh<Tdim>::generate_material_points(unsigned nquadratures,
 template <unsigned Tdim>
 bool mpm::Mesh<Tdim>::create_particles(
     const std::string& particle_type, const std::vector<VectorDim>& coordinates,
-    unsigned material_id, unsigned pset_id, bool check_duplicates) {
+    const std::vector<unsigned>& material_ids, unsigned pset_id,
+    bool check_duplicates) {
   bool status = true;
   try {
     // Particle ids
     std::vector<mpm::Index> pids;
     // Get material
-    auto material = materials_.at(material_id);
+    std::vector<std::shared_ptr<mpm::Material<Tdim>>> materials;
+    for (auto m_id : material_ids) materials.emplace_back(materials_.at(m_id));
     // Check if particle coordinates is empty
     if (coordinates.empty())
       throw std::runtime_error("List of coordinates is empty");
@@ -560,7 +564,8 @@ bool mpm::Mesh<Tdim>::create_particles(
 
       // If insertion is successful
       if (insert_status) {
-        map_particles_[pid]->assign_material(material);
+        for (unsigned phase = 0; phase < materials.size(); phase++)
+          map_particles_[pid]->assign_material(materials[phase], phase);
         pids.emplace_back(pid);
       } else
         throw std::runtime_error("Addition of particle to mesh failed!");
@@ -1644,13 +1649,19 @@ bool mpm::Mesh<Tdim>::generate_particles(const std::shared_ptr<mpm::IO>& io,
       auto particle_type =
           generator["particle_type"].template get<std::string>();
       // Material id
-      unsigned material_id = generator["material_id"].template get<unsigned>();
+      std::vector<unsigned> material_ids;
+      if (generator.at("material_id").is_array())
+        material_ids =
+            generator["material_id"].template get<std::vector<unsigned>>();
+      else
+        material_ids.emplace_back(
+            generator["material_id"].template get<unsigned>());
       // Cell set id
       int cset_id = generator["cset_id"].template get<int>();
       // Particle set id
       unsigned pset_id = generator["pset_id"].template get<unsigned>();
       status = this->generate_material_points(nparticles_dir, particle_type,
-                                              material_id, cset_id, pset_id);
+                                              material_ids, cset_id, pset_id);
     }
 
     // Generate material points at the Gauss location in all cells
@@ -1663,7 +1674,12 @@ bool mpm::Mesh<Tdim>::generate_particles(const std::shared_ptr<mpm::IO>& io,
       inject.particle_type =
           generator["particle_type"].template get<std::string>();
       // Material id
-      inject.material_id = generator["material_id"].template get<unsigned>();
+      if (generator.at("material_id").is_array())
+        inject.material_ids =
+            generator["material_id"].template get<std::vector<unsigned>>();
+      else
+        inject.material_ids.emplace_back(
+            generator["material_id"].template get<unsigned>());
       // Cell set id
       inject.cell_set_id = generator["cset_id"].template get<int>();
       // Duration of injection
@@ -1709,7 +1725,10 @@ void mpm::Mesh<Tdim>::inject_particles(double current_time) {
     unsigned pid = this->nparticles();
     bool checks = false;
     // Get material
-    auto material = materials_.at(injection.material_id);
+    std::vector<std::shared_ptr<mpm::Material<Tdim>>> materials;
+    for (auto m_id : injection.material_ids)
+      materials.emplace_back(materials_.at(m_id));
+
     // Check if duration is within the current time
     if (injection.start_time <= current_time &&
         injection.end_time > current_time) {
@@ -1742,7 +1761,8 @@ void mpm::Mesh<Tdim>::inject_particles(double current_time) {
             unsigned status = this->add_particle(particle, checks);
             if (status) {
               map_particles_[pid]->assign_cell(*citr);
-              map_particles_[pid]->assign_material(material);
+              for (unsigned phase = 0; phase < materials.size(); phase++)
+                map_particles_[pid]->assign_material(materials[phase], phase);
               ++pid;
               injected_particles.emplace_back(particle);
             }
@@ -1773,7 +1793,13 @@ bool mpm::Mesh<Tdim>::read_particles_file(const std::shared_ptr<mpm::IO>& io,
   bool check_duplicates = generator["check_duplicates"].template get<bool>();
 
   // Material id
-  unsigned material_id = generator["material_id"].template get<unsigned>();
+  std::vector<unsigned> material_ids;
+  if (generator.at("material_id").is_array())
+    material_ids =
+        generator["material_id"].template get<std::vector<unsigned>>();
+  else
+    material_ids.emplace_back(
+        generator["material_id"].template get<unsigned>());
 
   const std::string reader = generator["io_type"].template get<std::string>();
 
@@ -1784,7 +1810,7 @@ bool mpm::Mesh<Tdim>::read_particles_file(const std::shared_ptr<mpm::IO>& io,
   auto coords = particle_io->read_particles(file_loc);
 
   // Create particles from coordinates
-  bool status = this->create_particles(particle_type, coords, material_id,
+  bool status = this->create_particles(particle_type, coords, material_ids,
                                        pset_id, check_duplicates);
 
   if (!status) throw std::runtime_error("Addition of particles to mesh failed");

include/particles/particle.h

@@ -62,10 +62,12 @@ class Particle : public ParticleBase<Tdim> {
   //! Assign material history variables
   //! \param[in] state_vars State variables
   //! \param[in] material Material associated with the particle
+  //! \param[in] phase Index to indicate material phase
   //! \retval status Status of cloning HDF5 particle
   bool assign_material_state_vars(
       const mpm::dense_map& state_vars,
-      const std::shared_ptr<mpm::Material<Tdim>>& material) override;
+      const std::shared_ptr<mpm::Material<Tdim>>& material,
+      unsigned phase = mpm::ParticlePhase::Solid) override;
 
   //! Retrun particle data as HDF5
   //! \retval particle HDF5 data of the particle
@@ -157,8 +159,9 @@ class Particle : public ParticleBase<Tdim> {
 
   //! Assign material
   //! \param[in] material Pointer to a material
-  bool assign_material(
-      const std::shared_ptr<Material<Tdim>>& material) override;
+  //! \param[in] phase Index to indicate phase
+  bool assign_material(const std::shared_ptr<Material<Tdim>>& material,
+                       unsigned phase = mpm::ParticlePhase::Solid) override;
 
   //! Compute strain
   //! \param[in] dt Analysis time step
@@ -233,24 +236,31 @@ class Particle : public ParticleBase<Tdim> {
 
   //! Return a state variable
   //! \param[in] var State variable
+  //! \param[in] phase Index to indicate phase
   //! \retval Quantity of the state history variable
-  double state_variable(const std::string& var) const override {
-    return (state_variables_.find(var) != state_variables_.end())
-               ? state_variables_.at(var)
+  double state_variable(
+      const std::string& var,
+      unsigned phase = mpm::ParticlePhase::Solid) const override {
+    return (state_variables_[phase].find(var) != state_variables_[phase].end())
+               ? state_variables_[phase].at(var)
                : std::numeric_limits<double>::quiet_NaN();
   }
 
   //! Map particle pressure to nodes
-  bool map_pressure_to_nodes() noexcept override;
+  bool map_pressure_to_nodes(
+      unsigned phase = mpm::ParticlePhase::Solid) noexcept override;
 
   //! Compute pressure smoothing of the particle based on nodal pressure
   //! $$\hat{p}_p = \sum_{i = 1}^{n_n} N_i(x_p) p_i$$
-  bool compute_pressure_smoothing() noexcept override;
+  bool compute_pressure_smoothing(
+      unsigned phase = mpm::ParticlePhase::Solid) noexcept override;
 
   //! Return pressure of the particles
-  double pressure() const override {
-    return (state_variables_.find("pressure") != state_variables_.end())
-               ? state_variables_.at("pressure")
+  //! \param[in] phase Index to indicate phase
+  double pressure(unsigned phase = mpm::ParticlePhase::Solid) const override {
+    return (state_variables_[phase].find("pressure") !=
+            state_variables_[phase].end())
+               ? state_variables_[phase].at("pressure")
                : std::numeric_limits<double>::quiet_NaN();
   }
 
@@ -279,8 +289,20 @@ class Particle : public ParticleBase<Tdim> {
   //! Return neighbour ids
   std::vector<mpm::Index> neighbours() const override { return neighbours_; };
 
+ protected:
+  //! Initialise particle material container
+  //! \details This function allocate memory and initialise the material related
+  //! containers according to the particle phase, i.e. solid or fluid particle
+  //! has phase_size = 1, whereas two-phase (solid-fluid) or three-phase
+  //! (solid-water-air) particle have phase_size = 2 and 3, respectively.
+  //! \param[in] phase_size The material phase size
+  void initialise_material(unsigned phase_size = 1);
+
  private:
   //! Compute strain rate
+  //! \param[in] dn_dx The spatial gradient of shape function
+  //! \param[in] phase Index to indicate phase
+  //! \retval strain rate at particle inside a cell
   inline Eigen::Matrix<double, 6, 1> compute_strain_rate(
       const Eigen::MatrixXd& dn_dx, unsigned phase) noexcept;