PostProcContact.hpp

Go to the documentation of this file.

/**
 * \file PostProcContact.hpp
 *
 *
 * @copyright Copyright (c) 2023
 */
 
namespace ContactOps {
 
template <int DIM> struct PostProcEleByDim;
 
template <> struct PostProcEleByDim<2> {
  using PostProcEleDomain = PostProcBrokenMeshInMoabBaseCont<DomainEle>;
  using PostProcEleBdy = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;
  using SideEle = PipelineManager::ElementsAndOpsByDim<2>::FaceSideEle;
};
 
template <> struct PostProcEleByDim<3> {
  using PostProcEleDomain = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;
  using PostProcEleBdy = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;
  using SideEle = PipelineManager::ElementsAndOpsByDim<3>::FaceSideEle;
};
 
using PostProcEleDomain = PostProcEleByDim<SPACE_DIM>::PostProcEleDomain;
using SideEle = PostProcEleByDim<SPACE_DIM>::SideEle;
using PostProcEleBdy = PostProcEleByDim<SPACE_DIM>::PostProcEleBdy;
 
struct Monitor : public FEMethod {
 
  Monitor(SmartPetscObj<DM> &dm, double scale,
          boost::shared_ptr<MFrontInterface> mfront_interface_ptr = nullptr,
          bool is_axisymmetric = false, bool is_large_strain = false)
      : dM(dm), moabVertex(mbVertexPostproc), sTEP(0),
        mfrontInterfacePtr(mfront_interface_ptr),
        isLargeStrain(is_large_strain) {
 
    MoFEM::Interface *m_field_ptr;
    CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
 
    using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
    struct OpScale : public ForcesAndSourcesCore::UserDataOperator {
      OpScale(boost::shared_ptr<MatrixDouble> m_ptr, double s)
          : ForcesAndSourcesCore::UserDataOperator(NOSPACE, OPSPACE),
            mPtr(m_ptr), scale(s) {}
      MoFEMErrorCode doWork(int, EntityType, EntitiesFieldData::EntData &) {
        *mPtr *= 1. / scale;
        return 0;
      }
 
    private:
      boost::shared_ptr<MatrixDouble> mPtr;
      double scale;
    };
 
    auto push_domain_ops = [&](auto &pip) {
      CHK_THROW_MESSAGE((AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
                            pip, {H1, HDIV}, "GEOMETRY")),
                        "Apply base transform");
      // Define variant type that holds either HookeOps or HenckyOps
      // CommonData.
      using CommonDataVariant =
          std::variant<boost::shared_ptr<HookeOps::CommonData>,
                       boost::shared_ptr<HenckyOps::CommonData>>;
      auto contact_stress_ptr = boost::make_shared<MatrixDouble>();
      pip.push_back(new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
          "SIGMA", contact_stress_ptr));
 
      if (!isLargeStrain) { // use HookeOps
        auto hooke_common_data_ptr =
            HookeOps::commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
                *m_field_ptr, pip, "U", "MAT_ELASTIC", Sev::inform, scale);
 
