CrackPropagation.cpp

Go to the documentation of this file.

/** \file CrackPropagation.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/>. */
 
 
#define SINGULARITY
#ifdef WITH_TETGEN
 
#include <tetgen.h>
#ifdef REAL
#undef REAL
#endif
 
#endif
 
#include <moab/SpatialLocator.hpp>
#include <moab/ElemEvaluator.hpp>
 
#include <MoFEM.hpp>
using namespace MoFEM;
 
#include <cholesky.hpp>
 
#include <BasicFiniteElements.hpp>
#include <Mortar.hpp>
#include <Hooke.hpp>
#include <AnalyticalFun.hpp>
 
// TODO implementation in ComplexConstArea is obsolete, not follows naming
// convention and code is unnecessarily complicated and in places difficult to
// follow.
#include <ComplexConstArea.hpp>
#include <ConstantArea.hpp>
 
extern "C" {
#include <phg-quadrule/quad.h>
}
 
#include <NeoHookean.hpp>
#include <Hooke.hpp>
#include <Smoother.hpp>
#include <VolumeLengthQuality.hpp>
#include <SurfaceSlidingConstrains.hpp>
#include <Mortar.hpp>
 
#include <CrackFrontElement.hpp>
#include <MWLS.hpp>
#include <GriffithForceElement.hpp>
#include <CrackPropagation.hpp>
#include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/
 
#include <CPSolvers.hpp>
#include <CPMeshCut.hpp>
 
#include <thread>
 
static MoFEMErrorCode snes_monitor_fields(SNES snes, PetscInt its,
                                          PetscReal fgnorm,
                                          PetscViewerAndFormat *vf) {
  PetscViewer viewer = vf->viewer;
  Vec res;
  DM dm;
  PetscSection s;
  const PetscScalar *r;
  PetscReal *lnorms, *norms;
  PetscInt numFields, f, pStart, pEnd, p;
 
  MoFEMFunctionBegin;
  LogManager::setLog("PETSC");
  MOFEM_LOG_TAG("PETSC", "SnesProp");
 
  PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 4);
  CHKERR SNESGetFunction(snes, &res, 0, 0);
  CHKERR SNESGetDM(snes, &dm);
  CHKERR DMGetDefaultSection(dm, &s);
  CHKERR PetscSectionGetNumFields(s, &numFields);
  CHKERR PetscSectionGetChart(s, &pStart, &pEnd);
  CHKERR PetscCalloc2(numFields, &lnorms, numFields, &norms);
  CHKERR VecGetArrayRead(res, &r);
  for (p = pStart; p < pEnd; ++p) {
    for (f = 0; f < numFields; ++f) {
      PetscInt fdof, foff, d;
 
      CHKERR PetscSectionGetFieldDof(s, p, f, &fdof);
      CHKERR PetscSectionGetFieldOffset(s, p, f, &foff);
      for (d = 0; d < fdof; ++d)
        lnorms[f] += PetscRealPart(PetscSqr(r[foff + d]));
    }
  }
  CHKERR VecRestoreArrayRead(res, &r);
  CHKERR MPIU_Allreduce(lnorms, norms, numFields, MPIU_REAL, MPIU_SUM,
                        PetscObjectComm((PetscObject)dm));
  CHKERR PetscViewerPushFormat(viewer, vf->format);
  CHKERR PetscViewerASCIIAddTab(viewer, ((PetscObject)snes)->tablevel);
  CHKERR PetscViewerASCIIPrintf(viewer, "%3D SNES Function norm %14.12e\n", its,
                                (double)fgnorm);
  for (f = 0; f < numFields; ++f) {
    const char *field_name;
    CHKERR PetscSectionGetFieldName(s, f, &field_name);
    CHKERR PetscViewerASCIIPrintf(viewer, "\t%30.30s\t%14.12e\n", field_name,
                                  (double)PetscSqrtReal(norms[f]));
  }
  CHKERR PetscViewerASCIISubtractTab(viewer, ((PetscObject)snes)->tablevel);
  CHKERR PetscViewerPopFormat(viewer);
  CHKERR PetscFree2(lnorms, norms);
 
  LogManager::setLog("PETSC");
  MoFEMFunctionReturn(0);
}
 
namespace po = boost::program_options;
 
namespace FractureMechanics {
 
MoFEMErrorCode clean_pcomms(moab::Interface &moab,
                            boost::shared_ptr<WrapMPIComm> moab_comm_wrap) {
  MoFEMFunctionBegin;
  std::vector<ParallelComm *> list_pcomms;
  ParallelComm::get_all_pcomm(&moab, list_pcomms);
  if (list_pcomms[MYPCOMM_INDEX]) {
    CHKERR moab.tag_delete(list_pcomms[MYPCOMM_INDEX]->pstatus_tag());
    CHKERR moab.tag_delete(list_pcomms[MYPCOMM_INDEX]->sharedps_tag());
    CHKERR moab.tag_delete(list_pcomms[MYPCOMM_INDEX]->sharedhs_tag());
    CHKERR moab.tag_delete(list_pcomms[MYPCOMM_INDEX]->sharedp_tag());
    CHKERR moab.tag_delete(list_pcomms[MYPCOMM_INDEX]->sharedh_tag());
  }
  for (auto p : list_pcomms) {
    delete p;
  }
  int pcomm_id = MYPCOMM_INDEX;
  list_pcomms[MYPCOMM_INDEX] =
      new ParallelComm(&moab, moab_comm_wrap->get_comm(), &pcomm_id);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode broadcast_entities(moab::Interface &moab,
                                  moab::Interface &moab_tmp,
                                  const int from_proc, Range &entities,
                                  const bool adjacencies, const bool tags) {
  ErrorCode result = MB_SUCCESS;
  int success;
  int buff_size;
 
  auto add_verts = [&](Range &sent_ents) {
    // Get the verts adj to these entities, since we'll have to send those
    // too
 
    // First check sets
    std::pair<Range::const_iterator, Range::const_iterator> set_range =
        sent_ents.equal_range(MBENTITYSET);
    ErrorCode result = MB_SUCCESS, tmp_result;
    for (Range::const_iterator rit = set_range.first; rit != set_range.second;
         ++rit) {
      tmp_result = moab_tmp.get_entities_by_type(*rit, MBVERTEX, sent_ents);
      CHK_MOAB_THROW(tmp_result, "Failed to get contained verts");
    }
 
    // Now non-sets
    Range tmp_ents;
    std::copy(sent_ents.begin(), set_range.first, range_inserter(tmp_ents));
    result = moab_tmp.get_adjacencies(tmp_ents, 0, false, sent_ents,
                                      moab::Interface::UNION);
    CHK_MOAB_THROW(result, "Failed to get vertices adj to ghosted ents");
 
    return result;
  };
 
  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab_tmp, MYPCOMM_INDEX);
  auto &procConfig = pcomm->proc_config();
  const int MAX_BCAST_SIZE = (1 << 28);
 
  ParallelComm::Buffer buff(ParallelComm::INITIAL_BUFF_SIZE);
  buff.reset_ptr(sizeof(int));
  if ((int)procConfig.proc_rank() == from_proc) {
 
    result = add_verts(entities);
    CHK_MOAB_THROW(result, "Failed to add adj vertices");
 
    buff.reset_ptr(sizeof(int));
    result = pcomm->pack_buffer(entities, adjacencies, tags, false, -1, &buff);
    CHK_MOAB_THROW(result,
                   "Failed to compute buffer size in broadcast_entities");
    buff.set_stored_size();
    buff_size = buff.buff_ptr - buff.mem_ptr;
  }
 
  success =
      MPI_Bcast(&buff_size, 1, MPI_INT, from_proc, procConfig.proc_comm());
  if (MPI_SUCCESS != success) {
    THROW_MESSAGE("MPI_Bcast of buffer size failed");
  }
 
  if (!buff_size) // No data
    return MB_SUCCESS;
 
  if ((int)procConfig.proc_rank() != from_proc)
    buff.reserve(buff_size);
 
  size_t offset = 0;
  while (buff_size) {
    int sz = std::min(buff_size, MAX_BCAST_SIZE);
    success = MPI_Bcast(buff.mem_ptr + offset, sz, MPI_UNSIGNED_CHAR, from_proc,
                        procConfig.proc_comm());
    if (MPI_SUCCESS != success) {
      THROW_MESSAGE("MPI_Bcast of buffer failed");
    }
 
    offset += sz;
    buff_size -= sz;
  }
 
  std::vector<std::vector<EntityHandle>> dum1a, dum1b;
  std::vector<std::vector<int>> dum1p;
  std::vector<EntityHandle> dum2, dum4;
  std::vector<unsigned int> dum3;
  buff.reset_ptr(sizeof(int));
  if ((int)procConfig.proc_rank() == from_proc) {
    ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
    result = pcomm->unpack_buffer(buff.buff_ptr, false, from_proc, -1, dum1a,
                                  dum1b, dum1p, dum2, dum2, dum3, dum4);
  } else {
    result = pcomm->unpack_buffer(buff.buff_ptr, false, from_proc, -1, dum1a,
                                  dum1b, dum1p, dum2, dum2, dum3, dum4);
  }
  CHK_MOAB_THROW(result, "Failed to unpack buffer in broadcast_entities");
  std::copy(dum4.begin(), dum4.end(), range_inserter(entities));
 
  return MB_SUCCESS;
};
 
struct GetSmoothingElementsSkin {
  CrackPropagation &cP;
  double gcFixThreshold;
  GetSmoothingElementsSkin(CrackPropagation &cp) : cP(cp), gcFixThreshold(0.5) {
    ierr = getOptions();
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
  }
 
  MoFEMErrorCode getOptions() {
    MoFEMFunctionBegin;
    ierr =
        PetscOptionsBegin(PETSC_COMM_WORLD, "", "Fixing nodes options", "none");
    CHKERRG(ierr);
    CHKERR PetscOptionsScalar("-gc_fix_threshold",
                              "fix nodes below given threshold", "",
                              gcFixThreshold, &gcFixThreshold, PETSC_NULL);
    CHKERR PetscPrintf(PETSC_COMM_WORLD,
                       "### Input parameter: -gc_fix_threshold %6.4e \n",
                       gcFixThreshold);
 
    ierr = PetscOptionsEnd();
    CHKERRG(ierr);
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode operator()(Range &fix_material_nodes,
                            const bool fix_small_g = true,
                            const bool debug = false) const {
    MoFEM::Interface &m_field = cP.mField;
    MoFEMFunctionBegin;
 
    Skinner skin(&m_field.get_moab());
    Range tets_level, prisms_level;
    CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
        cP.mapBitLevel["material_domain"], BitRefLevel().set(), MBTET,
        tets_level);
    CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
        cP.mapBitLevel["material_domain"], BitRefLevel().set(), MBPRISM,
        prisms_level);
 
    prisms_level = intersect(prisms_level, cP.contactElements);
    tets_level.merge(prisms_level);
 
    Range mesh_skin;
    CHKERR skin.find_skin(0, tets_level, false, mesh_skin);
 
    mesh_skin = mesh_skin.subset_by_type(MBTRI);
 
    Range faces_to_remove, contact_faces;
    faces_to_remove = unite(cP.crackFaces, cP.bodySkin);
    contact_faces = unite(cP.contactMasterFaces, cP.contactSlaveFaces);
    faces_to_remove = subtract(faces_to_remove, contact_faces);
 
    mesh_skin = subtract(mesh_skin, faces_to_remove);
 
    CHKERR m_field.get_moab().get_connectivity(mesh_skin, fix_material_nodes);
    // fix fixed edges nodes
    Range fixed_edges = cP.getInterface<CPMeshCut>()->getFixedEdges();
    Range fixed_edges_nodes;
    CHKERR m_field.get_moab().get_connectivity(fixed_edges, fixed_edges_nodes,
                                               true);
 
    Range crack_faces_and_fixed_edges_nodes;
    CHKERR m_field.get_moab().get_connectivity(
        cP.crackFaces, crack_faces_and_fixed_edges_nodes, true);
    crack_faces_and_fixed_edges_nodes =
        intersect(crack_faces_and_fixed_edges_nodes, fixed_edges_nodes);
 
    Range crack_faces_edges;
    CHKERR m_field.get_moab().get_adjacencies(
        cP.crackFaces, 1, false, crack_faces_edges, moab::Interface::UNION);
    Range not_fixed_edges = intersect(crack_faces_edges, fixed_edges);
    Range not_fixed_edges_nodes;
    CHKERR m_field.get_moab().get_connectivity(not_fixed_edges,
                                               not_fixed_edges_nodes, true);
    crack_faces_and_fixed_edges_nodes =
        subtract(crack_faces_and_fixed_edges_nodes, not_fixed_edges_nodes);
    fix_material_nodes.merge(crack_faces_and_fixed_edges_nodes);
 
    Range not_fixed_edges_skin;
    CHKERR skin.find_skin(0, not_fixed_edges, false, not_fixed_edges_skin);
    fix_material_nodes.merge(
        subtract(not_fixed_edges_skin, cP.crackFrontNodes));
    if (debug && !m_field.get_comm_rank()) {
      EntityHandle meshset;
      CHKERR m_field.get_moab().create_meshset(MESHSET_SET, meshset);
      CHKERR m_field.get_moab().add_entities(meshset, mesh_skin);
      CHKERR m_field.get_moab().write_file("fixed_faces.vtk", "VTK", "",
                                           &meshset, 1);
      CHKERR m_field.get_moab().delete_entities(&meshset, 1);
    }
 
    // fix corner nodes
    fix_material_nodes.merge(cP.getInterface<CPMeshCut>()->getCornerNodes());
 
    // fix nodes on kinematic and static boundary conditions
    auto bc_fix_nodes = [&m_field,
                         &fix_material_nodes](const EntityHandle meshset) {
      MoFEMFunctionBegin;
      Range ents;
      CHKERR m_field.get_moab().get_entities_by_handle(meshset, ents, true);
      Range nodes;
      CHKERR m_field.get_moab().get_connectivity(ents, nodes);
      fix_material_nodes.merge(nodes);
      MoFEMFunctionReturn(0);
    };
    for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
             m_field, NODESET | DISPLACEMENTSET, it)) {
      CHKERR bc_fix_nodes(it->getMeshset());
    }
 
    if (!cP.isSurfaceForceAle) {
      for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field,
                                                      NODESET | FORCESET, it)) {
        CHKERR bc_fix_nodes(it->getMeshset());
      }
    }
 
    if (!cP.isPressureAle) {
      // fix nodes on surfaces under pressure
      for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
               m_field, SIDESET | PRESSURESET, it)) {
        CHKERR bc_fix_nodes(it->getMeshset());
      }
    }
 
    // fix nodes on contact surfaces
    if (cP.fixContactNodes) {
      Range contact_nodes;
      CHKERR m_field.get_moab().get_connectivity(cP.contactElements,
                                                 contact_nodes);
      fix_material_nodes.merge(contact_nodes);
    }
 
    if (!cP.areSpringsAle) {
      // fix nodes on spring nodes
      for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, bit)) {
        if (bit->getName().compare(0, 9, "SPRING_BC") == 0) {
          CHKERR bc_fix_nodes(bit->meshset);
        }
      }
    }
 
    {
      Range contact_nodes;
      CHKERR m_field.get_moab().get_connectivity(cP.mortarContactElements,
                                                 contact_nodes);
      fix_material_nodes.merge(contact_nodes);
    }
 
    // fix nodes on the constrained interface
    if (!cP.constrainedInterface.empty()) {
      Range interface_nodes;
      CHKERR m_field.get_moab().get_connectivity(cP.constrainedInterface,
                                                 interface_nodes);
      fix_material_nodes.merge(interface_nodes);
    }
 
    auto add_nodes_on_opposite_side_of_prism = [&m_field, &fix_material_nodes,
                                                this](bool only_contact_prisms =
                                                          false) {
      MoFEMFunctionBegin;
      Range adj_prisms;
      CHKERR m_field.get_moab().get_adjacencies(
          fix_material_nodes, 3, false, adj_prisms, moab::Interface::UNION);
      adj_prisms = adj_prisms.subset_by_type(MBPRISM);
      if (only_contact_prisms) {
        adj_prisms = intersect(adj_prisms, cP.contactElements);
      }
      for (auto p : adj_prisms) {
        int num_nodes;
        const EntityHandle *conn;
        CHKERR m_field.get_moab().get_connectivity(p, conn, num_nodes, false);
        Range add_nodes;
        for (int n = 0; n != 3; ++n) {
          if (fix_material_nodes.find(conn[n]) != fix_material_nodes.end()) {
            add_nodes.insert(conn[3 + n]);
          }
          if (fix_material_nodes.find(conn[3 + n]) !=
              fix_material_nodes.end()) {
            add_nodes.insert(conn[n]);
          }
        }
        fix_material_nodes.merge(add_nodes);
      }
      MoFEMFunctionReturn(0);
    };
 
    CHKERR add_nodes_on_opposite_side_of_prism();
 
    cP.fractionOfFixedNodes = 0;
    if (fix_small_g) {
 
      if (!cP.mField.get_comm_rank()) {
        auto copy_map = [&fix_material_nodes](const auto &map) {
          std::map<EntityHandle, double> map_copy;
          for (auto &m : map) {
            if (fix_material_nodes.find(m.first) == fix_material_nodes.end()) {
              map_copy[m.first] = m.second;
            }
          }
          return map_copy;
        };
 
        const std::map<EntityHandle, double> map_j = copy_map(cP.mapJ);
        const std::map<EntityHandle, double> map_g1 = copy_map(cP.mapG1);
 
        auto find_max = [](const auto &map) {
          double max = 0;
          for (auto &m : map)
            max = fmax(m.second, max);
          return max;
        };
 
        const double max_j = find_max(map_j);
        const double max_g1 = find_max(map_g1);
 
        for (auto &m : map_j) {
          if ((m.second < gcFixThreshold * max_j) &&
              (map_g1.at(m.first) < gcFixThreshold * max_g1)) {
            fix_material_nodes.insert(m.first);
          } else {
            cP.fractionOfFixedNodes += 1;
          }
        }
 
        if (map_j.empty() || map_g1.empty() || cP.fractionOfFixedNodes == 0) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "No points to move on crack front");
        }
 
        cP.fractionOfFixedNodes /= map_j.size();
      }
 
      Vec vec;
      CHKERR VecCreateMPI(cP.mField.get_comm(), PETSC_DETERMINE, 1, &vec);
      if (!cP.mField.get_comm_rank()) {
        CHKERR VecSetValue(vec, cP.mField.get_comm_rank(),
                           cP.fractionOfFixedNodes, INSERT_VALUES);
      }
      CHKERR VecMax(vec, PETSC_NULL, &cP.fractionOfFixedNodes);
      CHKERR VecDestroy(&vec);
    }
 
    CHKERR add_nodes_on_opposite_side_of_prism(true);
 
    MoFEMFunctionReturn(0);
  }
};
 
const char *materials_list[] = {"HOOKE", "KIRCHHOFF", "NEOHOOKEAN",
                                "BONEHOOKE"};
 
MoFEMErrorCode
CrackPropagation::query_interface(boost::typeindex::type_index type_index,
                                  UnknownInterface **iface) const {
  MoFEMFunctionBeginHot;
  *iface = NULL;
 
  if (type_index == boost::typeindex::type_id<CrackPropagation>()) {
    *iface = const_cast<CrackPropagation *>(this);
    return 0;
  }
 
  if (type_index == boost::typeindex::type_id<CPSolvers>()) {
    CHKERR cpSolversPtr->query_interface(type_index, iface);
    return 0;
  }
 
  if (type_index == boost::typeindex::type_id<CPMeshCut>()) {
    CHKERR cpMeshCutPtr->query_interface(type_index, iface);
    return 0;
  }
  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "unknown interface");
  MoFEMFunctionReturnHot(0);
}
 
CrackPropagation::CrackPropagation(MoFEM::Interface &m_field,
                                   const int approx_order,
                                   const int geometry_order)
    : mField(m_field), contactPostProcMoab(contactPostProcCore),
      setSingularCoordinates(false), approxOrder(approx_order),
      geometryOrder(geometry_order), gC(0), betaGc(0),
      isGcHeterogeneous(PETSC_FALSE), isPartitioned(PETSC_FALSE),
      propagateCrack(PETSC_FALSE), refAtCrackTip(0), refOrderAtTip(0),
      residualStressBlock(-1), densityMapBlock(-1),
      addAnalyticalInternalStressOperators(PETSC_FALSE), postProcLevel(0),
      startStep(0), crackFrontLength(-1), crackSurfaceArea(-1), rHo0(1),
      nBone(0), cpSolversPtr(new CPSolvers(*this)),
      cpMeshCutPtr(new CPMeshCut(*this)) {
 
  if (!LogManager::checkIfChannelExist("CPWorld")) {
    auto core_log = logging::core::get();
 
    core_log->add_sink(
        LogManager::createSink(LogManager::getStrmWorld(), "CPWorld"));
    core_log->add_sink(
        LogManager::createSink(LogManager::getStrmSync(), "CPSync"));
    core_log->add_sink(
        LogManager::createSink(LogManager::getStrmSelf(), "CPSelf"));
 
    LogManager::setLog("CPWorld");
    LogManager::setLog("CPSync");
    LogManager::setLog("CPSelf");
 
    MOFEM_LOG_TAG("CPWorld", "CP");
    MOFEM_LOG_TAG("CPSync", "CP");
    MOFEM_LOG_TAG("CPSelf", "CP");
  }
 
  MOFEM_LOG("CPWorld", Sev::noisy) << "CPSolve created";
 
#ifdef GIT_UM_SHA1_NAME
  MOFEM_LOG_C("CPWorld", Sev::inform, "User module git commit id %s",
              GIT_UM_SHA1_NAME);
#endif
 
  Version version;
  getInterfaceVersion(version);
  MOFEM_LOG("CPWorld", Sev::inform)
      << "Fracture module version " << version.strVersion();
 
#ifdef GIT_FM_SHA1_NAME
  MOFEM_LOG_C("CPWorld", Sev::inform, "Fracture module git commit id %s",
              GIT_FM_SHA1_NAME);
#endif
 
  ierr = registerInterface<CrackPropagation>();
  CHKERRABORT(PETSC_COMM_SELF, ierr);
  ierr = registerInterface<CPSolvers>();
  CHKERRABORT(PETSC_COMM_SELF, ierr);
  ierr = registerInterface<CPMeshCut>();
  CHKERRABORT(PETSC_COMM_SELF, ierr);
 
  mapBitLevel["mesh_cut"] = BitRefLevel().set(0);
  mapBitLevel["spatial_domain"] = BitRefLevel().set(1);
  mapBitLevel["material_domain"] = BitRefLevel().set(2);
 
  if (!moabCommWorld)
    moabCommWorld = boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD, false);
}
 
CrackPropagation::~CrackPropagation() {}
 
