ADV-1: Contact problem

Note

Prerequisites of this tutorial include VEC-0: Linear elasticity

Note

Intended learning outcome:

• general structure of a program developed using MoFEM
• idea of Simple Interface in MoFEM and how to use it
• idea of Domain element in MoFEM and how to use it
• Use default Forms Integrals
• how to push the developed UDOs to the Pipeline
• utilisation of tools to convert outputs (MOAB) and visualise them (Paraview)

/**
 * \file contact.cpp
 * \CONTACT contact.cpp
 *
 * CONTACT of contact problem
 *
 * @copyright Copyright (c) 2023
 */
 
/* The above code is a preprocessor directive in C++ that checks if the macro
"EXECUTABLE_DIMENSION" has been defined. If it has not been defined, it is set
to 3" */
#ifndef EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION 3
#endif
 
#ifndef SCHUR_ASSEMBLE
#define SCHUR_ASSEMBLE 0
#endif
 
#include <MoFEM.hpp>
#include <MatrixFunction.hpp>
 
#ifdef PYTHON_SDF
#include <boost/python.hpp>
#include <boost/python/def.hpp>
#include <boost/python/numpy.hpp>
namespace bp = boost::python;
namespace np = boost::python::numpy;
#endif
 
using namespace MoFEM;
 
constexpr AssemblyType AT =
    (SCHUR_ASSEMBLE) ? AssemblyType::SCHUR
                     : AssemblyType::PETSC; //< selected assembly type
constexpr IntegrationType IT =
    IntegrationType::GAUSS; //< selected integration type
 
template <int DIM> struct ElementsAndOps;
 
template <> struct ElementsAndOps<2> : PipelineManager::ElementsAndOpsByDim<2> {
  static constexpr FieldSpace CONTACT_SPACE = HCURL;
};
 
template <> struct ElementsAndOps<3> : PipelineManager::ElementsAndOpsByDim<3> {
  static constexpr FieldSpace CONTACT_SPACE = HDIV;
};
 
constexpr FieldSpace ElementsAndOps<2>::CONTACT_SPACE;
constexpr FieldSpace ElementsAndOps<3>::CONTACT_SPACE;
 
constexpr int SPACE_DIM =
    EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
 
/* The above code is defining an alias `EntData` for the type
`EntitiesFieldData::EntData`. This is a C++ syntax for creating a new name for
an existing type. */
using EntData = EntitiesFieldData::EntData;
using DomainEle = ElementsAndOps<SPACE_DIM>::DomainEle;
using DomainEleOp = DomainEle::UserDataOperator;
using BoundaryEle = ElementsAndOps<SPACE_DIM>::BoundaryEle;
using BoundaryEleOp = BoundaryEle::UserDataOperator;
 
//! [Specialisation for assembly]
 
// Assemble to A matrix, by default, however, some terms are assembled only to
// preconditioning. 
 
template <>
typename MoFEM::OpBaseImpl<SCHUR, DomainEleOp>::MatSetValuesHook
    MoFEM::OpBaseImpl<SCHUR, DomainEleOp>::matSetValuesHook =
        [](ForcesAndSourcesCore::UserDataOperator *op_ptr,
           const EntitiesFieldData::EntData &row_data,
           const EntitiesFieldData::EntData &col_data, MatrixDouble &m) {
          return MatSetValues<AssemblyTypeSelector<SCHUR>>(
              op_ptr->getKSPA(), row_data, col_data, m, ADD_VALUES);
        };
 
template <>
typename MoFEM::OpBaseImpl<SCHUR, BoundaryEleOp>::MatSetValuesHook
    MoFEM::OpBaseImpl<SCHUR, BoundaryEleOp>::matSetValuesHook =
        [](ForcesAndSourcesCore::UserDataOperator *op_ptr,
           const EntitiesFieldData::EntData &row_data,
           const EntitiesFieldData::EntData &col_data, MatrixDouble &m) {
          return MatSetValues<AssemblyTypeSelector<SCHUR>>(
              op_ptr->getKSPA(), row_data, col_data, m, ADD_VALUES);
        };
 
/**
 * @brief Element used to specialise assembly
 * 
 */
struct BoundaryEleOpStab : public BoundaryEleOp {
  using BoundaryEleOp::BoundaryEleOp;
};
 
/**
 * @brief Specialise assembly for Stabilised matrix
 *
 * @tparam
 */
template <>
typename MoFEM::OpBaseImpl<SCHUR, BoundaryEleOpStab>::MatSetValuesHook
    MoFEM::OpBaseImpl<SCHUR, BoundaryEleOpStab>::matSetValuesHook =
        [](ForcesAndSourcesCore::UserDataOperator *op_ptr,
           const EntitiesFieldData::EntData &row_data,
           const EntitiesFieldData::EntData &col_data, MatrixDouble &m) {
          return MatSetValues<AssemblyTypeSelector<SCHUR>>(
              op_ptr->getKSPB(), row_data, col_data, m, ADD_VALUES);
        };
//! [Specialisation for assembly]
 
constexpr FieldSpace CONTACT_SPACE = ElementsAndOps<SPACE_DIM>::CONTACT_SPACE;
 
//! [Operators used for contact]
using OpSpringLhs = FormsIntegrators<BoundaryEleOp>::Assembly<AT>::BiLinearForm<
    IT>::OpMass<1, SPACE_DIM>;
using OpSpringRhs = FormsIntegrators<BoundaryEleOp>::Assembly<AT>::LinearForm<
    IT>::OpBaseTimesVector<1, SPACE_DIM, 1>;
//! [Operators used for contact]
 
PetscBool is_quasi_static = PETSC_TRUE;
 
int order = 2;
int contact_order = 2;
int geom_order = 1;
double young_modulus = 100;
double poisson_ratio = 0.25;
double rho = 0.0;
double spring_stiffness = 0.0;
double vis_spring_stiffness = 0.0;
double alpha_damping = 0;
 
double scale = 1.;
 
PetscBool is_axisymmetric = PETSC_FALSE; //< Axisymmetric model
 
int atom_test = 0;
 
namespace ContactOps {
double cn_contact = 0.1;
}; // namespace ContactOps
 
//#define HENCKY_SMALL_STRAIN
 
#include <HenckyOps.hpp>
using namespace HenckyOps;
#include <ContactOps.hpp>
#include <PostProcContact.hpp>
#include <ContactNaturalBC.hpp>
 
#ifdef WITH_MODULE_MFRONT_INTERFACE
#include <MFrontMoFEMInterface.hpp>
#endif
 
using DomainRhsBCs = NaturalBC<DomainEleOp>::Assembly<AT>::LinearForm<IT>;
using OpDomainRhsBCs =
    DomainRhsBCs::OpFlux<ContactOps::DomainBCs, 1, SPACE_DIM>;
using BoundaryRhsBCs = NaturalBC<BoundaryEleOp>::Assembly<AT>::LinearForm<IT>;
using OpBoundaryRhsBCs =
    BoundaryRhsBCs::OpFlux<ContactOps::BoundaryBCs, 1, SPACE_DIM>;
using BoundaryLhsBCs = NaturalBC<BoundaryEleOp>::Assembly<AT>::BiLinearForm<IT>;
using OpBoundaryLhsBCs =
    BoundaryLhsBCs::OpFlux<ContactOps::BoundaryBCs, 1, SPACE_DIM>;
 
using namespace ContactOps;
 
struct Contact {
 
  Contact(MoFEM::Interface &m_field) : mField(m_field) {}
 
  MoFEMErrorCode runProblem();
 
private:
  MoFEM::Interface &mField;
 
  MoFEMErrorCode setupProblem();
  MoFEMErrorCode createCommonData();
  MoFEMErrorCode bC();
  MoFEMErrorCode OPs();
  MoFEMErrorCode tsSolve();
  MoFEMErrorCode checkResults();
 
  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
 
  boost::shared_ptr<GenericElementInterface> mfrontInterface;
  boost::shared_ptr<Monitor> monitorPtr;
 
#ifdef PYTHON_SDF
  boost::shared_ptr<SDFPython> sdfPythonPtr;
#endif
 
  struct ScaledTimeScale : public MoFEM::TimeScale {
    using MoFEM::TimeScale::TimeScale;
    double getScale(const double time) {
      return scale * MoFEM::TimeScale::getScale(time);
    };
  };
};
 
//! [Run problem]
MoFEMErrorCode Contact::runProblem() {
  MoFEMFunctionBegin;
  CHKERR setupProblem();
  CHKERR createCommonData();
  CHKERR bC();
  CHKERR OPs();
  CHKERR tsSolve();
  CHKERR checkResults();
  MoFEMFunctionReturn(0);
}
//! [Run problem]
 
