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";
 
    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;
    }
 
    MoFEMFunctionReturn(0);
  }
 
  static MoFEMErrorCode postStepFun(TS ts) {
    MoFEMFunctionBegin;
    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);
    // double time;
    // CHKERR TSGetTime(ts, &time);
    // preProcRhs->ts_t = time;
    // CHKERR DMoFEMPreProcessFiniteElements(prjDM, preProcRhs.get());
    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;