RankTwoTensor.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 "MooseTypes.h"
#include "ADRankTwoTensorForward.h"
#include "ADRankFourTensorForward.h"
#include "ADRankThreeTensorForward.h"
#include "ADSymmetricRankTwoTensorForward.h"
#include "MooseUtils.h"
#include "MathUtils.h"

// Any requisite includes here
#include "libmesh/libmesh.h"
#include "libmesh/tensor_value.h"
#include "libmesh/vector_value.h"
#include "libmesh/int_range.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 TypeVector;
template <typename>
class TypeTensor;
template <typename>
class TensorValue;
namespace TensorTools
{
template <>
struct IncrementRank<RankTwoTensor>
{
  typedef RankThreeTensor type;
};
}
}

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

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

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

template <typename T>
class RankTwoTensorTempl : public libMesh::TensorValue<T>
{
public:

  static constexpr unsigned int N = Moose::dim;

  static constexpr unsigned int N2 = N * N;

  enum InitMethod
  {
    initNone,
    initIdentity
  };

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

  static void initRandom(unsigned int);

  [[nodiscard]] static MooseEnum fillMethodEnum();

  void fillRealTensor(libMesh::TensorValue<T> &);

  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;



  RankTwoTensorTempl();

  RankTwoTensorTempl(const InitMethod);

  RankTwoTensorTempl(const libMesh::TypeVector<T> & row1,
                     const libMesh::TypeVector<T> & row2,
                     const libMesh::TypeVector<T> & row3);

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

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

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

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

  RankTwoTensorTempl(const libMesh::TensorValue<T> & a) : libMesh::TensorValue<T>(a) {}

  RankTwoTensorTempl(const libMesh::TypeTensor<T> & a) : libMesh::TensorValue<T>(a) {}

  template <typename T2>
  RankTwoTensorTempl(const SymmetricRankTwoTensorTempl<T2> & a)
    : RankTwoTensorTempl<T>(a(0),
                            a(1),
                            a(2),
                            a(3) / MathUtils::sqrt2,
                            a(4) / MathUtils::sqrt2,
                            a(5) / MathUtils::sqrt2)
  {
  }

  template <typename T2>
  RankTwoTensorTempl(const RankTwoTensorTempl<T2> & a) : libMesh::TensorValue<T>(a)
  {
  }


  [[nodiscard]] static RankTwoTensorTempl initializeSymmetric(const libMesh::TypeVector<T> & v0,
                                                              const libMesh::TypeVector<T> & v1,
                                                              const libMesh::TypeVector<T> & v2);

  [[nodiscard]] static RankTwoTensorTempl initializeFromRows(const libMesh::TypeVector<T> & row0,
                                                             const libMesh::TypeVector<T> & row1,
                                                             const libMesh::TypeVector<T> & row2);

  [[nodiscard]] static RankTwoTensorTempl
  initializeFromColumns(const libMesh::TypeVector<T> & col0,
                        const libMesh::TypeVector<T> & col1,
                        const libMesh::TypeVector<T> & col2);

  [[nodiscard]] static RankTwoTensorTempl Identity() { return RankTwoTensorTempl(initIdentity); }

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

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

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

  [[nodiscard]] static RankTwoTensorTempl<T> outerProduct(const libMesh::TypeVector<T> &,
                                                          const libMesh::TypeVector<T> &);

  [[nodiscard]] static RankTwoTensorTempl<T> selfOuterProduct(const libMesh::TypeVector<T> &);

  [[nodiscard]] static RankTwoTensorTempl<T> genRandomTensor(T stddev, T mean);

  [[nodiscard]] static RankTwoTensorTempl<T> genRandomSymmTensor(T stddev, T mean);



  libMesh::VectorValue<T> column(const unsigned int i) const;



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

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

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

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

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

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

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

  RankTwoTensorTempl<T> & operator*=(const libMesh::TypeTensor<T> & a);

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

  void fillFromScalarVariable(const VariableValue & scalar_variable);

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

  void vectorOuterProduct(const libMesh::TypeVector<T> &, const libMesh::TypeVector<T> &);

  void fillRow(unsigned int r, const libMesh::TypeVector<T> & v);

  void fillColumn(unsigned int c, const libMesh::TypeVector<T> & v);

  void setToIdentity();


  void addIa(const T & a);

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



  RankTwoTensorTempl<T> square() const;

  RankTwoTensorTempl<T> rotated(const RankTwoTensorTempl<T> & R) const;

  RankTwoTensorTempl<T> rotateXyPlane(T a);

