RankFourTensor.h

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://www.mooseframework.org
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#pragma once

#include "Moose.h"
#include "MooseTypes.h"
#include "ADRankTwoTensorForward.h"
#include "ADRankFourTensorForward.h"
#include "ADRankThreeTensorForward.h"

#include "libmesh/libmesh.h"
#include "libmesh/tuple_of.h"
#include "libmesh/int_range.h"

#include "metaphysicl/raw_type.h"

#include <petscsys.h>

#include <Eigen/Core>
#include <Eigen/Dense>

using libMesh::Real;
using libMesh::tuple_of;
namespace libMesh
{
template <typename>
class TensorValue;
template <typename>
class TypeTensor;
template <typename>
class VectorValue;
}

// Forward declarations
class MooseEnum;

namespace MathUtils
{
template <typename T>
void mooseSetToZero(T & v);

template <>
void mooseSetToZero<RankFourTensor>(RankFourTensor & v);
template <>
void mooseSetToZero<ADRankFourTensor>(ADRankFourTensor & v);
}

template <typename T>
class RankFourTensorTempl
{
public:
  static constexpr unsigned int N = Moose::dim;
  static constexpr unsigned int N2 = N * N;
  static constexpr unsigned int N3 = N * N * N;
  static constexpr unsigned int N4 = N * N * N * N;

  typedef tuple_of<4, unsigned int> index_type;
  typedef T value_type;

  enum InitMethod
  {
    initNone,
    initIdentity,
    initIdentityFour,
    initIdentitySymmetricFour,
    initIdentityDeviatoric
  };

  enum FillMethod
  {
    antisymmetric,
    symmetric9,
    symmetric21,
    general_isotropic,
    symmetric_isotropic,
    symmetric_isotropic_E_nu,
    antisymmetric_isotropic,
    axisymmetric_rz,
    general,
    principal,
    orthotropic
  };

  template <template <typename> class Tensor, typename Scalar>
  struct TwoTensorMultTraits
  {
    static const bool value = false;
  };
  template <typename Scalar>
  struct TwoTensorMultTraits<RankTwoTensorTempl, Scalar>
  {
    static const bool value = ScalarTraits<Scalar>::value;
  };
  template <typename Scalar>
  struct TwoTensorMultTraits<TensorValue, Scalar>
  {
    static const bool value = ScalarTraits<Scalar>::value;
  };
  template <typename Scalar>
  struct TwoTensorMultTraits<TypeTensor, Scalar>
  {
    static const bool value = ScalarTraits<Scalar>::value;
  };

  RankFourTensorTempl();

  RankFourTensorTempl(const InitMethod);

  RankFourTensorTempl(const std::vector<T> &, FillMethod);

  RankFourTensorTempl(const RankFourTensorTempl<T> & a) = default;

  template <typename T2>
  RankFourTensorTempl(const RankFourTensorTempl<T2> & copy);

  // Named constructors
  static RankFourTensorTempl<T> Identity() { return RankFourTensorTempl<T>(initIdentity); }
  static RankFourTensorTempl<T> IdentityFour() { return RankFourTensorTempl<T>(initIdentityFour); };
  static RankFourTensorTempl<T> IdentityDeviatoric()
  {
    return RankFourTensorTempl<T>(initIdentityDeviatoric);
  };

  inline T & operator()(unsigned int i, unsigned int j, unsigned int k, unsigned int l)
  {
    return _vals[i * N3 + j * N2 + k * N + l];
  }

  inline const T & operator()(unsigned int i, unsigned int j, unsigned int k, unsigned int l) const
  {
    return _vals[i * N3 + j * N2 + k * N + l];
  }

  void zero();

  void print(std::ostream & stm = Moose::out) const;

  void printReal(std::ostream & stm = Moose::out) const;

  RankFourTensorTempl<T> & operator=(const RankFourTensorTempl<T> & a);

  template <typename Scalar>
  typename boostcopy::enable_if_c<ScalarTraits<Scalar>::value, RankFourTensorTempl &>::type
  operator=(const Scalar & libmesh_dbg_var(p))
  {
    libmesh_assert_equal_to(p, Scalar(0));
    this->zero();
    return *this;
  }

