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
 */
 
constexpr int SPACE_DIM = 3;
constexpr auto A = AssemblyType::BLOCK_SCHUR;
 
#ifndef __ESHELBIAN_PLASTICITY_HPP__
  #define __ESHELBIAN_PLASTICITY_HPP__
 
namespace EshelbianPlasticity {
 
struct CGGUserPolynomialBase : public TetPolynomialBase {
 
  using CachePhi = boost::tuple<int, int, MatrixDouble>;
 
  CGGUserPolynomialBase(boost::shared_ptr<CachePhi> cache_phi = nullptr);
  ~CGGUserPolynomialBase() = default;
 
  MoFEMErrorCode query_interface(boost::typeindex::type_index type_index,
                                 BaseFunctionUnknownInterface **iface) const;
  MoFEMErrorCode getValue(MatrixDouble &pts,
                          boost::shared_ptr<BaseFunctionCtx> ctx_ptr);
 
private:
 
  MatrixDouble shapeFun;
  boost::shared_ptr<CachePhi> cachePhiPtr;
 
  MoFEMErrorCode getValueHdivForCGGBubble(MatrixDouble &pts);
};
 
struct ContactTree;
 
enum RotSelector { SMALL_ROT, MODERATE_ROT, LARGE_ROT, NO_H1_CONFIGURATION };
enum StretchSelector { LINEAR, LOG };
 
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;
 
struct PhysicalEquations;
struct DataAtIntegrationPts
    : public boost::enable_shared_from_this<DataAtIntegrationPts> {
 
  MatrixDouble approxPAtPts;
  MatrixDouble approxSigmaAtPts;
  MatrixDouble divPAtPts;
  MatrixDouble divSigmaAtPts;
  MatrixDouble wL2AtPts;
  MatrixDouble wL2DotAtPts;
  MatrixDouble wL2DotDotAtPts;
  MatrixDouble logStretchTensorAtPts;
  MatrixDouble stretchTensorAtPts;
  VectorDouble jacobianAtPts;
 
  MatrixDouble diffStretchTensorAtPts;
  VectorDouble detStretchTensorAtPts;
  MatrixDouble detStretchTensorAtPts_du;
  MatrixDouble logStretchDotTensorAtPts;
  MatrixDouble rotAxisAtPts;
  MatrixDouble rotAxisDotAtPts;
  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 leviKirchhoffPAtPts;
  MatrixDouble leviKirchhoffOmegaAtPts;
  MatrixDouble adjointPdstretchAtPts;
  MatrixDouble adjointPdUAtPts;
  MatrixDouble adjointPdUdOmegaAtPts;
  MatrixDouble adjointPdUdPAtPts;
  MatrixDouble rotMatAtPts;
  MatrixDouble diffRotMatAtPts;
  MatrixDouble diff2RotMatAtPts;
  MatrixDouble PAtPts;
  MatrixDouble SigmaAtPts;
  VectorDouble energyAtPts; //< this is density of energy at integration points
  MatrixDouble flowL2AtPts;
 
  MatrixDouble P_du;
 
  MatrixDouble eigenVals;
  MatrixDouble eigenVecs;
  VectorInt nbUniq;
 
  MatrixDouble matD;
  MatrixDouble matAxiatorD;
  MatrixDouble matDeviatorD;
  MatrixDouble matInvD;
 
  double mu;
  double lambda;
 
  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 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 getRotAxisAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &rotAxisAtPts);
  }
 
  inline MatrixPtr getRotAxisDotAtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(),
                                           &rotAxisDotAtPts);
  }
 
  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);
  }
 
  boost::shared_ptr<PhysicalEquations> physicsPtr;
};
 
struct OpJacobian;
 
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 recordTape(const int tag, DTensor2Ptr *t_h) = 0;
 
  virtual OpJacobian *
  returnOpJacobian(const int tag, const bool eval_rhs, const bool eval_lhs,
                   boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                   boost::shared_ptr<PhysicalEquations> physics_ptr) = 0;
 
  virtual VolUserDataOperator *
  returnOpSpatialPhysical(const std::string &field_name,
                          boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                          const double alpha_u);
 
  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);
 
  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;
};
typedef std::vector<BcDisp> BcDispVec;
 
