EshelbianMonitor.cpp

Go to the documentation of this file.

/** @file
                @brief Contains definition of EshelbianMonitor class.
                @ingroup EshelbianPlasticty
*/
 
struct EshelbianMonitor : public FEMethod {
 
  using Ele = ForcesAndSourcesCore;
  using VolEle = VolumeElementForcesAndSourcesCore;
  using VolOp = VolumeElementForcesAndSourcesCore::UserDataOperator;
  using SetPtsData = FieldEvaluatorInterface::SetPtsData;
  using PtsHashMap = std::map<std::string, std::array<double, 3>>;
  using FaceSideEle =
      PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
 
  EshelbianCore &eP;
  boost::shared_ptr<SetPtsData> dataFieldEval;
  boost::shared_ptr<VolEle> volPostProcEnergy;
  boost::shared_ptr<double> gEnergy;
  PtsHashMap ptsHashMap;
  PetscBool writeRestart = PETSC_FALSE;
 
  EshelbianMonitor(EshelbianCore &ep)
      : eP(ep), dataFieldEval(ep.mField.getInterface<FieldEvaluatorInterface>()
                                  ->getData<VolEle>()),
        volPostProcEnergy(new VolEle(ep.mField)), gEnergy(new double) {
    CHK_THROW_MESSAGE(
        ep.mField.getInterface<FieldEvaluatorInterface>()->buildTree<SPACE_DIM>(
            dataFieldEval, "EP"),
        "build field evaluator tree");
 
    auto no_rule = [](int, int, int) { return -1; };
 
    auto set_element_for_field_eval = [&]() {
      MoFEMFunctionBegin;
      boost::shared_ptr<Ele> vol_ele(dataFieldEval->feMethodPtr.lock());
      vol_ele->getRuleHook = no_rule;
      vol_ele->getUserPolynomialBase() =
          boost::make_shared<CGGUserPolynomialBase>();
      EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
          vol_ele->getOpPtrVector(), {HDIV, H1, L2}, eP.materialH1Positions,
          ep.frontAdjEdges);
 
      auto piola_scale_ptr = boost::make_shared<double>(1.0);
      vol_ele->getOpPtrVector().push_back(
          eP.physicalEquations->returnOpSetScale(piola_scale_ptr,
                                                 eP.physicalEquations));
      vol_ele->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
          eP.piolaStress, eP.dataAtPts->getApproxPAtPts(), piola_scale_ptr));
      vol_ele->getOpPtrVector().push_back(
          new OpCalculateHTensorTensorField<3, 3>(
              eP.bubbleField, eP.dataAtPts->getApproxPAtPts(), piola_scale_ptr,
              SmartPetscObj<Vec>(), MBMAXTYPE));
      if (EshelbianCore::noStretch) {
        vol_ele->getOpPtrVector().push_back(
            eP.physicalEquations->returnOpCalculateStretchFromStress(
                eP.dataAtPts, eP.physicalEquations));
      } else {
        vol_ele->getOpPtrVector().push_back(
            new OpCalculateTensor2SymmetricFieldValues<3>(
                eP.stretchTensor, eP.dataAtPts->getLogStretchTensorAtPts(),
                MBTET));
      }
      vol_ele->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
          eP.rotAxis, eP.dataAtPts->getRotAxisAtPts(), MBTET));
      vol_ele->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
          eP.rotAxis, eP.dataAtPts->getRotAxis0AtPts(), eP.solTSStep, MBTET));
      vol_ele->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
          eP.spatialL2Disp, eP.dataAtPts->getSmallWL2AtPts(), MBTET));
 
      // H1 displacements
      vol_ele->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
          eP.spatialH1Disp, eP.dataAtPts->getSmallWH1AtPts()));
      vol_ele->getOpPtrVector().push_back(
          new OpCalculateVectorFieldGradient<3, 3>(
              eP.spatialH1Disp, eP.dataAtPts->getSmallWGradH1AtPts()));
 
      vol_ele->getOpPtrVector().push_back(
          new OpCalculateRotationAndSpatialGradient(eP.dataAtPts));
      MoFEMFunctionReturn(0);
    };
 
    auto set_element_for_post_process_energy = [&]() {
      MoFEMFunctionBegin;
 
      CHKERR eP.setBaseVolumeElementOps(
          1, false, false, false, SmartPetscObj<Vec>(), volPostProcEnergy);
 
      if (auto op = eP.physicalEquations->returnOpCalculateEnergy(eP.dataAtPts,
                                                                  gEnergy)) {
 
        if (eP.noStretch ||
            eP.materialModel == EshelbianCore::MaterialModel::Hencky) {
          // We have to use actual strains to evaluate J integral and energy,
          // in this case. Note actual stresses, and actual energy can only
          // drive crack growth
 
          // Note: Calling this bellow we overwrite the stretch calculation
          // from setBaseVolumeElementOps 
 
          volPostProcEnergy->getOpPtrVector().push_back(
              eP.physicalEquations->returnOpCalculateStretchFromStress(
                  eP.dataAtPts, eP.physicalEquations));
        }
 
        volPostProcEnergy->getOpPtrVector().push_back(op);
      }
 
      MoFEMFunctionReturn(0);
    };
 
    auto reads_post_proc_data = [](PtsHashMap &pts_hash_map) {
      MoFEMFunctionBegin;
      std::ifstream file(
          "points.txt"); // Open the file with the name "data.txt"
 
      if (!file.is_open()) {
        MoFEMFunctionReturnHot(0);
      }
 
      std::string line;
 
      while (std::getline(file, line)) {
        std::istringstream iss(line);
        std::string col1;
        double col2, col3, col4;
 
        if (iss >> col1 >> col2 >> col3 >> col4) {
          MOFEM_LOG("EP", Sev::verbose) << "Read: " << col1 << ", " << col2
                                        << ", " << col3 << ", " << col4;
          pts_hash_map[col1] = {col2, col3, col4};
        } else {
          MOFEM_LOG("EP", Sev::error) << "Error parsing line: " << line;
        }
      }
 
      file.close(); // Close the file
      MoFEMFunctionReturn(0);
    };
 
    CHK_THROW_MESSAGE(set_element_for_field_eval(), "set element for field");
    CHK_THROW_MESSAGE(set_element_for_post_process_energy(),
                      "set element for post energy");
 
    PetscBool test_cook_flg = PETSC_FALSE;
    CHK_THROW_MESSAGE(PetscOptionsGetBool(PETSC_NULLPTR, "", "-test_cook_pts",
                                          &test_cook_flg, PETSC_NULLPTR),
                      "get post proc points");
    if (test_cook_flg) {
      ptsHashMap["Point A"] = {48., 60., 4.999};
      ptsHashMap["Point B"] = {48. / 2., 44. + (60. - 44.) / 2., 0.};
      ptsHashMap["Point C"] = {48. / 2., (44. - 0.) / 2., 0.};
    }
 
    PetscInt atom_test = 0;
 
    CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-atom_test", &atom_test,
                              PETSC_NULLPTR);
    if (atom_test == 14) {
      // Points for atom test 14: check external strain
      ptsHashMap["Point (2.5, 0., 0.)"] = {2.5, 0, 0.};
    }
    CHK_MOAB_THROW(reads_post_proc_data(ptsHashMap), "read post proc points");
 
    CHK_THROW_MESSAGE(PetscOptionsGetBool(PETSC_NULLPTR, "", "-write_restart",
                                          &writeRestart, PETSC_NULLPTR),
                      "get write restart option");
 
    MOFEM_LOG("EP", Sev::inform)
            << "EshelbianMonitor writeRestart " << writeRestart
        ? " true"
        : " false";
  }
 
  MoFEMErrorCode preProcess() { return 0; }
 
  MoFEMErrorCode operator()() { return 0; }
 
  MoFEMErrorCode postProcess() {
    MoFEMFunctionBegin;
 
    MOFEM_LOG("EP", Sev::inform) << "Monitor postProcess";
 
    auto get_step = [](auto ts_step) {
      std::ostringstream ss;
      ss << boost::str(boost::format("%d") % static_cast<int>(ts_step));
      std::string s = ss.str();
      return s;
    };
 
    auto calculate_reaction_forces = [&]() {
      MoFEMFunctionBegin;
 
      // Get boundary faces marked in blocks name "SPATIAL_DISP_...".
      auto get_ents_on_mesh_skin = [&]() {
        std::map<std::string, Range> boundary_entities_vec;
 
        auto get_block_vec_impl = [&](auto block_name) {
          return eP.mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
              std::regex(
 
                  (boost::format("%s(.*)") % block_name).str()
 
                      ));
        };
 
        auto add_blockset_ents_impl = [&](auto &&vec) {
          MoFEMFunctionBegin;
          for (auto it : vec) {
            Range boundary_entities;
            CHKERR it->getMeshsetIdEntitiesByDimension(eP.mField.get_moab(), 2,
                                                       boundary_entities, true);
            std::string meshset_name = it->getName();
            boundary_entities_vec[meshset_name] = boundary_entities;
            boundary_entities.clear();
          }
          MoFEMFunctionReturn(0);
        };
 
        for (auto block_name : {"SPATIAL_DISP", "FIX", "CONTACT",
                                "SPATIAL_ROTATION", "NORMAL_DISPLACEMENT"}) {
          CHK_THROW_MESSAGE(
              add_blockset_ents_impl(get_block_vec_impl(block_name)),
              "add blockset entities");
        }
 
        return boundary_entities_vec;
      };
 
      auto boundary_entities_vec = get_ents_on_mesh_skin();
 
      std::vector<std::tuple<std::string, Range, std::array<double, 6>>>
          reactionForces;
      for (const auto &pair : boundary_entities_vec) {
        reactionForces.push_back(std::make_tuple(
            pair.first, pair.second,
            std::array<double, 6>{0.0, 0.0, 0.0, 0.0, 0.0, 0.0}));
      }
 
      auto integration_rule_face = [](int, int, int approx_order) {
        return 2 * approx_order + 1;
      };
      auto face_fe =
          boost::make_shared<FaceElementForcesAndSourcesCore>(eP.mField);
      auto no_rule = [](int, int, int) { return -1; };
      face_fe->getRuleHook = integration_rule_face;
      CHKERR
      EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
          face_fe->getOpPtrVector(), {L2}, eP.materialH1Positions,
          eP.frontAdjEdges);
 
      auto op_side =
          new OpLoopSide<FaceSideEle>(eP.mField, "EP", SPACE_DIM, Sev::noisy);
      face_fe->getOpPtrVector().push_back(op_side);
      auto side_fe_ptr = op_side->getSideFEPtr();
      auto base_ptr =
          boost::make_shared<EshelbianPlasticity::CGGUserPolynomialBase>();
      side_fe_ptr->getUserPolynomialBase() = base_ptr;
      CHKERR
      EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
          side_fe_ptr->getOpPtrVector(), {H1, HDIV, L2}, eP.materialH1Positions,
          eP.frontAdjEdges);
      auto piola_scale_ptr = boost::make_shared<double>(1.0);
      side_fe_ptr->getOpPtrVector().push_back(
          eP.physicalEquations->returnOpSetScale(piola_scale_ptr,
                                                 eP.physicalEquations));
      side_fe_ptr->getOpPtrVector().push_back(
          new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
              eP.piolaStress, eP.dataAtPts->getApproxPAtPts(),
              piola_scale_ptr));
      side_fe_ptr->getOpPtrVector().push_back(
          new OpCalculateHTensorTensorField<3, 3>(
              eP.bubbleField, eP.dataAtPts->getApproxPAtPts(), MBMAXTYPE));
      side_fe_ptr->getOpPtrVector().push_back(
          new OpCalculateTractionFromSideEl(eP.dataAtPts));
      side_fe_ptr->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<3>(
              eP.spatialL2Disp, eP.dataAtPts->getSmallWL2AtPts(), MBTET));
      for (auto &[name, ents, reaction_vec] : reactionForces) {
        face_fe->getOpPtrVector().push_back(new OpCalculateReactionForces(
            eP.dataAtPts, name, ents, reaction_vec));
      }
 
      CHKERR eP.mField.loop_finite_elements(problemPtr->getName(),
                                            eP.skinElement, *face_fe);
 
      for (auto &[name, ents, reaction_vec] : reactionForces) {
        std::array<double, 6> block_reaction_force{0.0, 0.0, 0.0,
                                                   0.0, 0.0, 0.0};
 
        MPI_Allreduce(reaction_vec.data(), &block_reaction_force, 6, MPI_DOUBLE,
                      MPI_SUM, eP.mField.get_comm());
 
        for (auto &force : block_reaction_force) {
          if (std::abs(force) < 1e-12) {
            force = 0.0;
          }
        }
        MOFEM_LOG_C("EP", Sev::inform,
                    "Step %d time %3.4g Block %s Reaction force [%3.6e, "
                    "%3.6e, %3.6e]",
                    ts_step, ts_t, name.c_str(), block_reaction_force[0],
                    block_reaction_force[1], block_reaction_force[2]);
        MOFEM_LOG_C("EP", Sev::inform,
                    "Step %d time %3.4g  Block %s Moment [%3.6e, %3.6e, %3.6e]",
                    ts_step, ts_t, name.c_str(), block_reaction_force[3],
                    block_reaction_force[4], block_reaction_force[5]);
      }
      MoFEMFunctionReturn(0);
    };
 
    auto save_restart_file = [&]() {
      MoFEMFunctionBegin;
 
      if (writeRestart && ts_u) {
        PetscViewer viewer;
        CHKERR PetscViewerBinaryOpen(
            PETSC_COMM_WORLD, ("restart_" + get_step(ts_step) + ".dat").c_str(),
            FILE_MODE_WRITE, &viewer);
        CHKERR PetscViewerBinarySetSkipInfo(viewer, PETSC_TRUE);
        CHKERR VecView(ts_u, viewer);
        CHKERR PetscViewerDestroy(&viewer);
      }
      MoFEMFunctionReturn(0);
    };
 
    auto post_process_skin = [&]() {
      MoFEMFunctionBegin;
      CHKERR eP.postProcessResults(1, "out_sol_elastic_" + get_step(ts_step) +
                                          ".h5m");
      MoFEMFunctionReturn(0);
    };
 
    auto post_process_material_forces = [&]() {
      MoFEMFunctionBegin;
      // Function to get material force tags
      CHKERR eP.calculateFaceMaterialForce(1, ts);
      MoFEMFunctionReturn(0);
    };
 
    auto post_process_skeleton_results = [&]() {
      MoFEMFunctionBegin;
      auto get_material_force_tags = [&]() {
        auto &moab = eP.mField.get_moab();
        std::vector<Tag> tag(2);
        CHK_MOAB_THROW(moab.tag_get_handle("MaterialForce", tag[0]),
                       "can't get tag");
        CHK_MOAB_THROW(moab.tag_get_handle("FacePressure", tag[1]),
                       "can't get tag");
        return tag;
      };
      // Post-process skeleton elements
      bool post_process_skeleton = true;
#ifndef NDEBUG
      post_process_skeleton = true;
#endif
      if (post_process_skeleton) {
        CHKERR eP.postProcessSkeletonResults(
            1, "out_skeleton_" + get_step(ts_step) + ".h5m", PETSC_NULLPTR,
            get_material_force_tags());
      }
      MoFEMFunctionReturn(0);
    };
 
    auto calculate_energy = [&]() {
      MoFEMFunctionBegin;
      // Loop boundary elements with traction boundary conditions
      *gEnergy = 0;
      TetPolynomialBase::switchCacheBaseOn<HDIV>({volPostProcEnergy.get()});
      CHKERR eP.mField.loop_finite_elements(problemPtr->getName(), "EP",
                                            *volPostProcEnergy);
      TetPolynomialBase::switchCacheBaseOff<HDIV>({volPostProcEnergy.get()});
 
      double body_energy = 0;
      MPI_Allreduce(gEnergy.get(), &body_energy, 1, MPI_DOUBLE, MPI_SUM,
                    eP.mField.get_comm());
      if (EshelbianCore::crackingOn) {
        auto crack_area_ptr = boost::make_shared<double>(0.0);
        CHKERR eP.calculateCrackArea(crack_area_ptr);
        MOFEM_LOG_C("EP", Sev::inform,
                    "Step %d time %3.4g strain energy %3.6e crack "
                    "area %3.6e",
                    ts_step, ts_t, body_energy, *crack_area_ptr);
        if (eP.mField.get_comm_rank() == 0) {
          auto crack_faces = *eP.crackFaces;
          if (crack_faces.empty()) {
            crack_faces = *eP.frontEdges;
          }
          CHKERR save_range(
              eP.mField.get_moab(),
              (boost::format("crack_faces_step_%d.vtk") % ts_step).str(),
              crack_faces);
          CHKERR save_range(
              eP.mField.get_moab(),
              (boost::format("front_edges_step_%d.vtk") % ts_step).str(),
              *eP.frontEdges);
        }
 
      } else {
        MOFEM_LOG_C("EP", Sev::inform, "Step %d time %3.4g strain energy %3.6e",
                    ts_step, ts_t, body_energy);
      }
      MoFEMFunctionReturn(0);
    };
 
    auto post_proc_at_points = [&](std::array<double, 3> point,
                                   std::string str) {
      MoFEMFunctionBegin;
 
      dataFieldEval->setEvalPoints(point.data(), point.size() / 3);
 
      struct OpPrint : public VolOp {
 
        EshelbianCore &eP;
        std::array<double, 3> point;
        std::string str;
 
        OpPrint(EshelbianCore &ep, std::array<double, 3> &point,
                std::string &str)
            : VolOp(ep.spatialL2Disp, VolOp::OPROW), eP(ep), point(point),
              str(str) {}
 
        MoFEMErrorCode doWork(int side, EntityType type,
                              EntitiesFieldData::EntData &data) {
          MoFEMFunctionBegin;
          if (type == MBTET) {
            if (getGaussPts().size2()) {
 
              auto t_h = getFTensor2FromMat<3, 3>(eP.dataAtPts->hAtPts);
              auto t_approx_P =
                  getFTensor2FromMat<3, 3>(eP.dataAtPts->approxPAtPts);
 
              FTensor::Index<'i', 3> i;
              FTensor::Index<'j', 3> j;
              FTensor::Index<'k', 3> k;
              const double jac = determinantTensor3by3(t_h);
              FTensor::Tensor2<double, 3, 3> t_cauchy;
              t_cauchy(i, j) = (1. / jac) * (t_approx_P(i, k) * t_h(j, k));
 
              auto add = [&]() {
                std::ostringstream s;
                s << str << " elem " << getFEEntityHandle() << " ";
                return s.str();
              };
 
              auto print_tensor = [](auto &t) {
                std::ostringstream s;
                s << t;
                return s.str();
              };
 
              std::ostringstream print;
              MOFEM_LOG("EPSYNC", Sev::inform)
                  << add() << "comm rank " << eP.mField.get_comm_rank();
              MOFEM_LOG("EPSYNC", Sev::inform)
                  << add() << "point " << getVectorAdaptor(point.data(), 3);
              MOFEM_LOG("EPSYNC", Sev::inform)
                  << add() << "coords at gauss pts " << getCoordsAtGaussPts();
              MOFEM_LOG("EPSYNC", Sev::inform)
                  << add() << "w " << eP.dataAtPts->wL2AtPts;
              MOFEM_LOG("EPSYNC", Sev::inform)
                  << add() << "Piola " << eP.dataAtPts->approxPAtPts;
              MOFEM_LOG("EPSYNC", Sev::inform)
                  << add() << "Cauchy " << print_tensor(t_cauchy);
            }
          }
          MoFEMFunctionReturn(0);
        }
      };
 
      if (auto fe_ptr = dataFieldEval->feMethodPtr.lock()) {
 
        fe_ptr->getOpPtrVector().push_back(new OpPrint(eP, point, str));
        CHKERR eP.mField.getInterface<FieldEvaluatorInterface>()
            ->evalFEAtThePoint<SPACE_DIM>(
                point.data(), 1e-12, problemPtr->getName(), "EP", dataFieldEval,
                eP.mField.get_comm_rank(), eP.mField.get_comm_rank(), nullptr,
                MF_EXIST, QUIET);
        fe_ptr->getOpPtrVector().pop_back();
      }
 
      MoFEMFunctionReturn(0);
    };
 
    auto check_external_strain = [&](std::array<double, 3> point,
                                     std::string str, PetscInt atom_test) {
      MoFEMFunctionBegin;
 
      dataFieldEval->setEvalPoints(point.data(), point.size() / 3);
      auto vec = createVectorMPI(eP.mField.get_comm(), 1, 1);
      eP.dataAtPts->wL2AtPts.resize(0, 0, false);
 
      if (auto fe_ptr = dataFieldEval->feMethodPtr.lock()) {
        CHKERR eP.mField.getInterface<FieldEvaluatorInterface>()
            ->evalFEAtThePoint<SPACE_DIM>(
                point.data(), 1e-12, problemPtr->getName(), "EP", dataFieldEval,
                eP.mField.get_comm_rank(), eP.mField.get_comm_rank(), nullptr,
                MF_EXIST, QUIET);
      }
      if (eP.dataAtPts->wL2AtPts.size2() == 0) {
        VecSetValue(vec, 0, 0.0, ADD_VALUES);
      } else if (eP.dataAtPts->wL2AtPts.size2() == 1) {
        auto add = [&]() {
          std::ostringstream s;
          s << str << " elem " << getFEEntityHandle() << " ";
          return s.str();
        };
        MOFEM_LOG("EPSYNC", Sev::inform)
            << add() << "comm rank " << eP.mField.get_comm_rank();
        MOFEM_LOG("EPSYNC", Sev::inform)
            << add() << "point " << getVectorAdaptor(point.data(), 3);
        MOFEM_LOG("EPSYNC", Sev::inform)
            << add() << "w " << eP.dataAtPts->wL2AtPts;
        double disp_at_point = eP.dataAtPts->wL2AtPts(0);
        VecSetValue(vec, 0, disp_at_point, ADD_VALUES);
      }
 
      MOFEM_LOG_SEVERITY_SYNC(eP.mField.get_comm(), Sev::inform);
      if (ts_t > 0.0) {
        CHKERR VecAssemblyBegin(vec);
        CHKERR VecAssemblyEnd(vec);
        double error;
        PetscInt idx = 0;
        if (eP.mField.get_comm_rank() == 0) {
          CHKERR VecGetValues(vec, 1, &idx, &error);
        }
        MPI_Bcast(&error, 1, MPI_DOUBLE, 0, PETSC_COMM_WORLD);
        // Check if the displacement is correct for the applied external
        // strain For the bar problem, we expect a displacement of 0.25 at
        // point A
        if (std::abs(error - 0.25) > 1e-5) {
          SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                  "Atom test %d failed: wrong displacement.", atom_test);
        }
      }
      MoFEMFunctionReturn(0);
    };
 
    CHKERR save_restart_file();
    CHKERR calculate_reaction_forces();
    CHKERR calculate_energy();
 
    CHKERR post_process_material_forces();
    CHKERR post_process_skin();
    CHKERR post_process_skeleton_results();
 
    PetscInt atom_test = 0;
    CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-atom_test", &atom_test,
                              PETSC_NULLPTR);
 
    if (atom_test == 14) {
      // Points for external strain
      for (auto &pts : ptsHashMap) {
        CHKERR check_external_strain(pts.second, pts.first, atom_test);
      }
    } else {
      // Points for Cook beam
      for (auto &pts : ptsHashMap) {
        CHKERR post_proc_at_points(pts.second, pts.first);
        MOFEM_LOG_SEVERITY_SYNC(eP.mField.get_comm(), Sev::inform);
      }
    }
 
    MoFEMFunctionReturn(0);
  }
};