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)) {
        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");
  }
 
  MoFEMErrorCode preProcess() { return 0; }
 
  MoFEMErrorCode operator()() { return 0; }
 
  MoFEMErrorCode postProcess() {
    MoFEMFunctionBegin;
 
    MOFEM_LOG("EP", Sev::inform) << "Monitor postProcess";
 
    // Get boundary faces marked in blocks name "SPATIAL_DISP_...".
    auto get_ents_on_mesh_skin = [&]() {
      
      std::map<std::string, Range> boundary_entities_vec;
      Range boundary_entities;
      for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(eP.mField, BLOCKSET, it)) {
        std::string meshset_name = it->getName();
        std::string block_name_disp = "SPATIAL_DISP";
        std::string block_name_fix = "FIX";
        if (meshset_name.compare(0, block_name_disp.size(), block_name_disp) == 0 ||
            meshset_name.compare(0, block_name_fix.size(), block_name_fix) == 0) {
          CHKERR it->getMeshsetIdEntitiesByDimension(eP.mField.get_moab(), 2,
                                                     boundary_entities, true);
          // Add vertices to boundary entities
          boundary_entities_vec[meshset_name] = boundary_entities;
          boundary_entities.clear();
        }
      }
      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]);
    }
 
    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;
    };
 
    if (writeRestart && ts_u) {
      PetscViewer viewer;
      CHKERR PetscViewerBinaryOpen(
          PETSC_COMM_WORLD, ("restart_" + get_step(ts_step) + ".dat").c_str(),
          FILE_MODE_WRITE, &viewer);
      CHKERR VecView(ts_u, viewer);
      CHKERR PetscViewerDestroy(&viewer);
    }
 
    CHKERR eP.postProcessResults(1, "out_sol_elastic_" + get_step(ts_step) +
                                        ".h5m");
 
    // Post-process skeleton elements
    bool post_process_skeleton = false;
    #ifndef NDEBUG
      post_process_skeleton = true;
    #endif
    if (post_process_skeleton) {
 
      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;
      };
 
      CHKERR eP.calculateFaceMaterialForce(1, ts);
      CHKERR eP.postProcessSkeletonResults(
          1, "out_skeleton_" + get_step(ts_step) + ".h5m", PETSC_NULLPTR,
          get_material_force_tags());
    }
 
    // 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);
    }
 
    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);
    };
 
    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);
      }
    }
 
    if (eP.crackingOn) {
    }
 
    MoFEMFunctionReturn(0);
  }
};