CohesiveInterfaceElement.hpp

Go to the documentation of this file.

/** \file CohesiveInterfaceElement.hpp
  \brief Implementation of linear interface element
 
*/
 
namespace CohesiveElement {
 
/** \brief Cohesive element implementation
 
  \bug Interface element not working with HO geometry.
*/
struct CohesiveInterfaceElement {
 
  struct CommonData {
    MatrixDouble gapGlob;
    MatrixDouble gapLoc;
    ublas::vector<MatrixDouble> R;
  };
  CommonData commonData;
 
  struct MyPrism : public MoFEM::FlatPrismElementForcesAndSourcesCore {
    MyPrism(MoFEM::Interface &m_field)
        : MoFEM::FlatPrismElementForcesAndSourcesCore(m_field) {
      meshPositionsFieldName = "GEOMETRY";
    }
    int getRule(int order) { return 2 * order; };
  };
  MyPrism feRhs;
  MyPrism feLhs;
  MyPrism feHistory;
 
  CohesiveInterfaceElement(MoFEM::Interface &m_field)
      : feRhs(m_field), feLhs(m_field), feHistory(m_field){};
 
  virtual ~CohesiveInterfaceElement() {}
 
  MyPrism &getFeRhs() { return feRhs; }
  MyPrism &getFeLhs() { return feLhs; }
  MyPrism &getFeHistory() { return feHistory; }
 
  /** \brief Constitutive (physical) equation for interface
 
    This is linear degradation model. Material parameters are: strength
    \f$f_t\f$, interface fracture energy \f$G_f\f$, elastic material stiffness
    \f$E\f$. Parameter \f$\beta\f$ controls how interface opening is calculated.
 
    Model parameter is interface penalty thickness \f$h\f$.
 
    */
  struct PhysicalEquation {
 
    MoFEM::Interface &mField;
    bool isInitialised;
    PhysicalEquation(MoFEM::Interface &m_field)
        : mField(m_field), isInitialised(false){};
 
    virtual ~PhysicalEquation() {}
 
    double h, youngModulus, beta, ft, Gf;
    Range pRisms;
    Tag thKappa, thDamagedPrism;
 
    double E0, g0, kappa1;
 
    /** \brief Initialize history variable data
 
    Create tag on the prism/interface to store damage history variable
 
    */
    MoFEMErrorCode iNitailise(const FEMethod *fe_method) {
      MoFEMFunctionBeginHot;
 
      double def_damaged = 0;
      CHKERR mField.get_moab().tag_get_handle(
          "DAMAGED_PRISM", 1, MB_TYPE_INTEGER, thDamagedPrism,
          MB_TAG_CREAT | MB_TAG_SPARSE, &def_damaged);
      const int def_len = 0;
      CHKERR mField.get_moab().tag_get_handle(
          "_KAPPA", def_len, MB_TYPE_DOUBLE, thKappa,
          MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN, nullptr);
      E0 = youngModulus / h;
      g0 = ft / E0;
      kappa1 = 2 * Gf / ft;
      MoFEMFunctionReturnHot(0);
    }
 
    /** \brief Calculate gap opening
 
    \f[
    g = \sqrt{ g_n^2 + \beta(g_{s1}^2 + g_{s2}^2)}
    \f]
 
    */
    double calcG(int gg, MatrixDouble gap_loc) {
      return sqrt(pow(gap_loc(gg, 0), 2) +
                  beta * (pow(gap_loc(gg, 1), 2) + pow(gap_loc(gg, 2), 2)));
    }
 
    double *kappaPtr;
    int kappaSize;
 
    /** \brief Get pointer from the mesh to histoy variables \f$\kappa\f$
    */
    MoFEMErrorCode getKappa(int nb_gauss_pts, const FEMethod *fe_method) {
      MoFEMFunctionBeginHot;
      EntityHandle ent = fe_method->numeredEntFiniteElementPtr->getEnt();
 
      rval = mField.get_moab().tag_get_by_ptr(
          thKappa, &ent, 1, (const void **)&kappaPtr, &kappaSize);
      if (rval != MB_SUCCESS || kappaSize != nb_gauss_pts) {
        VectorDouble kappa;
        kappa.resize(nb_gauss_pts);
        kappa.clear();
        int tag_size[1];
        tag_size[0] = nb_gauss_pts;
        void const *tag_data[] = {&kappa[0]};
        CHKERR mField.get_moab().tag_set_by_ptr(thKappa, &ent, 1, tag_data,
                                                tag_size);
        CHKERR mField.get_moab().tag_get_by_ptr(
            thKappa, &ent, 1, (const void **)&kappaPtr, &kappaSize);
      }
      MoFEMFunctionReturnHot(0);
    }
 