  template <template <typename> class Tensor, typename T2>
  auto operator*(const Tensor<T2> & a) const ->
      typename std::enable_if<TwoTensorMultTraits<Tensor, T2>::value,
                              RankTwoTensorTempl<decltype(T() * T2())>>::type;

  template <typename T2>
  auto operator*(const T2 & a) const ->
      typename std::enable_if<ScalarTraits<T2>::value,
                              RankFourTensorTempl<decltype(T() * T2())>>::type;

  RankFourTensorTempl<T> & operator*=(const T & a);

  template <typename T2>
  auto operator/(const T2 & a) const ->
      typename std::enable_if<ScalarTraits<T2>::value,
                              RankFourTensorTempl<decltype(T() / T2())>>::type;

  RankFourTensorTempl<T> & operator/=(const T & a);

  RankFourTensorTempl<T> & operator+=(const RankFourTensorTempl<T> & a);

  template <typename T2>
  auto operator+(const RankFourTensorTempl<T2> & a) const
      -> RankFourTensorTempl<decltype(T() + T2())>;

  RankFourTensorTempl<T> & operator-=(const RankFourTensorTempl<T> & a);

  template <typename T2>
  auto operator-(const RankFourTensorTempl<T2> & a) const
      -> RankFourTensorTempl<decltype(T() - T2())>;

  RankFourTensorTempl<T> operator-() const;

  template <typename T2>
  auto operator*(const RankFourTensorTempl<T2> & a) const
      -> RankFourTensorTempl<decltype(T() * T2())>;

  T L2norm() const;

  RankFourTensorTempl<T> invSymm() const;

  RankFourTensorTempl<T> inverse() const;

  void rotate(const TypeTensor<T> & R);

  RankFourTensorTempl<T> transposeMajor() const;

  RankFourTensorTempl<T> transposeIj() const;

  RankFourTensorTempl<T> transposeKl() const;

  template <int m>
  RankThreeTensorTempl<T> contraction(const VectorValue<T> & b) const;

  void surfaceFillFromInputVector(const std::vector<T> & input);

  static MooseEnum fillMethodEnum();

  void fillFromInputVector(const std::vector<T> & input, FillMethod fill_method);

  void fillGeneralIsotropic(const T & i0, const T & i1, const T & i2);
  void fillAntisymmetricIsotropic(const T & i0);
  void fillSymmetricIsotropic(const T & i0, const T & i1);
  void fillSymmetricIsotropicEandNu(const T & E, const T & nu);

  template <typename T2>
  void fillSymmetric9FromInputVector(const T2 & input);

  template <typename T2>
  void fillSymmetric21FromInputVector(const T2 & input);

  RankTwoTensorTempl<T> innerProductTranspose(const RankTwoTensorTempl<T> &) const;

  T contractionIj(unsigned int, unsigned int, const RankTwoTensorTempl<T> &) const;

  T contractionKl(unsigned int, unsigned int, const RankTwoTensorTempl<T> &) const;

  T sum3x3() const;

  VectorValue<T> sum3x1() const;

  RankFourTensorTempl<T> tripleProductJkl(const RankTwoTensorTempl<T> &,
                                          const RankTwoTensorTempl<T> &,
                                          const RankTwoTensorTempl<T> &) const;
  RankFourTensorTempl<T> tripleProductIkl(const RankTwoTensorTempl<T> &,
                                          const RankTwoTensorTempl<T> &,
                                          const RankTwoTensorTempl<T> &) const;
  RankFourTensorTempl<T> tripleProductIjl(const RankTwoTensorTempl<T> &,
                                          const RankTwoTensorTempl<T> &,
                                          const RankTwoTensorTempl<T> &) const;
  RankFourTensorTempl<T> tripleProductIjk(const RankTwoTensorTempl<T> &,
                                          const RankTwoTensorTempl<T> &,
                                          const RankTwoTensorTempl<T> &) const;

  RankFourTensorTempl<T> singleProductI(const RankTwoTensorTempl<T> &) const;
  RankFourTensorTempl<T> singleProductJ(const RankTwoTensorTempl<T> &) const;
  RankFourTensorTempl<T> singleProductK(const RankTwoTensorTempl<T> &) const;
  RankFourTensorTempl<T> singleProductL(const RankTwoTensorTempl<T> &) const;

