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
 *
 * \todo Implementation of plasticity
 */
 
constexpr int SPACE_DIM = 3;
constexpr AssemblyType A = AssemblyType::SCHUR;
 
#ifndef __ESHELBIAN_PLASTICITY_HPP__
#define __ESHELBIAN_PLASTICITY_HPP__
 
namespace EshelbianPlasticity {
 
struct ContactTree;
 
enum RotSelector { SMALL_ROT, MODERATE_ROT, LARGE_ROT };
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;
 
  MatrixDouble diffStretchTensorAtPts;
  VectorDouble detStretchTensorAtPts;
  MatrixDouble detStretchTensorAtPts_du;
  MatrixDouble logStretchDotTensorAtPts;
  MatrixDouble rotAxisAtPts;
  MatrixDouble rotAxisDotAtPts;
  MatrixDouble WAtPts;
  MatrixDouble W0AtPts;
  MatrixDouble GAtPts;
  MatrixDouble G0AtPts;
  MatrixDouble wH1AtPts;
  MatrixDouble contactL2AtPts;
  MatrixDouble wGradH1AtPts;
  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 phiAtPts;
  MatrixDouble flowL2AtPts;
 
  MatrixDouble P_du;
 
  MatrixDouble eigenVals;
  MatrixDouble eigenVecs;
  VectorInt nbUniq;
 
