Collaboration diagram for ArcLengthInterfaceExample:

[legend]

Public Member Functions
	ArcLengthInterfaceExample (moab::Interface &moab)

MoFEMErrorCode	runProblem ()
	[Run problem]

Protected Member Functions
MoFEMErrorCode	readMesh ()
	[Run problem]

MoFEMErrorCode	setupProblem ()
	[Read mesh]

MoFEMErrorCode	boundaryCondition ()
	[Set up problem]

MoFEMErrorCode	assembleSystem ()
	[Boundary condition]

MoFEMErrorCode	solveSystem ()
	[Solve]

MoFEMErrorCode	outputResults ()
	[Solve]

MoFEMErrorCode	checkResults ()
	[Postprocess results]

MoFEMErrorCode	evaluateFieldAtPoint (const int load_step)
	[Assemble system]

MoFEMErrorCode	evaluateReactions (const int load_step)

Protected Attributes
moab::Interface &	mOab

boost::shared_ptr< MoFEM::Core >	corePtr

MoFEM::Interface *	mField = nullptr

double *	stepSizePtr = nullptr

int *	stepPtr = nullptr

PetscScalar	stepSizeReduction = 1.

PetscInt	maxSteps = 5

PetscInt	itsD = 6

PetscInt	order = 2

const std::string	naturalBcFeName = "NATURAL_BC"

SmartPetscObj< Vec >	f

SmartPetscObj< Vec >	d

Vec	dRaw = nullptr

SmartPetscObj< Mat >	aij

SmartPetscObj< Mat >	shellAij

SmartPetscObj< SNES >	snes

boost::ptr_vector< CohesiveInterfaceElement::PhysicalEquation >	interfaceMaterials

boost::shared_ptr< ArcLengthCtx >	arcCtx

boost::shared_ptr< ArcLengthIntElemFEMethod >	arcMethodPtr

boost::shared_ptr< ArcLengthSnesCtx >	snesCtxPtr

boost::shared_ptr< AssembleRhsVectors >	prePostProcFePtr

boost::shared_ptr< CohesiveInterfaceElement >	cohesiveElementsPtr

boost::shared_ptr< DomainEle >	elasticRhsFe

boost::shared_ptr< DomainEle >	elasticLhsFe

boost::shared_ptr< BoundaryEle >	naturalBcFe

boost::shared_ptr< PostProcEle >	postProcPtr

boost::scoped_ptr< ArcLengthMatShell >	matCtx

boost::scoped_ptr< PCArcLengthCtx >	pcCtx

std::array< double, SPACE_DIM >	fieldEvalCoords {0., 0., 0.}

PetscBool	doEvalField = PETSC_FALSE

boost::shared_ptr< FieldEvaluatorInterface::SetPtsData >	fieldEvalData

boost::shared_ptr< MatrixDouble >	fieldEvalDispPtr

boost::shared_ptr< MatrixDouble >	fieldEvalStressPtr

double	gamma = 0.5

double	reduction = 1.

Detailed Description

Examples: mofem/tutorials/cor-12_cohesive_interface/arc_length_interface.cpp.

Definition at line 104 of file arc_length_interface.cpp.

Constructor & Destructor Documentation

◆ ArcLengthInterfaceExample()

ArcLengthInterfaceExample::ArcLengthInterfaceExample ( moab::Interface & moab )

inline

Definition at line 106 of file arc_length_interface.cpp.

106: mOab(moab) {}

ArcLengthInterfaceExample::mOab

moab::Interface & mOab

Definition arc_length_interface.cpp:111

Member Function Documentation

◆ assembleSystem()

MoFEMErrorCode ArcLengthInterfaceExample::assembleSystem ( )

protected

[Boundary condition]

[Assemble system]

Examples: mofem/tutorials/cor-12_cohesive_interface/arc_length_interface.cpp.

