SymmetricRankTwoTensor.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 "ADSymmetricRankTwoTensorForward.h"
#include "ADSymmetricRankFourTensorForward.h"
#include "MooseUtils.h"
#include "MathUtils.h"

#include "libmesh/libmesh.h"
#include "libmesh/tensor_value.h"

#include "metaphysicl/raw_type.h"

#include <petscsys.h>
#include <vector>

// Forward declarations
class MooseEnum;
template <typename T>
class MooseArray;
typedef MooseArray<Real> VariableValue;
template <typename>
class ColumnMajorMatrixTempl;
namespace libMesh
{
template <typename>
class VectorValue;
template <typename>
class TypeVector;
template <typename>
class TypeTensor;
template <typename>
class TensorValue;
}

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

template <>
void mooseSetToZero<SymmetricRankTwoTensor>(SymmetricRankTwoTensor & v);

template <>
void mooseSetToZero<ADSymmetricRankTwoTensor>(ADSymmetricRankTwoTensor & v);
}

template <typename T>
class SymmetricRankTwoTensorTempl
{
public:
  typedef T value_type;

  static constexpr unsigned int Ndim = Moose::dim;
  static constexpr unsigned int N = Ndim + (Ndim * (Ndim - 1)) / 2;

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

  // Reverse tensor indices in the Mandel/Voigt representation
  static constexpr unsigned int reverse_index[3][3] = {{0, 5, 4}, {5, 1, 3}, {4, 3, 2}};

  static constexpr Real mandelFactor(unsigned int i) { return i < Ndim ? 1.0 : MathUtils::sqrt2; }

  // Select initialization
  enum InitMethod
  {
    initNone,
    initIdentity
  };

  SymmetricRankTwoTensorTempl();

  SymmetricRankTwoTensorTempl(const InitMethod);

  enum FillMethod
  {
    autodetect = 0,
    isotropic1 = 1,
    diagonal3 = 3,
    symmetric6 = 6
  };

  SymmetricRankTwoTensorTempl(const std::vector<T> & input) { this->fillFromInputVector(input); };

  SymmetricRankTwoTensorTempl(
      const T & S11, const T & S22, const T & S33, const T & S23, const T & S13, const T & S12);

  // explicit cast to a full tensor
  explicit operator RankTwoTensorTempl<T>();

private:
  SymmetricRankTwoTensorTempl(const T & S11,
                              const T & S21,
                              const T & S31,
                              const T & S12,
                              const T & S22,
                              const T & S32,
                              const T & S13,
                              const T & S23,
                              const T & S33);

  // internal use named constructor that does not apply the sqrt(2) factors
  static SymmetricRankTwoTensorTempl fromRawComponents(
      const T & S11, const T & S22, const T & S33, const T & S23, const T & S13, const T & S12);

public:
  explicit SymmetricRankTwoTensorTempl(const TensorValue<T> & a);

  explicit SymmetricRankTwoTensorTempl(const TypeTensor<T> & a);

  template <typename T2>
  SymmetricRankTwoTensorTempl(const SymmetricRankTwoTensorTempl<T2> & a)
  {
    for (const auto i : make_range(N))
      _vals[i] = a(i);
  }

  // Named constructors
  static SymmetricRankTwoTensorTempl identity()
  {
    return SymmetricRankTwoTensorTempl(initIdentity);
  }

  static SymmetricRankTwoTensorTempl
  initializeSymmetric(const TypeVector<T> & v0, const TypeVector<T> & v1, const TypeVector<T> & v2);

  static MooseEnum fillMethodEnum();

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

  void fillFromScalarVariable(const VariableValue & scalar_variable);

  inline T & operator()(unsigned int i) { return _vals[i]; }

  inline T operator()(unsigned int i) const { return _vals[i]; }

  inline T operator()(unsigned int i, unsigned int j) const
  {
    const auto a = reverse_index[i][j];
    return _vals[a] / mandelFactor(a);
  }

  libMesh::VectorValue<T> row(const unsigned int n) const;

  libMesh::VectorValue<T> column(const unsigned int n) const { return row(n); }

  [[nodiscard]] static SymmetricRankTwoTensorTempl<T> timesTranspose(const RankTwoTensorTempl<T> &);
  [[nodiscard]] static SymmetricRankTwoTensorTempl<T>
  timesTranspose(const SymmetricRankTwoTensorTempl<T> &);

  [[nodiscard]] static SymmetricRankTwoTensorTempl<T> transposeTimes(const RankTwoTensorTempl<T> &);
  [[nodiscard]] static SymmetricRankTwoTensorTempl<T>
  transposeTimes(const SymmetricRankTwoTensorTempl<T> &);

  [[nodiscard]] static SymmetricRankTwoTensorTempl<T> plusTranspose(const RankTwoTensorTempl<T> &);
  [[nodiscard]] static SymmetricRankTwoTensorTempl<T>
  plusTranspose(const SymmetricRankTwoTensorTempl<T> &);

