ContactOperators.hpp

Go to the documentation of this file.

/* This file is part of MoFEM.
 * MoFEM is free software: you can redistribute it and/or modify it under
 * the terms of the GNU Lesser General Public License as published by the
 * Free Software Foundation, either version 3 of the License, or (at your
 * option) any later version.
 *
 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 * License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
 
namespace OpContactTools {
 
 
//! [Common data]
struct CommonData : public OpPlasticTools::CommonData {
  boost::shared_ptr<MatrixDouble> contactStressPtr;
  boost::shared_ptr<MatrixDouble> contactStressDivergencePtr;
  boost::shared_ptr<MatrixDouble> contactTractionPtr;
  boost::shared_ptr<MatrixDouble> contactNormalPtr;
  boost::shared_ptr<VectorDouble> contactGapPtr;
  boost::shared_ptr<MatrixDouble> contactGapDiffPtr;
  boost::shared_ptr<MatrixDouble> contactNormalDiffPtr;
  // boost::shared_ptr<MatrixDouble> mDispPtr;
  Tag rollerTag;
};
//! [Common data]
 
struct OpPostProcVertex : public BoundaryEleOp {
  OpPostProcVertex(const std::string field_name,
                   boost::shared_ptr<CommonData> common_data_ptr,
                   moab::Interface *moab_inst);
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
  moab::Interface *moab;
};
 
struct OpConstrainBoundaryRhs : public BoundaryEleOpAssembly {
  OpConstrainBoundaryRhs(const std::string field_name,
                         boost::shared_ptr<CommonData> common_data_ptr);
  MoFEMErrorCode iNtegrate(EntData &data);
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
};
 
struct OpConstrainBoundaryLhs_dU : public BoundaryEleOpAssembly {
  OpConstrainBoundaryLhs_dU(const std::string row_field_name,
                            const std::string col_field_name,
                            boost::shared_ptr<CommonData> common_data_ptr);
  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
 
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
};
 
struct OpConstrainBoundaryLhs_dTraction : public BoundaryEleOpAssembly {
  OpConstrainBoundaryLhs_dTraction(
      const std::string row_field_name, const std::string col_field_name,
      boost::shared_ptr<CommonData> common_data_ptr);
  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
 
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
};
 
struct OpGetClosestRigidBodyData : public BoundaryEleOp {
  OpGetClosestRigidBodyData(MoFEM::Interface &m_field,
                            const std::string field_name,
                            boost::shared_ptr<CommonData> common_data_ptr,
                            bool update);
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
private:
  MoFEM::Interface &mField;
  boost::shared_ptr<CommonData> commonDataPtr;
  bool uPdate;
};
 
 
template <typename T>
inline double gap0(double g, FTensor::Tensor1<T, 3> &t_disp,
                   FTensor::Tensor1<T, 3> &t_normal) {
  return -g + t_disp(i) * t_normal(i);
}
 
template <typename T1, typename T2>
inline double normal_traction(Tensor1<T1, 3> &t_traction,
                              Tensor1<T2, 3> &t_normal) {
  return -t_traction(i) * t_normal(i);
}
 
inline double sign(double x) {
  if (x == 0)
    return 0;
  else if (x > 0)
    return 1;
  else
    return -1;
};
 
inline double w(const double g, const double t) {
  return g - (*cache).cn_cont * t;
}
 
inline double constraint(double g0, double g, double &&t) {
  return (w(g, t) + std::abs(w(g, t))) / 2 + g0;
};
 
inline double diff_constraints_dtraction(double g, double &&t) {
  return -(*cache).cn_cont * (1 + sign(w(g, t))) / 2;
}
 
inline double diff_constraints_dgap(double g, double &&t) {
  return (sign(w(g, t)) - 1) / 2;
}
 
OpConstrainBoundaryRhs::OpConstrainBoundaryRhs(
    const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
    : BoundaryEleOpAssembly(field_name, field_name, BoundaryEleOp::OPROW),
      commonDataPtr(common_data_ptr) {}
 
MoFEMErrorCode OpConstrainBoundaryRhs::iNtegrate(EntData &data) {
  MoFEMFunctionBegin;
 
  const size_t nb_gauss_pts = getGaussPts().size2();
  const size_t nb_dofs = data.getIndices().size();
 
  if (nb_dofs) {
 
    auto t_normal = getFTensor1Normal();
    const double ls = sqrt(t_normal(i) * t_normal(i));
    t_normal(i) /= ls;
 
    auto t_w = getFTensor0IntegrationWeight();
    auto t_disp = getFTensor1FromMat<3>(*(commonDataPtr->mDispPtr));
    auto t_traction =
        getFTensor1FromMat<3>(*(commonDataPtr->contactTractionPtr));
    auto t_coords = getFTensor1CoordsAtGaussPts();
 
    size_t nb_base_functions = data.getN().size2() / 3;
    auto t_base = data.getFTensor1N<3>();
 
    auto t_contact_normal =
        getFTensor1FromMat<3>(*(commonDataPtr->contactNormalPtr));
    auto t_gap = getFTensor0FromVec(*(commonDataPtr->contactGapPtr));
 
    for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
 
      auto t_nf = getFTensor1FromArray<3, 3>(locF);
 
      const double alpha = t_w * getMeasure();
      if (getNormalsAtGaussPts().size1() == nb_gauss_pts) {
        VectorDouble n = getNormalsAtGaussPts(gg);
        auto t_n = getFTensor1FromPtr<3>(&*n.data().begin());
        t_normal(i) = t_n(i) / sqrt(t_n(j) * t_n(j));
      }
      Tensor2<double, 3, 3> t_P;
      t_P(i, j) = t_contact_normal(i) * t_contact_normal(j);
 
      Tensor2<double, 3, 3> t_Q;
      t_Q(i, j) = kronecker_delta(i, j) - t_P(i, j);
 
      const double gap_0 = gap0(t_gap, t_disp, t_contact_normal);
      Tensor1<double, 3> t_rhs_constrains;
      t_rhs_constrains(i) =
          t_contact_normal(i) *
          constraint(gap_0, t_gap,
                     normal_traction(t_traction, t_contact_normal));
 
      Tensor1<double, 3> t_rhs_tangent_disp, t_rhs_tangent_traction;
      t_rhs_tangent_disp(i) = t_Q(i, j) * t_disp(j);
      t_rhs_tangent_traction(i) = (*cache).cn_cont * t_Q(i, j) * t_traction(j);
 
      size_t bb = 0;
      for (; bb != nb_dofs / 3; ++bb) {
        const double beta = alpha * (t_base(i) * t_normal(i));
 
        t_nf(i) -= beta * t_rhs_constrains(i);
        t_nf(i) -= beta * t_rhs_tangent_disp(i);
        t_nf(i) += beta * t_rhs_tangent_traction(i);
 
        ++t_nf;
        ++t_base;
      }
      for (; bb < nb_base_functions; ++bb)
        ++t_base;
 
      ++t_contact_normal;
      ++t_gap;
      ++t_disp;
      ++t_traction;
      ++t_coords;
      ++t_w;
    }
 
  }
 
  MoFEMFunctionReturn(0);
}
 
OpPostProcVertex::OpPostProcVertex(
    const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr,
    moab::Interface *moab_inst)
    : BoundaryEleOp(field_name, field_name, BoundaryEleOp::OPROW),
      commonDataPtr(common_data_ptr), moab(moab_inst) {
  std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
  doEntities[MBTRI] = doEntities[MBTRI] = doEntities[MBQUAD] = true;
}
 
MoFEMErrorCode OpPostProcVertex::doWork(int side, EntityType type,
                                        EntData &data) {
  MoFEMFunctionBegin;
 
  const size_t nb_gauss_pts = getGaussPts().size2();
  const size_t nb_dofs = data.getIndices().size();
 
  bool empty = std::all_of(data.getIndices().begin(), data.getIndices().end(),
                           [](int i) { return i == -1; });
  if (empty)
    MoFEMFunctionReturnHot(0);
 
  EntityHandle ent = getFEEntityHandle();
  auto get_tag = [&](const std::string name, size_t size) {
    Tag th;
    double def[3] = {0, 0, 0};
    CHKERR moab->tag_get_handle(name.c_str(), size, MB_TYPE_DOUBLE, th,
                                MB_TAG_CREAT | MB_TAG_SPARSE, def);
    return th;
  };
 
  auto th_gap = get_tag("GAP", 1);
  auto th_cons = get_tag("CONSTRAINT", 1);
  auto th_roller = get_tag("ROLLER_ID", 1);
  auto th_traction = get_tag("TRACTION", 3);
  auto th_normal = get_tag("NORMAL", 3);
  auto th_normal_tri = get_tag("NORMAL_TRI", 3);
  auto th_disp = get_tag("DISPLACEMENT", 3);
 
  auto t_contact_normal =
      getFTensor1FromMat<3>(*(commonDataPtr->contactNormalPtr));
  auto t_gap = getFTensor0FromVec(*(commonDataPtr->contactGapPtr));
  auto t_traction = getFTensor1FromMat<3>(*(commonDataPtr->contactTractionPtr));
  auto t_disp = getFTensor1FromMat<3>(*(commonDataPtr->mDispPtr));
  auto t_coords = getFTensor1CoordsAtGaussPts();
 
  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
 
    double coords[] = {0, 0, 0};
    EntityHandle new_vertex;
    for (int dd = 0; dd != 3; dd++) {
      coords[dd] = getCoordsAtGaussPts()(gg, dd);
    }
    CHKERR moab->create_vertex(&coords[0], new_vertex);
 
    Tensor1<double, 3> trac(t_traction(0), t_traction(1), t_traction(2));
    Tensor1<double, 3> disp(t_disp(0), t_disp(1), t_disp(2));
    Tensor1<double, 3> norm(t_contact_normal(0), t_contact_normal(1),
                            t_contact_normal(2));
    trac(i) /= (*cache).young_modulus_inv;
    double fin_gap = t_gap;
    const double gap_0 = gap0(t_gap, t_disp, t_contact_normal);
    double constra =
        constraint(gap_0, t_gap, normal_traction(t_traction, t_contact_normal));
    CHKERR moab->tag_set_data(th_gap, &new_vertex, 1, &fin_gap);
    CHKERR moab->tag_set_data(th_cons, &new_vertex, 1, &constra);
 
    CHKERR moab->tag_set_data(th_traction, &new_vertex, 1, &trac(0));
    CHKERR moab->tag_set_data(th_normal, &new_vertex, 1, &norm(0));
    CHKERR moab->tag_set_data(th_disp, &new_vertex, 1, &disp(0));
    auto el_normal = getNormal();
    if (getNormalsAtGaussPts().size1() == nb_gauss_pts) {
      VectorDouble n = getNormalsAtGaussPts(gg);
      el_normal = n;
    }
    CHKERR moab->tag_set_data(th_normal_tri, &new_vertex, 1, &el_normal(0));
 
    double d_tag_data = (*cache).closest_roller_id;
    CHKERR moab->tag_set_data(th_roller, &new_vertex, 1, &d_tag_data);
 
    ++t_contact_normal;
    ++t_gap;
    ++t_disp;
    ++t_coords;
    ++t_traction;
  }
 
  MoFEMFunctionReturn(0);
}
 
OpConstrainBoundaryLhs_dU::OpConstrainBoundaryLhs_dU(
    const std::string row_field_name, const std::string col_field_name,
    boost::shared_ptr<CommonData> common_data_ptr)
    : BoundaryEleOpAssembly(row_field_name, col_field_name, OPROWCOL),
      commonDataPtr(common_data_ptr) {
  sYmm = false;
}
MoFEMErrorCode OpConstrainBoundaryLhs_dU::iNtegrate(EntData &row_data,
                                                    EntData &col_data) {
  MoFEMFunctionBegin;
 
  const size_t nb_gauss_pts = getGaussPts().size2();
  const size_t row_nb_dofs = row_data.getIndices().size();
  const size_t col_nb_dofs = col_data.getIndices().size();
 
  if (row_nb_dofs && col_nb_dofs) {
 
 
    auto t_normal = getFTensor1Normal();
    const double ls = sqrt(t_normal(i) * t_normal(i));
    t_normal(i) /= ls;
 
    auto t_disp = getFTensor1FromMat<3>(*(commonDataPtr->mDispPtr));
    auto t_traction =
        getFTensor1FromMat<3>(*(commonDataPtr->contactTractionPtr));
    auto t_coords = getFTensor1CoordsAtGaussPts();
 
    auto t_w = getFTensor0IntegrationWeight();
    auto t_row_base = row_data.getFTensor1N<3>();
    size_t nb_face_functions = row_data.getN().size2() / 3;
 
    auto t_contact_normal =
        getFTensor1FromMat<3>(*(commonDataPtr->contactNormalPtr));
    auto t_gap = getFTensor0FromVec(*(commonDataPtr->contactGapPtr));
    auto t_gap_diff =
        getFTensor1FromMat<3>(*(commonDataPtr->contactGapDiffPtr));
    auto t_diff_normal =
        getFTensor2FromMat<3, 3>(*(commonDataPtr->contactNormalDiffPtr));
    Tensor1<double, 3> t1;
    auto &cn = (*cache).cn_cont;
 
    for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
 
      const double alpha = t_w * getMeasure();
 
      if (getNormalsAtGaussPts().size1() == nb_gauss_pts) {
        VectorDouble n = getNormalsAtGaussPts(gg);
        auto t_n = getFTensor1FromPtr<3>(&*n.data().begin());
        t_normal(i) = t_n(i) / sqrt(t_n(j) * t_n(j));
      }
 
      Tensor3<double, 3, 3, 3> t_diff_contact_tangent_tensor;
      t_diff_contact_tangent_tensor(i, j, k) =
          -t_diff_normal(i, k) * t_contact_normal(j) -
          t_contact_normal(i) * t_diff_normal(j, k);
 
      Tensor2<double, 3, 3> t_P;
      t_P(i, j) = t_contact_normal(i) * t_contact_normal(j);
 
      Tensor2<double, 3, 3> t_Q;
      t_Q(i, j) = kronecker_delta(i, j) - t_P(i, j);
      const double gap_0 = gap0(t_gap, t_disp, t_contact_normal);
 
      t1(i) = t_contact_normal(i) - 0.5 * t_gap_diff(i) +
              0.5 * cn * t_traction(j) * t_diff_normal(i, j) +
              t_disp(j) * t_diff_normal(i, j) +
              0.5 * sign(t_gap + cn * t_traction(k) * t_contact_normal(k)) *
                  (t_gap_diff(i) + cn * t_traction(j) * t_diff_normal(i, j));
 
      auto constrain = constraint(
          gap_0, t_gap, normal_traction(t_traction, t_contact_normal));
      size_t rr = 0;
      for (; rr != row_nb_dofs / 3; ++rr) {
 
        auto t_mat = getFTensor2FromArray<3, 3, 3>(locMat, 3 * rr);
        const double row_base = t_row_base(i) * t_normal(i);
 
        auto t_col_base = col_data.getFTensor0N(gg, 0);
 
        for (size_t cc = 0; cc != col_nb_dofs / 3; ++cc) {
          const double beta = alpha * row_base * t_col_base;
 
          // linearisation of t_rhs_constrains(i)
          t_mat(i, j) -= beta * t_contact_normal(i) * t1(j);
          t_mat(i, j) -= (beta * constrain) * t_diff_normal(i, j);
 
          // // linearisation of t_rhs_tangent_disp(i)
          t_mat(i, j) -= beta * t_Q(i, j);
          t_mat(i, k) -=
              beta * t_diff_contact_tangent_tensor(i, j, k) * t_disp(j);
 
          // // linearisation of  t_rhs_tangent_traction(i);
          t_mat(i, k) += beta * t_diff_contact_tangent_tensor(i, j, k) *
                         (*cache).cn_cont * t_traction(j);
 
          ++t_col_base;
          ++t_mat;
        }
        ++t_row_base;
      }
      for (; rr < nb_face_functions; ++rr)
        ++t_row_base;
 
      ++t_contact_normal;
      ++t_gap;
      ++t_gap_diff;
      ++t_diff_normal;
      ++t_disp;
      ++t_traction;
      ++t_coords;
      ++t_w;
    }
 
  }
 
  MoFEMFunctionReturn(0);
}
 
OpConstrainBoundaryLhs_dTraction::OpConstrainBoundaryLhs_dTraction(
    const std::string row_field_name, const std::string col_field_name,
    boost::shared_ptr<CommonData> common_data_ptr)
    : BoundaryEleOpAssembly(row_field_name, col_field_name, OPROWCOL),
      commonDataPtr(common_data_ptr) {
  sYmm = false;
}
MoFEMErrorCode OpConstrainBoundaryLhs_dTraction::iNtegrate(
    EntData &row_data, EntData &col_data) {
  MoFEMFunctionBegin;
 
  const size_t nb_gauss_pts = getGaussPts().size2();
  const size_t row_nb_dofs = row_data.getIndices().size();
  const size_t col_nb_dofs = col_data.getIndices().size();
 
  if (row_nb_dofs && col_nb_dofs) {
 
    auto t_normal = getFTensor1Normal();
    const double ls = sqrt(t_normal(i) * t_normal(i));
    t_normal(i) /= ls;
 
    auto t_disp = getFTensor1FromMat<3>(*(commonDataPtr->mDispPtr));
    auto t_traction =
        getFTensor1FromMat<3>(*(commonDataPtr->contactTractionPtr));
    auto t_coords = getFTensor1CoordsAtGaussPts();
 
    auto t_w = getFTensor0IntegrationWeight();
    auto t_row_base = row_data.getFTensor1N<3>();
    size_t nb_face_functions = row_data.getN().size2() / 3;
 
    auto t_contact_normal =
        getFTensor1FromMat<3>(*(commonDataPtr->contactNormalPtr));
    auto t_gap = getFTensor0FromVec(*(commonDataPtr->contactGapPtr));
 
    for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
 
      const double alpha = t_w * getMeasure();
 
      if (getNormalsAtGaussPts().size1() == nb_gauss_pts) {
        VectorDouble n = getNormalsAtGaussPts(gg);
        auto t_n = getFTensor1FromPtr<3>(&*n.data().begin());
        t_normal(i) = t_n(i) / sqrt(t_n(j) * t_n(j));
      }
 
      // auto t_contact_normal = normal(t_coords, t_disp);
      Tensor2<double, 3, 3> t_P;
      t_P(i, j) = t_contact_normal(i) * t_contact_normal(j);
 
      Tensor2<double, 3, 3> t_Q;
      t_Q(i, j) = kronecker_delta(i, j) - t_P(i, j);
 
      const double diff_traction = diff_constraints_dtraction(
          t_gap, normal_traction(t_traction, t_contact_normal));
 
      size_t rr = 0;
      for (; rr != row_nb_dofs / 3; ++rr) {
        auto t_mat = getFTensor2FromArray<3, 3, 3>(locMat, 3 * rr);
        const double row_base = t_row_base(i) * t_normal(i);
 
        auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
 
        for (size_t cc = 0; cc != col_nb_dofs / 3; ++cc) {
          const double col_base = t_col_base(i) * t_normal(i);
          const double beta = alpha * row_base * col_base;
 
          t_mat(i, j) += (beta * diff_traction) * t_P(i, j);
          t_mat(i, j) += beta * (*cache).cn_cont * t_Q(i, j);
 
          ++t_col_base;
          ++t_mat;
        }
 
        ++t_row_base;
      }
      for (; rr < nb_face_functions; ++rr)
        ++t_row_base;
 
      ++t_contact_normal;
      ++t_gap;
      ++t_disp;
      ++t_traction;
      ++t_coords;
      ++t_w;
    }
 
  }
 
  MoFEMFunctionReturn(0);
}
 
OpGetClosestRigidBodyData::OpGetClosestRigidBodyData(
    MoFEM::Interface &m_field, const std::string field_name,
    boost::shared_ptr<CommonData> common_data_ptr, bool update)
    : mField(m_field),
      BoundaryEleOp(field_name, field_name, DomainEleOp::OPROW),
      commonDataPtr(common_data_ptr), uPdate(update) {
  std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
  doEntities[MBTRI] = doEntities[MBQUAD] = true;
}
 
MoFEMErrorCode OpGetClosestRigidBodyData::doWork(int side, EntityType type,
                                                 EntData &data) {
  MoFEMFunctionBegin;
 
  const size_t nb_gauss_pts = getGaussPts().size2();
  const size_t nb_dofs = data.getIndices().size();
 
  EntityHandle ent = getFEEntityHandle();
  auto t_coords = getFTensor1CoordsAtGaussPts();
  auto t_disp = getFTensor1FromMat<3>(*(commonDataPtr->mDispPtr));
 
  commonDataPtr->contactNormalPtr->resize(3, nb_gauss_pts, false);
  commonDataPtr->contactNormalPtr->clear();
  commonDataPtr->contactNormalDiffPtr->resize(9, nb_gauss_pts, false);
  commonDataPtr->contactNormalDiffPtr->clear();
  commonDataPtr->contactGapDiffPtr->resize(3, nb_gauss_pts, false);
  commonDataPtr->contactGapDiffPtr->clear();
  commonDataPtr->contactGapPtr->resize(nb_gauss_pts, false);
  commonDataPtr->contactGapPtr->clear();
 
  auto t_contact_normal =
      getFTensor1FromMat<3>(*(commonDataPtr->contactNormalPtr));
  auto t_diff_normal =
      getFTensor2FromMat<3, 3>(*(commonDataPtr->contactNormalDiffPtr));
  auto t_gap_diff = getFTensor1FromMat<3>(*(commonDataPtr->contactGapDiffPtr));
  auto t_gap = getFTensor0FromVec(*(commonDataPtr->contactGapPtr));
 
  if (uPdate) {
    // FIXME: takes first gauss pt
    int roll_id = (*cache).getClosestRollerID(t_coords, t_disp);
    CHKERR mField.get_moab().tag_set_data(commonDataPtr->rollerTag, &ent, 1,
                                          &roll_id);
  }
  int tag_data;
  CHKERR mField.get_moab().tag_get_data(commonDataPtr->rollerTag, &ent, 1,
                                        &tag_data);
  (*cache).closest_roller = &(*cache).rollerDataVec.at(tag_data);
  (*cache).closest_roller_id = tag_data;
 
  const bool is_lhs = (int)getTSCtx() == 3;
  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
 
    auto t_contact = (*cache).closest_roller->getNormal(t_coords, t_disp);
    t_contact_normal(i) = t_contact(i);
    t_gap = (*cache).closest_roller->getGap(t_coords);
 
    if (is_lhs) {
      auto t_dn_du = ((*cache).closest_roller)
                         ->getDiffNormal(t_coords, t_disp, t_contact_normal);
      t_diff_normal(i, j) = t_dn_du(i, j);
      auto t_d_gap =
          (*cache).closest_roller->getdGap(t_coords, t_contact_normal);
      t_gap_diff(i) = t_d_gap(i);
    }
 
    ++t_contact_normal;
    ++t_gap;
    ++t_gap_diff;
    ++t_coords;
    ++t_disp;
    ++t_diff_normal;
  }
 
  // } //dofs
 
  MoFEMFunctionReturn(0);
}
 
struct OpPostProcContact : public DomainEleOp {
  OpPostProcContact(const std::string field_name,
                    moab::Interface &post_proc_mesh,
                    std::vector<EntityHandle> &map_gauss_pts,
                    boost::shared_ptr<CommonData> common_data_ptr);
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
private:
  moab::Interface &postProcMesh;
  std::vector<EntityHandle> &mapGaussPts;
  boost::shared_ptr<CommonData> commonDataPtr;
};
//! [Operators definitions]
 
OpPostProcContact::OpPostProcContact(
    const std::string field_name, moab::Interface &post_proc_mesh,
    std::vector<EntityHandle> &map_gauss_pts,
    boost::shared_ptr<CommonData> common_data_ptr)
    : DomainEleOp(field_name, field_name, DomainEleOp::OPROW),
      postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
      commonDataPtr(common_data_ptr) {
  // Operator is only executed for vertices
  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
  doEntities[MBVERTEX] = true;
}
 
//! [Postprocessing]
MoFEMErrorCode OpPostProcContact::doWork(int side, EntityType type,
                                         EntData &data) {
  MoFEMFunctionBegin;
 
  auto get_tag_mat = [&](const std::string name) {
    std::array<double, 9> def;
    std::fill(def.begin(), def.end(), 0);
    Tag th;
    CHKERR postProcMesh.tag_get_handle(name.c_str(), 9, MB_TYPE_DOUBLE, th,
                                       MB_TAG_CREAT | MB_TAG_SPARSE,
                                       def.data());
    return th;
  };
 
  auto get_tag_vec = [&](const std::string name) {
    std::array<double, 3> def;
    std::fill(def.begin(), def.end(), 0);
    Tag th;
    CHKERR postProcMesh.tag_get_handle(name.c_str(), 3, MB_TYPE_DOUBLE, th,
                                       MB_TAG_CREAT | MB_TAG_SPARSE,
                                       def.data());
    return th;
  };
 
  MatrixDouble3by3 mat(3, 3);
 
  auto set_matrix = [&](auto &t) -> MatrixDouble3by3 & {
    mat.clear();
    for (size_t r = 0; r != 3; ++r)
      for (size_t c = 0; c != 3; ++c)
        mat(r, c) = t(r, c);
    return mat;
  };
 
  auto set_vector = [&](auto &t) -> MatrixDouble3by3 & {
    mat.clear();
    for (size_t r = 0; r != 3; ++r)
      mat(0, r) = t(r);
    return mat;
  };
 
  auto set_tag = [&](auto th, auto gg, MatrixDouble3by3 &mat) {
    return postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
                                     &*mat.data().begin());
  };
 
  auto th_stress = get_tag_mat("SIGMA");
  auto th_div = get_tag_vec("DIV_SIGMA");
 
  size_t nb_gauss_pts = getGaussPts().size2();
  auto t_stress = getFTensor2FromMat<3, 3>(*(commonDataPtr->contactStressPtr));
  auto t_div =
      getFTensor1FromMat<3>(*(commonDataPtr->contactStressDivergencePtr));
 
  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
 
    t_stress(i, j) /= (*cache).young_modulus_inv;
    t_div(i) /= (*cache).young_modulus_inv;
 
    // FIXME: add option for more postprocessing
    // CHKERR set_tag(th_stress, gg, set_matrix(t_stress));
    CHKERR set_tag(th_div, gg, set_vector(t_div));
    ++t_stress;
    ++t_div;
  }
 
  // Tag th_time;
  // const EntityHandle root_meshset = postProcMesh.get_root_set();
  // double def_t_val = 0;
  // rval = postProcMesh.tag_get_handle("_TS_TIME_", 1, MB_TYPE_DOUBLE, th_time,
  //                                  MB_TAG_CREAT | MB_TAG_EXCL | MB_TAG_MESH,
  //                                  &def_t_val);
  // if (rval == MB_ALREADY_ALLOCATED) {
  //   // MOAB_THROW(rval);
  // } else {
  //   def_t_val = getFEMethod()->ts_t;
  //   CHKERR postProcMesh.tag_set_data(th_time, &root_meshset, 1, &def_t_val);
  //   MOAB_THROW(rval);
  // }
 
  MoFEMFunctionReturn(0);
}
//! [Postprocessing]
 
struct OpGetGaussPtsContactState : public BoundaryEleOp {
 
  OpGetGaussPtsContactState(const std::string field_name,
                            boost::shared_ptr<CommonData> common_data_ptr)
      : BoundaryEleOp(field_name, field_name, BoundaryEleOp::OPROW),
        commonDataPtr(common_data_ptr) {
    std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
    doEntities[MBTRI] = doEntities[MBQUAD] = true;
  }
 
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
    MoFEMFunctionBegin;
    const size_t nb_dofs = data.getIndices().size();
    if (nb_dofs) {
      bool empty =
          std::all_of(data.getIndices().begin(), data.getIndices().end(),
                      [](int i) { return i == -1; });
      if (empty)
        MoFEMFunctionReturnHot(0);
      size_t nb_gauss_pts = getGaussPts().size2();
      auto t_traction =
          getFTensor1FromMat<3>(*(commonDataPtr->contactTractionPtr));
 
      auto &data_vec = commonDataPtr->stateArrayCont;
      auto t_contact_normal =
          getFTensor1FromMat<3>(*(commonDataPtr->contactNormalPtr));
      auto t_gap = getFTensor0FromVec(*(commonDataPtr->contactGapPtr));
      for (int gg = 0; gg != nb_gauss_pts; ++gg) {
 
        const double trac = normal_traction(t_traction, t_contact_normal);
 
        if (t_gap <= (*cache).cn_cont * trac)
          data_vec[0] += 1;
        data_vec[1] += 1;
 
        ++t_contact_normal;
        ++t_gap;
        ++t_traction;
      }
    }
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
};
 
struct OpGetContactArea : public BoundaryEleOp {
 
  OpGetContactArea(const std::string field_name,
                   boost::shared_ptr<CommonData> common_data_ptr)
      : BoundaryEleOp(field_name, field_name, BoundaryEleOp::OPROW),
        commonDataPtr(common_data_ptr) {
    std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
    doEntities[MBTRI] = doEntities[MBQUAD] = true;
  }
 
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
    MoFEMFunctionBegin;
    const size_t nb_dofs = data.getIndices().size();
    if (nb_dofs) {
      auto ent = getFEEntityHandle();
      bool empty =
          std::all_of(data.getIndices().begin(), data.getIndices().end(),
                      [](int i) { return i == -1; });
      if (empty)
        MoFEMFunctionReturnHot(0);
 
      size_t nb_gauss_pts = getGaussPts().size2();
      auto t_w = getFTensor0IntegrationWeight();
      auto t_traction =
          getFTensor1FromMat<3>(*(commonDataPtr->contactTractionPtr));
 
      auto &data_vec = commonDataPtr->stateArrayArea;
      auto t_contact_normal =
          getFTensor1FromMat<3>(*(commonDataPtr->contactNormalPtr));
      auto t_gap = getFTensor0FromVec(*(commonDataPtr->contactGapPtr));
 
      int id;
      CHKERR getPtrFE()->mField.get_moab().tag_get_data(
          commonDataPtr->rollerTag, &ent, 1, &id);
      auto &l_reactions = commonDataPtr->bodyReactions;
 
      for (int gg = 0; gg != nb_gauss_pts; ++gg) {
        const double alpha = t_w * getMeasure();
        const double trac = normal_traction(t_traction, t_contact_normal);
        const auto sign_w = sign(w(t_gap, trac));
 
        for (int dd = 0; dd != 3; ++dd)
          l_reactions[3 * id + dd] +=
              t_traction(dd) * alpha / (*cache).young_modulus_inv;
 
        if (sign_w < 0)
          data_vec[0] += alpha;
        data_vec[1] += alpha;
 
        ++t_contact_normal;
        ++t_gap;
        ++t_traction;
        ++t_w;
      }
    }
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
};
 
}; // namespace OpContactTools