1044 int nb_integration_pts = row_data.
getN().size1();
1045 int row_nb_dofs = row_data.
getIndices().size();
1046 int col_nb_dofs = col_data.
getIndices().size();
1050 &
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2), &
m(
r + 0,
c + 3),
1051 &
m(
r + 0,
c + 4), &
m(
r + 0,
c + 5),
1053 &
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2), &
m(
r + 1,
c + 3),
1054 &
m(
r + 1,
c + 4), &
m(
r + 1,
c + 5),
1056 &
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2), &
m(
r + 2,
c + 3),
1057 &
m(
r + 2,
c + 4), &
m(
r + 2,
c + 5),
1059 &
m(
r + 3,
c + 0), &
m(
r + 3,
c + 1), &
m(
r + 3,
c + 2), &
m(
r + 3,
c + 3),
1060 &
m(
r + 3,
c + 4), &
m(
r + 3,
c + 5),
1062 &
m(
r + 4,
c + 0), &
m(
r + 4,
c + 1), &
m(
r + 4,
c + 2), &
m(
r + 4,
c + 3),
1063 &
m(
r + 4,
c + 4), &
m(
r + 4,
c + 5),
1065 &
m(
r + 5,
c + 0), &
m(
r + 5,
c + 1), &
m(
r + 5,
c + 2), &
m(
r + 5,
c + 3),
1066 &
m(
r + 5,
c + 4), &
m(
r + 5,
c + 5)
1081 auto v = getVolume();
1082 auto t_w = getFTensor0IntegrationWeight();
1086 auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
1087 auto t_approx_P = getFTensor2FromMat<3, 3>(
dataAtPts->approxPAtPts);
1088 auto t_P = getFTensor2FromMat<3, 3>(
dataAtPts->PAtPts);
1090 getFTensor2SymmetricFromMat<3>(
dataAtPts->logStreachDotTensorAtPts);
1091 auto t_u = getFTensor2SymmetricFromMat<3>(
dataAtPts->streachTensorAtPts);
1093 getFTensor4DdgFromMat<3, 3, 1>(
dataAtPts->diffStreachTensorAtPts);
1094 auto t_eigen_vals = getFTensor1FromMat<3>(
dataAtPts->eigenVals);
1095 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(
dataAtPts->eigenVecs);
1100 t_one4(
i,
j,
k,
l) = t_kd(
j,
l) * t_kd(
i,
k);
1103 t_one4_symmetric(
i,
j,
k,
l) = 0;
1104 for (
auto ii : {0, 1, 2})
1105 for (
auto jj : {0, 1, 2})
1106 for (
auto kk : {0, 1, 2})
1107 for (
auto ll : {0, 1, 2}) {
1110 t_one4_symmetric(ii, jj, kk, ll) =
1111 t_one4(ii, jj, kk, ll) + t_one4(ii, jj, ll, kk);
1113 t_one4_symmetric(ii, jj, kk, ll) = t_one4(ii, jj, kk, ll);
1116 int row_nb_base_functions = row_data.
getN().size2();
1121 for (
auto ii : {0, 1, 2})
1124 const double ts_a = getTSa();
1125 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
1129 t_P_dh(
i,
j,
N0,
k) = t_P_dh0(
i,
j,
k);
1130 t_P_dh(
i,
j,
N1,
k) = t_P_dh1(
i,
j,
k);
1131 t_P_dh(
i,
j,
N2,
k) = t_P_dh2(
i,
j,
k);
1134 t_RTP(
i,
j) = t_R(
k,
i) * t_approx_P(
k,
j);
1136 t_deltaP(
i,
j) = t_RTP(
i,
j) - t_P(
i,
j);
1139 t_dot_u(
i,
j) = t_u(
i,
k) * t_dot_log_u(
k,
j);
1142 t_stress_diff(
i,
j,
k,
l) = 0;
1143 for (
int ii = 0; ii != 3; ++ii)
1144 for (
int jj = 0; jj != 3; ++jj)
1145 for (
int kk = 0; kk != 3; ++kk)
1146 for (
int ll = 0; ll != 3; ++ll)
1147 for (
int mm = 0; mm != 3; ++mm)
1148 for (
int nn = 0; nn != 3; ++nn)
1149 t_stress_diff(ii, jj, mm, nn) -=
1150 t_P_dh(ii, jj, kk, ll) * t_diff_u(kk, ll, mm, nn);
1153 t_viscous_stress(
i,
j) =
alphaU * t_dot_log_u(
i,
j);
1155 t_viscous_stress_diff(
i,
j,
k,
l) = t_one4_symmetric(
i,
j,
k,
l);
1156 t_viscous_stress_diff(
i,
j,
k,
l) *= (ts_a *
alphaU);
1159 t_deltaP2(
i,
j) = t_deltaP(
i,
j) - t_viscous_stress(
i,
j);
1162 t_eigen_vals, t_eigen_vecs,
f,
d_f,
dd_f, t_deltaP2, nbUniq[gg]);
1163 for (
int ii = 0; ii != 3; ++ii) {
1164 for (
int jj = 0; jj < ii; ++jj) {
1165 for (
int kk = 0; kk != 3; ++kk) {
1166 for (
int ll = 0; ll != kk; ++ll) {
1167 t_diff2_uP2(ii, jj, kk, ll) *= 2;
1172 for (
int ii = 0; ii != 3; ++ii) {
1173 for (
int jj = ii; jj != 3; ++jj) {
1174 for (
int kk = 0; kk != 3; ++kk) {
1175 for (
int ll = 0; ll < kk; ++ll) {
1176 t_diff2_uP2(ii, jj, kk, ll) *= 2;
1183 for (; rr != row_nb_dofs / 6; ++rr) {
1186 auto t_m = get_ftensor4(
K, 6 * rr, 0);
1188 for (
int cc = 0; cc != col_nb_dofs / 6; ++cc) {
1190 const double b = a * t_row_base_fun * t_col_base_fun;
1192 for (
int ii = 0; ii != 3; ++ii)
1193 for (
int jj = ii; jj != 3; ++jj)
1194 for (
int kk = 0; kk != 3; ++kk)
1195 for (
int ll = kk; ll != 3; ++ll)
1196 for (
int mm = 0; mm != 3; ++mm)
1197 for (
int nn = 0; nn != 3; ++nn)
1198 t_m(ii, jj, kk, ll) +=
1199 b * t_diff_u(mm, nn, ii, jj) *
1200 (t_stress_diff(mm, nn, kk, ll) -
1201 t_viscous_stress_diff(mm, nn, kk, ll));
1203 t_m(
i,
j,
k,
l) +=
b * t_diff2_uP2(
i,
j,
k,
l);
1211 for (; rr != row_nb_base_functions; ++rr)