EshelbianOperators.hpp

Go to the documentation of this file.

/**
 * @file EshelbianOperators.hpp
 * @author your name (you@domain.com)
 * @brief 
 * @version 0.1
 * @date 2024-12-31
 * 
 * @copyright Copyright (c) 2024
 * 
 */
 
 
 
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> {
  using OpAssembleBasic<VolUserDataOperator>::OpAssembleBasic;
};
 
struct OpAssembleFace : public OpAssembleBasic<FaceUserDataOperator> {
  using OpAssembleBasic<FaceUserDataOperator>::OpAssembleBasic;
};
 
struct OpAssembleVolumePositiveDefine : OpAssembleVolume {
  using OpAssembleVolume::OpAssembleVolume;
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
    MoFEMFunctionBegin;
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
            "OpAssembleVolumePositiveDefine not implemented");
    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);
 
    auto symmetrize = [](MatrixDouble &&mat) {
      return (mat + trans(mat)) / 2.0;
    };
    auto check_positive_definite = [](MatrixDouble &&mat) {
      int n = mat.size1();
      return (lapack_dpotrf('U', n, mat.data().data(), n) == 0);
    };
    auto frobenius_norm = [](MatrixDouble &mat) {
      double sum = 0.0;
      for (std::size_t i = 0; i < mat.size1(); ++i)
        for (std::size_t j = 0; j < mat.size2(); ++j)
          sum += mat(i, j) * mat(i, j);
      return std::sqrt(sum);
    };
 
    auto higham_method_one_pass = [&](MatrixDouble &mat) {
      int n = mat.size1();
      const int lda = n;
      const int size = (n + 2) * n;
      int lwork = size;
      double *work = new double[size];
 
      VectorDouble eig(n, 0.0);
      MatrixDouble eigen_vec = prod(mat, trans(mat));
      if (lapack_dsyev('V', 'U', n, eigen_vec.data().data(), lda,
                       eig.data().data(), work, lwork) > 0)
        CHK_THROW_MESSAGE(MOFEM_INVALID_DATA,
                          "The algorithm failed to compute eigenvalues.");
      delete[] work;
 
      double max_e = std::abs(eig[0]);
      for (auto e : eig) {
        max_e = std::max(std::abs(e), max_e);
      }
      constexpr auto eps = std::numeric_limits<float>::epsilon();
      auto min_e = std::sqrt(max_e) * eps;
      for (auto &e : eig) {
        e = std::sqrt(std::abs(e));
        e = std::max(e, min_e);
        e = std::min(e, static_cast<double>(std::numeric_limits<float>::max()));
      }
 
      auto diagonal_matrix =
          ublas::diagonal_matrix<double, ublas::row_major,
                                 VecAllocator<double>>(n, eig.data());
      MatrixDouble symm_mat;
      symm_mat = prod(trans(eigen_vec), diagonal_matrix);
      symm_mat = prod(symm_mat, eigen_vec);
 
      MatrixDouble hat_m = (mat + symm_mat) / 2.0;
 
      return hat_m;
    };
 
    auto hat_k = higham_method_one_pass(K);
 
    // MatrixDouble tmp_mat = K - hat_k;
    // auto nrm2 = frobenius_norm(tmp_mat);
    // MOFEM_LOG("EPSELF", Sev::inform) <<"Higham method one pass norm: " << nrm2;
    // SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
    //         "Higham method one pass not implemented");
 
    K.swap(hat_k);
 
    CHKERR assemble(row_side, col_side, row_type, col_type, row_data, col_data);
    MoFEMFunctionReturn(0);
  }
 
private:
};
 
// 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, double alpha_omega = 0)
      : VolUserDataOperator(NOSPACE, OPSPACE), dataAtPts(data_ptr),
        alphaOmega(alpha_omega) {}
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
private:
  boost::shared_ptr<DataAtIntegrationPts>
      dataAtPts; ///< data at integration pts
  double alphaOmega;
};
 
struct OpCalculateTractionFromSideEl : public SideEleOp {
  OpCalculateTractionFromSideEl(
      boost::shared_ptr<DataAtIntegrationPts> data_ptr)
      : SideEleOp(NOSPACE, SideEleOp::OPSPACE), dataAtPts(data_ptr) {
    std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
    for (auto t = moab::CN::TypeDimensionMap[SPACE_DIM].first;
         t <= moab::CN::TypeDimensionMap[SPACE_DIM].second; ++t)
      doEntities[t] = true;
    sYmm = false;
  }
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
private:
  boost::shared_ptr<DataAtIntegrationPts>
      dataAtPts; ///< data at integration pts
};
 
