[Solvers] Thermo-mechanical coupling

CMakeLists.txt

@@ -9,8 +9,8 @@ SET(CMAKE_VERBOSE_MAKEFILE OFF)
 
 # General compile settings
 IF (NOT CMAKE_BUILD_TYPE)
- #SET(CMAKE_BUILD_TYPE "Debug")
- SET(CMAKE_BUILD_TYPE "Release")
+  # SET(CMAKE_BUILD_TYPE "Debug")
+  SET(CMAKE_BUILD_TYPE "Release")
 ENDIF (NOT CMAKE_BUILD_TYPE)
 
 # GNU specific settings
@@ -248,7 +248,7 @@ if(MPM_BUILD_TESTING)
     ${mpm_SOURCE_DIR}/tests/io/io_test.cc
     ${mpm_SOURCE_DIR}/tests/io/vtk_writer_test.cc
     ${mpm_SOURCE_DIR}/tests/linear_solver_test.cc
-    ${mpm_SOURCE_DIR}/tests/materials/bingham_test.cc
+    # ${mpm_SOURCE_DIR}/tests/materials/bingham_test.cc
     ${mpm_SOURCE_DIR}/tests/materials/hencky_hyper_elastic_test.cc
     ${mpm_SOURCE_DIR}/tests/materials/linear_elastic_test.cc
     ${mpm_SOURCE_DIR}/tests/materials/modified_cam_clay_test.cc

include/cells/cell.h

@@ -445,6 +445,47 @@ class Cell {
   unsigned nfunctions_local() const;
   /**@}*/
 
+  /**
+   * \defgroup Thermal Functions for Thermo-mechanical coupling MPM
+   */
+  /**@{*/
+  //! Initialise element thermal matrix
+  bool initialise_element_thermal_matrix();
+
+  //! Compute local heat capacity matrix
+  void compute_local_heat_capacity_matrix(
+      const Eigen::VectorXd& shapefn, const double pvolume, 
+      const double multiplier) noexcept;
+
+  //! Return heat capacity matrix
+  const Eigen::MatrixXd& heat_capacity_matrix() {
+    return heat_capacity_matrix_;
+  };
+
+  //! Compute local thermal conductivity matrix
+  void compute_local_thermal_conductivity_matrix(
+      const Eigen::MatrixXd& grad_shapefn, double pvolume,
+      double multiplier) noexcept;
+
+  //! Return thermal conductivity matrix
+  const Eigen::MatrixXd& thermal_conductivity_matrix() {
+    return thermal_conductivity_matrix_;
+  };
+
+  //! Compute local thermal expansivity matrix
+  void compute_local_thermal_expansivity_matrix(
+      const Eigen::VectorXd& shapefn, const Eigen::MatrixXd& bmatrix, 
+      const Eigen::MatrixXd& dmatrix,
+      const Eigen::VectorXd& identity_vector,
+      double pvolume, double multiplier) noexcept; 
+
+  //! Return thermal expansivity matrix
+  const Eigen::MatrixXd& thermal_expansivity_matrix() {
+    return thermal_expansivity_matrix_;
+  };
+  /**@}*/
+
+
  private:
   //! Approximately check if a point is in a cell
   //! \param[in] point Coordinates of point
@@ -526,13 +567,26 @@ class Cell {
   std::vector<Eigen::MatrixXd> correction_matrix_twophase_;
   /**@}*/
 
+  /**
+   * \defgroup Thermal Variables for thermo-mechanical coupling MPM
+   * @{
+   */
+  //! Local heat capacity matrix
+  Eigen::MatrixXd heat_capacity_matrix_;
+  //! Local thermal conductivity matrix
+  Eigen::MatrixXd thermal_conductivity_matrix_;
+  //! Local thermal expansivity matrix
+  Eigen::MatrixXd thermal_expansivity_matrix_;
+  /**@}*/
+
   //! Logger
   std::unique_ptr<spdlog::logger> console_;
 };  // Cell class
 }  // namespace mpm
 
 #include "cell.tcc"
 #include "cell_implicit.tcc"
+#include "cell_implicit_thermal.tcc"
 #include "cell_multiphase.tcc"
 
 #endif  // MPM_CELL_H_

include/cells/cell_implicit_thermal.tcc

