EshelbianPlasticity.hpp

Go to the documentation of this file.

/**
 * \file EshelbianPlasticity.hpp
 * \brief Eshelbian plasticity interface
 *
 * \brief Problem implementation for mix element for large-strain elasticity
 *
 * For reference on mixed formulation see: \cite gopalakrishnan2012second and
 * \cite cockburn2010new
 *
 * \todo Implementation of plasticity
 */
 
// DO NOT DELETE MIGHT BE USEFULL
// #include <boost/multiprecision/mpfr.hpp>
// using mpfr50 = boost::multiprecision::number<
//     boost::multiprecision::mpfr_float_backend<50>>;
 
// inline mpfr50 exp_hi(const mpfr50 x) { return boost::multiprecision::exp(x); }
// inline mpfr50 log_hi(const mpfr50 x) { return boost::multiprecision::log(x); }
 
constexpr int SPACE_DIM = 3;
// constexpr auto A = AssemblyType::BLOCK_SCHUR;
constexpr auto A = AssemblyType::BLOCK_MAT;
 
#ifndef __ESHELBIAN_PLASTICITY_HPP__
  #define __ESHELBIAN_PLASTICITY_HPP__
 
  #ifdef ENABLE_PYTHON_BINDING
    #include <boost/python.hpp>
    #include <boost/python/def.hpp>
    #include <boost/python/numpy.hpp>
namespace bp = boost::python;
namespace np = boost::python::numpy;
  #endif