struct OpCalculateReactionForces : public FaceUserDataOperator {
  OpCalculateReactionForces(boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                            std::string block_name, Range block_entities,
                            std::array<double, 6> &reaction_vec)
      : FaceUserDataOperator(NOSPACE, OPSPACE), dataAtPts(data_ptr),
        blockName(block_name), blockEntities(block_entities),
        reactionVec(reaction_vec) {}
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
private:
  boost::shared_ptr<DataAtIntegrationPts>
      dataAtPts; ///< data at integration pts
  Range blockEntities;
  std::string blockName;
  std::array<double, 6> &reactionVec;
};
 
struct OpSpatialEquilibrium : public OpAssembleVolume {
  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);
 
private:
  const double alphaW;
  const double alphaRho;
 
};
 
struct OpSpatialRotation : public OpAssembleVolume {
  OpSpatialRotation(const std::string &field_name,
                    boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                    double alpha_omega)
      : OpAssembleVolume(field_name, data_ptr, OPROW), alphaOmega(alpha_omega) {
  }
  MoFEMErrorCode integrate(EntData &data);
 
private:
  double alphaOmega;
};
 
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 <int INTERP_ORDER>
struct OpGetInternalStress : public UserDataOperator {
  OpGetInternalStress(boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                      std::string tag_name = "")
      : UserDataOperator(H1, OPLAST), dataAtPts(data_ptr), tagName(tag_name) {
    std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
    doEntities[MBVERTEX] = true;
  }
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
private:
  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
  std::string tagName;
};
 
template <bool VOIGT> 
struct OpSpatialPhysicalInternalStress : public OpAssembleVolume {
  OpSpatialPhysicalInternalStress(
      const std::string &field_name,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr)
      : OpAssembleVolume(field_name, data_ptr, OPROW) {}
  MoFEMErrorCode integrate(EntData &row_data);
};
 
 
struct OpSpatialPhysicalExternalStrain : public OpAssembleVolume  {
 
  OpSpatialPhysicalExternalStrain(
      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,
      boost::shared_ptr<Range> ents_ptr = nullptr)
 
      : OpAssembleVolume(field_name, data_ptr, OPROW), dataAtPts{data_ptr}, externalStrainVecPtr(external_strain_vec_ptr), scalingMethodsMap{smv} {}
 
protected:
  MoFEMErrorCode integrate(EntData &row_data);
 
  boost::shared_ptr<DataAtIntegrationPts> dataAtPts; 
  
  boost::shared_ptr<ExternalStrainVec> externalStrainVecPtr;
  std::map<std::string, boost::shared_ptr<ScalingMethod>> scalingMethodsMap;
};
 
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::map<std::string, boost::shared_ptr<ScalingMethod>> smv,
      boost::shared_ptr<Range> ents_ptr = nullptr)
      : OP(broken_base_side_data, ents_ptr), bcDispPtr(bc_disp_ptr),
        scalingMethodsMap(smv) {}
 
protected:
  MoFEMErrorCode iNtegrate(EntData &data);
  boost::shared_ptr<BcDispVec> bcDispPtr;
  std::map<std::string, boost::shared_ptr<ScalingMethod>> scalingMethodsMap;
};
 
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::map<std::string, boost::shared_ptr<ScalingMethod>> smv,
      boost::shared_ptr<Range> ents_ptr = nullptr)
      : OP(broken_base_side_data, ents_ptr), bcRotPtr(bc_rot_ptr),
        scalingMethodsMap(smv) {}
 
protected:
  MoFEMErrorCode iNtegrate(EntData &data);
  boost::shared_ptr<BcRotVec> bcRotPtr;
  std::map<std::string, boost::shared_ptr<ScalingMethod>> scalingMethodsMap;
};
 
/**
 * @brief Apply rotation boundary condition
 *
 * The boundary condition is applied getting data from the SPATIAL_ROTATION_BC
 * blockset.
 *
 * Centre of rotation is given by the first three attributes of the blockset.
 * The fort parameter is angle of rotation
 *
 */
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,
      boost::shared_ptr<double> piola_scale_ptr,
      std::map<std::string, boost::shared_ptr<ScalingMethod>> smv)
      : FormsIntegrators<FaceUserDataOperator>::Assembly<PETSC>::OpBase(
            row_field, row_field, FaceUserDataOperator::OPROW),
        bcData(bc_data), piolaScalePtr(), scalingMethodsMap(smv) {}
 
  MoFEMErrorCode iNtegrate(EntData &data);
 
protected:
  boost::shared_ptr<TractionBcVec> bcData;
  boost::shared_ptr<double> piolaScalePtr;
  std::map<std::string, boost::shared_ptr<ScalingMethod>> scalingMethodsMap;
};
 
