EshelbianCore.hpp

Go to the documentation of this file.

/**
 * @file EshelbianCore.hpp
 * @author your name (you@domain.com)
 * @brief
 * @version 0.1
 * @date 2024-12-31
 *
 * @copyright Copyright (c) 2024
 *
 */
 
struct EshelbianCore : public MoFEM::UnknownInterface {
 
  static inline const char *listSolvers[] = {
      "time_solver",
      "dynamic_relaxation",
      "incremental_optimization",
      "load_factor",
      "shape_optimization",
      "test_topological_derivative",
      "test_equilibrated_mechanical_value",
      "test_incremental_optimization_layout",
      "test_incremental_optimization_transaction",
      "test_incremental_optimization_objective_derivative",
      "test_incremental_optimization_constraint_derivative"};
 
  enum SolverType {
    TimeSolver,
    DynamicRelaxation,
    IncrementalOptimization,
    LoadFactor,
    ShapeOptimization,
    TestTopologicalDerivative,
    TestEquilibratedMechanicalValue,
    TestIncrementalOptimizationLayout,
    TestIncrementalOptimizationTransaction,
    TestIncrementalOptimizationObjectiveDerivative,
    TestIncrementalOptimizationConstraintDerivative,
    LastSolver
  };
 
  // That exploits fact that full system is symmetric
  static inline constexpr enum SymmetrySelector symmetrySelector = SYMMETRIC;
 
  static inline enum SolverType solverType = TimeSolver;
  static inline enum RotSelector rotSelector = LARGE_ROT;
  static inline enum RotSelector gradApproximator = LARGE_ROT;
  static inline enum StretchSelector stretchSelector = LOG;
  static inline PetscBool setSingularity = PETSC_FALSE; //< set singularity
  static inline PetscBool physicalTimeFlg =
      PETSC_FALSE; //< switch on/off dynamic relaxation
  static inline PetscBool crackingOn = PETSC_FALSE; //< cracking on
  static inline PetscBool propagateUnderCompression =
      PETSC_TRUE; //< allow crack propagation under compression
  static inline double crackingStartTime = 0; //< time when crack starts to grow
  static inline double crackingAddTime = 0;   //< time to add new crack surface
  static inline bool noCrackExtension = false; //< flag for no crack extension
  static inline int nbStepsNoCrackExtension =
      0; //< number of steps with no crack extension
  static inline bool potentialCrackArrest =
      false; //< flag potential crack arrest
  static inline int nbJIntegralContours =
      0; //< number of contours for J integral evaluation
  static inline double finalPhysicalTime =
      0; //< Maximum step for non-time solver
  static inline double currentPhysicalTime =
      0;                               //< Current step for non-time solver
  static inline double physicalDt = 0; //< Step size for non-time solver
  static inline int physicalMaxSteps =
      20; //< maximum iterations for non-time solver
  static inline int physicalStepNumber = 0; //< Step number for non-time solver
  static inline PetscBool physicalH1Update =
      PETSC_FALSE; //< update H1 space at each non-time solver step
  static inline PetscBool l2UserBaseScale = PETSC_FALSE; //< scale L2 user base
  static inline int addCrackMeshsetId = 1000; //< add crack meshset id
  static inline double griffithEnergy = 1;    ///< Griffith energy
  static inline double crackingRtol = 1e-10;  ///< Cracking relative tolerance
  static inline double crackingAtol = 1e-12;  ///< Cracking absolute tolerance
  static inline enum EnergyReleaseSelector energyReleaseSelector =
      GRIFFITH_SKELETON; //< energy release selector
  static inline std::string internalStressTagName =
      ""; //< internal stress tag name
  static inline PetscBool internalStressVoigt =
      PETSC_FALSE; //< internal stress index notation
  static inline PetscBool interfaceCrack =
      PETSC_FALSE; //< interface crack tracking
  static inline PetscBool plasticVolume =
      PETSC_FALSE; //< restrict plasticity to the plastic volume block
 
