24      : 
eP(ep), 
dataFieldEval(ep.mField.getInterface<FieldEvaluatorInterface>()
 
   30        "build field evaluator tree");
 
   32    auto no_rule = [](int, int, int) { 
return -1; };
 
   34    auto set_element_for_field_eval = [&]() {
 
   36      boost::shared_ptr<Ele> vol_ele(
dataFieldEval->feMethodPtr.lock());
 
   37      vol_ele->getRuleHook = no_rule;
 
   38      vol_ele->getUserPolynomialBase() =
 
   39          boost::make_shared<CGGUserPolynomialBase>();
 
   40      EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
 
   44      auto piola_scale_ptr = boost::make_shared<double>(1.0);
 
   45      vol_ele->getOpPtrVector().push_back(
 
   48      vol_ele->getOpPtrVector().push_back(
new OpCalculateHVecTensorField<3, 3>(
 
   50      vol_ele->getOpPtrVector().push_back(
 
   51          new OpCalculateHTensorTensorField<3, 3>(
 
   53              SmartPetscObj<Vec>(), MBMAXTYPE));
 
   55        vol_ele->getOpPtrVector().push_back(
 
   59        vol_ele->getOpPtrVector().push_back(
 
   60            new OpCalculateTensor2SymmetricFieldValues<3>(
 
   64      vol_ele->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<3>(
 
   66      vol_ele->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<3>(
 
   68      vol_ele->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<3>(
 
   72      vol_ele->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<3>(
 
   74      vol_ele->getOpPtrVector().push_back(
 
   75          new OpCalculateVectorFieldGradient<3, 3>(
 
   78      vol_ele->getOpPtrVector().push_back(
 
   83    auto set_element_for_post_process_energy = [&]() {
 
   97    auto reads_post_proc_data = [](
PtsHashMap &pts_hash_map) {
 
  102      if (!file.is_open()) {
 
  108      while (std::getline(file, line)) {
 
  109        std::istringstream iss(line);
 
  111        double col2, col3, col4;
 
  113        if (iss >> col1 >> col2 >> col3 >> col4) {
 
  114          MOFEM_LOG(
"EP", Sev::verbose) << 
"Read: " << col1 << 
", " << col2
 
  115                                        << 
", " << col3 << 
", " << col4;
 
  116          pts_hash_map[col1] = {col2, col3, col4};
 
  118          MOFEM_LOG(
"EP", Sev::error) << 
"Error parsing line: " << line;
 
  128                      "set element for post energy");
 
  130    PetscBool test_cook_flg = PETSC_FALSE;
 
  132                                          &test_cook_flg, PETSC_NULLPTR),
 
  133                      "get post proc points");
 
  136      ptsHashMap[
"Point B"] = {48. / 2., 44. + (60. - 44.) / 2., 0.};
 
  137      ptsHashMap[
"Point C"] = {48. / 2., (44. - 0.) / 2., 0.};
 
  146      ptsHashMap[
"Point (2.5, 0., 0.)"] = {2.5, 0, 0.};
 
  152                      "get write restart option");
 
 
  162    MOFEM_LOG(
"EP", Sev::inform) << 
"Monitor postProcess";
 
  165    auto get_ents_on_mesh_skin = [&]() {
 
  167      std::map<std::string, Range> boundary_entities_vec;
 
  168      Range boundary_entities;
 
  170        std::string meshset_name = it->getName();
 
  171        std::string block_name_disp = 
"SPATIAL_DISP";
 
  172        std::string block_name_fix = 
"FIX";
 
  173        if (meshset_name.compare(0, block_name_disp.size(), block_name_disp) == 0 ||
 
  174            meshset_name.compare(0, block_name_fix.size(), block_name_fix) == 0) {
 
  176                                                     boundary_entities, 
true);
 
  178          boundary_entities_vec[meshset_name] = boundary_entities;
 
  179          boundary_entities.clear();
 
  182      return boundary_entities_vec;
 
  185    auto boundary_entities_vec = get_ents_on_mesh_skin();
 
  187    std::vector<std::tuple<std::string, Range, std::array<double, 6>>>
 
  189    for (
const auto &pair : boundary_entities_vec) {
 
  190      reactionForces.push_back(
 
  191          std::make_tuple(pair.first, pair.second,
 
  192                          std::array<double, 6>{0.0, 0.0, 0.0, 0.0, 0.0, 0.0}));
 
  195    auto integration_rule_face = [](int, int, 
int approx_order) {
 
  199        boost::make_shared<FaceElementForcesAndSourcesCore>(
eP.
mField);
 
  200    auto no_rule = [](int, int, int) { 
return -1; };
 
  201    face_fe->getRuleHook = integration_rule_face;
 
  203    EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
 
  209    face_fe->getOpPtrVector().push_back(op_side);
 
  210    auto side_fe_ptr = op_side->getSideFEPtr();
 
  212        boost::make_shared<EshelbianPlasticity::CGGUserPolynomialBase>();
 
  213    side_fe_ptr->getUserPolynomialBase() = base_ptr;
 
  214    CHKERR EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
 
  217    auto piola_scale_ptr = boost::make_shared<double>(1.0);
 
  218    side_fe_ptr->getOpPtrVector().push_back(
 
  221    side_fe_ptr->getOpPtrVector().push_back(
 
  222        new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
 
  224    side_fe_ptr->getOpPtrVector().push_back(
 
  225        new OpCalculateHTensorTensorField<3, 3>(
 
  227    side_fe_ptr->getOpPtrVector().push_back(
 
  229    side_fe_ptr->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<3>(
 
  231    for (
auto &[name, ents, reaction_vec] : reactionForces) {
 
  239    for (
auto &[name, ents, reaction_vec] : reactionForces) {
 
  240      std::array<double, 6> block_reaction_force{0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
 
  242      MPI_Allreduce(reaction_vec.data(), &block_reaction_force, 6, MPI_DOUBLE,
 
  245      for (
auto &force : block_reaction_force) {
 
  246        if (std::abs(force) < 1e-12) {
 
  252          "Step %d time %3.4g Block %s Reaction force [%3.6e, %3.6e, %3.6e]",
 
  253          ts_step, ts_t, name.c_str(), block_reaction_force[0],
 
  254          block_reaction_force[1], block_reaction_force[2]);
 
  256                  "Step %d time %3.4g  Block %s Moment [%3.6e, %3.6e, %3.6e]",
 
  257                  ts_step, ts_t, name.c_str(), block_reaction_force[3],
 
  258                  block_reaction_force[4], block_reaction_force[5]);
 
  261    auto get_step = [](
auto ts_step) {
 
  262      std::ostringstream ss;
 
  263      ss << boost::str(boost::format(
"%d") % 
static_cast<int>(ts_step));
 
  264      std::string s = ss.str();
 
  270      CHKERR PetscViewerBinaryOpen(
 
  271          PETSC_COMM_WORLD, (
"restart_" + get_step(ts_step) + 
".dat").c_str(),
 
  272          FILE_MODE_WRITE, &viewer);
 
  273      CHKERR VecView(ts_u, viewer);
 
  274      CHKERR PetscViewerDestroy(&viewer);
 
  281    bool post_process_skeleton = 
true;
 
  283      post_process_skeleton = 
true;
 
  285    if (post_process_skeleton) {
 
  287      auto get_material_force_tags = [&]() {
 
  289        std::vector<Tag> tag(2);
 
  299          1, 
"out_skeleton_" + get_step(ts_step) + 
".h5m", PETSC_NULLPTR,
 
  300          get_material_force_tags());
 
  310    double body_energy = 0;
 
  311    MPI_Allreduce(
gEnergy.get(), &body_energy, 1, MPI_DOUBLE, MPI_SUM,
 
  314      auto crack_area_ptr = boost::make_shared<double>(0.0);
 
  317                  "Step %d time %3.4g strain energy %3.6e crack " 
  319                  ts_step, ts_t, body_energy, *crack_area_ptr);
 
  322        if (crack_faces.empty()) {
 
  327            (boost::format(
"crack_faces_step_%d.vtk") % ts_step).str(),
 
  331            (boost::format(
"front_edges_step_%d.vtk") % ts_step).str(),
 
  336      MOFEM_LOG_C(
"EP", Sev::inform, 
"Step %d time %3.4g strain energy %3.6e",
 
  337                  ts_step, ts_t, body_energy);
 
  340    auto post_proc_at_points = [&](std::array<double, 3> point,
 
  344      dataFieldEval->setEvalPoints(point.data(), point.size() / 3);
 
  346      struct OpPrint : 
public VolOp {
 
  349        std::array<double, 3> point;
 
  357        MoFEMErrorCode doWork(
int side, EntityType type,
 
  358                              EntitiesFieldData::EntData &data) {
 
  361            if (getGaussPts().size2()) {
 
  363              auto t_h = getFTensor2FromMat<3, 3>(eP.
dataAtPts->hAtPts);
 
  365                  getFTensor2FromMat<3, 3>(eP.
dataAtPts->approxPAtPts);
 
  370              const double jac = determinantTensor3by3(t_h);
 
  372              t_cauchy(
i, 
j) = (1. / jac) * (t_approx_P(
i, 
k) * t_h(
j, 
k));
 
  375                std::ostringstream s;
 
  376                s << str << 
" elem " << getFEEntityHandle() << 
" ";
 
  380              auto print_tensor = [](
auto &
t) {
 
  381                std::ostringstream s;
 
  386              std::ostringstream print;
 
  390                  << add() << 
"point " << getVectorAdaptor(point.data(), 3);
 
  392                  << add() << 
"coords at gauss pts " << getCoordsAtGaussPts();
 
  394                  << add() << 
"w " << eP.
dataAtPts->wL2AtPts;
 
  396                  << add() << 
"Piola " << eP.
dataAtPts->approxPAtPts;
 
  398                  << add() << 
"Cauchy " << print_tensor(t_cauchy);
 
  407        fe_ptr->getOpPtrVector().push_back(
new OpPrint(
eP, point, str));
 
  409            ->evalFEAtThePoint<SPACE_DIM>(
 
  410                point.data(), 1e-12, problemPtr->getName(), 
"EP", 
dataFieldEval,
 
  413        fe_ptr->getOpPtrVector().pop_back();
 
  419    auto check_external_strain = [&](std::array<double, 3> point,
 
  423      dataFieldEval->setEvalPoints(point.data(), point.size() / 3);
 
  429            ->evalFEAtThePoint<SPACE_DIM>(
 
  430                point.data(), 1e-12, problemPtr->getName(), 
"EP", 
dataFieldEval,
 
  435        VecSetValue(vec, 0, 0.0, ADD_VALUES);
 
  438          std::ostringstream s;
 
  439          s << str << 
" elem " << getFEEntityHandle() << 
" ";
 
  445            << add() << 
"point " << getVectorAdaptor(point.data(), 3);
 
  449        VecSetValue(vec, 0, disp_at_point, ADD_VALUES);
 
  454        CHKERR VecAssemblyBegin(vec);
 
  455        CHKERR VecAssemblyEnd(vec);
 
  459          CHKERR VecGetValues(vec, 1, &idx, &error);
 
  461        MPI_Bcast(&error, 1, MPI_DOUBLE, 0, PETSC_COMM_WORLD);
 
  465        if (std::abs(error - 0.25) > 1e-5) {
 
  467                   "Atom test %d failed: wrong displacement.", 
atom_test);
 
  485        CHKERR post_proc_at_points(pts.second, pts.first);