struct BcRot {
  BcRot(std::string name, std::vector<double> &attr, Range &faces);
  std::string blockName;
  Range faces;
  VectorDouble3 vals;
  double theta;
};
typedef std::vector<BcRot> BcRotVec;
 
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;
};
typedef std::vector<TractionBc> TractionBcVec;
 
struct OpJacobian : public UserDataOperator {
 
  OpJacobian(const int tag, const bool eval_rhs, const bool eval_lhs,
             boost::shared_ptr<DataAtIntegrationPts> data_ptr,
             boost::shared_ptr<PhysicalEquations> physics_ptr)
      : UserDataOperator(NOSPACE, OPSPACE), tAg(tag), evalRhs(eval_rhs),
        evalLhs(eval_lhs), dataAtPts(data_ptr), physicsPtr(physics_ptr) {}
 
  virtual MoFEMErrorCode evaluateRhs(EntData &data) = 0;
  virtual MoFEMErrorCode evaluateLhs(EntData &data) = 0;
 
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
protected:
  OpJacobian() : UserDataOperator(NOSPACE, OPSPACE) {};
 
  int tAg = -1; ///< adol-c tape
  bool evalRhs = false;
  bool evalLhs = false;
 
  boost::shared_ptr<DataAtIntegrationPts>
      dataAtPts; ///< data at integration pts
  boost::shared_ptr<PhysicalEquations>
      physicsPtr; ///< material physical equations
};
 
template <typename T> struct OpAssembleBasic : public T {
 
  using ScaleOff = boost::function<double()>;
  const bool assembleSymmetry;
 
  boost::shared_ptr<DataAtIntegrationPts>
      dataAtPts; ///< data at integration pts
 
  OpAssembleBasic(const std::string &field_name,
                  boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                  const char type)
      : T(field_name, type), dataAtPts(data_ptr), assembleSymmetry(false) {}
 
  OpAssembleBasic(
      std::string row_field, std::string col_field,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr, const char type,
      const bool assemble_symmetry, ScaleOff scale_off = []() { return 1; })
      : T(row_field, col_field, type, false), dataAtPts(data_ptr),
        assembleSymmetry(assemble_symmetry), scaleOff(scale_off) {}
 
  VectorDouble nF; ///< local right hand side vector
  MatrixDouble K;  ///< local tangent matrix
  MatrixDouble transposeK;
 
  ScaleOff scaleOff;
 
  virtual MoFEMErrorCode integrate(EntData &data) {
    MoFEMFunctionBegin;
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
    MoFEMFunctionReturn(0);
  }
 
  virtual MoFEMErrorCode integrate(int row_side, EntityType row_type,
                                   EntData &data) {
    MoFEMFunctionBegin;
    CHKERR integrate(data);
    MoFEMFunctionReturn(0);
  }
 
  virtual MoFEMErrorCode assemble(EntData &data) {
    MoFEMFunctionBegin;
    double *vec_ptr = &*nF.begin();
    int nb_dofs = data.getIndices().size();
    int *ind_ptr = &*data.getIndices().begin();
    CHKERR VecSetValues(this->getTSf(), nb_dofs, ind_ptr, vec_ptr, ADD_VALUES);
    MoFEMFunctionReturn(0);
  }
 
  virtual MoFEMErrorCode assemble(int row_side, EntityType row_type,
                                  EntData &data) {
    MoFEMFunctionBegin;
    CHKERR assemble(data);
    MoFEMFunctionReturn(0);
  }
 
  virtual MoFEMErrorCode integrate(EntData &row_data, EntData &col_data) {
    MoFEMFunctionBegin;
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
    MoFEMFunctionReturn(0);
  }
 
  virtual MoFEMErrorCode assemble(int row_side, int col_side,
                                  EntityType row_type, EntityType col_type,
                                  EntData &row_data, EntData &col_data) {
    MoFEMFunctionBegin;
 
    if (assembleSymmetry) {
      const auto row_nb_dofs = row_data.getIndices().size();
      const auto col_nb_dofs = col_data.getIndices().size();
      transposeK.resize(col_nb_dofs, row_nb_dofs, false);
      noalias(transposeK) = trans(K);
      transposeK *= scaleOff();
    }
 
    CHKERR MatSetValues<AssemblyTypeSelector<A>>(this->getTSB(), row_data,
                                                 col_data, K, ADD_VALUES);
    if (assembleSymmetry) {
      CHKERR MatSetValues<AssemblyTypeSelector<A>>(
          this->getTSB(), col_data, row_data, transposeK, ADD_VALUES);
    }
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
    MoFEMFunctionBegin;
    if (data.getIndices().empty())
      MoFEMFunctionReturnHot(0);
    nF.resize(data.getIndices().size(), false);
    nF.clear();
    CHKERR integrate(side, type, data);
    CHKERR assemble(side, type, data);
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
                        EntityType col_type, EntData &row_data,
                        EntData &col_data) {
    MoFEMFunctionBegin;
    if (row_data.getIndices().empty())
      MoFEMFunctionReturnHot(0);
    if (col_data.getIndices().empty())
      MoFEMFunctionReturnHot(0);
    K.resize(row_data.getIndices().size(), col_data.getIndices().size(), false);
    K.clear();
    CHKERR integrate(row_data, col_data);
    CHKERR assemble(row_side, col_side, row_type, col_type, row_data, col_data);
    MoFEMFunctionReturn(0);
  }
};
 
