MiniTensor Version of the Day
Loading...
Searching...
No Matches
MiniTensor_Tensor4.h
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// MiniTensor Package
4//
5// Copyright 2016 NTESS and the MiniTensor contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#if !defined(MiniTensor_Tensor4_h)
11#define MiniTensor_Tensor4_h
12
13#include <iomanip>
14
15#include "MiniTensor_Tensor3.h"
16
17namespace minitensor {
18
21
25template<typename T, Index N>
26using tensor4_store = Storage<T, dimension_power<N, 4>::value>;
27
31template<typename T, Index N = DYNAMIC>
32class Tensor4 : public TensorBase<T, tensor4_store<T, N>>
33{
34public:
35
39 static constexpr
40 Index
41 ORDER = 4;
42
46 static constexpr
47 bool
48 IS_DYNAMIC = N == DYNAMIC;
49
53 using Store = tensor4_store<T, N>;
54
58 KOKKOS_INLINE_FUNCTION
59 static constexpr
60 Index
61 get_order()
62 {
63 return ORDER;
64 }
65
69 explicit
70 KOKKOS_INLINE_FUNCTION
71 Tensor4();
72
77 explicit
78 KOKKOS_INLINE_FUNCTION
79 Tensor4(Index const dimension);
80
85 explicit
86 KOKKOS_INLINE_FUNCTION
87 Tensor4(Filler const value);
88
94 explicit
95 KOKKOS_INLINE_FUNCTION
96 Tensor4(Index const dimension, Filler const value);
97
102 explicit
103 KOKKOS_INLINE_FUNCTION
104 Tensor4(T const * data_ptr);
105
111 explicit
112 KOKKOS_INLINE_FUNCTION
113 Tensor4(Index const dimension, T const * data_ptr);
114
119 KOKKOS_INLINE_FUNCTION
120 Tensor4(Tensor4<T, N> const & A);
121
125 KOKKOS_INLINE_FUNCTION
126 Tensor4(Tensor<T, dimension_square<N>::value> const & A);
127
131 virtual
132 KOKKOS_INLINE_FUNCTION
133 ~Tensor4();
134
142 KOKKOS_INLINE_FUNCTION
143 T const &
144 operator()(
145 Index const i,
146 Index const j,
147 Index const k,
148 Index const l) const;
149
157 KOKKOS_INLINE_FUNCTION
158 T &
159 operator()(
160 Index const i,
161 Index const j,
162 Index const k,
163 Index const l);
164
168 KOKKOS_INLINE_FUNCTION
169 Index
170 get_dimension() const;
171
175 KOKKOS_INLINE_FUNCTION
176 void
177 set_dimension(Index const dimension);
178
179};
180
187template<typename S, typename T, Index N>
188KOKKOS_INLINE_FUNCTION
189Tensor4<typename Promote<S, T>::type, N>
190operator+(Tensor4<S, N> const & A, Tensor4<T, N> const & B);
191
198template<typename S, typename T, Index N>
199KOKKOS_INLINE_FUNCTION
200Tensor4<typename Promote<S, T>::type, N>
201operator-(Tensor4<S, N> const & A, Tensor4<T, N> const & B);
202
207template<typename T, Index N>
208KOKKOS_INLINE_FUNCTION
209Tensor4<T, N>
210operator-(Tensor4<T, N> const & A);
211
216template<typename T, Index N>
217KOKKOS_INLINE_FUNCTION
218bool
219operator==(Tensor4<T, N> const & A, Tensor4<T, N> const & B);
220
225template<typename T, Index N>
226KOKKOS_INLINE_FUNCTION
227bool
228operator!=(Tensor4<T, N> const & A, Tensor4<T, N> const & B);
229
236template<typename S, typename T, Index N>
237KOKKOS_INLINE_FUNCTION
238typename lazy_disable_if< order_1234<S>, apply_tensor4< Promote<S,T>, N>>::type
239operator*(S const & s, Tensor4<T, N> const & A);
240
247template<typename S, typename T, Index N>
248KOKKOS_INLINE_FUNCTION
249typename lazy_disable_if< order_1234<S>, apply_tensor4< Promote<S,T>, N>>::type
250operator*(Tensor4<T, N> const & A, S const & s);
251
258template<typename S, typename T, Index N>
259KOKKOS_INLINE_FUNCTION
260Tensor4<typename Promote<S, T>::type, N>
261operator/(Tensor4<T, N> const & A, S const & s);
262
269template<typename S, typename T, Index N>
270KOKKOS_INLINE_FUNCTION
271Tensor4<typename Promote<S, T>::type, N>
272operator/(S const & s, Tensor4<T, N> const & A);
273
277template<typename T, Index N>
278KOKKOS_INLINE_FUNCTION
279Tensor4<T, N>
280transpose(Tensor4<T, N> const & A);
281
286template<typename T, Index N>
287KOKKOS_INLINE_FUNCTION
288Tensor4<T, N> const
289identity_1();
290
294template<typename T>
295KOKKOS_INLINE_FUNCTION
296Tensor4<T, DYNAMIC> const
297identity_1(Index const dimension);
298
299template<typename T, Index N>
300KOKKOS_INLINE_FUNCTION
301Tensor4<T, N> const
302identity_1(Index const dimension);
303
308template<typename T, Index N>
309KOKKOS_INLINE_FUNCTION
310Tensor4<T, N> const
311identity_2();
312
316template<typename T>
317KOKKOS_INLINE_FUNCTION
318Tensor4<T, DYNAMIC> const
319identity_2(Index const dimension);
320
321template<typename T, Index N>
322KOKKOS_INLINE_FUNCTION
323Tensor4<T, N> const
324identity_2(Index const dimension);
325
330template<typename T, Index N>
331KOKKOS_INLINE_FUNCTION
332Tensor4<T, N> const
333identity_3();
334
338template<typename T>
339KOKKOS_INLINE_FUNCTION
340Tensor4<T, DYNAMIC> const
341identity_3(Index const dimension);
342
343template<typename T, Index N>
344KOKKOS_INLINE_FUNCTION
345Tensor4<T, N> const
346identity_3(Index const dimension);
347
351template<typename T, Index N>
352KOKKOS_INLINE_FUNCTION
353Tensor4<T, N> const
354levi_civita_4();
355
356template<typename T>
357KOKKOS_INLINE_FUNCTION
358Tensor4<T, DYNAMIC> const
359levi_civita_4(Index const dimension);
360
361template<typename T, Index N>
362KOKKOS_INLINE_FUNCTION
363Tensor4<T, N> const
364levi_civita_4(Index const dimension);
365
369template<typename T, Index N>
370KOKKOS_INLINE_FUNCTION
371Tensor4<T, N> const
372permutation_4();
373
374template<typename T>
375KOKKOS_INLINE_FUNCTION
376Tensor4<T, DYNAMIC> const
377permutation_4(Index const dimension);
378
379template<typename T, Index N>
380KOKKOS_INLINE_FUNCTION
381Tensor4<T, N> const
382permutation_4(Index const dimension);
383
387template<typename T, Index N>
388KOKKOS_INLINE_FUNCTION
389Tensor4<T, N> const
390alternator_4();
391
392template<typename T>
393KOKKOS_INLINE_FUNCTION
394Tensor4<T, DYNAMIC> const
395alternator_4(Index const dimension);
396
397template<typename T, Index N>
398KOKKOS_INLINE_FUNCTION
399Tensor4<T, N> const
400alternator_4(Index const dimension);
401
406template<typename T, Index N>
407KOKKOS_INLINE_FUNCTION
408Tensor4<T, N>
409inverse(Tensor4<T, N> const & A);
410
417template<typename S, typename T, Index N>
418KOKKOS_INLINE_FUNCTION
419Tensor3<typename Promote<S, T>::type, N>
420dot(Tensor4<T, N> const & A, Vector<S, N> const & u);
421
428template<typename S, typename T, Index N>
429KOKKOS_INLINE_FUNCTION
430Tensor3<typename Promote<S, T>::type, N>
431dot(Vector<S, N> const & u, Tensor4<T, N> const & A);
432
439template<typename S, typename T, Index N>
440KOKKOS_INLINE_FUNCTION
441Tensor3<typename Promote<S, T>::type, N>
442dot2(Tensor4<T, N> const & A, Vector<S, N> const & u);
443
450template<typename S, typename T, Index N>
451KOKKOS_INLINE_FUNCTION
452Tensor3<typename Promote<S, T>::type, N>
453dot2(Vector<S, N> const & u, Tensor4<T, N> const & A);
454
461template<typename S, typename T, Index N>
462KOKKOS_INLINE_FUNCTION
463Tensor<typename Promote<S, T>::type, N>
464dotdot(Tensor4<T, N> const & A, Tensor<S, N> const & B);
465
472template<typename S, typename T, Index N>
473KOKKOS_INLINE_FUNCTION
474Tensor<typename Promote<S, T>::type, N>
475dotdot(Tensor<S, N> const & B, Tensor4<T, N> const & A);
476
483template<typename S, typename T, Index N>
484KOKKOS_INLINE_FUNCTION
485Tensor4<typename Promote<S, T>::type, N>
486dotdot(Tensor4<S, N> const & A, Tensor4<T, N> const & B);
487
494template<typename S, typename T, Index N>
495KOKKOS_INLINE_FUNCTION
496Tensor4<typename Promote<S, T>::type, N>
497tensor(Tensor<S, N> const & A, Tensor<T, N> const & B);
498
505template<typename S, typename T, Index N>
506KOKKOS_INLINE_FUNCTION
507Tensor4<typename Promote<S, T>::type, N>
508tensor2(Tensor<S, N> const & A, Tensor<T, N> const & B);
509
516template<typename S, typename T, Index N>
517KOKKOS_INLINE_FUNCTION
518Tensor4<typename Promote<S, T>::type, N>
519tensor3(Tensor<S, N> const & A, Tensor<T, N> const & B);
520
527template<typename S, typename T, Index N>
528KOKKOS_INLINE_FUNCTION
529Tensor4<typename Promote<S, T>::type, N>
530dot(Tensor4<T, N> const & A, Tensor<S, N> const & B);
531
538template<typename S, typename T, Index N>
539KOKKOS_INLINE_FUNCTION
540Tensor4<typename Promote<S, T>::type, N>
541dot_t(Tensor4<T, N> const & A, Tensor<S, N> const & B);
542
549template<typename S, typename T, Index N>
550KOKKOS_INLINE_FUNCTION
551Tensor4<typename Promote<S, T>::type, N>
552dot(Tensor<S> const & A, Tensor4<T, N> const & B);
553
560template<typename S, typename T, Index N>
561KOKKOS_INLINE_FUNCTION
562Tensor4<typename Promote<S, T>::type, N>
563t_dot(Tensor<S, N> const & A, Tensor4<T, N> const & B);
564
571template<typename S, typename T, Index N>
572KOKKOS_INLINE_FUNCTION
573Tensor4<typename Promote<S, T>::type, N>
574dot2(Tensor4<T, N> const & A, Tensor<S, N> const & B);
575
582template<typename S, typename T, Index N>
583KOKKOS_INLINE_FUNCTION
584Tensor4<typename Promote<S, T>::type, N>
585dot2_t(Tensor4<T, N> const & A, Tensor<S, N> const & B);
586
593template<typename S, typename T, Index N>
594KOKKOS_INLINE_FUNCTION
595Tensor4<typename Promote<S, T>::type, N>
596dot2(Tensor<S, N> const & A, Tensor4<T, N> const & B);
597
604template<typename S, typename T, Index N>
605KOKKOS_INLINE_FUNCTION
606Tensor4<typename Promote<S, T>::type, N>
607t_dot2(Tensor<S, N> const & A, Tensor4<T, N> const & B);
608
617template<typename S, typename T, Index N>
618KOKKOS_INLINE_FUNCTION
619Tensor4<typename Promote<S, T>::type, N>
620odot(Tensor<S, N> const & A, Tensor<T, N> const & B);
621
628template<typename S, typename T, Index N>
629KOKKOS_INLINE_FUNCTION
630Tensor4<typename Promote<S, T>::type, N>
631kronecker(Tensor<S, N> const & A, Tensor4<T, N> const & B);
632
639template<typename T, Index N>
640std::istream &
641operator>>(std::istream & is, Tensor4<T, N> & A);
642
649template<typename T, Index N>
650std::ostream &
651operator<<(std::ostream & os, Tensor4<T, N> const & A);
652
653} // namespace minitensor
654
655namespace minitensor
656{
657
658//
659// 4th-order tensor constructor with NaNs
660//
661template<typename T, Index N>
662KOKKOS_INLINE_FUNCTION
663Tensor4<T, N>::Tensor4() :
664TensorBase<T, Store>::TensorBase()
665{
666 set_dimension(N);
667 return;
668}
669
670template<typename T, Index N>
671KOKKOS_INLINE_FUNCTION
672Tensor4<T, N>::Tensor4(Index const dimension) :
673TensorBase<T, Store>::TensorBase(dimension, ORDER)
674{
675 return;
676}
677
678//
679// 4th-order tensor constructor with a specified value
680//
681template<typename T, Index N>
682KOKKOS_INLINE_FUNCTION
683Tensor4<T, N>::Tensor4(Filler const value) :
684TensorBase<T, Store>::TensorBase(N, ORDER, value)
685{
686 return;
687}
688
689template<typename T, Index N>
690KOKKOS_INLINE_FUNCTION
691Tensor4<T, N>::Tensor4(Index const dimension, Filler const value) :
692TensorBase<T, Store>::TensorBase(dimension, ORDER, value)
693{
694 return;
695}
696
697//
698// Create 4th-order tensor from array
699//
700template<typename T, Index N>
701KOKKOS_INLINE_FUNCTION
702Tensor4<T, N>::Tensor4(T const * data_ptr) :
703TensorBase<T, Store>::TensorBase(N, ORDER, data_ptr)
704{
705 return;
706}
707
708template<typename T, Index N>
709KOKKOS_INLINE_FUNCTION
710Tensor4<T, N>::Tensor4(Index const dimension, T const * data_ptr) :
711TensorBase<T, Store>::TensorBase(dimension, ORDER, data_ptr)
712{
713 return;
714}
715
716//
717// Copy constructor
718//
719template<typename T, Index N>
720KOKKOS_INLINE_FUNCTION
721Tensor4<T, N>::Tensor4(Tensor4<T, N> const & A) :
722TensorBase<T, Store>::TensorBase(A)
723{
724 return;
725}
726
727//
728// 4th-order tensor from 2nd-order tensor
729//
730
731namespace {
732
733KOKKOS_INLINE_FUNCTION
734Index
735second_to_fourth_dimension(Index const dimension_2nd)
736{
737 Index
738 dimension_4th = 0;
739
740 switch (dimension_2nd) {
741
742 default:
743 MT_ERROR_EXIT("Invalid dimension for 2nd-order tensor.");
744 break;
745
746 case 1:
747 dimension_4th = 1;
748 break;
749
750 case 4:
751 dimension_4th = 2;
752 break;
753
754 case 9:
755 dimension_4th = 3;
756 break;
757
758 case 16:
759 dimension_4th = 4;
760 break;
761
762 }
763
764 return dimension_4th;
765}
766
767} //anonymous namespace
768
769template<typename T, Index N>
770KOKKOS_INLINE_FUNCTION
771Tensor4<T, N>::Tensor4(Tensor<T, dimension_square<N>::value> const & A)
772{
773 Index const
774 dimension_2nd = A.get_dimension();
775
776 Index const
777 dimension_4th = second_to_fourth_dimension(dimension_2nd);
778
779 Tensor4<T, N> &
780 self = (*this);
781
782 self.set_dimension(dimension_4th);
783
784 Index const
785 number_components = dimension_2nd * dimension_2nd;
786
787 for (Index i = 0; i < number_components; ++i) {
788 self[i] = A[i];
789 }
790
791 return;
792}
793//
794// 4th-order tensor simple destructor
795//
796template<typename T, Index N>
797KOKKOS_INLINE_FUNCTION
798Tensor4<T, N>::~Tensor4()
799{
800 return;
801}
802
803//
804// Get dimension
805//
806template<typename T, Index N>
807KOKKOS_INLINE_FUNCTION
808Index
813
814//
815// Set dimension
816//
817template<typename T, Index N>
818KOKKOS_INLINE_FUNCTION
819void
820Tensor4<T, N>::set_dimension(Index const dimension)
821{
822 TensorBase<T, Store>::set_dimension(dimension, ORDER);
823 return;
824}
825
826//
827// 4th-order tensor addition
828//
829template<typename S, typename T, Index N>
830KOKKOS_INLINE_FUNCTION
831Tensor4<typename Promote<S, T>::type, N>
832operator+(Tensor4<S, N> const & A, Tensor4<T, N> const & B)
833{
834 Tensor4<typename Promote<S, T>::type, N>
835 C(A.get_dimension());
836
837 add(A, B, C);
838
839 return C;
840}
841
842//
843// 4th-order tensor subtraction
844//
845template<typename S, typename T, Index N>
846KOKKOS_INLINE_FUNCTION
847Tensor4<typename Promote<S, T>::type, N>
848operator-(Tensor4<S, N> const & A, Tensor4<T, N> const & B)
849{
850 Tensor4<typename Promote<S, T>::type, N>
851 C(A.get_dimension());
852
853 subtract(A, B, C);
854
855 return C;
856}
857
858//
859// 4th-order tensor minus
860//
861template<typename T, Index N>
862KOKKOS_INLINE_FUNCTION
863Tensor4<T, N>
864operator-(Tensor4<T, N> const & A)
865{
866 Tensor4<T, N>
867 B(A.get_dimension());
868
869 minus(A, B);
870
871 return B;
872}
873
874//
875// 4th-order equality
876//
877template<typename T, Index N>
878KOKKOS_INLINE_FUNCTION
879bool
880operator==(Tensor4<T, N> const & A, Tensor4<T, N> const & B)
881{
882 return equal(A, B);
883}
884
885//
886// 4th-order inequality
887//
888template<typename T, Index N>
889KOKKOS_INLINE_FUNCTION
890bool
891operator!=(Tensor4<T, N> const & A, Tensor4<T, N> const & B)
892{
893 return not_equal(A, B);
894}
895
896//
897// Scalar 4th-order tensor product
898//
899template<typename S, typename T, Index N>
900KOKKOS_INLINE_FUNCTION
901typename lazy_disable_if< order_1234<S>, apply_tensor4< Promote<S,T>, N>>::type
902operator*(S const & s, Tensor4<T, N> const & A)
903{
904 Tensor4<typename Promote<S, T>::type, N>
905 B(A.get_dimension());
906
907 scale(A, s, B);
908
909 return B;
910}
911
912//
913// 4th-order tensor scalar product
914//
915template<typename S, typename T, Index N>
916KOKKOS_INLINE_FUNCTION
917typename lazy_disable_if< order_1234<S>, apply_tensor4< Promote<S,T>, N>>::type
918operator*(Tensor4<T, N> const & A, S const & s)
919{
920 Tensor4<typename Promote<S, T>::type, N>
921 B(A.get_dimension());
922
923 scale(A, s, B);
924
925 return B;
926}
927
928//
929// 4th-order tensor scalar division
930//
931template<typename S, typename T, Index N>
932KOKKOS_INLINE_FUNCTION
933Tensor4<typename Promote<S, T>::type, N>
934operator/(Tensor4<T, N> const & A, S const & s)
935{
936 Tensor4<typename Promote<S, T>::type, N>
937 B(A.get_dimension());
938
939 divide(A, s, B);
940
941 return B;
942}
943
944//
945// 4th-order scalar tensor division
946//
947template<typename S, typename T, Index N>
948KOKKOS_INLINE_FUNCTION
949Tensor4<typename Promote<S, T>::type, N>
950operator/(S const & s, Tensor4<T, N> const & A)
951{
952 Tensor4<typename Promote<S, T>::type, N>
953 B(A.get_dimension());
954
955 split(A, s, B);
956
957 return B;
958}
959
960//
961// Indexing for constant 4th order tensor
962// \param i index
963// \param j index
964// \param k index
965// \param l index
966//
967template<typename T, Index N>
968KOKKOS_INLINE_FUNCTION
969T const &
970Tensor4<T, N>::operator()(
971 Index const i, Index const j, Index const k, Index const l) const
972{
973 Tensor4<T, N> const &
974 self = (*this);
975
976 Index const
977 dimension = self.get_dimension();
978
979 return self[((i * dimension + j) * dimension + k) * dimension + l];
980}
981
982//
983// 4th-order tensor indexing
984// \param i index
985// \param j index
986// \param k index
987// \param l index
988//
989template<typename T, Index N>
990KOKKOS_INLINE_FUNCTION
991T &
992Tensor4<T, N>::operator()(
993 Index const i, Index const j, Index const k, Index const l)
994{
995 Tensor4<T, N> &
996 self = (*this);
997
998 Index const
999 dimension = self.get_dimension();
1000
1001 return self[((i * dimension + j) * dimension + k) * dimension + l];
1002}
1003
1004//
1005// 4th-order inverse
1006//
1007template<typename T, Index N>
1008KOKKOS_INLINE_FUNCTION
1009Tensor4<T, N>
1014
1015} // namespace minitensor
1016namespace minitensor {
1017
1018namespace {
1019
1020template< typename T, Index N>
1021KOKKOS_INLINE_FUNCTION
1022void ones_in_ikjl(Tensor4<T, N> & A)
1023{
1024
1025 Index const
1026 dimension = A.get_dimension();
1027
1028 for (Index i = 0; i < dimension; ++i) {
1029 for (Index j = 0; j < dimension; ++j) {
1030 for (Index k = 0; k < dimension; ++k) {
1031 for (Index l = 0; l < dimension; ++l) {
1032 if (i == k && j == l) {
1033 A(i,j,k,l) = 1;
1034 }
1035 }
1036 }
1037 }
1038 }
1039
1040 return;
1041}
1042
1043template< typename T, Index N>
1044KOKKOS_INLINE_FUNCTION
1045void ones_in_iljk(Tensor4<T, N> & A)
1046{
1047
1048 Index const
1049 dimension = A.get_dimension();
1050
1051 for (Index i = 0; i < dimension; ++i) {
1052 for (Index j = 0; j < dimension; ++j) {
1053 for (Index k = 0; k < dimension; ++k) {
1054 for (Index l = 0; l < dimension; ++l) {
1055 if (i == l && j == k) {
1056 A(i,j,k,l) = 1;
1057 }
1058 }
1059 }
1060 }
1061 }
1062
1063 return;
1064}
1065
1066template< typename T, Index N>
1067KOKKOS_INLINE_FUNCTION
1068void ones_in_ijkl(Tensor4<T, N> & A)
1069{
1070
1071 Index const
1072 dimension = A.get_dimension();
1073
1074 for (Index i = 0; i < dimension; ++i) {
1075 for (Index j = 0; j < dimension; ++j) {
1076 for (Index k = 0; k < dimension; ++k) {
1077 for (Index l = 0; l < dimension; ++l) {
1078 if (i == j && k == l) {
1079 A(i,j,k,l) = 1;
1080 }
1081 }
1082 }
1083 }
1084 }
1085
1086 return;
1087}
1088
1089template<typename T, Index N>
1090KOKKOS_INLINE_FUNCTION
1091void fill_levi_civita(Tensor4<T, N> & A)
1092{
1093 Index const
1094 dimension = A.get_dimension();
1095
1096 for (Index i = 0; i < dimension; ++i) {
1097 for (Index j = 0; j < dimension; ++j) {
1098 for (Index k = 0; k < dimension; ++k) {
1099 for (Index l = 0; l < dimension; ++l) {
1100 A(i, j, k, l) = levi_civita<T>(i, j, k, l);
1101 }
1102 }
1103 }
1104 }
1105
1106 return;
1107}
1108
1109} // anonymous namespace
1110
1111//
1112// 4th-order identity I1
1113// \return \f$ \delta_{ik} \delta_{jl} \f$ such that \f$ A = I_1 A \f$
1114//
1115template<typename T, Index N>
1116KOKKOS_INLINE_FUNCTION
1117Tensor4<T, N> const
1118identity_1()
1119{
1120 Tensor4<T, N> I(N, Filler::ZEROS);
1121 ones_in_ikjl(I);
1122 return I;
1123}
1124
1125template<typename T>
1126KOKKOS_INLINE_FUNCTION
1127const Tensor4<T, DYNAMIC>
1128identity_1(Index const dimension)
1129{
1130 Tensor4<T, DYNAMIC> I(dimension, Filler::ZEROS);
1131 ones_in_ikjl(I);
1132 return I;
1133}
1134
1138template<typename T, Index N>
1139KOKKOS_INLINE_FUNCTION
1140Tensor4<T, N> const
1141identity_1(Index const dimension)
1142{
1143 if (N != DYNAMIC) assert(dimension == N);
1144
1145 Tensor4<T, N> I(dimension, Filler::ZEROS);
1146 ones_in_ikjl(I);
1147 return I;
1148}
1149
1150//
1151// 4th-order identity I2
1152// \return \f$ \delta_{il} \delta_{jk} \f$ such that \f$ A^T = I_2 A \f$
1153//
1154template<typename T, Index N>
1155KOKKOS_INLINE_FUNCTION
1156Tensor4<T, N> const
1157identity_2()
1158{
1159 Tensor4<T, N> I(N, Filler::ZEROS);
1160 ones_in_iljk(I);
1161 return I;
1162}
1163
1164template<typename T>
1165KOKKOS_INLINE_FUNCTION
1166const Tensor4<T, DYNAMIC>
1167identity_2(Index const dimension)
1168{
1169 Tensor4<T, DYNAMIC> I(dimension, Filler::ZEROS);
1170 ones_in_iljk(I);
1171 return I;
1172}
1173
1177template<typename T, Index N>
1178KOKKOS_INLINE_FUNCTION
1179Tensor4<T, N> const
1180identity_2(Index const dimension)
1181{
1182 if (N != DYNAMIC) assert(dimension == N);
1183
1184 Tensor4<T, N> I(dimension, Filler::ZEROS);
1185 ones_in_iljk(I);
1186 return I;
1187}
1188
1189//
1190// 4th-order identity I3
1191// \return \f$ \delta_{ij} \delta_{kl} \f$ such that \f$ I_A I = I_3 A \f$
1192//
1193template<typename T, Index N>
1194KOKKOS_INLINE_FUNCTION
1195Tensor4<T, N> const
1196identity_3()
1197{
1198 Tensor4<T, N> I(N, Filler::ZEROS);
1199 ones_in_ijkl(I);
1200 return I;
1201}
1202
1203template<typename T>
1204KOKKOS_INLINE_FUNCTION
1205const Tensor4<T, DYNAMIC>
1206identity_3(Index const dimension)
1207{
1208 Tensor4<T, DYNAMIC> I(dimension, Filler::ZEROS);
1209 ones_in_ijkl(I);
1210 return I;
1211}
1212
1216template<typename T, Index N>
1217KOKKOS_INLINE_FUNCTION
1218Tensor4<T, N> const
1219identity_3(Index const dimension)
1220{
1221 Tensor4<T, N> I(dimension, Filler::ZEROS);
1222 ones_in_ijkl(I);
1223 return I;
1224}
1225
1226//
1227// Levi-Civita symbol
1228//
1229template<typename T, Index N>
1230KOKKOS_INLINE_FUNCTION
1231Tensor4<T, N> const
1232levi_civita_4()
1233{
1234 Tensor4<T, N>
1235 A(N, Filler::ZEROS);
1236
1237 fill_levi_civita(A);
1238
1239 return A;
1240}
1241
1245template<typename T>
1246KOKKOS_INLINE_FUNCTION
1247Tensor4<T, DYNAMIC> const
1248levi_civita_4(Index const dimension)
1249{
1250 Tensor4<T, DYNAMIC>
1251 A(dimension, Filler::ZEROS);
1252
1253 fill_levi_civita(A);
1254
1255 return A;
1256}
1257
1261template<typename T, Index N>
1262KOKKOS_INLINE_FUNCTION
1263Tensor4<T, N> const
1264levi_civita_4(Index const dimension)
1265{
1266 if (N != DYNAMIC) assert(dimension == N);
1267
1268 Tensor4<T, DYNAMIC>
1269 A(dimension, Filler::ZEROS);
1270
1271 fill_levi_civita(A);
1272
1273 return A;
1274}
1275
1276//
1277// Permutation symbol
1278//
1279template<typename T, Index N>
1280KOKKOS_INLINE_FUNCTION
1281Tensor4<T, N> const
1282permutation_4()
1283{
1284 return levi_civita_4<T, N>();
1285}
1286
1290template<typename T>
1291KOKKOS_INLINE_FUNCTION
1292Tensor4<T, DYNAMIC> const
1293permutation_4(Index const dimension)
1294{
1295 return levi_civita_4<T>(dimension);
1296}
1297
1301template<typename T, Index N>
1302KOKKOS_INLINE_FUNCTION
1303Tensor4<T, N> const
1304permutation_4(Index const dimension)
1305{
1306 return levi_civita_4<T, N>(dimension);
1307}
1308
1309//
1310// Alternating symbol
1311//
1312template<typename T, Index N>
1313KOKKOS_INLINE_FUNCTION
1314Tensor4<T, N> const
1315alternator_4()
1316{
1317 return levi_civita_4<T, N>();
1318}
1319
1323template<typename T>
1324KOKKOS_INLINE_FUNCTION
1325Tensor4<T, DYNAMIC> const
1326alternator_4(Index const dimension)
1327{
1328 return levi_civita_4<T>(dimension);
1329}
1330
1331
1335template<typename T, Index N>
1336KOKKOS_INLINE_FUNCTION
1337Tensor4<T, N> const
1338alternator_4(Index const dimension)
1339{
1340 return levi_civita_4<T, N>(dimension);
1341}
1342
1343//
1344// 4th-order tensor transpose
1345// per Holzapfel 1.157
1346//
1347template<typename T, Index N>
1348KOKKOS_INLINE_FUNCTION
1349Tensor4<T, N>
1350transpose(Tensor4<T, N> const & A)
1351{
1352 Index const
1353 dimension = A.get_dimension();
1354
1355 Tensor4<T, N>
1356 B(dimension);
1357
1358 for (Index i = 0; i < dimension; ++i) {
1359 for (Index j = 0; j < dimension; ++j) {
1360 for (Index k = 0; k < dimension; ++k) {
1361 for (Index l = 0; l < dimension; ++l) {
1362 B(i, j, k, l) = A(k, l, i, j);
1363 }
1364 }
1365 }
1366 }
1367
1368 return B;
1369}
1370
1371//
1372// 4th-order tensor vector dot product
1373// \param A 4th-order tensor
1374// \param u vector
1375// \return 3rd-order tensor \f$ B = A \cdot u := B_{ijk}=A_{ijkp} u_{p} \f$
1376//
1377template<typename S, typename T, Index N>
1378KOKKOS_INLINE_FUNCTION
1379Tensor3<typename Promote<S, T>::type, N>
1380dot(Tensor4<T, N> const & A, Vector<S, N> const & u)
1381{
1382 Index const
1383 dimension= A.get_dimension();
1384
1385 assert(u.get_dimension() == dimension);
1386
1387 Tensor3<typename Promote<S, T>::type, N>
1388 B(dimension);
1389
1390 for (Index i = 0; i < dimension; ++i) {
1391 for (Index j = 0; j < dimension; ++j) {
1392 for (Index k = 0; k < dimension; ++k) {
1393
1394 typename Promote<S, T>::type
1395 s = 0.0;
1396
1397 for (Index p = 0; p < dimension; ++p) {
1398 s += A(i,j,k,p) * u(p);
1399 }
1400 B(i,j,k) = s;
1401 }
1402 }
1403 }
1404
1405 return B;
1406}
1407
1408//
1409// vector 4th-order tensor dot product
1410// \param A 4th-order tensor
1411// \param u vector
1412// \return 3rd-order tensor \f$ u dot A \f$ as \f$ B_{ijk}=u_{p} A_{pijk} \f$
1413//
1414template<typename S, typename T, Index N>
1415KOKKOS_INLINE_FUNCTION
1416Tensor3<typename Promote<S, T>::type, N>
1417dot(Vector<S, N> const & u, Tensor4<T, N> const & A)
1418{
1419 Index const
1420 dimension = A.get_dimension();
1421
1422 assert(u.get_dimension() == dimension);
1423
1424 Tensor3<typename Promote<S, T>::type, N>
1425 B(dimension);
1426
1427 for (Index i = 0; i < dimension; ++i) {
1428 for (Index j = 0; j < dimension; ++j) {
1429 for (Index k = 0; k < dimension; ++k) {
1430
1431 typename Promote<S, T>::type
1432 s = 0.0;
1433
1434 for (Index p = 0; p < dimension; ++p) {
1435 s += u(p) * A(p,i,j,k);
1436 }
1437 B(i,j,k) = s;
1438 }
1439 }
1440 }
1441
1442 return B;
1443}
1444
1445//
1446// 4th-order tensor vector dot2 product
1447// \param A 4th-order tensor
1448// \param u vector
1449// \return 3rd-order tensor \f$ B = A \cdot u := B_{ijk} = A_{ijpk} u_{p} \f$
1450//
1451template<typename S, typename T, Index N>
1452KOKKOS_INLINE_FUNCTION
1453Tensor3<typename Promote<S, T>::type, N>
1454dot2(Tensor4<T, N> const & A, Vector<S, N> const & u)
1455{
1456 Index const
1457 dimension = A.get_dimension();
1458
1459 assert(u.get_dimension() == dimension);
1460
1461 Tensor3<typename Promote<S, T>::type, N>
1462 B(dimension);
1463
1464 for (Index i = 0; i < dimension; ++i) {
1465 for (Index j = 0; j < dimension; ++j) {
1466 for (Index k = 0; k < dimension; ++k) {
1467
1468 typename Promote<S, T>::type
1469 s = 0.0;
1470
1471 for (Index p = 0; p < dimension; ++p) {
1472 s += A(i,j,p,k) * u(p);
1473 }
1474 B(i,j,k) = s;
1475 }
1476 }
1477 }
1478
1479 return B;
1480}
1481
1482//
1483// vector 4th-order tensor dot2 product
1484// \param A 4th-order tensor
1485// \param u vector
1486// \return 3rd-order tensor \f$ u dot2 A \f$ as \f$ B_{ijk}=u_{p} A_{ipjk} \f$
1487//
1488template<typename S, typename T, Index N>
1489KOKKOS_INLINE_FUNCTION
1490Tensor3<typename Promote<S, T>::type, N>
1491dot2(Vector<S, N> const & u, Tensor4<T, N> const & A)
1492{
1493 Index const
1494 dimension = A.get_dimension();
1495
1496 assert(u.get_dimension() == dimension);
1497
1498 Tensor3<typename Promote<S, T>::type, N>
1499 B(dimension);
1500
1501 for (Index i = 0; i < dimension; ++i) {
1502 for (Index j = 0; j < dimension; ++j) {
1503 for (Index k = 0; k < dimension; ++k) {
1504
1505 typename Promote<S, T>::type
1506 s = 0.0;
1507
1508 for (Index p = 0; p < dimension; ++p) {
1509 s += u(p) * A(i,p,j,k);
1510 }
1511 B(i,j,k) = s;
1512 }
1513 }
1514 }
1515
1516 return B;
1517}
1518
1519//
1520// \return 2nd-order tensor \f$ C = A : B := C_{ij} = A_{ijpq} B_{pq} \f$
1521//
1522template<typename S, typename T, Index N>
1523KOKKOS_INLINE_FUNCTION
1524Tensor<typename Promote<S, T>::type, N>
1525dotdot(Tensor4<T, N> const & A, Tensor<S, N> const & B)
1526{
1527 Index const
1528 dimension = A.get_dimension();
1529
1530 assert(B.get_dimension() == dimension);
1531
1532 Tensor<typename Promote<S, T>::type, N>
1533 C(dimension);
1534
1535 for (Index i = 0; i < dimension; ++i) {
1536 for (Index j = 0; j < dimension; ++j) {
1537
1538 typename Promote<S, T>::type
1539 s = 0.0;
1540
1541 for (Index p = 0; p < dimension; ++p) {
1542 for (Index q = 0; q < dimension; ++q) {
1543 s += A(i,j,p,q) * B(p,q);
1544 }
1545 }
1546 C(i,j) = s;
1547 }
1548 }
1549
1550 return C;
1551}
1552
1553//
1554// \return 2nd-order tensor \f$ C = B : A := C_{ij} = B_{pq} A_{pqij} \f$
1555//
1556template<typename S, typename T, Index N>
1557KOKKOS_INLINE_FUNCTION
1558Tensor<typename Promote<S, T>::type, N>
1559dotdot(Tensor<S, N> const & B, Tensor4<T, N> const & A)
1560{
1561 Index const
1562 dimension = A.get_dimension();
1563
1564 assert(B.get_dimension() == dimension);
1565
1566 Tensor<typename Promote<S, T>::type, N>
1567 C(dimension);
1568
1569 for (Index i = 0; i < dimension; ++i) {
1570 for (Index j = 0; j < dimension; ++j) {
1571
1572 typename Promote<S, T>::type
1573 s = 0.0;
1574
1575 for (Index p = 0; p < dimension; ++p) {
1576 for (Index q = 0; q < dimension; ++q) {
1577 s += B(p,q) * A(p,q,i,j);
1578 }
1579 }
1580 C(i,j) = s;
1581 }
1582 }
1583
1584 return C;
1585}
1586
1587//
1588// \return \f$ C = A : B := C_{ijkl} = A_{ijpq} B{pqkl} \f$
1589//
1590template<typename S, typename T, Index N>
1591KOKKOS_INLINE_FUNCTION
1592Tensor4<typename Promote<S, T>::type, N>
1593dotdot(Tensor4<S, N> const & A, Tensor4<T, N> const & B)
1594{
1595 Index const
1596 dimension = A.get_dimension();
1597
1598 assert(B.get_dimension() == dimension);
1599
1600 Tensor4<typename Promote<S, T>::type, N>
1601 C(dimension);
1602
1603 for (Index i = 0; i < dimension; ++i) {
1604 for (Index j = 0; j < dimension; ++j) {
1605 for (Index k = 0; k < dimension; ++k) {
1606 for (Index l = 0; l < dimension; ++l) {
1607
1608 typename Promote<S, T>::type
1609 s = 0.0;
1610
1611 for (Index p = 0; p < dimension; ++p) {
1612 for (Index q = 0; q < dimension; ++q) {
1613 s += A(i,j,p,q) * B(p,q,k,l);
1614 }
1615 }
1616 C(i,j,k,l) = s;
1617 }
1618 }
1619 }
1620 }
1621
1622 return C;
1623}
1624
1625//
1626// \return \f$ C = A \otimes B := C_{ijkl} = A_{ij} B_{kl} \f$
1627//
1628template<typename S, typename T, Index N>
1629KOKKOS_INLINE_FUNCTION
1630Tensor4<typename Promote<S, T>::type, N>
1631tensor(Tensor<S, N> const & A, Tensor<T, N> const & B)
1632{
1633 Index const
1634 dimension = A.get_dimension();
1635
1636 assert(B.get_dimension() == dimension);
1637
1638 Tensor4<typename Promote<S, T>::type, N>
1639 C(dimension);
1640
1641 for (Index i = 0; i < dimension; ++i) {
1642 for (Index j = 0; j < dimension; ++j) {
1643 for (Index k = 0; k < dimension; ++k) {
1644 for (Index l = 0; l < dimension; ++l) {
1645 C(i,j,k,l) = A(i,j) * B(k,l);
1646 }
1647 }
1648 }
1649 }
1650
1651 return C;
1652}
1653
1654//
1655// \return \f$ C_{ijkl} = A_{ik} B_{jl} \f$
1656//
1657template<typename S, typename T, Index N>
1658KOKKOS_INLINE_FUNCTION
1659Tensor4<typename Promote<S, T>::type, N>
1660tensor2(Tensor<S, N> const & A, Tensor<T, N> const & B)
1661{
1662 Index const
1663 dimension = A.get_dimension();
1664
1665 assert(B.get_dimension() == dimension);
1666
1667 Tensor4<typename Promote<S, T>::type, N>
1668 C(dimension);
1669
1670 for (Index i = 0; i < dimension; ++i) {
1671 for (Index j = 0; j < dimension; ++j) {
1672 for (Index k = 0; k < dimension; ++k) {
1673 for (Index l = 0; l < dimension; ++l) {
1674 C(i,j,k,l) = A(i,k) * B(j,l);
1675 }
1676 }
1677 }
1678 }
1679
1680 return C;
1681}
1682
1683//
1684// \return \f$ C_{ijkl} = A_{il} B_{kj} \f$
1685//
1686template<typename S, typename T, Index N>
1687KOKKOS_INLINE_FUNCTION
1688Tensor4<typename Promote<S, T>::type, N>
1689tensor3(Tensor<S, N> const & A, Tensor<T, N> const & B)
1690{
1691 Index const
1692 dimension = A.get_dimension();
1693
1694 assert(B.get_dimension() == dimension);
1695
1696 Tensor4<typename Promote<S, T>::type, N>
1697 C(dimension);
1698
1699 for (Index i = 0; i < dimension; ++i) {
1700 for (Index j = 0; j < dimension; ++j) {
1701 for (Index k = 0; k < dimension; ++k) {
1702 for (Index l = 0; l < dimension; ++l) {
1703 C(i,j,k,l) = A(i,l) * B(k,j);
1704 }
1705 }
1706 }
1707 }
1708
1709 return C;
1710}
1711
1712//
1713// \return \f$ C = A \cdot B := C_{ijkl} = A_{ijkp} B_{pl} \f$
1714//
1715template<typename S, typename T, Index N>
1716KOKKOS_INLINE_FUNCTION
1717Tensor4<typename Promote<S, T>::type, N>
1718dot(Tensor4<T, N> const & A, Tensor<S, N> const & B)
1719{
1720 Index const
1721 dimension = A.get_dimension();
1722
1723 assert(B.get_dimension() == dimension);
1724
1725 Tensor4<typename Promote<S, T>::type, N>
1726 C(dimension);
1727
1728 for (Index i = 0; i < dimension; ++i) {
1729 for (Index j = 0; j < dimension; ++j) {
1730 for (Index k = 0; k < dimension; ++k) {
1731 for (Index l = 0; l < dimension; ++l) {
1732
1733 typename Promote<S, T>::type
1734 s = 0.0;
1735
1736 for (Index p = 0; p < dimension; ++p) {
1737 s += A(i,j,k,p) * B(p,l);
1738 }
1739 C(i,j,k,l) = s;
1740 }
1741 }
1742 }
1743 }
1744
1745 return C;
1746}
1747
1748//
1749// \return \f$ C = A \cdot B^T := C_{ijkl} = A_{ijkp} B_{lp} \f$
1750//
1751template<typename S, typename T, Index N>
1752KOKKOS_INLINE_FUNCTION
1753Tensor4<typename Promote<S, T>::type, N>
1754dot_t(Tensor4<T, N> const & A, Tensor<S, N> const & B)
1755{
1756 Index const
1757 dimension = A.get_dimension();
1758
1759 assert(B.get_dimension() == dimension);
1760
1761 Tensor4<typename Promote<S, T>::type, N>
1762 C(dimension);
1763
1764 for (Index i = 0; i < dimension; ++i) {
1765 for (Index j = 0; j < dimension; ++j) {
1766 for (Index k = 0; k < dimension; ++k) {
1767 for (Index l = 0; l < dimension; ++l) {
1768
1769 typename Promote<S, T>::type
1770 s = 0.0;
1771
1772 for (Index p = 0; p < dimension; ++p) {
1773 s += A(i,j,k,p) * B(l,p);
1774 }
1775 C(i,j,k,l) = s;
1776 }
1777 }
1778 }
1779 }
1780
1781 return C;
1782}
1783
1784//
1785// \return \f$ C = A \cdot B := C_{ijkl} = A_{ip} B_{pjkl} \f$
1786//
1791template<typename S, typename T, Index N>
1792KOKKOS_INLINE_FUNCTION
1793Tensor4<typename Promote<S, T>::type, N>
1794dot(Tensor<S, N> const & A, Tensor4<T, N> const & B)
1795{
1796 Index const
1797 dimension = A.get_dimension();
1798
1799 assert(B.get_dimension() == dimension);
1800
1801 Tensor4<typename Promote<S, T>::type, N>
1802 C(dimension);
1803
1804 for (Index i = 0; i < dimension; ++i) {
1805 for (Index j = 0; j < dimension; ++j) {
1806 for (Index k = 0; k < dimension; ++k) {
1807 for (Index l = 0; l < dimension; ++l) {
1808
1809 typename Promote<S, T>::type
1810 s = 0.0;
1811
1812 for (Index p = 0; p < dimension; ++p) {
1813 s += A(i,p) * B(p,j,k,l);
1814 }
1815 C(i,j,k,l) = s;
1816 }
1817 }
1818 }
1819 }
1820
1821 return C;
1822}
1823
1824//
1825// \return \f$ C = A^T \cdot B := C_{ijkl} = A_{pi} B_{pjkl} \f$
1826//
1827template<typename S, typename T, Index N>
1828KOKKOS_INLINE_FUNCTION
1829Tensor4<typename Promote<S, T>::type, N>
1830t_dot(Tensor<S, N> const & A, Tensor4<T, N> const & B)
1831{
1832 Index const
1833 dimension = A.get_dimension();
1834
1835 assert(B.get_dimension() == dimension);
1836
1837 Tensor4<typename Promote<S, T>::type, N>
1838 C(dimension);
1839
1840 for (Index i = 0; i < dimension; ++i) {
1841 for (Index j = 0; j < dimension; ++j) {
1842 for (Index k = 0; k < dimension; ++k) {
1843 for (Index l = 0; l < dimension; ++l) {
1844
1845 typename Promote<S, T>::type
1846 s = 0.0;
1847
1848 for (Index p = 0; p < dimension; ++p) {
1849 s += A(p,i) * B(p,j,k,l);
1850 }
1851 C(i,j,k,l) = s;
1852 }
1853 }
1854 }
1855 }
1856
1857 return C;
1858}
1859
1860//
1861// \return \f$ C = A \cdot B := C_{ijkl} = A_{ijpl} B_{pk} \f$
1862//
1863template<typename S, typename T, Index N>
1864KOKKOS_INLINE_FUNCTION
1865Tensor4<typename Promote<S, T>::type, N>
1866dot2(Tensor4<T, N> const & A, Tensor<S, N> const & B)
1867{
1868 Index const
1869 dimension = A.get_dimension();
1870
1871 assert(B.get_dimension() == dimension);
1872
1873 Tensor4<typename Promote<S, T>::type, N>
1874 C(dimension);
1875
1876 for (Index i = 0; i < dimension; ++i) {
1877 for (Index j = 0; j < dimension; ++j) {
1878 for (Index k = 0; k < dimension; ++k) {
1879 for (Index l = 0; l < dimension; ++l) {
1880
1881 typename Promote<S, T>::type
1882 s = 0.0;
1883
1884 for (Index p = 0; p < dimension; ++p) {
1885 s += A(i,j,p,l) * B(p,k);
1886 }
1887 C(i,j,k,l) = s;
1888 }
1889 }
1890 }
1891 }
1892
1893 return C;
1894}
1895
1896//
1897// \return \f$ C = A \cdot B^T := C_{ijkl} = A_{ijpl} B_{kp} \f$
1898//
1899template<typename S, typename T, Index N>
1900KOKKOS_INLINE_FUNCTION
1901Tensor4<typename Promote<S, T>::type, N>
1902dot2_t(Tensor4<T, N> const & A, Tensor<S, N> const & B)
1903{
1904 Index const
1905 dimension = A.get_dimension();
1906
1907 assert(B.get_dimension() == dimension);
1908
1909 Tensor4<typename Promote<S, T>::type, N>
1910 C(dimension);
1911
1912 for (Index i = 0; i < dimension; ++i) {
1913 for (Index j = 0; j < dimension; ++j) {
1914 for (Index k = 0; k < dimension; ++k) {
1915 for (Index l = 0; l < dimension; ++l) {
1916
1917 typename Promote<S, T>::type
1918 s = 0.0;
1919
1920 for (Index p = 0; p < dimension; ++p) {
1921 s += A(i,j,p,l) * B(k,p);
1922 }
1923 C(i,j,k,l) = s;
1924 }
1925 }
1926 }
1927 }
1928
1929 return C;
1930}
1931
1932//
1933// \return \f$ C = A \cdot B := C_{ijkl} = A_{jp} B_{ipkl} \f$
1934//
1935template<typename S, typename T, Index N>
1936KOKKOS_INLINE_FUNCTION
1937Tensor4<typename Promote<S, T>::type, N>
1938dot2(Tensor<S, N> const & A, Tensor4<T, N> const & B)
1939{
1940 Index const
1941 dimension = A.get_dimension();
1942
1943 assert(B.get_dimension() == dimension);
1944
1945 Tensor4<typename Promote<S, T>::type, N>
1946 C(dimension);
1947
1948 for (Index i = 0; i < dimension; ++i) {
1949 for (Index j = 0; j < dimension; ++j) {
1950 for (Index k = 0; k < dimension; ++k) {
1951 for (Index l = 0; l < dimension; ++l) {
1952
1953 typename Promote<S, T>::type
1954 s = 0.0;
1955
1956 for (Index p = 0; p < dimension; ++p) {
1957 s += A(j,p) * B(i,p,k,l);
1958 }
1959 C(i,j,k,l) = s;
1960 }
1961 }
1962 }
1963 }
1964
1965 return C;
1966}
1967
1968//
1969// \return \f$ C = A^T \cdot B := C_{ijkl} = A_{pj} B_{ipkl} \f$
1970//
1971template<typename S, typename T, Index N>
1972KOKKOS_INLINE_FUNCTION
1973Tensor4<typename Promote<S, T>::type, N>
1974t_dot2(Tensor<S, N> const & A, Tensor4<T, N> const & B)
1975{
1976 Index const
1977 dimension = A.get_dimension();
1978
1979 assert(B.get_dimension() == dimension);
1980
1981 Tensor4<typename Promote<S, T>::type, N>
1982 C(dimension);
1983
1984 for (Index i = 0; i < dimension; ++i) {
1985 for (Index j = 0; j < dimension; ++j) {
1986 for (Index k = 0; k < dimension; ++k) {
1987 for (Index l = 0; l < dimension; ++l) {
1988
1989 typename Promote<S, T>::type
1990 s = 0.0;
1991
1992 for (Index p = 0; p < dimension; ++p) {
1993 s += A(p,j) * B(i,p,k,l);
1994 }
1995 C(i,j,k,l) = s;
1996 }
1997 }
1998 }
1999 }
2000
2001 return C;
2002}
2003
2004//
2005// odot operator useful for \f$ \frac{\partial A^{-1}}{\partial A} \f$
2006// see Holzapfel eqn 6.165
2007// \param A 2nd-order tensor
2008// \param B 2nd-order tensor
2009// \return \f$ A \odot B \f$ which is
2010// \f$ C_{ijkl} = \frac{1}{2}(A_{ik} B_{jl} + A_{il} B_{jk}) \f$
2011//
2012template<typename S, typename T, Index N>
2013KOKKOS_INLINE_FUNCTION
2014Tensor4<typename Promote<S, T>::type, N>
2015odot(Tensor<S, N> const & A, Tensor<T, N> const & B)
2016{
2017 Index const
2018 dimension = A.get_dimension();
2019
2020 assert(B.get_dimension() == dimension);
2021
2022 Tensor4<typename Promote<S, T>::type, N>
2023 C(dimension);
2024
2025 for (Index i = 0; i < dimension; ++i) {
2026 for (Index j = 0; j < dimension; ++j) {
2027 for (Index k = 0; k < dimension; ++k) {
2028 for (Index l = 0; l < dimension; ++l) {
2029 C(i,j,k,l) = 0.5 * (A(i,k) * B(j,l) + A(i,l) * B(j,k));
2030 }
2031 }
2032 }
2033 }
2034
2035 return C;
2036}
2037
2038//
2039// \return \f$ C'_{i'j'k'l'} = Q_{i'i} Q_{j'j} Q_{k'k} Q_{l'l} C_{ijkl} \f$
2040//
2041template<typename S, typename T, Index N>
2042KOKKOS_INLINE_FUNCTION
2043Tensor4<typename Promote<S, T>::type, N>
2044kronecker(Tensor<S, N> const & A, Tensor4<T, N> const & B)
2045{
2046 Index const
2047 dimension = A.get_dimension();
2048
2049 assert(B.get_dimension() == dimension);
2050
2051 Tensor4<typename Promote<S, T>::type, N>
2052 C(dimension);
2053
2054 for (Index i = 0; i < dimension; ++i) {
2055 for (Index j = 0; j < dimension; ++j) {
2056 for (Index k = 0; k < dimension; ++k) {
2057 for (Index l = 0; l < dimension; ++l) {
2058
2059 typename Promote<S, T>::type
2060 s = 0.0;
2061
2062 // we assume A is the direction cosine matrix
2063 // have s_{i}\be_{i} = s_{i'}\be_{i'}, want s_{i}
2064 // s_{i} = A_{i',i}s_{i'}
2065 for (Index p = 0; p < dimension; ++p) {
2066 for (Index q = 0; q < dimension; ++q) {
2067 for (Index m = 0; m < dimension; ++m) {
2068 for (Index n = 0; n < dimension; ++n) {
2069 s += A(p,i) * A(q,j) * A(m,k) * A(n,l) * B(p,q,m,n);
2070 }
2071 }
2072 }
2073 }
2074 C(i,j,k,l) = s;
2075 }
2076 }
2077 }
2078 }
2079
2080 return C;
2081}
2082
2083//
2084// 4th-order input
2085// \param A 4th-order tensor
2086// \param is input stream
2087// \return is input stream
2088//
2089template<typename T, Index N>
2090std::istream &
2091operator>>(std::istream & is, Tensor4<T, N> & A)
2092{
2093 Index const
2094 dimension = A.get_dimension();
2095
2096 for (Index i = 0; i < dimension; ++i) {
2097 for (Index j = 0; j < dimension; ++j) {
2098 for (Index k = 0; k < dimension; ++k) {
2099 for (Index l = 0; l < dimension; ++l) {
2100 is >> A(i,j,k,l);
2101 }
2102 }
2103 }
2104 }
2105
2106 return is;
2107}
2108
2109//
2110// 4th-order output
2111// \param A 4th-order tensor
2112// \param os output stream
2113// \return os output stream
2114//
2115template<typename T, Index N>
2116std::ostream &
2117operator<<(std::ostream & os, Tensor4<T, N> const & A)
2118{
2119 Index const
2120 dimension = A.get_dimension();
2121
2122 if (dimension == 0) {
2123 return os;
2124 }
2125
2126 os << std::scientific << std::setprecision(17);
2127
2128 for (Index i = 0; i < dimension; ++i) {
2129
2130 for (Index j = 0; j < dimension; ++j) {
2131
2132 for (Index k = 0; k < dimension; ++k) {
2133
2134 os << std::setw(24) << A(i,j,k,0);
2135
2136 for (Index l = 1; l < dimension; ++l) {
2137
2138 os << "," << std::setw(24) << A(i,j,k,l);
2139 }
2140
2141 os << std::endl;
2142
2143 }
2144
2145 os << std::endl;
2146 os << std::endl;
2147
2148 }
2149
2150 os << std::endl;
2151
2152 }
2153
2154 return os;
2155}
2156
2158} // namespace minitensor
2159
2160#endif //MiniTensor_Tensor4_h
MiniTensor_Tensor3.h
KOKKOS_INLINE_FUNCTION
#define KOKKOS_INLINE_FUNCTION
Definition MiniTensor_Traits.h:20
MT_ERROR_EXIT
#define MT_ERROR_EXIT(...)
Definition MiniTensor_Traits.h:27
minitensor::operator+
KOKKOS_INLINE_FUNCTION Matrix< typename Promote< S, T >::type, M, N > operator+(Matrix< S, M, N > const &A, Matrix< T, M, N > const &B)
Definition MiniTensor_Matrix.h:1250
minitensor::kronecker
KOKKOS_INLINE_FUNCTION Tensor4< typename Promote< S, T >::type, N > kronecker(Tensor< S, N > const &A, Tensor4< T, N > const &B)
Definition MiniTensor_Tensor4.h:2044
minitensor::minus
KOKKOS_INLINE_FUNCTION void minus(TensorBase< T, ST > const &A, TensorBase< T, ST > &B)
Definition MiniTensor_TensorBase.h:1740
minitensor::operator-
KOKKOS_INLINE_FUNCTION Matrix< typename Promote< S, T >::type, M, N > operator-(Matrix< S, M, N > const &A, Matrix< T, M, N > const &B)
Definition MiniTensor_Matrix.h:1266
minitensor::levi_civita_4
KOKKOS_INLINE_FUNCTION Tensor4< T, N > const levi_civita_4()
Definition MiniTensor_Tensor4.h:1232
minitensor::operator<<
std::ostream & operator<<(std::ostream &os, Matrix< T, M, N > const &A)
Definition MiniTensor_Matrix.h:1841
minitensor::add
KOKKOS_INLINE_FUNCTION void add(TensorBase< R, SR > const &A, TensorBase< S, SS > const &B, TensorBase< T, ST > &C)
Definition MiniTensor_TensorBase.h:1690
minitensor::Filler
Filler
Definition MiniTensor_TensorBase.h:29
minitensor::tensor2
KOKKOS_INLINE_FUNCTION Tensor4< typename Promote< S, T >::type, N > tensor2(Tensor< S, N > const &A, Tensor< T, N > const &B)
Definition MiniTensor_Tensor4.h:1660
minitensor::Tensor4::Tensor4
KOKKOS_INLINE_FUNCTION Tensor4(Index const dimension, T const *data_ptr)
Definition MiniTensor_Tensor4.h:710
minitensor::Tensor4::operator()
KOKKOS_INLINE_FUNCTION T const & operator()(Index const i, Index const j, Index const k, Index const l) const
Definition MiniTensor_Tensor4.h:970
minitensor::operator>>
std::istream & operator>>(std::istream &is, Matrix< T, M, N > &A)
Definition MiniTensor_Matrix.h:1819
minitensor::Tensor4::~Tensor4
virtual KOKKOS_INLINE_FUNCTION ~Tensor4()
Definition MiniTensor_Tensor4.h:798
minitensor::Tensor4::Tensor4
KOKKOS_INLINE_FUNCTION Tensor4(T const *data_ptr)
Definition MiniTensor_Tensor4.h:702
minitensor::Tensor4::Tensor4
KOKKOS_INLINE_FUNCTION Tensor4(Tensor4< T, N > const &A)
Definition MiniTensor_Tensor4.h:721
minitensor::Tensor4::Tensor4
KOKKOS_INLINE_FUNCTION Tensor4(Tensor< T, dimension_square< N >::value > const &A)
Definition MiniTensor_Tensor4.h:771
minitensor::Tensor4::Tensor4
KOKKOS_INLINE_FUNCTION Tensor4(Index const dimension, Filler const value)
Definition MiniTensor_Tensor4.h:691
minitensor::odot
KOKKOS_INLINE_FUNCTION Tensor4< typename Promote< S, T >::type, N > odot(Tensor< S, N > const &A, Tensor< T, N > const &B)
Definition MiniTensor_Tensor4.h:2015
minitensor::operator/
KOKKOS_INLINE_FUNCTION Matrix< typename Promote< S, T >::type, M, N > operator/(Matrix< T, M, N > const &A, S const &s)
Definition MiniTensor_Matrix.h:1354
minitensor::operator*
KOKKOS_INLINE_FUNCTION Vector< typename Promote< S, T >::type, M > operator*(Matrix< T, M, N > const &A, Vector< S, N > const &u)
Definition MiniTensor_Matrix.h:1386
minitensor::dotdot
KOKKOS_INLINE_FUNCTION Promote< S, T >::type dotdot(Matrix< S, M, N > const &A, Matrix< T, M, N > const &B)
Definition MiniTensor_Matrix.h:1732
minitensor::Vector::get_dimension
KOKKOS_INLINE_FUNCTION Index get_dimension() const
Definition MiniTensor_Vector.h:995
minitensor::tensor3
KOKKOS_INLINE_FUNCTION Tensor4< typename Promote< S, T >::type, N > tensor3(Tensor< S, N > const &A, Tensor< T, N > const &B)
Definition MiniTensor_Tensor4.h:1689
minitensor::equal
KOKKOS_INLINE_FUNCTION bool equal(TensorBase< T, ST > const &A, TensorBase< T, ST > const &B)
Definition MiniTensor_TensorBase.h:1760
minitensor::Tensor4::get_dimension
KOKKOS_INLINE_FUNCTION Index get_dimension() const
Definition MiniTensor_Tensor4.h:809
minitensor::Tensor::get_dimension
KOKKOS_INLINE_FUNCTION Index get_dimension() const
Definition MiniTensor_Tensor.h:1390
minitensor::TensorBase::get_dimension
KOKKOS_INLINE_FUNCTION Index get_dimension(Index const order) const
Definition MiniTensor_TensorBase.h:806
minitensor::Tensor4::set_dimension
KOKKOS_INLINE_FUNCTION void set_dimension(Index const dimension)
Definition MiniTensor_Tensor4.h:820
minitensor::dot_t
KOKKOS_INLINE_FUNCTION Matrix< typename Promote< S, T >::type, M, N > dot_t(Matrix< S, M, P > const &A, Matrix< T, N, P > const &B)
Definition MiniTensor_Matrix.h:1653
minitensor::scale
KOKKOS_INLINE_FUNCTION void scale(TensorBase< R, SR > const &A, S const &s, TensorBase< T, ST > &B)
Definition MiniTensor_TensorBase.h:1791
minitensor::operator!=
KOKKOS_INLINE_FUNCTION bool operator!=(Matrix< T, M, N > const &A, Matrix< T, M, N > const &B)
Definition MiniTensor_Matrix.h:1309
minitensor::Tensor4::Tensor4
KOKKOS_INLINE_FUNCTION Tensor4()
Definition MiniTensor_Tensor4.h:663
minitensor::permutation_4
KOKKOS_INLINE_FUNCTION Tensor4< T, N > const permutation_4()
Definition MiniTensor_Tensor4.h:1282
minitensor::dot2_t
KOKKOS_INLINE_FUNCTION Tensor4< typename Promote< S, T >::type, N > dot2_t(Tensor4< T, N > const &A, Tensor< S, N > const &B)
Definition MiniTensor_Tensor4.h:1902
minitensor::operator==
KOKKOS_INLINE_FUNCTION bool operator==(Matrix< T, M, N > const &A, Matrix< T, M, N > const &B)
Definition MiniTensor_Matrix.h:1298
minitensor::dot
KOKKOS_INLINE_FUNCTION Vector< typename Promote< S, T >::type, M > dot(Matrix< T, M, N > const &A, Vector< S, N > const &u)
Definition MiniTensor_Matrix.h:1441
minitensor::Tensor4::get_order
static KOKKOS_INLINE_FUNCTION constexpr Index get_order()
Definition MiniTensor_Tensor4.h:61
minitensor::alternator_4
KOKKOS_INLINE_FUNCTION Tensor4< T, N > const alternator_4()
Definition MiniTensor_Tensor4.h:1315
minitensor::split
KOKKOS_INLINE_FUNCTION void split(TensorBase< R, SR > const &A, S const &s, TensorBase< T, ST > &B)
Definition MiniTensor_TensorBase.h:1831
minitensor::Tensor4::ORDER
static constexpr Index ORDER
Definition MiniTensor_Tensor4.h:41
minitensor::Tensor4::IS_DYNAMIC
static constexpr bool IS_DYNAMIC
Definition MiniTensor_Tensor4.h:48
minitensor::t_dot
KOKKOS_INLINE_FUNCTION Matrix< typename Promote< S, T >::type, M, N > t_dot(Matrix< S, P, M > const &A, Matrix< T, P, N > const &B)
Definition MiniTensor_Matrix.h:1615
minitensor::tensor
KOKKOS_INLINE_FUNCTION Tensor< typename Promote< S, T >::type, N > tensor(Vector< S, N > const &u, Vector< T, N > const &v)
Definition MiniTensor_Tensor.h:2212
minitensor::identity_2
KOKKOS_INLINE_FUNCTION Tensor4< T, N > const identity_2()
Definition MiniTensor_Tensor4.h:1157
minitensor::identity_3
KOKKOS_INLINE_FUNCTION Tensor4< T, N > const identity_3()
Definition MiniTensor_Tensor4.h:1196
minitensor::not_equal
KOKKOS_INLINE_FUNCTION bool not_equal(TensorBase< T, ST > const &A, TensorBase< T, ST > const &B)
Definition MiniTensor_TensorBase.h:1780
minitensor::Tensor4::Tensor4
KOKKOS_INLINE_FUNCTION Tensor4(Filler const value)
Definition MiniTensor_Tensor4.h:683
minitensor::identity_1
KOKKOS_INLINE_FUNCTION Tensor4< T, N > const identity_1()
Definition MiniTensor_Tensor4.h:1118
minitensor::TensorBase::set_dimension
KOKKOS_INLINE_FUNCTION void set_dimension(Index const dimension, Index const order)
Definition MiniTensor_TensorBase.h:817
minitensor::Tensor4::operator()
KOKKOS_INLINE_FUNCTION T & operator()(Index const i, Index const j, Index const k, Index const l)
Definition MiniTensor_Tensor4.h:992
minitensor::Tensor4::Tensor4
KOKKOS_INLINE_FUNCTION Tensor4(Index const dimension)
Definition MiniTensor_Tensor4.h:672
minitensor::t_dot2
KOKKOS_INLINE_FUNCTION Tensor4< typename Promote< S, T >::type, N > t_dot2(Tensor< S, N > const &A, Tensor4< T, N > const &B)
Definition MiniTensor_Tensor4.h:1974
minitensor::transpose
KOKKOS_INLINE_FUNCTION Matrix< T, M, N > transpose(Matrix< T, N, M > const &A)
Definition MiniTensor_Matrix.h:1791
minitensor::divide
KOKKOS_INLINE_FUNCTION void divide(TensorBase< R, SR > const &A, S const &s, TensorBase< T, ST > &B)
Definition MiniTensor_TensorBase.h:1811
minitensor::dot2
KOKKOS_INLINE_FUNCTION Tensor< typename Promote< S, T >::type, N > dot2(Tensor3< T, N > const &A, Vector< S > const &u)
minitensor::subtract
KOKKOS_INLINE_FUNCTION void subtract(TensorBase< R, SR > const &A, TensorBase< S, SS > const &B, TensorBase< T, ST > &C)
Definition MiniTensor_TensorBase.h:1716
minitensor::Filler::ZEROS
@ ZEROS
minitensor::inverse
KOKKOS_INLINE_FUNCTION Tensor< T, N > inverse(Tensor< T, N > const &A)
Definition MiniTensor_Inverse.h:182
minitensor::Index
uint32_t Index
Indexing type.
Definition MiniTensor_Traits.h:50
minitensor::DYNAMIC
constexpr Index DYNAMIC
Indicator for dynamic storage.
Definition MiniTensor_Traits.h:82
minitensor::apply_tensor4
Definition MiniTensor_Traits.h:192
minitensor::dimension_square
Integer square for manipulations between 2nd and 4rd-order tensors.
Definition MiniTensor_Storage.h:136