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 bool eval_rhs, const bool eval_lhs,
             boost::shared_ptr<DataAtIntegrationPts> data_ptr,
             boost::shared_ptr<PhysicalEquations> physics_ptr)
      : UserDataOperator(NOSPACE, OPSPACE), 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(FieldSpace space, UserDataOperator::OpType op_type)
      : UserDataOperator(space, op_type) {}
 
  OpJacobian() : UserDataOperator(NOSPACE, OPSPACE) {};
 
  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), assembleSymmetry(false), dataAtPts(data_ptr) {}
 
  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),
        assembleSymmetry(assemble_symmetry), dataAtPts(data_ptr),
        scaleOff(scale_off) {}
 
  OpAssembleBasic(const FieldSpace space)
      : T(space, T::OPSPACE), assembleSymmetry(false) {}
 
  VectorDouble nF; ///< local right hand side vector
  MatrixDouble K;  ///< local tangent matrix
  MatrixDouble transposeK;
 
  ScaleOff scaleOff;
 
  virtual MoFEMErrorCode integrate(EntData &) {
    MoFEMFunctionBegin;
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
    MoFEMFunctionReturn(0);
  }
 
  virtual MoFEMErrorCode integrate(int, EntityType, EntData &data) {
    MoFEMFunctionBegin;
    CHKERR integrate(data);
    MoFEMFunctionReturn(0);
  }
 
  virtual MoFEMErrorCode assemble(EntData &data) {
    MoFEMFunctionBegin;
    double *vec_ptr = nF.data().data();
    int nb_dofs = data.getIndices().size();
    int *ind_ptr = data.getIndices().data().data();
    CHKERR VecSetValues(this->getTSf(), nb_dofs, ind_ptr, vec_ptr, ADD_VALUES);
    MoFEMFunctionReturn(0);
  }
 
  virtual MoFEMErrorCode assemble(int, EntityType, EntData &data) {
    MoFEMFunctionBegin;
    CHKERR assemble(data);
    MoFEMFunctionReturn(0);
  }
 
  virtual MoFEMErrorCode integrate(EntData &, EntData &) {
    MoFEMFunctionBegin;
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
    MoFEMFunctionReturn(0);
  }
 
  virtual MoFEMErrorCode assemble(int, int, EntityType, EntityType,
                                  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 OP = OpAssembleBasic<VolUserDataOperator>;
  using ScaleOff = typename OP::ScaleOff;
 
  using OP::OpAssembleBasic;
 
  MoFEMErrorCode assemble(int row_side, int col_side, EntityType row_type,
                          EntityType col_type, EntData &row_data,
                          EntData &col_data) {
    MoFEMFunctionBegin;
    CHKERR OP::assemble(row_side, col_side, row_type, col_type, row_data,
                        col_data);
    constexpr bool store_matrices = false;
    if constexpr (store_matrices) {
      auto &m = mapMatrix[{OP::rowFieldName, OP::colFieldName}];
      m.resize(K.size1(), K.size2(), false);
      noalias(m) = K;
      if (assembleSymmetry) {
        auto &m = mapMatrix[{OP::colFieldName, OP::rowFieldName}];
        m.resize(K.size2(), K.size1(), false);
        noalias(m) = transposeK;
      }
    }
    MoFEMFunctionReturn(0);
  };
 
protected:
  static inline std::map<std::pair<std::string, std::string>, MatrixDouble>
      mapMatrix;
};
 
struct OpAssembleFace : public OpAssembleBasic<FaceUserDataOperator> {
  using OpAssembleBasic<FaceUserDataOperator>::OpAssembleBasic;
};
 
struct OpAssembleVolumePositiveDefine : OpAssembleVolume {
  using OpAssembleVolume::OpAssembleVolume;
  MoFEMErrorCode doWork(int, EntityType, EntData &) {
    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);
 
    constexpr bool run_psd = true;
 
    if (run_psd) {
 
      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(trans(mat), 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;
 
        for (auto &e : eig) {
          constexpr double eps = std::numeric_limits<double>::epsilon();
          e = std::max(std::sqrt(std::abs(e)), eps);
        }
 
        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);
      K.swap(hat_k);
    }
 
    CHKERR assemble(row_side, col_side, row_type, col_type, row_data, col_data);
    MoFEMFunctionReturn(0);
  }
 
private:
};
 
