Levelset + Barrier Method

include/contacts/contact.h

@@ -20,6 +20,14 @@ class Contact {
   //! Compute contact forces
   virtual inline void compute_contact_forces(){};
 
+  //! Compute contact forces
+  //! \param[in] dt Analysis time step
+  //! \param[in] levelset_damping Levelset damping factor
+  //! \param[in] levelset_pic Particle in cell method bool for contact velocity
+  virtual inline void compute_contact_forces(const double levelset_damping,
+                                             const bool levelset_pic,
+                                             double dt){};
+
  protected:
   //! Mesh object
   std::shared_ptr<mpm::Mesh<Tdim>> mesh_;

include/contacts/contact_friction.h

@@ -15,10 +15,17 @@ class ContactFriction : public Contact<Tdim> {
   ContactFriction(const std::shared_ptr<mpm::Mesh<Tdim>>& mesh);
 
   //! Intialize
-  virtual inline void initialise() override;
+  void initialise() override;
 
   //! Compute contact forces
-  virtual inline void compute_contact_forces() override;
+  void compute_contact_forces() override;
+
+  //! Compute contact forces
+  //! \param[in] dt Analysis time step
+  //! \param[in] levelset_damping Levelset damping factor
+  //! \param[in] levelset_pic Particle in cell method bool for contact velocity
+  void compute_contact_forces(const double levelset_damping,
+                              const bool levelset_pic, double dt) override;
 
  protected:
   //! Mesh object

include/contacts/contact_friction.tcc

@@ -6,7 +6,7 @@
 
 //! Initialize nodal properties
 template <unsigned Tdim>
-inline void mpm::ContactFriction<Tdim>::initialise() {
+void mpm::ContactFriction<Tdim>::initialise() {
   // Initialise nodal properties
   mesh_->initialise_nodal_properties();
 
@@ -18,7 +18,43 @@
 
 //! Compute contact forces
 template <unsigned Tdim>
-inline void mpm::ContactFriction<Tdim>::compute_contact_forces() {
+void mpm::ContactFriction<Tdim>::compute_contact_forces() {
+
+  // Map multimaterial properties from particles to nodes
+  mesh_->iterate_over_particles(std::bind(
+      &mpm::ParticleBase<Tdim>::map_multimaterial_mass_momentum_to_nodes,
+      std::placeholders::_1));
+
+  // Map multimaterial displacements from particles to nodes
+  mesh_->iterate_over_particles(std::bind(
+      &mpm::ParticleBase<Tdim>::map_multimaterial_displacements_to_nodes,
+      std::placeholders::_1));
+
+  // Map multimaterial domain gradients from particles to nodes
+  mesh_->iterate_over_particles(std::bind(
+      &mpm::ParticleBase<Tdim>::map_multimaterial_domain_gradients_to_nodes,
+      std::placeholders::_1));
+
+  // Compute multimaterial change in momentum
+  mesh_->iterate_over_nodes(
+      std::bind(&mpm::NodeBase<Tdim>::compute_multimaterial_change_in_momentum,
+                std::placeholders::_1));
+
+  // Compute multimaterial separation vector
+  mesh_->iterate_over_nodes(
+      std::bind(&mpm::NodeBase<Tdim>::compute_multimaterial_separation_vector,
+                std::placeholders::_1));
+
+  // Compute multimaterial normal unit vector
+  mesh_->iterate_over_nodes(
+      std::bind(&mpm::NodeBase<Tdim>::compute_multimaterial_normal_unit_vector,
+                std::placeholders::_1));
+}
+
+//! Compute contact forces
+template <unsigned Tdim>
+void mpm::ContactFriction<Tdim>::compute_contact_forces(
+    const double levelset_damping, const bool levelset_pic, double dt) {
 
   // Map multimaterial properties from particles to nodes
   mesh_->iterate_over_particles(std::bind(

include/contacts/contact_levelset.h

@@ -0,0 +1,41 @@
+#ifndef MPM_CONTACT_LEVELSET_H_
+#define MPM_CONTACT_LEVELSET_H_
+
+#include "contact.h"
+#include "mesh_levelset.h"
+#include "particle_levelset.h"
+
+namespace mpm {
+
+//! ContactLevelset class
+//! \brief ContactLevelset base class
+//! \tparam Tdim Dimension
+template <unsigned Tdim>
+class ContactLevelset : public Contact<Tdim> {
+ public:
+  //! Default constructor with mesh class
+  ContactLevelset(const std::shared_ptr<mpm::Mesh<Tdim>>& mesh);
+
+  //! Intialize
+  void initialise() override;
+
+  //! Compute contact forces
+  //! \param[in] dt Analysis time step
+  //! \param[in] levelset_damping Levelset damping factor
+  //! \param[in] levelset_pic Particle in cell method bool for contact velocity
+  void compute_contact_forces(const double levelset_damping,
+                              const bool levelset_pic, double dt) override;
+
+  //! Mesh levelset object
+  std::shared_ptr<mpm::MeshLevelset<Tdim>> mesh_levelset_;
+
+ protected:
+  //! Mesh object
+  using mpm::Contact<Tdim>::mesh_;
+
+};  // ContactLevelset class
+}  // namespace mpm
+
+#include "contact_levelset.tcc"
+
+#endif  // MPM_CONTACT_LEVELSET_H_

include/contacts/contact_levelset.tcc

@@ -0,0 +1,26 @@
+//! Constructor of contact with mesh
+template <unsigned Tdim>
+mpm::ContactLevelset<Tdim>::ContactLevelset(
+    const std::shared_ptr<mpm::Mesh<Tdim>>& mesh)
+    : mpm::Contact<Tdim>(mesh) {
+  // Assign mesh
+  mesh_ = mesh;
+}
+
+//! Initialize nodal properties
+template <unsigned Tdim>
+void mpm::ContactLevelset<Tdim>::initialise() {
+  // Initialise nodal properties
+  mesh_->initialise_nodal_properties();
+}
+
+//! Compute contact forces
+template <unsigned Tdim>
+void mpm::ContactLevelset<Tdim>::compute_contact_forces(
+    const double levelset_damping, const bool levelset_pic, double dt) {
+
+  // Compute and map contact forces to nodes
+  mesh_->iterate_over_particles(
+      std::bind(&mpm::ParticleBase<Tdim>::map_particle_contact_force_to_nodes,
+                std::placeholders::_1, levelset_damping, levelset_pic, dt));
+}

include/io/io_mesh.h

@@ -64,11 +64,11 @@ class IOMesh {
   virtual std::vector<Eigen::Matrix<double, 6, 1>> read_particles_stresses(
       const std::string& particles_stresses) = 0;
 
-  //! Read particle scalar properties
+  //! Read scalar properties for particles or nodes
   //! \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;
+  //! \retval Vector of scalar properties for particles or nodes
+  virtual std::vector<std::tuple<mpm::Index, double>> read_scalar_properties(
+      const std::string& scalar_file) = 0;
 
   //! Read pressure constraints file
   //! \param[in] pressure_constraints_files file name with pressure constraints
@@ -123,6 +123,12 @@ class IOMesh {
       read_adhesion_constraints(
           const std::string& adhesion_constraints_file) = 0;
 
+  //! Read levelset file
+  //! \param[in] levelset_input_file file name with levelset values
+  virtual std::vector<
+      std::tuple<mpm::Index, double, double, double, double, double>>
+      read_levelset_input(const std::string& levelset_input_file) = 0;
+
   //! Read forces file
   //! \param[in] forces_file file name with nodal concentrated force
   virtual std::vector<std::tuple<mpm::Index, unsigned, double>> read_forces(

include/io/io_mesh_ascii.h

@@ -51,10 +51,10 @@ 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(
+  //! Read scalar properties for particles or nodes
+  //! \param[in] scalar_file file name with particle or node scalar properties
+  //! \retval Vector of scalar properties for particles or nodes
+  std::vector<std::tuple<mpm::Index, double>> read_scalar_properties(
       const std::string& scalar_file) override;
 
   //! Read pressure constraints file
@@ -110,6 +110,11 @@ class IOMeshAscii : public IOMesh<Tdim> {
       read_adhesion_constraints(
           const std::string& adhesion_constraints_file) override;
 
+  //! Read levelset file
+  //! \param[in] levelset_input_file file name with levelset values
+  std::vector<std::tuple<mpm::Index, double, double, double, double, double>>
+      read_levelset_input(const std::string& levelset_input_file) override;
+
   //! Read traction file
   //! \param[in] forces_file file name with nodal concentrated force
   std::vector<std::tuple<mpm::Index, unsigned, double>> read_forces(

include/io/io_mesh_ascii.tcc

@@ -254,13 +254,13 @@
   return stresses;
 }
 
-//! Return particles scalar properties
+//! Return scalar properties for particles or nodes
 template <unsigned Tdim>
 std::vector<std::tuple<mpm::Index, double>>
-    mpm::IOMeshAscii<Tdim>::read_particles_scalar_properties(
+    mpm::IOMeshAscii<Tdim>::read_scalar_properties(
         const std::string& scalar_file) {
 
-  // Particles scalar properties
+  // Scalar properties for particles or nodes
   std::vector<std::tuple<mpm::Index, double>> scalar_properties;
 
   // input file stream
@@ -293,7 +293,7 @@
     }
     file.close();
   } catch (std::exception& exception) {
-    console_->error("Read particle {} #{}: {}\n", __FILE__, __LINE__,
+    console_->error("Read particle/node {} #{}: {}\n", __FILE__, __LINE__,
                     exception.what());
     file.close();
   }
@@ -641,7 +641,7 @@
   return constraints;
 }
 
-//! Return friction constraints of particles
+//! Return friction constraints
 template <unsigned Tdim>
 std::vector<std::tuple<mpm::Index, unsigned, int, double>>
     mpm::IOMeshAscii<Tdim>::read_friction_constraints(
@@ -748,6 +748,64 @@
   return constraints;
 }
 
+//! Return nodal levelset information
+template <unsigned Tdim>
+std::vector<std::tuple<mpm::Index, double, double, double, double, double>>
+    mpm::IOMeshAscii<Tdim>::read_levelset_input(
+        const std::string& levelset_input_file) {
+
+  // Nodal levelset information
+  std::vector<std::tuple<mpm::Index, double, double, double, double, double>>
+      levelset_inputs;
+  levelset_inputs.clear();
+
+  // input file stream
+  std::fstream file;
+  file.open(levelset_input_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 != "")) {
+          // ID
+          mpm::Index id;
+          // Levelset value
+          double levelset;
+          // Friction
+          double levelset_mu;
+          // Adhesion coefficient
+          double levelset_alpha;
+          // Barrier stiffness
+          double barrier_stiffness;
+          // Slip threshold
+          double slip_threshold;
+          while (istream.good()) {
+            // Read stream
+            istream >> id >> levelset >> levelset_mu >> levelset_alpha >>
+                barrier_stiffness >> slip_threshold;
+            levelset_inputs.emplace_back(
+                std::make_tuple(id, levelset, levelset_mu, levelset_alpha,
+                                barrier_stiffness, slip_threshold));
+          }
+        }
+      }
+    } else {
+      throw std::runtime_error("File not open or not good!");
+    }
+    file.close();
+  } catch (std::exception& exception) {
+    console_->error("Read levelset inputs: {}", exception.what());
+    file.close();
+  }
+  return levelset_inputs;
+}
+
 //! Return particles force
 template <unsigned Tdim>
 std::vector<std::tuple<mpm::Index, unsigned, double>>

include/mesh/mesh.h

@@ -164,7 +164,7 @@
   //! \param[in] cells Node ids of cells
   //! \param[in] check_duplicates Parameter to check duplicates
   //! \retval status Create cells status
-  bool create_cells(mpm::Index gnid,
+  bool create_cells(mpm::Index gcid,
                     const std::shared_ptr<mpm::Element<Tdim>>& element,
                     const std::vector<std::vector<mpm::Index>>& cells,
                     bool check_duplicates = true);
@@ -263,7 +263,7 @@
   //! Locate particles in a cell
   //! Iterate over all cells in a mesh to find the cell in which particles
   //! are located.
-  //! \retval particles Particles which cannot be located in the mesh
+  //! \retval particles which cannot be located in the mesh
   std::vector<std::shared_ptr<mpm::ParticleBase<Tdim>>> locate_particles_mesh();
 
   //! Iterate over particles
@@ -475,7 +475,7 @@
   //! Create map of vector of particles in sets
   //! \param[in] map of particles ids in sets
   //! \param[in] check_duplicates Parameter to check duplicates
-  //! \retval status Status of  create particle sets
+  //! \retval status Status of create particle sets
   bool create_particle_sets(
       const tsl::robin_map<mpm::Index, std::vector<mpm::Index>>& particle_sets,
       bool check_duplicates);
@@ -518,11 +518,34 @@
   void inject_particles(double current_time);
 
   // Create the nodal properties' map
-  void create_nodal_properties();
+  virtual void create_nodal_properties();
 
   // Initialise the nodal properties' map
   void initialise_nodal_properties();
 
+  /**
+   * \defgroup Levelset Functions
+   */
+  /**@{*/
+
+  //! Assign nodal levelset values
+  //! \ingroup Levelset
+  //! \param[in] levelset Levelset value at the particle
+  //! \param[in] levelset_mu Levelset friction
+  //! \param[in] levelset_alpha Levelset adhesion coefficient
+  //! \param[in] barrier_stiffness Barrier stiffness
+  //! \param[in] slip_threshold Slip threshold
+  virtual bool assign_nodal_levelset_values(
+      const std::vector<std::tuple<mpm::Index, double, double, double, double,
+                                   double>>& levelset_input_file) {
+    throw std::runtime_error(
+        "Calling the base class function (assign_nodal_levelset_values) in "
+        "Mesh:: illegal operation!");
+    return false;
+  };
+
+  /**@}*/
+
   /**
    * \defgroup MultiPhase Functions dealing with multi-phase MPM
    */
@@ -673,10 +696,11 @@
                            const Json& generator, unsigned pset_id);
 
   // Locate a particle in mesh cells
+  //! \param[in] particle of interest
   bool locate_particle_cells(
       const std::shared_ptr<mpm::ParticleBase<Tdim>>& particle);
 
- private:
+ protected:
   //! mesh id
   unsigned id_{std::numeric_limits<unsigned>::max()};
   //! Isoparametric mesh