//! [Set up problem]
MoFEMErrorCode Contact::setupProblem() {
  MoFEMFunctionBegin;
  Simple *simple = mField.getInterface<Simple>();
  PetscBool use_mfront = PETSC_FALSE;
 
  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-contact_order", &contact_order,
                            PETSC_NULL);
  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-geom_order", &geom_order,
                            PETSC_NULL);
 
  MOFEM_LOG("CONTACT", Sev::inform) << "Order " << order;
  if (contact_order != order)
    MOFEM_LOG("CONTACT", Sev::inform) << "Contact order " << contact_order;
  MOFEM_LOG("CONTACT", Sev::inform) << "Geom order " << geom_order;
 
  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-use_mfront", &use_mfront,
                             PETSC_NULL);
  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_axisymmetric",
                             &is_axisymmetric, PETSC_NULL);
  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-atom_test", &atom_test, PETSC_NULL);
 
  // Select base
  enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
  const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
  PetscInt choice_base_value = AINSWORTH;
  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
                              LASBASETOPT, &choice_base_value, PETSC_NULL);
 
  FieldApproximationBase base;
  switch (choice_base_value) {
  case AINSWORTH:
    base = AINSWORTH_LEGENDRE_BASE;
    MOFEM_LOG("CONTACT", Sev::inform)
        << "Set AINSWORTH_LEGENDRE_BASE for displacements";
    break;
  case DEMKOWICZ:
    base = DEMKOWICZ_JACOBI_BASE;
    MOFEM_LOG("CONTACT", Sev::inform)
        << "Set DEMKOWICZ_JACOBI_BASE for displacements";
    break;
  default:
    base = LASTBASE;
    break;
  }
 
  // Note: For tets we have only H1 Ainsworth base, for Hex we have only H1
  // Demkowicz base. We need to implement Demkowicz H1 base on tet.
  CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
  CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
  CHKERR simple->addDomainField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
                                SPACE_DIM);
  CHKERR simple->addBoundaryField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
                                  SPACE_DIM);
  CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
 
  CHKERR simple->setFieldOrder("U", order);
  CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
 
  auto get_skin = [&]() {
    Range body_ents;
    CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
    Skinner skin(&mField.get_moab());
    Range skin_ents;
    CHKERR skin.find_skin(0, body_ents, false, skin_ents);
    return skin_ents;
  };
 
  auto filter_blocks = [&](auto skin) {
    bool is_contact_block = false;
    Range contact_range;
    for (auto m :
         mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
 
             (boost::format("%s(.*)") % "CONTACT").str()
 
                 ))
 
    ) {
      is_contact_block =
          true; ///< blocs interation is collective, so that is set irrespective
                ///< if there are entities in given rank or not in the block
      MOFEM_LOG("CONTACT", Sev::inform)
          << "Find contact block set:  " << m->getName();
      auto meshset = m->getMeshset();
      Range contact_meshset_range;
      CHKERR mField.get_moab().get_entities_by_dimension(
          meshset, SPACE_DIM - 1, contact_meshset_range, true);
 
      CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
          contact_meshset_range);
      contact_range.merge(contact_meshset_range);
    }
    if (is_contact_block) {
      MOFEM_LOG("SYNC", Sev::inform)
          << "Nb entities in contact surface: " << contact_range.size();
      MOFEM_LOG_SYNCHRONISE(mField.get_comm());
      skin = intersect(skin, contact_range);
    }
    return skin;
  };
 
  auto filter_true_skin = [&](auto skin) {
    Range boundary_ents;
    ParallelComm *pcomm =
        ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
    CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
                                 PSTATUS_NOT, -1, &boundary_ents);
    return boundary_ents;
  };
 
  auto boundary_ents = filter_true_skin(filter_blocks(get_skin()));
  CHKERR simple->setFieldOrder("SIGMA", 0);
  int sigma_order = std::max(order, contact_order) - 1;
  CHKERR simple->setFieldOrder("SIGMA", sigma_order, &boundary_ents);
 
  if (contact_order > order) {
    Range ho_ents;
    if constexpr (SPACE_DIM == 3) {
      CHKERR mField.get_moab().get_adjacencies(boundary_ents, 1, false, ho_ents,
                                               moab::Interface::UNION);
    } else {
      ho_ents = boundary_ents;
    }
    CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(ho_ents);
    CHKERR simple->setFieldOrder("U", contact_order, &ho_ents);
    CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities("U");
  }
 
  if (is_axisymmetric) {
    if (SPACE_DIM == 3) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Use executable contact_2d with axisymmetric model");
    } else {
      if (!use_mfront) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
                "Axisymmetric model is only available with MFront (set "
                "use_mfront to 1)");
      } else {
        MOFEM_LOG("CONTACT", Sev::inform) << "Using axisymmetric model";
      }
    }
  } else {
    if (SPACE_DIM == 2) {
      MOFEM_LOG("CONTACT", Sev::inform) << "Using plane strain model";
    }
  }
 
  if (!use_mfront) {
    CHKERR simple->setUp();
  } else {
#ifndef WITH_MODULE_MFRONT_INTERFACE
    SETERRQ(
        PETSC_COMM_SELF, MOFEM_NOT_FOUND,
        "MFrontInterface module was not found while use_mfront was set to 1");
#else
    if (SCHUR_ASSEMBLE) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
              "MFrontInterface module is not compatible with Schur assembly");
    }
    if (SPACE_DIM == 3) {
      mfrontInterface =
          boost::make_shared<MFrontMoFEMInterface<TRIDIMENSIONAL>>(
              mField, "U", "GEOMETRY", true, is_quasi_static);
    } else if (SPACE_DIM == 2) {
      if (is_axisymmetric) {
        mfrontInterface =
            boost::make_shared<MFrontMoFEMInterface<AXISYMMETRICAL>>(
                mField, "U", "GEOMETRY", true, is_quasi_static);
      } else {
        mfrontInterface = boost::make_shared<MFrontMoFEMInterface<PLANESTRAIN>>(
            mField, "U", "GEOMETRY", true, is_quasi_static);
      }
    }
#endif
 
    CHKERR mfrontInterface->getCommandLineParameters();
    CHKERR mfrontInterface->addElementFields();
    CHKERR mfrontInterface->createElements();
 
    CHKERR simple->defineFiniteElements();
    CHKERR simple->defineProblem(PETSC_TRUE);
    CHKERR simple->buildFields();
    CHKERR simple->buildFiniteElements();
 
    CHKERR mField.build_finite_elements("MFRONT_EL");
    CHKERR mfrontInterface->addElementsToDM(simple->getDM());
 
    CHKERR simple->buildProblem();
  }
 
  auto dm = simple->getDM();
  monitorPtr =
      boost::make_shared<Monitor>(dm, mfrontInterface, is_axisymmetric);
  if (use_mfront) {
    mfrontInterface->setMonitorPtr(monitorPtr);
  }
 
  auto project_ho_geometry = [&]() {
    Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
    return mField.loop_dofs("GEOMETRY", ent_method);
  };
  
  PetscBool project_geometry = PETSC_TRUE;
  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-project_geometry",
                               &project_geometry, PETSC_NULL);
  if (project_geometry){
    CHKERR project_ho_geometry();
  }
 
  MoFEMFunctionReturn(0);
} //! [Set up problem]
 
//! [Create common data]
MoFEMErrorCode Contact::createCommonData() {
  MoFEMFunctionBegin;
 
  auto get_options = [&]() {
    MoFEMFunctionBegin;
    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-scale", &scale, PETSC_NULL);
    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-young_modulus",
                                 &young_modulus, PETSC_NULL);
    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-poisson_ratio",
                                 &poisson_ratio, PETSC_NULL);
    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-rho", &rho, PETSC_NULL);
    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn", &cn_contact,
                                 PETSC_NULL);
    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-spring_stiffness",
                                 &spring_stiffness, PETSC_NULL);
    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-vis_spring_stiffness",
                                 &vis_spring_stiffness, PETSC_NULL);
    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-alpha_damping",
                                 &alpha_damping, PETSC_NULL);
                               
    if (!mfrontInterface) {
      MOFEM_LOG("CONTACT", Sev::inform) << "Young modulus " << young_modulus;
      MOFEM_LOG("CONTACT", Sev::inform) << "Poisson_ratio " << poisson_ratio;
    } else {
      MOFEM_LOG("CONTACT", Sev::inform) << "Using MFront for material model";
    }
 
    MOFEM_LOG("CONTACT", Sev::inform) << "Density " << rho;
    MOFEM_LOG("CONTACT", Sev::inform) << "cn_contact " << cn_contact;
    MOFEM_LOG("CONTACT", Sev::inform)
        << "Spring stiffness " << spring_stiffness;
    MOFEM_LOG("CONTACT", Sev::inform)
        << "Vis spring_stiffness " << vis_spring_stiffness;
 
    MOFEM_LOG("CONTACT", Sev::inform) << "alpha_damping " << alpha_damping;
 
    PetscBool use_scale = PETSC_FALSE;
    CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-use_scale", &use_scale,
                               PETSC_NULL);
    if (use_scale) {
      scale /= young_modulus;
    }
 
    MOFEM_LOG("CONTACT", Sev::inform) << "Scale " << scale;
 
    CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_quasi_static",
                               &is_quasi_static, PETSC_NULL);
    MOFEM_LOG("CONTACT", Sev::inform)
        << "Is quasi-static: " << (is_quasi_static ? "true" : "false");
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR get_options();
 
#ifdef PYTHON_SDF
  char sdf_file_name[255];
  CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-sdf_file",
                               sdf_file_name, 255, PETSC_NULL);
 
  sdfPythonPtr = boost::make_shared<SDFPython>();
  CHKERR sdfPythonPtr->sdfInit(sdf_file_name);
  sdfPythonWeakPtr = sdfPythonPtr;
#endif
 
  MoFEMFunctionReturn(0);
}
//! [Create common data]
 
//! [Boundary condition]
MoFEMErrorCode Contact::bC() {
  MoFEMFunctionBegin;
  auto bc_mng = mField.getInterface<BcManager>();
  auto simple = mField.getInterface<Simple>();
 
  for (auto f : {"U", "SIGMA"}) {
    CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
                                             "REMOVE_X", f, 0, 0);
    CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
                                             "REMOVE_Y", f, 1, 1);
    CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
                                             "REMOVE_Z", f, 2, 2);
    CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
                                             "REMOVE_ALL", f, 0, 3);
  }
 
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_X",
                                           "SIGMA", 0, 0, false, true);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Y",
                                           "SIGMA", 1, 1, false, true);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Z",
                                           "SIGMA", 2, 2, false, true);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_ALL",
                                           "SIGMA", 0, 3, false, true);
  CHKERR bc_mng->removeBlockDOFsOnEntities(
      simple->getProblemName(), "NO_CONTACT", "SIGMA", 0, 3, false, true);
 
  // Note remove has to be always before push. Then node marking will be
  // corrupted.
  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
      simple->getProblemName(), "U");
 
  MoFEMFunctionReturn(0);
}
//! [Boundary condition]
 
