Apply ground motion input as nodal acceleration constraint

CMakeLists.txt

@@ -185,6 +185,7 @@ target_link_libraries(mpm git)
 if(MPM_BUILD_TESTING)
   SET(test_src
     ${mpm_SOURCE_DIR}/tests/test_main.cc
+    ${mpm_SOURCE_DIR}/tests/acceleration_constraint_test.cc
     ${mpm_SOURCE_DIR}/tests/cell_test.cc
     ${mpm_SOURCE_DIR}/tests/cell_vector_test.cc
     ${mpm_SOURCE_DIR}/tests/contact_test.cc

external/csv_license.md

@@ -0,0 +1,27 @@
+Copyright (c) 2015, ben-strasser
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+* Redistributions of source code must retain the above copyright notice, this
+  list of conditions and the following disclaimer.
+
+* Redistributions in binary form must reproduce the above copyright notice,
+  this list of conditions and the following disclaimer in the documentation
+  and/or other materials provided with the distribution.
+
+* Neither the name of fast-cpp-csv-parser nor the names of its
+  contributors may be used to endorse or promote products derived from
+  this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

include/io/io_mesh.h

@@ -12,6 +12,7 @@
 #include "Eigen/Dense"
 #include <boost/algorithm/string.hpp>
 
+#include "csv/csv.h"
 #include "spdlog/spdlog.h"
 
 #include "data_types.h"
@@ -87,22 +88,33 @@ class IOMesh {
       const std::vector<std::array<mpm::Index, 2>>& particles_cells) = 0;
 
   //! Read velocity constraints file
-  //! \param[in] velocity_constraints_files file name with constraints
+  //! \param[in] velocity_constraints_file file name with constraints
   virtual std::vector<std::tuple<mpm::Index, unsigned, double>>
       read_velocity_constraints(
           const std::string& velocity_constraints_file) = 0;
 
+  //! Read acceleration constraints file
+  //! \param[in] acceleration_constraints_file file name with constraints
+  virtual std::vector<std::tuple<mpm::Index, unsigned, double>>
+      read_acceleration_constraints(
+          const std::string& acceleration_constraints_file) = 0;
+
   //! Read friction constraints file
-  //! \param[in] friction_constraints_files file name with frictions
+  //! \param[in] friction_constraints_file file name with frictions
   virtual std::vector<std::tuple<mpm::Index, unsigned, int, double>>
       read_friction_constraints(
           const std::string& friction_constraints_file) = 0;
 
   //! Read forces file
-  //! \param[in] forces_files file name with nodal concentrated force
+  //! \param[in] forces_file file name with nodal concentrated force
   virtual std::vector<std::tuple<mpm::Index, unsigned, double>> read_forces(
       const std::string& forces_file) = 0;
 
+  //! Read math function file
+  //! \param[in] function_file file name with linear math function entries
+  virtual std::array<std::vector<double>, 2> read_math_functions(
+      const std::string& math_file) = 0;
+
 };  // IOMesh class
 }  // namespace mpm
 

include/io/io_mesh_ascii.h

@@ -74,14 +74,20 @@ class IOMeshAscii : public IOMesh<Tdim> {
       const std::string& particles_cells_file,
       const std::vector<std::array<mpm::Index, 2>>& particles_cells) override;
 
-  //! Read constraints file
-  //! \param[in] velocity_constraints_files file name with constraints
+  //! Read velocity constraints file
+  //! \param[in] velocity_constraints_file file name with constraints
   std::vector<std::tuple<mpm::Index, unsigned, double>>
       read_velocity_constraints(
           const std::string& velocity_constraints_file) override;
 
+  //! Read acceleration constraints file
+  //! \param[in] acceleration_constraints_file file name with constraints
+  std::vector<std::tuple<mpm::Index, unsigned, double>>
+      read_acceleration_constraints(
+          const std::string& acceleration_constraints_file) override;
+
   //! Read friction constraints file
-  //! \param[in] friction_constraints_files file name with frictions
+  //! \param[in] friction_constraints_file file name with frictions
   std::vector<std::tuple<mpm::Index, unsigned, int, double>>
       read_friction_constraints(
           const std::string& friction_constraints_file) override;