struct OpAssembleVolume : public OpAssembleBasic<VolUserDataOperator> {
  OpAssembleVolume(const std::string &field,
                   boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                   const char type)
      : OpAssembleBasic<VolUserDataOperator>(field, data_ptr, type) {}
 
  OpAssembleVolume(std::string row_field, std::string col_field,
                   boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                   const char type, const bool assemble_symmetry)
      : OpAssembleBasic<VolUserDataOperator>(row_field, col_field, data_ptr,
                                             type, assemble_symmetry) {}
};
 
struct OpAssembleFace : public OpAssembleBasic<FaceUserDataOperator> {
  OpAssembleFace(const std::string &field,
                 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                 const char type)
      : OpAssembleBasic<FaceUserDataOperator>(field, data_ptr, type) {}
 
  OpAssembleFace(std::string row_field, std::string col_field,
                 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                 const char type, const bool assemble_symmetry)
      : OpAssembleBasic<FaceUserDataOperator>(row_field, col_field, data_ptr,
                                              type, assemble_symmetry) {}
};
 
// Add operator to calculate Eshelby tensor
 
struct OpCalculateEshelbyStress
    : public ForcesAndSourcesCore::UserDataOperator {
  OpCalculateEshelbyStress(boost::shared_ptr<DataAtIntegrationPts> data_ptr)
      : ForcesAndSourcesCore::UserDataOperator(NOSPACE, OPSPACE),
        dataAtPts(data_ptr) {}
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
private:
  boost::shared_ptr<DataAtIntegrationPts>
      dataAtPts; ///< data at integration pts
};
 
struct OpCalculateRotationAndSpatialGradient : public VolUserDataOperator {
 
  OpCalculateRotationAndSpatialGradient(
      boost::shared_ptr<DataAtIntegrationPts> data_ptr)
      : VolUserDataOperator(NOSPACE, OPSPACE), dataAtPts(data_ptr) {}
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
private:
  boost::shared_ptr<DataAtIntegrationPts>
      dataAtPts; ///< data at integration pts
};
 