Definition at line 583 of file arc_length_interface.cpp.

                                                         {
  MoFEMFunctionBegin;
  auto &m_field = *mField;
  auto &step = *stepPtr;
  auto &step_size = *stepSizePtr;
 
  // create matrices
  CHKERR m_field.getInterface<VecManager>()->vecCreateGhost("ELASTIC_MECHANICS",
                                                            COL, f);
  d = vectorDuplicate(f);
  dRaw = d;
  CHKERR m_field.getInterface<MatrixManager>()
      ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("ELASTIC_MECHANICS", aij);
 
  auto meshsets_manager = m_field.getInterface<MeshsetsManager>();
  const std::string interface_block_name = "MAT_INTERF";
  for (auto m : meshsets_manager->getCubitMeshsetPtr(std::regex(
           (boost::format("%s(.*)") % interface_block_name).str()))) {
    MOFEM_LOG("ARC_LENGTH", Sev::inform) << *m;
    Mat_Interf mydata;
    CHKERR m->getAttributeDataStructure(mydata);
    MOFEM_LOG("ARC_LENGTH", Sev::inform) << mydata;
 
    interfaceMaterials.push_back(
        new CohesiveInterfaceElement::PhysicalEquation(m_field));
    interfaceMaterials.back().h = 1;
    interfaceMaterials.back().youngModulus = mydata.data.alpha;
    interfaceMaterials.back().beta = mydata.data.beta;
    interfaceMaterials.back().ft = mydata.data.ft;
    interfaceMaterials.back().Gf = mydata.data.Gf;
 
    EntityHandle meshset = m->getMeshset();
    Range tris;
    CHKERR mOab.get_entities_by_type(meshset, MBTRI, tris, true);
    Range ents3d;
    CHKERR mOab.get_adjacencies(tris, 3, false, ents3d, moab::Interface::UNION);
    interfaceMaterials.back().pRisms = ents3d.subset_by_type(MBPRISM);
  }
 
  arcCtx = boost::make_shared<ArcLengthCtx>(m_field, "ELASTIC_MECHANICS");
  arcMethodPtr =
      boost::make_shared<ArcLengthIntElemFEMethod>(m_field.get_moab(), arcCtx);
  snesCtxPtr = boost::make_shared<ArcLengthSnesCtx>(
      m_field, "ELASTIC_MECHANICS", arcCtx);
  prePostProcFePtr = boost::make_shared<AssembleRhsVectors>(arcCtx);
  auto &snes_ctx = *snesCtxPtr;
  auto &pre_post_proc_fe = *prePostProcFePtr;
 
  auto essential_pre_proc_zero_rhs = boost::make_shared<FEMethod>();
  auto essential_pre_proc_lhs = boost::make_shared<FEMethod>();
  auto essential_post_proc_rhs = boost::make_shared<FEMethod>();
  auto essential_post_proc_lhs = boost::make_shared<FEMethod>();
 
  struct ArcScale : public ScalingMethod {
    boost::shared_ptr<ArcLengthCtx> arcCtx;
 
    ArcScale(boost::shared_ptr<ArcLengthCtx> arc_ctx) : arcCtx(arc_ctx) {}
 
    double getScale(const double time) override {
      return arcCtx->getFieldData();
    }
  };
  auto arc_scale = boost::make_shared<ArcScale>(arcCtx);
 
  struct ArcZero : public ScalingMethod {
    double getScale(const double time) override {
      return 0.0;
    }
  };
  auto arc_zero = boost::make_shared<ArcZero>();
 
  auto get_essential_pre_proc_zero =
      [&m_field, essential_pre_proc_zero_rhs, arc_zero]() {
    return EssentialPreProc<DisplacementCubitBcData>(
        m_field, essential_pre_proc_zero_rhs, {arc_zero});
  };
  essential_pre_proc_zero_rhs->preProcessHook =
      get_essential_pre_proc_zero();
  auto get_essential_pre_proc_lhs =
      [&m_field, essential_pre_proc_lhs, arc_scale]() {
    return EssentialPreProc<DisplacementCubitBcData>(
        m_field, essential_pre_proc_lhs, {arc_scale});
  };
  essential_pre_proc_lhs->preProcessHook = get_essential_pre_proc_lhs();
 
  essential_post_proc_rhs->postProcessHook =
      [&m_field, essential_post_proc_rhs, arc_scale]() {
    MoFEMFunctionBegin;
    CHKERR EssentialPreProc<DisplacementCubitBcData>(
        m_field, essential_post_proc_rhs, {arc_scale})();
    CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(
        m_field, essential_post_proc_rhs, 1.)();
    MoFEMFunctionReturn(0);
  };
 
  // ArcScale sets the reference prescribed field to lambda * F_lambda, so
  // the constrained residual is u_D - lambda * F_lambda,D.
  essential_post_proc_lhs->postProcessHook =
      EssentialPostProcLhs<DisplacementCubitBcData>(
          m_field, essential_post_proc_lhs, 1.);
 
  elasticRhsFe = boost::make_shared<DomainEle>(m_field);
  elasticLhsFe = boost::make_shared<DomainEle>(m_field);
  auto elastic_integration_rule = [](int, int, int approx_order) {
    return 2 * approx_order + 1;
  };
  elasticRhsFe->getRuleHook = elastic_integration_rule;
  elasticLhsFe->getRuleHook = elastic_integration_rule;
 
  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      elasticRhsFe->getOpPtrVector(), {H1}, "GEOMETRY");
  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      elasticLhsFe->getOpPtrVector(), {H1}, "GEOMETRY");
  CHKERR HookeOps::opFactoryDomainRhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
      m_field, elasticRhsFe->getOpPtrVector(), "DISPLACEMENT", "MAT_ELASTIC",
      Sev::verbose, true);
  CHKERR HookeOps::opFactoryDomainLhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
      m_field, elasticLhsFe->getOpPtrVector(), "DISPLACEMENT", "MAT_ELASTIC",
      Sev::verbose);
 
  cohesiveElementsPtr = boost::make_shared<CohesiveInterfaceElement>(m_field);
  auto &cohesive_elements = *cohesiveElementsPtr;
  CHKERR cohesive_elements.addOps("DISPLACEMENT", interfaceMaterials);
 
  PetscInt M, N;
  CHKERR MatGetSize(aij, &M, &N);
  PetscInt m, n;
  CHKERR MatGetLocalSize(aij, &m, &n);
  matCtx.reset(new ArcLengthMatShell(aij, arcCtx, "ELASTIC_MECHANICS"));
  Mat shell_aij;
  CHKERR MatCreateShell(PETSC_COMM_WORLD, m, n, M, N, (void *)matCtx.get(),
                        &shell_aij);
  shellAij = SmartPetscObj<Mat>(shell_aij);
  CHKERR MatShellSetOperation(shellAij, MATOP_MULT,
                              (void (*)(void))ArcLengthMatMultShellOp);
 
  using BoundaryRhsBCs = NaturalBC<BoundaryEle::UserDataOperator>::Assembly<
      PETSC>::LinearForm<GAUSS>;
  using OpBoundaryRhsBCs = BoundaryRhsBCs::OpFlux<BoundaryBCs, 1, SPACE_DIM>;
 
  naturalBcFe = boost::make_shared<BoundaryEle>(m_field);
  naturalBcFe->getRuleHook = [](int, int, int approx_order) {
    return 2 * approx_order + 1;
  };
  CHKERR AddHOOps<2, 3, 3>::add(naturalBcFe->getOpPtrVector(), {NOSPACE},
                                "GEOMETRY");
  CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(
      naturalBcFe->getOpPtrVector(), m_field, "DISPLACEMENT", Sev::inform);
 
  snes = createSNES(PETSC_COMM_WORLD);
  CHKERR SNESSetApplicationContext(snes, &snes_ctx);
  CHKERR SNESSetFunction(snes, f, SnesRhs, &snes_ctx);
  CHKERR SNESSetJacobian(snes, shellAij, aij, SnesMat, &snes_ctx);
  CHKERR SNESSetFromOptions(snes);
 
  KSP ksp;
  CHKERR SNESGetKSP(snes, &ksp);
  PC pc;
  CHKERR KSPGetPC(ksp, &pc);
  pcCtx.reset(new PCArcLengthCtx(shellAij, aij, arcCtx));
  CHKERR PCSetType(pc, PCSHELL);
  CHKERR PCShellSetContext(pc, pcCtx.get());
  CHKERR PCShellSetApply(pc, PCApplyArcLength);
  CHKERR PCShellSetSetUp(pc, PCSetupArcLength);
 
  // Rhs
  SnesCtx::FEMethodsSequence &loops_to_do_Rhs = snes_ctx.getComputeRhs();
  snes_ctx.getPreProcComputeRhs().push_back(essential_pre_proc_zero_rhs);
  snes_ctx.getPreProcComputeRhs().push_back(&pre_post_proc_fe);
  loops_to_do_Rhs.push_back(
      SnesCtx::PairNameFEMethodPtr("INTERFACE", &cohesive_elements.getFeRhs()));
  loops_to_do_Rhs.push_back(
      SnesCtx::PairNameFEMethodPtr("ELASTIC", elasticRhsFe.get()));
  loops_to_do_Rhs.push_back(
      SnesCtx::PairNameFEMethodPtr("ARC_LENGTH", arcMethodPtr.get()));
  // Add F_lambda before essential processing replaces the constrained rows
  // with u_D - lambda * F_lambda,D.
  snes_ctx.getPostProcComputeRhs().push_back(&pre_post_proc_fe);
  snes_ctx.getPostProcComputeRhs().push_back(essential_post_proc_rhs);
 
  // Mat
  SnesCtx::FEMethodsSequence &loops_to_do_Mat = snes_ctx.getSetOperators();
  snes_ctx.getPreProcSetOperators().push_back(essential_pre_proc_lhs);
  loops_to_do_Mat.push_back(
      SnesCtx::PairNameFEMethodPtr("INTERFACE", &cohesive_elements.getFeLhs()));
  loops_to_do_Mat.push_back(
      SnesCtx::PairNameFEMethodPtr("ELASTIC", elasticLhsFe.get()));
  loops_to_do_Mat.push_back(
      SnesCtx::PairNameFEMethodPtr("ARC_LENGTH", arcMethodPtr.get()));
  snes_ctx.getPostProcSetOperators().push_back(essential_post_proc_lhs);
 
  // step = 1;
  if (step == 1) {
    step_size = stepSizeReduction;
  } else {
    reduction = stepSizeReduction;
    step++;
  }
 
  auto f_lambda_on_essential_bc = boost::make_shared<FEMethod>();
  f_lambda_on_essential_bc->preProcessHook =
      EssentialPreProc<DisplacementCubitBcData>(
          m_field, f_lambda_on_essential_bc, {}, false);
  // EssentialPreProc stores the unit prescribed field. Its derivative in
  // u_D - lambda * F_lambda,D has the opposite sign, hence vDiag = -1.
  f_lambda_on_essential_bc->postProcessHook =
      EssentialPostProcRhs<DisplacementCubitBcData>(
          m_field, f_lambda_on_essential_bc, -1., arcCtx->F_lambda);
 
  CHKERR VecZeroEntries(arcCtx->F_lambda);
  CHKERR VecGhostUpdateBegin(arcCtx->F_lambda, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(arcCtx->F_lambda, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR m_field.problem_basic_method_preProcess("ELASTIC_MECHANICS",
                                                 *f_lambda_on_essential_bc);
  elasticRhsFe->ksp_f = arcCtx->F_lambda;
  CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", "ELASTIC",
                                      *elasticRhsFe);
  naturalBcFe->ksp_f = arcCtx->F_lambda;
  CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", naturalBcFeName,
                                      *naturalBcFe);
  CHKERR m_field.problem_basic_method_postProcess("ELASTIC_MECHANICS",
                                                  *f_lambda_on_essential_bc);
  CHKERR VecGhostUpdateBegin(arcCtx->F_lambda, ADD_VALUES, SCATTER_REVERSE);
  CHKERR VecGhostUpdateEnd(arcCtx->F_lambda, ADD_VALUES, SCATTER_REVERSE);
  CHKERR VecAssemblyBegin(arcCtx->F_lambda);
  CHKERR VecAssemblyEnd(arcCtx->F_lambda);
 
  // F_lambda2
  CHKERR VecDot(arcCtx->F_lambda, arcCtx->F_lambda, &arcCtx->F_lambda2);
  MOFEM_LOG_C("ARC_LENGTH", Sev::inform, "\tFlambda2 = %6.4e",
              arcCtx->F_lambda2);
 
  if (step > 1) {
    CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
        "ELASTIC_MECHANICS", COL, d, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR m_field.getInterface<VecManager>()->setOtherGlobalGhostVector(
        "ELASTIC_MECHANICS", "DISPLACEMENT", "X0_DISPLACEMENT", COL, arcCtx->x0,
        INSERT_VALUES, SCATTER_FORWARD);
    double x0_nrm;
    CHKERR VecNorm(arcCtx->x0, NORM_2, &x0_nrm);
    MOFEM_LOG_C("ARC_LENGTH", Sev::inform,
                "\tRead x0_nrm = %6.4e dlambda = %6.4e", x0_nrm,
                arcCtx->dLambda);
    CHKERR arcCtx->setAlphaBeta(1, 0);
  } else {
    CHKERR arcCtx->setS(0);
    CHKERR arcCtx->setAlphaBeta(0, 1);
  }
  CHKERR SnesRhs(snes, d, f, &snes_ctx);
 
  postProcPtr = boost::make_shared<PostProcEle>(m_field);
  auto &post_proc = *postProcPtr;
  auto disp_ptr = boost::make_shared<MatrixDouble>();
  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      post_proc.getOpPtrVector(), {H1}, "GEOMETRY");
  auto hooke_common_ptr =
      HookeOps::commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
          m_field, post_proc.getOpPtrVector(), "DISPLACEMENT", "MAT_ELASTIC",
          Sev::verbose);
  post_proc.getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<SPACE_DIM>("DISPLACEMENT", disp_ptr));
  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
  post_proc.getOpPtrVector().push_back(
      new OpPPMap(post_proc.getPostProcMesh(), post_proc.getMapGaussPts(), {},
                  {{"DISPLACEMENT", disp_ptr}},
                  {{"DISPLACEMENT_GRAD", hooke_common_ptr->matGradPtr}},
                  {{"STRAIN", hooke_common_ptr->getMatStrain()},
                   {"STRESS", hooke_common_ptr->getMatCauchyStress()}}));
 
  int coords_dim = SPACE_DIM;
  CHKERR PetscOptionsGetRealArray(PETSC_NULLPTR, PETSC_NULLPTR,
                                  "-field_eval_coords", fieldEvalCoords.data(),
                                  &coords_dim, &doEvalField);
  if (doEvalField && coords_dim != SPACE_DIM)
    SETERRQ(PETSC_COMM_WORLD, MOFEM_INVALID_DATA,
            "-field_eval_coords requires exactly three coordinates");
 
  if (doEvalField) {
    fieldEvalData =
        m_field.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
    CHKERR m_field.getInterface<FieldEvaluatorInterface>()
        ->buildTree<SPACE_DIM>(fieldEvalData, "ELASTIC");
    fieldEvalData->setEvalPoints(fieldEvalCoords.data(), 1);
 
    auto field_eval_fe = fieldEvalData->feMethodPtr;
    field_eval_fe->getRuleHook = [](int, int, int) { return -1; };
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
        field_eval_fe->getOpPtrVector(), {H1}, "GEOMETRY");
 
    auto field_eval_hooke_common_ptr =
        HookeOps::commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
            m_field, field_eval_fe->getOpPtrVector(), "DISPLACEMENT",
            "MAT_ELASTIC", Sev::verbose);
    fieldEvalStressPtr = field_eval_hooke_common_ptr->getMatCauchyStress();
    fieldEvalDispPtr = boost::make_shared<MatrixDouble>();
    field_eval_fe->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<SPACE_DIM>("DISPLACEMENT",
                                                    fieldEvalDispPtr));
  }
 
  MoFEMFunctionReturn(0);
}