/** @brief Calculate q = 3 K_ext epsilon_ext at integration points. */
struct OpCalculateExternalPressure : public VolUserDataOperator {
  OpCalculateExternalPressure(
      VectorPtr external_pressure_ptr,
      boost::shared_ptr<ExternalStrainVec> external_strain_vec_ptr,
      std::map<std::string, boost::shared_ptr<ScalingMethod>> smv);
 
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) override;
 
private:
  VectorPtr externalPressurePtr;
  boost::shared_ptr<ExternalStrainVec> externalStrainVecPtr;
  std::map<std::string, boost::shared_ptr<ScalingMethod>> scalingMethodsMap;
};
 
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 OpCalculateTractionFromSideEle : public SideEleOp {
  OpCalculateTractionFromSideEle(
      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),
        blockEntities(block_entities), blockName(block_name),
        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_r, double alpha_omega, double alpha_viscous_r,
                    double alpha_viscous_omega)
      : OpAssembleVolume(field_name, data_ptr, OPROW), alphaR(alpha_r),
        alphaOmega(alpha_omega), alphaViscousR(alpha_viscous_r),
        alphaViscousOmega(alpha_viscous_omega) {}
  MoFEMErrorCode integrate(EntData &data);
 
private:
  double alphaR;
  double alphaOmega;
  double alphaViscousR;
  double alphaViscousOmega;
};
 
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,
      boost::shared_ptr<ScalingMethod> scaling_method_ptr)
      : OpAssembleVolume(field_name, data_ptr, OPROW),
        scalingMethodPtr(scaling_method_ptr) {}
  MoFEMErrorCode integrate(EntData &row_data);
 
protected:
  boost::shared_ptr<ScalingMethod> scalingMethodPtr;
};
 
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(piola_scale_ptr),
        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,
      boost::shared_ptr<MatrixDouble> hybrid_grad_disp,
      std::map<std::string, boost::shared_ptr<ScalingMethod>> smv)
      : FormsIntegrators<FaceUserDataOperator>::Assembly<PETSC>::OpBase(
            row_field, row_field, FaceUserDataOperator::OPROW),
        bcData(bc_data), piolaScalePtr(piola_scale_ptr),
        hybridGradDispPtr(hybrid_grad_disp), scalingMethodsMap(smv) {}
 
  MoFEMErrorCode iNtegrate(EntData &data);
 
protected:
  boost::shared_ptr<PressureBcVec> bcData;
  boost::shared_ptr<double> piolaScalePtr;
  boost::shared_ptr<MatrixDouble> hybridGradDispPtr;
  std::map<std::string, boost::shared_ptr<ScalingMethod>> scalingMethodsMap;
};
 
template <AssemblyType A, IntegrationType I>
struct OpBrokenPressureBcLhsImpl_dU;
 
template <AssemblyType A>
struct OpBrokenPressureBcLhsImpl_dU<A, GAUSS>
    : public FormsIntegrators<FaceUserDataOperator>::Assembly<A>::OpBase {
 
  using OP =
      typename FormsIntegrators<FaceUserDataOperator>::Assembly<A>::OpBase;
 
  OpBrokenPressureBcLhsImpl_dU(
      std::string row_field, boost::shared_ptr<PressureBcVec> bc_data,
      boost::shared_ptr<MatrixDouble> hybrid_grad_disp,
      std::map<std::string, boost::shared_ptr<ScalingMethod>> smv)
      : OP(row_field, row_field, OP::OPROWCOL), bcData(bc_data),
        hybridGradDispPtr(hybrid_grad_disp), scalingMethodsMap(smv) {}
 
protected:
  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
  boost::shared_ptr<PressureBcVec> bcData;
  boost::shared_ptr<MatrixDouble> hybridGradDispPtr;
  std::map<std::string, boost::shared_ptr<ScalingMethod>> scalingMethodsMap;
};
 
struct OpBrokenPressureBcLhs_dU
    : public OpBrokenPressureBcLhsImpl_dU<A, GAUSS> {
  using OP = OpBrokenPressureBcLhsImpl_dU<A, GAUSS>;
  using OP::OpBrokenPressureBcLhsImpl_dU;
 
protected:
  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
};
 
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(piola_scale_ptr),
        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;
};
 