namespace EshelbianPlasticity {
 
struct ContactTree;
 
enum SymmetrySelector { SYMMETRIC, NOT_SYMMETRIC, ANIT_SYMMETRIC };
enum RotSelector { SMALL_ROT, MODERATE_ROT, LARGE_ROT, NO_H1_CONFIGURATION };
enum StretchSelector {
  LINEAR,
  LOG /*linar extenion*/,
  LOG_QUADRATIC /*quadratic extension*/,
  STRETCH_SELECTOR_LAST
};
enum EnergyReleaseSelector { GRIFFITH_FORCE, GRIFFITH_SKELETON };
 
using MatrixPtr = boost::shared_ptr<MatrixDouble>;
using VectorPtr = boost::shared_ptr<VectorDouble>;
 
using EntData = EntitiesFieldData::EntData;
using UserDataOperator = ForcesAndSourcesCore::UserDataOperator;
using VolUserDataOperator = VolumeElementForcesAndSourcesCore::UserDataOperator;
using FaceUserDataOperator = FaceElementForcesAndSourcesCore::UserDataOperator;
using EleOnSide = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
using SideEleOp = EleOnSide::UserDataOperator;
 
struct PhysicalEquations;
 
struct AnalyticalExprPython;
 
struct DataAtIntegrationPts
    : public boost::enable_shared_from_this<DataAtIntegrationPts> {
 
  MatrixDouble approxPAtPts;
  MatrixDouble approxSigmaAtPts;
  MatrixDouble divPAtPts;
  MatrixDouble gradPAtPts;
  MatrixDouble divSigmaAtPts;
  MatrixDouble wL2AtPts;
  MatrixDouble wL2DotAtPts;
  MatrixDouble wL2DotDotAtPts;
  MatrixDouble logStretchTensorAtPts;
  MatrixDouble stretchTensorAtPts;
  VectorDouble jacobianAtPts;
  MatrixDouble tractionAtPts;
 
  MatrixDouble diffStretchTensorAtPts;
  VectorDouble detStretchTensorAtPts;
  MatrixDouble detStretchTensorAtPts_du;
  MatrixDouble logStretchDotTensorAtPts;
  MatrixDouble gradLogStretchDotTensorAtPts;
  MatrixDouble rotAxisAtPts;
  MatrixDouble rotAxisDotAtPts;
  MatrixDouble rotAxisGradDotAtPts;
  MatrixDouble rotAxis0AtPts;
  MatrixDouble WAtPts;
  MatrixDouble W0AtPts;
  MatrixDouble GAtPts;
  MatrixDouble G0AtPts;
  MatrixDouble wH1AtPts;
  MatrixDouble XH1AtPts;
  MatrixDouble contactL2AtPts;
  MatrixDouble wGradH1AtPts;
  MatrixDouble logStretch2H1AtPts;
  MatrixDouble logStretchTotalTensorAtPts;
 
  MatrixDouble hAtPts;
  MatrixDouble hdOmegaAtPts;
  MatrixDouble hdLogStretchAtPts;
  MatrixDouble leviKirchhoffAtPts;
  MatrixDouble leviKirchhoffdOmegaAtPts;
  MatrixDouble leviKirchhoffdLogStreatchAtPts;
  MatrixDouble leviKirchhoffPAtPts;
  MatrixDouble adjointPdstretchAtPts;
  MatrixDouble adjointPdUAtPts;
  MatrixDouble adjointPdUdOmegaAtPts;
  MatrixDouble adjointPdUdPAtPts;
  MatrixDouble rotMatAtPts;
  MatrixDouble PAtPts;
  MatrixDouble SigmaAtPts;
  VectorDouble energyAtPts; //< this is density of energy at integration points
  MatrixDouble flowL2AtPts;
 
  MatrixDouble varRotAxis;
  MatrixDouble varGradRotAxis;
  MatrixDouble varLogStreach;
  MatrixDouble varPiola;
  MatrixDouble varDivPiola;
  MatrixDouble varWL2;
  MatrixDouble varHybridDispAtPts;
 
  MatrixDouble P_du;
 
  MatrixDouble eigenVals;
  MatrixDouble eigenVecs;
  VectorInt nbUniq;
  MatrixDouble eigenValsC;
  MatrixDouble eigenVecsC;
  VectorInt nbUniqC;
 
  MatrixDouble matD;
  MatrixDouble matAxiatorD;
  MatrixDouble matDeviatorD;
  MatrixDouble matInvD;
 
  MatrixDouble internalStressAtPts;
 
  MatrixDouble facePiolaAtPts;
  MatrixDouble hybridDispAtPts;
  MatrixDouble hybridDispDotAtPts;
  MatrixDouble gradHybridDispAtPts;
  MatrixDouble faceMaterialForceAtPts;
  VectorDouble normalPressureAtPts;
 
  double mu;
  double lambda;
  double piolaScale = 1.;
  double faceEnergy = 0.;
 
  inline auto getPiolaScalePtr() {
    return boost::shared_ptr<double>(shared_from_this(), &piolaScale);
  }
 
  inline MatrixPtr getApproxSigmaAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(),
                                           &approxSigmaAtPts);
  }
  inline MatrixPtr getApproxPAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &approxPAtPts);
  }
 
  inline MatrixPtr getDivPAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &divPAtPts);
  }
 
  inline MatrixPtr getGradPAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &gradPAtPts);
  }
 
  inline MatrixPtr getDivSigmaAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &divSigmaAtPts);
  }
 
  inline MatrixPtr getSmallWL2AtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &wL2AtPts);
  }
 
  inline MatrixPtr getSmallWL2DotAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &wL2DotAtPts);
  }
 
  inline MatrixPtr getSmallWL2DotDotAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &wL2DotDotAtPts);
  }
 
  inline MatrixPtr getLogStretchTensorAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(),
                                           &logStretchTensorAtPts);
  }
 
  inline MatrixPtr getStretchTensorAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(),
                                           &stretchTensorAtPts);
  }
 
  inline MatrixPtr getLogStretchDotTensorAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(),
                                           &logStretchDotTensorAtPts);
  }
 
  inline MatrixPtr getGradLogStretchDotTensorAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(),
                                           &gradLogStretchDotTensorAtPts);
  }
 
  inline MatrixPtr getRotAxisAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &rotAxisAtPts);
  }
 
  inline MatrixPtr getRotAxis0AtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &rotAxis0AtPts);
  }
 
  inline MatrixPtr getRotAxisDotAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(),
                                           &rotAxisDotAtPts);
  }
 
  inline MatrixPtr getRotAxisGradDotAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(),
                                           &rotAxisGradDotAtPts);
  }
 
  inline MatrixPtr getBigGAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &GAtPts);
  }
 
  inline MatrixPtr getBigG0AtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &GAtPts);
  }
 
  inline MatrixPtr getMatDPtr() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matD);
  }
 
  inline MatrixPtr getMatInvDPtr() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matInvD);
  }
 
  inline MatrixPtr getMatAxiatorDPtr() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matAxiatorD);
  }
 
  inline MatrixPtr getMatDeviatorDPtr() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matDeviatorD);
  }
 
  inline MatrixPtr getSmallWH1AtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &wH1AtPts);
  }
 
  inline MatrixPtr getLargeXH1AtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &XH1AtPts);
  }
 
  inline MatrixPtr getContactL2AtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &contactL2AtPts);
  }
 
  inline MatrixPtr getSmallWGradH1AtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &wGradH1AtPts);
  }
 
  inline VectorPtr getJacobianAtPts() {
    return boost::shared_ptr<VectorDouble>(shared_from_this(), &jacobianAtPts);
  }
 
  inline MatrixPtr getLeviKirchhoffAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(),
                                           &leviKirchhoffAtPts);
  };
 
  inline MatrixPtr getVarRotAxisPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &varRotAxis);
  }
 
  inline MatrixPtr getVarGradRotAxisPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(),
                                           &varGradRotAxis);
  }
 
  inline MatrixPtr getVarLogStreachPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &varLogStreach);
  }
 
  inline MatrixPtr getVarPiolaPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &varPiola);
  }
 
  inline MatrixPtr getDivVarPiolaPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &varDivPiola);
  }
 
  inline MatrixPtr getVarWL2Pts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &varWL2);
  }
 
  inline MatrixPtr getVarHybridDispAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(),
                                           &varHybridDispAtPts);
  }
 
  inline MatrixPtr getInternalStressAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(),
                                           &internalStressAtPts);
  }
 
  inline MatrixPtr getFacePiolaAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &facePiolaAtPts);
  }
 
  inline MatrixPtr getHybridDispAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(),
                                           &hybridDispAtPts);
  }
 
  inline MatrixPtr getHybridDispDotAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(),
                                           &hybridDispDotAtPts);
  }
 
  inline MatrixPtr getGradHybridDispAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(),
                                           &gradHybridDispAtPts);
  }
 
  inline MatrixPtr getEigenValsAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &eigenVals);
  } 
 
  boost::shared_ptr<PhysicalEquations> physicsPtr;
};
 
