1379 {
1381
1388
1389 auto nb_integration_pts =
stressPtr->size1();
1390#ifndef NDEBUG
1391 if (nb_integration_pts != getGaussPts().size2()) {
1393 "inconsistent number of integration points");
1394 }
1395#endif
1396
1399
1400 auto get_strain =
1403 auto t_strain = get_strain();
1405 auto t_inv_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, STRIDEMATD>(
1407#ifndef NDEBUG
1408 auto t_D =
1409 getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, STRIDEMATD>(
dataAtPts->matD);
1410#endif
1411
1413
1414
1415
1416 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
1417 t_strain(
i,
j) = t_inv_D(
i,
j,
k,
l) * t_stress(
k,
l);
1418
1419#ifndef NDEBUG
1421 t_stress_symm_debug(
i,
j) = (t_stress(
i,
j) || t_stress(
j,
i)) / 2;
1423 t_stress_symm_debug_diff(
i,
j) =
1424 t_D(
i,
j,
k,
l) * t_strain(
k,
l) - t_stress_symm_debug(
i,
j);
1425 double nrm =
1426 t_stress_symm_debug_diff(
i,
j) * t_stress_symm_debug_diff(
i,
j);
1427 double nrm0 = t_stress_symm_debug(
i,
j) * t_stress_symm_debug(
i,
j) +
1428 std::numeric_limits<double>::epsilon();
1429 constexpr double eps = 1e-10;
1430 if (std::fabs(std::sqrt(nrm / nrm0)) >
eps) {
1432 << "Stress symmetry check failed: " << std::endl
1433 << t_stress_symm_debug_diff << std::endl
1434 << t_stress;
1436 "Norm is too big: " + std::to_string(nrm / nrm0));
1437 }
1438 ++t_D;
1439#endif
1440
1441 ++t_strain;
1442 ++t_stress;
1443 ++t_inv_D;
1444 }
1445
1447}
#define FTENSOR_INDEX(DIM, I)
Kronecker Delta class symmetric.
#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