struct OpSpringRhsBc
    : public FormsIntegrators<FaceUserDataOperator>::Assembly<PETSC>::OpBase {
 
  using OP =
      typename FormsIntegrators<FaceUserDataOperator>::Assembly<PETSC>::OpBase;
 
  OpSpringRhsBc(
      std::string row_field, boost::shared_ptr<MatrixDouble> hybrid_disp_ptr,
      boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_side_data_ptr,
      boost::shared_ptr<SpringBcVec> bc_spring_ptr)
      : OP(row_field, row_field, OP::OPROW), hybridDispPtr(hybrid_disp_ptr),
        brokenBaseSideDataPtr(broken_side_data_ptr),
        bcSpringPtr(bc_spring_ptr) {}
 
  MoFEMErrorCode iNtegrate(EntData &data);
 
protected:
  boost::shared_ptr<MatrixDouble> hybridDispPtr;
  boost::shared_ptr<std::vector<BrokenBaseSideData>> brokenBaseSideDataPtr;
  boost::shared_ptr<SpringBcVec> bcSpringPtr;
};
 
struct OpSpringLhsBc_dU
    : public FormsIntegrators<FaceUserDataOperator>::Assembly<A>::OpBase {
 
  using OP =
      typename FormsIntegrators<FaceUserDataOperator>::Assembly<A>::OpBase;
 
  OpSpringLhsBc_dU(std::string row_field,
                   boost::shared_ptr<SpringBcVec> bc_spring_ptr)
      : OP(row_field, row_field, OP::OPROWCOL), bcSpringPtr(bc_spring_ptr) {}
 
  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
 
protected:
  boost::shared_ptr<SpringBcVec> bcSpringPtr;
};
 
struct OpSpringLhsBc_dP
    : public FormsIntegrators<FaceUserDataOperator>::Assembly<A>::OpBrokenBase {
 
  using OP = typename FormsIntegrators<FaceUserDataOperator>::Assembly<
      A>::OpBrokenBase;
 
  OpSpringLhsBc_dP(
      std::string row_field,
      boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
      boost::shared_ptr<SpringBcVec> bc_spring_ptr)
      : OP(row_field, broken_base_side_data, false, false, nullptr),
        bcSpringPtr(bc_spring_ptr) {}
 
  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
 
protected:
  boost::shared_ptr<SpringBcVec> bcSpringPtr;
};
 
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<std::vector<BrokenBaseSideData>> broken_side_data_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),
        brokenBaseSideDataPtr(
            broken_side_data_ptr), /*piolaStressPtr(piola_stress_ptr),*/
        bcDispPtr(bc_disp_ptr), scalingMethodsMap(smv) {}
 
protected:
  MoFEMErrorCode iNtegrate(EntData &data);
  boost::shared_ptr<MatrixDouble> hybridDispPtr;
  boost::shared_ptr<std::vector<BrokenBaseSideData>> brokenBaseSideDataPtr;
  // 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,
                              [[maybe_unused]] 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_r,
      double alpha_omega, double alpha_viscous_r, double alpha_viscous_omega)
      : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL, false),
        alphaR(alpha_r), alphaOmega(alpha_omega),
        alphaViscousR(alpha_viscous_r), alphaViscousOmega(alpha_viscous_omega) {
    sYmm = false;
  }
  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
 
private:
  double alphaR;
  double alphaOmega;
  double alphaViscousR;
  double alphaViscousOmega;
};
 
/// Add an intermediate-configuration rotation-gradient block to Schur PC
/// storage.
struct OpSpatialRotationPrec
    : public FormsIntegrators<VolUserDataOperator>::Assembly<PETSC>::OpBase {
  using OP = FormsIntegrators<VolUserDataOperator>::Assembly<PETSC>::OpBase;
 
  OpSpatialRotationPrec(const std::string &field_name,
                        boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                        const double alpha_omega_prec)
      : OP(field_name, field_name, OPROWCOL), dataAtPts(data_ptr),
        alphaOmegaPrec(alpha_omega_prec) {
    sYmm = false;
  }
 
protected:
  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data) override;
  MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data,
                          const bool) override;
 
private:
  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
  const double alphaOmegaPrec;
};
 
struct OpSpatialConsistency_dP_dP : public OpAssembleVolume {
  OpSpatialConsistency_dP_dP(std::string row_field, std::string col_field,
                             boost::shared_ptr<DataAtIntegrationPts> data_ptr,
                             const bool has_nonhomogeneous_mat_block)
      : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL, false),
        hasNonhomogeneousMatBlock(has_nonhomogeneous_mat_block) {
    sYmm = false;
  }
  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
 