    MatrixDouble Dglob, Dloc;
 
    /** \brief Calculate stiffness material matrix
 
    \f[
    \mathbf{D}_\textrm{loc} = (1-\Omega) \mathbf{I} E_0
    \f]
    where \f$E_0\f$ is initial interface penalty stiffness
 
    \f[
    \mathbf{D}_\textrm{glob} = \mathbf{R}^\textrm{T} \mathbf{D}_\textrm{loc}\mathbf{R}
    \f]
 
    */
    MoFEMErrorCode calcDglob(const double omega, MatrixDouble &R) {
      MoFEMFunctionBeginHot;
      Dglob.resize(3, 3);
      Dloc.resize(3, 3);
      Dloc.clear();
      double E = (1 - omega) * E0;
      Dloc(0, 0) = E;
      Dloc(1, 1) = E;
      Dloc(2, 2) = E;
      Dglob = prod(Dloc, R);
      Dglob = prod(trans(R), Dglob);
      MoFEMFunctionReturnHot(0);
    }
 
    /** \brief Calculate damage
 
    \f[
    \Omega = \frac{1}{2} \frac{(2 G_f E_0+f_t^2)\kappa}{(ft+E_0 \kappa)G_f}
    \f]
 
    */
    MoFEMErrorCode calcOmega(const double kappa, double &omega) {
      MoFEMFunctionBeginHot;
      omega = 0;
      if (kappa >= kappa1) {
        omega = 1;
        MoFEMFunctionReturnHot(0);
      } else if (kappa > 0) {
        double a = (2.0 * Gf * E0 + ft * ft) * kappa;
        double b = (ft + E0 * kappa) * Gf;
        omega = 0.5 * a / b;
      }
      MoFEMFunctionReturnHot(0);
    }
 
    /** \brief Calculate tangent material stiffness
    */
    MoFEMErrorCode calcTangetDglob(const double omega, double g,
                                   const VectorDouble &gap_loc,
                                   MatrixDouble &R) {
      MoFEMFunctionBeginHot;
      Dglob.resize(3, 3);
      Dloc.resize(3, 3);
      double domega =
          0.5 * (2 * Gf * E0 + ft * ft) / ((ft + (g - ft / E0) * E0) * Gf) -
          0.5 * ((g - ft / E0) * (2 * Gf * E0 + ft * ft) * E0) /
              (pow(ft + (g - ft / E0) * E0, 2) * Gf);
      Dloc.resize(3, 3);
      //r0
      Dloc(0, 0) = (1 - omega) * E0 - domega * E0 * gap_loc[0] * gap_loc[0] / g;
      Dloc(0, 1) = -domega * E0 * gap_loc[0] * beta * gap_loc[1] / g;
      Dloc(0, 2) = -domega * E0 * gap_loc[0] * beta * gap_loc[2] / g;
      //r1
      Dloc(1, 0) = -domega * E0 * gap_loc[1] * gap_loc[0] / g;
      Dloc(1, 1) =
          (1 - omega) * E0 - domega * E0 * gap_loc[1] * beta * gap_loc[1] / g;
      Dloc(1, 2) = -domega * E0 * gap_loc[1] * beta * gap_loc[2] / g;
      //r2
      Dloc(2, 0) = -domega * E0 * gap_loc[2] * gap_loc[0] / g;
      Dloc(2, 1) = -domega * E0 * gap_loc[2] * beta * gap_loc[1] / g;
      Dloc(2, 2) =
          (1 - omega) * E0 - domega * E0 * gap_loc[2] * beta * gap_loc[2] / g;
      Dglob = prod(Dloc, R);
      Dglob = prod(trans(R), Dglob);
      MoFEMFunctionReturnHot(0);
    }
 
    /** \brief Calculate tractions
 
    \f[
    \mathbf{t} = \mathbf{D}_\textrm{glob}\mathbf{g}
    \f]
 
    */
    virtual MoFEMErrorCode calculateTraction(VectorDouble &traction, int gg,
                                             CommonData &common_data,
                                             const FEMethod *fe_method) {
      MoFEMFunctionBeginHot;
 
      if (!isInitialised) {
        CHKERR iNitailise(fe_method);
        isInitialised = true;
      }
      if (gg == 0) {
        CHKERR getKappa(common_data.gapGlob.size1(), fe_method);
      }
      double g = calcG(gg, common_data.gapLoc);
      double kappa = fmax(g - g0, kappaPtr[gg]);
      double omega = 0;
      CHKERR calcOmega(kappa, omega);
      CHKERR calcDglob(omega, common_data.R[gg]);
      traction.resize(3);
      ublas::matrix_row<MatrixDouble> gap_glob(common_data.gapGlob, gg);
      noalias(traction) = prod(Dglob, gap_glob);
      MoFEMFunctionReturnHot(0);
    }
 
