CrackPropagation.hpp

Go to the documentation of this file.

/** \file CrackPropagation.hpp
  \brief Main class for crack propagation
*/
 
/* 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/>. */
 
#ifndef __CRACK_PROPAGATION_HPP__
#define __CRACK_PROPAGATION_HPP__
 
namespace FractureMechanics {
 
MoFEMErrorCode clean_pcomms(moab::Interface &moab,
                            boost::shared_ptr<WrapMPIComm> moab_comm_wrap);
 
MoFEMErrorCode broadcast_entities(moab::Interface &moab,
                                  moab::Interface &moab_tmp,
                                  const int from_proc, Range &entities,
                                  const bool adjacencies, const bool tags);
 
/**
 * @brief Tapes numbers used by ADOL-C.
 *
 */
enum AdolcTags {
  ELASTIC_TAG = 1,
  MATERIAL_TAG,
  GRIFFITH_FORCE_TAG,
  GRIFFITH_CONSTRAINS_TAG,
  CONSTANT_AREA_TAG,
  ARC_LENGTH_TAG,
  FRONT_TANGENT,
  EXTERIOR_DERIVATIVE_TAG,
  SMOOTHING_TAG,
  SURFACE_SLIDING_TAG,
  EDGE_SLIDING_TAG
};
 
/**
 * @brief Names of tangent matrices tests
 *
 */
enum tangent_tests {
  MWLS_STRESS_TAN,
  MWLS_DENSITY_TAN,
  MWLS_GRIFFITH_TAN,
  LASTTAN
};
/**
 * Crack propagation data and methods
 * \nosubgrouping
 *
 * - Functions strarting with build, for example buildElastic,
 * add and build fields.
 *
 * - Functions starting with declare, for example declareElasticFE,
 * add and build finite elements.
 *
 * - Functions starting with assemble, for example assembleElasticDM,
 * create finite element instances and operators.
 *
 * - Functions starting with create, for example createElasticDM, create
 * problems.
 *
 */
struct CrackPropagation : public MoFEM::UnknownInterface {
 
  static bool debug;
 
  static bool parallelReadAndBroadcast; ///< That is unstable development, for
                                        ///< some meshses (propably generated by
                                        ///< cubit) this is not working. Error
                                        ///< can be attributed to bug in MOAB.
 
  Version fileVersion;
 
  MoFEMErrorCode getInterfaceVersion(Version &version) const {
    MoFEMFunctionBeginHot;
    version = Version(FM_VERSION_MAJOR, FM_VERSION_MINOR, FM_VERSION_BUILD);
    MoFEMFunctionReturnHot(0);
  }
 
  /**
   * \brief Getting interface of core database
   * @param  uuid  unique ID of interface that can be either for
   * CrackPropagation, CPSolvers or CPMeshCut interface
   * @param  iface returned pointer to interface
   * @return       error code
   */
  MoFEMErrorCode query_interface(boost::typeindex::type_index type_index,
                                 UnknownInterface **iface) const;
 
  boost::shared_ptr<WrapMPIComm> moabCommWorld;
 
  MoFEM::Interface &mField;
  bool setSingularCoordinates;
  int approxOrder;
  int geometryOrder;
 
  double gC;        ///< Griffith energy
  double griffithE; ///< Griffith stability parameter
  double griffithR; ///< Griffith regularisation parameter
 
  double betaGc;               ///< heterogeneous Griffith energy exponent
  PetscBool isGcHeterogeneous; ///< flag for heterogeneous gc
 
  PetscBool isPartitioned;
  PetscBool propagateCrack; ///< If true crack propagation is calculated
  int refAtCrackTip;
  int refOrderAtTip;
  PetscBool addSingularity;
  PetscBool otherSideConstrains;
 
  PetscBool isPressureAle; ///< If true surface pressure is considered in ALE
  PetscBool areSpringsAle; ///< If true surface spring is considered in ALE
  PetscBool
      isSurfaceForceAle; ///< If true surface pressure is considered in ALE
  PetscBool
      ignoreMaterialForce; ///< If true surface pressure is considered in ALE
 
  std::string mwlsApproxFile;         ///< Name of file with internal stresses
  std::string mwlsStressTagName;      ///< Name of tag with internal stresses
  std::string mwlsEigenStressTagName; ///< Name of tag with eigen stresses
  std::string mwlsRhoTagName;         ///< Name of tag with density
  int residualStressBlock; ///< Block on which residual stress is applied
  int densityMapBlock;     ///< Block on which density is applied (Note: if used
                           ///< with internal stress, it has to be the same ID!)
  std::string configFile;
  bool isStressTagSavedOnNodes;
 
  double
      partitioningWeightPower; ///< Weight controling load balance. Elements at
                               ///< the crack front are higher order, also have
                               ///< more fildes, associated with ALE
                               ///< formulation. Higher number put more weight
                               ///< on those elements while partitioning.
 
  PetscBool resetMWLSCoeffsEveryPropagationStep; ///< If true, MWLS coefficients
                                                 ///< are recalulated every
                                                 ///< propagaton step.
 
  static FTensor::Tensor2_symmetric<double, 3> analyticalStrainFunction(
 
      FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3> &t_coords
 
  );
 
  PetscBool
      addAnalyticalInternalStressOperators; ///< Add analytical internal stress
                                            ///< operators for MWLS test
 
  PetscBool solveEigenStressProblem;        ///< Solve eigen problem
  PetscBool useEigenPositionsSimpleContact; ///< Use eigen positions for
                                            ///< matching meshes contact
 
  int postProcLevel; ///< level of postprocessing (amount of output files)
  int startStep;
 
  double crackFrontLength;
  double crackSurfaceArea;
 
  /**
   * Map associating string literals to bit refinement levels
   * The following convention is currently used:
   * "mesh_cut"        : level 0
   * "spatial_domain"  : level 1
   * "material_domain" : level 2
   * */
  std::map<std::string, BitRefLevel> mapBitLevel;
 
  /**
   * \brief Constructor of Crack Propagation, prints fracture module version and
   registers the three interfaces
   * used to solve fracture problem, i.e. CrackPropagation, CPSolvers and
   * CPMeshCut
   * @param  m_field interface to MoAB database
   * @param  approx_order order of approximation space, default is 2
   * @param  geometry_order order of geometry, default is 1 (can be up to 2)
   */
  CrackPropagation(MoFEM::Interface &m_field, const int approx_order = 2,
                   const int geometry_order = 1);
 