struct OpJacobian;
 
// Forward declarations
struct ExternalStrain;
using ExternalStrainVec = std::vector<ExternalStrain>;
 
struct PhysicalEquations {
 
  typedef FTensor::Tensor1<adouble, 3> ATensor1;
  typedef FTensor::Tensor2<adouble, 3, 3> ATensor2;
  typedef FTensor::Tensor3<adouble, 3, 3, 3> ATensor3;
  typedef FTensor::Tensor1<double, 3> DTensor1;
  typedef FTensor::Tensor2<double, 3, 3> DTensor2;
  typedef FTensor::Tensor3<double, 3, 3, 3> DTensor3;
 
  typedef FTensor::Tensor0<FTensor::PackPtr<double *, 1>> DTensor0Ptr;
 
  typedef FTensor::Tensor2<FTensor::PackPtr<double *, 1>, 3, 3> DTensor2Ptr;
  typedef FTensor::Tensor3<FTensor::PackPtr<double *, 1>, 3, 3, 3> DTensor3Ptr;
 
  PhysicalEquations() = delete;
  PhysicalEquations(const int size_active, const int size_dependent)
      : activeVariables(size_active, 0),
        dependentVariablesPiola(size_dependent, 0),
        dependentVariablesPiolaDirevatives(size_dependent * size_active, 0) {}
  virtual ~PhysicalEquations() = default;
 
  virtual MoFEMErrorCode setParams(const int tag, EntityHandle ent, int gg) {
    return 0;
  }
  virtual MoFEMErrorCode recordTape(const int tag, DTensor2Ptr *t_h) = 0;
 
  virtual UserDataOperator *
  returnOpJacobian(const int tag, const bool eval_rhs, const bool eval_lhs,
                   boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                   boost::shared_ptr<PhysicalEquations> physics_ptr);
 
  virtual VolUserDataOperator *
  returnOpSpatialPhysical(const std::string &field_name,
                          boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                          const double alpha_u);
 
  virtual VolUserDataOperator *
    returnOpSpatialPhysicalExternalStrain(
      const std::string &field_name,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      boost::shared_ptr<ExternalStrainVec> external_strain_vec_ptr,
      std::map<std::string, boost::shared_ptr<ScalingMethod>> smv);
 
  virtual VolUserDataOperator *returnOpSpatialPhysical_du_du(
      std::string row_field, std::string col_field,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr, const double alpha);
 