//! [Push operators to pip]
MoFEMErrorCode Contact::OPs() {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  auto *pip_mng = mField.getInterface<PipelineManager>();
  auto bc_mng = mField.getInterface<BcManager>();
  auto time_scale = boost::make_shared<ScaledTimeScale>();
  auto body_force_time_scale =
      boost::make_shared<ScaledTimeScale>("body_force_hist.txt");
 
  auto integration_rule_vol = [](int, int, int approx_order) {
    return 2 * approx_order + geom_order - 1;
  };
  auto integration_rule_boundary = [](int, int, int approx_order) {
    return 2 * approx_order + geom_order - 1;
  };
 
  auto add_domain_base_ops = [&](auto &pip) {
    MoFEMFunctionBegin;
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, HDIV},
                                                          "GEOMETRY");
    MoFEMFunctionReturn(0);
  };
 
  auto henky_common_data_ptr = boost::make_shared<HenckyOps::CommonData>();
  henky_common_data_ptr->matDPtr = boost::make_shared<MatrixDouble>();
  henky_common_data_ptr->matGradPtr = boost::make_shared<MatrixDouble>();
 
  auto add_domain_ops_lhs = [&](auto &pip) {
    MoFEMFunctionBegin;
 
    //! [Only used for dynamics]
    using OpMass = FormsIntegrators<DomainEleOp>::Assembly<AT>::BiLinearForm<
        GAUSS>::OpMass<1, SPACE_DIM>;
    //! [Only used for dynamics]
    if (is_quasi_static == PETSC_FALSE) {
 
      auto *pip_mng = mField.getInterface<PipelineManager>();
      auto fe_domain_lhs = pip_mng->getDomainLhsFE();
 
      auto get_inertia_and_mass_damping =
          [this, fe_domain_lhs](const double, const double, const double) {
            return (rho * scale) * fe_domain_lhs->ts_aa +
                   (alpha_damping * scale) * fe_domain_lhs->ts_a;
          };
      pip.push_back(new OpMass("U", "U", get_inertia_and_mass_damping));
    } else {
 
      auto *pip_mng = mField.getInterface<PipelineManager>();
      auto fe_domain_lhs = pip_mng->getDomainLhsFE();
 
      auto get_inertia_and_mass_damping =
          [this, fe_domain_lhs](const double, const double, const double) {
            return (alpha_damping * scale) * fe_domain_lhs->ts_a;
          };
      pip.push_back(new OpMass("U", "U", get_inertia_and_mass_damping));
 
    }
 
    if (!mfrontInterface) {
      CHKERR HenckyOps::opFactoryDomainLhs<SPACE_DIM, AT, IT, DomainEleOp>(
          mField, pip, "U", "MAT_ELASTIC", Sev::verbose, scale);
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto add_domain_ops_rhs = [&](auto &pip) {
    MoFEMFunctionBegin;
 
    CHKERR DomainRhsBCs::AddFluxToPipeline<OpDomainRhsBCs>::add(
        pip, mField, "U", {body_force_time_scale}, Sev::inform);
 
    //! [Only used for dynamics]
    using OpInertiaForce = FormsIntegrators<DomainEleOp>::Assembly<
        AT>::LinearForm<IT>::OpBaseTimesVector<1, SPACE_DIM, 1>;
    //! [Only used for dynamics]
 
    // only in case of dynamics
    if (is_quasi_static == PETSC_FALSE) {
      auto mat_acceleration = boost::make_shared<MatrixDouble>();
      pip.push_back(new OpCalculateVectorFieldValuesDotDot<SPACE_DIM>(
          "U", mat_acceleration));
      pip.push_back(
          new OpInertiaForce("U", mat_acceleration, [](double, double, double) {
            return rho * scale;
          }));
 
    }
 
    // only in case of viscosity
    if (alpha_damping > 0) {
      auto mat_velocity = boost::make_shared<MatrixDouble>();
      pip.push_back(
          new OpCalculateVectorFieldValuesDot<SPACE_DIM>("U", mat_velocity));
      pip.push_back(
          new OpInertiaForce("U", mat_velocity, [](double, double, double) {
            return alpha_damping * scale;
          }));
    }
 
    if (!mfrontInterface) {
      CHKERR HenckyOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
          mField, pip, "U", "MAT_ELASTIC", Sev::inform, scale);
    }
 
    CHKERR ContactOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
        pip, "SIGMA", "U", is_axisymmetric);
 
    MoFEMFunctionReturn(0);
  };
 
  auto add_boundary_base_ops = [&](auto &pip) {
    MoFEMFunctionBegin;
    CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {HDIV},
                                                              "GEOMETRY");
    // We have to integrate on curved face geometry, thus integration weight
    // have to adjusted.
    pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
    MoFEMFunctionReturn(0);
  };
 
  auto add_boundary_ops_lhs = [&](auto &pip) {
    MoFEMFunctionBegin;
 
    //! [Operators used for contact]
    using OpSpringLhs = FormsIntegrators<BoundaryEleOp>::Assembly<
        AT>::BiLinearForm<IT>::OpMass<1, SPACE_DIM>;
    //! [Operators used for contact]
 
    // Add Natural BCs to LHS
    CHKERR BoundaryLhsBCs::AddFluxToPipeline<OpBoundaryLhsBCs>::add(
        pip, mField, "U", Sev::inform);
 
    if (spring_stiffness > 0 || vis_spring_stiffness > 0) {
 
      auto *pip_mng = mField.getInterface<PipelineManager>();
      auto fe_boundary_lhs = pip_mng->getBoundaryLhsFE();
 
      pip.push_back(new OpSpringLhs(
          "U", "U",
 
          [this, fe_boundary_lhs](double, double, double) {
            return spring_stiffness * scale +
                   (vis_spring_stiffness * scale) * fe_boundary_lhs->ts_a;
          }
 
          ));
    }
 
    CHKERR
    ContactOps::opFactoryBoundaryLhs<SPACE_DIM, AT, GAUSS, BoundaryEleOp>(
        pip, "SIGMA", "U", is_axisymmetric);
    CHKERR ContactOps::opFactoryBoundaryToDomainLhs<SPACE_DIM, AT, GAUSS,
                                                    DomainEle>(
        mField, pip, simple->getDomainFEName(), "SIGMA", "U", "GEOMETRY",
        integration_rule_vol, is_axisymmetric);
 
    MoFEMFunctionReturn(0);
  };
 
  auto add_boundary_ops_rhs = [&](auto &pip) {
    MoFEMFunctionBegin;
 
    //! [Operators used for contact]
    using OpSpringRhs = FormsIntegrators<BoundaryEleOp>::Assembly<
        AT>::LinearForm<IT>::OpBaseTimesVector<1, SPACE_DIM, 1>;
    //! [Operators used for contact]
    
    // Add Natural BCs to RHS
    CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(
        pip, mField, "U", {time_scale}, Sev::inform);
 
    if (spring_stiffness > 0 || vis_spring_stiffness > 0) {
      auto u_disp = boost::make_shared<MatrixDouble>();
      auto dot_u_disp = boost::make_shared<MatrixDouble>();
      pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_disp));
      pip.push_back(
          new OpCalculateVectorFieldValuesDot<SPACE_DIM>("U", dot_u_disp));
      pip.push_back(
          new OpSpringRhs("U", u_disp, [this](double, double, double) {
            return spring_stiffness * scale;
          }));
      pip.push_back(
          new OpSpringRhs("U", dot_u_disp, [this](double, double, double) {
            return vis_spring_stiffness * scale;
          }));
    }
 
    CHKERR
    ContactOps::opFactoryBoundaryRhs<SPACE_DIM, AT, GAUSS, BoundaryEleOp>(
        pip, "SIGMA", "U", is_axisymmetric);
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR add_domain_base_ops(pip_mng->getOpDomainLhsPipeline());
  CHKERR add_domain_base_ops(pip_mng->getOpDomainRhsPipeline());
  CHKERR add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
  CHKERR add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
 
  CHKERR add_boundary_base_ops(pip_mng->getOpBoundaryLhsPipeline());
  CHKERR add_boundary_base_ops(pip_mng->getOpBoundaryRhsPipeline());
  CHKERR add_boundary_ops_lhs(pip_mng->getOpBoundaryLhsPipeline());
  CHKERR add_boundary_ops_rhs(pip_mng->getOpBoundaryRhsPipeline());
 
  if (mfrontInterface) {
    auto t_type = GenericElementInterface::IM2;
    if (is_quasi_static == PETSC_TRUE)
      t_type = GenericElementInterface::IM;
 
    CHKERR mfrontInterface->setOperators();
    CHKERR mfrontInterface->setupSolverFunctionTS(t_type);
    CHKERR mfrontInterface->setupSolverJacobianTS(t_type);
  }
 
  CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule_vol);
  CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule_vol);
  CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule_boundary);
  CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule_boundary);
 
  MoFEMFunctionReturn(0);
}
//! [Push operators to pip]
 
//! [Solve]
struct SetUpSchur {
  static boost::shared_ptr<SetUpSchur>
  createSetUpSchur(MoFEM::Interface &m_field);
  virtual MoFEMErrorCode setUp(SmartPetscObj<TS> solver) = 0;
 
protected:
  SetUpSchur() = default;
};
 
MoFEMErrorCode Contact::tsSolve() {
  MoFEMFunctionBegin;
 
  Simple *simple = mField.getInterface<Simple>();
  PipelineManager *pip_mng = mField.getInterface<PipelineManager>();
  ISManager *is_manager = mField.getInterface<ISManager>();
 
  auto set_section_monitor = [&](auto solver) {
    MoFEMFunctionBegin;
    SNES snes;
    CHKERR TSGetSNES(solver, &snes);
    PetscViewerAndFormat *vf;
    CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
                                      PETSC_VIEWER_DEFAULT, &vf);
    CHKERR SNESMonitorSet(
        snes,
        (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
        vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
    MoFEMFunctionReturn(0);
  };
 
  auto scatter_create = [&](auto D, auto coeff) {
    SmartPetscObj<IS> is;
    CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
                                                   ROW, "U", coeff, coeff, is);
    int loc_size;
    CHKERR ISGetLocalSize(is, &loc_size);
    Vec v;
    CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
    VecScatter scatter;
    CHKERR VecScatterCreate(D, is, v, PETSC_NULL, &scatter);
    return std::make_tuple(SmartPetscObj<Vec>(v),
                           SmartPetscObj<VecScatter>(scatter));
  };
 
  auto set_time_monitor = [&](auto dm, auto solver) {
    MoFEMFunctionBegin;
    monitorPtr->setScatterVectors(uXScatter, uYScatter, uZScatter);
    boost::shared_ptr<ForcesAndSourcesCore> null;
    CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
                               monitorPtr, null, null);
    MoFEMFunctionReturn(0);
  };
 
  auto set_essential_bc = [&]() {
    MoFEMFunctionBegin;
    // This is low level pushing finite elements (pipelines) to solver
    auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
    auto pre_proc_ptr = boost::make_shared<FEMethod>();
    auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
    auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
 
    // Add boundary condition scaling
    auto time_scale = boost::make_shared<TimeScale>();
 
    auto get_bc_hook_rhs = [&]() {
      EssentialPreProc<DisplacementCubitBcData> hook(mField, pre_proc_ptr,
                                                     {time_scale}, false);
      return hook;
    };
    pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
 
    auto get_post_proc_hook_rhs = [&]() {
      return EssentialPostProcRhs<DisplacementCubitBcData>(
          mField, post_proc_rhs_ptr, 1.);
    };
    auto get_post_proc_hook_lhs = [&]() {
      return EssentialPostProcLhs<DisplacementCubitBcData>(
          mField, post_proc_lhs_ptr, 1.);
    };
    post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs();
 
    ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
    ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
    ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
    if (AT != AssemblyType::SCHUR) {
      post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
      ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
    }
    MoFEMFunctionReturn(0);
  };
 
  auto set_schur_pc = [&](auto solver) {
    boost::shared_ptr<SetUpSchur> schur_ptr;
    if (AT == AssemblyType::SCHUR) {
      schur_ptr = SetUpSchur::createSetUpSchur(mField);
      CHK_MOAB_THROW(schur_ptr->setUp(solver), "SetUpSchur::setUp");
    }
    return schur_ptr;
  };
 
  auto dm = simple->getDM();
  auto D = createDMVector(dm);
 
  ContactOps::CommonData::createTotalTraction(mField);
 
  uXScatter = scatter_create(D, 0);
  uYScatter = scatter_create(D, 1);
  if (SPACE_DIM == 3)
    uZScatter = scatter_create(D, 2);
 
  // Add extra finite elements to SNES solver pipelines to resolve essential
  // boundary conditions
  CHKERR set_essential_bc();
 
  if (is_quasi_static == PETSC_TRUE) {
    auto solver = pip_mng->createTSIM();
    CHKERR TSSetFromOptions(solver);
 
    auto D = createDMVector(dm);
    auto schur_pc_ptr = set_schur_pc(solver);
 
    CHKERR set_section_monitor(solver);
    CHKERR set_time_monitor(dm, solver);
    CHKERR TSSetSolution(solver, D);
    CHKERR TSSetUp(solver);
    CHKERR TSSolve(solver, NULL);
  } else {
    auto solver = pip_mng->createTSIM2();
    CHKERR TSSetFromOptions(solver);
 
    auto dm = simple->getDM();
    auto D = createDMVector(dm);
    auto DD = vectorDuplicate(D);
    auto schur_pc_ptr = set_schur_pc(solver);
 
    CHKERR set_section_monitor(solver);
    CHKERR set_time_monitor(dm, solver);
    CHKERR TS2SetSolution(solver, D, DD);
    CHKERR TSSetUp(solver);
    CHKERR TSSolve(solver, NULL);
  }
 
  MoFEMFunctionReturn(0);
}
//! [Solve]
 
