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;
 
  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");
 
    CHK_THROW_MESSAGE(setElementForFieldEval(), "set element for field");
    CHK_MOAB_THROW(readPtsFromFile(), "read post proc points");
    CHK_THROW_MESSAGE(setElementForPostProcessEnergy(),
                      "set element for post energy");
 
    CHK_THROW_MESSAGE(PetscOptionsGetBool(PETSC_NULLPTR, "", "-write_restart",
                                          &writeRestart, PETSC_NULLPTR),
                      "get write restart option");
    CHK_THROW_MESSAGE(PetscOptionsGetInt(PETSC_NULLPTR, "",
                                         "-write_restart_every",
                                         &writeRestartEvery, PETSC_NULLPTR),
                      "get write restart option");
 
    MOFEM_LOG("EP", Sev::inform) << "EshelbianMonitor writeRestart "
                                 << (writeRestart ? "true" : "false");
 
    MOFEM_LOG("EP", Sev::inform) << "EshelbianMonitor writeRestartEvery "
                                 << writeRestartEvery << " steps";
  }
 
  virtual ~EshelbianMonitor() {
    adjointGradientVector.reset();
  }
 
  MoFEMErrorCode preProcess() { return 0; }
 
  MoFEMErrorCode operator()() { return 0; }
 
  MoFEMErrorCode postProcess() {
    MoFEMFunctionBegin;
 
    MOFEM_LOG("EP", Sev::inform) << "Monitor postProcess";
 
    CHKERR saveRestartFile();
    CHKERR calculateReactionForces();
    CHKERR calculateEnergy();
 
    CHKERR postProcessMaterialForces();
    CHKERR postProcessSkin();
    CHKERR postProcessSkeletonResults();
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode setElementForPostProcessEnergy();
  std::string getStep();
  MoFEMErrorCode calculateReactionForces();
  MoFEMErrorCode calculateEnergy();
  MoFEMErrorCode saveRestartFile();
  MoFEMErrorCode postProcessSkin();
  MoFEMErrorCode postProcessMaterialForces();
  MoFEMErrorCode postProcessSkeletonResults();
  MoFEMErrorCode readPtsFromFile();
  MoFEMErrorCode setElementForFieldEval();
  MoFEMErrorCode postProcAtPoints(std::array<double, 3> point, std::string str);
 
  std::map<std::string, std::array<double, 6>> &getReactionMap() {
    return reactionForcesMap;
  }
 
  EshelbianCore &getEpCore() { return eP; }
 
  boost::shared_ptr<SetPtsData> &getDataField() { return dataFieldEval; }
 
  PtsHashMap &getHashMap() { return ptsHashMap; }
 
  boost::shared_ptr<double> &getEnergy() { return gEnergy; }
 
  inline static SmartPetscObj<Vec> adjointGradientVector = SmartPetscObj<Vec>();
 
protected:
  std::map<std::string, std::array<double, 6>> reactionForcesMap;
  EshelbianCore &eP;
  boost::shared_ptr<SetPtsData> dataFieldEval;
  PtsHashMap ptsHashMap;
  boost::shared_ptr<double> gEnergy;
  boost::shared_ptr<VolEle> volPostProcEnergy;
  PetscBool writeRestart = PETSC_FALSE;
  PetscInt writeRestartEvery = 1;
};
 
MoFEMErrorCode EshelbianMonitor::setElementForPostProcessEnergy() {
  MoFEMFunctionBegin;
  // set time to post processing elements
  volPostProcEnergy->data_ctx |= PetscData::CTX_SET_TIME;
  CHKERR eP.setBaseVolumeElementOps(1, false, false, false, volPostProcEnergy);
 
  if (auto op = eP.physicalEquations->returnOpCalculateEnergy(eP.dataAtPts,
                                                              gEnergy)) {
 
    if (EshelbianCore::stretchHandling == EshelbianCore::NO_STREACH ||
        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);
};
 
std::string EshelbianMonitor::getStep() {
  std::ostringstream ss;
  ss << boost::str(boost::format("%d") % static_cast<int>(ts_step));
  return ss.str();
};
 
MoFEMErrorCode EshelbianMonitor::calculateReactionForces() {
  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(
      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;
      }
    }
    // Store the reaction forces for later access
    reactionForcesMap[name] = block_reaction_force;
    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);
};
 