private:
  template <int S>
  MoFEMErrorCode integrateImpl(EntData &row_data, EntData &col_data);
  const bool hasNonhomogeneousMatBlock;
};
 
struct OpSpatialConsistency_dBubble_dBubble : public OpAssembleVolume {
  OpSpatialConsistency_dBubble_dBubble(
      std::string row_field, std::string col_field,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      const bool has_nonhomogeneous_mat_block)
      : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL, false),
        hasNonhomogeneousMatBlock(has_nonhomogeneous_mat_block) {
    sYmm = false;
  }
  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
 
private:
  template <int S>
  MoFEMErrorCode integrateImpl(EntData &row_data, EntData &col_data);
  const bool hasNonhomogeneousMatBlock;
};
 
struct OpSpatialConsistency_dBubble_dP : public OpAssembleVolume {
  OpSpatialConsistency_dBubble_dP(
      std::string row_field, std::string col_field,
      boost::shared_ptr<DataAtIntegrationPts> data_ptr,
      const bool has_nonhomogeneous_mat_block)
      : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL, true),
        hasNonhomogeneousMatBlock(has_nonhomogeneous_mat_block) {
    sYmm = false;
  }
  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
 
private:
  template <int S>
  MoFEMErrorCode integrateImpl(EntData &row_data, EntData &col_data);
  const bool hasNonhomogeneousMatBlock;
};
 
using OpStressGram_dP_dP = FormsIntegrators<VolUserDataOperator>::Assembly<
    A>::BiLinearForm<GAUSS>::OpMass<3, 9>;
 
using OpStressGram_dBubble_dBubble = FormsIntegrators<
    VolUserDataOperator>::Assembly<A>::BiLinearForm<GAUSS>::OpMass<9, 9>;
 
struct OpStressGram_dBubble_dP : public OpAssembleVolume {
  OpStressGram_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);
};
 
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 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,
      boost::shared_ptr<MatrixDouble> jac = nullptr,
      boost::shared_ptr<VectorDouble> det = nullptr,
      boost::shared_ptr<MatrixDouble> inv_jac = nullptr);
};
 
using OpBaseTimesVectorFace =
    FormsIntegrators<FaceElementForcesAndSourcesCore::UserDataOperator>::
        Assembly<A>::LinearForm<GAUSS>::OpBaseTimesVector<1, SPACE_DIM,
                                                          SPACE_DIM>;
 
template <typename OP>
struct OpStabBrokenBaseImpl : public OpBrokenBaseImpl<OP> {
 
  using BASE = OpBrokenBaseImpl<OP>;
  using BASE::BASE;
 
  MoFEMErrorCode doWork(int row_side, EntityType row_type,
                        EntitiesFieldData::EntData &row_data) {
    MoFEMFunctionBegin;
 
    if (OP::entsPtr) {
      if (OP::entsPtr->find(this->getFEEntityHandle()) == OP::entsPtr->end())
        MoFEMFunctionReturnHot(0);
    }
 
#ifndef NDEBUG
    if (!BASE::brokenBaseSideData) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "space not set");
    }
#endif // NDEBUG
 
    auto do_work_rhs = [this](int, EntityType,
                              EntitiesFieldData::EntData &row_data) {
      MoFEMFunctionBegin;
      // get number of dofs on row
      OP::nbRows = row_data.getIndices().size();
      if (!OP::nbRows)
        MoFEMFunctionReturnHot(0);
      // get number of integration points
      OP::nbIntegrationPts = OP::getGaussPts().size2();
      // get row base functions
      OP::nbRowBaseFunctions = OP::getNbOfBaseFunctions(row_data);
      // resize and clear the right hand side vector
      OP::locF.resize(OP::nbRows, false);
      OP::locF.clear();
      // integrate local vector
      CHKERR this->iNtegrate(row_data);
      // assemble local vector
      CHKERR this->aSsemble(row_data);
      MoFEMFunctionReturn(0);
    };
 
    auto do_work_lhs = [this](int row_side, int col_side, EntityType row_type,
                              EntityType col_type,
                              EntitiesFieldData::EntData &row_data,
                              EntitiesFieldData::EntData &col_data) {
      MoFEMFunctionBegin;
 
      auto check_if_assemble_transpose = [&] {
        if (this->sYmm) {
          if (OP::rowSide != OP::colSide || OP::rowType != OP::colType)
            return true;
          else
            return false;
        } else if (OP::assembleTranspose) {
          return true;
        }
        return false;
      };
 
      OP::rowSide = row_side;
      OP::rowType = row_type;
      OP::colSide = col_side;
      OP::colType = col_type;
      OP::nbCols = col_data.getIndices().size();
      OP::locMat.resize(OP::nbRows, OP::nbCols, false);
      OP::locMat.clear();
      CHKERR this->iNtegrate(row_data, col_data);
      CHKERR this->aSsemble(row_data, col_data, check_if_assemble_transpose());
      MoFEMFunctionReturn(0);
    };
 
    switch (OP::opType) {
    case OP::OPROW:
 
      OP::nbRows = row_data.getIndices().size();
      if (!OP::nbRows)
        MoFEMFunctionReturnHot(0);
      OP::nbIntegrationPts = OP::getGaussPts().size2();
      OP::nbRowBaseFunctions = OP::getNbOfBaseFunctions(row_data);
 
      if (!OP::nbRows)
        MoFEMFunctionReturnHot(0);
 
      for (auto &bd : *BASE::brokenBaseSideData) {
 
#ifndef NDEBUG
        if (!bd.getData().getNSharedPtr(bd.getData().getBase())) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "base functions not set");
        }