@@ -91,6 +97,11 @@ class IOMeshAscii : public IOMesh<Tdim> {
   std::vector<std::tuple<mpm::Index, unsigned, double>> read_forces(
       const std::string& forces_file) override;
 
+  //! Read math function file
+  //! \param[in] function_file file name with linear math function entries
+  std::array<std::vector<double>, 2> read_math_functions(
+      const std::string& math_file) override;
+
  private:
   //! Logger
   std::shared_ptr<spdlog::logger> console_;

include/io/io_mesh_ascii.tcc

@@ -270,12 +270,12 @@ std::map<mpm::Index, Eigen::Matrix<double, Tdim, 1>>
         // ignore comment lines (# or !) or blank lines
         if ((line.find('#') == std::string::npos) &&
             (line.find('!') == std::string::npos) && (line != "")) {
+          // ID and read stream
+          mpm::Index id;
+          // Angles and ream stream
+          Eigen::Matrix<double, Tdim, 1> angles;
           while (istream.good()) {
-            // ID and read stream
-            mpm::Index id;
             istream >> id;
-            // Angles and ream stream
-            Eigen::Matrix<double, Tdim, 1> angles;
             for (unsigned i = 0; i < Tdim; ++i) istream >> angles[i];
             euler_angles.emplace(std::make_pair(id, angles));
           }
@@ -314,11 +314,11 @@ std::vector<std::tuple<mpm::Index, double>>
         // ignore comment lines (# or !) or blank lines
         if ((line.find('#') == std::string::npos) &&
             (line.find('!') == std::string::npos) && (line != "")) {
+          // ID
+          mpm::Index id;
+          // Volume
+          double volume;
           while (istream.good()) {
-            // ID
-            mpm::Index id;
-            // Volume
-            double volume;
             // Read stream
             istream >> id >> volume;
             volumes.emplace_back(std::make_tuple(id, volume));
@@ -358,9 +358,9 @@ std::vector<std::array<mpm::Index, 2>>
         // ignore comment lines (# or !) or blank lines
         if ((line.find('#') == std::string::npos) &&
             (line.find('!') == std::string::npos) && (line != "")) {
+          // ID
+          mpm::Index pid, cid;
           while (istream.good()) {
-            // ID
-            mpm::Index pid, cid;
             // Read stream
             istream >> pid >> cid;
             particles_cells.emplace_back(std::array<mpm::Index, 2>({pid, cid}));
@@ -416,13 +416,13 @@ std::vector<std::tuple<mpm::Index, unsigned, double>>
         // ignore comment lines (# or !) or blank lines
         if ((line.find('#') == std::string::npos) &&
             (line.find('!') == std::string::npos) && (line != "")) {
+          // ID
+          mpm::Index id;
+          // Direction
+          unsigned dir;
+          // Velocity
+          double velocity;
           while (istream.good()) {
-            // ID
-            mpm::Index id;
-            // Direction
-            unsigned dir;
-            // Velocity
-            double velocity;
             // Read stream
             istream >> id >> dir >> velocity;
             constraints.emplace_back(std::make_tuple(id, dir, velocity));
@@ -438,6 +438,52 @@ std::vector<std::tuple<mpm::Index, unsigned, double>>
   return constraints;
 }
 
+//! Return acceleration constraints of nodes or particles
+template <unsigned Tdim>
+std::vector<std::tuple<mpm::Index, unsigned, double>>
+    mpm::IOMeshAscii<Tdim>::read_acceleration_constraints(
+        const std::string& acceleration_constraints_file) {
+
+  // Nodal or particle acceleration constraints
+  std::vector<std::tuple<mpm::Index, unsigned, double>> constraints;
+  constraints.clear();
+
+  // input file stream
+  std::fstream file;
+  file.open(acceleration_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 != "")) {
+          // ID
+          mpm::Index id;
+          // Direction
+          unsigned dir;
+          // Velocity
+          double acceleration;
+          while (istream.good()) {
+            // Read stream
+            istream >> id >> dir >> acceleration;
+            constraints.emplace_back(std::make_tuple(id, dir, acceleration));
+          }
+        }
+      }
+    }
+    file.close();
+  } catch (std::exception& exception) {
+    console_->error("Read acceleration constraints: {}", exception.what());
+    file.close();
+  }
+  return constraints;
+}
+
 //! Return friction constraints of particles
 template <unsigned Tdim>
 std::vector<std::tuple<mpm::Index, unsigned, int, double>>
@@ -462,15 +508,15 @@ std::vector<std::tuple<mpm::Index, unsigned, int, double>>
         // ignore comment lines (# or !) or blank lines
         if ((line.find('#') == std::string::npos) &&
             (line.find('!') == std::string::npos) && (line != "")) {
+          // ID
+          mpm::Index id;
+          // Direction
+          unsigned dir;
+          // Sign
+          int sign;
+          // Friction
+          double friction;
           while (istream.good()) {
-            // ID
-            mpm::Index id;
-            // Direction
-            unsigned dir;
-            // Sign
-            int sign;
-            // Friction
-            double friction;
             // Read stream
             istream >> id >> dir >> sign >> friction;
             constraints.emplace_back(std::make_tuple(id, dir, sign, friction));
@@ -509,13 +555,13 @@ std::vector<std::tuple<mpm::Index, unsigned, double>>
         // ignore comment lines (# or !) or blank lines
         if ((line.find('#') == std::string::npos) &&
             (line.find('!') == std::string::npos) && (line != "")) {
+          // ID
+          mpm::Index id;
+          // Direction
+          unsigned dir;
+          // Force
+          double force;
           while (istream.good()) {
-            // ID
-            mpm::Index id;
-            // Direction
-            unsigned dir;
-            // Force
-            double force;
             // Read stream
             istream >> id >> dir >> force;
             forces.emplace_back(std::make_tuple(id, dir, force));
@@ -530,3 +576,25 @@ std::vector<std::tuple<mpm::Index, unsigned, double>>
   }
   return forces;
 }
+
+// Return array with math function entries
+template <unsigned Tdim>
+std::array<std::vector<double>, 2> mpm::IOMeshAscii<Tdim>::read_math_functions(
+    const std::string& math_file) {
+  // Initialise vector with 2 empty vectors
+  std::array<std::vector<double>, 2> xfx_values;
+
+  // Read from csv file
+  try {
+    io::CSVReader<2> in(math_file.c_str());
+    double x_value, fx_value;
+    while (in.read_row(x_value, fx_value)) {
+      xfx_values[0].push_back(x_value);
+      xfx_values[1].push_back(fx_value);
+    }
+  } catch (std::exception& exception) {
+    console_->error("Read math functions: {}", exception.what());
+  }
+
+  return xfx_values;
+}

include/loads_bcs/acceleration_constraint.h

@@ -0,0 +1,56 @@
+#ifndef MPM_ACCELERATION_CONSTRAINT_H_
+#define MPM_ACCELERATION_CONSTRAINT_H_
+
+//! Alias for JSON
+#include "json.hpp"
+using Json = nlohmann::json;
+
+#include "function_base.h"
+
+namespace mpm {
+
+//! AccelerationConstraint class to store acceleration constraint on a set
+//! \brief AccelerationConstraint class to store a constraint on a set
+//! \details AccelerationConstraint stores the constraint as a static value
+class AccelerationConstraint {
+ public:
+  // Constructor
+  //! \param[in] setid  set id
+  //! \param[in] acceleration_fn Math function if defined
+  //! \param[in] dir Direction of constraint load
+  //! \param[in] acceleration Constraint acceleration
+  AccelerationConstraint(
+      int setid, const std::shared_ptr<mpm::FunctionBase>& acceleration_fn,
+      unsigned dir, double acceleration)
+      : setid_{setid},
+        acceleration_fn_{acceleration_fn},
+        dir_{dir},
+        acceleration_{acceleration} {};
+
+  // Set id
+  int setid() const { return setid_; }
+
+  // Direction
+  unsigned dir() const { return dir_; }
+
+  // Return acceleration
+  double acceleration(double current_time) const {
+    // Static load when no math function is defined
+    double scalar = (this->acceleration_fn_ != nullptr)
+                        ? (this->acceleration_fn_)->value(current_time)
+                        : 1.0;
+    return acceleration_ * scalar;
+  }
+
+ private:
+  // ID
+  int setid_;
+  // Math function
+  std::shared_ptr<mpm::FunctionBase> acceleration_fn_;
+  // Direction
+  unsigned dir_;
+  // Acceleration
+  double acceleration_;
+};
+}  // namespace mpm
+#endif  // MPM_ACCELERATION_CONSTRAINT_H_

include/loads_bcs/constraints.h

@@ -3,6 +3,7 @@
 
 #include <memory>
 
+#include "acceleration_constraint.h"
 #include "friction_constraint.h"
 #include "logger.h"
 #include "mesh.h"
@@ -34,6 +35,21 @@ class Constraints {
       const std::vector<std::tuple<mpm::Index, unsigned, double>>&
           velocity_constraints);
 
+  //! Assign nodal acceleration constraints
+  //! \param[in] setid Node set id
+  //! \param[in] acceleration_constraints Accelerartion constraint at node, dir,
+  //! acceleration
+  bool assign_nodal_acceleration_constraint(
+      int set_id,
+      const std::shared_ptr<mpm::AccelerationConstraint>& constraint);
+
+  //! Assign acceleartion constraints to nodes
+  //! \param[in] acceleration_constraints Constraint at node, dir, and
+  //! acceleration
+  bool assign_nodal_acceleration_constraints(
+      const std::vector<std::tuple<mpm::Index, unsigned, double>>&
+          acceleration_constraints);
+
   //! Assign nodal frictional constraints
   //! \param[in] setid Node set id
   //! \param[in] friction_constraints Constraint at node, dir, sign, friction