  SymmetricRankTwoTensorTempl<T> square() const;

  T tr() const { return _vals[0] + _vals[1] + _vals[2]; }

  void zero();

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

  SymmetricRankTwoTensorTempl<T> transpose() const;

  template <typename T2>
  SymmetricRankTwoTensorTempl<T> & operator=(const SymmetricRankTwoTensorTempl<T2> & a);

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

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

  template <typename T2>
  SymmetricRankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
  operator+(const SymmetricRankTwoTensorTempl<T2> & a) const;

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

  template <typename T2>
  SymmetricRankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
  operator-(const SymmetricRankTwoTensorTempl<T2> & a) const;

  SymmetricRankTwoTensorTempl<T> operator-() const;

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

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

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

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

  template <typename T2>
  TypeVector<typename libMesh::CompareTypes<T, T2>::supertype>
  operator*(const TypeVector<T2> & a) const;

  template <typename T2>
  bool operator==(const SymmetricRankTwoTensorTempl<T2> & a) const;

  bool isSymmetric() const { return true; }

  template <typename T2>
  bool operator!=(const SymmetricRankTwoTensorTempl<T2> & a) const;

  SymmetricRankTwoTensorTempl<T> & operator=(const ColumnMajorMatrixTempl<T> & a);

  T doubleContraction(const SymmetricRankTwoTensorTempl<T> & a) const;

  SymmetricRankFourTensorTempl<T> outerProduct(const SymmetricRankTwoTensorTempl<T> & a) const;

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

  SymmetricRankTwoTensorTempl<T> deviatoric() const;

  T trace() const;

  SymmetricRankTwoTensorTempl<T> inverse() const;

  SymmetricRankTwoTensorTempl<T> dtrace() const;

  T generalSecondInvariant() const;

  T secondInvariant() const;

  SymmetricRankTwoTensorTempl<T> dsecondInvariant() const;

  SymmetricRankFourTensorTempl<T> d2secondInvariant() const;

  T sin3Lode(const T & r0, const T & r0_value) const;

  SymmetricRankTwoTensorTempl<T> dsin3Lode(const T & r0) const;

  T thirdInvariant() const;

  SymmetricRankTwoTensorTempl<T> dthirdInvariant() const;

  SymmetricRankFourTensorTempl<T> d2thirdInvariant() const;

  // determinant of the tensor
  T det() const;

  SymmetricRankTwoTensorTempl<T> ddet() const;

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

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

  void printADReal(unsigned int nDual, std::ostream & stm = Moose::out) const;

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

  void addIa(const T & a);

  T L2norm() const;

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

  void symmetricEigenvalues(std::vector<T> & eigvals) const;

  void symmetricEigenvaluesEigenvectors(std::vector<T> & eigvals,
                                        RankTwoTensorTempl<T> & eigvecs) const;

  static void initRandom(unsigned int);

  [[nodiscard]] static SymmetricRankTwoTensorTempl<T> genRandomSymmTensor(T, T);

  [[nodiscard]] static SymmetricRankTwoTensorTempl<T> selfOuterProduct(const TypeVector<T> &);

  SymmetricRankTwoTensorTempl<T>
  initialContraction(const SymmetricRankFourTensorTempl<T> & b) const;

  void setToIdentity();

protected:
  void syev(const char * calculation_type, std::vector<T> & eigvals, std::vector<T> & a) const;

private:
  static constexpr std::array<Real, N> identityCoords = {{1, 1, 1, 0, 0, 0}};

  // tensor components
  std::array<T, N> _vals;

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

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

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

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

    return ret;
  }
};
}

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

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

template <typename T>
template <typename T2>
TypeVector<typename libMesh::CompareTypes<T, T2>::supertype>
SymmetricRankTwoTensorTempl<T>::operator*(const TypeVector<T2> & a) const
{
  TypeVector<typename libMesh::CompareTypes<T, T2>::supertype> ret;
  ret(0) = a(0) * _vals[0] + a(1) * _vals[5] + a(2) * _vals[4];
  ret(1) = a(0) * _vals[5] + a(1) * _vals[1] + a(2) * _vals[3];
  ret(2) = a(0) * _vals[4] + a(1) * _vals[3] + a(2) * _vals[2];
}

template <typename T>
template <typename T2>
bool
SymmetricRankTwoTensorTempl<T>::operator==(const SymmetricRankTwoTensorTempl<T2> & a) const
{
  for (std::size_t i = 0; i < N; ++i)
    if (_vals[i] != a._vals[i])
      return false;
  return true;
}

template <typename T>
template <typename T2>
bool
SymmetricRankTwoTensorTempl<T>::operator!=(const SymmetricRankTwoTensorTempl<T2> & a) const
{
  for (std::size_t i = 0; i < N; ++i)
    if (_vals[i] != a._vals[i])
      return true;
  return false;
}