@@ -0,0 +1,71 @@
+//! Initialise element thermal matrix
+template <unsigned Tdim>
+bool mpm::Cell<Tdim>::initialise_element_thermal_matrix() {
+  bool status = true;
+  if (this->status()) {
+    try {
+      // Initialse heat capacity matrix N x N
+      heat_capacity_matrix_.resize(nnodes_, nnodes_);
+      heat_capacity_matrix_.setZero();      
+      // Initialse thermal conductivity matrix N x N
+      thermal_conductivity_matrix_.resize(nnodes_, nnodes_);
+      thermal_conductivity_matrix_.setZero();
+      // Initialse thermal expansivity matrix (N * Tdim) x N
+      thermal_expansivity_matrix_.resize(nnodes_ * Tdim, nnodes_);
+      thermal_expansivity_matrix_.setZero();
+    } catch (std::exception& exception) {
+      console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what());
+      status = false;
+    }
+  }
+  return status;
+}
+
+//! Compute local heat capacity stiffness matrix
+template <unsigned Tdim>
+void mpm::Cell<Tdim>::compute_local_heat_capacity_matrix(
+    const Eigen::VectorXd& shapefn, const double pvolume, 
+    const double multiplier) noexcept {
+
+  std::lock_guard<std::mutex> guard(cell_mutex_);
+  // Heat capacity matrix N x N
+  // Consistent matrix
+  // thermal_conductivity_matrix_ +=
+  //     shapefn * (shapefn.transpose()) * pvolume * multiplier;
+  // Lumped matrix
+  for (unsigned i = 0; i < this->nnodes_; ++i) 
+    thermal_conductivity_matrix_(i, i) += shapefn(i) * multiplier * pvolume;
+}
+
+//! Compute local thermal conductivity matrix
+template <unsigned Tdim>
+void mpm::Cell<Tdim>::compute_local_thermal_conductivity_matrix(
+    const Eigen::MatrixXd& grad_shapefn, double pvolume,
+    double multiplier) noexcept {
+
+  std::lock_guard<std::mutex> guard(cell_mutex_);
+  // Thermal conductivity matrix N x N
+  thermal_conductivity_matrix_ +=
+    grad_shapefn * (grad_shapefn.transpose()) * pvolume * multiplier;;
+}
+
+//! Compute local thermal expansivity matrix
+template <unsigned Tdim>
+void mpm::Cell<Tdim>::compute_local_thermal_expansivity_matrix(
+    const Eigen::VectorXd& shapefn, const Eigen::MatrixXd& bmatrix, 
+    const Eigen::MatrixXd& dmatrix,
+    const Eigen::VectorXd& identity_vector,
+    double pvolume, double multiplier) noexcept {
+
+  std::lock_guard<std::mutex> guard(cell_mutex_);
+  // Thermal expansivity matrix (N * Tdim) * N
+  // Tdim * (Tdim + 1) / 2 = M, M=1, 3, 6 for 1D, 2D, 3D
+  // bmatrix: M * (N * Tdim), bmatrix.transpose(): (N * Tdim) * M
+  // dmatrix: M * M
+  // identity_vector: M * 1
+  // shapefn: N * 1, shapefn.transpose(): 1 * N
+  thermal_expansivity_matrix_.noalias() +=
+      -bmatrix.transpose() * dmatrix * identity_vector * shapefn.transpose() * 
+                                                   multiplier * pvolume;
+}
+

include/io/io_mesh.h

@@ -146,6 +146,16 @@ class IOMesh {
           const std::string& displacement_constraints_file) = 0;
   /**@}*/
 
+  /**
+   * \defgroup Thermal Functions for Thermo-mechanical MPM
+   */
+  /**@{*/
+  //! Read temperature constraints file
+  virtual std::vector<std::tuple<mpm::Index, double>>
+      read_temperature_constraints(
+          const std::string& temperature_constraints_file) = 0;
+  /**@}*/  
+
 };  // IOMesh class
 }  // namespace mpm
 

include/io/io_mesh_ascii.h

@@ -133,6 +133,16 @@ class IOMeshAscii : public IOMesh<Tdim> {
           const std::string& displacement_constraints_file) override;
   /**@}*/
 
+  /**
+   * \defgroup Thermal Functions for Thermo-mechanical MPM
+   */
+  /**@{*/
+  //! Read temperature constraints file
+  std::vector<std::tuple<mpm::Index, double>>
+      read_temperature_constraints(
+          const std::string& temperature_constraints_file) override;
+  /**@}*/
+
  private:
   //! Logger
   std::shared_ptr<spdlog::logger> console_;

include/io/io_mesh_ascii.tcc

@@ -815,4 +815,47 @@ std::array<std::vector<double>, 2> mpm::IOMeshAscii<Tdim>::read_math_functions(
   }
 
   return xfx_values;
+}
+
+//! Read temperature constraints file
+template <unsigned Tdim>
+std::vector<std::tuple<mpm::Index, double>>
+    mpm::IOMeshAscii<Tdim>::read_temperature_constraints(
+        const std::string& temperature_constraints_file) {
+  // Particle temperature constraints
+  std::vector<std::tuple<mpm::Index, double>> constraints;
+  constraints.clear();
+
+  // input file stream
+  std::fstream file;
+  file.open(temperature_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;
+            // temperature
+            double temperature;
+            // Read stream
+            istream >> id >> temperature;
+            constraints.emplace_back(std::make_tuple(id, temperature));
+          }
+        }
+      }
+    }
+    file.close();
+  } catch (std::exception& exception) {
+    console_->error("Read temperature constraints: {}", exception.what());
+    file.close();
+  }
+  return constraints;
 }