◆ boundaryCondition()

MoFEMErrorCode ArcLengthInterfaceExample::boundaryCondition ( )

protected

[Set up problem]

[Boundary condition]

Examples: mofem/tutorials/cor-12_cohesive_interface/arc_length_interface.cpp.

Definition at line 562 of file arc_length_interface.cpp.

                                                            {
  MoFEMFunctionBegin;
  auto &m_field = *mField;
 
  auto bc_mng = m_field.getInterface<BcManager>();
  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
      "ELASTIC_MECHANICS", "DISPLACEMENT");
 
  // print bcs
  MeshsetsManager *mmanager_ptr;
  CHKERR m_field.getInterface(mmanager_ptr);
  CHKERR mmanager_ptr->printDisplacementSet();
  CHKERR mmanager_ptr->printForceSet();
  // print block sets with materials
  CHKERR mmanager_ptr->printMaterialsSet();
 
  MoFEMFunctionReturn(0);
}

◆ checkResults()

MoFEMErrorCode ArcLengthInterfaceExample::checkResults ( )

protected

[Postprocess results]

[Check]

Examples: mofem/tutorials/cor-12_cohesive_interface/arc_length_interface.cpp.

Definition at line 1113 of file arc_length_interface.cpp.

                                                       {
  MoFEMFunctionBegin;
  MoFEMFunctionReturn(0);
}