MoFEMErrorCode CrackPropagation::getOptions() {
  MoFEMFunctionBegin;
  ierr = PetscOptionsBegin(mField.get_comm(), "", "Fracture options", "none");
  CHKERRQ(ierr);
  {
    CHKERR PetscOptionsBool("-my_propagate_crack",
                            "true if crack is propagated", "", propagateCrack,
                            &propagateCrack, NULL);
    CHKERR PetscOptionsInt("-my_order", "approximation order", "", approxOrder,
                           &approxOrder, PETSC_NULL);
    CHKERR PetscOptionsInt("-my_geom_order", "approximation geometry order", "",
                           geometryOrder, &geometryOrder, PETSC_NULL);
    CHKERR PetscOptionsScalar("-my_gc", "release energy", "", gC, &gC,
                              PETSC_NULL);
    CHKERR PetscOptionsScalar("-beta_gc", "heterogeneous gc beta coefficient",
                              "", betaGc, &betaGc, &isGcHeterogeneous);
    CHKERR PetscOptionsBool("-my_is_partitioned", "true if mesh is partitioned",
                            "", isPartitioned, &isPartitioned, NULL);
    CHKERR PetscOptionsInt("-my_ref", "crack tip mesh refinement level", "",
                           refAtCrackTip, &refAtCrackTip, PETSC_NULL);
    CHKERR PetscOptionsInt("-my_ref_order",
                           "crack tip refinement approximation order level", "",
                           refOrderAtTip, &refOrderAtTip, PETSC_NULL);
 
    MOFEM_LOG_C("CPWorld", Sev::inform, "### Input parameter: -my_order %d",
                approxOrder);
    MOFEM_LOG_C("CPWorld", Sev::inform,
                "### Input parameter: -my_geom_order %d", geometryOrder);
    MOFEM_LOG_C("CPWorld", Sev::inform, "### Input parameter: -my_gc %6.4e",
                gC);
    MOFEM_LOG_C("CPWorld", Sev::inform,
                "### Input parameter: -my_propagate_crack %d", propagateCrack);
    MOFEM_LOG_C("CPWorld", Sev::inform, "### Input parameter: -my_ref %d",
                refAtCrackTip);
    MOFEM_LOG_C("CPWorld", Sev::inform, "### Input parameter: -my_ref_order %d",
                refOrderAtTip);
 
    onlyHookeFromOptions = PETSC_TRUE;
    CHKERR PetscOptionsBool(
        "-my_hook_elastic",
        "if false force use of nonlinear element for hooke material", "",
        onlyHookeFromOptions, &onlyHookeFromOptions, NULL);
 
    isPressureAle = PETSC_TRUE;
    CHKERR PetscOptionsBool("-my_add_pressure_ale",
                            "if set surface pressure is considered in ALE", "",
                            isPressureAle, &isPressureAle, NULL);
 
    areSpringsAle = PETSC_TRUE;
    CHKERR PetscOptionsBool("-my_add_springs_ale",
                            "if set surface springs is considered in ALE", "",
                            areSpringsAle, &areSpringsAle, NULL);
 
    isSurfaceForceAle = PETSC_TRUE;
    CHKERR PetscOptionsBool("-my_add_surface_force_ale",
                            "if set surface force is considered in ALE", "",
                            isSurfaceForceAle, &isSurfaceForceAle, NULL);
 
    ignoreMaterialForce = PETSC_FALSE;
    CHKERR PetscOptionsBool(
        "-my_ignore_material_surface_force",
        "if set material forces arising from surface force are ignored", "",
        ignoreMaterialForce, &ignoreMaterialForce, NULL);
 
    addSingularity = PETSC_TRUE;
    CHKERR PetscOptionsBool("-my_add_singularity",
                            "if set singularity added to crack front", "",
                            addSingularity, &addSingularity, NULL);
    MOFEM_LOG_C("CPWorld", Sev::inform,
                "### Input parameter: -my_add_singularity %d", addSingularity);
    {
      PetscBool flg;
      char file_name[255] = "mwls.med";
      CHKERR PetscOptionsString(
          "-my_mwls_approx_file",
          "file with data from med file (code-aster) with radiation data", "",
          file_name, file_name, 255, &flg);
      if (flg == PETSC_TRUE) {
        mwlsApproxFile = std::string(file_name);
        MOFEM_LOG_C("CPWorld", Sev::inform,
                    "### Input parameter: -my_mwls_approx_file %s", file_name);
      }
      char tag_name[255] = "SIGP";
      CHKERR PetscOptionsString("-my_internal_stress_name",
                                "name of internal stress tag", "", tag_name,
                                tag_name, 255, &flg);
      MOFEM_LOG_C("CPWorld", Sev::inform,
                  "### Input parameter: -my_internal_stress_name %s", tag_name);
      mwlsStressTagName = "MED_" + std::string(tag_name);
      CHKERR PetscOptionsString("-my_eigen_stress_name",
                                "name of eigen stress tag", "", tag_name,
                                tag_name, 255, &flg);
      MOFEM_LOG_C("CPWorld", Sev::inform,
                  "### Input parameter: -my_eigen_stress_name %s", tag_name);
      mwlsEigenStressTagName = "MED_" + std::string(tag_name);
    }
    CHKERR PetscOptionsInt(
        "-my_residual_stress_block",
        "block  in  mechanical (i.e. cracked) mesh  to which residual stress "
        "and density mapping are applied",
        "", residualStressBlock, &residualStressBlock, PETSC_NULL);
    MOFEM_LOG_C("CPWorld", Sev::inform,
                "### Input parameter: -my_residual_stress_block %d",
                residualStressBlock);
 
    CHKERR PetscOptionsInt("-my_density_block",
                           "block to which density is mapped", "",
                           densityMapBlock, &densityMapBlock, PETSC_NULL);
    MOFEM_LOG_C("CPWorld", Sev::inform,
                "### Input parameter: -my_density_block %d", densityMapBlock);
 
    char density_tag_name[255] = "RHO";
    CHKERR PetscOptionsString("-my_density_tag_name",
                              "name of density tag (RHO)", "", density_tag_name,
                              density_tag_name, 255, PETSC_NULL);
    MOFEM_LOG_C("CPWorld", Sev::inform,
                "### Input parameter: -my_density_tag_name %s",
                density_tag_name);
    mwlsRhoTagName = std::string(density_tag_name);
 
    {
      PetscBool flg;
      char file_name[255] = "add_cubit_meshsets.in";
      CHKERR PetscOptionsString("-meshsets_config",
                                "meshsets config  file name", "", file_name,
                                file_name, 255, &flg);
      if (flg == PETSC_TRUE) {
        ifstream f(file_name);
        if (!f.good()) {
          SETERRQ1(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
                   "File configuring meshsets ( %s ) cannot be open\n",
                   file_name);
        }
        configFile = file_name;
      }
    }
    {
      PetscBool flg;
      defaultMaterial = HOOKE;
      CHKERR PetscOptionsEList("-material", "Material type", "", materials_list,
                               LASTOP, materials_list[HOOKE], &defaultMaterial,
                               &flg);
      MOFEM_LOG_C("CPWorld", Sev::inform, "### Input parameter: -material %s",
                  materials_list[defaultMaterial]);
    }
    // Set paramaters with bone density like material
    {
      CHKERR PetscOptionsScalar("-my_rho0", "release energy", "", rHo0, &rHo0,
                                PETSC_NULL);
      MOFEM_LOG_C("CPWorld", Sev::inform, "### Input parameter: -my_rho0 %6.4e",
                  rHo0);
 
      CHKERR PetscOptionsScalar("-my_n_bone", "release energy", "", nBone,
                                &nBone, PETSC_NULL);
      MOFEM_LOG_C("CPWorld", Sev::inform,
                  "### Input parameter: -my_n_bone %6.4e", nBone);
    }
    {
      nbLoadSteps = 1;
      CHKERR PetscOptionsInt("-nb_load_steps", "number of load steps", "",
                             nbLoadSteps, &nbLoadSteps, PETSC_NULL);
      nbCutSteps = 0;
      CHKERR PetscOptionsInt("-nb_cut_steps", "number of cut mesh steps", "",
                             nbCutSteps, &nbCutSteps, PETSC_NULL);
      loadScale = 1;
      CHKERR PetscOptionsScalar("-load_scale", "load scale", "", loadScale,
                                &loadScale, PETSC_NULL);
      MOFEM_LOG_C("CPWorld", Sev::inform,
                  "### Input parameter: -load_scale %6.4e", loadScale);
    }
    {
      otherSideConstrains = PETSC_FALSE;
      CHKERR PetscOptionsBool(
          "-other_side_constrain",
          "Set surface constrain on other side of crack surface", "",
          otherSideConstrains, &otherSideConstrains, PETSC_NULL);
    }
    // Parameters for smoothing
    {
      smootherAlpha = 1;
      smootherGamma = 0;
      CHKERR PetscOptionsReal("-smoother_alpha", "Control mesh smoother", "",
                              smootherAlpha, &smootherAlpha, PETSC_NULL);
      CHKERR PetscOptionsReal("-my_gamma", "Controls mesh smoother", "",
                              smootherGamma, &smootherGamma, PETSC_NULL);
 
      MOFEM_LOG_C("CPWorld", Sev::inform,
                  "### Input parameter: -smoother_alpha %6.4e", smootherAlpha);
    }
    // Get Parameters for arc length
    {
      arcAlpha = 1;
      arcBeta = 0;
      arcS = 0;
      CHKERR PetscOptionsReal("-arc_alpha", "Arc length alpha", "", arcAlpha,
                              &arcAlpha, PETSC_NULL);
      CHKERR PetscOptionsReal("-arc_beta", "Arc length beta", "", arcBeta,
                              &arcBeta, PETSC_NULL);
      CHKERR PetscOptionsReal("-arc_s", "Arc length step size", "", arcS, &arcS,
                              PETSC_NULL);
      MOFEM_LOG_C("CPWorld", Sev::inform,
                  "### Input parameter: -arc_alpha %6.4e", arcAlpha);
      MOFEM_LOG_C("CPWorld", Sev::inform,
                  "### Input parameter: -arc_beta %6.4e", arcBeta);
      MOFEM_LOG_C("CPWorld", Sev::inform, "### Input parameter: -arc_s %6.4e",
                  arcS);
    }
    // Get parameters for contact
    {
      contactOrder = 2;
      contactLambdaOrder = 1;
      rValue = 1;
      cnValue = 1;
      ignoreContact = PETSC_FALSE;
      fixContactNodes = PETSC_FALSE;
      printContactState = PETSC_FALSE;
      contactOutputIntegPts = PETSC_TRUE;
      CHKERR PetscOptionsBool("-my_ignore_contact", "If true ignore contact",
                              "", ignoreContact, &ignoreContact, NULL);
      CHKERR PetscOptionsBool("-my_fix_contact_nodes",
                              "If true fix contact nodes", "", fixContactNodes,
                              &fixContactNodes, NULL);
      CHKERR PetscOptionsBool("-my_print_contact_state",
                              "If true print contact state", "",
                              printContactState, &printContactState, NULL);
      CHKERR PetscOptionsBool(
          "-my_contact_output_integ_pts",
          "If true output data at contact integration points", "",
          contactOutputIntegPts, &contactOutputIntegPts, NULL);
      CHKERR PetscOptionsReal("-my_r_value", "Contact regularisation parameter",
                              "", rValue, &rValue, PETSC_NULL);
      CHKERR PetscOptionsReal("-my_cn_value", "Contact augmentation parameter",
                              "", cnValue, &cnValue, PETSC_NULL);
      CHKERR PetscOptionsInt("-my_contact_order", "contact approximation order",
                             "", contactOrder, &contactOrder, PETSC_NULL);
      CHKERR PetscOptionsInt(
          "-my_contact_lambda_order", "contact Lagrange multipliers order", "",
          contactLambdaOrder, &contactLambdaOrder, PETSC_NULL);
      MOFEM_LOG_C("CPWorld", Sev::inform,
                  "### Input parameter: -my_cn_value %6.4e", cnValue);
      MOFEM_LOG_C("CPWorld", Sev::inform,
                  "### Input parameter: -my_contact_order %d", contactOrder);
      MOFEM_LOG_C("CPWorld", Sev::inform,
                  "### Input parameter: -my_cn_value %6.4e", cnValue);
      MOFEM_LOG_C("CPWorld", Sev::inform,
                  "### Input parameter: -my_contact_order %d", contactOrder);
    }
  }
  // Get Regine/Split/Cut options
  {
    doCutMesh = PETSC_FALSE;
    CHKERR PetscOptionsBool("-cut_mesh", "If true mesh is cut by surface", "",
                            doCutMesh, &doCutMesh, NULL);
    crackAccelerationFactor = 1;
    CHKERR PetscOptionsReal("-cut_factor", "Crack acceleration factor", "",
                            crackAccelerationFactor, &crackAccelerationFactor,
                            PETSC_NULL);
    doElasticWithoutCrack = PETSC_FALSE;
    CHKERR PetscOptionsBool(
        "-run_elastic_without_crack", "If true mesh is cut by surface", "",
        doElasticWithoutCrack, &doElasticWithoutCrack, NULL);
 
    MOFEM_LOG_C("CPWorld", Sev::inform, "### Input parameter: -cut_mesh %d",
                doCutMesh);
    MOFEM_LOG_C("CPWorld", Sev::inform,
                "### Input parameter: -cut_factor %6.4e",
                crackAccelerationFactor);
    MOFEM_LOG_C("CPWorld", Sev::inform,
                "### Input parameter: -run_elastic_without_crack %d",
                doElasticWithoutCrack);
  }
  CHKERR PetscOptionsInt("-post_proc_level", "level of output files", "",
                         postProcLevel, &postProcLevel, PETSC_NULL);
  MOFEM_LOG_C("CPWorld", Sev::verbose,
              "### Input parameter: -post_proc_level %6.4e", postProcLevel);
 
  addAnalyticalInternalStressOperators = PETSC_FALSE;
  CHKERR PetscOptionsBool(
      "-add_analytical_internal_stress_operators",
      "If true add analytical internal stress operators for mwls test", "",
      addAnalyticalInternalStressOperators,
      &addAnalyticalInternalStressOperators, NULL);
 
  solveEigenStressProblem = PETSC_FALSE;
  CHKERR PetscOptionsBool("-solve_eigen_problem", "If true solve eigen problem",
                          "", solveEigenStressProblem, &solveEigenStressProblem,
                          NULL);
  useEigenPositionsSimpleContact = PETSC_TRUE;
  CHKERR PetscOptionsBool(
      "-use_eigen_pos_simple_contact",
      "If true use eigen positions for matching meshes contact", "",
      useEigenPositionsSimpleContact, &useEigenPositionsSimpleContact, NULL);
 
  MOFEM_LOG_C("CPWorld", Sev::inform,
              "### Input parameter: -solve_eigen_problem %d",
              solveEigenStressProblem);
  MOFEM_LOG_C("CPWorld", Sev::inform,
              "### Input parameter: -use_eigen_pos_simple_contact %d",
              useEigenPositionsSimpleContact);
 
  // Partitioning
  partitioningWeightPower = 1;
  CHKERR PetscOptionsScalar("-part_weight_power", "Partitioning weight power",
                            "", partitioningWeightPower,
                            &partitioningWeightPower, PETSC_NULL);
 
  resetMWLSCoeffsEveryPropagationStep = PETSC_FALSE;
  CHKERR PetscOptionsBool(
      "-calc_mwls_coeffs_every_propagation_step",
      "If true recaulculate MWLS coefficients every propagation step", "",
      resetMWLSCoeffsEveryPropagationStep, &resetMWLSCoeffsEveryPropagationStep,
      NULL);
 
  ierr = PetscOptionsEnd();
  CHKERRQ(ierr);
 
  initialSmootherAlpha = smootherAlpha;
 
  CHKERRQ(ierr);
 
  // adding mums options
  char mumps_options[] =
      "-mat_mumps_icntl_14 800 -mat_mumps_icntl_24 1 -mat_mumps_icntl_13 1 "
      "-propagation_fieldsplit_0_mat_mumps_icntl_14 800 "
      "-propagation_fieldsplit_0_mat_mumps_icntl_24 1 "
      "-propagation_fieldsplit_0_mat_mumps_icntl_13 1 "
      "-propagation_fieldsplit_1_mat_mumps_icntl_14 800 "
      "-propagation_fieldsplit_1_mat_mumps_icntl_24 1 "
      "-propagation_fieldsplit_1_mat_mumps_icntl_13 1 "
      "-mat_mumps_icntl_20 0 "
      "-propagation_fieldsplit_0_mat_mumps_icntl_20 0 "
      "-propagation_fieldsplit_1_mat_mumps_icntl_20 0";
  CHKERR PetscOptionsInsertString(NULL, mumps_options);
 
  // Get options for interfaces
  CHKERR getInterface<CPMeshCut>()->getOptions();
  CHKERR getInterface<CPSolvers>()->getOptions();
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::tetsingReleaseEnergyCalculation() {
  PetscBool flg;
  enum tests {
    ANALYTICAL_MODE_1,
    CIRCULAR_PLATE_IN_INFINITE_BODY,
    ANGLE_ANALYTICAL_MODE_1,
    MWLS_INTERNAL_STRESS,
    MWLS_INTERNAL_STRESS_ANALYTICAL,
    CUT_CIRCLE_PLATE,
    NOTCH_WITH_DENSITY,
    LASTOP
  };
  const char *list[] = {"analytical_mode_1",
                        "circular_plate_in_infinite_body",
                        "angle_analytical_mode_1",
                        "mwls_internal_stress",
                        "mwls_internal_stress_analytical",
                        "cut_circle_plate",
                        "notch_with_density"};
 
  int choice_value = ANALYTICAL_MODE_1;
  double test_tol = 0;
  MoFEMFunctionBegin;
  CHKERR PetscOptionsBegin(mField.get_comm(), "", "Testing code", "none");
  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-test", list, LASTOP,
                              &choice_value, &flg);
  CHKERR PetscOptionsScalar("-test_tol", "test tolerance", "", test_tol,
                            &test_tol, PETSC_NULL);
  ierr = PetscOptionsEnd();
  CHKERRQ(ierr);
 
  if (flg == PETSC_TRUE) {
 
    if (mField.get_comm_rank() == 0) {
 
      switch (choice_value) {
      case ANALYTICAL_MODE_1: {
        double ave_g1 = 0;
        for (map<EntityHandle, double>::iterator mit = mapG1.begin();
             mit != mapG1.end(); mit++) {
          ave_g1 += mit->second;
        }
        ave_g1 /= mapG1.size();
        // Exact solution
        const double gc = 1e-5;
        // Tolerance
        const double error = fabs(ave_g1 - gc) / gc;
        CHKERR PetscPrintf(
            PETSC_COMM_SELF,
            "ave_gc = %6.4e gc = %6.4e tolerance = %6.4e error = %6.4e\n",
            ave_g1, gc, test_tol, error);
        if (error > test_tol || error != error) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                  "ANALYTICAL_MODE_1 test failed");
        }
      } break;
      case CIRCULAR_PLATE_IN_INFINITE_BODY: {
        const double gc = 1.1672e+01;
        double max_error = 0;
        for (map<EntityHandle, double>::iterator mit = mapG1.begin();
             mit != mapG1.end(); mit++) {
          const double g = mit->second;
          const double error = fabs(g - gc) / gc;
          max_error = (max_error < error) ? error : max_error;
          CHKERR PetscPrintf(
              PETSC_COMM_SELF,
              "g = %6.4e gc = %6.4e tolerance = %6.4e error = %6.4e\n", g, gc,
              test_tol, error);
        }
        if (max_error > test_tol || max_error != max_error)
          SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                  "CIRCULAR_PLATE_IN_INFINITE_BODY test failed");
 
      } break;
      case ANGLE_ANALYTICAL_MODE_1: {
        const double gc = 1e-5;
        double max_error = 0;
        for (map<EntityHandle, double>::iterator mit = mapG1.begin();
             mit != mapG1.end(); mit++) {
          const double g = mit->second;
          const double error = fabs(g - gc) / gc;
          max_error = (max_error < error) ? error : max_error;
          CHKERR PetscPrintf(
              PETSC_COMM_SELF,
              "g = %6.4e gc = %6.4e tolerance = %6.4e error = %6.4e\n", g, gc,
              test_tol, error);
        }
        if (max_error > test_tol || max_error != max_error) {
          SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                   "ANGLE_ANALYTICAL_MODE_1 test failed %6.4e > %6.4e",
                   max_error, test_tol);
        }
      } break;
      case MWLS_INTERNAL_STRESS: {
        const double gc =
            1.5527e+06; // this is max g1 for big problem (not extact solution)
        double max_error = 0;
        for (map<EntityHandle, double>::iterator mit = mapG1.begin();
             mit != mapG1.end(); mit++) {
          const double g = mit->second;
          const double error = fabs(g - gc) / gc;
          max_error = (max_error < error) ? error : max_error;
          CHKERR PetscPrintf(
              PETSC_COMM_SELF,
              "g = %6.4e gc = %6.4e tolerance = %6.4e error = %6.4e\n", g, gc,
              test_tol, error);
        }
        if (max_error > test_tol || max_error != max_error) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                  "MWLS_INTERNAL_STRESS test failed");
        }
      } break;
      case MWLS_INTERNAL_STRESS_ANALYTICAL: {
        const double gc =
            1.5698e+06; // this is max g1 for big problem (not extact solution)
        double max_error = 0;
        for (map<EntityHandle, double>::iterator mit = mapG1.begin();
             mit != mapG1.end(); mit++) {
          const double g = mit->second;
          const double error = fabs(g - gc) / gc;
          max_error = (max_error < error) ? error : max_error;
          CHKERR PetscPrintf(
              PETSC_COMM_SELF,
              "g = %6.4e gc = %6.4e tolerance = %6.4e error = %6.4e\n", g, gc,
              test_tol, error);
        }
        if (max_error > test_tol || max_error != max_error) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                  "MWLS_INTERNAL_STRESS_ANALYTICAL test failed");
        }
      } break;
      case CUT_CIRCLE_PLATE: {
        // This is simple regression test
        double max_g = 0;
        for (auto &m : mapG1) {
          const double g = m.second;
          max_g = fmax(g, max_g);
        }
        const double gc =
            7.5856e-02; // this is max g1 for big problem (not extact solution)
        double error = fabs(max_g - gc) / gc;
        if (error > test_tol || error != error) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                  "CUT_CIRCLE_PLATE test failed");
        }
      } break;
      case NOTCH_WITH_DENSITY: {
        // This is simple regression test
        double max_g = 0;
        for (map<EntityHandle, double>::iterator mit = mapG1.begin();
             mit != mapG1.end(); mit++) {
          const double g = mit->second;
          max_g = fmax(g, max_g);
        }
        const double gc = 1.6642e-04; // this is max g1 for a small problem (not
                                      // extact solution)
        double error = fabs(max_g - gc) / gc;
        if (error > test_tol || error != error) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                  "NOTCH_WITH_DENSITY test failed");
        }
      } break;
      }
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::saveEachPart(const std::string prefix,
                                              const Range &ents) {
  MoFEMFunctionBegin;
  std::ostringstream file;
  file << prefix << mField.get_comm_rank() << ".vtk";
  EntityHandle meshset_out;
  CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset_out);
  CHKERR mField.get_moab().add_entities(meshset_out, ents);
  CHKERR mField.get_moab().write_file(file.str().c_str(), "VTK", "",
                                      &meshset_out, 1, NULL, 0);
  CHKERR mField.get_moab().delete_entities(&meshset_out, 1);
  MoFEMFunctionReturn(0);
};
 