    /** \brief Calculate tangent stiffness
    */
    virtual MoFEMErrorCode
    calculateTangentStiffeness(MatrixDouble &tangent_matrix, int gg,
                               CommonData &common_data,
                               const FEMethod *fe_method) {
      MoFEMFunctionBeginHot;
 
      try {
        if (!isInitialised) {
          CHKERR iNitailise(fe_method);
          isInitialised = true;
        }
        if (gg == 0) {
          CHKERR getKappa(common_data.gapGlob.size1(), fe_method);
        }
        double g = calcG(gg, common_data.gapLoc);
        double kappa = fmax(g - g0, kappaPtr[gg]);
        double omega = 0;
        CHKERR calcOmega(kappa, omega);
        int iter;
        CHKERR SNESGetIterationNumber(fe_method->snes, &iter);
        if ((kappa <= kappaPtr[gg]) || (kappa >= kappa1) || (iter <= 1)) {
          CHKERR calcDglob(omega, common_data.R[gg]);
        } else {
          ublas::matrix_row<MatrixDouble> g_loc(common_data.gapLoc, gg);
          CHKERR calcTangetDglob(omega, g, g_loc, common_data.R[gg]);
        }
        tangent_matrix.resize(3, 3);
        noalias(tangent_matrix) = Dglob;
      } catch (const std::exception &ex) {
        std::ostringstream ss;
        ss << "throw in method: " << ex.what() << std::endl;
        SETERRQ(PETSC_COMM_SELF, MOFEM_STD_EXCEPTION_THROW, "%s",
                ss.str().c_str());
      }
      MoFEMFunctionReturnHot(0);
    }
 
    /** \brief Update history variables when converged
    */
    virtual MoFEMErrorCode updateHistory(CommonData &common_data,
                                         const FEMethod *fe_method) {
      MoFEMFunctionBeginHot;
 
      if (!isInitialised) {
        CHKERR iNitailise(fe_method);
        isInitialised = true;
      }
      CHKERR getKappa(common_data.gapGlob.size1(), fe_method);
      bool all_gauss_pts_damaged = true;
      for (unsigned int gg = 0; gg < common_data.gapGlob.size1(); gg++) {
        double omega = 0;
        double g = calcG(gg, common_data.gapLoc);
        double kappa = fmax(g - g0, kappaPtr[gg]);
        kappaPtr[gg] = kappa;
        CHKERR calcOmega(kappa, omega);
        //if(omega < 1.) {
        all_gauss_pts_damaged = false;
        //}
      }
      if (all_gauss_pts_damaged) {
        EntityHandle ent = fe_method->numeredEntFiniteElementPtr->getEnt();
        int set_prism_as_demaged = 1;
        CHKERR mField.get_moab().tag_set_data(thDamagedPrism, &ent, 1,
                                              &set_prism_as_demaged);
      }
      MoFEMFunctionReturnHot(0);
    }
  };
 
  /** \brief Set negative sign to shape functions on face 4
    */
  struct OpSetSignToShapeFunctions
      : public FlatPrismElementForcesAndSourcesCore::UserDataOperator {
 
    OpSetSignToShapeFunctions(const std::string field_name)
        : FlatPrismElementForcesAndSourcesCore::UserDataOperator(
              field_name, ForcesAndSourcesCore::UserDataOperator::OPROW) {}
 
    MoFEMErrorCode doWork(int side, EntityType type,
                          EntitiesFieldData::EntData &data) {
      MoFEMFunctionBeginHot;
      if (data.getN().size1() == 0)
        MoFEMFunctionReturnHot(0);
      if (data.getN().size2() == 0)
        MoFEMFunctionReturnHot(0);
      switch (type) {
      case MBVERTEX:
        for (unsigned int gg = 0; gg < data.getN().size1(); gg++) {
          for (int nn = 3; nn < 6; nn++) {
            data.getN()(gg, nn) *= -1;
          }
        }
        break;
      case MBEDGE:
        if (side < 3)
          MoFEMFunctionReturnHot(0);
        data.getN() *= -1;
        break;
      case MBTRI:
        if (side == 3)
          MoFEMFunctionReturnHot(0);
        data.getN() *= -1;
        break;
      default:
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Unexpected entity type when setting interface shape-function signs");
      }
      MoFEMFunctionReturnHot(0);
    }
  };
 
