EshelbianAux.hpp

Go to the documentation of this file.

/**
 * @file EshelbianAux.hpp
 * @brief Auxilary functions for Eshelbian plasticity
 * @date 2024-08-30
 *
 * @copyright Copyright (c) 2024
 *
 */
 
#include <Lie.hpp>
 
namespace EshelbianPlasticity {
 
inline auto diff_deviator(FTensor::Ddg<double, 3, 3> &&t_diff_stress) {
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
 
  FTensor::Ddg<double, 3, 3> t_diff_deviator;
  t_diff_deviator(i, j, k, l) = t_diff_stress(i, j, k, l);
 
  constexpr double third = boost::math::constants::third<double>();
 
  t_diff_deviator(0, 0, 0, 0) -= third;
  t_diff_deviator(0, 0, 1, 1) -= third;
 
  t_diff_deviator(1, 1, 0, 0) -= third;
  t_diff_deviator(1, 1, 1, 1) -= third;
 
  t_diff_deviator(2, 2, 0, 0) -= third;
  t_diff_deviator(2, 2, 1, 1) -= third;
 
  t_diff_deviator(0, 0, 2, 2) -= third;
  t_diff_deviator(1, 1, 2, 2) -= third;
  t_diff_deviator(2, 2, 2, 2) -= third;
 
  return t_diff_deviator;
}
 
constexpr static auto size_symm = (3 * (3 + 1)) / 2;
 
inline auto diff_tensor() {
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
  FTensor::Ddg<double, 3, 3> t_diff;
  constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
  t_diff(i, j, k, l) = (t_kd(i, k) ^ t_kd(j, l)) / 4.;
  return t_diff;
};
 
inline auto symm_L_tensor() {
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'L', size_symm> L;
  FTensor::Dg<double, 3, size_symm> t_L;
  t_L(i, j, L) = 0;
  t_L(0, 0, 0) = 1;
  t_L(1, 1, 3) = 1;
  t_L(2, 2, 5) = 1;
  t_L(1, 0, 1) = 1;
  t_L(2, 0, 2) = 1;
  t_L(2, 1, 4) = 1;
  return t_L;
}
 
inline auto diff_symmetrize() {
 
  FTensor::Index<'i', SPACE_DIM> i;
  FTensor::Index<'j', SPACE_DIM> j;
  FTensor::Index<'k', SPACE_DIM> k;
  FTensor::Index<'l', SPACE_DIM> l;
 
  FTensor::Ddg<double, SPACE_DIM, SPACE_DIM> t_diff;
 
  t_diff(i, j, k, l) = 0;
  t_diff(0, 0, 0, 0) = 1;
  t_diff(1, 1, 1, 1) = 1;
 
  t_diff(1, 0, 1, 0) = 0.5;
  t_diff(1, 0, 0, 1) = 0.5;
 
  t_diff(0, 1, 0, 1) = 0.5;
  t_diff(0, 1, 1, 0) = 0.5;
 
  if constexpr (SPACE_DIM == 3) {
    t_diff(2, 2, 2, 2) = 1;
 
    t_diff(2, 0, 2, 0) = 0.5;
    t_diff(2, 0, 0, 2) = 0.5;
    t_diff(0, 2, 0, 2) = 0.5;
    t_diff(0, 2, 2, 0) = 0.5;
 
    t_diff(2, 1, 2, 1) = 0.5;
    t_diff(2, 1, 1, 2) = 0.5;
    t_diff(1, 2, 1, 2) = 0.5;
    t_diff(1, 2, 2, 1) = 0.5;
  }
 
  return t_diff;
}
 
template <class T> struct TensorTypeExtractor {
  typedef typename std::remove_pointer<T>::type Type;
};
template <class T, int I> struct TensorTypeExtractor<FTensor::PackPtr<T, I>> {
  typedef typename std::remove_pointer<T>::type Type;
};
 
struct isEq {
  static inline auto check(const double &a, const double &b) {
    constexpr double eps = std::numeric_limits<double>::epsilon();
    return std::abs(a - b) / absMax < eps;
  }
  static inline double absMax = 1.;
};
 
inline auto is_eq(const double &a, const double &b) {
  return isEq::check(a, b);
};
 
template <int DIM> inline auto get_uniq_nb(double *ptr) {
  std::array<double, DIM> tmp;
  std::copy(ptr, &ptr[DIM], tmp.begin());
  std::sort(tmp.begin(), tmp.end());
  isEq::absMax = std::max(std::abs(tmp[0]), std::abs(tmp[DIM - 1]));
  constexpr double eps = std::numeric_limits<double>::epsilon();
  isEq::absMax = std::max(isEq::absMax, static_cast<double>(eps));
  return std::distance(tmp.begin(), std::unique(tmp.begin(), tmp.end(), is_eq));
}
 
template <typename A, typename B>
inline auto sort_eigen_vals(A &eig, B &eigen_vec) {
 
  constexpr int dim = 3;
 
  isEq::absMax =
      std::max(std::max(std::abs(eig(0)), std::abs(eig(1))), std::abs(eig(2)));
  constexpr double eps = std::numeric_limits<double>::epsilon();
  isEq::absMax = std::max(isEq::absMax, static_cast<double>(eps));
 
  int i = 0, j = 1, k = 2;
 
  if (is_eq(eig(0), eig(1))) {
    i = 0;
    j = 2;
    k = 1;
  } else if (is_eq(eig(0), eig(2))) {
    i = 0;
    j = 1;
    k = 2;
  } else if (is_eq(eig(1), eig(2))) {
    i = 1;
    j = 0;
    k = 2;
  }
 
  FTensor::Tensor2<double, 3, 3> eigen_vec_c{
      eigen_vec(i, 0), eigen_vec(i, 1), eigen_vec(i, 2),
 
      eigen_vec(j, 0), eigen_vec(j, 1), eigen_vec(j, 2),
 
      eigen_vec(k, 0), eigen_vec(k, 1), eigen_vec(k, 2)};
 
  FTensor::Tensor1<double, 3> eig_c{eig(i), eig(j), eig(k)};
 
  {
    FTENSOR_INDEX(dim, i);
    FTENSOR_INDEX(dim, j);
    eig(i) = eig_c(i);
    eigen_vec(i, j) = eigen_vec_c(i, j);
  }
}
}