MoFEMErrorCode CrackPropagation::partitionMesh(BitRefLevel bit1,
                                               BitRefLevel bit2, int verb,
                                               const bool debug) {
  MoFEMFunctionBegin;
 
  moab::Interface &moab = mField.get_moab();
  CHKERR FractureMechanics::clean_pcomms(moab, moabCommWorld);
 
  // create map with layers with different weights for partitioning
  map<int, Range> layers;
  Range ref_nodes = crackFrontNodes;
  Range tets;
  CHKERR moab.get_adjacencies(ref_nodes, 3, false, tets,
                              moab::Interface::UNION);
  CHKERR mField.getInterface<BitRefManager>()->filterEntitiesByRefLevel(
      bit1, BitRefLevel().set(), tets);
  layers[0] = tets;
  for (int ll = 0; ll != refOrderAtTip; ll++) {
 
    CHKERR moab.get_connectivity(tets, ref_nodes, false);
    CHKERR moab.get_adjacencies(ref_nodes, 3, false, tets,
                                moab::Interface::UNION);
    CHKERR mField.getInterface<BitRefManager>()->filterEntitiesByRefLevel(
        bit1, BitRefLevel().set(), tets);
 
    layers[ll + 1] = subtract(tets, layers[ll]);
  }
  Range tets_bit_2;
  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
      bit2, BitRefLevel().set(), MBTET, tets_bit_2);
  int last = layers.size();
  auto rest_elements_in_the_bubble = subtract(tets_bit_2, layers[last - 1]);
  if (rest_elements_in_the_bubble.size())
    layers[last] = rest_elements_in_the_bubble;
 
  // create tag with weights
  Tag th_tet_weight;
  rval = moab.tag_get_handle("TETS_WEIGHT", th_tet_weight);
  if (rval == MB_SUCCESS) {
    CHKERR moab.tag_delete(th_tet_weight);
  }
  int def_val = 1;
  CHKERR moab.tag_get_handle("TETS_WEIGHT", 1, MB_TYPE_INTEGER, th_tet_weight,
                             MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
 
  for (unsigned int ll = 0; ll != layers.size(); ll++) {
    int weight = pow(layers.size() + 2 - ll, partitioningWeightPower);
    CHKERR moab.tag_clear_data(th_tet_weight, layers[ll], &weight);
  }
 
  // partition mesh
  Range ents;
  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
      bit1, BitRefLevel().set(), MBTET, ents);
  if (debug) {
    CHKERR saveEachPart("out_ents_to_partition_", ents);
  }
  CHKERR mField.getInterface<ProblemsManager>()->partitionMesh(
      ents, 3, 0, mField.get_comm_size(), &th_tet_weight, NULL, NULL, verb,
      false);
  CHKERR moab.tag_delete(th_tet_weight);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::resolveShared(const Range &tets,
                                               Range &proc_ents, const int verb,
                                               const bool debug) {
  moab::Interface &moab = mField.get_moab();
  Skinner skin(&mField.get_moab());
  MoFEMFunctionBegin;
 
  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
 
  proc_ents.clear();
  Range all_proc_ents;
 
  // get entities on processor
  Tag part_tag = pcomm->part_tag();
  Range tagged_sets;
  CHKERR mField.get_moab().get_entities_by_type_and_tag(
      0, MBENTITYSET, &part_tag, NULL, 1, tagged_sets, moab::Interface::UNION);
  for (Range::iterator mit = tagged_sets.begin(); mit != tagged_sets.end();
       mit++) {
    int part;
    CHKERR moab.tag_get_data(part_tag, &*mit, 1, &part);
    if (part == mField.get_comm_rank()) {
      CHKERR moab.get_entities_by_dimension(*mit, 3, proc_ents, true);
      CHKERR moab.get_entities_by_handle(*mit, all_proc_ents, true);
    }
  }
  proc_ents = intersect(proc_ents, tets);
 
  // get body skin
  Range tets_skin;
  CHKERR skin.find_skin(0, tets, false, tets_skin);
 
  // get skin on processor
  Range proc_ents_skin[4];
 
  Range contact_tets;
  if (!contactElements.empty() && !ignoreContact && !fixContactNodes) {
    BitRefLevel bit2_again = mapBitLevel["material_domain"];
    Range prisms_level;
    CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
        bit2_again, BitRefLevel().set(), MBPRISM, prisms_level);
 
    Range contact_prisms = intersect(prisms_level, contactElements);
 
    Range contact_prisms_to_remove =
        intersect(prisms_level, mortarContactElements);
    prisms_level = subtract(prisms_level, contact_prisms_to_remove);
 
    Range face_on_prism;
    CHKERR moab.get_adjacencies(contact_prisms, 2, false, face_on_prism,
                                moab::Interface::UNION);
 
    Range tris_on_prism = face_on_prism.subset_by_type(MBTRI);
 
    Range adj_vols_to_prisms;
    CHKERR moab.get_adjacencies(tris_on_prism, 3, false, adj_vols_to_prisms,
                                moab::Interface::UNION);
 
    Range tet_level;
    CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
        bit2_again, BitRefLevel().set(), MBTET, tet_level);
 
    Range contact_tets_from_vols = adj_vols_to_prisms.subset_by_type(MBTET);
 
    contact_tets = intersect(tet_level, contact_tets_from_vols);
 
    moab::Interface &moab = mField.get_moab();
    ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
    Range contact_tets_container = contact_tets;
    // FIXME: can we delete this?
    // Distribute prisms elements
    // {
 
    //   Range tris_level;
 
    //       //Find tris of tets
    //   CHKERR moab.get_adjacencies(contact_tets, 2, false, tris_level,
    //                               moab::Interface::UNION);
    //   for (Range::iterator mit = tagged_sets.begin(); mit !=
    //   tagged_sets.end();
    //        mit++) {
    //     Range part_tris;
    //     CHKERR moab.get_entities_by_type(*mit, MBTRI, part_tris);
    //     part_tris = intersect(part_tris, tris_level);
    //     Range adj_vols;
    //     CHKERR moab.get_adjacencies(part_tris, 3, false, adj_vols,
    //                                 moab::Interface::UNION);
    //     adj_vols = intersect(adj_vols, contact_tets_container);
    //     contact_tets_container = subtract(contact_tets_container, adj_vols);
    //     int part;
    //     CHKERR moab.tag_get_data(part_tag, &*mit, 1, &part);
    //     std::vector<int> tag(adj_vols.size(), part);
    //     CHKERR moab.tag_set_data(part_tag, adj_vols, &*tag.begin());
    //     CHKERR moab.add_entities(*mit, adj_vols);
    //   }
    // }
    contactTets = contact_tets;
    proc_ents.merge(contact_tets);
 
    Range contact_tets_vert;
    CHKERR moab.get_adjacencies(contact_tets, 0, false, contact_tets_vert,
                                moab::Interface::UNION);
    Range contact_tets_edges;
    CHKERR moab.get_adjacencies(contact_tets, 1, false, contact_tets_edges,
                                moab::Interface::UNION);
    Range contact_tets_faces;
    CHKERR moab.get_adjacencies(contact_tets, 2, false, contact_tets_faces,
                                moab::Interface::UNION);
    proc_ents_skin[2].merge(contact_tets_faces);
    proc_ents_skin[1].merge(contact_tets_edges);
    proc_ents_skin[0].merge(contact_tets_vert);
  }
 
  proc_ents_skin[3] = proc_ents;
  CHKERR skin.find_skin(0, proc_ents, false, proc_ents_skin[2]);
  // and get shared entities
  proc_ents_skin[2] = subtract(proc_ents_skin[2], tets_skin);
  // Note that all crack faces are shared
  proc_ents_skin[2].merge(crackFaces);
  proc_ents_skin[2].merge(contactSlaveFaces);
  proc_ents_skin[2].merge(contactMasterFaces);
  proc_ents_skin[2].merge(mortarContactSlaveFaces);
  proc_ents_skin[2].merge(mortarContactMasterFaces);
 
  CHKERR moab.get_adjacencies(proc_ents_skin[2], 1, false, proc_ents_skin[1],
                              moab::Interface::UNION);
  CHKERR moab.get_connectivity(proc_ents_skin[1], proc_ents_skin[0], true);
  proc_ents_skin[1].merge(crackFront);
  Range crack_faces_nodes;
  CHKERR moab.get_connectivity(crackFaces, crack_faces_nodes, true);
  proc_ents_skin[0].merge(crack_faces_nodes);
  // proc_ents_skin[0].merge(crackFrontNodes);
 
  if (mField.check_field("LAMBDA_ARC_LENGTH")) {
    EntityHandle lambda_meshset = mField.get_field_meshset("LAMBDA_ARC_LENGTH");
    Range arc_length_vertices;
    CHKERR moab.get_entities_by_type(lambda_meshset, MBVERTEX,
                                     arc_length_vertices, true);
    if (arc_length_vertices.size() != 1) {
      SETERRQ1(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
               "Should be one vertex in <LAMBDA_ARC_LENGTH> field but is %d",
               arc_length_vertices.size());
    }
    arcLengthVertex = arc_length_vertices[0];
    proc_ents_skin[0].merge(arc_length_vertices);
  } else {
    const double coords[] = {0, 0, 0};
    CHKERR moab.create_vertex(coords, arcLengthVertex);
    proc_ents_skin[0].insert(arcLengthVertex);
  }
  {
    Tag th_gid;
    CHKERR mField.get_moab().tag_get_handle(GLOBAL_ID_TAG_NAME, th_gid);
    int gid;
    ierr = mField.get_moab().get_number_entities_by_type(0, MBVERTEX, gid);
    CHKERR mField.get_moab().tag_set_data(th_gid, &arcLengthVertex, 1, &gid);
  }
 
  Range all_ents_no_on_part;
  CHKERR moab.get_entities_by_handle(0, all_ents_no_on_part, true);
  all_ents_no_on_part = subtract(all_ents_no_on_part, all_proc_ents);
  for (int dd = 0; dd != 4; ++dd) {
    all_ents_no_on_part = subtract(all_ents_no_on_part, proc_ents_skin[dd]);
  }
  all_ents_no_on_part = subtract(
      all_ents_no_on_part, all_ents_no_on_part.subset_by_type(MBENTITYSET));
  CHKERR mField.remove_ents_from_finite_element(all_ents_no_on_part);
  CHKERR mField.remove_ents_from_field(all_ents_no_on_part);
 
  {
    Range all_ents;
    CHKERR moab.get_entities_by_handle(0, all_ents);
    std::vector<unsigned char> pstat_tag(all_ents.size(), 0);
    CHKERR moab.tag_set_data(pcomm->pstatus_tag(), all_ents,
                             &*pstat_tag.begin());
  }
 
  CHKERR pcomm->resolve_shared_ents(0, proc_ents, 3, -1, proc_ents_skin);
 
  const RefEntity_multiIndex *refined_ents_ptr;
  ierr = mField.get_ref_ents(&refined_ents_ptr);
  if (debug) {
    std::ostringstream file_skin;
    file_skin << "out_skin_" << mField.get_comm_rank() << ".vtk";
    EntityHandle meshset_skin;
    CHKERR moab.create_meshset(MESHSET_SET, meshset_skin);
    CHKERR moab.add_entities(meshset_skin, proc_ents_skin[2]);
    CHKERR moab.add_entities(meshset_skin, proc_ents_skin[1]);
    CHKERR moab.add_entities(meshset_skin, proc_ents_skin[0]);
 
    CHKERR moab.write_file(file_skin.str().c_str(), "VTK", "", &meshset_skin,
                           1);
    std::ostringstream file_owned;
    file_owned << "out_owned_" << mField.get_comm_rank() << ".vtk";
    EntityHandle meshset_owned;
    CHKERR moab.create_meshset(MESHSET_SET, meshset_owned);
    CHKERR moab.add_entities(meshset_owned, proc_ents);
    CHKERR moab.write_file(file_owned.str().c_str(), "VTK", "", &meshset_owned,
                           1);
    CHKERR moab.delete_entities(&meshset_owned, 1);
  }
 
  if (debug) {
    Range shared_ents;
    // Get entities shared with all other processors
    CHKERR pcomm->get_shared_entities(-1, shared_ents);
    std::ostringstream file_shared_owned;
    file_shared_owned << "out_shared_owned_" << mField.get_comm_rank()
                      << ".vtk";
    EntityHandle meshset_shared_owned;
    CHKERR moab.create_meshset(MESHSET_SET, meshset_shared_owned);
    CHKERR moab.add_entities(meshset_shared_owned, shared_ents);
    CHKERR moab.write_file(file_shared_owned.str().c_str(), "VTK", "",
                           &meshset_shared_owned, 1);
    CHKERR moab.delete_entities(&meshset_shared_owned, 1);
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::resolveSharedBitRefLevel(const BitRefLevel bit,
                                                          const int verb,
                                                          const bool debug) {
  moab::Interface &moab = mField.get_moab();
  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
  MoFEMFunctionBegin;
  bitEnts.clear();
  bitProcEnts.clear();
  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
      bit, BitRefLevel().set(), MBTET, bitEnts);
 
  if (debug) {
    std::ostringstream file_owned;
    file_owned << "out_bit_" << mField.get_comm_rank() << ".vtk";
    EntityHandle meshset;
    CHKERR moab.create_meshset(MESHSET_SET, meshset);
    CHKERR moab.add_entities(meshset, bitEnts.subset_by_dimension(3));
    CHKERR moab.write_file(file_owned.str().c_str(), "VTK", "", &meshset, 1);
    CHKERR moab.delete_entities(&meshset, 1);
  }
  // Distribute prisms elements
  {
    Tag part_tag = pcomm->part_tag();
    Range tagged_sets;
    CHKERR mField.get_moab().get_entities_by_type_and_tag(
        0, MBENTITYSET, &part_tag, NULL, 1, tagged_sets,
        moab::Interface::UNION);
    Range tris_level;
    CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
        bit, BitRefLevel().set(), MBTRI, tris_level);
    Range prisms_sides;
    prisms_sides.merge(oneSideCrackFaces);
    prisms_sides.merge(contactSlaveFaces);
    prisms_sides.merge(mortarContactSlaveFaces);
 
    tris_level = intersect(tris_level, prisms_sides);
    Range prisms_level;
    CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
        bit, BitRefLevel().set(), MBPRISM, prisms_level);
 
    for (Range::iterator mit = tagged_sets.begin(); mit != tagged_sets.end();
         mit++) {
      Range part_tris;
      CHKERR moab.get_entities_by_type(*mit, MBTRI, part_tris);
      part_tris = intersect(part_tris, tris_level);
      Range adj_prisms;
      CHKERR moab.get_adjacencies(part_tris, 3, false, adj_prisms,
                                  moab::Interface::UNION);
      adj_prisms = intersect(adj_prisms, prisms_level);
      prisms_level = subtract(prisms_level, adj_prisms);
 
      int part;
      CHKERR moab.tag_get_data(part_tag, &*mit, 1, &part);
      std::vector<int> tag(adj_prisms.size(), part);
      CHKERR moab.tag_set_data(part_tag, adj_prisms, &*tag.begin());
      CHKERR moab.add_entities(*mit, adj_prisms);
    }
  }
 
  CHKERR resolveShared(bitEnts, bitProcEnts, verb, debug);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::buildElasticFields(
    const BitRefLevel bit, const BitRefLevel mask, const bool proc_only,
    const bool build_fields, const int verb, const bool debug) {
  moab::Interface &moab = mField.get_moab();
  MoFEMFunctionBegin;
 
  auto write_file = [&](const std::string file_name, Range ents) {
    MoFEMFunctionBegin;
    EntityHandle meshset;
    CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset);
    CHKERR mField.get_moab().add_entities(meshset, ents);
    if (mField.get_comm_rank() == 0) {
      CHKERR mField.get_moab().write_file(file_name.c_str(), "VTK", "",
                                          &meshset, 1);
    }
    CHKERR mField.get_moab().delete_entities(&meshset, 1);
    MoFEMFunctionReturn(0);
  };
 
  Range tets_level, tets_level_all;
  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
      bit, mask, MBTET, tets_level_all);
  if (proc_only) {
    tets_level = intersect(bitProcEnts, tets_level_all);
  }
  // That is ok, since all crack faces are shared
  Range prisms_level;
  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
      bit, mask, MBPRISM, prisms_level);
  if (proc_only) {
    prisms_level = intersect(prisms_level, bitProcEnts);
  }
  Range prisms_level_tris;
  rval = mField.get_moab().get_adjacencies(
      prisms_level, 2, false, prisms_level_tris, moab::Interface::UNION);
  // MoAB will report problems at crack front since quads are not well defined.
  if (rval != MB_SUCCESS && rval != MB_MULTIPLE_ENTITIES_FOUND) {
    CHKERRG(rval);
  } else {
    rval = MB_SUCCESS;
  }
  prisms_level_tris = intersect(prisms_level_tris, crackFaces);
  Range proc_ents = tets_level;
  proc_ents.merge(prisms_level_tris);
 
  CHKERR mField.add_field("SPATIAL_POSITION", H1, AINSWORTH_LOBATTO_BASE, 3,
                          MB_TAG_SPARSE, MF_ZERO, verb);
  if (solveEigenStressProblem) {
    CHKERR mField.add_field("EIGEN_SPATIAL_POSITIONS", H1,
                            AINSWORTH_LEGENDRE_BASE, 3, MB_TAG_SPARSE, MF_ZERO,
                            verb);
  }
  CHKERR mField.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE, 3,
                          MB_TAG_SPARSE, MF_ZERO, verb);
 
  // Set approx order
  Range dofs_ents;
  dofs_ents.merge(proc_ents);
  for (int dd = 1; dd != 3; ++dd) {
    CHKERR moab.get_adjacencies(proc_ents, dd, false, dofs_ents,
                                moab::Interface::UNION);
  }
 
  Range dofs_nodes;
  CHKERR moab.get_connectivity(proc_ents, dofs_nodes, true);
  dofs_nodes.merge(crackFrontNodes);
  dofs_ents.merge(crackFront);
 
  EntityHandle spatial_meshset = mField.get_field_meshset("SPATIAL_POSITION");
  Range
      ents_to_remove_from_field; // Entities which are no part of problem any
                                 // more. Very likely partition has been chaned.
  CHKERR mField.get_moab().get_entities_by_handle(spatial_meshset,
                                                  ents_to_remove_from_field);
  ents_to_remove_from_field = subtract(ents_to_remove_from_field, dofs_nodes);
  ents_to_remove_from_field = subtract(ents_to_remove_from_field, dofs_ents);
  CHKERR mField.remove_ents_from_field("SPATIAL_POSITION",
                                       ents_to_remove_from_field, verb);
  CHKERR mField.remove_ents_from_field("MESH_NODE_POSITIONS",
                                       ents_to_remove_from_field, verb);
  if (solveEigenStressProblem)
    CHKERR mField.remove_ents_from_field("EIGEN_SPATIAL_POSITIONS",
                                         ents_to_remove_from_field, verb);
 
  // Add entities to field
  auto add_spatial_field_ents = [&](const auto field_name) {
    MoFEMFunctionBegin;
    CHKERR mField.add_ents_to_field_by_type(tets_level, MBTET, field_name,
                                            verb);
    CHKERR mField.add_ents_to_field_by_type(prisms_level_tris, MBTRI,
                                            field_name, verb);
    CHKERR mField.add_ents_to_field_by_type(crackFrontNodes, MBVERTEX,
                                            field_name, verb);
    CHKERR mField.add_ents_to_field_by_type(crackFront, MBEDGE, field_name,
                                            verb);
    CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities(
        field_name, verb);
    MoFEMFunctionReturn(0);
  };
  CHKERR add_spatial_field_ents("SPATIAL_POSITION");
  if (solveEigenStressProblem) {
    CHKERR add_spatial_field_ents("EIGEN_SPATIAL_POSITIONS");
  }
 
  Range current_field_ents;
  CHKERR mField.get_moab().get_entities_by_handle(spatial_meshset,
                                                  current_field_ents);
  dofs_nodes = current_field_ents.subset_by_type(MBVERTEX);
  dofs_ents = subtract(current_field_ents, dofs_nodes);
  CHKERR mField.add_ents_to_field_by_type(dofs_nodes, MBVERTEX,
                                          "MESH_NODE_POSITIONS", verb);
 
  auto set_spatial_field_order = [&](const auto field_name) {
    MoFEMFunctionBegin;
    CHKERR mField.set_field_order(dofs_ents, field_name, approxOrder, verb);
    CHKERR mField.set_field_order(dofs_nodes, field_name, 1, verb);
    MoFEMFunctionReturn(0);
  };
  CHKERR set_spatial_field_order("SPATIAL_POSITION");
  if (solveEigenStressProblem) {
    CHKERR set_spatial_field_order("EIGEN_SPATIAL_POSITIONS");
  }
 
  CHKERR mField.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", geometryOrder,
                                verb);
  CHKERR mField.set_field_order(0, MBTRI, "MESH_NODE_POSITIONS", geometryOrder,
                                verb);
  CHKERR mField.set_field_order(0, MBTET, "MESH_NODE_POSITIONS", geometryOrder,
                                verb);
  CHKERR mField.set_field_order(dofs_nodes, "MESH_NODE_POSITIONS", 1, verb);
 
  if (debug) {
    CHKERR write_file("global_order_" + std::to_string(approxOrder) + ".vtk",
                      dofs_ents);
  }
 
  Range non_ref_ents;
  non_ref_ents.merge(dofs_ents);
 
  Range contact_faces;
  contact_faces.merge(contactSlaveFaces);
  contact_faces.merge(contactMasterFaces);
  contact_faces.merge(mortarContactMasterFaces);
  contact_faces.merge(mortarContactSlaveFaces);
 
  // Set ho approx. at crack front
  if (refOrderAtTip > 0) {
    Range proc_ents;
    CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
        bit, mask, MBTET, proc_ents);
    map<int, Range> layers;
    Range ref_nodes = crackFrontNodes;
    for (int ll = 0; ll != refOrderAtTip; ll++) {
      Range tets;
      CHKERR moab.get_adjacencies(ref_nodes, 3, false, tets,
                                  moab::Interface::UNION);
      tets = intersect(tets_level_all, tets);
      layers[ll] = tets;
      Range ref_order_ents;
      CHKERR moab.get_adjacencies(tets, 2, false, ref_order_ents,
                                  moab::Interface::UNION);
      CHKERR moab.get_adjacencies(tets, 1, false, ref_order_ents,
                                  moab::Interface::UNION);
      layers[ll].merge(ref_order_ents);
      CHKERR moab.get_connectivity(tets, ref_nodes, false);
 
      if (!contact_faces.empty()) {
        Range contact_layer = intersect(layers[ll], contact_faces);
        for (Range::iterator tit = contact_layer.begin();
             tit != contact_layer.end(); tit++) {
          Range contact_prisms, contact_tris;
          CHKERR moab.get_adjacencies(&*tit, 1, 3, false, contact_prisms,
                                      moab::Interface::UNION);
          contact_prisms = contact_prisms.subset_by_type(MBPRISM);
 
          Range contact_elements;
          if (!contactElements.empty() || !mortarContactElements.empty()) {
            contact_elements.merge(contactElements);
            contact_elements.merge(mortarContactElements);
            contact_elements = intersect(contact_prisms, contact_elements);
            if (contact_elements.empty()) {
              SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                      "Expecting adjacent contact prism(s), but none found");
            }
            CHKERR moab.get_adjacencies(contact_elements, 2, false,
                                        contact_tris, moab::Interface::UNION);
            contact_tris = contact_tris.subset_by_type(MBTRI);
            if (contact_tris.size() < 2) {
              SETERRQ1(
                  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "Expecting 2 or more adjacent contact tris, but %d found",
                  contact_tris.size());
            }
          }
 
          layers[ll].merge(contact_tris);
          Range contact_edges, contact_nodes;
          CHKERR moab.get_adjacencies(contact_tris, 1, false, contact_edges,
                                      moab::Interface::UNION);
          layers[ll].merge(contact_edges);
          CHKERR moab.get_connectivity(contact_tris, contact_nodes, false);
          ref_nodes.merge(contact_nodes);
        }
      }
    }
    int poo = 1;
    for (map<int, Range>::reverse_iterator mit = layers.rbegin();
         mit != layers.rend(); mit++, poo++) {
      mit->second = intersect(mit->second, dofs_ents);
      CHKERR mField.set_field_order(mit->second, "SPATIAL_POSITION",
                                    approxOrder + poo, verb);
      non_ref_ents = subtract(non_ref_ents, mit->second);
 
      if (debug) {
        CHKERR write_file("layer_order_" + std::to_string(approxOrder + poo) +
                              ".vtk",
                          mit->second);
      }
    }
  }
 
  if (!contactElements.empty() || !mortarContactElements.empty()) {
    if (contactOrder > approxOrder) {
      Range cont_ents, cont_adj_ent;
      cont_ents.merge(contactSlaveFaces);
      cont_ents.merge(contactMasterFaces);
      cont_ents.merge(mortarContactMasterFaces);
      cont_ents.merge(mortarContactSlaveFaces);
 
      CHKERR moab.get_adjacencies(cont_ents, 1, false, cont_adj_ent,
                                  moab::Interface::UNION);
      cont_ents.merge(cont_adj_ent);
      cont_ents = intersect(cont_ents, non_ref_ents);
      CHKERR mField.set_field_order(cont_ents, "SPATIAL_POSITION", contactOrder,
                                    verb);
      if (debug) {
        CHKERR write_file("contact_order_" + std::to_string(contactOrder) +
                              ".vtk",
                          cont_ents);
      }
    }
 
    // add Lagrange multipliers field for contact constraints on slave faces
    CHKERR mField.add_field("LAMBDA_CONTACT", H1, AINSWORTH_LEGENDRE_BASE, 1,
                            MB_TAG_SPARSE, MF_ZERO);
 
    if (!contactSlaveFaces.empty())
      CHKERR mField.add_ents_to_field_by_type(contactSlaveFaces, MBTRI,
                                              "LAMBDA_CONTACT");
 
    if (!mortarContactSlaveFaces.empty())
      CHKERR mField.add_ents_to_field_by_type(mortarContactSlaveFaces, MBTRI,
                                              "LAMBDA_CONTACT");
 
    CHKERR mField.set_field_order(0, MBTRI, "LAMBDA_CONTACT",
                                  contactLambdaOrder);
    CHKERR mField.set_field_order(0, MBEDGE, "LAMBDA_CONTACT",
                                  contactLambdaOrder);
    CHKERR mField.set_field_order(0, MBVERTEX, "LAMBDA_CONTACT", 1);
  }
 
  CHKERR mField.build_fields(verb);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::buildArcLengthField(const BitRefLevel bit,
                                                     const bool build_fields,
                                                     const int verb) {
  const RefEntity_multiIndex *refined_ents_ptr;
  const FieldEntity_multiIndex *field_ents_ptr;
  MoFEMFunctionBegin;
  CHKERR mField.get_ref_ents(&refined_ents_ptr);
  CHKERR mField.get_field_ents(&field_ents_ptr);
  if (mField.check_field("LAMBDA_ARC_LENGTH")) {
    auto vit = refined_ents_ptr->get<Ent_mi_tag>().find(arcLengthVertex);
    if (vit == refined_ents_ptr->get<Ent_mi_tag>().end()) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Vertex of arc-length field not found");
    }
    *(const_cast<RefEntity *>(vit->get())->getBitRefLevelPtr()) |= bit;
    auto fit = field_ents_ptr->get<Unique_mi_tag>().find(
        FieldEntity::getLocalUniqueIdCalculate(
            mField.get_field_bit_number("LAMBDA_ARC_LENGTH"), arcLengthVertex));
    if (fit == field_ents_ptr->get<Unique_mi_tag>().end())
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Arc length field entity not found");
 
    MOFEM_LOG("CPWorld", Sev::noisy) << **fit;
    MOFEM_LOG_C("CPWorld", Sev::inform, "Arc-Length field lambda = %3.4e",
                (*fit)->getEntFieldData()[0]);
    MoFEMFunctionReturnHot(0);
  }
 
  CHKERR mField.add_field("LAMBDA_ARC_LENGTH", NOFIELD, NOBASE, 1,
                          MB_TAG_SPARSE, MF_ZERO, verb);
  EntityHandle lambda_meshset = mField.get_field_meshset("LAMBDA_ARC_LENGTH");
  CHKERR mField.get_moab().add_entities(lambda_meshset, &arcLengthVertex, 1);
  // Add vertex to MoFEM database
  std::pair<RefEntity_multiIndex::iterator, bool> p_ent =
      const_cast<RefEntity_multiIndex *>(refined_ents_ptr)
          ->insert(boost::make_shared<RefEntity>(
              mField.get_basic_entity_data_ptr(), arcLengthVertex));
  *(const_cast<RefEntity *>(p_ent.first->get())->getBitRefLevelPtr()) |=
      (bit | BitRefLevel().set(BITREFLEVEL_SIZE - 2));
  if (build_fields) {
    CHKERR mField.build_fields(verb);
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::buildSurfaceFields(const BitRefLevel bit,
                                                    const bool proc_only,
                                                    const bool build_fields,
                                                    const int verb,
                                                    const bool debug) {
  moab::Interface &moab = mField.get_moab();
  MeshsetsManager *meshset_manager_ptr;
  MoFEMFunctionBegin;
  CHKERR mField.getInterface(meshset_manager_ptr);
 
  Range level_tris;
  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
      bit, BitRefLevel().set(), MBTRI, level_tris);
 
  Range body_skin_proc_tris;
  body_skin_proc_tris = intersect(bodySkin, level_tris);
 
  if (proc_only) {
    Range proc_tris;
    CHKERR moab.get_adjacencies(bitProcEnts, 2, false, proc_tris,
                                moab::Interface::UNION);
    body_skin_proc_tris = intersect(body_skin_proc_tris, proc_tris);
  }
 
  if (debug) {
    EntityHandle meshset;
    CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset);
    CHKERR mField.get_moab().add_entities(meshset, body_skin_proc_tris);
    std::string field_name;
    field_name = "out_proc_body_skin_" +
                 boost::lexical_cast<std::string>(mField.get_comm_rank()) +
                 ".vtk";
    CHKERR mField.get_moab().write_file(field_name.c_str(), "VTK", "", &meshset,
                                        1);
    CHKERR mField.get_moab().delete_entities(&meshset, 1);
  }
 
  std::string field_name =
      "LAMBDA_SURFACE" + boost::lexical_cast<std::string>(
                             getInterface<CPMeshCut>()->getSkinOfTheBodyId());
  Range dofs_nodes;
  CHKERR moab.get_connectivity(body_skin_proc_tris, dofs_nodes, true);
  if (proc_only)
    CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(dofs_nodes,
                                                                     QUIET);
 
  CHKERR mField.add_field(field_name, H1, AINSWORTH_LEGENDRE_BASE, 1,
                          MB_TAG_SPARSE, MF_ZERO, verb);
 
  Range ents_to_remove;
  CHKERR mField.get_field_entities_by_handle(field_name, ents_to_remove);
  ents_to_remove = subtract(ents_to_remove, dofs_nodes);
  CHKERR mField.remove_ents_from_field(field_name, ents_to_remove, verb);
 
  CHKERR mField.add_ents_to_field_by_type(dofs_nodes, MBVERTEX, field_name,
                                          verb);
  CHKERR mField.set_field_order(dofs_nodes, field_name, 1, verb);
 
  if (build_fields) {
    CHKERR mField.build_fields(verb);
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::buildEdgeFields(const BitRefLevel bit,
                                                 const bool proc_only,
                                                 const bool build_fields,
                                                 const int verb,
                                                 const bool debug) {
  moab::Interface &moab = mField.get_moab();
  auto *meshset_manager_ptr = mField.getInterface<MeshsetsManager>();
  auto *bit_ref_manager_ptr = mField.getInterface<BitRefManager>();
  MoFEMFunctionBegin;
 
  Range level_edges;
  CHKERR bit_ref_manager_ptr->getEntitiesByTypeAndRefLevel(
      bit, BitRefLevel().set(), MBEDGE, level_edges);
  Range edges_ents;
  if (meshset_manager_ptr->checkMeshset(
          getInterface<CPMeshCut>()->getEdgesBlockSet(), BLOCKSET)) {
    Range meshset_ents;
    CHKERR meshset_manager_ptr->getEntitiesByDimension(
        getInterface<CPMeshCut>()->getEdgesBlockSet(), BLOCKSET, 1,
        meshset_ents, true);
    edges_ents = intersect(meshset_ents, level_edges);
  }
 
  if (proc_only) {
    Range proc_edges;
    CHKERR moab.get_adjacencies(bitProcEnts, 1, false, proc_edges,
                                moab::Interface::UNION);
    edges_ents = intersect(edges_ents, proc_edges);
  }
 
  // Get nodes on other side faces
  Range other_side_crack_faces_nodes;
  rval = mField.get_moab().get_connectivity(otherSideCrackFaces,
                                            other_side_crack_faces_nodes, true);
  other_side_crack_faces_nodes =
      subtract(other_side_crack_faces_nodes, crackFrontNodes);
 
  Range dofs_nodes;
  CHKERR moab.get_connectivity(edges_ents, dofs_nodes, true);
  dofs_nodes = subtract(dofs_nodes, other_side_crack_faces_nodes);
  if (proc_only)
    CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(dofs_nodes,
                                                                     QUIET);
 
  std::string field_name = "LAMBDA_EDGE";
  CHKERR mField.add_field(field_name, H1, AINSWORTH_LEGENDRE_BASE, 2,
                          MB_TAG_SPARSE, MF_ZERO, verb);
 
  Range ents_to_remove;
  CHKERR mField.get_field_entities_by_handle(field_name, ents_to_remove);
  ents_to_remove = subtract(ents_to_remove, dofs_nodes);
  CHKERR mField.remove_ents_from_field(field_name, ents_to_remove, verb);
 
  CHKERR mField.add_ents_to_field_by_type(dofs_nodes, MBVERTEX, field_name,
                                          verb);
  CHKERR mField.set_field_order(dofs_nodes, field_name, 1, verb);
 
  auto get_fields_multi_index = [this]() {
    const Field_multiIndex *fields_ptr;
    CHKERR mField.get_fields(&fields_ptr);
    return *fields_ptr;
  };
  for (auto field : get_fields_multi_index()) {
    if (field->getName().compare(0, 14, "LAMBDA_SURFACE") == 0) {
      CHKERR mField.remove_ents_from_field(field->getName(), dofs_nodes, verb);
    }
  }
 
  if (build_fields) {
    CHKERR mField.build_fields(verb);
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::buildCrackSurfaceFieldId(
    const BitRefLevel bit, const bool proc_only, const bool build_fields,
    const int verb, const bool debug) {
  moab::Interface &moab = mField.get_moab();
  MeshsetsManager *meshset_manager_ptr;
  MoFEMFunctionBegin;
 
  CHKERR mField.getInterface(meshset_manager_ptr);
  Range crack_surface_ents = oneSideCrackFaces;
  Range level_tris;
  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
      bit, BitRefLevel().set(), MBTRI, level_tris);
  crack_surface_ents = intersect(crack_surface_ents, level_tris);
 
  if (proc_only) {
    Range proc_tris;
    CHKERR moab.get_adjacencies(bitProcEnts, 2, false, proc_tris,
                                moab::Interface::UNION);
    crack_surface_ents = intersect(crack_surface_ents, proc_tris);
  }
 
  if (debug) {
    EntityHandle meshset;
    CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset);
    CHKERR mField.get_moab().add_entities(meshset, crack_surface_ents);
    std::string field_name;
    field_name = "out_proc_crack_" +
                 boost::lexical_cast<std::string>(mField.get_comm_rank()) +
                 ".vtk";
    CHKERR mField.get_moab().write_file(field_name.c_str(), "VTK", "", &meshset,
                                        1);
    CHKERR mField.get_moab().delete_entities(&meshset, 1);
  }
 
  std::string field_name =
      "LAMBDA_SURFACE" + boost::lexical_cast<std::string>(
                             getInterface<CPMeshCut>()->getCrackSurfaceId());
  Range dofs_nodes;
  CHKERR moab.get_connectivity(crack_surface_ents, dofs_nodes, true);
  dofs_nodes = subtract(dofs_nodes, crackFrontNodes);
  if (proc_only)
    CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(dofs_nodes,
                                                                     QUIET);
 
  CHKERR mField.add_field(field_name, H1, AINSWORTH_LEGENDRE_BASE, 1,
                          MB_TAG_SPARSE, MF_ZERO, verb);
  Range ents_to_remove;
  CHKERR mField.get_field_entities_by_handle(field_name, ents_to_remove);
  ents_to_remove = subtract(ents_to_remove, dofs_nodes);
  CHKERR mField.remove_ents_from_field(field_name, ents_to_remove, verb);
 
  CHKERR mField.add_ents_to_field_by_type(dofs_nodes, MBVERTEX, field_name,
                                          verb);
  CHKERR mField.set_field_order(dofs_nodes, field_name, 1, verb);
 
  if (build_fields) {
    CHKERR mField.build_fields(verb);
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::buildBothSidesFieldId(
    const BitRefLevel bit_spatial, const BitRefLevel bit_material,
    const bool proc_only, const bool build_fields, const int verb,
    const bool debug) {
  MoFEMFunctionBegin;
 
  auto add_both_sides_field = [&](const std::string lambda_field_name,
                                  Range master_nodes, Range nodes_to_remove,
                                  bool add_ho) {
    MoFEMFunctionBegin;
    // This field set constrains that both side of crack surface have the same
    // node positions
    CHKERR mField.add_field(lambda_field_name, H1, AINSWORTH_LEGENDRE_BASE, 3,
                            MB_TAG_SPARSE, MF_ZERO);
 
    if (debug) {
      std::ostringstream file;
      file << lambda_field_name << ".vtk";
      EntityHandle meshset;
      CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset);
      CHKERR mField.get_moab().add_entities(meshset, master_nodes);
      CHKERR mField.get_moab().write_file(file.str().c_str(), "VTK", "",
                                          &meshset, 1);
      CHKERR mField.get_moab().delete_entities(&meshset, 1);
    }
 
    // remove old ents from field which are no longer part of the field
    Range ents_to_remove;
    CHKERR mField.get_field_entities_by_handle(lambda_field_name,
                                               ents_to_remove);
    ents_to_remove = subtract(ents_to_remove, master_nodes);
    CHKERR mField.remove_ents_from_field(lambda_field_name, ents_to_remove,
                                         verb);
 
    CHKERR mField.add_ents_to_field_by_type(master_nodes, MBVERTEX,
                                            lambda_field_name);
    if (approxOrder > 1 && add_ho) {
      CHKERR mField.add_ents_to_field_by_type(master_nodes, MBEDGE,
                                              lambda_field_name);
      CHKERR mField.add_ents_to_field_by_type(master_nodes, MBTRI,
                                              lambda_field_name);
    }
    CHKERR mField.set_field_order(master_nodes.subset_by_dimension(0),
                                  lambda_field_name, 1);
    if (approxOrder > 1 && add_ho) {
      CHKERR mField.set_field_order(master_nodes.subset_by_dimension(1),
                                    lambda_field_name, approxOrder);
      CHKERR mField.set_field_order(master_nodes.subset_by_dimension(2),
                                    lambda_field_name, approxOrder);
    }
 
    // Since we state constrains that mesh nodes, in material space, have to
    // move the same on the top and bottom side of crack surface, we have to
    // remove constraint from the all surface constrains fields, except
    // "LAMBDA_SURFACE", otherwise system will be over-constrained.
    const Field_multiIndex *fields_ptr;
    CHKERR mField.get_fields(&fields_ptr);
    for (auto const &field : (*fields_ptr)) {
      const std::string field_name = field->getName();
      if (field_name.compare(0, 14, "LAMBDA_SURFACE") == 0) {
        CHKERR mField.remove_ents_from_field(field_name, nodes_to_remove, verb);
      }
      if (field_name.compare(0, 11, "LAMBDA_EDGE") == 0) {
        CHKERR mField.remove_ents_from_field(field_name, nodes_to_remove, verb);
      }
    }
    MoFEMFunctionReturn(0);
  };
 
  auto write_file = [&](const std::string file_name, Range ents) {
    MoFEMFunctionBegin;
    EntityHandle meshset;
    CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset);
    CHKERR mField.get_moab().add_entities(meshset, ents);
    if (mField.get_comm_rank() == 0) {
      CHKERR mField.get_moab().write_file(file_name.c_str(), "VTK", "",
                                          &meshset, 1);
    }
    CHKERR mField.get_moab().delete_entities(&meshset, 1);
    MoFEMFunctionReturn(0);
  };
 
  auto get_prisms_level_nodes = [&](auto bit) {
    Range prisms_level;
    CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
        bit, BitRefLevel().set(), MBPRISM, prisms_level);
    if (proc_only) {
      prisms_level = intersect(bitProcEnts, prisms_level);
    }
    Range prisms_level_nodes;
    CHKERR mField.get_moab().get_connectivity(prisms_level, prisms_level_nodes,
                                              true);
    return prisms_level_nodes;
  };
 
  auto get_other_side_crack_faces_nodes = [&](auto prisms_level_nodes) {
    Range other_side_crack_faces_nodes;
    CHKERR mField.get_moab().get_connectivity(
        otherSideCrackFaces, other_side_crack_faces_nodes, true);
    other_side_crack_faces_nodes =
        intersect(other_side_crack_faces_nodes, prisms_level_nodes);
    other_side_crack_faces_nodes =
        subtract(other_side_crack_faces_nodes, crackFrontNodes);
 
    return other_side_crack_faces_nodes;
  };
 
  auto prisms_level_nodes = get_prisms_level_nodes(bit_material);
  auto other_side_crack_faces_nodes =
      get_other_side_crack_faces_nodes(prisms_level_nodes);
 
  CHKERR add_both_sides_field("LAMBDA_BOTH_SIDES", other_side_crack_faces_nodes,
                              other_side_crack_faces_nodes, false);
 
  if (solveEigenStressProblem) {
    Range other_side_crack_faces_edges;
    if (approxOrder > 1) {
      CHKERR mField.get_moab().get_adjacencies(otherSideCrackFaces, 1, false,
                                               other_side_crack_faces_edges,
                                               moab::Interface::UNION);
      other_side_crack_faces_edges =
          subtract(other_side_crack_faces_edges, crackFront);
      other_side_crack_faces_edges.merge(otherSideCrackFaces);
    }
    other_side_crack_faces_edges.merge(
        get_other_side_crack_faces_nodes(get_prisms_level_nodes(bit_spatial)));
    CHKERR add_both_sides_field("LAMBDA_CLOSE_CRACK",
                                other_side_crack_faces_edges, Range(), true);
  }
 
  if (!contactElements.empty() && !ignoreContact && !fixContactNodes) {
    contactBothSidesMasterNodes.clear();
    CHKERR mField.get_moab().get_connectivity(
        contactBothSidesMasterFaces, contactBothSidesMasterNodes, true);
    contactBothSidesMasterNodes =
        intersect(contactBothSidesMasterNodes, prisms_level_nodes);
    contactBothSidesMasterNodes =
        subtract(contactBothSidesMasterNodes, other_side_crack_faces_nodes);
    Range contact_both_sides_slave_nodes;
    for (Range::iterator nit = contactBothSidesMasterNodes.begin();
         nit != contactBothSidesMasterNodes.end(); nit++) {
      Range contact_prisms;
      CHKERR mField.get_moab().get_adjacencies(
          &*nit, 1, 3, false, contact_prisms, moab::Interface::UNION);
      contact_prisms = intersect(contact_prisms, contactElements);
      EntityHandle prism = contact_prisms.front();
      const EntityHandle *conn;
      int side_number, other_side_number, sense, offset, number_nodes = 0;
      CHKERR mField.get_moab().get_connectivity(prism, conn, number_nodes);
      CHKERR mField.get_moab().side_number(prism, *nit, side_number, sense,
                                           offset);
      if (side_number < 3)
        contact_both_sides_slave_nodes.insert(conn[side_number + 3]);
      else
        contact_both_sides_slave_nodes.insert(conn[side_number - 3]);
    }
    contact_both_sides_slave_nodes =
        intersect(contact_both_sides_slave_nodes, prisms_level_nodes);
 
    if (debug) {
      CHKERR write_file("contact_both_sides_slave_nodes.vtk",
                        contact_both_sides_slave_nodes);
    }
 
    CHKERR add_both_sides_field("LAMBDA_BOTH_SIDES_CONTACT",
                                contactBothSidesMasterNodes,
                                contact_both_sides_slave_nodes, false);
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::buildCrackFrontFieldId(const BitRefLevel bit,
                                                        const bool build_fields,
                                                        const int verb,
                                                        const bool debug) {
  moab::Interface &moab = mField.get_moab();
  MoFEMFunctionBegin;
 
  Range level_edges;
  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
      bit, BitRefLevel().set(), MBEDGE, level_edges);
  Range crack_front_edges = intersect(crackFront, level_edges);
  Range crack_front_nodes;
  CHKERR moab.get_connectivity(crack_front_edges, crack_front_nodes, true);
 
  CHKERR mField.add_field("LAMBDA_CRACKFRONT_AREA", H1, AINSWORTH_LEGENDRE_BASE,
                          1, MB_TAG_SPARSE, MF_ZERO, verb);
  // Add dofs to LAMBDA_CRACKFRONT_AREA field
  {
    Range ents_to_remove;
    CHKERR mField.get_field_entities_by_handle("LAMBDA_CRACKFRONT_AREA",
                                               ents_to_remove);
    ents_to_remove = subtract(ents_to_remove, crack_front_nodes);
    CHKERR mField.remove_ents_from_field("LAMBDA_CRACKFRONT_AREA",
                                         ents_to_remove, verb);
    CHKERR mField.add_ents_to_field_by_type(crack_front_nodes, MBVERTEX,
                                            "LAMBDA_CRACKFRONT_AREA", verb);
    CHKERR mField.set_field_order(crack_front_nodes, "LAMBDA_CRACKFRONT_AREA",
                                  1, verb);
  }
 
  CHKERR mField.add_field("LAMBDA_CRACKFRONT_AREA_TANGENT", H1,
                          AINSWORTH_LEGENDRE_BASE, 1, MB_TAG_SPARSE, MF_ZERO,
                          verb);
  // Add dofs to LAMBDA_CRACKFRONT_AREA_TANGENT field
  {
    // Remove dofs on crack front ends
    Skinner skin(&mField.get_moab());
    Range skin_crack_front;
    rval = skin.find_skin(0, crack_front_edges, false, skin_crack_front);
    Range field_ents = subtract(crack_front_nodes, skin_crack_front);
    // remove old ents from field which are no longer part of the field
    Range ents_to_remove;
    CHKERR mField.get_field_entities_by_handle("LAMBDA_CRACKFRONT_AREA_TANGENT",
                                               ents_to_remove);
    ents_to_remove = subtract(ents_to_remove, field_ents);
    CHKERR mField.remove_ents_from_field("LAMBDA_CRACKFRONT_AREA_TANGENT",
                                         ents_to_remove, verb);
    CHKERR mField.add_ents_to_field_by_type(
        field_ents, MBVERTEX, "LAMBDA_CRACKFRONT_AREA_TANGENT", verb);
    CHKERR mField.set_field_order(field_ents, "LAMBDA_CRACKFRONT_AREA_TANGENT",
                                  1, verb);
  }
 
  if (build_fields) {
    CHKERR mField.build_fields(verb);
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::declareElasticFE(
    const BitRefLevel bit1, const BitRefLevel mask1, const BitRefLevel bit2,
    const BitRefLevel mask2, const bool add_forces, const bool proc_only,
    const int verb) {
  MoFEMFunctionBegin;
 
  Range tets_level1;
  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
      bit1, mask1, MBTET, tets_level1);
  if (proc_only) {
    tets_level1 = intersect(bitProcEnts, tets_level1);
  }
  Range tets_not_level2;
  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
      bit2, mask2, MBTET, tets_not_level2);
  tets_not_level2 = subtract(tets_level1, tets_not_level2);
 
  // Add finite elements
  CHKERR mField.add_finite_element("ELASTIC", MF_ZERO, verb);
  CHKERR mField.modify_finite_element_add_field_row("ELASTIC",
                                                    "SPATIAL_POSITION");
  CHKERR mField.modify_finite_element_add_field_col("ELASTIC",
                                                    "SPATIAL_POSITION");
  CHKERR mField.modify_finite_element_add_field_data("ELASTIC",
                                                     "SPATIAL_POSITION");
  CHKERR mField.modify_finite_element_add_field_data("ELASTIC",
                                                     "MESH_NODE_POSITIONS");
  // Add ents and build finite elements. Only build elements which has been
  // added.
  {
    Range current_ents_with_fe;
    CHKERR mField.get_finite_element_entities_by_handle("ELASTIC",
                                                        current_ents_with_fe);
    Range ents_to_remove;
    ents_to_remove = subtract(current_ents_with_fe, tets_level1);
 
    CHKERR mField.remove_ents_from_finite_element("ELASTIC", ents_to_remove);
    CHKERR mField.add_ents_to_finite_element_by_type(tets_level1, MBTET,
                                                     "ELASTIC");
    CHKERR mField.build_finite_elements("ELASTIC", &tets_level1, VERBOSE);
  }
 
  // Add finite elements not on crack front
  CHKERR mField.add_finite_element("ELASTIC_NOT_ALE", MF_ZERO, verb);
  CHKERR mField.modify_finite_element_add_field_row("ELASTIC_NOT_ALE",
                                                    "SPATIAL_POSITION");
  CHKERR mField.modify_finite_element_add_field_col("ELASTIC_NOT_ALE",
                                                    "SPATIAL_POSITION");
  CHKERR mField.modify_finite_element_add_field_data("ELASTIC_NOT_ALE",
                                                     "SPATIAL_POSITION");
  CHKERR mField.modify_finite_element_add_field_data("ELASTIC_NOT_ALE",
                                                     "MESH_NODE_POSITIONS");
 
  // Add ents and build finite elements
  {
    Range current_ents_with_fe;
    CHKERR mField.get_finite_element_entities_by_handle("ELASTIC_NOT_ALE",
                                                        current_ents_with_fe);
    Range ents_to_remove;
    ents_to_remove = subtract(current_ents_with_fe, tets_not_level2);
 
    CHKERR mField.remove_ents_from_finite_element("ELASTIC_NOT_ALE",
                                                  ents_to_remove);
    CHKERR mField.add_ents_to_finite_element_by_type(tets_not_level2, MBTET,
                                                     "ELASTIC_NOT_ALE");
    CHKERR mField.build_finite_elements("ELASTIC_NOT_ALE", &tets_not_level2,
                                        verb);
  }
 
  // Add skin elements for faster postprocessing
  {
    CHKERR mField.add_finite_element("SKIN", MF_ZERO, verb);
    Skinner skin(&mField.get_moab());
    Range proc_tris, skin_tris, tets_level;
    CHKERR mField.get_moab().get_adjacencies(bitProcEnts, 2, false, proc_tris,
                                             moab::Interface::UNION);
 
    CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
        bit1, mask1, MBTET, tets_level);
 
    CHKERR skin.find_skin(0, tets_level, false, skin_tris);
    skin_tris = intersect(skin_tris, proc_tris);
 
    CHKERR mField.add_ents_to_finite_element_by_type(skin_tris, MBTRI, "SKIN");
 
    CHKERR mField.modify_finite_element_add_field_row("SKIN",
                                                      "SPATIAL_POSITION");
    CHKERR mField.modify_finite_element_add_field_col("SKIN",
                                                      "SPATIAL_POSITION");
    CHKERR mField.modify_finite_element_add_field_data("SKIN",
                                                       "SPATIAL_POSITION");
    CHKERR mField.modify_finite_element_add_field_data("SKIN",
                                                       "MESH_NODE_POSITIONS");
    if (solveEigenStressProblem)
      CHKERR mField.modify_finite_element_add_field_data(
          "SKIN", "EIGEN_SPATIAL_POSITIONS");
    CHKERR mField.build_finite_elements("SKIN");
  }
  // Add spring boundary condition applied on surfaces.
  // This is only declaration not implementation.
  CHKERR MetaSpringBC::addSpringElements(mField, "SPATIAL_POSITION",
                                         "MESH_NODE_POSITIONS");
  CHKERR mField.build_finite_elements("SPRING");
 
  auto cn_value_ptr = boost::make_shared<double>(cnValue);
 
  contactProblem = boost::shared_ptr<SimpleContactProblem>(
      new SimpleContactProblem(mField, cn_value_ptr));
  // add fields to the global matrix by adding the element
 
  CHKERR contactProblem->addContactElement(
      "CONTACT", "SPATIAL_POSITION", "LAMBDA_CONTACT", contactElements,
      solveEigenStressProblem, "EIGEN_SPATIAL_POSITIONS");
 
  CHKERR mField.build_finite_elements("CONTACT");
 
  mortarContactProblemPtr =
      boost::shared_ptr<MortarContactProblem>(new MortarContactProblem(
          mField, contactSearchMultiIndexPtr, cn_value_ptr));
 
  CHKERR mortarContactProblemPtr->addMortarContactElement(
      "MORTAR_CONTACT", "SPATIAL_POSITION", "LAMBDA_CONTACT",
      mortarContactElements, "MESH_NODE_POSITIONS", solveEigenStressProblem,
      "EIGEN_SPATIAL_POSITIONS");
 
  CHKERR mField.build_finite_elements("MORTAR_CONTACT");
 
  Range all_slave_tris;
  all_slave_tris.merge(mortarContactSlaveFaces);
  all_slave_tris.merge(contactSlaveFaces);
 
  CHKERR contactProblem->addPostProcContactElement(
      "CONTACT_POST_PROC", "SPATIAL_POSITION", "LAMBDA_CONTACT",
      "MESH_NODE_POSITIONS", all_slave_tris);
  CHKERR mField.build_finite_elements("CONTACT_POST_PROC");
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::declareArcLengthFE(const BitRefLevel bit,
                                                    const int verb) {
  const EntFiniteElement_multiIndex *fe_ent_ptr;
  MoFEMFunctionBegin;
  if (mField.check_finite_element("ARC_LENGTH")) {
    Range arc_fe_meshsets;
    for (auto fit = mField.get_fe_by_name_begin("ARC_LENGTH");
         fit != mField.get_fe_by_name_end("ARC_LENGTH"); ++fit) {
      arc_fe_meshsets.insert(fit->get()->getEnt());
      if (auto ptr = fit->get()->getPartProcPtr())
        *ptr = 0;
    }
    CHKERR mField.getInterface<BitRefManager>()->addBitRefLevel(arc_fe_meshsets,
                                                                bit);
    CHKERR mField.build_finite_elements("ARC_LENGTH", NULL, verb);
    MoFEMFunctionReturnHot(0);
  }
 
  CHKERR mField.add_finite_element("ARC_LENGTH", MF_ZERO, verb);
  CHKERR mField.modify_finite_element_add_field_row("ARC_LENGTH",
                                                    "LAMBDA_ARC_LENGTH");
  CHKERR mField.modify_finite_element_add_field_col("ARC_LENGTH",
                                                    "LAMBDA_ARC_LENGTH");
  // elem data
  CHKERR mField.modify_finite_element_add_field_data("ARC_LENGTH",
                                                     "LAMBDA_ARC_LENGTH");
 
  // finally add created meshset to the ARC_LENGTH finite element
  EntityHandle meshset_fe_arc_length;
  {
    CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset_fe_arc_length);
    CHKERR mField.get_moab().add_entities(meshset_fe_arc_length,
                                          &arcLengthVertex, 1);
    CHKERR mField.getInterface<BitRefManager>()->setBitLevelToMeshset(
        meshset_fe_arc_length, bit | BitRefLevel().set(BITREFLEVEL_SIZE - 2));
  }
  CHKERR mField.add_ents_to_finite_element_by_MESHSET(meshset_fe_arc_length,
                                                      "ARC_LENGTH", false);
  CHKERR mField.build_finite_elements("ARC_LENGTH", NULL, verb);
  
  for (auto fit = mField.get_fe_by_name_begin("ARC_LENGTH");
       fit != mField.get_fe_by_name_end("ARC_LENGTH"); ++fit) {
    if (auto ptr = fit->get()->getPartProcPtr())
      *ptr = 0;
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::declareExternalForcesFE(const BitRefLevel bit,
                                                         const BitRefLevel mask,
                                                         const bool proc_only) {
  moab::Interface &moab = mField.get_moab();
  MoFEMFunctionBegin;
 
  Range tets_level;
  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
      bit, mask, MBTET, tets_level);
  if (proc_only) {
    tets_level = intersect(bitProcEnts, tets_level);
  }
 
  // Set approx order
  Range dofs_ents;
  dofs_ents.merge(tets_level);
  for (int dd = 1; dd != 3; dd++) {
    CHKERR moab.get_adjacencies(tets_level, dd, false, dofs_ents,
                                moab::Interface::UNION);
  }
  Range dofs_nodes;
  rval = moab.get_connectivity(tets_level, dofs_nodes, true);
  CHKERRQ_MOAB(rval);
  Range lower_dim_ents = unite(dofs_ents, dofs_nodes);
 
#ifdef __ANALITICAL_TRACTION__
  {
    CHKERR mField.add_finite_element("ANALITICAL_TRACTION", MF_ZERO);
    CHKERR mField.modify_finite_element_add_field_row("ANALITICAL_TRACTION",
                                                      "SPATIAL_POSITION");
    CHKERR mField.modify_finite_element_add_field_col("ANALITICAL_TRACTION",
                                                      "SPATIAL_POSITION");
    CHKERR mField.modify_finite_element_add_field_data("ANALITICAL_TRACTION",
                                                       "SPATIAL_POSITION");
    CHKERR mField.modify_finite_element_add_field_data("ANALITICAL_TRACTION",
                                                       "MESH_NODE_POSITIONS");
    MeshsetsManager *meshset_manager_ptr;
    CHKERR mField.getInterface(meshset_manager_ptr);
    if (meshset_manager_ptr->checkMeshset("ANALITICAL_TRACTION")) {
      const CubitMeshSets *cubit_meshset_ptr;
      meshset_manager_ptr->getCubitMeshsetPtr("ANALITICAL_TRACTION",
                                              &cubit_meshset_ptr);
      Range tris;
      CHKERR meshset_manager_ptr->getEntitiesByDimension(
          cubit_meshset_ptr->getMeshsetId(), BLOCKSET, 2, tris, true);
      tris = intersect(tris, lower_dim_ents);
      CHKERR mField.add_ents_to_finite_element_by_type(tris, MBTRI,
                                                       "ANALITICAL_TRACTION");
    }
    CHKERR mField.build_finite_elements("ANALITICAL_TRACTION", &lower_dim_ents);
  }
#endif //__ANALITICAL_TRACTION__
 
  CHKERR mField.add_finite_element("FORCE_FE", MF_ZERO);
  CHKERR mField.modify_finite_element_add_field_row("FORCE_FE",
                                                    "SPATIAL_POSITION");
  CHKERR mField.modify_finite_element_add_field_col("FORCE_FE",
                                                    "SPATIAL_POSITION");
  CHKERR mField.modify_finite_element_add_field_data("FORCE_FE",
                                                     "SPATIAL_POSITION");
  CHKERR mField.modify_finite_element_add_field_data("FORCE_FE",
                                                     "MESH_NODE_POSITIONS");
  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(mField, NODESET | FORCESET,
                                                  it)) {
    Range tris;
    CHKERR mField.get_moab().get_entities_by_type(it->meshset, MBTRI, tris,
                                                  true);
    Range edges;
    CHKERR mField.get_moab().get_entities_by_type(it->meshset, MBEDGE, edges,
                                                  true);
    Range tris_edges;
    CHKERR mField.get_moab().get_adjacencies(tris, 1, false, tris_edges,
                                             moab::Interface::UNION);
    Range tris_nodes;
    CHKERR mField.get_moab().get_connectivity(tris, tris_nodes);
    Range edges_nodes;
    CHKERR mField.get_moab().get_connectivity(edges, edges_nodes);
    Range nodes;
    CHKERR mField.get_moab().get_entities_by_type(it->meshset, MBVERTEX, nodes,
                                                  true);
    nodes = subtract(nodes, tris_nodes);
    nodes = subtract(nodes, edges_nodes);
    nodes = intersect(lower_dim_ents, nodes);
    edges = subtract(edges, tris_edges);
    edges = intersect(lower_dim_ents, edges);
    tris = intersect(lower_dim_ents, tris);
    CHKERR mField.add_ents_to_finite_element_by_type(tris, MBTRI, "FORCE_FE");
    CHKERR mField.add_ents_to_finite_element_by_type(edges, MBEDGE, "FORCE_FE");
    CHKERR mField.add_ents_to_finite_element_by_type(nodes, MBVERTEX,
                                                     "FORCE_FE");
  }
 
  // Reading forces from BLOCKSET (for bone meshes)
  const string block_set_force_name("FORCE");
  // search for block named FORCE and add its attributes to FORCE_FE element
  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
    if (it->getName().compare(0, block_set_force_name.length(),
                              block_set_force_name) == 0) {
      std::vector<double> mydata;
      CHKERR it->getAttributes(mydata);
 
      VectorDouble force(mydata.size());
      for (unsigned int ii = 0; ii < mydata.size(); ii++) {
        force[ii] = mydata[ii];
      }
      if (force.size() != 3) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Force vector not given");
      }
      Range tris;
      CHKERR mField.get_moab().get_entities_by_type(it->meshset, MBTRI, tris,
                                                    true);
      CHKERR mField.add_ents_to_finite_element_by_type(tris, MBTRI, "FORCE_FE");
    }
  }
 
  // CHKERR MetaNodalForces::addElement(mField,"SPATIAL_POSITION");
  // CHKERR MetaEdgeForces::addElement(mField,"SPATIAL_POSITION");
  CHKERR mField.build_finite_elements("FORCE_FE", &lower_dim_ents);
  CHKERR mField.add_finite_element("PRESSURE_FE", MF_ZERO);
  CHKERR mField.modify_finite_element_add_field_row("PRESSURE_FE",
                                                    "SPATIAL_POSITION");
  CHKERR mField.modify_finite_element_add_field_col("PRESSURE_FE",
                                                    "SPATIAL_POSITION");
  CHKERR mField.modify_finite_element_add_field_data("PRESSURE_FE",
                                                     "SPATIAL_POSITION");
  CHKERR mField.modify_finite_element_add_field_data("PRESSURE_FE",
                                                     "MESH_NODE_POSITIONS");
  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(mField, SIDESET | PRESSURESET,
                                                  it)) {
    Range tris;
    CHKERR mField.get_moab().get_entities_by_type(it->meshset, MBTRI, tris,
                                                  true);
    tris = intersect(lower_dim_ents, tris);
    CHKERR mField.add_ents_to_finite_element_by_type(tris, MBTRI,
                                                     "PRESSURE_FE");
  }
  const string block_set_linear_pressure_name("LINEAR_PRESSURE");
  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
    if (it->getName().compare(0, block_set_linear_pressure_name.length(),
                              block_set_linear_pressure_name) == 0) {
      Range tris;
      CHKERR mField.get_moab().get_entities_by_type(it->meshset, MBTRI, tris,
                                                    true);
      tris = intersect(lower_dim_ents, tris);
      CHKERR mField.add_ents_to_finite_element_by_type(tris, MBTRI,
                                                       "PRESSURE_FE");
    }
  }
 
  // search for block named PRESSURE and add its attributes to PRESSURE_FE
  // element (for bone meshes)
  const string block_set_pressure_name("PRESSURE");
  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
    if (it->getName().compare(0, block_set_pressure_name.length(),
                              block_set_pressure_name) == 0) {
      std::vector<double> mydata;
      CHKERR it->getAttributes(mydata);
      VectorDouble pressure(mydata.size());
      for (unsigned int ii = 0; ii < mydata.size(); ii++) {
        pressure[ii] = mydata[ii];
      }
      if (pressure.empty()) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Pressure not given");
      }
      Range tris;
      CHKERR mField.get_moab().get_entities_by_type(it->meshset, MBTRI, tris,
                                                    true);
 
      CHKERR
      mField.add_ents_to_finite_element_by_type(tris, MBTRI, "PRESSURE_FE");
    }
  }
 
  CHKERR mField.build_finite_elements("PRESSURE_FE", &lower_dim_ents);
 
