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.}
 *
 */
 
#define SINGULARITY
#include <MoFEM.hpp>
using namespace MoFEM;
 
#include <EshelbianPlasticity.hpp>
 
#include <boost/math/constants/constants.hpp>
#include <boost/math/special_functions/lambert_w.hpp>
 
#include <TSElasticPostStep.hpp>
 
namespace EshelbianPlasticity {
 
MoFEMErrorCode TSElasticPostStep::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);
};
 
MoFEMErrorCode TSElasticPostStep::postStepDestroy() {
  MoFEMFunctionBegin;
  prjKsp.reset();
  prjD.reset();
  prjF.reset();
  prjDM.reset();
  preProcRhs.reset();
  MoFEMFunctionReturn(0);
};
 
MoFEMErrorCode TSElasticPostStep::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();
    CHKERR epPtr->projectInternalStress();
  }
 
  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);
}
 
MoFEMErrorCode TSElasticPostStep::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);
    }
  }
 
  CHKERR VecZeroEntries(epPtr->solTSStep);
 
  MoFEMFunctionReturn(0);
};
}