◆ evaluateFieldAtPoint()

MoFEMErrorCode ArcLengthInterfaceExample::evaluateFieldAtPoint ( const int load_step )

protected

[Assemble system]

Examples: mofem/tutorials/cor-12_cohesive_interface/arc_length_interface.cpp.

Definition at line 888 of file arc_length_interface.cpp.

                                                                   {
  MoFEMFunctionBegin;
  auto &m_field = *mField;
 
  if (!doEvalField)
    MoFEMFunctionReturnHot(0);
 
  CHKERR m_field.getInterface<FieldEvaluatorInterface>()
      ->evalFEAtThePoint<SPACE_DIM>(
          fieldEvalCoords.data(), 1e-12, "ELASTIC_MECHANICS", "ELASTIC",
          fieldEvalData, m_field.get_comm_rank(), m_field.get_comm_rank(),
          nullptr, MF_EXIST, QUIET);
 
  int point_found = fieldEvalDispPtr->size1() && fieldEvalStressPtr->size1();
  int global_point_found = 0;
  MPI_Allreduce(&point_found, &global_point_found, 1, MPI_INT, MPI_SUM,
                m_field.get_comm());
 
  if (point_found) {
    auto t_disp = getFTensor1FromMat<SPACE_DIM>(*fieldEvalDispPtr);
    auto t_stress = getFTensor2SymmetricFromMat<SPACE_DIM>(*fieldEvalStressPtr);
    MOFEM_LOG_C("SYNC", Sev::inform,
                "FieldEvaluator step %d lambda %6.4e point [%6.4e, %6.4e, "
                "%6.4e]",
                load_step, arcCtx->getFieldData(), fieldEvalCoords[0],
                fieldEvalCoords[1], fieldEvalCoords[2]);
    MOFEM_LOG_C("SYNC", Sev::inform,
                "FieldEvaluator displacement [Ux, Uy, Uz] = [%6.4e, "
                "%6.4e, %6.4e]",
                t_disp(0), t_disp(1), t_disp(2));
    MOFEM_LOG_C("SYNC", Sev::inform,
                "FieldEvaluator Cauchy stress [Sxx, Syy, Szz, Sxy, Syz, Sxz] "
                "= [%6.4e, %6.4e, %6.4e, %6.4e, %6.4e, %6.4e]",
                t_stress(0, 0), t_stress(1, 1), t_stress(2, 2), t_stress(0, 1),
                t_stress(1, 2), t_stress(0, 2));
  } else if (!global_point_found && !m_field.get_comm_rank()) {
    MOFEM_LOG_C("ARC_LENGTH", Sev::warning,
                "FieldEvaluator did not find point [%6.4e, %6.4e, %6.4e] in "
                "the ELASTIC domain",
                fieldEvalCoords[0], fieldEvalCoords[1], fieldEvalCoords[2]);
  }
  MOFEM_LOG_SYNCHRONISE(m_field.get_comm());
 
  MoFEMFunctionReturn(0);
}

◆ evaluateReactions()

MoFEMErrorCode ArcLengthInterfaceExample::evaluateReactions ( const int load_step )

protected

Examples: mofem/tutorials/cor-12_cohesive_interface/arc_length_interface.cpp.

Definition at line 935 of file arc_length_interface.cpp.

                                                                {
  MoFEMFunctionBegin;
  auto &m_field = *mField;
 
  auto domain_reaction_fe = boost::make_shared<DomainEle>(m_field);
  domain_reaction_fe->getRuleHook = [](int, int, int approx_order) {
    return 2 * approx_order + 1;
  };
 
  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      domain_reaction_fe->getOpPtrVector(), {H1}, "GEOMETRY");
  CHKERR HookeOps::opFactoryDomainRhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
      m_field, domain_reaction_fe->getOpPtrVector(), "DISPLACEMENT",
      "MAT_ELASTIC", Sev::verbose, true);
 
  auto reactions = vectorDuplicate(f);
  CHKERR VecZeroEntries(reactions);
  domain_reaction_fe->f = reactions;
  domain_reaction_fe->vecAssembleSwitch =
      boost::movelib::make_unique<bool>(true);
  domain_reaction_fe->postProcessHook =
      EssentialPreProcReaction<DisplacementCubitBcData>(m_field,
                                                        domain_reaction_fe);
 
  MOFEM_LOG_C("ARC_LENGTH", Sev::inform,
              "Reactions step %d lambda %6.4e", load_step,
              arcCtx->getFieldData());
  CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", "ELASTIC",
                                      *domain_reaction_fe);
 
  MoFEMFunctionReturn(0);
}

