[Feature] State variable vtk output for different phases

include/mesh.h

@@ -266,7 +266,8 @@ class Mesh {
 
   //! Return particles scalar data
   //! \param[in] attribute Name of the scalar data attribute
-  std::vector<double> particles_statevars_data(const std::string& attribute);
+  std::vector<double> particles_statevars_data(
+      const std::string& attribute, unsigned phase = mpm::ParticlePhase::Solid);
 
   //! Compute and assign rotation matrix to nodes
   //! \param[in] euler_angles Map of node number and respective euler_angles

include/mesh.tcc

@@ -1104,12 +1104,12 @@ std::vector<Eigen::Matrix<double, Tsize, 1>>
 //! Return particle scalar data
 template <unsigned Tdim>
 std::vector<double> mpm::Mesh<Tdim>::particles_statevars_data(
-    const std::string& attribute) {
+    const std::string& attribute, unsigned phase) {
   std::vector<double> scalar_data;
   scalar_data.reserve(particles_.size());
   // Iterate over particles and add scalar value to data
   for (auto pitr = particles_.cbegin(); pitr != particles_.cend(); ++pitr)
-    scalar_data.emplace_back((*pitr)->state_variable(attribute));
+    scalar_data.emplace_back((*pitr)->state_variable(attribute, phase));
   return scalar_data;
 }
 

include/solvers/mpm_base.h

@@ -193,7 +193,7 @@ class MPMBase : public MPM {
   //! Mathematical functions
   std::map<unsigned, std::shared_ptr<mpm::FunctionBase>> math_functions_;
   //! VTK state variables
-  std::vector<std::string> vtk_statevars_;
+  std::map<unsigned, std::vector<std::string>> vtk_statevars_;
   //! Set node concentrated force
   bool set_node_concentrated_force_{false};
   //! Damping type

include/solvers/mpm_base.tcc

@@ -96,10 +96,44 @@ mpm::MPMBase<Tdim>::MPMBase(const std::shared_ptr<IO>& io) : mpm::MPM(io) {
   try {
     if (post_process_.at("vtk_statevars").is_array() &&
         post_process_.at("vtk_statevars").size() > 0) {
-      for (unsigned i = 0; i < post_process_.at("vtk_statevars").size(); ++i) {
-        std::string attribute =
-            post_process_["vtk_statevars"][i].template get<std::string>();
-        vtk_statevars_.emplace_back(attribute);
+      // Iterate over state_vars
+      for (const auto& svars : post_process_["vtk_statevars"]) {
+        // Initial phase id and state vars
+        unsigned phase_id = 0;
+        std::vector<std::string> state_var;
+
+        // If state variables are arranged as string
+        if (svars.is_string()) {
+          std::string attribute = svars.template get<std::string>();
+          state_var.emplace_back(attribute);
+        }
+        // If state variables are arranged as json object
+        else if (svars.is_object()) {
+          // Phase id
+          if (svars.contains("phase_id"))
+            phase_id = svars.at("phase_id").template get<unsigned>();
+
+          // State variables
+          if (svars.at("statevars").is_array() &&
+              svars.at("statevars").size() > 0) {
+            for (unsigned i = 0; i < svars.at("statevars").size(); ++i) {
+              std::string attribute =
+                  svars["statevars"][i].template get<std::string>();
+              state_var.emplace_back(attribute);
+            }
+          } else {
+            throw std::runtime_error(
+                "No VTK statevariable were specified, none will be generated");
+          }
+        }
+
+        // Check if vtk_statevars_ has the designated phase_id
+        if (vtk_statevars_.find(phase_id) != vtk_statevars_.end())
+          vtk_statevars_.at(phase_id).insert(vtk_statevars_.at(phase_id).end(),
+                                             state_var.begin(),
+                                             state_var.end());
+        else
+          vtk_statevars_.insert(std::make_pair(phase_id, state_var));
       }
     } else {
       throw std::runtime_error(
@@ -554,23 +588,29 @@ void mpm::MPMBase<Tdim>::write_vtk(mpm::Index step, mpm::Index max_steps) {
   }
 
   // VTK state variables
-  for (const auto& attribute : vtk_statevars_) {
-    // Write state variables
-    auto file =
-        io_->output_file(attribute, extension, uuid_, step, max_steps).string();
-    vtk_writer->write_scalar_point_data(
-        file, mesh_->particles_statevars_data(attribute), attribute);
-    // Write a parallel MPI VTK container file
+  for (auto const& vtk_statevar : vtk_statevars_) {
+    unsigned phase_id = vtk_statevar.first;
+    for (const auto& attribute : vtk_statevar.second) {
+      std::string phase_attribute = "P" + std::to_string(phase_id) + attribute;
+      // Write state variables
+      auto file =
+          io_->output_file(phase_attribute, extension, uuid_, step, max_steps)
+              .string();
+      vtk_writer->write_scalar_point_data(
+          file, mesh_->particles_statevars_data(attribute, phase_id),
+          phase_attribute);
+      // Write a parallel MPI VTK container file
 #ifdef USE_MPI
-    if (mpi_rank == 0 && mpi_size > 1) {
-      auto parallel_file = io_->output_file(attribute, ".pvtp", uuid_, step,
-                                            max_steps, write_mpi_rank)
-                               .string();
-      unsigned ncomponents = 1;
-      vtk_writer->write_parallel_vtk(parallel_file, attribute, mpi_size, step,
-                                     max_steps, ncomponents);
-    }
+      if (mpi_rank == 0 && mpi_size > 1) {
+        auto parallel_file = io_->output_file(phase_attribute, ".pvtp", uuid_,
+                                              step, max_steps, write_mpi_rank)
+                                 .string();
+        unsigned ncomponents = 1;
+        vtk_writer->write_parallel_vtk(parallel_file, phase_attribute, mpi_size,
+                                       step, max_steps, ncomponents);
+      }
 #endif
+    }
   }
 }
 #endif