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 std::enable_if<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;
}