#endif
 
        CHKERR do_work_lhs(
 
            // side
            row_side, bd.getSide(),
 
            // type
            row_type, bd.getType(),
 
            // row_data
            row_data, bd.getData()
 
        );
      }
 
      break;
    case OP::OPSPACE:
      for (auto &bd : *BASE::brokenBaseSideData) {
        fluxMatPtr = boost::shared_ptr<MatrixDouble>(BASE::brokenBaseSideData,
                                                     &bd.getFlux());
        CHKERR do_work_rhs(bd.getSide(), bd.getType(), bd.getData());
      }
      break;
    default:
      CHK_MOAB_THROW(
          MOFEM_IMPOSSIBLE_CASE,
          (std::string("wrong op type ") +
           OpBaseDerivativesBase::OpTypeNames[static_cast<unsigned char>(
               OP::opType)])
              .c_str());
    }
 
    MoFEMFunctionReturn(0);
  }
 
protected:
  boost::weak_ptr<MatrixDouble> fluxMatPtr;
};
 
struct OpBrokenBaseTimesHybridDisp
    : public OpStabBrokenBaseImpl<OpBaseTimesVectorFace> {
 
  using OP = OpStabBrokenBaseImpl<OpBaseTimesVectorFace>;
  OpBrokenBaseTimesHybridDisp(
      boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
      boost::shared_ptr<MatrixDouble> vec,
      ScalarFun beta_coeff = [](double, double, double) constexpr { return 1; },
      boost::shared_ptr<Range> ents_ptr = nullptr);
};
 
struct OpBrokenBaseTimesBrokenDisp
    : public OpStabBrokenBaseImpl<OpBaseTimesVectorFace> {
 
  using OP = OpStabBrokenBaseImpl<OpBaseTimesVectorFace>;
  OpBrokenBaseTimesBrokenDisp(
      boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
      ScalarFun beta_coeff = [](double, double, double) constexpr { return 1; },
      boost::shared_ptr<Range> ents_ptr = nullptr);
 
protected:
  MoFEMErrorCode doWork(int row_side, EntityType row_type,
                        EntitiesFieldData::EntData &row_data);
};
 
struct OpHybridBaseTimesBrokenDisp : public OpBaseTimesVectorFace {
 
  using OP = OpBaseTimesVectorFace;
  OpHybridBaseTimesBrokenDisp(
      const std::string row_field,
      boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
      ScalarFun beta_coeff, boost::shared_ptr<Range> ents_ptr = nullptr);
 
protected:
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
  boost::shared_ptr<std::vector<BrokenBaseSideData>> brokenBaseSideDataPtr;
};
 
using OpMassVectorFace =
    FormsIntegrators<FaceElementForcesAndSourcesCore::UserDataOperator>::
        Assembly<A>::BiLinearForm<GAUSS>::OpMass<1, SPACE_DIM>;
 
struct OpHyrbridBaseBrokenBase : public OpStabBrokenBaseImpl<OpMassVectorFace> {
 
  using OP = OpStabBrokenBaseImpl<OpMassVectorFace>;
  OpHyrbridBaseBrokenBase(
      std::string row_field,
      boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
      ScalarFun beta, const bool assmb_transpose, const bool only_transpose,
      boost::shared_ptr<Range> ents_ptr = nullptr);
};
 
