Kelvin-Voigt Boundary Implementation for #692

include/loads_bcs/absorbing_constraint.h

@@ -0,0 +1,68 @@
+#ifndef MPM_ABSORBING_CONSTRAINT_H_
+#define MPM_ABSORBING_CONSTRAINT_H_
+
+namespace mpm {
+
+//! AbsorbingConstraint class to store friction constraint on a set
+//! \brief AbsorbingConstraint class to store a constraint on a set
+//! \details AbsorbingConstraint stores the constraint as a static value
+class AbsorbingConstraint {
+ public:
+  // Constructor
+  //! \param[in] setid  set id
+  //! \param[in] dir Normal direction of boundary application
+  //! \param[in] delta Virtual viscous layer thickness
+  //! \param[in] h_min Cell height
+  //! \param[in] a Equation constant
+  //! \param[in] b Equation constant
+  //! \param[in] position Nodal position along boundary
+  AbsorbingConstraint(int setid, unsigned dir, double delta, double h_min,
+                      double a = 1., double b = 1.,
+                      std::string position = "null")
+      : setid_{setid},
+        dir_{dir},
+        delta_{delta},
+        h_min_{h_min},
+        a_{a},
+        b_{b},
+        position_{position} {};
+
+  // Set id
+  int setid() const { return setid_; }
+
+  // Direction
+  unsigned dir() const { return dir_; }
+
+  // Virtual viscous layer thickness
+  double delta() const { return delta_; }
+
+  // Cell height
+  double h_min() const { return h_min_; }
+
+  // Return a
+  double a() const { return a_; }
+
+  // Return b
+  double b() const { return b_; }
+
+  // Return position
+  std::string position() const { return position_; }
+
+ private:
+  // ID
+  int setid_;
+  // Direction
+  unsigned dir_;
+  // Delta
+  double delta_;
+  // Cell height
+  double h_min_;
+  // a
+  double a_;
+  // b
+  double b_;
+  // Node position
+  std::string position_;
+};
+}  // namespace mpm
+#endif  // MPM_ABSORBING_CONSTRAINT_H_

include/loads_bcs/constraints.h

@@ -3,6 +3,7 @@
 
 #include <memory>
 
+#include "absorbing_constraint.h"
 #include "friction_constraint.h"
 #include "logger.h"
 #include "mesh.h"
@@ -48,6 +49,22 @@ class Constraints {
       const std::vector<std::tuple<mpm::Index, unsigned, int, double>>&
           friction_constraints);
 
+  //! Assign nodal absorbing constraints
+  //! \param[in] setid Node set id
+  //! \param[in] absorbing_constraints Constraint at node, dir, delta, h_min, a,
+  //! b
+  bool assign_nodal_absorbing_constraint(
+      int nset_id,
+      const std::shared_ptr<mpm::AbsorbingConstraint>& absorbing_constraints);
+
+  //! Assign absorbing constraints to nodes
+  //! \param[in] absorbing_constraints Constraint at node, dir, delta, h_min, a,
+  //! b, and position
+  bool assign_nodal_absorbing_constraints(
+      const std::vector<std::tuple<mpm::Index, unsigned, double, double, double,
+                                   double, std::string>>&
+          absorbing_constraints);
+
  private:
   //! Mesh object
   std::shared_ptr<mpm::Mesh<Tdim>> mesh_;

include/loads_bcs/constraints.tcc

@@ -105,3 +105,81 @@ bool mpm::Constraints<Tdim>::assign_nodal_friction_constraints(
   }
   return status;
 }
