Get second directive of matrix.
\[
LS_{klmn} =
S_{ij} \frac{\partial^2 B_{ij}}{\partial A_{kl} \partial A_{mn} }
\]
940 {
941
942 for (auto aa = 0; aa != Dim; ++aa)
944
945 for (auto aa = 0; aa != Dim; ++aa)
947
948 for (auto aa = 0; aa != Dim; ++aa)
950
951 for (auto aa = 0; aa != Dim; ++aa)
952 for (auto bb = 0; bb != aa; ++bb) {
954 aF(bb, aa) = -
aF(aa, bb);
955 aF2(aa, bb) =
aF(aa, bb) *
aF(aa, bb);
956 }
957
962
963 for (auto aa = 0; aa != Dim; ++aa) {
964 for (auto bb = 0; bb != Dim; ++bb) {
965 if (aa != bb) {
967 const auto &
M =
aM[aa];
969 for (auto mm = 0; mm != Dim; ++mm) {
970 for (auto nn = mm; nn != Dim; ++nn) {
972 S(
i,
j,
k,
l) *
M(mm, nn);
973 }
974 }
975 }
976 }
977 }
978
979 for (auto aa = 0; aa != Dim; ++aa) {
980 for (auto mm = 0; mm != Dim; ++mm) {
981 for (auto nn = mm; nn != Dim; ++nn) {
983 }
984 }
985 }
986
987 if constexpr (NB == 3)
988 for (auto aa = 0; aa != Dim; ++aa) {
989 for (auto bb = 0; bb != Dim; ++bb) {
990 if (aa != bb) {
992 for (auto mm = 0; mm != Dim; ++mm) {
993 for (auto nn = mm; nn != Dim; ++nn) {
997 }
998 }
999 }
1000 }
1001 }
1002
1003 if constexpr (NB == 2)
1004 for (auto aa = 0; aa != Dim; ++aa) {
1005 for (auto bb = 0; bb != Dim; ++bb) {
1006 if (aa != bb) {
1008 if (aa == 1 || bb == 1)
1010 else
1012 for (auto mm = 0; mm != Dim; ++mm) {
1013 for (auto nn = mm; nn != Dim; ++nn) {
1017 }
1018 }
1019 }
1020 }
1021 }
1022
1023 if constexpr (NB == 1)
1024 for (auto aa = 0; aa != Dim; ++aa) {
1025 for (auto bb = 0; bb != Dim; ++bb) {
1026 if (aa != bb) {
1028 for (auto mm = 0; mm != Dim; ++mm) {
1029 for (auto nn = mm; nn != Dim; ++nn) {
1033 }
1034 }
1035 }
1036 }
1037 }
1038
1039 for (auto aa = 0; aa != Dim; ++aa) {
1040 for (auto mm = 0; mm != Dim; ++mm) {
1041 for (auto nn = 0; nn != Dim; ++nn) {
1042 if (nn != mm)
1044 }
1045 }
1046 }
1047
1048 if constexpr (NB == 3)
1049 for (auto aa = 0; aa != Dim; ++aa) {
1050 for (auto bb = 0; bb != Dim; ++bb) {
1051 if (aa != bb) {
1053 const auto v0 =
aF(aa, bb);
1054 for (auto cc = 0; cc != Dim; ++cc) {
1055 for (
auto dd = 0;
dd != Dim; ++
dd) {
1056 if (cc != dd) {
1057 const double v1 =
fVal(cc) *
aF(cc, dd);
1059 (v1 * v0) * S(
i,
j,
k,
l);
1060 }
1061 }
1062 }
1063 }
1064 }
1065 }
1066
1067 if constexpr (NB == 2)
1068 for (auto aa = 0; aa != Dim; ++aa) {
1069 for (auto bb = 0; bb != Dim; ++bb) {
1070 if (aa != bb)
1071 for (auto cc = 0; cc != Dim; ++cc) {
1072 for (
auto dd = 0;
dd != Dim; ++
dd)
1073 if (cc != dd) {
1074
1076
1077 if ((cc == 1 || dd == 1) && (aa == 1 || bb == 1))
1079 else if (cc != 1 && dd != 1 && aa != 1 && bb != 1) {
1080
1081 if ((aa != bb && bb != dd) && (aa != dd && bb != cc))
1083 else
1085
1086 } else if ((cc != 1 && dd != 1) && (aa == 1 || bb == 1))
1088 else if ((cc == 1 || dd == 1) && (aa != 1 && bb != 1)) {
1089
1090 if ((cc == 2 && dd == 1) || (cc == 2 && dd == 1))
1092
1094
1096
1097 ) *
1099
1100 else
1102
1103 } else
1105
1106 if (r)
1109 }
1110 }
1111 }
1112 }
1113
1114 if constexpr (NB == 1)
1115 for (auto aa = 0; aa != Dim; ++aa) {
1116 for (auto bb = 0; bb != Dim; ++bb) {
1117 if (aa != bb) {
1118 for (auto cc = 0; cc != Dim; ++cc) {
1119 for (
auto dd = 0;
dd != Dim; ++
dd) {
1120 if (cc != dd) {
1121 if ((bb != dd) && (aa !=
dd && bb != cc)) {
1122 const double r =
ddfVal(cc) / 4;
1125 }
1126 }
1127 }
1128 }
1129 }
1130 }
1131 }
1132
1133 using THIS = EigenMatrixImp<T1, T2, NB, Dim>;
1135
1136 T3 t_diff_A;
1137 GetDiffDiffMatImpl<THIS, V, T3, T>(*this, t_diff_A, t_S)
1138 .set(Number<Dim>(), Number<Dim>(), Number<Dim>(), Number<Dim>());
1139 return t_diff_A;
1140 }
const double v
phase velocity of light in medium (cm/ns)
auto get_nodiag_index(const Number< N1 > &, const Number< N2 > &, const Number< Dim > &)
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
FTensor::Tensor1< V, Dim > fVal
FTensor::Ddg< V, Dim, Dim > d2MType0[Dim][(Dim *(Dim+1))/2]
FTensor::Tensor2_symmetric< V, Dim > aF2
FTensor::Ddg< V, Dim, Dim > d2MType1[Dim][(Dim *(Dim+1))/2]
FTensor::Tensor1< V, Dim > ddfVal
FTensor::Ddg< V, Dim, Dim > aSM[(Dim - 1) *Dim][(Dim *(Dim+1))/2]
FTensor::Tensor2< V, Dim, Dim > aF
FTensor::Tensor1< V, Dim > dfVal