DiffDeviator.hpp

Go to the documentation of this file.

/**
 * @file DiffDeviator.hpp
 *
 */
 
#pragma once
 
namespace FTensor {
 
/**
 * @brief Fourth-order differential deviator tensor
 *
 * @tparam T
 * @tparam Dim
 */
template <class T = int, int Dim = 3> class DiffDeviator {
public:
  using value_type = typename promote<T, double>::V;
 
  constexpr value_type operator()(const int N0, const int N1, const int N2,
                                  const int N3) const {
    const auto diff = DiffTensor<T>()(N0, N1, N2, N3);
    const auto delta_ij = (N0 == N1) ? value_type(1) : value_type(0);
    const auto delta_kl = (N2 == N3) ? value_type(1) : value_type(0);
    return diff - (delta_ij * delta_kl) / value_type(3);
  }
 
  template <char i, char j, char k, char l, int Dim01, int Dim23>
  typename std::enable_if<
      (Dim01 == 3 && Dim23 == Dim),
      Ddg_Expr<DiffDeviator<T, Dim>, value_type, Dim01, Dim23, i, j, k,
               l>>::type
  operator()(const Index<i, Dim01> &, const Index<j, Dim01> &,
             const Index<k, Dim23> &, const Index<l, Dim23> &) const {
    return Ddg_Expr<DiffDeviator<T, Dim>, value_type, Dim01, Dim23, i, j, k,
                    l>(*this);
  };
};
 
template <class T = int, int Dim> constexpr auto diff_deviator() {
  return DiffDeviator<T, Dim>();
}
 
template <class T = int, int Dim>
constexpr auto diff_deviator(const Number<Dim> &) {
  return DiffDeviator<T, Dim>();
}
 
template <class T = int, int Dim, char i, char j, char k, char l, int Dim01,
          int Dim23>
typename std::enable_if<
    (Dim01 == 3 && Dim23 == Dim),
    Ddg_Expr<DiffDeviator<T, Dim>,
             typename DiffDeviator<T, Dim>::value_type, Dim01, Dim23, i, j, k,
             l>>::type
diff_deviator(const Index<i, Dim01> &, const Index<j, Dim01> &,
              const Index<k, Dim23> &, const Index<l, Dim23> &) {
  return DiffDeviator<T, Dim>()(Index<i, Dim01>(), Index<j, Dim01>(),
                                Index<k, Dim23>(), Index<l, Dim23>());
}
 
} // namespace FTensor