+
+//! Apply absorbing constraints to nodes
+template <unsigned Tdim>
+bool mpm::Constraints<Tdim>::assign_nodal_absorbing_constraint(
+    int nset_id,
+    const std::shared_ptr<mpm::AbsorbingConstraint>& absorbing_constraint) {
+  bool status = true;
+  try {
+    // int set_id = absorbing_constraint->setid();
+    int set_id = nset_id;
+    auto nset = mesh_->nodes(set_id);
+    if (nset.size() == 0)
+      throw std::runtime_error(
+          "Node set is empty for application of absorbing constraints");
+    unsigned dir = absorbing_constraint->dir();
+    double delta = absorbing_constraint->delta();
+    double h_min = absorbing_constraint->h_min();
+    double a = absorbing_constraint->a();
+    double b = absorbing_constraint->b();
+    std::string position = absorbing_constraint->position();
+    if (delta >= h_min / (2 * a) and delta >= h_min / (2 * b))
+      for (auto nitr = nset.cbegin(); nitr != nset.cend(); ++nitr) {
+        if (!(*nitr)->apply_absorbing_constraint(dir, delta, h_min, a, b,
+                                                 position))
+          throw std::runtime_error(
+              "Failed to apply absorbing constraint at node");
+      }
+    else
+      throw std::runtime_error("Invalid value for delta");
+  } catch (std::exception& exception) {
+    console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what());
+    status = false;
+  }
+  return status;
+}
+
+//! Assign absorbing constraints to nodes
+template <unsigned Tdim>
+bool mpm::Constraints<Tdim>::assign_nodal_absorbing_constraints(
+    const std::vector<std::tuple<mpm::Index, unsigned, double, double, double,
+                                 double, std::string>>& absorbing_constraints) {
+  bool status = true;
+  try {
+    for (const auto& absorbing_constraint : absorbing_constraints) {
+      // Node id
+      mpm::Index nid = std::get<0>(absorbing_constraint);
+      // Direction
+      unsigned dir = std::get<1>(absorbing_constraint);
+      // Delta
+      double delta = std::get<2>(absorbing_constraint);
+      // h_min
+      double h_min = std::get<3>(absorbing_constraint);
+      // a
+      double a = std::get<4>(absorbing_constraint);
+      // b
+      double b = std::get<5>(absorbing_constraint);
+      // Position
+      std::string position = std::get<6>(absorbing_constraint);
+      // delta check
+      if (delta >= h_min / (2 * a) and delta >= h_min / (2 * b)) {
+        if (position == "side" or position == "bottom" or
+            position == "corner") {
+          // Apply constraint
+          if (!mesh_->node(nid)->apply_absorbing_constraint(dir, delta, h_min,
+                                                            a, b, position))
+            throw std::runtime_error(
+                "Nodal absorbing constraints assignment failed");
+        } else
+          throw std::runtime_error("Invalid position");
+      } else
+        throw std::runtime_error("Invalid value for delta");
+    }
+  } catch (std::exception& exception) {
+    console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what());
+    status = false;
+  }
+  return status;
+}

include/mesh.h

@@ -24,6 +24,7 @@
 #include "json.hpp"
 using Json = nlohmann::json;
 
+#include "absorbing_constraint.h"
 #include "cell.h"
 #include "factory.h"
 #include "friction_constraint.h"

include/mesh.tcc

@@ -1991,6 +1991,10 @@ void mpm::Mesh<Tdim>::create_nodal_properties() {
                                        materials_.size());
     nodal_properties_->create_property("normal_unit_vectors", nrows,
                                        materials_.size());
+    nodal_properties_->create_property("wave_velocities", nrows,
+                                       materials_.size());
+    nodal_properties_->create_property("density", nodes_.size(),
+                                       materials_.size());
 
     // 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

include/node.h

