MortarContactProblem.cpp

Go to the documentation of this file.

/* This file is part of MoFEM.
 * MoFEM is free software: you can redistribute it and/or modify it under
 * the terms of the GNU Lesser General Public License as published by the
 * Free Software Foundation, either version 3 of the License, or (at your
 * option) any later version.
 *
 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 * License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
 
MoFEMErrorCode
MortarContactProblem::MortarContactElement::setGaussPts(int order) {
  MoFEMFunctionBegin;
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  auto get_tensor_vec = [](VectorDouble &v) {
    return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&v(0), &v(1),
                                                              &v(2));
  };
 
  order *= 2; // multiply by 2 due to the integrand of NTN (twice the
              // approximation)
 
  // Defined integration points for only 1 integrated tris
  int nb_gauss_pts_1tri;
  MatrixDouble gaussPts_1tri;
  nb_gauss_pts_1tri = QUAD_2D_TABLE[order]->npoints;
 
  gaussPts_1tri.resize(3, nb_gauss_pts_1tri, false);
  cblas_dcopy(nb_gauss_pts_1tri, &QUAD_2D_TABLE[order]->points[1], 3,
              &gaussPts_1tri(0, 0), 1);
  cblas_dcopy(nb_gauss_pts_1tri, &QUAD_2D_TABLE[order]->points[2], 3,
              &gaussPts_1tri(1, 0), 1);
  cblas_dcopy(nb_gauss_pts_1tri, QUAD_2D_TABLE[order]->weights, 1,
              &gaussPts_1tri(2, 0), 1);
  dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts_1tri, 3,
                                                         false);
  double *shape_ptr =
      &*dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
  cblas_dcopy(3 * nb_gauss_pts_1tri, QUAD_2D_TABLE[order]->points, 1, shape_ptr,
              1);
 
  // get the entity (prism in this case)
  EntityHandle prism_entity =
      numeredEntFiniteElementPtr->getEnt(); // we will get prism element here
  typedef ContactSearchKdTree::ContactCommonData_multiIndex::index<
      ContactSearchKdTree::Prism_tag>::type::iterator ItMultIndexPrism;
  ItMultIndexPrism it_mult_index_prism =
      contactCommondataMultiIndex->get<ContactSearchKdTree::Prism_tag>().find(
          prism_entity);
 
  // get integration tris of each prism
  Range range_poly_tris;
  range_poly_tris.clear();
  range_poly_tris = it_mult_index_prism->get()->commonIntegratedTriangle;
 
  // Number of integration points = number of Gauss points in each common
  // tris * no of common tris * 2 (2 here is due to different integratio
  // rule
  // of master and slave tris)
  int nb_gauss_pts = nb_gauss_pts_1tri * range_poly_tris.size();
  gaussPtsMaster.resize(3, nb_gauss_pts, false);
  gaussPtsSlave.resize(3, nb_gauss_pts, false);
 
  const EntityHandle *conn_slave = NULL;
  int num_nodes_prism = 0;
  rval = mField.get_moab().get_connectivity(prism_entity, conn_slave,
                                            num_nodes_prism);
  CHKERRQ_MOAB(rval);
  // for prism element bottom 3 nodes belong to slave tri and top 3 nodes
  // belong to master tri
  VectorDouble coords_prism; // 6*3=18 coordinates for 6 nodes in 3D
  coords_prism.resize(18, false);
  rval = mField.get_moab().get_coords(conn_slave, num_nodes_prism,
                                      &*coords_prism.data().begin());
  CHKERRQ_MOAB(rval);
 
  VectorDouble v_elem_coords_master, v_elem_coords_slave;
  v_elem_coords_master.resize(9, false);
  v_elem_coords_master = subrange(coords_prism, 0, 9);
  v_elem_coords_slave.resize(9, false);
  v_elem_coords_slave = subrange(coords_prism, 9, 18);
  auto t_elem_coords_master = get_tensor_vec(v_elem_coords_master);
  auto t_elem_coords_slave = get_tensor_vec(v_elem_coords_slave);
  // Here we need to calculate the local ara of integration triangles in
  // both master and slave surfaces (we need this for weight calculation of
  // the Gauss points)
 
  // To do this, first convert the master and slave triangles to 2D case
  // (or z=0)  as both master and slave triangles are oriented in 3D space
 
  // tansfermation(rotation) matrix (use the same matrix as defined in
  // ContactSearchKdTree)
  MatrixDouble m_rot;
  m_rot.resize(3, 3, false);
 
  FTensor::Tensor2<double *, 3, 3> t_m_rot{
      &m_rot(0, 0), &m_rot(0, 1), &m_rot(0, 2), &m_rot(1, 0), &m_rot(1, 1),
      &m_rot(1, 2), &m_rot(2, 0), &m_rot(2, 1), &m_rot(2, 2)};
 
  ContactSearchKdTree contact_problem_mulit_index(mField);
 
  CHKERR contact_problem_mulit_index.rotationMatrix(m_rot, v_elem_coords_slave);
 
  VectorDouble v_elem_coords_master_new, v_elem_coords_slave_new;
  v_elem_coords_master_new.resize(9, false);
  v_elem_coords_slave_new.resize(9, false);
  auto t_elem_coords_master_new = get_tensor_vec(v_elem_coords_master_new);
  auto t_elem_coords_slave_new = get_tensor_vec(v_elem_coords_slave_new);
 
  // master and slave tri elemnets (coord in 2D)
  VectorDouble v_elem_coords_master_2d, v_elem_coords_slave_2d;
  v_elem_coords_master_2d.resize(6, false);
  v_elem_coords_slave_2d.resize(6, false);
  FTensor::Tensor1<double *, 2> t_elem_coords_master_2d{
      &v_elem_coords_master_2d(0), &v_elem_coords_master_2d(1), 2};
 
  FTensor::Tensor1<double *, 2> t_elem_coords_slave_2d{
      &v_elem_coords_slave_2d(0), &v_elem_coords_slave_2d(1), 2};
  // transfer the slave and master tri coordinates (n=3) one by one to z
  // plane
  for (int ii = 0; ii != 3; ++ii) {
    t_elem_coords_slave_new(i) = t_m_rot(i, j) * t_elem_coords_slave(j);
    t_elem_coords_master_new(i) = t_m_rot(i, j) * t_elem_coords_master(j);
 
    t_elem_coords_master_2d(0) = t_elem_coords_master_new(0);
    t_elem_coords_master_2d(1) = t_elem_coords_master_new(1);
 
    t_elem_coords_slave_2d(0) = t_elem_coords_slave_new(0);
    t_elem_coords_slave_2d(1) = t_elem_coords_slave_new(1);
 
    ++t_elem_coords_slave_new;
    ++t_elem_coords_slave;
    ++t_elem_coords_master_new;
    ++t_elem_coords_master;
 
    ++t_elem_coords_master_2d;
    ++t_elem_coords_slave_2d;
  }
 
  // get nodal coordinates
  VectorDouble v_coords_integration_tri;
  v_coords_integration_tri.resize(9, false);
 
  // transfer coord to z plane
  VectorDouble v_coords_integration_tri_new;
  v_coords_integration_tri_new.resize(9, false);
 
  VectorDouble loc_coords_master; // starting point
  loc_coords_master.resize(2, false);
  double n_input[3] = {1, 0, 0}; // shape function at starting point
 
  // function to calculate the local coordinate of element based on its
  // global coordinates  master is the bottom surface
  VectorDouble loc_coords_slave;
  loc_coords_slave.resize(2, false);
 
  // global coordinate of each Gauss point in 2D (z plane)
  // x and y global coordinates
  VectorDouble v_glob_coords;
  v_glob_coords.resize(2, false);
 
  auto t_gaussPtsMaster = getFTensor1FromMat<3>(gaussPtsMaster);
  auto t_gaussPtsSlave = getFTensor1FromMat<3>(gaussPtsSlave);
  // for each prism loop over all the integration tris
  for (Range::iterator it_tri = range_poly_tris.begin();
       it_tri != range_poly_tris.end(); ++it_tri) {
    const EntityHandle *conn_face;
    int num_nodes_tri;
    // get nodes attached to the tri
    rval = mField.get_moab().get_connectivity(*it_tri, conn_face, num_nodes_tri,
                                              true);
    CHKERRQ_MOAB(rval);
    // v_coords_integration_tri.clear(); //[x1 y1 z1 x2 y2 z2 .......]
    rval = mField.get_moab().get_coords(
        conn_face, num_nodes_tri, &*v_coords_integration_tri.data().begin());
    CHKERRQ_MOAB(rval);
 
    auto t_coords_integration_tri = get_tensor_vec(v_coords_integration_tri);
    auto t_coords_integration_tri_new =
        get_tensor_vec(v_coords_integration_tri_new);
 
    for (int ii = 0; ii != 3; ++ii) {
      t_coords_integration_tri_new(i) =
          t_m_rot(i, j) * t_coords_integration_tri(j);
      t_coords_integration_tri_new(2) = 0;
      ++t_coords_integration_tri;
      ++t_coords_integration_tri_new;
    }
 
    // shape function derivative for tri elements (these are constant)
    // diff_n_tri = [dN1/dxi, dN1/deta, dN2/dxi, dN2/deta, dN3/dxi,
    // dN3/deta] = [-1 -1 1 0 0 1]
    double diff_n_tri[6];
    CHKERR ShapeDiffMBTRI(diff_n_tri);
 
    // calculate local coordinates of each integration tri
    // double n_input[3] = {1, 0, 0}; // shape function at starting point
    // function to calculate the local coordinate of element based on its
    // global coordinates  element local coordinates of nodes of integration
    // tri
    // double coords_integration_tri_loc_master[9],
    //     coords_integration_tri_loc_slave[9];
 
    VectorDouble coords_integration_tri_loc_master,
        coords_integration_tri_loc_slave;
    coords_integration_tri_loc_master.resize(6, false);
    coords_integration_tri_loc_slave.resize(6, false);
    FTensor::Tensor1<double *, 2> t_coords_integration_tri_loc_master{
        &coords_integration_tri_loc_master(0),
        &coords_integration_tri_loc_master(1), 2};
    FTensor::Tensor1<double *, 2> t_coords_integration_tri_loc_slave{
        &coords_integration_tri_loc_slave(0),
        &coords_integration_tri_loc_slave(1), 2};
 
    auto t_coords_integration_tri_new_2 =
        get_tensor_vec(v_coords_integration_tri_new);
 
    for (int ii = 0; ii != 3; ++ii) {
 
      double glob_coords_tri[2];
      glob_coords_tri[0] = t_coords_integration_tri_new_2(0);
      glob_coords_tri[1] = t_coords_integration_tri_new_2(1);
 
      // local coordinates of integration tri in master element
      double loc_coords_tri_master[2] = {0, 0}; // starting point
      CHKERR ShapeMBTRI_inverse(n_input, diff_n_tri,
                                &v_elem_coords_master_2d(0), glob_coords_tri,
                                loc_coords_tri_master);
 
      t_coords_integration_tri_loc_master(0) = loc_coords_tri_master[0];
      t_coords_integration_tri_loc_master(1) = loc_coords_tri_master[1];
 
      // local coordinates of integration tri in slave element
      double loc_coords_tri_slave[2] = {0, 0}; // starting point
      CHKERR ShapeMBTRI_inverse(n_input, diff_n_tri, &v_elem_coords_slave_2d(0),
                                glob_coords_tri, loc_coords_tri_slave);
 
      t_coords_integration_tri_loc_slave(0) = loc_coords_tri_slave[0];
      t_coords_integration_tri_loc_slave(1) = loc_coords_tri_slave[1];
 
      ++t_coords_integration_tri_loc_master;
      ++t_coords_integration_tri_loc_slave;
      ++t_coords_integration_tri_new_2;
    }
 
    // local (not global) area or jacobian of integration tri  in both
    // master and slave triangles (can be +ve or -ve based on surface
    // orientation)
    double area_integration_tri_master_loc, area_integration_tri_slave_loc;
    area_integration_tri_master_loc =
        std::abs(area2D(&coords_integration_tri_loc_master(0),
                        &coords_integration_tri_loc_master(2),
                        &coords_integration_tri_loc_master(4)));
 
    area_integration_tri_slave_loc =
        std::abs(area2D(&coords_integration_tri_loc_slave(0),
                        &coords_integration_tri_loc_slave(2),
                        &coords_integration_tri_loc_slave(4)));
 
    auto t_gaussPts_1tri = getFTensor1FromMat<3>(gaussPts_1tri);
 
    // for each integration tri loop over all its Gauss points
    // calculate global coordinates of each integration point and then
    // calculate the local coordinates  of this integration point in each
    // master and slave surface
    for (int gg = 0; gg != nb_gauss_pts_1tri; ++gg) {
 
      t_gaussPtsMaster(2) =
          t_gaussPts_1tri(2) * area_integration_tri_master_loc *
          2; // 2 here is due to as for ref tri A=1/2 and w=1  or w=2*A
      t_gaussPtsSlave(2) =
          t_gaussPts_1tri(2) * area_integration_tri_slave_loc * 2;
 
      // shape function for each Gauss point
      MatrixDouble N_tri;
      N_tri.resize(1, 3, false);
 
      CHKERR ShapeMBTRI(&N_tri(0, 0), &t_gaussPts_1tri(0), &t_gaussPts_1tri(1),
                        1);
 
      FTensor::Tensor1<double *, 3> t_N_tri{&N_tri(0, 0), &N_tri(0, 1),
                                            &N_tri(0, 2)};
 
      FTensor::Tensor1<double *, 3> t_coords_integration_tri_new{
          &v_coords_integration_tri_new(0), &v_coords_integration_tri_new(3),
          &v_coords_integration_tri_new(6), 1};
 
      v_glob_coords(0) = t_N_tri(i) * t_coords_integration_tri_new(i);
      ++t_coords_integration_tri_new;
      v_glob_coords(1) = t_N_tri(i) * t_coords_integration_tri_new(i);
 
      // calculate local coordinates of each Gauss point in both slave
      // and master tris
      CHKERR ShapeMBTRI_inverse(
          n_input, diff_n_tri, &v_elem_coords_master_2d(0),
          &*v_glob_coords.data().begin(), &loc_coords_master(0));
 
      t_gaussPtsMaster(0) = loc_coords_master(0);
      t_gaussPtsMaster(1) = loc_coords_master(1);
 
      // slave is the top surface
      CHKERR ShapeMBTRI_inverse(n_input, diff_n_tri, &v_elem_coords_slave_2d(0),
                                &*v_glob_coords.data().begin(),
                                &loc_coords_slave(0));
 
      t_gaussPtsSlave(0) = loc_coords_slave[0];
      t_gaussPtsSlave(1) = loc_coords_slave[1];
 
      ++t_gaussPts_1tri;
      ++t_gaussPtsSlave;
      ++t_gaussPtsMaster;
    }
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
MortarContactProblem::MortarConvectMasterContactElement::setGaussPts(
    int order) {
  MoFEMFunctionBegin;
  CHKERR MortarContactElement::setGaussPts(order);
  CHKERR convectPtr->convectSlaveIntegrationPts<true>();
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
MortarContactProblem::MortarConvectSlaveContactElement::setGaussPts(int order) {
  MoFEMFunctionBegin;
  CHKERR MortarContactElement::setGaussPts(order);
  CHKERR convectPtr->convectSlaveIntegrationPts<false>();
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode MortarContactProblem::setContactOperatorsRhs(
    boost::shared_ptr<MortarContactElement> fe_rhs_mortar_contact,
    boost::shared_ptr<CommonDataMortarContact> common_data_mortar_contact,
    string field_name, string lagrange_field_name, bool is_alm,
    bool is_eigen_pos_field, string eigen_pos_field_name,
    bool use_reference_coordinates) {
  MoFEMFunctionBegin;
 
  fe_rhs_mortar_contact->getOpPtrVector().push_back(new OpGetNormalSlaveALE(
      "MESH_NODE_POSITIONS", common_data_mortar_contact));
 
  fe_rhs_mortar_contact->getOpPtrVector().push_back(new OpGetNormalMasterALE(
      "MESH_NODE_POSITIONS", common_data_mortar_contact));
 
  fe_rhs_mortar_contact->getOpPtrVector().push_back(
      new OpGetPositionAtGaussPtsMaster(field_name,
                                        common_data_mortar_contact));
 
  fe_rhs_mortar_contact->getOpPtrVector().push_back(
      new OpGetPositionAtGaussPtsSlave(field_name, common_data_mortar_contact));
 
  if (is_eigen_pos_field) {
    fe_rhs_mortar_contact->getOpPtrVector().push_back(
        new OpGetDeformationFieldForDisplAtGaussPtsMaster(
            eigen_pos_field_name, common_data_mortar_contact));
 
    fe_rhs_mortar_contact->getOpPtrVector().push_back(
        new OpGetMatPosForDisplAtGaussPtsMaster("MESH_NODE_POSITIONS",
                                                common_data_mortar_contact));
 
    fe_rhs_mortar_contact->getOpPtrVector().push_back(
        new OpGetMatPosForDisplAtGaussPtsSlave("MESH_NODE_POSITIONS",
                                               common_data_mortar_contact));
 
    fe_rhs_mortar_contact->getOpPtrVector().push_back(
        new OpGetDeformationFieldForDisplAtGaussPtsSlave(
            eigen_pos_field_name, common_data_mortar_contact));
  }
 
  if (use_reference_coordinates) {
    fe_rhs_mortar_contact->getOpPtrVector().push_back(
        new OpGetDeformationFieldForDisplAtGaussPtsMaster(
            "MESH_NODE_POSITIONS", common_data_mortar_contact));
 
    fe_rhs_mortar_contact->getOpPtrVector().push_back(
        new OpGetDeformationFieldForDisplAtGaussPtsSlave(
            "MESH_NODE_POSITIONS", common_data_mortar_contact));
  }
 
  fe_rhs_mortar_contact->getOpPtrVector().push_back(
      new OpGetGapSlave(field_name, common_data_mortar_contact));
 
  fe_rhs_mortar_contact->getOpPtrVector().push_back(
      new OpGetLagMulAtGaussPtsSlave(lagrange_field_name,
                                     common_data_mortar_contact));
 
  fe_rhs_mortar_contact->getOpPtrVector().push_back(new OpGetGaussPtsState(
      lagrange_field_name, common_data_mortar_contact, cnValue, is_alm));
 
  if (!is_alm) {
    fe_rhs_mortar_contact->getOpPtrVector().push_back(
        new OpCalContactTractionOnSlave(field_name,
                                        common_data_mortar_contact));
 
    fe_rhs_mortar_contact->getOpPtrVector().push_back(new OpCalIntCompFunSlave(
        lagrange_field_name, common_data_mortar_contact, cnValuePtr));
  } else {
    fe_rhs_mortar_contact->getOpPtrVector().push_back(
        new OpGetAugmentedLambdaSlave(field_name, common_data_mortar_contact,
                                      cnValue));
 
    fe_rhs_mortar_contact->getOpPtrVector().push_back(
        new OpCalAugmentedTractionRhsSlave(field_name,
                                           common_data_mortar_contact));
 
    fe_rhs_mortar_contact->getOpPtrVector().push_back(
        new OpGapConstraintAugmentedRhs(lagrange_field_name,
                                        common_data_mortar_contact, cnValue));
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode MortarContactProblem::setMasterForceOperatorsRhs(
    boost::shared_ptr<MortarContactElement> fe_rhs_mortar_contact,
    boost::shared_ptr<CommonDataMortarContact> common_data_mortar_contact,
    string field_name, string lagrange_field_name, bool is_alm,
    bool is_eigen_pos_field, string eigen_pos_field_name,
    bool use_reference_coordinates) {
  MoFEMFunctionBegin;
 
  fe_rhs_mortar_contact->getOpPtrVector().push_back(new OpGetNormalSlaveALE(
      "MESH_NODE_POSITIONS", common_data_mortar_contact));
 
  fe_rhs_mortar_contact->getOpPtrVector().push_back(new OpGetNormalMasterALE(
      "MESH_NODE_POSITIONS", common_data_mortar_contact));
 
  fe_rhs_mortar_contact->getOpPtrVector().push_back(
      new OpGetLagMulAtGaussPtsSlave(lagrange_field_name,
                                     common_data_mortar_contact));
 
  //   fe_rhs_mortar_contact->getOpPtrVector().push_back(
  //       new OpGetPositionAtGaussPtsMaster(field_name,
  //       common_data_mortar_contact));
 
  //   fe_rhs_mortar_contact->getOpPtrVector().push_back(
  //       new OpGetPositionAtGaussPtsSlave(field_name,
  //       common_data_mortar_contact));
 
  if (!is_alm) {
    fe_rhs_mortar_contact->getOpPtrVector().push_back(
        new OpCalContactTractionOnMaster(field_name,
                                         common_data_mortar_contact));
  } else {
 
    fe_rhs_mortar_contact->getOpPtrVector().push_back(
        new OpGetPositionAtGaussPtsMaster(field_name,
                                          common_data_mortar_contact));
 
    fe_rhs_mortar_contact->getOpPtrVector().push_back(
        new OpGetPositionAtGaussPtsSlave(field_name,
                                         common_data_mortar_contact));
 
    if (is_eigen_pos_field) {
      fe_rhs_mortar_contact->getOpPtrVector().push_back(
          new OpGetDeformationFieldForDisplAtGaussPtsMaster(
              eigen_pos_field_name, common_data_mortar_contact));
 
      fe_rhs_mortar_contact->getOpPtrVector().push_back(
          new OpGetMatPosForDisplAtGaussPtsMaster("MESH_NODE_POSITIONS",
                                                  common_data_mortar_contact));
 
      fe_rhs_mortar_contact->getOpPtrVector().push_back(
          new OpGetMatPosForDisplAtGaussPtsSlave("MESH_NODE_POSITIONS",
                                                 common_data_mortar_contact));
 
      fe_rhs_mortar_contact->getOpPtrVector().push_back(
          new OpGetDeformationFieldForDisplAtGaussPtsSlave(
              eigen_pos_field_name, common_data_mortar_contact));
    }
 
    if (use_reference_coordinates) {
      fe_rhs_mortar_contact->getOpPtrVector().push_back(
          new OpGetDeformationFieldForDisplAtGaussPtsMaster(
              "MESH_NODE_POSITIONS", common_data_mortar_contact));
 
      fe_rhs_mortar_contact->getOpPtrVector().push_back(
          new OpGetDeformationFieldForDisplAtGaussPtsSlave(
              "MESH_NODE_POSITIONS", common_data_mortar_contact));
    }
 
    fe_rhs_mortar_contact->getOpPtrVector().push_back(
        new OpGetGapSlave(field_name, common_data_mortar_contact));
 
    fe_rhs_mortar_contact->getOpPtrVector().push_back(
        new OpGetAugmentedLambdaSlave(field_name, common_data_mortar_contact,
                                      cnValue));
 
    fe_rhs_mortar_contact->getOpPtrVector().push_back(
        new OpCalAugmentedTractionRhsMaster(field_name,
                                            common_data_mortar_contact));
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode MortarContactProblem::setContactOperatorsLhs(
    boost::shared_ptr<MortarContactElement> fe_lhs_mortar_contact,
    boost::shared_ptr<CommonDataMortarContact> common_data_mortar_contact,
    string field_name, string lagrange_field_name, bool is_alm,
    bool is_eigen_pos_field, string eigen_pos_field_name,
    bool use_reference_coordinates) {
  MoFEMFunctionBegin;
 
  fe_lhs_mortar_contact->getOpPtrVector().push_back(
      new OpGetPositionAtGaussPtsMaster(field_name,
                                        common_data_mortar_contact));
 
  fe_lhs_mortar_contact->getOpPtrVector().push_back(new OpGetNormalSlaveALE(
      "MESH_NODE_POSITIONS", common_data_mortar_contact));
 
  fe_lhs_mortar_contact->getOpPtrVector().push_back(new OpGetNormalMasterALE(
      "MESH_NODE_POSITIONS", common_data_mortar_contact));
 
  fe_lhs_mortar_contact->getOpPtrVector().push_back(
      new OpGetPositionAtGaussPtsSlave(field_name, common_data_mortar_contact));
 
  if (is_eigen_pos_field) {
    fe_lhs_mortar_contact->getOpPtrVector().push_back(
        new OpGetDeformationFieldForDisplAtGaussPtsMaster(
            eigen_pos_field_name, common_data_mortar_contact));
 
    fe_lhs_mortar_contact->getOpPtrVector().push_back(
        new OpGetMatPosForDisplAtGaussPtsMaster("MESH_NODE_POSITIONS",
                                                common_data_mortar_contact));
 
    fe_lhs_mortar_contact->getOpPtrVector().push_back(
        new OpGetMatPosForDisplAtGaussPtsSlave("MESH_NODE_POSITIONS",
                                               common_data_mortar_contact));
 
    fe_lhs_mortar_contact->getOpPtrVector().push_back(
        new OpGetDeformationFieldForDisplAtGaussPtsSlave(
            eigen_pos_field_name, common_data_mortar_contact));
  }
 
  if (use_reference_coordinates) {
    fe_lhs_mortar_contact->getOpPtrVector().push_back(
        new OpGetDeformationFieldForDisplAtGaussPtsMaster(
            "MESH_NODE_POSITIONS", common_data_mortar_contact));
 
    fe_lhs_mortar_contact->getOpPtrVector().push_back(
        new OpGetDeformationFieldForDisplAtGaussPtsSlave(
            "MESH_NODE_POSITIONS", common_data_mortar_contact));
  }
 
  fe_lhs_mortar_contact->getOpPtrVector().push_back(
      new OpGetGapSlave(field_name, common_data_mortar_contact));
 
  fe_lhs_mortar_contact->getOpPtrVector().push_back(
      new OpGetLagMulAtGaussPtsSlave(lagrange_field_name,
                                     common_data_mortar_contact));
  if (!is_alm) {
    fe_lhs_mortar_contact->getOpPtrVector().push_back(
        new OpCalContactTractionOverLambdaSlaveSlave(
            field_name, lagrange_field_name, common_data_mortar_contact));
 
    fe_lhs_mortar_contact->getOpPtrVector().push_back(
        new OpCalDerIntCompFunOverLambdaSlaveSlave(
            lagrange_field_name, common_data_mortar_contact, cnValuePtr));
 
    fe_lhs_mortar_contact->getOpPtrVector().push_back(
        new OpCalDerIntCompFunOverSpatPosSlaveMaster(
            lagrange_field_name, field_name, common_data_mortar_contact,
            cnValuePtr));
 
    fe_lhs_mortar_contact->getOpPtrVector().push_back(
        new OpCalDerIntCompFunOverSpatPosSlaveSlave(
            lagrange_field_name, field_name, common_data_mortar_contact,
            cnValuePtr));
  } else {
    fe_lhs_mortar_contact->getOpPtrVector().push_back(
        new OpGetAugmentedLambdaSlave(field_name, common_data_mortar_contact,
                                      cnValue));
 
    fe_lhs_mortar_contact->getOpPtrVector().push_back(
        new OpCalContactAugmentedTractionOverLambdaSlaveSlave(
            field_name, lagrange_field_name, common_data_mortar_contact));
 
    fe_lhs_mortar_contact->getOpPtrVector().push_back(
        new OpCalContactAugmentedTractionOverSpatialSlaveSlave(
            field_name, field_name, cnValue, common_data_mortar_contact));
 
    fe_lhs_mortar_contact->getOpPtrVector().push_back(
        new OpCalContactAugmentedTractionOverSpatialSlaveMaster(
            field_name, field_name, cnValue, common_data_mortar_contact));
 
    fe_lhs_mortar_contact->getOpPtrVector().push_back(
        new OpGapConstraintAugmentedOverLambda(
            lagrange_field_name, common_data_mortar_contact, cnValue));
 
    fe_lhs_mortar_contact->getOpPtrVector().push_back(
        new OpGapConstraintAugmentedOverSpatialMaster(
            field_name, lagrange_field_name, common_data_mortar_contact,
            cnValue));
 
    fe_lhs_mortar_contact->getOpPtrVector().push_back(
        new OpGapConstraintAugmentedOverSpatialSlave(
            field_name, lagrange_field_name, common_data_mortar_contact,
            cnValue));
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode MortarContactProblem::setContactOperatorsLhs(
    boost::shared_ptr<MortarConvectMasterContactElement> fe_lhs_mortar_contact,
    boost::shared_ptr<CommonDataMortarContact> common_data_mortar_contact,
    string field_name, string lagrange_field_name, bool is_alm,
    bool is_eigen_pos_field, string eigen_pos_field_name,
    bool use_reference_coordinates) {
  MoFEMFunctionBegin;
 
  CHKERR MortarContactProblem::setContactOperatorsLhs(
      boost::dynamic_pointer_cast<MortarContactElement>(fe_lhs_mortar_contact),
      common_data_mortar_contact, field_name, lagrange_field_name, is_alm,
      is_eigen_pos_field, eigen_pos_field_name, use_reference_coordinates);
 
  //   CHKERR setContactOperatorsLhs(
  //       boost::dynamic_pointer_cast<MortarContactElement>(
  //           fe_lhs_mortar_contact),
  //       common_data_mortar_contact, field_name, lagrange_field_name);
 
  fe_lhs_mortar_contact->getOpPtrVector().push_back(
      new OpCalculateGradPositionXi(field_name, common_data_mortar_contact));
 
  fe_lhs_mortar_contact->getOpPtrVector().push_back(
      new OpLhsConvectIntegrationPtsConstrainMasterGap(
          lagrange_field_name, field_name, common_data_mortar_contact,
          cnValuePtr, ContactOp::FACESLAVESLAVE,
          fe_lhs_mortar_contact->getConvectPtr()->getDiffKsiSpatialSlave()));
 
  fe_lhs_mortar_contact->getOpPtrVector().push_back(
      new OpLhsConvectIntegrationPtsConstrainMasterGap(
          lagrange_field_name, field_name, common_data_mortar_contact,
          cnValuePtr, ContactOp::FACESLAVEMASTER,
          fe_lhs_mortar_contact->getConvectPtr()->getDiffKsiSpatialMaster()));
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode MortarContactProblem::setMasterForceOperatorsLhs(
    boost::shared_ptr<MortarContactElement> fe_lhs_mortar_contact,
    boost::shared_ptr<CommonDataMortarContact> common_data_mortar_contact,
    string field_name, string lagrange_field_name, bool is_alm,
    bool is_eigen_pos_field, string eigen_pos_field_name,
    bool use_reference_coordinates) {
  MoFEMFunctionBegin;
 
  fe_lhs_mortar_contact->getOpPtrVector().push_back(new OpGetNormalSlaveALE(
      "MESH_NODE_POSITIONS", common_data_mortar_contact));
 
  fe_lhs_mortar_contact->getOpPtrVector().push_back(new OpGetNormalMasterALE(
      "MESH_NODE_POSITIONS", common_data_mortar_contact));
 
  fe_lhs_mortar_contact->getOpPtrVector().push_back(
      new OpGetLagMulAtGaussPtsSlave(lagrange_field_name,
                                     common_data_mortar_contact));
  if (!is_alm) {
    fe_lhs_mortar_contact->getOpPtrVector().push_back(
        new OpCalContactTractionOverLambdaMasterSlave(
            field_name, lagrange_field_name, common_data_mortar_contact));
  } else {
 
    fe_lhs_mortar_contact->getOpPtrVector().push_back(
        new OpGetPositionAtGaussPtsMaster(field_name,
                                          common_data_mortar_contact));
 
    fe_lhs_mortar_contact->getOpPtrVector().push_back(
        new OpGetPositionAtGaussPtsSlave(field_name,
                                         common_data_mortar_contact));
 
    if (is_eigen_pos_field) {
      fe_lhs_mortar_contact->getOpPtrVector().push_back(
          new OpGetDeformationFieldForDisplAtGaussPtsMaster(
              eigen_pos_field_name, common_data_mortar_contact));
 
      fe_lhs_mortar_contact->getOpPtrVector().push_back(
          new OpGetMatPosForDisplAtGaussPtsMaster("MESH_NODE_POSITIONS",
                                                  common_data_mortar_contact));
 
      fe_lhs_mortar_contact->getOpPtrVector().push_back(
          new OpGetMatPosForDisplAtGaussPtsSlave("MESH_NODE_POSITIONS",
                                                 common_data_mortar_contact));
 
      fe_lhs_mortar_contact->getOpPtrVector().push_back(
          new OpGetDeformationFieldForDisplAtGaussPtsSlave(
              eigen_pos_field_name, common_data_mortar_contact));
    }
 
    if (use_reference_coordinates) {
      fe_lhs_mortar_contact->getOpPtrVector().push_back(
          new OpGetDeformationFieldForDisplAtGaussPtsMaster(
              "MESH_NODE_POSITIONS", common_data_mortar_contact));
 
      fe_lhs_mortar_contact->getOpPtrVector().push_back(
          new OpGetDeformationFieldForDisplAtGaussPtsSlave(
              "MESH_NODE_POSITIONS", common_data_mortar_contact));
    }
 
    fe_lhs_mortar_contact->getOpPtrVector().push_back(
        new OpGetGapSlave(field_name, common_data_mortar_contact));
 
    fe_lhs_mortar_contact->getOpPtrVector().push_back(
        new OpGetAugmentedLambdaSlave(field_name, common_data_mortar_contact,
                                      cnValue));
 
    fe_lhs_mortar_contact->getOpPtrVector().push_back(
        new OpCalContactAugmentedTractionOverLambdaMasterSlave(
            field_name, lagrange_field_name, common_data_mortar_contact));
 
    fe_lhs_mortar_contact->getOpPtrVector().push_back(
        new OpCalContactAugmentedTractionOverSpatialMasterMaster(
            field_name, field_name, cnValue, common_data_mortar_contact));
 
    fe_lhs_mortar_contact->getOpPtrVector().push_back(
        new OpCalContactAugmentedTractionOverSpatialMasterSlave(
            field_name, field_name, cnValue, common_data_mortar_contact));
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode MortarContactProblem::setMasterForceOperatorsLhs(
    boost::shared_ptr<MortarConvectSlaveContactElement> fe_lhs_mortar_contact,
    boost::shared_ptr<CommonDataMortarContact> common_data_mortar_contact,
    string field_name, string lagrange_field_name, bool is_alm,
    bool is_eigen_pos_field, string eigen_pos_field_name,
    bool use_reference_coordinates) {
  MoFEMFunctionBegin;
 
  CHKERR MortarContactProblem::setMasterForceOperatorsLhs(
      boost::dynamic_pointer_cast<MortarContactElement>(fe_lhs_mortar_contact),
      common_data_mortar_contact, field_name, lagrange_field_name, is_alm,
      is_eigen_pos_field, eigen_pos_field_name, use_reference_coordinates);
 
  fe_lhs_mortar_contact->getOpPtrVector().push_back(new OpCalculateGradLambdaXi(
      lagrange_field_name, common_data_mortar_contact));
 
  fe_lhs_mortar_contact->getOpPtrVector().push_back(
      new OpLhsConvectIntegrationPtsContactTraction(
          field_name, field_name, common_data_mortar_contact,
          ContactOp::FACEMASTERSLAVE,
          fe_lhs_mortar_contact->getConvectPtr()->getDiffKsiSpatialSlave()));
 
  fe_lhs_mortar_contact->getOpPtrVector().push_back(
      new OpLhsConvectIntegrationPtsContactTraction(
          field_name, field_name, common_data_mortar_contact,
          ContactOp::FACEMASTERMASTER,
          fe_lhs_mortar_contact->getConvectPtr()->getDiffKsiSpatialMaster()));
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode MortarContactProblem::setContactOperatorsForPostProc(
    boost::shared_ptr<MortarContactElement> fe_post_proc_mortar_contact,
    boost::shared_ptr<CommonDataMortarContact> common_data_mortar_contact,
    MoFEM::Interface &m_field, string field_name, string lagrange_field_name,
    moab::Interface &moab_out, bool alm_flag, bool is_eigen_pos_field,
    string eigen_pos_field_name, bool is_displacements,
    Range post_proc_surface, double post_proc_gap_tol) {
  MoFEMFunctionBegin;
 
  fe_post_proc_mortar_contact->getOpPtrVector().push_back(
      new OpGetNormalMasterALE("MESH_NODE_POSITIONS",
                               common_data_mortar_contact));
 
  fe_post_proc_mortar_contact->getOpPtrVector().push_back(
      new OpGetNormalSlaveALE("MESH_NODE_POSITIONS",
                              common_data_mortar_contact));
 
  fe_post_proc_mortar_contact->getOpPtrVector().push_back(
      new OpGetPositionAtGaussPtsMaster(field_name,
                                        common_data_mortar_contact));
 
  fe_post_proc_mortar_contact->getOpPtrVector().push_back(
      new OpGetPositionAtGaussPtsSlave(field_name, common_data_mortar_contact));
 
  if (is_eigen_pos_field) {
    fe_post_proc_mortar_contact->getOpPtrVector().push_back(
        new OpGetDeformationFieldForDisplAtGaussPtsMaster(
            eigen_pos_field_name, common_data_mortar_contact));
 
    fe_post_proc_mortar_contact->getOpPtrVector().push_back(
        new OpGetMatPosForDisplAtGaussPtsMaster("MESH_NODE_POSITIONS",
                                                common_data_mortar_contact));
 
    fe_post_proc_mortar_contact->getOpPtrVector().push_back(
        new OpGetMatPosForDisplAtGaussPtsSlave("MESH_NODE_POSITIONS",
                                               common_data_mortar_contact));
 
    fe_post_proc_mortar_contact->getOpPtrVector().push_back(
        new OpGetDeformationFieldForDisplAtGaussPtsSlave(
            eigen_pos_field_name, common_data_mortar_contact));
  }
 
  if (is_displacements && !is_eigen_pos_field) {
    fe_post_proc_mortar_contact->getOpPtrVector().push_back(
        new OpGetDeformationFieldForDisplAtGaussPtsMaster(
            "MESH_NODE_POSITIONS", common_data_mortar_contact));
 
    fe_post_proc_mortar_contact->getOpPtrVector().push_back(
        new OpGetDeformationFieldForDisplAtGaussPtsSlave(
            "MESH_NODE_POSITIONS", common_data_mortar_contact));
  }
  fe_post_proc_mortar_contact->getOpPtrVector().push_back(
      new OpGetGapSlave(field_name, common_data_mortar_contact));
 
  fe_post_proc_mortar_contact->getOpPtrVector().push_back(
      new OpGetLagMulAtGaussPtsSlave(lagrange_field_name,
                                     common_data_mortar_contact));
 
  fe_post_proc_mortar_contact->getOpPtrVector().push_back(
      new OpLagGapProdGaussPtsSlave(lagrange_field_name,
                                    common_data_mortar_contact));
 
  fe_post_proc_mortar_contact->getOpPtrVector().push_back(
      new OpGetGaussPtsState(lagrange_field_name, common_data_mortar_contact,
                             cnValue, alm_flag));
 
  fe_post_proc_mortar_contact->getOpPtrVector().push_back(
      new OpMakeVtkSlave(m_field, field_name, common_data_mortar_contact,
                         moab_out, NO_TAG, post_proc_surface));
 
  fe_post_proc_mortar_contact->getOpPtrVector().push_back(new OpGetContactArea(
      lagrange_field_name, common_data_mortar_contact, cnValue, alm_flag,
      post_proc_surface, post_proc_gap_tol));
 
  MoFEMFunctionReturn(0);
}
 
// setup operators for calculation of active set
MoFEMErrorCode MortarContactProblem::setContactOperatorsRhsALEMaterial(
    boost::shared_ptr<MortarContactElement> fe_rhs_mortar_contact_ale,
    boost::shared_ptr<CommonDataMortarContact> common_data_mortar_contact,
    const string field_name, const string mesh_node_field_name,
    const string lagrange_field_name, const string side_fe_name) {
  MoFEMFunctionBegin;
 
  boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnContactPrismSide>
      fe_mat_side_rhs_master = boost::make_shared<
          VolumeElementForcesAndSourcesCoreOnContactPrismSide>(mField);
 
  fe_mat_side_rhs_master->getOpPtrVector().push_back(
      new OpCalculateVectorFieldGradient<3, 3>(
          mesh_node_field_name, common_data_mortar_contact->HMat));
  fe_mat_side_rhs_master->getOpPtrVector().push_back(
      new OpCalculateVectorFieldGradient<3, 3>(
          field_name, common_data_mortar_contact->hMat));
 
  boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnContactPrismSide>
      fe_mat_side_rhs_slave = boost::make_shared<
          VolumeElementForcesAndSourcesCoreOnContactPrismSide>(mField);
 
  fe_mat_side_rhs_slave->getOpPtrVector().push_back(
      new OpCalculateVectorFieldGradient<3, 3>(
          mesh_node_field_name, common_data_mortar_contact->HMat));
  fe_mat_side_rhs_slave->getOpPtrVector().push_back(
      new OpCalculateVectorFieldGradient<3, 3>(
          field_name, common_data_mortar_contact->hMat));
 
  fe_rhs_mortar_contact_ale->getOpPtrVector().push_back(new OpGetNormalSlaveALE(
      "MESH_NODE_POSITIONS", common_data_mortar_contact));
 
  fe_rhs_mortar_contact_ale->getOpPtrVector().push_back(
      new OpGetNormalMasterALE("MESH_NODE_POSITIONS",
                               common_data_mortar_contact));
 
  fe_rhs_mortar_contact_ale->getOpPtrVector().push_back(
      new OpGetPositionAtGaussPtsMaster(field_name,
                                        common_data_mortar_contact));
 
  fe_rhs_mortar_contact_ale->getOpPtrVector().push_back(
      new OpGetPositionAtGaussPtsSlave(field_name, common_data_mortar_contact));
 
  fe_rhs_mortar_contact_ale->getOpPtrVector().push_back(
      new OpGetGapSlave(field_name, common_data_mortar_contact));
 
  fe_rhs_mortar_contact_ale->getOpPtrVector().push_back(
      new OpGetLagMulAtGaussPtsSlave(lagrange_field_name,
                                     common_data_mortar_contact));
 
  fe_mat_side_rhs_master->getOpPtrVector().push_back(new OpCalculateDeformation(
      mesh_node_field_name, common_data_mortar_contact, false));
 
  fe_rhs_mortar_contact_ale->getOpPtrVector().push_back(
      new OpLoopForSideOfContactPrism(mesh_node_field_name,
                                      fe_mat_side_rhs_master, side_fe_name,
                                      ContactOp::FACEMASTER));
 
  fe_rhs_mortar_contact_ale->getOpPtrVector().push_back(
      new OpCalMatForcesALEMaster(mesh_node_field_name,
                                  common_data_mortar_contact));
 
  fe_mat_side_rhs_slave->getOpPtrVector().push_back(new OpCalculateDeformation(
      mesh_node_field_name, common_data_mortar_contact, false));
 
  fe_rhs_mortar_contact_ale->getOpPtrVector().push_back(
      new OpLoopForSideOfContactPrism(mesh_node_field_name,
                                      fe_mat_side_rhs_slave, side_fe_name,
                                      ContactOp::FACESLAVE));
 
  fe_rhs_mortar_contact_ale->getOpPtrVector().push_back(
      new OpCalMatForcesALESlave(mesh_node_field_name,
                                 common_data_mortar_contact));
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode MortarContactProblem::setContactOperatorsLhsALEMaterial(
    boost::shared_ptr<MortarContactElement> fe_lhs_mortar_contact_ale,
    boost::shared_ptr<CommonDataMortarContact> common_data_mortar_contact,
    const string field_name, const string mesh_node_field_name,
    const string lagrange_field_name, const string side_fe_name) {
  MoFEMFunctionBegin;
 
  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(new OpGetNormalSlaveALE(
      mesh_node_field_name, common_data_mortar_contact));
 
  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
      new OpGetNormalMasterALE(mesh_node_field_name,
                               common_data_mortar_contact));
 
  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
      new OpGetLagMulAtGaussPtsSlave(lagrange_field_name,
                                     common_data_mortar_contact));
 
  boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnContactPrismSide>
      feMatSideLhs_dx = boost::make_shared<
          VolumeElementForcesAndSourcesCoreOnContactPrismSide>(mField);
 
  feMatSideLhs_dx->getOpPtrVector().push_back(
      new OpCalculateVectorFieldGradient<3, 3>(
          mesh_node_field_name, common_data_mortar_contact->HMat));
 
  feMatSideLhs_dx->getOpPtrVector().push_back(
      new OpCalculateVectorFieldGradient<3, 3>(
          field_name, common_data_mortar_contact->hMat));
 
  //   // Master derivative over spatial
  feMatSideLhs_dx->getOpPtrVector().push_back(new OpCalculateDeformation(
      mesh_node_field_name, common_data_mortar_contact, false));
 
  feMatSideLhs_dx->getOpPtrVector().push_back(
      new OpContactMaterialVolOnSideLhs_dX_dX(
          mesh_node_field_name, mesh_node_field_name,
          common_data_mortar_contact, true));
 
  feMatSideLhs_dx->getOpPtrVector().push_back(
      new OpContactMaterialVolOnSideLhs_dX_dx(
          mesh_node_field_name, field_name, common_data_mortar_contact, true));
 
  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
      new OpLoopForSideOfContactPrism(mesh_node_field_name, feMatSideLhs_dx,
                                      side_fe_name, ContactOp::FACEMASTER));
 
  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
      new OpContactMaterialMasterOnFaceLhs_dX_dX(
          mesh_node_field_name, mesh_node_field_name,
          common_data_mortar_contact, POSITION_RANK, POSITION_RANK));
 
  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
      new OpContactMaterialMasterSlaveLhs_dX_dLagmult(
          mesh_node_field_name, lagrange_field_name, common_data_mortar_contact,
          POSITION_RANK, LAGRANGE_RANK));
 
  boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnContactPrismSide>
      feMatSideLhsSlave_dx = boost::make_shared<
          VolumeElementForcesAndSourcesCoreOnContactPrismSide>(mField);
 
  feMatSideLhsSlave_dx->getOpPtrVector().push_back(
      new OpCalculateVectorFieldGradient<3, 3>(
          mesh_node_field_name, common_data_mortar_contact->HMat));
 
  feMatSideLhsSlave_dx->getOpPtrVector().push_back(
      new OpCalculateVectorFieldGradient<3, 3>(
          field_name, common_data_mortar_contact->hMat));
 
  feMatSideLhsSlave_dx->getOpPtrVector().push_back(new OpCalculateDeformation(
      mesh_node_field_name, common_data_mortar_contact, false));
 
  feMatSideLhsSlave_dx->getOpPtrVector().push_back(
      new OpContactMaterialVolOnSideLhs_dX_dX(
          mesh_node_field_name, mesh_node_field_name,
          common_data_mortar_contact, false));
 
  feMatSideLhsSlave_dx->getOpPtrVector().push_back(
      new OpContactMaterialVolOnSideLhs_dX_dx(
          mesh_node_field_name, field_name, common_data_mortar_contact, false));
 
  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
      new OpLoopForSideOfContactPrism(mesh_node_field_name,
                                      feMatSideLhsSlave_dx, side_fe_name,
                                      ContactOp::FACESLAVE));
 
  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
      new OpContactMaterialSlaveOnFaceLhs_dX_dX(
          mesh_node_field_name, mesh_node_field_name,
          common_data_mortar_contact, POSITION_RANK, POSITION_RANK));
 
  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
      new OpContactMaterialSlaveSlaveLhs_dX_dLagmult(
          mesh_node_field_name, lagrange_field_name, common_data_mortar_contact,
          POSITION_RANK, LAGRANGE_RANK));
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode MortarContactProblem::setContactOperatorsLhsALE(
    boost::shared_ptr<MortarContactElement> fe_lhs_mortar_contact_ale,
    boost::shared_ptr<CommonDataMortarContact> common_data_mortar_contact,
    const string field_name, const string mesh_node_field_name,
    const string lagrange_field_name) {
  MoFEMFunctionBegin;
 
  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(new OpGetNormalSlaveALE(
      mesh_node_field_name, common_data_mortar_contact));
 
  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
      new OpGetNormalMasterALE(mesh_node_field_name,
                               common_data_mortar_contact));
 
  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
      new OpGetPositionAtGaussPtsMaster(field_name,
                                        common_data_mortar_contact));
 
  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
      new OpGetPositionAtGaussPtsSlave(field_name, common_data_mortar_contact));
 
  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
      new OpGetGapSlave(field_name, common_data_mortar_contact));
 
  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
      new OpGetLagMulAtGaussPtsSlave(lagrange_field_name,
                                     common_data_mortar_contact));
 
  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
      new OpContactTractionSlaveSlave_dX(field_name, mesh_node_field_name,
                                         common_data_mortar_contact,
                                         POSITION_RANK, POSITION_RANK));
 
  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
      new OpContactTractionMasterSlave_dX(field_name, mesh_node_field_name,
                                          common_data_mortar_contact,
                                          POSITION_RANK, POSITION_RANK));
 
  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
      new OpContactTractionMasterMaster_dX(field_name, mesh_node_field_name,
                                           common_data_mortar_contact,
                                           POSITION_RANK, POSITION_RANK));
 
  fe_lhs_mortar_contact_ale->getOpPtrVector().push_back(
      new OpCalDerIntCompFunSlaveSlave_dX(
          lagrange_field_name, mesh_node_field_name, cnValuePtr,
          common_data_mortar_contact, LAGRANGE_RANK, POSITION_RANK));
  MoFEMFunctionReturn(0);
}