  RankTwoTensorTempl<T> transpose() const;

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

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

  RankTwoTensorTempl<T> operator-() const;

  template <typename T2, typename std::enable_if<libMesh::ScalarTraits<T2>::value, int>::type = 0>
  RankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
  operator*(const T2 & a) const;

  template <typename T2, typename std::enable_if<libMesh::ScalarTraits<T2>::value, int>::type = 0>
  RankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
  operator/(const T2 & a) const;

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

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

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

  template <int n, int o, int p, int q>
  RankFourTensorTempl<T> times(const RankTwoTensorTempl<T> & b) const;

  template <int n, int o, int p, int q, int r, int s>
  RankFourTensorTempl<T> times(const RankFourTensorTempl<T> & b) const;

  RankFourTensorTempl<T> outerProduct(const RankTwoTensorTempl<T> & b) const
  {
    usingTensorIndices(i_, j_, k_, l_);
    return times<i_, j_, k_, l_>(b);
  }

  RankThreeTensorTempl<T> contraction(const RankThreeTensorTempl<T> & b) const;

  RankThreeTensorTempl<T> mixedProductJkI(const libMesh::VectorValue<T> & b) const;

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

  RankTwoTensorTempl<T> deviatoric() const;

  T trace() const;

  RankTwoTensorTempl<T> dtrace() const;

  RankTwoTensorTempl<T> inverse() const;

  T generalSecondInvariant() const;

  T secondInvariant() const;

  RankTwoTensorTempl<T> dsecondInvariant() const;

  RankFourTensorTempl<T> d2secondInvariant() const;

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

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

  RankFourTensorTempl<T> d2sin3Lode(const T & r0) const;

  T thirdInvariant() const;

  RankTwoTensorTempl<T> dthirdInvariant() const;

  RankFourTensorTempl<T> d2thirdInvariant() const;

  RankTwoTensorTempl<T> ddet() const;

  T L2norm() const;

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

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

  void dsymmetricEigenvalues(std::vector<T> & eigvals,
                             std::vector<RankTwoTensorTempl<T>> & deigvals) const;

  void d2symmetricEigenvalues(std::vector<RankFourTensorTempl<T>> & deriv) const;

  void getRUDecompositionRotation(RankTwoTensorTempl<T> & rot) const;

  RankTwoTensorTempl<T> initialContraction(const RankFourTensorTempl<T> & b) const;



  bool operator==(const RankTwoTensorTempl<T> & a) const;

  bool isSymmetric() const;


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

private:
  static constexpr Real identityCoords[N2] = {1, 0, 0, 0, 1, 0, 0, 0, 1};

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

  using libMesh::TensorValue<T>::_coords;

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

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

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

    return ret;
  }
};
}

template <typename T>
template <typename T2>
RankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
RankTwoTensorTempl<T>::operator+(const libMesh::TypeTensor<T2> & b) const
{
  return libMesh::TensorValue<T>::operator+(b);
}

template <typename T>
template <typename T2>
RankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
RankTwoTensorTempl<T>::operator-(const libMesh::TypeTensor<T2> & b) const
{
  return libMesh::TensorValue<T>::operator-(b);
}

template <typename T>
template <typename T2, typename std::enable_if<libMesh::ScalarTraits<T2>::value, int>::type>
RankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
RankTwoTensorTempl<T>::operator*(const T2 & b) const
{
  return libMesh::TensorValue<T>::operator*(b);
}

template <typename T>
template <typename T2>
libMesh::TypeVector<typename libMesh::CompareTypes<T, T2>::supertype>
RankTwoTensorTempl<T>::operator*(const libMesh::TypeVector<T2> & b) const
{
  return libMesh::TensorValue<T>::operator*(b);
}

template <typename T>
template <typename T2>
RankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
RankTwoTensorTempl<T>::operator*(const libMesh::TypeTensor<T2> & b) const
{
  return libMesh::TensorValue<T>::operator*(b);
}

template <typename T>
template <typename T2, typename std::enable_if<libMesh::ScalarTraits<T2>::value, int>::type>
RankTwoTensorTempl<typename libMesh::CompareTypes<T, T2>::supertype>
RankTwoTensorTempl<T>::operator/(const T2 & b) const
{
  return libMesh::TensorValue<T>::operator/(b);
}

