1397 {
1399
1400 if (type != MBVERTEX) {
1402 }
1403
1404 auto tensor_to_tensor = [](const auto &t1, auto &t2) {
1405 t2(0, 0) = t1(0, 0);
1406 t2(1, 1) = t1(1, 1);
1407 t2(2, 2) = t1(2, 2);
1408 t2(0, 1) = t2(1, 0) = t1(1, 0);
1409 t2(0, 2) = t2(2, 0) = t1(2, 0);
1410 t2(1, 2) = t2(2, 1) = t1(2, 1);
1411 };
1412
1413 std::array<double, 9> def_val;
1414 def_val.fill(0);
1415
1416 auto make_tag = [&](auto name, auto size) {
1419 MB_TAG_CREAT | MB_TAG_SPARSE,
1420 def_val.data());
1422 };
1423
1424 auto th_stress = make_tag("STRESS", 9);
1425 auto th_psi = make_tag("ENERGY", 1);
1426
1427 const int nb_integration_pts =
mapGaussPts.size();
1428
1433
1436
1437 dataAtPts->stiffnessMat->resize(36, 1,
false);
1440
1444 if (type == MBTRI || type == MBQUAD) {
1446 auto &m_field = this->getPtrFE()->mField;
1447 CHKERR m_field.get_moab().get_adjacencies(&ent, 1, 3,
false, ents,
1448 moab::Interface::UNION);
1449#ifndef NDEBUG
1450 if (ents.empty())
1452 "Could not find a 3D element adjacent to a given face element");
1453#endif
1454 ent_3d = ents.front();
1455 }
1456
1457 bool found_block = false;
1458 int block_id = -1;
1460 if (
m.second.tEts.find(ent_3d) !=
m.second.tEts.end()) {
1461 const double young =
m.second.E;
1462 const double poisson =
m.second.PoissonRatio;
1463 const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
1464 block_id =
m.second.iD;
1465
1466 t_D(
i,
j,
k,
l) = 0.;
1467 t_D(0, 0, 0, 0) = t_D(1, 1, 1, 1) = t_D(2, 2, 2, 2) = 1 - poisson;
1468 t_D(0, 1, 0, 1) = t_D(0, 2, 0, 2) = t_D(1, 2, 1, 2) =
1469 0.5 * (1 - 2 * poisson);
1470 t_D(0, 0, 1, 1) = t_D(1, 1, 0, 0) = t_D(0, 0, 2, 2) = t_D(2, 2, 0, 0) =
1471 t_D(1, 1, 2, 2) = t_D(2, 2, 1, 1) = poisson;
1472 t_D(
i,
j,
k,
l) *= coefficient;
1473
1474 found_block = true;
1475 break;
1476 }
1477 }
1478 if (!found_block)
1480 "Element not found in any of material blocksets");
1481
1482 int def_val_int = 0;
1485 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val_int);
1486 double detH = 0.;
1493
1494 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1495
1497 t_h(0, 0) += 1;
1498 t_h(1, 1) += 1;
1499 t_h(2, 2) += 1;
1500 }
1501
1503 t_small_strain_symm(
i,
j) = (t_h(
i,
j) || t_h(
j,
i)) / 2.;
1504 } else {
1507 t_F(
i,
j) = t_h(
i,
k) * t_invH(
k,
j);
1508 t_small_strain_symm(
i,
j) = (t_F(
i,
j) || t_F(
j,
i)) / 2.;
1509 ++t_H;
1510 }
1511
1512 t_small_strain_symm(0, 0) -= 1;
1513 t_small_strain_symm(1, 1) -= 1;
1514 t_small_strain_symm(2, 2) -= 1;
1515
1516
1517 t_stress_symm(
i,
j) = t_D(
i,
j,
k,
l) * t_small_strain_symm(
k,
l);
1518 tensor_to_tensor(t_stress_symm, t_stress);
1519
1520 const double psi = 0.5 * t_stress_symm(
i,
j) * t_small_strain_symm(
i,
j);
1521
1524 &t_stress(0, 0));
1526 &block_id);
1527
1528 ++t_h;
1529 }
1530
1532}
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#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.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
auto type_from_handle(const EntityHandle h)
get type from entity handle
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, Tensor_Dim1, Tensor_Dim2 > getFTensor2FromMat(MatrixDouble &data)
Get tensor rank 2 (matrix) form data matrix.
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
FTensor::Index< 'm', 3 > m