488 {
490
491 if (type != MBVERTEX) {
493 }
494
495 auto tensor_to_tensor = [](const auto &t1, auto &t2) {
496 t2(0, 0) = t1(0, 0);
497 t2(1, 1) = t1(1, 1);
498 t2(2, 2) = t1(2, 2);
499 t2(0, 1) = t2(1, 0) = t1(1, 0);
500 t2(0, 2) = t2(2, 0) = t1(2, 0);
501 t2(1, 2) = t2(2, 1) = t1(2, 1);
502 };
503
504 auto get_tag_handle = [&](auto name, auto size) {
506 std::vector<double> def_vals(size, 0.0);
508 MB_TAG_CREAT | MB_TAG_SPARSE,
509 def_vals.data());
511 };
512
513 auto th_stress = get_tag_handle("STRESS", 9);
514 auto th_strain = get_tag_handle("STRAIN", 9);
515 auto th_psi = get_tag_handle("ENERGY", 1);
516 auto th_disp = get_tag_handle("DISPLACEMENT", 3);
517
519
524
525 auto t_h = getFTensor2FromMat<3, 3>(*
dataAtPts->hMat);
526 auto t_H = getFTensor2FromMat<3, 3>(*
dataAtPts->HMat);
527
528 dataAtPts->stiffnessMat->resize(36, 1,
false);
532 const double young =
m.second.E;
533 const double poisson =
m.second.PoissonRatio;
534
535 const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
536
537 t_D(
i,
j,
k,
l) = 0.;
538 t_D(0, 0, 0, 0) = t_D(1, 1, 1, 1) = t_D(2, 2, 2, 2) = 1 - poisson;
539 t_D(0, 1, 0, 1) = t_D(0, 2, 0, 2) = t_D(1, 2, 1, 2) =
540 0.5 * (1 - 2 * poisson);
541 t_D(0, 0, 1, 1) = t_D(1, 1, 0, 0) = t_D(0, 0, 2, 2) = t_D(2, 2, 0, 0) =
542 t_D(1, 1, 2, 2) = t_D(2, 2, 1, 1) = poisson;
543 t_D(
i,
j,
k,
l) *= coefficient;
544
545 break;
546 }
547
548 double detH = 0.;
553
556
557 auto t_spat_pos = getFTensor1FromMat<3>(*
dataAtPts->spatPosMat);
558 auto t_mesh_node_pos = getFTensor1FromMat<3>(*
dataAtPts->meshNodePosMat);
560
561 for (int gg = 0; gg != nb_integration_pts; ++gg) {
562
564 t_h(0, 0) += 1;
565 t_h(1, 1) += 1;
566 t_h(2, 2) += 1;
567 }
568
570 t_small_strain_symm(
i,
j) = (t_h(
i,
j) || t_h(
j,
i)) / 2.;
571 } else {
574 t_F(
i,
j) = t_h(
i,
k) * t_invH(
k,
j);
575 t_small_strain_symm(
i,
j) = (t_F(
i,
j) || t_F(
j,
i)) / 2.;
576 ++t_H;
577 }
578
579 t_small_strain_symm(0, 0) -= 1;
580 t_small_strain_symm(1, 1) -= 1;
581 t_small_strain_symm(2, 2) -= 1;
582
583
584 t_stress_symm(
i,
j) = t_D(
i,
j,
k,
l) * t_small_strain_symm(
k,
l);
585 tensor_to_tensor(t_stress_symm, t_stress);
586 tensor_to_tensor(t_small_strain_symm, t_small_strain);
587
588 t_disp(
i) = t_spat_pos(
i) - t_mesh_node_pos(
i);
589
590 const double psi = 0.5 * t_stress_symm(
i,
j) * t_small_strain_symm(
i,
j);
591
594 &t_stress(0, 0));
596 &t_small_strain(0, 0));
598
599 ++t_h;
600 ++t_spat_pos;
601 ++t_mesh_node_pos;
602 }
603
605}
#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 ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FTensor::Index< 'm', SPACE_DIM > m
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
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.
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.