  bool isSymmetric() const;

  bool isIsotropic() const;

protected:
  T _vals[N4];

  void fillAntisymmetricFromInputVector(const std::vector<T> & input);

  void fillGeneralIsotropicFromInputVector(const std::vector<T> & input);

  void fillAntisymmetricIsotropicFromInputVector(const std::vector<T> & input);

  void fillSymmetricIsotropicFromInputVector(const std::vector<T> & input);

  void fillSymmetricIsotropicEandNuFromInputVector(const std::vector<T> & input);

  void fillAxisymmetricRZFromInputVector(const std::vector<T> & input);

  void fillGeneralFromInputVector(const std::vector<T> & input);

  void fillPrincipalFromInputVector(const std::vector<T> & input);

  void fillGeneralOrthotropicFromInputVector(const std::vector<T> & input);

  template <class T2>
  friend void dataStore(std::ostream &, RankFourTensorTempl<T2> &, void *);

  template <class T2>
  friend void dataLoad(std::istream &, RankFourTensorTempl<T2> &, void *);

  template <typename T2>
  friend class RankTwoTensorTempl;
  template <typename T2>
  friend class RankFourTensorTempl;
  template <typename T2>
  friend class RankThreeTensorTempl;
};

namespace MetaPhysicL
{
template <typename T>
struct RawType<RankFourTensorTempl<T>>
{
  typedef RankFourTensorTempl<typename RawType<T>::value_type> value_type;

  static value_type value(const RankFourTensorTempl<T> & in)
  {
    constexpr auto N = RankFourTensorTempl<T>::N;
    value_type ret;
    for (auto i : make_range(N))
      for (auto j : make_range(N))
        for (auto k : make_range(N))
          for (auto l : make_range(N))
            ret(i, j, k, l) = raw_value(in(i, j, k, l));

    return ret;
  }
};
}

template <typename T1, typename T2>
inline auto
operator*(const T1 & a, const RankFourTensorTempl<T2> & b) ->
    typename std::enable_if<ScalarTraits<T1>::value,
                            RankFourTensorTempl<decltype(T1() * T2())>>::type
{
  return b * a;
}

template <typename T>
template <typename T2>
RankFourTensorTempl<T>::RankFourTensorTempl(const RankFourTensorTempl<T2> & copy)
{
  for (auto i : make_range(N4))
    _vals[i] = copy._vals[i];
}

template <typename T>
template <typename T2>
auto
RankFourTensorTempl<T>::operator*(const T2 & b) const ->
    typename std::enable_if<ScalarTraits<T2>::value,
                            RankFourTensorTempl<decltype(T() * T2())>>::type
{
  typedef decltype(T() * T2()) ValueType;
  RankFourTensorTempl<ValueType> result;

  for (auto i : make_range(N4))
    result._vals[i] = _vals[i] * b;

  return result;
}

template <typename T>
template <typename T2>
auto
RankFourTensorTempl<T>::operator/(const T2 & b) const ->
    typename std::enable_if<ScalarTraits<T2>::value,
                            RankFourTensorTempl<decltype(T() / T2())>>::type
{
  RankFourTensorTempl<decltype(T() / T2())> result;
  for (auto i : make_range(N4))
    result._vals[i] = _vals[i] / b;
  return result;
}

