3088                                                              {
 3090 
 3091  if (
tagSense != getSkeletonSense())
 
 3093 
 3094  auto get_tag = [&](auto name) {
 3095    auto &mob = getPtrFE()->mField.get_moab();
 3098    return tag;
 3099  };
 3100 
 3101  auto get_tag_value = [&](auto &&tag, int dim) {
 3102    auto &mob = getPtrFE()->mField.get_moab();
 3103    auto face = getSidePtrFE()->getFEEntityHandle();
 3104    std::vector<double> value(dim);
 3105    CHK_MOAB_THROW(mob.tag_get_data(tag, &face, 1, value.data()), 
"set tag");
 
 3106    return value;
 3107  };
 3108 
 3109  auto create_tag = [this](const std::string tag_name, const int size) {
 3110    double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
 3113                                       th, MB_TAG_CREAT | MB_TAG_SPARSE,
 
 3114                                       def_VAL);
 3116  };
 3117 
 3118  Tag th_cauchy_streess = create_tag(
"CauchyStress", 9);
 
 3119  Tag th_detF = create_tag(
"detF", 1);
 
 3120  Tag th_traction = create_tag(
"traction", 3);
 
 3121  Tag th_disp_error = create_tag(
"DisplacementError", 1);
 
 3122 
 3123  Tag th_energy = create_tag(
"Energy", 1);
 
 3124 
 3125  auto t_w = getFTensor1FromMat<3>(
dataAtPts->wL2AtPts);
 
 3126  auto t_h = getFTensor2FromMat<3, 3>(
dataAtPts->hAtPts);
 
 3127  auto t_approx_P = getFTensor2FromMat<3, 3>(
dataAtPts->approxPAtPts);
 
 3128 
 3129  auto t_normal = getFTensor1NormalsAtGaussPts();
 3130  auto t_disp = getFTensor1FromMat<3>(
dataAtPts->wH1AtPts);
 
 3131 
 3132  auto next = [&]() {
 3133    ++t_w;
 3134    ++t_h;
 3135    ++t_approx_P;
 3136    ++t_normal;
 3137    ++t_disp;
 3138  };
 3139 
 3140  if (
dataAtPts->energyAtPts.size() == 0) {
 
 3141    
 3142    dataAtPts->energyAtPts.resize(getGaussPts().size2());
 
 3144  }
 3146 
 3151 
 3152  auto set_float_precision = [](const double x) {
 3153    if (std::abs(x) < std::numeric_limits<float>::epsilon())
 3154      return 0.;
 3155    else
 3156      return x;
 3157  };
 3158 
 3159  
 3160  auto save_scal_tag = [&](
auto &
th, 
auto v, 
const int gg) {
 
 3162    v = set_float_precision(
v);
 
 3165  };
 3166 
 3167  
 3170  auto save_vec_tag = [&](
auto &
th, 
auto &t_d, 
const int gg) {
 
 3173    for (
auto &
a : 
v.data())
 
 3174      a = set_float_precision(
a);
 
 3176                                     &*
v.data().begin());
 
 3178  };
 3179 
 3180  
 3181 
 3184      &
m(0, 0), &
m(0, 1), &
m(0, 2),
 
 3185 
 3186      &
m(1, 0), &
m(1, 1), &
m(1, 2),
 
 3187 
 3188      &
m(2, 0), &
m(2, 1), &
m(2, 2));
 
 3189 
 3190  auto save_mat_tag = [&](
auto &
th, 
auto &t_d, 
const int gg) {
 
 3192    t_m(
i, 
j) = t_d(
i, 
j);
 
 3193    for (
auto &
v : 
m.data())
 
 3194      v = set_float_precision(
v);
 
 3196                                     &*
m.data().begin());
 
 3198  };
 3199 
 3200  const auto nb_gauss_pts = getGaussPts().size2();
 3201  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
 3202 
 3204    t_traction(
i) = t_approx_P(
i, 
j) * t_normal(
j) / t_normal.
l2();
 
 3205    
 3206    CHKERR save_vec_tag(th_traction, t_traction, gg);
 
 3207 
 3208    double u_error = sqrt((t_disp(
i) - t_w(
i)) * (t_disp(
i) - t_w(
i)));
 
 3209    CHKERR save_scal_tag(th_disp_error, u_error, gg);
 
 3210    CHKERR save_scal_tag(th_energy, t_energy, gg);
 
 3211 
 3214    t_cauchy(
i, 
j) = (1. / jac) * (t_approx_P(
i, 
k) * t_h(
j, 
k));
 
 3215    CHKERR save_mat_tag(th_cauchy_streess, t_cauchy, gg);
 
 3217 
 3218    next();
 3219  }
 3220 
 3222}
#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(ublas::vector< T, A > &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