TSElasticPostStep.cpp

Go to the documentation of this file.

/** @file TSElasticPostStep.cpp
 * @brief
 * @date 2023-05-13
 *
 * @license{This project is released under the MIT License.}
 *
 */
 
struct TSElasticPostStep {
 
  static MoFEMErrorCode postStepInitialise(EshelbianCore *ep_ptr) {
    MoFEMFunctionBegin;
 
    epPtr = ep_ptr;
 
    auto create_post_step_ksp = [&]() {
      auto ksp = createKSP(epPtr->mField.get_comm());
 
      auto set_up = [&]() {
        MoFEMFunctionBegin;
        using DomainEle = VolumeElementForcesAndSourcesCore;
        using DomainEleOp = DomainEle::UserDataOperator;
        using OpDomainMass = FormsIntegrators<DomainEleOp>::Assembly<
            PETSC>::BiLinearForm<GAUSS>::OpMass<1, 3>;
        using OpRhs = FormsIntegrators<DomainEleOp>::Assembly<
            PETSC>::LinearForm<GAUSS>::OpBaseTimesVector<1, 3, 1>;
        auto fe_lhs = boost::make_shared<DomainEle>(ep_ptr->mField);
        auto fe_rhs = boost::make_shared<DomainEle>(ep_ptr->mField);
 
        fe_lhs->getUserPolynomialBase() =
            boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
        fe_rhs->getUserPolynomialBase() =
            boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
        CHKERR
        EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
            fe_lhs->getOpPtrVector(), {L2}, ep_ptr->materialH1Positions,
            ep_ptr->frontAdjEdges);
        CHKERR
        EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
            fe_rhs->getOpPtrVector(), {L2}, ep_ptr->materialH1Positions,
            ep_ptr->frontAdjEdges);
 
        fe_lhs->getOpPtrVector().push_back(
            new OpDomainMass(ep_ptr->spatialH1Disp, ep_ptr->spatialH1Disp));
        auto w_ptr = boost::make_shared<MatrixDouble>();
        fe_rhs->getOpPtrVector().push_back(
            new OpCalculateVectorFieldValues<3>(ep_ptr->spatialL2Disp, w_ptr));
        fe_rhs->getOpPtrVector().push_back(
            new OpRhs(ep_ptr->spatialH1Disp, w_ptr));
 
        CHKERR DMMoFEMKSPSetComputeOperators(ep_ptr->dmPrjSpatial,
                                             ep_ptr->elementVolumeName, fe_lhs,
                                             nullptr, nullptr);
        CHKERR DMMoFEMKSPSetComputeRHS(ep_ptr->dmPrjSpatial,
                                       ep_ptr->elementVolumeName, fe_rhs,
                                       nullptr, nullptr);
 
        // preProcRhs = boost::make_shared<FEMethod>();
        // struct MinusOne : public ScalingMethod {
        //   double getScale(const double time) { return -time; }
        // };
        // preProcRhs->preProcessHook = EssentialPreProc<DisplacementCubitBcData>(
        //     ep_ptr->mField, preProcRhs, {boost::make_shared<MinusOne>()});
 
        CHKERR KSPAppendOptionsPrefix(ksp, "prjspatial_");
        CHKERR KSPSetFromOptions(ksp);
        CHKERR KSPSetDM(ksp, ep_ptr->dmPrjSpatial);
        CHKERR KSPSetUp(ksp);
        MoFEMFunctionReturn(0);
      };
 
      CHK_THROW_MESSAGE(set_up(), "set up");
 