  virtual ~CrackPropagation();
 
  /**
   * \brief Get options form command line
   * @return error code
   */
  MoFEMErrorCode getOptions();
 
  /**
   * \brief This is run with ctest
   * @return error code
   */
  MoFEMErrorCode tetsingReleaseEnergyCalculation();
 
  inline int &getNbLoadSteps() { return nbLoadSteps; }
  inline int getNbLoadSteps() const { return nbLoadSteps; }
  inline double &getLoadScale() { return loadScale; }
  inline double getLoadScale() const { return loadScale; }
  inline int &getNbCutSteps() { return nbCutSteps; }
  inline int getNbCutSteps() const { return nbCutSteps; }
 
  /// \brief read mesh file
  MoFEMErrorCode readMedFile();
 
  /**
   * @brief Save entities on ech processor.
   *
   * Usually used for debugging, to check if meshes are consistent on all
   * processors.
   *
   * @param prefix file name prefix
   * @param ents entities to save
   * @return MoFEMErrorCode error code
   */
  MoFEMErrorCode saveEachPart(const std::string prefix, const Range &ents);
 
  /** \name  Resolve shared entities and partition mesh */
 
  /**@{*/
 
  /// \brief partotion mesh
  MoFEMErrorCode partitionMesh(BitRefLevel bit1, BitRefLevel bit2,
                               int verb = QUIET, const bool debug = false);
 
  /// \brief resolve shared entities
  MoFEMErrorCode resolveShared(const Range &tets, Range &proc_ents,
                               const int verb = QUIET,
                               const bool debug = false);
 
  /// \brief resole shared entities by bit level
  MoFEMErrorCode resolveSharedBitRefLevel(const BitRefLevel bit,
                                          const int verb = QUIET,
                                          const bool debug = false);
 
  /**@}*/
 
  /** \name  Set bit level for material problem */
 
  /**@{*/
 
  /**
   * @brief Set bit ref level for entities adjacent to crack front
   *
   * @param from_bit
   * @param bit
   * @param nb_levels adjacency bridge levels
   * @return MoFEMErrorCode
   */
  MoFEMErrorCode setCrackFrontBitLevel(BitRefLevel from_bit, BitRefLevel bit,
                                       const int nb_levels = 2,
                                       const bool debug = false);
 
  /**@}*/
 
  /** \name Build problem */
 
  /**@{*/
 
  /**
   * \brief Build problem fields
   * @param  bit1  bit level for computational domain
   * @param  mask1 bit mask
   * @param  bit2  bit level for elements near crack front
   * @param  mask2 bit mask 2
   * @param  verb  verbosity level
   * @param  debug
   * @return       error code
   *
   */
  MoFEMErrorCode buildProblemFields(const BitRefLevel &bit1,
                                    const BitRefLevel &mask1,
                                    const BitRefLevel &bit2,
                                    const int verb = QUIET,
                                    const bool debug = false);
 
  /**
   * \brief Build problem finite elements
   * @param  bit1  bit level for computational domain
   * @param  bit2  bit level for elements near crack front
   * @param  surface_ids meshsets Ids on body skin
   * @param  verb  verbosity level
   * @param  debug
   * @return       error code
   *
   */
  MoFEMErrorCode buildProblemFiniteElements(BitRefLevel bit1, BitRefLevel bit2,
                                            std::vector<int> surface_ids,
                                            const int verb = QUIET,
                                            const bool debug = false);
 
  /**
   * \brief Construct problem data structures
   * @param  surface_ids surface_ids meshsets Ids on body skin
   * @param  verb        verbosity level
   * @param  debug
   * @return             error code
   */
  MoFEMErrorCode createProblemDataStructures(const std::vector<int> surface_ids,
                                             const int verb = QUIET,
                                             const bool debug = false);
 
  /**
   * \brief Crate DMs for all problems and sub problems
   * @param  dm_elastic            elastic DM
   * @param  dm_eiegn_elastci      eigen_elastic DM to solve proble wih eigen
   * stresses
   * @param  dm_material           material DM
   * @param  dm_crack_propagation  composition of elastic DM and material DM
   * @param  dm_material_forces    used to calculate material forces, sub-dm of
   * dm_crack_propagation
   * @param  dm_surface_projection dm to project material forces on surfaces,
   * sub-dm of dm_crack_propagation
   * @param  dm_crack_srf_area     dm to project material forces on crack
   * surfaces, sub-dm of dm_crack_propagation
   * @param  surface_ids           IDs of surface meshsets
   * @param  fe_surf_proj_list     name of elements on surface elements
   * @return                       error code
   */
  MoFEMErrorCode createDMs(SmartPetscObj<DM> &dm_elastic,
                           SmartPetscObj<DM> &dm_eigen_elastic,
                           SmartPetscObj<DM> &dm_material,
                           SmartPetscObj<DM> &dm_crack_propagation,
                           SmartPetscObj<DM> &dm_material_forces,
                           SmartPetscObj<DM> &dm_surface_projection,
                           SmartPetscObj<DM> &dm_crack_srf_area,
                           std::vector<int> surface_ids,
                           std::vector<std::string> fe_surf_proj_list);
 
  MoFEMErrorCode deleteEntities(const int verb = QUIET,
                                const bool debug = false);
 
  /**@}*/
 
  /** \name Build fields */
 
  /**@{*/
 
  /**
   * \brief Declate fields for elastic analysis
   * @param  bit bit refinement level
   * @return     error code
   */
  MoFEMErrorCode buildElasticFields(
      const BitRefLevel bit, const BitRefLevel mask = BitRefLevel().set(),
      const bool proc_only = true, const bool build_fields = true,
      const int verb = QUIET, const bool debug = false);
 
  /**
   * \brief Declate field for arc-length
   * @param  bit bit refinement level
   * @return     error code
   */
  MoFEMErrorCode buildArcLengthField(const BitRefLevel bit,
                                     const bool build_fields = true,
                                     const int verb = QUIET);
 
  /// \brief build fields with Lagrange multipliers to constrain surfaces
  MoFEMErrorCode buildSurfaceFields(const BitRefLevel bit,
                                    const bool proc_only = true,
                                    const bool build_fields = true,
                                    const int verb = QUIET,
                                    const bool debug = false);
 
  /// \brief build fields with Lagrange multipliers to constrain edges
  MoFEMErrorCode buildEdgeFields(const BitRefLevel bit,
                                 const bool proc_only = true,
                                 const bool build_fields = true,
                                 const int verb = QUIET,
                                 const bool debug = false);
 
