[Solver] Two-Phase One-Point Explicit

.gitignore

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

CMakeLists.txt

@@ -71,6 +71,7 @@ if (HDF5_FOUND)
   include_directories(${HDF5_INCLUDE_DIRS})
   link_libraries(${HDF5_LIBRARIES} ${HDF5_HL_LIBRARIES} ${HDF5_CXX_HL_LIBRARIES})
   add_definitions(${HDF5_DEFINITIONS})
+  SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-invalid-offsetof")
 endif()
 
 # OpenMP
@@ -142,8 +143,10 @@ include_directories(BEFORE
   ${mpm_SOURCE_DIR}/include/loads_bcs/
   ${mpm_SOURCE_DIR}/include/materials/
   ${mpm_SOURCE_DIR}/include/particles/
+  ${mpm_SOURCE_DIR}/include/particles/pod_particles
   ${mpm_SOURCE_DIR}/include/solvers/
   ${mpm_SOURCE_DIR}/include/solvers/mpm_scheme/
+  ${mpm_SOURCE_DIR}/include/utilities/
   ${mpm_SOURCE_DIR}/external/
   ${mpm_SOURCE_DIR}/tests/include/
 )
@@ -157,7 +160,8 @@ SET(mpm_src
   ${mpm_SOURCE_DIR}/src/functions/linear_function.cc
   ${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/pod_particle.cc
+  ${mpm_SOURCE_DIR}/src/pod_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
@@ -191,6 +195,7 @@ if(MPM_BUILD_TESTING)
     ${mpm_SOURCE_DIR}/tests/elements/triangle_quadrature_test.cc
     ${mpm_SOURCE_DIR}/tests/functions/linear_function_test.cc
     ${mpm_SOURCE_DIR}/tests/functions/sin_function_test.cc
+    ${mpm_SOURCE_DIR}/tests/functions/radial_basis_function_test.cc
     ${mpm_SOURCE_DIR}/tests/graph_test.cc
     ${mpm_SOURCE_DIR}/tests/interface_test.cc
     ${mpm_SOURCE_DIR}/tests/io/io_mesh_ascii_test.cc
@@ -205,25 +210,33 @@ if(MPM_BUILD_TESTING)
     ${mpm_SOURCE_DIR}/tests/materials/newtonian_test.cc
     ${mpm_SOURCE_DIR}/tests/materials/norsand_test.cc
     ${mpm_SOURCE_DIR}/tests/materials/material_utility_test.cc
+    ${mpm_SOURCE_DIR}/tests/mesh_free_surface_test.cc
     ${mpm_SOURCE_DIR}/tests/mesh_neighbours_test.cc
     ${mpm_SOURCE_DIR}/tests/mesh_test_2d.cc
     ${mpm_SOURCE_DIR}/tests/mesh_test_3d.cc
     ${mpm_SOURCE_DIR}/tests/mpi_transfer_particle_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/particles/particle_serialize_deserialize_test.cc
+    ${mpm_SOURCE_DIR}/tests/particles/particle_serialize_deserialize_twophase_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
     ${mpm_SOURCE_DIR}/tests/solvers/mpm_scheme_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_serialize_deserialize_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/point_in_cell_test.cc
   )
   add_executable(mpmtest ${mpm_src} ${test_src})
   add_test(NAME mpmtest COMMAND $<TARGET_FILE:mpmtest>)

include/cell.h

@@ -217,6 +217,50 @@ class Cell {
   //! Return previous mpi rank
   unsigned previous_mpirank() const;
 
+  /**
+   * \defgroup MultiPhase Functions dealing with multi-phase MPM
+   */
+  /**@{*/
+
+  //! Assign solving status
+  //! \ingroup MultiPhase
+  //! \param[in] status Cell solving status for parallel free-surface detection
+  void assign_solving_status(bool status) { solving_status_ = status; }
+
+  //! Return solving status of a cell: active (if a particle is present in all
+  //! MPI rank)
+  //! \ingroup MultiPhase
+  bool solving_status() const { return solving_status_; }
+
+  //! Assign free surface
+  //! \ingroup MultiPhase
+  //! \param[in] free_surface boolean indicating free surface cell
+  void assign_free_surface(bool free_surface) { free_surface_ = free_surface; };
+
+  //! Return free surface bool
+  //! \ingroup MultiPhase
+  //! \retval free_surface_ indicating free surface cell
+  bool free_surface() const { return free_surface_; };
+
+  //! Assign volume traction to node
+  //! \ingroup MultiPhase
+  //! \param[in] volume_fraction cell volume fraction
+  void assign_volume_fraction(double volume_fraction) {
+    volume_fraction_ = volume_fraction;
+  };
+
+  //! Return cell volume fraction
+  //! \ingroup MultiPhase
+  //! \retval volume_fraction_ cell volume fraction
+  double volume_fraction() const { return volume_fraction_; };
+
+  //! Map cell volume to the nodes
+  //! \ingroup MultiPhase
+  //! \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,11 +308,17 @@ class Cell {
   //! Normal of face
   //! first-> face_id, second->vector of the normal
   std::map<unsigned, Eigen::VectorXd> face_normals_;
+  //! Solving status
+  bool solving_status_{false};
+  //! Free surface bool
+  bool free_surface_{false};
+  //! Volume fraction
+  double volume_fraction_{0.0};
   //! Logger
   std::unique_ptr<spdlog::logger> console_;
 };  // Cell class
 }  // namespace mpm
 
 #include "cell.tcc"
 
-#endif  // MPM_CELL_H_
+#endif  // MPM_CELL_H_

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,17 @@ template <unsigned Tdim>
 inline unsigned mpm::Cell<Tdim>::previous_mpirank() const {
   return this->previous_mpirank_;
 }
+
+//! Map cell volume to nodes
+template <unsigned Tdim>
+void mpm::Cell<Tdim>::map_cell_volume_to_nodes(unsigned phase) {
+  if (this->status()) {
+    // 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/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