#ifdef __ANALITICAL_DISPLACEMENT__
  {
    MeshsetsManager *meshset_manager_ptr;
    CHKERR mField.getInterface(meshset_manager_ptr);
    if (meshset_manager_ptr->checkMeshset("ANALITICAL_DISP")) {
      const CubitMeshSets *cubit_meshset_ptr;
      meshset_manager_ptr->getCubitMeshsetPtr("ANALITICAL_DISP",
                                              &cubit_meshset_ptr);
      Range tris;
      CHKERR meshset_manager_ptr->getEntitiesByDimension(
          cubit_meshset_ptr->getMeshsetId(), BLOCKSET, 2, tris, true);
      if (!tris.empty()) {
        CHKERR mField.add_finite_element("BC_FE", MF_ZERO);
        CHKERR mField.modify_finite_element_add_field_row("BC_FE",
                                                          "SPATIAL_POSITION");
        CHKERR mField.modify_finite_element_add_field_col("BC_FE",
                                                          "SPATIAL_POSITION");
        CHKERR mField.modify_finite_element_add_field_data("BC_FE",
                                                           "SPATIAL_POSITION");
        CHKERR mField.modify_finite_element_add_field_data(
            "BC_FE", "MESH_NODE_POSITIONS");
        tris = intersect(lower_dim_ents, tris);
        CHKERR mField.add_ents_to_finite_element_by_type(tris, MBTRI, "BC_FE");
        CHKERR mField.build_finite_elements("BC_FE", &tris);
      }
    }
  }