@@ -213,6 +213,18 @@ class Node : public NodeBase<Tdim> {
   //! \param[in] dt Time-step
   void apply_friction_constraints(double dt) override;
 
+  //! Apply absorbing constraint
+  //! Directions can take values between 0 and Dim * Nphases
+  //! \param[in] dir Direction of absorbing constraint
+  //! \param[in] delta Virtual viscous layer thickness
+  //! \param[in] h_min Characteristic Length
+  //! \param[in] a Equation constant
+  //! \param[in] a Equation constant
+  //! \param[in] position Nodal position along boundary
+  bool apply_absorbing_constraint(unsigned dir, double delta, double h_min,
+                                  double a, double b,
+                                  std::string position) override;
+
   //! Assign rotation matrix
   //! \param[in] rotation_matrix Rotation matrix of the node
   void assign_rotation_matrix(

include/node.tcc

@@ -346,6 +346,80 @@ void mpm::Node<Tdim, Tdof, Tnphases>::apply_velocity_constraints() {
   }
 }
 
+// !Apply absorbing constraints
+template <unsigned Tdim, unsigned Tdof, unsigned Tnphases>
+bool mpm::Node<Tdim, Tdof, Tnphases>::apply_absorbing_constraint(
+    unsigned dir, double delta, double h_min, double a, double b,
+    std::string position) {
+  bool status = true;
+  try {
+
+    if (dir >= Tdim) {
+      throw std::runtime_error("Direction is out of bounds");
+      status = false;
+    }
+    if (material_ids_.size() != 0) {
+      // Get material id
+      auto mat_id = material_ids_.begin();
+
+      // Extract material properties and displacements
+      double pwave_v = this->property_handle_->property(
+          "wave_velocities", prop_id_, *mat_id, Tdim)(0);
+      double swave_v = this->property_handle_->property(
+          "wave_velocities", prop_id_, *mat_id, Tdim)(1);
+      double density =
+          this->property_handle_->property("density", prop_id_, *mat_id)(0);
+      Eigen::Matrix<double, Tdim, 1> material_displacement =
+          this->property_handle_->property("displacements", prop_id_, *mat_id,
+                                           Tdim);
+      material_displacement /= this->mass(*mat_id);
+
+      // Wave velocity Eigen Matrix
+      Eigen::Matrix<double, Tdim, 1> wave_velocity =
+          Eigen::MatrixXd::Constant(Tdim, 1, b * swave_v);
+      wave_velocity(dir, 0) = a * pwave_v;
+      // Spring constant Eigen matrix
+      double k_s = (density * pow(swave_v, 2)) / delta;
+      double k_p = (density * pow(pwave_v, 2)) / delta;
+      Eigen::Matrix<double, Tdim, 1> spring_constant =
+          Eigen::MatrixXd::Constant(Tdim, 1, k_s);
+      spring_constant(dir, 0) = k_p;
+
+      // Iterate through each phase
+      for (unsigned phase = 0; phase < Tnphases; ++phase) {
+        // Calculate Aborbing Traction
+        Eigen::Matrix<double, Tdim, 1> absorbing_traction_ =
+            this->velocity_.col(phase).cwiseProduct(wave_velocity) * density +
+            material_displacement.cwiseProduct(spring_constant);
+
+        // Update external force
+        if (Tdim == 2) {
+          Eigen::Matrix<double, Tdim, 1> absorbing_force_;
+          if (position == "corner") {
+            absorbing_force_ = 0.5 * h_min * absorbing_traction_;
+            this->update_external_force(true, phase, -absorbing_force_);
+          } else if (position == "bottom") {
+            double dir0 = 0.5 * h_min * absorbing_traction_(0);
+            double dir1 = 1.0 * h_min * absorbing_traction_(1);
+            absorbing_force_ << dir0, dir1;
+            this->update_external_force(true, phase, -absorbing_force_);
+          } else if (position == "side") {
+            double dir0 = 1.0 * h_min * absorbing_traction_(0);
+            double dir1 = 0.5 * h_min * absorbing_traction_(1);
+            absorbing_force_ << dir0, dir1;
+            this->update_external_force(true, phase, -absorbing_force_);
+          } else
+            throw std::runtime_error("Invalid absorbing boundary position");
+        }
+      }
+    }
+  } catch (std::exception& exception) {
+    console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what());
+    status = false;
+  }
+  return status;
+}
+
 //! Assign friction constraint
 //! Constrain directions can take values between 0 and Dim * Nphases
 template <unsigned Tdim, unsigned Tdof, unsigned Tnphases>
