1120 {
1122
1129
1130 auto nb_integration_pts =
stressPtr->size1();
1131#ifndef NDEBUG
1132 if (nb_integration_pts != getGaussPts().size2()) {
1134 "inconsistent number of integration points");
1135 }
1136#endif
1137
1140
1141 auto get_strain =
1144 auto t_strain = get_strain();
1146 auto t_inv_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, STRIDEMATD>(
1148#ifndef NDEBUG
1149 auto t_D =
1150 getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, STRIDEMATD>(
dataAtPts->matD);
1151#endif
1152
1153
1154
1155 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
1156 t_strain(
i,
j) = t_inv_D(
i,
j,
k,
l) * t_stress(
k,
l);
1157
1158#ifndef NDEBUG
1160 t_stress_symm_debug(
i,
j) = (t_stress(
i,
j) || t_stress(
j,
i)) / 2;
1162 t_stress_symm_debug_diff(
i,
j) =
1163 t_D(
i,
j,
k,
l) * t_strain(
k,
l) - t_stress_symm_debug(
i,
j);
1164 double nrm =
1165 t_stress_symm_debug_diff(
i,
j) * t_stress_symm_debug_diff(
i,
j);
1166 double nrm0 = t_stress_symm_debug(
i,
j) * t_stress_symm_debug(
i,
j) +
1167 std::numeric_limits<double>::epsilon();
1168 constexpr double eps = 1e-10;
1169 if (std::fabs(std::sqrt(nrm / nrm0)) >
eps) {
1171 << "Stress symmetry check failed: " << std::endl
1172 << t_stress_symm_debug_diff << std::endl
1173 << t_stress;
1175 "Norm is too big: " + std::to_string(nrm / nrm0));
1176 }
1177 ++t_D;
1178#endif
1179
1180 ++t_strain;
1181 ++t_stress;
1182 ++t_inv_D;
1183 }
1184
1186}
#define FTENSOR_INDEX(DIM, I)
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MOFEM_LOG(channel, severity)
Log.
FTensor::Index< 'i', SPACE_DIM > i
const double n
refractive index of diffusive medium
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
DataLayoutTraits< DataLayout::GaussByCoeffs > DL
decltype(GetFTensor2SymmetricFromMatImpl< Tensor_Dim, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor2SymmetricFromMatType
auto getFTensor2FromMat(M &data)
Get tensor rank 2 (matrix) form data matrix.
FTensor::Index< 'm', 3 > m