  static inline int interfaceRemoveLevel =
      0; //< number of levels of elements to remove around interface crack
  static inline std::string heterogeneousYoungModTagName =
      ""; //< heterogenous young's modulus
  static inline std::string meshTransferSourceMeshFileName =
      ""; //< source mesh file name for projection material tags
  static inline int meshTransferInterpOrder =
      1; //< interpolation order for projection material tags
  static inline PetscBool meshTransferSourceMeshFileSpecified =
      PETSC_FALSE; //< flag to check if source mesh file is specified for
                   // projection material tags
  static inline PetscBool meshTransferHybridInterp =
      PETSC_TRUE; //< flag to use hybrid interpolation for projection material
                  // tags
  static inline std::vector<std::string>
      listTagsToProject; // list of tags to project from source mesh to target
                         // mesh
  static inline FieldApproximationBase brokenHdivBase =
      DEMKOWICZ_JACOBI_BASE; //< approximation base for broken HDIV stress
 
  template <FieldApproximationBase HdivBase> struct FieldOrders;
 
  template <typename Op> MoFEMErrorCode withFieldOrders(Op &&op) const {
    MoFEMFunctionBegin;
    switch (EshelbianCore::brokenHdivBase) {
    case DEMKOWICZ_JACOBI_BASE:
      CHKERR op.template operator()<DEMKOWICZ_JACOBI_BASE>();
      break;
    case AINSWORTH_LEGENDRE_BASE:
      CHKERR op.template operator()<AINSWORTH_LEGENDRE_BASE>();
      break;
    default:
      SETERRQ(mField.get_comm(), MOFEM_NOT_IMPLEMENTED,
              "Broken HDIV base not implemented");
    }
    MoFEMFunctionReturn(0);
  }
 
  static inline double maxCrackExtension =
      50; //< maximum crack extension in one physical step
 
  static boost::function<double(const double)> f;
  static boost::function<double(const double)> d_f;
  static boost::function<double(const double)> dd_f;
  static boost::function<double(const double)> inv_f;
  static boost::function<double(const double)> inv_d_f;
  static boost::function<double(const double)> inv_dd_f;
 
  static inline constexpr double v_max = 24;
 
  static double f_log_e_quadratic(const double v) {
    if (v > v_max) {
      double e = static_cast<double>(std::exp(v_max));
      double dv = v - v_max;
      return 0.5 * e * dv * dv + e * dv + e;
    } else {
      return static_cast<double>(std::exp(v));
    }
  }
 
  static double d_f_log_e_quadratic(const double v) {
    if (v > v_max) {
      double e = static_cast<double>(std::exp(v_max));
      double dv = v - v_max;
      return e * dv + e;
    } else {
      return static_cast<double>(std::exp(v));
    }
  }
 
  static double dd_f_log_e_quadratic(const double v) {
    if (v > v_max) {
      return static_cast<double>(std::exp(v_max));
    } else {
      return static_cast<double>(std::exp(v));
    }
  }
 
  static double inv_f_log_e_quadratic(const double stretch) {
    const double transition_stretch = std::exp(v_max);
    if (stretch <= transition_stretch) {
      return std::log(stretch);
    }
    return v_max - 1. + std::sqrt(2. * stretch / transition_stretch - 1.);
  }
 
  static double inv_d_f_log_e_quadratic(const double stretch) {
    const double transition_stretch = std::exp(v_max);
    if (stretch <= transition_stretch) {
      return 1. / stretch;
    }
    const double root = std::sqrt(2. * stretch / transition_stretch - 1.);
    return 1. / (transition_stretch * root);
  }
 
  static double inv_dd_f_log_e_quadratic(const double stretch) {
    const double transition_stretch = std::exp(v_max);
    if (stretch <= transition_stretch) {
      return -1. / (stretch * stretch);
    }
    const double root = std::sqrt(2. * stretch / transition_stretch - 1.);
    return -1. / (transition_stretch * transition_stretch * root * root * root);
  }
 
