RankThreeTensor.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 "ADRankThreeTensorForward.h"
#include "ADRankFourTensorForward.h"
 
#include "libmesh/libmesh.h"
#include "libmesh/int_range.h"
 
#include "metaphysicl/raw_type.h"
 
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<RankThreeTensorTempl<Real>>(RankThreeTensorTempl<Real> & v);
template <>
void mooseSetToZero<RankThreeTensorTempl<ADReal>>(RankThreeTensorTempl<ADReal> & v);
}
 
template <typename T>
class RankThreeTensorTempl
{
public:
  static constexpr unsigned int N = Moose::dim;
  static constexpr unsigned int N2 = N * N;
  static constexpr unsigned int N3 = N * N * N;
 
  enum InitMethod
  {
    initNone
  };
 
  enum FillMethod
  {
    automatic,
    general,
    plane_normal
  };
 
  RankThreeTensorTempl();
 
  RankThreeTensorTempl(const RankThreeTensorTempl<T> & a) = default;
 
  template <typename T2>
  RankThreeTensorTempl(const RankThreeTensorTempl<T2> & copy);
 
  RankThreeTensorTempl(const InitMethod);
 
  RankThreeTensorTempl(const std::vector<T> &, FillMethod method = automatic);
 
  inline T & operator()(unsigned int i, unsigned int j, unsigned int k)
  {
    return _vals[((i * N + j) * N + k)];
  }
 
  inline T operator()(unsigned int i, unsigned int j, unsigned int k) const
  {
    return _vals[((i * N + j) * N + k)];
  }
 
  RankThreeTensorTempl<T> & operator=(const T & value);
 
  void zero();
 
  void print(std::ostream & stm = Moose::out) const;
 
  friend std::ostream & operator<<(std::ostream & os, const RankThreeTensorTempl<T> & t)
  {
    t.print(os);
    return os;
  }
 
  RankThreeTensorTempl<T> & operator=(const RankThreeTensorTempl<T> & a);
 
  template <typename T2>
  RankThreeTensorTempl<T> & operator=(const RankThreeTensorTempl<T2> & a);
 
  libMesh::VectorValue<T> operator*(const RankTwoTensorTempl<T> & a) const;
 
  RankTwoTensorTempl<T> operator*(const libMesh::VectorValue<T> & a) const;
 
  RankThreeTensorTempl<T> operator*(const T a) const;
 
  RankThreeTensorTempl<T> & operator*=(const T a);
 
  RankThreeTensorTempl<T> operator/(const T a) const;
 
  RankThreeTensorTempl<T> & operator/=(const T a);
 
  RankThreeTensorTempl<T> & operator+=(const RankThreeTensorTempl<T> & a);
 
  RankThreeTensorTempl<T> operator+(const RankThreeTensorTempl<T> & a) const;
 
  RankThreeTensorTempl<T> & operator-=(const RankThreeTensorTempl<T> & a);
 
  RankThreeTensorTempl<T> operator-(const RankThreeTensorTempl<T> & a) const;
 
  RankThreeTensorTempl<T> operator-() const;
 
  T L2norm() const;
 
  template <class T2>
  void rotate(const T2 & R);
 
  void rotate(const libMesh::TensorValue<T> & R);
 
  static MooseEnum fillMethodEnum();
 
  void fillFromInputVector(const std::vector<T> & input, FillMethod fill_method = automatic);
 
  void fillFromPlaneNormal(const libMesh::VectorValue<T> & input);
 
  RankFourTensorTempl<T> mixedProductRankFour(const RankTwoTensorTempl<T> & a) const;
 
  libMesh::VectorValue<T> doubleContraction(const RankTwoTensorTempl<T> & b) const;
 
protected:
  T _vals[N3];
 
  void fillGeneralFromInputVector(const std::vector<T> & input);
 
  template <class T2>
  friend void dataStore(std::ostream &, RankThreeTensorTempl<T2> &, void *);
 
  template <class T2>
  friend void dataLoad(std::istream &, RankThreeTensorTempl<T2> &, void *);
 
  template <class T2>
  friend class RankTwoTensorTempl;
 
  template <class T2>
  friend class RankThreeTensorTempl;
 
  template <class T2>
  friend class RankFourTensorTempl;
};
 
namespace MetaPhysicL
{
template <typename T>
struct RawType<RankThreeTensorTempl<T>>
{
  typedef RankThreeTensorTempl<typename RawType<T>::value_type> value_type;
 
  static value_type value(const RankThreeTensorTempl<T> & in)
  {
    value_type ret;
    for (const auto i : libMesh::make_range(RankThreeTensorTempl<T>::N))
      for (const auto j : libMesh::make_range(RankThreeTensorTempl<T>::N))
        for (const auto k : libMesh::make_range(RankThreeTensorTempl<T>::N))
          ret(i, j, k) = raw_value(in(i, j, k));
 
    return ret;
  }
};
}
 
template <typename T>
template <typename T2>
RankThreeTensorTempl<T>::RankThreeTensorTempl(const RankThreeTensorTempl<T2> & copy)
{
  for (const auto i : libMesh::make_range(N3))
    _vals[i] = copy._vals[i];
}
 
template <typename T>
template <class T2>
void
RankThreeTensorTempl<T>::rotate(const T2 & R)
{
  unsigned int index = 0;
  for (const auto i : libMesh::make_range(N))
    for (const auto j : libMesh::make_range(N))
      for (const auto k : libMesh::make_range(N))
      {
        unsigned int index2 = 0;
        T sum = 0.0;
        for (const auto m : libMesh::make_range(N))
        {
          T a = R(i, m);
          for (const auto n : libMesh::make_range(N))
          {
            T ab = a * R(j, n);
            for (const auto o : libMesh::make_range(N))
              sum += ab * R(k, o) * _vals[index2++];
          }
        }
        _vals[index++] = sum;
      }
}
 
template <typename T>
RankThreeTensorTempl<T>
operator*(T a, const RankThreeTensorTempl<T> & b)
{
  return b * a;
}
 
template <typename T>
RankTwoTensorTempl<T>
operator*(const libMesh::VectorValue<T> & p, const RankThreeTensorTempl<T> & b)
{
  static_assert(RankThreeTensorTempl<T>::N == RankTwoTensorTempl<T>::N,
                "RankTwoTensor and RankThreeTensor have to have the same dimension N.");
  RankTwoTensorTempl<T> result;
 
  for (const auto i : libMesh::make_range(RankThreeTensorTempl<T>::N))
    for (const auto j : libMesh::make_range(RankThreeTensorTempl<T>::N))
      for (const auto k : libMesh::make_range(RankThreeTensorTempl<T>::N))
        result(i, j) += p(k) * b(k, i, j);
 
  return result;
}
 
template <typename T>
template <typename T2>
RankThreeTensorTempl<T> &
RankThreeTensorTempl<T>::operator=(const RankThreeTensorTempl<T2> & a)
{
  for (const auto i : libMesh::make_range(N))
    for (const auto j : libMesh::make_range(N))
      for (const auto k : libMesh::make_range(N))
        (*this)(i, j, k) = a(i, j, k);
 
  return *this;
}