3912 {
3914
3915 if (
tagSense != getSkeletonSense())
3917
3918 auto get_tag = [&](auto name) {
3919 auto &mob = getPtrFE()->mField.get_moab();
3922 return tag;
3923 };
3924
3925 auto get_tag_value = [&](auto &&tag, int dim) {
3926 auto &mob = getPtrFE()->mField.get_moab();
3927 auto face = getSidePtrFE()->getFEEntityHandle();
3928 std::vector<double> value(dim);
3929 CHK_MOAB_THROW(mob.tag_get_data(tag, &face, 1, value.data()),
"set tag");
3930 return value;
3931 };
3932
3933 auto create_tag = [this](const std::string tag_name, const int size) {
3934 double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
3937 th, MB_TAG_CREAT | MB_TAG_SPARSE,
3938 def_VAL);
3940 };
3941
3942 Tag th_cauchy_streess = create_tag(
"CauchyStress", 9);
3943 Tag th_detF = create_tag(
"detF", 1);
3944 Tag th_traction = create_tag(
"traction", 3);
3945 Tag th_disp_error = create_tag(
"DisplacementError", 1);
3946
3947 Tag th_energy = create_tag(
"Energy", 1);
3948 Tag th_young_modulus = create_tag(
"YoungModulus", 1);
3949
3950 const auto nb_gauss_pts = getGaussPts().size2();
3951 auto t_w =
dataAtPts->getFTensorSmallWL2(nb_gauss_pts);
3952 auto t_h =
dataAtPts->getFTensorSmallH(nb_gauss_pts);
3953 auto t_approx_P =
dataAtPts->getFTensorApproxP(nb_gauss_pts);
3954
3955 auto t_normal = getFTensor1NormalsAtGaussPts();
3956 auto t_disp =
dataAtPts->getFTensorSmallWH1(nb_gauss_pts);
3957
3958
3959
3960 if (
dataAtPts->energyAtPts.size() == 0) {
3961
3962 dataAtPts->energyAtPts.resize(nb_gauss_pts);
3964 }
3967
3968 auto next = [&]() {
3969 ++t_w;
3970 ++t_h;
3971 ++t_approx_P;
3972 ++t_normal;
3973 ++t_disp;
3974 ++t_youngs_modulus;
3975 ++t_energy;
3976 };
3977
3982
3983 auto set_float_precision = [](const double x) {
3984 if (std::abs(x) < std::numeric_limits<float>::epsilon())
3985 return 0.;
3986 else
3987 return x;
3988 };
3989
3990
3991 auto save_scal_tag = [&](
auto &
th,
auto v,
const int gg) {
3993 v = set_float_precision(
v);
3996 };
3997
3998
4001 auto save_vec_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
4004 for (
auto &
a :
v.data())
4005 a = set_float_precision(
a);
4007 &*
v.data().begin());
4009 };
4010
4011
4012
4015 &
m(0, 0), &
m(0, 1), &
m(0, 2),
4016
4017 &
m(1, 0), &
m(1, 1), &
m(1, 2),
4018
4019 &
m(2, 0), &
m(2, 1), &
m(2, 2));
4020
4021 auto save_mat_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
4023 t_m(
i,
j) = t_d(
i,
j);
4024 for (
auto &
v :
m.data())
4025 v = set_float_precision(
v);
4027 &*
m.data().begin());
4029 };
4030
4031 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
4032
4034 t_traction(
i) = t_approx_P(
i,
j) * t_normal(
j) / t_normal.
l2();
4035
4037 CHKERR save_vec_tag(th_traction, t_traction, gg);
4038
4039 double u_error = sqrt((t_disp(
i) - t_w(
i)) * (t_disp(
i) - t_w(
i)));
4040 if (!std::isfinite(u_error))
4041 u_error = -1.;
4042 CHKERR save_scal_tag(th_disp_error, u_error, gg);
4043 CHKERR save_scal_tag(th_energy, t_energy, gg);
4045 CHKERR save_scal_tag(th_young_modulus, t_youngs_modulus, gg);
4046
4049 t_cauchy(
i,
j) = (1. / jac) * (t_approx_P(
i,
k) * t_h(
j,
k));
4050 CHKERR save_mat_tag(th_cauchy_streess, t_cauchy, gg);
4052
4053 next();
4054 }
4055
4057}
#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 CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
#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
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
VectorBoundedArray< double, 3 > VectorDouble3
static auto getFTensor0FromVec(V &data)
Get tensor rank 0 (scalar) form data vector.
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
FTensor::Index< 'm', 3 > m
static bool hasNonHomogeneousMaterialBlock