struct OpSpatialEquilibrium : public OpAssembleVolume {
  const double alphaW;
  const double alphaRho;
  OpSpatialEquilibrium(const std::string &field_name,
                       boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                       const double alpha, const double rho)
      : OpAssembleVolume(field_name, data_ptr, OPROW), alphaW(alpha),
        alphaRho(rho) {}
  MoFEMErrorCode integrate(EntData &data);
};
 
struct OpSpatialRotation : public OpAssembleVolume {
  OpSpatialRotation(const std::string &field_name,
                    boost::shared_ptr<DataAtIntegrationPts> data_ptr)
      : OpAssembleVolume(field_name, data_ptr, OPROW) {}
  MoFEMErrorCode integrate(EntData &data);
};
 
struct OpSpatialConsistencyP : public OpAssembleVolume {
  OpSpatialConsistencyP(const std::string &field_name,
                        boost::shared_ptr<DataAtIntegrationPts> data_ptr)
      : OpAssembleVolume(field_name, data_ptr, OPROW) {}
  MoFEMErrorCode integrate(EntData &data);
 
private:
};
 
struct OpSpatialConsistencyBubble : public OpAssembleVolume {
  OpSpatialConsistencyBubble(const std::string &field_name,
                             boost::shared_ptr<DataAtIntegrationPts> data_ptr)
      : OpAssembleVolume(field_name, data_ptr, OPROW) {}
  MoFEMErrorCode integrate(EntData &data);
 
private:
};
 
struct OpSpatialConsistencyDivTerm : public OpAssembleVolume {
  OpSpatialConsistencyDivTerm(const std::string &field_name,
                              boost::shared_ptr<DataAtIntegrationPts> data_ptr)
      : OpAssembleVolume(field_name, data_ptr, OPROW) {}
  MoFEMErrorCode integrate(EntData &data);
};
 
template <AssemblyType A, IntegrationType I> struct OpDispBcImpl;
 
template <AssemblyType A>
struct OpDispBcImpl<A, GAUSS>
    : public FormsIntegrators<FaceUserDataOperator>::Assembly<A>::OpBrokenBase {
 
  using OP = typename FormsIntegrators<FaceUserDataOperator>::Assembly<
      A>::OpBrokenBase;
 
  OpDispBcImpl(
      boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
      boost::shared_ptr<BcDispVec> &bc_disp_ptr,
      std::vector<boost::shared_ptr<ScalingMethod>> smv,
      boost::shared_ptr<Range> ents_ptr = nullptr)
      : OP(broken_base_side_data, ents_ptr), bcDispPtr(bc_disp_ptr),
        scalingMethodsVec(smv) {}
 
protected:
  MoFEMErrorCode iNtegrate(EntData &data);
  boost::shared_ptr<BcDispVec> bcDispPtr;
  std::vector<boost::shared_ptr<ScalingMethod>> scalingMethodsVec;
};
 
struct OpDispBc : public OpDispBcImpl<PETSC, GAUSS> {
  using OP = OpDispBcImpl<PETSC, GAUSS>;
  using OP::OpDispBcImpl;
 
protected:
  MoFEMErrorCode iNtegrate(EntData &data);
};
 
template <AssemblyType A, IntegrationType I> struct OpRotationBcImpl;
 
template <AssemblyType A>
struct OpRotationBcImpl<A, GAUSS>
    : public FormsIntegrators<FaceUserDataOperator>::Assembly<A>::OpBrokenBase {
 
  using OP = typename FormsIntegrators<FaceUserDataOperator>::Assembly<
      A>::OpBrokenBase;
 
  OpRotationBcImpl(
      boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
      boost::shared_ptr<BcRotVec> &bc_rot_ptr,
      std::vector<boost::shared_ptr<ScalingMethod>> smv,
      boost::shared_ptr<Range> ents_ptr = nullptr)
      : OP(broken_base_side_data, ents_ptr), bcRotPtr(bc_rot_ptr),
        scalingMethodsVec(smv) {}
 
protected:
  MoFEMErrorCode iNtegrate(EntData &data);
  boost::shared_ptr<BcRotVec> bcRotPtr;
  std::vector<boost::shared_ptr<ScalingMethod>> scalingMethodsVec;
};
 