  static double f_log_e(const double v) { return std::exp(v); }
  static double d_f_log_e(const double v) { return std::exp(v); }
  static double dd_f_log_e(const double v) { return std::exp(v); }
  static double inv_f_log_e(const double v) { return std::log(v); }
  static double inv_d_f_log_e(const double v) { return 1. / v; }
  static double inv_dd_f_log_e(const double v) { return -1. / (v * v); }
 
  static double f_linear(const double v) { return v + 1; }
  static double d_f_linear(const double) { return 1; }
  static double dd_f_linear(const double) { return 0; }
 
  static double inv_f_linear(const double v) { return v - 1; }
  static double inv_d_f_linear(const double) { return 1; }
  static double inv_dd_f_linear(const double) { return 0; }
 
  /**
   * \brief Getting interface of core database
   * @param  uuid  unique ID of interface
   * @param  iface returned pointer to interface
   * @return       error code
   */
  MoFEMErrorCode query_interface(boost::typeindex::type_index type_index,
                                 UnknownInterface **iface) const;
 
  MoFEM::Interface &mField;
 
  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
  boost::shared_ptr<PhysicalEquations> physicalEquations;
  boost::shared_ptr<AnalyticalExprPython> AnalyticalExprPythonPtr;
 
  boost::shared_ptr<VolumeElementForcesAndSourcesCore> elasticFeRhs;
  boost::shared_ptr<VolumeElementForcesAndSourcesCore> elasticFeLhs;
  boost::shared_ptr<FaceElementForcesAndSourcesCore> elasticBcLhs;
  boost::shared_ptr<FaceElementForcesAndSourcesCore> elasticBcRhs;
  boost::shared_ptr<ForcesAndSourcesCore>
      contactTreeRhs; ///< Make a contact tree
 
  SmartPetscObj<DM> dM;           ///< Coupled problem all fields
  SmartPetscObj<DM> dmElastic;    ///< Elastic problem
  SmartPetscObj<DM> dmMaterial;   ///< Material problem
  SmartPetscObj<DM> dmPrjSpatial; ///< Projection spatial displacement
  SmartPetscObj<DM>
      dmIncrementalOptimization; ///< Incremental-optimization control problem
  /** Trial control read directly by fixed-control FE pipelines. */
  SmartPetscObj<Vec> incrementalTrialControl;
 
  const std::string piolaStress = "P";
  const std::string spatialL2Disp = "wL2";
  const std::string spatialH1Disp = "wH1";
  const std::string materialH1Positions = "XH1";
  const std::string hybridSpatialDisp = "hybridSpatialDisp";
  const std::string contactDisp = "contactDisp";
  const std::string stretchTensor = "u";
  const std::string logDeviator = "D";
  const std::string logJacobian = "theta";
  const std::string auxiliaryLogStress = "Td";
  const std::string rotAxis = "omega";
  const std::string bubbleField = "bubble";
  const std::string plasticFlowField = "plasticFlow";
  const std::string plasticKappaField = "plasticKappa";
  const std::string plasticHField = "plasticH";
 
  const std::string elementVolumeName = "EP";
  const std::string naturalBcElement = "NATURAL_BC";
  const std::string skinElement = "SKIN";
  const std::string skeletonElement = "SKELETON";
  const std::string contactElement = "CONTACT";
 
  EshelbianCore(MoFEM::Interface &m_field);
  virtual ~EshelbianCore();
 
  int spaceOrder = 2;
  int spaceH1Order = -1;
  int materialH1Order = 1;
  double alphaU = 0;
  double alphaW = 0;
  double alphaOmega = 0;
  double alphaR = 0;
  double alphaViscousOmega = 0;
  double alphaViscousR = 0;
  double alphaRho = 0;
  double alphaTau = 0;
  double alphaTauLin = 0;
  double alphaTauBcDisp = 0;
  double alphaTauBcDisp0 = 0;
  double dynamicAtol = 0;
  double dynamicRtol = 0;
  double dynamicInitialResidual = -1;
 
  int contactRefinementLevels = 1; //< refinement levels for contact integration
  int frontLayers = 3;             //< front layers with material element
  double loadFactor = 1.0;         //< load factor
  double avgGriffithsEnergy = -1;
 