  /** \brief Operator calculate gap, normal vector and rotation matrix
  */
  struct OpCalculateGapGlobal
      : public FlatPrismElementForcesAndSourcesCore::UserDataOperator {
 
    CommonData &commonData;
    OpCalculateGapGlobal(const std::string field_name, CommonData &common_data)
        : FlatPrismElementForcesAndSourcesCore::UserDataOperator(
              field_name, ForcesAndSourcesCore::UserDataOperator::OPROW),
          commonData(common_data) {}
 
    MoFEMErrorCode doWork(int side, EntityType type,
                          EntitiesFieldData::EntData &data) {
      MoFEMFunctionBeginHot;
      try {
        int nb_dofs = data.getIndices().size();
        if (nb_dofs == 0)
          MoFEMFunctionReturnHot(0);
        int nb_gauss_pts = data.getN().size1();
        if (type == MBVERTEX) {
          commonData.R.resize(nb_gauss_pts);
          for (int gg = 0; gg < nb_gauss_pts; gg++) {
            commonData.R[gg].resize(3, 3);
            double nrm2_normal = 0;
            double nrm2_tangent1 = 0;
            double nrm2_tangent2 = 0;
            for (int dd = 0; dd < 3; dd++) {
              nrm2_normal += pow(getNormalsAtGaussPtsF3()(gg, dd), 2);
              nrm2_tangent1 += pow(getTangent1AtGaussPtF3()(gg, dd), 2);
              nrm2_tangent2 += pow(getTangent2AtGaussPtF3()(gg, dd), 2);
            }
            nrm2_normal = sqrt(nrm2_normal);
            nrm2_tangent1 = sqrt(nrm2_tangent1);
            nrm2_tangent2 = sqrt(nrm2_tangent2);
            for (int dd = 0; dd < 3; dd++) {
              commonData.R[gg](0, dd) =
                  getNormalsAtGaussPtsF3()(gg, dd) / nrm2_normal;
              commonData.R[gg](1, dd) =
                  getTangent1AtGaussPtF3()(gg, dd) / nrm2_tangent1;
              commonData.R[gg](2, dd) =
                  getTangent2AtGaussPtF3()(gg, dd) / nrm2_tangent2;
            }
          }
        }
        if (type == MBVERTEX) {
          commonData.gapGlob.resize(nb_gauss_pts, 3);
          commonData.gapGlob.clear();
        }
        for (int gg = 0; gg < nb_gauss_pts; gg++) {
          for (int dd = 0; dd < 3; dd++) {
            commonData.gapGlob(gg, dd) += cblas_ddot(
                nb_dofs / 3, &data.getN(gg)[0], 1, &data.getFieldData()[dd], 3);
          }
        }
      } catch (const std::exception &ex) {
        std::ostringstream ss;
        ss << "throw in method: " << ex.what() << std::endl;
        SETERRQ(PETSC_COMM_SELF, MOFEM_STD_EXCEPTION_THROW, "%s",
                ss.str().c_str());
      }
      MoFEMFunctionReturnHot(0);
    }
  };
 
  /** \brief Operator calculate gap in local coordinate system
  */
  struct OpCalculateGapLocal
      : public FlatPrismElementForcesAndSourcesCore::UserDataOperator {
 
    CommonData &commonData;
    OpCalculateGapLocal(const std::string field_name, CommonData &common_data)
        : FlatPrismElementForcesAndSourcesCore::UserDataOperator(
              field_name, ForcesAndSourcesCore::UserDataOperator::OPROW),
          commonData(common_data) {}
 
    MoFEMErrorCode doWork(int side, EntityType type,
                          EntitiesFieldData::EntData &data) {
      MoFEMFunctionBeginHot;
      try {
        if (type == MBVERTEX) {
          int nb_gauss_pts = data.getN().size1();
          commonData.gapLoc.resize(nb_gauss_pts, 3);
          for (int gg = 0; gg < nb_gauss_pts; gg++) {
            ublas::matrix_row<MatrixDouble> gap_glob(commonData.gapGlob, gg);
            ublas::matrix_row<MatrixDouble> gap_loc(commonData.gapLoc, gg);
            gap_loc = prod(commonData.R[gg], gap_glob);
          }
        }
      } catch (const std::exception &ex) {
        std::ostringstream ss;
        ss << "throw in method: " << ex.what() << std::endl;
        SETERRQ(PETSC_COMM_SELF, MOFEM_STD_EXCEPTION_THROW, "%s",
                ss.str().c_str());
      }
      MoFEMFunctionReturnHot(0);
    }
  };
 