  virtual VolUserDataOperator *
  returnOpCalculateEnergy(boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                          boost::shared_ptr<double> total_energy_ptr);
 
  virtual VolUserDataOperator *returnOpCalculateStretchFromStress(
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      boost::shared_ptr<PhysicalEquations> physics_ptr);
 
  virtual VolUserDataOperator *returnOpCalculateVarStretchFromStress(
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      boost::shared_ptr<PhysicalEquations> physics_ptr);
 
  virtual VolUserDataOperator *
  returnOpSetScale(boost::shared_ptr<double> scale_ptr,
                   boost::shared_ptr<PhysicalEquations> physics_ptr);
 
  std::vector<double> activeVariables;
  std::vector<double> dependentVariablesPiola;
  std::vector<double> dependentVariablesPiolaDirevatives;
 
  /** \name Active variables */
 
  /**@{*/
 
  template <int S>
  inline static DTensor2Ptr get_VecTensor2(std::vector<double> &v) {
    return DTensor2Ptr(&v[S + 0], &v[S + 1], &v[S + 2], &v[S + 3], &v[S + 4],
                       &v[S + 5], &v[S + 6], &v[S + 7], &v[S + 8]);
  }
 
  template <int S>
  inline static DTensor0Ptr get_VecTensor0(std::vector<double> &v) {
    return DTensor0Ptr(&v[S + 0]);
  }
 
  template <int S0>
  inline static DTensor3Ptr get_vecTensor3(std::vector<double> &v,
                                           const int nba) {
 
    const int A00 = nba * 0 + S0;
    const int A01 = nba * 1 + S0;
    const int A02 = nba * 2 + S0;
    const int A10 = nba * 3 + S0;
    const int A11 = nba * 4 + S0;
    const int A12 = nba * 5 + S0;
    const int A20 = nba * 6 + S0;
    const int A21 = nba * 7 + S0;
    const int A22 = nba * 8 + S0;
 
    return DTensor3Ptr(
 
        &v[A00 + 0], &v[A00 + 1], &v[A00 + 2], &v[A01 + 0], &v[A01 + 1],
        &v[A01 + 2], &v[A02 + 0], &v[A02 + 1], &v[A02 + 2],
 
        &v[A10 + 0], &v[A10 + 1], &v[A10 + 2], &v[A11 + 0], &v[A11 + 1],
        &v[A11 + 2], &v[A12 + 0], &v[A12 + 1], &v[A12 + 2],
 
        &v[A20 + 0], &v[A20 + 1], &v[A20 + 2], &v[A21 + 0], &v[A21 + 1],
        &v[A21 + 2], &v[A22 + 0], &v[A22 + 1], &v[A22 + 2]
 
    );
  }
 
  /**@}*/
 
  /** \name Active variables */
 
  /**@{*/
 
  inline DTensor2Ptr get_h() { return get_VecTensor2<0>(activeVariables); }
 
  /**@}*/
 
  /** \name Dependent variables */
 
  /**@{*/
 
  inline DTensor2Ptr get_P() {
    return get_VecTensor2<0>(dependentVariablesPiola);
  }
 
  /**@}*/
 
  /** \name Derivatives of dependent variables */
 
  /**@{*/
 
  inline DTensor3Ptr get_P_dh0() {
    return get_vecTensor3<0>(dependentVariablesPiolaDirevatives,
                             activeVariables.size());
  }
  inline DTensor3Ptr get_P_dh1() {
    return get_vecTensor3<3>(dependentVariablesPiolaDirevatives,
                             activeVariables.size());
  }
  inline DTensor3Ptr get_P_dh2() {
    return get_vecTensor3<6>(dependentVariablesPiolaDirevatives,
                             activeVariables.size());
  }
 
  /**@}*/
};
 
struct BcDisp {
  BcDisp(std::string name, std::vector<double> attr, Range faces);
  std::string blockName;
  Range faces;
  VectorDouble3 vals;
  VectorInt3 flags;
};
using BcDispVec = std::vector<BcDisp>;
 
struct BcRot {
  BcRot(std::string name, std::vector<double> attr, Range faces);
  std::string blockName;
  Range faces;
  VectorDouble vals;
  double theta;
};
using BcRotVec = std::vector<BcRot>;
 
typedef std::vector<Range> TractionFreeBc;
 
