[Solver] Two-Phase One-Point Explicit

.gitignore

@@ -76,3 +76,4 @@ configure.scan
 *missing
 *stamp-h1
 benchmarks/
+bin/

CMakeLists.txt

@@ -142,6 +142,7 @@ include_directories(BEFORE
   ${mpm_SOURCE_DIR}/include/materials/
   ${mpm_SOURCE_DIR}/include/particles/
   ${mpm_SOURCE_DIR}/include/solvers/
+  ${mpm_SOURCE_DIR}/include/utilities/
   ${mpm_SOURCE_DIR}/external/
   ${mpm_SOURCE_DIR}/tests/include/
 )
@@ -156,6 +157,7 @@ SET(mpm_src
   ${mpm_SOURCE_DIR}/src/functions/sin_function.cc
   ${mpm_SOURCE_DIR}/src/geometry.cc
   ${mpm_SOURCE_DIR}/src/hdf5_particle.cc
+  ${mpm_SOURCE_DIR}/src/hdf5_particle_twophase.cc
   ${mpm_SOURCE_DIR}/src/io/io.cc
   ${mpm_SOURCE_DIR}/src/io/io_mesh.cc
   ${mpm_SOURCE_DIR}/src/io/logger.cc
@@ -206,20 +208,27 @@ if(MPM_BUILD_TESTING)
     ${mpm_SOURCE_DIR}/tests/mesh_test_2d.cc
     ${mpm_SOURCE_DIR}/tests/mesh_test_3d.cc
     ${mpm_SOURCE_DIR}/tests/mpi_particle_test.cc
+    ${mpm_SOURCE_DIR}/tests/mpi_particle_twophase_test.cc
     ${mpm_SOURCE_DIR}/tests/mpi_transfer_particle_test.cc
-    ${mpm_SOURCE_DIR}/tests/mpm_explicit_usf_test.cc
-    ${mpm_SOURCE_DIR}/tests/mpm_explicit_usf_unitcell_test.cc
-    ${mpm_SOURCE_DIR}/tests/mpm_explicit_usl_test.cc
-    ${mpm_SOURCE_DIR}/tests/mpm_explicit_usl_unitcell_test.cc
-    ${mpm_SOURCE_DIR}/tests/nodal_properties_test.cc
-    ${mpm_SOURCE_DIR}/tests/node_map_test.cc
-    ${mpm_SOURCE_DIR}/tests/node_test.cc
-    ${mpm_SOURCE_DIR}/tests/node_vector_test.cc
-    ${mpm_SOURCE_DIR}/tests/particle_cell_crossing_test.cc
-    ${mpm_SOURCE_DIR}/tests/particle_test.cc
-    ${mpm_SOURCE_DIR}/tests/particle_traction_test.cc
-    ${mpm_SOURCE_DIR}/tests/particle_vector_test.cc
+    ${mpm_SOURCE_DIR}/tests/nodes/nodal_properties_test.cc
+    ${mpm_SOURCE_DIR}/tests/nodes/node_map_test.cc
+    ${mpm_SOURCE_DIR}/tests/nodes/node_test.cc
+    ${mpm_SOURCE_DIR}/tests/nodes/node_twophase_test.cc
+    ${mpm_SOURCE_DIR}/tests/nodes/node_vector_test.cc
+    ${mpm_SOURCE_DIR}/tests/particles/particle_cell_crossing_test.cc
+    ${mpm_SOURCE_DIR}/tests/particles/particle_test.cc
+    ${mpm_SOURCE_DIR}/tests/particles/particle_twophase_test.cc
+    ${mpm_SOURCE_DIR}/tests/particles/particle_traction_test.cc
+    ${mpm_SOURCE_DIR}/tests/particles/particle_vector_test.cc
     ${mpm_SOURCE_DIR}/tests/point_in_cell_test.cc
+    ${mpm_SOURCE_DIR}/tests/solvers/mpm_explicit_usf_test.cc
+    ${mpm_SOURCE_DIR}/tests/solvers/mpm_explicit_usf_unitcell_test.cc
+    ${mpm_SOURCE_DIR}/tests/solvers/mpm_explicit_usl_test.cc
+    ${mpm_SOURCE_DIR}/tests/solvers/mpm_explicit_usl_unitcell_test.cc
+    ${mpm_SOURCE_DIR}/tests/solvers/mpm_explicit_twophase_usf_test.cc
+    ${mpm_SOURCE_DIR}/tests/solvers/mpm_explicit_twophase_usf_unitcell_test.cc
+    ${mpm_SOURCE_DIR}/tests/solvers/mpm_explicit_twophase_usl_test.cc
+    ${mpm_SOURCE_DIR}/tests/solvers/mpm_explicit_twophase_usl_unitcell_test.cc
   )
   add_executable(mpmtest ${mpm_src} ${test_src})
   add_test(NAME mpmtest COMMAND $<TARGET_FILE:mpmtest>)

