DiffTensor.hpp

Go to the documentation of this file.

/**
 * @file DiffTensor.hpp
 *
 */
 
#pragma once
 
namespace FTensor {
 
/**
 * @brief Fourth-order differential tensor symmetric in both index pairs
 *
 * @tparam T
 */
template <class T = int> class DiffTensor {
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 first =
        (N0 == N2 && N1 == N3) ? value_type(1) : value_type(0);
    const auto second =
        (N0 == N3 && N1 == N2) ? value_type(1) : value_type(0);
    return (first + second) / value_type(2);
  }
 
  template <char i, char j, char k, char l, int Dim01, int Dim23>
  Ddg_Expr<DiffTensor<T>, value_type, Dim01, Dim23, i, j, k, l>
  operator()(const Index<i, Dim01> &, const Index<j, Dim01> &,
             const Index<k, Dim23> &, const Index<l, Dim23> &) const {
    return Ddg_Expr<DiffTensor<T>, value_type, Dim01, Dim23, i, j, k, l>(
        *this);
  };
};
 
template <class T = int> constexpr auto diff_tensor() {
  return DiffTensor<T>();
}
 
template <class T = int, char i, char j, char k, char l, int Dim01, int Dim23>
Ddg_Expr<DiffTensor<T>, typename DiffTensor<T>::value_type, Dim01, Dim23, i, j,
         k, l>
diff_tensor(const Index<i, Dim01> &, const Index<j, Dim01> &,
            const Index<k, Dim23> &, const Index<l, Dim23> &) {
  return DiffTensor<T>()(Index<i, Dim01>(), Index<j, Dim01>(),
                         Index<k, Dim23>(), Index<l, Dim23>());
}
 
} // namespace FTensor