include/io/logger.h

@@ -61,6 +61,15 @@ struct Logger {
   // Create a logger for MPM Semi-implicit Two Phase
   static const std::shared_ptr<spdlog::logger>
       mpm_semi_implicit_two_phase_logger;
+
+  // Create a logger for Thermo-Mechanical MPM explicit
+  static const std::shared_ptr<spdlog::logger>
+      mpm_explicit_thermal_logger;
+
+  // Create a logger for Thermo-Mechanical MPM mplicit
+  static const std::shared_ptr<spdlog::logger>
+      mpm_implicit_thermal_logger;      
+
 };
 
 }  // namespace mpm

include/linear_solvers/assemblers/assembler_base.h

@@ -173,6 +173,116 @@ class AssemblerBase {
 
   /**@{*/
 
+  /**
+   * \defgroup Thermal Functions forthermo-mechancial coupling MPM
+   */
+  /**@{*/
+  //! Assemble thermal expansivity matrix
+  virtual bool assemble_thermal_expansivity_matrix() {
+    throw std::runtime_error(
+        "Calling the base class function (assemble_thermal_expansivity_matrix) "
+        "in AssemblerBase:: illegal operation!");
+    return 0;   
+  };
+
+  //! Assemble thermal conductivity matrix
+  virtual bool assemble_thermal_conductivity_matrix() {
+    throw std::runtime_error(
+        "Calling the base class function (aassemble_thermal_conductivity_"
+        "matrix) in AssemblerBase:: illegal operation!");
+    return 0;   
+  };
+
+  // Assemble residual heat right vector
+  virtual bool assemble_residual_heat_right() {
+    throw std::runtime_error(
+        "Calling the base class function (assemble_residual_heat_right) in "
+        "AssemblerBase:: illegal operation!");
+    return 0;   
+  };
+
+  //! Assemble global stiffness matrix
+  virtual bool assemble_global_stiffness_matrix() {
+    throw std::runtime_error(
+        "Calling the base class function (assemble_global_stiffness_matrix) in "
+        "AssemblerBase:: illegal operation!");
+    return 0;   
+  };
+
+  //! Return global stiffness matrix
+  virtual Eigen::SparseMatrix<double>& global_stiffness_matrix() {
+    throw std::runtime_error(
+        "Calling the base class function (global_stiffness_matrix) in "
+        "AssemblerBase:: illegal operation!");
+  };
+
+  //! Assemble global residual matrix
+  virtual bool assemble_global_residual_right() {
+    throw std::runtime_error(
+        "Calling the base class function (assemble_global_residual_right) in "
+        "AssemblerBase:: illegal operation!");   
+  };
+
+  //! Return global residual rhs vector
+  virtual Eigen::VectorXd& global_residual_rhs_vector() {
+    throw std::runtime_error(
+        "Calling the base class function (global_residual_rhs_vector) in "
+        "AssemblerBase:: illegal operation!");
+  };
+
+  //! Assign temperature constraints
+  virtual bool assign_temperature_constraints(double current_time) {
+    throw std::runtime_error(
+        "Calling the base class function (assign_temperature_constraints) in "
+        "AssemblerBase:: illegal operation!");
+    return 0;   
+  };
+
+  //! Apply temperature increment constraints vector
+  virtual void apply_temperature_increment_constraints() {
+    throw std::runtime_error(
+        "Calling the base class function (apply_temperature_constraints) in "
+        "AssemblerBase:: illegal operation!");
+  };
+
+  //! Apply temperature constraints vector
+  virtual void apply_coupling_constraints() {
+    throw std::runtime_error(
+        "Calling the base class function (apply_coupling_constraints) in "
+        "AssemblerBase:: illegal operation!");
+  };
+
+  //! Assign solution increment
+  virtual void assign_solution_increment(
+      const Eigen::VectorXd& solution_increment) {
+    throw std::runtime_error(
+        "Calling the base class function (assign_solution_increment) in "
+        "AssemblerBase:: illegal operation!");
+  };
+
+  //! Return solution increment
+  virtual Eigen::VectorXd& solution_increment() {
+    throw std::runtime_error(
+        "Calling the base class function (solution_increment) in "
+        "AssemblerBase:: illegal operation!");
+  };
+
+  //! Assign temperature increment
+  virtual void assign_temperature_increment(
+      const Eigen::VectorXd& temperature_increment) {
+    throw std::runtime_error(
+        "Calling the base class function (assign_temperature_increment) in "
+        "AssemblerBase:: illegal operation!");
+  };
+
+  //! Return temperature increment
+  virtual Eigen::VectorXd& temperature_increment() {
+    throw std::runtime_error(
+        "Calling the base class function (temperature_increment) in "
+        "AssemblerBase:: illegal operation!");
+  };  
+  /**@{*/
+
   //! Navier-Stokes
   //! functions-------------------------------------------------------- Return
   //! laplacian matrix