struct OpBrokenBaseBrokenBase : public OpStabBrokenBaseImpl<OpMassVectorFace> {
 
  using OP = OpStabBrokenBaseImpl<OpMassVectorFace>;
  OpBrokenBaseBrokenBase(
      boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
      ScalarFun beta, boost::shared_ptr<Range> ents_ptr = nullptr);
 
  MoFEMErrorCode doWork(int row_side, EntityType row_type,
                        EntitiesFieldData::EntData &row_data);
};
 
struct OpTauStabilizationDispRhsBc : public OpBrokenBaseTimesBrokenDisp {
 
  using OP = OpBrokenBaseTimesBrokenDisp;
 
  OpTauStabilizationDispRhsBc(
      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,
      ScalarFun tau_coeff, boost::shared_ptr<Range> ents_ptr = nullptr)
      : OpBrokenBaseTimesBrokenDisp(broken_base_side_data, tau_coeff, 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 OpTauStabilizationDispLhsBc : public OpBrokenBaseBrokenBase {
 
  using OP = OpBrokenBaseBrokenBase;
 
  OpTauStabilizationDispLhsBc(
      boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
      boost::shared_ptr<BcDispVec> &bc_disp_ptr, ScalarFun tau_coeff,
      boost::shared_ptr<Range> ents_ptr = nullptr)
      : OP(broken_base_side_data, tau_coeff, ents_ptr), bcDispPtr(bc_disp_ptr) {
  }
 
protected:
  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
  boost::shared_ptr<BcDispVec> bcDispPtr;
};
 
struct OpTauStabilizationOpAnalyticalDispBc
    : public OpBrokenBaseTimesBrokenDisp {
 
  using OP = OpBrokenBaseTimesBrokenDisp;
 
  OpTauStabilizationOpAnalyticalDispBc(
      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,
      ScalarFun tau_coeff, boost::shared_ptr<Range> ents_ptr = nullptr)
      : OpBrokenBaseTimesBrokenDisp(broken_base_side_data, tau_coeff, 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 OpTauStabilizationOpAnalyticalDispLhsBc : public OpBrokenBaseBrokenBase {
 
  using OP = OpBrokenBaseBrokenBase;
 
  OpTauStabilizationOpAnalyticalDispLhsBc(
      boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
      boost::shared_ptr<AnalyticalDisplacementBcVec> &bc_disp_ptr,
      ScalarFun tau_coeff, boost::shared_ptr<Range> ents_ptr = nullptr)
      : OP(broken_base_side_data, tau_coeff, ents_ptr), bcDispPtr(bc_disp_ptr) {
  }
 
protected:
  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
  boost::shared_ptr<AnalyticalDisplacementBcVec> bcDispPtr;
};
 
struct OpTauStabilizationRotationRhsBc : public OpBrokenBaseTimesBrokenDisp {
 
  using OP = OpBrokenBaseTimesBrokenDisp;
 
  OpTauStabilizationRotationRhsBc(
      boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
      boost::shared_ptr<BcRotVec> &bc_ptr,
      std::map<std::string, boost::shared_ptr<ScalingMethod>> smv,
      ScalarFun tau_coeff, boost::shared_ptr<Range> ents_ptr = nullptr)
      : OpBrokenBaseTimesBrokenDisp(broken_base_side_data, tau_coeff, ents_ptr),
        bcRotPtr(bc_ptr), scalingMethodsMap(smv) {}
 
protected:
  MoFEMErrorCode iNtegrate(EntData &data);
  boost::shared_ptr<BcRotVec> bcRotPtr;
  std::map<std::string, boost::shared_ptr<ScalingMethod>> scalingMethodsMap;
};
 
struct OpTauStabilizationRotationLhsBc : public OpBrokenBaseBrokenBase {
 
  using OP = OpBrokenBaseBrokenBase;
 
  OpTauStabilizationRotationLhsBc(
      boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
      boost::shared_ptr<BcRotVec> &bc_ptr, ScalarFun tau_coeff,
      boost::shared_ptr<Range> ents_ptr = nullptr)
      : OP(broken_base_side_data, tau_coeff, ents_ptr), bcRotPtr(bc_ptr) {}
 
protected:
  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
  boost::shared_ptr<BcRotVec> bcRotPtr;
};