template <typename T>
SymmetricRankFourTensorTempl<T>
SymmetricRankTwoTensorTempl<T>::positiveProjectionEigenDecomposition(
    std::vector<T> & eigval, RankTwoTensorTempl<T> & eigvec) const
{
  using std::abs;

  // The calculate of projection tensor follows
  // C. Miehe and M. Lambrecht, Commun. Numer. Meth. Engng 2001; 17:337~353

  // Compute eigenvectors and eigenvalues of this tensor
  if (MooseUtils::IsLikeReal<T>::value)
  {
    this->symmetricEigenvaluesEigenvectors(eigval, eigvec);

    // Separate out positive and negative eigen values
    std::array<T, N> epos;
    std::array<T, N> d;
    for (unsigned int i = 0; i < Ndim; ++i)
    {
      epos[i] = (abs(eigval[i]) + eigval[i]) / 2.0;
      d[i] = 0 < eigval[i] ? 1.0 : 0.0;
    }

    // projection tensor
    SymmetricRankFourTensorTempl<T> proj_pos;

    for (unsigned int a = 0; a < Ndim; ++a)
    {
      const auto Ma = SymmetricRankTwoTensorTempl<T>::selfOuterProduct(eigvec.column(a));
      proj_pos += d[a] * Ma.outerProduct(Ma);
    }

    for (const auto a : make_range(Ndim))
      for (const auto b : make_range(a))
      {
        const auto Ma = SymmetricRankTwoTensorTempl<T>::selfOuterProduct(eigvec.column(a));
        const auto Mb = SymmetricRankTwoTensorTempl<T>::selfOuterProduct(eigvec.column(b));

        SymmetricRankFourTensorTempl<T> Gabba(SymmetricRankFourTensorTempl<T>::initNone);
        for (const auto aa : make_range(N))
          for (const auto bb : make_range(N))
          {
            const auto i = SymmetricRankFourTensorTempl<T>::full_index[aa][bb][0];
            const auto j = SymmetricRankFourTensorTempl<T>::full_index[aa][bb][1];
            const auto k = SymmetricRankFourTensorTempl<T>::full_index[aa][bb][2];
            const auto l = SymmetricRankFourTensorTempl<T>::full_index[aa][bb][3];

            Gabba(aa, bb) = (Ma(i, k) * Mb(j, l) + Ma(i, l) * Mb(j, k) + Ma(j, l) * Mb(i, k) +
                             Ma(j, k) * Mb(i, l)) *
                            SymmetricRankFourTensorTempl<T>::mandelFactor(aa, bb);
          }

        T theta_ab;
        if (!MooseUtils::relativeFuzzyEqual(eigval[a], eigval[b]))
          theta_ab = 0.5 * (epos[a] - epos[b]) / (eigval[a] - eigval[b]);
        else
          theta_ab = 0.25 * (d[a] + d[b]);

        proj_pos += theta_ab * Gabba;
      }
    return proj_pos;
  }
  else
    mooseError("positiveProjectionEigenDecomposition is only available for ordered tensor "
               "component types");
}

template <typename T>
T
SymmetricRankTwoTensorTempl<T>::sin3Lode(const T & r0, const T & r0_value) const
{
  using std::max, std::min, std::sqrt, std::pow;
  if (MooseUtils::IsLikeReal<T>::value)
  {
    T bar = secondInvariant();
    if (bar <= r0)
      // in this case the Lode angle is not defined
      return r0_value;
    else
      // the min and max here gaurd against precision-loss when bar is tiny but nonzero.
      return max(min(thirdInvariant() * -1.5 * sqrt(3.0) / pow(bar, 1.5), 1.0), -1.0);
  }
  else
    mooseError("sin3Lode is only available for ordered tensor component types");
}

template <typename T>
SymmetricRankTwoTensorTempl<T>
SymmetricRankTwoTensorTempl<T>::dsin3Lode(const T & r0) const
{
  using std::sqrt, std::pow;
  if (MooseUtils::IsLikeReal<T>::value)
  {
    T bar = secondInvariant();
    if (bar <= r0)
      return SymmetricRankTwoTensorTempl<T>();
    else
      return -1.5 * sqrt(3.0) *
             (dthirdInvariant() / pow(bar, 1.5) -
              1.5 * dsecondInvariant() * thirdInvariant() / pow(bar, 2.5));
  }
  else
    mooseError("dsin3Lode is only available for ordered tensor component types");
}

template <typename T>
template <typename T2>
SymmetricRankTwoTensorTempl<T> &
SymmetricRankTwoTensorTempl<T>::operator=(const SymmetricRankTwoTensorTempl<T2> & a)
{
  for (const auto i : make_range(N))
    (*this)(i) = a(i);
  return *this;
}

// non-member operators

template <typename T, typename Scalar>
inline typename std::enable_if_t<
    libMesh::ScalarTraits<Scalar>::value,
    SymmetricRankTwoTensorTempl<typename libMesh::CompareTypes<T, Scalar>::supertype>>
operator*(const Scalar & factor, const SymmetricRankTwoTensorTempl<T> & t)
{
  return t * factor;
}