  /// \brief declare crack surface files
  MoFEMErrorCode buildCrackSurfaceFieldId(const BitRefLevel bit,
                                          const bool proc_only = true,
                                          const bool build_fields = true,
                                          const int verb = QUIET,
                                          const bool debug = false);
 
  /// \brief declare crack surface files
  MoFEMErrorCode buildCrackFrontFieldId(const BitRefLevel bit,
                                        const bool build_fields = true,
                                        const int verb = QUIET,
                                        const bool debug = false);
 
  /**
   * \brief Lagrange multipliers field which constrains material displacements
   *
   * Nodes on both sides are constrainded to have the same material
   * displacements
   *
   */
  MoFEMErrorCode buildBothSidesFieldId(const BitRefLevel bit_spatial,
                                       const BitRefLevel bit_material,
                                       const bool proc_only = false,
                                       const bool build_fields = true,
                                       const int verb = QUIET,
                                       const bool debug = false);
 
  /** \brief Zero fields with lagrange multipliers
   */
  MoFEMErrorCode zeroLambdaFields();
 
  /**@}*/
 
  /** \name  Build finite elements */
 
  /**@{*/
 
  /// \brief declare elastic finite elements
  MoFEMErrorCode
  declareElasticFE(const BitRefLevel bit1, const BitRefLevel mask1,
                   const BitRefLevel bit2, const BitRefLevel mask2,
                   const bool add_forces = true, const bool proc_only = true,
                   const int verb = QUIET);
 
  /// \brief create arc-length element entity and declare elemets
  MoFEMErrorCode declareArcLengthFE(const BitRefLevel bits,
                                    const int verb = QUIET);
 
  MoFEMErrorCode
  declareExternalForcesFE(const BitRefLevel bit,
                          const BitRefLevel mask = BitRefLevel().set(),
                          const bool proc_only = true);
 
  /** \brief Declare FE for pressure BC in ALE formulation (in material domain)
   * @param  bit            material domain bit ref level
   * @param  mask           bit ref level mask for the problem
   * @return                error code
   */
  MoFEMErrorCode
  declareSurfaceForceAleFE(const BitRefLevel bit,
                           const BitRefLevel mask = BitRefLevel().set(),
                           const bool proc_only = true);
 
  /** \brief Declare FE for pressure BC in ALE formulation (in material domain)
   * @param  bit            material domain bit ref level
   * @param  mask           bit ref level mask for the problem
   * @return                error code
   */
  MoFEMErrorCode
  declarePressureAleFE(const BitRefLevel bit,
                       const BitRefLevel mask = BitRefLevel().set(),
                       const bool proc_only = true);
 
  /** \brief Declare FE for spring BC in ALE formulation (in material domain)
   * @param  bit            material domain bit ref level
   * @param  mask           bit ref level mask for the problem
   * @return                error code
   */
  MoFEMErrorCode
  declareSpringsAleFE(const BitRefLevel bit,
                      const BitRefLevel mask = BitRefLevel().set(),
                      const bool proc_only = true);
 
  /** \brief Declare FE for pressure BC in ALE formulation (in material domain)
   * @param  bit            material domain bit ref level
   * @param  mask           bit ref level mask for the problem
   * @return                error code
   */
  MoFEMErrorCode
  declareSimpleContactAleFE(const BitRefLevel bit,
                            const BitRefLevel mask = BitRefLevel().set(),
                            const bool proc_only = true);
 
  /// \brief declare material finite elements
  MoFEMErrorCode declareMaterialFE(const BitRefLevel bit,
                                   const BitRefLevel mask = BitRefLevel().set(),
                                   const bool proc_only = true,
                                   const bool verb = QUIET);
 
  MoFEMErrorCode declareFrontFE(const BitRefLevel bit,
                                const BitRefLevel mask = BitRefLevel().set(),
                                const bool proc_only = true,
                                const bool verb = QUIET);
 
  MoFEMErrorCode
  declareBothSidesFE(const BitRefLevel bit_spatial,
                     const BitRefLevel bit_material,
                     const BitRefLevel mask = BitRefLevel().set(),
                     const bool proc_only = true, const bool verb = QUIET);
 
  /// \brief declare mesh smoothing finite elements
  MoFEMErrorCode
  declareSmoothingFE(const BitRefLevel bit,
                     const BitRefLevel mask = BitRefLevel().set(),
                     const bool proc_only = true, const bool verb = QUIET);
 
  /// \brief declare surface sliding elements
  MoFEMErrorCode declareSurfaceFE(std::string fe_name, const BitRefLevel bit,
                                  const BitRefLevel mask,
                                  const std::vector<int> &ids,
                                  const bool proc_only = true,
                                  const int verb = QUIET,
                                  const bool debug = false);
 
  // \brief declare edge sliding elements
  MoFEMErrorCode declareEdgeFE(std::string fe_name, const BitRefLevel bit,
                               const BitRefLevel mask,
                               const bool proc_only = true,
                               const int verb = QUIET,
                               const bool debug = false);
 
  MoFEMErrorCode
  declareCrackSurfaceFE(std::string fe_name, const BitRefLevel bit,
                        const BitRefLevel mask, const bool proc_only = true,
                        const int verb = QUIET, const bool debug = false);
 
  /**@}*/
 
  /** \name  Create DMs & problems */
 
  /**@{*/
 
  /** \brief Create DM by composition of elastic DM and material DM
   * @param  dm                    composition of elastic DM and material DM
   * @param  prb_name              problem name
   * @param  dm_elastic            elastic DM
   * @param  dm_material           material DM
   * @param  bit                   domain bit ref level
   * @param  mask                  mask ref level for problem
   * @return                       error code
   */
  MoFEMErrorCode
  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);
 
  /** \brief Create elastic problem DM
   * @param  dm                    elastic DM
   * @param  prb_name              problem name
   * @param  bit                   domain bit ref level
   * @param  mask                  mask ref level for problem
   * @return                       error code
   */
  MoFEMErrorCode createElasticDM(SmartPetscObj<DM> &dm,
                                 const std::string prb_name,
                                 const BitRefLevel bit,
                                 const BitRefLevel mask = BitRefLevel().set());
 
