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 {
 
  enum MaterialModel {
    StVenantKirchhoff,
    MooneyRivlin,
    Hencky,
    Neohookean,
    LastMaterial
  };
 
  enum SolverType { TimeSolver, DynamicRelaxation, Cohesive, LastSolver };
 
  static inline enum SolverType solverType = TimeSolver;
  static inline enum MaterialModel materialModel = MooneyRivlin;
  static inline enum SymmetrySelector symmetrySelector = NOT_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_FALSE; //< set singularity
  static inline PetscBool dynamicRelaxation =
      PETSC_FALSE; //< switch on/off dynamic relaxation
  static inline PetscBool crackingOn = PETSC_FALSE; //< cracking on
  static inline double crackingStartTime = 0; //< time when crack starts to grow
  static inline int nbJIntegralContours =
      0; //< number of contours for J integral evaluation
  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_TRUE; //< scale L2 user base
  static inline int addCrackMeshsetId = 1000;           //< add crack meshset id
  static inline double griffithEnergy = 1;              ///< Griffith energy
  static inline enum EnergyReleaseSelector energyReleaseSelector =
      GRIFFITH_SKELETON; //< energy release selector
  static inline std::string internalStressTagName =
      ""; //< internal stress tag name
  static inline int internalStressInterpOrder =
      1; //< internal stress interpolation order
  static inline PetscBool internalStressVoigt =
      PETSC_FALSE; //< internal stress index notation
  static inline PetscBool intefaceCrack =
      PETSC_FALSE; //< interface crack tracking
 
  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 inline constexpr bool use_quadratic_exp = true;
  static inline constexpr double v_max = 24;
  static inline constexpr double v_min = -v_max;
 
  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 f_log_e(const double v) {
    if constexpr(use_quadratic_exp) {
      return f_log_e_quadratic(v);
    } else {
      if (v > v_max)
        // y = exp(v_max) * v + exp(v_max) * (1 - v_max);
        // y/exp(v_max) = v + (1 - v_max);
        // y/exp(v_max) - (1 - v_max) = v;
        return std::exp(v_max) * v + std::exp(v_max) * (1 - v_max);
      else 
        return std::exp(v);
    }
  }
  static double d_f_log_e(const double v) {
    if constexpr (use_quadratic_exp) {
      return d_f_log_e_quadratic(v);
    } else {
      if (v > v_max)
        return std::exp(v_max);
      else 
        return std::exp(v);
    }
  }
  static double dd_f_log_e(const double v) {
    if constexpr (use_quadratic_exp) {
      return dd_f_log_e_quadratic(v);
    } else {
      if (v > v_max)
        return 0.;
      else 
        return std::exp(v);
    }
  }
  static double inv_f_log_e(const double v) {
    if (v > std::exp(v_min))
      return std::log(v);
    else
      // y = exp(v_min) * v + exp(v_min) * (1 - v_min);
      // y/exp(v_min) = v + (1 - v_min);
      // y/exp(v_min) - (1 - v_min) = v;
      return v / exp(v_min) - (1 - v_min);
  }
  static double inv_d_f_log_e(const double v) {
    if (v > std::exp(v_min))
      return 1. / v;
    else
      return 1. / exp(v_min);
  }
  static double inv_dd_f_log_e(const double v) {
    if (v > std::exp(v_min))
      return -1. / (v * v);
    else
      return 0.;
  }
 
  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<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
 
  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";
 
  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;
  double alphaTau = 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> bcSpatialTractionVecPtr;
  boost::shared_ptr<TractionFreeBc> bcSpatialFreeTractionVecPtr;
  boost::shared_ptr<NormalDisplacementBcVec> bcSpatialNormalDisplacementVecPtr;
  boost::shared_ptr<AnalyticalDisplacementBcVec>
      bcSpatialAnalyticalDisplacementVecPtr;
  boost::shared_ptr<AnalyticalTractionBcVec> bcSpatialAnalyticalTractionVecPtr;
  boost::shared_ptr<PressureBcVec> bcSpatialPressureVecPtr;
  boost::shared_ptr<ExternalStrainVec> externalStrainVecPtr;
 
  std::map<std::string, boost::shared_ptr<ScalingMethod>> timeScaleMap;
 
  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) {
        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);
    }
    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);
 
    auto ts_rotation =
        boost::make_shared<DynamicRelaxationTimeScale>("rotation_history.txt");
    for (auto &bc : *bcSpatialRotationVecPtr) {
      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 addFields(const EntityHandle meshset = 0);
  MoFEMErrorCode projectGeometry(const EntityHandle meshset = 0);
  MoFEMErrorCode projectInternalStress(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 setFaceInterfaceOps(
      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);
 
  friend struct solve_elastic_set_up;
  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 cohesive crack growth problem
   * 
   * @param ts 
   * @param x 
   * @return * MoFEMErrorCode 
   */
  MoFEMErrorCode solveCohesiveCrackGrowth(TS ts, Vec x, int start_step,
                                          double start_time);
 
  MoFEMErrorCode setBlockTagsOnSkin();
 
  MoFEMErrorCode postProcessResults(const int tag, const std::string file,
                                    Vec f_residual = PETSC_NULLPTR,
                                    Vec var_vec = PETSC_NULLPTR,
                                    std::vector<Tag> tags_to_transfer = {});
  MoFEMErrorCode
  postProcessSkeletonResults(const int tag, const std::string file,
                             Vec f_residual = PETSC_NULLPTR,
                             std::vector<Tag> tags_to_transfer = {});
 
  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 {
      switch (EshelbianCore::solverType) {
      case EshelbianCore::Cohesive:
      case EshelbianCore::DynamicRelaxation:
        return TimeScale::getScale(EshelbianCore::dynamicTime);
      default:
        return TimeScale::getScale(time);
      }
    }
  };
 
  MoFEMErrorCode calculateFaceMaterialForce(const int tag, TS ts);
  MoFEMErrorCode calculateOrientation(const int tag, bool set_orientation);
  MoFEMErrorCode setNewFrontCoordinates();
  MoFEMErrorCode addCrackSurfaces(const bool debug = false);
  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> maxMovedFaces;
  boost::shared_ptr<Range> interfaceFaces;
 
  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_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
};