#endif //__ANALITICAL_DISPLACEMENT__
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::declarePressureAleFE(const BitRefLevel bit,
                                                      const BitRefLevel mask,
                                                      const bool proc_only) {
 
  moab::Interface &moab = mField.get_moab();
  MoFEMFunctionBegin;
 
  Range tets_level;
  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
      bit, mask, MBTET, tets_level);
  if (proc_only) {
    tets_level = intersect(bitProcEnts, tets_level);
  }
 
  Range dofs_ents;
  dofs_ents.merge(tets_level);
  for (int dd = 1; dd != 3; dd++) {
    CHKERR moab.get_adjacencies(tets_level, dd, false, dofs_ents,
                                moab::Interface::UNION);
  }
  Range dofs_nodes;
  rval = moab.get_connectivity(tets_level, dofs_nodes, true);
  CHKERRQ_MOAB(rval);
  Range lower_dim_ents = unite(dofs_ents, dofs_nodes);
 
  CHKERR mField.add_finite_element("PRESSURE_ALE", MF_ZERO);
  CHKERR mField.modify_finite_element_add_field_row("PRESSURE_ALE",
                                                    "SPATIAL_POSITION");
  CHKERR mField.modify_finite_element_add_field_col("PRESSURE_ALE",
                                                    "SPATIAL_POSITION");
  CHKERR mField.modify_finite_element_add_field_data("PRESSURE_ALE",
                                                     "SPATIAL_POSITION");
  CHKERR mField.modify_finite_element_add_field_row("PRESSURE_ALE",
                                                    "MESH_NODE_POSITIONS");
  CHKERR mField.modify_finite_element_add_field_col("PRESSURE_ALE",
                                                    "MESH_NODE_POSITIONS");
  CHKERR mField.modify_finite_element_add_field_data("PRESSURE_ALE",
                                                     "MESH_NODE_POSITIONS");
 
  Range ents_to_add;
  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(mField, SIDESET | PRESSURESET,
                                                  it)) {
    Range tris;
    CHKERR mField.get_moab().get_entities_by_type(it->meshset, MBTRI, tris,
                                                  true);
    tris = intersect(lower_dim_ents, tris);
    ents_to_add.merge(tris);
  }
 
  Range current_ents_with_fe;
  CHKERR mField.get_finite_element_entities_by_handle("PRESSURE_ALE",
                                                      current_ents_with_fe);
  Range ents_to_remove;
  ents_to_remove = subtract(current_ents_with_fe, ents_to_add);
  CHKERR mField.remove_ents_from_finite_element("PRESSURE_ALE", ents_to_remove);
  CHKERR mField.add_ents_to_finite_element_by_type(ents_to_add, MBTRI,
                                                   "PRESSURE_ALE");
 
  CHKERR mField.build_finite_elements("PRESSURE_ALE", &lower_dim_ents);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::declareSurfaceForceAleFE(
    const BitRefLevel bit, const BitRefLevel mask, const bool proc_only) {
 
  moab::Interface &moab = mField.get_moab();
  MoFEMFunctionBegin;
 
  Range tets_level;
  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
      bit, mask, MBTET, tets_level);
  if (proc_only) {
    tets_level = intersect(bitProcEnts, tets_level);
  }
 
  Range dofs_ents;
  dofs_ents.merge(tets_level);
  for (int dd = 1; dd != 3; dd++) {
    CHKERR moab.get_adjacencies(tets_level, dd, false, dofs_ents,
                                moab::Interface::UNION);
  }
  Range dofs_nodes;
  rval = moab.get_connectivity(tets_level, dofs_nodes, true);
  CHKERRQ_MOAB(rval);
  Range lower_dim_ents = unite(dofs_ents, dofs_nodes);
 
  CHKERR mField.add_finite_element("FORCE_FE_ALE", MF_ZERO);
  CHKERR mField.modify_finite_element_add_field_row("FORCE_FE_ALE",
                                                    "SPATIAL_POSITION");
  CHKERR mField.modify_finite_element_add_field_col("FORCE_FE_ALE",
                                                    "SPATIAL_POSITION");
  CHKERR mField.modify_finite_element_add_field_data("FORCE_FE_ALE",
                                                     "SPATIAL_POSITION");
  CHKERR mField.modify_finite_element_add_field_row("FORCE_FE_ALE",
                                                    "MESH_NODE_POSITIONS");
  CHKERR mField.modify_finite_element_add_field_col("FORCE_FE_ALE",
                                                    "MESH_NODE_POSITIONS");
  CHKERR mField.modify_finite_element_add_field_data("FORCE_FE_ALE",
                                                     "MESH_NODE_POSITIONS");
 
  Range ents_to_add;
  for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(mField, NODESET | FORCESET,
                                                  it)) {
    Range tris;
    CHKERR mField.get_moab().get_entities_by_type(it->meshset, MBTRI, tris,
                                                  true);
    tris = intersect(lower_dim_ents, tris);
    ents_to_add.merge(tris);
  }
 
  const string block_set_force_name("FORCE");
  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
    if (it->getName().compare(0, block_set_force_name.length(),
                              block_set_force_name) == 0) {
      Range tris;
      CHKERR mField.get_moab().get_entities_by_type(it->meshset, MBTRI, tris,
                                                    true);
      tris = intersect(lower_dim_ents, tris);
      ents_to_add.merge(tris);
    }
  }
 
  Range current_ents_with_fe;
  CHKERR mField.get_finite_element_entities_by_handle("FORCE_FE_ALE",
                                                      current_ents_with_fe);
  Range ents_to_remove;
  ents_to_remove = subtract(current_ents_with_fe, ents_to_add);
  CHKERR mField.remove_ents_from_finite_element("FORCE_FE_ALE", ents_to_remove);
  CHKERR mField.add_ents_to_finite_element_by_type(ents_to_add, MBTRI,
                                                   "FORCE_FE_ALE");
 
  CHKERR mField.build_finite_elements("FORCE_FE_ALE", &lower_dim_ents);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::declareSpringsAleFE(const BitRefLevel bit,
                                                     const BitRefLevel mask,
                                                     const bool proc_only) {
 
  moab::Interface &moab = mField.get_moab();
  MoFEMFunctionBegin;
 
  Range tets_level;
  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
      bit, mask, MBTET, tets_level);
  if (proc_only) {
    tets_level = intersect(bitProcEnts, tets_level);
  }
 
  Range dofs_ents;
  dofs_ents.merge(tets_level);
  for (int dd = 1; dd != 3; dd++) {
    CHKERR moab.get_adjacencies(tets_level, dd, false, dofs_ents,
                                moab::Interface::UNION);
  }
  Range dofs_nodes;
  rval = moab.get_connectivity(tets_level, dofs_nodes, true);
  CHKERRQ_MOAB(rval);
  Range lower_dim_ents = unite(dofs_ents, dofs_nodes);
 
  CHKERR mField.add_finite_element("SPRING_ALE", MF_ZERO);
  CHKERR mField.modify_finite_element_add_field_row("SPRING_ALE",
                                                    "SPATIAL_POSITION");
  CHKERR mField.modify_finite_element_add_field_col("SPRING_ALE",
                                                    "SPATIAL_POSITION");
  CHKERR mField.modify_finite_element_add_field_data("SPRING_ALE",
                                                     "SPATIAL_POSITION");
  CHKERR mField.modify_finite_element_add_field_row("SPRING_ALE",
                                                    "MESH_NODE_POSITIONS");
  CHKERR mField.modify_finite_element_add_field_col("SPRING_ALE",
                                                    "MESH_NODE_POSITIONS");
  CHKERR mField.modify_finite_element_add_field_data("SPRING_ALE",
                                                     "MESH_NODE_POSITIONS");
 
  Range ents_to_add;
  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, bit)) {
    if (bit->getName().compare(0, 9, "SPRING_BC") == 0) {
      Range tris;
      CHKERR mField.get_moab().get_entities_by_type(bit->meshset, MBTRI, tris,
                                                    true);
      tris = intersect(lower_dim_ents, tris);
      ents_to_add.merge(tris);
    }
  }
 
  Range current_ents_with_fe;
  CHKERR mField.get_finite_element_entities_by_handle("SPRING_ALE",
                                                      current_ents_with_fe);
  Range ents_to_remove;
  ents_to_remove = subtract(current_ents_with_fe, ents_to_add);
  CHKERR mField.remove_ents_from_finite_element("SPRING_ALE", ents_to_remove);
  CHKERR mField.add_ents_to_finite_element_by_type(ents_to_add, MBTRI,
                                                   "SPRING_ALE");
 
  CHKERR mField.build_finite_elements("SPRING_ALE", &lower_dim_ents);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::declareSimpleContactAleFE(
    const BitRefLevel bit, const BitRefLevel mask, const bool proc_only) {
 
  moab::Interface &moab = mField.get_moab();
  MoFEMFunctionBegin;
 
  Range prisms_level;
  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
      bit, mask, MBPRISM, prisms_level);
 
  if (proc_only) {
    prisms_level = intersect(bitProcEnts, prisms_level);
  }
 
  Range contact_prisms = intersect(prisms_level, contactElements);
 
  CHKERR contactProblem->addContactElementALE(
      "SIMPLE_CONTACT_ALE", "SPATIAL_POSITION", "LAMBDA_CONTACT",
      "MESH_NODE_POSITIONS", contact_prisms, solveEigenStressProblem,
      "EIGEN_SPATIAL_POSITIONS");
 
  CHKERR mField.build_finite_elements("SIMPLE_CONTACT_ALE", &contact_prisms);
 
  Range face_on_prism;
  CHKERR moab.get_adjacencies(contact_prisms, 2, false, face_on_prism,
                              moab::Interface::UNION);
 
  Range tris_on_prism = face_on_prism.subset_by_type(MBTRI);
 
  Range check_tris =
      intersect(tris_on_prism, unite(contactSlaveFaces, contactMasterFaces));
 
  Range adj_vols_to_prisms;
  CHKERR moab.get_adjacencies(check_tris, 3, false, adj_vols_to_prisms,
                              moab::Interface::UNION);
 
  Range tet_level;
  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
      bit, mask, MBTET, tet_level);
 
  Range contact_tets_from_vols = adj_vols_to_prisms.subset_by_type(MBTET);
 
  Range contact_tets = intersect(tet_level, contact_tets_from_vols);
 
  CHKERR mField.add_finite_element("MAT_CONTACT", MF_ZERO);
  CHKERR mField.modify_finite_element_add_field_row("MAT_CONTACT",
                                                    "SPATIAL_POSITION");
  CHKERR mField.modify_finite_element_add_field_col("MAT_CONTACT",
                                                    "SPATIAL_POSITION");
  CHKERR mField.modify_finite_element_add_field_row("MAT_CONTACT",
                                                    "MESH_NODE_POSITIONS");
  CHKERR mField.modify_finite_element_add_field_col("MAT_CONTACT",
                                                    "MESH_NODE_POSITIONS");
  CHKERR mField.modify_finite_element_add_field_data("MAT_CONTACT",
                                                     "SPATIAL_POSITION");
  CHKERR mField.modify_finite_element_add_field_data("MAT_CONTACT",
                                                     "MESH_NODE_POSITIONS");
 
  Range current_ale_tets;
  CHKERR mField.get_finite_element_entities_by_handle("MAT_CONTACT",
                                                      current_ale_tets);
  Range ale_tets_to_remove;
  ale_tets_to_remove = subtract(current_ale_tets, contact_tets);
  CHKERR mField.remove_ents_from_finite_element("MAT_CONTACT",
                                                ale_tets_to_remove);
 
  CHKERR mField.add_ents_to_finite_element_by_type(contact_tets, MBTET,
                                                   "MAT_CONTACT");
  CHKERR mField.build_finite_elements("MAT_CONTACT", &contact_tets);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::declareMaterialFE(const BitRefLevel bit,
                                                   const BitRefLevel mask,
                                                   const bool proc_only,
                                                   const bool verb) {
  MoFEMFunctionBegin;
 
  Range tets_level;
  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
      bit, mask, MBTET, tets_level);
  if (proc_only) {
    tets_level = intersect(bitProcEnts, tets_level);
  }
 
  // Add finite elements
  CHKERR mField.add_finite_element("MATERIAL", MF_ZERO, verb);
  CHKERR mField.modify_finite_element_add_field_row("MATERIAL",
                                                    "SPATIAL_POSITION");
  CHKERR mField.modify_finite_element_add_field_col("MATERIAL",
                                                    "SPATIAL_POSITION");
  CHKERR mField.modify_finite_element_add_field_row("MATERIAL",
                                                    "MESH_NODE_POSITIONS");
  CHKERR mField.modify_finite_element_add_field_col("MATERIAL",
                                                    "MESH_NODE_POSITIONS");
  CHKERR mField.modify_finite_element_add_field_data("MATERIAL",
                                                     "SPATIAL_POSITION");
  CHKERR mField.modify_finite_element_add_field_data("MATERIAL",
                                                     "MESH_NODE_POSITIONS");
 
  {
    Range current_ents_with_fe;
    CHKERR mField.get_finite_element_entities_by_handle("MATERIAL",
                                                        current_ents_with_fe);
    Range ents_to_remove;
    ents_to_remove = subtract(current_ents_with_fe, tets_level);
    CHKERR mField.remove_ents_from_finite_element("MATERIAL", ents_to_remove);
    CHKERR mField.add_ents_to_finite_element_by_type(tets_level, MBTET,
                                                     "MATERIAL");
    // Range ents_to_build = subtract(tets_level, current_ents_with_fe);
    // CHKERR mField.build_finite_elements("MATERIAL", &ents_to_build);
    CHKERR mField.build_finite_elements("MATERIAL", &tets_level, verb);
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::declareFrontFE(const BitRefLevel bit,
                                                const BitRefLevel mask,
                                                const bool proc_only,
                                                const bool verb) {
  moab::Interface &moab = mField.get_moab();
  MoFEMFunctionBegin;
 
  CHKERR mField.add_finite_element("GRIFFITH_FORCE_ELEMENT", MF_ZERO, verb);
  CHKERR mField.modify_finite_element_add_field_row("GRIFFITH_FORCE_ELEMENT",
                                                    "MESH_NODE_POSITIONS");
  CHKERR mField.modify_finite_element_add_field_col("GRIFFITH_FORCE_ELEMENT",
                                                    "MESH_NODE_POSITIONS");
  CHKERR mField.modify_finite_element_add_field_data("GRIFFITH_FORCE_ELEMENT",
                                                     "MESH_NODE_POSITIONS");
  Range crack_front_adj_tris;
  CHKERR mField.get_moab().get_adjacencies(
      crackFrontNodes, 2, false, crack_front_adj_tris, moab::Interface::UNION);
  crack_front_adj_tris = intersect(crack_front_adj_tris, crackFaces);
  Range tris_level;
  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
      bit, mask, MBTRI, tris_level);
  crack_front_adj_tris = intersect(crack_front_adj_tris, tris_level);
  crackFrontNodesTris = crack_front_adj_tris;
  Range ents_to_remove;
  CHKERR mField.get_finite_element_entities_by_handle("GRIFFITH_FORCE_ELEMENT",
                                                      ents_to_remove);
  ents_to_remove = subtract(ents_to_remove, crack_front_adj_tris);
  CHKERR mField.remove_ents_from_finite_element("GRIFFITH_FORCE_ELEMENT",
                                                ents_to_remove, verb);
 
  CHKERR mField.add_ents_to_finite_element_by_type(crack_front_adj_tris, MBTRI,
                                                   "GRIFFITH_FORCE_ELEMENT");
  CHKERR mField.build_finite_elements("GRIFFITH_FORCE_ELEMENT",
                                      &crack_front_adj_tris, verb);
 
#ifdef __ANALITICAL_TRACTION__
  {
    Range tets_level;
    CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
        bit, mask, MBTET, tets_level);
    if (proc_only) {
      tets_level = intersect(bitProcEnts, tets_level);
    }
 
    Range dofs_ents;
    dofs_ents.merge(tets_level);
    for (int dd = 1; dd != 3; dd++) {
      CHKERR moab.get_adjacencies(tets_level, dd, false, dofs_ents,
                                  moab::Interface::UNION);
    }
    Range dofs_nodes;
    CHKERR moab.get_connectivity(tets_level, dofs_nodes, true);
    Range lower_dim_ents = unite(dofs_ents, dofs_nodes);
 
    CHKERR mField.add_finite_element("ANALITICAL_METERIAL_TRACTION", MF_ZERO);
    CHKERR mField.modify_finite_element_add_field_row(
        "ANALITICAL_METERIAL_TRACTION", "MESH_NODE_POSITIONS");
    CHKERR mField.modify_finite_element_add_field_col(
        "ANALITICAL_METERIAL_TRACTION", "MESH_NODE_POSITIONS");
    CHKERR mField.modify_finite_element_add_field_data(
        "ANALITICAL_METERIAL_TRACTION", "SPATIAL_POSITION");
    CHKERR mField.modify_finite_element_add_field_data(
        "ANALITICAL_METERIAL_TRACTION", "MESH_NODE_POSITIONS");
    MeshsetsManager *meshset_manager_ptr;
    ierr = mField.getInterface(meshset_manager_ptr);
    CHKERRQ(ierr);
    if (meshset_manager_ptr->checkMeshset("ANALITICAL_TRACTION")) {
      const CubitMeshSets *cubit_meshset_ptr;
      meshset_manager_ptr->getCubitMeshsetPtr("ANALITICAL_TRACTION",
                                              &cubit_meshset_ptr);
      Range tris;
      CHKERR meshset_manager_ptr->getEntitiesByDimension(
          cubit_meshset_ptr->getMeshsetId(), BLOCKSET, 2, tris, true);
      tris = intersect(tris, lower_dim_ents);
      CHKERR mField.add_ents_to_finite_element_by_type(
          tris, MBTRI, "ANALITICAL_METERIAL_TRACTION");
    }
    CHKERR mField.build_finite_elements("ANALITICAL_METERIAL_TRACTION",
                                        &lower_dim_ents);
  }