        pip.push_back(new OpScale(contact_stress_ptr, scale));
        return std::make_tuple(CommonDataVariant(hooke_common_data_ptr),
                               contact_stress_ptr);
      } else { // use HenckyOps
        auto hencky_common_data_ptr =
            HenckyOps::commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
                *m_field_ptr, pip, "U", "MAT_ELASTIC", Sev::inform, scale);
        pip.push_back(new OpScale(contact_stress_ptr, scale));
        return std::make_tuple(CommonDataVariant(hencky_common_data_ptr),
                               contact_stress_ptr);
      }
    };
    auto push_bdy_ops = [&](auto &pip) {
      // evaluate traction
      auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
      pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
          "U", common_data_ptr->contactDispPtr()));
      pip.push_back(new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
          "SIGMA", common_data_ptr->contactTractionPtr()));
      pip.push_back(new OpScale(common_data_ptr->contactTractionPtr(), scale));
      using C = ContactIntegrators<BoundaryEleOp>;
      pip.push_back(new typename C::template OpEvaluateSDF<SPACE_DIM, GAUSS>(
          common_data_ptr));
      return common_data_ptr;
    };
 
    auto get_domain_pip = [&](auto &pip)
        -> boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> & {
      if constexpr (SPACE_DIM == 3) {
        auto op_loop_side = new OpLoopSide<SideEle>(
            *m_field_ptr, "dFE", SPACE_DIM, Sev::noisy,
            boost::make_shared<
                ForcesAndSourcesCore::UserDataOperator::AdjCache>());
        pip.push_back(op_loop_side);
        return op_loop_side->getOpPtrVector();
      } else {
        return pip;
      }
    };
 
    auto get_post_proc_domain_fe = [&]() {
      auto post_proc_fe =
          boost::make_shared<PostProcEleDomain>(*m_field_ptr, postProcMesh);
      auto &pip = post_proc_fe->getOpPtrVector();
      // Get common data pointer (either Hooke or Hencky)
      auto [hencky_or_hooke_common_data_ptr, contact_stress_ptr] =
          push_domain_ops(get_domain_pip(pip));
 
      auto u_ptr = boost::make_shared<MatrixDouble>();
      pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
      auto X_ptr = boost::make_shared<MatrixDouble>();
      pip.push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", X_ptr));
      /**
       * @brief  Adds operators for evaluating stress and strain based on the
       * material model: i.e. Visit the variant to handle either
       * HookeOps::CommonData or HenckyOps::CommonData. i.e.
       *   - For `HookeOps::CommonData`, maps Cauchy stress and small strain.
       *   - For `HenckyOps::CommonData`, maps first Piola-Kirchhoff stress and
       * Hencky strain.
       * @return A commod data pointer to the configured post-processing
       * hencky_or_hooke_common_data_ptr.
       */
      std::visit(
          [&](auto &common_data_ptr) {
            if constexpr (std::is_same_v<
                              std::decay_t<decltype(common_data_ptr)>,
                              boost::shared_ptr<HookeOps::CommonData>>) {
              pip.push_back(new OpPPMap(
                  post_proc_fe->getPostProcMesh(),
                  post_proc_fe->getMapGaussPts(), {},
                  {{"U", u_ptr}, {"GEOMETRY", X_ptr}},
                  {{"SIGMA", contact_stress_ptr},
                   {"G", common_data_ptr->matGradPtr}}, // Displacement gradient
                  {{"STRESS", common_data_ptr->getMatCauchyStress()},
                   {"STRAIN", common_data_ptr->getMatStrain()}}));
            }
            // Handle HenckyOps::CommonData
            else if constexpr (std::is_same_v<
                                   std::decay_t<decltype(common_data_ptr)>,
                                   boost::shared_ptr<HenckyOps::CommonData>>) {
 
              pip.push_back(new OpPPMap(
                  post_proc_fe->getPostProcMesh(),
                  post_proc_fe->getMapGaussPts(), {},
                  {{"U", u_ptr}, {"GEOMETRY", X_ptr}},
                  {{"SIGMA", contact_stress_ptr},
                   {"G", common_data_ptr->matGradPtr},
                   {"PK1", common_data_ptr->getMatFirstPiolaStress()}
 
                  },
                  {{"HENCKY_STRAIN", common_data_ptr->getMatLogC()}}));
            }
          },
          hencky_or_hooke_common_data_ptr);
 
      if (SPACE_DIM == 3) {
 
        CHK_THROW_MESSAGE((AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
                              pip, {HDIV}, "GEOMETRY")),
                          "Apply transform");
        auto common_data_ptr = push_bdy_ops(pip);
 
        pip.push_back(
 
            new OpPPMap(
 
                post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
 
                {{"SDF", common_data_ptr->sdfPtr()},
                 {"CONSTRAINT_CONTACT", common_data_ptr->constraintPtr()}},
 
                {
 
                    {"TRACTION_CONTACT", common_data_ptr->contactTractionPtr()},
                    {"GRAD_SDF", common_data_ptr->gradSdfPtr()}
 
                },
 
                {},
 
                {{"HESS_SDF", common_data_ptr->hessSdfPtr()}}
 
                )
 
        );
      }
 
      return post_proc_fe;
    };
 
    auto get_post_proc_bdy_fe = [&]() {
      auto post_proc_fe =
          boost::make_shared<PostProcEleBdy>(*m_field_ptr, postProcMesh);
      auto &pip = post_proc_fe->getOpPtrVector();
 
      CHK_THROW_MESSAGE((AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
                            pip, {HDIV}, "GEOMETRY")),
                        "Apply transform");
      auto common_data_ptr = push_bdy_ops(pip);
 
      // create OP which run element on side
      auto op_loop_side = new OpLoopSide<SideEle>(
          *m_field_ptr, "dFE", SPACE_DIM, Sev::noisy,
          boost::make_shared<
              ForcesAndSourcesCore::UserDataOperator::AdjCache>());
      pip.push_back(op_loop_side);
 
      auto [hencky_or_hooke_common_data_ptr, contact_stress_ptr] =
          push_domain_ops(op_loop_side->getOpPtrVector());
 
      auto X_ptr = boost::make_shared<MatrixDouble>();
      pip.push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", X_ptr));
 
      pip.push_back(
 
          new OpPPMap(
 
              post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
 
              {{"SDF", common_data_ptr->sdfPtr()},
               {"CONSTRAINT_CONTACT", common_data_ptr->constraintPtr()}},
 
              {{"U", common_data_ptr->contactDispPtr()},
               {"GEOMETRY", X_ptr},
               {"TRACTION_CONTACT", common_data_ptr->contactTractionPtr()},
               {"GRAD_SDF", common_data_ptr->gradSdfPtr()}
 
              },
 
              {},
 
              {{"HESS_SDF", common_data_ptr->hessSdfPtr()}}
 
              )
 
      );
 
      return post_proc_fe;
    };
 
    auto get_integrate_traction = [&]() {
      auto integrate_traction = boost::make_shared<BoundaryEle>(*m_field_ptr);
      CHK_THROW_MESSAGE(
          (AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
              integrate_traction->getOpPtrVector(), {HDIV}, "GEOMETRY")),
          "Apply transform");
      // We have to integrate on curved face geometry, thus integration weight
      // have to adjusted.
      integrate_traction->getOpPtrVector().push_back(
          new OpSetHOWeightsOnSubDim<SPACE_DIM>());
      integrate_traction->getRuleHook = [](int, int, int approx_order) {
        return 2 * approx_order + geom_order - 1;
      };
 
      CHK_THROW_MESSAGE(
          (opFactoryCalculateTraction<SPACE_DIM, GAUSS, BoundaryEleOp>(
              integrate_traction->getOpPtrVector(), "SIGMA", is_axisymmetric)),
          "push operators to calculate traction");
 
      return integrate_traction;
    };
 
    auto get_integrate_area = [&]() {
      auto integrate_area = boost::make_shared<BoundaryEle>(*m_field_ptr);
 
      CHK_THROW_MESSAGE(
          (AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
              integrate_area->getOpPtrVector(), {HDIV}, "GEOMETRY")),
          "Apply transform");
      // We have to integrate on curved face geometry, thus integration weight
      // have to adjusted.
      integrate_area->getOpPtrVector().push_back(
          new OpSetHOWeightsOnSubDim<SPACE_DIM>());
      integrate_area->getRuleHook = [](int, int, int approx_order) {
        return 2 * approx_order + geom_order - 1;
      };
      Range contact_range;
      for (auto m :
           m_field_ptr->getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
               std::regex((boost::format("%s(.*)") % "CONTACT").str()))) {
        auto meshset = m->getMeshset();
        Range contact_meshset_range;
        CHKERR m_field_ptr->get_moab().get_entities_by_dimension(
            meshset, SPACE_DIM - 1, contact_meshset_range, true);
 
        CHKERR m_field_ptr->getInterface<CommInterface>()->synchroniseEntities(
            contact_meshset_range);
        contact_range.merge(contact_meshset_range);
      }
 
      auto contact_range_ptr = boost::make_shared<Range>(contact_range);
 
      auto op_loop_side = new OpLoopSide<SideEle>(
          *m_field_ptr, m_field_ptr->getInterface<Simple>()->getDomainFEName(),
          SPACE_DIM);
      CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
          op_loop_side->getOpPtrVector(), {H1}, "GEOMETRY");
 
      CHK_THROW_MESSAGE(
          (opFactoryCalculateArea<SPACE_DIM, GAUSS, BoundaryEleOp>(
              integrate_area->getOpPtrVector(), op_loop_side, "SIGMA", "U",
              is_axisymmetric, contact_range_ptr)),
          "push operators to calculate area");
 
      return integrate_area;
    };
 
    postProcDomainFe = get_post_proc_domain_fe();
    if constexpr (SPACE_DIM == 2)
      postProcBdyFe = get_post_proc_bdy_fe();
 
    integrateTraction = get_integrate_traction();
    integrateArea = get_integrate_area();
  }
 
  MoFEMErrorCode preProcess() { return 0; }
  MoFEMErrorCode operator()() { return 0; }
 
  MoFEMErrorCode postProcess() {
    MoFEMFunctionBegin;
    MoFEM::Interface *m_field_ptr;
    CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
 
    auto post_proc = [&]() {
      MoFEMFunctionBegin;
 
      if (!mfrontInterfacePtr) {
        auto post_proc_begin =
            boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(*m_field_ptr,
                                                                  postProcMesh);
        auto post_proc_end =
            boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(*m_field_ptr,
                                                                postProcMesh);
 
        CHKERR DMoFEMPreProcessFiniteElements(dM,
                                              post_proc_begin->getFEMethod());
        if (!postProcBdyFe) {
          postProcDomainFe->copyTs(*this); // this here is a Monitor
          CHKERR DMoFEMLoopFiniteElements(dM, "bFE", postProcDomainFe);
        } else {
          postProcDomainFe->copyTs(*this); // this here is a Monitor
          postProcBdyFe->copyTs(*this);
          CHKERR DMoFEMLoopFiniteElements(dM, "dFE", postProcDomainFe);
          CHKERR DMoFEMLoopFiniteElements(dM, "bFE", postProcBdyFe);
        }
        CHKERR DMoFEMPostProcessFiniteElements(dM,
                                               post_proc_end->getFEMethod());
 
        CHKERR post_proc_end->writeFile(
            "out_contact_" + boost::lexical_cast<std::string>(sTEP) + ".h5m");
      } else {
        CHKERR mfrontInterfacePtr->postProcess(
            ts_step, dM,
            m_field_ptr->getInterface<Simple>()->getDomainFEName());
      }
 
      MoFEMFunctionReturn(0);
    };
 
    auto calculate_force = [&] {
      MoFEMFunctionBegin;
      CHKERR VecZeroEntries(CommonData::totalTractionAndArea);
      CHKERR DMoFEMLoopFiniteElements(dM, "bFE", integrateTraction);
      CHKERR VecAssemblyBegin(CommonData::totalTractionAndArea);
      CHKERR VecAssemblyEnd(CommonData::totalTractionAndArea);
      MoFEMFunctionReturn(0);
    };
 
    auto calculate_area = [&] {
      MoFEMFunctionBegin;
      integrateArea->copyTs(*this);
      CHKERR DMoFEMLoopFiniteElements(dM, "bFE", integrateArea);
      CHKERR VecAssemblyBegin(CommonData::totalTractionAndArea);
      CHKERR VecAssemblyEnd(CommonData::totalTractionAndArea);
      MoFEMFunctionReturn(0);
    };
 
    auto calculate_reactions = [&]() {
      MoFEMFunctionBegin;
 
      auto res = createDMVector(dM);
 
      auto assemble_domain = [&]() {
        MoFEMFunctionBegin;
        auto fe_rhs = boost::make_shared<DomainEle>(*m_field_ptr);
        auto &pip = fe_rhs->getOpPtrVector();
        fe_rhs->f = res;
 
        auto integration_rule = [](int, int, int approx_order) {
          return 2 * approx_order + geom_order - 1;
        };
        fe_rhs->getRuleHook = integration_rule;
 
        CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1},
                                                              "GEOMETRY");
        CHKERR
        ContactOps::opFactoryDomainRhs<SPACE_DIM, PETSC, IT, DomainEleOp>(
            pip, "SIGMA", "U", is_axisymmetric);
 
        if (!mfrontInterfacePtr) {
          CHKERR
          HenckyOps::opFactoryDomainRhs<SPACE_DIM, PETSC, IT, DomainEleOp>(
              *m_field_ptr, pip, "U", "MAT_ELASTIC", Sev::inform, scale);
        } else {
          CHKERR mfrontInterfacePtr->opFactoryDomainRhs(pip, "U");
        }
        CHKERR DMoFEMLoopFiniteElements(dM, "dFE", fe_rhs);
 
        CHKERR VecAssemblyBegin(res);
        CHKERR VecAssemblyEnd(res);
        CHKERR VecGhostUpdateBegin(res, ADD_VALUES, SCATTER_REVERSE);
        CHKERR VecGhostUpdateEnd(res, ADD_VALUES, SCATTER_REVERSE);
 
        MoFEMFunctionReturn(0);
      };
 
      auto assemble_boundary = [&]() {
        MoFEMFunctionBegin;
        auto fe_rhs = boost::make_shared<BoundaryEle>(*m_field_ptr);
        auto &pip = fe_rhs->getOpPtrVector();
        fe_rhs->f = res;
 
        auto integration_rule = [](int, int, int approx_order) {
          return 2 * approx_order + geom_order - 1;
        };
        fe_rhs->getRuleHook = integration_rule;
 
        CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {},
                                                                  "GEOMETRY");
        // We have to integrate on curved face geometry, thus integration weight
        // have to adjusted.
        pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
 
        auto u_disp = boost::make_shared<MatrixDouble>();
        pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_disp));
        pip.push_back(
            new OpSpringRhs("U", u_disp, [this](double, double, double) {
              return spring_stiffness;
            }));
 
        CHKERR DMoFEMLoopFiniteElements(dM, "bFE", fe_rhs);
 
        MoFEMFunctionReturn(0);
      };
 
      CHKERR assemble_domain();
      CHKERR assemble_boundary();
 
      auto fe_post_proc_ptr = boost::make_shared<FEMethod>();
      auto get_post_proc_hook_rhs = [this, fe_post_proc_ptr, res,
                                     m_field_ptr]() {
        MoFEMFunctionBegin;
        CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
            *m_field_ptr, fe_post_proc_ptr, res)();
        MoFEMFunctionReturn(0);
      };
      fe_post_proc_ptr->postProcessHook = get_post_proc_hook_rhs;
      CHKERR DMoFEMPostProcessFiniteElements(dM, fe_post_proc_ptr.get());
 
      MoFEMFunctionReturn(0);
    };
 
    auto print_max_min = [&](auto &tuple, const std::string msg) {
      MoFEMFunctionBegin;
      CHKERR VecScatterBegin(std::get<1>(tuple), ts_u, std::get<0>(tuple),
                             INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecScatterEnd(std::get<1>(tuple), ts_u, std::get<0>(tuple),
                           INSERT_VALUES, SCATTER_FORWARD);
      double max, min;
      CHKERR VecMax(std::get<0>(tuple), PETSC_NULLPTR, &max);
      CHKERR VecMin(std::get<0>(tuple), PETSC_NULLPTR, &min);
      MOFEM_LOG_C("CONTACT", Sev::inform, "%s time %6.4e min %6.4e max %6.4e",
                  msg.c_str(), ts_t, min, max);
      MoFEMFunctionReturn(0);
    };
 
    auto print_force_and_area = [&]() {
      MoFEMFunctionBegin;
      MoFEM::Interface *m_field_ptr;
      CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
      if (!m_field_ptr->get_comm_rank()) {
        const double *t_ptr;
        CHKERR VecGetArrayRead(CommonData::totalTractionAndArea, &t_ptr);
        MOFEM_LOG_C("CONTACT", Sev::inform,
                    "Contact force: time %6.3e Fx: %6.6e Fy: %6.6e Fz: %6.6e",
                    ts_t, t_ptr[TX], t_ptr[TY], t_ptr[TZ]);
        MOFEM_LOG_C("CONTACT", Sev::inform,
                    "Contact area: time %6.3e Active: %6.6e Potential: %6.6e",
                    ts_t, t_ptr[AREA_ACTIVE], t_ptr[AREA_POTENTIAL]);
        CHKERR VecRestoreArrayRead(CommonData::totalTractionAndArea, &t_ptr);
      }
      MoFEMFunctionReturn(0);
    };
 
    if (mfrontInterfacePtr) {
      CHKERR mfrontInterfacePtr->updateInternalVariables(
          dM, m_field_ptr->getInterface<Simple>()->getDomainFEName());
    }
 
    int se = 1;
    CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-save_every", &se,
                              PETSC_NULLPTR);
 
    if (!(ts_step % se)) {
      MOFEM_LOG("CONTACT", Sev::inform)
          << "Write file at time " << ts_t << " write step " << sTEP;
      CHKERR post_proc();
    }
    CHKERR calculate_force();
    CHKERR calculate_area();
 
    CHKERR calculate_reactions();
 
    CHKERR print_max_min(uXScatter, "Ux");
    CHKERR print_max_min(uYScatter, "Uy");
    if (SPACE_DIM == 3)
      CHKERR print_max_min(uZScatter, "Uz");
    CHKERR print_force_and_area();
    ++sTEP;
 
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode setScatterVectors(
      std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> ux_scatter,
      std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uy_scatter,
      std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uz_scatter) {
    MoFEMFunctionBegin;
    uXScatter = ux_scatter;
    uYScatter = uy_scatter;
    uZScatter = uz_scatter;
    MoFEMFunctionReturn(0);
  }
 
private:
  SmartPetscObj<DM> dM;
  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
 
  boost::shared_ptr<moab::Core> postProcMesh = boost::make_shared<moab::Core>();
 
  boost::shared_ptr<PostProcEleDomain> postProcDomainFe;
  boost::shared_ptr<PostProcEleBdy> postProcBdyFe;
 
  boost::shared_ptr<BoundaryEle> integrateTraction;
  boost::shared_ptr<BoundaryEle> integrateArea;
 
  moab::Core mbVertexPostproc;
  moab::Interface &moabVertex;
 
  double lastTime;
  double deltaTime;
  int sTEP;
  bool isLargeStrain;
 
  boost::shared_ptr<MFrontInterface> mfrontInterfacePtr;
};
 
} // namespace ContactOps