include/cell.h

@@ -217,6 +217,26 @@ class Cell {
   //! Return previous mpi rank
   unsigned previous_mpirank() const;
 
+  //! Assign free surface
+  //! \param[in] free_surface boolean indicating free surface cell
+  void assign_free_surface(bool free_surface) { free_surface_ = free_surface; };
+
+  //! Return free surface bool
+  //! \retval free_surface_ indicating free surface cell
+  bool free_surface() const { return free_surface_; };
+
+  //! Assign volume traction to node
+  //! \param[in] volume_fraction cell volume fraction
+  void assign_volume_fraction(double volume_fraction);
+
+  //! Return cell volume fraction
+  //! \retval volume_fraction_ cell volume fraction
+  double volume_fraction() const { return volume_fraction_; };
+
+  //! Map cell volume to the nodes
+  //! \param[in] phase to map volume
+  void map_cell_volume_to_nodes(unsigned phase);
+
  private:
   //! Approximately check if a point is in a cell
   //! \param[in] point Coordinates of point
@@ -264,6 +284,10 @@ class Cell {
   //! Normal of face
   //! first-> face_id, second->vector of the normal
   std::map<unsigned, Eigen::VectorXd> face_normals_;
+  //! Free surface bool
+  bool free_surface_{false};
+  //! Volume fraction
+  double volume_fraction_{0.0};
   //! Logger
   std::unique_ptr<spdlog::logger> console_;
 };  // Cell class

include/cell.tcc

@@ -412,7 +412,7 @@ inline Eigen::Matrix<double, 2, 1> mpm::Cell<2>::local_coordinates_point(
     if (indices.size() == 3) {
       //   2 0
       //     |\
-      //     | \  
+      //     | \
       //   c |  \ b
       //     |   \
       //     |    \
@@ -844,3 +844,20 @@ template <unsigned Tdim>
 inline unsigned mpm::Cell<Tdim>::previous_mpirank() const {
   return this->previous_mpirank_;
 }
+
+//! Assign volume traction to node
+template <unsigned Tdim>
+void mpm::Cell<Tdim>::assign_volume_fraction(double volume_fraction) {
+  volume_fraction_ = volume_fraction;
+}
+
+//! Map cell volume to nodes
+template <unsigned Tdim>
+void mpm::Cell<Tdim>::map_cell_volume_to_nodes(unsigned phase) {
+  // Check if cell volume is set
+  if (volume_ == std::numeric_limits<double>::lowest()) this->compute_volume();
+
+  for (unsigned i = 0; i < nodes_.size(); ++i) {
+    nodes_[i]->update_volume(true, phase, volume_ / nnodes_);
+  }
+}

include/hdf5_particle_twophase.h

@@ -0,0 +1,49 @@
+#ifndef MPM_HDF5_TWOPHASE_H_
+#define MPM_HDF5_TWOPHASE_H_
+
+// HDF5 Particle
+#include "hdf5_particle.h"
+
+namespace mpm {
+// Define a struct of particle
+typedef struct HDF5ParticleTwoPhase : HDF5Particle {
+  // Liquid Mass
+  double liquid_mass;
+  // Liquid Velocity
+  double liquid_velocity_x, liquid_velocity_y, liquid_velocity_z;
+  // Porosity
+  double porosity;
+  // Liquid Saturation
+  double liquid_saturation;
+  // Material id
+  unsigned liquid_material_id;
+  // Number of state variables
+  unsigned nliquid_state_vars;
+  // State variables (init to zero)
+  double liquid_svars[5] = {0};
+} HDF5ParticleTwoPhase;
+
+namespace hdf5 {
+namespace particletwophase {
+const hsize_t NFIELDS = 66;
+
+const size_t dst_size = sizeof(HDF5ParticleTwoPhase);
+
+// Destination offset
+extern const size_t dst_offset[NFIELDS];
+
+// Destination size
+extern const size_t dst_sizes[NFIELDS];
+
+// Define particle field information
+extern const char* field_names[NFIELDS];
+
+// Initialize field types
+extern const hid_t field_type[NFIELDS];
+
+}  // namespace particletwophase
+}  // namespace hdf5
+
+}  // namespace mpm
+
+#endif  // MPM_HDF5_TWOPHASE_H_

include/io/io_mesh.h

@@ -63,6 +63,17 @@ class IOMesh {
   virtual std::vector<Eigen::Matrix<double, 6, 1>> read_particles_stresses(
       const std::string& particles_stresses) = 0;
 
+  //! Read particle scalar properties
+  //! \param[in] scalar_file file name with particle scalar properties
+  //! \retval Vector of particles scalar properties
+  virtual std::vector<std::tuple<mpm::Index, double>>
+      read_particles_scalar_properties(const std::string& scalar_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(
+      const std::string& _pressure_constraints_file) = 0;
+
   //! Read nodal euler angles file
   //! \param[in] nodal_euler_angles_file file name with nodal id and respective
   //! euler angles

include/io/io_mesh_ascii.h

@@ -51,6 +51,18 @@ class IOMeshAscii : public IOMesh<Tdim> {
   std::vector<Eigen::Matrix<double, 6, 1>> read_particles_stresses(
       const std::string& particles_stresses) override;
 
+  //! Read particle scalar properties
+  //! \param[in] scalar_file file name with particle scalar properties
+  //! \retval Vector of particles scalar properties
+  std::vector<std::tuple<mpm::Index, double>> read_particles_scalar_properties(
+      const std::string& scalar_file) override;
+
+  //! Read pressure constraints file
+  //! \param[in] pressure_constraints_files file name with pressure
+  //! constraints
+  std::vector<std::tuple<mpm::Index, double>> read_pressure_constraints(
+      const std::string& pressure_constraints_file) override;
+
   //! Read nodal euler angles file
   //! \param[in] nodal_euler_angles_file file name with nodal id and respective
   //! euler angles

include/io/io_mesh_ascii.tcc

@@ -246,6 +246,93 @@ std::vector<Eigen::Matrix<double, 6, 1>>
   return stresses;
 }
 
+//! Return particles scalar properties
+template <unsigned Tdim>
+std::vector<std::tuple<mpm::Index, double>>
+    mpm::IOMeshAscii<Tdim>::read_particles_scalar_properties(
+        const std::string& scalar_file) {
+
+  // Particles scalar properties
+  std::vector<std::tuple<mpm::Index, double>> scalar_properties;
+
+  // input file stream
+  std::fstream file;
+  file.open(scalar_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 != "")) {
+          while (istream.good()) {
+            // ID
+            mpm::Index id;
+            // Scalar
+            double scalar;
+            // Read stream
+            istream >> id >> scalar;
+            scalar_properties.emplace_back(std::make_tuple(id, scalar));
+          }
+        }
+      }
+      file.close();
+    }
+  } catch (std::exception& exception) {
+    console_->error("Read particle {} #{}: {}\n", __FILE__, __LINE__,
+                    exception.what());
+    file.close();
+  }
+  return scalar_properties;
+}
+
+//! Read pressure constraints file
+template <unsigned Tdim>
+std::vector<std::tuple<mpm::Index, double>>
+    mpm::IOMeshAscii<Tdim>::read_pressure_constraints(
+        const std::string& pressure_constraints_file) {
+  // Particle pressure constraints
+  std::vector<std::tuple<mpm::Index, double>> constraints;
+  constraints.clear();
+
+  // input file stream
+  std::fstream file;
+  file.open(pressure_constraints_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 != "")) {
+          while (istream.good()) {
+            // ID
+            mpm::Index id;
+            // Pressure
+            double pressure;
+            // Read stream
+            istream >> id >> pressure;
+            constraints.emplace_back(std::make_tuple(id, pressure));
+          }
+        }
+      }
+    }
+    file.close();
+  } catch (std::exception& exception) {
+    console_->error("Read pressure constraints: {}", exception.what());
+    file.close();
+  }
+  return constraints;
+}
+
 //! Return euler angles of nodes
 template <unsigned Tdim>
 std::map<mpm::Index, Eigen::Matrix<double, Tdim, 1>>