MoFEMErrorCode EshelbianMonitor::calculateEnergy() {
  MoFEMFunctionBegin;
  // Loop boundary elements with traction boundary conditions
  *gEnergy = 0;
  // set time to post proc energy
  const bool has_ts = (ts != PETSC_NULLPTR);
  if (has_ts) {
    CHKERR TSGetTime(ts, &(volPostProcEnergy->ts_t));
    CHKERR TSGetTimeStep(ts, &(volPostProcEnergy->ts_dt));
  }
  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());
  *gEnergy = body_energy;
  if (EshelbianCore::crackingOn) {
    CHKERR eP.calculateCrackArea(eP.currentCrackAreaPtr);
    // store strain energy
    eP.strainEnergy = body_energy;
    MOFEM_LOG_C("EP", Sev::inform,
                "Step %d time %3.4g strain energy %3.12e crack "
                "area %3.12e",
                ts_step, ts_t, body_energy, *eP.currentCrackAreaPtr);
    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);
};
 
MoFEMErrorCode EshelbianMonitor::saveRestartFile() {
  MoFEMFunctionBegin;
 
  if (writeRestart && ts_u && (ts_step % writeRestartEvery == 0)) {
    PetscViewer viewer;
    CHKERR PetscViewerBinaryOpen(PETSC_COMM_WORLD,
                                 ("restart_" + getStep() + ".dat").c_str(),
                                 FILE_MODE_WRITE, &viewer);
    CHKERR PetscViewerBinarySetSkipInfo(viewer, PETSC_TRUE);
    CHKERR VecView(ts_u, viewer);
    CHKERR PetscViewerDestroy(&viewer);
    CHKERR eP.postProcessRestartMesh(1, "out_restart_mesh_" +
                                            getStep() + ".h5m");
  }
  MoFEMFunctionReturn(0);
};
 
MoFEMErrorCode EshelbianMonitor::postProcessSkin() {
  MoFEMFunctionBegin;
  CHKERR eP.postProcessResults(1, "out_sol_elastic_" + getStep() + ".h5m",
                               PETSC_NULLPTR, PETSC_NULLPTR,
                               adjointGradientVector, {}, ts);
  MoFEMFunctionReturn(0);
};
 
MoFEMErrorCode EshelbianMonitor::postProcessMaterialForces() {
  MoFEMFunctionBegin;
  // Function to get material force tags
  if (!EshelbianCore::crackingOn) {
    int nb_front_faces = eP.frontEdges->size();
    MPI_Allreduce(MPI_IN_PLACE, &nb_front_faces, 1, MPI_INT, MPI_SUM,
                  eP.mField.get_comm());
    if (nb_front_faces > 0) {
      CHKERR eP.calculateFaceMaterialForce(1, ts, &adjointGradientVector);
    }
  }
  MoFEMFunctionReturn(0);
};
 
MoFEMErrorCode EshelbianMonitor::postProcessSkeletonResults() {
  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 = false;
#ifndef NDEBUG
  post_process_skeleton = true;
#endif
  if (post_process_skeleton) {
    CHKERR eP.postProcessSkeletonResults(
        1, "out_skeleton_" + getStep() + ".h5m", PETSC_NULLPTR,
        get_material_force_tags());
  }
  MoFEMFunctionReturn(0);
};
 
MoFEMErrorCode EshelbianMonitor::readPtsFromFile() {
  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;
      ptsHashMap[col1] = {col2, col3, col4};
    } else {
      MOFEM_LOG("EP", Sev::error) << "Error parsing line: " << line;
    }
  }
 
  file.close(); // Close the file
  MoFEMFunctionReturn(0);
};
 
MoFEMErrorCode EshelbianMonitor::setElementForFieldEval() {
  MoFEMFunctionBegin;
  boost::shared_ptr<Ele> vol_ele(dataFieldEval->feMethodPtr);
  auto no_rule = [](int, int, int) { return -1; };
  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(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::stretchHandling == EshelbianCore::NO_STREACH) {
    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);
};
 
MoFEMErrorCode EshelbianMonitor::postProcAtPoints(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) {
        const int nb_gauss_pts = getGaussPts().size2();
        if (nb_gauss_pts) {
 
          auto t_h = eP.dataAtPts->getFTensorSmallH(nb_gauss_pts);
          auto t_approx_P = eP.dataAtPts->getFTensorApproxP(nb_gauss_pts);
 
          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->getSmallWL2AtPts());
          MOFEM_LOG("EPSYNC", Sev::inform)
              << add() << "Piola " << *(eP.dataAtPts->getApproxPAtPts());
          MOFEM_LOG("EPSYNC", Sev::inform)
              << add() << "Cauchy " << print_tensor(t_cauchy);
        }
      }
      MoFEMFunctionReturn(0);
    }
  };
 
  if (auto fe_ptr = dataFieldEval->feMethodPtr) {
 
    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);
}