diff --git a/include/generators/injection.h b/include/generators/injection.h
index 633fb861a..6efc09988 100644
--- a/include/generators/injection.h
+++ b/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
diff --git a/include/mesh.h b/include/mesh.h
index 78496a607..2de4f49e9 100644
--- a/include/mesh.h
+++ b/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
diff --git a/include/mesh.tcc b/include/mesh.tcc
index a59c8b046..c603a2fea 100644
--- a/include/mesh.tcc
+++ b/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");
diff --git a/include/particles/particle.h b/include/particles/particle.h
index f5f02925f..849b72f56 100644
--- a/include/particles/particle.h
+++ b/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;
 
diff --git a/include/particles/particle.tcc b/include/particles/particle.tcc
index 0fa0121ef..4b0f32e65 100644
--- a/include/particles/particle.tcc
+++ b/include/particles/particle.tcc
@@ -7,8 +7,8 @@ mpm::Particle<Tdim>::Particle(Index id, const VectorDim& coord)
   cell_ = nullptr;
   // Nodes
   nodes_.clear();
-  // Set material pointer to null
-  material_ = nullptr;
+  // Set material containers
+  this->initialise_material(1);
   // Logger
   std::string logger =
       "particle" + std::to_string(Tdim) + "d::" + std::to_string(id);
@@ -22,7 +22,8 @@ mpm::Particle<Tdim>::Particle(Index id, const VectorDim& coord, bool status)
   this->initialise();
   cell_ = nullptr;
   nodes_.clear();
-  material_ = nullptr;
+  // Set material containers
+  this->initialise_material(1);
   //! Logger
   std::string logger =
       "particle" + std::to_string(Tdim) + "d::" + std::to_string(id);
@@ -96,7 +97,7 @@ bool mpm::Particle<Tdim>::initialise_particle(const HDF5Particle& particle) {
   this->nodes_.clear();
 
   // Material id
-  this->material_id_ = particle.material_id;
+  this->material_id_[mpm::ParticlePhase::Solid] = particle.material_id;
 
   return true;
 }