◆ outputResults()

MoFEMErrorCode ArcLengthInterfaceExample::outputResults ( )

protected

[Solve]

[Postprocess results]

Examples: mofem/tutorials/cor-12_cohesive_interface/arc_length_interface.cpp.

Definition at line 1106 of file arc_length_interface.cpp.

                                                        {
  MoFEMFunctionBegin;
  MoFEMFunctionReturn(0);
}

◆ readMesh()

MoFEMErrorCode ArcLengthInterfaceExample::readMesh ( )

protected

[Run problem]

[Read mesh]

Examples: mofem/tutorials/cor-12_cohesive_interface/arc_length_interface.cpp.

Definition at line 179 of file arc_length_interface.cpp.

                                                   {
  MoFEMFunctionBegin;
 
  CHKERR JsonConfigManager::setMeshFileFromJson();
 
  // Read the mesh file before constructing MoFEM::Core so that the MOAB
  // database can be populated first. JsonConfigManager::setMeshFileFromJson()
  // maps the primary JSON mesh to -file_name.
  PetscBool flg = PETSC_TRUE;
  char mesh_file_name[PETSC_MAX_PATH_LEN] = "";
  CHKERR PetscOptionsGetString(PETSC_NULLPTR, PETSC_NULLPTR, "-my_file",
                               mesh_file_name, sizeof(mesh_file_name), &flg);
  if (flg != PETSC_TRUE) {
    CHKERR PetscOptionsGetString(PETSC_NULLPTR, PETSC_NULLPTR, "-file_name",
                                 mesh_file_name, sizeof(mesh_file_name), &flg);
    if (flg != PETSC_TRUE)
      SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
              "*** ERROR -my_file or -file_name (MESH FILE NEEDED)");
  }
 
  // Read mesh to MOAB
  const char *option;
  option = "";
  CHKERR mOab.load_file(mesh_file_name, 0, option);
 
  // Data stored on mesh for restart
  Tag th_step_size, th_step;
  double def_step_size = 1;
  rval = mOab.tag_get_handle("_STEPSIZE", 1, MB_TYPE_DOUBLE, th_step_size,
                             MB_TAG_CREAT | MB_TAG_MESH, &def_step_size);
  if (rval == MB_ALREADY_ALLOCATED)
    rval = MB_SUCCESS;
  CHKERR rval;
  int def_step = 1;
  rval = mOab.tag_get_handle("_STEP", 1, MB_TYPE_INTEGER, th_step,
                             MB_TAG_CREAT | MB_TAG_MESH, &def_step);
  if (rval == MB_ALREADY_ALLOCATED)
    rval = MB_SUCCESS;
  CHKERR rval;
  const void *tag_data_step_size[1];
  EntityHandle root = 0;
  CHKERR mOab.tag_get_by_ptr(th_step_size, &root, 1, tag_data_step_size);
  double &step_size = *(double *)tag_data_step_size[0];
  const void *tag_data_step[1];
  CHKERR mOab.tag_get_by_ptr(th_step, &root, 1, tag_data_step);
  int &step = *(int *)tag_data_step[0];
  // end of data stored for restart
  MOFEM_LOG_C("ARC_LENGTH", Sev::inform, "Start step %d and step_size = %6.4e",
              step, step_size);
 
  stepSizePtr = (double *)tag_data_step_size[0];
  stepPtr = (int *)tag_data_step[0];
 
  // Create MoFEM database after loading the mesh. This preserves loading of
  // JSON options by the Core constructor.
  corePtr = boost::make_shared<MoFEM::Core>(mOab);
  mField = corePtr.get();
 
  MoFEMFunctionReturn(0);
}

◆ runProblem()

MoFEMErrorCode ArcLengthInterfaceExample::runProblem ( )

[Run problem]

Examples: mofem/tutorials/cor-12_cohesive_interface/arc_length_interface.cpp.

Definition at line 165 of file arc_length_interface.cpp.

                                                     {
  MoFEMFunctionBegin;
  CHKERR readMesh();
  CHKERR setupProblem();
  CHKERR boundaryCondition();
  CHKERR assembleSystem();
  CHKERR solveSystem();
  CHKERR outputResults();
  CHKERR checkResults();
  MoFEMFunctionReturn(0);
}

◆ setupProblem()

MoFEMErrorCode ArcLengthInterfaceExample::setupProblem ( )

protected

[Read mesh]

[Set up problem]

Examples: mofem/tutorials/cor-12_cohesive_interface/arc_length_interface.cpp.