  MoFEMErrorCode getOptions();
  MoFEMErrorCode applyTestSolverMonitorOptions(TS ts);
  MoFEMErrorCode applyProjectionSolverMonitorOptions();
 
  boost::shared_ptr<BcDispVec> bcSpatialDispVecPtr;
  boost::shared_ptr<BcRotVec> bcSpatialRotationVecPtr;
  boost::shared_ptr<TractionBcVec> bcSpatialTractionVecPtr;
  boost::shared_ptr<TractionFreeBc> bcSpatialFreeTractionVecPtr;
  boost::shared_ptr<NormalDisplacementBcVec> bcSpatialNormalDisplacementVecPtr;
  boost::shared_ptr<SpringBcVec> bcSpatialSpringVecPtr;
  boost::shared_ptr<AnalyticalDisplacementBcVec>
      bcSpatialAnalyticalDisplacementVecPtr;
  boost::shared_ptr<AnalyticalTractionBcVec> bcSpatialAnalyticalTractionVecPtr;
  boost::shared_ptr<PressureBcVec> bcSpatialPressureVecPtr;
  boost::shared_ptr<ExternalStrainVec> externalStrainVecPtr;
  double oldCrackArea = 0.;
  double oldStrainEnergy = 0.;
  double oldLoadFactor = 1.0;
  double strainEnergy = 0.;
  boost::shared_ptr<double> currentCrackAreaPtr;
 
  std::map<std::string, boost::shared_ptr<ScalingMethod>> timeScaleMap;
 
  std::string getStringArgumentFromJsonBlockset(const std::string &type_name,
                                                const int meshset_id,
                                                const std::string &param_name) {
    const auto string_params =
        mField.getInterface<JsonConfigManager>()->getStringParamsFromBlockset(
            type_name, meshset_id);
    if (const auto it = string_params.find(param_name);
        it != string_params.end()) {
      return it->second;
    }
    return "";
  }
 