@@ -108,17 +109,18 @@ bool mpm::Particle<Tdim>::initialise_particle(
     const std::shared_ptr<mpm::Material<Tdim>>& material) {
   bool status = this->initialise_particle(particle);
   if (material != nullptr) {
-    if (this->material_id_ == material->id() ||
-        this->material_id_ == std::numeric_limits<unsigned>::max()) {
+    if (this->material_id() == material->id() ||
+        this->material_id() == std::numeric_limits<unsigned>::max()) {
       bool assign_mat = this->assign_material(material);
       if (!assign_mat) throw std::runtime_error("Material assignment failed");
       // Reinitialize state variables
-      auto mat_state_vars = material_->initialise_state_variables();
+      auto mat_state_vars = (this->material())->initialise_state_variables();
       if (mat_state_vars.size() == particle.nstate_vars) {
         unsigned i = 0;
-        auto state_variables = material_->state_variables();
+        auto state_variables = (this->material())->state_variables();
         for (const auto& state_var : state_variables) {
-          this->state_variables_.at(state_var) = particle.svars[i];
+          this->state_variables_[mpm::ParticlePhase::Solid].at(state_var) =
+              particle.svars[i];
           ++i;
         }
       }
@@ -163,8 +165,9 @@ mpm::HDF5Particle mpm::Particle<Tdim>::hdf5() const {
   particle_data.mass = this->mass();
   particle_data.volume = this->volume();
   particle_data.pressure =
-      (state_variables_.find("pressure") != state_variables_.end())
-          ? state_variables_.at("pressure")
+      (state_variables_[mpm::ParticlePhase::Solid].find("pressure") !=
+       state_variables_[mpm::ParticlePhase::Solid].end())
+          ? state_variables_[mpm::ParticlePhase::Solid].at("pressure")
           : 0.;
 
   particle_data.coord_x = coordinates[0];
@@ -206,14 +209,16 @@ mpm::HDF5Particle mpm::Particle<Tdim>::hdf5() const {
   particle_data.material_id = this->material_id();
 
   // Write state variables
-  if (material_ != nullptr) {
-    particle_data.nstate_vars = state_variables_.size();
-    if (state_variables_.size() > 20)
+  if (this->material() != nullptr) {
+    particle_data.nstate_vars =
+        state_variables_[mpm::ParticlePhase::Solid].size();
+    if (state_variables_[mpm::ParticlePhase::Solid].size() > 20)
       throw std::runtime_error("# of state variables cannot be more than 20");
     unsigned i = 0;
-    auto state_variables = material_->state_variables();
+    auto state_variables = (this->material())->state_variables();
     for (const auto& state_var : state_variables) {
-      particle_data.svars[i] = state_variables_.at(state_var);
+      particle_data.svars[i] =
+          state_variables_[mpm::ParticlePhase::Solid].at(state_var);
       ++i;
     }
   }
@@ -245,19 +250,31 @@ void mpm::Particle<Tdim>::initialise() {
   this->properties_["displacements"] = [&]() { return displacement(); };
 }
 
+//! Initialise particle material container
+template <unsigned Tdim>
+void mpm::Particle<Tdim>::initialise_material(unsigned phase_size) {
+  material_.resize(phase_size);
+  material_id_.resize(phase_size);
+  state_variables_.resize(phase_size);
+  std::fill(material_.begin(), material_.end(), nullptr);
+  std::fill(material_id_.begin(), material_id_.end(),
+            std::numeric_limits<unsigned>::max());
+  std::fill(state_variables_.begin(), state_variables_.end(), mpm::dense_map());
+}
+
 //! Assign material state variables from neighbour particle
 template <unsigned Tdim>
 bool mpm::Particle<Tdim>::assign_material_state_vars(
     const mpm::dense_map& state_vars,
-    const std::shared_ptr<mpm::Material<Tdim>>& material) {
+    const std::shared_ptr<mpm::Material<Tdim>>& material, unsigned phase) {
   bool status = false;
-  if (material != nullptr && this->material_ != nullptr &&
-      this->material_id_ == material->id()) {
+  if (material != nullptr && this->material() != nullptr &&
+      this->material_id() == material->id()) {
     // Clone state variables
-    auto mat_state_vars = material_->initialise_state_variables();
-    if (this->state_variables_.size() == state_vars.size() &&
+    auto mat_state_vars = (this->material())->initialise_state_variables();
+    if (state_variables_[phase].size() == state_vars.size() &&
         mat_state_vars.size() == state_vars.size()) {
-      this->state_variables_ = state_vars;
+      this->state_variables_[phase] = state_vars;
       status = true;
     }
   }
@@ -373,14 +390,15 @@ void mpm::Particle<Tdim>::remove_cell() {
 // Assign a material to particle
 template <unsigned Tdim>
 bool mpm::Particle<Tdim>::assign_material(
-    const std::shared_ptr<Material<Tdim>>& material) {
+    const std::shared_ptr<Material<Tdim>>& material, unsigned phase) {
   bool status = false;
   try {
     // Check if material is valid and properties are set
     if (material != nullptr) {
-      material_ = material;
-      material_id_ = material_->id();
-      state_variables_ = material_->initialise_state_variables();
+      material_.at(phase) = material;
+      material_id_.at(phase) = material_[phase]->id();
+      state_variables_.at(phase) =
+          material_[phase]->initialise_state_variables();
       status = true;
     } else {
       throw std::runtime_error("Material is undefined!");
@@ -482,10 +500,11 @@ void mpm::Particle<Tdim>::update_volume() noexcept {
 template <unsigned Tdim>
 void mpm::Particle<Tdim>::compute_mass() noexcept {
   // Check if particle volume is set and material ptr is valid
-  assert(volume_ != std::numeric_limits<double>::max() && material_ != nullptr);
+  assert(volume_ != std::numeric_limits<double>::max() &&
+         this->material() != nullptr);
   // Mass = volume of particle * mass_density
   this->mass_density_ =
-      material_->template property<double>(std::string("density"));
+      (this->material())->template property<double>(std::string("density"));
   this->mass_ = volume_ * mass_density_;
 }
 
@@ -516,9 +535,10 @@ void mpm::Particle<Tdim>::map_multimaterial_mass_momentum_to_nodes() noexcept {
   // Map mass and momentum to nodal property taking into account the material id
   for (unsigned i = 0; i < nodes_.size(); ++i) {
     nodal_mass(0, 0) = mass_ * shapefn_[i];
-    nodes_[i]->update_property(true, "masses", nodal_mass, material_id_, 1);
+    nodes_[i]->update_property(true, "masses", nodal_mass, this->material_id(),
+                               1);
     nodes_[i]->update_property(true, "momenta", velocity_ * nodal_mass,
-                               material_id_, Tdim);
+                               this->material_id(), Tdim);
   }
 }
 
@@ -533,7 +553,7 @@ void mpm::Particle<Tdim>::map_multimaterial_displacements_to_nodes() noexcept {
   for (unsigned i = 0; i < nodes_.size(); ++i) {
     const auto& displacement = mass_ * shapefn_[i] * displacement_;
     nodes_[i]->update_property(true, "displacements", displacement,
-                               material_id_, Tdim);
+                               this->material_id(), Tdim);
   }
 }
 
@@ -549,8 +569,8 @@ void mpm::Particle<
   for (unsigned i = 0; i < nodes_.size(); ++i) {
     Eigen::Matrix<double, Tdim, 1> gradient;
     for (unsigned j = 0; j < Tdim; ++j) gradient[j] = volume_ * dn_dx_(i, j);
-    nodes_[i]->update_property(true, "domain_gradients", gradient, material_id_,
-                               Tdim);
+    nodes_[i]->update_property(true, "domain_gradients", gradient,
+                               this->material_id(), Tdim);
   }
 }
 
@@ -636,10 +656,12 @@ void mpm::Particle<Tdim>::compute_strain(double dt) noexcept {
 template <unsigned Tdim>
 void mpm::Particle<Tdim>::compute_stress() noexcept {
   // Check if material ptr is valid
-  assert(material_ != nullptr);
+  assert(this->material() != nullptr);
   // Calculate stress
   this->stress_ =
-      material_->compute_stress(stress_, dstrain_, this, &state_variables_);
+      (this->material())
+          ->compute_stress(stress_, dstrain_, this,
+                           &state_variables_[mpm::ParticlePhase::Solid]);
 }
 
 //! Map body force
@@ -781,19 +803,19 @@ void mpm::Particle<Tdim>::compute_updated_position(
 
 //! Map particle pressure to nodes
 template <unsigned Tdim>
-bool mpm::Particle<Tdim>::map_pressure_to_nodes() noexcept {
+bool mpm::Particle<Tdim>::map_pressure_to_nodes(unsigned phase) noexcept {
   // Mass is initialized
   assert(mass_ != std::numeric_limits<double>::max());
 
   bool status = false;
   // Check if particle mass is set and state variable pressure is found
   if (mass_ != std::numeric_limits<double>::max() &&
-      (state_variables_.find("pressure") != state_variables_.end())) {
+      (state_variables_[phase].find("pressure") !=
+       state_variables_[phase].end())) {
     // Map particle pressure to nodes
     for (unsigned i = 0; i < nodes_.size(); ++i)
       nodes_[i]->update_mass_pressure(
-          mpm::ParticlePhase::Solid,
-          shapefn_[i] * mass_ * state_variables_["pressure"]);
+          phase, shapefn_[i] * mass_ * state_variables_[phase]["pressure"]);
 
     status = true;
   }
@@ -802,21 +824,21 @@ bool mpm::Particle<Tdim>::map_pressure_to_nodes() noexcept {
 
 // Compute pressure smoothing of the particle based on nodal pressure
 template <unsigned Tdim>
-bool mpm::Particle<Tdim>::compute_pressure_smoothing() noexcept {
+bool mpm::Particle<Tdim>::compute_pressure_smoothing(unsigned phase) noexcept {
   // Assert
   assert(cell_ != nullptr);
 
   bool status = false;
   // Check if particle has a valid cell ptr
-  if (cell_ != nullptr &&
-      (state_variables_.find("pressure") != state_variables_.end())) {
+  if (cell_ != nullptr && (state_variables_[phase].find("pressure") !=
+                           state_variables_[phase].end())) {
 
     double pressure = 0.;
     // Update particle pressure to interpolated nodal pressure
     for (unsigned i = 0; i < this->nodes_.size(); ++i)
-      pressure += shapefn_[i] * nodes_[i]->pressure(mpm::ParticlePhase::Solid);
+      pressure += shapefn_[i] * nodes_[i]->pressure(phase);
 
-    state_variables_["pressure"] = pressure;
+    state_variables_[phase]["pressure"] = pressure;
     status = true;
   }
   return status;
@@ -840,7 +862,7 @@ Eigen::VectorXd mpm::Particle<Tdim>::tensor_data(const std::string& property) {
 template <unsigned Tdim>
 void mpm::Particle<Tdim>::append_material_id_to_nodes() const {
   for (unsigned i = 0; i < nodes_.size(); ++i)
-    nodes_[i]->append_material_id(material_id_);
+    nodes_[i]->append_material_id(this->material_id());
 }
 
 //! Assign neighbour particles
diff --git a/include/particles/particle_base.h b/include/particles/particle_base.h
index d8e006fb3..da0d22a99 100644
--- a/include/particles/particle_base.h
+++ b/include/particles/particle_base.h
@@ -67,10 +67,12 @@ class ParticleBase {
   //! 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
   virtual bool assign_material_state_vars(
       const mpm::dense_map& state_vars,
-      const std::shared_ptr<mpm::Material<Tdim>>& material) = 0;
+      const std::shared_ptr<mpm::Material<Tdim>>& material,
+      unsigned phase = mpm::ParticlePhase::Solid) = 0;
 
   //! Retrun particle data as HDF5
   //! \retval particle HDF5 data of the particle
@@ -149,11 +151,28 @@ class ParticleBase {
   virtual void map_multimaterial_domain_gradients_to_nodes() noexcept = 0;
 
   //! Assign material
-  virtual bool assign_material(
-      const std::shared_ptr<Material<Tdim>>& material) = 0;
+  virtual bool assign_material(const std::shared_ptr<Material<Tdim>>& material,
+                               unsigned phase = mpm::ParticlePhase::Solid) = 0;
+
+  //! Return material of particle
+  //! \param[in] phase Index to indicate material phase
+  std::shared_ptr<Material<Tdim>> material(
+      unsigned phase = mpm::ParticlePhase::Solid) const {
+    return material_[phase];
+  }
 
   //! Return material id
-  unsigned material_id() const { return material_id_; }
+  //! \param[in] phase Index to indicate material phase
+  unsigned material_id(unsigned phase = mpm::ParticlePhase::Solid) const {
+    return material_id_[phase];
+  }
+
+  //! Return state variables
+  //! \param[in] phase Index to indicate material phase
+  mpm::dense_map state_variables(
+      unsigned phase = mpm::ParticlePhase::Solid) const {
+    return state_variables_[phase];
+  }
 
   //! Assign status
   void assign_status(bool status) { status_ = status; }
@@ -171,7 +190,7 @@ class ParticleBase {
   virtual double mass() const = 0;
 
   //! Return pressure
-  virtual double pressure() const = 0;
+  virtual double pressure(unsigned phase = mpm::ParticlePhase::Solid) const = 0;
 
   //! Compute strain
   virtual void compute_strain(double dt) noexcept = 0;
@@ -204,10 +223,12 @@ class ParticleBase {
   virtual void map_internal_force() noexcept = 0;
 
   //! Map particle pressure to nodes
-  virtual bool map_pressure_to_nodes() noexcept = 0;
+  virtual bool map_pressure_to_nodes(
+      unsigned phase = mpm::ParticlePhase::Solid) noexcept = 0;
 
   //! Compute pressure smoothing of the particle based on nodal pressure
-  virtual bool compute_pressure_smoothing() noexcept = 0;
+  virtual bool compute_pressure_smoothing(
+      unsigned phase = mpm::ParticlePhase::Solid) noexcept = 0;
 
   //! Assign velocity
   virtual bool assign_velocity(const VectorDim& velocity) = 0;
@@ -232,7 +253,9 @@ class ParticleBase {
       double dt, bool velocity_update = false) noexcept = 0;
 
   //! Return a state variable
-  virtual double state_variable(const std::string& var) const = 0;
+  virtual double state_variable(
+      const std::string& var,
+      unsigned phase = mpm::ParticlePhase::Solid) const = 0;
 
   //! Return tensor data of particles
   //! \param[in] property Property string
@@ -275,11 +298,11 @@ class ParticleBase {
   //! Vector of nodal pointers
   std::vector<std::shared_ptr<NodeBase<Tdim>>> nodes_;
   //! Material
-  std::shared_ptr<Material<Tdim>> material_;
+  std::vector<std::shared_ptr<Material<Tdim>>> material_;
   //! Unsigned material id
-  unsigned material_id_{std::numeric_limits<unsigned>::max()};
+  std::vector<unsigned> material_id_;
   //! Material state history variables
-  mpm::dense_map state_variables_;
+  std::vector<mpm::dense_map> state_variables_;
   //! Vector of particle neighbour ids
   std::vector<mpm::Index> neighbours_;
 };  // ParticleBase class
diff --git a/include/solvers/mpm_base.h b/include/solvers/mpm_base.h
index 162271ac1..ca017f1fc 100644
--- a/include/solvers/mpm_base.h
+++ b/include/solvers/mpm_base.h
@@ -84,6 +84,10 @@ class MPMBase : public MPM {
   //! \param[in] initial_step Start of simulation or later steps
   void mpi_domain_decompose(bool initial_step = false) override;
 
+  //! Pressure smoothing
+  //! \param[in] phase Phase to smooth pressure
+  void pressure_smoothing(unsigned phase);
+
  private:
   //! Return if a mesh will be isoparametric or not
   //! \retval isoparametric Status of mesh type
diff --git a/include/solvers/mpm_base.tcc b/include/solvers/mpm_base.tcc
index c236cbd58..8426ccf4f 100644
--- a/include/solvers/mpm_base.tcc
+++ b/include/solvers/mpm_base.tcc
@@ -336,17 +336,20 @@ bool mpm::MPMBase<Tdim>::initialise_particles() {
         for (const auto& material_set : material_sets) {
           unsigned material_id =
               material_set["material_id"].template get<unsigned>();
+          unsigned phase_id = mpm::ParticlePhase::Solid;
+          if (material_set.contains("phase_id"))
+            phase_id = material_set["phase_id"].template get<unsigned>();
           unsigned pset_id = material_set["pset_id"].template get<unsigned>();
           // Update material_id for particles in each pset
           mesh_->iterate_over_particle_set(
-              pset_id,
-              std::bind(&mpm::ParticleBase<Tdim>::assign_material,
-                        std::placeholders::_1, materials_.at(material_id)));
+              pset_id, std::bind(&mpm::ParticleBase<Tdim>::assign_material,
+                                 std::placeholders::_1,
+                                 materials_.at(material_id), phase_id));
         }
       }
     } catch (std::exception& exception) {
-      console_->warn("{} #{}: Material sets are not specified", __FILE__, __LINE__,
-                     exception.what());
+      console_->warn("{} #{}: Material sets are not specified", __FILE__,
+                     __LINE__, exception.what());
     }
 
   } catch (std::exception& exception) {
@@ -1157,3 +1160,33 @@ void mpm::MPMBase<Tdim>::mpi_domain_decompose(bool initial_step) {
   }
 #endif  // MPI
 }
+
+//! MPM pressure smoothing
+template <unsigned Tdim>
+void mpm::MPMBase<Tdim>::pressure_smoothing(unsigned phase) {
+  // Assign pressure to nodes
+  mesh_->iterate_over_particles(
+      std::bind(&mpm::ParticleBase<Tdim>::map_pressure_to_nodes,
+                std::placeholders::_1, phase));
+
+#ifdef USE_MPI
+  int mpi_size = 1;
+
+  // Get number of MPI ranks
+  MPI_Comm_size(MPI_COMM_WORLD, &mpi_size);
+
+  // Run if there is more than a single MPI task
+  if (mpi_size > 1) {
+    // MPI all reduce nodal pressure
+    mesh_->template nodal_halo_exchange<double, 1>(
+        std::bind(&mpm::NodeBase<Tdim>::pressure, std::placeholders::_1, phase),
+        std::bind(&mpm::NodeBase<Tdim>::assign_pressure, std::placeholders::_1,
+                  phase, std::placeholders::_2));
+  }
+#endif
+
+  // Smooth pressure over particles
+  mesh_->iterate_over_particles(
+      std::bind(&mpm::ParticleBase<Tdim>::compute_pressure_smoothing,
+                std::placeholders::_1, phase));
+}
diff --git a/include/solvers/mpm_explicit.h b/include/solvers/mpm_explicit.h
index 60ee53040..47ae93446 100644
--- a/include/solvers/mpm_explicit.h
+++ b/include/solvers/mpm_explicit.h
@@ -22,10 +22,6 @@ class MPMExplicit : public MPMBase<Tdim> {
   //! Solve
   bool solve() override;
 
-  //! Pressure smoothing
-  //! \param[in] phase Phase to smooth pressure
-  void pressure_smoothing(unsigned phase);
-
   //! Compute stress strain
   //! \param[in] phase Phase to smooth pressure
   void compute_stress_strain(unsigned phase);
diff --git a/include/solvers/mpm_explicit.tcc b/include/solvers/mpm_explicit.tcc
index 05ddf36c2..67261d041 100644
--- a/include/solvers/mpm_explicit.tcc
+++ b/include/solvers/mpm_explicit.tcc
@@ -6,35 +6,6 @@ mpm::MPMExplicit<Tdim>::MPMExplicit(const std::shared_ptr<IO>& io)
   console_ = spdlog::get("MPMExplicit");
 }
 
-//! MPM Explicit pressure smoothing
-template <unsigned Tdim>
-void mpm::MPMExplicit<Tdim>::pressure_smoothing(unsigned phase) {
-  // Assign pressure to nodes
-  mesh_->iterate_over_particles(std::bind(
-      &mpm::ParticleBase<Tdim>::map_pressure_to_nodes, std::placeholders::_1));
-
-#ifdef USE_MPI
-  int mpi_size = 1;
-
-  // Get number of MPI ranks
-  MPI_Comm_size(MPI_COMM_WORLD, &mpi_size);
-
-  // Run if there is more than a single MPI task
-  if (mpi_size > 1) {
-    // MPI all reduce nodal pressure
-    mesh_->template nodal_halo_exchange<double, 1>(
-        std::bind(&mpm::NodeBase<Tdim>::pressure, std::placeholders::_1, phase),
-        std::bind(&mpm::NodeBase<Tdim>::assign_pressure, std::placeholders::_1,
-                  phase, std::placeholders::_2));
-  }
-#endif
-
-  // Smooth pressure over particles
-  mesh_->iterate_over_particles(
-      std::bind(&mpm::ParticleBase<Tdim>::compute_pressure_smoothing,
-                std::placeholders::_1));
-}
-
 //! MPM Explicit compute stress strain
 template <unsigned Tdim>
 void mpm::MPMExplicit<Tdim>::compute_stress_strain(unsigned phase) {
diff --git a/tests/io/write_mesh_particles.cc b/tests/io/write_mesh_particles.cc
index e5e6b0c51..555d5abb5 100644
--- a/tests/io/write_mesh_particles.cc
+++ b/tests/io/write_mesh_particles.cc
@@ -65,7 +65,8 @@ bool write_json(unsigned dim, bool resume, const std::string& analysis,
          {"density", 2300.},
          {"youngs_modulus", 1.5E+6},
          {"poisson_ratio", 0.25}}}},
-      {"material_sets", {{{"material_id", 1}, {"pset_id", 2}}}},
+      {"material_sets",
+       {{{"material_id", 1}, {"phase_id", 0}, {"pset_id", 2}}}},
       {"external_loading_conditions",
        {{"gravity", gravity},
         {"particle_surface_traction",
diff --git a/tests/io/write_mesh_particles_unitcell.cc b/tests/io/write_mesh_particles_unitcell.cc
index d3f970889..3767f52b1 100644
--- a/tests/io/write_mesh_particles_unitcell.cc
+++ b/tests/io/write_mesh_particles_unitcell.cc
@@ -14,7 +14,7 @@ bool write_json_unitcell(unsigned dim, const std::string& analysis,
   auto cell_type = "ED2Q4";
   auto io_type = "Ascii2D";
   std::string material = "LinearElastic2D";
-  unsigned material_id = 1;
+  std::vector<unsigned> material_id{{1}};
   std::vector<double> gravity{{0., -9.81}};
   std::vector<double> xvalues{{0.0, 0.5, 1.0}};
   std::vector<double> fxvalues{{0.0, 1.0, 1.0}};
diff --git a/tests/mesh_test_2d.cc b/tests/mesh_test_2d.cc
index 7b564bd94..c82b04f2e 100644
--- a/tests/mesh_test_2d.cc
+++ b/tests/mesh_test_2d.cc
@@ -51,6 +51,7 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") {
 
   // Assign material
   unsigned mid = 0;
+  std::vector<unsigned> mids(1, mid);
   // Initialise material
   Json jmaterial;
   jmaterial["density"] = 1000.;
@@ -530,7 +531,7 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") {
     // Generate material points in cell
     REQUIRE(mesh->nparticles() == 0);
 
-    REQUIRE(mesh->generate_material_points(1, particle_type, mid, -1, 0) ==
+    REQUIRE(mesh->generate_material_points(1, particle_type, mids, -1, 0) ==
             false);
     REQUIRE(mesh->nparticles() == 0);
 
@@ -539,19 +540,19 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") {
 
     SECTION("Check generating 1 particle / cell") {
       // Generate material points in cell
-      REQUIRE(mesh->generate_material_points(1, particle_type, mid, -1, 0) ==
+      REQUIRE(mesh->generate_material_points(1, particle_type, mids, -1, 0) ==
               true);
       REQUIRE(mesh->nparticles() == 1);
     }
 
     SECTION("Check generating 2 particle / cell") {
-      REQUIRE(mesh->generate_material_points(2, particle_type, mid, -1, 0) ==
+      REQUIRE(mesh->generate_material_points(2, particle_type, mids, -1, 0) ==
               true);
       REQUIRE(mesh->nparticles() == 4);
     }
 
     SECTION("Check generating 3 particle / cell") {
-      REQUIRE(mesh->generate_material_points(3, particle_type, mid, -1, 0) ==
+      REQUIRE(mesh->generate_material_points(3, particle_type, mids, -1, 0) ==
               true);
       REQUIRE(mesh->nparticles() == 9);
     }
@@ -577,7 +578,7 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") {
         // Gauss point generation
         Json jgen;
         jgen["type"] = "gauss";
-        jgen["material_id"] = mid;
+        jgen["material_id"] = {0};
         jgen["cset_id"] = 1;
         jgen["particle_type"] = "P2D";
         jgen["check_duplicates"] = false;
@@ -594,7 +595,7 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") {
         // Gauss point generation
         Json jgen;
         jgen["type"] = "inject";
-        jgen["material_id"] = mid;
+        jgen["material_id"] = {0};
         jgen["cset_id"] = 1;
         jgen["particle_type"] = "P2D";
         jgen["check_duplicates"] = false;
@@ -781,7 +782,7 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") {
             // Particle type 2D
             const std::string particle_type = "P2D";
             // Create particles from file
-            mesh->create_particles(particle_type, coordinates, mid, 0, false);
+            mesh->create_particles(particle_type, coordinates, mids, 0, false);
             // Check if mesh has added particles
             REQUIRE(mesh->nparticles() == coordinates.size());
 
@@ -840,7 +841,7 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") {
             unsigned nparticles = coordinates.size();
             coordinates.clear();
             // This fails with empty list error in particle creation
-            mesh->create_particles(particle_type, coordinates, mid, 0, false);
+            mesh->create_particles(particle_type, coordinates, mids, 0, false);
             REQUIRE(mesh->nparticles() == nparticles);
 
             const unsigned phase = 0;
diff --git a/tests/mesh_test_3d.cc b/tests/mesh_test_3d.cc
index 56d05b503..a6c1db5bb 100644
--- a/tests/mesh_test_3d.cc
+++ b/tests/mesh_test_3d.cc
@@ -43,6 +43,7 @@ TEST_CASE("Mesh is checked for 3D case", "[mesh][3D]") {
 
   // Assign material
   unsigned mid = 0;
+  std::vector<unsigned> mids(1, mid);
   // Initialise material
   Json jmaterial;
   jmaterial["density"] = 1000.;
@@ -587,7 +588,7 @@ TEST_CASE("Mesh is checked for 3D case", "[mesh][3D]") {
 
     REQUIRE(mesh->nparticles() == 0);
     // Generate material points in cell
-    REQUIRE(mesh->generate_material_points(1, particle_type, mid, -1, 0) ==
+    REQUIRE(mesh->generate_material_points(1, particle_type, mids, -1, 0) ==
             false);
     REQUIRE(mesh->nparticles() == 0);
 
@@ -596,19 +597,19 @@ TEST_CASE("Mesh is checked for 3D case", "[mesh][3D]") {
 
     SECTION("Check generating 1 particle / cell") {
       // Generate material points in cell
-      REQUIRE(mesh->generate_material_points(1, particle_type, mid, -1, 0) ==
+      REQUIRE(mesh->generate_material_points(1, particle_type, mids, -1, 0) ==
               true);
       REQUIRE(mesh->nparticles() == 1);
     }
 
     SECTION("Check generating 2 particle / cell") {
-      REQUIRE(mesh->generate_material_points(2, particle_type, mid, -1, 0) ==
+      REQUIRE(mesh->generate_material_points(2, particle_type, mids, -1, 0) ==
               true);
       REQUIRE(mesh->nparticles() == 8);
     }
 
     SECTION("Check generating 3 particle / cell") {
-      REQUIRE(mesh->generate_material_points(3, particle_type, mid, -1, 0) ==
+      REQUIRE(mesh->generate_material_points(3, particle_type, mids, -1, 0) ==
               true);
       REQUIRE(mesh->nparticles() == 27);
     }
@@ -877,7 +878,7 @@ TEST_CASE("Mesh is checked for 3D case", "[mesh][3D]") {
             // Particle type 3D
             const std::string particle_type = "P3D";
             // Create particles from file
-            mesh->create_particles(particle_type, coordinates, mid, 0, false);
+            mesh->create_particles(particle_type, coordinates, mids, 0, false);
             // Check if mesh has added particles
             REQUIRE(mesh->nparticles() == coordinates.size());
             // Clear coordinates and try creating a list of particles with an
@@ -885,7 +886,7 @@ TEST_CASE("Mesh is checked for 3D case", "[mesh][3D]") {
             unsigned nparticles = coordinates.size();
             coordinates.clear();
             // This fails with empty list error in particle creation
-            mesh->create_particles(particle_type, coordinates, mid, 1, false);
+            mesh->create_particles(particle_type, coordinates, mids, 1, false);
             REQUIRE(mesh->nparticles() == nparticles);
 
             // Test assign particles cells again should fail
diff --git a/tests/mpm_explicit_usf_test.cc b/tests/mpm_explicit_usf_test.cc
index 2256f7482..4670037c4 100644
--- a/tests/mpm_explicit_usf_test.cc
+++ b/tests/mpm_explicit_usf_test.cc
@@ -88,6 +88,15 @@ TEST_CASE("MPM 2D Explicit implementation is checked",
     // Solve
     REQUIRE(mpm->solve() == true);
   }
+
+  SECTION("Check pressure smoothing") {
+    // Create an IO object
+    auto io = std::make_unique<mpm::IO>(argc, argv);
+    // Run explicit MPM
+    auto mpm = std::make_unique<mpm::MPMExplicit<Dim>>(std::move(io));
+    // Pressure smoothing
+    REQUIRE_NOTHROW(mpm->pressure_smoothing(0));
+  }
 }
 
 // Check MPM Explicit
@@ -168,4 +177,13 @@ TEST_CASE("MPM 3D Explicit implementation is checked",
     // Solve
     REQUIRE(mpm->solve() == true);
   }
+
+  SECTION("Check pressure smoothing") {
+    // Create an IO object
+    auto io = std::make_unique<mpm::IO>(argc, argv);
+    // Run explicit MPM
+    auto mpm = std::make_unique<mpm::MPMExplicit<Dim>>(std::move(io));
+    // Pressure smoothing
+    REQUIRE_NOTHROW(mpm->pressure_smoothing(0));
+  }
 }
diff --git a/tests/particle_cell_crossing_test.cc b/tests/particle_cell_crossing_test.cc
index 05ca405fe..ea3a1eb35 100644
--- a/tests/particle_cell_crossing_test.cc
+++ b/tests/particle_cell_crossing_test.cc
@@ -162,7 +162,7 @@ TEST_CASE("Particle cell crossing is checked for 2D case",
   // Iterate over each particle to assign material
   mesh->iterate_over_particles(
       std::bind(&mpm::ParticleBase<Dim>::assign_material, std::placeholders::_1,
-                material));
+                material, mpm::ParticlePhase::Solid));
 
   // Compute volume
   mesh->iterate_over_particles(std::bind(
@@ -418,7 +418,7 @@ TEST_CASE("Particle cell crossing is checked for 3D case",
   // Iterate over each particle to assign material
   mesh->iterate_over_particles(
       std::bind(&mpm::ParticleBase<Dim>::assign_material, std::placeholders::_1,
-                material));
+                material, mpm::ParticlePhase::Solid));
 
   // Compute volume
   mesh->iterate_over_particles(std::bind(