Definition at line 242 of file arc_length_interface.cpp.

                                                       {
  MoFEMFunctionBegin;
  auto &m_field = *mField;
  auto &step = *stepPtr;
  PetscBool flg = PETSC_FALSE;
 
  // Constructing MoFEM::Core loads -json_config and inserts its "petsc" and
  // "mofem" sections into the PETSc options database. Read application
  // controls only after that point so they can be supplied by JSON.
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR, "-my_sr",
                             &stepSizeReduction, &flg);
  if (flg != PETSC_TRUE) {
    stepSizeReduction = 1.;
  }
 
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, PETSC_NULLPTR, "-my_ms", &maxSteps,
                            &flg);
  if (flg != PETSC_TRUE) {
    maxSteps = 5;
  }
 
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, PETSC_NULLPTR, "-my_its_d", &itsD,
                            &flg);
  if (flg != PETSC_TRUE) {
    itsD = 6;
  }
 
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, PETSC_NULLPTR, "-my_order", &order,
                            &flg);
  if (flg != PETSC_TRUE) {
    order = 2;
  }
 
  MeshsetsManager *meshsets_manager_ptr;
  CHKERR m_field.getInterface(meshsets_manager_ptr);
  CHKERR meshsets_manager_ptr->setMeshsetFromFile();
 
  PrismInterface *interface_ptr;
  CHKERR m_field.getInterface(interface_ptr);
 
  Tag th_my_ref_level;
  BitRefLevel def_bit_level = 0;
  CHKERR m_field.get_moab().tag_get_handle(
      "_MY_REFINEMENT_LEVEL", sizeof(BitRefLevel), MB_TYPE_OPAQUE,
      th_my_ref_level, MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_BYTES,
      &def_bit_level);
  const EntityHandle root_meshset = m_field.get_moab().get_root_set();
  BitRefLevel *ptr_bit_level0;
  CHKERR m_field.get_moab().tag_get_by_ptr(th_my_ref_level, &root_meshset, 1,
                                           (const void **)&ptr_bit_level0);
  BitRefLevel &bit_level0 = *ptr_bit_level0;
  BitRefLevel problem_bit_level = bit_level0;
 
  if (step == 1) {
 
    // ref meshset ref level 0
    CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
        0, 3, BitRefLevel().set(0));
 
    std::vector<BitRefLevel> bit_levels;
    bit_levels.push_back(BitRefLevel().set(0));
 
    int ll = 1;
 
    auto split_interface = [&](const EntityHandle cubit_meshset,
                               const int meshset_id) -> MoFEMErrorCode {
      MoFEMFunctionBegin;
      MOFEM_LOG_C("ARC_LENGTH", Sev::inform, "Insert Interface %d", meshset_id);
      {
        // get tet entities form back bit_level
        EntityHandle ref_level_meshset = 0;
        CHKERR mOab.create_meshset(MESHSET_SET, ref_level_meshset);
        CHKERR m_field.getInterface<BitRefManager>()
            ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
                                           BitRefLevel().set(), MBTET,
                                           ref_level_meshset);
        CHKERR m_field.getInterface<BitRefManager>()
            ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
                                           BitRefLevel().set(), MBPRISM,
                                           ref_level_meshset);
        Range ref_level_tets;
        CHKERR mOab.get_entities_by_handle(ref_level_meshset, ref_level_tets,
                                           true);
        // get faces and test to split
        CHKERR interface_ptr->getSides(cubit_meshset, bit_levels.back(), true,
                                       0);
        // set new bit level
        bit_levels.push_back(BitRefLevel().set(ll++));
        // split faces and
        CHKERR interface_ptr->splitSides(ref_level_meshset, bit_levels.back(),
                                         cubit_meshset, true, true, 0);
        // clean meshsets
        CHKERR mOab.delete_entities(&ref_level_meshset, 1);
      }
      // Update cubit meshsets
      for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, ciit)) {
        EntityHandle cubit_meshset = ciit->meshset;
        CHKERR m_field.getInterface<BitRefManager>()
            ->updateMeshsetByEntitiesChildren(cubit_meshset, bit_levels.back(),
                                              cubit_meshset, MBMAXTYPE, true);
      }
      MoFEMFunctionReturn(0);
    };
 
    // Prefer Cubit INTERFACESET side sets. Use MAT_INTERF blocksets only as
    // a fallback for meshes configured through JsonConfigManager, since a
    // legacy mesh can store the same interface under both set types.
    std::set<int> split_interface_ids;
    bool split_from_sideset = false;
    for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
             m_field, SIDESET | INTERFACESET, cit)) {
      if (split_interface_ids.insert(cit->getMeshsetId()).second) {
        CHKERR split_interface(cit->getMeshset(), cit->getMeshsetId());
        split_from_sideset = true;
      }
    }
    if (!split_from_sideset) {
      const std::string block_name = "MAT_INTERF";
      for (auto m : m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
               std::regex((boost::format("%s(.*)") % block_name).str()))) {
        if (split_interface_ids.insert(m->getMeshsetId()).second)
          CHKERR split_interface(m->getMeshset(), m->getMeshsetId());
      }
    }
 
    bit_level0 = bit_levels.back();
    problem_bit_level = bit_level0;
 
    /***/
    // Define problem
 
    // Fields
    CHKERR m_field.add_field("DISPLACEMENT", H1, AINSWORTH_LEGENDRE_BASE, 3);
    CHKERR
    m_field.add_field("GEOMETRY", H1, AINSWORTH_LEGENDRE_BASE, 3);
 
    CHKERR m_field.add_field("LAMBDA", NOFIELD, NOBASE, 1);
 
    // Field for ArcLength
    CHKERR
    m_field.add_field("X0_DISPLACEMENT", H1, AINSWORTH_LEGENDRE_BASE, 3);
 
    // FE
    CHKERR m_field.add_finite_element("ELASTIC");
 
    // Define rows/cols and element data
    CHKERR m_field.modify_finite_element_add_field_row("ELASTIC",
                                                       "DISPLACEMENT");
    CHKERR m_field.modify_finite_element_add_field_col("ELASTIC",
                                                       "DISPLACEMENT");
    CHKERR m_field.modify_finite_element_add_field_data("ELASTIC",
                                                        "DISPLACEMENT");
    CHKERR
    m_field.modify_finite_element_add_field_data("ELASTIC", "GEOMETRY");
    CHKERR m_field.modify_finite_element_add_field_row("ELASTIC", "LAMBDA");
    CHKERR m_field.modify_finite_element_add_field_col("ELASTIC", "LAMBDA");
    // this is for paremtis
    CHKERR m_field.modify_finite_element_add_field_data("ELASTIC", "LAMBDA");
 
    // FE Interface
    CHKERR m_field.add_finite_element("INTERFACE");
    CHKERR m_field.modify_finite_element_add_field_row("INTERFACE",
                                                       "DISPLACEMENT");
    CHKERR m_field.modify_finite_element_add_field_col("INTERFACE",
                                                       "DISPLACEMENT");
    CHKERR m_field.modify_finite_element_add_field_data("INTERFACE",
                                                        "DISPLACEMENT");
    CHKERR
    m_field.modify_finite_element_add_field_data("INTERFACE", "GEOMETRY");
 
    // FE ArcLength
    CHKERR m_field.add_finite_element("ARC_LENGTH");
 
    // Define rows/cols and element data
    CHKERR m_field.modify_finite_element_add_field_row("ARC_LENGTH", "LAMBDA");
    CHKERR m_field.modify_finite_element_add_field_col("ARC_LENGTH", "LAMBDA");
 
    // elem data
    CHKERR
    m_field.modify_finite_element_add_field_data("ARC_LENGTH", "LAMBDA");
 
    // define problems
    CHKERR m_field.add_problem("ELASTIC_MECHANICS");
 
    // set finite elements for problem
    CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
                                                     "ELASTIC");
    CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
                                                     "INTERFACE");
    CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
                                                     "ARC_LENGTH");
    // set refinement level for problem
    CHKERR m_field.modify_problem_ref_level_add_bit("ELASTIC_MECHANICS",
                                                    problem_bit_level);
 
    /***/
    // Declare problem
 
    // add entities (by tets) to the field
    CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "DISPLACEMENT");
    CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "X0_DISPLACEMENT");
    CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "GEOMETRY");
 
    // add finite elements entities
    CHKERR m_field.add_ents_to_finite_element_by_bit_ref(
        problem_bit_level, BitRefLevel().set(), "ELASTIC", MBTET);
    CHKERR m_field.add_ents_to_finite_element_by_bit_ref(
        problem_bit_level, BitRefLevel().set(), "INTERFACE", MBPRISM);
 
    // Setting up LAMBDA field and ARC_LENGTH interface
    {
      // Add dummy no-field vertex
      EntityHandle no_field_vertex;
      {
        const double coords[] = {0, 0, 0};
        CHKERR m_field.get_moab().create_vertex(coords, no_field_vertex);
        Range range_no_field_vertex;
        range_no_field_vertex.insert(no_field_vertex);
        CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(
            range_no_field_vertex, BitRefLevel().set());
 
        EntityHandle lambda_meshset = m_field.get_field_meshset("LAMBDA");
        CHKERR m_field.get_moab().add_entities(lambda_meshset,
                                               range_no_field_vertex);
      }
      // this entity will carray data for this finite element
      EntityHandle meshset_fe_arc_length;
      {
        CHKERR mOab.create_meshset(MESHSET_SET, meshset_fe_arc_length);
        CHKERR mOab.add_entities(meshset_fe_arc_length, &no_field_vertex, 1);
        CHKERR m_field.getInterface<BitRefManager>()->setBitLevelToMeshset(
            meshset_fe_arc_length, BitRefLevel().set());
      }
      // finally add created meshset to the ARC_LENGTH finite element
      CHKERR m_field.add_ents_to_finite_element_by_MESHSET(
          meshset_fe_arc_length, "ARC_LENGTH", false);
    }
 
    // set app. order
    // see Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes
    // (Mark Ainsworth & Joe Coyle)
    CHKERR m_field.set_field_order(0, MBTET, "DISPLACEMENT", order);
    CHKERR m_field.set_field_order(0, MBTRI, "DISPLACEMENT", order);
    CHKERR m_field.set_field_order(0, MBEDGE, "DISPLACEMENT", order);
    CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", 1);
 
    CHKERR m_field.set_field_order(0, MBTET, "X0_DISPLACEMENT", order);
    CHKERR m_field.set_field_order(0, MBTRI, "X0_DISPLACEMENT", order);
    CHKERR m_field.set_field_order(0, MBEDGE, "X0_DISPLACEMENT", order);
    CHKERR m_field.set_field_order(0, MBVERTEX, "X0_DISPLACEMENT", 1);
 
    CHKERR m_field.set_field_order(0, MBTET, "GEOMETRY", 2);
    CHKERR m_field.set_field_order(0, MBTRI, "GEOMETRY", 2);
    CHKERR m_field.set_field_order(0, MBEDGE, "GEOMETRY", 2);
    CHKERR m_field.set_field_order(0, MBVERTEX, "GEOMETRY", 1);
  }
 
  if (!m_field.check_finite_element(naturalBcFeName)) {
    CHKERR m_field.add_finite_element(naturalBcFeName);
    CHKERR m_field.modify_finite_element_add_field_row(naturalBcFeName,
                                                       "DISPLACEMENT");
    CHKERR m_field.modify_finite_element_add_field_col(naturalBcFeName,
                                                       "DISPLACEMENT");
    CHKERR m_field.modify_finite_element_add_field_data(naturalBcFeName,
                                                        "DISPLACEMENT");
    CHKERR m_field.modify_finite_element_add_field_data(naturalBcFeName,
                                                        "GEOMETRY");
    CHKERR m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS",
                                                     naturalBcFeName);
 
    Range natural_bc_faces;
    auto meshsets_manager = m_field.getInterface<MeshsetsManager>();
    for (auto m : meshsets_manager->getCubitMeshsetPtr(NODESET | FORCESET)) {
      CHKERR mOab.get_entities_by_type(m->getMeshset(), MBTRI, natural_bc_faces,
                                       true);
    }
    for (auto m : meshsets_manager->getCubitMeshsetPtr(SIDESET | PRESSURESET)) {
      CHKERR mOab.get_entities_by_type(m->getMeshset(), MBTRI, natural_bc_faces,
                                       true);
    }
    for (auto m : meshsets_manager->getCubitMeshsetPtr(
             std::regex("(FORCE|PRESSURE)(.*)"))) {
      CHKERR mOab.get_entities_by_type(m->getMeshset(), MBTRI, natural_bc_faces,
                                       true);
    }
    CHKERR m_field.add_ents_to_finite_element_by_type(natural_bc_faces, MBTRI,
                                                      naturalBcFeName);
  }
 
  /****/
  // build database
 
  // build field
  CHKERR m_field.build_fields();
  Projection10NodeCoordsOnField ent_method_material(m_field, "GEOMETRY");
  CHKERR m_field.loop_dofs("GEOMETRY", ent_method_material);
 
  // build finite elemnts
  CHKERR m_field.build_finite_elements();
 
  // build adjacencies
  CHKERR m_field.build_adjacencies(problem_bit_level);
 
  /****/
  ProblemsManager *prb_mng_ptr;
  CHKERR m_field.getInterface(prb_mng_ptr);
  // build problem
  CHKERR prb_mng_ptr->buildProblem("ELASTIC_MECHANICS", true);
  // partition
  CHKERR prb_mng_ptr->partitionProblem("ELASTIC_MECHANICS");
  CHKERR prb_mng_ptr->partitionFiniteElements("ELASTIC_MECHANICS", false, 0,
                                              m_field.get_comm_size());
  // what are ghost nodes, see Petsc Manual
  CHKERR prb_mng_ptr->partitionGhostDofs("ELASTIC_MECHANICS");
 
  MoFEMFunctionReturn(0);
}

