SymmetricRankFourTensor.h

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* 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 "ADRankTwoTensorForward.h"
#include "ADRankFourTensorForward.h"
#include "ADSymmetricRankTwoTensorForward.h"
#include "ADSymmetricRankFourTensorForward.h"
#include "MooseUtils.h"
#include "MathUtils.h"

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

#include "metaphysicl/raw_type.h"

#include <petscsys.h>

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

#include <array>

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 <>
void mooseSetToZero<SymmetricRankFourTensor>(SymmetricRankFourTensor & v);
template <>
void mooseSetToZero<ADSymmetricRankFourTensor>(ADSymmetricRankFourTensor & v);
}

template <typename T>
class SymmetricRankFourTensorTempl
{
public:
  static constexpr unsigned int Ndim = LIBMESH_DIM;
  static constexpr unsigned int N = Ndim + Ndim * (Ndim - 1) / 2;
  static constexpr unsigned int N2 = N * N;

  // Full tensor indices in the Mandel/Voigt representation
  static constexpr unsigned int full_index[6][6][4] = {
      {{0, 0, 0, 0}, {0, 0, 1, 1}, {0, 0, 2, 2}, {0, 0, 1, 2}, {0, 0, 0, 2}, {0, 0, 0, 1}},
      {{1, 1, 0, 0}, {1, 1, 1, 1}, {1, 1, 2, 2}, {1, 1, 1, 2}, {1, 1, 0, 2}, {1, 1, 0, 1}},
      {{2, 2, 0, 0}, {2, 2, 1, 1}, {2, 2, 2, 2}, {2, 2, 1, 2}, {2, 2, 0, 2}, {2, 2, 0, 1}},
      {{1, 2, 0, 0}, {1, 2, 1, 1}, {1, 2, 2, 2}, {1, 2, 1, 2}, {1, 2, 0, 2}, {1, 2, 0, 1}},
      {{0, 2, 0, 0}, {0, 2, 1, 1}, {0, 2, 2, 2}, {0, 2, 1, 2}, {0, 2, 0, 2}, {0, 2, 0, 1}},
      {{0, 1, 0, 0}, {0, 1, 1, 1}, {0, 1, 2, 2}, {0, 1, 1, 2}, {0, 1, 0, 2}, {0, 1, 0, 1}}};

  static constexpr Real mandelFactor(unsigned int i, unsigned int j)
  {
    return SymmetricRankTwoTensorTempl<T>::mandelFactor(i) *
           SymmetricRankTwoTensorTempl<T>::mandelFactor(j);
  }

  enum InitMethod
  {
    initNone,
    initIdentity,
    initIdentitySymmetricFour
  };

  enum FillMethod
  {
    symmetric9,               // fillSymmetric9FromInputVector
    symmetric21,              // fillSymmetric21FromInputVector
    symmetric_isotropic,      // fillSymmetricIsotropicFromInputVector
    symmetric_isotropic_E_nu, // fillSymmetricIsotropicEandNu
    axisymmetric_rz,          // fillAxisymmetricRZFromInputVector
    principal,                // fillPrincipalFromInputVector
    orthotropic               // fillGeneralOrthotropicFromInputVector
  };

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

  SymmetricRankFourTensorTempl();

  SymmetricRankFourTensorTempl(const InitMethod);

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

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

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

  explicit SymmetricRankFourTensorTempl(const RankFourTensorTempl<T> & a);

  explicit operator RankFourTensorTempl<T>();

  // Named constructors
  static SymmetricRankFourTensorTempl<T> identity()
  {
    return SymmetricRankFourTensorTempl<T>(initIdentity);
  }
  static SymmetricRankFourTensorTempl<T> identitySymmetricFour()
  {
    return SymmetricRankFourTensorTempl<T>(initIdentitySymmetricFour);
  };

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

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

  void zero();

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

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

  friend std::ostream & operator<<(std::ostream & os, const SymmetricRankFourTensorTempl<T> & t)
  {
    t.print(os);
    return os;
  }

  SymmetricRankFourTensorTempl<T> & operator=(const SymmetricRankFourTensorTempl<T> & a) = default;

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

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

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

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

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

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

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

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

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

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

  SymmetricRankFourTensorTempl<T> operator-() const;

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

  T L2norm() const;

  SymmetricRankFourTensorTempl<T> invSymm() const;

  static SymmetricRankFourTensorTempl<T> rotationMatrix(const TypeTensor<T> & R);

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