template <typename T>
template <typename T2>
void
RankFourTensorTempl<T>::fillSymmetric9FromInputVector(const T2 & input)
{
  mooseAssert(input.size() == 9,
              "To use fillSymmetric9FromInputVector, your input must have size 9.");
  zero();

  (*this)(0, 0, 0, 0) = input[0]; // C1111
  (*this)(1, 1, 1, 1) = input[3]; // C2222
  (*this)(2, 2, 2, 2) = input[5]; // C3333

  (*this)(0, 0, 1, 1) = input[1]; // C1122
  (*this)(1, 1, 0, 0) = input[1];

  (*this)(0, 0, 2, 2) = input[2]; // C1133
  (*this)(2, 2, 0, 0) = input[2];

  (*this)(1, 1, 2, 2) = input[4]; // C2233
  (*this)(2, 2, 1, 1) = input[4];

  (*this)(1, 2, 1, 2) = input[6]; // C2323
  (*this)(2, 1, 2, 1) = input[6];
  (*this)(2, 1, 1, 2) = input[6];
  (*this)(1, 2, 2, 1) = input[6];

  (*this)(0, 2, 0, 2) = input[7]; // C1313
  (*this)(2, 0, 2, 0) = input[7];
  (*this)(2, 0, 0, 2) = input[7];
  (*this)(0, 2, 2, 0) = input[7];

  (*this)(0, 1, 0, 1) = input[8]; // C1212
  (*this)(1, 0, 1, 0) = input[8];
  (*this)(1, 0, 0, 1) = input[8];
  (*this)(0, 1, 1, 0) = input[8];
}
template <typename T>
template <typename T2>
void
RankFourTensorTempl<T>::fillSymmetric21FromInputVector(const T2 & input)
{
  mooseAssert(input.size() == 21,
              "To use fillSymmetric21FromInputVector, your input must have size 21.");

  (*this)(0, 0, 0, 0) = input[0];  // C1111
  (*this)(1, 1, 1, 1) = input[6];  // C2222
  (*this)(2, 2, 2, 2) = input[11]; // C3333

  (*this)(0, 0, 1, 1) = input[1]; // C1122
  (*this)(1, 1, 0, 0) = input[1];

  (*this)(0, 0, 2, 2) = input[2]; // C1133
  (*this)(2, 2, 0, 0) = input[2];

  (*this)(1, 1, 2, 2) = input[7]; // C2233
  (*this)(2, 2, 1, 1) = input[7];

  (*this)(0, 0, 0, 2) = input[4]; // C1113
  (*this)(0, 0, 2, 0) = input[4];
  (*this)(0, 2, 0, 0) = input[4];
  (*this)(2, 0, 0, 0) = input[4];

  (*this)(0, 0, 0, 1) = input[5]; // C1112
  (*this)(0, 0, 1, 0) = input[5];
  (*this)(0, 1, 0, 0) = input[5];
  (*this)(1, 0, 0, 0) = input[5];

  (*this)(1, 1, 1, 2) = input[8]; // C2223
  (*this)(1, 1, 2, 1) = input[8];
  (*this)(1, 2, 1, 1) = input[8];
  (*this)(2, 1, 1, 1) = input[8];

  (*this)(1, 1, 1, 0) = input[10];
  (*this)(1, 1, 0, 1) = input[10];
  (*this)(1, 0, 1, 1) = input[10];
  (*this)(0, 1, 1, 1) = input[10]; // C2212 //flipped for filling purposes

  (*this)(2, 2, 2, 1) = input[12];
  (*this)(2, 2, 1, 2) = input[12];
  (*this)(2, 1, 2, 2) = input[12];
  (*this)(1, 2, 2, 2) = input[12]; // C3323 //flipped for filling purposes

  (*this)(2, 2, 2, 0) = input[13];
  (*this)(2, 2, 0, 2) = input[13];
  (*this)(2, 0, 2, 2) = input[13];
  (*this)(0, 2, 2, 2) = input[13]; // C3313 //flipped for filling purposes

  (*this)(0, 0, 1, 2) = input[3]; // C1123
  (*this)(0, 0, 2, 1) = input[3];
  (*this)(1, 2, 0, 0) = input[3];
  (*this)(2, 1, 0, 0) = input[3];

  (*this)(1, 1, 0, 2) = input[9];
  (*this)(1, 1, 2, 0) = input[9];
  (*this)(0, 2, 1, 1) = input[9]; // C2213  //flipped for filling purposes
  (*this)(2, 0, 1, 1) = input[9];

  (*this)(2, 2, 0, 1) = input[14];
  (*this)(2, 2, 1, 0) = input[14];
  (*this)(0, 1, 2, 2) = input[14]; // C3312 //flipped for filling purposes
  (*this)(1, 0, 2, 2) = input[14];

  (*this)(1, 2, 1, 2) = input[15]; // C2323
  (*this)(2, 1, 2, 1) = input[15];
  (*this)(2, 1, 1, 2) = input[15];
  (*this)(1, 2, 2, 1) = input[15];

  (*this)(0, 2, 0, 2) = input[18]; // C1313
  (*this)(2, 0, 2, 0) = input[18];
  (*this)(2, 0, 0, 2) = input[18];
  (*this)(0, 2, 2, 0) = input[18];

  (*this)(0, 1, 0, 1) = input[20]; // C1212
  (*this)(1, 0, 1, 0) = input[20];
  (*this)(1, 0, 0, 1) = input[20];
  (*this)(0, 1, 1, 0) = input[20];

  (*this)(1, 2, 0, 2) = input[16];
  (*this)(0, 2, 1, 2) = input[16]; // C2313 //flipped for filling purposes
  (*this)(2, 1, 0, 2) = input[16];
  (*this)(1, 2, 2, 0) = input[16];
  (*this)(2, 0, 1, 2) = input[16];
  (*this)(0, 2, 2, 1) = input[16];
  (*this)(2, 1, 2, 0) = input[16];
  (*this)(2, 0, 2, 1) = input[16];

  (*this)(1, 2, 0, 1) = input[17];
  (*this)(0, 1, 1, 2) = input[17]; // C2312 //flipped for filling purposes
  (*this)(2, 1, 0, 1) = input[17];
  (*this)(1, 2, 1, 0) = input[17];
  (*this)(1, 0, 1, 2) = input[17];
  (*this)(0, 1, 2, 1) = input[17];
  (*this)(2, 1, 1, 0) = input[17];
  (*this)(1, 0, 2, 1) = input[17];

  (*this)(0, 2, 0, 1) = input[19];
  (*this)(0, 1, 0, 2) = input[19]; // C1312 //flipped for filling purposes
  (*this)(2, 0, 0, 1) = input[19];
  (*this)(0, 2, 1, 0) = input[19];
  (*this)(1, 0, 0, 2) = input[19];
  (*this)(0, 1, 2, 0) = input[19];
  (*this)(2, 0, 1, 0) = input[19];
  (*this)(1, 0, 2, 0) = input[19];
}