template <typename T>
RankFourTensorTempl<T>
RankTwoTensorTempl<T>::positiveProjectionEigenDecomposition(std::vector<T> & eigval,
                                                            RankTwoTensorTempl<T> & eigvec) const
{
  using std::abs;
  if constexpr (MooseUtils::IsLikeReal<T>::value)
  {
    // Compute eigenvectors and eigenvalues of this tensor
    this->symmetricEigenvaluesEigenvectors(eigval, eigvec);

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

    // projection tensor
    RankFourTensorTempl<T> proj_pos;
    RankFourTensorTempl<T> Gab, Gba;

    for (auto a : libMesh::make_range(N))
    {
      const auto Ma = RankTwoTensorTempl<T>::selfOuterProduct(eigvec.column(a));
      proj_pos += d[a] * Ma.outerProduct(Ma);
    }

    usingTensorIndices(i_, j_, k_, l_);
    for (const auto a : libMesh::make_range(N))
      for (const auto b : libMesh::make_range(a))
      {
        const auto Ma = RankTwoTensorTempl<T>::selfOuterProduct(eigvec.column(a));
        const auto Mb = RankTwoTensorTempl<T>::selfOuterProduct(eigvec.column(b));

        Gab = Ma.template times<i_, k_, j_, l_>(Mb) + Ma.template times<i_, l_, j_, k_>(Mb);
        Gba = Mb.template times<i_, k_, j_, l_>(Ma) + Mb.template times<i_, l_, j_, k_>(Ma);

        T theta_ab;
        if (!MooseUtils::absoluteFuzzyEqual(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 * (Gab + Gba);
      }
    return proj_pos;
  }
  else
    mooseError("positiveProjectionEigenDecomposition is only available for ordered tensor "
               "component types");
}

template <typename T>
T
RankTwoTensorTempl<T>::sin3Lode(const T & r0, const T & r0_value) const
{
  using std::pow, std::max, std::sqrt, std::min;
  if constexpr (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(-1.5 * sqrt(3.0) * thirdInvariant() / pow(bar, 1.5), 1.0), -1.0);
  }
  else
    mooseError("sin3Lode is only available for ordered tensor component types");
}

template <typename T>
RankTwoTensorTempl<T>
RankTwoTensorTempl<T>::dsin3Lode(const T & r0) const
{
  using std::sqrt, std::pow;
  if constexpr (MooseUtils::IsLikeReal<T>::value)
  {
    T bar = secondInvariant();
    if (bar <= r0)
      return RankTwoTensorTempl<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>
RankFourTensorTempl<T>
RankTwoTensorTempl<T>::d2sin3Lode(const T & r0) const
{
  using std::pow, std::sqrt;
  if constexpr (MooseUtils::IsLikeReal<T>::value)
  {
    T bar = secondInvariant();
    if (bar <= r0)
      return RankFourTensorTempl<T>();

    T J3 = thirdInvariant();
    RankTwoTensorTempl<T> dII = dsecondInvariant();
    RankTwoTensorTempl<T> dIII = dthirdInvariant();
    RankFourTensorTempl<T> deriv =
        d2thirdInvariant() / pow(bar, 1.5) - 1.5 * d2secondInvariant() * J3 / pow(bar, 2.5);

    for (unsigned i = 0; i < N; ++i)
      for (unsigned j = 0; j < N; ++j)
        for (unsigned k = 0; k < N; ++k)
          for (unsigned l = 0; l < N; ++l)
            deriv(i, j, k, l) +=
                (-1.5 * dII(i, j) * dIII(k, l) - 1.5 * dIII(i, j) * dII(k, l)) / pow(bar, 2.5) +
                1.5 * 2.5 * dII(i, j) * dII(k, l) * J3 / pow(bar, 3.5);

    deriv *= -1.5 * sqrt(3.0);
    return deriv;
  }
  else
    mooseError("d2sin3Lode is only available for ordered tensor component types");
}

template <typename T>
template <int n, int o, int p, int q>
RankFourTensorTempl<T>
RankTwoTensorTempl<T>::times(const RankTwoTensorTempl<T> & b) const
{
  RankFourTensorTempl<T> result;
  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[0], x[1], x[2], x[3]) = (*this)(x[n], x[o]) * b(x[p], x[q]);

  return result;
}

template <typename T>
template <int n, int o, int p, int q, int r, int s>
RankFourTensorTempl<T>
RankTwoTensorTempl<T>::times(const RankFourTensorTempl<T> & b) const
{
  RankFourTensorTempl<T> result;
  std::size_t x[5];
  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])
          for (x[4] = 0; x[4] < N; ++x[4])
            result(x[0], x[1], x[2], x[3]) += (*this)(x[n], x[o]) * b(x[p], x[q], x[r], x[s]);

  return result;
}