#endif //__ANALITICAL_TRACTION__
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::declareBothSidesFE(
    const BitRefLevel bit_spatial, const BitRefLevel bit_material,
    const BitRefLevel mask, const bool proc_only, const bool verb) {
  MoFEMFunctionBegin;
 
  auto get_prisms_level = [&](auto bit) {
    Range prisms_level;
    CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
        bit, mask, MBPRISM, prisms_level);
    if (proc_only) {
      prisms_level = intersect(bitProcEnts, prisms_level);
    }
    return prisms_level;
  };
 
  auto get_crack_prisms = [&](auto prisms_level) {
    Range crack_prisms = subtract(prisms_level, contactElements);
    crack_prisms = subtract(crack_prisms, mortarContactElements);
    return crack_prisms;
  };
 
  auto add_both_sides_fe = [&](const std::string fe_name,
                               const std::string lambda_field_name,
                               const std::string spatial_field, Range ents) {
    MoFEMFunctionBegin;
 
    CHKERR mField.add_finite_element(fe_name, MF_ZERO);
 
    if (!ents.empty()) {
      CHKERR mField.modify_finite_element_add_field_row(fe_name,
                                                        lambda_field_name);
      CHKERR mField.modify_finite_element_add_field_row(fe_name, spatial_field);
 
      CHKERR mField.modify_finite_element_add_field_col(fe_name,
                                                        lambda_field_name);
      CHKERR mField.modify_finite_element_add_field_col(fe_name, spatial_field);
 
      CHKERR mField.modify_finite_element_add_field_data(fe_name,
                                                         lambda_field_name);
      CHKERR mField.modify_finite_element_add_field_data(fe_name,
                                                         spatial_field);
 
      CHKERR mField.add_ents_to_finite_element_by_type(ents, MBPRISM, fe_name);
    }
 
    CHKERR mField.build_finite_elements(fe_name, &ents);
 
    MoFEMFunctionReturn(0);
  };
 
  if (solveEigenStressProblem) {
    CHKERR add_both_sides_fe("CLOSE_CRACK", "LAMBDA_CLOSE_CRACK",
                             "SPATIAL_POSITION",
                             get_crack_prisms(get_prisms_level(bit_spatial)));
  }
  CHKERR add_both_sides_fe("BOTH_SIDES_OF_CRACK", "LAMBDA_BOTH_SIDES",
                           "MESH_NODE_POSITIONS",
                           get_crack_prisms(get_prisms_level(bit_material)));
  {
    Range contact_prisms =
        intersect(get_prisms_level(bit_material), contactElements);
    CHKERR add_both_sides_fe("BOTH_SIDES_OF_CONTACT",
                             "LAMBDA_BOTH_SIDES_CONTACT", "MESH_NODE_POSITIONS",
                             contact_prisms);
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::declareSmoothingFE(const BitRefLevel bit,
                                                    const BitRefLevel mask,
                                                    const bool proc_only,
                                                    const bool verb) {
  MoFEMFunctionBegin;
  Range tets_level;
  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
      bit, mask, MBTET, tets_level);
  if (proc_only) {
    tets_level = intersect(bitProcEnts, tets_level);
  }
 
  // Add finite elements
  CHKERR mField.add_finite_element("SMOOTHING", MF_ZERO, verb);
  CHKERR mField.modify_finite_element_add_field_row("SMOOTHING",
                                                    "MESH_NODE_POSITIONS");
  CHKERR mField.modify_finite_element_add_field_col("SMOOTHING",
                                                    "MESH_NODE_POSITIONS");
  CHKERR mField.modify_finite_element_add_field_data("SMOOTHING",
                                                     "MESH_NODE_POSITIONS");
  CHKERR mField.modify_finite_element_add_field_data("SMOOTHING",
                                                     "SPATIAL_POSITION");
  CHKERR mField.modify_finite_element_add_field_row(
      "SMOOTHING", "LAMBDA_CRACKFRONT_AREA_TANGENT");
  CHKERR mField.modify_finite_element_add_field_col(
      "SMOOTHING", "LAMBDA_CRACKFRONT_AREA_TANGENT");
 
  // remove old ents no longer used as a this element
  Range ents_to_remove;
  CHKERR mField.get_finite_element_entities_by_handle("SMOOTHING",
                                                      ents_to_remove);
  ents_to_remove = subtract(ents_to_remove, tets_level);
  CHKERR mField.remove_ents_from_finite_element("SMOOTHING", ents_to_remove);
 
  CHKERR mField.add_ents_to_finite_element_by_type(tets_level, MBTET,
                                                   "SMOOTHING");
  CHKERR mField.build_finite_elements("SMOOTHING", &tets_level, verb);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::declareSurfaceFE(
    std::string fe_name, const BitRefLevel bit, const BitRefLevel mask,
    const std::vector<int> &ids, const bool proc_only, const int verb,
    const bool debug) {
  moab::Interface &moab = mField.get_moab();
  MeshsetsManager *meshset_manager_ptr;
  MoFEMFunctionBegin;
  CHKERR mField.getInterface(meshset_manager_ptr);
 
  Range level_tris;
  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
      bit, BitRefLevel().set(), MBTRI, level_tris);
  Range surface_ents = intersect(bodySkin, level_tris);
 
  if (proc_only) {
    // Note that some elements are on this proc, but not owned by that proc,
    // since
    // some tetrahedral share lower dimension entities
    Range proc_tris;
    CHKERR moab.get_adjacencies(bitProcEnts, 2, false, proc_tris,
                                moab::Interface::UNION);
    surface_ents = intersect(surface_ents, proc_tris);
  }
 
  CHKERR mField.add_finite_element(fe_name.c_str(), MF_ZERO, verb);
  CHKERR mField.modify_finite_element_add_field_row(fe_name,
                                                    "MESH_NODE_POSITIONS");
  CHKERR mField.modify_finite_element_add_field_col(fe_name,
                                                    "MESH_NODE_POSITIONS");
  CHKERR mField.modify_finite_element_add_field_data(fe_name,
                                                     "MESH_NODE_POSITIONS");
 
  for (std::vector<int>::const_iterator it = ids.begin(); it != ids.end();
       it++) {
    std::string field_name =
        "LAMBDA_SURFACE" + boost::lexical_cast<std::string>(*it);
    for (_IT_GET_ENT_FIELD_BY_NAME_FOR_LOOP_(mField, field_name, dit)) {
      if (dit->get()->getOwnerProc() == mField.get_comm_rank()) {
        EntityHandle ent = dit->get()->getEnt();
        Range adj_tris;
        CHKERR moab.get_adjacencies(&ent, 1, 2, false, adj_tris);
        adj_tris = intersect(adj_tris, bodySkin);
        surface_ents.merge(adj_tris);
      }
    }
    // Add finite elements
    CHKERR mField.modify_finite_element_add_field_row(fe_name,
                                                      field_name.c_str());
    CHKERR mField.modify_finite_element_add_field_col(fe_name,
                                                      field_name.c_str());
    CHKERR mField.modify_finite_element_add_field_data(fe_name,
                                                       field_name.c_str());
 
    // remove old ents no longer used as a this element
    Range ents_to_remove;
    CHKERR mField.get_finite_element_entities_by_handle(fe_name,
                                                        ents_to_remove);
    ents_to_remove = subtract(ents_to_remove, surface_ents);
    CHKERR mField.remove_ents_from_finite_element(fe_name, ents_to_remove);
    CHKERR mField.add_ents_to_finite_element_by_type(surface_ents, MBTRI,
                                                     fe_name);
  }
 
  CHKERR mField.build_finite_elements(fe_name, &surface_ents, verb);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
CrackPropagation::declareEdgeFE(std::string fe_name, const BitRefLevel bit,
                                const BitRefLevel mask, const bool proc_only,
                                const int verb, const bool debug) {
  auto meshset_manager_ptr = mField.getInterface<MeshsetsManager>();
  auto bit_ref_manager_ptr = mField.getInterface<BitRefManager>();
  MoFEMFunctionBegin;
 
  auto get_edges_block_ents = [this, meshset_manager_ptr]() {
    EntityHandle edge_block_set = getInterface<CPMeshCut>()->getEdgesBlockSet();
    Range meshset_ents;
    if (meshset_manager_ptr->checkMeshset(edge_block_set, BLOCKSET)) {
      CHKERR meshset_manager_ptr->getEntitiesByDimension(
          edge_block_set, BLOCKSET, 1, meshset_ents, true);
    }
    return meshset_ents;
  };
 
  auto get_skin = [this, bit_ref_manager_ptr]() {
    Range tets;
    CHKERR bit_ref_manager_ptr->getEntitiesByTypeAndRefLevel(
        mapBitLevel["mesh_cut"], BitRefLevel().set(), MBTET, tets);
    Skinner skin(&mField.get_moab());
    Range skin_faces;
    CHKERR skin.find_skin(0, tets, false, skin_faces);
    return skin_faces;
  };
 
  auto intersect_with_skin_edges = [this, get_edges_block_ents, get_skin]() {
    Range adj_edges;
    CHKERR mField.get_moab().get_adjacencies(get_skin(), 1, false, adj_edges,
                                             moab::Interface::UNION);
    return intersect(get_edges_block_ents(), adj_edges);
  };
  Range edges_ents = intersect_with_skin_edges();
 
  Range contact_faces;
  contact_faces.merge(contactSlaveFaces);
  contact_faces.merge(contactMasterFaces);
 
  contactOrientation =
      boost::shared_ptr<SurfaceSlidingConstrains::DriverElementOrientation>(
          new FaceOrientation(false, contact_faces, mapBitLevel["mesh_cut"]));
 
  CHKERR EdgeSlidingConstrains::CalculateEdgeBase::numberSurfaces(
      mField.get_moab(), edges_ents, get_skin());
  CHKERR EdgeSlidingConstrains::CalculateEdgeBase::setTags(
      mField.get_moab(), edges_ents, get_skin(), false, contactOrientation,
      &mField);
 
  if (debug) {
    Range faces = get_skin();
    CHKERR EdgeSlidingConstrains::CalculateEdgeBase::saveEdges(
        mField.get_moab(), "edges_normals.vtk", edges_ents, &faces);
  }
 
  auto intersect_with_bit_level_edges = [bit_ref_manager_ptr, &bit,
                                         &edges_ents]() {
    Range level_edges;
    CHKERR bit_ref_manager_ptr->getEntitiesByTypeAndRefLevel(
        bit, BitRefLevel().set(), MBEDGE, level_edges);
    return intersect(edges_ents, level_edges);
  };
  edges_ents = intersect_with_bit_level_edges();
 
  auto intersect_with_proc_edges = [this, &edges_ents, proc_only]() {
    if (proc_only) {
      Range proc_edges;
      CHKERR mField.get_moab().get_adjacencies(
          bitProcEnts, 1, false, proc_edges, moab::Interface::UNION);
      return intersect(edges_ents, proc_edges);
    }
    return edges_ents;
  };
  edges_ents = intersect_with_proc_edges();
  if (proc_only)
    CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(edges_ents,
                                                                     QUIET);
 
  CHKERR mField.add_finite_element(fe_name.c_str(), MF_ZERO, verb);
 
  auto add_field = [this, &fe_name](std::string &&field_name) {
    MoFEMFunctionBegin;
    CHKERR mField.modify_finite_element_add_field_row(fe_name, field_name);
    CHKERR mField.modify_finite_element_add_field_col(fe_name, field_name);
    CHKERR mField.modify_finite_element_add_field_data(fe_name, field_name);
    MoFEMFunctionReturn(0);
  };
  CHKERR add_field("MESH_NODE_POSITIONS");
  CHKERR add_field("LAMBDA_EDGE");
 
  if (debug) {
    CHKERR EdgeSlidingConstrains::CalculateEdgeBase::saveEdges(
        mField.get_moab(),
        "edges_normals_on_proc_" +
            boost::lexical_cast<std::string>(mField.get_comm_rank()) + ".vtk",
        edges_ents);
  }
 
  // remove old ents no longer used as a this element
  Range ents_to_remove;
  CHKERR mField.get_finite_element_entities_by_handle(fe_name, ents_to_remove);
  ents_to_remove = subtract(ents_to_remove, edges_ents);
  CHKERR mField.remove_ents_from_finite_element(fe_name, ents_to_remove);
  CHKERR mField.add_ents_to_finite_element_by_type(edges_ents, MBEDGE, fe_name);
 
  CHKERR mField.build_finite_elements(fe_name, &edges_ents, verb);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::declareCrackSurfaceFE(
    std::string fe_name, const BitRefLevel bit, const BitRefLevel mask,
    const bool proc_only, const int verb, const bool debug) {
  moab::Interface &moab = mField.get_moab();
  MeshsetsManager *meshset_manager_ptr;
  MoFEMFunctionBegin;
  CHKERR mField.getInterface(meshset_manager_ptr);
 
  Range level_tris;
  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
      bit, mask, MBTRI, level_tris);
  Range crack_surface;
  crack_surface = intersect(oneSideCrackFaces, level_tris);
  if (proc_only) {
    Range proc_tris;
    CHKERR moab.get_adjacencies(bitProcEnts, 2, false, proc_tris,
                                moab::Interface::UNION);
    crack_surface = intersect(crack_surface, proc_tris);
  }
  Range crack_front_adj_tris;
  CHKERR mField.get_moab().get_adjacencies(
      crackFrontNodes, 2, false, crack_front_adj_tris, moab::Interface::UNION);
  crack_front_adj_tris = intersect(crack_front_adj_tris, crackFaces);
  crack_front_adj_tris = intersect(crack_front_adj_tris, level_tris);
 
  std::string field_name =
      "LAMBDA_SURFACE" + boost::lexical_cast<std::string>(
                             getInterface<CPMeshCut>()->getCrackSurfaceId());
 
  // Add finite elements
  CHKERR mField.add_finite_element(fe_name.c_str(), MF_ZERO, verb);
  CHKERR mField.modify_finite_element_add_field_row(fe_name,
                                                    field_name.c_str());
  CHKERR mField.modify_finite_element_add_field_row(fe_name,
                                                    "MESH_NODE_POSITIONS");
  CHKERR mField.modify_finite_element_add_field_col(fe_name,
                                                    field_name.c_str());
  CHKERR mField.modify_finite_element_add_field_col(fe_name,
                                                    "MESH_NODE_POSITIONS");
  CHKERR mField.modify_finite_element_add_field_data(fe_name,
                                                     field_name.c_str());
  CHKERR mField.modify_finite_element_add_field_data(fe_name,
                                                     "MESH_NODE_POSITIONS");
  if (proc_only)
    CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
        crack_surface, QUIET);
 
  {
    // remove ents no longer part of this element
    Range ents;
    CHKERR mField.get_finite_element_entities_by_handle(fe_name, ents);
    ents = subtract(ents, crack_surface);
    if (!ents.empty()) {
      CHKERR mField.remove_ents_from_finite_element(fe_name, ents);
    }
  }
 
  CHKERR mField.add_ents_to_finite_element_by_type(crack_surface, MBTRI,
                                                   fe_name);
  CHKERR mField.build_finite_elements(fe_name, &crack_surface, verb);
 
  if (debug) {
    std::ostringstream file;
    file << "out_crack_faces_" << mField.get_comm_rank() << ".vtk";
    EntityHandle meshset;
    CHKERR moab.create_meshset(MESHSET_SET, meshset);
    CHKERR moab.add_entities(meshset, crack_surface);
    CHKERR moab.write_file(file.str().c_str(), "VTK", "", &meshset, 1);
    CHKERR moab.delete_entities(&meshset, 1);
  }
 
  // Add CRACKFRONT_AREA_ELEM
  CHKERR mField.add_finite_element("CRACKFRONT_AREA_ELEM", MF_ZERO, verb);
  CHKERR mField.modify_finite_element_add_field_row("CRACKFRONT_AREA_ELEM",
                                                    "LAMBDA_CRACKFRONT_AREA");
  CHKERR mField.modify_finite_element_add_field_row("CRACKFRONT_AREA_ELEM",
                                                    "MESH_NODE_POSITIONS");
  CHKERR mField.modify_finite_element_add_field_col("CRACKFRONT_AREA_ELEM",
                                                    "LAMBDA_CRACKFRONT_AREA");
  CHKERR mField.modify_finite_element_add_field_col("CRACKFRONT_AREA_ELEM",
                                                    "MESH_NODE_POSITIONS");
  CHKERR mField.modify_finite_element_add_field_data("CRACKFRONT_AREA_ELEM",
                                                     "MESH_NODE_POSITIONS");
  CHKERR mField.modify_finite_element_add_field_data("CRACKFRONT_AREA_ELEM",
                                                     "LAMBDA_CRACKFRONT_AREA");
  // Add CRACKFRONT_AREA_TANGENT_ELEM
  CHKERR mField.add_finite_element("CRACKFRONT_AREA_TANGENT_ELEM", MF_ZERO,
                                   verb);
  CHKERR mField.modify_finite_element_add_field_row(
      "CRACKFRONT_AREA_TANGENT_ELEM", "LAMBDA_CRACKFRONT_AREA_TANGENT");
  CHKERR mField.modify_finite_element_add_field_row(
      "CRACKFRONT_AREA_TANGENT_ELEM", "MESH_NODE_POSITIONS");
  CHKERR mField.modify_finite_element_add_field_col(
      "CRACKFRONT_AREA_TANGENT_ELEM", "LAMBDA_CRACKFRONT_AREA_TANGENT");
  CHKERR mField.modify_finite_element_add_field_col(
      "CRACKFRONT_AREA_TANGENT_ELEM", "MESH_NODE_POSITIONS");
  CHKERR mField.modify_finite_element_add_field_data(
      "CRACKFRONT_AREA_TANGENT_ELEM", "MESH_NODE_POSITIONS");
  CHKERR mField.modify_finite_element_add_field_data(
      "CRACKFRONT_AREA_TANGENT_ELEM", "LAMBDA_CRACKFRONT_AREA_TANGENT");
 
  if (proc_only)
    CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
        crack_front_adj_tris, QUIET);
 
  auto remeve_ents_from_fe =
      [this, &crack_front_adj_tris](const std::string fe_name) {
        MoFEMFunctionBegin;
        Range ents;
        CHKERR mField.get_finite_element_entities_by_handle(fe_name, ents);
        ents = subtract(ents, crack_front_adj_tris);
        if (!ents.empty()) {
          CHKERR mField.remove_ents_from_finite_element(fe_name, ents);
        }
        MoFEMFunctionReturn(0);
      };
  CHKERR remeve_ents_from_fe("CRACKFRONT_AREA_ELEM");
  CHKERR remeve_ents_from_fe("CRACKFRONT_AREA_TANGENT_ELEM");
 
  CHKERR mField.add_ents_to_finite_element_by_type(crack_front_adj_tris, MBTRI,
                                                   "CRACKFRONT_AREA_ELEM");
  CHKERR mField.add_ents_to_finite_element_by_type(
      crack_front_adj_tris, MBTRI, "CRACKFRONT_AREA_TANGENT_ELEM");
 
  if (debug) {
    std::ostringstream file;
    file << "out_crack_front_adj_tris_" << mField.get_comm_rank() << ".vtk";
    EntityHandle meshset;
    CHKERR moab.create_meshset(MESHSET_SET, meshset);
    CHKERR moab.add_entities(meshset, crack_front_adj_tris);
    CHKERR moab.write_file(file.str().c_str(), "VTK", "", &meshset, 1);
    CHKERR moab.delete_entities(&meshset, 1);
  }
 
  CHKERR mField.build_finite_elements("CRACKFRONT_AREA_ELEM",
                                      &crack_front_adj_tris, verb);
  CHKERR mField.build_finite_elements("CRACKFRONT_AREA_TANGENT_ELEM",
                                      &crack_front_adj_tris, verb);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::createCrackPropagationDM(
    SmartPetscObj<DM> &dm, const std::string prb_name,
    SmartPetscObj<DM> dm_elastic, SmartPetscObj<DM> dm_material,
    const BitRefLevel bit, const BitRefLevel mask,
    const std::vector<std::string> fe_list) {
  const MoFEM::Problem *prb_elastic_ptr;
  const MoFEM::Problem *prb_material_ptr;
  MoFEMFunctionBegin;
  CHKERR DMMoFEMGetProblemPtr(dm_elastic, &prb_elastic_ptr);
  CHKERR DMMoFEMGetProblemPtr(dm_material, &prb_material_ptr);
  dm = createSmartDM(PETSC_COMM_WORLD, "MOFEM");
  CHKERR DMMoFEMCreateMoFEM(dm, &mField, prb_name.c_str(), bit, mask);
 
  CHKERR DMMoFEMAddElement(dm, "ELASTIC");
  CHKERR DMMoFEMAddElement(dm, "ELASTIC_NOT_ALE");
  CHKERR DMMoFEMAddElement(dm, "SPRING");
  CHKERR DMMoFEMAddElement(dm, "CONTACT");
  CHKERR DMMoFEMAddElement(dm, "MORTAR_CONTACT");
  CHKERR DMMoFEMAddElement(dm, "CONTACT_POST_PROC");
  CHKERR DMMoFEMAddElement(dm, "SKIN"); // for postprocessing
  CHKERR DMMoFEMAddElement(dm, "ARC_LENGTH");
  CHKERR DMMoFEMAddElement(dm, "FORCE_FE");
  CHKERR DMMoFEMAddElement(dm, "PRESSURE_FE");
 
  if (isSurfaceForceAle) {
    CHKERR DMMoFEMAddElement(dm, "FORCE_FE_ALE");
  }
  if (isPressureAle) {
    CHKERR DMMoFEMAddElement(dm, "PRESSURE_ALE");
  }
  if (areSpringsAle) {
    CHKERR DMMoFEMAddElement(dm, "SPRING_ALE");
  }
  if (!contactElements.empty() && !ignoreContact && !fixContactNodes) {
    CHKERR DMMoFEMAddElement(dm, "SIMPLE_CONTACT_ALE");
    CHKERR DMMoFEMAddElement(dm, "MAT_CONTACT");
  }
  CHKERR DMMoFEMAddElement(dm, "MATERIAL");
 
  CHKERR DMMoFEMAddElement(dm, "GRIFFITH_FORCE_ELEMENT");
  for (std::vector<std::string>::const_iterator it = fe_list.begin();
       it != fe_list.end(); it++)
    CHKERR DMMoFEMAddElement(dm, it->c_str());
 
  CHKERR DMMoFEMAddElement(dm, "CRACK_SURFACE");
  CHKERR DMMoFEMAddElement(dm, "CRACKFRONT_AREA_ELEM");
  CHKERR DMMoFEMAddElement(dm, "CRACKFRONT_AREA_TANGENT_ELEM");
#ifdef __ANALITICAL_TRACTION__
  CHKERR DMMoFEMAddElement(dm, "ANALITICAL_METERIAL_TRACTION");
#endif
  if (mField.check_finite_element("SMOOTHING"))
    CHKERR DMMoFEMAddElement(dm, "SMOOTHING");
 
  if (mField.check_finite_element("EDGE"))
    CHKERR DMMoFEMAddElement(dm, "EDGE");
 
  if (mField.check_finite_element("BOTH_SIDES_OF_CRACK"))
    CHKERR DMMoFEMAddElement(dm, "BOTH_SIDES_OF_CRACK");
 
  if (mField.check_finite_element("BOTH_SIDES_OF_CONTACT"))
    CHKERR DMMoFEMAddElement(dm, "BOTH_SIDES_OF_CONTACT");
 
  CHKERR DMMoFEMAddRowCompositeProblem(dm, prb_elastic_ptr->getName().c_str());
  CHKERR DMMoFEMAddRowCompositeProblem(dm, prb_material_ptr->getName().c_str());
  CHKERR DMMoFEMSetIsPartitioned(dm, PETSC_TRUE);
  CHKERR DMSetUp(dm);
 
  // Create material section. This is used by field split solver and snes
  // monitor
  PetscSection section;
  CHKERR mField.getInterface<ISManager>()->sectionCreate(prb_name, &section);
  CHKERR DMSetDefaultSection(dm, section);
  CHKERR DMSetDefaultGlobalSection(dm, section);
  // CHKERR PetscSectionView(section,PETSC_VIEWER_STDOUT_WORLD);
  CHKERR PetscSectionDestroy(&section);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::createElasticDM(SmartPetscObj<DM> &dm,
                                                 const std::string prb_name,
                                                 const BitRefLevel bit,
                                                 const BitRefLevel mask) {
  ProblemsManager *prb_mng_ptr;
  MoFEMFunctionBegin;
  CHKERR mField.getInterface(prb_mng_ptr);
  CHKERR mField.build_adjacencies(bit, QUIET);
  dm = createSmartDM(PETSC_COMM_WORLD, "MOFEM");
  CHKERR DMMoFEMCreateMoFEM(dm, &mField, prb_name.c_str(), bit, mask);
  CHKERR DMMoFEMAddElement(dm, "ELASTIC");
  CHKERR DMMoFEMAddElement(dm, "SPRING");
  CHKERR DMMoFEMAddElement(dm, "CONTACT");
  CHKERR DMMoFEMAddElement(dm, "MORTAR_CONTACT");
  CHKERR DMMoFEMAddElement(dm, "CONTACT_POST_PROC");
  CHKERR DMMoFEMAddElement(dm, "SKIN");
  CHKERR DMMoFEMAddElement(dm, "ARC_LENGTH");
  CHKERR DMMoFEMAddElement(dm, "FORCE_FE");
  CHKERR DMMoFEMAddElement(dm, "PRESSURE_FE");
#ifdef __ANALITICAL_TRACTION__
  CHKERR DMMoFEMAddElement(dm, "ANALITICAL_TRACTION");
#endif
 
  CHKERR DMMoFEMSetIsPartitioned(dm, PETSC_TRUE);
  mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
 
  CHKERR DMSetUp(dm);
 
  // Create material section. This is used by field split solver and snes
  // monitor
  PetscSection section;
  CHKERR mField.getInterface<ISManager>()->sectionCreate(prb_name, &section);
  CHKERR DMSetDefaultSection(dm, section);
  CHKERR DMSetDefaultGlobalSection(dm, section);
  // CHKERR PetscSectionView(section,PETSC_VIEWER_STDOUT_WORLD);
  CHKERR PetscSectionDestroy(&section);
 
  mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_FALSE;
 
  // Setting arc-length control
  const MoFEM::Problem *problem_ptr;
  CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
  arcCtx = boost::make_shared<ArcLengthCtx>(mField, problem_ptr->getName(),
                                            "LAMBDA_ARC_LENGTH");
  CHKERR arcCtx->setAlphaBeta(0, 1);
  auto arc_snes_ctx = boost::make_shared<ArcLengthSnesCtx>(
      mField, problem_ptr->getName(), arcCtx);
  CHKERR DMMoFEMSetSnesCtx(dm, arc_snes_ctx);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::createEigenElasticDM(
    SmartPetscObj<DM> &dm, const std::string prb_name, const BitRefLevel bit,
    const BitRefLevel mask) {
  ProblemsManager *prb_mng_ptr;
  MoFEMFunctionBegin;
  CHKERR mField.getInterface(prb_mng_ptr);
  CHKERR mField.build_adjacencies(bit, QUIET);
  dm = createSmartDM(PETSC_COMM_WORLD, "MOFEM");
  CHKERR DMMoFEMCreateMoFEM(dm, &mField, prb_name.c_str(), bit, mask);
  CHKERR DMMoFEMAddElement(dm, "ELASTIC");
  CHKERR DMMoFEMAddElement(dm, "SPRING");
  CHKERR DMMoFEMAddElement(dm, "CONTACT");
  CHKERR DMMoFEMAddElement(dm, "MORTAR_CONTACT");
  CHKERR DMMoFEMAddElement(dm, "CONTACT_POST_PROC");
  CHKERR DMMoFEMAddElement(dm, "SKIN");
  CHKERR DMMoFEMAddElement(dm, "ARC_LENGTH");
  CHKERR DMMoFEMAddElement(dm, "FORCE_FE");
  CHKERR DMMoFEMAddElement(dm, "PRESSURE_FE");
  CHKERR DMMoFEMAddElement(dm, "CLOSE_CRACK");
#ifdef __ANALITICAL_TRACTION__
  CHKERR DMMoFEMAddElement(dm, "ANALITICAL_TRACTION");
#endif
 
  CHKERR DMMoFEMSetIsPartitioned(dm, PETSC_TRUE);
  mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
 
  CHKERR DMSetUp(dm);
 
  auto remove_higher_dofs_than_approx_dofs_from_eigen_elastic = [&]() {
    MoFEMFunctionBegin;
    Range ents;
    CHKERR mField.get_moab().get_entities_by_handle(0, ents, true);
    auto prb_mng = mField.getInterface<ProblemsManager>();
    CHKERR prb_mng->removeDofsOnEntities("EIGEN_ELASTIC", "SPATIAL_POSITION",
                                         ents, 0, 3, approxOrder + 1, 100);
    MoFEMFunctionReturn(0);
  };
  CHKERR remove_higher_dofs_than_approx_dofs_from_eigen_elastic();
 
  // Create material section. This is used by field split solver and snes
  // monitor
  PetscSection section;
  CHKERR mField.getInterface<ISManager>()->sectionCreate(prb_name, &section);
  CHKERR DMSetDefaultSection(dm, section);
  CHKERR DMSetDefaultGlobalSection(dm, section);
  // CHKERR PetscSectionView(section,PETSC_VIEWER_STDOUT_WORLD);
  CHKERR PetscSectionDestroy(&section);
 
  mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_FALSE;
 
  // Setting arc-length control
  const MoFEM::Problem *problem_ptr;
  CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
  arcEigenCtx = boost::make_shared<ArcLengthCtx>(mField, problem_ptr->getName(),
                                                 "LAMBDA_ARC_LENGTH");
  CHKERR arcEigenCtx->setAlphaBeta(0, 1);
  auto arc_snes_ctx = boost::make_shared<ArcLengthSnesCtx>(
      mField, problem_ptr->getName(), arcEigenCtx);
  CHKERR DMMoFEMSetSnesCtx(dm, arc_snes_ctx);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::createBcDM(SmartPetscObj<DM> &dm,
                                            const std::string prb_name,
                                            const BitRefLevel bit,
                                            const BitRefLevel mask) {
  ProblemsManager *prb_mng_ptr;
  MeshsetsManager *meshset_manager_ptr;
  MoFEMFunctionBegin;
  CHKERR mField.getInterface(prb_mng_ptr);
  CHKERR mField.getInterface(meshset_manager_ptr);
  CHKERR mField.build_adjacencies(bit, QUIET);
#ifdef __ANALITICAL_DISPLACEMENT__
  if (meshset_manager_ptr->checkMeshset("ANALITICAL_DISP")) {
    const CubitMeshSets *cubit_meshset_ptr;
    meshset_manager_ptr->getCubitMeshsetPtr("ANALITICAL_DISP",
                                            &cubit_meshset_ptr);
    Range tris;
    CHKERR meshset_manager_ptr->getEntitiesByDimension(
        cubit_meshset_ptr->getMeshsetId(), BLOCKSET, 2, tris, true);
    if (!tris.empty()) {
      dm = createSmartDM(PETSC_COMM_WORLD, "MOFEM");
      CHKERR DMMoFEMCreateMoFEM(dm, &mField, prb_name.c_str(), bit, mask);
      CHKERR DMMoFEMAddElement(dm, "BC_FE");
      CHKERR DMMoFEMSetIsPartitioned(dm, PETSC_TRUE);
      mField.getInterface<ProblemsManager>()->synchroniseProblemEntities =
          PETSC_TRUE;
      CHKERR DMSetUp(dm);
      mField.getInterface<ProblemsManager>()->synchroniseProblemEntities =
          PETSC_FALSE;
    }
  }
#endif //__ANALITICAL_DISPLACEMENT__
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::createMaterialDM(
    SmartPetscObj<DM> &dm, const std::string prb_name, const BitRefLevel bit,
    const BitRefLevel mask, const std::vector<std::string> fe_list,
    const bool debug) {
  MoFEMFunctionBegin;
  CHKERR mField.build_adjacencies(bit, QUIET);
  dm = createSmartDM(PETSC_COMM_WORLD, "MOFEM");
  CHKERR DMMoFEMCreateMoFEM(dm, &mField, prb_name.c_str(), bit, mask);
  CHKERR DMMoFEMAddElement(dm, "GRIFFITH_FORCE_ELEMENT");
 
  for (auto &fe : fe_list)
    CHKERR DMMoFEMAddElement(dm, fe.c_str());
 
  CHKERR DMMoFEMAddElement(dm, "CRACK_SURFACE");
  CHKERR DMMoFEMAddElement(dm, "CRACKFRONT_AREA_ELEM");
  CHKERR DMMoFEMAddElement(dm, "CRACKFRONT_AREA_TANGENT_ELEM");
#ifdef __ANALITICAL_TRACTION__
  CHKERR DMMoFEMAddElement(dm, "ANALITICAL_METERIAL_TRACTION");
#endif
  if (mField.check_finite_element("SMOOTHING"))
    CHKERR DMMoFEMAddElement(dm, "SMOOTHING");
 
  if (mField.check_finite_element("EDGE"))
    CHKERR DMMoFEMAddElement(dm, "EDGE");
 
  if (mField.check_finite_element("BOTH_SIDES_OF_CRACK"))
    CHKERR DMMoFEMAddElement(dm, "BOTH_SIDES_OF_CRACK");
 
  if (mField.check_finite_element("BOTH_SIDES_OF_CONTACT"))
    CHKERR DMMoFEMAddElement(dm, "BOTH_SIDES_OF_CONTACT");
 
  CHKERR DMMoFEMSetIsPartitioned(dm, PETSC_TRUE);
  CHKERR DMSetUp(dm);
 
  // Create crack propagation section
  PetscSection section;
  CHKERR mField.getInterface<ISManager>()->sectionCreate(prb_name, &section);
  CHKERR DMSetDefaultSection(dm, section);
  CHKERR DMSetDefaultGlobalSection(dm, section);
  // CHKERR PetscSectionView(section,PETSC_VIEWER_STDOUT_WORLD);
  CHKERR PetscSectionDestroy(&section);
 
  if (debug) {
    std::vector<std::string> fe_lists;
    const MoFEM::Problem *problem_ptr;
    CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
    const FiniteElement_multiIndex *fe_ptr;
    CHKERR mField.get_finite_elements(&fe_ptr);
    for (FiniteElement_multiIndex::iterator fit = fe_ptr->begin();
         fit != fe_ptr->end(); fit++) {
      if ((fit->get()->getId() & problem_ptr->getBitFEId()).any()) {
        std::string fe_name = fit->get()->getName();
        EntityHandle meshset;
        CHKERR mField.getInterface<ProblemsManager>()->getFEMeshset(
            prb_name, fe_name, &meshset);
        CHKERR mField.get_moab().write_file(
            (prb_name + "_" + fe_name + ".h5m").c_str(), "MOAB",
            "PARALLEL=WRITE_PART", &meshset, 1);
        CHKERR mField.get_moab().delete_entities(&meshset, 1);
      }
    }
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::createMaterialForcesDM(
    SmartPetscObj<DM> &dm, SmartPetscObj<DM> dm_material,
    const std::string prb_name, const int verb) {
  MoFEMFunctionBegin;
  dm = createSmartDM(PETSC_COMM_WORLD, "MOFEM");
  CHKERR DMMoFEMSetVerbosity(dm, verb);
  CHKERR DMMoFEMCreateSubDM(dm, dm_material, prb_name.c_str());
  CHKERR DMMoFEMAddSubFieldRow(dm, "MESH_NODE_POSITIONS");
  CHKERR DMMoFEMAddSubFieldCol(dm, "MESH_NODE_POSITIONS");
  CHKERR DMMoFEMAddElement(dm, "MATERIAL");
  if (mField.check_finite_element("SMOOTHING")) {
    CHKERR DMMoFEMAddElement(dm, "SMOOTHING");
  }
  CHKERR DMMoFEMAddElement(dm, "GRIFFITH_FORCE_ELEMENT");
#ifdef __ANALITICAL_TRACTION__
  CHKERR DMMoFEMAddElement(dm, "ANALITICAL_METERIAL_TRACTION");
#endif
  CHKERR DMMoFEMSetSquareProblem(dm, PETSC_TRUE);
  CHKERR DMMoFEMSetIsPartitioned(dm, PETSC_TRUE);
  CHKERR DMSetUp(dm);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::createSurfaceProjectionDM(
    SmartPetscObj<DM> &dm, SmartPetscObj<DM> dm_material,
    const std::string prb_name, const std::vector<int> surface_ids,
    std::vector<std::string> fe_list, const int verb) {
  MoFEMFunctionBegin;
  dm = createSmartDM(PETSC_COMM_WORLD, "MOFEM");
  CHKERR DMMoFEMSetVerbosity(dm, verb);
  CHKERR DMMoFEMCreateSubDM(dm, dm_material, prb_name.c_str());
  CHKERR DMMoFEMSetSquareProblem(dm, PETSC_FALSE);
  for (auto id : surface_ids) {
    CHKERR DMMoFEMAddSubFieldRow(
        dm, ("LAMBDA_SURFACE" + boost::lexical_cast<std::string>(id)).c_str());
  }
  CHKERR DMMoFEMAddSubFieldCol(dm, "MESH_NODE_POSITIONS");
  for (auto fe : fe_list) {
    CHKERR DMMoFEMAddElement(dm, fe.c_str());
  }
  auto get_edges_block_set = [this]() {
    return getInterface<CPMeshCut>()->getEdgesBlockSet();
  };
  if (get_edges_block_set()) {
    CHKERR DMMoFEMAddSubFieldRow(dm, "LAMBDA_EDGE");
    CHKERR DMMoFEMAddElement(dm, "EDGE");
  }
  CHKERR DMMoFEMSetIsPartitioned(dm, PETSC_TRUE);
  CHKERR DMSetUp(dm);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode CrackPropagation::createCrackFrontAreaDM(
    SmartPetscObj<DM> &dm, SmartPetscObj<DM> dm_material,
    const std::string prb_name, const bool verb) {
  MoFEMFunctionBegin;
  dm = createSmartDM(PETSC_COMM_WORLD, "MOFEM");
  CHKERR DMMoFEMSetVerbosity(dm, verb);
  CHKERR DMMoFEMCreateSubDM(dm, dm_material, prb_name.c_str());
  CHKERR DMMoFEMAddSubFieldRow(dm, "LAMBDA_CRACKFRONT_AREA");
  CHKERR DMMoFEMAddSubFieldRow(dm, "LAMBDA_CRACKFRONT_AREA_TANGENT");
  CHKERR DMMoFEMAddSubFieldCol(dm, "MESH_NODE_POSITIONS");
  CHKERR DMMoFEMAddElement(dm, "CRACKFRONT_AREA_ELEM");
  CHKERR DMMoFEMAddElement(dm, "CRACKFRONT_AREA_TANGENT_ELEM");
  CHKERR DMMoFEMSetSquareProblem(dm, PETSC_FALSE);
  CHKERR DMMoFEMSetIsPartitioned(dm, PETSC_TRUE);
  CHKERR DMSetUp(dm);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
CrackPropagation::assembleElasticDM(const std::string mwls_stress_tag_name,
                                    const int verb, const bool debug) {
  MoFEMFunctionBegin;
 
  // Create elastic element data structure
  elasticFe = boost::make_shared<NonlinearElasticElement>(mField, ELASTIC_TAG);
 
  // Create elastic element finite element instance for residuals. Note that
  // this element approx. singularity at crack front.
  feRhs = boost::make_shared<CrackFrontElement>(
      mField, setSingularCoordinates, crackFrontNodes, crackFrontNodesEdges,
      crackFrontElements, addSingularity);
 
  // Create finite element instance for assembly of tangent matrix
  feLhs = boost::make_shared<CrackFrontElement>(
      mField, setSingularCoordinates, crackFrontNodes, crackFrontNodesEdges,
      crackFrontElements, addSingularity);
 
  // Arbitrary lagrangian formulation, mesh node positions are taken into
  // account by operators.
  feRhs->meshPositionsFieldName = "NONE";
  feLhs->meshPositionsFieldName = "NONE";
  feRhs->addToRule = 0;
  feLhs->addToRule = 0;
 
  // Create operators to calculate finite element matrices for elastic element
  onlyHooke = true;
  using BlockData = NonlinearElasticElement::BlockData;
  boost::shared_ptr<map<int, BlockData>> block_sets_ptr(
      elasticFe, &elasticFe->setOfBlocks);
  {
    for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
             mField, BLOCKSET | MAT_ELASTICSET, it)) {
 
      // Get data from block
      Mat_Elastic mydata;
      CHKERR it->getAttributeDataStructure(mydata);
      int id = it->getMeshsetId();
      EntityHandle meshset = it->getMeshset();
      CHKERR mField.get_moab().get_entities_by_type(
          meshset, MBTET, elasticFe->setOfBlocks[id].tEts, true);
      elasticFe->setOfBlocks[id].iD = id;
      elasticFe->setOfBlocks[id].E = mydata.data.Young;
      elasticFe->setOfBlocks[id].PoissonRatio = mydata.data.Poisson;
 
      // Print material properties of blocks after being read
      CHKERR PetscPrintf(PETSC_COMM_WORLD, "\nMaterial block %d \n", id);
      CHKERR PetscPrintf(PETSC_COMM_WORLD, "\tYoung's modulus %6.4e \n",
                         mydata.data.Young);
      CHKERR PetscPrintf(PETSC_COMM_WORLD, "\tPoisson's ratio %6.4e \n\n",
                         mydata.data.Poisson);
 
      // Set default material to elastic blocks. Create material instances
      switch (defaultMaterial) {
      case HOOKE:
        elasticFe->setOfBlocks[id].materialDoublePtr =
            boost::make_shared<Hooke<double>>();
        elasticFe->setOfBlocks[id].materialAdoublePtr =
            boost::make_shared<Hooke<adouble>>();
        break;
      case KIRCHHOFF:
        elasticFe->setOfBlocks[id].materialDoublePtr = boost::make_shared<
            NonlinearElasticElement::FunctionsToCalculatePiolaKirchhoffI<
                double>>();
        elasticFe->setOfBlocks[id].materialAdoublePtr = boost::make_shared<
            NonlinearElasticElement::FunctionsToCalculatePiolaKirchhoffI<
                adouble>>();
        onlyHooke = false;
        break;
      case NEOHOOKEAN:
        elasticFe->setOfBlocks[id].materialDoublePtr =
            boost::make_shared<NeoHookean<double>>();
        elasticFe->setOfBlocks[id].materialAdoublePtr =
            boost::make_shared<NeoHookean<adouble>>();
        onlyHooke = false;
        break;
      case BONEHOOKE:
        elasticFe->setOfBlocks[id].materialDoublePtr =
            boost::make_shared<Hooke<double>>();
        elasticFe->setOfBlocks[id].materialAdoublePtr =
            boost::make_shared<Hooke<adouble>>();
        break;
      default:
        SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
                "Material type not yet implemented");
        break;
      }
    }
    if (onlyHookeFromOptions == PETSC_FALSE)
      onlyHooke = false;
 
    elasticFe->commonData.spatialPositions = "SPATIAL_POSITION";
    elasticFe->commonData.meshPositions = "MESH_NODE_POSITIONS";
 
    auto data_hooke_element_at_pts =
        boost::make_shared<HookeElement::DataAtIntegrationPts>();
 
    // calculate material positions at integration points
    auto mat_pos_at_pts_ptr = boost::make_shared<MatrixDouble>();
    feRhs->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
        "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr));
    boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr;
    if (residualStressBlock != -1 || densityMapBlock != -1) {
      mat_grad_pos_at_pts_ptr = boost::make_shared<MatrixDouble>();
      feRhs->getOpPtrVector().push_back(
          new OpCalculateVectorFieldGradient<3, 3>("MESH_NODE_POSITIONS",
                                                   mat_grad_pos_at_pts_ptr));
    }
 
    /// Rhs
    if (!onlyHooke) {
      // Calculate spatial positions and gradients (of deformation) at
      // integration pts. This operator takes into account singularity at crack
      // front
      feRhs->getOpPtrVector().push_back(new OpGetCrackFrontCommonDataAtGaussPts(
          "SPATIAL_POSITION", elasticFe->commonData, feRhs->singularElement,
          feRhs->invSJac));
      // Calculate material positions and gradients at integration pts.
      feRhs->getOpPtrVector().push_back(
          new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
              "MESH_NODE_POSITIONS", elasticFe->commonData));
      // Transfrom base functions to create singularity at crack front
      feRhs->getOpPtrVector().push_back(new OpTransfromSingularBaseFunctions(
          feRhs->singularElement, feRhs->detS, feRhs->invSJac));
      // Iterate over material blocks
      map<int, NonlinearElasticElement::BlockData>::iterator sit =
          elasticFe->setOfBlocks.begin();
      for (; sit != elasticFe->setOfBlocks.end(); sit++) {
        // Evaluate stress at integration pts.
        feRhs->getOpPtrVector().push_back(
            new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
                "SPATIAL_POSITION", sit->second, elasticFe->commonData,
                ELASTIC_TAG, false, true, false));
        // Assemble internal forces for right hand side
        feRhs->getOpPtrVector().push_back(
            new NonlinearElasticElement::OpRhsPiolaKirchhoff(
                "SPATIAL_POSITION", sit->second, elasticFe->commonData));
      }
    } else {
      feRhs->getOpPtrVector().push_back(
          new OpCalculateVectorFieldGradient<3, 3>(
              "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
      mat_grad_pos_at_pts_ptr = data_hooke_element_at_pts->HMat;
 
      // Evaluate density at integration points fro material like bone
      if (defaultMaterial == BONEHOOKE) {
        if (!mwlsApprox) {
          // Create instance for moving least square approximation
          mwlsApprox = boost::make_shared<MWLSApprox>(mField);
          // Load mesh with stresses
          CHKERR mwlsApprox->loadMWLSMesh(mwlsApproxFile.c_str());
          MeshsetsManager *block_manager_ptr;
          CHKERR mField.getInterface(block_manager_ptr);
          CHKERR block_manager_ptr->getEntitiesByDimension(
              densityMapBlock, BLOCKSET, 3, mwlsApprox->tetsInBlock, true);
          Tag th_rho;
          if (mwlsApprox->mwlsMoab.tag_get_handle(mwlsRhoTagName.c_str(),
                                                  th_rho) != MB_SUCCESS) {
            SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                     "Density tag name %s cannot be found. Please "
                     "provide the correct name.",
                     mwlsRhoTagName.c_str());
          }
          CHKERR mwlsApprox->mwlsMoab.tag_get_handle(mwlsRhoTagName.c_str(),
                                                     th_rho);
          CHKERR mwlsApprox->getValuesToNodes(th_rho);
        } else {
          mwlsApprox->tetsInBlock.clear();
          CHKERR mField.getInterface<MeshsetsManager>()->getEntitiesByDimension(
              densityMapBlock, BLOCKSET, 3, mwlsApprox->tetsInBlock, true);
        }
        // Calculate spatial positions and gradients (of deformation) at
        // integration pts. This operator takes into account singularity at
        // crack front
        feRhs->getOpPtrVector().push_back(
            new OpGetCrackFrontDataGradientAtGaussPts(
                "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
                feRhs->singularElement, feRhs->invSJac));
        // Evaluate density at integration pts.
        if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
          feRhs->getOpPtrVector().push_back(new MWLSApprox::OpMWLSRhoAtGaussPts(
              mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feRhs, mwlsApprox,
              mwlsRhoTagName, true, true));
        else {
          feRhs->getOpPtrVector().push_back(
              new MWLSApprox::OpMWLSCalculateBaseCoeffcientsAtGaussPts(
                  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feRhs,
                  mwlsApprox));
          feRhs->getOpPtrVector().push_back(
              new MWLSApprox::OpMWLSRhoAtGaussUsingPrecalulatedCoeffs(
                  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feRhs,
                  mwlsApprox, mwlsRhoTagName, true, false));
        }
 
        feRhs->getOpPtrVector().push_back(
            new HookeElement::OpCalculateStiffnessScaledByDensityField(
                "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
                data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts, nBone,
                rHo0));
        feRhs->getOpPtrVector().push_back(
            new HookeElement::OpCalculateStrainAle("SPATIAL_POSITION",
                                                   "SPATIAL_POSITION",
                                                   data_hooke_element_at_pts));
        feRhs->getOpPtrVector().push_back(
            new HookeElement::OpCalculateStress<1>("SPATIAL_POSITION",
                                                   "SPATIAL_POSITION",
                                                   data_hooke_element_at_pts));
      } else {
        feRhs->getOpPtrVector().push_back(
            new HookeElement::OpCalculateHomogeneousStiffness<0>(
                "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
                data_hooke_element_at_pts));
        // Calculate spatial positions and gradients (of deformation) at
        // integration pts. This operator takes into account singularity at
        // crack front
        feRhs->getOpPtrVector().push_back(
            new OpGetCrackFrontDataGradientAtGaussPts(
                "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
                feRhs->singularElement, feRhs->invSJac));
        feRhs->getOpPtrVector().push_back(
            new HookeElement::OpCalculateStrainAle("SPATIAL_POSITION",
                                                   "SPATIAL_POSITION",
                                                   data_hooke_element_at_pts));
        feRhs->getOpPtrVector().push_back(
            new HookeElement::OpCalculateStress<0>("SPATIAL_POSITION",
                                                   "SPATIAL_POSITION",
                                                   data_hooke_element_at_pts));
      }
 
      feRhs->getOpPtrVector().push_back(new OpTransfromSingularBaseFunctions(
          feRhs->singularElement, feRhs->detS, feRhs->invSJac));
      feRhs->getOpPtrVector().push_back(new HookeElement::OpAleRhs_dx(
          "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts));
    }
 
    // Add operators if internal stresses are present
    if (residualStressBlock != -1) {
      if (!mwlsApprox) {
        // Create instance for moving least square approximation
        mwlsApprox = boost::make_shared<MWLSApprox>(mField);
        // Load mesh with stresses
        CHKERR mwlsApprox->loadMWLSMesh(mwlsApproxFile.c_str());
        MeshsetsManager *block_manager_ptr;
        CHKERR mField.getInterface(block_manager_ptr);
        CHKERR block_manager_ptr->getEntitiesByDimension(
            residualStressBlock, BLOCKSET, 3, mwlsApprox->tetsInBlock, true);
        Tag th;
        if (mwlsApprox->mwlsMoab.tag_get_handle(mwls_stress_tag_name.c_str(),
                                                th) != MB_SUCCESS) {
          SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                   "Internal stress tag name %s cannot be found. Please "
                   "provide the correct name.",
                   mwls_stress_tag_name.substr(4).c_str());
        }
      }
 
      Tag th;
      if (mwlsApprox->mwlsMoab.tag_get_handle(mwls_stress_tag_name.c_str(),
                                              th) != MB_SUCCESS) {
        SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                 "Internal stress tag name %s cannot be found. Please "
                 "provide the correct name.",
                 mwls_stress_tag_name.substr(4).c_str());
      }
      CHKERR mwlsApprox->mwlsMoab.tag_get_handle(mwls_stress_tag_name.c_str(),
                                                 th);
      CHKERR mwlsApprox->getValuesToNodes(th);
 
      // Evaluate internal stresses at integration pts.
      if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
        feRhs->getOpPtrVector().push_back(
            new MWLSApprox::OpMWLSStressAtGaussPts(
                mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feRhs, mwlsApprox,
                mwls_stress_tag_name, false));
      else {
        feRhs->getOpPtrVector().push_back(
            new MWLSApprox::OpMWLSCalculateBaseCoeffcientsAtGaussPts(
                mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feRhs,
                mwlsApprox));
        feRhs->getOpPtrVector().push_back(
            new MWLSApprox::OpMWLSStressAtGaussUsingPrecalulatedCoeffs(
                mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feRhs, mwlsApprox,
                mwls_stress_tag_name, false));
      }
 
      // Assemble internall stresses forces
      feRhs->getOpPtrVector().push_back(new MWLSApprox::OpMWLSSpatialStressRhs(
          mat_grad_pos_at_pts_ptr, mwlsApprox));
      if (addAnalyticalInternalStressOperators) {
        feRhs->getOpPtrVector().push_back(
            new HookeElement::OpAnalyticalInternalAleStrain_dx<0>(
                "SPATIAL_POSITION", data_hooke_element_at_pts,
                analyticalStrainFunction,
                boost::shared_ptr<MatrixDouble>(
                    mwlsApprox, &mwlsApprox->mwlsMaterialPositions)));
      }
    }
 
    /// Lhs
    if (!onlyHooke) {
      // Calculate spatial positions and gradients (of deformation) at
      // integration pts. This operator takes into account singularity at crack
      // front
      feLhs->getOpPtrVector().push_back(new OpGetCrackFrontCommonDataAtGaussPts(
          "SPATIAL_POSITION", elasticFe->commonData, feLhs->singularElement,
          feLhs->invSJac));
      // Calculate material positions and gradients at integration pts.
      feLhs->getOpPtrVector().push_back(
          new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
              "MESH_NODE_POSITIONS", elasticFe->commonData));
      // Transform base functions to get singular base functions at crack front
      feLhs->getOpPtrVector().push_back(new OpTransfromSingularBaseFunctions(
          feLhs->singularElement, feLhs->detS, feLhs->invSJac));
      // Iterate over material blocks
      map<int, NonlinearElasticElement::BlockData>::iterator sit =
          elasticFe->setOfBlocks.begin();
      for (; sit != elasticFe->setOfBlocks.end(); sit++) {
        // Evaluate stress at integration pts
        feLhs->getOpPtrVector().push_back(
            new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
                "SPATIAL_POSITION", sit->second, elasticFe->commonData,
                ELASTIC_TAG, true, true, false));
        // Assemble tanget matrix, derivative of spatial poisotions
        feLhs->getOpPtrVector().push_back(
            new NonlinearElasticElement::OpLhsPiolaKirchhoff_dx(
                "SPATIAL_POSITION", "SPATIAL_POSITION", sit->second,
                elasticFe->commonData));
      }
    } else {
      // calculate material positions at integration points
      feLhs->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
          "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr));
      boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr;
      if (residualStressBlock != -1 || densityMapBlock != -1) {
        mat_grad_pos_at_pts_ptr = boost::make_shared<MatrixDouble>();
        feLhs->getOpPtrVector().push_back(
            new OpCalculateVectorFieldGradient<3, 3>("MESH_NODE_POSITIONS",
                                                     mat_grad_pos_at_pts_ptr));
      }
 
      feLhs->getOpPtrVector().push_back(
          new OpCalculateVectorFieldGradient<3, 3>(
              "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
      mat_grad_pos_at_pts_ptr = data_hooke_element_at_pts->HMat;
 
      if (defaultMaterial == BONEHOOKE) {
        if (!mwlsApprox) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "mwlsApprox not allocated");
        } else {
          mwlsApprox->tetsInBlock.clear();
          CHKERR mField.getInterface<MeshsetsManager>()->getEntitiesByDimension(
              densityMapBlock, BLOCKSET, 3, mwlsApprox->tetsInBlock,
              true); // FIXME: do we still need this?
        }
 
        // Calculate spatial positions and gradients (of deformation) at
        // integration pts. This operator takes into account singularity at
        // crack front
        feLhs->getOpPtrVector().push_back(
            new OpGetCrackFrontDataGradientAtGaussPts(
                "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
                feLhs->singularElement, feLhs->invSJac));
 
        // Evaluate density at integration pts.
        if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
          feLhs->getOpPtrVector().push_back(new MWLSApprox::OpMWLSRhoAtGaussPts(
              mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feLhs, mwlsApprox,
              mwlsRhoTagName, true, true));
        else {
          feLhs->getOpPtrVector().push_back(
              new MWLSApprox::OpMWLSCalculateBaseCoeffcientsAtGaussPts(
                  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feLhs,
                  mwlsApprox));
          feLhs->getOpPtrVector().push_back(
              new MWLSApprox::OpMWLSRhoAtGaussUsingPrecalulatedCoeffs(
                  mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feLhs,
                  mwlsApprox, mwlsRhoTagName, true, false));
        }
 
        feLhs->getOpPtrVector().push_back(
            new HookeElement::OpCalculateStiffnessScaledByDensityField(
                "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
                data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts, nBone,
                rHo0));
        feLhs->getOpPtrVector().push_back(
            new HookeElement::OpCalculateStrainAle("SPATIAL_POSITION",
                                                   "SPATIAL_POSITION",
                                                   data_hooke_element_at_pts));
        feLhs->getOpPtrVector().push_back(
            new HookeElement::OpCalculateStress<1>("SPATIAL_POSITION",
                                                   "SPATIAL_POSITION",
                                                   data_hooke_element_at_pts));
 
        feLhs->getOpPtrVector().push_back(new OpTransfromSingularBaseFunctions(
            feLhs->singularElement, feLhs->detS, feLhs->invSJac));
        feLhs->getOpPtrVector().push_back(new HookeElement::OpAleLhs_dx_dx<1>(
            "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts));
 
      } else {
        feLhs->getOpPtrVector().push_back(
            new HookeElement::OpCalculateHomogeneousStiffness<0>(
                "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
                data_hooke_element_at_pts));
        // Calculate spatial positions and gradients (of deformation) at
        // integration pts. This operator takes into account singularity at
        // crack front
        feLhs->getOpPtrVector().push_back(
            new OpGetCrackFrontDataGradientAtGaussPts(
                "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
                feLhs->singularElement, feLhs->invSJac));
        feLhs->getOpPtrVector().push_back(
            new HookeElement::OpCalculateStrainAle("SPATIAL_POSITION",
                                                   "SPATIAL_POSITION",
                                                   data_hooke_element_at_pts));
        feLhs->getOpPtrVector().push_back(
            new HookeElement::OpCalculateStress<0>("SPATIAL_POSITION",
                                                   "SPATIAL_POSITION",
                                                   data_hooke_element_at_pts));
        // Transform base functions to get singular base functions at crack
        // front
        feLhs->getOpPtrVector().push_back(new OpTransfromSingularBaseFunctions(
            feLhs->singularElement, feLhs->detS, feLhs->invSJac));
        feLhs->getOpPtrVector().push_back(new HookeElement::OpAleLhs_dx_dx<0>(
            "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts));
      }
    }
  }
 
  // Assembly forces
 
  surfaceForces =
      boost::make_shared<boost::ptr_map<string, NeumannForcesSurface>>();
  // Surface forces
  {
    string fe_name_str = "FORCE_FE";
    surfaceForces->insert(fe_name_str, new NeumannForcesSurface(mField));
    auto &nf = surfaceForces->at(fe_name_str);
    auto &fe = nf.getLoopFe();
 
    CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", fe, false, false);
 
    for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(mField, NODESET | FORCESET,
                                                    it)) {
      CHKERR surfaceForces->at(fe_name_str)
          .addForce("SPATIAL_POSITION", PETSC_NULL, it->getMeshsetId(), true);
    }
 
    const string block_set_force_name("FORCE");
    // search for block named FORCE and add its attributes to FORCE_FE element
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
      if (it->getName().compare(0, block_set_force_name.length(),
                                block_set_force_name) == 0) {
        CHKERR surfaceForces->at(fe_name_str)
            .addForce("SPATIAL_POSITION", PETSC_NULL, (it->getMeshsetId()),
                      true, true);
      }
    }
  }
 
  if (isSurfaceForceAle) {
    surfaceForceAle =
        boost::make_shared<boost::ptr_map<string, NeumannForcesSurface>>();
    commonDataSurfaceForceAle =
        boost::make_shared<NeumannForcesSurface::DataAtIntegrationPts>();
    commonDataSurfaceForceAle->forcesOnlyOnEntitiesRow = crackFrontNodes;
 
    string fe_name_str = "FORCE_FE_ALE";
    surfaceForceAle->insert(fe_name_str, new NeumannForcesSurface(mField));
 
    auto &nf = surfaceForceAle->at(fe_name_str);
    auto &fe_lhs = nf.getLoopFeLhs();
    auto &fe_mat_rhs = nf.getLoopFeMatRhs();
    auto &fe_mat_lhs = nf.getLoopFeMatLhs();
 
    CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", fe_lhs, false, false);
    CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", fe_mat_rhs, false, false);
    CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", fe_mat_lhs, false, false);
 
    for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(mField, NODESET | FORCESET,
                                                    bit)) {
      CHKERR surfaceForceAle->at(fe_name_str)
          .addForceAle("SPATIAL_POSITION", "MESH_NODE_POSITIONS",
                       commonDataSurfaceForceAle, "MATERIAL", PETSC_NULL,
                       PETSC_NULL, bit->getMeshsetId(), true, false,
                       ignoreMaterialForce);
    }
 
    const string block_set_force_name("FORCE");
    // search for block named FORCE and add its attributes to FORCE_FE element
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
      if (it->getName().compare(0, block_set_force_name.length(),
                                block_set_force_name) == 0) {
        CHKERR surfaceForceAle->at(fe_name_str)
            .addForceAle("SPATIAL_POSITION", "MESH_NODE_POSITIONS",
                         commonDataSurfaceForceAle, "MATERIAL", PETSC_NULL,
                         PETSC_NULL, it->getMeshsetId(), true, true,
                         ignoreMaterialForce);
      }
    }
  }
 
  // assemble surface pressure
  surfacePressure =
      boost::make_shared<boost::ptr_map<string, NeumannForcesSurface>>();
  {
    string fe_name_str = "PRESSURE_FE";
    surfacePressure->insert(fe_name_str, new NeumannForcesSurface(mField));
 
    auto &nf = surfacePressure->at(fe_name_str);
    auto &fe = nf.getLoopFe();
 
    CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", fe, false, false);
 
    // iterate over sidestep where pressure is applied
    for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
             mField, SIDESET | PRESSURESET, it)) {
      CHKERR surfacePressure->at(fe_name_str)
          .addPressure("SPATIAL_POSITION", PETSC_NULL, it->getMeshsetId(),
                       true);
    }
 
    const string block_set_pressure_name("PRESSURE");
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
      if (it->getName().compare(0, block_set_pressure_name.length(),
                                block_set_pressure_name) == 0) {
        CHKERR surfacePressure->at(fe_name_str)
            .addPressure("SPATIAL_POSITION", PETSC_NULL, it->getMeshsetId(),
                         true, true);
      }
    }
 
    const string block_set_linear_pressure_name("LINEAR_PRESSURE");
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
      if (it->getName().compare(0, block_set_linear_pressure_name.length(),
                                block_set_linear_pressure_name) == 0) {
        CHKERR surfacePressure->at(fe_name_str)
            .addLinearPressure("SPATIAL_POSITION", PETSC_NULL,
                               it->getMeshsetId(), true);
      }
    }
  }
  // assemble surface pressure (ALE)
  if (isPressureAle) {
    surfacePressureAle =
        boost::make_shared<boost::ptr_map<string, NeumannForcesSurface>>();
    commonDataSurfacePressureAle =
        boost::make_shared<NeumannForcesSurface::DataAtIntegrationPts>();
    commonDataSurfacePressureAle->forcesOnlyOnEntitiesRow = crackFrontNodes;
    string fe_name_str = "PRESSURE_ALE";
    surfacePressureAle->insert(fe_name_str, new NeumannForcesSurface(mField));
 
    auto &nf = surfacePressureAle->at(fe_name_str);
    auto &fe_lhs = nf.getLoopFeLhs();
    auto &fe_mat_rhs = nf.getLoopFeMatRhs();
    auto &fe_mat_lhs = nf.getLoopFeMatLhs();
 
    CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", fe_lhs, false, false);
    CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", fe_mat_rhs, false, false);
    CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", fe_mat_lhs, false, false);
 
    // iterate over sidesets where pressure is applied
    for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(
             mField, SIDESET | PRESSURESET, it)) {
      CHKERR surfacePressureAle->at(fe_name_str)
          .addPressureAle("SPATIAL_POSITION", "MESH_NODE_POSITIONS",
                          commonDataSurfacePressureAle, "MATERIAL", PETSC_NULL,
                          PETSC_NULL, it->getMeshsetId(), true);
    }
  }
 
  if (areSpringsAle) {
    feRhsSpringALEMaterial =
        boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
 
    feLhsSpringALE =
        boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
    feLhsSpringALEMaterial =
        boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
    commonDataSpringsALE =
        boost::make_shared<MetaSpringBC::DataAtIntegrationPtsSprings>(mField);
    commonDataSpringsALE->forcesOnlyOnEntitiesRow = crackFrontNodes;
 
    CHKERR MetaSpringBC::setSpringOperatorsMaterial(
        mField, feLhsSpringALE, feLhsSpringALEMaterial, feRhsSpringALEMaterial,
        commonDataSpringsALE, "SPATIAL_POSITION", "MESH_NODE_POSITIONS",
        "MATERIAL");
  }
 
  // Implementation of spring element
  // Create new instances of face elements for springs
  feSpringLhsPtr = boost::shared_ptr<FaceElementForcesAndSourcesCore>(
      new FaceElementForcesAndSourcesCore(mField));
  feSpringRhsPtr = boost::shared_ptr<FaceElementForcesAndSourcesCore>(
      new FaceElementForcesAndSourcesCore(mField));
 
  CHKERR MetaSpringBC::setSpringOperators(mField, feSpringLhsPtr,
                                          feSpringRhsPtr, "SPATIAL_POSITION",
                                          "MESH_NODE_POSITIONS");
 
  bool use_eigen_pos =
      solveEigenStressProblem && (mwls_stress_tag_name == mwlsStressTagName);
 
  bool use_eigen_pos_simple_contact =
      use_eigen_pos && useEigenPositionsSimpleContact;
 
  // Set contact operators
  feRhsSimpleContact =
      boost::make_shared<SimpleContactProblem::SimpleContactElement>(mField);
  feLhsSimpleContact =
      boost::make_shared<SimpleContactProblem::SimpleContactElement>(mField);
  commonDataSimpleContact =
      boost::make_shared<SimpleContactProblem::CommonDataSimpleContact>(mField);
 
  if (!ignoreContact) {
    if (printContactState) {
      feRhsSimpleContact->contactStateVec =
          commonDataSimpleContact->gaussPtsStateVec;
    }
    contactProblem->setContactOperatorsRhs(
        feRhsSimpleContact, commonDataSimpleContact, "SPATIAL_POSITION",
        "LAMBDA_CONTACT", false, use_eigen_pos_simple_contact,
        "EIGEN_SPATIAL_POSITIONS");
    contactProblem->setMasterForceOperatorsRhs(
        feRhsSimpleContact, commonDataSimpleContact, "SPATIAL_POSITION",
        "LAMBDA_CONTACT", false, use_eigen_pos_simple_contact,
        "EIGEN_SPATIAL_POSITIONS");
 
    contactProblem->setContactOperatorsLhs(
        feLhsSimpleContact, commonDataSimpleContact, "SPATIAL_POSITION",
        "LAMBDA_CONTACT", false, use_eigen_pos_simple_contact,
        "EIGEN_SPATIAL_POSITIONS");
 
    contactProblem->setMasterForceOperatorsLhs(
        feLhsSimpleContact, commonDataSimpleContact, "SPATIAL_POSITION",
        "LAMBDA_CONTACT", false, use_eigen_pos_simple_contact,
        "EIGEN_SPATIAL_POSITIONS");
  }
 
  // Close crack constrains
  if (solveEigenStressProblem) {
    closeCrackConstrains =
        boost::shared_ptr<BothSurfaceConstrains>(new BothSurfaceConstrains(
            mField, "LAMBDA_CLOSE_CRACK", "SPATIAL_POSITION"));
    closeCrackConstrains->aLpha = 1;
  }
 
  // Set contact operators
 
  if (!contactElements.empty() && !ignoreContact && !fixContactNodes) {
    // Set contact operators
    feRhsSimpleContactALEMaterial =
        boost::make_shared<SimpleContactProblem::SimpleContactElement>(mField);
    feLhsSimpleContactALEMaterial =
        boost::make_shared<SimpleContactProblem::SimpleContactElement>(mField);
    feLhsSimpleContactALE =
        boost::make_shared<SimpleContactProblem::SimpleContactElement>(mField);
    commonDataSimpleContactALE =
        boost::make_shared<SimpleContactProblem::CommonDataSimpleContact>(
            mField);
 
    commonDataSimpleContactALE->forcesOnlyOnEntitiesRow = crackFrontNodes;
 
    contactProblem->setContactOperatorsRhsALEMaterial(
        feRhsSimpleContactALEMaterial, commonDataSimpleContactALE,
        "SPATIAL_POSITION", "MESH_NODE_POSITIONS", "LAMBDA_CONTACT",
        "MAT_CONTACT");
 
    contactProblem->setContactOperatorsLhsALEMaterial(
        feLhsSimpleContactALEMaterial, commonDataSimpleContactALE,
        "SPATIAL_POSITION", "MESH_NODE_POSITIONS", "LAMBDA_CONTACT",
        "MAT_CONTACT");
 
    contactProblem->setContactOperatorsLhsALE(
        feLhsSimpleContactALE, commonDataSimpleContactALE, "SPATIAL_POSITION",
        "MESH_NODE_POSITIONS", "LAMBDA_CONTACT", use_eigen_pos_simple_contact,
        "EIGEN_SPATIAL_POSITIONS");
  }
 
  feRhsMortarContact =
      boost::make_shared<MortarContactProblem::MortarContactElement>(
          mField, contactSearchMultiIndexPtr, "SPATIAL_POSITION",
          "MESH_NODE_POSITIONS");
  feLhsMortarContact =
      boost::make_shared<MortarContactProblem::MortarContactElement>(
          mField, contactSearchMultiIndexPtr, "SPATIAL_POSITION",
          "MESH_NODE_POSITIONS");
  commonDataMortarContact =
      boost::make_shared<MortarContactProblem::CommonDataMortarContact>(mField);
 
  if (!ignoreContact) {
    if (printContactState) {
      feRhsMortarContact->contactStateVec =
          commonDataMortarContact->gaussPtsStateVec;
    }
    mortarContactProblemPtr->setContactOperatorsRhs(
        feRhsMortarContact, commonDataMortarContact, "SPATIAL_POSITION",
        "LAMBDA_CONTACT", false, use_eigen_pos, "EIGEN_SPATIAL_POSITIONS");
 
    mortarContactProblemPtr->setMasterForceOperatorsRhs(
        feRhsMortarContact, commonDataMortarContact, "SPATIAL_POSITION",
        "LAMBDA_CONTACT", false, use_eigen_pos, "EIGEN_SPATIAL_POSITIONS");
 
    mortarContactProblemPtr->setContactOperatorsLhs(
        feLhsMortarContact, commonDataMortarContact, "SPATIAL_POSITION",
        "LAMBDA_CONTACT", false, use_eigen_pos, "EIGEN_SPATIAL_POSITIONS");
 
    mortarContactProblemPtr->setMasterForceOperatorsLhs(
        feLhsMortarContact, commonDataMortarContact, "SPATIAL_POSITION",
        "LAMBDA_CONTACT", false, use_eigen_pos, "EIGEN_SPATIAL_POSITIONS");
  }
 
  // set contact post proc operators
 
  fePostProcSimpleContact =
      boost::make_shared<SimpleContactProblem::SimpleContactElement>(mField);
 
  fePostProcMortarContact =
      boost::make_shared<MortarContactProblem::MortarContactElement>(
          mField, contactSearchMultiIndexPtr, "SPATIAL_POSITION",
          "MESH_NODE_POSITIONS");
 
  if (!ignoreContact) {
    CHKERR contactProblem->setContactOperatorsForPostProc(
        fePostProcSimpleContact, commonDataSimpleContact, mField,
        "SPATIAL_POSITION", "LAMBDA_CONTACT", contactPostProcMoab, false,
        use_eigen_pos_simple_contact, "EIGEN_SPATIAL_POSITIONS");
 
    CHKERR mortarContactProblemPtr->setContactOperatorsForPostProc(
        fePostProcMortarContact, commonDataMortarContact, mField,
        "SPATIAL_POSITION", "LAMBDA_CONTACT", contactPostProcMoab, false,
        use_eigen_pos, "EIGEN_SPATIAL_POSITIONS");
  }
 
  // assemble nodal forces
  nodalForces = boost::make_shared<boost::ptr_map<string, NodalForce>>();
  {
    string fe_name_str = "FORCE_FE";
    nodalForces->insert(fe_name_str, new NodalForce(mField));
    for (_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(mField, NODESET | FORCESET,
                                                    it)) {
      CHKERR nodalForces->at(fe_name_str)
          .addForce("SPATIAL_POSITION", PETSC_NULL, it->getMeshsetId(), false);
    }
  }
  // assemble edge forces
  edgeForces = boost::make_shared<boost::ptr_map<string, EdgeForce>>();