1045 {
1047
1048 int nb_integration_pts = row_data.
getN().size1();
1049 int row_nb_dofs = row_data.
getIndices().size();
1050 int col_nb_dofs = col_data.
getIndices().size();
1053
1054 &
m(r + 0,
c + 0), &
m(r + 0,
c + 1), &
m(r + 0,
c + 2), &
m(r + 0,
c + 3),
1055 &
m(r + 0,
c + 4), &
m(r + 0,
c + 5),
1056
1057 &
m(r + 1,
c + 0), &
m(r + 1,
c + 1), &
m(r + 1,
c + 2), &
m(r + 1,
c + 3),
1058 &
m(r + 1,
c + 4), &
m(r + 1,
c + 5),
1059
1060 &
m(r + 2,
c + 0), &
m(r + 2,
c + 1), &
m(r + 2,
c + 2), &
m(r + 2,
c + 3),
1061 &
m(r + 2,
c + 4), &
m(r + 2,
c + 5),
1062
1063 &
m(r + 3,
c + 0), &
m(r + 3,
c + 1), &
m(r + 3,
c + 2), &
m(r + 3,
c + 3),
1064 &
m(r + 3,
c + 4), &
m(r + 3,
c + 5),
1065
1066 &
m(r + 4,
c + 0), &
m(r + 4,
c + 1), &
m(r + 4,
c + 2), &
m(r + 4,
c + 3),
1067 &
m(r + 4,
c + 4), &
m(r + 4,
c + 5),
1068
1069 &
m(r + 5,
c + 0), &
m(r + 5,
c + 1), &
m(r + 5,
c + 2), &
m(r + 5,
c + 3),
1070 &
m(r + 5,
c + 4), &
m(r + 5,
c + 5)
1071
1072 );
1073 };
1080
1084
1085 auto v = getVolume();
1086 auto t_w = getFTensor0IntegrationWeight();
1090 auto t_R = getFTensor2FromMat<3, 3>(
dataAtPts->rotMatAtPts);
1091 auto t_approx_P = getFTensor2FromMat<3, 3>(
dataAtPts->approxPAtPts);
1092 auto t_P = getFTensor2FromMat<3, 3>(
dataAtPts->PAtPts);
1093 auto t_dot_log_u =
1094 getFTensor2SymmetricFromMat<3>(
dataAtPts->logStreachDotTensorAtPts);
1095 auto t_u = getFTensor2SymmetricFromMat<3>(
dataAtPts->streachTensorAtPts);
1096 auto t_diff_u =
1097 getFTensor4DdgFromMat<3, 3, 1>(
dataAtPts->diffStreachTensorAtPts);
1098 auto t_eigen_vals = getFTensor1FromMat<3>(
dataAtPts->eigenVals);
1099 auto t_eigen_vecs = getFTensor2FromMat<3, 3>(
dataAtPts->eigenVecs);
1101
1105
1107 t_one4_symmetric(
i,
j,
k,
l) = 0;
1108 for (auto ii : {0, 1, 2})
1109 for (auto jj : {0, 1, 2})
1110 for (auto kk : {0, 1, 2})
1111 for (auto ll : {0, 1, 2}) {
1112
1113 if (ll != kk)
1114 t_one4_symmetric(ii, jj, kk, ll) =
1115 t_one4(ii, jj, kk, ll) + t_one4(ii, jj, ll, kk);
1116 else
1117 t_one4_symmetric(ii, jj, kk, ll) = t_one4(ii, jj, kk, ll);
1118 }
1119
1120 int row_nb_base_functions = row_data.
getN().size2();
1122
1125 for (auto ii : {0, 1, 2})
1126 t_one(ii, ii) = 1;
1127
1128 const double ts_a = getTSa();
1129 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1131
1133 t_P_dh(
i,
j,
N0,
k) = t_P_dh0(
i,
j,
k);
1134 t_P_dh(
i,
j,
N1,
k) = t_P_dh1(
i,
j,
k);
1135 t_P_dh(
i,
j,
N2,
k) = t_P_dh2(
i,
j,
k);
1136
1138 t_RTP(
i,
j) = t_R(
k,
i) * t_approx_P(
k,
j);
1140 t_deltaP(
i,
j) = t_RTP(
i,
j) - t_P(
i,
j);
1141
1143 t_dot_u(
i,
j) = t_u(
i,
k) * t_dot_log_u(
k,
j);
1144
1146 t_stress_diff(
i,
j,
k,
l) = 0;
1147 for (int ii = 0; ii != 3; ++ii)
1148 for (int jj = 0; jj != 3; ++jj)
1149 for (int kk = 0; kk != 3; ++kk)
1150 for (int ll = 0; ll != 3; ++ll)
1151 for (int mm = 0; mm != 3; ++mm)
1152 for (int nn = 0; nn != 3; ++nn)
1153 t_stress_diff(ii, jj, mm, nn) -=
1154 t_P_dh(ii, jj, kk, ll) * t_diff_u(kk, ll, mm, nn);
1155
1157 t_viscous_stress(
i,
j) =
alphaU * t_dot_log_u(
i,
j);
1159 t_viscous_stress_diff(
i,
j,
k,
l) = t_one4_symmetric(
i,
j,
k,
l);
1160 t_viscous_stress_diff(
i,
j,
k,
l) *= (ts_a *
alphaU);
1161
1163 t_deltaP2(
i,
j) = t_deltaP(
i,
j) - t_viscous_stress(
i,
j);
1164
1166 t_eigen_vals, t_eigen_vecs,
f,
d_f,
dd_f, t_deltaP2, nbUniq[gg]);
1167 for (int ii = 0; ii != 3; ++ii) {
1168 for (int jj = 0; jj < ii; ++jj) {
1169 for (int kk = 0; kk != 3; ++kk) {
1170 for (int ll = 0; ll != kk; ++ll) {
1171 t_diff2_uP2(ii, jj, kk, ll) *= 2;
1172 }
1173 }
1174 }
1175 }
1176 for (int ii = 0; ii != 3; ++ii) {
1177 for (int jj = ii; jj != 3; ++jj) {
1178 for (int kk = 0; kk != 3; ++kk) {
1179 for (int ll = 0; ll < kk; ++ll) {
1180 t_diff2_uP2(ii, jj, kk, ll) *= 2;
1181 }
1182 }
1183 }
1184 }
1185
1186 int rr = 0;
1187 for (; rr != row_nb_dofs / 6; ++rr) {
1188
1190 auto t_m = get_ftensor4(
K, 6 * rr, 0);
1191
1192 for (int cc = 0; cc != col_nb_dofs / 6; ++cc) {
1193
1194 const double b =
a * t_row_base_fun * t_col_base_fun;
1195
1196 for (int ii = 0; ii != 3; ++ii)
1197 for (int jj = ii; jj != 3; ++jj)
1198 for (int kk = 0; kk != 3; ++kk)
1199 for (int ll = kk; ll != 3; ++ll)
1200 for (int mm = 0; mm != 3; ++mm)
1201 for (int nn = 0; nn != 3; ++nn)
1202 t_m(ii, jj, kk, ll) +=
1203 b * t_diff_u(mm, nn, ii, jj) *
1204 (t_stress_diff(mm, nn, kk, ll) -
1205 t_viscous_stress_diff(mm, nn, kk, ll));
1206
1207 t_m(
i,
j,
k,
l) +=
b * t_diff2_uP2(
i,
j,
k,
l);
1208
1209 ++t_m;
1210 ++t_col_base_fun;
1211 }
1212
1213 ++t_row_base_fun;
1214 }
1215 for (; rr != row_nb_base_functions; ++rr)
1216 ++t_row_base_fun;
1217
1218 ++t_w;
1219 ++t_P_dh0;
1220 ++t_P_dh1;
1221 ++t_P_dh2;
1222 ++t_R;
1223 ++t_approx_P;
1224 ++t_P;
1225 ++t_dot_log_u;
1226 ++t_u;
1227 ++t_diff_u;
1228 ++t_eigen_vals;
1229 ++t_eigen_vecs;
1230 }
1232}
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
FTensor::Ddg< double, 3, 3 > getDiffDiffMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f, Fun< double > d_f, Fun< double > dd_f, FTensor::Tensor2< double, 3, 3 > &t_S, const int nb)
FTensor::Tensor3< FTensor::PackPtr< double *, 1 >, 3, 3, 3 > getFTensor3FromMat(MatrixDouble &m)
double dd_f(const double v)
double d_f(const double v)
MatrixDouble K
local tangent matrix
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of dofs on entity.