  SymmetricRankFourTensorTempl<T> transposeMajor() const;

  SymmetricRankFourTensorTempl<T> transposeIj() const { return *this; }

  static MooseEnum fillMethodEnum();

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

  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);

  T sum3x3() const;

  libMesh::VectorValue<T> sum3x1() const;

  bool isSymmetric() const;

  bool isIsotropic() const;

protected:
  std::array<T, N2> _vals;

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

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

  void fillAxisymmetricRZFromInputVector(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 &, SymmetricRankFourTensorTempl<T2> &, void *);

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

  template <typename T2>
  friend class SymmetricRankTwoTensorTempl;
  template <typename T2>
  friend class SymmetricRankFourTensorTempl;
  template <typename T2>
  friend class RankThreeTensorTempl;
};

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

  static value_type value(const SymmetricRankFourTensorTempl<T> & in)
  {
    value_type ret;
    for (unsigned int i = 0; i < SymmetricRankFourTensorTempl<T>::N; ++i)
      for (unsigned int j = 0; j < SymmetricRankFourTensorTempl<T>::N; ++j)
        ret(i, j) = raw_value(in(i, j));

    return ret;
  }
};
}

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

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

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

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

  return result;
}

template <typename T>
template <typename T2>
auto
SymmetricRankFourTensorTempl<T>::operator*(const SymmetricRankTwoTensorTempl<T2> & b) const
    -> SymmetricRankTwoTensorTempl<decltype(T() * T2())>
{
  typedef decltype(T() * T2()) ValueType;
  SymmetricRankTwoTensorTempl<ValueType> result;

  std::size_t index = 0;
  for (const auto i : make_range(N))
  {
    ValueType tmp = 0.0;
    for (const auto j : make_range(N))
      tmp += _vals[index++] * b._vals[j];
    result._vals[i] = tmp;
  }

  return result;
}

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

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

  _vals[0] = input[0];  // C1111
  _vals[7] = input[3];  // C2222
  _vals[14] = input[5]; // C3333

  _vals[1] = _vals[6] = input[1];  // C1122
  _vals[2] = _vals[12] = input[2]; // C1133
  _vals[8] = _vals[13] = input[4]; // C2233

  static constexpr std::size_t C2323 = 21;
  static constexpr std::size_t C1313 = 28;
  static constexpr std::size_t C1212 = 35;

  _vals[C2323] = 2.0 * input[6];
  _vals[C1313] = 2.0 * input[7];
  _vals[C1212] = 2.0 * input[8];
}

template <typename T>
template <typename T2>
void
SymmetricRankFourTensorTempl<T>::fillSymmetric21FromInputVector(const T2 & input)
{
  // C1111 C1122 C1133 C1123 C1113 C1112
  //       C2222 C2233 C2223 C2213 C2212
  //             C3333 C3323 C3313 C3312
  //                   C2323 C2313 C2312
  //                         C1313 C1312
  //                               C1212

  mooseAssert(LIBMESH_DIM == 3, "This method assumes LIBMESH_DIM == 3");
  mooseAssert(input.size() == 21,
              "To use fillSymmetric21FromInputVector, your input must have size 21.");
  std::size_t index = 0;
  for (const auto i : make_range(N))
    for (const auto j : make_range(i, N))
    {
      _vals[i + N * j] = mandelFactor(i, j) * input[index];
      _vals[j + N * i] = mandelFactor(j, i) * input[index];
      index++;
    }
}

template <typename T>
SymmetricRankFourTensorTempl<T>
SymmetricRankFourTensorTempl<T>::invSymm() const
{
  SymmetricRankFourTensorTempl<T> result(initNone);

  if constexpr (SymmetricRankFourTensorTempl<T>::N2 * sizeof(T) > EIGEN_STACK_ALLOCATION_LIMIT)
  {
    Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> mat(N, N);
    for (const auto i : make_range(N))
      for (const auto j : make_range(N))
        mat(i, j) = (*this)(i, j);
    mat = mat.inverse();
    for (const auto i : make_range(N))
      for (const auto j : make_range(N))
        result(i, j) = mat(i, j);
  }
  else
  {
    const Eigen::Map<const Eigen::Matrix<T, N, N, Eigen::RowMajor>> mat(&_vals[0]);
    Eigen::Map<Eigen::Matrix<T, N, N, Eigen::RowMajor>> res(&result._vals[0]);
    res = mat.inverse();
  }

  return result;
}