struct OpRotationBc : public OpRotationBcImpl<PETSC, GAUSS> {
  using OP = OpRotationBcImpl<PETSC, GAUSS>;
  using OP::OpRotationBcImpl;
 
protected:
  MoFEMErrorCode iNtegrate(EntData &data);
};
 
struct OpBrokenTractionBc
    : public FormsIntegrators<FaceUserDataOperator>::Assembly<PETSC>::OpBase {
  OpBrokenTractionBc(const std::string row_field,
                     boost::shared_ptr<TractionBcVec> bc_data)
      : FormsIntegrators<FaceUserDataOperator>::Assembly<PETSC>::OpBase(
            row_field, row_field, FaceUserDataOperator::OPROW),
        bcData(bc_data) {}
 
  MoFEMErrorCode iNtegrate(EntData &data);
 
protected:
  boost::shared_ptr<TractionBcVec> bcData;
};
 
struct OpSpatialEquilibrium_dw_dP : public OpAssembleVolume {
  OpSpatialEquilibrium_dw_dP(std::string row_field, std::string col_field,
                             boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                             const bool assemble_off = false)
      : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
                         assemble_off) {
    sYmm = false;
  }
  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
};
 
struct OpSpatialEquilibrium_dw_dw : public OpAssembleVolume {
  const double alphaW;
  const double alphaRho;
  OpSpatialEquilibrium_dw_dw(std::string row_field, std::string col_field,
                             boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                             const double alpha, const double rho)
      : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL, false),
        alphaW(alpha), alphaRho(rho) {
    sYmm = true;
  }
  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
};
 
struct OpSpatialPhysical_du_dP : public OpAssembleVolume {
  OpSpatialPhysical_du_dP(std::string row_field, std::string col_field,
                          boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                          const bool assemble_off = false)
      : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
                         assemble_off) {
    sYmm = false;
  }
 
  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
};
 
struct OpSpatialPhysical_du_dBubble : public OpAssembleVolume {
  OpSpatialPhysical_du_dBubble(std::string row_field, std::string col_field,
                               boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                               const bool assemble_off = false)
      : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
                         assemble_off) {
    sYmm = false;
  }
 
  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
};
 
struct OpSpatialPhysical_du_domega : public OpAssembleVolume {
  OpSpatialPhysical_du_domega(std::string row_field, std::string col_field,
                              boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                              const bool assemble_off = false)
      : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
                         assemble_off) {
    sYmm = false;
  }
  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
};
 
struct OpSpatialRotation_domega_dP : public OpAssembleVolume {
  OpSpatialRotation_domega_dP(std::string row_field, std::string col_field,
                              boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                              const bool assemble_off)
      : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
                         assemble_off) {
    sYmm = false;
  }
  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
};
 
struct OpSpatialRotation_domega_dBubble : public OpAssembleVolume {
  OpSpatialRotation_domega_dBubble(
      std::string row_field, std::string col_field,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr, const bool assemble_off)
      : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
                         assemble_off) {
    sYmm = false;
  }
  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
};
 
struct OpSpatialRotation_domega_domega : public OpAssembleVolume {
  OpSpatialRotation_domega_domega(
      std::string row_field, std::string col_field,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr)
      : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL, false) {
    sYmm = false;
  }
  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
};
 
struct OpSpatialConsistency_dP_dP : public OpAssembleVolume {
  OpSpatialConsistency_dP_dP(std::string row_field, std::string col_field,
                             boost::shared_ptr<DataAtIntegrationPts> data_ptr)
      : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL, false) {
    sYmm = false;
  }
  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
 
private:
  template <int S>
  MoFEMErrorCode integrateImpl(EntData &row_data, EntData &col_data);
};
 
struct OpSpatialConsistency_dBubble_dBubble : public OpAssembleVolume {
  OpSpatialConsistency_dBubble_dBubble(
      std::string row_field, std::string col_field,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr)
      : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL, false) {
    sYmm = false;
  }
  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
 