struct TractionBc {
  TractionBc(std::string name, std::vector<double> attr, Range faces);
  std::string blockName;
  Range faces;
  VectorDouble3 vals;
  VectorInt3 flags;
};
using TractionBcVec = std::vector<TractionBc>;
 
struct NormalDisplacementBc {
  NormalDisplacementBc(std::string name, std::vector<double> attr, Range faces);
  std::string blockName;
  Range faces;
  double val;
};
using NormalDisplacementBcVec = std::vector<NormalDisplacementBc>;
 
struct AnalyticalDisplacementBc {
  AnalyticalDisplacementBc(std::string name, std::vector<double> attr,
                           Range faces);
  std::string blockName;
  Range faces;
  VectorInt3 flags;
};
using AnalyticalDisplacementBcVec = std::vector<AnalyticalDisplacementBc>;
 
struct AnalyticalTractionBc {
  AnalyticalTractionBc(std::string name, std::vector<double> attr, Range faces);
  std::string blockName;
  Range faces;
  VectorInt3 flags;
};
using AnalyticalTractionBcVec = std::vector<AnalyticalTractionBc>;
 
struct PressureBc {
  PressureBc(std::string name, std::vector<double> attr, Range faces);
  std::string blockName;
  Range faces;
  double val;
};
using PressureBcVec = std::vector<PressureBc>;
 
struct ExternalStrain {
  ExternalStrain(std::string name, std::vector<double> attr, Range ents);
  std::string blockName;
  Range ents;
  double val;
  double bulkModulusK;
};
 
template <typename OP_PTR>
std::tuple<std::string, MatrixDouble>
getAnalyticalExpr(OP_PTR op_ptr, MatrixDouble &analytical_expr,
                  const std::string block_name);
 
template <typename OP_PTR>
std::tuple<std::string, VectorDouble>
getAnalyticalExternalStrain(OP_PTR op_ptr, VectorDouble &analytical_expr,
                            const std::string block_name);
 
MatrixDouble analytical_expr_function(double delta_t, double t,
                                      int nb_gauss_pts,
                                      MatrixDouble &m_ref_coords,
                                      MatrixDouble &m_ref_normals,
                                      const std::string block_name);
 
                                      
VectorDouble analytical_externalstrain_function(double delta_t, double t,
                                             int nb_gauss_pts,
                                             MatrixDouble &m_ref_coords,
                                             const std::string block_name);
 
  #ifdef ENABLE_PYTHON_BINDING
struct AnalyticalExprPython {
  AnalyticalExprPython() = default;
  virtual ~AnalyticalExprPython() = default;
 
  MoFEMErrorCode analyticalExprInit(const std::string py_file);
  MoFEMErrorCode evalAnalyticalDisp(double delta_t, double t, np::ndarray x,
                                    np::ndarray y, np::ndarray z,
                                    np::ndarray nx, np::ndarray ny,
                                    np::ndarray nz,
                                    const std::string &block_name,
                                    np::ndarray &analytical_expr);
  MoFEMErrorCode evalAnalyticalTraction(double delta_t, double t, np::ndarray x,
                                        np::ndarray y, np::ndarray z,
                                        np::ndarray nx, np::ndarray ny,
                                        np::ndarray nz,
                                        const std::string &block_name,
                                        np::ndarray &analytical_expr);
 
  MoFEMErrorCode evalAnalyticalExternalStrain(double delta_t, double t,
                                              np::ndarray x, np::ndarray y,
                                              np::ndarray z,
                                              const std::string &block_name,
                                              np::ndarray &analytical_expr);
  template <typename T>
  inline std::vector<T>
  py_list_to_std_vector(const boost::python::object &iterable) {
    return std::vector<T>(boost::python::stl_input_iterator<T>(iterable),
                          boost::python::stl_input_iterator<T>());
  }
 
private:
  bp::object mainNamespace;
  bp::object analyticalDispFun;
  bp::object analyticalTractionFun;
  bp::object analyticalExternalStrainFun;
};
 
extern boost::weak_ptr<AnalyticalExprPython> AnalyticalExprPythonWeakPtr;
  #endif
 
  #include "CGGUserPolynomialBase.hpp"
  #include "EshelbianCore.hpp"
  #include "EshelbianOperators.hpp"
 
} // namespace EshelbianPlasticity
 
#endif //__ESHELBIAN_PLASTICITY_HPP__