include/io/logger.h

@@ -36,6 +36,9 @@ struct Logger {
   // Create a logger for MPM Explicit
   static const std::shared_ptr<spdlog::logger> mpm_explicit_logger;
 
+  // Create a logger for MPM Explicit Two Phase
+  static const std::shared_ptr<spdlog::logger> mpm_explicit_two_phase_logger;
+
   // Create a logger for MPM Explicit USF
   static const std::shared_ptr<spdlog::logger> mpm_explicit_usf_logger;
 

include/loads_bcs/constraints.h

@@ -48,6 +48,21 @@ class Constraints {
       const std::vector<std::tuple<mpm::Index, unsigned, int, double>>&
           friction_constraints);
 
+  //! Assign nodal velocity constraints
+  //! \param[in] mfunction Math function if defined
+  //! \param[in] setid Node set id
+  //! \param[in] phase Index corresponding to the phase
+  //! \param[in] pconstraint Pressure constraint at node
+  bool assign_nodal_pressure_constraint(
+      const std::shared_ptr<FunctionBase>& mfunction, int set_id,
+      const unsigned phase, const unsigned pconstraint);
+
+  //! Assign nodal pressure constraints to nodes
+  //! \param[in] pressure_constraints Constraint at node, pressure
+  bool assign_nodal_pressure_constraints(
+      const unsigned phase,
+      const std::vector<std::tuple<mpm::Index, double>>& pressure_constraints);
+
  private:
   //! Mesh object
   std::shared_ptr<mpm::Mesh<Tdim>> mesh_;

include/loads_bcs/constraints.tcc

@@ -105,3 +105,53 @@ bool mpm::Constraints<Tdim>::assign_nodal_friction_constraints(
   }
   return status;
 }