  MatrixDouble matD;
  MatrixDouble matAxiatorD;
  MatrixDouble matDeviatorD;
 
  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 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 getContactL2AtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &contactL2AtPts);
  }
 
  inline MatrixPtr getSmallWGradH1AtPts() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(),
                                           &wGradH1AtPts);
  } 
 
  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;
 
  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(
      const std::string &row_field, const 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<SchurL2Mats>(this->getTSA(), row_data, col_data, K,
                                     ADD_VALUES);
    if (assembleSymmetry) {
      CHKERR MatSetValues<SchurL2Mats>(this->getTSA(), 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(const std::string &row_field, const 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(const std::string &row_field, const 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) {}
};
 
struct OpCalculateRotationAndSpatialGradient : public VolUserDataOperator {
  boost::shared_ptr<DataAtIntegrationPts>
      dataAtPts; ///< data at integration pts
  OpCalculateRotationAndSpatialGradient(
      boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
      : VolUserDataOperator(NOSPACE, OPSPACE), dataAtPts(data_ptr) {}
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
};
 
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 OpSpatialPhysical : public OpAssembleVolume {
 
  OpSpatialPhysical(const std::string &field_name,
                    boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
                    const double alpha_u)
      : OpAssembleVolume(field_name, data_ptr, OPROW), alphaU(alpha_u) {
 
    CHK_MOAB_THROW(PetscOptionsGetBool(PETSC_NULL, "", "-poly_convex",
                                       &polyConvex, PETSC_NULL),
                   "get ployconvex option failed");
  }
 
  MoFEMErrorCode integrate(EntData &data);
 
  MoFEMErrorCode integratePiola(EntData &data);
  MoFEMErrorCode integrateHencky(EntData &data);
  MoFEMErrorCode integratePolyconvexHencky(EntData &data);
 
private:
  const double alphaU;
  PetscBool polyConvex = PETSC_FALSE;
};
 
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);
};
 
struct OpDispBc : public OpAssembleFace {
  boost::shared_ptr<BcDispVec> bcDispPtr;
  OpDispBc(const std::string &field_name,
           boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
           boost::shared_ptr<BcDispVec> &bc_disp_ptr,
           std::vector<boost::shared_ptr<ScalingMethod>> smv)
      : OpAssembleFace(field_name, data_ptr, OPROW), bcDispPtr(bc_disp_ptr),
        scalingMethodsVec(smv) {}
  MoFEMErrorCode integrate(EntData &data);
 
protected:
  std::vector<boost::shared_ptr<ScalingMethod>> scalingMethodsVec;
};
 
struct OpDispBc_dx : public OpAssembleFace {
  boost::shared_ptr<BcDispVec> bcDispPtr;
  OpDispBc_dx(const std::string &row_field_name,
              const std::string &col_field_name,
              boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
              boost::shared_ptr<BcDispVec> &bc_disp_ptr)
      : OpAssembleFace(row_field_name, col_field_name, data_ptr, OPROWCOL,
                       false),
        bcDispPtr(bc_disp_ptr) {}
  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
};
 
struct OpRotationBc : public OpAssembleFace {
  boost::shared_ptr<BcRotVec> bcRotPtr;
  OpRotationBc(const std::string &field_name,
               boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
               boost::shared_ptr<BcRotVec> &bc_rot_ptr)
      : OpAssembleFace(field_name, data_ptr, OPROW), bcRotPtr(bc_rot_ptr) {}
  MoFEMErrorCode integrate(EntData &data);
};
 
struct OpRotationBc_dx : public OpAssembleFace {
  boost::shared_ptr<BcRotVec> bcRotPtr;
  OpRotationBc_dx(const std::string &row_field_name,
                  const std::string &col_field_name,
                  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
                  boost::shared_ptr<BcRotVec> &bc_rot_ptr)
      : OpAssembleFace(row_field_name, col_field_name, data_ptr, OPROWCOL,
                       false),
        bcRotPtr(bc_rot_ptr) {}
  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
};
 
struct FeTractionBc;
 
struct OpTractionBc : public FaceUserDataOperator {
  OpTractionBc(std::string row_field, FeTractionBc &bc_fe)
      : FaceUserDataOperator(row_field, FaceUserDataOperator::OPROW),
        bcFe(bc_fe) {}
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
protected:
  FeTractionBc &bcFe;
  MatrixDouble matPP;
  MatrixDouble vecPv;
};
 
struct FeTractionBc : public FEMethod {
  int nbConstrainsDofs;
  FeTractionBc(MoFEM::Interface &m_field, const std::string field_name,
               boost::shared_ptr<TractionBcVec> &bc)
      : FEMethod(), mField(m_field), bcData(bc), fieldName(field_name) {}
  MoFEMErrorCode preProcess();
  MoFEMErrorCode postProcess() { return 0; }
 
  friend MoFEMErrorCode OpTractionBc::doWork(int side, EntityType type,
                                             EntData &data);
 
protected:
  MoFEM::Interface &mField;
  boost::shared_ptr<TractionBcVec> bcData;
  std::string fieldName;
};
 
struct OpSpatialEquilibrium_dw_dP : public OpAssembleVolume {
  OpSpatialEquilibrium_dw_dP(const std::string &row_field,
                             const 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(const std::string &row_field,
                             const 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_du : public OpAssembleVolume {
  const double alphaU;
  OpSpatialPhysical_du_du(const std::string &row_field,
                          const std::string &col_field,
                          boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
                          const double alpha)
      : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL, false),
        alphaU(alpha) {
    sYmm = false;
 
    CHK_MOAB_THROW(PetscOptionsGetBool(PETSC_NULL, "", "-poly_convex",
                                       &polyConvex, PETSC_NULL),
                   "get ployconvex option failed");
  }
  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
  MoFEMErrorCode integratePiola(EntData &row_data, EntData &col_data);
  MoFEMErrorCode integrateHencky(EntData &row_data, EntData &col_data);
  MoFEMErrorCode integratePolyconvexHencky(EntData &row_data,
                                           EntData &col_data);
private:
  PetscBool polyConvex = PETSC_FALSE;
 
  MatrixDouble dP;
};
 
struct OpSpatialPhysical_du_dP : public OpAssembleVolume {
  OpSpatialPhysical_du_dP(const std::string &row_field,
                          const 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(
      const std::string &row_field, const 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(const std::string &row_field,
                              const 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_dx : public OpAssembleVolume {
  OpSpatialPhysical_du_dx(const std::string &row_field,
                          const 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(const std::string &row_field,
                              const 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(
      const std::string &row_field, const 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(
      const std::string &row_field, const 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_domega : public OpAssembleVolume {
  OpSpatialConsistency_dP_domega(
      const std::string &row_field, const 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(
      const std::string &row_field, const 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 FaceUserDataOperator {
 
  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)
      : FaceUserDataOperator(NOSPACE, UserDataOperator::OPSPACE),
        postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
        dataAtPts(data_ptr) {}
 
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
};
 
struct OpSpatialPrj : public OpAssembleVolume {
  OpSpatialPrj(const 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(const std::string &row_field, const 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(const std::string &row_field, const 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 OpCalculateStrainEnergy : public VolUserDataOperator {
 
  OpCalculateStrainEnergy(const std::string field_name,
                          boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
                          boost::shared_ptr<double> &e)
      : VolUserDataOperator(field_name, VolUserDataOperator::OPROW),
        dataAtPts(data_at_pts), energy(e) {}
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data);
 
private:
  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
  boost::shared_ptr<double> energy;
};
 
struct EshelbianCore : public MoFEM::UnknownInterface {
 
  static enum RotSelector rotSelector;
  static enum RotSelector gradApperoximator;
  static enum StretchSelector stretchSelector;
 
  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(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<FaceElementForcesAndSourcesCore> contactRhs;
  boost::shared_ptr<ContactTree> contactTreeRhs; ///< Make a contact tree
 
  SmartPetscObj<DM> dM;           ///< Coupled problem all fields
  SmartPetscObj<DM> dmElastic;    ///< Elastic problem
  SmartPetscObj<DM> dmPrjSpatial; ///< Projection spatial displacement
 
  const std::string piolaStress;
  const std::string eshelbyStress;
  const std::string spatialL2Disp;
  const std::string spatialH1Disp;
  const std::string contactDisp;
  const std::string materialL2Disp;
  const std::string stretchTensor;
  const std::string rotAxis;
  const std::string materialGradient;
  const std::string tauField;
  const std::string lambdaField;
  const std::string bubbleField;
 
  const std::string elementVolumeName;
  const std::string naturalBcElement;
  const std::string essentialBcElement;
  const std::string skinElement;
  const std::string contactElement;
 
  EshelbianCore(MoFEM::Interface &m_field);
  virtual ~EshelbianCore();
 
  int spaceOrder;
  int spaceH1Order;
  int materialOrder;
  double alphaU;
  double alphaW;
  double alphaRho;
  double precEps;
  double precEpsOmega;
  double precEpsW;
  double precEpsContactDisp;
 
  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 addVolumeFiniteElement(const EntityHandle meshset = 0);
  MoFEMErrorCode addBoundaryFiniteElement(const EntityHandle meshset = 0);
  MoFEMErrorCode addDMs(const BitRefLevel bit = BitRefLevel().set(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_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 setContactElementOps(
 
      boost::shared_ptr<ContactTree> &fe_contact_tree,
      boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_rhs,
      boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_lhs
 
  );
 
  MoFEMErrorCode setElasticElementOps(const int tag);
  MoFEMErrorCode setElasticElementToTs(DM dm);
 
  MoFEMErrorCode solveElastic(TS ts, Mat m, Vec f, 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, SmartPetscObj<Mat> m, SmartPetscObj<Mat> p,
        EshelbianCore *ep_core_ptr
 
    );
    virtual MoFEMErrorCode setUp(KSP solver) = 0;
 
    virtual MoFEMErrorCode preProc() = 0;
    virtual MoFEMErrorCode postProc() = 0;
 
  protected:
    SetUpSchur() = default;
  };
};
 
} // namespace EshelbianPlasticity
 
#endif //__ESHELBIAN_PLASTICITY_HPP__