MortarContactInterface.cpp

Go to the documentation of this file.

 
 
#ifdef __cplusplus
extern "C" {
#endif
#include <phg-quadrule/quad.h>
#ifdef __cplusplus
}
#endif
 
#include <Mortar.hpp>
 
std::vector<BitRefLevel> &MortarContactInterface::getBitRefLevelVector() {
  return bitLevels;
}
 
boost::ptr_deque<MoFEM::ForcesAndSourcesCore::UserDataOperator> &
MortarContactInterface::getPostProcSimpleContactPipeline() {
  return fePostProcSimpleContact->getOpPtrVector();
}
 
MoFEMErrorCode MortarContactInterface::getCommandLineParameters() {
  MoFEMFunctionBeginHot;
 
  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Mortar Config", "none");
 
  CHKERR PetscOptionsInt(
      "-order_contact",
      "approximation order of spatial positions in contact interface", "", 1,
      &orderContact, PETSC_NULL);
  CHKERR PetscOptionsInt("-ho_levels_num", "number of higher order levels", "",
                         0, &nbHoLevels, PETSC_NULL);
  CHKERR PetscOptionsInt("-order_lambda",
                         "approximation order of Lagrange multipliers", "", 1,
                         &orderLambda, PETSC_NULL);
  isPartitioned = PETSC_FALSE;
  CHKERR PetscOptionsBool("-is_partitioned",
                          "set if mesh is partitioned (this result that each "
                          "process keeps only part of the mes",
                          "", PETSC_FALSE, &isPartitioned, PETSC_NULL);
 
  CHKERR PetscOptionsReal("-cn_value", "default regularisation cn value", "",
                          1., &cnValue, PETSC_NULL);
 
  CHKERR PetscOptionsReal("-post_proc_gap_tol",
                          "gap tolerance for post processing", "", 0.,
                          &postProcGapTol, PETSC_NULL);
 
  CHKERR PetscOptionsBool("-is_newton_cotes",
                          "set if Newton-Cotes quadrature rules are used", "",
                          PETSC_FALSE, &isNewtonCotes, PETSC_NULL);
  CHKERR PetscOptionsBool("-convect", "set to convect integration pts", "",
                          PETSC_FALSE, &convectPts, PETSC_NULL);
  CHKERR PetscOptionsBool("-print_contact_state",
                          "output number of active gp at every iteration", "",
                          PETSC_FALSE, &printContactState, PETSC_NULL);
  CHKERR PetscOptionsBool("-alm_flag", "if set use ALM, if not use C-function",
                          "", PETSC_FALSE, &almFlag, PETSC_NULL);
  CHKERR PetscOptionsBool("-debug", "output debug files", "", PETSC_FALSE,
                          &dEbug, PETSC_NULL);
 
  int order = 1;
  CHKERR PetscOptionsInt("-order", "approximation order of spatial positions",
                         "", 1, &order, PETSC_NULL);
  CHKERR this->setPositionFieldOrder(order);
 
  ierr = PetscOptionsEnd();
  CHKERRQ(ierr);
 
  MoFEMFunctionReturnHot(0);
}
 
