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