  MoFEMErrorCode
  getStringArgumentFromJsonBlocksets(const std::string &type_name,
                                     const std::string &param_name,
                                     std::string &param_value) {
    MoFEMFunctionBegin;
    param_value.clear();
    for (auto it : mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
             std::regex((boost::format("%s(.*)") % type_name).str()))) {
      const auto block_param = getStringArgumentFromJsonBlockset(
          type_name, it->getMeshsetId(), param_name);
      if (block_param.empty()) {
        continue;
      }
      if (!param_value.empty() && param_value != block_param) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "JSON %s blocksets specify different '%s' values",
                type_name.c_str(), param_name.c_str());
      }
      param_value = block_param;
    }
    MoFEMFunctionReturn(0);
  }
 
  template <typename BC>
  MoFEMErrorCode getBc(boost::shared_ptr<BC> &bc_vec_ptr,
                       const std::string block_name, const int nb_attributes) {
    MoFEMFunctionBegin;
    for (auto it :
         mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
 
             (boost::format("%s(.*)") % block_name).str()
 
                 ))
 
    ) {
      std::vector<double> block_attributes;
      CHKERR it->getAttributes(block_attributes);
      if (block_attributes.size() < static_cast<size_t>(nb_attributes)) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "In block %s expected %d attributes, but given %ld",
                it->getName().c_str(), nb_attributes, block_attributes.size());
      }
      Range faces;
      CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
                                                 true);
      bc_vec_ptr->emplace_back(
          it->getName(), block_attributes, faces,
          getStringArgumentFromJsonBlockset(block_name, it->getMeshsetId(),
                                            "load_history"));
    }
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode getSpatialDispBc();
 
  inline MoFEMErrorCode getSpatialRotationBc() {
    MoFEMFunctionBegin;
    bcSpatialRotationVecPtr = boost::make_shared<BcRotVec>();
    CHKERR getBc(bcSpatialRotationVecPtr, "SPATIAL_ROTATION_BC", 4);
    CHKERR getBc(bcSpatialRotationVecPtr, "SPATIAL_ROTATION_AXIS_BC", 7);
 
    for (auto &bc : *bcSpatialRotationVecPtr) {
      MOFEM_LOG("EP", Sev::inform)
          << "Found spatial rotation BC on block " << bc.blockName;
      MOFEM_LOG("EP", Sev::inform) << " with attributes: " << bc.vals;
      MOFEM_LOG("EP", Sev::inform) << " and rotation angle: " << bc.theta;
      MOFEM_LOG("EP", Sev::inform) << " and nb of faces: " << bc.faces.size();
    }
 
    auto ts_rotation =
        boost::make_shared<DynamicRelaxationTimeScale>("rotation_history.txt");
    for (auto &bc : *bcSpatialRotationVecPtr) {
      if (!bc.loadHistoryFile.empty()) {
        MOFEM_LOG("EP", Sev::inform)
            << "Rotation load history from JSON for " << bc.blockName << ": "
            << bc.loadHistoryFile;
        timeScaleMap[bc.blockName] =
            boost::make_shared<DynamicRelaxationTimeScale>(bc.loadHistoryFile);
      } else {
        timeScaleMap[bc.blockName] =
            GetBlockScalingMethod<DynamicRelaxationTimeScale>::get(
                ts_rotation, "rotation_history", ".txt", bc.blockName);
      }
    }
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode getSpatialTractionBc();
 
  /**
   * @brief Remove all, but entities where kinematic constrains are applied.
   *
   * @param meshset
   * @param bc_ptr
   * @param disp_block_set_name
   * @param rot_block_set_name
   * @param contact_set_name
   * @return MoFEMErrorCode
   */
  MoFEMErrorCode getTractionFreeBc(const EntityHandle meshset,
                                   boost::shared_ptr<TractionFreeBc> &bc_ptr,
                                   const std::string contact_set_name);
 
  inline MoFEMErrorCode
  getSpatialTractionFreeBc(const EntityHandle meshset = 0) {
    bcSpatialFreeTractionVecPtr =
        boost::shared_ptr<TractionFreeBc>(new TractionFreeBc());
    return getTractionFreeBc(meshset, bcSpatialFreeTractionVecPtr, "CONTACT");
  }
 
  MoFEMErrorCode getExternalStrain();
 
  MoFEMErrorCode createExchangeVectors(Sev sev);
 
  MoFEMErrorCode resolveDissipationEntities(const EntityHandle meshset = 0);
  MoFEMErrorCode addFields(const EntityHandle meshset = 0,
                           const bool add_bubble = true);
  MoFEMErrorCode projectGeometry(const EntityHandle meshset = 0,
                                 double time = 0);
  MoFEMErrorCode projectMaterialTags(const EntityHandle meshset = 0);
 
  MoFEMErrorCode addVolumeFiniteElement(const EntityHandle meshset = 0,
                                        const bool add_bubble = true);
  MoFEMErrorCode addBoundaryFiniteElement(const EntityHandle meshset = 0);
  MoFEMErrorCode addDMs(const BitRefLevel bit = BitRefLevel().set(0),
                        const EntityHandle meshset = 0);
 
  MoFEMErrorCode setBaseVolumeElementOps(
      const int tag, const bool do_rhs, const bool do_lhs,
      const bool calc_rates,
      boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe,
      const bool add_bubble = true);
 
  MoFEMErrorCode setVolumeElementOps(
      const int tag, const bool add_elastic, const bool add_material,
      boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_rhs,
      boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_lhs);
 
  MoFEMErrorCode
  pushVolumeA00Ops(boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_lhs);
 
  MoFEMErrorCode pushStressGramOps(
      boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_lhs);
 
  MoFEMErrorCode pushPiolaStressGramOps(
      boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_lhs);
 
  MoFEMErrorCode
  setFaceElementOps(const bool add_elastic, const bool add_material,
                    boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_rhs,
                    boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_lhs);
 
  MoFEMErrorCode setContactElementRhsOps(
 
      boost::shared_ptr<ForcesAndSourcesCore> &fe_contact_tree
 
  );
 
  MoFEMErrorCode setElasticElementOps(const int tag);
  MoFEMErrorCode setElasticElementToTs(DM dm);
 
  /** \brief Add debug to model
   *
   * That prints information every SNES step
   */
  MoFEMErrorCode addDebugModel(TS ts);
 
  MoFEMErrorCode solveElastic(TS ts, Vec x);
 
  /**
   * @brief Solve problem using dynamic relaxation method
   *
   * @param ts solver time stepper
   * @param x solution vector
   * @param start_step starting step number
   * @param start_time starting time
   * @return MoFEMErrorCode
   */
  MoFEMErrorCode solveDynamicRelaxation(TS ts, Vec x, int start_step,
                                        double start_time);
 
  /**
   * @brief Solve the incremental constitutive optimization with TAO
   *
   * @param ts
   * @param x
   * @return * MoFEMErrorCode
   */
  MoFEMErrorCode solveIncrementalOptimizationTAO(TS ts, Vec x, int start_step,
                                                 double start_time);
 
protected:
  MoFEMErrorCode runIncrementalOptimizationTAO(TS ts, Vec x);
 
public:
  /**
   * @brief   Solve load factor crack growth problem
   *
   * @param ts
   * @param x
   * @return * MoFEMErrorCode
   */
  MoFEMErrorCode solveLoadFactor(TS ts, Vec x, int start_step,
                                 double start_time);
 
  /**
   * @brief   Solve shape optimisation problem
   *
   * @param ts
   * @param x
   * @return * MoFEMErrorCode
   */
  MoFEMErrorCode solveSchapeOptimisation(TS ts, Vec x, int start_step,
                                         double start_time);
 
  MoFEMErrorCode solveTestTopologicalDerivative(TS ts, Vec x, int start_step,
                                                double start_time);
 
  MoFEMErrorCode solveTestEquilibratedMechanicalValue(TS ts, Vec x,
                                                      int start_step,
                                                      double start_time);
 
  MoFEMErrorCode solveTestIncrementalOptimizationLayout(TS ts, Vec x,
                                                        int start_step,
                                                        double start_time);
 
  MoFEMErrorCode solveTestIncrementalOptimizationTransaction(TS ts, Vec x,
                                                             int start_step,
                                                             double start_time);
 
  MoFEMErrorCode solveTestIncrementalOptimizationObjectiveDerivative(
      TS ts, Vec x, int start_step, double start_time);
 
  MoFEMErrorCode solveTestIncrementalOptimizationConstraintDerivative(
      TS ts, Vec x, int start_step, double start_time);
 
  MoFEMErrorCode setBlockTagsOnSkin();
 
  MoFEMErrorCode postProcessRestartMesh(const int tag, const std::string file,
                                        std::vector<Tag> tags_to_transfer = {});
 
  MoFEMErrorCode postProcessResults(const int tag, const std::string file,
                                    Vec f_residual = PETSC_NULLPTR,
                                    Vec var_vec = PETSC_NULLPTR,
                                    Vec gradient = PETSC_NULLPTR,
                                    std::vector<Tag> tags_to_transfer = {},
                                    TS ts = PETSC_NULLPTR);
  MoFEMErrorCode postProcessSkeletonResults(
      const int tag, const std::string file, Vec f_residual = PETSC_NULLPTR,
      std::vector<Tag> tags_to_transfer = {}, TS ts = PETSC_NULLPTR);
 
  MoFEMErrorCode calculateCrackArea(
      boost::shared_ptr<double> &area_ptr); // calculate crack area
 
  MoFEMErrorCode gettingNorms();
 
  struct SetUpSchur {
    static boost::shared_ptr<SetUpSchur> createSetUpSchur(
 
        MoFEM::Interface &m_field, EshelbianCore *ep_core_ptr
 
    );
    virtual MoFEMErrorCode setUp(TS) = 0;
 
  protected:
    SetUpSchur() = default;
  };
 
  struct DynamicRelaxationTimeScale : public TimeScale {
 
    using TimeScale::TimeScale;
 
    double getScale(const double time) override {
      if (EshelbianCore::physicalTimeFlg)
        return TimeScale::getScale(EshelbianCore::currentPhysicalTime);
      else
        return TimeScale::getScale(time);
    }
  };
 
  MoFEMErrorCode calculateFaceMaterialForce(
      const int tag, TS ts,
      SmartPetscObj<Vec> *adjoint_gradient_vector = nullptr);
  MoFEMErrorCode calculateOrientation(const int tag, bool set_orientation);
  MoFEMErrorCode setNewFrontCoordinates();
  MoFEMErrorCode addCrackSurfaces(const bool debug = false);
  MoFEMErrorCode createCrackSurfaceMeshset();
 
  boost::shared_ptr<Range> contactFaces;
  boost::shared_ptr<Range> crackFaces;
  boost::shared_ptr<Range> frontEdges;
  boost::shared_ptr<Range> frontAdjEdges;
  boost::shared_ptr<Range> frontVertices;
  boost::shared_ptr<Range> skeletonFaces;
  boost::shared_ptr<Range> maxMovedFaces;
  boost::shared_ptr<Range> interfaceFaces;
  boost::shared_ptr<Range> plasticVolumes;
 
  boost::shared_ptr<ParentFiniteElementAdjacencyFunctionSkeleton<2>>
      parentAdjSkeletonFunctionDim2;
 
  BitRefLevel bitAdjParent = BitRefLevel().set(); ///< bit ref level for parent
  BitRefLevel bitAdjParentMask =
      BitRefLevel().set(); ///< bit ref level for parent parent
  BitRefLevel bitAdjEnt = BitRefLevel().set(); ///< bit ref level for parent
  BitRefLevel bitAdjEntMask =
      BitRefLevel().set(); ///< bit ref level for parent parent
 
  SmartPetscObj<Vec> solTSStep;
  PetscBool loadFactorTSSolveExecuted = PETSC_FALSE;
 
  CommInterface::EntitiesPetscVector volumeExchange;
  CommInterface::EntitiesPetscVector faceExchange;
  CommInterface::EntitiesPetscVector edgeExchange;
  CommInterface::EntitiesPetscVector vertexExchange;
 
  std::vector<Tag>
      listTagsToTransfer; ///< list of tags to transfer to postprocessor
 
  Mat S = PETSC_NULLPTR;                 //< Schur complement matrix
  AO aoS = PETSC_NULLPTR;                //< AO for Schur complement matrix
  SmartPetscObj<IS> crackHybridIs;       //< IS for crack hybrid field
  std::vector<std::string> a00FieldList; //< list of fields for Schur complement
  std::vector<boost::shared_ptr<Range>>
      a00RangeList; //< list of ranges for Schur complement
 
  int nbCrackFaces = 0; //< number of crack faces
};
 
template <> struct EshelbianCore::FieldOrders<DEMKOWICZ_JACOBI_BASE> {
  static inline int stress(const int o) { return o; }
  static inline int bubble(const int o) { return o; }
  static inline int disp(const int o) { return o - 1; }
  static inline int rot(const int o) { return o - 1; }
  static inline int stretch(const int o) { return o; }
  static inline int logDeviator(const int o) { return stretch(o); }
  static inline int logJacobian(const int o) { return stretch(o); }
  static inline int auxiliaryLogStress(const int o) { return stretch(o); }
  static inline int hybrid(const int o) { return o - 1; }
};
 
template <> struct EshelbianCore::FieldOrders<AINSWORTH_LEGENDRE_BASE> {
  static inline int stress(const int o) { return o; }
  static inline int bubble(const int o) { return o + 1; }
  static inline int disp(const int o) { return o - 1; }
  static inline int rot(const int o) { return o; }
  static inline int stretch(const int o) { return o + 1; }
  static inline int logDeviator(const int o) { return stretch(o); }
  static inline int logJacobian(const int o) { return stretch(o); }
  static inline int auxiliaryLogStress(const int o) { return stretch(o); }
  static inline int hybrid(const int o) { return o; }
};