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 {
 
namespace {
 
MoFEMErrorCode resetIncrementalKinematics(EshelbianCore *ep_ptr,
                                          Vec solution) {
  MoFEMFunctionBegin;
 
  auto zero_field = [&](const std::string &field_name) {
    MoFEMFunctionBegin;
    auto is_mng = ep_ptr->mField.getInterface<ISManager>();
    SmartPetscObj<IS> is;
    CHKERR is_mng->isCreateProblemFieldAndRankLocal(
        "ELASTIC_PROBLEM", ROW, field_name, 0, MAX_DOFS_ON_ENTITY, is);
    const int *index_ptr;
    CHKERR ISGetIndices(is, &index_ptr);
    int size;
    CHKERR ISGetLocalSize(is, &size);
    double *array;
    CHKERR VecGetArray(solution, &array);
    for (int ii = 0; ii != size; ++ii) {
      array[index_ptr[ii]] = 0;
    }
    CHKERR VecRestoreArray(solution, &array);
    CHKERR ISRestoreIndices(is, &index_ptr);
    MoFEMFunctionReturn(0);
  };
 
  switch (EshelbianCore::gradApproximator) {
  case NO_H1_CONFIGURATION:
    break;
  case LARGE_ROT:
  case MODERATE_ROT:
  case SMALL_ROT:
    MOFEM_LOG("EP", Sev::verbose)
        << "Resetting incremental stretch and rotation after H1 transfer";
    CHKERR zero_field(ep_ptr->stretchTensor);
    CHKERR zero_field(ep_ptr->rotAxis);
    CHKERR VecGhostUpdateBegin(solution, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(solution, INSERT_VALUES, SCATTER_FORWARD);
    break;
  }
 
  MoFEMFunctionReturn(0);
}
 
} // namespace
 
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, 3>;
      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(nullptr, true));
      fe_rhs->getUserPolynomialBase() = boost::shared_ptr<BaseFunction>(
          new CGGUserPolynomialBase(nullptr, true));
      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) {
 
    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 append_meshsets_from_block = [](MoFEM::Interface &m_field,
                                             const std::string &block_name,
                                             const int dim,
                                             std::vector<EntityHandle> &r) {
          MoFEMFunctionBegin;
          auto mesh_mng = m_field.getInterface<MeshsetsManager>();
          auto bcs = mesh_mng->getCubitMeshsetPtr(
 
              std::regex((boost::format("%s(.*)") % block_name).str())
 
          );
 
          for (auto bc : bcs) {
            Range ents;
            CHKERR m_field.get_moab().get_entities_by_dimension(
                bc->getMeshset(), dim, ents, true);
            if (!ents.empty()) {
              r.push_back(bc->getMeshset());
            }
          }
 
          MoFEMFunctionReturn(0);
        };
 
        if (!epPtr->mField.get_comm_rank()) {
          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);
 
          CHKERR append_meshsets_from_block(epPtr->mField, "CRACK", 2,
                                            meshsets);
          CHKERR append_meshsets_from_block(epPtr->mField, "FRONT", 1,
                                            meshsets);
          CHKERR append_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(),
                                                     static_cast<int>(meshsets.size()));
        }
      }
 
      MoFEMFunctionReturn(0);
    };
 
      if (EshelbianCore::crackingStartTime < time) {
        CHKERR epPtr->setNewFrontCoordinates();
        CHKERR debug_crack();
        CHKERR epPtr->addCrackSurfaces();
      }
      CHKERR epPtr->projectGeometry(0, time);
      CHKERR epPtr->projectMaterialTags();
    }
 
  Vec T;
  CHKERR TSGetSolution(ts, &T);
 
  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);
 
  // The projected H1 field now carries the accepted deformation. Reset the
  // incremental stretch and rotation before crack calculations and TS
  // monitors so that the accepted increment is not applied a second time.
  Vec T;
  CHKERR TSGetSolution(ts, &T);
  CHKERR resetIncrementalKinematics(epPtr, T);
 
  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 (time - EshelbianCore::crackingStartTime >
        std::numeric_limits<double>::epsilon()) {
      MOFEM_LOG("EP", Sev::inform) << "Calculate faces orientation";
      CHKERR epPtr->calculateOrientation(1, true);
    }
    CHKERR epPtr->calculateCrackArea(epPtr->currentCrackAreaPtr);
  }
 
  CHKERR VecZeroEntries(epPtr->solTSStep);
 
  PetscBool local_stop = PETSC_FALSE;
  PetscBool global_stop = PETSC_FALSE;
  switch (EshelbianCore::solverType) {
  case EshelbianCore::LoadFactor:
    if (!epPtr->mField.get_comm_rank() && EshelbianCore::crackingOn) {
      const double area_increase =
          *epPtr->currentCrackAreaPtr - epPtr->oldCrackArea;
      if (area_increase > EshelbianCore::maxCrackExtension ||
          epPtr->avgGriffithsEnergy < epPtr->griffithEnergy) {
        MOFEM_LOG_C("EP", Sev::warning,
                    "Stopping crack growth. Area increase: %3.12e, "
                    "average Griffiths energy: %3.12e, critical Griffiths "
                    "energy: %3.12e",
                    area_increase, epPtr->avgGriffithsEnergy,
                    epPtr->griffithEnergy);
        local_stop = PETSC_TRUE;
      }
    }
 
    MPI_Bcast(&local_stop, 1, MPIU_BOOL, 0, PETSC_COMM_WORLD);
    global_stop = local_stop;
 
    if (global_stop) {
      TSSetConvergedReason(ts, TS_CONVERGED_USER);
    }
    break;
  default:
    break;
  }
 
  MoFEMFunctionReturn(0);
};
} // namespace EshelbianPlasticity