//! [Check]
MoFEMErrorCode Contact::checkResults() {
  MoFEMFunctionBegin;
  if (atom_test == 1 && !mField.get_comm_rank()) {
    const double *t_ptr;
    CHKERR VecGetArrayRead(ContactOps::CommonData::totalTraction, &t_ptr);
    double hertz_tract = 158.73;
    double tol = 4e-3;
    if (fabs(t_ptr[1] - hertz_tract) / hertz_tract > tol) {
      SETERRQ3(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
               "atom test %d diverged! %3.4e != %3.4e", atom_test, t_ptr[1],
               hertz_tract);
    }
    CHKERR VecRestoreArrayRead(ContactOps::CommonData::totalTraction, &t_ptr);
  }
 
  ContactOps::CommonData::totalTraction.reset();
 
  MoFEMFunctionReturn(0);
}
//! [Check]
 
static char help[] = "...\n\n";
 
int main(int argc, char *argv[]) {
 
#ifdef PYTHON_SDF
  Py_Initialize();
  np::initialize();
#endif
 
  // Initialisation of MoFEM/PETSc and MOAB data structures
  const char param_file[] = "param_file.petsc";
  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
 
  // Add logging channel for CONTACT
  auto core_log = logging::core::get();
  core_log->add_sink(
      LogManager::createSink(LogManager::getStrmWorld(), "CONTACT"));
  LogManager::setLog("CONTACT");
  MOFEM_LOG_TAG("CONTACT", "Indent");
 
  try {
 
    //! [Register MoFEM discrete manager in PETSc]
    DMType dm_name = "DMMOFEM";
    CHKERR DMRegister_MoFEM(dm_name);
    //! [Register MoFEM discrete manager in PETSc
 
    //! [Create MoAB]
    moab::Core mb_instance;              ///< mesh database
    moab::Interface &moab = mb_instance; ///< mesh database interface
    //! [Create MoAB]
 
    //! [Create MoFEM]
    MoFEM::Core core(moab);           ///< finite element database
    MoFEM::Interface &m_field = core; ///< finite element database interface
    //! [Create MoFEM]
 
    //! [Load mesh]
    Simple *simple = m_field.getInterface<Simple>();
    CHKERR simple->getOptions();
    CHKERR simple->loadFile("");
    //! [Load mesh]
 
    //! [CONTACT]
    Contact ex(m_field);
    CHKERR ex.runProblem();
    //! [CONTACT]
  }
  CATCH_ERRORS;
 
  CHKERR MoFEM::Core::Finalize();
 
#ifdef PYTHON_SDF
  if (Py_FinalizeEx() < 0) {
    exit(120);
  }
#endif
 
  return 0;
}
 
struct SetUpSchurImpl : public SetUpSchur {
 
  SetUpSchurImpl(MoFEM::Interface &m_field) : SetUpSchur(), mField(m_field) {}
 
  virtual ~SetUpSchurImpl() {
    A.reset();
    P.reset();
    S.reset();
  }
 
  MoFEMErrorCode setUp(SmartPetscObj<TS> solver);
 
private:
  MoFEMErrorCode setEntities();
  MoFEMErrorCode setOperator();
  MoFEMErrorCode setPC(PC pc);
 
  SmartPetscObj<DM> createSubDM();
 
  SmartPetscObj<Mat> A;
  SmartPetscObj<Mat> P;
  SmartPetscObj<Mat> S;
 
  MoFEM::Interface &mField;
 
  SmartPetscObj<DM> subDM;
};
 
MoFEMErrorCode SetUpSchurImpl::setUp(SmartPetscObj<TS> solver) {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  auto pip = mField.getInterface<PipelineManager>();
  auto dm = simple->getDM();
 
  SNES snes;
  CHKERR TSGetSNES(solver, &snes);
  KSP ksp;
  CHKERR SNESGetKSP(snes, &ksp);
  CHKERR KSPSetFromOptions(ksp);
 
  PC pc;
  CHKERR KSPGetPC(ksp, &pc);
 
  PetscBool is_pcfs = PETSC_FALSE;
  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
  if (is_pcfs) {
 
    MOFEM_LOG("CONTACT", Sev::inform) << "Setup Schur pc";
 
    if (A || P || S) {
      CHK_THROW_MESSAGE(
          MOFEM_DATA_INCONSISTENCY,
          "It is expected that Schur matrix is not allocated. This is "
          "possible only if PC is set up twice");
    }
 
    A = createDMMatrix(dm);
    P = matDuplicate(A, MAT_DO_NOT_COPY_VALUES);
    subDM = createSubDM();
    S = createDMMatrix(subDM);
    CHKERR MatSetBlockSize(S, SPACE_DIM);
 
    auto ts_ctx_ptr = getDMTsCtx(dm);
    CHKERR TSSetIJacobian(solver, A, P, TsSetIJacobian, ts_ctx_ptr.get());
 
    CHKERR setOperator();
    CHKERR setPC(pc);
 
  } else {
    MOFEM_LOG("CONTACT", Sev::inform) << "No Schur pc";
 
    pip->getOpBoundaryLhsPipeline().push_front(new OpSchurAssembleBegin());
    pip->getOpBoundaryLhsPipeline().push_back(
        new OpSchurAssembleEnd<SCHUR_DGESV>({}, {}, {}, {}, {}));
    pip->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
    pip->getOpDomainLhsPipeline().push_back(
        new OpSchurAssembleEnd<SCHUR_DGESV>({}, {}, {}, {}, {}));
 
    auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
    post_proc_schur_lhs_ptr->postProcessHook = [this,
                                                post_proc_schur_lhs_ptr]() {
      MoFEMFunctionBegin;
      CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(
          mField, post_proc_schur_lhs_ptr, 1.)();
      MoFEMFunctionReturn(0);
    };
    auto ts_ctx_ptr = getDMTsCtx(dm);
    ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_schur_lhs_ptr);
  }
  MoFEMFunctionReturn(0);
}
 
SmartPetscObj<DM> SetUpSchurImpl::createSubDM() {
  auto simple = mField.getInterface<Simple>();
  auto sub_dm = createDM(mField.get_comm(), "DMMOFEM");
  auto set_up = [&]() {
    MoFEMFunctionBegin;
    CHKERR DMMoFEMCreateSubDM(sub_dm, simple->getDM(), "SUB");
    CHKERR DMMoFEMSetSquareProblem(sub_dm, PETSC_TRUE);
    CHKERR DMMoFEMAddElement(sub_dm, simple->getDomainFEName());
    CHKERR DMMoFEMAddSubFieldRow(sub_dm, "U");
    CHKERR DMSetUp(sub_dm);
    MoFEMFunctionReturn(0);
  };
  CHK_THROW_MESSAGE(set_up(), "sey up dm");
  return sub_dm;
}
 
MoFEMErrorCode SetUpSchurImpl::setOperator() {
  MoFEMFunctionBegin;
 
  double eps_stab = 1e-4;
  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-eps_stab", &eps_stab,
                               PETSC_NULL);
 
  using B =
      FormsIntegrators<BoundaryEleOpStab>::Assembly<SCHUR>::BiLinearForm<IT>;
  using OpMassStab = B::OpMass<3, SPACE_DIM * SPACE_DIM>;
 
  auto pip = mField.getInterface<PipelineManager>();
  // Boundary
  auto dm_is = getDMSubData(subDM)->getSmartRowIs();
  auto ao_up = createAOMappingIS(dm_is, PETSC_NULL);
 
  pip->getOpBoundaryLhsPipeline().push_front(new OpSchurAssembleBegin());
  pip->getOpBoundaryLhsPipeline().push_back(
      new OpMassStab("SIGMA", "SIGMA",
                     [eps_stab](double, double, double) { return eps_stab; }));
  pip->getOpBoundaryLhsPipeline().push_back(new OpSchurAssembleEnd<SCHUR_DGESV>(
      {"SIGMA"}, {nullptr}, {ao_up}, {S}, {false}, false));
 
  // Domain
  pip->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
  pip->getOpDomainLhsPipeline().push_back(new OpSchurAssembleEnd<SCHUR_DGESV>(
      {"SIGMA"}, {nullptr}, {ao_up}, {S}, {false}, false));
 
  auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
  auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
 
  pre_proc_schur_lhs_ptr->preProcessHook = [this]() {
    MoFEMFunctionBegin;
    CHKERR MatZeroEntries(A);
    CHKERR MatZeroEntries(P);
    CHKERR MatZeroEntries(S);
    MOFEM_LOG("CONTACT", Sev::verbose) << "Lhs Assemble Begin";
    MoFEMFunctionReturn(0);
  };
 
  post_proc_schur_lhs_ptr->postProcessHook = [this, ao_up,
                                              post_proc_schur_lhs_ptr]() {
    MoFEMFunctionBegin;
    MOFEM_LOG("CONTACT", Sev::verbose) << "Lhs Assemble End";
 
    *post_proc_schur_lhs_ptr->matAssembleSwitch = false;
 
    auto print_mat_norm = [this](auto a, std::string prefix) {
      MoFEMFunctionBegin;
      double nrm;
      CHKERR MatNorm(a, NORM_FROBENIUS, &nrm);
      MOFEM_LOG("CONTACT", Sev::noisy) << prefix << " norm = " << nrm;
      MoFEMFunctionReturn(0);
    };
 
    CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
    CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
    CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(
        mField, post_proc_schur_lhs_ptr, 1, A)();
 
    CHKERR MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY);
    CHKERR MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY);
    CHKERR MatAXPY(P, 1, A, SAME_NONZERO_PATTERN);
 
    CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
    CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
 
    CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(
        mField, post_proc_schur_lhs_ptr, 1, S, ao_up)();
 
#ifndef NDEBUG
    CHKERR print_mat_norm(A, "A");
    CHKERR print_mat_norm(P, "P");
    CHKERR print_mat_norm(S, "S");
#endif // NDEBUG
 
    MOFEM_LOG("CONTACT", Sev::verbose) << "Lhs Assemble Finish";
    MoFEMFunctionReturn(0);
  };
 
  auto simple = mField.getInterface<Simple>();
  auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
  ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_schur_lhs_ptr);
  ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_schur_lhs_ptr);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode SetUpSchurImpl::setPC(PC pc) {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  SmartPetscObj<IS> is;
  mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
      simple->getProblemName(), ROW, "SIGMA", 0, SPACE_DIM, is);
  CHKERR PCFieldSplitSetIS(pc, NULL, is);
  CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
  MoFEMFunctionReturn(0);
}
 
boost::shared_ptr<SetUpSchur>
SetUpSchur::createSetUpSchur(MoFEM::Interface &m_field) {
  return boost::shared_ptr<SetUpSchur>(new SetUpSchurImpl(m_field));
}

 
 
/**
 * \file ContactOps.hpp
 * \example ContactOps.hpp
 */
 
#ifndef __CONTACTOPS_HPP__
#define __CONTACTOPS_HPP__
 
