1456 {
1458
1465
1466 auto nb_integration_pts =
stressPtr->size1();
1467#ifndef NDEBUG
1468 if (nb_integration_pts != getGaussPts().size2()) {
1470 "inconsistent number of integration points");
1471 }
1472#endif
1473
1476
1477 auto get_strain =
1480 auto t_strain = get_strain();
1482 auto t_inv_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, STRIDEMATD>(
1484#ifndef NDEBUG
1485 auto t_D =
1486 getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, STRIDEMATD>(
dataAtPts->matD);
1487#endif
1488
1490
1491
1492
1493 for (auto gg = 0; gg != nb_integration_pts; ++gg) {
1494 t_strain(
i,
j) = t_inv_D(
i,
j,
k,
l) * t_stress(
k,
l);
1495
1496#ifndef NDEBUG
1498 t_stress_symm_debug(
i,
j) = (t_stress(
i,
j) || t_stress(
j,
i)) / 2;
1500 t_stress_symm_debug_diff(
i,
j) =
1501 t_D(
i,
j,
k,
l) * t_strain(
k,
l) - t_stress_symm_debug(
i,
j);
1502 double nrm =
1503 t_stress_symm_debug_diff(
i,
j) * t_stress_symm_debug_diff(
i,
j);
1504 double nrm0 = t_stress_symm_debug(
i,
j) * t_stress_symm_debug(
i,
j) +
1505 std::numeric_limits<double>::epsilon();
1506 constexpr double eps = 1e-10;
1507 if (std::fabs(std::sqrt(nrm / nrm0)) >
eps) {
1509 << "Stress symmetry check failed: " << std::endl
1510 << t_stress_symm_debug_diff << std::endl
1511 << t_stress;
1513 "Norm is too big: " + std::to_string(nrm / nrm0));
1514 }
1515 ++t_D;
1516#endif
1517
1518 ++t_strain;
1519 ++t_stress;
1520 ++t_inv_D;
1521 }
1522
1524}
#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