+
+//! Assign nodal pressure constraints
+template <unsigned Tdim>
+bool mpm::Constraints<Tdim>::assign_nodal_pressure_constraint(
+    const std::shared_ptr<FunctionBase>& mfunction, int set_id,
+    const unsigned phase, const unsigned pconstraint) {
+  bool status = true;
+  try {
+    auto nset = mesh_->nodes(set_id);
+    if (nset.size() == 0)
+      throw std::runtime_error(
+          "Node set is empty for assignment of pressure constraints");
+
+#pragma omp parallel for schedule(runtime)
+    for (auto nitr = nset.cbegin(); nitr != nset.cend(); ++nitr) {
+      if (!(*nitr)->assign_pressure_constraint(phase, pconstraint, mfunction))
+        throw std::runtime_error("Setting pressure constraint failed");
+    }
+  } catch (std::exception& exception) {
+    console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what());
+    status = false;
+  }
+  return status;
+}
+
+//! Assign nodal pressure constraints to nodes
+template <unsigned Tdim>
+bool mpm::Constraints<Tdim>::assign_nodal_pressure_constraints(
+    const unsigned phase,
+    const std::vector<std::tuple<mpm::Index, double>>& pressure_constraints) {
+  bool status = true;
+  try {
+    for (const auto& pressure_constraint : pressure_constraints) {
+      // Node id
+      mpm::Index nid = std::get<0>(pressure_constraint);
+      // Pressure
+      double pressure = std::get<1>(pressure_constraint);
+
+      // Apply constraint
+      if (!mesh_->node(nid)->assign_pressure_constraint(phase, pressure,
+                                                        nullptr))
+        throw std::runtime_error(
+            "Nodal pressure constraints assignment failed");
+    }
+  } catch (std::exception& exception) {
+    console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what());
+    status = false;
+  }
+  return status;
+}