namespace ContactOps {
 
//! [Common data]
struct CommonData : public boost::enable_shared_from_this<CommonData> {
  // MatrixDouble contactStress;
  MatrixDouble contactTraction;
  MatrixDouble contactDisp;
 
  VectorDouble sdfVals;  ///< size is equal to number of gauss points on element
  MatrixDouble gradsSdf; ///< nb of rows is equals to dimension, and nb of cols
                         ///< is equals to number of gauss points on element
  MatrixDouble hessSdf;  ///< nb of rows is equals to nb of element of symmetric
                        ///< matrix, and nb of cols is equals to number of gauss
                        ///< points on element
  VectorDouble constraintVals;
 
  static SmartPetscObj<Vec>
      totalTraction; // User have to release and create vector when appropiate.
 
  static auto createTotalTraction(MoFEM::Interface &m_field) {
    constexpr int ghosts[] = {0, 1, 2};
    totalTraction =
        createGhostVector(m_field.get_comm(),
 
                          (m_field.get_comm_rank() == 0) ? 3 : 0, 3,
 
                          (m_field.get_comm_rank() == 0) ? 0 : 3, ghosts);
    return totalTraction;
  }
 
  static auto getFTensor1TotalTraction() {
    if (CommonData::totalTraction) {
      const double *t_ptr;
      CHK_THROW_MESSAGE(VecGetArrayRead(CommonData::totalTraction, &t_ptr),
                        "get array");
      FTensor::Tensor1<double, 3> t{t_ptr[0], t_ptr[1], t_ptr[2]};
      CHK_THROW_MESSAGE(VecRestoreArrayRead(CommonData::totalTraction, &t_ptr),
                        "restore array");
      return t;
    } else {
      return FTensor::Tensor1<double, 3>{0., 0., 0.};
    }
  }
 
  inline auto contactTractionPtr() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(),
                                           &contactTraction);
  }
 
  inline auto contactDispPtr() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &contactDisp);
  }
 
  inline auto sdfPtr() {
    return boost::shared_ptr<VectorDouble>(shared_from_this(), &sdfVals);
  }
 
  inline auto gradSdfPtr() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &gradsSdf);
  }
 
  inline auto hessSdfPtr() {
    return boost::shared_ptr<MatrixDouble>(shared_from_this(), &hessSdf);
  }
 
  inline auto constraintPtr() {
    return boost::shared_ptr<VectorDouble>(shared_from_this(), &constraintVals);
  }
};
 
SmartPetscObj<Vec> CommonData::totalTraction;
 
//! [Common data]
 
//! [Surface distance function from python]
#ifdef PYTHON_SDF
struct SDFPython {
  SDFPython() = default;
  virtual ~SDFPython() = default;
 
  MoFEMErrorCode sdfInit(const std::string py_file) {
    MoFEMFunctionBegin;
    try {
 
      // create main module
      auto main_module = bp::import("__main__");
      mainNamespace = main_module.attr("__dict__");
      bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
      // create a reference to python function
      sdfFun = mainNamespace["sdf"];
      sdfGradFun = mainNamespace["grad_sdf"];
      sdfHessFun = mainNamespace["hess_sdf"];
 
    } catch (bp::error_already_set const &) {
      // print all other errors to stderr
      PyErr_Print();
      CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
    }
    MoFEMFunctionReturn(0);
  };
 
  template <typename T>
  inline std::vector<T>
  py_list_to_std_vector(const boost::python::object &iterable) {
    return std::vector<T>(boost::python::stl_input_iterator<T>(iterable),
                          boost::python::stl_input_iterator<T>());
  }
 
  MoFEMErrorCode evalSdf(
 
      double delta_t, double t, np::ndarray x, np::ndarray y, np::ndarray z,
      np::ndarray tx, np::ndarray ty, np::ndarray tz, int block_id,
      np::ndarray &sdf
 
  ) {
    MoFEMFunctionBegin;
    try {
 
      // call python function
      sdf = bp::extract<np::ndarray>(
          sdfFun(delta_t, t, x, y, z, tx, ty, tz, block_id));
 
    } catch (bp::error_already_set const &) {
      // print all other errors to stderr
      PyErr_Print();
      CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
    }
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode evalGradSdf(
 
      double delta_t, double t, np::ndarray x, np::ndarray y, np::ndarray z,
      np::ndarray tx, np::ndarray ty, np::ndarray tz, int block_id,
      np::ndarray &grad_sdf
 
  ) {
    MoFEMFunctionBegin;
    try {
 
      // call python function
      grad_sdf = bp::extract<np::ndarray>(
          sdfGradFun(delta_t, t, x, y, z, tx, ty, tz, block_id));
 
    } catch (bp::error_already_set const &) {
      // print all other errors to stderr
      PyErr_Print();
      CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
    }
    MoFEMFunctionReturn(0);
  }
 
  MoFEMErrorCode evalHessSdf(
 
      double delta_t, double t, np::ndarray x, np::ndarray y, np::ndarray z,
      np::ndarray tx, np::ndarray ty, np::ndarray tz, int block_id,
      np::ndarray &hess_sdf
 
  ) {
    MoFEMFunctionBegin;
    try {
 
      // call python function
      hess_sdf = bp::extract<np::ndarray>(
          sdfHessFun(delta_t, t, x, y, z, tx, ty, tz, block_id));
 
    } catch (bp::error_already_set const &) {
      // print all other errors to stderr
      PyErr_Print();
      CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
    }
    MoFEMFunctionReturn(0);
  }
 
private:
  bp::object mainNamespace;
  bp::object sdfFun;
  bp::object sdfGradFun;
  bp::object sdfHessFun;
};
 
static boost::weak_ptr<SDFPython> sdfPythonWeakPtr;
 
inline np::ndarray convert_to_numpy(VectorDouble &data, int nb_gauss_pts,
                                    int id) {
  auto dtype = np::dtype::get_builtin<double>();
  auto size = bp::make_tuple(nb_gauss_pts);
  auto stride = bp::make_tuple(3 * sizeof(double));
  return (np::from_data(&data[id], dtype, size, stride, bp::object()));
};
#endif
//! [Surface distance function from python]
 
using SurfaceDistanceFunction = boost::function<VectorDouble(
    double delta_t, double t, int nb_gauss_pts, MatrixDouble &spatial_coords,
    MatrixDouble &normals_at_pts, int block_id)>;
 
using GradSurfaceDistanceFunction = boost::function<MatrixDouble(
    double delta_t, double t, int nb_gauss_pts, MatrixDouble &spatial_coords,
    MatrixDouble &normals_at_pts, int block_id)>;
 
using HessSurfaceDistanceFunction = boost::function<MatrixDouble(
    double delta_t, double t, int nb_gauss_pts, MatrixDouble &spatial_coords,
    MatrixDouble &normals_at_pts, int block_id)>;
 
inline VectorDouble surface_distance_function(double delta_t, double t,
                                              int nb_gauss_pts,
                                              MatrixDouble &m_spatial_coords,
                                              MatrixDouble &m_normals_at_pts,
                                              int block_id) {
 
#ifdef PYTHON_SDF
  if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
 
    VectorDouble v_spatial_coords = m_spatial_coords.data();
    VectorDouble v_normal_at_pts = m_normals_at_pts.data();
 
    bp::list python_coords;
    bp::list python_normals;
 
    for (int idx = 0; idx < 3; ++idx) {
      python_coords.append(
          convert_to_numpy(v_spatial_coords, nb_gauss_pts, idx));
      python_normals.append(
          convert_to_numpy(v_normal_at_pts, nb_gauss_pts, idx));
    }
 
    np::ndarray np_sdf = np::empty(bp::make_tuple(nb_gauss_pts),
                                   np::dtype::get_builtin<double>());
    CHK_MOAB_THROW(sdf_ptr->evalSdf(delta_t, t,
                                    bp::extract<np::ndarray>(python_coords[0]),
                                    bp::extract<np::ndarray>(python_coords[1]),
                                    bp::extract<np::ndarray>(python_coords[2]),
                                    bp::extract<np::ndarray>(python_normals[0]),
                                    bp::extract<np::ndarray>(python_normals[1]),
                                    bp::extract<np::ndarray>(python_normals[2]),
                                    block_id, np_sdf),
                   "Failed python call");
 
    double *sdf_val_ptr = reinterpret_cast<double *>(np_sdf.get_data());
 
    VectorDouble v_sdf;
    v_sdf.resize(nb_gauss_pts, false);
 
    for (size_t gg = 0; gg < nb_gauss_pts; ++gg)
      v_sdf[gg] = *(sdf_val_ptr + gg);
 
    return v_sdf;
  }
#endif
  VectorDouble v_sdf;
  v_sdf.resize(nb_gauss_pts, false);
  m_spatial_coords.resize(3, nb_gauss_pts);
  auto t_coords = getFTensor1FromMat<3>(m_spatial_coords);
 
  for (size_t gg = 0; gg < nb_gauss_pts; ++gg)
    v_sdf[gg] = t_coords(1) + 0.5;
 
  return v_sdf;
}
 
inline MatrixDouble
grad_surface_distance_function(double delta_t, double t, int nb_gauss_pts,
                               MatrixDouble &m_spatial_coords,
                               MatrixDouble &m_normals_at_pts, int block_id) {
#ifdef PYTHON_SDF
  if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
 
    VectorDouble v_spatial_coords = m_spatial_coords.data();
    VectorDouble v_normal_at_pts = m_normals_at_pts.data();
 
    bp::list python_coords;
    bp::list python_normals;
 
    for (int idx = 0; idx < 3; ++idx) {
      python_coords.append(
          convert_to_numpy(v_spatial_coords, nb_gauss_pts, idx));
      python_normals.append(
          convert_to_numpy(v_normal_at_pts, nb_gauss_pts, idx));
    }
 
    np::ndarray np_grad_sdf = np::empty(bp::make_tuple(nb_gauss_pts, 3),
                                        np::dtype::get_builtin<double>());
    CHK_MOAB_THROW(sdf_ptr->evalGradSdf(
                       delta_t, t, bp::extract<np::ndarray>(python_coords[0]),
                       bp::extract<np::ndarray>(python_coords[1]),
                       bp::extract<np::ndarray>(python_coords[2]),
                       bp::extract<np::ndarray>(python_normals[0]),
                       bp::extract<np::ndarray>(python_normals[1]),
                       bp::extract<np::ndarray>(python_normals[2]), block_id,
                       np_grad_sdf),
                   "Failed python call");
 
    double *grad_ptr = reinterpret_cast<double *>(np_grad_sdf.get_data());
 
    MatrixDouble m_grad_sdf;
    m_grad_sdf.resize(3, nb_gauss_pts, false);
    for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
      for (int idx = 0; idx < 3; ++idx)
        m_grad_sdf(idx, gg) = *(grad_ptr + (3 * gg + idx));
    }
    return m_grad_sdf;
  }