struct OpBrokenPressureBc
    : public FormsIntegrators<FaceUserDataOperator>::Assembly<PETSC>::OpBase {
  OpBrokenPressureBc(
      const std::string row_field, boost::shared_ptr<PressureBcVec> bc_data,
      boost::shared_ptr<double> piola_scale_ptr,
      std::map<std::string, boost::shared_ptr<ScalingMethod>> smv)
      : FormsIntegrators<FaceUserDataOperator>::Assembly<PETSC>::OpBase(
            row_field, row_field, FaceUserDataOperator::OPROW),
        bcData(bc_data), piolaScalePtr(), scalingMethodsMap(smv) {}
 
  MoFEMErrorCode iNtegrate(EntData &data);
 
protected:
  boost::shared_ptr<PressureBcVec> bcData;
  boost::shared_ptr<double> piolaScalePtr;
  std::map<std::string, boost::shared_ptr<ScalingMethod>> scalingMethodsMap;
};
 
struct OpBrokenAnalyticalTractionBc
    : public FormsIntegrators<FaceUserDataOperator>::Assembly<PETSC>::OpBase {
  OpBrokenAnalyticalTractionBc(
      const std::string row_field, boost::shared_ptr<AnalyticalTractionBcVec> bc_data,
      boost::shared_ptr<double> piola_scale_ptr,
      std::map<std::string, boost::shared_ptr<ScalingMethod>> smv)
      : FormsIntegrators<FaceUserDataOperator>::Assembly<PETSC>::OpBase(
            row_field, row_field, FaceUserDataOperator::OPROW),
        bcData(bc_data), piolaScalePtr(), scalingMethodsMap(smv) {}
 
  MoFEMErrorCode iNtegrate(EntData &data);
 
protected:
  boost::shared_ptr<AnalyticalTractionBcVec> bcData;
  boost::shared_ptr<double> piolaScalePtr;
  std::map<std::string, boost::shared_ptr<ScalingMethod>> scalingMethodsMap;
};
 
template <AssemblyType A, IntegrationType I> struct OpNormalDispBcRhsImpl;
template <AssemblyType A, IntegrationType I> struct OpNormalDispBcLhsImpl_dU;
template <AssemblyType A, IntegrationType I> struct OpNormalDispBcLhsImpl_dP;
 
template <AssemblyType A>
struct OpNormalDispBcRhsImpl<A, GAUSS>
    : public FormsIntegrators<FaceUserDataOperator>::Assembly<PETSC>::OpBase {
 
  using OP =
      typename FormsIntegrators<FaceUserDataOperator>::Assembly<PETSC>::OpBase;
 
  OpNormalDispBcRhsImpl(
      std::string row_field, boost::shared_ptr<MatrixDouble> hybrid_disp_ptr,
      boost::shared_ptr<MatrixDouble> piola_stress_ptr,
      boost::shared_ptr<NormalDisplacementBcVec> &bc_disp_ptr,
      std::map<std::string, boost::shared_ptr<ScalingMethod>> smv,
      boost::shared_ptr<Range> ents_ptr = nullptr)
      : FormsIntegrators<FaceUserDataOperator>::Assembly<PETSC>::OpBase(
            row_field, row_field, OP::OPROW),
        hybridDispPtr(hybrid_disp_ptr), piolaStressPtr(piola_stress_ptr),
        bcDispPtr(bc_disp_ptr), scalingMethodsMap(smv) {}
 
protected:
  MoFEMErrorCode iNtegrate(EntData &data);
  boost::shared_ptr<MatrixDouble> hybridDispPtr;
  boost::shared_ptr<MatrixDouble> piolaStressPtr;
  boost::shared_ptr<NormalDisplacementBcVec> bcDispPtr;
  std::map<std::string, boost::shared_ptr<ScalingMethod>> scalingMethodsMap;
};
 
template <AssemblyType A>
struct OpNormalDispBcLhsImpl_dU<A, GAUSS>
    : public FormsIntegrators<FaceUserDataOperator>::Assembly<A>::OpBase {
 
  using OP =
      typename FormsIntegrators<FaceUserDataOperator>::Assembly<A>::OpBase;
 
  OpNormalDispBcLhsImpl_dU(
      std::string row_field,
      boost::shared_ptr<NormalDisplacementBcVec> &bc_disp_ptr,
      std::map<std::string, boost::shared_ptr<ScalingMethod>> smv,
      boost::shared_ptr<Range> ents_ptr = nullptr)
      : OP(row_field, row_field, OP::OPROWCOL), bcDispPtr(bc_disp_ptr),
        scalingMethodsMap(smv) {}
 
protected:
  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
  boost::shared_ptr<NormalDisplacementBcVec> bcDispPtr;
  std::map<std::string, boost::shared_ptr<ScalingMethod>> scalingMethodsMap;
};
 