  /** \brief Operator calculate right hand side vector
  */
  struct OpRhs : public FlatPrismElementForcesAndSourcesCore::UserDataOperator {
 
    CommonData &commonData;
    PhysicalEquation &physicalEqations;
    OpRhs(const std::string field_name, CommonData &common_data,
          PhysicalEquation &physical_eqations)
        : FlatPrismElementForcesAndSourcesCore::UserDataOperator(
              field_name, ForcesAndSourcesCore::UserDataOperator::OPROW),
          commonData(common_data), physicalEqations(physical_eqations) {}
 
    VectorDouble traction, Nf;
    MoFEMErrorCode doWork(int side, EntityType type,
                          EntitiesFieldData::EntData &data) {
      MoFEMFunctionBeginHot;
 
      try {
        int nb_dofs = data.getIndices().size();
        if (nb_dofs == 0)
          MoFEMFunctionReturnHot(0);
        if (physicalEqations.pRisms.find(
                getNumeredEntFiniteElementPtr()->getEnt()) ==
            physicalEqations.pRisms.end()) {
          MoFEMFunctionReturnHot(0);
        }
        Nf.resize(nb_dofs);
        Nf.clear();
        int nb_gauss_pts = data.getN().size1();
        for (int gg = 0; gg < nb_gauss_pts; gg++) {
          CHKERR physicalEqations.calculateTraction(traction, gg, commonData,
                                                    getFEMethod());
          double w = getGaussPts()(2, gg) *
                     cblas_dnrm2(3, &getNormalsAtGaussPtsF3()(gg, 0), 1) * 0.5;
          for (int nn = 0; nn < nb_dofs / 3; nn++) {
            for (int dd = 0; dd < 3; dd++) {
              Nf[3 * nn + dd] += w * data.getN(gg)[nn] * traction[dd];
            }
          }
        }
        CHKERR VecSetValues(getFEMethod()->snes_f, data.getIndices().size(),
                            &data.getIndices()[0], &Nf[0], ADD_VALUES);
      } catch (const std::exception &ex) {
        std::ostringstream ss;
        ss << "throw in method: " << ex.what() << std::endl;
        SETERRQ(PETSC_COMM_SELF, MOFEM_STD_EXCEPTION_THROW, "%s",
                ss.str().c_str());
      }
      MoFEMFunctionReturnHot(0);
    }
  };
 
  /** \brief Operator calculate element stiffens matrix
  */
  struct OpLhs : public FlatPrismElementForcesAndSourcesCore::UserDataOperator {
 
    CommonData &commonData;
    PhysicalEquation &physicalEqations;
    OpLhs(const std::string field_name, CommonData &common_data,
          PhysicalEquation &physical_eqations)
        : FlatPrismElementForcesAndSourcesCore::UserDataOperator(
              field_name, ForcesAndSourcesCore::UserDataOperator::OPROWCOL),
          commonData(common_data), physicalEqations(physical_eqations) {
      sYmm = false;
    }
 