#endif
  MatrixDouble m_grad_sdf;
  m_grad_sdf.resize(3, nb_gauss_pts, false);
  FTensor::Index<'i', 3> i;
  FTensor::Tensor1<double, 3> t_grad_sdf_set{0.0, 0.0, 1.0};
  auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
 
  for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
    t_grad_sdf(i) = t_grad_sdf_set(i);
    ++t_grad_sdf;
  }
 
  return m_grad_sdf;
}
 
inline MatrixDouble
hess_surface_distance_function(double delta_t, double t, int nb_gauss_pts,
                               MatrixDouble &m_spatial_coords,
                               MatrixDouble &m_normals_at_pts, int block_id) {
#ifdef PYTHON_SDF
  if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
 
    VectorDouble v_spatial_coords = m_spatial_coords.data();
    VectorDouble v_normal_at_pts = m_normals_at_pts.data();
 
    bp::list python_coords;
    bp::list python_normals;
 
    for (int idx = 0; idx < 3; ++idx) {
      python_coords.append(
          convert_to_numpy(v_spatial_coords, nb_gauss_pts, idx));
      python_normals.append(
          convert_to_numpy(v_normal_at_pts, nb_gauss_pts, idx));
    };
 
    np::ndarray np_hess_sdf = np::empty(bp::make_tuple(nb_gauss_pts, 6),
                                        np::dtype::get_builtin<double>());
    CHK_MOAB_THROW(sdf_ptr->evalHessSdf(
                       delta_t, t, bp::extract<np::ndarray>(python_coords[0]),
                       bp::extract<np::ndarray>(python_coords[1]),
                       bp::extract<np::ndarray>(python_coords[2]),
                       bp::extract<np::ndarray>(python_normals[0]),
                       bp::extract<np::ndarray>(python_normals[1]),
                       bp::extract<np::ndarray>(python_normals[2]), block_id,
                       np_hess_sdf),
                   "Failed python call");
 
    double *hess_ptr = reinterpret_cast<double *>(np_hess_sdf.get_data());
 
    MatrixDouble m_hess_sdf;
    m_hess_sdf.resize(6, nb_gauss_pts, false);
    for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
      for (int idx = 0; idx < 6; ++idx)
        m_hess_sdf(idx, gg) =
            *(hess_ptr + (6 * gg + idx));
    }
    return m_hess_sdf;
  }
#endif
  MatrixDouble m_hess_sdf;
  m_hess_sdf.resize(6, nb_gauss_pts, false);
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Tensor2_symmetric<double, 3> t_hess_sdf_set{0., 0., 0., 0., 0., 0.};
  auto t_hess_sdf = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
 
  for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
    t_hess_sdf(i, j) = t_hess_sdf_set(i, j);
    ++t_hess_sdf;
  }
  return m_hess_sdf;
}
 
 
 
template <int DIM, IntegrationType I, typename BoundaryEleOp>
struct OpAssembleTotalContactTractionImpl;
 
template <int DIM, IntegrationType I, typename BoundaryEleOp>
struct OpEvaluateSDFImpl;
 
template <int DIM, IntegrationType I, typename AssemblyBoundaryEleOp>
struct OpConstrainBoundaryRhsImpl;
 
template <int DIM, IntegrationType I, typename AssemblyBoundaryEleOp>
struct OpConstrainBoundaryLhs_dUImpl;
 
template <int DIM, IntegrationType I, typename AssemblyBoundaryEleOp>
struct OpConstrainBoundaryLhs_dTractionImpl;
 
template <typename T1, typename T2, int DIM1, int DIM2>
inline auto get_spatial_coords(FTensor::Tensor1<T1, DIM1> &&t_coords,
                               FTensor::Tensor1<T2, DIM2> &&t_disp,
                               size_t nb_gauss_pts) {
  MatrixDouble m_spatial_coords(nb_gauss_pts, 3);
  m_spatial_coords.clear();
  auto t_spatial_coords = getFTensor1FromPtr<3>(&m_spatial_coords(0, 0));
  FTensor::Index<'i', DIM2> i;
  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
    t_spatial_coords(i) = t_coords(i) + t_disp(i);
    ++t_spatial_coords;
    ++t_coords;
    ++t_disp;
  }
  return m_spatial_coords;
}
 
template <typename T1, int DIM1>
inline auto get_normalize_normals(FTensor::Tensor1<T1, DIM1> &&t_normal_at_pts,
                                  size_t nb_gauss_pts) {
  MatrixDouble m_normals_at_pts(3, nb_gauss_pts);
  m_normals_at_pts.clear();
  FTensor::Index<'i', DIM1> i;
  auto t_set_normal = getFTensor1FromMat<3>(m_normals_at_pts);
  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
    t_set_normal(i) = t_normal_at_pts(i) / t_normal_at_pts.l2();
    ++t_set_normal;
    ++t_normal_at_pts;
  }
  return m_normals_at_pts;
}
 
template <int DIM, typename BoundaryEleOp>
struct OpAssembleTotalContactTractionImpl<DIM, GAUSS, BoundaryEleOp>
    : public BoundaryEleOp {
  OpAssembleTotalContactTractionImpl(
      boost::shared_ptr<CommonData> common_data_ptr, double scale = 1,
      bool is_axisymmetric = false);
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
  const double scaleTraction;
  bool isAxisymmetric;
};
 
template <int DIM, typename BoundaryEleOp>
struct OpEvaluateSDFImpl<DIM, GAUSS, BoundaryEleOp> : public BoundaryEleOp {
  OpEvaluateSDFImpl(boost::shared_ptr<CommonData> common_data_ptr);
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
 
  SurfaceDistanceFunction surfaceDistanceFunction = surface_distance_function;
  GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
      grad_surface_distance_function;
  HessSurfaceDistanceFunction hessSurfaceDistanceFunction =
      hess_surface_distance_function;
};
 
template <int DIM, typename AssemblyBoundaryEleOp>
struct OpConstrainBoundaryRhsImpl<DIM, GAUSS, AssemblyBoundaryEleOp>
    : public AssemblyBoundaryEleOp {
  OpConstrainBoundaryRhsImpl(const std::string field_name,
                             boost::shared_ptr<CommonData> common_data_ptr,
                             bool is_axisymmetric = false);
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
 
  SurfaceDistanceFunction surfaceDistanceFunction = surface_distance_function;
  GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
      grad_surface_distance_function;
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
  bool isAxisymmetric;
};
 
template <int DIM, typename AssemblyBoundaryEleOp>
struct OpConstrainBoundaryLhs_dUImpl<DIM, GAUSS, AssemblyBoundaryEleOp>
    : public AssemblyBoundaryEleOp {
  OpConstrainBoundaryLhs_dUImpl(const std::string row_field_name,
                                const std::string col_field_name,
                                boost::shared_ptr<CommonData> common_data_ptr,
                                bool is_axisymmetric = false);
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
                           EntitiesFieldData::EntData &col_data);
 
  SurfaceDistanceFunction surfaceDistanceFunction = surface_distance_function;
  GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
      grad_surface_distance_function;
  HessSurfaceDistanceFunction hessSurfaceDistanceFunction =
      hess_surface_distance_function;
 
  boost::shared_ptr<CommonData> commonDataPtr;
  bool isAxisymmetric;
};
 
template <int DIM, typename AssemblyBoundaryEleOp>
struct OpConstrainBoundaryLhs_dTractionImpl<DIM, GAUSS, AssemblyBoundaryEleOp>
    : public AssemblyBoundaryEleOp {
  OpConstrainBoundaryLhs_dTractionImpl(
      const std::string row_field_name, const std::string col_field_name,
      boost::shared_ptr<CommonData> common_data_ptr,
      bool is_axisymmetric = false);
  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
                           EntitiesFieldData::EntData &col_data);
 
  SurfaceDistanceFunction surfaceDistanceFunction = surface_distance_function;
  GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
      grad_surface_distance_function;
 
private:
  boost::shared_ptr<CommonData> commonDataPtr;
  bool isAxisymmetric;
};
 
template <typename BoundaryEleOp> struct ContactIntegrators {
 
  template <int DIM, IntegrationType I>
  using OpAssembleTotalContactTraction =
      OpAssembleTotalContactTractionImpl<DIM, I, BoundaryEleOp>;
 
  template <int DIM, IntegrationType I>
  using OpEvaluateSDF = OpEvaluateSDFImpl<DIM, I, BoundaryEleOp>;
 
  template <AssemblyType A> struct Assembly {
 
    using AssemblyBoundaryEleOp =
        typename FormsIntegrators<BoundaryEleOp>::template Assembly<A>::OpBase;
 
    template <int DIM, IntegrationType I>
    using OpConstrainBoundaryRhs =
        OpConstrainBoundaryRhsImpl<DIM, I, AssemblyBoundaryEleOp>;
 
    template <int DIM, IntegrationType I>
    using OpConstrainBoundaryLhs_dU =
        OpConstrainBoundaryLhs_dUImpl<DIM, I, AssemblyBoundaryEleOp>;
 
    template <int DIM, IntegrationType I>
    using OpConstrainBoundaryLhs_dTraction =
        OpConstrainBoundaryLhs_dTractionImpl<DIM, I, AssemblyBoundaryEleOp>;
  };
};
 
inline double sign(double x) {
  constexpr auto eps = std::numeric_limits<float>::epsilon();
  if (std::abs(x) < eps)
    return 0;
  else if (x > eps)
    return 1;
  else
    return -1;
};
 
inline double w(const double sdf, const double tn) {
  return sdf - cn_contact * tn;
}
 
/**
 * @brief constrain function
 *
 * return 1 if negative sdf or positive tn
 *
 * @param sdf signed distance
 * @param tn traction
 * @return double
 */
inline double constrain(double sdf, double tn) {
  const auto s = sign(w(sdf, tn));
  return (1 - s) / 2;
}
 
template <int DIM, typename BoundaryEleOp>
OpAssembleTotalContactTractionImpl<DIM, GAUSS, BoundaryEleOp>::
    OpAssembleTotalContactTractionImpl(
        boost::shared_ptr<CommonData> common_data_ptr, double scale,
        bool is_axisymmetric)
    : BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE),
      commonDataPtr(common_data_ptr), scaleTraction(scale),
      isAxisymmetric(is_axisymmetric) {}
 
template <int DIM, typename BoundaryEleOp>
MoFEMErrorCode
OpAssembleTotalContactTractionImpl<DIM, GAUSS, BoundaryEleOp>::doWork(
    int side, EntityType type, EntData &data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'i', DIM> i;
  FTensor::Tensor1<double, 3> t_sum_t{0., 0., 0.};
 
  auto t_w = BoundaryEleOp::getFTensor0IntegrationWeight();
  auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
  auto t_coords = BoundaryEleOp::getFTensor1CoordsAtGaussPts();
 
  const auto nb_gauss_pts = BoundaryEleOp::getGaussPts().size2();
  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
    double jacobian = 1.;
    if (isAxisymmetric) {
      jacobian = 2. * M_PI * t_coords(0);
    }
    const double alpha = t_w * jacobian * BoundaryEleOp::getMeasure();
    t_sum_t(i) += alpha * t_traction(i);
    ++t_w;
    ++t_traction;
    ++t_coords;
  }
 
  t_sum_t(i) *= scaleTraction;
 
  constexpr int ind[] = {0, 1, 2};
  CHKERR VecSetValues(commonDataPtr->totalTraction, 3, ind, &t_sum_t(0),
                      ADD_VALUES);
 
  MoFEMFunctionReturn(0);
}
 
