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 enum SymmetrySelector symmetrySelector = SYMMETRIC;
  static inline enum RotSelector rotSelector = LARGE_ROT;
  static inline enum RotSelector gradApproximator = LARGE_ROT;
  static inline enum StretchSelector stretchSelector = LOG;
  static inline PetscBool noStretch = PETSC_FALSE;     //< no stretch field
  static inline PetscBool setSingularity = PETSC_TRUE; //< set singularity
  static inline PetscBool dynamicRelaxation =
      PETSC_FALSE; //< switch on/off dynamic relaxation
  static inline PetscBool crackingOn = PETSC_FALSE; //< cracking on
  static inline double dynamicTime =
      0; //< true time used when dynamic relaxation is used.
  static inline int dynamicStep =
      0; //< true time used when dynamic relaxation is used.
  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 double exponentBase;
  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 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_log(const double v) {
    return pow(EshelbianCore::exponentBase, v);
  }
  static double d_f_log(const double v) {
    return pow(exponentBase, v) * log(EshelbianCore::exponentBase);
  }
  static double dd_f_log(const double v) {
    return pow(EshelbianCore::exponentBase, v) *
           log(EshelbianCore::exponentBase) * log(EshelbianCore::exponentBase);
  }
 
  static double inv_f_log(const double v) {
    return log(v) / log(EshelbianCore::exponentBase);
  }
  static double inv_d_f_log(const double v) {
    return (1. / v) / log(EshelbianCore::exponentBase);
  }
  static double inv_dd_f_log(const double v) {
    return -(1. / (v * v)) / log(EshelbianCore::exponentBase);
  }
 
  static double f_linear(const double v) { return v + 1; }
  static double d_f_linear(const double v) { return 1; }
  static double dd_f_linear(const double v) { return 0; }
 
  static double inv_f_linear(const double v) { return v - 1; }
  static double inv_d_f_linear(const double v) { return 0; }
  static double inv_dd_f_linear(const double v) { 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<VolumeElementForcesAndSourcesCore> elasticFeRhs;
  boost::shared_ptr<VolumeElementForcesAndSourcesCore> elasticFeLhs;
  boost::shared_ptr<FaceElementForcesAndSourcesCore> elasticBcLhs;
  boost::shared_ptr<FaceElementForcesAndSourcesCore> elasticBcRhs;
  boost::shared_ptr<ContactTree> 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
 
  const std::string piolaStress = "P";
  const std::string eshelbyStress = "S";
  const std::string spatialL2Disp = "wL2";
  const std::string materialL2Disp = "WL2";
  const std::string spatialH1Disp = "wH1";
  const std::string materialH1Positions = "XH1";
  const std::string hybridSpatialDisp = "hybridSpatialDisp";
  const std::string hybridMaterialDisp = "hybridMaterialDisp";
  const std::string contactDisp = "contactDisp";
  const std::string stretchTensor = "u";
  const std::string rotAxis = "omega";
  const std::string bubbleField = "bubble";
 
  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";
  const std::string materialVolumeElement = "MEP";
  const std::string materialSkeletonElement = "MATERIAL_SKELETON";
 
  EshelbianCore(MoFEM::Interface &m_field);
  virtual ~EshelbianCore();
 
  int spaceOrder = 2;
  int spaceH1Order = -1;
  int materialOrder = 1;
  double alphaU = 0;
  double alphaW = 0;
  double alphaOmega = 0;
  double alphaRho = 0;
 
  int contactRefinementLevels = 1; //< refinement levels for contact integration
  int frontLayers = 3;             //< front layers with material element
 
  MoFEMErrorCode getOptions();
 
  boost::shared_ptr<BcDispVec> bcSpatialDispVecPtr;
  boost::shared_ptr<BcRotVec> bcSpatialRotationVecPtr;
  boost::shared_ptr<TractionBcVec> bcSpatialTraction;
  boost::shared_ptr<TractionFreeBc> bcSpatialFreeTraction;
 
  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() < nb_attributes) {
        SETERRQ3(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                 "In block %s expected %d attributes, but given %d",
                 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);
    }
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode getSpatialDispBc();
 
  inline MoFEMErrorCode getSpatialRotationBc() {
    bcSpatialRotationVecPtr = boost::make_shared<BcRotVec>();
    return getBc(bcSpatialRotationVecPtr, "SPATIAL_ROTATION_BC", 4);
  }
 
  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) {
    bcSpatialFreeTraction =
        boost::shared_ptr<TractionFreeBc>(new TractionFreeBc());
    return getTractionFreeBc(meshset, bcSpatialFreeTraction, "CONTACT");
  }
 
  MoFEMErrorCode createExchangeVectors(Sev sev);
 
  MoFEMErrorCode addFields(const EntityHandle meshset = 0);
  MoFEMErrorCode projectGeometry(const EntityHandle meshset = 0);
 
  MoFEMErrorCode addVolumeFiniteElement(const EntityHandle meshset = 0);
  MoFEMErrorCode addBoundaryFiniteElement(const EntityHandle meshset = 0);
  MoFEMErrorCode addDMs(const BitRefLevel bit = BitRefLevel().set(0),
                        const EntityHandle meshset = 0);
 
  MoFEMErrorCode addMaterial_HMHHStVenantKirchhoff(const int tape,
                                                   const double lambda,
                                                   const double mu,
                                                   const double sigma_y);
 
  MoFEMErrorCode addMaterial_HMHMooneyRivlin(const int tape, const double alpha,
                                             const double beta,
                                             const double lambda,
                                             const double sigma_y);
 
  MoFEMErrorCode addMaterial_HMHNeohookean(const int tape, const double c10,
                                           const double K);
 
  MoFEMErrorCode addMaterial_Hencky(double E, double nu);
 
  MoFEMErrorCode setBaseVolumeElementOps(
      const int tag, const bool do_rhs, const bool do_lhs,
      const bool calc_rates, SmartPetscObj<Vec> ver_vec,
      boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe);
 
  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
  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<ContactTree> &fe_contact_tree
 
  );
 
  MoFEMErrorCode setElasticElementOps(const int tag);
  MoFEMErrorCode setElasticElementToTs(DM dm);
 
  friend struct solve_elastic_set_up;
  MoFEMErrorCode solveElastic(TS ts, Vec x);
  MoFEMErrorCode solveDynamicRelaxation(TS ts, Vec x);
 
  MoFEMErrorCode setBlockTagsOnSkin();
 
  MoFEMErrorCode postProcessResults(const int tag, const std::string file,
                                    Vec f_residual = PETSC_NULL,
                                    Vec var_vec = PETSC_NULL,
                                    std::vector<Tag> tags_to_transfer = {});
  MoFEMErrorCode postProcessSkeletonResults(const int tag,
                                            const std::string file,
                                            Vec f_residual = PETSC_NULL);
 
  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;
  };
 
  MoFEMErrorCode calculateEnergyRelease(const int tag, TS ts);
  MoFEMErrorCode calculateOrientation(const int tag, bool set_orientation);
  MoFEMErrorCode setNewFrontCoordinates();
  MoFEMErrorCode addCrackSurfaces();
  MoFEMErrorCode saveOrgCoords();
  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> materialSkeletonFaces;
  boost::shared_ptr<Range> maxMovedFaces;
 
  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;
 
  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_NULL; //< Schur complement matrix
  AO aoS = PETSC_NULL; //< 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
 
 
};