  /** \brief Create elastic problem DM
   * @param  dm                    elastic DM
   * @param  prb_name              problem name
   * @param  bit                   domain bit ref level
   * @param  mask                  mask ref level for problem
   * @return                       error code
   */
  MoFEMErrorCode
  createEigenElasticDM(SmartPetscObj<DM> &dm, const std::string prb_name,
                       const BitRefLevel bit,
                       const BitRefLevel mask = BitRefLevel().set());
 
  /** \brief Create problem to calculate boundary conditions
   * @param  dm                    newly created boundary condition DM
   * @param  prb_name              problem name
   * @param  bit                   auxiliary bit ref level
   * @param  mask                  mask ref level for problem
   * @return                       error code
   * \note BcDM is only useed for testing, to enforce boundary conditions
staisifing anlylitical solution to which numerical solution is compared.
   */
  MoFEMErrorCode createBcDM(SmartPetscObj<DM> &dm, const std::string prb_name,
                            const BitRefLevel bit,
                            const BitRefLevel mask = BitRefLevel().set());
 
  /** \brief Create DM fto calculate material problem
   * @param  dm                    used to calculate material forces, sub-dm of
   * dm_crack_propagation
   * @param  prb_name              problem name
   * @param  bit                   material domain
   * @param  mask                  dof mask ref level for problem
   * @param  fe_list               name of elements on surface elements
   * @param  debug                 flag for debugging
   * @return                       error code
   */
  MoFEMErrorCode 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 = false);
 
  /** \brief Create DM for calculation of material forces (sub DM of DM
   * material)
   * @param  dm                    used to calculate material forces, sub-dm of
   * dm_crack_propagation
   * @param  prb_name              problem name
   * @param  verb                  compilation parameter determining the amount
   * of information printed
   * @return                       error code
   */
  MoFEMErrorCode createMaterialForcesDM(SmartPetscObj<DM> &dm,
                                        SmartPetscObj<DM> dm_material,
                                        const std::string prb_name,
                                        const int verb = QUIET);
 
  /** \brief create DM to calculate projection matrices (sub DM of DM
   * material)
   * @param  dm                    DM to project material forces on surfaces,
   * sub-dm of dm_crack_propagation
   * @param  dm_material           material DM
   * @param  prb_name              problem name
   * @param  surface_ids           IDs of surface meshsets
   * @param  fe_list               name of elements on surface elements
   * @param  verb                  compilation parameter determining the amount
   * of information printed
   * @return                       error code
   */
  MoFEMErrorCode createSurfaceProjectionDM(
      SmartPetscObj<DM> &dm, SmartPetscObj<DM> dm_material,
      const std::string prb_name, const std::vector<int> surface_ids,
      const std::vector<std::string> fe_list, const int verb = QUIET);
 
  /** \brief create DM to calculate Griffith energy
   * @param  dm                    dm to project material forces on crack
   * surfaces, sub-dm of dm_crack_propagation
   * @param  dm_material           material DM
   * @param  prb_name              problem name
   * @param  verb                  compilation parameter determining the amount
   * of information printed
   * @return                       error code
   */
  MoFEMErrorCode createCrackFrontAreaDM(SmartPetscObj<DM> &dm,
                                        SmartPetscObj<DM> dm_material,
                                        const std::string prb_name,
                                        const bool verb = QUIET);
 
  /**@}*/
 
  /** \name  Create finite element instances */
 
  /**@{*/
 
  /** \brief create elastic finite element instance for spatial assembly
   * @param  mwls_stress_tag_name  Name of internal stress tag
   * @param  close_crack           if true, crack surface is closed
   * @param  verb                  compilation parameter determining the amount
   * of information printed (not in use here)
   * @param  debug                 flag for debugging
   * @return                       error code
   */
  MoFEMErrorCode assembleElasticDM(const std::string mwls_stress_tag_name,
                                   const int verb = QUIET,
                                   const bool debug = false);
 
  /** \brief create elastic finite element instance for spatial assembly
   * @param  verb                   compilation parameter determining the amount
   * of information printed (not in use here)
   * @param  debug                 flag for debugging
   * @return                       error code
   */
  MoFEMErrorCode assembleElasticDM(const int verb = QUIET,
                                   const bool debug = false);
 
  // \brief add elastic element instances to SNES
  MoFEMErrorCode addElasticFEInstancesToSnes(
      DM dm, Mat m, Vec q, Vec f,
      boost::shared_ptr<FEMethod> arc_method = boost::shared_ptr<FEMethod>(),
      boost::shared_ptr<ArcLengthCtx> arc_ctx = nullptr, const int verb = QUIET,
      const bool debug = false);
 
  /// \brief create material element instance
  MoFEMErrorCode assembleMaterialForcesDM(DM dm, const int verb = QUIET,
                                          const bool debug = false);
 
  /// \brief create smoothing element instance
  MoFEMErrorCode assembleSmootherForcesDM(DM dm, const std::vector<int> ids,
                                          const int verb = QUIET,
                                          const bool debug = false);
 
  /// \brief assemble coupling element instances
  MoFEMErrorCode assembleCouplingForcesDM(DM dm, const int verb = QUIET,
                                          const bool debug = false);
 
  /// \brief add material elements instances to SNES
  MoFEMErrorCode addMaterialFEInstancesToSnes(DM dm, const bool fix_crack_front,
                                              const int verb = QUIET,
                                              const bool debug = false);
 
  /// \brief add softening elements instances to SNES
  MoFEMErrorCode addSmoothingFEInstancesToSnes(DM dm,
                                               const bool fix_crack_front,
                                               const int verb = QUIET,
                                               const bool debug = false);
 
  /// \brief add finite element to SNES for crack propagation problem
  MoFEMErrorCode
  addPropagationFEInstancesToSnes(DM dm, boost::shared_ptr<FEMethod> arc_method,
                                  boost::shared_ptr<ArcLengthCtx> arc_ctx,
                                  const std::vector<int> &surface_ids,
                                  const int verb = QUIET,
                                  const bool debug = false);
 
  /**  \brief test LHS Jacobians
   * @param  bit                     bit level
   * @param  mask                    bitref mask
   * @param  test_nb                 name of the test
   */
  MoFEMErrorCode testJacobians(const BitRefLevel bit, const BitRefLevel mask,
                               tangent_tests test);
 
  MoFEMErrorCode
  addMWLSStressOperators(boost::shared_ptr<CrackFrontElement> &fe_rhs,
                         boost::shared_ptr<CrackFrontElement> &fe_lhs);
 
  MoFEMErrorCode
  addMWLSDensityOperators(boost::shared_ptr<CrackFrontElement> &fe_rhs,
                          boost::shared_ptr<CrackFrontElement> &fe_lhs);
 
  /// \brief Update fixed nodes
  MoFEMErrorCode updateMaterialFixedNode(const bool fix_front,
                                         const bool fix_small_g,
                                         const bool debug = false);
 
  /**@}*/
 
  /** \name  Assemble and calculate */
 
  /**@{*/
 
  /// \brief assemble elastic part of matrix, by running elastic finite element
  /// instance
  MoFEMErrorCode calculateElasticEnergy(DM dm, const std::string msg = "");
 
  /** \brief assemble material forces, by running material finite element
   * instance
   * @param  dm                    dm used to calculate material forces, sub-dm
   * of dm_crack_propagation
   * @param  q                     DoF sub-vector for crack propagation DM
   * @param  f                     Rhs sub-vector for crack propagation DM
   * @param  verb                   compilation parameter determining the amount
   * of information printed
   * @param  debug                 flag for debugging
   * @return                       error code
   */
  MoFEMErrorCode calculateMaterialForcesDM(DM dm, Vec q, Vec f,
                                           const int verb = QUIET,
                                           const bool debug = false);
 
  /** \brief assemble smoothing forces, by running material finite element
   * instance
   * @param  dm                    crack propagation dm that is composition of
   * elastic DM and material DM
   * @param  q                     DoF sub-vector for crack propagation DM
   * @param  f                     Rhs sub-vector for crack propagation DM
   * @param  verb                  compilation parameter determining the amount
   * of information printed
   * @param  debug                 flag for debugging
   * @return                       error code
   */
  MoFEMErrorCode calculateSmoothingForcesDM(DM dm, Vec q, Vec f,
                                            const int verb = QUIET,
                                            const bool debug = false);
 
  MoFEMErrorCode calculateSmoothingForceFactor(const int verb = QUIET,
                                               const bool debug = true);
 
  /// \brief assemble projection matrices
  MoFEMErrorCode calculateSurfaceProjectionMatrix(DM dm_front, DM dm_project,
                                                  const std::vector<int> &ids,
                                                  const int verb = QUIET,
                                                  const bool debug = false);
 
  /// \brief assemble crack front projection matrix (that constrains crack area
  /// growth)
  MoFEMErrorCode calculateFrontProjectionMatrix(DM dm_surface, DM dm_project,
                                                const int verb = QUIET,
                                                const bool debug = false);
 
  /** \brief project material forces along the crack elongation direction
   * @param  dm_project            material forces DM, sub-dm of
   * dm_crack_propagation
   * @param  f_proj                Rhs sub-vector of projected material forces
   * for material forces DM
   * @param  f                     Rhs sub-vector for material forces DM
   * @param  verb                  compilation parameter determining the amount
   * of information printed
   * @param  debug                 flag for debugging
   * @return                       error code
   */
  MoFEMErrorCode projectMaterialForcesDM(DM dm_project, Vec f, Vec f_proj,
                                         const int verb = QUIET,
                                         const bool debug = false);
 
  /// \brief calculate Griffith (driving) force
  MoFEMErrorCode calculateGriffithForce(DM dm, const double gc, Vec f_griffith,
                                        const int verb = QUIET,
                                        const bool debug = false);
 
  /// \brief project Griffith forces
  MoFEMErrorCode projectGriffithForce(DM dm, Vec f_griffith,
                                      Vec f_griffith_proj,
                                      const int verb = QUIET,
                                      const bool debug = false);
 
  /// \brief calculate release energy
  MoFEMErrorCode calculateReleaseEnergy(DM dm, Vec f_material_proj,
                                        Vec f_griffith_proj, Vec f_lambda,
                                        const double gc, const int verb = QUIET,
                                        const bool debug = true);
 
  /** \name  Solvers */
 
  /**@{*/
 
  /** \brief solve elastic problem
   * @param  dm                    elastic DM
   * @param  snes                  snes solver for elastic DM
   * @param  m                     stiffness sub-matrix for elastic DM
   * @param  q                     DoF sub-vector for elastic DM
   * @param  f                     Rhs sub-vector for material DM
   * @param  snes_set_up           boolean flag wheather to determin snes solver
   * @param  shell_m               stiffness sub-matrix for elastic DM
   * @return                       error code
   */
  MoFEMErrorCode solveElasticDM(DM dm, SNES snes, Mat m, Vec q, Vec f,
                                bool snes_set_up, Mat *shell_m);
 
  /** \brief solve crack propagation problem
   * @param  dm                    crack propagation DM that is composed of
   * elastic DM and material DM
   * @param  dm_elastic            elastic DM
   * @param  snes                  snes solver for crack propagation DM
   * @param  m                     stiffness sub-matrix for crack propagation DM
   * @param  q                     DoF sub-vector for crack propagation DM
   * @param  f                     Rhs sub-vector for crack propagation DM
   * @param  snes_set_up           boolean flag wheather to determin snes solver
   * @return                       error code
   */
  MoFEMErrorCode solvePropagationDM(DM dm, DM dm_elastic, SNES snes, Mat m,
                                    Vec q, Vec f);
 
  /**@}*/
 
  /**@}*/
 
  /** \name  Postprocess results */
 
  /**@{*/
 
  MoFEMErrorCode savePositionsOnCrackFrontDM(DM dm, Vec q,
                                             const int verb = QUIET,
                                             const bool debug = false);
 
  MoFEMErrorCode writeCrackFont(const BitRefLevel &bit, const int step);
 
  /// \brief post-process results for elastic solution
  MoFEMErrorCode postProcessDM(DM dm, const int step, const std::string fe_name,
                               const bool approx_internal_stress);
 
  /**@}*/
 
  /** \name  Set data from mesh to vectors and back */
 
  /**@{*/
 
  /// \brief set field from node positions
  MoFEMErrorCode setFieldFromCoords(const std::string field_name);
 
  /// \brief set material field from nodes
  MoFEMErrorCode setMaterialPositionFromCoords();
 
  /// \brief set spatial field from nodes
  MoFEMErrorCode setSpatialPositionFromCoords();
 
  /// \brief set coordinates from field
  MoFEMErrorCode
  setCoordsFromField(const std::string field_name = "MESH_NODE_POSITIONS");
 
  /// \brief set singular dofs (on edges adjacent to crack front) from geometry
  MoFEMErrorCode setSingularDofs(const string field_name,
                                 const int verb = QUIET);
 
  /// \brief set singular dofs of material field (on edges adjacent to crack
  /// front) from geometry
  MoFEMErrorCode setSingularElementMatrialPositions(const int verb = QUIET);
 
  /// \brief remove singularity
  MoFEMErrorCode unsetSingularElementMatrialPositions();
 
  /// \brief set maetrial field order to one
  MoFEMErrorCode cleanSingularElementMatrialPositions();
 
  /**@}*/
 
  /** \name  Auxiliary method */
 
  /**@{*/
 
  boost::shared_ptr<DofEntity> arcLengthDof;
 
  /// \brief set pointer to arc-length DOF
  MoFEMErrorCode getArcLengthDof();
 
  inline boost::shared_ptr<ArcLengthCtx> &getArcCtx() { return arcCtx; }
  inline boost::shared_ptr<ArcLengthCtx> &getEigenArcCtx() {
    return arcEigenCtx;
  }
 
  inline boost::shared_ptr<FEMethod>
  getFrontArcLengthControl(boost::shared_ptr<ArcLengthCtx> arc_ctx) {
    boost::shared_ptr<FEMethod> arc_method(
        new GriffithForceElement::FrontArcLengthControl(
            ARC_LENGTH_TAG, griffithForceElement->blockData[0],
            griffithForceElement->commonData, "LAMBDA_CRACKFRONT_AREA",
            arc_ctx));
    return arc_method;
  }
 
  inline boost::shared_ptr<AnalyticalDirichletBC::DirichletBC>
  getAnalyticalDirichletBc() {
    return analyticalDirichletBc;
  }
 
  boost::shared_ptr<MWLSApprox> &getMWLSApprox() { return mwlsApprox; }
 
  /**@}*/
 
  /**
   * \brief operator to post-process results on crack front nodes
   */
  struct PostProcVertexMethod : public MoFEM::DofMethod {
 
    MoFEM::Interface &mField;
    Vec F;
    std::string tagName;
    PostProcVertexMethod(MoFEM::Interface &m_field, Vec f, std::string tag_name)
        : mField(m_field), F(f), tagName(tag_name) {}
 
    Tag tH;
    double *aRray;
    MoFEMErrorCode preProcess();
 
    MoFEMErrorCode postProcess();
 
    MoFEMErrorCode operator()();
  };
 
  /**
   * \brief Constrains material displacement on both sides
   *
   * Both sides of crack surface have the same material displacements.
   *
   */
  struct BothSurfaceConstrains : public FEMethod {
 
    /**
     * @brief Construct a new Both Surface Constrains object
     *
     * @param m_field
     * @param lambda_field_name lagrange multiplayers
     * @param spatial_field_name constrained spatial field
     */
    BothSurfaceConstrains(
        MoFEM::Interface &m_field,
        const std::string lambda_field_name = "LAMBDA_BOTH_SIDES",
        const std::string spatial_field_name = "MESH_NODE_POSITIONS")
        : mField(m_field), lambdaFieldName(lambda_field_name),
          spatialFieldName(spatial_field_name), aLpha(1) {
      ierr = getOptions();
      CHKERRABORT(PETSC_COMM_WORLD, ierr);
    }
 
    MoFEMErrorCode getOptions() {
      MoFEMFunctionBegin;
      ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "",
                               "Get Booth surface constrains element scaling",
                               "none");
      CHKERRQ(ierr);
      CHKERR PetscOptionsScalar("-booth_side_alpha", "scaling parameter", "",
                                aLpha, &aLpha, PETSC_NULL);
      ierr = PetscOptionsEnd();
      CHKERRQ(ierr);
      MoFEMFunctionReturn(0);
    }
 
    MoFEMErrorCode preProcess() { return 0; }
    MoFEMErrorCode operator()();
    MoFEMErrorCode postProcess() { return 0; }
 
    double aLpha;      ///! Scaling parameter
    Range masterNodes; ///! Node on which lagrange multiplier is set
 
  private:
    MoFEM::Interface &mField;
    const std::string lambdaFieldName;
    const std::string spatialFieldName;
  };
 
  /// \brief Determine face orientation
  struct FaceOrientation
      : public SurfaceSlidingConstrains::DriverElementOrientation {
    bool useProjectionFromCrackFront;
    Range contactFaces;
    BitRefLevel defaultBitLevel;
    FaceOrientation(bool crack_front, Range contact_faces,
                    BitRefLevel bit_level)
        : useProjectionFromCrackFront(crack_front), contactFaces(contact_faces),
          defaultBitLevel(bit_level) {}
    virtual ~FaceOrientation() {}
    MoFEMErrorCode getElementOrientation(MoFEM::Interface &m_field,
                                         const EntityHandle face,
                                         const BitRefLevel bit);
    MoFEMErrorCode getElementOrientation(MoFEM::Interface &m_field,
                                         const FEMethod *fe_method_ptr) {
      MoFEMFunctionBegin;
      EntityHandle face = fe_method_ptr->numeredEntFiniteElementPtr->getEnt();
      BitRefLevel bit = fe_method_ptr->problemPtr->getBitRefLevel();
      CHKERR getElementOrientation(m_field, face, bit);
      MoFEMFunctionReturn(0);
    }
    MoFEMErrorCode getElementOrientation(MoFEM::Interface &m_field,
                                         const EntityHandle face) {
      MoFEMFunctionBegin;
      CHKERR getElementOrientation(m_field, face, defaultBitLevel);
      MoFEMFunctionReturn(0);
    }
  };
 
  struct ArcLengthSnesCtx : public SnesCtx {
 
    ArcLengthSnesCtx(MoFEM::Interface &m_field, const std::string problem_name,
                     boost::shared_ptr<ArcLengthCtx> &arc_ptr)
        : SnesCtx(m_field, problem_name), arcPtr(arc_ptr.get()) {
      MOFEM_LOG("CPWorld", Sev::noisy) << "ArcLengthSnesCtx create";
    }
    virtual ~ArcLengthSnesCtx() {
      MOFEM_LOG("CPWorld", Sev::noisy) << "ArcLengthSnesCtx destroyed";
    };
 
    inline auto getArcPtr() { return arcPtr; }
    inline void setVecDiagM(Vec diag_m) { diagM = diag_m; }
    inline auto getVecDiagM() { return diagM; }
 
  private:
    ArcLengthCtx *arcPtr;
    Vec diagM;
  };
 
  // private:
 
  /** \name  Pointers to finite element instances */
 
  /**@{*/
 
  bool onlyHooke;                 ///< True if only Hooke material is applied
  PetscBool onlyHookeFromOptions; ///< True if only Hooke material is applied
 
  /// Integrate spatial stresses, matrix and vector (not only element but all
  /// data structure)
  boost::shared_ptr<NonlinearElasticElement> elasticFe;
  boost::shared_ptr<NonlinearElasticElement> materialFe;
  boost::shared_ptr<CrackFrontElement> feLhs; ///< Integrate elastic FE
  boost::shared_ptr<CrackFrontElement> feRhs; ///< Integrate elastic FE
  boost::shared_ptr<CrackFrontElement>
      feMaterialRhs; ///< Integrate material stresses, assemble vector
  boost::shared_ptr<CrackFrontElement>
      feMaterialLhs; ///< Integrate material stresses, assemble matrix
  boost::shared_ptr<CrackFrontElement> feEnergy; ///< Integrate energy
  boost::shared_ptr<ConstrainMatrixCtx>
      projSurfaceCtx; ///< Data structure to project on the body surface
  boost::shared_ptr<ConstrainMatrixCtx>
      projFrontCtx; ///< Data structure to project on crack front
  boost::shared_ptr<MWLSApprox>
      mwlsApprox; ///< Data struture for MWLS approximation
  /// Surface elment to calculate tractions in material space
  boost::shared_ptr<NeumannForcesSurface::MyTriangleFE>
      feMaterialAnaliticalTraction;
  boost::shared_ptr<DirichletSpatialPositionsBc>
      spatialDirichletBc; ///< apply Dirichlet BC to sparial positions
  boost::shared_ptr<AnalyticalDirichletBC::DirichletBC>
      analyticalDirichletBc; ///< apply analytical Dirichlet BC to sparial
                             ///< positions
 
  boost::shared_ptr<boost::ptr_map<string, NeumannForcesSurface>>
      surfaceForces; ///< assemble surface forces
  boost::shared_ptr<boost::ptr_map<string, NeumannForcesSurface>>
      surfaceForceAle; ///< assemble surface pressure (ALE)
  boost::shared_ptr<NeumannForcesSurface::DataAtIntegrationPts>
      commonDataSurfaceForceAle; ///< common data at integration points (ALE)
  boost::shared_ptr<boost::ptr_map<string, NeumannForcesSurface>>
      surfacePressure; ///< assemble surface pressure
  boost::shared_ptr<boost::ptr_map<string, NeumannForcesSurface>>
      surfacePressureAle; ///< assemble surface pressure (ALE)
  boost::shared_ptr<NeumannForcesSurface::DataAtIntegrationPts>
      commonDataSurfacePressureAle; ///< common data at integration points (ALE)
  boost::shared_ptr<boost::ptr_map<string, EdgeForce>>
      edgeForces; ///< assemble edge forces
  boost::shared_ptr<boost::ptr_map<string, NodalForce>>
      nodalForces;                             ///< assemble nodal forces
  boost::shared_ptr<FEMethod> assembleFlambda; ///< assemble F_lambda vector
  boost::shared_ptr<FEMethod> zeroFlambda;     ///< assemble F_lambda vector
 
  boost::shared_ptr<CrackFrontElement>
      feCouplingElasticLhs; ///< FE instance to assemble coupling terms
  boost::shared_ptr<CrackFrontElement>
      feCouplingMaterialLhs; ///< FE instance to assemble coupling terms
 
  boost::shared_ptr<Smoother> smootherFe;
  boost::shared_ptr<Smoother::MyVolumeFE>
      feSmootherRhs; ///< Integrate smoothing operators
  boost::shared_ptr<Smoother::MyVolumeFE>
      feSmootherLhs; ///< Integrate smoothing operators
  boost::shared_ptr<VolumeLengthQuality<double>> volumeLengthDouble;
  boost::shared_ptr<VolumeLengthQuality<adouble>> volumeLengthAdouble;
 
  boost::shared_ptr<
      ObosleteUsersModules::TangentWithMeshSmoothingFrontConstrain>
      tangentConstrains; ///< Constrains crack front in tangent directiona
  boost::shared_ptr<SurfaceSlidingConstrains::DriverElementOrientation>
      skinOrientation;
  boost::shared_ptr<SurfaceSlidingConstrains::DriverElementOrientation>
      crackOrientation;
  boost::shared_ptr<SurfaceSlidingConstrains::DriverElementOrientation>
      contactOrientation;
  map<int, boost::shared_ptr<SurfaceSlidingConstrains>> surfaceConstrain;
  map<int, boost::shared_ptr<EdgeSlidingConstrains>> edgeConstrains;
 
  boost::shared_ptr<BothSurfaceConstrains> bothSidesConstrains;
  boost::shared_ptr<BothSurfaceConstrains> closeCrackConstrains;
  boost::shared_ptr<BothSurfaceConstrains> bothSidesContactConstrains;
  boost::shared_ptr<DirichletFixFieldAtEntitiesBc> fixMaterialEnts;
 
  boost::shared_ptr<GriffithForceElement> griffithForceElement;
  boost::shared_ptr<GriffithForceElement::MyTriangleFE> feGriffithForceRhs;
  boost::shared_ptr<GriffithForceElement::MyTriangleFE> feGriffithForceLhs;
  boost::shared_ptr<GriffithForceElement::MyTriangleFEConstrainsDelta>
      feGriffithConstrainsDelta;
  boost::shared_ptr<GriffithForceElement::MyTriangleFEConstrains>
      feGriffithConstrainsRhs;
  boost::shared_ptr<GriffithForceElement::MyTriangleFEConstrains>
      feGriffithConstrainsLhs;
 
  boost::shared_ptr<boost::ptr_map<string, NeumannForcesSurface>>
      analiticalSurfaceElement; ///< Finite element to integrate tractions given
                                ///< by analytical formula
  boost::shared_ptr<FaceElementForcesAndSourcesCore> feSpringLhsPtr;
  boost::shared_ptr<FaceElementForcesAndSourcesCore> feSpringRhsPtr;
  /**@}*/
 
  /** \name Entity ranges */
 
  /**@{*/
 
  Range bitEnts;
  Range bitProcEnts;
  Range bodySkin;
  Range crackFaces;
  Range oneSideCrackFaces;
  Range otherSideCrackFaces;
  Range crackFront;
  Range crackFrontNodes;
  Range crackFrontNodesEdges;
  Range crackFrontNodesTris;
  Range crackFrontElements;
 
  EntityHandle
      arcLengthVertex; ///< Vertex associated with dof for arc-length control
 
  /**@}*/
 
  /** \name  Auxiliary data */
 
  /**@{*/
 
  int defaultMaterial;
 
  int nbLoadSteps;
  int nbCutSteps;
  double loadScale;
 
  double maxG1;
  double maxJ;
  double fractionOfFixedNodes;
  map<EntityHandle, double> mapG1; ///< hashmap of g1 - release energy at nodes
  map<EntityHandle, double> mapG3; ///< hashmap of g3 - release energy at nodes
  map<EntityHandle, double> mapJ;  ///< hashmap of J - release energy at nodes
  map<EntityHandle, VectorDouble3>
      mapMatForce; ///< hashmap of material force at nodes
  map<EntityHandle, VectorDouble3>
      mapGriffith; ///< hashmap of Griffith energy at nodes
 
  map<EntityHandle, double> mapSmoothingForceFactor;
 
  bool doSurfaceProjection; ///< If crack surface is immersed in the body do not
                            ///< projection
 
  double smootherAlpha;        ///< Controls mesh smoothing
  double initialSmootherAlpha; ///< Controls mesh smoothing, set at start of
                               ///< analysis
  double smootherGamma;        ///< Controls mesh smoothing
 
  /**@}*/
 
  /** \name  arc-length options */
 
  /**@{*/
 
  double arcAlpha;
  double arcBeta;
  double arcS;
 
  boost::shared_ptr<ArcLengthCtx> arcCtx;
  boost::shared_ptr<ArcLengthCtx> arcEigenCtx;
 
  /**@}*/
 
  /** \name  Contact parameters, entity ranges and data structures */
 
  /**@{*/
 
  int contactOrder;       ///< Approximation order for spatial positions
                          ///< used in contact elements
  int contactLambdaOrder; ///< Approximation order for field of Lagrange
                          ///< multipliers used to enforce contact
  double rValue;          ///< Parameter for regularizing absolute value
                          ///< function in the complementarity function
  double cnValue; ///< Augmentation parameter in the complementarity function
                  ///< used to enforce KKT constraints
 
  PetscBool ignoreContact;   ///< If set true, contact interfaces are ignored
                             ///< and prism interfaces are not created
  PetscBool fixContactNodes; ///< If set true, contact nodes are fixed in
                             ///< the material configuration
  PetscBool printContactState;
 
  Range contactMeshsetFaces;
 
  Range contactElements;
  Range contactMasterFaces;
  Range contactSlaveFaces;
 
  Range mortarContactElements;
  Range mortarContactMasterFaces;
  Range mortarContactSlaveFaces;
 
  boost::shared_ptr<ContactSearchKdTree::ContactCommonData_multiIndex>
      contactSearchMultiIndexPtr;
  Range contactBothSidesMasterFaces;
  Range contactBothSidesSlaveFaces;
  Range contactBothSidesMasterNodes;
  Range contactTets;
 
  EntityHandle meshsetFaces;
 
  boost::shared_ptr<FaceElementForcesAndSourcesCore> feLhsSpringALE;
  boost::shared_ptr<FaceElementForcesAndSourcesCore> feLhsSpringALEMaterial;
  boost::shared_ptr<FaceElementForcesAndSourcesCore> feRhsSpringALEMaterial;
  boost::shared_ptr<FaceElementForcesAndSourcesCore>
      feLhsSurfaceForceALEMaterial;
  boost::shared_ptr<FaceElementForcesAndSourcesCore>
      feRhsSurfaceForceALEMaterial;
  boost::shared_ptr<MetaSpringBC::DataAtIntegrationPtsSprings>
      commonDataSpringsALE;
 
  boost::shared_ptr<SimpleContactProblem> contactProblem;
  boost::shared_ptr<SimpleContactProblem::SimpleContactElement>
      feRhsSimpleContact;
  boost::shared_ptr<SimpleContactProblem::SimpleContactElement>
      feLhsSimpleContact;
  boost::shared_ptr<SimpleContactProblem::CommonDataSimpleContact>
      commonDataSimpleContact;
 
  boost::shared_ptr<MortarContactProblem> mortarContactProblemPtr;
 
  boost::shared_ptr<MortarContactProblem::MortarContactElement>
      feRhsMortarContact;
  boost::shared_ptr<MortarContactProblem::MortarContactElement>
      feLhsMortarContact;
  boost::shared_ptr<MortarContactProblem::CommonDataMortarContact>
      commonDataMortarContact;
  boost::shared_ptr<SimpleContactProblem::SimpleContactElement>
      feRhsSimpleContactALEMaterial;
  boost::shared_ptr<SimpleContactProblem::SimpleContactElement>
      feLhsSimpleContactALEMaterial;
  boost::shared_ptr<SimpleContactProblem::SimpleContactElement>
      feLhsSimpleContactALE;
  boost::shared_ptr<SimpleContactProblem::CommonDataSimpleContact>
      commonDataSimpleContactALE;
 
  PetscBool contactOutputIntegPts;
  boost::shared_ptr<SimpleContactProblem::SimpleContactElement>
      fePostProcSimpleContact;
  boost::shared_ptr<MortarContactProblem::MortarContactElement>
      fePostProcMortarContact;
  moab::Core contactPostProcCore;
  moab::Interface &contactPostProcMoab;
 
  Range constrainedInterface; ///< Range of faces on the constrained interface
 
  /**@}*/
 
  /** \name  Data for bone */
 
  /**@{*/
 
  double rHo0;  ///< Reference density if bone is analyzed
  double nBone; ///< Exponent parameter in bone density
 
  /**@}*/
 
  /** \name  Interfaces */
 
  /**@{*/
 
  const boost::scoped_ptr<UnknownInterface> cpSolversPtr;
  const boost::scoped_ptr<UnknownInterface> cpMeshCutPtr;
 
  /**@}*/
 
  /** \name Cut/Refine/Split options */
 
  /**@{*/
 
  PetscBool doCutMesh;
  double crackAccelerationFactor;
  PetscBool doElasticWithoutCrack;
 
  /**@}*/
};
 
} // namespace FractureMechanics
 
#endif // __CRACK_PROPAGATION_HPP__