@@ -543,10 +617,15 @@ void mpm::Node<Tdim, Tdof, Tnphases>::update_property(
     bool update, const std::string& property,
     const Eigen::MatrixXd& property_value, unsigned mat_id,
     unsigned nprops) noexcept {
+
   // Update/assign property
   node_mutex_.lock();
-  property_handle_->update_property(property, prop_id_, mat_id, property_value,
-                                    nprops);
+  if (update)
+    property_handle_->update_property(property, prop_id_, mat_id,
+                                      property_value, nprops);
+  else
+    property_handle_->assign_property(property, prop_id_, mat_id,
+                                      property_value, nprops);
   node_mutex_.unlock();
 }
 

include/node_base.h

@@ -205,6 +205,18 @@ class NodeBase {
   //! \param[in] dt Time-step
   virtual void apply_friction_constraints(double dt) = 0;
 
+  //! Apply absorbing constraint
+  //! Directions can take values between 0 and Dim * Nphases
+  //! \param[in] dir Direction of absorbing constraint
+  //! \param[in] delta Virtual viscous layer thickness
+  //! \param[in] h_min Characteristic Length
+  //! \param[in] a Equation constant
+  //! \param[in] a Equation constant
+  //! \param[in] position Nodal position along boundary
+  virtual bool apply_absorbing_constraint(unsigned dir, double delta,
+                                          double h_min, double a, double b,
+                                          std::string position) = 0;
+
   //! Assign rotation matrix
   //! \param[in] rotation_matrix Rotation matrix of the node
   virtual void assign_rotation_matrix(

include/particles/particle.h

@@ -149,6 +149,9 @@ class Particle : public ParticleBase<Tdim> {
   //! Map multimaterial domain gradients to nodes
   void map_multimaterial_domain_gradients_to_nodes() noexcept override;
 
+  // ! Map linear elastic wave velocities to nodes
+  void map_wave_velocities_to_nodes() noexcept override;
+
   //! Assign nodal mass to particles
   //! \param[in] mass Mass from the particles in a cell
   //! \retval status Assignment status

include/particles/particle.tcc

@@ -582,6 +582,30 @@ void mpm::Particle<
   }
 }
 
+//! Map linear elastic wave velocities to nodes
+template <unsigned Tdim>
+void mpm::Particle<Tdim>::map_wave_velocities_to_nodes() noexcept {
+  // if (std::isnan(this->material())->template
+  // property<double>(std::string("pwave_velocity")) == false) {
+  Eigen::Matrix<double, 1, 1> density;
+  density(0) =
+      (this->material())->template property<double>(std::string("density"));
+  Eigen::Matrix<double, Tdim, 1> wave_velocities;
+  wave_velocities(0) =
+      (this->material())
+          ->template property<double>(std::string("pwave_velocity"));
+  wave_velocities(1) =
+      (this->material())
+          ->template property<double>(std::string("swave_velocity"));
+  for (unsigned i = 0; i < nodes_.size(); ++i) {
+    nodes_[i]->update_property(false, "wave_velocities", wave_velocities,
+                               this->material_id(), Tdim);
+    nodes_[i]->update_property(false, "density", density, this->material_id(),
+                               1);
+    //}
+  }
+}
+
 // Compute strain rate of the particle
 template <>
 inline Eigen::Matrix<double, 6, 1> mpm::Particle<1>::compute_strain_rate(

include/particles/particle_base.h

@@ -159,6 +159,9 @@ class ParticleBase {
   //! Map multimaterial domain gradients to nodes
   virtual void map_multimaterial_domain_gradients_to_nodes() noexcept = 0;
 
+  // ! Map linear elastic wave velocities to nodes
+  virtual void map_wave_velocities_to_nodes() noexcept = 0;
+
   //! Assign material
   virtual bool assign_material(const std::shared_ptr<Material<Tdim>>& material,
                                unsigned phase = mpm::ParticlePhase::Solid) = 0;