template <int DIM, typename BoundaryEleOp>
OpEvaluateSDFImpl<DIM, GAUSS, BoundaryEleOp>::OpEvaluateSDFImpl(
    boost::shared_ptr<CommonData> common_data_ptr)
    : BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE),
      commonDataPtr(common_data_ptr) {}
 
template <int DIM, typename BoundaryEleOp>
MoFEMErrorCode
OpEvaluateSDFImpl<DIM, GAUSS, BoundaryEleOp>::doWork(int side, EntityType type,
                                                     EntData &data) {
  MoFEMFunctionBegin;
 
  const auto nb_gauss_pts = BoundaryEleOp::getGaussPts().size2();
  auto &sdf_vec = commonDataPtr->sdfVals;
  auto &grad_mat = commonDataPtr->gradsSdf;
  auto &hess_mat = commonDataPtr->hessSdf;
  auto &constraint_vec = commonDataPtr->constraintVals;
  auto &contactTraction_mat = commonDataPtr->contactTraction;
 
  sdf_vec.resize(nb_gauss_pts, false);
  grad_mat.resize(DIM, nb_gauss_pts, false);
  hess_mat.resize((DIM * (DIM + 1)) / 2, nb_gauss_pts, false);
  constraint_vec.resize(nb_gauss_pts, false);
 
  auto t_traction = getFTensor1FromMat<DIM>(contactTraction_mat);
 
  auto t_sdf = getFTensor0FromVec(sdf_vec);
  auto t_grad_sdf = getFTensor1FromMat<DIM>(grad_mat);
  auto t_hess_sdf = getFTensor2SymmetricFromMat<DIM>(hess_mat);
  auto t_constraint = getFTensor0FromVec(constraint_vec);
 
  auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
  auto t_coords = BoundaryEleOp::getFTensor1CoordsAtGaussPts();
  auto t_normal_at_pts = BoundaryEleOp::getFTensor1NormalsAtGaussPts();
 
  FTensor::Index<'i', DIM> i;
  FTensor::Index<'j', DIM> j;
 
  auto ts_time = BoundaryEleOp::getTStime();
  auto ts_time_step = BoundaryEleOp::getTStimeStep();
 
  auto m_spatial_coords = get_spatial_coords(
      BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
      getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
  auto m_normals_at_pts = get_normalize_normals(
      BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
 
  // placeholder to pass boundary block id to python
  int block_id = 0;
 
  auto v_sdf =
      surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
                              m_spatial_coords, m_normals_at_pts, block_id);
 
  auto m_grad_sdf =
      gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
                                  m_spatial_coords, m_normals_at_pts, block_id);
 
  auto m_hess_sdf =
      hessSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
                                  m_spatial_coords, m_normals_at_pts, block_id);
 
  auto t_sdf_v = getFTensor0FromVec(v_sdf);
  auto t_grad_sdf_v = getFTensor1FromMat<3>(m_grad_sdf);
  auto t_hess_sdf_v = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
 
  auto next = [&]() {
    ++t_sdf;
    ++t_sdf_v;
    ++t_grad_sdf;
    ++t_grad_sdf_v;
    ++t_hess_sdf;
    ++t_hess_sdf_v;
    ++t_disp;
    ++t_traction;
    ++t_constraint;
  };
 
  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
 
    auto tn = -t_traction(i) * t_grad_sdf_v(i);
    auto c = constrain(t_sdf_v, tn);
 
    t_sdf = t_sdf_v;
    t_grad_sdf(i) = t_grad_sdf_v(i);
    t_hess_sdf(i, j) = t_hess_sdf_v(i, j);
    t_constraint = c;
 
    next();
  }
 
  MoFEMFunctionReturn(0);
}
 
template <int DIM, typename AssemblyBoundaryEleOp>
OpConstrainBoundaryRhsImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::
    OpConstrainBoundaryRhsImpl(const std::string field_name,
                               boost::shared_ptr<CommonData> common_data_ptr,
                               bool is_axisymmetric)
    : AssemblyBoundaryEleOp(field_name, field_name,
                            AssemblyBoundaryEleOp::OPROW),
      commonDataPtr(common_data_ptr), isAxisymmetric(is_axisymmetric) {}
 
template <int DIM, typename AssemblyBoundaryEleOp>
MoFEMErrorCode
OpConstrainBoundaryRhsImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::iNtegrate(
    EntitiesFieldData::EntData &data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'i', DIM> i;
  FTensor::Index<'j', DIM> j;
  FTensor::Index<'k', DIM> k;
  FTensor::Index<'l', DIM> l;
 
  const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
 
  auto &nf = AssemblyBoundaryEleOp::locF;
 
  auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
 
  auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
  auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
  auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
  auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
 
  size_t nb_base_functions = data.getN().size2() / 3;
  auto t_base = data.getFTensor1N<3>();
 
  auto m_spatial_coords = get_spatial_coords(
      BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
      getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
  auto m_normals_at_pts = get_normalize_normals(
      BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
 
  auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
 
  auto ts_time = AssemblyBoundaryEleOp::getTStime();
  auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
 
  // placeholder to pass boundary block id to python
  int block_id = 0;
 
  auto v_sdf =
      surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
                              m_spatial_coords, m_normals_at_pts, block_id);
 
  auto m_grad_sdf =
      gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
                                  m_spatial_coords, m_normals_at_pts, block_id);
 
  auto t_sdf = getFTensor0FromVec(v_sdf);
  auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
 
  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
 
    auto t_nf = getFTensor1FromPtr<DIM>(&nf[0]);
 
    double jacobian = 1.;
    if (isAxisymmetric) {
      jacobian = 2. * M_PI * t_coords(0);
    }
    const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
 
    auto tn = -t_traction(i) * t_grad_sdf(i);
    auto c = constrain(t_sdf, tn);
 
    FTensor::Tensor2<double, DIM, DIM> t_cP;
    t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
    FTensor::Tensor2<double, DIM, DIM> t_cQ;
    t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
 
    FTensor::Tensor1<double, DIM> t_rhs;
    t_rhs(i) =
 
        t_cQ(i, j) * (t_disp(j) - cn_contact * t_traction(j))
 
        +
 
        t_cP(i, j) * t_disp(j) +
        c * (t_sdf * t_grad_sdf(i)); // add gap0 displacements
 
    size_t bb = 0;
    for (; bb != AssemblyBoundaryEleOp::nbRows / DIM; ++bb) {
      const double beta = alpha * (t_base(i) * t_normal(i));
      t_nf(i) -= beta * t_rhs(i);
 
      ++t_nf;
      ++t_base;
    }
    for (; bb < nb_base_functions; ++bb)
      ++t_base;
 
    ++t_disp;
    ++t_traction;
    ++t_coords;
    ++t_w;
    ++t_normal;
    ++t_sdf;
    ++t_grad_sdf;
  }
 
  MoFEMFunctionReturn(0);
}
 
template <int DIM, typename AssemblyBoundaryEleOp>
OpConstrainBoundaryLhs_dUImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::
    OpConstrainBoundaryLhs_dUImpl(const std::string row_field_name,
                                  const std::string col_field_name,
                                  boost::shared_ptr<CommonData> common_data_ptr,
                                  bool is_axisymmetric)
    : AssemblyBoundaryEleOp(row_field_name, col_field_name,
                            AssemblyBoundaryEleOp::OPROWCOL),
      commonDataPtr(common_data_ptr), isAxisymmetric(is_axisymmetric) {
  AssemblyBoundaryEleOp::sYmm = false;
}
 
template <int DIM, typename AssemblyBoundaryEleOp>
MoFEMErrorCode
OpConstrainBoundaryLhs_dUImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::iNtegrate(
    EntitiesFieldData::EntData &row_data,
    EntitiesFieldData::EntData &col_data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'i', DIM> i;
  FTensor::Index<'j', DIM> j;
  FTensor::Index<'k', DIM> k;
 
  const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
  auto &locMat = AssemblyBoundaryEleOp::locMat;
 
  auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
  auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
  auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
 
  auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
  auto t_row_base = row_data.getFTensor1N<3>();
  size_t nb_face_functions = row_data.getN().size2() / 3;
 
  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
 
  auto m_spatial_coords = get_spatial_coords(
      BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
      getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
  auto m_normals_at_pts = get_normalize_normals(
      BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
 
  auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
 
  auto ts_time = AssemblyBoundaryEleOp::getTStime();
  auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
 
  // placeholder to pass boundary block id to python
  int block_id = 0;
 
  auto v_sdf =
      surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
                              m_spatial_coords, m_normals_at_pts, block_id);
 
  auto m_grad_sdf =
      gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
                                  m_spatial_coords, m_normals_at_pts, block_id);
 
  auto m_hess_sdf =
      hessSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
                                  m_spatial_coords, m_normals_at_pts, block_id);
 
  auto t_sdf = getFTensor0FromVec(v_sdf);
  auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
  auto t_hess_sdf = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
 
  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
 
    double jacobian = 1.;
    if (isAxisymmetric) {
      jacobian = 2. * M_PI * t_coords(0);
    }
    const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
 
    auto tn = -t_traction(i) * t_grad_sdf(i);
    auto c = constrain(t_sdf, tn);
 
    FTensor::Tensor2<double, DIM, DIM> t_cP;
    t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
    FTensor::Tensor2<double, DIM, DIM> t_cQ;
    t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
 
    FTensor::Tensor2<double, DIM, DIM> t_res_dU;
    t_res_dU(i, j) = kronecker_delta(i, j) + t_cP(i, j);
 
    if (c > 0) {
      t_res_dU(i, j) +=
          (c * cn_contact) *
              (t_hess_sdf(i, j) * (t_grad_sdf(k) * t_traction(k)) +
               t_grad_sdf(i) * t_hess_sdf(k, j) * t_traction(k)) +
          c * t_sdf * t_hess_sdf(i, j);
    }
 
    size_t rr = 0;
    for (; rr != AssemblyBoundaryEleOp::nbRows / DIM; ++rr) {
 
      auto t_mat = getFTensor2FromArray<DIM, DIM, DIM>(locMat, DIM * 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 != AssemblyBoundaryEleOp::nbCols / DIM; ++cc) {
        const double beta = alpha * row_base * t_col_base;
 
        t_mat(i, j) -= beta * t_res_dU(i, j);
 
        ++t_col_base;
        ++t_mat;
      }
 
      ++t_row_base;
    }
    for (; rr < nb_face_functions; ++rr)
      ++t_row_base;
 
    ++t_traction;
    ++t_coords;
    ++t_w;
    ++t_normal;
    ++t_sdf;
    ++t_grad_sdf;
    ++t_hess_sdf;
  }
 
  MoFEMFunctionReturn(0);
}
 