MoFEMErrorCode MortarContactInterface::setupElementEntities() {
  MoFEMFunctionBegin;
  auto &moab = mField.get_moab();
 
  bitLevels.push_back(BitRefLevel().set(0));
  CHKERR mField.getInterface<BitRefManager>()->setBitRefLevelByDim(
      0, 3, bitLevels.back());
 
  Range meshset_level0;
  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByRefLevel(
      bitLevels.back(), BitRefLevel().set(), meshset_level0);
  MOFEM_LOG("WORLD", Sev::verbose)
      << "Meshset_level0 " << meshset_level0.size();
 
  contactCommondataMultiIndex = boost::make_shared<contactMIndex>();
 
  CHKERR findContactSurfacePairs();
 
  for (auto &contact_surface_pair : contactSurfacePairs) {
    slaveTris.merge(contact_surface_pair.first);
    masterTris.merge(contact_surface_pair.second);
 
    ContactSearchKdTree contact_search_kd_tree(mField);
 
    CHKERR contact_search_kd_tree.buildTree(contact_surface_pair.second);
 
    CHKERR contact_search_kd_tree.contactSearchAlgorithm(
        contact_surface_pair.second, contact_surface_pair.first,
        contactCommondataMultiIndex, contactPrisms);
  }
 
  if (!masterTris.empty() && isPartitioned)
    SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
            "Partitioned mesh is not supported with mortar contact");
 
  EntityHandle meshset_slave_master_prisms;
  CHKERR moab.create_meshset(MESHSET_SET, meshset_slave_master_prisms);
  CHKERR
  moab.add_entities(meshset_slave_master_prisms, contactPrisms);
  CHKERR mField.getInterface<BitRefManager>()->setBitRefLevelByDim(
      meshset_slave_master_prisms, 3, bitLevels.back());
 
  if (dEbug)
    CHKERR moab.write_mesh("slave_master_prisms.vtk",
                           &meshset_slave_master_prisms, 1);
 
  Range tris_from_prism;
  // find actual masters and slave used
  CHKERR mField.get_moab().get_adjacencies(
      contactPrisms, 2, true, tris_from_prism, moab::Interface::UNION);
 
  tris_from_prism = tris_from_prism.subset_by_type(MBTRI);
  slaveTris = intersect(tris_from_prism, slaveTris);
 
  if (dEbug) {
    EntityHandle mshs_surf_slave;
    CHKERR mField.get_moab().create_meshset(MESHSET_SET, mshs_surf_slave);
    CHKERR mField.get_moab().add_entities(mshs_surf_slave, slaveTris);
    CHKERR mField.get_moab().write_mesh("surf_slave.vtk", &mshs_surf_slave, 1);
    // EntityHandle meshset_tri_slave;
    // CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset_tri_slave);
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode MortarContactInterface::addElementFields() {
  MoFEMFunctionBegin;
  auto &moab = mField.get_moab();
 
  if (!contactCommondataMultiIndex)
    CHKERR this->setupElementEntities();
 
  if (!mField.check_field(meshNodeField)) {
 
    CHKERR mField.add_field(meshNodeField, H1, AINSWORTH_LEGENDRE_BASE, 3,
                            MB_TAG_SPARSE, MF_ZERO);
    CHKERR mField.add_ents_to_field_by_type(0, MBTET, meshNodeField);
    CHKERR mField.set_field_order(0, MBTET, meshNodeField, 1);
    CHKERR mField.set_field_order(0, MBTRI, meshNodeField, 1);
    CHKERR mField.set_field_order(0, MBEDGE, meshNodeField, 1);
    CHKERR mField.set_field_order(0, MBVERTEX, meshNodeField, 1);
  }
 
  if (!mField.check_field(positionField)) {
    // FIXME: TODO: we need to discuss how to
    if (oRder < 0)
      SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
              "set appropriate position field order, "
              "setPositionFieldOrder(int order)");
 
    // Declare problem add entities (by tets) to the field
    CHKERR mField.add_field(positionField, H1, AINSWORTH_LEGENDRE_BASE, 3,
                            MB_TAG_SPARSE, MF_ZERO);
    CHKERR mField.add_ents_to_field_by_type(0, MBTET, positionField);
    CHKERR mField.set_field_order(0, MBTET, positionField, oRder);
    CHKERR mField.set_field_order(0, MBTRI, positionField, oRder);
    CHKERR mField.set_field_order(0, MBEDGE, positionField, oRder);
    CHKERR mField.set_field_order(0, MBVERTEX, positionField, 1);
  }
 
  CHKERR mField.add_field("LAGMULT", H1, AINSWORTH_LEGENDRE_BASE, 1,
                          MB_TAG_SPARSE, MF_ZERO);
  CHKERR mField.add_ents_to_field_by_type(slaveTris, MBTRI, "LAGMULT");
  CHKERR mField.set_field_order(0, MBTRI, "LAGMULT", orderLambda);
  CHKERR mField.set_field_order(0, MBEDGE, "LAGMULT", orderLambda);
  CHKERR mField.set_field_order(0, MBVERTEX, "LAGMULT", 1);
 
  auto set_contact_order = [&]() {
    MoFEMFunctionBegin;
    Range contact_tris, contact_edges;
    CHKERR moab.get_adjacencies(contactPrisms, 2, false, contact_tris,
                                moab::Interface::UNION);
    contact_tris = contact_tris.subset_by_type(MBTRI);
    CHKERR moab.get_adjacencies(contact_tris, 1, false, contact_edges,
                                moab::Interface::UNION);
    Range ho_ents;
    ho_ents.merge(contact_tris);
    ho_ents.merge(contact_edges);
    for (int ll = 0; ll != nbHoLevels; ll++) {
      Range ents, verts, tets;
      CHKERR moab.get_connectivity(ho_ents, verts, true);
      CHKERR moab.get_adjacencies(verts, 3, false, tets,
                                  moab::Interface::UNION);
      tets = tets.subset_by_type(MBTET);
      for (auto d : {1, 2}) {
        CHKERR moab.get_adjacencies(tets, d, false, ents,
                                    moab::Interface::UNION);
      }
      ho_ents = unite(ho_ents, ents);
      ho_ents = unite(ho_ents, tets);
    }
 
    CHKERR mField.set_field_order(ho_ents, positionField, orderContact);
 
    MoFEMFunctionReturn(0);
  };
 
  if (orderContact > oRder)
    CHKERR set_contact_order();
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode MortarContactInterface::setOperators() { return 0; }
 
MoFEMErrorCode MortarContactInterface::createElements() {
  MoFEMFunctionBegin;
 
  contactProblemPtr = boost::make_shared<MortarContactProblem>(
      mField, contactCommondataMultiIndex, boost::make_shared<double>(cnValue),
      isNewtonCotes);
 
  // add fields to the global matrix by adding the element
  contactProblemPtr->addMortarContactElement(
      "CONTACT_ELEM", positionField, "LAGMULT", contactPrisms, meshNodeField);
  contactProblemPtr->addPostProcContactElement(
      "CONTACT_POST_PROC", positionField, "LAGMULT", meshNodeField, slaveTris);
 
  MoFEMFunctionReturn(0);
}
 
BitRefLevel MortarContactInterface::getBitRefLevel() {
  return this->getBitRefLevelVector().back();
};
 
MoFEMErrorCode MortarContactInterface::addElementsToDM(SmartPetscObj<DM> dm) {
  MoFEMFunctionBeginHot;
 
  this->dM = dm;
 
  CHKERR DMMoFEMAddElement(this->dM, "CONTACT_ELEM");
  CHKERR DMMoFEMAddElement(this->dM, "CONTACT_POST_PROC");
 
  auto simple = mField.getInterface<Simple>();
  simple->getOtherFiniteElements().push_back("CONTACT_ELEM");
  simple->getOtherFiniteElements().push_back("CONTACT_POST_PROC");
 
  MoFEMFunctionReturnHot(0);
}
 
template <typename T, bool RHS>
MoFEMErrorCode MortarContactInterface::setupSolverImpl(const TSType type) {
  MoFEMFunctionBegin;
 
  if (!dM)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "DM not defined, addElementsToDM(dm) first");
 
  auto &moab = mField.get_moab();
 
  auto make_contact_element = [&]() {
    return boost::make_shared<MortarContactProblem::MortarContactElement>(
        mField, contactCommondataMultiIndex, positionField, meshNodeField);
  };
 
  auto make_convective_master_element = [&]() {
    return boost::make_shared<
        MortarContactProblem::MortarConvectMasterContactElement>(
        mField, contactCommondataMultiIndex, positionField, meshNodeField);
  };
 
  auto make_convective_slave_element = [&]() {
    return boost::make_shared<
        MortarContactProblem::MortarConvectSlaveContactElement>(
        mField, contactCommondataMultiIndex, positionField, meshNodeField);
  };
 
  auto make_contact_common_data = [&]() {
    return boost::make_shared<MortarContactProblem::CommonDataMortarContact>(
        mField);
  };
 
  auto get_contact_rhs = [&](auto make_element, bool is_alm = false,
                             bool is_displacements = false) {
    auto fe_rhs_simple_contact = make_element();
    auto common_data_mortar_contact = make_contact_common_data();
    if (printContactState) {
      fe_rhs_simple_contact->contactStateVec =
          common_data_mortar_contact->gaussPtsStateVec;
    }
    fe_rhs_simple_contact->getOpPtrVector().push_back(
        new OpSetBc(positionField, true, mBoundaryMarker));
    contactProblemPtr->setContactOperatorsRhs(
        fe_rhs_simple_contact, common_data_mortar_contact, positionField,
        "LAGMULT", is_alm, false, "", is_displacements);
    fe_rhs_simple_contact->getOpPtrVector().push_back(
        new OpUnSetBc(positionField));
    return fe_rhs_simple_contact;
  };
 
  auto get_master_traction_rhs = [&](auto make_element, bool is_alm = false,
                                     bool is_displacements = false) {
    auto fe_rhs_simple_contact = make_element();
    fe_rhs_simple_contact->getOpPtrVector().push_back(
        new OpSetBc(positionField, true, mBoundaryMarker));
    auto common_data_mortar_contact = make_contact_common_data();
    contactProblemPtr->setMasterForceOperatorsRhs(
        fe_rhs_simple_contact, common_data_mortar_contact, positionField,
        "LAGMULT", is_alm, false, "", is_displacements);
    fe_rhs_simple_contact->getOpPtrVector().push_back(
        new OpUnSetBc(positionField));
    return fe_rhs_simple_contact;
  };
 
  auto get_master_traction_lhs = [&](auto make_element, bool is_alm = false,
                                     bool is_displacements = false) {
    auto fe_lhs_simple_contact = make_element();
    auto common_data_mortar_contact = make_contact_common_data();
    fe_lhs_simple_contact->getOpPtrVector().push_back(
        new OpSetBc(positionField, true, mBoundaryMarker));
    contactProblemPtr->setMasterForceOperatorsLhs(
        fe_lhs_simple_contact, common_data_mortar_contact, positionField,
        "LAGMULT", is_alm, false, "", is_displacements);
    fe_lhs_simple_contact->getOpPtrVector().push_back(
        new OpUnSetBc(positionField));
    return fe_lhs_simple_contact;
  };
 
  auto get_master_help_traction_lhs = [&](auto make_element,
                                          bool is_alm = false,
                                          bool is_displacements = false) {
    auto fe_lhs_simple_contact = make_element();
    auto common_data_mortar_contact = make_contact_common_data();
    fe_lhs_simple_contact->getOpPtrVector().push_back(
        new OpSetBc(positionField, true, mBoundaryMarker));
    contactProblemPtr->setMasterForceOperatorsLhs(
        fe_lhs_simple_contact, common_data_mortar_contact, positionField,
        "LAGMULT", is_alm, false, "", is_displacements);
    fe_lhs_simple_contact->getOpPtrVector().push_back(
        new OpUnSetBc(positionField));
    return fe_lhs_simple_contact;
  };
 
  auto get_contact_lhs = [&](auto make_element, bool is_alm = false,
                             bool is_displacements = false) {
    auto fe_lhs_simple_contact = make_element();
    auto common_data_mortar_contact = make_contact_common_data();
    fe_lhs_simple_contact->getOpPtrVector().push_back(
        new OpSetBc(positionField, true, mBoundaryMarker));
    contactProblemPtr->setContactOperatorsLhs(
        fe_lhs_simple_contact, common_data_mortar_contact, positionField,
        "LAGMULT", is_alm, false, "", is_displacements);
    fe_lhs_simple_contact->getOpPtrVector().push_back(
        new OpUnSetBc(positionField));
    return fe_lhs_simple_contact;
  };
 
  auto get_contact_help_lhs = [&](auto make_element, bool is_alm = false,
                                  bool is_displacements = false) {
    auto fe_lhs_simple_contact = make_element();
    auto common_data_mortar_contact = make_contact_common_data();
    fe_lhs_simple_contact->getOpPtrVector().push_back(
        new OpSetBc(positionField, true, mBoundaryMarker));
    contactProblemPtr->setContactOperatorsLhs(
        fe_lhs_simple_contact, common_data_mortar_contact, positionField,
        "LAGMULT", is_alm, false, "", is_displacements);
    fe_lhs_simple_contact->getOpPtrVector().push_back(
        new OpUnSetBc(positionField));
    return fe_lhs_simple_contact;
  };
 
  auto set_solver_pipelines = [&](PetscErrorCode (*function)(
                                      DM, const std::string,
                                      boost::shared_ptr<MoFEM::FEMethod>,
                                      boost::shared_ptr<MoFEM::BasicMethod>,
                                      boost::shared_ptr<MoFEM::BasicMethod>),
                                  PetscErrorCode (*jacobian)(
                                      DM, const std::string,
                                      boost::shared_ptr<MoFEM::FEMethod>,
                                      boost::shared_ptr<MoFEM::BasicMethod>,
                                      boost::shared_ptr<MoFEM::BasicMethod>)) {
    MoFEMFunctionBeginHot;
 
    if (convectPts) {
      if (RHS) {
 
        CHKERR function(this->dM, "CONTACT_ELEM",
                        get_contact_rhs(make_convective_master_element, almFlag,
                                        isDisplacementField),
                        PETSC_NULL, PETSC_NULL);
        CHKERR function(this->dM, "CONTACT_ELEM",
                        get_master_traction_rhs(make_convective_slave_element,
                                                almFlag, isDisplacementField),
                        PETSC_NULL, PETSC_NULL);
      } else {
        if (jacobian) {
 
          CHKERR jacobian(this->dM, "CONTACT_ELEM",
                          get_contact_help_lhs(make_convective_master_element,
                                               almFlag, isDisplacementField),
                          PETSC_NULL, PETSC_NULL);
 
          CHKERR jacobian(
              this->dM, "CONTACT_ELEM",
              get_master_help_traction_lhs(make_convective_slave_element),
              PETSC_NULL, PETSC_NULL);
        }
      }
    } else {
      if (RHS) {
 
        CHKERR function(
            this->dM, "CONTACT_ELEM",
            get_contact_rhs(make_contact_element, almFlag, isDisplacementField),
            PETSC_NULL, PETSC_NULL);
        CHKERR function(this->dM, "CONTACT_ELEM",
                        get_master_traction_rhs(make_contact_element, almFlag,
                                                isDisplacementField),
                        PETSC_NULL, PETSC_NULL);
      } else {
        if (jacobian) {
 
          CHKERR jacobian(this->dM, "CONTACT_ELEM",
                          get_contact_lhs(make_contact_element, almFlag,
                                          isDisplacementField),
                          PETSC_NULL, PETSC_NULL);
          CHKERR jacobian(this->dM, "CONTACT_ELEM",
                          get_master_traction_lhs(make_contact_element, almFlag,
                                                  isDisplacementField),
                          PETSC_NULL, PETSC_NULL);
        }
      }
    }
 
    MoFEMFunctionReturnHot(0);
  };
 
  if constexpr (std::is_same_v<T, SNES>) {
 
    CHKERR set_solver_pipelines(&DMMoFEMSNESSetFunction,
                                &DMMoFEMSNESSetJacobian);
  } else if (std::is_same_v<T, TS>) {
 
    switch (type) {
    case IM:
      CHKERR set_solver_pipelines(&DMMoFEMTSSetIFunction,
                                  &DMMoFEMTSSetIJacobian);
      break;
    case IM2:
      CHKERR set_solver_pipelines(&DMMoFEMTSSetI2Function,
                                  &DMMoFEMTSSetI2Jacobian);
      break;
    case EX:
      CHKERR set_solver_pipelines(&DMMoFEMTSSetRHSFunction,
                                  nullptr);
      break;
    default:
      SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
              "This TS is not yet implemented for contact");
      break;
    }
 
  } else
    static_assert(!std::is_same_v<T, KSP>,
                  "this solver has not been implemented for contact yet");
 
  if (!commonDataPostProc) {
    commonDataPostProc = make_contact_common_data();
    if (convectPts)
      fePostProcSimpleContact = make_convective_master_element();
    else
      fePostProcSimpleContact = make_contact_element();
 
    CHKERR findPostProcSurface();
 
    contactProblemPtr->setContactOperatorsForPostProc(
        fePostProcSimpleContact, commonDataPostProc, mField, positionField,
        "LAGMULT", moabPostProcMesh, almFlag, false, "", false, postProcSurface,
        postProcGapTol);
  }
 
  MoFEMFunctionReturn(0);
}
 