template <AssemblyType A>
struct OpNormalDispBcLhsImpl_dP<A, GAUSS>
    : public FormsIntegrators<FaceUserDataOperator>::Assembly<A>::OpBrokenBase {
 
  using OP = typename FormsIntegrators<FaceUserDataOperator>::Assembly<
      A>::OpBrokenBase;
 
  OpNormalDispBcLhsImpl_dP(
      std::string row_field,
      boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
      boost::shared_ptr<NormalDisplacementBcVec> &bc_disp_ptr,
      std::map<std::string, boost::shared_ptr<ScalingMethod>> smv,
      boost::shared_ptr<Range> ents_ptr = nullptr)
      : OP(row_field, broken_base_side_data, false, false, ents_ptr),
        bcDispPtr(bc_disp_ptr), scalingMethodsMap(smv) {}
 
protected:
  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
  boost::shared_ptr<NormalDisplacementBcVec> bcDispPtr;
  std::map<std::string, boost::shared_ptr<ScalingMethod>> scalingMethodsMap;
};
 
struct OpNormalDispRhsBc : public OpNormalDispBcRhsImpl<PETSC, GAUSS> {
  using OP = OpNormalDispBcRhsImpl<PETSC, GAUSS>;
  using OP::OpNormalDispBcRhsImpl;
 
protected:
  MoFEMErrorCode iNtegrate(EntData &data);
};
 
struct OpNormalDispLhsBc_dU : public OpNormalDispBcLhsImpl_dU<A, GAUSS> {
  using OP = OpNormalDispBcLhsImpl_dU<A, GAUSS>;
  using OP::OpNormalDispBcLhsImpl_dU;
 
protected:
  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
};
 
struct OpNormalDispLhsBc_dP : public OpNormalDispBcLhsImpl_dP<A, GAUSS> {
  using OP = OpNormalDispBcLhsImpl_dP<A, GAUSS>;
  using OP::OpNormalDispBcLhsImpl_dP;
 
protected:
  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
};
 
template <AssemblyType A, IntegrationType I> struct OpAnalyticalDispBcImpl;
 
template <AssemblyType A>
struct OpAnalyticalDispBcImpl<A, GAUSS>
    : public FormsIntegrators<FaceUserDataOperator>::Assembly<A>::OpBrokenBase {
 
  using OP = typename FormsIntegrators<FaceUserDataOperator>::Assembly<
      A>::OpBrokenBase;
 
  OpAnalyticalDispBcImpl(
      boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
      boost::shared_ptr<AnalyticalDisplacementBcVec> &bc_disp_ptr,
      std::map<std::string, boost::shared_ptr<ScalingMethod>> smv,
      boost::shared_ptr<Range> ents_ptr = nullptr)
      : OP(broken_base_side_data, ents_ptr), bcDispPtr(bc_disp_ptr),
        scalingMethodsMap(smv) {}
 
protected:
  MoFEMErrorCode iNtegrate(EntData &data);
  boost::shared_ptr<AnalyticalDisplacementBcVec> bcDispPtr;
  std::map<std::string, boost::shared_ptr<ScalingMethod>> scalingMethodsMap;
};
 
struct OpAnalyticalDispBc : public OpAnalyticalDispBcImpl<PETSC, GAUSS> {
  using OP = OpAnalyticalDispBcImpl<PETSC, GAUSS>;
  using OP::OpAnalyticalDispBcImpl;
 
protected:
  MoFEMErrorCode iNtegrate(EntData &data);
};
 
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 {
  using OpAssembleVolume::OpAssembleVolume;
  OpSpatialPhysical_du_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);
};
 
struct OpSpatialRotation_domega_du : public OpAssembleVolume {
  OpSpatialRotation_domega_du(std::string row_field, std::string col_field,
                              boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                              double alpha_omega)
      : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL, false) {
    sYmm = false;
  }
  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
private:
};
 
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, double alpha_omega)
      : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL, false),
        alphaOmega(alpha_omega) {
    sYmm = false;
  }
  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
 
private:
  double alphaOmega;
};
 
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 OpFaceSideMaterialForce
    : public VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator {
 
  // @note This is not on the side, since, integrate energy in the volume
  using OP = VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator;
 
  OpFaceSideMaterialForce(boost::shared_ptr<DataAtIntegrationPts> data_ptr)
      : OP(NOSPACE, OPSPACE), dataAtPts(data_ptr) {}
 
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
private:
  boost::shared_ptr<DataAtIntegrationPts>
      dataAtPts; ///< data at integration pts
};
 
struct OpFaceMaterialForce
    : public FaceElementForcesAndSourcesCore::UserDataOperator {
 
  // @note This is not on the side, since, integrate energy in the volume
  using OP = FaceElementForcesAndSourcesCore::UserDataOperator;
 
  OpFaceMaterialForce(boost::shared_ptr<DataAtIntegrationPts> data_ptr)
      : OP(NOSPACE, OPSPACE), dataAtPts(data_ptr) {}
 
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
private:
  boost::shared_ptr<DataAtIntegrationPts>
      dataAtPts; ///< data at integration pts
};
 
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<2, 2, 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);
};