◆ solveSystem()

MoFEMErrorCode ArcLengthInterfaceExample::solveSystem ( )

protected

[Solve]

Examples: mofem/tutorials/cor-12_cohesive_interface/arc_length_interface.cpp.

Definition at line 969 of file arc_length_interface.cpp.

                                                      {
  MoFEMFunctionBegin;
 
  double min_arc_length_step = std::numeric_limits<double>::epsilon();
  double max_arc_length_step = std::numeric_limits<double>::max();
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR, "-gamma", &gamma,
                             PETSC_NULLPTR);
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR,
                             "-min_arc_length_step", &min_arc_length_step,
                             PETSC_NULLPTR);
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, PETSC_NULLPTR,
                             "-max_arc_length_step", &max_arc_length_step,
                             PETSC_NULLPTR);
 
  if (min_arc_length_step <= 0.)
    SETERRQ(PETSC_COMM_WORLD, MOFEM_INVALID_DATA,
            "-min_arc_length_step must be positive");
  if (max_arc_length_step <= 0.)
    SETERRQ(PETSC_COMM_WORLD, MOFEM_INVALID_DATA,
            "-max_arc_length_step must be positive");
  if (min_arc_length_step > max_arc_length_step)
    SETERRQ(PETSC_COMM_WORLD, MOFEM_INVALID_DATA,
            "-min_arc_length_step must not exceed -max_arc_length_step");
 
  auto &m_field = *mField;
  auto &step = *stepPtr;
  auto &step_size = *stepSizePtr;
 
  auto &cohesive_elements = *cohesiveElementsPtr;
  auto &post_proc = *postProcPtr;
 
  bool converged_state = false;
  for (; step < maxSteps; step++) {
 
    if (step == 1) {
      MOFEM_LOG_C("ARC_LENGTH", Sev::inform, "Load Step %d step_size = %6.4e",
                  step, step_size);
      CHKERR arcCtx->setS(step_size);
      CHKERR arcCtx->setAlphaBeta(0, 1);
      CHKERR VecCopy(d, arcCtx->x0);
      double dlambda;
      CHKERR arcMethodPtr->calculate_init_dlambda(&dlambda);
      CHKERR arcMethodPtr->set_dlambda_to_x(dRaw, dlambda);
    } else if (step == 2) {
      CHKERR arcCtx->setAlphaBeta(1, 0);
      CHKERR arcMethodPtr->calculate_dx_and_dlambda(dRaw);
      CHKERR arcMethodPtr->calculate_lambda_int(step_size);
      CHKERR arcCtx->setS(step_size);
      double dlambda = arcCtx->dLambda;
      double dx_nrm;
      CHKERR VecNorm(arcCtx->dx, NORM_2, &dx_nrm);
      MOFEM_LOG_C("ARC_LENGTH", Sev::inform,
                  "Load Step %d step_size = %6.4e dlambda0 = %6.4e "
                  "dx_nrm = %6.4e dx2 = %6.4e",
                  step, step_size, dlambda, dx_nrm, arcCtx->dx2);
      CHKERR VecCopy(d, arcCtx->x0);
      CHKERR VecAXPY(d, 1., arcCtx->dx);
      CHKERR arcMethodPtr->set_dlambda_to_x(dRaw, dlambda);
    } else {
      CHKERR arcMethodPtr->calculate_dx_and_dlambda(dRaw);
      CHKERR arcMethodPtr->calculate_lambda_int(step_size);
      // step_size0_1/step_size0 = step_stize1/step_size
      // step_size0_1 = step_size0*(step_stize1/step_size)
      step_size *= reduction;
      CHKERR arcCtx->setS(step_size);
      double dlambda = reduction * arcCtx->dLambda;
      CHKERR VecScale(arcCtx->dx, reduction);
      double dx_nrm;
      CHKERR VecNorm(arcCtx->dx, NORM_2, &dx_nrm);
      MOFEM_LOG_C("ARC_LENGTH", Sev::inform,
                  "Load Step %d step_size = %6.4e dlambda0 = %6.4e "
                  "dx_nrm = %6.4e dx2 = %6.4e",
                  step, step_size, dlambda, dx_nrm, arcCtx->dx2);
      CHKERR VecCopy(d, arcCtx->x0);
      CHKERR VecAXPY(d, 1., arcCtx->dx);
      CHKERR arcMethodPtr->set_dlambda_to_x(dRaw, dlambda);
    }
 
    CHKERR SNESSolve(snes, PETSC_NULLPTR, d);
 
    int its;
    CHKERR SNESGetIterationNumber(snes, &its);
    MOFEM_LOG_C("ARC_LENGTH", Sev::inform, "number of Newton iterations = %d",
                its);
 
    SNESConvergedReason reason;
    CHKERR SNESGetConvergedReason(snes, &reason);
 
    if (reason <= 0) {
      CHKERR arcCtx->setAlphaBeta(1, 0);
      reduction = 0.1;
      converged_state = false;
      continue;
    } else {
      if (step > 1 && converged_state) {
        reduction = pow((double)itsD / (double)(its + 1), gamma);
        reduction =
            std::clamp(reduction, min_arc_length_step / std::abs(step_size),
                       max_arc_length_step / std::abs(step_size));
        MOFEM_LOG_C("ARC_LENGTH", Sev::inform, "reduction step_size = %6.4e",
                    reduction);
      }
 
      // Save data on mesh
      CHKERR m_field.getInterface<VecManager>()->setGlobalGhostVector(
          "ELASTIC_MECHANICS", COL, d, INSERT_VALUES, SCATTER_REVERSE);
      // Commit cohesive history only after SNES convergence.
      CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", "INTERFACE",
                                          cohesive_elements.getFeHistory(), 0,
                                          m_field.get_comm_size());
      CHKERR arcMethodPtr->remove_damaged_prisms_nodes();
      CHKERR m_field.getInterface<VecManager>()->setOtherGlobalGhostVector(
          "ELASTIC_MECHANICS", "DISPLACEMENT", "X0_DISPLACEMENT", COL,
          arcCtx->x0, INSERT_VALUES, SCATTER_REVERSE);
      converged_state = true;
    }
    //
    if (reason > 0) {
      CHKERR evaluateFieldAtPoint(step);
      CHKERR evaluateReactions(step);
    }
 
    if (step % 1 == 0) {
 
      CHKERR m_field.loop_finite_elements("ELASTIC_MECHANICS", "ELASTIC",
                                          post_proc);
      std::ostringstream ss;
      ss << "out_values_" << step << ".h5m";
      CHKERR post_proc.writeFile(ss.str().c_str());
    }
  }
 
  MoFEMFunctionReturn(0);
}