      return ksp;
    };
 
    prjKsp = create_post_step_ksp();
    prjD = createDMVector(ep_ptr->dmPrjSpatial);
    prjF = vectorDuplicate(prjD);
    prjDM = ep_ptr->dmPrjSpatial;
 
    MoFEMFunctionReturn(0);
  };
 
  static MoFEMErrorCode postStepDestroy() {
    MoFEMFunctionBegin;
    prjKsp.reset();
    prjD.reset();
    prjF.reset();
    prjDM.reset();
    preProcRhs.reset();
    MoFEMFunctionReturn(0);
  };
 
  static MoFEMErrorCode preStepFun(TS ts) {
    MoFEMFunctionBegin;
    MOFEM_LOG("EP", Sev::inform) << "Pre step";
 
    double time;
    CHKERR TSGetTime(ts, &time);
 
    if (EshelbianCore::crackingOn && EshelbianCore::crackingStartTime < time) {
 
      auto debug_crack = [&]() {
        MoFEMFunctionBegin;
 
        PetscBool debug_crack_mesh = PETSC_FALSE;
        CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-debug_crack_mesh",
                                   &debug_crack_mesh, PETSC_NULLPTR);
        if (debug_crack_mesh) {
 
          auto get_meshsets_from_block =
              [](MoFEM::Interface &m_field, std::string block_name, int dim,
                 std::vector<EntityHandle> r) {
 
                auto mesh_mng = m_field.getInterface<MeshsetsManager>();
                auto bcs = mesh_mng->getCubitMeshsetPtr(
 
                    std::regex((boost::format("%s(.*)") % block_name).str())
 
                );
 
                for (auto bc : bcs) {
                  r.push_back(bc->getMeshset());
                }
 
                return r;
              };
 
          if (epPtr->mField.get_comm_rank() == 0) {
 
            if (epPtr->maxMovedFaces) {
              CHKERR epPtr->mField.getInterface<MeshsetsManager>()
                  ->addEntitiesToMeshset(BLOCKSET, epPtr->addCrackMeshsetId,
                                         *epPtr->maxMovedFaces);
            }
 
            auto meshset_ptr = get_temp_meshset_ptr(epPtr->mField.get_moab());
            Range tets;
            CHKERR epPtr->mField.get_moab().get_entities_by_dimension(
                *meshset_ptr, 3, tets);
            CHKERR epPtr->mField.get_moab().add_entities(*meshset_ptr, tets);
 
            std::vector<EntityHandle> meshsets;
            meshsets.push_back(*meshset_ptr);
 
            meshsets =
                get_meshsets_from_block(epPtr->mField, "CRACK", 2, meshsets);
            meshsets =
                get_meshsets_from_block(epPtr->mField, "FRONT", 1, meshsets);
            meshsets =
                get_meshsets_from_block(epPtr->mField, "EDGE", 1, meshsets);
 
            int time_step = 0;
            CHKERR TSGetStepNumber(ts, &time_step);
            std::string file_name =
                "crack_meshsets_" + std::to_string(time_step) + ".h5m";
            CHKERR epPtr->mField.get_moab().write_file(
                file_name.c_str(), "MOAB", nullptr, meshsets.data());
          }
        }
 
        MoFEMFunctionReturn(0);
      };
 
      CHKERR epPtr->setNewFrontCoordinates();
      CHKERR debug_crack();
      CHKERR epPtr->addCrackSurfaces();
      CHKERR epPtr->projectGeometry();
    }
 
    Vec T;
    CHKERR TSGetSolution(ts, &T);
 
    auto zero_filled = [&](auto name) {
      MoFEMFunctionBegin;
      auto is_mng = epPtr->mField.getInterface<ISManager>();
      SmartPetscObj<IS> is;
      CHKERR is_mng->isCreateProblemFieldAndRankLocal(
          "ELASTIC_PROBLEM", ROW, name, 0, MAX_DOFS_ON_ENTITY, is);
      const int *index_ptr;
      CHKERR ISGetIndices(is, &index_ptr);
      int size;
      CHKERR ISGetLocalSize(is, &size);
      double *a;
      CHKERR VecGetArray(T, &a);
      for (auto i = 0; i != size; i++) {
        a[index_ptr[i]] = 0;
      }
      CHKERR VecRestoreArray(T, &a);
      CHKERR ISRestoreIndices(is, &index_ptr);
      MoFEMFunctionReturn(0);
    };
 
    switch (EshelbianCore::gradApproximator) {
    case NO_H1_CONFIGURATION:
      break;
    case LARGE_ROT:
      MOFEM_LOG("EP", Sev::verbose) << "Zeroing (large) fields";
      CHKERR zero_filled(epPtr->stretchTensor);
      CHKERR zero_filled(epPtr->rotAxis);
      break;
    case MODERATE_ROT:
      MOFEM_LOG("EP", Sev::verbose) << "Zeroing (moderate) fields";
      CHKERR zero_filled(epPtr->stretchTensor);
      CHKERR zero_filled(epPtr->rotAxis);
      break;
    case SMALL_ROT:
      MOFEM_LOG("EP", Sev::verbose) << "Zeroing (small) fields";
      CHKERR zero_filled(epPtr->stretchTensor);
      CHKERR zero_filled(epPtr->rotAxis);
      break;
    }
 
    CHKERR VecCopy(T, epPtr->solTSStep);
    CHKERR VecGhostUpdateBegin(epPtr->solTSStep, INSERT_VALUES,
                               SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(epPtr->solTSStep, INSERT_VALUES, SCATTER_FORWARD);
 
    MoFEMFunctionReturn(0);
  }
 
  static MoFEMErrorCode postStepFun(TS ts) {
    MoFEMFunctionBegin;
 
    double time;
    CHKERR TSGetTime(ts, &time);
 
    MOFEM_LOG("EP", Sev::inform) << "Solve H1 post-step";
    CHKERR VecZeroEntries(prjF);
    CHKERR KSPSolve(prjKsp, prjF, prjD);
    CHKERR VecGhostUpdateBegin(prjD, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(prjD, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR DMoFEMMeshToLocalVector(prjDM, prjD, INSERT_VALUES, SCATTER_REVERSE);
 
    if (EshelbianCore::crackingOn) {
      switch (EshelbianCore::energyReleaseSelector) {
      case GRIFFITH_FORCE:
      case GRIFFITH_SKELETON:
        MOFEM_LOG("EP", Sev::inform) << "Calculate Griffith force";
        CHKERR epPtr->calculateFaceMaterialForce(1, ts);
        break;
      default:
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Energy release selector not handled");
      };
      if (EshelbianCore::crackingStartTime < time) {
        MOFEM_LOG("EP", Sev::inform) << "Calculate faces orientation";
 
        CHKERR epPtr->calculateOrientation(1, true);
      }
    }
 
 
    // double time;
    // CHKERR TSGetTime(ts, &time);
    // preProcRhs->ts_t = time;
    // CHKERR DMoFEMPreProcessFiniteElements(prjDM, preProcRhs.get());
    CHKERR VecZeroEntries(epPtr->solTSStep);
 
    MoFEMFunctionReturn(0);
  };
 
private:
  static SmartPetscObj<KSP> prjKsp;
  static SmartPetscObj<Vec> prjD;
  static SmartPetscObj<Vec> prjF;
  static SmartPetscObj<DM> prjDM;
  static EshelbianCore *epPtr;
 
  static boost::shared_ptr<FEMethod> preProcRhs;
};
 
SmartPetscObj<KSP> TSElasticPostStep::prjKsp;
SmartPetscObj<Vec> TSElasticPostStep::prjD;
SmartPetscObj<Vec> TSElasticPostStep::prjF;
SmartPetscObj<DM> TSElasticPostStep::prjDM;
EshelbianCore *TSElasticPostStep::epPtr;
boost::shared_ptr<FEMethod> TSElasticPostStep::preProcRhs;