    MatrixDouble K, D, ND;
    MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
                          EntityType col_type,
                          EntitiesFieldData::EntData &row_data,
                          EntitiesFieldData::EntData &col_data) {
      MoFEMFunctionBeginHot;
 
      try {
        int nb_row = row_data.getIndices().size();
        if (nb_row == 0)
          MoFEMFunctionReturnHot(0);
        int nb_col = col_data.getIndices().size();
        if (nb_col == 0)
          MoFEMFunctionReturnHot(0);
        if (physicalEqations.pRisms.find(
                getNumeredEntFiniteElementPtr()->getEnt()) ==
            physicalEqations.pRisms.end()) {
          MoFEMFunctionReturnHot(0);
        }
        ND.resize(nb_row, 3);
        K.resize(nb_row, nb_col);
        K.clear();
        int nb_gauss_pts = row_data.getN().size1();
        for (int gg = 0; gg < nb_gauss_pts; gg++) {
          CHKERR physicalEqations.calculateTangentStiffeness(D, gg, commonData,
                                                             getFEMethod());
          double w = getGaussPts()(2, gg) *
                     cblas_dnrm2(3, &getNormalsAtGaussPtsF3()(gg, 0), 1) * 0.5;
          ND.clear();
          for (int nn = 0; nn < nb_row / 3; nn++) {
            for (int dd = 0; dd < 3; dd++) {
              for (int DD = 0; DD < 3; DD++) {
                ND(3 * nn + dd, DD) += row_data.getN(gg)[nn] * D(dd, DD);
              }
            }
          }
          for (int nn = 0; nn < nb_row / 3; nn++) {
            for (int dd = 0; dd < 3; dd++) {
              for (int NN = 0; NN < nb_col / 3; NN++) {
                for (int DD = 0; DD < 3; DD++) {
                  K(3 * nn + dd, 3 * NN + DD) +=
                      w * ND(3 * nn + dd, DD) * col_data.getN(gg)[NN];
                }
              }
            }
          }
        }
        CHKERR MatSetValues(getFEMethod()->snes_B, nb_row,
                            &row_data.getIndices()[0], nb_col,
                            &col_data.getIndices()[0], &K(0, 0), ADD_VALUES);
      } catch (const std::exception &ex) {
        std::ostringstream ss;
        ss << "throw in method: " << ex.what() << std::endl;
        SETERRQ(PETSC_COMM_SELF, MOFEM_STD_EXCEPTION_THROW, "%s",
                ss.str().c_str());
      }
      MoFEMFunctionReturnHot(0);
    }
  };
 
  /** \brief Operator update history variables
  */
  struct OpHistory
      : public FlatPrismElementForcesAndSourcesCore::UserDataOperator {
 
    CommonData &commonData;
    PhysicalEquation &physicalEqations;
    OpHistory(const std::string field_name, CommonData &common_data,
              PhysicalEquation &physical_eqations)
        : FlatPrismElementForcesAndSourcesCore::UserDataOperator(
              field_name, ForcesAndSourcesCore::UserDataOperator::OPROW),
          commonData(common_data), physicalEqations(physical_eqations) {}
 
    MoFEMErrorCode doWork(int side, EntityType type,
                          EntitiesFieldData::EntData &data) {
      MoFEMFunctionBeginHot;
 
      if (type != MBVERTEX)
        MoFEMFunctionReturnHot(0);
      if (physicalEqations.pRisms.find(
              getNumeredEntFiniteElementPtr()->getEnt()) ==
          physicalEqations.pRisms.end()) {
        MoFEMFunctionReturnHot(0);
      }
      CHKERR physicalEqations.updateHistory(commonData, getFEMethod());
      MoFEMFunctionReturnHot(0);
    }
  };
 
  /** \brief Driver function settting all operators needed for interface element
  */
  MoFEMErrorCode
  addOps(const std::string field_name,
         boost::ptr_vector<CohesiveInterfaceElement::PhysicalEquation>
             &interfaces) {
    MoFEMFunctionBeginHot;
 
    //Rhs
    feRhs.getOpPtrVector().push_back(new OpSetSignToShapeFunctions(field_name));
    feRhs.getOpPtrVector().push_back(
        new OpCalculateGapGlobal(field_name, commonData));
    feRhs.getOpPtrVector().push_back(
        new OpCalculateGapLocal(field_name, commonData));
    //Lhs
    feLhs.getOpPtrVector().push_back(new OpSetSignToShapeFunctions(field_name));
    feLhs.getOpPtrVector().push_back(
        new OpCalculateGapGlobal(field_name, commonData));
    feLhs.getOpPtrVector().push_back(
        new OpCalculateGapLocal(field_name, commonData));
    //History
    feHistory.getOpPtrVector().push_back(
        new OpSetSignToShapeFunctions(field_name));
    feHistory.getOpPtrVector().push_back(
        new OpCalculateGapGlobal(field_name, commonData));
    feHistory.getOpPtrVector().push_back(
        new OpCalculateGapLocal(field_name, commonData));
 
    //add equations/data for physical interfaces
    boost::ptr_vector<CohesiveInterfaceElement::PhysicalEquation>::iterator pit;
    for (pit = interfaces.begin(); pit != interfaces.end(); pit++) {
      feRhs.getOpPtrVector().push_back(new OpRhs(field_name, commonData, *pit));
      feLhs.getOpPtrVector().push_back(new OpLhs(field_name, commonData, *pit));
      feHistory.getOpPtrVector().push_back(
          new OpHistory(field_name, commonData, *pit));
    }
 
    MoFEMFunctionReturnHot(0);
  }
};
 
} // namespace CohesiveElement