template <int DIM, typename AssemblyBoundaryEleOp>
OpConstrainBoundaryLhs_dTractionImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::
    OpConstrainBoundaryLhs_dTractionImpl(
        const std::string row_field_name, const std::string col_field_name,
        boost::shared_ptr<CommonData> common_data_ptr, bool is_axisymmetric)
    : AssemblyBoundaryEleOp(row_field_name, col_field_name,
                            AssemblyBoundaryEleOp::OPROWCOL),
      commonDataPtr(common_data_ptr), isAxisymmetric(is_axisymmetric) {
  AssemblyBoundaryEleOp::sYmm = false;
}
 
template <int DIM, typename AssemblyBoundaryEleOp>
MoFEMErrorCode
OpConstrainBoundaryLhs_dTractionImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::
    iNtegrate(EntitiesFieldData::EntData &row_data,
              EntitiesFieldData::EntData &col_data) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'i', DIM> i;
  FTensor::Index<'j', DIM> j;
  FTensor::Index<'k', DIM> k;
 
  const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
  auto &locMat = AssemblyBoundaryEleOp::locMat;
 
  auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
  auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
  auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
 
  auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
  auto t_row_base = row_data.getFTensor1N<3>();
  size_t nb_face_functions = row_data.getN().size2() / 3;
 
  auto m_spatial_coords = get_spatial_coords(
      BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
      getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
  auto m_normals_at_pts = get_normalize_normals(
      BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
 
  auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
 
  auto ts_time = AssemblyBoundaryEleOp::getTStime();
  auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
 
  // placeholder to pass boundary block id to python
  int block_id = 0;
 
  auto v_sdf =
      surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
                              m_spatial_coords, m_normals_at_pts, block_id);
 
  auto m_grad_sdf =
      gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
                                  m_spatial_coords, m_normals_at_pts, block_id);
 
  auto t_sdf = getFTensor0FromVec(v_sdf);
  auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
 
  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
 
    double jacobian = 1.;
    if (isAxisymmetric) {
      jacobian = 2. * M_PI * t_coords(0);
    }
    const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
 
    auto tn = -t_traction(i) * t_grad_sdf(i);
    auto c = constrain(t_sdf, tn);
 
    FTensor::Tensor2<double, DIM, DIM> t_cP;
    t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
    FTensor::Tensor2<double, DIM, DIM> t_cQ;
    t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
 
    FTensor::Tensor2<double, DIM, DIM> t_res_dt;
    t_res_dt(i, j) = -cn_contact * t_cQ(i, j);
 
    size_t rr = 0;
    for (; rr != AssemblyBoundaryEleOp::nbRows / DIM; ++rr) {
 
      auto t_mat = getFTensor2FromArray<DIM, DIM, DIM>(locMat, DIM * 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 != AssemblyBoundaryEleOp::nbCols / DIM; ++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 * t_res_dt(i, j);
 
        ++t_col_base;
        ++t_mat;
      }
 
      ++t_row_base;
    }
    for (; rr < nb_face_functions; ++rr)
      ++t_row_base;
 
    ++t_traction;
    ++t_coords;
    ++t_w;
    ++t_normal;
    ++t_sdf;
    ++t_grad_sdf;
  }
 
  MoFEMFunctionReturn(0);
}
 
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
MoFEMErrorCode opFactoryDomainRhs(
    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
    std::string sigma, std::string u, bool is_axisymmetric = false) {
  MoFEMFunctionBegin;
 
  using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
      A>::template LinearForm<I>;
  using OpMixDivURhs = typename B::template OpMixDivTimesU<3, DIM, DIM>;
  using OpMixDivUCylRhs =
      typename B::template OpMixDivTimesU<3, DIM, DIM, CYLINDRICAL>;
 
  using OpMixLambdaGradURhs = typename B::template OpMixTensorTimesGradU<DIM>;
  using OpMixUTimesDivLambdaRhs =
      typename B::template OpMixVecTimesDivLambda<SPACE_DIM>;
  using OpMixUTimesLambdaRhs =
      typename B::template OpGradTimesTensor<1, DIM, DIM>;
 
  auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
  auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
  auto div_stress_ptr = boost::make_shared<MatrixDouble>();
  auto contact_stress_ptr = boost::make_shared<MatrixDouble>();
 
  auto jacobian = [is_axisymmetric](const double r, const double,
                                    const double) {
    if (is_axisymmetric)
      return 2. * M_PI * r;
    else
      return 1.;
  };
 
  pip.push_back(new OpCalculateVectorFieldValues<DIM>(
      u, common_data_ptr->contactDispPtr()));
  pip.push_back(
      new OpCalculateHVecTensorField<DIM, DIM>(sigma, contact_stress_ptr));
 
  if (!is_axisymmetric) {
    pip.push_back(
        new OpCalculateHVecTensorDivergence<DIM, DIM>(sigma, div_stress_ptr));
  } else {
    pip.push_back(new OpCalculateHVecTensorDivergence<DIM, DIM, CYLINDRICAL>(
        sigma, div_stress_ptr));
  }
 
  pip.push_back(new OpCalculateVectorFieldGradient<DIM, DIM>(u, mat_grad_ptr));
 
  if (!is_axisymmetric) {
    pip.push_back(
        new OpMixDivURhs(sigma, common_data_ptr->contactDispPtr(), jacobian));
  } else {
    pip.push_back(new OpMixDivUCylRhs(sigma, common_data_ptr->contactDispPtr(),
                                      jacobian));
  }
 
  pip.push_back(new OpMixLambdaGradURhs(sigma, mat_grad_ptr, jacobian));
  pip.push_back(new OpMixUTimesDivLambdaRhs(u, div_stress_ptr, jacobian));
  pip.push_back(new OpMixUTimesLambdaRhs(u, contact_stress_ptr, jacobian));
 
  MoFEMFunctionReturn(0);
}
 
template <typename OpMixLhs> struct OpMixLhsSide : public OpMixLhs {
  using OpMixLhs::OpMixLhs;
  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
                        EntityType col_type,
                        EntitiesFieldData::EntData &row_data,
                        EntitiesFieldData::EntData &col_data) {
    MoFEMFunctionBegin;
    auto side_fe_entity = OpMixLhs::getSidePtrFE()->getFEEntityHandle();
    auto side_fe_data = OpMixLhs::getSideEntity(row_side, row_type);
    // Only assemble side which correspond to edge entity on boundary
    if (side_fe_entity == side_fe_data) {
      CHKERR OpMixLhs::doWork(row_side, col_side, row_type, col_type, row_data,
                              col_data);
    }
    MoFEMFunctionReturn(0);
  }
};
 
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEle>
MoFEMErrorCode opFactoryBoundaryToDomainLhs(
    MoFEM::Interface &m_field,
    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
    std::string fe_domain_name, std::string sigma, std::string u,
    std::string geom, ForcesAndSourcesCore::RuleHookFun rule,
    bool is_axisymmetric = false) {
  MoFEMFunctionBegin;
 
  using DomainEleOp = typename DomainEle::UserDataOperator;
 
  auto op_loop_side = new OpLoopSide<DomainEle>(
      m_field, fe_domain_name, DIM, Sev::noisy,
      boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>());
  pip.push_back(op_loop_side);
 
  CHKERR AddHOOps<DIM, DIM, DIM>::add(op_loop_side->getOpPtrVector(),
                                      {H1, HDIV}, geom);
 
  using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
      A>::template BiLinearForm<I>;
 
  using OpMixDivULhs = typename B::template OpMixDivTimesVec<DIM>;
  using OpMixDivUCylLhs =
      typename B::template OpMixDivTimesVec<DIM, CYLINDRICAL>;
  using OpLambdaGraULhs = typename B::template OpMixTensorTimesGrad<DIM>;
 
  using OpMixDivULhsSide = OpMixLhsSide<OpMixDivULhs>;
  using OpMixDivUCylLhsSide = OpMixLhsSide<OpMixDivUCylLhs>;
  using OpLambdaGraULhsSide = OpMixLhsSide<OpLambdaGraULhs>;
 
  auto unity = []() { return 1; };
  auto jacobian = [is_axisymmetric](const double r, const double,
                                    const double) {
    if (is_axisymmetric)
      return 2. * M_PI * r;
    else
      return 1.;
  };
 
  if (!is_axisymmetric) {
    op_loop_side->getOpPtrVector().push_back(
        new OpMixDivULhsSide(sigma, u, unity, jacobian, true));
  } else {
    op_loop_side->getOpPtrVector().push_back(
        new OpMixDivUCylLhsSide(sigma, u, unity, jacobian, true));
  }
  op_loop_side->getOpPtrVector().push_back(
      new OpLambdaGraULhsSide(sigma, u, unity, jacobian, true));
 
  op_loop_side->getSideFEPtr()->getRuleHook = rule;
  MoFEMFunctionReturn(0);
}
 
template <int DIM, AssemblyType A, IntegrationType I, typename BoundaryEleOp>
MoFEMErrorCode opFactoryBoundaryLhs(
    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
    std::string sigma, std::string u, bool is_axisymmetric = false) {
  MoFEMFunctionBegin;
 
  using C = ContactIntegrators<BoundaryEleOp>;
 
  auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
 
  pip.push_back(new OpCalculateVectorFieldValues<DIM>(
      u, common_data_ptr->contactDispPtr()));
  pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
      sigma, common_data_ptr->contactTractionPtr()));
  pip.push_back(
      new typename C::template Assembly<A>::template OpConstrainBoundaryLhs_dU<
          DIM, GAUSS>(sigma, u, common_data_ptr, is_axisymmetric));
  pip.push_back(new typename C::template Assembly<A>::
                    template OpConstrainBoundaryLhs_dTraction<DIM, GAUSS>(
                        sigma, sigma, common_data_ptr, is_axisymmetric));
 
  MoFEMFunctionReturn(0);
}
 
template <int DIM, AssemblyType A, IntegrationType I, typename BoundaryEleOp>
MoFEMErrorCode opFactoryBoundaryRhs(
    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
    std::string sigma, std::string u, bool is_axisymmetric = false) {
  MoFEMFunctionBegin;
 
  using C = ContactIntegrators<BoundaryEleOp>;
 
  auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
 
  pip.push_back(new OpCalculateVectorFieldValues<DIM>(
      u, common_data_ptr->contactDispPtr()));
  pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
      sigma, common_data_ptr->contactTractionPtr()));
  pip.push_back(
      new typename C::template Assembly<A>::template OpConstrainBoundaryRhs<
          DIM, GAUSS>(sigma, common_data_ptr, is_axisymmetric));
 
  MoFEMFunctionReturn(0);
}
 
template <int DIM, IntegrationType I, typename BoundaryEleOp>
MoFEMErrorCode opFactoryCalculateTraction(
    boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
    std::string sigma, bool is_axisymmetric = false) {
  MoFEMFunctionBegin;
 
  using C = ContactIntegrators<BoundaryEleOp>;
 
  auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
  pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
      sigma, common_data_ptr->contactTractionPtr()));
  pip.push_back(new typename C::template OpAssembleTotalContactTraction<DIM, I>(
      common_data_ptr, 1. / scale, is_axisymmetric));
 
  MoFEMFunctionReturn(0);
}
 
}; // namespace ContactOps
 
#endif // __CONTACTOPS_HPP__