mortar_contact_thermal.cpp

Implementation of mortar contact between surfaces with non-matching meshes taking into account internal stress resulting from the thermal expansion

/** \file mortar_contact_thermal.cpp
 * \example mortar_contact_thermal.cpp
 *
 * Implementation of mortar contact between surfaces with non-matching meshes
 * taking into account internal stress resulting from the thermal expansion
 *
 **/
 
/* 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/>.
 */
 
#include <Mortar.hpp>
#include <HookeInternalStressElement.hpp>
 
static char help[] = "\n";
double MortarContactProblem::LoadScale::lAmbda = 1;
using BlockData = NonlinearElasticElement::BlockData;
using VolSideFe = VolumeElementForcesAndSourcesCoreOnSide;
 
struct VolRule {
  int operator()(int, int, int order) const { return 2 * (order - 1); }
};
 
int main(int argc, char *argv[]) {
 
  const string default_options = "-ksp_type fgmres \n"
                                 "-pc_type lu \n"
                                 "-pc_factor_mat_solver_type mumps \n"
                                 "-mat_mumps_icntl_13 1 \n"
                                 "-mat_mumps_icntl_14 800 \n"
                                 "-mat_mumps_icntl_20 0 \n"
                                 "-mat_mumps_icntl_24 1 \n"
                                 "-snes_type newtonls \n"
                                 "-snes_linesearch_type basic \n"
                                 "-snes_divergence_tolerance 0 \n"
                                 "-snes_max_it 50 \n"
                                 "-snes_atol 1e-8 \n"
                                 "-snes_rtol 1e-10 \n"
                                 "-snes_monitor \n"
                                 "-ksp_monitor \n"
                                 "-snes_converged_reason \n"
                                 "-order 1 \n"
                                 "-order_lambda 1 \n"
                                 "-order_contact 2 \n"
                                 "-ho_levels_num 1 \n"
                                 "-step_num 1 \n"
                                 "-cn_value 1. \n"
                                 "-r_value 1. \n"
                                 "-alm_flag 0 \n"
                                 "-eigen_pos_flag 0 \n";
 
  string param_file = "param_file.petsc";
  if (!static_cast<bool>(ifstream(param_file))) {
    std::ofstream file(param_file.c_str(), std::ios::ate);
    if (file.is_open()) {
      file << default_options;
      file.close();
    }
  }
 
  // Initialize MoFEM
  MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
 
  // Create mesh database
  moab::Core mb_instance;              // create database
  moab::Interface &moab = mb_instance; // create interface to database
 
  try {
    PetscBool flg_file;
    PetscBool flg_file_out;
 
    char mesh_file_name[255];
    char output_mesh_name[255];
    PetscInt order = 1;
    PetscInt order_contact = 1;
    PetscInt nb_ho_levels = 0;
    PetscInt order_lambda = 1;
    PetscReal r_value = 1.;
    PetscReal cn_value = -1;
    PetscInt nb_sub_steps = 1;
    PetscBool is_partitioned = PETSC_FALSE;
    PetscBool is_newton_cotes = PETSC_FALSE;
    PetscInt test_num = 0;
    PetscBool convect_pts = PETSC_FALSE;
    PetscBool print_contact_state = PETSC_FALSE;
    PetscBool alm_flag = PETSC_FALSE;
    PetscBool eigen_pos_flag = PETSC_FALSE;
 
    PetscReal thermal_expansion_coef = 1e-5;
    PetscReal init_temp = 250.0;
    PetscReal scale_factor = 1.0;
    PetscBool analytical_input = PETSC_TRUE;
    char stress_tag_name[255];
    PetscBool flg_tag_name;
    PetscBool save_mean_stress = PETSC_TRUE;
    PetscBool ignore_contact = PETSC_FALSE;
    PetscBool ignore_pressure = PETSC_FALSE;
 
    CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Elastic Config", "none");
 
    CHKERR PetscOptionsString("-file_name", "mesh file name", "", "mesh.h5m",
                              mesh_file_name, 255, &flg_file);
    CHKERR PetscOptionsString("-output_mesh_file", "output mesh file name", "",
                              "mesh.h5m", output_mesh_name, 255, &flg_file_out);
 
    CHKERR PetscOptionsInt("-order", "approximation order of spatial positions",
                           "", 1, &order, PETSC_NULL);
    CHKERR PetscOptionsInt(
        "-order_contact",
        "approximation order of spatial positions in contact interface", "", 1,
        &order_contact, PETSC_NULL);
    CHKERR PetscOptionsInt("-ho_levels_num", "number of higher order levels",
                           "", 0, &nb_ho_levels, PETSC_NULL);
    CHKERR PetscOptionsInt("-order_lambda",
                           "approximation order of Lagrange multipliers", "", 1,
                           &order_lambda, PETSC_NULL);
 
    CHKERR PetscOptionsInt("-step_num", "number of steps", "", nb_sub_steps,
                           &nb_sub_steps, PETSC_NULL);
 
    CHKERR PetscOptionsBool("-is_partitioned",
                            "set if mesh is partitioned (this result that each "
                            "process keeps only part of the mes",
                            "", PETSC_FALSE, &is_partitioned, PETSC_NULL);
 
    CHKERR PetscOptionsReal("-cn_value", "default regularisation cn value", "",
                            1., &cn_value, PETSC_NULL);
 
    CHKERR PetscOptionsBool("-is_newton_cotes",
                            "set if Newton-Cotes quadrature rules are used", "",
                            PETSC_FALSE, &is_newton_cotes, PETSC_NULL);
    CHKERR PetscOptionsInt("-test_num", "test number", "", 0, &test_num,
                           PETSC_NULL);
    CHKERR PetscOptionsBool("-convect", "set to convect integration pts", "",
                            PETSC_FALSE, &convect_pts, PETSC_NULL);
    CHKERR PetscOptionsBool("-print_contact_state",
                            "output number of active gp at every iteration", "",
                            PETSC_FALSE, &print_contact_state, PETSC_NULL);
    CHKERR PetscOptionsBool("-alm_flag",
                            "if set use ALM, if not use C-function", "",
                            PETSC_FALSE, &alm_flag, PETSC_NULL);
    CHKERR PetscOptionsBool("-eigen_pos_flag",
                            "if set use eigen spatial positions are taken into "
                            "account for predeformed configuration",
                            "", PETSC_FALSE, &eigen_pos_flag, PETSC_NULL);
 
    CHKERR PetscOptionsReal("-scale_factor", "scale factor", "", scale_factor,
                            &scale_factor, PETSC_NULL);
 
    CHKERR PetscOptionsBool("-ignore_contact", "if set true, ignore contact",
                            "", PETSC_FALSE, &ignore_contact, PETSC_NULL);
    CHKERR PetscOptionsBool("-ignore_pressure", "if set true, ignore pressure",
                            "", PETSC_FALSE, &ignore_pressure, PETSC_NULL);
    CHKERR PetscOptionsBool("-analytical_input",
                            "if set true, use analytical strain", "",
                            PETSC_TRUE, &analytical_input, PETSC_NULL);
    CHKERR PetscOptionsBool("-save_mean_stress",
                            "if set true, save mean stress", "", PETSC_TRUE,
                            &save_mean_stress, PETSC_NULL);
    CHKERR PetscOptionsString("-stress_tag_name", "stress tag name file name",
                              "", "INTERNAL_STRESS", stress_tag_name, 255,
                              &flg_tag_name);
    CHKERR PetscOptionsReal(
        "-thermal_expansion_coef", "thermal expansion coef ", "",
        thermal_expansion_coef, &thermal_expansion_coef, PETSC_NULL);
    CHKERR PetscOptionsReal("-init_temp", "init_temp ", "", init_temp,
                            &init_temp, PETSC_NULL);
 
    ierr = PetscOptionsEnd();
    CHKERRQ(ierr);
 
    // Check if mesh file was provided
    if (flg_file != PETSC_TRUE) {
      SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -file_name (MESH FILE NEEDED)");
    }
 
    if (is_partitioned) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
              "Partitioned mesh is not supported");
    }
 
    const char *option;
    option = "";
    CHKERR moab.load_file(mesh_file_name, 0, option);
    Version file_ver;
    CHKERR MoFEM::Interface::getFileVersion(moab, file_ver);
    MOFEM_LOG("WORLD", Sev::inform) << "File version " << file_ver.strVersion();
 
    // Create MoFEM database and link it to MoAB
    MoFEM::Core core(moab);
    MoFEM::Interface &m_field = core;
 
    std::vector<BitRefLevel> bit_levels;
    bit_levels.push_back(BitRefLevel().set(0));
    auto bit_ref_manager = m_field.getInterface<BitRefManager>();
    CHKERR bit_ref_manager->setBitRefLevelByDim(0, 3, bit_levels.back());
 
    auto add_prism_interface = [&](std::vector<BitRefLevel> &bit_levels) {
      MoFEMFunctionBegin;
      auto prism_interface = m_field.getInterface<PrismInterface>();
 
      for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, cit)) {
        if (cit->getName().compare(0, 11, "INT_CONTACT") == 0) {
          CHKERR PetscPrintf(PETSC_COMM_WORLD, "Insert %s (id: %d)\n",
                             cit->getName().c_str(), cit->getMeshsetId());
          EntityHandle cubit_meshset = cit->getMeshset();
 
          // get tet entities from back bit_level
          EntityHandle ref_level_meshset;
          CHKERR moab.create_meshset(MESHSET_SET, ref_level_meshset);
          CHKERR bit_ref_manager->getEntitiesByTypeAndRefLevel(
              bit_levels.back(), BitRefLevel().set(), MBTET, ref_level_meshset);
          CHKERR bit_ref_manager->getEntitiesByTypeAndRefLevel(
              bit_levels.back(), BitRefLevel().set(), MBPRISM,
              ref_level_meshset);
 
          // get faces and tets to split
          CHKERR prism_interface->getSides(cubit_meshset, bit_levels.back(),
                                           true, 0);
          // set new bit level
          bit_levels.push_back(BitRefLevel().set(bit_levels.size()));
          // split faces and tets
          CHKERR prism_interface->splitSides(ref_level_meshset,
                                             bit_levels.back(), cubit_meshset,
                                             true, true, 0);
          // clean meshsets
          CHKERR moab.delete_entities(&ref_level_meshset, 1);
 
          // update cubit meshsets
          for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, ciit)) {
            EntityHandle cubit_meshset = ciit->meshset;
            CHKERR bit_ref_manager->updateMeshsetByEntitiesChildren(
                cubit_meshset, bit_levels.back(), cubit_meshset, MBVERTEX,
                true);
            CHKERR bit_ref_manager->updateMeshsetByEntitiesChildren(
                cubit_meshset, bit_levels.back(), cubit_meshset, MBEDGE, true);
            CHKERR bit_ref_manager->updateMeshsetByEntitiesChildren(
                cubit_meshset, bit_levels.back(), cubit_meshset, MBTRI, true);
            CHKERR bit_ref_manager->updateMeshsetByEntitiesChildren(
                cubit_meshset, bit_levels.back(), cubit_meshset, MBTET, true);
          }
 
          CHKERR bit_ref_manager->shiftRightBitRef(1);
          bit_levels.pop_back();
        }
      }
 
      MoFEMFunctionReturn(0);
    };
 
    auto find_contact_prisms = [&](std::vector<BitRefLevel> &bit_levels,
                                   Range &simple_contact_prisms,
                                   Range &simple_master_tris,
                                   Range &simple_slave_tris) {
      MoFEMFunctionBegin;
 
      EntityHandle meshset_prisms;
      CHKERR moab.create_meshset(MESHSET_SET, meshset_prisms);
      CHKERR bit_ref_manager->getEntitiesByTypeAndRefLevel(
          bit_levels.back(), BitRefLevel().set(), MBPRISM, meshset_prisms);
      CHKERR moab.get_entities_by_handle(meshset_prisms, simple_contact_prisms);
      CHKERR moab.delete_entities(&meshset_prisms, 1);
 
      EntityHandle tri;
      for (Range::iterator pit = simple_contact_prisms.begin();
           pit != simple_contact_prisms.end(); pit++) {
        CHKERR moab.side_element(*pit, 2, 3, tri);
        simple_master_tris.insert(tri);
        CHKERR moab.side_element(*pit, 2, 4, tri);
        simple_slave_tris.insert(tri);
      }
 
      MoFEMFunctionReturn(0);
    };
 
    Range simple_contact_prisms, simple_master_tris, simple_slave_tris;
 
    if (!ignore_contact) {
      CHKERR add_prism_interface(bit_levels);
      CHKERR find_contact_prisms(bit_levels, simple_contact_prisms,
                                 simple_master_tris, simple_slave_tris);
    }
 
    Range mortar_contact_prisms, mortar_master_tris, mortar_slave_tris;
 
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
      if (it->getName().compare(0, 13, "MORTAR_MASTER") == 0) {
        CHKERR m_field.get_moab().get_entities_by_type(
            it->meshset, MBTRI, mortar_master_tris, true);
      }
    }
 
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
      if (it->getName().compare(0, 12, "MORTAR_SLAVE") == 0) {
        CHKERR m_field.get_moab().get_entities_by_type(it->meshset, MBTRI,
                                                       mortar_slave_tris, true);
      }
    }
 
    ContactSearchKdTree contact_search_kd_tree(m_field);
 
    boost::shared_ptr<ContactSearchKdTree::ContactCommonData_multiIndex>
        contact_commondata_multi_index;
 
    contact_commondata_multi_index =
        boost::shared_ptr<ContactSearchKdTree::ContactCommonData_multiIndex>(
            new ContactSearchKdTree::ContactCommonData_multiIndex());
 
    CHKERR contact_search_kd_tree.buildTree(mortar_master_tris);
 
    CHKERR contact_search_kd_tree.contactSearchAlgorithm(
        mortar_master_tris, mortar_slave_tris, contact_commondata_multi_index,
        mortar_contact_prisms);
 
    EntityHandle meshset_slave_master_prisms;
    CHKERR moab.create_meshset(MESHSET_SET, meshset_slave_master_prisms);
 
    CHKERR
    moab.add_entities(meshset_slave_master_prisms, mortar_contact_prisms);
 
    CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
        meshset_slave_master_prisms, 3, bit_levels.back());
 
    Range tris_from_prism;
    // find actual masters and slave used
    CHKERR m_field.get_moab().get_adjacencies(mortar_contact_prisms, 2, true,
                                              tris_from_prism,
                                              moab::Interface::UNION);
 
    tris_from_prism = tris_from_prism.subset_by_type(MBTRI);
    mortar_slave_tris = intersect(tris_from_prism, mortar_slave_tris);
 
    Range all_contact_prisms, all_slave_tris;
    all_contact_prisms = unite(simple_contact_prisms, mortar_contact_prisms);
    all_slave_tris = unite(simple_slave_tris, mortar_slave_tris);
 
    CHKERR m_field.add_field("SPATIAL_POSITION", H1, AINSWORTH_LEGENDRE_BASE, 3,
                             MB_TAG_SPARSE, MF_ZERO);
 
    CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
                             3, MB_TAG_SPARSE, MF_ZERO);
 
    if (!eigen_pos_flag) {
      CHKERR m_field.add_field("EIGEN_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
                               3, MB_TAG_SPARSE, MF_ZERO);
    }
 
    CHKERR m_field.add_field("LAGMULT", H1, AINSWORTH_LEGENDRE_BASE, 1,
                             MB_TAG_SPARSE, MF_ZERO);
 
    CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
    CHKERR m_field.set_field_order(0, MBTET, "MESH_NODE_POSITIONS", 1);
    CHKERR m_field.set_field_order(0, MBTRI, "MESH_NODE_POSITIONS", 1);
    CHKERR m_field.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 1);
    CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
 
    // Declare problem add entities (by tets) to the field
    CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "SPATIAL_POSITION");
    CHKERR m_field.set_field_order(0, MBTET, "SPATIAL_POSITION", order);
    CHKERR m_field.set_field_order(0, MBTRI, "SPATIAL_POSITION", order);
    CHKERR m_field.set_field_order(0, MBEDGE, "SPATIAL_POSITION", order);
    CHKERR m_field.set_field_order(0, MBVERTEX, "SPATIAL_POSITION", 1);
 
    if (!eigen_pos_flag) {
      CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "EIGEN_POSITIONS");
      CHKERR m_field.set_field_order(0, MBTET, "EIGEN_POSITIONS", order);
      CHKERR m_field.set_field_order(0, MBTRI, "EIGEN_POSITIONS", order);
      CHKERR m_field.set_field_order(0, MBEDGE, "EIGEN_POSITIONS", order);
      CHKERR m_field.set_field_order(0, MBVERTEX, "EIGEN_POSITIONS", 1);
    }
 
    CHKERR m_field.add_ents_to_field_by_type(all_slave_tris, MBTRI, "LAGMULT");
 
    CHKERR m_field.set_field_order(0, MBTRI, "LAGMULT", order_lambda);
    CHKERR m_field.set_field_order(0, MBEDGE, "LAGMULT", order_lambda);
    CHKERR m_field.set_field_order(0, MBVERTEX, "LAGMULT", 1);
 
    auto set_contact_order = [&](Range &contact_prisms, int order_contact,
                                 int nb_ho_levels) {
      MoFEMFunctionBegin;
      Range contact_tris, contact_edges;
      CHKERR moab.get_adjacencies(contact_prisms, 2, false, contact_tris,
                                  moab::Interface::UNION);
      contact_tris = contact_tris.subset_by_type(MBTRI);
      CHKERR moab.get_adjacencies(contact_tris, 1, false, contact_edges,
                                  moab::Interface::UNION);
      Range ho_ents;
      ho_ents.merge(contact_tris);
      ho_ents.merge(contact_edges);
      for (int ll = 0; ll < nb_ho_levels; ll++) {
        Range ents, verts, tets;
        CHKERR moab.get_connectivity(ho_ents, verts, true);
        CHKERR moab.get_adjacencies(verts, 3, false, tets,
                                    moab::Interface::UNION);
        tets = tets.subset_by_type(MBTET);
        for (auto d : {1, 2}) {
          CHKERR moab.get_adjacencies(tets, d, false, ents,
                                      moab::Interface::UNION);
        }
        ho_ents = unite(ho_ents, ents);
        ho_ents = unite(ho_ents, tets);
      }
 
      CHKERR m_field.set_field_order(ho_ents, "SPATIAL_POSITION",
                                     order_contact);
 
      MoFEMFunctionReturn(0);
    };
 
    if (!ignore_contact && order_contact > order) {
      CHKERR set_contact_order(all_contact_prisms, order_contact, nb_ho_levels);
    }
 
    // build field
    CHKERR m_field.build_fields();
 
    // Projection on "x" field
    {
      Projection10NodeCoordsOnField ent_method(m_field, "SPATIAL_POSITION");
      CHKERR m_field.loop_dofs("SPATIAL_POSITION", ent_method);
    }
    // Projection on "X" field
    {
      Projection10NodeCoordsOnField ent_method(m_field, "MESH_NODE_POSITIONS");
      CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method);
    }
 
    Range slave_verts;
    CHKERR moab.get_connectivity(all_slave_tris, slave_verts, true);
    CHKERR m_field.getInterface<FieldBlas>()->setField(0.0, MBVERTEX,
                                                       slave_verts, "LAGMULT");
 
    // Add elastic element
    boost::shared_ptr<std::map<int, BlockData>> block_sets_ptr =
        boost::make_shared<std::map<int, BlockData>>();
    CHKERR HookeElement::setBlocks(m_field, block_sets_ptr);
 
    boost::shared_ptr<ForcesAndSourcesCore> fe_elastic_lhs_ptr(
        new VolumeElementForcesAndSourcesCore(m_field));
    boost::shared_ptr<ForcesAndSourcesCore> fe_elastic_rhs_ptr(
        new VolumeElementForcesAndSourcesCore(m_field));
    fe_elastic_lhs_ptr->getRuleHook = VolRule();
    fe_elastic_rhs_ptr->getRuleHook = VolRule();
 
    CHKERR HookeElement::addElasticElement(m_field, block_sets_ptr, "ELASTIC",
                                           "SPATIAL_POSITION",
                                           "MESH_NODE_POSITIONS", false);
 
    auto data_hooke_element_at_pts =
        boost::make_shared<HookeInternalStressElement::DataAtIntegrationPts>();
    CHKERR HookeElement::setOperators(fe_elastic_lhs_ptr, fe_elastic_rhs_ptr,
                                      block_sets_ptr, "SPATIAL_POSITION",
                                      "MESH_NODE_POSITIONS", false, false,
                                      MBTET, data_hooke_element_at_pts);
    auto thermal_strain =
        [&thermal_expansion_coef, &init_temp, &test_num](
            FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> &t_coords) {
          FTensor::Tensor2_symmetric<double, 3> t_thermal_strain;
          FTensor::Index<'i', 3> i;
          FTensor::Index<'j', 3> j;
 
          constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
          double temp;
 
          t_thermal_strain(i, j) = 0.0;
 
          switch (test_num) {
          case 0:
            // Put here analytical formula which may depend on coordinates
            temp = init_temp - 1.0;
            t_thermal_strain(i, j) =
                thermal_expansion_coef * (temp - init_temp) * t_kd(i, j);
            break;
          case 1:
          case 2:
            t_thermal_strain(i, j) = -thermal_expansion_coef * t_kd(i, j);
            break;
          case 3:
            t_thermal_strain(2, 2) = thermal_expansion_coef;
            break;
          case 4:
            t_thermal_strain(i, j) = thermal_expansion_coef * t_kd(i, j);
            break;
          default:
            break;
          }
 
          return t_thermal_strain;
        };
 
    if (analytical_input) {
      fe_elastic_rhs_ptr->getOpPtrVector().push_back(
          new HookeElement::OpAnalyticalInternalStrain_dx<0>(
              "SPATIAL_POSITION", data_hooke_element_at_pts, thermal_strain));
      fe_elastic_rhs_ptr->getOpPtrVector().push_back(
          new HookeInternalStressElement::OpGetAnalyticalInternalStress<0>(
              "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts,
              thermal_strain));
    } else {
      fe_elastic_rhs_ptr->getOpPtrVector().push_back(
          new HookeInternalStressElement::OpGetInternalStress(
              "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts,
              moab, stress_tag_name));
      fe_elastic_rhs_ptr->getOpPtrVector().push_back(
          new HookeInternalStressElement::OpInternalStrain_dx(
              "SPATIAL_POSITION", data_hooke_element_at_pts));
    }
 
    fe_elastic_rhs_ptr->getOpPtrVector().push_back(
        new OpCalculateVectorFieldGradient<3, 3>(
            "SPATIAL_POSITION", data_hooke_element_at_pts->hMat));
    fe_elastic_rhs_ptr->getOpPtrVector().push_back(
        new OpCalculateVectorFieldGradient<3, 3>(
            "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
    fe_elastic_rhs_ptr->getOpPtrVector().push_back(
        new HookeInternalStressElement::OpSaveStress(
            "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts,
            *block_sets_ptr.get(), moab, scale_factor, save_mean_stress, false,
            false));
 
    Range all_tets;
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, bit)) {
      if (bit->getName().compare(0, 11, "MAT_ELASTIC") == 0) {
        Range tets;
        CHKERR moab.get_entities_by_dimension(bit->getMeshset(), 3, tets, true);
        all_tets.merge(tets);
      }
    }
    Skinner skinner(&moab);
    Range skin_tris;
    CHKERR skinner.find_skin(0, all_tets, false, skin_tris);
 
    CHKERR m_field.add_finite_element("SKIN", MF_ZERO);
    CHKERR m_field.add_ents_to_finite_element_by_type(skin_tris, MBTRI, "SKIN");
    CHKERR m_field.modify_finite_element_add_field_row("SKIN",
                                                       "SPATIAL_POSITION");
    CHKERR m_field.modify_finite_element_add_field_col("SKIN",
                                                       "SPATIAL_POSITION");
    CHKERR m_field.modify_finite_element_add_field_data("SKIN",
                                                        "SPATIAL_POSITION");
    CHKERR m_field.modify_finite_element_add_field_data("SKIN",
                                                        "MESH_NODE_POSITIONS");
    CHKERR m_field.modify_finite_element_add_field_data("SKIN",
                                                        "EIGEN_POSITIONS");
 
    auto make_mortar_contact_element = [&]() {
      return boost::make_shared<MortarContactProblem::MortarContactElement>(
          m_field, contact_commondata_multi_index, "SPATIAL_POSITION",
          "MESH_NODE_POSITIONS");
    };
 
    auto make_convective_mortar_master_element = [&]() {
      return boost::make_shared<
          MortarContactProblem::MortarConvectMasterContactElement>(
          m_field, contact_commondata_multi_index, "SPATIAL_POSITION",
          "MESH_NODE_POSITIONS");
    };
 
    auto make_convective_mortar_slave_element = [&]() {
      return boost::make_shared<
          MortarContactProblem::MortarConvectSlaveContactElement>(
          m_field, contact_commondata_multi_index, "SPATIAL_POSITION",
          "MESH_NODE_POSITIONS");
    };
 
    auto make_mortar_contact_common_data = [&]() {
      return boost::make_shared<MortarContactProblem::CommonDataMortarContact>(
          m_field);
    };
 
    auto get_mortar_contact_rhs = [&](auto contact_problem, auto make_element,
                                      bool is_alm = false) {
      auto fe_rhs_mortar_contact = make_element();
      auto common_data_mortar_contact = make_mortar_contact_common_data();
      if (print_contact_state) {
        fe_rhs_mortar_contact->contactStateVec =
            common_data_mortar_contact->gaussPtsStateVec;
      }
      contact_problem->setContactOperatorsRhs(
          fe_rhs_mortar_contact, common_data_mortar_contact, "SPATIAL_POSITION",
          "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
      return fe_rhs_mortar_contact;
    };
 
    auto get_mortar_master_traction_rhs = [&](auto contact_problem,
                                              auto make_element,
                                              bool is_alm = false) {
      auto fe_rhs_mortar_contact = make_element();
      auto common_data_mortar_contact = make_mortar_contact_common_data();
      contact_problem->setMasterForceOperatorsRhs(
          fe_rhs_mortar_contact, common_data_mortar_contact, "SPATIAL_POSITION",
          "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
      return fe_rhs_mortar_contact;
    };
 
    auto get_mortar_master_traction_lhs = [&](auto contact_problem,
                                              auto make_element,
                                              bool is_alm = false) {
      auto fe_lhs_mortar_contact = make_element();
      auto common_data_mortar_contact = make_mortar_contact_common_data();
      contact_problem->setMasterForceOperatorsLhs(
          fe_lhs_mortar_contact, common_data_mortar_contact, "SPATIAL_POSITION",
          "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
      return fe_lhs_mortar_contact;
    };
 
    auto get_mortar_master_help_traction_lhs =
        [&](auto contact_problem, auto make_element, bool is_alm = false) {
          auto fe_lhs_mortar_contact = make_element();
          auto common_data_mortar_contact = make_mortar_contact_common_data();
          contact_problem->setMasterForceOperatorsLhs(
              fe_lhs_mortar_contact, common_data_mortar_contact,
              "SPATIAL_POSITION", "LAGMULT", is_alm);
          return fe_lhs_mortar_contact;
        };
 
    auto get_mortar_contact_lhs = [&](auto contact_problem, auto make_element,
                                      bool is_alm = false) {
      auto fe_lhs_mortar_contact = make_element();
      auto common_data_mortar_contact = make_mortar_contact_common_data();
      contact_problem->setContactOperatorsLhs(
          fe_lhs_mortar_contact, common_data_mortar_contact, "SPATIAL_POSITION",
          "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
      return fe_lhs_mortar_contact;
    };
 
    auto get_mortar_contact_help_lhs =
        [&](auto contact_problem, auto make_element, bool is_alm = false) {
          auto fe_lhs_mortar_contact = make_element();
          auto common_data_mortar_contact = make_mortar_contact_common_data();
          contact_problem->setContactOperatorsLhs(
              fe_lhs_mortar_contact, common_data_mortar_contact,
              "SPATIAL_POSITION", "LAGMULT", is_alm);
          return fe_lhs_mortar_contact;
        };
 
    auto cn_value_ptr = boost::make_shared<double>(cn_value);
    auto mortar_contact_problem = boost::make_shared<MortarContactProblem>(
        m_field, contact_commondata_multi_index, cn_value_ptr, is_newton_cotes);
 
    // add fields to the global matrix by adding the element
 
    if (!eigen_pos_flag)
      mortar_contact_problem->addMortarContactElement(
          "MORTAR_CONTACT_ELEM", "SPATIAL_POSITION", "LAGMULT",
          mortar_contact_prisms);
    else
      mortar_contact_problem->addMortarContactElement(
          "MORTAR_CONTACT_ELEM", "SPATIAL_POSITION", "LAGMULT",
          mortar_contact_prisms, "MESH_NODE_POSITIONS", eigen_pos_flag,
          "EIGEN_POSITIONS");
 
    auto make_simple_contact_element = [&]() {
      return boost::make_shared<SimpleContactProblem::SimpleContactElement>(
          m_field);
    };
 
    auto make_convective_simple_master_element = [&]() {
      return boost::make_shared<
          SimpleContactProblem::ConvectMasterContactElement>(
          m_field, "SPATIAL_POSITION", "MESH_NODE_POSITIONS");
    };
 
    auto make_convective_simple_slave_element = [&]() {
      return boost::make_shared<
          SimpleContactProblem::ConvectSlaveContactElement>(
          m_field, "SPATIAL_POSITION", "MESH_NODE_POSITIONS");
    };
 
    auto make_simple_contact_common_data = [&]() {
      return boost::make_shared<SimpleContactProblem::CommonDataSimpleContact>(
          m_field);
    };
 
    auto get_simple_contact_rhs = [&](auto contact_problem, auto make_element,
                                      bool is_alm = false) {
      auto fe_rhs_simple_contact = make_element();
      auto common_data_simple_contact = make_simple_contact_common_data();
      if (print_contact_state) {
        fe_rhs_simple_contact->contactStateVec =
            common_data_simple_contact->gaussPtsStateVec;
      }
      contact_problem->setContactOperatorsRhs(
          fe_rhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
          "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
      return fe_rhs_simple_contact;
    };
 
    auto get_simple_master_traction_rhs = [&](auto contact_problem,
                                              auto make_element,
                                              bool is_alm = false) {
      auto fe_rhs_simple_contact = make_element();
      auto common_data_simple_contact = make_simple_contact_common_data();
      contact_problem->setMasterForceOperatorsRhs(
          fe_rhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
          "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
      return fe_rhs_simple_contact;
    };
 
    auto get_simple_master_traction_lhs =
        [&](auto contact_problem, auto make_element, bool is_alm = false) {
          auto fe_lhs_simple_contact = make_element();
          auto common_data_simple_contact = make_simple_contact_common_data();
          contact_problem->setMasterForceOperatorsLhs(
              fe_lhs_simple_contact, common_data_simple_contact,
              "SPATIAL_POSITION", "LAGMULT", is_alm);
          return fe_lhs_simple_contact;
        };
 
    auto get_contact_lhs = [&](auto contact_problem, auto make_element,
                               bool is_alm = false) {
      auto fe_lhs_simple_contact = make_element();
      auto common_data_simple_contact = make_simple_contact_common_data();
      contact_problem->setContactOperatorsLhs(
          fe_lhs_simple_contact, common_data_simple_contact, "SPATIAL_POSITION",
          "LAGMULT", is_alm, eigen_pos_flag, "EIGEN_POSITIONS");
      return fe_lhs_simple_contact;
    };
 
    auto simple_contact_problem = boost::make_shared<SimpleContactProblem>(
        m_field, cn_value_ptr, is_newton_cotes);
    simple_contact_problem->addContactElement("SIMPLE_CONTACT_ELEM",
                                              "SPATIAL_POSITION", "LAGMULT",
                                              simple_contact_prisms);
 
    simple_contact_problem->addPostProcContactElement(
        "CONTACT_POST_PROC", "SPATIAL_POSITION", "LAGMULT",
        "MESH_NODE_POSITIONS", all_slave_tris);
 
    CHKERR MetaNeumannForces::addNeumannBCElements(m_field, "SPATIAL_POSITION");
 
    // Add spring boundary condition applied on surfaces.
    CHKERR MetaSpringBC::addSpringElements(m_field, "SPATIAL_POSITION",
                                           "MESH_NODE_POSITIONS");
 
    // build finite elemnts
    CHKERR m_field.build_finite_elements();
 
    // build adjacencies
    CHKERR m_field.build_adjacencies(bit_levels.back());
 
    // define problems
    CHKERR m_field.add_problem("CONTACT_PROB", MF_ZERO);
 
    // set refinement level for problem
    CHKERR m_field.modify_problem_ref_level_add_bit("CONTACT_PROB",
                                                    bit_levels.back());
 
    DMType dm_name = "DMMOFEM";
    CHKERR DMRegister_MoFEM(dm_name);
 
    SmartPetscObj<DM> dm;
    dm = createSmartDM(m_field.get_comm(), dm_name);
 
    // create dm instance
    CHKERR DMSetType(dm, dm_name);
 
    // set dm datastruture which created mofem datastructures
    CHKERR DMMoFEMCreateMoFEM(dm, &m_field, "CONTACT_PROB", bit_levels.back());
    CHKERR DMSetFromOptions(dm);
    CHKERR DMMoFEMSetIsPartitioned(dm, is_partitioned);
    // add elements to dm
    CHKERR DMMoFEMAddElement(dm, "MORTAR_CONTACT_ELEM");
    CHKERR DMMoFEMAddElement(dm, "SIMPLE_CONTACT_ELEM");
    CHKERR DMMoFEMAddElement(dm, "ELASTIC");
    CHKERR DMMoFEMAddElement(dm, "PRESSURE_FE");
    CHKERR DMMoFEMAddElement(dm, "SPRING");
    CHKERR DMMoFEMAddElement(dm, "CONTACT_POST_PROC");
    CHKERR DMMoFEMAddElement(dm, "SKIN");
 
    CHKERR DMSetUp(dm);
 
    // Vector of DOFs and the RHS
    auto D = smartCreateDMVector(dm);
    auto F = smartVectorDuplicate(D);
 
    // Stiffness matrix
    auto Aij = smartCreateDMMatrix(dm);
 
    CHKERR VecZeroEntries(D);
    CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
 
    CHKERR VecZeroEntries(F);
    CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
 
    CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
    CHKERR MatZeroEntries(Aij);
 
    fe_elastic_rhs_ptr->snes_f = F;
    fe_elastic_lhs_ptr->snes_B = Aij;
 
    // Dirichlet BC
    boost::shared_ptr<FEMethod> dirichlet_bc_ptr =
        boost::shared_ptr<FEMethod>(new DirichletSpatialPositionsBc(
            m_field, "SPATIAL_POSITION", Aij, D, F));
 
    dirichlet_bc_ptr->snes_ctx = SnesMethod::CTX_SNESNONE;
    dirichlet_bc_ptr->snes_x = D;
 
    // Assemble pressure and traction forces
    boost::ptr_map<std::string, NeumannForcesSurface> neumann_forces;
    CHKERR MetaNeumannForces::setMomentumFluxOperators(
        m_field, neumann_forces, NULL, "SPATIAL_POSITION");
 
    boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
        neumann_forces.begin();
    for (; mit != neumann_forces.end(); mit++) {
      mit->second->methodsOp.push_back(new MortarContactProblem::LoadScale());
      CHKERR DMMoFEMSNESSetFunction(dm, mit->first.c_str(),
                                    &mit->second->getLoopFe(), NULL, NULL);
    }
 
    // Implementation of spring element
    // Create new instances of face elements for springs
    boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr(
        new FaceElementForcesAndSourcesCore(m_field));
    boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
        new FaceElementForcesAndSourcesCore(m_field));
 
    CHKERR MetaSpringBC::setSpringOperators(
        m_field, fe_spring_lhs_ptr, fe_spring_rhs_ptr, "SPATIAL_POSITION",
        "MESH_NODE_POSITIONS");
 
    CHKERR DMoFEMPreProcessFiniteElements(dm, dirichlet_bc_ptr.get());
    CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
    CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, NULL,
                                  dirichlet_bc_ptr.get(), NULL);
    if (convect_pts == PETSC_TRUE) {
      CHKERR DMMoFEMSNESSetFunction(
          dm, "MORTAR_CONTACT_ELEM",
          get_mortar_contact_rhs(mortar_contact_problem,
                                 make_convective_mortar_master_element),
          PETSC_NULL, PETSC_NULL);
      CHKERR DMMoFEMSNESSetFunction(
          dm, "MORTAR_CONTACT_ELEM",
          get_mortar_master_traction_rhs(mortar_contact_problem,
                                         make_convective_mortar_slave_element),
          PETSC_NULL, PETSC_NULL);
    } else {
      CHKERR DMMoFEMSNESSetFunction(
          dm, "MORTAR_CONTACT_ELEM",
          get_mortar_contact_rhs(mortar_contact_problem,
                                 make_mortar_contact_element, alm_flag),
          PETSC_NULL, PETSC_NULL);
      CHKERR DMMoFEMSNESSetFunction(
          dm, "MORTAR_CONTACT_ELEM",
          get_mortar_master_traction_rhs(mortar_contact_problem,
                                         make_mortar_contact_element, alm_flag),
          PETSC_NULL, PETSC_NULL);
    }
 
    if (convect_pts == PETSC_TRUE) {
      CHKERR DMMoFEMSNESSetFunction(
          dm, "SIMPLE_CONTACT_ELEM",
          get_simple_contact_rhs(simple_contact_problem,
                                 make_convective_simple_master_element),
          PETSC_NULL, PETSC_NULL);
      CHKERR DMMoFEMSNESSetFunction(
          dm, "SIMPLE_CONTACT_ELEM",
          get_simple_master_traction_rhs(simple_contact_problem,
                                         make_convective_simple_slave_element),
          PETSC_NULL, PETSC_NULL);
    } else {
      CHKERR DMMoFEMSNESSetFunction(
          dm, "SIMPLE_CONTACT_ELEM",
          get_simple_contact_rhs(simple_contact_problem,
                                 make_simple_contact_element, alm_flag),
          PETSC_NULL, PETSC_NULL);
      CHKERR DMMoFEMSNESSetFunction(
          dm, "SIMPLE_CONTACT_ELEM",
          get_simple_master_traction_rhs(simple_contact_problem,
                                         make_simple_contact_element, alm_flag),
          PETSC_NULL, PETSC_NULL);
    }
 
    CHKERR DMMoFEMSNESSetFunction(dm, "ELASTIC", fe_elastic_rhs_ptr, PETSC_NULL,
                                  PETSC_NULL);
 
    CHKERR DMMoFEMSNESSetFunction(dm, "SPRING", fe_spring_rhs_ptr, PETSC_NULL,
                                  PETSC_NULL);
    CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, PETSC_NULL, PETSC_NULL,
                                  dirichlet_bc_ptr.get());
 
    boost::shared_ptr<FEMethod> fe_null;
    CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, fe_null, dirichlet_bc_ptr,
                                  fe_null);
    if (convect_pts == PETSC_TRUE) {
      CHKERR DMMoFEMSNESSetJacobian(
          dm, "MORTAR_CONTACT_ELEM",
          get_mortar_contact_help_lhs(mortar_contact_problem,
                                      make_convective_mortar_master_element),
          PETSC_NULL, PETSC_NULL);
      CHKERR DMMoFEMSNESSetJacobian(
          dm, "MORTAR_CONTACT_ELEM",
          get_mortar_master_help_traction_lhs(
              mortar_contact_problem, make_convective_mortar_slave_element),
          PETSC_NULL, PETSC_NULL);
    } else {
      CHKERR DMMoFEMSNESSetJacobian(
          dm, "MORTAR_CONTACT_ELEM",
          get_mortar_contact_lhs(mortar_contact_problem,
                                 make_mortar_contact_element, alm_flag),
          PETSC_NULL, PETSC_NULL);
      CHKERR DMMoFEMSNESSetJacobian(
          dm, "MORTAR_CONTACT_ELEM",
          get_mortar_master_traction_lhs(mortar_contact_problem,
                                         make_mortar_contact_element, alm_flag),
          PETSC_NULL, PETSC_NULL);
    }
 
    if (convect_pts == PETSC_TRUE) {
      CHKERR DMMoFEMSNESSetJacobian(
          dm, "SIMPLE_CONTACT_ELEM",
          get_contact_lhs(simple_contact_problem,
                          make_convective_simple_master_element),
          PETSC_NULL, PETSC_NULL);
      CHKERR DMMoFEMSNESSetJacobian(
          dm, "SIMPLE_CONTACT_ELEM",
          get_simple_master_traction_lhs(simple_contact_problem,
                                         make_convective_simple_slave_element),
          PETSC_NULL, PETSC_NULL);
    } else {
      CHKERR DMMoFEMSNESSetJacobian(dm, "SIMPLE_CONTACT_ELEM",
                                    get_contact_lhs(simple_contact_problem,
                                                    make_simple_contact_element,
                                                    alm_flag),
                                    PETSC_NULL, PETSC_NULL);
      CHKERR DMMoFEMSNESSetJacobian(
          dm, "SIMPLE_CONTACT_ELEM",
          get_simple_master_traction_lhs(simple_contact_problem,
                                         make_simple_contact_element, alm_flag),
          PETSC_NULL, PETSC_NULL);
    }
    CHKERR DMMoFEMSNESSetJacobian(dm, "ELASTIC", fe_elastic_lhs_ptr, PETSC_NULL,
                                  PETSC_NULL);
    CHKERR DMMoFEMSNESSetJacobian(dm, "SPRING", fe_spring_lhs_ptr, PETSC_NULL,
                                  PETSC_NULL);
    CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, fe_null, fe_null,
                                  dirichlet_bc_ptr);
 
    if (test_num) {
      char testing_options[] = "-ksp_type fgmres "
                               "-pc_type lu "
                               "-pc_factor_mat_solver_type mumps "
                               "-snes_type newtonls "
                               "-snes_linesearch_type basic "
                               "-snes_max_it 10 "
                               "-snes_atol 1e-8 "
                               "-snes_rtol 1e-8 ";
      CHKERR PetscOptionsInsertString(PETSC_NULL, testing_options);
    }
 
    auto snes = MoFEM::createSNES(m_field.get_comm());
    CHKERR SNESSetDM(snes, dm);
    SnesCtx *snes_ctx;
    // create snes nonlinear solver
    {
      CHKERR SNESSetDM(snes, dm);
      CHKERR DMMoFEMGetSnesCtx(dm, &snes_ctx);
      CHKERR SNESSetFunction(snes, F, SnesRhs, snes_ctx);
      CHKERR SNESSetJacobian(snes, Aij, Aij, SnesMat, snes_ctx);
      CHKERR SNESSetFromOptions(snes);
    }
 
    /// Post proc on the skin
    PostProcFaceOnRefinedMesh post_proc_skin(m_field);
    CHKERR post_proc_skin.generateReferenceElementMesh();
    CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", post_proc_skin, false, false);
    CHKERR post_proc_skin.addFieldValuesPostProc("SPATIAL_POSITION");
    CHKERR post_proc_skin.addFieldValuesPostProc("MESH_NODE_POSITIONS");
    CHKERR post_proc_skin.addFieldValuesPostProc("EIGEN_POSITIONS");
 
    struct OpGetFieldGradientValuesOnSkin
        : public FaceElementForcesAndSourcesCore::UserDataOperator {
 
      const std::string feVolName;
      boost::shared_ptr<VolSideFe> sideOpFe;
 
      OpGetFieldGradientValuesOnSkin(const std::string field_name,
                                     const std::string vol_fe_name,
                                     boost::shared_ptr<VolSideFe> side_fe)
          : FaceElementForcesAndSourcesCore::UserDataOperator(
                field_name, UserDataOperator::OPROW),
            feVolName(vol_fe_name), sideOpFe(side_fe) {}
 
      MoFEMErrorCode doWork(int side, EntityType type,
                            EntitiesFieldData::EntData &data) {
        MoFEMFunctionBegin;
        if (type != MBVERTEX)
          MoFEMFunctionReturnHot(0);
        CHKERR loopSideVolumes(feVolName, *sideOpFe);
        MoFEMFunctionReturn(0);
      }
    };
 
    boost::shared_ptr<VolSideFe> vol_side_fe_ptr =
        boost::make_shared<VolSideFe>(m_field);
    vol_side_fe_ptr->getOpPtrVector().push_back(
        new OpCalculateVectorFieldGradient<3, 3>(
            "SPATIAL_POSITION", data_hooke_element_at_pts->hMat));
    vol_side_fe_ptr->getOpPtrVector().push_back(
        new OpCalculateVectorFieldGradient<3, 3>(
            "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
 
    post_proc_skin.getOpPtrVector().push_back(
        new OpGetFieldGradientValuesOnSkin("SPATIAL_POSITION", "ELASTIC",
                                           vol_side_fe_ptr));
    post_proc_skin.getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<3>(
            "SPATIAL_POSITION", data_hooke_element_at_pts->spatPosMat));
    post_proc_skin.getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<3>(
            "MESH_NODE_POSITIONS", data_hooke_element_at_pts->meshNodePosMat));
 
    post_proc_skin.getOpPtrVector().push_back(
        new HookeInternalStressElement::OpPostProcHookeElement<
            FaceElementForcesAndSourcesCore>(
            "SPATIAL_POSITION", data_hooke_element_at_pts,
            *block_sets_ptr.get(), post_proc_skin.postProcMesh,
            post_proc_skin.mapGaussPts, false, false));
 
    for (int ss = 0; ss != nb_sub_steps; ++ss) {
      if (!ignore_pressure) {
        MortarContactProblem::LoadScale::lAmbda = (ss + 1.0) / nb_sub_steps;
      } else {
        MortarContactProblem::LoadScale::lAmbda = 0;
        CHKERR PetscPrintf(PETSC_COMM_WORLD, "Ignoring pressure...\n");
      }
      CHKERR PetscPrintf(PETSC_COMM_WORLD, "Load scale: %6.4e\n",
                         MortarContactProblem::LoadScale::lAmbda);
 
      CHKERR SNESSolve(snes, PETSC_NULL, D);
 
      CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
    }
 
    // save on mesh
    CHKERR DMoFEMMeshToGlobalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
 
    SmartPetscObj<Vec> v_energy;
    CHKERR HookeElement::calculateEnergy(dm, block_sets_ptr, "SPATIAL_POSITION",
                                         "MESH_NODE_POSITIONS", false, false,
                                         v_energy);
    double energy;
    CHKERR VecSum(v_energy, &energy);
    // Print elastic energy
    MOFEM_LOG_C("WORLD", Sev::inform, "Elastic energy: %8.8e", energy);
 
    // moab_instance
    moab::Core mb_post;                   // create database
    moab::Interface &moab_proc = mb_post; // create interface to database
 
    auto common_data_mortar_contact = make_mortar_contact_common_data();
 
    boost::shared_ptr<MortarContactProblem::MortarContactElement>
        fe_post_proc_mortar_contact;
    if (convect_pts == PETSC_TRUE) {
      fe_post_proc_mortar_contact = make_convective_mortar_master_element();
    } else {
      fe_post_proc_mortar_contact = make_mortar_contact_element();
    }
 
    std::array<double, 2> nb_gauss_pts;
    std::array<double, 2> contact_area;
 
    auto get_contact_data = [&](auto vec, std::array<double, 2> &data) {
      MoFEMFunctionBegin;
      CHKERR VecAssemblyBegin(vec);
      CHKERR VecAssemblyEnd(vec);
      const double *array;
      CHKERR VecGetArrayRead(vec, &array);
      if (m_field.get_comm_rank() == 0) {
        for (int i : {0, 1})
          data[i] = array[i];
      }
      CHKERR VecRestoreArrayRead(vec, &array);
      MoFEMFunctionReturn(0);
    };
 
    mb_post.delete_mesh();
 
    if (!mortar_contact_prisms.empty()) {
      mortar_contact_problem->setContactOperatorsForPostProc(
          fe_post_proc_mortar_contact, common_data_mortar_contact, m_field,
          "SPATIAL_POSITION", "LAGMULT", mb_post, alm_flag, eigen_pos_flag,
          "EIGEN_POSITIONS");
 
      CHKERR VecZeroEntries(common_data_mortar_contact->gaussPtsStateVec);
      CHKERR VecZeroEntries(common_data_mortar_contact->contactAreaVec);
 
      CHKERR DMoFEMLoopFiniteElements(dm, "MORTAR_CONTACT_ELEM",
                                      fe_post_proc_mortar_contact);
 
      CHKERR get_contact_data(common_data_mortar_contact->gaussPtsStateVec,
                              nb_gauss_pts);
      CHKERR get_contact_data(common_data_mortar_contact->contactAreaVec,
                              contact_area);
 
      if (m_field.get_comm_rank() == 0) {
        PetscPrintf(PETSC_COMM_SELF,
                    "Active gauss pts (mortar): %d out of %d\n",
                    (int)nb_gauss_pts[0], (int)nb_gauss_pts[1]);
 
        PetscPrintf(
            PETSC_COMM_SELF,
            "Active contact area (mortar): %8.8f out of %8.8f (%8.8f% %)\n",
            contact_area[0], contact_area[1],
            contact_area[0] / contact_area[1] * 100.);
      }
    }
 
    auto common_data_simple_contact = make_simple_contact_common_data();
 
    boost::shared_ptr<SimpleContactProblem::SimpleContactElement>
        fe_post_proc_simple_contact;
    if (convect_pts == PETSC_TRUE) {
      fe_post_proc_simple_contact = make_convective_simple_master_element();
    } else {
      fe_post_proc_simple_contact = make_simple_contact_element();
    }
 
    if (!simple_contact_prisms.empty()) {
      simple_contact_problem->setContactOperatorsForPostProc(
          fe_post_proc_simple_contact, common_data_simple_contact, m_field,
          "SPATIAL_POSITION", "LAGMULT", mb_post, alm_flag, eigen_pos_flag,
          "EIGEN_POSITIONS");
 
      CHKERR VecZeroEntries(common_data_simple_contact->gaussPtsStateVec);
      CHKERR VecZeroEntries(common_data_simple_contact->contactAreaVec);
 
      CHKERR DMoFEMLoopFiniteElements(dm, "SIMPLE_CONTACT_ELEM",
                                      fe_post_proc_simple_contact);
 
      CHKERR get_contact_data(common_data_simple_contact->gaussPtsStateVec,
                              nb_gauss_pts);
      CHKERR get_contact_data(common_data_simple_contact->contactAreaVec,
                              contact_area);
 
      if (m_field.get_comm_rank() == 0) {
        MOFEM_LOG_C("SELF", Sev::inform,
                    "Active gauss pts (matching): %d out of %d",
                    (int)nb_gauss_pts[0], (int)nb_gauss_pts[1]);
        MOFEM_LOG_C(
            "SELF", Sev::inform,
            "Active contact area (matching): %8.8f out of %8.8f (%8.8f% %)",
            contact_area[0], contact_area[1],
            contact_area[0] / contact_area[1] * 100.);
      }
    }
 
    if (!ignore_contact) {
      string out_file_name;
      std::ostringstream strm;
      strm << "out_contact_integ_pts"
           << ".h5m";
      out_file_name = strm.str();
      MOFEM_LOG_C("WORLD", Sev::inform, "Write file %s\n",
                  out_file_name.c_str());
      CHKERR mb_post.write_file(out_file_name.c_str(), "MOAB",
                                "PARALLEL=WRITE_PART");
    }
 
    boost::shared_ptr<PostProcFaceOnRefinedMesh> post_proc_contact_ptr(
        new PostProcFaceOnRefinedMesh(m_field));
 
    CHKERR post_proc_contact_ptr->generateReferenceElementMesh();
    CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", *post_proc_contact_ptr, false,
                          false);
    CHKERR post_proc_contact_ptr->addFieldValuesPostProc("LAGMULT");
    CHKERR post_proc_contact_ptr->addFieldValuesPostProc("SPATIAL_POSITION");
    CHKERR post_proc_contact_ptr->addFieldValuesPostProc("MESH_NODE_POSITIONS");
 
    if (!all_slave_tris.empty()) {
      CHKERR DMoFEMLoopFiniteElements(dm, "CONTACT_POST_PROC",
                                      post_proc_contact_ptr);
      string out_file_name;
      std::ostringstream stm;
      stm << "out_contact_surface"
          << ".h5m";
      out_file_name = stm.str();
      MOFEM_LOG_C("WORLD", Sev::inform, "Write file %s",
                  out_file_name.c_str());
      CHKERR post_proc_contact_ptr->postProcMesh.write_file(
          out_file_name.c_str(), "MOAB", "PARALLEL=WRITE_PART");
    }
 
    CHKERR m_field.getInterface<FieldBlas>()->fieldCopy(1., "SPATIAL_POSITION",
                                                        "EIGEN_POSITIONS");
 
    {
      PetscPrintf(PETSC_COMM_WORLD, "Loop post proc on the skin\n");
      CHKERR DMoFEMLoopFiniteElements(dm, "SKIN", &post_proc_skin);
      ostringstream stm;
      string out_file_name;
      stm << "out_skin.h5m";
      out_file_name = stm.str();
      MOFEM_LOG_C("WORLD", Sev::inform, "Write file %s",
                  out_file_name.c_str());
      CHKERR post_proc_skin.writeFile(stm.str());
    }
 
    const int n_parts = m_field.get_comm_size();
    if (m_field.get_comm_rank() == 0) {
      CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(
          dm, "ELASTIC", fe_elastic_rhs_ptr, 0, n_parts);
 
      Range faces, tris, quads, tris_edges, quads_edges, ents_to_delete;
 
      CHKERR moab.get_adjacencies(all_contact_prisms, 2, true, faces,
                                  moab::Interface::UNION);
      tris = faces.subset_by_type(MBTRI);
      quads = faces.subset_by_type(MBQUAD);
      CHKERR moab.get_adjacencies(tris, 1, true, tris_edges,
                                  moab::Interface::UNION);
      CHKERR moab.get_adjacencies(quads, 1, true, quads_edges,
                                  moab::Interface::UNION);
 
      ents_to_delete.merge(all_contact_prisms);
      ents_to_delete.merge(quads);
      ents_to_delete.merge(subtract(quads_edges, tris_edges));
 
      CHKERR moab.delete_entities(ents_to_delete);
 
      MOFEM_LOG_C("WORLD", Sev::inform, "Write file %s", output_mesh_name);
      CHKERR moab.write_file(output_mesh_name, "MOAB");
 
      auto get_tag_handle = [&](auto name, auto size) {
        Tag th;
        std::vector<double> def_vals(size, 0.0);
        CHKERR moab.tag_get_handle(name, size, MB_TYPE_DOUBLE, th,
                                   MB_TAG_CREAT | MB_TAG_SPARSE,
                                   def_vals.data());
        return th;
      };
 
      if (test_num) {
        Range tets;
        CHKERR moab.get_entities_by_dimension(0, 3, tets);
        EntityHandle tet = tets.front();
        std::array<double, 9> internal_stress, actual_stress;
        std::array<double, 9> internal_stress_ref, actual_stress_ref;
        std::array<double, 2> nb_gauss_pts_ref, contact_area_ref;
        switch (test_num) {
        case 3:
          actual_stress_ref = {0., 0., -100., 0., 0., 0., 0., 0., 0.};
          if (strcmp(stress_tag_name, "INTERNAL_STRESS") == 0)
            internal_stress_ref = {0., 0., -200., 0., 0., 0., 0., 0., 0.};
          else
            internal_stress_ref = {0., 0., -100., 0., 0., 0., 0., 0., 0.};
          break;
        case 4:
          nb_gauss_pts_ref = {216, 492};
          contact_area_ref = {0.125, 0.25};
          break;
        default:
          SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "Test number %d not found",
                   test_num);
        }
 
        auto th_internal_stress = get_tag_handle("MED_INTERNAL_STRESS", 9);
        auto th_actual_stress = get_tag_handle("MED_ACTUAL_STRESS", 9);
        CHKERR moab.tag_get_data(th_internal_stress, &tet, 1,
                                 internal_stress.data());
        CHKERR moab.tag_get_data(th_actual_stress, &tet, 1,
                                 actual_stress.data());
        if (test_num == 4) {
          const double eps = 1e-3;
          if (std::abs(nb_gauss_pts_ref[0] - nb_gauss_pts[0]) > eps) {
            SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                     "Wrong number of active contact gauss pts: should be %d "
                     "but is %d",
                     (int)nb_gauss_pts_ref[0], (int)nb_gauss_pts[0]);
          }
          if (std::abs(nb_gauss_pts_ref[1] - nb_gauss_pts[1]) > eps) {
            SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                     "Wrong total number of contact gauss pts: should be %d "
                     "but is %d",
                     (int)nb_gauss_pts_ref[1], (int)nb_gauss_pts[1]);
          }
          if (std::abs(contact_area_ref[0] - contact_area[0]) > eps) {
            SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                     "Wrong active contact area: should be %g "
                     "but is %g",
                     contact_area_ref[0], contact_area[0]);
          }
          if (std::abs(contact_area_ref[1] - contact_area[1]) > eps) {
            SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                     "Wrong potential contact area: should be %g "
                     "but is %g",
                     contact_area_ref[1], contact_area[1]);
          }
        } else {
          if (save_mean_stress) {
            const double eps = 1e-10;
            for (int i = 0; i < 9; i++) {
              if (std::abs(internal_stress[i] - internal_stress_ref[i]) > eps) {
                SETERRQ3(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                         "Wrong component %d of internal stress: should be %g "
                         "but is %g",
                         i, internal_stress_ref[i], internal_stress[i]);
              }
              if (std::abs(actual_stress[i] - actual_stress_ref[i]) > eps) {
                SETERRQ3(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                         "Wrong component %d of actual stress: should be %g "
                         "but is %g",
                         i, actual_stress_ref[i], actual_stress[i]);
              }
            }
          }
        }
      }
    }
  }
  CATCH_ERRORS;
 
  // finish work cleaning memory, getting statistics, etc
  MoFEM::Core::Finalize();
 
  return 0;
}