private:
  template <int S>
  MoFEMErrorCode integrateImpl(EntData &row_data, EntData &col_data);
};
 
struct OpSpatialConsistency_dBubble_dP : public OpAssembleVolume {
  OpSpatialConsistency_dBubble_dP(
      std::string row_field, std::string col_field,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr)
      : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL, true) {
    sYmm = false;
  }
  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
 
private:
  template <int S>
  MoFEMErrorCode integrateImpl(EntData &row_data, EntData &col_data);
};
 
struct OpSpatialConsistency_dP_domega : public OpAssembleVolume {
  OpSpatialConsistency_dP_domega(
      std::string row_field, std::string col_field,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr, const bool assemble_off)
      : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
                         assemble_off) {
    sYmm = false;
  }
  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
 
private:
};
 
struct OpSpatialConsistency_dBubble_domega : public OpAssembleVolume {
  OpSpatialConsistency_dBubble_domega(
      std::string row_field, std::string col_field,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr, const bool assemble_off)
      : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
                         assemble_off) {
    sYmm = false;
  }
  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
 
private:
};
 
struct OpPostProcDataStructure
    : public VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator {
 
  using OP = VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator;
 
  moab::Interface &postProcMesh;
  std::vector<EntityHandle> &mapGaussPts;
  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
 
  OpPostProcDataStructure(moab::Interface &post_proc_mesh,
                          std::vector<EntityHandle> &map_gauss_pts,
                          boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                          int sense)
      : OP(NOSPACE, UserDataOperator::OPSPACE), postProcMesh(post_proc_mesh),
        mapGaussPts(map_gauss_pts), dataAtPts(data_ptr), tagSense(sense) {}
 
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
private:
  int tagSense;
};
 
struct OpSpatialPrj : public OpAssembleVolume {
  OpSpatialPrj(std::string row_field,
               boost::shared_ptr<DataAtIntegrationPts> data_ptr)
      : OpAssembleVolume(row_field, data_ptr, OPROW) {}
  MoFEMErrorCode integrate(EntData &row_data);
};
 
struct OpSpatialPrj_dx_dx : public OpAssembleVolume {
  OpSpatialPrj_dx_dx(std::string row_field, std::string col_field,
                     boost::shared_ptr<DataAtIntegrationPts> data_ptr)
      : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL, false) {
    // FIXME: That is symmetric
    sYmm = false;
  }
  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
};
 
struct OpSpatialPrj_dx_dw : public OpAssembleVolume {
  OpSpatialPrj_dx_dw(std::string row_field, std::string col_field,
                     boost::shared_ptr<DataAtIntegrationPts> data_ptr)
      : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL, false) {
    sYmm = false;
  }
  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
};
 
struct EshelbianCore : public MoFEM::UnknownInterface {
 
  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 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 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 crackElement = "CRACK";
  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 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 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,
      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);
 
  MoFEMErrorCode solveElastic(TS ts, Vec x);
 
  MoFEMErrorCode postProcessResults(const int tag, const std::string file);
  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;
  };
 
  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<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
};
 
template <int FE_DIM, int PROBLEM_DIM, int SPACE_DIM> struct AddHOOps;
 
template <> struct AddHOOps<2, 3, 3> {
  AddHOOps() = delete;
  static MoFEMErrorCode
  add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
      std::vector<FieldSpace> space, std::string geom_field_name,
      boost::shared_ptr<Range> crack_front_edges_ptr);
};
 
template <> struct AddHOOps<3, 3, 3> {
  AddHOOps() = delete;
  static MoFEMErrorCode
  add(boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
      std::vector<FieldSpace> space, std::string geom_field_name,
      boost::shared_ptr<Range> crack_front_edges_ptr);
};
 
} // namespace EshelbianPlasticity
 
#endif //__ESHELBIAN_PLASTICITY_HPP__