19 #include "libmesh/libmesh.h" 20 #include "libmesh/tuple_of.h" 21 #include "libmesh/int_range.h" 23 #include "metaphysicl/raw_type.h" 28 #include <Eigen/Dense> 70 static constexpr
unsigned int N2 =
N *
N;
71 static constexpr
unsigned int N3 =
N *
N *
N;
72 static constexpr
unsigned int N4 =
N *
N *
N *
N;
108 template <
template <
typename>
class Tensor,
typename Scalar>
113 template <
typename Scalar>
118 template <
typename Scalar>
123 template <
typename Scalar>
144 template <
typename T2>
150 template <
typename T2>
157 const auto i =
idx[0];
158 const auto j =
idx[1];
159 const auto k =
idx[2];
160 const auto l =
idx[3];
180 inline T &
operator()(
unsigned int i,
unsigned int j,
unsigned int k,
unsigned int l)
189 inline const T &
operator()(
unsigned int i,
unsigned int j,
unsigned int k,
unsigned int l)
const 198 void print(std::ostream & stm = Moose::out)
const;
200 friend std::ostream & operator<<(std::ostream & os, const RankFourTensorTempl<T> & t)
207 void printReal(std::ostream & stm = Moose::out)
const;
217 template <
typename Scalar>
222 libmesh_assert_equal_to(p, Scalar(0));
228 template <
template <
typename>
class Tensor,
typename T2>
230 typename std::enable_if<TwoTensorMultTraits<Tensor, T2>::value,
234 template <
typename T2>
236 typename std::enable_if<libMesh::ScalarTraits<T2>::value,
243 template <
typename T2>
245 typename std::enable_if<libMesh::ScalarTraits<T2>::value,
255 template <
typename T2>
263 template <
typename T2>
271 template <
typename T2>
294 void rotate(
const TypeTensor<T> & R);
372 template <
typename T2>
383 template <
typename T2>
530 template <
typename T2>
532 template <
typename T2>
534 template <
typename T2>
540 template <
typename T>
553 ret(i, j, k, l) =
raw_value(in(i, j, k, l));
560 template <
typename T1,
typename T2>
563 typename std::enable_if<libMesh::ScalarTraits<T1>::value,
569 template <
typename T>
570 template <
typename T2>
574 _vals[i] = copy.
_vals[i];
577 template <
typename T>
578 template <
typename T2>
581 typename std::enable_if<libMesh::ScalarTraits<T2>::value,
584 typedef decltype(T() * T2()) ValueType;
588 result.
_vals[i] = _vals[i] * b;
593 template <
typename T>
594 template <
typename T2>
597 typename std::enable_if<libMesh::ScalarTraits<T2>::value,
602 result.
_vals[i] = _vals[i] / b;
606 template <
typename T>
607 template <
typename T2>
611 mooseAssert(input.size() == 9,
612 "To use fillSymmetric9FromInputVector, your input must have size 9.");
615 (*this)(0, 0, 0, 0) = input[0];
616 (*this)(1, 1, 1, 1) = input[3];
617 (*this)(2, 2, 2, 2) = input[5];
619 (*this)(0, 0, 1, 1) = input[1];
620 (*this)(1, 1, 0, 0) = input[1];
622 (*this)(0, 0, 2, 2) = input[2];
623 (*this)(2, 2, 0, 0) = input[2];
625 (*this)(1, 1, 2, 2) = input[4];
626 (*this)(2, 2, 1, 1) = input[4];
628 (*this)(1, 2, 1, 2) = input[6];
629 (*this)(2, 1, 2, 1) = input[6];
630 (*this)(2, 1, 1, 2) = input[6];
631 (*this)(1, 2, 2, 1) = input[6];
633 (*this)(0, 2, 0, 2) = input[7];
634 (*this)(2, 0, 2, 0) = input[7];
635 (*this)(2, 0, 0, 2) = input[7];
636 (*this)(0, 2, 2, 0) = input[7];
638 (*this)(0, 1, 0, 1) = input[8];
639 (*this)(1, 0, 1, 0) = input[8];
640 (*this)(1, 0, 0, 1) = input[8];
641 (*this)(0, 1, 1, 0) = input[8];
643 template <
typename T>
644 template <
typename T2>
648 mooseAssert(input.size() == 21,
649 "To use fillSymmetric21FromInputVector, your input must have size 21.");
651 (*this)(0, 0, 0, 0) = input[0];
652 (*this)(1, 1, 1, 1) = input[6];
653 (*this)(2, 2, 2, 2) = input[11];
655 (*this)(0, 0, 1, 1) = input[1];
656 (*this)(1, 1, 0, 0) = input[1];
658 (*this)(0, 0, 2, 2) = input[2];
659 (*this)(2, 2, 0, 0) = input[2];
661 (*this)(1, 1, 2, 2) = input[7];
662 (*this)(2, 2, 1, 1) = input[7];
664 (*this)(0, 0, 0, 2) = input[4];
665 (*this)(0, 0, 2, 0) = input[4];
666 (*this)(0, 2, 0, 0) = input[4];
667 (*this)(2, 0, 0, 0) = input[4];
669 (*this)(0, 0, 0, 1) = input[5];
670 (*this)(0, 0, 1, 0) = input[5];
671 (*this)(0, 1, 0, 0) = input[5];
672 (*this)(1, 0, 0, 0) = input[5];
674 (*this)(1, 1, 1, 2) = input[8];
675 (*this)(1, 1, 2, 1) = input[8];
676 (*this)(1, 2, 1, 1) = input[8];
677 (*this)(2, 1, 1, 1) = input[8];
679 (*this)(1, 1, 1, 0) = input[10];
680 (*this)(1, 1, 0, 1) = input[10];
681 (*this)(1, 0, 1, 1) = input[10];
682 (*this)(0, 1, 1, 1) = input[10];
684 (*this)(2, 2, 2, 1) = input[12];
685 (*this)(2, 2, 1, 2) = input[12];
686 (*this)(2, 1, 2, 2) = input[12];
687 (*this)(1, 2, 2, 2) = input[12];
689 (*this)(2, 2, 2, 0) = input[13];
690 (*this)(2, 2, 0, 2) = input[13];
691 (*this)(2, 0, 2, 2) = input[13];
692 (*this)(0, 2, 2, 2) = input[13];
694 (*this)(0, 0, 1, 2) = input[3];
695 (*this)(0, 0, 2, 1) = input[3];
696 (*this)(1, 2, 0, 0) = input[3];
697 (*this)(2, 1, 0, 0) = input[3];
699 (*this)(1, 1, 0, 2) = input[9];
700 (*this)(1, 1, 2, 0) = input[9];
701 (*this)(0, 2, 1, 1) = input[9];
702 (*this)(2, 0, 1, 1) = input[9];
704 (*this)(2, 2, 0, 1) = input[14];
705 (*this)(2, 2, 1, 0) = input[14];
706 (*this)(0, 1, 2, 2) = input[14];
707 (*this)(1, 0, 2, 2) = input[14];
709 (*this)(1, 2, 1, 2) = input[15];
710 (*this)(2, 1, 2, 1) = input[15];
711 (*this)(2, 1, 1, 2) = input[15];
712 (*this)(1, 2, 2, 1) = input[15];
714 (*this)(0, 2, 0, 2) = input[18];
715 (*this)(2, 0, 2, 0) = input[18];
716 (*this)(2, 0, 0, 2) = input[18];
717 (*this)(0, 2, 2, 0) = input[18];
719 (*this)(0, 1, 0, 1) = input[20];
720 (*this)(1, 0, 1, 0) = input[20];
721 (*this)(1, 0, 0, 1) = input[20];
722 (*this)(0, 1, 1, 0) = input[20];
724 (*this)(1, 2, 0, 2) = input[16];
725 (*this)(0, 2, 1, 2) = input[16];
726 (*this)(2, 1, 0, 2) = input[16];
727 (*this)(1, 2, 2, 0) = input[16];
728 (*this)(2, 0, 1, 2) = input[16];
729 (*this)(0, 2, 2, 1) = input[16];
730 (*this)(2, 1, 2, 0) = input[16];
731 (*this)(2, 0, 2, 1) = input[16];
733 (*this)(1, 2, 0, 1) = input[17];
734 (*this)(0, 1, 1, 2) = input[17];
735 (*this)(2, 1, 0, 1) = input[17];
736 (*this)(1, 2, 1, 0) = input[17];
737 (*this)(1, 0, 1, 2) = input[17];
738 (*this)(0, 1, 2, 1) = input[17];
739 (*this)(2, 1, 1, 0) = input[17];
740 (*this)(1, 0, 2, 1) = input[17];
742 (*this)(0, 2, 0, 1) = input[19];
743 (*this)(0, 1, 0, 2) = input[19];
744 (*this)(2, 0, 0, 1) = input[19];
745 (*this)(0, 2, 1, 0) = input[19];
746 (*this)(1, 0, 0, 2) = input[19];
747 (*this)(0, 1, 2, 0) = input[19];
748 (*this)(2, 0, 1, 0) = input[19];
749 (*this)(1, 0, 2, 0) = input[19];
752 template <
typename T>
767 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> mat(9, 9);
774 result.
_vals[i] = mat(i);
779 const Eigen::Map<const Eigen::Matrix<T, 9, 9, Eigen::RowMajor>> mat(&_vals[0]);
780 Eigen::Map<Eigen::Matrix<T, 9, 9, Eigen::RowMajor>> res(&result.
_vals[0]);
787 template <
typename T>
793 static constexpr std::size_t z[4][3] = {{1, 2, 3}, {0, 2, 3}, {0, 1, 3}, {0, 1, 2}};
795 for (x[0] = 0; x[0] <
N; ++x[0])
796 for (x[1] = 0; x[1] <
N; ++x[1])
797 for (x[2] = 0; x[2] <
N; ++x[2])
798 for (x[3] = 0; x[3] <
N; ++x[3])
799 result(x[z[m][0]], x[z[m][1]], x[z[m][2]]) += (*this)(x[0], x[1], x[2], x[3]) * b(x[m]);
RankFourTensorTempl is designed to handle any N-dimensional fourth order tensor, C.
RankFourTensorTempl< T > inverse() const
This returns A_ijkl such that C_ijkl*A_klmn = de_im de_jn i.e.
RankFourTensorTempl< T > singleProductJ(const RankTwoTensorTempl< T > &) const
Calculates C_imkl A_jm.
void mooseSetToZero(T &v)
Helper function templates to set a variable to zero.
void fillSymmetric9FromInputVector(const T2 &input)
fillSymmetric9FromInputVector takes 9 inputs to fill in the Rank-4 tensor with the appropriate crysta...
void fillGeneralFromInputVector(const std::vector< T > &input)
fillGeneralFromInputVector takes 81 inputs to fill the Rank-4 tensor No symmetries are explicitly mai...
void fillGeneralOrthotropicFromInputVector(const std::vector< T > &input)
fillGeneralOrhotropicFromInputVector takes 10 inputs to fill the Rank-4 tensor It defines a general o...
void fillAntisymmetricFromInputVector(const std::vector< T > &input)
fillAntisymmetricFromInputVector takes 6 inputs to fill the the antisymmetric Rank-4 tensor with the ...
RankFourTensorTempl< T > singleProductL(const RankTwoTensorTempl< T > &) const
Calculates C_ijkm A_lm.
RankFourTensorTempl< T > & operator=(const RankFourTensorTempl< T > &a)=default
copies values from a into this tensor
friend void dataStore(std::ostream &, RankFourTensorTempl< T2 > &, void *)
static RankFourTensorTempl< T > IdentityFour()
RankFourTensorTempl< T > tripleProductIjl(const RankTwoTensorTempl< T > &, const RankTwoTensorTempl< T > &, const RankTwoTensorTempl< T > &) const
Calculates C_mnkt A_im B_jn C_lt.
RankFourTensorTempl< T > & operator*=(const T &a)
C_ijkl *= a.
void fillSymmetric21FromInputVector(const T2 &input)
fillSymmetric21FromInputVector takes either 21 inputs to fill in the Rank-4 tensor with the appropria...
RankTwoTensorTempl< T > innerProductTranspose(const RankTwoTensorTempl< T > &) const
Inner product of the major transposed tensor with a rank two tensor.
static constexpr unsigned int N3
tuple_of< 4, unsigned int > index_type
static RankFourTensorTempl< T > Identity()
auto operator*(const T1 &a, const RankFourTensorTempl< T2 > &b) -> typename std::enable_if< libMesh::ScalarTraits< T1 >::value, RankFourTensorTempl< decltype(T1() *T2())>>::type
friend void dataLoad(std::istream &, RankFourTensorTempl< T2 > &, void *)
static constexpr std::size_t dim
This is the dimension of all vector and tensor datastructures used in MOOSE.
void fillPrincipalFromInputVector(const std::vector< T > &input)
fillPrincipalFromInputVector takes 9 inputs to fill a Rank-4 tensor C1111 = input0 C1122 = input1 C11...
RankFourTensorTempl(const SymmetricRankFourTensorTempl< T2 > &t)
The conversion operator from a SymmetricRankFourTensorTempl
static constexpr Real mandelFactor(unsigned int i, unsigned int j)
returns the 1, sqrt(2), or 2 prefactor in the Mandel notation for the indices i,j ranging from 0-5...
void mooseSetToZero< ADRankFourTensor >(ADRankFourTensor &v)
The following methods are specializations for using the libMesh::Parallel::packed_range_* routines fo...
std::basic_ostream< charT, traits > * os
void fillFromInputVector(const std::vector< T > &input, FillMethod fill_method)
fillFromInputVector takes some number of inputs to fill the Rank-4 tensor.
static constexpr unsigned int N2
void surfaceFillFromInputVector(const std::vector< T > &input)
Fills the tensor entries ignoring the last dimension (ie, C_ijkl=0 if any of i, j, k, or l = 3).
static constexpr unsigned int N4
T contractionIj(unsigned int, unsigned int, const RankTwoTensorTempl< T > &) const
Sum C_ijkl M_kl for a given i,j.
void print(std::ostream &stm=Moose::out) const
Print the rank four tensor.
RankFourTensorTempl< T > & operator/=(const T &a)
C_ijkl /= a for all i, j, k, l.
RankThreeTensor is designed to handle any N-dimensional third order tensor, r.
auto operator*(const Tensor< T2 > &a) const -> typename std::enable_if< TwoTensorMultTraits< Tensor, T2 >::value, RankTwoTensorTempl< decltype(T() *T2())>>::type
C_ijkl*a_kl.
RankFourTensorTempl< T > tripleProductJkl(const RankTwoTensorTempl< T > &, const RankTwoTensorTempl< T > &, const RankTwoTensorTempl< T > &) const
Calculates C_imnt A_jm B_kn C_lt.
RankFourTensorTempl< T > tripleProductIjk(const RankTwoTensorTempl< T > &, const RankTwoTensorTempl< T > &, const RankTwoTensorTempl< T > &) const
Calculates C_mntl A_im B_jn C_kt.
void fillSymmetricIsotropic(const T &i0, const T &i1)
auto operator/(const T2 &a) const -> typename std::enable_if< libMesh::ScalarTraits< T2 >::value, RankFourTensorTempl< decltype(T()/T2())>>::type
C_ijkl/a.
auto operator+(const RankFourTensorTempl< T2 > &a) const -> RankFourTensorTempl< decltype(T()+T2())>
C_ijkl + a_ijkl.
void fillGeneralIsotropicFromInputVector(const std::vector< T > &input)
fillGeneralIsotropicFromInputVector takes 3 inputs to fill the Rank-4 tensor with symmetries C_ijkl =...
void zero()
Zeros out the tensor.
typename tuple_n< Index, T >::template type<> tuple_of
FillMethod
To fill up the 81 entries in the 4th-order tensor, fillFromInputVector is called with one of the foll...
T L2norm() const
sqrt(C_ijkl*C_ijkl)
void mooseSetToZero< RankFourTensor >(RankFourTensor &v)
Helper function template specialization to set an object to zero.
void fillSymmetricIsotropicEandNu(const T &E, const T &nu)
bool isSymmetric() const
checks if the tensor is symmetric
void fillGeneralIsotropic(const T &i0, const T &i1, const T &i2)
Vector-less fill API functions. See docs of the corresponding ...FromInputVector methods.
RankFourTensorTempl()
Default constructor; fills to zero.
This is a "smart" enum class intended to replace many of the shortcomings in the C++ enum type It sho...
void printReal(std::ostream &stm=Moose::out) const
Print the values of the rank four tensor.
T sum3x3() const
Calculates the sum of Ciijj for i and j varying from 0 to 2.
RankFourTensorTempl< T > transposeMajor() const
Transpose the tensor by swapping the first pair with the second pair of indices.
bool isIsotropic() const
checks if the tensor is isotropic
T _vals[N4]
The values of the rank-four tensor stored by index=(((i * LIBMESH_DIM + j) * LIBMESH_DIM + k) * LIBME...
void fillAxisymmetricRZFromInputVector(const std::vector< T > &input)
fillAxisymmetricRZFromInputVector takes 5 inputs to fill the axisymmetric Rank-4 tensor with the appr...
void rotate(const TypeTensor< T > &R)
Rotate the tensor using C_ijkl = R_im R_jn R_ko R_lp C_mnop.
void fillSymmetricIsotropicEandNuFromInputVector(const std::vector< T > &input)
fillSymmetricIsotropicEandNuFromInputVector is a variation of the fillSymmetricIsotropicFromInputVect...
libMesh::VectorValue< T > sum3x1() const
Calculates the vector a[i] = sum over j Ciijj for i and j varying from 0 to 2.
T contractionKl(unsigned int, unsigned int, const RankTwoTensorTempl< T > &) const
Sum M_ij C_ijkl for a given k,l.
T & operator()(unsigned int i, unsigned int j, unsigned int k, unsigned int l)
Gets the value for the indices specified. Takes indices ranging from 0-2 for i, j, k, and l.
libMesh::boostcopy::enable_if_c< libMesh::ScalarTraits< Scalar >::value, RankFourTensorTempl & >::type operator=(const Scalar &libmesh_dbg_var(p))
Assignment-from-scalar operator.
RankFourTensorTempl< T > transposeIj() const
Transpose the tensor by swapping the first two indeces.
void fillSymmetricIsotropicFromInputVector(const std::vector< T > &input)
fillSymmetricIsotropicFromInputVector takes 2 inputs to fill the the symmetric Rank-4 tensor with the...
RankFourTensorTempl< T > singleProductI(const RankTwoTensorTempl< T > &) const
Calculates C_mjkl A_im.
RankTwoTensorTempl is designed to handle the Stress or Strain Tensor for a fully anisotropic material...
RankThreeTensorTempl< T > contraction(const libMesh::VectorValue< T > &b) const
single contraction of a RankFourTensor with a vector over index m
InitMethod
Initialization method.
SymmetricRankFourTensorTempl is designed to handle an N-dimensional fourth order tensor with minor sy...
RankFourTensorTempl< T > operator-() const
-C_ijkl
IntRange< T > make_range(T beg, T end)
void fillAntisymmetricIsotropicFromInputVector(const std::vector< T > &input)
fillAntisymmetricIsotropicFromInputVector takes 1 input to fill the the antisymmetric Rank-4 tensor w...
static MooseEnum fillMethodEnum()
Static method for use in validParams for getting the "fill_method".
RankFourTensorTempl< T > singleProductK(const RankTwoTensorTempl< T > &) const
Calculates C_ijml A_km.
static constexpr unsigned int N
tensor dimension and powers of the dimension
RankFourTensorTempl< T > & operator+=(const RankFourTensorTempl< T > &a)
C_ijkl += a_ijkl for all i, j, k, l.
const T & operator()(unsigned int i, unsigned int j, unsigned int k, unsigned int l) const
Gets the value for the indices specified.
void fillAntisymmetricIsotropic(const T &i0)
RankFourTensorTempl< T > tripleProductIkl(const RankTwoTensorTempl< T > &, const RankTwoTensorTempl< T > &, const RankTwoTensorTempl< T > &) const
Calculates C_mjnt A_im B_kn C_lt.
static RankFourTensorTempl< T > IdentityDeviatoric()
Identity of type {ik} {jl} - {ij} {kl} / 3.
RankFourTensorTempl< T > transposeKl() const
Transpose the tensor by swapping the last two indeces.
RankFourTensorTempl< T > invSymm() const
This returns A_ijkl such that C_ijkl*A_klmn = 0.5*(de_im de_jn + de_in de_jm) This routine assumes th...
RankFourTensorTempl< T > & operator-=(const RankFourTensorTempl< T > &a)
C_ijkl -= a_ijkl.