CPSolvers.cpp

Go to the documentation of this file.

/** \file CPSolvers.cpp
 */
 
/* This file is part of MoFEM.
 * MoFEM is free software: you can redistribute it and/or modify it under
 * the terms of the GNU Lesser General Public License as published by the
 * Free Software Foundation, either version 3 of the License, or (at your
 * option) any later version.
 *
 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 * License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
 
#ifdef WITH_TETGEN
 
#include <tetgen.h>
#ifdef REAL
#undef REAL
#endif
 
#endif
 
#include <MoFEM.hpp>
using namespace MoFEM;
#include <BasicFiniteElements.hpp>
#include <Mortar.hpp>
#include <NeoHookean.hpp>
#include <Hooke.hpp>
#include <CrackFrontElement.hpp>
#include <ComplexConstArea.hpp>
#include <MWLS.hpp>
#include <GriffithForceElement.hpp>
#include <VolumeLengthQuality.hpp>
#include <CrackPropagation.hpp>
#include <AnalyticalFun.hpp>
#include <CPMeshCut.hpp>
#include <CPSolvers.hpp>
 
#define SolveFunctionBegin                                                     \
  MoFEMFunctionBegin;                                                          \
  MOFEM_LOG_FUNCTION();
 
namespace FractureMechanics {
 
MoFEMErrorCode
CPSolvers::query_interface(boost::typeindex::type_index type_index,
                           UnknownInterface **iface) const {
  *iface = const_cast<CPSolvers *>(this);
  return 0;
}
 
CPSolvers::CPSolvers(CrackPropagation &cp) : cP(cp) {
 
  if (!LogManager::checkIfChannelExist("CPSolverWorld")) {
    auto core_log = logging::core::get();
 
    core_log->add_sink(
        LogManager::createSink(LogManager::getStrmWorld(), "CPSolverWorld"));
    core_log->add_sink(
        LogManager::createSink(LogManager::getStrmSync(), "CPSolverSync"));
    core_log->add_sink(
        LogManager::createSink(LogManager::getStrmSelf(), "CPSolverSelf"));
 
    LogManager::setLog("CPSolverWorld");
    LogManager::setLog("CPSolverSync");
    LogManager::setLog("CPSolverSelf");
 
    MOFEM_LOG_TAG("CPSolverWorld", "CPSolve");
    MOFEM_LOG_TAG("CPSolverSync", "CPSolve");
    MOFEM_LOG_TAG("CPSolverSelf", "CPSolve");
  }
 
  MOFEM_LOG("CPSolverWorld", Sev::noisy) << "CPSolve created";
}
 
MoFEMErrorCode CPSolvers::getOptions() {
  MoFEMFunctionBegin;
 
  // CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Solvers options", "none");
  // ierr = PetscOptionsEnd();
  // CHKERRG(ierr);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CPSolvers::solveElastic(DM dm_elastic, Vec init,
                                       const double load_factor) {
  MoFEM::Interface &m_field = cP.mField;
  MoFEMFunctionBegin;
 
  // Create matrix
  auto m_elastic = smartCreateDMMatrix(dm_elastic);
  // Create vector
  auto q_elastic = smartCreateDMVector(dm_elastic);
  auto f_elastic = smartVectorDuplicate(q_elastic);
 
  const MoFEM::Problem *problem_ptr;
  CHKERR DMMoFEMGetProblemPtr(dm_elastic, &problem_ptr);
 
  auto solve_elastic_problem = [&](auto init_x) {
    MoFEMFunctionBegin;
 
    /// Clear finite elements added SNES solver
    smartGetDMSnesCtx(dm_elastic)->clearLoops();
 
    boost::shared_ptr<ArcLengthCtx> arc_ctx;
    if (problem_ptr->getName() == "ELASTIC")
      arc_ctx = cP.getArcCtx();
    if (problem_ptr->getName() == "EIGEN_ELASTIC")
      arc_ctx = cP.getEigenArcCtx();
 
    // Set operators for elasticity analysis
    boost::shared_ptr<SimpleArcLengthControl> arc_method(
        new SimpleArcLengthControl(arc_ctx));
 
    // Add operators to DM SNES
    CHKERR cP.addElasticFEInstancesToSnes(dm_elastic, PETSC_NULL, PETSC_NULL,
                                          PETSC_NULL, arc_method, arc_ctx,
                                          VERBOSE, false);
 
    CHKERR arc_ctx->setAlphaBeta(0, 1);
    CHKERR arc_ctx->setS(0);
    arc_ctx->getFieldData() = load_factor;
 
    CHKERR cP.getArcLengthDof();
 
    if (m_field.check_problem("BC_PROBLEM")) {
      // solve problem to find dofs for analytical displacements
      cP.getAnalyticalDirichletBc()->snes_B = m_elastic;
      cP.getAnalyticalDirichletBc()->snes_x = q_elastic;
      cP.getAnalyticalDirichletBc()->snes_f = f_elastic;
      // solve for Dirichlet bc dofs
      AnalyticalDirichletBC analytical_bc(m_field);
      CHKERR analytical_bc.setUpProblem(m_field, "BC_PROBLEM");
      boost::shared_ptr<AnalyticalDisp> testing_function =
          boost::shared_ptr<AnalyticalDisp>(new AnalyticalDisp());
      // EntityHandle fe_meshset  =
      // m_field.get_finite_element_meshset("BC_FE");
      CHKERR analytical_bc.setApproxOps(m_field, "SPATIAL_POSITION",
                                        testing_function, 0,
                                        "MESH_NODE_POSITIONS");
 
      analytical_bc.approxField.getLoopFeApprox().addToRule = 1;
      CHKERR analytical_bc.solveProblem(m_field, "BC_PROBLEM", "BC_FE",
                                        *cP.getAnalyticalDirichletBc());
 
      CHKERR analytical_bc.destroyProblem();
    }
 
    // set vector values from field data
    CHKERR DMoFEMMeshToLocalVector(dm_elastic, q_elastic, INSERT_VALUES,
                                   SCATTER_FORWARD);
 
    auto snes = createSNES(m_field.get_comm());
    Mat shell_m;
 
    // for(int ll = 0;ll!=cP.getNbLoadSteps();ll++) {
    CHKERR VecCopy(q_elastic, init_x);
    CHKERR cP.solveElasticDM(dm_elastic, snes, m_elastic, q_elastic, f_elastic,
                             true, &shell_m);
    CHKERR DMoFEMMeshToLocalVector(dm_elastic, q_elastic, INSERT_VALUES,
                                   SCATTER_REVERSE);
 
    CHKERR MatDestroy(&shell_m);
 
    MoFEMFunctionReturn(0);
  };
 
  // If some flag is set, we solve eigen problem, i.e. calculate eigen
  // displacements
  if (problem_ptr->getName() == "EIGEN_ELASTIC")
    CHKERR cP.assembleElasticDM(cP.mwlsEigenStressTagName);
  if (problem_ptr->getName() == "ELASTIC")
    CHKERR cP.assembleElasticDM(cP.mwlsStressTagName);
 
  // + and ops which will close crack
  CHKERR solve_elastic_problem(init);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CPSolvers::calculateGc(DM dm_material_forces,
                                      DM dm_surface_projection,
                                      DM dm_crack_srf_area,
                                      const std::vector<int> &surface_ids) {
  MoFEM::Interface &m_field = cP.mField;
  MoFEMFunctionBegin;
 
  // Create vector of material forces and set operators
  Vec q_material, f_material;
  CHKERR DMCreateGlobalVector(dm_material_forces, &q_material);
  CHKERR VecDuplicate(q_material, &f_material);
  Vec f_material_proj, f_griffith, f_griffith_proj;
  CHKERR VecDuplicate(f_material, &f_material_proj);
  CHKERR VecDuplicate(f_material, &f_griffith);
  CHKERR VecDuplicate(f_griffith, &f_griffith_proj);
  Vec f_lambda;
  {
    const MoFEM::Problem *prb_ptr;
    CHKERR DMMoFEMGetProblemPtr(dm_crack_srf_area, &prb_ptr);
    // Need to use directly MoFEM interface to create vector, matrix for non
    // square matrix
    CHKERR m_field.getInterface<VecManager>()->vecCreateGhost(
        prb_ptr->getName(), ROW, &f_lambda);
  }
 
  // Calculate projection matrices
  CHKERR cP.calculateSurfaceProjectionMatrix(
      dm_surface_projection, dm_material_forces, surface_ids, VERBOSE, false);
  CHKERR cP.calculateFrontProjectionMatrix(dm_crack_srf_area,
                                           dm_material_forces, VERBOSE, false);
  // Calculate griffith force vector (material resistance)
  CHKERR cP.calculateGriffithForce(dm_material_forces, cP.gC ? cP.gC : 1,
                                   f_griffith, VERBOSE, false);
  CHKERR cP.projectGriffithForce(dm_material_forces, f_griffith,
                                 f_griffith_proj, VERBOSE, false);
  // Calculate material forces
  CHKERR DMoFEMMeshToLocalVector(dm_material_forces, q_material, INSERT_VALUES,
                                 SCATTER_FORWARD);
 
  CHKERR cP.calculateMaterialForcesDM(dm_material_forces, q_material,
                                      f_material, VERBOSE, false);
 
  // project material forces
  CHKERR cP.projectMaterialForcesDM(dm_material_forces, f_material,
                                    f_material_proj, VERBOSE, false);
  // calculate griffith energy
  CHKERR cP.calculateReleaseEnergy(dm_crack_srf_area, f_material_proj,
                                   f_griffith_proj, f_lambda, cP.gC, VERBOSE,
                                   false);
 
  CHKERR VecDestroy(&q_material);
  CHKERR VecDestroy(&f_material);
  CHKERR VecDestroy(&f_material_proj);
  CHKERR VecDestroy(&f_griffith);
  CHKERR VecDestroy(&f_griffith_proj);
  CHKERR VecDestroy(&f_lambda);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
CPSolvers::solvePropagation(DM dm_crack_propagation, DM dm_elastic,
                            DM dm_material, DM dm_material_forces,
                            DM dm_surface_projection, DM dm_crack_srf_area,
                            const std::vector<int> &surface_ids,
                            int cut_mesh_it, const bool set_cut_surface) {
  MoFEM::Interface &m_field = cP.mField;
  MoFEMFunctionBegin;
 
  const MoFEM::Problem *problem_elastic_ptr;
  CHKERR DMMoFEMGetProblemPtr(dm_elastic, &problem_elastic_ptr);
  cP.arcCtx = boost::make_shared<ArcLengthCtx>(
      cP.mField, problem_elastic_ptr->getName(), "LAMBDA_ARC_LENGTH");
  auto arc_snes_ctx = boost::make_shared<CrackPropagation::ArcLengthSnesCtx>(
      cP.mField, problem_elastic_ptr->getName(), cP.arcCtx);
  CHKERR DMMoFEMSetSnesCtx(dm_elastic, arc_snes_ctx);
 
  // set finite element instances and user data operators on instances
  CHKERR cP.assembleElasticDM(VERBOSE, false);
  MOFEM_LOG("CPSolverWorld", Sev::noisy)
      << "Solve ELASTIC problem in solvePropagation";
  CHKERR solveElastic(dm_elastic, cP.getArcCtx()->x0,
                      cP.getArcCtx()->getFieldData());
 
  CHKERR cP.assembleMaterialForcesDM(PETSC_NULL);
  CHKERR cP.assembleSmootherForcesDM(PETSC_NULL, surface_ids, VERBOSE, false);
  CHKERR cP.assembleCouplingForcesDM(PETSC_NULL, VERBOSE, false);
 
  struct StepAdaptivity {
 
    CrackPropagation &cP;
 
    int itsD;
    double gAmma;
    double minStep;
    double maxStep;
    PetscBool minStepFlg;
    PetscBool maxStepFlg;
 
    StepAdaptivity(CrackPropagation &cp, int its_d, double gamma)
        : cP(cp), itsD(its_d), gAmma(gamma), minStep(0), maxStep(0) {
      ierr = getOptions();
      CHKERRABORT(PETSC_COMM_WORLD, ierr);
    }
 
    MoFEMErrorCode getOptions() {
      SolveFunctionBegin;
      ierr =
          PetscOptionsBegin(PETSC_COMM_WORLD, "", "Fracture options", "none");
      CHKERRQ(ierr);
      {
        CHKERR PetscOptionsInt("-adapt_step_its_d",
                               "number of desired iterations", "", itsD, &itsD,
                               PETSC_NULL);
        CHKERR PetscOptionsScalar("-adapt_step_min_s", "minimal setp", "",
                                  minStep, &minStep, &minStepFlg);
        CHKERR PetscOptionsScalar("-adapt_step_max_s", "maximal setp", "",
                                  maxStep, &maxStep, &maxStepFlg);
      }
      ierr = PetscOptionsEnd();
      CHKERRQ(ierr);
 
      MOFEM_LOG_C("CPSolverWorld", Sev::inform,
                  "### Input parameter: -adapt_step_its_d %d", itsD);
      MOFEM_LOG_C("CPSolverWorld", Sev::inform,
                  "### Input parameter: -adapt_step_min_s %6.4e", minStep);
      MOFEM_LOG_C("CPSolverWorld", Sev::inform,
                  "### Input parameter: -adapt_step_max_s %6.4e", maxStep);
 
      if (minStepFlg == PETSC_TRUE && maxStepFlg == PETSC_TRUE)
        if (minStep > maxStep) {
          SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
                  "Minimal step size cannot be bigger than maximal step size");
        }
      MoFEMFunctionReturn(0);
    }
 
    MoFEMErrorCode operator()(DM dm, SNES snes,
                              boost::shared_ptr<ArcLengthCtx> arc_ctx,
                              double &s) const {
      SolveFunctionBegin;
 
      auto min_max = [this](const double s) {
        double new_s = s;
        if (minStepFlg == PETSC_TRUE) {
          new_s = std::max(new_s, minStep);
        }
        if (maxStepFlg == PETSC_TRUE) {
          new_s = std::min(new_s, maxStep);
        }
        return new_s;
      };
 
      SNESConvergedReason reason;
      CHKERR SNESGetConvergedReason(snes, &reason);
      if (reason < 0) {
        s = min_max(0.25 * s);
        MOFEM_LOG_C("CPSolverWorld", Sev::warning,
                    "*** Diverged ***"
                    "Setting reduced load step arc_s = %3.4g",
                    s);
        CHKERR DMoFEMMeshToLocalVector(dm, arc_ctx->x0, INSERT_VALUES,
                                       SCATTER_REVERSE);
 
      } else {
        int its;
        CHKERR SNESGetIterationNumber(snes, &its);
        s = min_max(s * pow((double)itsD / (double)(its + 1), gAmma));
      }
      MoFEMFunctionReturn(0);
    }
  };
 
  struct SmoothingAlphaAdaptivity {
 
    CrackPropagation &cP;
    double gAmma;
    double dAlpha;
    double minAlpha;
    double maxAlpha;
    double incrementAlpha;
 
    SmoothingAlphaAdaptivity(CrackPropagation &cp, const double gamma,
                             const double d_alpha, const double m_alpha)
        : cP(cp), gAmma(gamma), dAlpha(d_alpha), minAlpha(m_alpha) {
      ierr = getOptions();
      CHKERRABORT(PETSC_COMM_WORLD, ierr);
    }
 
    MoFEMErrorCode getOptions() {
      MoFEMFunctionBegin;
      ierr =
          PetscOptionsBegin(PETSC_COMM_WORLD, "", "Fracture options", "none");
      CHKERRQ(ierr);
      CHKERR PetscOptionsScalar("-adapt_min_smoother_alpha",
                                "minimal smoother alpha", "", minAlpha,
                                &minAlpha, PETSC_NULL);
      maxAlpha = 1e3;
      CHKERR PetscOptionsScalar("-adapt_max_smoother_alpha",
                                "minimal smoother alpha", "", maxAlpha,
                                &maxAlpha, PETSC_NULL);
      incrementAlpha = 10;
      CHKERR PetscOptionsScalar("-adapt_incitement_smoother_alpha",
                                "minimal smoother alpha", "", incrementAlpha,
                                &incrementAlpha, PETSC_NULL);
 
      CHKERR PetscOptionsScalar(
          "-adapt_desired_alpha",
          "Desired alpha (allowed error for gc calculation)", "", dAlpha,
          &dAlpha, PETSC_NULL);
 
      MOFEM_LOG_C("CPSolverWorld", Sev::inform,
                  "### Input parameter: -adapt_min_smoother_alpha %6.4e",
                  minAlpha);
      MOFEM_LOG_C("CPSolverWorld", Sev::inform,
                  "### Input parameter: -adapt_max_smoother_alpha %6.4e",
                  maxAlpha);
      MOFEM_LOG_C("CPSolverWorld", Sev::inform,
                  "### Input parameter: -adapt_desired_alpha %6.4e", dAlpha);
 
      ierr = PetscOptionsEnd();
 
      MoFEMFunctionReturn(0);
    }
 
    MoFEMErrorCode operator()(DM dm, SNES snes, Vec q) const {
      MoFEMFunctionBegin;
 
      if (cP.smootherFe->smootherData.sTabilised) {
        SNESConvergedReason reason;
        CHKERR SNESGetConvergedReason(snes, &reason);
        if (reason < 0) {
          const double old_alpha = cP.smootherAlpha;
          cP.smootherAlpha =
              std::min(incrementAlpha * cP.smootherAlpha, maxAlpha);
          MOFEM_LOG_C("CPSolverWorld", Sev::inform,
                      "Increase smoothing alpha = %6.5g (old %6.5e)",
                      cP.smootherAlpha, old_alpha);
        } else {
          auto f_smoothing = smartVectorDuplicate(q);
          auto f_griffith = smartVectorDuplicate(q);
 
          CHKERR cP.calculateGriffithForce(dm, cP.gC ? cP.gC : 1, f_griffith,
                                           VERBOSE, false);
          // Calculate material forces
          CHKERR DMoFEMMeshToLocalVector(dm, q, INSERT_VALUES, SCATTER_FORWARD);
          CHKERR cP.calculateSmoothingForcesDM(dm, q, f_smoothing, VERBOSE,
                                               false);
          CHKERR cP.calculateSmoothingForceFactor(VERBOSE, false);
 
          double max_alpha = 0;
          for (const auto &kv : cP.mapSmoothingForceFactor)
            max_alpha = std::max(max_alpha, kv.second);
 
          const double old_alpha = cP.smootherAlpha;
          cP.smootherAlpha *= pow(dAlpha / max_alpha, gAmma);
          cP.smootherAlpha = std::max(cP.smootherAlpha, minAlpha);
          MOFEM_LOG_C("CPSolverWorld", Sev::inform,
                      "Set smoothing alpha = %6.5g (old %6.5e)",
                      cP.smootherAlpha, old_alpha);
        }
 
        cP.volumeLengthAdouble->aLpha = cP.smootherAlpha;
        cP.volumeLengthDouble->aLpha = cP.smootherAlpha;
        for (auto &kv : cP.surfaceConstrain)
          kv.second->aLpha = cP.smootherAlpha;
 
        for (auto &kv : cP.edgeConstrains)
          kv.second->aLpha = cP.smootherAlpha;
 
        cP.bothSidesConstrains->aLpha = cP.smootherAlpha;
        cP.bothSidesContactConstrains->aLpha = cP.smootherAlpha;
      }
 
      MoFEMFunctionReturn(0);
    }
  };
 
  const MoFEM::Problem *problem_cp_ptr;
  CHKERR DMMoFEMGetProblemPtr(dm_crack_propagation, &problem_cp_ptr);
 
  auto cp_arc_ctx = boost::make_shared<ArcLengthCtx>(
      m_field, problem_cp_ptr->getName(), "LAMBDA_ARC_LENGTH");
 
  CHKERR cP.mField.getInterface<VecManager>()->vecCreateGhost(
      problem_cp_ptr->getName(), COL, cp_arc_ctx->x0);
  cp_arc_ctx->dx = smartVectorDuplicate(cp_arc_ctx->x0);
  arc_snes_ctx = boost::make_shared<CrackPropagation::ArcLengthSnesCtx>(
      m_field, problem_cp_ptr->getName(), cp_arc_ctx);
  CHKERR DMMoFEMSetSnesCtx(dm_crack_propagation, arc_snes_ctx);
 
  boost::shared_ptr<FEMethod> arc_method =
      cP.getFrontArcLengthControl(cp_arc_ctx);
 
  // Add operators to DM SNES
 
  CHKERR cP.addPropagationFEInstancesToSnes(dm_crack_propagation, arc_method,
                                            cp_arc_ctx, surface_ids, 1, false);
 
  {
    // create vectors and matrix
    auto m = smartCreateDMMatrix(dm_crack_propagation);
    auto f = smartCreateDMVector(dm_crack_propagation);
    auto q = smartVectorDuplicate(f);
 
    // Create snes
    auto snes_crack_propagation = createSNES(m_field.get_comm());
    CHKERR SNESAppendOptionsPrefix(snes_crack_propagation, "propagation_");
 
    for (int ll = 0; ll != cP.getNbLoadSteps(); ll++) {
 
      int run_step = cP.getNbLoadSteps() * cut_mesh_it + ll;
      int step = cP.startStep++;
 
      CHKERR PetscPrintf(PETSC_COMM_WORLD, "\n\nLoad step %d ( %d ) [ %d ]\n\n",
                         ll, run_step, step);
 
      // set vector values from mesh data
      CHKERR cP.zeroLambdaFields();
      CHKERR DMoFEMMeshToLocalVector(dm_crack_propagation, q, INSERT_VALUES,
                                     SCATTER_FORWARD);
 
      CHKERR VecCopy(q, cp_arc_ctx->x0);
      CHKERR cp_arc_ctx->setAlphaBeta(cP.arcAlpha, cP.arcBeta);
 
      auto solve_problem = [&]() {
        MoFEMFunctionBegin;
 
        const double arc_s = cP.arcS;
 
        cP.arcS *= cP.fractionOfFixedNodes;
        CHKERR cp_arc_ctx->setS(cP.arcS);
 
        CHKERR cP.solvePropagationDM(dm_crack_propagation, dm_elastic,
                                     snes_crack_propagation, m, q, f);
        cP.arcS = arc_s;
        CHKERR cp_arc_ctx->setS(cP.arcS);
 
        // Adapt step
        CHKERR StepAdaptivity(cP, 6, 0.5)(
            dm_crack_propagation, snes_crack_propagation, cp_arc_ctx, cP.arcS);
        CHKERR SmoothingAlphaAdaptivity(cP, 0.25, 1e-2,
                                        cP.initialSmootherAlpha * 1e-3)(
            dm_crack_propagation, snes_crack_propagation, q);
        MoFEMFunctionReturn(0);
      };
 
      CHKERR solve_problem();
 
      auto if_not_converged = [&](Vec q) {
        MoFEMFunctionBegin;
        SNESConvergedReason reason;
        CHKERR SNESGetConvergedReason(snes_crack_propagation, &reason);
        if (reason < 0) {
          CHKERR VecCopy(cp_arc_ctx->x0, q);
          CHKERR DMoFEMMeshToLocalVector(dm_crack_propagation, q, INSERT_VALUES,
                                         SCATTER_FORWARD);
        }
        MoFEMFunctionReturn(0);
      };
      CHKERR if_not_converged(q);
 
      auto scatter_forward_vec = [](Vec q) {
        MoFEMFunctionBegin;
        CHKERR VecGhostUpdateBegin(q, INSERT_VALUES, SCATTER_FORWARD);
        CHKERR VecGhostUpdateEnd(q, INSERT_VALUES, SCATTER_FORWARD);
        MoFEMFunctionReturn(0);
      };
      CHKERR scatter_forward_vec(q);
 
      // Save material displacements on the mesh
      CHKERR VecAYPX(cp_arc_ctx->x0, -1, q);
      CHKERR scatter_forward_vec(cP.getArcCtx()->x0);
 
      CrackPropagation::PostProcVertexMethod ent_front_disp(
          cP.mField, cp_arc_ctx->x0, "W");
      CHKERR cP.mField.loop_dofs(problem_cp_ptr->getName(),
                                 "MESH_NODE_POSITIONS", ROW, ent_front_disp, 0,
                                 cP.mField.get_comm_rank());
 
      // Displace mesh
      CHKERR cP.setCoordsFromField("MESH_NODE_POSITIONS");
 
      // calculate energy
      SNESConvergedReason reason;
      CHKERR SNESGetConvergedReason(snes_crack_propagation, &reason);
      std::string step_msg;
      if (reason < 0) {
        step_msg = "Not Converged Propagation step " +
                   boost::lexical_cast<std::string>(step) + " ( " +
                   boost::lexical_cast<std::string>(ll) + " )";
      } else {
        step_msg = "Propagation step " +
                   boost::lexical_cast<std::string>(step) + " ( " +
                   boost::lexical_cast<std::string>(ll) + " )";
      }
      CHKERR cP.calculateElasticEnergy(dm_elastic, step_msg);
      // calculate griffith energy
      CHKERR calculateGc(dm_material_forces, dm_surface_projection,
                         dm_crack_srf_area, surface_ids);
      CHKERR cP.updateMaterialFixedNode(false, true, false);
 
      // post-processing
      CHKERR cP.postProcessDM(dm_elastic, step, "ELASTIC",
                              false || run_step == 0);
      // save crack front
      CHKERR cP.savePositionsOnCrackFrontDM(dm_elastic, q, QUIET, false);
 
      // Save positions
      auto save_positions_on_mesh_tags_from_field =
          [this, problem_cp_ptr, &q](const std::string field_name) {
            MoFEMFunctionBegin;
            CrackPropagation::PostProcVertexMethod ent_post_proc(cP.mField, q,
                                                                 field_name);
            CHKERR cP.mField.loop_dofs(problem_cp_ptr->getName(), field_name,
                                       ROW, ent_post_proc, 0,
                                       cP.mField.get_comm_size());
            MoFEMFunctionReturn(0);
          };
      CHKERR save_positions_on_mesh_tags_from_field("MESH_NODE_POSITIONS");
      CHKERR save_positions_on_mesh_tags_from_field("SPATIAL_POSITION");
 
      // write state
      CHKERR cP.writeCrackFont(cP.mapBitLevel["spatial_domain"], step);
 
      auto save_restart_file = [&]() {
        auto bit_interface = cP.mField.getInterface<BitRefManager>();
 
        MoFEMFunctionBegin;
 
        auto bit = BitRefLevel().set(BITREFLEVEL_SIZE - 1);
        bit.set(BITREFLEVEL_SIZE - 2);
        auto meshset_ptr = get_temp_meshset_ptr(m_field.get_moab());
 
        Range ents_on_last_bits;
        CHKERR bit_interface->getEntitiesByRefLevel(bit, BitRefLevel().set(),
                                                    ents_on_last_bits);
        CHKERR cP.mField.get_moab().add_entities(*meshset_ptr,
                                                 ents_on_last_bits);
 
        Range ents_spatial_domain;
        CHKERR bit_interface->getEntitiesByRefLevel(
            cP.mapBitLevel.at("spatial_domain"), BitRefLevel().set(),
            ents_spatial_domain);
 
        Range entities_on_skin;
        entities_on_skin.merge(cP.bodySkin);
        entities_on_skin.merge(cP.crackFaces);
        entities_on_skin.merge(cP.crackFront);
 
        entities_on_skin = intersect(ents_spatial_domain, entities_on_skin);
        CHKERR cP.mField.get_moab().add_entities(*meshset_ptr,
                                                 entities_on_skin);
        CHKERR cP.mField.get_moab().add_entities(*meshset_ptr,
                                                 ents_spatial_domain);
 
        Range entities_bit;
        entities_bit.merge(ents_on_last_bits);
        entities_bit.merge(ents_spatial_domain);
 
        auto add_field = [&](auto field_name, auto &tags_list) {
          MoFEMFunctionBegin;
          auto field_ptr = cP.mField.get_field_structure(field_name);
          auto meshset = field_ptr->getMeshset();
          CHKERR cP.mField.get_moab().add_entities(*meshset_ptr, &meshset, 1);
 
          Tag th_field_id;
          CHKERR cP.mField.get_moab().tag_get_handle("_FieldId", th_field_id);
          Tag th_field_base;
          CHKERR cP.mField.get_moab().tag_get_handle("_FieldBase",
                                                     th_field_base);
          Tag th_field_space;
          CHKERR cP.mField.get_moab().tag_get_handle("_FieldSpace",
                                                     th_field_space);
          Tag th_field_name;
          CHKERR cP.mField.get_moab().tag_get_handle("_FieldName",
                                                     th_field_name);
          Tag th_field_name_data_name_prefix;
          CHKERR cP.mField.get_moab().tag_get_handle(
              "_FieldName_DataNamePrefix", th_field_name_data_name_prefix);
 
          tags_list.push_back(th_field_id);
          tags_list.push_back(th_field_base);
          tags_list.push_back(th_field_space);
          tags_list.push_back(th_field_name);
          tags_list.push_back(th_field_name_data_name_prefix);
          tags_list.push_back(field_ptr->th_FieldDataVerts);
          tags_list.push_back(field_ptr->th_FieldData);
          tags_list.push_back(field_ptr->th_AppOrder);
          tags_list.push_back(field_ptr->th_FieldRank);
          MoFEMFunctionReturn(0);
        };
 
        auto add_fe = [&](auto fe_name, auto &tags_list) {
          MoFEMFunctionBegin;
          auto arc_length_finite_element =
              cP.mField.get_finite_element_meshset(fe_name);
          CHKERR cP.mField.get_moab().add_entities(
              *meshset_ptr, &arc_length_finite_element, 1);
          Tag th_fe_id;
          CHKERR cP.mField.get_moab().tag_get_handle("_FEId", th_fe_id);
          Tag th_fe_name;
          CHKERR cP.mField.get_moab().tag_get_handle("_FEName", th_fe_name);
          Tag th_fe_id_col, th_fe_id_row, th_fe_id_data;
          CHKERR cP.mField.get_moab().tag_get_handle("_FEIdCol", th_fe_id_col);
          CHKERR cP.mField.get_moab().tag_get_handle("_FEIdRow", th_fe_id_row);
          CHKERR cP.mField.get_moab().tag_get_handle("_FEIdData",
                                                     th_fe_id_data);
          tags_list.push_back(th_fe_id);
          tags_list.push_back(th_fe_name);
          tags_list.push_back(th_fe_id_col);
          tags_list.push_back(th_fe_id_row);
          tags_list.push_back(th_fe_id_data);
          MoFEMFunctionReturn(0);
        };
 
        auto arc_field_meshset =
            cP.mField.get_field_meshset("LAMBDA_ARC_LENGTH");
        auto arc_length_finite_element =
            cP.mField.get_finite_element_meshset("ARC_LENGTH");
        CHKERR cP.mField.get_moab().add_entities(*meshset_ptr,
                                                 &arc_field_meshset, 1);
        CHKERR cP.mField.get_moab().add_entities(*meshset_ptr,
                                                 &arc_length_finite_element, 1);
 
        auto get_tags_list = [&]() {
          Tag th_griffith_force;
          CHKERR cP.mField.get_moab().tag_get_handle("GRIFFITH_FORCE",
                                                     th_griffith_force);
          Tag th_griffith_force_projected;
          CHKERR cP.mField.get_moab().tag_get_handle(
              "GRIFFITH_FORCE_PROJECTED", th_griffith_force_projected);
          Tag th_w;
          CHKERR cP.mField.get_moab().tag_get_handle("W", th_w);
          Tag th_material_force;
          CHKERR cP.mField.get_moab().tag_get_handle("MATERIAL_FORCE",
                                                     th_material_force);
          Tag th_interface_side;
          CHKERR cP.mField.get_moab().tag_get_handle("INTERFACE_SIDE",
                                                     th_interface_side);
          Tag th_org_pos;
          CHKERR cP.mField.get_moab().tag_get_handle("ORG_POSITION",
                                                     th_org_pos);
          Tag th_version;
          CHKERR cP.mField.get_moab().tag_get_handle("MOFEM_VERSION",
                                                     th_version);
 
          std::vector<Tag> tags_list;
          tags_list.push_back(bit_interface->get_th_RefBitLevel());
 
          tags_list.push_back(th_version);
          tags_list.push_back(th_griffith_force);
          tags_list.push_back(th_griffith_force_projected);
          tags_list.push_back(th_w);
          tags_list.push_back(th_material_force);
          tags_list.push_back(th_interface_side);
          tags_list.push_back(th_org_pos);
          return tags_list;
        };
 
        auto add_meshsets = [&](Range bit, auto &tags_list) {
          auto meshsets_mng = cP.mField.getInterface<MeshsetsManager>();
          std::vector<boost::shared_ptr<TempMeshset>> meshses_tmp_list;
          auto &list = meshsets_mng->getMeshsetsMultindex();
          for (auto &m : list) {
            meshses_tmp_list.push_back(
                get_temp_meshset_ptr(cP.mField.get_moab()));
            EntityHandle new_meshset = *meshses_tmp_list.back();
            auto meshset = m.getMeshset();
            std::vector<Tag> tmp_tags_list;
            CHKERR cP.mField.get_moab().tag_get_tags_on_entity(meshset,
                                                               tmp_tags_list);
            Range ents;
            CHKERR cP.mField.get_moab().get_entities_by_handle(meshset, ents,
                                                               true);
            ents = intersect(ents, bit);
            CHKERR cP.mField.get_moab().add_entities(new_meshset, ents);
            for (auto t : tmp_tags_list) {
              void *tag_vals[1];
              int tag_size[1];
              CHKERR cP.mField.get_moab().tag_get_by_ptr(
                  t, &meshset, 1, (const void **)tag_vals, tag_size);
              CHKERR cP.mField.get_moab().tag_set_by_ptr(t, &new_meshset, 1,
                                                         tag_vals, tag_size);
            }
 
            std::vector<std::string> remove_tags;
            remove_tags.push_back("CATEGORY");
            remove_tags.push_back("OBB");
            remove_tags.push_back("ROOTSETSURF");
            remove_tags.push_back("__PARALLEL_");
 
            for (auto t : tmp_tags_list) {
 
              std::string tag_name;
              CHKERR cP.mField.get_moab().tag_get_name(t, tag_name);
              bool add = true;
 
              for (auto &p : remove_tags) {
                if (tag_name.compare(0, p.size(), p) == 0) {
                  add = false;
                  break;
                }
              }
 
              if (add)
                tags_list.push_back(t);
            }
          }
          return meshses_tmp_list;
        };
 
        auto tags_list = get_tags_list();
        CHKERR add_field("LAMBDA_ARC_LENGTH", tags_list);
        CHKERR add_fe("ARC_LENGTH", tags_list);
 
        auto meshses_tmp_list = add_meshsets(entities_bit, tags_list);
        for (auto &m_ptr : meshses_tmp_list) {
          EntityHandle m = *m_ptr;
          CHKERR cP.mField.get_moab().add_entities(*meshset_ptr, &m, 1);
        }
 
        {
          Tag gid_tag = cP.mField.get_moab().globalId_tag();
          int negative = -1;
          CHKERR cP.mField.get_moab().tag_clear_data(gid_tag, entities_bit,
                                                     &negative);
          Range verts_ents = entities_bit.subset_by_type(MBVERTEX);
          int gid = 1; // moab indexing from 1
          for (auto v : verts_ents) {
            CHKERR cP.mField.get_moab().tag_set_data(gid_tag, &v, 1, &gid);
            gid++;
          }
        }
 
        std::sort(tags_list.begin(), tags_list.end());
        auto new_end = std::unique(tags_list.begin(), tags_list.end());
        tags_list.resize(std::distance(tags_list.begin(), new_end));
 
        for (auto t : tags_list) {
          std::string tag_name;
          CHKERR cP.mField.get_moab().tag_get_name(t, tag_name);
          MOFEM_LOG("CPSolverWorld", Sev::noisy)
              << "Add tag to restart file: " << tag_name;
        }
 
        EntityHandle save_meshset = *meshset_ptr;
        const std::string file_name =
            "restart_" + boost::lexical_cast<std::string>(step) + ".h5m";
        CHKERR cP.mField.get_moab().write_file(file_name.c_str(), "MOAB", "",
                                               &save_meshset, 1, &tags_list[0],
                                               tags_list.size());
 
        MoFEMFunctionReturn(0);
      };
 
      if (cP.mField.get_comm_rank() == 0)
        CHKERR save_restart_file();
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CPSolvers::solveCutMesh(SmartPetscObj<DM> &dm_crack_propagation,
                                       SmartPetscObj<DM> &dm_elastic,
                                       SmartPetscObj<DM> &dm_eigen_elastic,
                                       SmartPetscObj<DM> &dm_material,
                                       SmartPetscObj<DM> &dm_material_forces,
                                       SmartPetscObj<DM> &dm_surface_projection,
                                       SmartPetscObj<DM> &dm_crack_srf_area,
                                       std::vector<int> &surface_ids,
                                       const bool debug) {
  SolveFunctionBegin;
 
  auto calculate_load_factor = [this]() {
    double a = fabs(cP.maxG1) / pow(cP.getArcCtx()->getFieldData(), 2);
    double load_factor =
        copysign(sqrt(cP.gC / a), cP.getArcCtx()->getFieldData());
    return load_factor;
  };
  cP.getArcCtx()->getFieldData() = calculate_load_factor();
  MOFEM_LOG_C("CPSolverWorld", Sev::inform, "Calculated load factor %g",
              cP.getArcCtx()->getFieldData());
  MOFEM_LOG_C("CPSolverSync", Sev::noisy, "Calculated load factor %g",
              cP.getArcCtx()->getFieldData());
  MOFEM_LOG_SYNCHRONISE(cP.mField.get_comm());
 
  auto save_position_on_mesh_tags_from_coords =
      [this](const std::string tag_name) {
        Tag th;
        MoFEMFunctionBegin;
        rval = cP.mField.get_moab().tag_get_handle(tag_name.c_str(), th);
        if (rval != MB_SUCCESS) {
          double def[] = {0, 0, 0};
          CHKERR cP.mField.get_moab().tag_get_handle(
              tag_name.c_str(), 3, MB_TYPE_DOUBLE, th,
              MB_TAG_CREAT | MB_TAG_SPARSE, def);
        }
        Range verts;
        CHKERR cP.mField.getInterface<BitRefManager>()
            ->getEntitiesByTypeAndRefLevel(cP.mapBitLevel["spatial_domain"],
                                           BitRefLevel().set(), MBVERTEX,
                                           verts);
        std::vector<double> coords(verts.size() * 3);
        CHKERR cP.mField.get_moab().get_coords(verts, &*coords.begin());
        CHKERR cP.mField.get_moab().tag_set_data(th, verts, &*coords.begin());
        MoFEMFunctionReturn(0);
      };
 
  int nb_cuts = 0;
  for (;;) {
 
    CHKERR save_position_on_mesh_tags_from_coords("MESH_NODE_POSITIONS");
    CHKERR solvePropagation(dm_crack_propagation, dm_elastic, dm_material,
                            dm_material_forces, dm_surface_projection,
                            dm_crack_srf_area, surface_ids, nb_cuts, true);
 
    // Destroy DMs
    auto reset_dms = [&]() {
      dm_crack_propagation.reset();
      dm_elastic.reset();
      dm_eigen_elastic.reset();
      dm_material.reset();
      dm_material_forces.reset();
      dm_surface_projection.reset();
      dm_crack_srf_area.reset();
    };
 
    auto reset_fes = [&]() {
      std::vector<std::pair<boost::weak_ptr<FEMethod>, std::string>> vec_fe;
      std::vector<
          std::pair<boost::weak_ptr<NonlinearElasticElement>, std::string>>
          v_elastic_ele_str;
 
      auto push_fes = [&]() {
        v_elastic_ele_str.emplace_back(
            std::make_pair(cP.elasticFe, "elasticFe"));
        v_elastic_ele_str.emplace_back(
            std::make_pair(cP.materialFe, "materialFe"));
        vec_fe.emplace_back(std::make_pair(cP.feLhs, "feLhs"));
        vec_fe.emplace_back(std::make_pair(cP.feRhs, "feRhs"));
        vec_fe.emplace_back(std::make_pair(cP.feMaterialRhs, "feMaterialRhs"));
        vec_fe.emplace_back(std::make_pair(cP.feMaterialLhs, "feMaterialLhs"));
        vec_fe.emplace_back(std::make_pair(cP.feEnergy, "feEnergy"));
        vec_fe.emplace_back(
            std::make_pair(cP.feCouplingElasticLhs, "feCouplingElasticLhs"));
        vec_fe.emplace_back(
            std::make_pair(cP.feCouplingMaterialLhs, "feCouplingMaterialLhs"));
        vec_fe.emplace_back(std::make_pair(cP.feSmootherRhs, "feSmootherRhs"));
        vec_fe.emplace_back(std::make_pair(cP.feSmootherLhs, "feSmootherLhs"));
        vec_fe.emplace_back(
            std::make_pair(cP.feGriffithForceRhs, "feGriffithForceRhs"));
        vec_fe.emplace_back(std::make_pair(cP.feGriffithConstrainsDelta,
                                           "feGriffithConstrainsDelta"));
        vec_fe.emplace_back(std::make_pair(cP.feGriffithConstrainsRhs,
                                           "feGriffithConstrainsRhs"));
        vec_fe.emplace_back(std::make_pair(cP.feGriffithConstrainsLhs,
                                           "feGriffithConstrainsLhs"));
        vec_fe.emplace_back(
            std::make_pair(cP.feSpringLhsPtr, "feSpringLhsPtr"));
        vec_fe.emplace_back(
            std::make_pair(cP.feSpringRhsPtr, "feSpringRhsPtr"));
        vec_fe.emplace_back(
            std::make_pair(cP.feRhsSimpleContact, "feRhsSimpleContact"));
        vec_fe.emplace_back(
            std::make_pair(cP.feLhsSimpleContact, "feLhsSimpleContact"));
        vec_fe.emplace_back(std::make_pair(cP.feMaterialAnaliticalTraction,
                                           "feMaterialAnaliticalTraction"));
 
        vec_fe.emplace_back(
            std::make_pair(cP.bothSidesConstrains, "bothSidesConstrains"));
        vec_fe.emplace_back(
            std::make_pair(cP.closeCrackConstrains, "closeCrackConstrains"));
        vec_fe.emplace_back(std::make_pair(cP.bothSidesContactConstrains,
                                           "bothSidesContactConstrains"));
        vec_fe.emplace_back(
            std::make_pair(cP.fixMaterialEnts, "fixMaterialEnts"));
        vec_fe.emplace_back(
            std::make_pair(cP.feSpringLhsPtr, "feSpringLhsPtr"));
        vec_fe.emplace_back(
            std::make_pair(cP.feSpringRhsPtr, "feSpringRhsPtr"));
        vec_fe.emplace_back(
            std::make_pair(cP.assembleFlambda, "assembleFlambda"));
        vec_fe.emplace_back(std::make_pair(cP.zeroFlambda, "zeroFlambda"));
      };
 
      auto reset = [&]() {
        cP.smootherFe.reset();
        cP.feSmootherRhs.reset();
        cP.feSmootherLhs.reset();
        cP.volumeLengthDouble.reset();
        cP.volumeLengthAdouble.reset();
        cP.griffithForceElement.reset();
 
        cP.tangentConstrains.reset();
        cP.skinOrientation.reset();
        cP.crackOrientation.reset();
        cP.contactOrientation.reset();
        for (auto &m : cP.surfaceConstrain)
          m.second.reset();
        for (auto &m : cP.edgeConstrains)
          m.second.reset();
 
        cP.projSurfaceCtx.reset();
        cP.projFrontCtx.reset();
 
        cP.elasticFe.reset();
        cP.materialFe.reset();
        cP.feLhs.reset();
        cP.feRhs.reset();
        cP.feMaterialRhs.reset();
        cP.feMaterialLhs.reset();
        cP.feEnergy.reset();
        cP.feCouplingElasticLhs.reset();
        cP.feCouplingMaterialLhs.reset();
        cP.feSmootherRhs.reset();
        cP.feSmootherLhs.reset();
        cP.feGriffithForceRhs.reset();
        cP.feGriffithConstrainsDelta.reset();
        cP.feGriffithConstrainsRhs.reset();
        cP.feGriffithConstrainsLhs.reset();
        cP.feSpringLhsPtr.reset();
        cP.feSpringRhsPtr.reset();
        cP.feRhsSimpleContact.reset();
        cP.feLhsSimpleContact.reset();
        cP.feMaterialAnaliticalTraction.reset();
        cP.bothSidesConstrains.reset();
        cP.bothSidesContactConstrains.reset();
        cP.fixMaterialEnts.reset();
        cP.assembleFlambda.reset();
        cP.zeroFlambda.reset();
        cP.closeCrackConstrains.reset();
      };
 
      auto check = [&]() {
        for (auto v : vec_fe)
          if (auto a = v.first.lock()) {
            MOFEM_LOG("CPSolverWorld", Sev::error)
                << "fe " << v.second << " not destroyed " << a.use_count();
          }
 
        for (auto v : v_elastic_ele_str)
          if (auto a = v.first.lock()) {
            MOFEM_LOG("CPSolverWorld", Sev::error)
                << "fe elastic " << v.second << " not destroyed "
                << a.use_count();
          }
 
        if (dm_crack_propagation.use_count())
          MOFEM_LOG("CPSolverWorld", Sev::error)
              << "dm_crack_propagation not destroyed "
              << dm_crack_propagation.use_count();
        if (dm_elastic.use_count())
          MOFEM_LOG("CPSolverWorld", Sev::error)
              << "dm_elastic not destroyed" << dm_elastic.use_count();
        if (dm_material_forces.use_count())
          MOFEM_LOG("CPSolverWorld", Sev::error)
              << "dm_material_forces not destroyed "
              << dm_material_forces.use_count();
        if (dm_surface_projection.use_count())
          MOFEM_LOG("CPSolverWorld", Sev::error)
              << "dm_surface_projection not destroyed "
              << dm_surface_projection.use_count();
        if (dm_crack_srf_area.use_count())
          MOFEM_LOG("CPSolverWorld", Sev::error)
              << "dm_crack_srf_area not destroyed "
              << dm_crack_srf_area.use_count();
 
        if (cP.tangentConstrains.use_count())
          MOFEM_LOG("CPSolverWorld", Sev::error)
              << "tangentConstrains not destroyed";
        if (cP.skinOrientation.use_count())
          MOFEM_LOG("CPSolverWorld", Sev::error)
              << "skinOrientation not destroyed";
        if (cP.crackOrientation.use_count())
          MOFEM_LOG("CPSolverWorld", Sev::error)
              << "crackOrientation not destroyed";
        if (cP.contactOrientation.use_count())
          MOFEM_LOG("CPSolverWorld", Sev::error)
              << "contactOrientation not destroyed";
        for (auto &m : cP.surfaceConstrain)
          if (m.second.use_count())
            MOFEM_LOG("CPSolverWorld", Sev::error)
                << "surfaceConstrain not destroyed : " << m.first;
        for (auto &m : cP.edgeConstrains)
          if (m.second.use_count())
            MOFEM_LOG("CPSolverWorld", Sev::error)
                << "edgeConstrains not destroyed : " << m.first;
        cP.surfaceConstrain.clear();
        cP.edgeConstrains.clear();
 
        if (cP.projSurfaceCtx.use_count())
          MOFEM_LOG("CPSolverWorld", Sev::error)
              << "projSurfaceCtx not destroyed";
        if (cP.projFrontCtx.use_count())
          MOFEM_LOG("CPSolverWorld", Sev::error)
              << "projFrontCtx not destroyed";
      };
 
      push_fes();
      reset();
      check();
    };
 
    reset_dms();
    reset_fes();
 
    if (cP.mwlsApprox) {
      MOFEM_LOG("MWLSWorld", Sev::inform)
          << "Resest MWLS approximation coefficients. Coefficients will be "
             "recalulated for current material positions.";
      cP.mwlsApprox->invABMap.clear();
      cP.mwlsApprox->influenceNodesMap.clear();
      cP.mwlsApprox->dmNodesMap.clear();
    }
 
    ++nb_cuts;
    if (nb_cuts <= cP.getNbCutSteps()) {
 
      // clear database
      CHKERR cP.mField.clear_database(QUIET);
      CHKERR cP.getInterface<CPMeshCut>()->clearData();
 
      auto tag_gid = cP.mField.get_moab().globalId_tag();
      std::vector<Tag> tag_handles;
      CHKERR cP.mField.get_moab().tag_get_tags(tag_handles);
      for (auto th : tag_handles)
        if (th != tag_gid)
          CHKERR cP.mField.get_moab().tag_delete(th);
      CHKERR cP.mField.get_moab().delete_mesh();
 
      // reaload mesh
      CHKERR PetscBarrier(PETSC_NULL);
      const std::string file_name =
          "restart_" + boost::lexical_cast<std::string>(cP.startStep - 1) +
          ".h5m";
 
      auto add_bits_tets = [&]() {
        MoFEMFunctionBegin;
        CHKERR cP.mField.getInterface<BitRefManager>()
            ->addToDatabaseBitRefLevelByType(MBTET, BitRefLevel().set(),
                                             BitRefLevel().set());
        MoFEMFunctionReturn(0);
      };
 
      std::vector<int> surface_ids;
      surface_ids.push_back(cP.getInterface<CPMeshCut>()->getSkinOfTheBodyId());
      std::vector<std::string> fe_surf_proj_list;
      fe_surf_proj_list.push_back("SURFACE");
 
      ParallelComm *pcomm =
          ParallelComm::get_pcomm(&cP.mField.get_moab(), MYPCOMM_INDEX);
      if (pcomm == NULL)
        pcomm = new ParallelComm(&cP.mField.get_moab(),
                                 cP.moabCommWorld->get_comm());
 
      if (cP.mField.get_comm_size() == 1 ||
          !CrackPropagation::parallelReadAndBroadcast) {
        const char *option = "";
        CHKERR cP.mField.get_moab().load_file(file_name.c_str(), 0, option);
      } else {
        if (cP.mField.get_comm_rank() == 0) {
          // Read mesh file
          moab::Core mb_instance_read;
          moab::Interface &moab_read = mb_instance_read;
          ParallelComm *pcomm_read =
              ParallelComm::get_pcomm(&moab_read, MYPCOMM_INDEX);
          if (pcomm_read == NULL)
            pcomm_read =
                new ParallelComm(&moab_read, cP.moabCommWorld->get_comm());
 
          const char *option = "";
          CHKERR moab_read.load_file(file_name.c_str(), 0, option);
          Range ents;
          CHKERR moab_read.get_entities_by_handle(0, ents, false);
          CHKERR moab_read.get_entities_by_type(0, MBENTITYSET, ents, false);
          CHKERR FractureMechanics::broadcast_entities(
              cP.mField.get_moab(), moab_read, 0, ents, false, true);
        } else {
          Range ents;
          CHKERR FractureMechanics::broadcast_entities(
              cP.mField.get_moab(), cP.mField.get_moab(), 0, ents, false, true);
        }
      }
 
      CHKERR MoFEM::Interface::setFileVersion(cP.mField.get_moab(),
                                              cP.fileVersion);
 
      CHKERR cP.mField.set_moab_interface(cP.mField.get_moab(), QUIET);
      CHKERR cP.getInterface<CPMeshCut>()->getInterfacesPtr();
      CHKERR add_bits_tets();
 
      Range ents;
      CHKERR cP.mField.getInterface<BitRefManager>()
          ->getAllEntitiesNotInDatabase(ents);
      Range meshsets;
      CHKERR cP.mField.get_moab().get_entities_by_type(0, MBENTITYSET, meshsets,
                                                       false);
      for (auto m : meshsets)
        CHKERR cP.mField.get_moab().remove_entities(m, ents);
      CHKERR cP.mField.get_moab().delete_entities(ents);
 
      string cutting_surface_name =
          "cutting_surface_" + boost::lexical_cast<std::string>(cP.startStep) +
          ".vtk";
 
      if (cP.mField.get_comm_rank() == 0)
        CHKERR cP.getInterface<CPMeshCut>()->rebuildCrackSurface(
            cP.crackAccelerationFactor, cutting_surface_name, QUIET, debug);
      else
        CHKERR cP.getInterface<CPMeshCut>()->rebuildCrackSurface(
            cP.crackAccelerationFactor, "", QUIET, debug);
 
      CHKERR cP.getInterface<CPMeshCut>()->refineMesh(nullptr, false, QUIET,
                                                      debug);
      CHKERR cP.getInterface<CPMeshCut>()->cutRefineAndSplit(QUIET, debug);
 
      // Create crack propagation fields and elements
      CHKERR cP.mField.build_field("LAMBDA_ARC_LENGTH");
      CHKERR cP.mField.build_finite_elements("ARC_LENGTH");
 
      CHKERR cP.createProblemDataStructures(surface_ids, QUIET, debug);
 
      // set inital coordinates
      CHKERR cP.setMaterialPositionFromCoords();
      CHKERR cP.setSpatialPositionFromCoords();
      if (cP.addSingularity == PETSC_TRUE)
        CHKERR cP.setSingularDofs("SPATIAL_POSITION");
 
      CHKERR cP.createDMs(dm_elastic, dm_eigen_elastic, dm_material,
                          dm_crack_propagation, dm_material_forces,
                          dm_surface_projection, dm_crack_srf_area, surface_ids,
                          fe_surf_proj_list);
 
      auto set_up_arc_length = [&]() {
        MoFEMFunctionBegin;
        const MoFEM::Problem *problem_ptr;
        CHKERR DMMoFEMGetProblemPtr(dm_elastic, &problem_ptr);
        cP.arcCtx = boost::shared_ptr<ArcLengthCtx>(new ArcLengthCtx(
            cP.mField, problem_ptr->getName(), "LAMBDA_ARC_LENGTH"));
        CHKERR cP.arcCtx->setAlphaBeta(0, 1);
        auto arc_snes_ctx =
            boost::make_shared<CrackPropagation::ArcLengthSnesCtx>(
                cP.mField, problem_ptr->getName(), cP.arcCtx);
        CHKERR DMMoFEMSetSnesCtx(dm_elastic, arc_snes_ctx);
        MoFEMFunctionReturn(0);
      };
 
      CHKERR set_up_arc_length();
 
      CHKERR cP.setSpatialPositionFromCoords();
      if (cP.addSingularity == PETSC_TRUE)
        CHKERR cP.setSingularDofs("SPATIAL_POSITION");
 
      const auto load_factor = std::abs(cP.getArcCtx()->getFieldData());
      MOFEM_LOG("CPSolverSync", Sev::noisy) << "Lambda factor " << load_factor;
      MOFEM_LOG_SYNCHRONISE(cP.mField.get_comm());
      if (cP.solveEigenStressProblem) {
        MOFEM_LOG("CPSolverWorld", Sev::verbose)
            << "Solve Eigen ELASTIC problem";
        CHKERR cP.getInterface<CPSolvers>()->solveElastic(
            dm_eigen_elastic, cP.getEigenArcCtx()->x0, 1);
      }
      cP.getLoadScale() = load_factor;
      MOFEM_LOG("CPSolverSync", Sev::noisy)
          << "Reset lambda factor " << cP.getLoadScale();
      MOFEM_LOG("CPSolverWorld", Sev::noisy)
          << "Solve ELASTIC problem and calculate g";
      CHKERR cP.getInterface<CPSolvers>()->solveElastic(
          dm_elastic, cP.getArcCtx()->x0, load_factor);
 
      // set finite element instances and user data operators on instances
      CHKERR cP.assembleElasticDM(VERBOSE, false);
      CHKERR cP.assembleMaterialForcesDM(PETSC_NULL);
      CHKERR cP.assembleSmootherForcesDM(PETSC_NULL, surface_ids, VERBOSE,
                                         false);
      CHKERR cP.assembleCouplingForcesDM(PETSC_NULL, VERBOSE, false);
 
      CHKERR cP.calculateElasticEnergy(dm_elastic);
      CHKERR cP.getInterface<CPSolvers>()->calculateGc(
          dm_material_forces, dm_surface_projection, dm_crack_srf_area,
          surface_ids);
      cP.getArcCtx()->getFieldData() = calculate_load_factor();
      MOFEM_LOG_C("CPSolverWorld", Sev::inform, "Calculated load factor %g",
                  cP.getArcCtx()->getFieldData());
 
    } else
      break;
  }
 
  MoFEMFunctionReturn(0);
}
} // namespace FractureMechanics