template <typename T>
RankFourTensorTempl<T>
RankFourTensorTempl<T>::inverse() const
{

  // The inverse of a 3x3x3x3 in the C_ijkl*A_klmn = de_im de_jn sense is
  // simply the inverse of the 9x9 matrix of the tensor entries.
  // So all we need to do is inverse _vals (with the appropriate row-major
  // storage)

  RankFourTensorTempl<T> result;

  if constexpr (RankFourTensorTempl<T>::N4 * sizeof(T) > EIGEN_STACK_ALLOCATION_LIMIT)
  {
    // Allocate on the heap if you're going to exceed the stack size limit
    Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> mat(9, 9);
    for (auto i : make_range(9 * 9))
      mat(i) = _vals[i];

    mat = mat.inverse();

    for (auto i : make_range(9 * 9))
      result._vals[i] = mat(i);
  }
  else
  {
    // Allocate on the stack if small enough
    const Eigen::Map<const Eigen::Matrix<T, 9, 9, Eigen::RowMajor>> mat(&_vals[0]);
    Eigen::Map<Eigen::Matrix<T, 9, 9, Eigen::RowMajor>> res(&result._vals[0]);
    res = mat.inverse();
  }

  return result;
}

template <typename T>
template <int m>
RankThreeTensorTempl<T>
RankFourTensorTempl<T>::contraction(const VectorValue<T> & b) const
{
  RankThreeTensorTempl<T> result;
  static constexpr std::size_t z[4][3] = {{1, 2, 3}, {0, 2, 3}, {0, 1, 3}, {0, 1, 2}};
  std::size_t x[4];
  for (x[0] = 0; x[0] < N; ++x[0])
    for (x[1] = 0; x[1] < N; ++x[1])
      for (x[2] = 0; x[2] < N; ++x[2])
        for (x[3] = 0; x[3] < N; ++x[3])
          result(x[z[m][0]], x[z[m][1]], x[z[m][2]]) += (*this)(x[0], x[1], x[2], x[3]) * b(x[m]);

  return result;
}