🔨modify the solver scheme

include/mesh.h

@@ -462,14 +462,11 @@ class Mesh {
   // Initialise the nodal properties' map
   void initialise_nodal_properties();
 
-  //! Set particles lsm values
-  void assign_particle_levelset(const std::vector<double>& phi_list);
-
   // Create the nodal properties' map for discontinuity
   void create_nodal_properties_discontinuity();
 
   //! Initialise discontinuities
-  //! \param[in] discontinuities List of discontinuities 
+  //! \param[in] discontinuities List of discontinuities
   void initialise_discontinuities(
       const std::map<unsigned, std::shared_ptr<mpm::DiscontinuityBase<Tdim>>>&
           discontinuities) {
@@ -489,6 +486,9 @@ class Mesh {
   //! compute the normal vector of enriched nodes at the discontinuity
   void compute_normal_vector_discontinuity();
 
+  //! Initialise the level set function values
+  void initialise_levelset_discontinuity();
+
  private:
   // Read particles from file
   //! \param[in] pset_id Set ID of the particles

include/mesh.tcc

@@ -1975,8 +1975,7 @@ void mpm::Mesh<Tdim>::create_nodal_properties_discontinuity() {
     nodal_properties_->create_property("external_force_enrich", nrows, 1);
     nodal_properties_->create_property("normal_unit_vectors_discontinuity",
                                        nrows, 1);
-    nodal_properties_->create_property("friction_coef",
-                                       nodes_.size(), 1);
+    nodal_properties_->create_property("friction_coef", nodes_.size(), 1);
     // Iterate over all nodes to initialise the property handle in each node
     // and assign its node id as the prop id in the nodal property data pool
     for (auto nitr = nodes_.cbegin(); nitr != nodes_.cend(); ++nitr)
@@ -1994,13 +1993,6 @@ void mpm::Mesh<Tdim>::initialise_nodal_properties() {
   nodal_properties_->initialise_nodal_properties();
 }
 
-//! Set particles lsm values
-template <unsigned Tdim>
-void mpm::Mesh<Tdim>::assign_particle_levelset(const std::vector<double>& phi_list) {
-  for (mpm::Index i = 0; i < nparticles(); ++i)
-    particles_[i]->assign_levelsetphi(phi_list[i]);
-}
-
 //! Locate points in a cell
 template <unsigned Tdim>
 void mpm::Mesh<Tdim>::locate_discontinuity() {
@@ -2045,3 +2037,17 @@ void mpm::Mesh<Tdim>::compute_normal_vector_discontinuity() {
                                        normal, Tdim);
   }
 }
+
+// Initialise level set values particles
+template <unsigned Tdim>
+void mpm::Mesh<Tdim>::initialise_levelset_discontinuity() {
+
+  double phi_particle;
+  for (unsigned i = 0; i < discontinuities_.size(); ++i) {
+    for (mpm::Index j = 0; j < nparticles(); ++j) {
+      discontinuities_[i]->compute_levelset(particles_[j]->coordinates(),
+                                            phi_particle);
+      particles_[j]->assign_levelsetphi(phi_particle);
+    }
+  }
+}

include/node.h

@@ -279,7 +279,7 @@ class Node : public NodeBase<Tdim> {
   void assign_discontinuity_enrich(bool discontinuity) {
     discontinuity_enrich_ = discontinuity;
   };
-  
+
   //! Return whether the node is enriched
   bool discontinuity_enrich() const { return discontinuity_enrich_; };
 

include/node.tcc

@@ -34,7 +34,7 @@ void mpm::Node<Tdim, Tdof, Tnphases>::initialise() noexcept {
   acceleration_.setZero();
   status_ = false;
   material_ids_.clear();
-  //mark: need to check
+  // mark: need to check
   discontinuity_enrich_ = true;
 }
 

include/node_xmpm.tcc

@@ -121,8 +121,8 @@ void mpm::Node<Tdim, Tdof, Tnphases>::self_contact_discontinuity(
   //  single phase for solid
   unsigned phase = 0;
   const double tolerance = 1.0E-15;
-  
-  //obtain the enriched values of enriched nodes
+
+  // obtain the enriched values of enriched nodes
   Eigen::Matrix<double, 1, 1> mass_enrich =
       property_handle_->property("mass_enrich", discontinuity_prop_id_, 0, 1);
   Eigen::Matrix<double, Tdim, 1> momenta_enrich = property_handle_->property(
@@ -156,7 +156,8 @@ void mpm::Node<Tdim, Tdof, Tnphases>::self_contact_discontinuity(
 
   // friction_coef < 0: move together without slide
   // need to be done
-  double friction_coef = property_handle_->property("friction_coef", discontinuity_prop_id_, 0, 1)(0,0);
+  double friction_coef = property_handle_->property(
+      "friction_coef", discontinuity_prop_id_, 0, 1)(0, 0);
 
   if (friction_coef < 0) {
     property_handle_->update_property("momenta_enrich", discontinuity_prop_id_,
@@ -165,25 +166,27 @@ void mpm::Node<Tdim, Tdof, Tnphases>::self_contact_discontinuity(
                                       discontinuity_prop_id_, 0,
                                       force_contact.col(phase), Tdim);
   } else {
-    //the contact momentum, force value for sticking contact at normal direction
+    // the contact momentum, force value for sticking contact at normal
+    // direction
     double momentum_contact_norm =
         momentum_contact.col(phase).dot(normal_vector);
     double force_contact_norm = momentum_contact_norm / dt;
 
     // the maximum friction contact force
     double max_friction_force = friction_coef * abs(force_contact_norm);
 
-  // the contact momentum, force vector for sticking contact at tangential direction
+    // the contact momentum, force vector for sticking contact at tangential
+    // direction
     auto momentum_tangential =
         momentum_contact.col(phase) - momentum_contact_norm * normal_vector;
     auto force_tangential = momentum_tangential / dt;
 
-    //the friction force magnitude
+    // the friction force magnitude
     double force_tangential_value = force_tangential.norm();
 
     double force_friction = force_tangential_value < max_friction_force
-                                  ? force_tangential_value
-                                  : max_friction_force;
+                                ? force_tangential_value
+                                : max_friction_force;
 
     // adjust the momentum and force
     property_handle_->update_property(

include/particles/particle.tcc

@@ -977,8 +977,8 @@ std::vector<uint8_t> mpm::Particle<Tdim>::serialize() {
   unsigned nmaterials = material_id_.size();
   MPI_Pack(&nmaterials, 1, MPI_UNSIGNED, data_ptr, data.size(), &position,
            MPI_COMM_WORLD);
-  MPI_Pack(&material_id_[mpm::ParticlePhase::Solid], 1, MPI_UNSIGNED, data_ptr, data.size(), &position,
-           MPI_COMM_WORLD);
+  MPI_Pack(&material_id_[mpm::ParticlePhase::Solid], 1, MPI_UNSIGNED, data_ptr,
+           data.size(), &position, MPI_COMM_WORLD);
 
   // ID
   MPI_Pack(&id_, 1, MPI_UNSIGNED_LONG_LONG, data_ptr, data.size(), &position,

include/particles/particle_xmpm.h

@@ -68,6 +68,7 @@ class ParticleXMPM : public Particle<Tdim> {
   //! \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;
+
  private:
   //! Assign the level set function values
   //! \param[in] phivalue The level set values
@@ -142,8 +143,14 @@ class ParticleXMPM : public Particle<Tdim> {
   using Particle<Tdim>::dn_dx_centroid_;
   //! Logger
   using Particle<Tdim>::console_;
+  //! Map of scalar properties
+  using Particle<Tdim>::scalar_properties_;
   //! Map of vector properties
-  using Particle<Tdim>::properties_;
+  using Particle<Tdim>::vector_properties_;
+  //! Map of tensor properties
+  using Particle<Tdim>::tensor_properties_;
+  //! Pack size
+  using Particle<Tdim>::pack_size_;
 
  private:
   //! level set value： phi for discontinuity

include/solvers/xmpm_explicit.h

@@ -28,9 +28,8 @@ class XMPMExplicit : public MPMBase<Tdim> {
   void compute_stress_strain(unsigned phase);
 
   //! Initialise discontinuities
-  bool initialise_discontinuities();
-  //! Initialise the level set function values
-  bool initialise_levelset();
+  void initialise_discontinuities();
+
   // return the number of discontinuities
   mpm::Index ndiscontinuities() { return discontinuities_.size(); };
 
@@ -55,6 +54,12 @@ class XMPMExplicit : public MPMBase<Tdim> {
   using mpm::MPMBase<Tdim>::post_process_;
   //! Logger
   using mpm::MPMBase<Tdim>::console_;
+  //! MPM Scheme
+  using mpm::MPMBase<Tdim>::mpm_scheme_;
+  //! Stress update method
+  using mpm::MPMBase<Tdim>::stress_update_;
+  //! Interface scheme
+  using mpm::MPMBase<Tdim>::contact_;
 
 #ifdef USE_GRAPH_PARTITIONING
   //! Graph