template <>
MoFEMErrorCode
MortarContactInterface::setupSolverFunction<TS>(const TSType type) {
  MoFEMFunctionBeginHot;
  CHKERR setupSolverImpl<TS, true>(type);
  MoFEMFunctionReturnHot(0);
}
template <>
MoFEMErrorCode
MortarContactInterface::setupSolverFunction<SNES>(const TSType type) {
  MoFEMFunctionBeginHot;
  CHKERR setupSolverImpl<SNES, true>(type);
  MoFEMFunctionReturnHot(0);
}
template <>
MoFEMErrorCode
MortarContactInterface::setupSolverJacobian<SNES>(const TSType type) {
  MoFEMFunctionBeginHot;
  CHKERR setupSolverImpl<SNES, false>(type);
  MoFEMFunctionReturnHot(0);
}
template <>
MoFEMErrorCode
MortarContactInterface::setupSolverJacobian<TS>(const TSType type) {
  MoFEMFunctionBeginHot;
  CHKERR setupSolverImpl<TS, false>(type);
  MoFEMFunctionReturnHot(0);
}
 
MoFEMErrorCode MortarContactInterface::setupSolverFunctionSNES() {
  MoFEMFunctionBegin;
  CHKERR this->setupSolverFunction<SNES>();
  MoFEMFunctionReturn(0);
}
MoFEMErrorCode MortarContactInterface::setupSolverJacobianSNES() {
  MoFEMFunctionBegin;
  CHKERR this->setupSolverJacobian<SNES>();
  MoFEMFunctionReturn(0);
}
MoFEMErrorCode
MortarContactInterface::setupSolverFunctionTS(const TSType type) {
  MoFEMFunctionBegin;
  CHKERR this->setupSolverFunction<TS>(type);
  MoFEMFunctionReturn(0);
}
MoFEMErrorCode
MortarContactInterface::setupSolverJacobianTS(const TSType type) {
  MoFEMFunctionBegin;
  CHKERR this->setupSolverJacobian<TS>(type);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode MortarContactInterface::postProcessContactState(int step) {
  MoFEMFunctionBegin;
 
  moabPostProcMesh.delete_mesh();
  if (!commonDataPostProc)
    SETERRQ(
        PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
        "you need to setup the solver first to use postProcessContactState");
  CHKERR VecZeroEntries(commonDataPostProc->gaussPtsStateVec);
  CHKERR VecZeroEntries(commonDataPostProc->contactAreaVec);
  CHKERR DMoFEMLoopFiniteElements(dM, "CONTACT_ELEM", fePostProcSimpleContact);
 
  auto get_contact_data = [&](auto vec, std::array<double, 2> &data) {
    MoFEMFunctionBegin;
    CHKERR VecAssemblyBegin(vec);
    CHKERR VecAssemblyEnd(vec);
    const double *array;
    CHKERR VecGetArrayRead(vec, &array);
    if (mField.get_comm_rank() == 0) {
      for (int i : {0, 1})
        data[i] = array[i];
    }
    CHKERR VecRestoreArrayRead(vec, &array);
    MoFEMFunctionReturn(0);
  };
 
  CHKERR get_contact_data(commonDataPostProc->gaussPtsStateVec, nbGaussPts);
  CHKERR get_contact_data(commonDataPostProc->contactAreaVec, contactArea);
 
  if (mField.get_comm_rank() == 0) {
    MOFEM_LOG_C("WORLD", Sev::verbose, "Active gauss pts: %d out of %d",
                (int)nbGaussPts[0], (int)nbGaussPts[1]);
         
    if (!postProcSurface.empty())
      MOFEM_LOG("WORLD", Sev::inform)
          << "Computing contact area only on the post process surface";
    if (postProcGapTol > std::numeric_limits<double>::epsilon())
      MOFEM_LOG_C("WORLD", Sev::inform, "Using gap tolerance %3.3f",
                  postProcGapTol);
 
    MOFEM_LOG_C("WORLD", Sev::inform, "Active contact area: %3.3f out of %3.3f",
                contactArea[0], contactArea[1]);
  }
 
  auto out_file_name = "out_contact_integ_pts_" + to_string(step) + ".h5m";
  MOFEM_LOG_C("WORLD", Sev::verbose, "out file %s", out_file_name.c_str());
  CHKERR moabPostProcMesh.write_file(out_file_name.c_str(), "MOAB",
                                     "PARALLEL=WRITE_PART");
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode MortarContactInterface::postProcessElement(int step) {
  MoFEMFunctionBegin;
 
  if (masterTris.empty())
    MoFEMFunctionReturnHot(0);
 
  this->postProcessContactState(step);
 
  if (!postProcContactPtr) {
    postProcContactPtr = boost::make_shared<PostProcFaceOnRefinedMesh>(mField);
    CHKERR postProcContactPtr->generateReferenceElementMesh();
    CHKERR addHOOpsFace3D(meshNodeField, *postProcContactPtr, false, false);
    CHKERR postProcContactPtr->addFieldValuesPostProc("LAGMULT");
    CHKERR postProcContactPtr->addFieldValuesPostProc(positionField);
    CHKERR postProcContactPtr->addFieldValuesPostProc(meshNodeField);
  }
  CHKERR DMoFEMLoopFiniteElements(this->dM, "CONTACT_POST_PROC",
                                  postProcContactPtr);
 
  auto out_file_name = "out_contact_" + to_string(step) + ".h5m";
  MOFEM_LOG_C("WORLD", Sev::verbose, "out file %s", out_file_name.c_str());
  CHKERR postProcContactPtr->writeFile(out_file_name.c_str());
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode MortarContactInterface::setPositionFieldOrder(int order) {
  MoFEMFunctionBeginHot;
  oRder = order;
  MoFEMFunctionReturnHot(0);
}
 
MoFEMErrorCode MortarContactInterface::updateElementVariables() {
  MoFEMFunctionBeginHot;
  MoFEMFunctionReturnHot(0);
}
 
MoFEMErrorCode MortarContactInterface::findContactSurfacePairs() {
  MoFEMFunctionBegin;
 
  std::map<std::string, Range> slave_surface_map, master_surface_map;
 
  auto fill_surface_map = [&](auto &surface_map, std::string surface_name) {
    MoFEMFunctionBegin;
 
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
      if (it->getName().compare(0, surface_name.size(), surface_name) == 0) {
        Range tris;
        CHKERR mField.get_moab().get_entities_by_type(it->meshset, MBTRI, tris,
                                                      true);
        std::string key = it->getName().substr(surface_name.size());
        auto res = surface_map.insert(std::pair<std::string, Range>(key, tris));
        if (!res.second) {
          SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                   "Contact surface name %s used more than once",
                   it->getName().c_str());
        }
      }
    }
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR fill_surface_map(slave_surface_map, "MORTAR_SLAVE");
  CHKERR fill_surface_map(master_surface_map, "MORTAR_MASTER");
 
  for (const auto &[key, slave_range] : slave_surface_map) {
    auto master_search = master_surface_map.find(key);
    if (master_search != master_surface_map.end()) {
      contactSurfacePairs.push_back(
          std::pair<Range, Range>(slave_range, master_search->second));
      master_surface_map.erase(master_search);
    } else {
      SETERRQ2(
          PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
          "Cannot find pairing contact surface MORTAR_MASTER%s for existing "
          "surface MORTAR_SLAVE%s",
          key.c_str(), key.c_str());
    }
  }
 
  if (!master_surface_map.empty()) {
    std::string key = master_surface_map.begin()->first;
    SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
             "Cannot find pairing contact surface MORTAR_SLAVE%s for existing "
             "surface MORTAR_MASTER%s",
             key.c_str(), key.c_str());
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode MortarContactInterface::findPostProcSurface() {
  MoFEMFunctionBegin;
 
  std::string surface_name = "MORTAR_POST_PROC";
 
  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
    if (it->getName().compare(0, surface_name.size(), surface_name) == 0) {
      Range tris;
      CHKERR mField.get_moab().get_entities_by_type(it->meshset, MBTRI, tris,
                                                    true);
      postProcSurface.merge(tris);
    }
  }
 
  MoFEMFunctionReturn(0);
}