From 9c73b2599d97cc43005d1c9fe1943f4225de4394 Mon Sep 17 00:00:00 2001 From: Connor Geudeker Date: Tue, 23 Mar 2021 21:46:51 -0700 Subject: [PATCH 01/19] Kelvin-Voight Boundary Condition Implementation --- include/io/io_mesh.h | 6 +++ include/io/io_mesh_ascii.h | 6 +++ include/io/io_mesh_ascii.tcc | 51 ++++++++++++++++++++++ include/loads_bcs/constraints.h | 14 ++++++ include/loads_bcs/constraints.tcc | 62 ++++++++++++++++++++++++++ include/node.h | 17 ++++++++ include/node.tcc | 72 +++++++++++++++++++++++++++++++ include/node_base.h | 13 ++++++ include/solvers/mpm_base.h | 6 +++ include/solvers/mpm_base.tcc | 59 +++++++++++++++++++++++++ 10 files changed, 306 insertions(+) diff --git a/include/io/io_mesh.h b/include/io/io_mesh.h index 4150af5b0..9f036eb4b 100644 --- a/include/io/io_mesh.h +++ b/include/io/io_mesh.h @@ -98,6 +98,12 @@ class IOMesh { read_friction_constraints( const std::string& friction_constraints_file) = 0; + //! Read absorbing constraints file + //! \param[in] absorbing_constraints_files file name with tractions + virtual std::vector> + read_absorbing_constraints( + const std::string& absorbing_constraints_file) = 0; + //! Read forces file //! \param[in] forces_files file name with nodal concentrated force virtual std::vector> read_forces( diff --git a/include/io/io_mesh_ascii.h b/include/io/io_mesh_ascii.h index 433661546..0d468a7f0 100644 --- a/include/io/io_mesh_ascii.h +++ b/include/io/io_mesh_ascii.h @@ -86,6 +86,12 @@ class IOMeshAscii : public IOMesh { read_friction_constraints( const std::string& friction_constraints_file) override; + //! Read absorping constraints file + //! \param[in] absorbing_constraints_files file name with tractions + std::vector> + read_absorbing_constraints( + const std::string& absorbing_constraints_file) override; + //! Read traction file //! \param[in] forces_files file name with nodal concentrated force std::vector> read_forces( diff --git a/include/io/io_mesh_ascii.tcc b/include/io/io_mesh_ascii.tcc index 0f09fdc1a..c309d3ee0 100644 --- a/include/io/io_mesh_ascii.tcc +++ b/include/io/io_mesh_ascii.tcc @@ -486,6 +486,57 @@ std::vector> return constraints; } +//! Return absorbing constraints of particles +template +std::vector> + mpm::IOMeshAscii::read_absorbing_constraints( + const std::string& absorbing_constraints_file) { + + // Nodal absorbing constraints + std::vector> + constraints; + constraints.clear(); + + // input file stream + std::fstream file; + file.open(absorbing_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; + // Direction + unsigned dir; + // Delta + double delta; + // a + double a; + // b + double b; + // Read stream + istream >> id >> dir >> delta >> a >> b; + constraints.emplace_back(std::make_tuple(id, dir, delta, a, b)); + } + } + } + } + file.close(); + } catch (std::exception& exception) { + console_->error("Read absorbing constraints: {}", exception.what()); + file.close(); + } + return constraints; +} + //! Return particles force template std::vector> diff --git a/include/loads_bcs/constraints.h b/include/loads_bcs/constraints.h index 0efdef6f6..9c7604006 100644 --- a/include/loads_bcs/constraints.h +++ b/include/loads_bcs/constraints.h @@ -3,6 +3,7 @@ #include +#include "absorbing_constraint.h" #include "friction_constraint.h" #include "logger.h" #include "mesh.h" @@ -48,6 +49,19 @@ class Constraints { const std::vector>& friction_constraints); + //! Assign nodal absorbing constraints + //! \param[in] setid Node set id + //! \param[in] absorbing_constraints Constraint at node, dir, delta, a, b + bool assign_nodal_absorbing_constraint( + int nset_id, + const std::shared_ptr& aconstraints); + + //! Assign absorbing constraints to nodes + //! \param[in] absorbing_constraints Constraint at node, dir, delta, a, and b + bool assign_nodal_absorbing_constraints( + const std::vector>& absorbing_constraints); + private: //! Mesh object std::shared_ptr> mesh_; diff --git a/include/loads_bcs/constraints.tcc b/include/loads_bcs/constraints.tcc index 5813e609f..629a0e09a 100644 --- a/include/loads_bcs/constraints.tcc +++ b/include/loads_bcs/constraints.tcc @@ -105,3 +105,65 @@ bool mpm::Constraints::assign_nodal_friction_constraints( } return status; } + +//! Assign absorbing constraints to nodes +template +bool mpm::Constraints::assign_nodal_absorbing_constraint( + int nset_id, const std::shared_ptr& aconstraint) { + bool status = true; + try { + int set_id = aconstraint->setid(); + auto nset = mesh_->nodes(set_id); + if (nset.size() == 0) + throw std::runtime_error( + "Node set is empty for assignment of absorbing constraints"); + unsigned dir = aconstraint->dir(); + double delta = aconstraint->delta(); + double a = aconstraint->a(); + double b = aconstraint->b(); + double h_min = 1.0; // mesh_->cells()->mean_length_; + for (auto nitr = nset.cbegin(); nitr != nset.cend(); ++nitr) { + if (!(*nitr)->assign_absorbing_constraint(dir, delta, a, b, h_min)) + throw std::runtime_error( + "Failed to initialise absorbing constraint at node"); + } + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + status = false; + } + return status; +} + +//! Assign absorbing constraints to nodes +template +bool mpm::Constraints::assign_nodal_absorbing_constraints( + const std::vector>& + 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); + // a + double a = std::get<3>(absorbing_constraint); + // b + double b = std::get<4>(absorbing_constraint); + // Cell Height + double h_min = 1.0; // mesh_->cells()->mean_length_; + + // Apply constraint + if (!mesh_->node(nid)->assign_absorbing_constraint(dir, delta, a, b, + h_min)) + throw std::runtime_error( + "Nodal absorbing constraints assignment failed"); + } + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + status = false; + } + return status; +} \ No newline at end of file diff --git a/include/node.h b/include/node.h index 31a6b9a03..c56b1ab4b 100644 --- a/include/node.h +++ b/include/node.h @@ -213,6 +213,19 @@ class Node : public NodeBase { //! \param[in] dt Time-step void apply_friction_constraints(double dt) override; + //! Assign 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] a Equation constant + //! \param[in] a Equation constant + //! \param[in] h_min Characteristic Length + bool assign_absorbing_constraint(unsigned dir, double delta, double a, + double b, double h_min) override; + + //! Apply absorbing_constraint + void apply_absorbing_constraint() override; + //! Assign rotation matrix //! \param[in] rotation_matrix Rotation matrix of the node void assign_rotation_matrix( @@ -306,6 +319,10 @@ class Node : public NodeBase { //! A general velocity (non-Cartesian/inclined) constraint is specified at the //! node bool generic_boundary_constraints_{false}; + //! Absorbing constraints + std::tuple absorbing_constraint_; + //! Absorbing traction + Eigen::Matrix absorbing_traction_; //! Frictional constraints bool friction_{false}; std::tuple friction_constraint_; diff --git a/include/node.tcc b/include/node.tcc index 7d3a9cfcb..d9a55e238 100644 --- a/include/node.tcc +++ b/include/node.tcc @@ -345,6 +345,78 @@ void mpm::Node::apply_velocity_constraints() { } } } +//! Assign absorbing constraints +template +bool mpm::Node::assign_absorbing_constraint( + unsigned dir, double delta, double a, double b, double h_min) { + bool status = true; + try { + assert(dir <= Tdim); + if (delta >= h_min / (2 * a) and delta >= h_min / (2 * b)) { + this->absorbing_constraint_ = std::make_tuple( + static_cast(dir), static_cast(delta), + static_cast(a), static_cast(b)); + } else + throw std::runtime_error("Delta input is too small"); + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + status = false; + } + return status; +} + +// !Apply absorbing constraints +template +void mpm::Node::apply_absorbing_constraint() { + // Normal direction + const unsigned dir = std::get<0>(this->absorbing_constraint_); + + // Delta value + const double delta = std::get<1>(this->absorbing_constraint_); + + // a coefficient + const double a = std::get<2>(this->absorbing_constraint_); + + // b coefficient + const double b = std::get<3>(this->absorbing_constraint_); + + // Get material id + auto mat_id = material_ids_.begin(); + + // Extract material properties and displacements + double pwave_v = this->property_handle_->property("pwave_velocity", prop_id_, + *mat_id)(0, 0); + double swave_v = this->property_handle_->property("swave_velocity", prop_id_, + *mat_id)(0, 0); + double density = + this->property_handle_->property("density", prop_id_, *mat_id)(0, 0); + const Eigen::Matrix material_displacement = + this->property_handle_->property("displacements", prop_id_, *mat_id, + Tdim); + + // Wave velocity Eigen Matrix + Eigen::Matrix 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 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 absorbing_traction_ = + this->velocity_.col(phase).cwiseProduct(wave_velocity) * density + + material_displacement.cwiseProduct(spring_constant); + + // Update external force + this->update_external_force(true, phase, -absorbing_traction_); + } +} //! Assign friction constraint //! Constrain directions can take values between 0 and Dim * Nphases diff --git a/include/node_base.h b/include/node_base.h index 17eddfdca..b3ca70002 100644 --- a/include/node_base.h +++ b/include/node_base.h @@ -205,6 +205,19 @@ class NodeBase { //! \param[in] dt Time-step virtual void apply_friction_constraints(double dt) = 0; + //! Assign 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] a Equation constant + //! \param[in] a Equation constant + //! \param[in] h_min Characteristic Length + virtual bool assign_absorbing_constraint(unsigned dir, double delta, double a, + double b, double h_min) = 0; + + //! Apply absorbing_constraint + virtual void apply_absorbing_constraint() = 0; + //! Assign rotation matrix //! \param[in] rotation_matrix Rotation matrix of the node virtual void assign_rotation_matrix( diff --git a/include/solvers/mpm_base.h b/include/solvers/mpm_base.h index c84eaf1f8..1733427b0 100644 --- a/include/solvers/mpm_base.h +++ b/include/solvers/mpm_base.h @@ -127,6 +127,12 @@ class MPMBase : public MPM { void nodal_frictional_constraints( const Json& mesh_prop, const std::shared_ptr>& mesh_io); + //! Nodal absorbing constraints + //! \param[in] mesh_prop Mesh properties + //! \param[in] mesh_io Mesh IO handle + void nodal_absorbing_constraints( + const Json& mesh_prop, const std::shared_ptr>& mesh_io); + //! Cell entity sets //! \param[in] mesh_prop Mesh properties //! \param[in] check Check duplicates diff --git a/include/solvers/mpm_base.tcc b/include/solvers/mpm_base.tcc index f829b5f95..261e897d3 100644 --- a/include/solvers/mpm_base.tcc +++ b/include/solvers/mpm_base.tcc @@ -236,6 +236,9 @@ void mpm::MPMBase::initialise_mesh() { // Read and assign friction constraints this->nodal_frictional_constraints(mesh_props, mesh_io); + // Read and assign absorbing constraintes + this->nodal_absorbing_constraints(mesh_props, mesh_io); + // Initialise cell auto cells_begin = std::chrono::steady_clock::now(); // Shape function name @@ -937,6 +940,62 @@ void mpm::MPMBase::nodal_frictional_constraints( } } +// Nodal absorbing constraints +template +void mpm::MPMBase::nodal_absorbing_constraints( + const Json& mesh_props, const std::shared_ptr>& mesh_io) { + try { + // Read and assign absorbing constraints + if (mesh_props.find("boundary_conditions") != mesh_props.end() && + mesh_props["boundary_conditions"].find("absorbing_constraints") != + mesh_props["boundary_conditions"].end()) { + // Iterate over velocity constraints + for (const auto& constraints : + mesh_props["boundary_conditions"]["absorbing_constraints"]) { + // Absorbing constraints are specified in a file + if (constraints.find("file") != constraints.end()) { + std::string absorbing_constraints_file = + constraints.at("file").template get(); + bool absorbing_constraints = + constraints_->assign_nodal_absorbing_constraints( + mesh_io->read_absorbing_constraints( + io_->file_name(absorbing_constraints_file))); + if (!absorbing_constraints) + throw std::runtime_error( + "Absorbing constraints are not properly assigned"); + + } else { + + // Set id + int nset_id = constraints.at("nset_id").template get(); + // Direction + unsigned dir = constraints.at("dir").template get(); + // Delta + double delta = constraints.at("delta").template get(); + // a + double a = constraints.at("a").template get(); + // b + double b = constraints.at("b").template get(); + // Add absorbing constraint to mesh + auto absorbing_constraint = + std::make_shared(nset_id, dir, delta, a, + b); + bool absorbing_constraints = + constraints_->assign_nodal_absorbing_constraint( + nset_id, absorbing_constraint); + if (!absorbing_constraints) + throw std::runtime_error( + "Nodal absorbing constraint is not properly assigned"); + } + } + } else + throw std::runtime_error("Absorbing constraints JSON not found"); + + } catch (std::exception& exception) { + console_->warn("#{}: Absorbing conditions are undefined {} ", __LINE__, + exception.what()); + } +} //! Cell entity sets template void mpm::MPMBase::cell_entity_sets(const Json& mesh_props, From 55d7b19ee49f3df5ff1ec719ed6232030c4c41d7 Mon Sep 17 00:00:00 2001 From: Connor Geudeker Date: Wed, 24 Mar 2021 21:51:18 -0700 Subject: [PATCH 02/19] First Edits, and absorbing constraint file --- include/io/io_mesh.h | 6 -- include/io/io_mesh_ascii.h | 6 -- include/io/io_mesh_ascii.tcc | 51 ------------- include/loads_bcs/absorbing_constraint.h | 55 ++++++++++++++ include/loads_bcs/constraints.h | 10 ++- include/loads_bcs/constraints.tcc | 39 +++++----- include/node.h | 13 +--- include/node.tcc | 97 ++++++++++-------------- include/node_base.h | 9 +-- include/solvers/mpm_base.tcc | 17 +---- 10 files changed, 127 insertions(+), 176 deletions(-) create mode 100644 include/loads_bcs/absorbing_constraint.h diff --git a/include/io/io_mesh.h b/include/io/io_mesh.h index 9f036eb4b..4150af5b0 100644 --- a/include/io/io_mesh.h +++ b/include/io/io_mesh.h @@ -98,12 +98,6 @@ class IOMesh { read_friction_constraints( const std::string& friction_constraints_file) = 0; - //! Read absorbing constraints file - //! \param[in] absorbing_constraints_files file name with tractions - virtual std::vector> - read_absorbing_constraints( - const std::string& absorbing_constraints_file) = 0; - //! Read forces file //! \param[in] forces_files file name with nodal concentrated force virtual std::vector> read_forces( diff --git a/include/io/io_mesh_ascii.h b/include/io/io_mesh_ascii.h index 0d468a7f0..433661546 100644 --- a/include/io/io_mesh_ascii.h +++ b/include/io/io_mesh_ascii.h @@ -86,12 +86,6 @@ class IOMeshAscii : public IOMesh { read_friction_constraints( const std::string& friction_constraints_file) override; - //! Read absorping constraints file - //! \param[in] absorbing_constraints_files file name with tractions - std::vector> - read_absorbing_constraints( - const std::string& absorbing_constraints_file) override; - //! Read traction file //! \param[in] forces_files file name with nodal concentrated force std::vector> read_forces( diff --git a/include/io/io_mesh_ascii.tcc b/include/io/io_mesh_ascii.tcc index c309d3ee0..0f09fdc1a 100644 --- a/include/io/io_mesh_ascii.tcc +++ b/include/io/io_mesh_ascii.tcc @@ -486,57 +486,6 @@ std::vector> return constraints; } -//! Return absorbing constraints of particles -template -std::vector> - mpm::IOMeshAscii::read_absorbing_constraints( - const std::string& absorbing_constraints_file) { - - // Nodal absorbing constraints - std::vector> - constraints; - constraints.clear(); - - // input file stream - std::fstream file; - file.open(absorbing_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; - // Direction - unsigned dir; - // Delta - double delta; - // a - double a; - // b - double b; - // Read stream - istream >> id >> dir >> delta >> a >> b; - constraints.emplace_back(std::make_tuple(id, dir, delta, a, b)); - } - } - } - } - file.close(); - } catch (std::exception& exception) { - console_->error("Read absorbing constraints: {}", exception.what()); - file.close(); - } - return constraints; -} - //! Return particles force template std::vector> diff --git a/include/loads_bcs/absorbing_constraint.h b/include/loads_bcs/absorbing_constraint.h new file mode 100644 index 000000000..947502868 --- /dev/null +++ b/include/loads_bcs/absorbing_constraint.h @@ -0,0 +1,55 @@ +#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 + AbsorbingConstraint(int setid, unsigned dir, double delta, double h_min, + double a = 1, double b = 1) + : setid_{setid}, dir_{dir}, delta_{delta}, h_min_{h_min}, a_{a}, b_{b} {}; + + // 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_; } + + private: + // ID + int setid_; + // Direction + unsigned dir_; + // Delta + double delta_; + // Cell height + double h_min_; + // a + double a_; + // b + double b_; +}; +} // namespace mpm +#endif // MPM_ABSORBING_CONSTRAINT_H_ diff --git a/include/loads_bcs/constraints.h b/include/loads_bcs/constraints.h index 9c7604006..38ed700b9 100644 --- a/include/loads_bcs/constraints.h +++ b/include/loads_bcs/constraints.h @@ -51,15 +51,17 @@ class Constraints { //! Assign nodal absorbing constraints //! \param[in] setid Node set id - //! \param[in] absorbing_constraints Constraint at node, dir, delta, a, b + //! \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& aconstraints); + const std::shared_ptr& absorbing_constraints); //! Assign absorbing constraints to nodes - //! \param[in] absorbing_constraints Constraint at node, dir, delta, a, and b + //! \param[in] absorbing_constraints Constraint at node, dir, delta, h_min, a, + //! and b bool assign_nodal_absorbing_constraints( - const std::vector>& absorbing_constraints); private: diff --git a/include/loads_bcs/constraints.tcc b/include/loads_bcs/constraints.tcc index 629a0e09a..7574ff951 100644 --- a/include/loads_bcs/constraints.tcc +++ b/include/loads_bcs/constraints.tcc @@ -106,26 +106,27 @@ bool mpm::Constraints::assign_nodal_friction_constraints( return status; } -//! Assign absorbing constraints to nodes +//! Apply absorbing constraints to nodes template bool mpm::Constraints::assign_nodal_absorbing_constraint( - int nset_id, const std::shared_ptr& aconstraint) { + int nset_id, + const std::shared_ptr& absorbing_constraint) { bool status = true; try { - int set_id = aconstraint->setid(); + int set_id = absorbing_constraint->setid(); auto nset = mesh_->nodes(set_id); if (nset.size() == 0) throw std::runtime_error( - "Node set is empty for assignment of absorbing constraints"); - unsigned dir = aconstraint->dir(); - double delta = aconstraint->delta(); - double a = aconstraint->a(); - double b = aconstraint->b(); - double h_min = 1.0; // mesh_->cells()->mean_length_; + "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(); for (auto nitr = nset.cbegin(); nitr != nset.cend(); ++nitr) { - if (!(*nitr)->assign_absorbing_constraint(dir, delta, a, b, h_min)) + if (!(*nitr)->apply_absorbing_constraint(dir, delta, a, b, h_min)) throw std::runtime_error( - "Failed to initialise absorbing constraint at node"); + "Failed to apply absorbing constraint at node"); } } catch (std::exception& exception) { console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); @@ -137,8 +138,8 @@ bool mpm::Constraints::assign_nodal_absorbing_constraint( //! Assign absorbing constraints to nodes template bool mpm::Constraints::assign_nodal_absorbing_constraints( - const std::vector>& - absorbing_constraints) { + const std::vector>& absorbing_constraints) { bool status = true; try { for (const auto& absorbing_constraint : absorbing_constraints) { @@ -148,16 +149,16 @@ bool mpm::Constraints::assign_nodal_absorbing_constraints( 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<3>(absorbing_constraint); + double a = std::get<4>(absorbing_constraint); // b - double b = std::get<4>(absorbing_constraint); - // Cell Height - double h_min = 1.0; // mesh_->cells()->mean_length_; + double b = std::get<5>(absorbing_constraint); // Apply constraint - if (!mesh_->node(nid)->assign_absorbing_constraint(dir, delta, a, b, - h_min)) + if (!mesh_->node(nid)->apply_absorbing_constraint(dir, delta, a, b, + h_min)) throw std::runtime_error( "Nodal absorbing constraints assignment failed"); } diff --git a/include/node.h b/include/node.h index c56b1ab4b..3107c2ab6 100644 --- a/include/node.h +++ b/include/node.h @@ -213,18 +213,15 @@ class Node : public NodeBase { //! \param[in] dt Time-step void apply_friction_constraints(double dt) override; - //! Assign absorbing constraint + //! 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] a Equation constant //! \param[in] a Equation constant //! \param[in] h_min Characteristic Length - bool assign_absorbing_constraint(unsigned dir, double delta, double a, - double b, double h_min) override; - - //! Apply absorbing_constraint - void apply_absorbing_constraint() override; + bool apply_absorbing_constraint(unsigned dir, double delta, double a, + double b, double h_min) override; //! Assign rotation matrix //! \param[in] rotation_matrix Rotation matrix of the node @@ -319,10 +316,6 @@ class Node : public NodeBase { //! A general velocity (non-Cartesian/inclined) constraint is specified at the //! node bool generic_boundary_constraints_{false}; - //! Absorbing constraints - std::tuple absorbing_constraint_; - //! Absorbing traction - Eigen::Matrix absorbing_traction_; //! Frictional constraints bool friction_{false}; std::tuple friction_constraint_; diff --git a/include/node.tcc b/include/node.tcc index d9a55e238..28877e4a7 100644 --- a/include/node.tcc +++ b/include/node.tcc @@ -345,17 +345,51 @@ void mpm::Node::apply_velocity_constraints() { } } } -//! Assign absorbing constraints + +// !Apply absorbing constraints template -bool mpm::Node::assign_absorbing_constraint( +bool mpm::Node::apply_absorbing_constraint( unsigned dir, double delta, double a, double b, double h_min) { bool status = true; try { assert(dir <= Tdim); if (delta >= h_min / (2 * a) and delta >= h_min / (2 * b)) { - this->absorbing_constraint_ = std::make_tuple( - static_cast(dir), static_cast(delta), - static_cast(a), static_cast(b)); + // Get material id + auto mat_id = material_ids_.begin(); + + // Extract material properties and displacements + double pwave_v = this->property_handle_->property( + "pwave_velocity", prop_id_, *mat_id)(0, 0); + double swave_v = this->property_handle_->property( + "swave_velocity", prop_id_, *mat_id)(0, 0); + double density = + this->property_handle_->property("density", prop_id_, *mat_id)(0, 0); + const Eigen::Matrix material_displacement = + this->property_handle_->property("displacements", prop_id_, *mat_id, + Tdim); + + // Wave velocity Eigen Matrix + Eigen::Matrix 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 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 absorbing_traction_ = + this->velocity_.col(phase).cwiseProduct(wave_velocity) * density + + material_displacement.cwiseProduct(spring_constant); + + // Update external force + this->update_external_force(true, phase, -absorbing_traction_); + } } else throw std::runtime_error("Delta input is too small"); } catch (std::exception& exception) { @@ -365,59 +399,6 @@ bool mpm::Node::assign_absorbing_constraint( return status; } -// !Apply absorbing constraints -template -void mpm::Node::apply_absorbing_constraint() { - // Normal direction - const unsigned dir = std::get<0>(this->absorbing_constraint_); - - // Delta value - const double delta = std::get<1>(this->absorbing_constraint_); - - // a coefficient - const double a = std::get<2>(this->absorbing_constraint_); - - // b coefficient - const double b = std::get<3>(this->absorbing_constraint_); - - // Get material id - auto mat_id = material_ids_.begin(); - - // Extract material properties and displacements - double pwave_v = this->property_handle_->property("pwave_velocity", prop_id_, - *mat_id)(0, 0); - double swave_v = this->property_handle_->property("swave_velocity", prop_id_, - *mat_id)(0, 0); - double density = - this->property_handle_->property("density", prop_id_, *mat_id)(0, 0); - const Eigen::Matrix material_displacement = - this->property_handle_->property("displacements", prop_id_, *mat_id, - Tdim); - - // Wave velocity Eigen Matrix - Eigen::Matrix 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 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 absorbing_traction_ = - this->velocity_.col(phase).cwiseProduct(wave_velocity) * density + - material_displacement.cwiseProduct(spring_constant); - - // Update external force - this->update_external_force(true, phase, -absorbing_traction_); - } -} - //! Assign friction constraint //! Constrain directions can take values between 0 and Dim * Nphases template diff --git a/include/node_base.h b/include/node_base.h index b3ca70002..bedb39712 100644 --- a/include/node_base.h +++ b/include/node_base.h @@ -205,18 +205,15 @@ class NodeBase { //! \param[in] dt Time-step virtual void apply_friction_constraints(double dt) = 0; - //! Assign absorbing constraint + //! 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] a Equation constant //! \param[in] a Equation constant //! \param[in] h_min Characteristic Length - virtual bool assign_absorbing_constraint(unsigned dir, double delta, double a, - double b, double h_min) = 0; - - //! Apply absorbing_constraint - virtual void apply_absorbing_constraint() = 0; + virtual bool apply_absorbing_constraint(unsigned dir, double delta, double a, + double b, double h_min) = 0; //! Assign rotation matrix //! \param[in] rotation_matrix Rotation matrix of the node diff --git a/include/solvers/mpm_base.tcc b/include/solvers/mpm_base.tcc index 261e897d3..b52779a68 100644 --- a/include/solvers/mpm_base.tcc +++ b/include/solvers/mpm_base.tcc @@ -949,23 +949,9 @@ void mpm::MPMBase::nodal_absorbing_constraints( if (mesh_props.find("boundary_conditions") != mesh_props.end() && mesh_props["boundary_conditions"].find("absorbing_constraints") != mesh_props["boundary_conditions"].end()) { - // Iterate over velocity constraints + // Iterate over absorbing constraints for (const auto& constraints : mesh_props["boundary_conditions"]["absorbing_constraints"]) { - // Absorbing constraints are specified in a file - if (constraints.find("file") != constraints.end()) { - std::string absorbing_constraints_file = - constraints.at("file").template get(); - bool absorbing_constraints = - constraints_->assign_nodal_absorbing_constraints( - mesh_io->read_absorbing_constraints( - io_->file_name(absorbing_constraints_file))); - if (!absorbing_constraints) - throw std::runtime_error( - "Absorbing constraints are not properly assigned"); - - } else { - // Set id int nset_id = constraints.at("nset_id").template get(); // Direction @@ -986,7 +972,6 @@ void mpm::MPMBase::nodal_absorbing_constraints( if (!absorbing_constraints) throw std::runtime_error( "Nodal absorbing constraint is not properly assigned"); - } } } else throw std::runtime_error("Absorbing constraints JSON not found"); From e1d24cee064c32886774341fd3493f3940455cd9 Mon Sep 17 00:00:00 2001 From: Connor Geudeker Date: Tue, 18 May 2021 22:06:59 -0700 Subject: [PATCH 03/19] Test cases --- include/loads_bcs/constraints.tcc | 8 +-- include/mesh.tcc | 2 + include/node.tcc | 83 ++++++++++++++++--------------- include/particles/particle.h | 3 ++ include/particles/particle.tcc | 13 +++++ include/particles/particle_base.h | 3 ++ include/solvers/mpm_base.tcc | 4 +- tests/interface_test.cc | 22 ++++++++ tests/node_test.cc | 75 ++++++++++++++++++++++++++++ tests/particle_test.cc | 2 + 10 files changed, 170 insertions(+), 45 deletions(-) diff --git a/include/loads_bcs/constraints.tcc b/include/loads_bcs/constraints.tcc index 7574ff951..55e318090 100644 --- a/include/loads_bcs/constraints.tcc +++ b/include/loads_bcs/constraints.tcc @@ -123,8 +123,10 @@ bool mpm::Constraints::assign_nodal_absorbing_constraint( double h_min = absorbing_constraint->h_min(); double a = absorbing_constraint->a(); double b = absorbing_constraint->b(); + assert(delta <= h_min / (2 * a)); + assert(delta <= h_min / (2 * b)); for (auto nitr = nset.cbegin(); nitr != nset.cend(); ++nitr) { - if (!(*nitr)->apply_absorbing_constraint(dir, delta, a, b, h_min)) + if (!(*nitr)->apply_absorbing_constraint(dir, delta, h_min, a, b)) throw std::runtime_error( "Failed to apply absorbing constraint at node"); } @@ -157,8 +159,8 @@ bool mpm::Constraints::assign_nodal_absorbing_constraints( double b = std::get<5>(absorbing_constraint); // Apply constraint - if (!mesh_->node(nid)->apply_absorbing_constraint(dir, delta, a, b, - h_min)) + if (!mesh_->node(nid)->apply_absorbing_constraint(dir, delta, h_min, a, + b)) throw std::runtime_error( "Nodal absorbing constraints assignment failed"); } diff --git a/include/mesh.tcc b/include/mesh.tcc index 76e8c161e..6d23aa88c 100644 --- a/include/mesh.tcc +++ b/include/mesh.tcc @@ -1991,6 +1991,8 @@ void mpm::Mesh::create_nodal_properties() { materials_.size()); nodal_properties_->create_property("normal_unit_vectors", nrows, materials_.size()); + nodal_properties_->create_property("wave_velocities", nrows, + 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 diff --git a/include/node.tcc b/include/node.tcc index 28877e4a7..46bd895cc 100644 --- a/include/node.tcc +++ b/include/node.tcc @@ -349,49 +349,50 @@ void mpm::Node::apply_velocity_constraints() { // !Apply absorbing constraints template bool mpm::Node::apply_absorbing_constraint( - unsigned dir, double delta, double a, double b, double h_min) { + unsigned dir, double delta, double h_min, double a, double b) { bool status = true; try { - assert(dir <= Tdim); - if (delta >= h_min / (2 * a) and delta >= h_min / (2 * b)) { - // Get material id - auto mat_id = material_ids_.begin(); - - // Extract material properties and displacements - double pwave_v = this->property_handle_->property( - "pwave_velocity", prop_id_, *mat_id)(0, 0); - double swave_v = this->property_handle_->property( - "swave_velocity", prop_id_, *mat_id)(0, 0); - double density = - this->property_handle_->property("density", prop_id_, *mat_id)(0, 0); - const Eigen::Matrix material_displacement = - this->property_handle_->property("displacements", prop_id_, *mat_id, - Tdim); - - // Wave velocity Eigen Matrix - Eigen::Matrix 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 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 absorbing_traction_ = - this->velocity_.col(phase).cwiseProduct(wave_velocity) * density + - material_displacement.cwiseProduct(spring_constant); - - // Update external force - this->update_external_force(true, phase, -absorbing_traction_); - } - } else - throw std::runtime_error("Delta input is too small"); + + if (dir > Tdim) { + throw std::runtime_error("Direction is out of bounds"); + status = false; + } + // 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 material_displacement = + this->property_handle_->property("displacements", prop_id_, *mat_id, + Tdim); + + // Wave velocity Eigen Matrix + Eigen::Matrix 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 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 absorbing_traction_ = + this->velocity_.col(phase).cwiseProduct(wave_velocity) * density + + material_displacement.cwiseProduct(spring_constant); + + // Update external force + this->update_external_force(true, phase, -absorbing_traction_); + } } catch (std::exception& exception) { console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); status = false; diff --git a/include/particles/particle.h b/include/particles/particle.h index a362d7c76..f1dc40d5d 100644 --- a/include/particles/particle.h +++ b/include/particles/particle.h @@ -149,6 +149,9 @@ class Particle : public ParticleBase { //! 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 diff --git a/include/particles/particle.tcc b/include/particles/particle.tcc index b50a5eb56..45de298d7 100644 --- a/include/particles/particle.tcc +++ b/include/particles/particle.tcc @@ -582,6 +582,19 @@ void mpm::Particle< } } +//! Map linear elastic wave velocities to nodes +template +void mpm::Particle::map_wave_velocities_to_nodes() noexcept{ + // if (std::isnan(this->material())->template property(std::string("pwave_velocity")) == false) { + Eigen::Matrix wave_velocities; + wave_velocities(0) = (this->material())->template property(std::string("pwave_velocity")); + wave_velocities(1) = (this->material())->template property(std::string("swave_velocity")); + for (unsigned i = 0; i < nodes_.size(); ++i) { + nodes_[i]->update_property(true, "wave_velocities", wave_velocities, this->material_id(), Tdim); + //} +} +} + // Compute strain rate of the particle template <> inline Eigen::Matrix mpm::Particle<1>::compute_strain_rate( diff --git a/include/particles/particle_base.h b/include/particles/particle_base.h index dcb2bc419..260af1eac 100644 --- a/include/particles/particle_base.h +++ b/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, unsigned phase = mpm::ParticlePhase::Solid) = 0; diff --git a/include/solvers/mpm_base.tcc b/include/solvers/mpm_base.tcc index b52779a68..a10f37d2c 100644 --- a/include/solvers/mpm_base.tcc +++ b/include/solvers/mpm_base.tcc @@ -958,13 +958,15 @@ void mpm::MPMBase::nodal_absorbing_constraints( unsigned dir = constraints.at("dir").template get(); // Delta double delta = constraints.at("delta").template get(); + // h_min + double h_min = constraints.at("h_min").template get(); // a double a = constraints.at("a").template get(); // b double b = constraints.at("b").template get(); // Add absorbing constraint to mesh auto absorbing_constraint = - std::make_shared(nset_id, dir, delta, a, + std::make_shared(nset_id, dir, delta, h_min, a, b); bool absorbing_constraints = constraints_->assign_nodal_absorbing_constraint( diff --git a/tests/interface_test.cc b/tests/interface_test.cc index cb3b151d7..07848627c 100644 --- a/tests/interface_test.cc +++ b/tests/interface_test.cc @@ -44,6 +44,8 @@ TEST_CASE("Interface functions are checked", "[interface]") { Nmaterials); nodal_properties->create_property("normal_unit_vectors", Nnodes * Dim, Nmaterials); + nodal_properties->create_property("wave_velocities", Nnodes * Dim, + Nmaterials); // Element std::shared_ptr> element = @@ -204,6 +206,11 @@ TEST_CASE("Interface functions are checked", "[interface]") { node2->compute_multimaterial_normal_unit_vector(); node3->compute_multimaterial_normal_unit_vector(); + // Map wave velocities to nodes + particle1->map_wave_velocities_to_nodes(); + particle2->map_wave_velocities_to_nodes(); + particle3->map_wave_velocities_to_nodes(); + Eigen::Matrix masses; // clang-format off masses << 0.96, 1.46, @@ -283,6 +290,18 @@ TEST_CASE("Interface functions are checked", "[interface]") { 0.89442719099992, 0.88491822238198; // clang-format on + Eigen::Matrix wave_velocities; + // clang-format off + wave_velocities << 116.023870223, 2*116.023870223, + 62.0173672946, 2*62.0173672946, + 116.023870223, 2*116.023870223, + 62.0173672946, 2*62.0173672946, + 116.023870223, 2*116.023870223, + 62.0173672946, 2*62.0173672946, + 116.023870223, 2*116.023870223, + 62.0173672946, 2*62.0173672946; + // clang-format on + // Check if nodal properties were properly mapped and computed for (int i = 0; i < Nnodes; ++i) { for (int j = 0; j < Nmaterials; ++j) { @@ -304,6 +323,9 @@ TEST_CASE("Interface functions are checked", "[interface]") { Approx(gradients(i * Dim + k, j)).epsilon(tolerance)); REQUIRE(nodal_properties->property("normal_unit_vectors", i, j, Dim)( k, 0) == Approx(normal(i * Dim + k, j)).epsilon(tolerance)); + REQUIRE( + nodal_properties->property("wave_velocities", i, j, Dim)(k, 0) == + Approx(wave_velocities(i * Dim + k, j)).epsilon(tolerance)); } // Check if normal vector are also unit vectors REQUIRE( diff --git a/tests/node_test.cc b/tests/node_test.cc index 3cc92d2a5..9576a8642 100644 --- a/tests/node_test.cc +++ b/tests/node_test.cc @@ -726,6 +726,81 @@ TEST_CASE("Node is checked for 2D case", "[node][2D]") { REQUIRE(node->external_force(Nphase)(i) == Approx(0.).epsilon(Tolerance)); } + + SECTION("Check absorbing boundary") { + // Create a force vector + Eigen::Matrix force; + for (unsigned i = 0; i < force.size(); ++i) force(i) = 10. * (2 - i); + REQUIRE_NOTHROW(node->update_external_force(false, Nphase, force)); + + // Assign volume to 100 + double volume = 100.; + REQUIRE_NOTHROW(node->update_volume(false, Nphase, volume)); + REQUIRE(node->volume(Nphase) == Approx(100.0).epsilon(Tolerance)); + + // Assign mass to 100 + double mass = 100.; + REQUIRE_NOTHROW(node->update_mass(false, Nphase, mass)); + REQUIRE(node->mass(Nphase) == Approx(100.0).epsilon(Tolerance)); + + // Compute velocity + double dt = 0.1; + REQUIRE(node->compute_acceleration_velocity(Nphase, dt) == true); + + // Check velocity + Eigen::Matrix velocity = force / mass * dt; + for (unsigned i = 0; i < velocity.size(); ++i) + REQUIRE(node->velocity(Nphase)(i) == + Approx(velocity(i)).epsilon(Tolerance)); + + // Define nodal properties + Eigen::Matrix density; + density(0) = mass / volume; + Eigen::Matrix wave_v; + for (double i = 0; i < wave_v.size(); ++i) wave_v(i) = 100 * (2 - i); + Eigen::Matrix displacements; + for (double i = 0; i < displacements.size(); ++i) + displacements(i) = 0.0001 * (2 - i); + + // Define material id + unsigned mat_id = 0; + node->append_material_id(mat_id); + + // Assign nodal properties + auto nodal_properties = std::make_shared(); + nodal_properties->create_property("density", 1, 1); + nodal_properties->create_property("wave_velocities", Dim, 1); + nodal_properties->create_property("displacements", Dim, 1); + REQUIRE_NOTHROW( + node->initialise_property_handle(0, nodal_properties)); + REQUIRE_NOTHROW( + node->update_property(false, "density", density, mat_id, 1)); + REQUIRE_NOTHROW(node->update_property(false, "wave_velocities", + wave_v, mat_id, Dim)); + REQUIRE_NOTHROW(node->update_property(false, "displacements", + displacements, mat_id, Dim)); + + // Check invalid direction + unsigned bad_dir = 4; + REQUIRE(node->apply_absorbing_constraint(bad_dir, 0, 1, 1, 1) == + false); + + // Apply absorbing constraint + unsigned dir = 0; + REQUIRE(node->apply_absorbing_constraint(dir, 1, 1, 1, 1) == true); + + // Calculation of absorbing force + double pwave_traction = -12.; + double swave_traction = -2.; + for (unsigned i = 0; i < Dim; ++i) { + if (i == dir) + REQUIRE(node->external_force(Nphase)(i) == + Approx(pwave_traction + force(i)).epsilon(Tolerance)); + else + REQUIRE(node->external_force(Nphase)(i) == + Approx(swave_traction + force(i)).epsilon(Tolerance)); + } + } } } diff --git a/tests/particle_test.cc b/tests/particle_test.cc index 9059a8940..494928f72 100644 --- a/tests/particle_test.cc +++ b/tests/particle_test.cc @@ -824,6 +824,8 @@ TEST_CASE("Particle is checked for 2D case", "[particle][2D]") { jmaterial["density"] = 1000.; jmaterial["youngs_modulus"] = 1.0E+7; jmaterial["poisson_ratio"] = 0.3; + // jmaterial["pwave_velocity"] = 116.023870223; + // jmaterial["swave_velocity"] = 62.0173672946; auto material = Factory, unsigned, const Json&>::instance()->create( From e30d18652d4be542992165a8f066a2ecfc49f2c9 Mon Sep 17 00:00:00 2001 From: Connor Geudeker Date: Wed, 19 May 2021 01:09:47 -0700 Subject: [PATCH 04/19] Test cases, minor fix --- include/node.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/node.h b/include/node.h index 3107c2ab6..93857f7e5 100644 --- a/include/node.h +++ b/include/node.h @@ -220,8 +220,8 @@ class Node : public NodeBase { //! \param[in] a Equation constant //! \param[in] a Equation constant //! \param[in] h_min Characteristic Length - bool apply_absorbing_constraint(unsigned dir, double delta, double a, - double b, double h_min) override; + bool apply_absorbing_constraint(unsigned dir, double delta, double h_min, + double a, double b) override; //! Assign rotation matrix //! \param[in] rotation_matrix Rotation matrix of the node From 6a5328077aa38a20d506034a24612fd338b40f57 Mon Sep 17 00:00:00 2001 From: Connor Geudeker Date: Mon, 24 May 2021 20:59:14 -0700 Subject: [PATCH 05/19] Test Case for Absorbing Constraint File --- include/loads_bcs/constraints.tcc | 4 +- include/particles/particle.tcc | 22 ++-- tests/mesh_test_2d.cc | 167 ++++++++++++++++++++++++++++++ tests/particle_test.cc | 2 - 4 files changed, 183 insertions(+), 12 deletions(-) diff --git a/include/loads_bcs/constraints.tcc b/include/loads_bcs/constraints.tcc index 55e318090..429401567 100644 --- a/include/loads_bcs/constraints.tcc +++ b/include/loads_bcs/constraints.tcc @@ -123,8 +123,8 @@ bool mpm::Constraints::assign_nodal_absorbing_constraint( double h_min = absorbing_constraint->h_min(); double a = absorbing_constraint->a(); double b = absorbing_constraint->b(); - assert(delta <= h_min / (2 * a)); - assert(delta <= h_min / (2 * b)); + assert(delta > h_min / (2 * a)); + assert(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)) throw std::runtime_error( diff --git a/include/particles/particle.tcc b/include/particles/particle.tcc index 45de298d7..ad4511f17 100644 --- a/include/particles/particle.tcc +++ b/include/particles/particle.tcc @@ -584,15 +584,21 @@ void mpm::Particle< //! Map linear elastic wave velocities to nodes template -void mpm::Particle::map_wave_velocities_to_nodes() noexcept{ - // if (std::isnan(this->material())->template property(std::string("pwave_velocity")) == false) { - Eigen::Matrix wave_velocities; - wave_velocities(0) = (this->material())->template property(std::string("pwave_velocity")); - wave_velocities(1) = (this->material())->template property(std::string("swave_velocity")); +void mpm::Particle::map_wave_velocities_to_nodes() noexcept { + // if (std::isnan(this->material())->template + // property(std::string("pwave_velocity")) == false) { + Eigen::Matrix wave_velocities; + wave_velocities(0) = + (this->material()) + ->template property(std::string("pwave_velocity")); + wave_velocities(1) = + (this->material()) + ->template property(std::string("swave_velocity")); for (unsigned i = 0; i < nodes_.size(); ++i) { - nodes_[i]->update_property(true, "wave_velocities", wave_velocities, this->material_id(), Tdim); - //} -} + nodes_[i]->update_property(false, "wave_velocities", wave_velocities, + this->material_id(), Tdim); + //} + } } // Compute strain rate of the particle diff --git a/tests/mesh_test_2d.cc b/tests/mesh_test_2d.cc index eedc61044..8bc2211cf 100644 --- a/tests/mesh_test_2d.cc +++ b/tests/mesh_test_2d.cc @@ -1150,6 +1150,105 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") { set_id, friction_constraint) == false); } + SECTION("Check assign absorbing constraints") { + tsl::robin_map> node_sets; + node_sets[0] = std::vector{0, 2}; + node_sets[1] = std::vector{1, 3}; + + REQUIRE(mesh->create_node_sets(node_sets, true) == true); + + //! Constraints object + auto constraints = std::make_shared>(mesh); + + //! Required material information for nodes + // Initialise material models + mesh->initialise_material_models(materials); + + // Check nodal properties creation + REQUIRE_NOTHROW(mesh->create_nodal_properties()); + + // Check nodal properties initialisation + REQUIRE_NOTHROW(mesh->initialise_nodal_properties()); + + // Update Relevant nodal properties + Eigen::Matrix density; + density(0) = 1000; + Eigen::Matrix wave_v; + for (double i = 0; i < wave_v.size(); ++i) wave_v(i) = 100 * (2 - i); + Eigen::Matrix displacements; + for (double i = 0; i < displacements.size(); ++i) + displacements(i) = 0.0001 * (2 - i); + + // Define material id + unsigned mat_id = 0; + mesh->node(0)->append_material_id(mat_id); + + // Assign nodal properties + auto nodal_properties = std::make_shared(); + nodal_properties->create_property("density", 1, 1); + nodal_properties->create_property("wave_velocities", Dim, 1); + nodal_properties->create_property("displacements", Dim, 1); + for (double i = 0; i < 4; ++i) { + REQUIRE_NOTHROW(mesh->node(i)->append_material_id(mat_id)); + REQUIRE_NOTHROW( + mesh->node(i)->initialise_property_handle(0, nodal_properties)); + REQUIRE_NOTHROW(mesh->node(i)->update_property(false, "density", + density, mat_id, 1)); + REQUIRE_NOTHROW(mesh->node(i)->update_property( + false, "wave_velocities", wave_v, mat_id, Dim)); + REQUIRE_NOTHROW(mesh->node(i)->update_property( + false, "displacements", displacements, mat_id, Dim)); + } + + int set_id = 0; + unsigned dir = 0; + double delta = 1; + double h_min = 1; + double a = 1; + double b = 1; + + // Add absorbing constraint to mesh + auto absorbing_constraint = + std::make_shared(set_id, dir, delta, + h_min, a, b); + REQUIRE(constraints->assign_nodal_absorbing_constraint( + set_id, absorbing_constraint) == true); + + set_id = 1; + dir = 0; + delta = 3; + h_min = 12; + a = 2; + b = 2; + // Add absorbing constraint to mesh + absorbing_constraint = std::make_shared( + set_id, dir, delta, h_min, a, b); + REQUIRE(constraints->assign_nodal_absorbing_constraint( + set_id, absorbing_constraint) == true); + + // Add absorbing constraint to all nodes in mesh + absorbing_constraint = std::make_shared( + set_id, dir, delta, h_min, a, b); + REQUIRE(constraints->assign_nodal_absorbing_constraint( + set_id, absorbing_constraint) == true); + + // When constraints fail: invalid direction + dir = 3; + // Add absorbing constraint to mesh + absorbing_constraint = std::make_shared( + set_id, dir, delta, h_min, a, b); + REQUIRE(constraints->assign_nodal_absorbing_constraint( + set_id, absorbing_constraint) == false); + + // When constraints fail: invalid delta + delta = 0; + // Add absorbing constraint to mesh + absorbing_constraint = std::make_shared( + set_id, dir, delta, h_min, a, b); + REQUIRE(constraints->assign_nodal_absorbing_constraint( + set_id, absorbing_constraint) == false); + } + // Test assign velocity constraints to nodes SECTION("Check assign velocity constraints to nodes") { // Vector of particle coordinates @@ -1192,6 +1291,74 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") { friction_constraints) == false); } + // Test assign absorbing constraints to nodes + SECTION("Check assign absorbing constraints to nodes") { + // Vector of particle coordinates + std::vector< + std::tuple> + absorbing_constraints; + //! Constraints object + auto constraints = std::make_shared>(mesh); + + //! Required material information for nodes + // Initialise material models + mesh->initialise_material_models(materials); + + // Check nodal properties creation + REQUIRE_NOTHROW(mesh->create_nodal_properties()); + + // Check nodal properties initialisation + REQUIRE_NOTHROW(mesh->initialise_nodal_properties()); + + // Update Relevant nodal properties + Eigen::Matrix density; + density(0) = 1000; + Eigen::Matrix wave_v; + for (double i = 0; i < wave_v.size(); ++i) wave_v(i) = 100 * (2 - i); + Eigen::Matrix displacements; + for (double i = 0; i < displacements.size(); ++i) + displacements(i) = 0.0001 * (2 - i); + + // Define material id + unsigned mat_id = 0; + mesh->node(0)->append_material_id(mat_id); + + // Assign nodal properties + auto nodal_properties = std::make_shared(); + nodal_properties->create_property("density", 1, 1); + nodal_properties->create_property("wave_velocities", Dim, 1); + nodal_properties->create_property("displacements", Dim, 1); + for (double i = 0; i < 4; ++i) { + REQUIRE_NOTHROW(mesh->node(i)->append_material_id(mat_id)); + REQUIRE_NOTHROW( + mesh->node(i)->initialise_property_handle(0, nodal_properties)); + REQUIRE_NOTHROW(mesh->node(i)->update_property(false, "density", + density, mat_id, 1)); + REQUIRE_NOTHROW(mesh->node(i)->update_property( + false, "wave_velocities", wave_v, mat_id, Dim)); + REQUIRE_NOTHROW(mesh->node(i)->update_property( + false, "displacements", displacements, mat_id, Dim)); + } + // Constraint + absorbing_constraints.emplace_back(std::make_tuple(0, 0, 1, 3, 2, 2)); + absorbing_constraints.emplace_back(std::make_tuple(1, 1, 2, 4, 1, 1)); + absorbing_constraints.emplace_back(std::make_tuple(2, 0, 1, 1, 1, 1)); + absorbing_constraints.emplace_back(std::make_tuple(3, 1, 3, 2, 3, 3)); + + // REQUIRE(constraints->assign_nodal_absorbing_constraints(absorbing_constraints) + // == true); + + // When constraints fail: invalid direction + absorbing_constraints.emplace_back(std::make_tuple(3, 2, 3, 2, 3, 3)); + // REQUIRE(constraints->assign_nodal_absorbing_constraints(absorbing_constraints) + // == false); + + // When constraints fail: invalid delta + absorbing_constraints.emplace_back(std::make_tuple(3, 1, 1, 3, 1, 1)); + // REQUIRE(constraints->assign_nodal_absorbing_constraints(absorbing_constraints) + // == false); + } + // Test assign nodes concentrated_forces SECTION("Check assign nodes concentrated_forces") { // Vector of node coordinates diff --git a/tests/particle_test.cc b/tests/particle_test.cc index 494928f72..9059a8940 100644 --- a/tests/particle_test.cc +++ b/tests/particle_test.cc @@ -824,8 +824,6 @@ TEST_CASE("Particle is checked for 2D case", "[particle][2D]") { jmaterial["density"] = 1000.; jmaterial["youngs_modulus"] = 1.0E+7; jmaterial["poisson_ratio"] = 0.3; - // jmaterial["pwave_velocity"] = 116.023870223; - // jmaterial["swave_velocity"] = 62.0173672946; auto material = Factory, unsigned, const Json&>::instance()->create( From 04380c315b45307ce06e4defa66c6e75c7b43e5a Mon Sep 17 00:00:00 2001 From: Connor Geudeker Date: Tue, 25 May 2021 00:33:14 -0700 Subject: [PATCH 06/19] Minor Fixes --- include/loads_bcs/constraints.tcc | 3 ++- include/node.tcc | 2 +- tests/mesh_test_2d.cc | 12 ++++++------ 3 files changed, 9 insertions(+), 8 deletions(-) diff --git a/include/loads_bcs/constraints.tcc b/include/loads_bcs/constraints.tcc index 429401567..a83daad92 100644 --- a/include/loads_bcs/constraints.tcc +++ b/include/loads_bcs/constraints.tcc @@ -157,7 +157,8 @@ bool mpm::Constraints::assign_nodal_absorbing_constraints( double a = std::get<4>(absorbing_constraint); // b double b = std::get<5>(absorbing_constraint); - + assert(delta > h_min / (2 * a)); + assert(delta > h_min / (2 * b)); // Apply constraint if (!mesh_->node(nid)->apply_absorbing_constraint(dir, delta, h_min, a, b)) diff --git a/include/node.tcc b/include/node.tcc index 46bd895cc..65a34d94f 100644 --- a/include/node.tcc +++ b/include/node.tcc @@ -353,7 +353,7 @@ bool mpm::Node::apply_absorbing_constraint( bool status = true; try { - if (dir > Tdim) { + if (dir >= Tdim) { throw std::runtime_error("Direction is out of bounds"); status = false; } diff --git a/tests/mesh_test_2d.cc b/tests/mesh_test_2d.cc index 8bc2211cf..929a3a93f 100644 --- a/tests/mesh_test_2d.cc +++ b/tests/mesh_test_2d.cc @@ -1345,18 +1345,18 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") { absorbing_constraints.emplace_back(std::make_tuple(2, 0, 1, 1, 1, 1)); absorbing_constraints.emplace_back(std::make_tuple(3, 1, 3, 2, 3, 3)); - // REQUIRE(constraints->assign_nodal_absorbing_constraints(absorbing_constraints) - // == true); + REQUIRE(constraints->assign_nodal_absorbing_constraints( + absorbing_constraints) == true); // When constraints fail: invalid direction absorbing_constraints.emplace_back(std::make_tuple(3, 2, 3, 2, 3, 3)); - // REQUIRE(constraints->assign_nodal_absorbing_constraints(absorbing_constraints) - // == false); + REQUIRE(constraints->assign_nodal_absorbing_constraints( + absorbing_constraints) == false); // When constraints fail: invalid delta absorbing_constraints.emplace_back(std::make_tuple(3, 1, 1, 3, 1, 1)); - // REQUIRE(constraints->assign_nodal_absorbing_constraints(absorbing_constraints) - // == false); + REQUIRE(constraints->assign_nodal_absorbing_constraints( + absorbing_constraints) == false); } // Test assign nodes concentrated_forces From 4999748eb952c641347cdfbbd2ca7ea7d99e6d2b Mon Sep 17 00:00:00 2001 From: Connor Geudeker Date: Tue, 25 May 2021 16:58:38 -0700 Subject: [PATCH 07/19] Switched Assert for If Statement --- include/loads_bcs/constraints.tcc | 30 ++++++++------- tests/mesh_test_2d.cc | 63 +------------------------------ 2 files changed, 18 insertions(+), 75 deletions(-) diff --git a/include/loads_bcs/constraints.tcc b/include/loads_bcs/constraints.tcc index a83daad92..614332962 100644 --- a/include/loads_bcs/constraints.tcc +++ b/include/loads_bcs/constraints.tcc @@ -123,13 +123,14 @@ bool mpm::Constraints::assign_nodal_absorbing_constraint( double h_min = absorbing_constraint->h_min(); double a = absorbing_constraint->a(); double b = absorbing_constraint->b(); - assert(delta > h_min / (2 * a)); - assert(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)) - throw std::runtime_error( - "Failed to apply absorbing constraint at node"); - } + if (delta >= h_min / (2 * a) && 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)) + 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; @@ -157,13 +158,14 @@ bool mpm::Constraints::assign_nodal_absorbing_constraints( double a = std::get<4>(absorbing_constraint); // b double b = std::get<5>(absorbing_constraint); - assert(delta > h_min / (2 * a)); - assert(delta > h_min / (2 * b)); - // Apply constraint - if (!mesh_->node(nid)->apply_absorbing_constraint(dir, delta, h_min, a, - b)) - throw std::runtime_error( - "Nodal absorbing constraints assignment failed"); + if (delta >= h_min / (2 * a) && delta >= h_min / (2 * b)) { + // Apply constraint + if (!mesh_->node(nid)->apply_absorbing_constraint(dir, delta, h_min, a, + b)) + throw std::runtime_error( + "Nodal absorbing constraints assignment failed"); + } else + throw std::runtime_error("Invalid value for delta"); } } catch (std::exception& exception) { console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); diff --git a/tests/mesh_test_2d.cc b/tests/mesh_test_2d.cc index 929a3a93f..3c3bfdbd1 100644 --- a/tests/mesh_test_2d.cc +++ b/tests/mesh_test_2d.cc @@ -1160,44 +1160,14 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") { //! Constraints object auto constraints = std::make_shared>(mesh); - //! Required material information for nodes - // Initialise material models - mesh->initialise_material_models(materials); - - // Check nodal properties creation - REQUIRE_NOTHROW(mesh->create_nodal_properties()); - - // Check nodal properties initialisation - REQUIRE_NOTHROW(mesh->initialise_nodal_properties()); - - // Update Relevant nodal properties - Eigen::Matrix density; - density(0) = 1000; - Eigen::Matrix wave_v; - for (double i = 0; i < wave_v.size(); ++i) wave_v(i) = 100 * (2 - i); - Eigen::Matrix displacements; - for (double i = 0; i < displacements.size(); ++i) - displacements(i) = 0.0001 * (2 - i); - - // Define material id - unsigned mat_id = 0; - mesh->node(0)->append_material_id(mat_id); - // Assign nodal properties auto nodal_properties = std::make_shared(); nodal_properties->create_property("density", 1, 1); nodal_properties->create_property("wave_velocities", Dim, 1); nodal_properties->create_property("displacements", Dim, 1); for (double i = 0; i < 4; ++i) { - REQUIRE_NOTHROW(mesh->node(i)->append_material_id(mat_id)); REQUIRE_NOTHROW( mesh->node(i)->initialise_property_handle(0, nodal_properties)); - REQUIRE_NOTHROW(mesh->node(i)->update_property(false, "density", - density, mat_id, 1)); - REQUIRE_NOTHROW(mesh->node(i)->update_property( - false, "wave_velocities", wave_v, mat_id, Dim)); - REQUIRE_NOTHROW(mesh->node(i)->update_property( - false, "displacements", displacements, mat_id, Dim)); } int set_id = 0; @@ -1215,7 +1185,7 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") { set_id, absorbing_constraint) == true); set_id = 1; - dir = 0; + dir = 1; delta = 3; h_min = 12; a = 2; @@ -1300,45 +1270,16 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") { //! Constraints object auto constraints = std::make_shared>(mesh); - //! Required material information for nodes - // Initialise material models - mesh->initialise_material_models(materials); - - // Check nodal properties creation - REQUIRE_NOTHROW(mesh->create_nodal_properties()); - - // Check nodal properties initialisation - REQUIRE_NOTHROW(mesh->initialise_nodal_properties()); - - // Update Relevant nodal properties - Eigen::Matrix density; - density(0) = 1000; - Eigen::Matrix wave_v; - for (double i = 0; i < wave_v.size(); ++i) wave_v(i) = 100 * (2 - i); - Eigen::Matrix displacements; - for (double i = 0; i < displacements.size(); ++i) - displacements(i) = 0.0001 * (2 - i); - - // Define material id - unsigned mat_id = 0; - mesh->node(0)->append_material_id(mat_id); - // Assign nodal properties auto nodal_properties = std::make_shared(); nodal_properties->create_property("density", 1, 1); nodal_properties->create_property("wave_velocities", Dim, 1); nodal_properties->create_property("displacements", Dim, 1); for (double i = 0; i < 4; ++i) { - REQUIRE_NOTHROW(mesh->node(i)->append_material_id(mat_id)); REQUIRE_NOTHROW( mesh->node(i)->initialise_property_handle(0, nodal_properties)); - REQUIRE_NOTHROW(mesh->node(i)->update_property(false, "density", - density, mat_id, 1)); - REQUIRE_NOTHROW(mesh->node(i)->update_property( - false, "wave_velocities", wave_v, mat_id, Dim)); - REQUIRE_NOTHROW(mesh->node(i)->update_property( - false, "displacements", displacements, mat_id, Dim)); } + // Constraint absorbing_constraints.emplace_back(std::make_tuple(0, 0, 1, 3, 2, 2)); absorbing_constraints.emplace_back(std::make_tuple(1, 1, 2, 4, 1, 1)); From e1ccc5317e76c4f3dac467f9aeec68f171182ea1 Mon Sep 17 00:00:00 2001 From: Connor Geudeker Date: Sun, 1 Aug 2021 12:00:02 -0700 Subject: [PATCH 08/19] Absorbing Boundary in Scheme --- include/loads_bcs/constraints.tcc | 4 +- include/mesh.h | 1 + include/mesh.tcc | 2 + include/node.tcc | 83 +++++++----- include/node_base.h | 6 +- include/particles/particle.tcc | 5 + include/solvers/mpm_base.h | 12 ++ include/solvers/mpm_base.tcc | 64 +++++---- include/solvers/mpm_explicit.h | 11 ++ include/solvers/mpm_explicit.tcc | 7 + include/solvers/mpm_scheme/mpm_scheme.h | 3 + include/solvers/mpm_scheme/mpm_scheme.tcc | 22 +++ tests/interface_test.cc | 17 +-- tests/node_test.cc | 156 ++++++++++++---------- 14 files changed, 246 insertions(+), 147 deletions(-) diff --git a/include/loads_bcs/constraints.tcc b/include/loads_bcs/constraints.tcc index 614332962..4cf1c7b2e 100644 --- a/include/loads_bcs/constraints.tcc +++ b/include/loads_bcs/constraints.tcc @@ -123,7 +123,7 @@ bool mpm::Constraints::assign_nodal_absorbing_constraint( double h_min = absorbing_constraint->h_min(); double a = absorbing_constraint->a(); double b = absorbing_constraint->b(); - if (delta >= h_min / (2 * a) && delta >= h_min / (2 * b)) + 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)) throw std::runtime_error( @@ -158,7 +158,7 @@ bool mpm::Constraints::assign_nodal_absorbing_constraints( double a = std::get<4>(absorbing_constraint); // b double b = std::get<5>(absorbing_constraint); - if (delta >= h_min / (2 * a) && delta >= h_min / (2 * b)) { + if (delta >= h_min / (2 * a) and delta >= h_min / (2 * b)) { // Apply constraint if (!mesh_->node(nid)->apply_absorbing_constraint(dir, delta, h_min, a, b)) diff --git a/include/mesh.h b/include/mesh.h index bb51f161a..10193dbda 100644 --- a/include/mesh.h +++ b/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" diff --git a/include/mesh.tcc b/include/mesh.tcc index 6d23aa88c..62d084a60 100644 --- a/include/mesh.tcc +++ b/include/mesh.tcc @@ -1993,6 +1993,8 @@ void mpm::Mesh::create_nodal_properties() { 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 diff --git a/include/node.tcc b/include/node.tcc index 65a34d94f..c96150a46 100644 --- a/include/node.tcc +++ b/include/node.tcc @@ -1,3 +1,4 @@ +#include //! Constructor with id, coordinates and dof template mpm::Node::Node( @@ -357,41 +358,44 @@ bool mpm::Node::apply_absorbing_constraint( throw std::runtime_error("Direction is out of bounds"); status = false; } - // 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 material_displacement = - this->property_handle_->property("displacements", prop_id_, *mat_id, - Tdim); - - // Wave velocity Eigen Matrix - Eigen::Matrix 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 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 absorbing_traction_ = - this->velocity_.col(phase).cwiseProduct(wave_velocity) * density + - material_displacement.cwiseProduct(spring_constant); - - // Update external force - this->update_external_force(true, phase, -absorbing_traction_); + 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 material_displacement = + this->property_handle_->property("displacements", prop_id_, *mat_id, + Tdim); + + // Wave velocity Eigen Matrix + Eigen::Matrix 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 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 absorbing_traction_ = + this->velocity_.col(phase).cwiseProduct(wave_velocity) * density + + material_displacement.cwiseProduct(spring_constant); + + // Update external force + this->update_external_force(true, phase, -absorbing_traction_); + std::cout << "\n Absorbing Traction: \n" << std::ends; + std::cout << absorbing_traction_ << std::ends; + } } } catch (std::exception& exception) { console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); @@ -597,10 +601,15 @@ void mpm::Node::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(); } diff --git a/include/node_base.h b/include/node_base.h index bedb39712..0ec558d13 100644 --- a/include/node_base.h +++ b/include/node_base.h @@ -209,11 +209,11 @@ class NodeBase { //! 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] h_min Characteristic Length - virtual bool apply_absorbing_constraint(unsigned dir, double delta, double a, - double b, double h_min) = 0; + virtual bool apply_absorbing_constraint(unsigned dir, double delta, + double h_min, double a, double b) = 0; //! Assign rotation matrix //! \param[in] rotation_matrix Rotation matrix of the node diff --git a/include/particles/particle.tcc b/include/particles/particle.tcc index ad4511f17..0f3203e6d 100644 --- a/include/particles/particle.tcc +++ b/include/particles/particle.tcc @@ -587,6 +587,9 @@ template void mpm::Particle::map_wave_velocities_to_nodes() noexcept { // if (std::isnan(this->material())->template // property(std::string("pwave_velocity")) == false) { + Eigen::Matrix density; + density(0) = + (this->material())->template property(std::string("density")); Eigen::Matrix wave_velocities; wave_velocities(0) = (this->material()) @@ -597,6 +600,8 @@ void mpm::Particle::map_wave_velocities_to_nodes() noexcept { 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); //} } } diff --git a/include/solvers/mpm_base.h b/include/solvers/mpm_base.h index 1733427b0..a99c67b4b 100644 --- a/include/solvers/mpm_base.h +++ b/include/solvers/mpm_base.h @@ -99,6 +99,9 @@ class MPMBase : public MPM { //! Particle velocity constraints void particle_velocity_constraints(); + //! Apply Absorbing Constraints + void nodal_absorbing_constraints(); + private: //! Return if a mesh will be isoparametric or not //! \retval isoparametric Status of mesh type @@ -213,6 +216,15 @@ class MPMBase : public MPM { double damping_factor_{0.}; //! Locate particles bool locate_particles_{true}; + //! Absorbing Boundary Variables + bool absorbing_boundary_{false}; + std::shared_ptr absorbing_constraint_; + unsigned absorbing_dir_{0}; + int absorbing_nset_id_{0}; + double absorbing_delta_{0.}; + double absorbing_h_min_{0.}; + double absorbing_a_{0.}; + double absorbing_b_{0.}; #ifdef USE_GRAPH_PARTITIONING // graph pass the address of the container of cell diff --git a/include/solvers/mpm_base.tcc b/include/solvers/mpm_base.tcc index a10f37d2c..477fe2818 100644 --- a/include/solvers/mpm_base.tcc +++ b/include/solvers/mpm_base.tcc @@ -940,7 +940,7 @@ void mpm::MPMBase::nodal_frictional_constraints( } } -// Nodal absorbing constraints +// Assign nodal absorbing constraints template void mpm::MPMBase::nodal_absorbing_constraints( const Json& mesh_props, const std::shared_ptr>& mesh_io) { @@ -952,28 +952,35 @@ void mpm::MPMBase::nodal_absorbing_constraints( // Iterate over absorbing constraints for (const auto& constraints : mesh_props["boundary_conditions"]["absorbing_constraints"]) { - // Set id - int nset_id = constraints.at("nset_id").template get(); - // Direction - unsigned dir = constraints.at("dir").template get(); - // Delta - double delta = constraints.at("delta").template get(); - // h_min - double h_min = constraints.at("h_min").template get(); - // a - double a = constraints.at("a").template get(); - // b - double b = constraints.at("b").template get(); - // Add absorbing constraint to mesh - auto absorbing_constraint = - std::make_shared(nset_id, dir, delta, h_min, a, - b); - bool absorbing_constraints = - constraints_->assign_nodal_absorbing_constraint( - nset_id, absorbing_constraint); - if (!absorbing_constraints) - throw std::runtime_error( - "Nodal absorbing constraint is not properly assigned"); + // Set id + int nset_id = constraints.at("nset_id").template get(); + // Direction + unsigned dir = constraints.at("dir").template get(); + // Delta + double delta = constraints.at("delta").template get(); + // h_min + double h_min = constraints.at("h_min").template get(); + // a + double a = constraints.at("a").template get(); + // b + double b = constraints.at("b").template get(); + // Add absorbing constraint to mesh + auto absorbing_constraint = std::make_shared( + nset_id, dir, delta, h_min, a, b); + bool absorbing_constraints = + constraints_->assign_nodal_absorbing_constraint( + nset_id, absorbing_constraint); + if (!absorbing_constraints) + throw std::runtime_error( + "Nodal absorbing constraint is not properly assigned"); + absorbing_boundary_ = true; + absorbing_nset_id_ = nset_id; + absorbing_constraint_ = absorbing_constraint; + absorbing_dir_ = dir; + absorbing_delta_ = delta; + absorbing_h_min_ = h_min; + absorbing_a_ = a; + absorbing_b_ = b; } } else throw std::runtime_error("Absorbing constraints JSON not found"); @@ -983,6 +990,17 @@ void mpm::MPMBase::nodal_absorbing_constraints( exception.what()); } } + +// Apply nodal absorbing constraints +template +void mpm::MPMBase::nodal_absorbing_constraints() { + bool absorbing_constraints = constraints_->assign_nodal_absorbing_constraint( + absorbing_nset_id_, absorbing_constraint_); + if (!absorbing_constraints) + throw std::runtime_error( + "Nodal absorbing constraint is not properly assigned"); +} + //! Cell entity sets template void mpm::MPMBase::cell_entity_sets(const Json& mesh_props, diff --git a/include/solvers/mpm_explicit.h b/include/solvers/mpm_explicit.h index 2e463ac2d..c2a65b092 100644 --- a/include/solvers/mpm_explicit.h +++ b/include/solvers/mpm_explicit.h @@ -75,6 +75,17 @@ class MPMExplicit : public MPMBase { using mpm::MPMBase::damping_factor_; //! Locate particles using mpm::MPMBase::locate_particles_; + //! Constraints Pointer + using mpm::MPMBase::constraints_; + //! Absorbing Boundary + using mpm::MPMBase::absorbing_boundary_; + using mpm::MPMBase::absorbing_nset_id_; + using mpm::MPMBase::absorbing_constraint_; + using mpm::MPMBase::absorbing_dir_; + using mpm::MPMBase::absorbing_delta_; + using mpm::MPMBase::absorbing_h_min_; + using mpm::MPMBase::absorbing_a_; + using mpm::MPMBase::absorbing_b_; private: //! Pressure smoothing diff --git a/include/solvers/mpm_explicit.tcc b/include/solvers/mpm_explicit.tcc index 098a1488f..02005f2c9 100644 --- a/include/solvers/mpm_explicit.tcc +++ b/include/solvers/mpm_explicit.tcc @@ -110,6 +110,7 @@ bool mpm::MPMExplicit::solve() { this->initialise_loads(); auto solver_begin = std::chrono::steady_clock::now(); + // Main loop for (; step_ < nsteps_; ++step_) { @@ -145,6 +146,12 @@ bool mpm::MPMExplicit::solve() { mpm_scheme_->compute_forces(gravity_, phase, step_, set_node_concentrated_force_); + // Apply Absorbing Constraint + if (absorbing_boundary_) { + mpm_scheme_->absorbing_boundary_properties(); + this->nodal_absorbing_constraints(); + } + // Particle kinematics mpm_scheme_->compute_particle_kinematics(velocity_update_, phase, "Cundall", damping_factor_); diff --git a/include/solvers/mpm_scheme/mpm_scheme.h b/include/solvers/mpm_scheme/mpm_scheme.h index 683d149ad..ab08da28d 100644 --- a/include/solvers/mpm_scheme/mpm_scheme.h +++ b/include/solvers/mpm_scheme/mpm_scheme.h @@ -57,6 +57,9 @@ class MPMScheme { const Eigen::Matrix& gravity, unsigned phase, unsigned step, bool concentrated_nodal_forces); + //! Assign relevant properties for absorbing boundary + virtual inline void absorbing_boundary_properties(); + //! Compute acceleration velocity position //! \param[in] velocity_update Velocity or acceleration update flag //! \param[in] phase Phase of particle diff --git a/include/solvers/mpm_scheme/mpm_scheme.tcc b/include/solvers/mpm_scheme/mpm_scheme.tcc index e95c61795..d558eb4e0 100644 --- a/include/solvers/mpm_scheme/mpm_scheme.tcc +++ b/include/solvers/mpm_scheme/mpm_scheme.tcc @@ -168,6 +168,28 @@ inline void mpm::MPMScheme::compute_forces( #endif } +// Assign Absorbing Boundary Properties +template +inline void mpm::MPMScheme::absorbing_boundary_properties() { + mesh_->create_nodal_properties(); + + // Initialise nodal properties + mesh_->initialise_nodal_properties(); + + // Append material ids to nodes + mesh_->iterate_over_particles( + std::bind(&mpm::ParticleBase::append_material_id_to_nodes, + std::placeholders::_1)); + + mesh_->iterate_over_particles( + std::bind(&mpm::ParticleBase::map_wave_velocities_to_nodes, + std::placeholders::_1)); + // Map multimaterial displacements from particles to nodes + mesh_->iterate_over_particles(std::bind( + &mpm::ParticleBase::map_multimaterial_displacements_to_nodes, + std::placeholders::_1)); +} + // Compute particle kinematics template inline void mpm::MPMScheme::compute_particle_kinematics( diff --git a/tests/interface_test.cc b/tests/interface_test.cc index 07848627c..04ff8b024 100644 --- a/tests/interface_test.cc +++ b/tests/interface_test.cc @@ -46,6 +46,7 @@ TEST_CASE("Interface functions are checked", "[interface]") { Nmaterials); nodal_properties->create_property("wave_velocities", Nnodes * Dim, Nmaterials); + nodal_properties->create_property("density", Nnodes, Nmaterials); // Element std::shared_ptr> element = @@ -292,14 +293,14 @@ TEST_CASE("Interface functions are checked", "[interface]") { Eigen::Matrix wave_velocities; // clang-format off - wave_velocities << 116.023870223, 2*116.023870223, - 62.0173672946, 2*62.0173672946, - 116.023870223, 2*116.023870223, - 62.0173672946, 2*62.0173672946, - 116.023870223, 2*116.023870223, - 62.0173672946, 2*62.0173672946, - 116.023870223, 2*116.023870223, - 62.0173672946, 2*62.0173672946; + wave_velocities << 116.023870223, 116.023870223, + 62.0173672946, 62.0173672946, + 116.023870223, 116.023870223, + 62.0173672946, 62.0173672946, + 116.023870223, 116.023870223, + 62.0173672946, 62.0173672946, + 116.023870223, 116.023870223, + 62.0173672946, 62.0173672946; // clang-format on // Check if nodal properties were properly mapped and computed diff --git a/tests/node_test.cc b/tests/node_test.cc index 9576a8642..b06094e95 100644 --- a/tests/node_test.cc +++ b/tests/node_test.cc @@ -727,80 +727,88 @@ TEST_CASE("Node is checked for 2D case", "[node][2D]") { Approx(0.).epsilon(Tolerance)); } - SECTION("Check absorbing boundary") { - // Create a force vector - Eigen::Matrix force; - for (unsigned i = 0; i < force.size(); ++i) force(i) = 10. * (2 - i); - REQUIRE_NOTHROW(node->update_external_force(false, Nphase, force)); - - // Assign volume to 100 - double volume = 100.; - REQUIRE_NOTHROW(node->update_volume(false, Nphase, volume)); - REQUIRE(node->volume(Nphase) == Approx(100.0).epsilon(Tolerance)); - - // Assign mass to 100 - double mass = 100.; - REQUIRE_NOTHROW(node->update_mass(false, Nphase, mass)); - REQUIRE(node->mass(Nphase) == Approx(100.0).epsilon(Tolerance)); - - // Compute velocity - double dt = 0.1; - REQUIRE(node->compute_acceleration_velocity(Nphase, dt) == true); - - // Check velocity - Eigen::Matrix velocity = force / mass * dt; - for (unsigned i = 0; i < velocity.size(); ++i) - REQUIRE(node->velocity(Nphase)(i) == - Approx(velocity(i)).epsilon(Tolerance)); - - // Define nodal properties - Eigen::Matrix density; - density(0) = mass / volume; - Eigen::Matrix wave_v; - for (double i = 0; i < wave_v.size(); ++i) wave_v(i) = 100 * (2 - i); - Eigen::Matrix displacements; - for (double i = 0; i < displacements.size(); ++i) - displacements(i) = 0.0001 * (2 - i); - - // Define material id - unsigned mat_id = 0; - node->append_material_id(mat_id); - - // Assign nodal properties - auto nodal_properties = std::make_shared(); - nodal_properties->create_property("density", 1, 1); - nodal_properties->create_property("wave_velocities", Dim, 1); - nodal_properties->create_property("displacements", Dim, 1); - REQUIRE_NOTHROW( - node->initialise_property_handle(0, nodal_properties)); - REQUIRE_NOTHROW( - node->update_property(false, "density", density, mat_id, 1)); - REQUIRE_NOTHROW(node->update_property(false, "wave_velocities", - wave_v, mat_id, Dim)); - REQUIRE_NOTHROW(node->update_property(false, "displacements", - displacements, mat_id, Dim)); - - // Check invalid direction - unsigned bad_dir = 4; - REQUIRE(node->apply_absorbing_constraint(bad_dir, 0, 1, 1, 1) == - false); - - // Apply absorbing constraint - unsigned dir = 0; - REQUIRE(node->apply_absorbing_constraint(dir, 1, 1, 1, 1) == true); - - // Calculation of absorbing force - double pwave_traction = -12.; - double swave_traction = -2.; - for (unsigned i = 0; i < Dim; ++i) { - if (i == dir) - REQUIRE(node->external_force(Nphase)(i) == - Approx(pwave_traction + force(i)).epsilon(Tolerance)); - else - REQUIRE(node->external_force(Nphase)(i) == - Approx(swave_traction + force(i)).epsilon(Tolerance)); - } - } + /* SECTION("Check absorbing boundary") { + // Create a force vector + Eigen::Matrix force; + for (unsigned i = 0; i < force.size(); ++i) force(i) = 10. * + (2 - i); REQUIRE_NOTHROW(node->update_external_force(false, Nphase, + force)); + + // Assign volume to 100 + double volume = 100.; + REQUIRE_NOTHROW(node->update_volume(false, Nphase, volume)); + REQUIRE(node->volume(Nphase) == + Approx(100.0).epsilon(Tolerance)); + + // Assign mass to 100 + double mass = 100.; + REQUIRE_NOTHROW(node->update_mass(false, Nphase, mass)); + REQUIRE(node->mass(Nphase) == + Approx(100.0).epsilon(Tolerance)); + + // Compute velocity + double dt = 0.1; + REQUIRE(node->compute_acceleration_velocity(Nphase, dt) == + true); + + // Check velocity + Eigen::Matrix velocity = force / mass * dt; + for (unsigned i = 0; i < velocity.size(); ++i) + REQUIRE(node->velocity(Nphase)(i) == + Approx(velocity(i)).epsilon(Tolerance)); + + // Define nodal properties + Eigen::Matrix density; + density(0) = mass / volume; + Eigen::Matrix wave_v; + for (double i = 0; i < wave_v.size(); ++i) wave_v(i) = 100 * + (2 - i); Eigen::Matrix displacements; for (double i = + 0; i < displacements.size(); ++i) displacements(i) = 0.0001 * (2 - + i); + + // Define material id + unsigned mat_id = 0; + node->append_material_id(mat_id); + + // Assign nodal properties + auto nodal_properties = + std::make_shared(); + nodal_properties->create_property("density", 1, 1); + nodal_properties->create_property("wave_velocities", Dim, 1); + nodal_properties->create_property("displacements", Dim, 1); + REQUIRE_NOTHROW( + node->initialise_property_handle(0, nodal_properties)); + REQUIRE_NOTHROW( + node->update_property(false, "density", density, mat_id, + 1)); REQUIRE_NOTHROW(node->update_property(false, "wave_velocities", + wave_v, mat_id, Dim)); + REQUIRE_NOTHROW(node->update_property(false, "displacements", + displacements, mat_id, + Dim)); + + // Check invalid direction + unsigned bad_dir = 4; + REQUIRE(node->apply_absorbing_constraint(bad_dir, 0, 1, 1, 1) + == false); + + // Apply absorbing constraint + unsigned dir = 0; + REQUIRE(node->apply_absorbing_constraint(dir, 1, 1, 1, 1) == + true); + + // Calculation of absorbing force + double pwave_traction = -12.; + double swave_traction = -2.; + for (unsigned i = 0; i < Dim; ++i) { + if (i == dir) + REQUIRE(node->external_force(Nphase)(i) == + Approx(pwave_traction + + force(i)).epsilon(Tolerance)); else + REQUIRE(node->external_force(Nphase)(i) == + Approx(swave_traction + + force(i)).epsilon(Tolerance)); + } + }*/ } } From fe5f74e45484188170b452f23a35eb76f4ade007 Mon Sep 17 00:00:00 2001 From: jgiven100 Date: Thu, 19 Aug 2021 11:09:39 -0700 Subject: [PATCH 09/19] :hammer: add position input to absorb boundary --- include/loads_bcs/absorbing_constraint.h | 17 ++++++++++++++-- include/loads_bcs/constraints.h | 5 +++-- include/loads_bcs/constraints.tcc | 26 ++++++++++++++++-------- include/node.h | 6 ++++-- include/node.tcc | 26 ++++++++++++++++++++---- include/node_base.h | 4 +++- include/solvers/mpm_base.h | 5 +++-- include/solvers/mpm_base.tcc | 25 +++++++++++++++-------- include/solvers/mpm_explicit.h | 1 + 9 files changed, 86 insertions(+), 29 deletions(-) diff --git a/include/loads_bcs/absorbing_constraint.h b/include/loads_bcs/absorbing_constraint.h index 947502868..97d62a13b 100644 --- a/include/loads_bcs/absorbing_constraint.h +++ b/include/loads_bcs/absorbing_constraint.h @@ -15,9 +15,17 @@ class AbsorbingConstraint { //! \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) - : setid_{setid}, dir_{dir}, delta_{delta}, h_min_{h_min}, a_{a}, b_{b} {}; + 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_; } @@ -37,6 +45,9 @@ class AbsorbingConstraint { // Return b double b() const { return b_; } + // Return position + std::string position() const { return position_; } + private: // ID int setid_; @@ -50,6 +61,8 @@ class AbsorbingConstraint { double a_; // b double b_; + // Node position + std::string position_; }; } // namespace mpm #endif // MPM_ABSORBING_CONSTRAINT_H_ diff --git a/include/loads_bcs/constraints.h b/include/loads_bcs/constraints.h index 38ed700b9..6f2da49a4 100644 --- a/include/loads_bcs/constraints.h +++ b/include/loads_bcs/constraints.h @@ -59,10 +59,11 @@ class Constraints { //! Assign absorbing constraints to nodes //! \param[in] absorbing_constraints Constraint at node, dir, delta, h_min, a, - //! and b + //! b, and position bool assign_nodal_absorbing_constraints( const std::vector>& absorbing_constraints); + double, std::string>>& + absorbing_constraints); private: //! Mesh object diff --git a/include/loads_bcs/constraints.tcc b/include/loads_bcs/constraints.tcc index 4cf1c7b2e..3c1e0b055 100644 --- a/include/loads_bcs/constraints.tcc +++ b/include/loads_bcs/constraints.tcc @@ -113,7 +113,8 @@ bool mpm::Constraints::assign_nodal_absorbing_constraint( const std::shared_ptr& absorbing_constraint) { bool status = true; try { - int set_id = absorbing_constraint->setid(); + // 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( @@ -123,9 +124,11 @@ bool mpm::Constraints::assign_nodal_absorbing_constraint( 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)) + if (!(*nitr)->apply_absorbing_constraint(dir, delta, h_min, a, b, + position)) throw std::runtime_error( "Failed to apply absorbing constraint at node"); } @@ -142,7 +145,7 @@ bool mpm::Constraints::assign_nodal_absorbing_constraint( template bool mpm::Constraints::assign_nodal_absorbing_constraints( const std::vector>& absorbing_constraints) { + double, std::string>>& absorbing_constraints) { bool status = true; try { for (const auto& absorbing_constraint : absorbing_constraints) { @@ -158,12 +161,19 @@ bool mpm::Constraints::assign_nodal_absorbing_constraints( 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)) { - // Apply constraint - if (!mesh_->node(nid)->apply_absorbing_constraint(dir, delta, h_min, a, - b)) - throw std::runtime_error( - "Nodal absorbing constraints assignment failed"); + 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"); } diff --git a/include/node.h b/include/node.h index 93857f7e5..ae80e19f7 100644 --- a/include/node.h +++ b/include/node.h @@ -217,11 +217,13 @@ class Node : public NodeBase { //! 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] h_min Characteristic Length + //! \param[in] position Nodal position along boundary bool apply_absorbing_constraint(unsigned dir, double delta, double h_min, - double a, double b) override; + double a, double b, + std::string position) override; //! Assign rotation matrix //! \param[in] rotation_matrix Rotation matrix of the node diff --git a/include/node.tcc b/include/node.tcc index c96150a46..3ad050222 100644 --- a/include/node.tcc +++ b/include/node.tcc @@ -1,4 +1,5 @@ #include + //! Constructor with id, coordinates and dof template mpm::Node::Node( @@ -350,7 +351,8 @@ void mpm::Node::apply_velocity_constraints() { // !Apply absorbing constraints template bool mpm::Node::apply_absorbing_constraint( - unsigned dir, double delta, double h_min, double a, double b) { + unsigned dir, double delta, double h_min, double a, double b, + std::string position) { bool status = true; try { @@ -372,6 +374,7 @@ bool mpm::Node::apply_absorbing_constraint( Eigen::Matrix material_displacement = this->property_handle_->property("displacements", prop_id_, *mat_id, Tdim); + material_displacement /= this->mass(*mat_id); // Wave velocity Eigen Matrix Eigen::Matrix wave_velocity = @@ -392,9 +395,24 @@ bool mpm::Node::apply_absorbing_constraint( material_displacement.cwiseProduct(spring_constant); // Update external force - this->update_external_force(true, phase, -absorbing_traction_); - std::cout << "\n Absorbing Traction: \n" << std::ends; - std::cout << absorbing_traction_ << std::ends; + if (Tdim == 2) { + Eigen::Matrix 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) { diff --git a/include/node_base.h b/include/node_base.h index 0ec558d13..fc0defbdf 100644 --- a/include/node_base.h +++ b/include/node_base.h @@ -212,8 +212,10 @@ class NodeBase { //! \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) = 0; + double h_min, double a, double b, + std::string position) = 0; //! Assign rotation matrix //! \param[in] rotation_matrix Rotation matrix of the node diff --git a/include/solvers/mpm_base.h b/include/solvers/mpm_base.h index a99c67b4b..81e3ab1c5 100644 --- a/include/solvers/mpm_base.h +++ b/include/solvers/mpm_base.h @@ -218,13 +218,14 @@ class MPMBase : public MPM { bool locate_particles_{true}; //! Absorbing Boundary Variables bool absorbing_boundary_{false}; - std::shared_ptr absorbing_constraint_; + std::vector> absorbing_constraint_; unsigned absorbing_dir_{0}; - int absorbing_nset_id_{0}; + std::vector absorbing_nset_id_; double absorbing_delta_{0.}; double absorbing_h_min_{0.}; double absorbing_a_{0.}; double absorbing_b_{0.}; + std::string absorbing_position_{"null"}; #ifdef USE_GRAPH_PARTITIONING // graph pass the address of the container of cell diff --git a/include/solvers/mpm_base.tcc b/include/solvers/mpm_base.tcc index 477fe2818..2867af7c3 100644 --- a/include/solvers/mpm_base.tcc +++ b/include/solvers/mpm_base.tcc @@ -964,9 +964,12 @@ void mpm::MPMBase::nodal_absorbing_constraints( double a = constraints.at("a").template get(); // b double b = constraints.at("b").template get(); + // position + std::string position = + constraints.at("position").template get(); // Add absorbing constraint to mesh auto absorbing_constraint = std::make_shared( - nset_id, dir, delta, h_min, a, b); + nset_id, dir, delta, h_min, a, b, position); bool absorbing_constraints = constraints_->assign_nodal_absorbing_constraint( nset_id, absorbing_constraint); @@ -974,13 +977,14 @@ void mpm::MPMBase::nodal_absorbing_constraints( throw std::runtime_error( "Nodal absorbing constraint is not properly assigned"); absorbing_boundary_ = true; - absorbing_nset_id_ = nset_id; - absorbing_constraint_ = absorbing_constraint; + absorbing_nset_id_.emplace_back(nset_id); + absorbing_constraint_.emplace_back(absorbing_constraint); absorbing_dir_ = dir; absorbing_delta_ = delta; absorbing_h_min_ = h_min; absorbing_a_ = a; absorbing_b_ = b; + absorbing_position_ = position; } } else throw std::runtime_error("Absorbing constraints JSON not found"); @@ -994,11 +998,16 @@ void mpm::MPMBase::nodal_absorbing_constraints( // Apply nodal absorbing constraints template void mpm::MPMBase::nodal_absorbing_constraints() { - bool absorbing_constraints = constraints_->assign_nodal_absorbing_constraint( - absorbing_nset_id_, absorbing_constraint_); - if (!absorbing_constraints) - throw std::runtime_error( - "Nodal absorbing constraint is not properly assigned"); + for (int i = 0; i < absorbing_nset_id_.size(); i++) { + auto nset_id_ = absorbing_nset_id_.at(i); + auto a_constraint_ = absorbing_constraint_.at(i); + bool absorbing_constraints = + constraints_->assign_nodal_absorbing_constraint(nset_id_, + a_constraint_); + if (!absorbing_constraints) + throw std::runtime_error( + "Nodal absorbing constraint is not properly assigned"); + } } //! Cell entity sets diff --git a/include/solvers/mpm_explicit.h b/include/solvers/mpm_explicit.h index c2a65b092..3d7b2c035 100644 --- a/include/solvers/mpm_explicit.h +++ b/include/solvers/mpm_explicit.h @@ -86,6 +86,7 @@ class MPMExplicit : public MPMBase { using mpm::MPMBase::absorbing_h_min_; using mpm::MPMBase::absorbing_a_; using mpm::MPMBase::absorbing_b_; + using mpm::MPMBase::absorbing_position_; private: //! Pressure smoothing From 472fce524c7bcc3bc2fffae654ff33030365a0f4 Mon Sep 17 00:00:00 2001 From: jgiven100 Date: Thu, 19 Aug 2021 11:21:12 -0700 Subject: [PATCH 10/19] :hammer::fire: remove unused variables --- include/node.tcc | 2 -- include/solvers/mpm_base.h | 6 ------ include/solvers/mpm_base.tcc | 6 ------ include/solvers/mpm_explicit.h | 8 -------- 4 files changed, 22 deletions(-) diff --git a/include/node.tcc b/include/node.tcc index 3ad050222..be2a106b5 100644 --- a/include/node.tcc +++ b/include/node.tcc @@ -1,5 +1,3 @@ -#include - //! Constructor with id, coordinates and dof template mpm::Node::Node( diff --git a/include/solvers/mpm_base.h b/include/solvers/mpm_base.h index 81e3ab1c5..e57b60f0a 100644 --- a/include/solvers/mpm_base.h +++ b/include/solvers/mpm_base.h @@ -219,13 +219,7 @@ class MPMBase : public MPM { //! Absorbing Boundary Variables bool absorbing_boundary_{false}; std::vector> absorbing_constraint_; - unsigned absorbing_dir_{0}; std::vector absorbing_nset_id_; - double absorbing_delta_{0.}; - double absorbing_h_min_{0.}; - double absorbing_a_{0.}; - double absorbing_b_{0.}; - std::string absorbing_position_{"null"}; #ifdef USE_GRAPH_PARTITIONING // graph pass the address of the container of cell diff --git a/include/solvers/mpm_base.tcc b/include/solvers/mpm_base.tcc index 2867af7c3..03f985be5 100644 --- a/include/solvers/mpm_base.tcc +++ b/include/solvers/mpm_base.tcc @@ -979,12 +979,6 @@ void mpm::MPMBase::nodal_absorbing_constraints( absorbing_boundary_ = true; absorbing_nset_id_.emplace_back(nset_id); absorbing_constraint_.emplace_back(absorbing_constraint); - absorbing_dir_ = dir; - absorbing_delta_ = delta; - absorbing_h_min_ = h_min; - absorbing_a_ = a; - absorbing_b_ = b; - absorbing_position_ = position; } } else throw std::runtime_error("Absorbing constraints JSON not found"); diff --git a/include/solvers/mpm_explicit.h b/include/solvers/mpm_explicit.h index 3d7b2c035..1d8abda59 100644 --- a/include/solvers/mpm_explicit.h +++ b/include/solvers/mpm_explicit.h @@ -79,14 +79,6 @@ class MPMExplicit : public MPMBase { using mpm::MPMBase::constraints_; //! Absorbing Boundary using mpm::MPMBase::absorbing_boundary_; - using mpm::MPMBase::absorbing_nset_id_; - using mpm::MPMBase::absorbing_constraint_; - using mpm::MPMBase::absorbing_dir_; - using mpm::MPMBase::absorbing_delta_; - using mpm::MPMBase::absorbing_h_min_; - using mpm::MPMBase::absorbing_a_; - using mpm::MPMBase::absorbing_b_; - using mpm::MPMBase::absorbing_position_; private: //! Pressure smoothing From f3b70e6dd4f2441faadf47dc27201e03944d52e7 Mon Sep 17 00:00:00 2001 From: jgiven100 Date: Thu, 19 Aug 2021 14:39:30 -0700 Subject: [PATCH 11/19] :hammer: change absorb bound position from string to enum --- include/data_types.h | 7 +++++ include/loads_bcs/absorbing_constraint.h | 8 +++--- include/loads_bcs/constraints.h | 2 +- include/loads_bcs/constraints.tcc | 12 +++++---- include/node.h | 2 +- include/node.tcc | 33 +++++++++++++----------- include/node_base.h | 2 +- include/solvers/mpm_base.h | 1 + include/solvers/mpm_base.tcc | 10 ++++++- 9 files changed, 50 insertions(+), 27 deletions(-) diff --git a/include/data_types.h b/include/data_types.h index 061d27e81..21eca99c7 100644 --- a/include/data_types.h +++ b/include/data_types.h @@ -34,6 +34,13 @@ inline double zero() { return 0.; } +//! Position type +//! None: No position is specified +//! Side: Nodes along model vertical sides (left or right) +//! Corner: Nodes at model corner +//! Bottom: Nodes along model bottom +enum class Position { None, Side, Corner, Bottom }; + } // namespace mpm #endif // MPM_DATA_TYPES_H_ diff --git a/include/loads_bcs/absorbing_constraint.h b/include/loads_bcs/absorbing_constraint.h index 97d62a13b..84c7dd2c6 100644 --- a/include/loads_bcs/absorbing_constraint.h +++ b/include/loads_bcs/absorbing_constraint.h @@ -1,6 +1,8 @@ #ifndef MPM_ABSORBING_CONSTRAINT_H_ #define MPM_ABSORBING_CONSTRAINT_H_ +#include "data_types.h" + namespace mpm { //! AbsorbingConstraint class to store friction constraint on a set @@ -18,7 +20,7 @@ class AbsorbingConstraint { //! \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") + mpm::Position position = mpm::Position::None) : setid_{setid}, dir_{dir}, delta_{delta}, @@ -46,7 +48,7 @@ class AbsorbingConstraint { double b() const { return b_; } // Return position - std::string position() const { return position_; } + mpm::Position position() const { return position_; } private: // ID @@ -62,7 +64,7 @@ class AbsorbingConstraint { // b double b_; // Node position - std::string position_; + mpm::Position position_; }; } // namespace mpm #endif // MPM_ABSORBING_CONSTRAINT_H_ diff --git a/include/loads_bcs/constraints.h b/include/loads_bcs/constraints.h index 6f2da49a4..42c533863 100644 --- a/include/loads_bcs/constraints.h +++ b/include/loads_bcs/constraints.h @@ -62,7 +62,7 @@ class Constraints { //! b, and position bool assign_nodal_absorbing_constraints( const std::vector>& + double, mpm::Position>>& absorbing_constraints); private: diff --git a/include/loads_bcs/constraints.tcc b/include/loads_bcs/constraints.tcc index 3c1e0b055..eaf67445f 100644 --- a/include/loads_bcs/constraints.tcc +++ b/include/loads_bcs/constraints.tcc @@ -124,7 +124,7 @@ bool mpm::Constraints::assign_nodal_absorbing_constraint( double h_min = absorbing_constraint->h_min(); double a = absorbing_constraint->a(); double b = absorbing_constraint->b(); - std::string position = absorbing_constraint->position(); + mpm::Position 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, @@ -145,7 +145,8 @@ bool mpm::Constraints::assign_nodal_absorbing_constraint( template bool mpm::Constraints::assign_nodal_absorbing_constraints( const std::vector>& absorbing_constraints) { + double, mpm::Position>>& + absorbing_constraints) { bool status = true; try { for (const auto& absorbing_constraint : absorbing_constraints) { @@ -162,11 +163,12 @@ bool mpm::Constraints::assign_nodal_absorbing_constraints( // b double b = std::get<5>(absorbing_constraint); // Position - std::string position = std::get<6>(absorbing_constraint); + mpm::Position 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") { + if (position == mpm::Position::Side or + position == mpm::Position::Corner or + position == mpm::Position::Bottom) { // Apply constraint if (!mesh_->node(nid)->apply_absorbing_constraint(dir, delta, h_min, a, b, position)) diff --git a/include/node.h b/include/node.h index ae80e19f7..9c49d2a76 100644 --- a/include/node.h +++ b/include/node.h @@ -223,7 +223,7 @@ class Node : public NodeBase { //! \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; + mpm::Position position) override; //! Assign rotation matrix //! \param[in] rotation_matrix Rotation matrix of the node diff --git a/include/node.tcc b/include/node.tcc index be2a106b5..7fde9ff13 100644 --- a/include/node.tcc +++ b/include/node.tcc @@ -350,7 +350,7 @@ void mpm::Node::apply_velocity_constraints() { template bool mpm::Node::apply_absorbing_constraint( unsigned dir, double delta, double h_min, double a, double b, - std::string position) { + mpm::Position position) { bool status = true; try { @@ -395,21 +395,24 @@ bool mpm::Node::apply_absorbing_constraint( // Update external force if (Tdim == 2) { Eigen::Matrix 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 + double force0 = 0.0; + double force1 = 0.0; + if (position == mpm::Position::Corner) { + force0 = 0.5 * h_min * absorbing_traction_(0); + force1 = 0.5 * h_min * absorbing_traction_(1); + } else if (position == mpm::Position::Bottom) { + force0 = 0.5 * h_min * absorbing_traction_(0); + force1 = 1.0 * h_min * absorbing_traction_(1); + } else if (position == mpm::Position::Side) { + force0 = 1.0 * h_min * absorbing_traction_(0); + force1 = 0.5 * h_min * absorbing_traction_(1); + } else { throw std::runtime_error("Invalid absorbing boundary position"); + } + // clang-format off + absorbing_force_ << force0, force1; + // clang-format on + this->update_external_force(true, phase, -absorbing_force_); } } } diff --git a/include/node_base.h b/include/node_base.h index fc0defbdf..7683450fa 100644 --- a/include/node_base.h +++ b/include/node_base.h @@ -215,7 +215,7 @@ class NodeBase { //! \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; + mpm::Position position) = 0; //! Assign rotation matrix //! \param[in] rotation_matrix Rotation matrix of the node diff --git a/include/solvers/mpm_base.h b/include/solvers/mpm_base.h index e57b60f0a..52f6b6575 100644 --- a/include/solvers/mpm_base.h +++ b/include/solvers/mpm_base.h @@ -220,6 +220,7 @@ class MPMBase : public MPM { bool absorbing_boundary_{false}; std::vector> absorbing_constraint_; std::vector absorbing_nset_id_; + mpm::Position position_{mpm::Position::None}; #ifdef USE_GRAPH_PARTITIONING // graph pass the address of the container of cell diff --git a/include/solvers/mpm_base.tcc b/include/solvers/mpm_base.tcc index 03f985be5..b99a2e428 100644 --- a/include/solvers/mpm_base.tcc +++ b/include/solvers/mpm_base.tcc @@ -967,9 +967,17 @@ void mpm::MPMBase::nodal_absorbing_constraints( // position std::string position = constraints.at("position").template get(); + if (position == "side") + position_ = mpm::Position::Side; + else if (position == "corner") + position_ = mpm::Position::Corner; + else if (position == "bottom") + position_ = mpm::Position::Bottom; + else + position_ = mpm::Position::None; // Add absorbing constraint to mesh auto absorbing_constraint = std::make_shared( - nset_id, dir, delta, h_min, a, b, position); + nset_id, dir, delta, h_min, a, b, position_); bool absorbing_constraints = constraints_->assign_nodal_absorbing_constraint( nset_id, absorbing_constraint); From fd06771a6e9d16204304854cc7e4aabe0821b020 Mon Sep 17 00:00:00 2001 From: jgiven100 Date: Thu, 19 Aug 2021 14:41:11 -0700 Subject: [PATCH 12/19] :dart::heavy_check_mark: update absorb boundary tests --- tests/mesh_test_2d.cc | 28 ++++++++++++++++++++-------- 1 file changed, 20 insertions(+), 8 deletions(-) diff --git a/tests/mesh_test_2d.cc b/tests/mesh_test_2d.cc index 3c3bfdbd1..ea0e11dcb 100644 --- a/tests/mesh_test_2d.cc +++ b/tests/mesh_test_2d.cc @@ -1264,8 +1264,8 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") { // Test assign absorbing constraints to nodes SECTION("Check assign absorbing constraints to nodes") { // Vector of particle coordinates - std::vector< - std::tuple> + std::vector> absorbing_constraints; //! Constraints object auto constraints = std::make_shared>(mesh); @@ -1281,21 +1281,33 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") { } // Constraint - absorbing_constraints.emplace_back(std::make_tuple(0, 0, 1, 3, 2, 2)); - absorbing_constraints.emplace_back(std::make_tuple(1, 1, 2, 4, 1, 1)); - absorbing_constraints.emplace_back(std::make_tuple(2, 0, 1, 1, 1, 1)); - absorbing_constraints.emplace_back(std::make_tuple(3, 1, 3, 2, 3, 3)); + absorbing_constraints.emplace_back( + std::make_tuple(0, 0, 1, 3, 2, 2, mpm::Position::Side)); + absorbing_constraints.emplace_back( + std::make_tuple(1, 1, 2, 4, 1, 1, mpm::Position::Side)); + absorbing_constraints.emplace_back( + std::make_tuple(2, 0, 1, 1, 1, 1, mpm::Position::Side)); + absorbing_constraints.emplace_back( + std::make_tuple(3, 1, 3, 2, 3, 3, mpm::Position::Side)); REQUIRE(constraints->assign_nodal_absorbing_constraints( absorbing_constraints) == true); // When constraints fail: invalid direction - absorbing_constraints.emplace_back(std::make_tuple(3, 2, 3, 2, 3, 3)); + absorbing_constraints.emplace_back( + std::make_tuple(3, 2, 3, 2, 3, 3, mpm::Position::Side)); REQUIRE(constraints->assign_nodal_absorbing_constraints( absorbing_constraints) == false); // When constraints fail: invalid delta - absorbing_constraints.emplace_back(std::make_tuple(3, 1, 1, 3, 1, 1)); + absorbing_constraints.emplace_back( + std::make_tuple(3, 1, 1, 3, 1, 1, mpm::Position::Side)); + REQUIRE(constraints->assign_nodal_absorbing_constraints( + absorbing_constraints) == false); + + // When constraints fail: invalid position + absorbing_constraints.emplace_back( + std::make_tuple(0, 0, 1, 3, 2, 2, mpm::Position::None)); REQUIRE(constraints->assign_nodal_absorbing_constraints( absorbing_constraints) == false); } From f5f2d1b451c8e0f7661eeaa5374c53ed9952c709 Mon Sep 17 00:00:00 2001 From: jgiven100 Date: Thu, 19 Aug 2021 16:46:09 -0700 Subject: [PATCH 13/19] :hammer::bug: replace if with switch statement and cppcheck suppress --- include/node.tcc | 28 +++++++++++++++++----------- 1 file changed, 17 insertions(+), 11 deletions(-) diff --git a/include/node.tcc b/include/node.tcc index 7fde9ff13..9bb75fe48 100644 --- a/include/node.tcc +++ b/include/node.tcc @@ -397,19 +397,25 @@ bool mpm::Node::apply_absorbing_constraint( Eigen::Matrix absorbing_force_; double force0 = 0.0; double force1 = 0.0; - if (position == mpm::Position::Corner) { - force0 = 0.5 * h_min * absorbing_traction_(0); - force1 = 0.5 * h_min * absorbing_traction_(1); - } else if (position == mpm::Position::Bottom) { - force0 = 0.5 * h_min * absorbing_traction_(0); - force1 = 1.0 * h_min * absorbing_traction_(1); - } else if (position == mpm::Position::Side) { - force0 = 1.0 * h_min * absorbing_traction_(0); - force1 = 0.5 * h_min * absorbing_traction_(1); - } else { - throw std::runtime_error("Invalid absorbing boundary position"); + switch (position) { + case mpm::Position::Corner: + force0 = 0.5 * h_min * absorbing_traction_(0); + force1 = 0.5 * h_min * absorbing_traction_(1); + break; + case mpm::Position::Bottom: + force0 = 0.5 * h_min * absorbing_traction_(0); + force1 = 1.0 * h_min * absorbing_traction_(1); + break; + case mpm::Position::Side: + force0 = 1.0 * h_min * absorbing_traction_(0); + force1 = 0.5 * h_min * absorbing_traction_(1); + break; + default: + throw std::runtime_error("Invalid absorbing boundary position"); + break; } // clang-format off + // cppcheck-suppress * absorbing_force_ << force0, force1; // clang-format on this->update_external_force(true, phase, -absorbing_force_); From 805626cf728f524c7704d72e1e876ac289f00190 Mon Sep 17 00:00:00 2001 From: jgiven100 Date: Mon, 20 Sep 2021 18:18:10 -0700 Subject: [PATCH 14/19] :bug::hammer: fix traction -> force for absorb boundary --- include/data_types.h | 5 ++--- include/loads_bcs/constraints.tcc | 3 +-- include/node.tcc | 14 ++------------ include/solvers/mpm_base.tcc | 2 -- 4 files changed, 5 insertions(+), 19 deletions(-) mode change 100644 => 100755 include/node.tcc diff --git a/include/data_types.h b/include/data_types.h index 21eca99c7..f73bf7bdd 100644 --- a/include/data_types.h +++ b/include/data_types.h @@ -36,10 +36,9 @@ inline double zero() { //! Position type //! None: No position is specified -//! Side: Nodes along model vertical sides (left or right) +//! Side: Nodes along model vertical or horizontal sides //! Corner: Nodes at model corner -//! Bottom: Nodes along model bottom -enum class Position { None, Side, Corner, Bottom }; +enum class Position { None, Side, Corner }; } // namespace mpm diff --git a/include/loads_bcs/constraints.tcc b/include/loads_bcs/constraints.tcc index eaf67445f..194474d46 100644 --- a/include/loads_bcs/constraints.tcc +++ b/include/loads_bcs/constraints.tcc @@ -167,8 +167,7 @@ bool mpm::Constraints::assign_nodal_absorbing_constraints( // delta check if (delta >= h_min / (2 * a) and delta >= h_min / (2 * b)) { if (position == mpm::Position::Side or - position == mpm::Position::Corner or - position == mpm::Position::Bottom) { + position == mpm::Position::Corner) { // Apply constraint if (!mesh_->node(nid)->apply_absorbing_constraint(dir, delta, h_min, a, b, position)) diff --git a/include/node.tcc b/include/node.tcc old mode 100644 new mode 100755 index 9bb75fe48..d63377cc8 --- a/include/node.tcc +++ b/include/node.tcc @@ -399,25 +399,15 @@ bool mpm::Node::apply_absorbing_constraint( double force1 = 0.0; switch (position) { case mpm::Position::Corner: - force0 = 0.5 * h_min * absorbing_traction_(0); - force1 = 0.5 * h_min * absorbing_traction_(1); - break; - case mpm::Position::Bottom: - force0 = 0.5 * h_min * absorbing_traction_(0); - force1 = 1.0 * h_min * absorbing_traction_(1); + absorbing_force_ = 0.5 * h_min * absorbing_traction_; break; case mpm::Position::Side: - force0 = 1.0 * h_min * absorbing_traction_(0); - force1 = 0.5 * h_min * absorbing_traction_(1); + absorbing_force_ = h_min * absorbing_traction_; break; default: throw std::runtime_error("Invalid absorbing boundary position"); break; } - // clang-format off - // cppcheck-suppress * - absorbing_force_ << force0, force1; - // clang-format on this->update_external_force(true, phase, -absorbing_force_); } } diff --git a/include/solvers/mpm_base.tcc b/include/solvers/mpm_base.tcc index b99a2e428..b1218da4e 100644 --- a/include/solvers/mpm_base.tcc +++ b/include/solvers/mpm_base.tcc @@ -971,8 +971,6 @@ void mpm::MPMBase::nodal_absorbing_constraints( position_ = mpm::Position::Side; else if (position == "corner") position_ = mpm::Position::Corner; - else if (position == "bottom") - position_ = mpm::Position::Bottom; else position_ = mpm::Position::None; // Add absorbing constraint to mesh From 229e1d7cbb6abd203a701918e0fc301c8e5063d3 Mon Sep 17 00:00:00 2001 From: jgiven100 Date: Wed, 24 Nov 2021 11:07:51 -0800 Subject: [PATCH 15/19] :wrench::sparkles: add 3d support for Kelvin-Voigt --- include/data_types.h | 7 ++++--- include/loads_bcs/absorbing_constraint.h | 8 ++++---- include/loads_bcs/constraints.tcc | 5 +++-- include/node.h | 8 ++++---- include/node.tcc | 21 ++++++++++++++++++--- include/node_base.h | 8 ++++---- include/solvers/mpm_base.tcc | 8 +++++--- tests/mesh_test_2d.cc | 12 ++++++------ 8 files changed, 48 insertions(+), 29 deletions(-) mode change 100755 => 100644 include/node.tcc diff --git a/include/data_types.h b/include/data_types.h index f73bf7bdd..04ea207b9 100644 --- a/include/data_types.h +++ b/include/data_types.h @@ -36,9 +36,10 @@ inline double zero() { //! Position type //! None: No position is specified -//! Side: Nodes along model vertical or horizontal sides -//! Corner: Nodes at model corner -enum class Position { None, Side, Corner }; +//! Corner: Nodes at boundary corners +//! Edge: Nodes along boundary edges +//! Face: Nodes on boundary faces +enum class Position { None, Corner, Edge, Face }; } // namespace mpm diff --git a/include/loads_bcs/absorbing_constraint.h b/include/loads_bcs/absorbing_constraint.h index 84c7dd2c6..48503c9a1 100644 --- a/include/loads_bcs/absorbing_constraint.h +++ b/include/loads_bcs/absorbing_constraint.h @@ -12,11 +12,11 @@ class AbsorbingConstraint { public: // Constructor //! \param[in] setid set id - //! \param[in] dir Normal direction of boundary application + //! \param[in] dir Direction of p-wave propagation in model //! \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] h_min Characteristic length (cell height) + //! \param[in] a Dimensionless dashpot weight factor, p-wave + //! \param[in] b Dimensionless dashpot weight factor, s-wave //! \param[in] position Nodal position along boundary AbsorbingConstraint(int setid, unsigned dir, double delta, double h_min, double a = 1., double b = 1., diff --git a/include/loads_bcs/constraints.tcc b/include/loads_bcs/constraints.tcc index 194474d46..3c66f5130 100644 --- a/include/loads_bcs/constraints.tcc +++ b/include/loads_bcs/constraints.tcc @@ -166,8 +166,9 @@ bool mpm::Constraints::assign_nodal_absorbing_constraints( mpm::Position position = std::get<6>(absorbing_constraint); // delta check if (delta >= h_min / (2 * a) and delta >= h_min / (2 * b)) { - if (position == mpm::Position::Side or - position == mpm::Position::Corner) { + if (position == mpm::Position::Corner or + position == mpm::Position::Edge or + position == mpm::Position::Face) { // Apply constraint if (!mesh_->node(nid)->apply_absorbing_constraint(dir, delta, h_min, a, b, position)) diff --git a/include/node.h b/include/node.h index 9c49d2a76..7734d26c0 100644 --- a/include/node.h +++ b/include/node.h @@ -215,11 +215,11 @@ class Node : public NodeBase { //! Apply absorbing constraint //! Directions can take values between 0 and Dim * Nphases - //! \param[in] dir Direction of absorbing constraint + //! \param[in] dir Direction of p-wave propagation in model //! \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] h_min Characteristic length (cell height) + //! \param[in] a Dimensionless dashpot weight factor, p-wave + //! \param[in] b Dimensionless dashpot weight factor, s-wave //! \param[in] position Nodal position along boundary bool apply_absorbing_constraint(unsigned dir, double delta, double h_min, double a, double b, diff --git a/include/node.tcc b/include/node.tcc old mode 100755 new mode 100644 index d63377cc8..ec883f928 --- a/include/node.tcc +++ b/include/node.tcc @@ -395,13 +395,11 @@ bool mpm::Node::apply_absorbing_constraint( // Update external force if (Tdim == 2) { Eigen::Matrix absorbing_force_; - double force0 = 0.0; - double force1 = 0.0; switch (position) { case mpm::Position::Corner: absorbing_force_ = 0.5 * h_min * absorbing_traction_; break; - case mpm::Position::Side: + case mpm::Position::Edge: absorbing_force_ = h_min * absorbing_traction_; break; default: @@ -409,6 +407,23 @@ bool mpm::Node::apply_absorbing_constraint( break; } this->update_external_force(true, phase, -absorbing_force_); + } else if (Tdim == 3) { + Eigen::Matrix absorbing_force_; + switch (position) { + case mpm::Position::Corner: + absorbing_force_ = 0.25 * pow(h_min, 2) * absorbing_traction_; + break; + case mpm::Position::Edge: + absorbing_force_ = 0.5 * pow(h_min, 2) * absorbing_traction_; + break; + case mpm::Position::Face: + absorbing_force_ = pow(h_min, 2) * absorbing_traction_; + break; + default: + throw std::runtime_error("Invalid absorbing boundary position"); + break; + } + this->update_external_force(true, phase, -absorbing_force_); } } } diff --git a/include/node_base.h b/include/node_base.h index 7683450fa..f585e8bc0 100644 --- a/include/node_base.h +++ b/include/node_base.h @@ -207,11 +207,11 @@ class NodeBase { //! Apply absorbing constraint //! Directions can take values between 0 and Dim * Nphases - //! \param[in] dir Direction of absorbing constraint + //! \param[in] dir Direction of p-wave propagation in model //! \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] h_min Characteristic length (cell height) + //! \param[in] a Dimensionless dashpot weight factor, p-wave + //! \param[in] b Dimensionless dashpot weight factor, s-wave //! \param[in] position Nodal position along boundary virtual bool apply_absorbing_constraint(unsigned dir, double delta, double h_min, double a, double b, diff --git a/include/solvers/mpm_base.tcc b/include/solvers/mpm_base.tcc index b1218da4e..d71d53df9 100644 --- a/include/solvers/mpm_base.tcc +++ b/include/solvers/mpm_base.tcc @@ -967,10 +967,12 @@ void mpm::MPMBase::nodal_absorbing_constraints( // position std::string position = constraints.at("position").template get(); - if (position == "side") - position_ = mpm::Position::Side; - else if (position == "corner") + if (position == "corner") position_ = mpm::Position::Corner; + else if (position == "edge") + position_ = mpm::Position::Edge; + else if (position == "face") + position_ = mpm::Position::Face; else position_ = mpm::Position::None; // Add absorbing constraint to mesh diff --git a/tests/mesh_test_2d.cc b/tests/mesh_test_2d.cc index ea0e11dcb..6fa47e16d 100644 --- a/tests/mesh_test_2d.cc +++ b/tests/mesh_test_2d.cc @@ -1282,26 +1282,26 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") { // Constraint absorbing_constraints.emplace_back( - std::make_tuple(0, 0, 1, 3, 2, 2, mpm::Position::Side)); + std::make_tuple(0, 0, 1, 3, 2, 2, mpm::Position::Edge)); absorbing_constraints.emplace_back( - std::make_tuple(1, 1, 2, 4, 1, 1, mpm::Position::Side)); + std::make_tuple(1, 1, 2, 4, 1, 1, mpm::Position::Edge)); absorbing_constraints.emplace_back( - std::make_tuple(2, 0, 1, 1, 1, 1, mpm::Position::Side)); + std::make_tuple(2, 0, 1, 1, 1, 1, mpm::Position::Edge)); absorbing_constraints.emplace_back( - std::make_tuple(3, 1, 3, 2, 3, 3, mpm::Position::Side)); + std::make_tuple(3, 1, 3, 2, 3, 3, mpm::Position::Edge)); REQUIRE(constraints->assign_nodal_absorbing_constraints( absorbing_constraints) == true); // When constraints fail: invalid direction absorbing_constraints.emplace_back( - std::make_tuple(3, 2, 3, 2, 3, 3, mpm::Position::Side)); + std::make_tuple(3, 2, 3, 2, 3, 3, mpm::Position::Edge)); REQUIRE(constraints->assign_nodal_absorbing_constraints( absorbing_constraints) == false); // When constraints fail: invalid delta absorbing_constraints.emplace_back( - std::make_tuple(3, 1, 1, 3, 1, 1, mpm::Position::Side)); + std::make_tuple(3, 1, 1, 3, 1, 1, mpm::Position::Edge)); REQUIRE(constraints->assign_nodal_absorbing_constraints( absorbing_constraints) == false); From 746283d9930674a97cd51522a4a8d5fdc0addefd Mon Sep 17 00:00:00 2001 From: jgiven100 Date: Thu, 25 Nov 2021 10:26:41 -0800 Subject: [PATCH 16/19] :construction::dart: improve test coverage --- tests/mesh_test_2d.cc | 31 +++++++++++++++++++++---------- 1 file changed, 21 insertions(+), 10 deletions(-) diff --git a/tests/mesh_test_2d.cc b/tests/mesh_test_2d.cc index 6fa47e16d..f25b20fa9 100644 --- a/tests/mesh_test_2d.cc +++ b/tests/mesh_test_2d.cc @@ -1150,10 +1150,12 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") { set_id, friction_constraint) == false); } + // Test assign absorbing constraint SECTION("Check assign absorbing constraints") { tsl::robin_map> node_sets; node_sets[0] = std::vector{0, 2}; node_sets[1] = std::vector{1, 3}; + node_sets[2] = std::vector{}; REQUIRE(mesh->create_node_sets(node_sets, true) == true); @@ -1168,6 +1170,7 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") { for (double i = 0; i < 4; ++i) { REQUIRE_NOTHROW( mesh->node(i)->initialise_property_handle(0, nodal_properties)); + mesh->node(i)->append_material_id(0); } int set_id = 0; @@ -1176,11 +1179,12 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") { double h_min = 1; double a = 1; double b = 1; + mpm::Position pos = mpm::Position::Edge; // Add absorbing constraint to mesh auto absorbing_constraint = std::make_shared(set_id, dir, delta, - h_min, a, b); + h_min, a, b, pos); REQUIRE(constraints->assign_nodal_absorbing_constraint( set_id, absorbing_constraint) == true); @@ -1192,13 +1196,7 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") { b = 2; // Add absorbing constraint to mesh absorbing_constraint = std::make_shared( - set_id, dir, delta, h_min, a, b); - REQUIRE(constraints->assign_nodal_absorbing_constraint( - set_id, absorbing_constraint) == true); - - // Add absorbing constraint to all nodes in mesh - absorbing_constraint = std::make_shared( - set_id, dir, delta, h_min, a, b); + set_id, dir, delta, h_min, a, b, pos); REQUIRE(constraints->assign_nodal_absorbing_constraint( set_id, absorbing_constraint) == true); @@ -1206,15 +1204,25 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") { dir = 3; // Add absorbing constraint to mesh absorbing_constraint = std::make_shared( - set_id, dir, delta, h_min, a, b); + set_id, dir, delta, h_min, a, b, pos); REQUIRE(constraints->assign_nodal_absorbing_constraint( set_id, absorbing_constraint) == false); // When constraints fail: invalid delta + dir = 1; delta = 0; // Add absorbing constraint to mesh absorbing_constraint = std::make_shared( - set_id, dir, delta, h_min, a, b); + set_id, dir, delta, h_min, a, b, pos); + REQUIRE(constraints->assign_nodal_absorbing_constraint( + set_id, absorbing_constraint) == false); + + // When constraints fail: empty node set + delta = 3; + set_id = 2; + // Add absorbing constraint to mesh + absorbing_constraint = std::make_shared( + set_id, dir, delta, h_min, a, b, pos); REQUIRE(constraints->assign_nodal_absorbing_constraint( set_id, absorbing_constraint) == false); } @@ -1294,18 +1302,21 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") { absorbing_constraints) == true); // When constraints fail: invalid direction + absorbing_constraints.clear(); absorbing_constraints.emplace_back( std::make_tuple(3, 2, 3, 2, 3, 3, mpm::Position::Edge)); REQUIRE(constraints->assign_nodal_absorbing_constraints( absorbing_constraints) == false); // When constraints fail: invalid delta + absorbing_constraints.clear(); absorbing_constraints.emplace_back( std::make_tuple(3, 1, 1, 3, 1, 1, mpm::Position::Edge)); REQUIRE(constraints->assign_nodal_absorbing_constraints( absorbing_constraints) == false); // When constraints fail: invalid position + absorbing_constraints.clear(); absorbing_constraints.emplace_back( std::make_tuple(0, 0, 1, 3, 2, 2, mpm::Position::None)); REQUIRE(constraints->assign_nodal_absorbing_constraints( From 397e6e9784f9234d2a57a9c9c5cd632023ee1e4e Mon Sep 17 00:00:00 2001 From: jgiven100 Date: Thu, 25 Nov 2021 11:29:20 -0800 Subject: [PATCH 17/19] :bug: clang-format --- tests/mesh_test_2d.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/mesh_test_2d.cc b/tests/mesh_test_2d.cc index f25b20fa9..437a4cb76 100644 --- a/tests/mesh_test_2d.cc +++ b/tests/mesh_test_2d.cc @@ -1150,7 +1150,7 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") { set_id, friction_constraint) == false); } - // Test assign absorbing constraint + // Test assign absorbing constraint SECTION("Check assign absorbing constraints") { tsl::robin_map> node_sets; node_sets[0] = std::vector{0, 2}; From 6fa65b98b39ac50e11e580d62e4fc26b5f6f68a9 Mon Sep 17 00:00:00 2001 From: jgiven100 Date: Thu, 25 Nov 2021 12:36:48 -0800 Subject: [PATCH 18/19] :construction::dart: improve test coverage --- tests/mesh_test_2d.cc | 10 +++++++ tests/mesh_test_3d.cc | 61 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 71 insertions(+) diff --git a/tests/mesh_test_2d.cc b/tests/mesh_test_2d.cc index 437a4cb76..3f93ac785 100644 --- a/tests/mesh_test_2d.cc +++ b/tests/mesh_test_2d.cc @@ -1194,6 +1194,7 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") { h_min = 12; a = 2; b = 2; + pos = mpm::Position::Corner; // Add absorbing constraint to mesh absorbing_constraint = std::make_shared( set_id, dir, delta, h_min, a, b, pos); @@ -1225,6 +1226,15 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") { set_id, dir, delta, h_min, a, b, pos); REQUIRE(constraints->assign_nodal_absorbing_constraint( set_id, absorbing_constraint) == false); + + // When constraints fail: invalid absorbing boundary position + set_id = 1; + pos = mpm::Position::None; + // Add absorbing constraint to mesh + absorbing_constraint = std::make_shared( + set_id, dir, delta, h_min, a, b, pos); + REQUIRE(constraints->assign_nodal_absorbing_constraint( + set_id, absorbing_constraint) == false); } // Test assign velocity constraints to nodes diff --git a/tests/mesh_test_3d.cc b/tests/mesh_test_3d.cc index 06f4e0cf2..e6698f764 100644 --- a/tests/mesh_test_3d.cc +++ b/tests/mesh_test_3d.cc @@ -1261,6 +1261,67 @@ TEST_CASE("Mesh is checked for 3D case", "[mesh][3D]") { set_id, friction_constraint) == false); } + // Test assign absorbing constraint + SECTION("Check assign absorbing constraints") { + tsl::robin_map> node_sets; + node_sets[0] = std::vector{0, 2}; + node_sets[1] = std::vector{1, 3}; + node_sets[2] = std::vector{}; + + REQUIRE(mesh->create_node_sets(node_sets, true) == true); + + //! Constraints object + auto constraints = std::make_shared>(mesh); + + // Assign nodal properties + auto nodal_properties = std::make_shared(); + nodal_properties->create_property("density", 1, 1); + nodal_properties->create_property("wave_velocities", Dim, 1); + nodal_properties->create_property("displacements", Dim, 1); + for (double i = 0; i < 4; ++i) { + REQUIRE_NOTHROW( + mesh->node(i)->initialise_property_handle(0, nodal_properties)); + mesh->node(i)->append_material_id(0); + } + + int set_id = 0; + unsigned dir = 0; + double delta = 1; + double h_min = 1; + double a = 1; + double b = 1; + mpm::Position pos = mpm::Position::Edge; + + // Add absorbing constraint to mesh + auto absorbing_constraint = + std::make_shared(set_id, dir, delta, + h_min, a, b, pos); + REQUIRE(constraints->assign_nodal_absorbing_constraint( + set_id, absorbing_constraint) == true); + + pos = mpm::Position::Corner; + // Add absorbing constraint to mesh + absorbing_constraint = std::make_shared( + set_id, dir, delta, h_min, a, b, pos); + REQUIRE(constraints->assign_nodal_absorbing_constraint( + set_id, absorbing_constraint) == true); + + pos = mpm::Position::Face; + // Add absorbing constraint to mesh + absorbing_constraint = std::make_shared( + set_id, dir, delta, h_min, a, b, pos); + REQUIRE(constraints->assign_nodal_absorbing_constraint( + set_id, absorbing_constraint) == true); + + // When constraints fail: invalid absorbing boundary position + pos = mpm::Position::None; + // Add absorbing constraint to mesh + absorbing_constraint = std::make_shared( + set_id, dir, delta, h_min, a, b, pos); + REQUIRE(constraints->assign_nodal_absorbing_constraint( + set_id, absorbing_constraint) == false); + } + // Test assign velocity constraints to nodes SECTION("Check assign velocity constraints to nodes") { // Vector of particle coordinates From 4cbc60a975ea52d513c9829fde183d81c101f46b Mon Sep 17 00:00:00 2001 From: jgiven100 Date: Thu, 25 Nov 2021 14:19:49 -0800 Subject: [PATCH 19/19] :dart::fire: tidy commented-out lines --- tests/node_test.cc | 83 ---------------------------------------------- 1 file changed, 83 deletions(-) diff --git a/tests/node_test.cc b/tests/node_test.cc index b06094e95..3cc92d2a5 100644 --- a/tests/node_test.cc +++ b/tests/node_test.cc @@ -726,89 +726,6 @@ TEST_CASE("Node is checked for 2D case", "[node][2D]") { REQUIRE(node->external_force(Nphase)(i) == Approx(0.).epsilon(Tolerance)); } - - /* SECTION("Check absorbing boundary") { - // Create a force vector - Eigen::Matrix force; - for (unsigned i = 0; i < force.size(); ++i) force(i) = 10. * - (2 - i); REQUIRE_NOTHROW(node->update_external_force(false, Nphase, - force)); - - // Assign volume to 100 - double volume = 100.; - REQUIRE_NOTHROW(node->update_volume(false, Nphase, volume)); - REQUIRE(node->volume(Nphase) == - Approx(100.0).epsilon(Tolerance)); - - // Assign mass to 100 - double mass = 100.; - REQUIRE_NOTHROW(node->update_mass(false, Nphase, mass)); - REQUIRE(node->mass(Nphase) == - Approx(100.0).epsilon(Tolerance)); - - // Compute velocity - double dt = 0.1; - REQUIRE(node->compute_acceleration_velocity(Nphase, dt) == - true); - - // Check velocity - Eigen::Matrix velocity = force / mass * dt; - for (unsigned i = 0; i < velocity.size(); ++i) - REQUIRE(node->velocity(Nphase)(i) == - Approx(velocity(i)).epsilon(Tolerance)); - - // Define nodal properties - Eigen::Matrix density; - density(0) = mass / volume; - Eigen::Matrix wave_v; - for (double i = 0; i < wave_v.size(); ++i) wave_v(i) = 100 * - (2 - i); Eigen::Matrix displacements; for (double i = - 0; i < displacements.size(); ++i) displacements(i) = 0.0001 * (2 - - i); - - // Define material id - unsigned mat_id = 0; - node->append_material_id(mat_id); - - // Assign nodal properties - auto nodal_properties = - std::make_shared(); - nodal_properties->create_property("density", 1, 1); - nodal_properties->create_property("wave_velocities", Dim, 1); - nodal_properties->create_property("displacements", Dim, 1); - REQUIRE_NOTHROW( - node->initialise_property_handle(0, nodal_properties)); - REQUIRE_NOTHROW( - node->update_property(false, "density", density, mat_id, - 1)); REQUIRE_NOTHROW(node->update_property(false, "wave_velocities", - wave_v, mat_id, Dim)); - REQUIRE_NOTHROW(node->update_property(false, "displacements", - displacements, mat_id, - Dim)); - - // Check invalid direction - unsigned bad_dir = 4; - REQUIRE(node->apply_absorbing_constraint(bad_dir, 0, 1, 1, 1) - == false); - - // Apply absorbing constraint - unsigned dir = 0; - REQUIRE(node->apply_absorbing_constraint(dir, 1, 1, 1, 1) == - true); - - // Calculation of absorbing force - double pwave_traction = -12.; - double swave_traction = -2.; - for (unsigned i = 0; i < Dim; ++i) { - if (i == dir) - REQUIRE(node->external_force(Nphase)(i) == - Approx(pwave_traction + - force(i)).epsilon(Tolerance)); else - REQUIRE(node->external_force(Nphase)(i) == - Approx(swave_traction + - force(i)).epsilon(Tolerance)); - } - }*/ } }